Session L21: Computational Fluid Dynamics and Related Methods

Spatial Grand Canonical Monte Carlo Algorithms for Fluid Simulation

Strict detailed balance is essentially unnecessary for Markov chain Monte Carlo simulations to converge to equilibrium. Recently, we proposed a general Monte Carlo algorithm based on sequential updating that only satisfies the weaker balance condition. We have shown analytically that the new algorithm identifies the correct equilibrium distribution of states. Analysis of the diagonal elements of the transition matrices shows that the new algorithm is more mobile than the conventional Metropolis algorithm. Monte Carlo simulations of the Ising model and the lattice gas show that the new algorithm reduces autocorrelation time and thus improves the statistical quality of sampling. By exploiting the equivalence of the Ising model and the lattice gas, we demonstrate that the new method can also be applied to continuum systems, such as Lennard-Jones, in the grand canonical ensemble. Potential applications of the new algorithm are simulations of 1st and 2nd order phase transitions. Any massively parallel Monte Carlo simulation based on spatial decomposition involves simultaneous moves of atoms/molecules on multiple CPUs. Parallel Monte Carlo simulations based on spatial decomposition suffer from loss of precision due to the periodic switching of active domains. On the other hand, our algorithms are intrinsically sequential, and can be parallelized easily without compromising precision.   

Solving Parabolic Equations on a Random Grid Using a Generalized Finite Difference Method

A novel, generalized finite differencing scheme for solving parabolic initial value PDEs on a random grid will be described and results from its application to the diffusion equation will be presented. For a given number of points, $N,$ in a computational star, we parameterize the ($N-$6)-dimensional space of all possible approximations to the Laplacian that are accurate to first order. We have generalized von Neumann stability analysis to the random grid and we use simulated annealing to search parameter space to find a Laplacian that gives stable time evolution of the system. The creation of a stable Laplacian is moderately computationally intensive, but using it to evolve the PDE in time is of the same order as standard finite difference schemes on regular grids. We will present simulations using a Gaussian initial profile on a 10,000 point, annealed random grid in 2D. We also show that the same Laplacian also allows stable time evolution of the stiff, nonlinear, Cahn-Hilliard equation on the same grid. .   

A multiscale algorithm for pulled fronts in reaction-diffusion

I will present a numerical scheme for simulating front propagation in the discrete $A{\rightarrow}2A$ reaction-diffusion problem. The FKPP equation describes the continuum approximation to the particle model. However it has been suggested that the dynamics of the particle model converges very slowly to the continuum dynamics. As such, it is infeasible to probe the approach to continuum with a purely particle-based method. Instead, we have created a hybrid model which properly treats important fluctuations where needed.   

Numerical Model for Hydrovolcanic Explosions.

A hydrovolcanic explosion is generated by the interaction of hot magma with ground water. It is called Surtseyan after the 1963 explosive eruption off Iceland. The water flashes to steam and expands explosively. Liquid water becomes water gas at constant volume and generates pressures of about 3GPa. The Krakatoa hydrovolcanic explosion was modeled using the full Navier-Stokes AMR Eulerian compressible hydrodynamic code called SAGE [1] which includes the high pressure physics of explosions. The water in the hydrovolcanic explosion was described as liquid water heated by magma to 1100 K. The high temperature water is treated as an explosive with the hot liquid water going to water gas. The BKW [2] steady state detonation state has a peak pressure of 8.9 GPa, a propagation velocity of 5900 meters/sec and the water is compressed to 1.33 g/cc. \newline \newline [1] Numerical Modeling of Water Waves, Second Edition, Charles L. Mader, CRC Press 2004. \newline [2] Numerical Modeling of Explosions and Propellants, Charles L. Mader, CRC Press 1998.   

Efficient Nonlinear Atomization Model for Thin 3D Free Liquid Films

Reviewed is a nonlinear reduced-dimension thin-film model developed by the author and aimed at the prediction of spray formation from thin films such as those found in gas-turbine engines (e.g., prefilming air-blast atomizers), heavy-fuel-oil burners (e.g., rotary-cup atomizers) and in the paint industry (e.g., flat-fan atomizers). Various implementations of the model focusing on different model-aspects, i.e., effect of film geometry, surface tension, liquid viscosity, coupling with surrounding gas-phase flow, influence of long-range intermolecular forces during film rupture are reviewed together with a validation of the nonlinear wave propagation characteristics predicted by the model for inviscid planar films using a two-dimensional vortex- method. An extension and generalization of the current nonlinear film model for implementation into a commercial flow- solver is outlined.   

Discrete Bubble Modeling for Cavitation Bubbles

\textsc{Dynaflow, Inc.} has conducted extensive studies on non-spherical bubble dynamics and interactions with solid and free boundaries, vortical flow structures, and other bubbles. From these studies, emerged a simplified Surface Averaged Pressure (SAP) spherical bubble dynamics model and a Lagrangian bubble tracking scheme. In this SAP scheme, the pressure and velocity of the surrounding flow field are averaged on the bubble surface, and then used for the bubble motion and volume dynamics calculations. This model is implemented using the Fluent User Defined Function (UDF) as Discrete Bubble Model (DBM). The Bubble dynamics portion can be solved using an incompressible liquid modified Rayleigh-Plesset equation or a compressible liquid modified Gilmore equation. The Discrete Bubble Model is a very suitable tool for the studies on cavitation inception of foils and turbo machinery, bubble nuclei effects, noise from the bubbles, and can be used in many practical problems in industrial and naval applications associated with flows in pipes, jets, pumps, propellers, ships, and the ocean. Applications to propeller cavitation, wake signatures of waterjet propelled ships, bubble-wake interactions, modeling of cavitating jets, and bubble entrainments around a ship will be presented.   

Experimental and Numerical Studies on Multibubble Interaction

The behavior of multiple interacting bubbles can vary substantially from that of single bubble dynamics. We have conducted experiments and developed different levels of numerical tools for studying multiple bubble dynamics effects. In the experiments, multiple bubbles were generated simultaneously by spark and visualized using high speed video camera. The experimental results were then used to validate the numerical simulations obtained by 3\textsc{DynaFS}$^{\copyright }$\textsc{ }and\textsc{ PhantomCloud}$^{\copyright }$. 3\textsc{DynaFS}$^{\copyright }$\textsc{ }uses the boundary element method, which is capable of predicting non-spherical bubbled deformation as well as multibubble interaction. \textsc{PhantomCloud}$^{\copyright }$ is based on our earlier asymptotic expansion studies in which all bubbles are replaced with sources/sinks and dipoles whose intensities are determined by the underlying flow field and the presence of the other bubbles and their dynamics. The dynamics of each bubble is then recovered by a modified Rayleigh-Plesset equation and a center of the bubble equation of motion. In this approach, bubbles, which are punctual singularities, can interpenetrate without resulting in code failure, thus the ``phantom'' naming.   

Very-high temperature molecular dynamics of dense plasmas.

Finite temperature Orbital Free DFT coupled consistently with moleculer dynamics is applied to the hot-dense plasma regime up to 1000~eV and 100~$\mathrm{g\,cm^{-3}}$. Results obtained on dense iron and boron plasmas are compared with all-electron Quantum MD and effective classical theories like OCP and Yukawa OCP. A prescription on ionization for the classical model is made through the structural properties. Emphasis is also done on a comparison between Kubo-Greenwood and Ziman theories on dc conductivity in the very dense regime.   

Particle Beam Waist Location in Plasma Wakefield Acceleration

The role of beam waist location in interactions between a plasma and a particle beam is not yet fully understood. Nonlinear effects within the plasma make an analysis of such interactions difficult. I present five simulations in which I vary the waist location of a beam of ultra-relativistic electrons propagating through one meter of self-ionized lithium plasma. The simulation parameters are chosen to model the recent experiment 167 at the Stanford Linear Accelerator, relevant to the design of future plasma wakefield accelerating afterburners. I find that beams focused near the point of entry into the plasma propagate further into the plasma and accelerates witness particles to a greater maximum energy before disintegrating. These results could indicate that ion channel formation is dependent on the drive beam waist location and that the plasma accelerating medium can have an observable effect on the focusing of the drive beam.   

Distance of closest approach of two hard ellipses

The distance of closest approach of hard particles is a key parameter in their interaction and plays an important role in the resulting phase behavior. The distance of closest approach of the centers of hard spheres in 3-D or of hard circles in 2-D is the diameter. For non-spherical particles, the distance depends on orientation, and its calculation is surprisingly difficult. Although overlap criteria have been developed(1,2) for use in computer simulations, no analytic solutions have been obtained for ellipsoids in 3-D, or, until now, for ellipses in 2-D. We have succeeded in deriving an analytic expression for the distance of closest approach of the centers of two arbitrary hard ellipses as function of their orientation relative to the line joining their centers. We describe our method for solving this problem, illustrate our result, and discuss its usefulness in modeling and for simulating systems of anisometric particles such as liquid crystals.   

Discontinuous Molecular Dynamics for Rigid Bodies: Applications.

Event-driven molecular dynamics simulations are carried out on two rigid body systems which differ in the symmetry of their molecular mass distributions. First, simulations of methane in which the molecules interact via discontinuous potentials are compared with simulations in which the molecules interact through standard continuous Lennard-Jones potentials. It is shown that under similar conditions of temperature and pressure, the rigid discontinuous molecular dynamics method reproduces the essential dynamical and structural features found in continuous-potential simulations at both gas and liquid densities. Moreover, the discontinuous molecular dynamics approach is demonstrated to be between 3 to 100 times more efficient than the standard molecular dynamics method, depending on the specific conditions of the simulation. The rigid discontinuous molecular dynamics method is also applied to a discontinuous-potential model of a liquid composed of rigid benzene molecules, and equilibrium and dynamical properties are shown to be in qualitative agreement with more detailed continuous-potential models of benzene. The few qualitative differences in the angular dynamics of the two models are related to the relatively crude treatment of variations in the repulsive interactions as one benzene molecule rotates by another.\newline *This work was supported by grants from the National Sciences and Engineering Research Council of Canada (NSERC).   

Wave Phenomena during Lattice Boltzmann Simulation of Interdiffusion for Species having Unequal Masses

We discuss wave propagation phenomena in a Lattice Boltzmann (LB) model of a binary diffusion couple for species having unequal masses. LB simulations reveal oscillations in the position of the global center of mass of the couple as it moves toward the center of the couple from its initial position located toward the more massive species. These oscillations are related to waves in the total number density and barycentric velocity. Waves are generated at the initial discontinuity in composition and propagate toward the ends of the couple, from which they reflect. For a small difference in the masses of the diffusing species, we use a perturbation expansion to obtain driven wave equations for the total number density and the barycentric velocity. For sufficiently long samples these waves have a negligible effect on the composition versus distance profiles during interdiffusion; however, for microfluidic devices whose length is comparable to the diffusivity divided by the wave speed, the tails of the composition profiles get cut off. Periodic boundary conditions (PBC) were used in the direction perpendicular to the axis of the diffusion couple. If these PBC are replaced by no-slip walls, the waves are heavily damped.   

Potential-based Reduced Newton Algorithm for Nonlinear Multiphase Flow in Porous Media

We present a phase-based potential ordering for the finite volume discretization of the multiphase porous media flow equations. This ordering is an extension of the Cascade ordering introduced by Appleyard and Cheshire. The extension is valid for both two-phase and three-phase flow, and it can handle countercurrent flow due to gravity and/or capillarity. We show how this ordering can be used to reduce the nonlinear algebraic system that arises from the fully-implicit method (FIM) into one with only pressure dependence. The potential-based reduced Newton algorithm is then obtained by applying Newton's method to this reduced-order system. Numerical evidence shows that our potential-based reduced Newton solver is able to converge for time steps that are much larger than what the standard Newton's method can handle. In addition, when standard Newton does converge, the reduced Newton algorithm also converges, and often at a faster rate than standard Newton. Applications of the potential ordering to linear preconditioning will also be discussed.