Thesis Abstracts in Computational Applied Mathematics An Improved Mass Conserving Front Tracking Scheme by Kou-kung Alex Chang, Advisor: Brent Lindquist May, 1993 The front-tracking method of Glimm et al. is a general purpose hydrodynamics code for solving fluid flow problems, governed mainly by a hyperbolic system of equations, in which co-dimension one jump discontinuities in the solution are important. The codimension one grids track the position of the jump discontinuity surfaces in the solution. The original front-tracking algorithms related to moving the co-dimension one grids did not enforce the conservation of mass properties of the hyperbolic system on any grid of finite discretization size. In this thesis, five items are identified to contribute mass conservation errors, they are: (i) curve propagation, (ii) curve redistribution, (iii) curve tangling, (iv) barrier and domain boundary tangling, and (v) curve deletions. We have developed replacement algorithms for these five items, achieving substantially improved mass conservation properties. Domain Decomposition Approach to Mixed Finite Element Solution of Elliptic Problems by Dragan Mirkovic, Advisor: Brent Lindquist May, 1993 Efficient parallel algorithms for numerical solution of the second order elliptic problem are considered in which Raviart-Thomas mixed finite element spaces of the lowest order are used for a saddle-point variational formulation of the original problem. The resulting indefinite linear system is reformulated either as a symmetric or non-symmetric positive definite system. In the symmetric reformulation, preconditioned conjugate gradient method is used to solve linear system; in non-symmetric reformulation generalized conjugate-gradient- like methods are used. We construct efficient solution procedures for preconditioner evaluation, and give an analysis of the convergence. A parallel version of the algorithm was obtained using the nonoverlapping domain decomposition method proposed by Glowinski and Wheeler. For this method we propose an quasi-optimal preconditioner, for the subdomain interface problem. Numerical experiments from a parallel implementation on an Intel iPSC/860 hypercube computer are presented which show the method to be both efficient and scalable, The numerical examples also confirm the theoretical bounds of the condition number of the interface operator. Parallel Multigrid Algorithms for Inverse Problems of General 3D Elastic Wave Equation by Xiang Xiang Zhong, Advisor: Yung Ming Chen August, 1993 The parallel multigrid algorithm is a versatile and efficient iterative numerical algorithm for solving multi-parameter inverse problems of a system of nonlinear partial differential equations. Parallelism is introduced into the hierarchy of inversion algorithm to improve its efficiency for solving the multi-parameter inverse problems of general 3D elastic wave equation with a number of unknown parameters. A hierarchical multi-grid strategy is also incorporated into this parallel inversion algorithm in combination with the existing re-structuralization algorithm. Complexity analysis is carried out for parallel inversion algorithm, Implementation of this multigrid parallel algorithm on a network of workstations is quite suitable due to little inter- processor communication overhead of the algorithm. A speedup timing study is given on each stage of the algorithm. The algorithm has been tested on a cluster of 5 workstations for material parameter identification for the 3D elastic wave equation. As a result, the calculated parameters match the actual parameters within a percentage range with almost linear speedup when the number of processors matches that of Laplace parameters and load cases. Computer Determination of Interior Temperature Distribution from Surface Temperature Measurements in Microwave Processing of Materials by Xiaoan Chen, Advisor: Yung Ming Chen May, 1994 The Generalized Pulse-Spectrum Technique (GPST) is applied to the simulation of microwave processing of materials. Microwave heating has many advantages over conventional heating; however, there are problems. Among them is sample cracking caused by undesirable thermal gradients. Knowledge of temperature distribution is crucial in acquiring a comprehensive understanding of these problems. In recent years there has been an intensified effort in the analysis of microwave heating problems. Most of these analyses are based on the coupled electromagnetic and thermal system, which is very difficult to solve. In this dissertation, based only on the heat equation, the GPST inversion algorithm is used to reconstruct the interior transient temperature profiles from surface temperature measurements. A 2-D version of a computer code is developed, which is capable of determining not only the interior transient temperature profiles but also the internal equivalent heat source of an object under microwave heating. Numerical examples demonstrate that the GPST inversion algorithm is very effective in the analysis of microwave heating problems. A Numerical Investigation of the Richtmyer-Meshkov Instability Using Front Tracking by Richard Lansing Holmes, Advisor: John Grove August, 1994 The method of front tracking is used to simulate the shock-tube experiments of Meshkov and Benjamin to predict perturbation amplitude growth rates. The results of the simulations are in much better agreement with experiment than previous simulations or theory. This improved agreement is explained by showing that simplifications in the theories, such as assumptions of incompressibility or linearity, miss crucial aspects of the experiments and thus the theories give incorrect growth rates. In addition, it is shown that thus the theories give incorrect growth rates. In addition, it is shown that for simulations the correct time period for measurement of the instability is crucial in getting agreement with experiment. The long-time behavior of the growth rates is discussed and the simulations are compared to an incompressible potential flow model for bubble velocities due to Hecht, et al. It is determined that the model adequately describes late-time behavior if modifications are made to take into account early time compressibility effects. Finally, the simulations are validated through comparison with a linearized theory which is valid in the limit of small amplitudes. Mathematical Methods for Realistic Neuron Modeling by Steven Martin Rosenthal, Advisor: Reginald Tewarson August, 1994 This dissertation deals with the development and application of mathematical methods tailored to the problem of realistic neuron modeling. We have concentrated on modeling the movement of free calcium (Ca++). The intracellular concentration of second messengers such as calcium ([Ca++]) has a critical role in shaping neuronal information processing. The mechanisms that regulate the distribution of these substances include nonlinear processes with kinetics that span a wide range of spatial and temporal scales. Our models involved the numerical solution of nonlinear ordinary and partial differential, and algebraic equations. Several possible techniques to make this computationally more tractable are described. We develop algorithms for the numerical solution of large and sparse systems of equations which arise in realistic modeling applications. We describe how the structure of the system can be used to develop efficient algorithms. W e also make use of separate time steps for those parts of the system which change more "slowly" relative to the other parts of the system. Computational results from models employing these methods are analyzed. We provide some background on adaptive gridding techniques and examine the problem of maintaining conservation when regridding. We detail a method which we have found to be successful. Further experimental advances will require the use of models that take into account the discrete nature of the ionic channels and their nonuniform distribution. Towards this end, we describe how one can handle multi- dimensional diffusion. We discuss the numerical solution of a model where the channels are distributed in clusters. Two Phase Flow Analysis of Turbulent Mixing in the Rayleigh-Taylor Instability by Yupin Chen, Advisor: James Glimm December, 1994 We study the compressible two phase mixing problem by the Rayleigh-Taylor instability, an acceleration driven instability of an interface between two fluids having different densities. The computational data for the Rayleigh- Taylor unstable random interface mixing problem are derived from front tracking simulations. The dynamics of the interior portion of the mixing zone is simplified by the use of scaling variables. The size of the mixing zone suggests fixed-point behavior. The effect of compressibility on the mixing layer growth rate is also examined. It is found that the growth rates increase significantly with increasing compressibility. We compare two formulations of the statistical moments, one based on two phase flow and the other on turbulence models. These formulations describe different aspects of the mixing process. For the Rayleigh-Taylor instability problem considered here, the two phase flow moments appear to be preferable, in that they subsume the turbulence moments but not conversely. Effective dynamic equations are developed, in which no undetermined parameters are introduced in the interior of the mixing region. The data required to validate this model comes from numerical simulations of the two fluid Euler equations by the front tracking method. The equations change structure across a boundary separating one phase from two phase flow, and an additional closure condition is needed at this boundary. Front Tracking of Complex Wave Interactions by Brian Keith Boston, Advisor: John Grove August, 1995 The method of front tracking has been applied to several complex problems in gas dynamics. All interactions between tracked waves should be handled automatically and robustly, and this provides a challenge both analytically and from a programming standpoint. Such issues have formed the basis for the present research, and a description of the relevant algorithms and ideas is included in this thesis. Three applications are used to produce interesting interactions for study. The first is the bifurcation from regular to Mach reflection of a shock wave at a reflecting surface. The results indicate that, in an idealized limit neglecting many physical effects, front tracking prefers the mechanical equilibrium point as the transition criterion for pseudo-steady flows. The next application is the Richtmyer-Meshkov instability of an interface separating two gases. The refraction of a shock through the interface causes any perturbations to grow over time, becoming unstable and leading eventually to fingering and mixing of the gases. The focus has been on larger amplitude perturbations, leading to more complex and irregular refraction structures. The simulations show that the overall growth rate of the interface and mixing zone is continuous across the bifurcation to irregular configurations. The last application is the non-ideal blast phenomenon, which involves the interaction of an expanding blast wave with a layer of warmed air near a reflecting boundary. The thermal layer channels the energy and holds it close to the ground, leading to a cascading effect with many complex wave interactions. Numerous algorithms are discussed for dealing with such phenomena, including the automatic untracking of waves during a simulation. Since it is not possible to foresee every eventuality, the code should be robust enough to do something reasonable when it encounters an unexpected situation. Also needed is the ability to leave sufficiently weak scattered waves untracked during the installation of a complicated structure. All three applications also involved the implementation of overlapping subdomain boundary conditions, appropriate for periodic, reflecting and parallel processing artificial boundaries. Issues pertinent to the implementation of these algorithms are discussed. The Effects of Viscous Terms on Riemann Problem Solutions by Jane Marie Hurley, Advisor: Bradley Plohr August, 1995 In this thesis, we study systems of two conservation laws with quadratic flux functions. In particular, we are concerned with the Riemann problem, an initial-value problem where the initial data consist of two constant states on either side of a jump discontinuity. Because we allow certain discontinuous solutions (shock waves), we must impose an additional admissibility criterion in order to guarantee uniqueness. The criterion we use is the viscous profile criterion. The use of this criterion results in the presence of some nonclassical shock waves in the Riemann problem solutions. The most important of these non-classical shock waves are the transitional shock waves. The presence of transitional shock waves is highly dependent on the form of the viscosity matrix. Research to date using the viscous profile criterion has primarily set the viscosity matrix to be a multiple of the identity matrix. Though this choice allows some mathematical simplifications, it also leads to some non- generic behavior in the Riemann problem solutions. The goal of this thesis is to investigate how the structure of the solutions changes when the viscosity matrix is perturbed away from the identity. The results presented here are a combination of analytical and numerical work. We have worked in the framework of the fundamental wave manifold to derive a condition that, when satisfied, guarantees the presence of transitional shock waves. Using this condition, we are able to describe the exact regions in the wave manifold that correspond to transi-tional shock waves. Also, we determine the boundaries in the parameter plane that separate models with varying numbers of transitional regions. Finally, we use a computer program to generate the solution diagrams for several different models and viscosity matrices. Dispersion of Tracer Slugs for Flow in Porous Media by Ana Cristina da Silva Grossi, Advisor: James Glimm December, 1995 The dispersion of tracer slugs for flow in porous media is studied through numerical solutions of the governing system of partial differential equations. These solutions are obtained by a numerical method introduced in this thesis. This is a robust method based on the front tracking scheme, which main feature is the elimination of the numerical diffusion, and based on the analytic solution of the one-dimensional heat equation. The primary results of the analysis of tracer flows concern the determination of preferential flow directions in reservoirs and whether or not communication between wells exists. The method of tracer output interpretation discussed here may also be used to determine sweep efficiencies, reservoir permeabilities and permeability distribution. Geometric and Stochastic Analysis of Porous Media Structure by Sang-Moon Lee, Advisor: Brent Lindquist December, 1995 We investigate two techniques for studying stochastic properties of the micro-scale geometry of porous media from three dimensional images provided by X-ray computed microtomography. These two techniques are the 2-point correlation function and the medial axis transform. Analysis of the 2-point correlation function suggests a limited practicality of this application for measurements of porous media. Results are effectively restricted to computation of the porosity, the specific surface area (internal surface area per unit volume), and an estimate of the average grain size. In addition 2- point correlation functions are computationally expensive to achieve in three dimensions; practical speeds are obtained by using parallel architecture. The medial axis transform is considered a tool in the analysis of the geometric structure of void (pore) and grain space in porous media. From the analysis of the medial axis of the void space, the following useful measurements are obtained: a discrete version of the pore-size distribution, the distribution of disconnected void medial axis segments, tortuosity distributions between parallel faces, and the instantaneous maximum flow volume through the medial axis of the void structure. The same analysis is applied to the grain structure in the media. A new segmented algorithm is implemented that includes edge detection capabilities, such that for our analysis, we use the edge segmented data rather than simple threshold segmented data. Projective Block Lanczos Algorithm and Its Parallelization on Distributed- Memory Parallel Processors by Gen-Ching Lo, Advisor: Frank Webster (Chemistry) December, 1995 Lanczos methods are well-known algorithms for computing the eigenvalues and eigenvectors of large matrices. We develop in this thesis an improved version of the Lanczos algorithm which uses a minimal amount of memory and is tailored for parallel implementation. The development of our sequential projective block Lanczos algorithm is mainly guided by the following considerations. We should use operator algebra rather than matrix algebra in representing the eigensystem problem; the algorithm should be able to compute several hundred of the smallest magnitude eigenpairs; and Krylov vectors should not be stored in main memory. Our algorithmic optimization is constrained to large, non-sparse, Hermitian eigensystems where a significant fraction of the spectrum must be obtained reliably and completely. We have also designed a parallel version of the projective block Lanczos algorithm. Major operations of the sequential algorithm are identified and parallelized. Both the sequential and parallel algorithms have been implemented and were applied to solve large, non-sparse, Hermitian eigensystems where a significant fraction of the extreme eigenvalues and eigenvectors is sought. Numerical experiments are performed on clusters of magnetic dipoles, a quantum many-body system with direct coupling of all states and without pure spin states. Experimental results indicate that we have achieved both of the goals of better performance and being able to solve larger problems. Parallel Algorithms for Molecular Dynamics with Applications to Thin-Film Sputter Deposition by Roger Alan McCoy, Advisor: Yuefan Deng December, 1995 Material coatings of less than 1 micron in thickness are termed "thin-films". The miniaturization of microelectronic devices has created a need to understand atomic-scale behavior during thin film fabrication procedures such as sputter deposition. To study the dynamics of atomic systems and observe microstructure development during thin-film sputter deposition, an algorithm for parallel computers which performs large-scale, short-range molecular dynamics simulations has been developed. Our simulation package incorporates advanced computational techniques such as asynchronous communication, data caching, and a novel dynamic load balancing algorithm. We demonstrate the unique capability of efficient particle simulations in non-uniform, three- dimensional spatial regions on parallel computers. We have used our simulation package to obtain quantitative results regarding the sticking and penetration behavior of atoms deposited at low-energies for various incident angles. Shadowing and film-coverage behavior have been obtained which agree with laboratory observations. Also, a specialized localization procedure has been developed which substantially reduces the computational expense of deposition simulations. Analysis of the Equations of Conservative Elastoplasticity and Numerical Solution to the Riemann Problem by Steven Arno Tael, Advisor: John Grove December 1995 An analysis of the equations of conservative elastoplasticity, fixed point algorithms for solving the corresponding Riemann problem, and an application of Riemann solutions to front tracking is presented. The underlying system of equations, a novel strain-based system, is written as a first order hyperbolic form. An analysis of the smooth and discontinuous flow regimes is presented. It is shown that the familiar equations from gas dynamics can be derived from the system of equations for elastic flow. Furthermore. it is demonstrated how to differentiate the Cauchy stress tensor with respect to the underlying state variables. in a mathematically consistent manner. The equivalence of the jump relations for the continuity equation between the Lagrangian and Eulerian frames. for piecewise continuous flows, is also demonstrated. General constitutive modeling principles, and their application to modeling flow in metals (assuming a hyperelastic equation of state) are discussed. Details of calculating important tensorial derivatives of a specific model for the internal energy are presented. Details of several fixed point iterations for solving the Riemann problem are discussed. Comparison with the solution of the linearized pressure - shear dynamical equations for various sets of initial data are presented. In most cases, the agreement between the solution generated by the exact solver and an independently obtained solution is quite good. Furthermore, an application of the Riemann solver to front tracking is also presented. Numerical Study of Single-Electron Effects in Systems of Small Tunnel Junctions by Leonardo Fonseca Advisor: Konstantin Likharev (Physics) May, 1996 We describe a new and efficient method for the numerical study of the dynamics and statistics of single-electron systems presenting arbitrary combinations of small tunnel junctions, capacitances, and voltage sources. The method is based on numerical solution of a linear matrix equation for the vector of probabilities of the most important electric charge states of the system, with iterative account of new states. The method is able to describe very small deviations from the ideal behavior of a system, due to finite speed of applied signals, thermal activation, and macroscopic quantum tunneling of charge (cotunneling). The code is portable across a number of UNIX based platforms, including the Intel parallel computers. In computer time consuming applications, the code's performance can be harnessed by combining resources of typically heterogeneous networked computing platforms linked by PVM. Several single-electron systems were analyzed, showing that the method is computationally efficient and robust. We have studied the probability of dynamic and static errors in a particular single-electron memory cell, the single-electron trap. We have also studied several single- electron standards of dc current. where a new device and a new and more efficient rf drive waveform have been proposed. Finally, we have performed an optimization of those parameters of a single-electron pump which may influence the accuracy of this device as a standard of dc current. The common result among these studies is that the theoretical limitation on the performance of single-electron devices is much less restrictive than what has been obtained experimentally. Intersection Detection and Bifurcation of Triangulated Surfaces in Three Dimensions by Qian Li, Advisor: Brent Lindquist May, 1996 We describe a three dimensional interface intersection detection and untangling algorithm in the context of a front tracking method for fluid flow problems. Front tracking is a numerical method for solving systems of partial differential equations whose solutions, as functions of space and time, are known to have discontinuities with respect to space variables. Front tracking method is characterized by the explicit declaration and tracking of an interface that represents the locations where solution discontinuities occur. Interfaces are assumed to be consisting of hypersurfaces of codimension one. An interface decomposes a computational region into a number of components separated by its surfaces. The solution sought varies smoothly within each component but changes abruptly across their boundaries. The propagation of interfaces, performed at every time step, is an integral part of advancing the numerical solution in front tracking method. An interface in three dimensions is basically a set of surfaces. The surfaces may be connected at their boundaries which are declared curves, but they do not intersect one another or with themselves otherwise. This is so that each surface separates a unique pair of space components. Dynamic evolution of an interface sometimes causes its surfaces to cross one another despite the nonintersection requirement, and changes the topology of the system. When this occurs, certain topological surgery on the interface becomes necessary to: (1) make sure that the topology of the new interface correctly reflects the state of the underlying system; (2) make sure the new interface still satisfies the nonintersection condition. We call this process of topological surgery bifurcation. Two real life examples of interface intersection and bifurcation are shock wave bifurcation and bubble merger. Bifurcation of an interface depends on the geometry of its interface, as well as the physics of the underlying system. In this dissertation we study the physics independent aspects of bifurcation. Topics covered include intersection detection, surface retriangulation and intersection untangling. A Front Tracking Method for Regularization-Sensitive Shock Waves by Hyun-Cheol Hwang, Advisor: Bradley Plohr August, 1996 The regularization-sensitive waves, which often appear in non-strictly hyperbolic and mixed elliptic-hyperbolic systems of conservation laws, are investigated. These weak solutions do not satisfy the Lax entropy criterion, and yet they are admissible under the viscous profile criterion; however they are sensitively dependent on the parabolic regularization process. This sensitivity creates computational difficulties because of the artificial regularization, i.e., numerical diffusion, that is present in most numerical schemes. A new version of the front-tracking algorithm is developed to circumvent these difficulties. To propagate the discontinuous wave accurately, the algorithm utilizes a Riemann problem solver which directly accounts for the internal wave structure. Specifically, this Riemann solver uses the multiple shooting method to compute the saddle-to-saddle connection in the corresponding dynamical system. The results of front-tracking simulations of regularization-sensitive waves are compared with the results of similar computations using common shock capturing schemes such as Lax-Friedrichs, Lax-Wendroff, and higher- order Godunov. Two types of regularization-sensitive waves, of great theoretical and practical interest, are considered: transitional shock waves in quadratic flux models, and deflagration waves in combustion. It is found that only the front-tracking algorithm correctly and accurately resolves these waves. Nonequilibrium Condensation behind a Shock Wave in Retrograde Vapor by David Marc Saltz, Advisor: Bradley Plohr August, 1996 I have investigated the behavior of a one-dimensional shock wave propagating through a superheated retrograde vapor. A forerunner shock wave supersaturates the vapor, which condenses along a trailing condensation wave. A numerical simulation of the time-dependent flow reveals the nature of the approach to a steady-state partial liquefaction shock wave. I have identified four regimes of qualitatively different flow, characterized with respect to the laboratory observation time and most conveniently expressed by the quasi steady nucleation rates behind the frozen (initial) and equilibrium (steady state) forerunner shock waves. The equation of state of the vapor phase is based on a phenomenological cluster expansion, utilizing a cluster free energy model valid at small cluster size. The homogeneous nucleation rate is derived from the same free energy model with the conventional assumption of detailed balancing. Droplets are nucleated at larger-than-critical size and a size-dependent nucleation induction time is used to make the numerical solution insensitive to the arbitrary choice of initial droplet size. Droplet growth, occurring in the transitional and continuum Knudsen regimes, is modeled using Fuchs' interpolation formula, modified to include condensation and evaporation events involving dimers and larger clusters. The conservation laws of one-dimensional inviscid gas dynamics, supplemented by three rate equations describing the nonequilibrium condensation, are solved numerically using the front-tracking and operator- splitting methods. Front tracking is needed to accurately determine the supersaturated vapor state behind the shock front, while operator splitting is required for computational expedience because the condensation variables evolve on time scales that can be much smaller than the time scale of the macroscopic fluid motion. Because of the high heat capacity of the fluid, the dominant effect of the condensation is to lower the pressure behind the forerunner front, causing the forerunner wave to decelerate over time as the condensation wave widens. The system approaches a steady state consisting of a single thick shock wave. Only at high Mach number, when the nucleation rates behind the frozen and steady-state forerunner waves are sufficiently high, does the system approach a steady state rapidly with respect to the laboratory observation time. Contaminant Transport and In-Situ Bioremediation in Porous Media by William J. Thistleton, Advisor: Brent Lindquist August, 1996 Single phase fluid flow through heterogeneous porous media is studied. PDE's coupled through nonlinear reaction terms are used to model the in-situ bioremediation of a contaminant and are solved numerically. The numerical scheme is incorporated into the Front Tracking code of Glimm, et al. A prediction for the rate of progress of a reaction front between resident substrate and invading oxygenated water is developed and validated for the case of an instantaneous reaction. Numerical simulations in two spatial dimensions are performed and presented for the reacting system modeled with Monod kinetics. A mesh refinement technique has been developed and implemented in two spatial dimensions in order to reduce the time required for numerical solution. A particle tracking method has been developed and implemented to study tracer flow in correlated and uncorrelated heterogeneous porous media. This method dramatically reduces the time required for numerical solution compared to conventional front tracking. The effects on growth of a mixing layer between resident and invading fluid due to increasing importance of fluid flow transverse to the predominant direction of flow are studied. A modified particle tracking method has been developed and implemented. The result concerning the rate of progress of the reaction front between invading oxygenated water and ambient hydrocarbon is extended to two spatial dimensions. Simulations of fluid flow through families of correlated and uncorrelated porous media are presented and discussed. Parallel Simulated Annealing and Applications by Chung-Chiang Chou, Advisor: Yuefan Deng December 1996 We address three issues: (1) design of novel cooling schedule leading to faster simulated annealing algorithms, (2) efficient parallelization of the algorithms on distributed-memory MIMD systems, and (3) applications of these algorithms to problems such as Ginzburg-Landau equation which is solved by a new approach of renormalized energy minimization by the simulated annealing algorithm and molecular confirmation by finding low energy states of molecular clusters of identical atoms. A Numerical Study of the Richtmyer-Meshkov Instability in Cylindrical Geometry by Mary Jane Graham, Advisor: Qiang Zhang December, 1996 As an incident shock wave collides with a material interface between two fluids of different densities, the interface becomes unstable. Small disturbances at the interface start to grow. This interfacial instability is known as a Richtmyer-Meshkov instability and is the subject of this work. So far theoretical, numerical and experimental studies of the instability have been performed in plane geometry. In most physical applications-such as inertial confinement fusion and supernova-the Richtmyer-Meshkov instability occurs in curved geometry. The method of front tracking has been used to perform a systematic numerical study of the Richtmyer-Meshkov instability in cylindrical geometry for both the imploding and the exploding cases, including the cases of the incident shock wave imploding (exploding) from a heavy to light fluid phase or a light to heavy fluid phase. The curved geometry complicates the system considerably. For example, the unperturbed system does not have an analytic solution; the unperturbed system in plane geometry does; the occurrence of re-acceleration or reshock at the material interface caused by the waves reflecting back from the origin is unavoidable in curved geometry. This paper addresses those issues. In addition, a detailed study is presented of small amplitude perturbation solutions, nonlinear solutions, the effect due to the number of fingers at initialization, the effect of the initial perturbation amplitude, the phenomenon of phase inversion and twice phase inversions, the growth rate of the spikes and bubbles formed by the instability (including the overall growth rate), and the similarity scaling law for Richtmyer-Meshkov instabilities by shocks of large Mach number. The Analysis of Continuous Interfaces Using Compressible Multi-component Flow with Front Tracking by Andrea Lynn Pass, Advisor: John Grove December 1996 This dissertation develops numerical algorithms for the computation of multi- component compressible fluids. We use front tracking to isolate regions of initial high concentration gradient; the tracked boundaries eliminate excessive numerical diffusion. These algorithms serve two purposes. The first is to study physical situations in which the above hypotheses are important features, and the second is to analyze in controlled numerical analysis experiments the influence of tracking on shock contact interactions. Numerical simulations have consistently overestimated the growth of the mixing region in comparison to the experiments. The multi-component method allows for an untracked curve to separate fluids in a diffusive manner. The untracked interface models a slightly slower growth rate, leading to the belief that numerical diffusion may influence the shape of the mixing region. Analytical Nonlinear Theory of Unstable Fluid Mixing Driven by a Shock Wave by Sungik Son, Advisor: Qiang Zhang December, 1996 We present a nonlinear theory of the growth of fingers at interface between two fluids of different densities driven by a shock wave. This interfacial instability is known as the Richtmyer-Meshkov (RM) instability. Previous theoretical work focused on linear regime and failed to give a quantitatively correct prediction for the growth rate of Richtmyer-Meshkov unstabie interface in the nonlinear regime. Spikes and bubbles are formed at the unstable interface. A spike (bubble) is a portion of heavy (light) fluid penetrating into light (heavy) fluid. The theory presented in this dissertation provides analytical, explicit expressions for quantitative predictions of the overall growth rate as well as the growth rate of the spike and bubble of the unstable interface between fluids of arbitrary density ratios in two and three dimensions. Our theory contains no adjustable parameters. The theoretical predictions of our nonlinear theory are in excellent agreement with results of full nonlinear numerical simulations and the experimental data in two dimensions from the linear (small amplitude) to the nonlinear regimes. The predictions of linear theories are qualitatively incorrect at late times. Our theory has identified that the RM unstable system goes through a transition from a compressible and linear one at early times to a nonlinear and incompressible one at later times. The nonlinear theory has been extended to fluids in three dimensions. Our results show that the growth rates for different orientation angle q with fixed total wave number k are qualitatively similar, but quantitatively different. As a part of the theoretical derivations, we include the com- pressible linear theory in three dimensions for reflected shock and reflected rarefaction wave cases. Diffusion Problems in Fluid Flow Models by Kent Sun, Advisor: Reginald Tewarson December, 1996 In this dissertation, accurate numerical schemes for solving a system of nonlinear differential-algebraic equations are presented. Such problems arise in segmental models of the proximal tubule of the mammalian kidney. Specifically, methods for solving a system of nonlinear differential- algebraic equations of the form dy/dx +f(y,z) = 0, A(y,z) = 0, and dU(y,z,e,z')/dx' + V(y,z) = 0 are considered. It is shown that. for e = O the system of equations is relatively easy to solve. For c ? 0, however, the solution process is quite difficult not only due to the increase in the order of the system, but also from the presence of the small coefficient in the second derivative terms. An error analysis is given and the need for singular perturbation methods is discussed. Parallel Simulation and Optimization for Inventory Systems and Biosystems by Yuan Wang, Advisor: Yuefan Deng December 1996 The first part of the thesis studies two stochastic multivariate inventory models of multi-location and multi-item. A parallel heuristic algorithm based on the genetic algorithm is developed for solving the multi-location model while a parallel algorithm coupling the SUMT method, function-assembling, and pattern search technique for the multi-item model. The second part analyzes DNA sequences binded by protein in water solution. Our results demonstrated the previously unknown roles of water molecules in protein binding to DNA. Small Amplitude Theory of Richtmyer-Meshkov Instability in Cylindrical and Spherical Geometries by Ju-Hong Kim, Advisor: Qiang Zhang August, 1997 When an incident shock wave collides with an interface between two different materials, small perturbations of this interface grow into nonlinear structures. These interfacial perturbations grow in the form of bubbles and spikes. Bubbles (or spikes) are the portions of light (or heavy) fluid penetrating into heavy (or light) fluid. This interfacial instability is known as Richtmyer-Meshkov instability. It plays important role in the study of inertial confinement fusion and supernova. Most of theoretical, numerical, and experimental studies have been performed in planar geometry. However, in both inertial confinement fusion and supernova applications, the Richtmyer-Meshkov instability occurs in spherical geometry. Changing the geometry from planar to cylindrical or spherical results in a much more complicated system. Even for compressible fluid with cylindrically or spherically symmetries, which is a one dimensional problem, the analytical solutions have not been found. We present theoretical formulation for compressible fluids in cylindrical and spherical geometries with small perturbation at the material interface. Furthermore, small time solutions for the cylin-drically or spherically symmetric flow and for the perturbed system are determined analytically. The small amplitude solution, formulated in this thesis will be important for calibration of results from full numerical simulation of compressible fluids in cylindrical and spherical geometries. Analytical Approximate Solution for the Prices of American Exotic Options. by Tatiana S. Taksar, Advisor: Qiang Zhang August, 1997 In this thesis we will present analytical approximate solutions to the values of American barrier, lookback and Asian options. Both "out" barrier and "in" barrier options have been considered. At an "out" ("in") barrier, an option becomes worthless (starts to exist). A path-dependent option is an option whose payoff at exercise or maturity depends on the past history of the price of the underlying asset. Lookback options are one of the most common types of path-dependent options. They are path-dependent options whose payoff depends on the maximum or the minimum realized price of the underlying asset over the life of the option. Asian options are path-dependent options whose payoff depends on the average stock price observed during the life of option. So far, the analytical solutions to the American options have not been found for vanilla options, the simplest form of American options. American options with barriers and path-dependent options contain additional complications. The analytical approximate solution derived in this thesis can be used to price the values of barrier, lookback and Asian call and put options on stocks, commodities and commodity futures and currency exchanges. The theory presented in this thesis provides simple analytical formulas for evaluation of some types of American exotic options. To validate results, we compared the predictions of analytical approximate solutions with those of numerical finite-difference methods. We considered not only the values of options, but also predictions for critical stock price, or free boundary, since this quantity is crucial for the estimation of option prices. The predictions of our analytical approximate solutions are in excellent agreement with the numerical solutions for cases with constant parameters, such as interest rate, dividends rate and volatility. For the case with time dependent parameters, namely Asian options case, the approximate solutions were less accurate for some sets of parameters. As a part of validation tests, we include a systematic validation study of quadratic approximation over the range of parameters. Numerical Solutions of Differential Equations Arising in Flow Networks. by William Toro, Advisor: Reginald Tewarson December, 1997 In this dissertation we investigate and give efficient and accurate numerical methods for multipoint boundary value problems that arise in flow network models of the counter current multiplier of the mammalian kidney. The algorithms presented are based on centered implicit multistep schemes to solve iteratively the nonlinear equations. The solutions obtained are O(h5) locally and O(h4) globally without needing extra storage or doing too many operations. Extra operations are required only in the solution of the first four levels of the algebraic equations. After that these algorithms use the same storage amount and almost the same number of operations as in the methods based on the trapezoidal rule, which is O(h^3) locally and O(h^2) globally. Inverse Problem Algorithms and Applications in Renal Concentrating Mechanism Models by Mariano Marcano, Advisor: Reginald Tewarson May, 1998 An algorithm to solve the inverse problem for the renal concentrating mechanism is given. The main part of the inverse algorithm is an iterative process which leads to a constrained minimization problem. A method to solve this problem by a generalized inverse algorithm is discussed. Computational results are presented showing that the proposed iterative process is much faster and more accurate than other optimization methods. The inverse algorithm is applied to a central core model. Computational results are given showing that this model would generate a significant concentration gradient if nearly fifty percent of the parameters have half the values of their experimental lower bounds. Periodic Structure in Riemann Solutions for Hamilton-Jacobi Equations by John Pinezich, Advisor: James Glimm August, 1998 This dissertation investigates the structure of two-dimensional Riemann problems for Hamilton-Jacobi equations. We show that it is possible for the Riemann solution to have closed characteristic orbits, enclosing furthermore a periodic sonic structure, which in turn encloses a parabolic structure. This investigation was prompted by a numerical construction of a Riemann solution for a particular example displaying an even richer internal structure. Since there are no (known) robust methods to numerically construct two-dimensional Riemann solutions, this dissertation shows that examples of Riemann problems with Riemann solutions of comparable complexity do in fact exist. A Theoretical Study of Trace Flow Through Heterogeneous Porous Media in Curve Geometry by Won-Suck Lee, Advisor: Qiang Zhang August, 1998 We present fluid mixing behavior and mixing zone growth patterns under radially diverging flow conditions, caused by the constant injection of fluid form a source well which is located within a porous medium. In order to study the nature of the fluid flow near the well we consider a flow in a curved geometry. Interfacial fluid mixing induced by heterogeneous porous media will be analyzed in both cylindrical and spherical geometries and compared with the plane geometry case. To concentrate our study specifically upon the geological difficulty, we perform a tracer flow study, in which only geometrical complexities occur and non-linear flow features are absent. The non-uniform, in the mean, velocity field is formed by the radially diverging fluid from the injection point. This flow feature makes the velocity correction function non-homogeneous and thus contributes a non-local dispersion term to the governing equation of tracer flow in random porous media. We derive the functional form of the second moment of the fluctuation of the velocity field and deduce the approximate governing equation by using Corrsin's hypothesis. our governing equation, a convection-dispersion integro-differential equation with flow history dependence, and the solution we obtain show that the dispersion of fluids should not be considered locally. The stronger the heterogeneity of the porous media, the non-local consideration of the dispersion term is more essential, regardless of the geometry. The finite volume method is successfully utilized as a numerical scheme for our history dependent integro-differential equation. In place geometry, homogeneity of the velocity correction gives a further approximation of the dispersivity. It can be approximated by a diffusion term through neglect of the fourth order term with respect to dn, when ||dn|| < 1. However, it is not possible to make the approximation in curved geometry because of the non-homogeneity of the velocity correlation function. The lower order dn term, which is not present in the planar case, contributes to the dispersion term in the curved geometry case. It acts as a stress which yields, from the physical point of view, expansion of the zone with increasing time around the point where the value of concentration is one half of the injected concentration level. Late Time Phenomena of Single Mode Rayleigh-Taylor Instability by An-Der Lin, Advisor: James Glimm December, 2000 The development of single mode Rayleigh-Taylor (RT) instability, consists of three different stages, described as the linear, free fall and terminal velocity regimes. This instability and its three regimes have been studied in detail \cite{Zha91a}. The central result of this thesis is the discovery of a fourth regime of oscillatory or nonuniform behavior after an apparent terminal velocity has been reached, especially for fluids whose density constrast is not too large. This result is based on numerical investigation by the Front Tracking method. The discovery documented here requires a modification to common beliefs regarding the late time evolution of a single mode disturbance. After a short time period of spike (bubble) pseudo terminal velocity plateau, the velocity resumes its increase reaches a new peak, and then decreases. The overall velocity development is sensitive to the Atwood number of the fluids. We show a relation between the spike (bubble) velocity and the tip curvature evolution. A linear relation between the spike (bubble) velocity and a pressure difference at the tips is also found. This pressure difference is the difference between the pressure at the spike (bubble) tip and the ambient pressure at the tip. The pressure difference reflects the pressure gradient at the tip. Numerical evidence shows that the pressure difference is strongly correlated to the spike (bubble) velocity development. The purpose of this work is to report these newly observed phenomena. Automated Recognition Algorithms for Neural Studies by Ingrid Koh Advisor: Brent Lindquist May 2001 Morphological characterization in the biological sciences is important for the study of structure-function relationships. As new imaging approaches progress and the speed of computation increases, it has become possible to acquire large sets of two- or three-dimensional high resolution images. However, this advance in data acquisition technique has not been accompanied by comparable advances in data analysis techniques. This thesis describes a number of geometric-based automated image analysis algorithms that are applied to neurobiological studies, including the tracing of neurons and the detection of neuronal dendrites and their spines and the quantification of the escape behavior of larval zebrafish from high speed camera recordings. Neuron morphology is characterized by quantitative measurements obtained from the tree model of a neuron. Morphological plasticity in dendritic spines is quantified via detection and measurements of spine length, volume, density and shape classification for both static and time-lapse images. Accuracy of the spine detection algorithms is addressed by comparing the automated and manual spine length and density measurements. Normal escape behavior of larval zebrafish is quantified via constructing the midline of the fish. The approaches presented here are generalizable to other aspects of neuronal morphology. A Posteriori Error Estimate for Front-Tracking by Marc Laforest Advisor: James Glimm Dec 2001 We demonstrate an a posteriori error estimate in the $L^1$ norm for front-tracking approximate solutions to hyperbolic systems of nonlinear conservation laws. Extending the $L^1$-stability result of Bressan, Liu, and Yang we use their $L^1$-equivalent functional for pairs of front-tracking approximations and identify the leading order contribution to the numerical error. This leading term is closely related to the residual of the approximation and determines an a posteriori bound of the error. We demonstrate the estimate for the front-tracking approximations of Risebro, which are extensions to systems of Dafermos' polygonal approximations. The Incompressible Limit of Compressible Multiphase Flow Equations by Hyeonseong Jin, Advisor: James Glimm December, 2001 We study the motion of the slightly compressible multi-phase flow model proposed by Chen, Glimm, Sharp and Zhang. The multi-phase flow model provides an averaged, or coarse gained, description of the mixing layer formed when an unstable interface between two fluids is driven by acceleration. The interface velocity and constitutive law are analyzed by derivation of the exact quantity. A zero parameter model is presented with the velocity constitutive law given implicitly in terms of the solution. The equations of compressible multi-phase flow in appropriate nondimensional form are a nonlinear hyperbolic system depending on a large parameter $\lambda$, the reciprocal of the Mach number. The incompressible limit of the compressible multi-phase flow equations is a time-singular and layer-type problem which requires advanced techniques in asymptotics. We discuss the limiting behavior of the solutions of the compressible multi-phase equations as $\lambda \to \infty$. Using singular perturbation theory, we derive asymptotic expansions of the solutions of the compressible equations, describing the incompressible limit process, $\lambda \to \infty$, which are formally uniformly valid on a bou nded time interval. A necessary and sufficient condition for convergence of compressible pressures in the second order of the asymptotic expansion to the incompressible pressures is derived. An additional degree of freedom exists in the pressures for the incompressible limit. This condition fixes a degree of freedom related to the relative compressibility of the two fluids in the incompressible pressures. There exists a substantial interaction between an essential incompressible mean flow and small amplitude acoustical waves. These waves are highly oscillatory in time and appear as solutions to wave equations in a first order approximation of the velocities and in a second order approximation of the pressures. Our purpose here is to prove rigorously that the se waves are negligible in the incompressible approximation to compressible flow. Linearized Analysis of the Richtmyer-Meshkov Instability for Elastic Materials by Jee-Yeon Nam Advisor: James Glimm December 2001 We present a study of the Richtmyer-Meshkov instability for elastic materials. This instability, which is a shock-driven interfacial instability, was originally investigated for fluids. We consider two elastic materials in contact along a frictionless surface. The governing system of equations comprises conservation laws supplemented by constitutive equations. To analyze it, we linearize the equations around a one-dimensional background solution under the assumption that the perturbation is small. The background problem defines a Riemann problem that is solved numerically; its solution contains transmitted and reflected waves in the longitudinal modes. The linearized Rankine-Hugoniot condition provides the interface conditions at the longitudinal and shear waves, while linearization of the frictionless material interface condition supplies the material interface conditions. The resulting equations, a linear system of PDEs, is solved numerically using front tracking. To validate the numerical code, we reproduce the result for the fluid case analyzed by Richtmyer. The main result is a calculation of the growth rate of the perturbation amplitude. Trends concerning the effect of the nonzero shear modulus on the growth rate are shown and discussed. Optimal Parallelization of Simulated Annealing by State Mixing by King-Wai Chu Advisor: John Reinitz and Dec 2001 Yuefan Deng This thesis describes a new, efficient, and general purpose parallel simulated annealing algorithm. The algorithm is based on periodic mixing steps, in which favorable states reproduce and unfavorable ones are destroyed. It runs on a distributed memory Multiple Instructions Multiple Data architecture parallel computer. Parallel efficiency is controlled by the interval between mixing steps. In this thesis, it is shown that for certain values of this interval found by exhaustive search, the algorithm can give up to 100\% parallel efficiency on up to 50 processors and 80\% parallel efficiency on 100 processors. Moreover, for a given number of processors, there is a range of mixing interval which gives high parallel efficiency. In this thesis, two efficient statistical estimators, namely, the cross-correlation and variance among processors are defined for finding efficient mixing intervals are constructed which give parallel efficiency of 75\% without exhaustive search. This is done by tracking the two statistical estimators right after communication, so as to obtain the lower and upper bounds for the optimal mixing interval. The algorithm was tested on two problems. One was an inverse problem on gene expression dynamics in developmental biology, and the other was the Lennard-Jones potential problem. A Numerical and Experimental Study of a Shock-Accelerated Heavy Gas Cylinder by Cindy Anne Zoldi Advisors: James Glimm and May 2002 David Sharp (LANL) We study the evolution of an SF6 gas cylinder surrounded by air when accelerated by a planar Mach 1.2 shock wave. Vorticity generated by the interaction of the shock wave's pressure gradient with the density gradient at the air/SF6 interface drives the evolution of the cylinder into a vortex pair. This thesis contains two interrelated parts; the first part concerns the acquisition of experimental data and the second part concerns the use of this data to benchmark simulations of the experiment with the hydrodynamics code RAGE. RAGE, an adaptive-mesh Eulerian code, has had previous success in simulating shocked interfaces. Improved experimental diagnostics were used to acquire data, which allowed us to perform a more stringent test of the code's capabilities. From each shock tube experiment, we obtained one image of the experimental initial conditions and six images of the time evolution of the cylinder. Moreover, the implementation of Particle Image Velocimetry (PIV) also allowed us to determine the velocity field at the last experimental time. This thesis is the first code validation study of a shocked flow to use two- dimensional velocity field data for comparison. Simulations incorporating the two-dimensional image of the experimental initial conditions led to good qualitative agreement with the experimental images. Comparing length measurements of the evolving cylinder and velocity vectors at the last experimental time led to quantitative differences, particularly between the measured and computed velocity magnitudes. The computational study carried out in this thesis showed that agreement between the measured and computed velocities could be achieved by decreasing the peak SF6 concentration and diffusing the air/SF6 interface in the experimental initial conditions. These modifications are consistent with the observation that the SF6 gas diffuses faster than the glycol droplets used to track the gas. Further experimental improvements, including a better characterization of the experimental initial conditions, are discussed. This thesis demonstrates that quantitative measurements, in addition to qualitative images, should be examined when comparing experimental data and computational results. The quantitative differences between the measured and computed velocity magnitudes led to the conclusion that the experimental initial conditions were not well characterized. A Parallelized Point-Shifted Tetrahedral Grid for the Finite Element Method by Wei Guo, Advisor: Brent W. LIndquist May 2002 FronTier is a software package developed for the solution of partial differential equations using the front tracking technique. FronTier tracks a discontinuity interface and generates grids that conform to the interface. For parabolic and elliptic problems, it uses the so-called point-shifted triangular grid for the finite element method in two dimensions. In this thesis, we extend the method to three dimensions. The point-shifted triangular grid is a hybrid structured-unstructured grid scheme, providing local unstructuring within the confines of a global rectangular-based structure. At the first stage, a structured rectangular grid is constructed and point-shifted. At the second stage, each grid block is tetrahedralized in an unstructured way. The ideas, algorithms, theories, and implementations including data structures for the parallelized point-shifted tetrahedral grid construction will be discussed. The main focus is on the theoretical grounds for the tetrahedralization and the algorithm for the parallelization. Parallel Computation of Radiative Heat Transfer in an Axisymmetric Closed Chamber: An Application to Crystal Growth by William Garber Advisor: James Glimm Aug 2002 Radiative heat transfer plays an important part in the Czochralski crystal growth process. This dissertation discusses the computation of radiative heat exchange in the crystal and the gas above the melt. Parallel computation of view factors for grey diffuse surfaces is studied using a combination of methods developed by the author and other researchers, exploiting the strengths and advantages of the respective methods. Kobayashi and Miyahara introduced the use of a path integral to calculate the view factor and this method is accurate and robust, especially near the corners which are typically problematic. Dupret et. al. have a less accurate means of calculating the view factor integral but their visibility or blocking algorithm allows analysis of completely arbitrary axisymmetric closed geometries. In this dissertation we introduce optimizations to the Dupret blocking algorithm. The computation of view factors and blocking requires more computer time than solving for the flux, and the parallel code developed exhibits excellent scalability. The code was applied to the design of a heat shield with a controlled temperature. The growth velocity and temperature gradients are measured and studied in terms of the constraint imposed by the equation for the formation of an oxidation induced stacking fault (OSF) ring, Vg/Gs > Ccrit. Protein-DNA binding simulation on parallel computers by Alexandre Korobka Advisor: Yuefan Deng Aug 2002 The central part of this thesis is application of algorithms for analysis of binding specificity of protein-DNA complexes to complexes containing water-mediated bridges. We show that inclusion of water contacts improves the quality of prediction; examined data sets yield correct results for all base pairs that have two or more hydrogen bonds. Our prediction algorithm uses square-well, quadratic, and Lennard-Jones approximations of hydrogen bond potential in point set pattern matching, quadratic energy minimization and very fast simulated reannealing stages respectvely. The following assumptions regarding protein-DNA binding are made: (1) specificity depends on direct hydrogen bonding; (2) electrostatic forces contribute strongly to stabilization and weakly to specificity; (3) most of Van der Waals and electrostatic interactions can be neglected. A data decomposition algorithm for faster evaluation of electrostatic potential of proteins in rigid body simulation is proposed. We show that it produces expected results with precision suitable for modeling of biologic systems including protein-DNA complexes. We present a parallel computer system built and used in calculations for this work. A Throat Finding Algorithm for Medial Axis Analysis of 3D Images of Vesiculated Basalts by Hyunkyung Shin Advisor: Brent Lindquist Aug 2002 We present a high resolution study of the void space geometry of vesiculated basalts (porosities in the range 61\% to 80\%) from three dimensional digitized images obtained by synchrotron X-ray tomography. As the void space is composed of vesicles, the solidified remnants of expanding gas bubbles, the contact surface between coalescent vesicles represents a pore throat. We discuss a throat finding algorithm used in digital analysis to locate such contact surfaces. The algorithm requires a segmented image from which a medial axis skeleton has been computed. We present results for the distribution of vesicle volumes, contact surface areas, vesicle coordination number (the number of neighbors with which a vesicle has coalesced), and correlations between these observable. The distribution for contact surface area is log-normal (down to the resolution of an individual voxel); the distribution for coordination number is exponential; and the distribution of vesicle volumes is not simply characterizable. There are strong correlations between vesicle size/coordination number and vesicle size/contact surface area. A Study of Predictability in Stochastic Reservoirs and Scale-Up by Seung Yeon Cho Advisor: Brent Lindquist Aug 2002 For a reservoir model consisting of purely stochastic, a priori, small scale geological information we explore the a posteriori distribution of reservoir oil production behavior. In particular we are interested in the inherent limits for ensemble based prediction. Our model contains fluid and geometrical complexity representative of real reservoir behavior but is sufficiently simple to allow rapid solution and thus enable deeper analysis of the mapping from geology to production. Our results show that production uncertainty has a strong dependence on both the heterogeneity strength and the spatial correlation length of the geology and is only weakly related to the correlation between porosity and permeability. We estimate the minimum number of realizations needed to achieve a given level of accuracy of ensemble based means. We quantify the improvement of accuracy obtained from history matched ensembles. We show that the prediction of reservoir production has a Markovian property. This implies the pessimistic though not unexpected conclusion that no single realization, even one that is history matched, is likely to retain the same level of predictive accuracy for an extended period of time. In part 2, scale-up computation in reservoir engineering is also studied. Front tracking scheme is applied in order to reduce numerical dispersion. We implimented two different scale-up algorithms to {\it FronTier}. For validation, fine scale geologies for realistic reservoir cross sections are assigned on a $100\times 100$ grid, and coarse grid simulations are performed on the upscaled grids corresponding $5$ discretizations varing from $2\times 2$ to $20\times 20$. Fluid Flow in a Rock Fractures using Finite Difference Lattice Boltzmann Method by Imbunm Kim Advisor: Brent Lindquist Aug 2002 Fluid flow through a rock fracture is often approximated by the parallel plate model for which flow rate is proportional to the cube of the average or effective aperture width. It has been of longstanding interestto determine how actual fracture surface roughness and irregularity of contact modify the predictions of the parallel plate model. In this thesis, fluid flow through a mechanical fracture in Harcourt Granite is studied. We made numerical definitions of percolation threshold, fracture aperture, flux and permeability. Numerical results are compared with laboratory measurements in some fractures. For the numerical flow studies, the lattice Boltzmann method (LBM) has been used. To accommodate the vastly different aspect ratios of a fracture, we utilize the finite difference lattice Boltzmann method (FDLBM) in order to apply nonuniform grids. The FDLBM implementation was validated on simple fluid flow simulations: Poiseuille flow between two parallel plates with and without a constrictive neck and decaying Taylor vortex flow. For the simulation of fluid flow in a real fracture, the velocity field and fracture permeability were obtained for different values of mean aperture(confining pressure). The numerical values of the permeability as a function of confining pressure agree qualitatively with experimental results. The FDLBM computations are CPU intensive. This was addressed using parallel computation on a Beowulf class cluster of processor Fluid Mixing: Rayleigh-Taylor Mixing Rate and Diesel Fuel Injection Spray by Andrea Marchese Advisor: James Glimm Dec 2002 We present a Rayleigh-Taylor mixing rate simulation with an acceleration rate falling within the range of experiments. The simulation uses front tracking to prevent interfacial mass diffusion. We present evidence to support the assertion that the lower acceleration rate found in untracked simulations is caused, at least to a large extent, by a reduced buoyancy force due to numerical interfacial mass diffusion. Quantitative evidence includes results from a time dependent Atwood number analysis of the diffusive simulation, which yields a renormalized mixing rate coefficient for the diffusive simulation in agreement with experiment. Also reported in this thesis are the preliminary results on the numerical simulation of diesel fuel injection spray. The breakup of injected fuel into spray is of key interest to the design of a fuel efficient, nonpolluting diesel engine. Our simulation design is set to match experiments at ANL, and our present agreement is semi-quantitative. Future efforts will include mesh refinement studies, which will better model the turbulent flow. Magnetically driven Rayleigh-Taylor instability by Cahrels J. Ju Advisor: Bradley Plohr Dec 2002 In this thesis, we study the collapse of a cylindrically shaped fluid shell under Lorentz forces created with large axial electric currents. The collapse of this shell suffers a breakdown due to Rayleigh-Taylor (RT) instability. We are interested in the effects of magnetic flux on the growth of this instability. The magnetic diffusion into the fluid is expected to slow the rate of growth of the instability. The analytical formulation involves Maxwell's equations of electrodynamics and Euler's equations of fluid dynamics. Theoretical analysis yielded a set of perturbed first order equations together with analytical boundary conditions. The boundary conditions were derived by applying the Biot-Savart law at perturbed surface and with perturbed current, thus avoiding the evaluation of certain vacuum boundary conditions. The governing equations and boundary equations were discretized and remains to be implemented in the computer code. We propose a numerical experiment based on the results of this analysis. Estimation of Computational Simulation Errors in Gas Dynamics by Yunghee Kang Advisor: John Grove May 2003 Statistical models in numerical solution errors for shock physics problems are discussed. The statistical Riemann problem is a key issue in this paper, and in particular the numerical solution of wave interactions. Our goal is to develop models for the probability distribution of solution errors in numerical simulations of gas dynamical flows. Randomness in solution error occurs for a variety of reasons, including finite resolution of the numerical solution, and the dependency of the solution error on uncertain input parameters. Basically, we focus on how the initial conditions, the numerical method, and computational parameters such as mesh size, affect solution errors for the variables (\eg, pressure, density, or specific internal energy) that describe the flow state and observables such as wave arrival times and wave locations. These initial conditions include wave strengths, initial densities and velocities, numerical methods including second order Godunov type schemes, the Lax-Wendroff and Lax-Friedrichs methods. Computational parameters include mesh refinement as well as limiting or artificial viscosity parameters. Errors in Numerical Solutions of Shock Physics Problems by Yan Yu Advisor: James Glimm Dec 2004 We seek error models for shock physics simulations which are robust and understandable. We propose statistical models of uncertainty and error in numerical solutions. To represent errors efficiently in shock physics simulations we propose a composition law. The law allows us to estimate errors in the solutions of composite problems in terms of the errors from simpler ones. We formulate and validate this composition law for shock interactions in planar geometry. We also explore complications introduced by spherical flow in the analysis of errors in the numerical solutions. We illustrate that idea in a very simple context. For shock interactions in spherical geometry, we conduct a detailed analysis of the errors. One of our goals is to understand the relative magnitude of the input uncertainty vs. the errors created within the numerical solution. In more detail, we wish to understand the contribution of each wave interaction to the errors observed at the end of the simulation. Fluid Displacement in Rock Cores: A Study Based on Three Dimensional X-ray Microtomography Images by Masa Prodanovic Advisor: Brent Lindquist August 2005 The study of multiphase flow in porous media is of major interest in oil recovery and contaminant remediation in aquifiers. Porous media are modeled as networks of pores and throats (channels). These network flow models depend on a number of parameters supplied by experiments. Interpretation of experimental measurements of a porous medium property often rely heavily on idealized models and equations. In contrast, the development of robust techniques for analysis of X-ray microtomography images can provide direct quantification of geometric parameters for the pore space and fluids resident therein. Such results then serve either as input to network models via parameter estimation or for model verification. In moderately high porosity (40% to 50%) samples, successful throat finding and pore partitioning has proved particularly challenging. We present improved throat finding methods and introduce a robust pore partitioning scheme. We present results on pore characterization in Berea sandstone and polyethylene cores in the form of parameter distributions required as input to network models. The parameters include volume, surface area, shape factor, principal direction diameters and permeability for both pores and throats. Lattice Boltzmann computations of absolute permeabilities for individual pores/throats are based on exact (to the resolution of the image) configurations and can replace current permeability estimates produced from analytic models based upon geometrically simplified shapes. In polyethylene core experiments, fluid displacement patterns were obtained both after fixed volumetric injections of invading fluid, and after steady state conditions were reached at fixed fractional flow injection. We concentrate our fluid analysis on overall fluid phase connectivity, interphase contact area, as well as saturation distribution in pores of different sizes. We also report on the ability to differentiate three fluid phases (oil, water, air) in the pore space. In extensive studies of the Berea core it was observed that introducing water-based gels in the displacement process reduces permeability to water more than to oil. We focus our fluid analysis capability on the effects gelation has on residual fluid and we provide supporting evidence for the involvement of gel compaction (dehydration) and oil trapping, while discounting gel blockage in throats, as mechanisms contributing to this effect. Direct Numerical Simulation of Bubbly Flows and Interfacial Dynamics of Phase Transitions by Tianshi Lu Advisor: James Glimm August 2005 We studied the propagation of acoustic and shock waves in bubbly fluids using the Front Tracking hydrodynamic simulation code FronTier for axisymmetric flows. We compared the simulation results with the theoretical predictions and the experimental data. The method was applied to an engineering problem on the mitigation of cavitation erosion in the container of the Spallation Neutron Source liquid mercury target. The simulation of the pressure wave in the container and the subsequent analysis on the collapse of the cavitation bubbles confirmed the effectiveness of the non-condensable gas bubble injection method on reducing cavitation damage. Then we analyzed the interfacial dynamics of liquid-vapor phase transitions and the wave equations for immiscible thermal conductive fluids. The phase transition rate is associated by the kinetic theory with the deviation of the vapor pressure from the saturated pressure. Analytical solutions to the linearized equations have been explored. The adiabatic and the isothermal limits have been investigated for both the linearized and the nonlinear equations, for latter the method of travelling wave solutions has been used. The wave structure of the solution to the problem with Riemann data has been discussed. We also implemented a numerical scheme for solving the Euler equations with thermal conduction and phase transitions in the frame of front tracking. Heat conduction has been added to the interior state update with second order accuracy. Phase boundary propagation has been handled according to the interfacial dynamics. A numerical technique has been introduced to account for the thermal layer thinner than a grid cell. The scheme has been validated, extended to multi-dimension, adapted for cylindrical and spherical symmetry, and applied to the simulation of condensing and cavitating processes. Simulation of Magnetohydrodynamic Multiphase Flow By: Jian Du Advisor: James Glimm May 2007 Computational Magnetohydrodynamics(MHD), greatly inspired over the last decades by magnetic confinement fusion and astrophysics problems, has achieved significant results.However, the existence of moving free material interfaces with complex geometries in many important MHD problems creates major complications for numerical algorithms. We proposed a numerical algorithm for the study of Magnetohydrodynamics (MHD) of free surface flows at low magnetic Reynolds numbers. The numerical algorithm employs the method of front tracking and the Riemann problem for material interfaces, second order Godunov-type hyperbolic solvers, and the embedded boundary method for the elliptic problem in complex domains. The code is applied for the numerical simulations of free surface conductive liquid undergoing phase transitions. Examples include the simulations of mercury jet distortion in non-uniform magnetic field, liquid metal target for the Neutrino Factory/Muon Collider, and laser ablation plasma plume in the process of nanotube synthesis. Enhanced 3D Front Tracking Method with Locally Grid Based Interface Tracking By: Yuanhua Li Advisor: Xiaolin Li August 2007 We present a new interface tracking algorithm for 3D Front Tracking called Locally Grid Based tracking (LGB), which is demonstrated to be a significant improvement to the existing Front tracking method. It combines the best features of two previous 3D interface tracking algorithms. To be specific, it combines the robustness of Grid Based tracking with the accuracy of Grid Free tracking. We report the implementation of this algorithm and the comparison study with publicly distributed interface codes (the level set method), with published performance results (VOF and other methods) and with previous versions of front tracking. We also explore the application of this algorithm in the study of mean curvature flow and 3D chaotic fluid mixing problems. The Effects of Temperature Equilibrium in Mixed Cell Hydrodynamics By: Thomas Masser Advisor: John W. Grove December 2007 Eulerian numerical schemes work well for smooth single-component hydrodynamic flow modeling. For multicomponent flows in compressible regimes, the emergence of large gradients in flow variables and the presence of material contacts degrade the quality of traditional Eulerian finite difference methods. Although higher order Godunov methods adequately capture shocks, articial diffusion at contacts leads to poorly resolved interfaces. Many techniques, such as adaptive mesh refinement, interface preservers and explicit interface treatments, have been proposed to address these problems. We study two Eulerian codes that adopt different techniques to overcome the difficulties of modeling multicomponent compressible hydrodynamics. RAGE, a code developed for the U.S. Department of Energy's Advanced Simulations and Computing program, implements adaptive mesh refinement (AMR) to improve resolution of shocks and contacts. FronTier, a code developed at Stony Brook, resolves contacts and prevents numerical diffusion by tracking material interfaces and maintaining the jump in states across interfaces. The use of mixed cells distinguishes RAGE from FronTier. In RAGE, interfaces exist implicitly in mixed cells, where the code assumes pressure, temperature and velocity equilibrium among the multiple material components present, while FronTier maintains an explicit interface across material contacts. Although the equilibrium assumption of RAGE provides desirable mathematical properties and a relatively simple numerical implementation, the model depends on a microphysical assumption of that may not be valid for some shock interactions. By studying and comparing test problem results from FronTier and RAGE, we investigate how the mixed cell equilibrium assumption in influences flow hydrodynamics. In particular we consider several one and two- dimensional test problems. For the one-dimensional test problems and a two-dimensional test problem with similar materials, the two codes produce qualitatively similar results. However, in two-dimensional Richtmyer-Meshkov instabilities with strong shocks and high contrast materials, we find that the pressure-temperature-velocity equilibrium assumption leads to large scale differences in the inteface structure and relevant thermodynamic variables. Thermonuclear Flame Studies in Rectangular Geometry by Paul Allen LaVergne Advisor: Yongmin Zhang December 2007 We consider four aspects of the modeling of flame fronts associated with Type Ia supernova. In our study we 1) numerically approximate the value of the fire polishing wavelength, 2) show that for increased laminar flame speeds smaller surface areas result, 3) show that as fuel densities are decreased, a slower less stable flame with a larger surface area results, and 4) show that the turbulent flame speed becomes independent of the laminar flame speed. We discuss both front tracking and front capturing methods that have been used in similar studies and in particular, we discuss our implementation of a reactive front tracking model in a rectangular geometry for use in a attaining the above results.