3
short title
Research Signpost
37/661 (2), Fort P.O., Trivandrum-695 023, Kerala, India
Recent Res. Devel. Immunology, 3(2001): 145-159 ISBN: 81-7736-073-6
Modeling and Simulation of Turbulent MixingHyunkyung Lim, Yan Yu, James Glimm* and Xiaolin Li
Department of Applied Mathematics and Statistics, Stony Brook University, Stony Brook NY 11794-3600 USA
Thomas Masser, Baolian Cheng, John Grove, and David H. Sharp
Los Alamos National Laboratory, Los Alamos NM, 87545, USA
Abstract
Principal results of the authors and coworkers for the modeling and simulation of turbulent mixing are presented. We present simulations and theories which agree with experiment in several macroscopic flow variables for 3D Rayleigh-Taylor unstable flows. The simulations depend on advanced numerical techniques and on accurate physical modeling. Sensitivity to the choice of numerical method is demonstrated. Richtmyer-Meshkov unstable flows (in 2D) show dependence on physical modeling and on the numerical method for temperature and concentration probability density functions.
Correspondence/Reprint request: Prof. James Glimm, Department of Applied Mathematics and Statistics, Stony Brook University, Stony Brook NY 11794-3600, USA. E-mail: .
Introduction
Rayleigh-Taylor (RT) and Richtmyer-Meshkov (RM) unstable flows have been a source of continued interest [1:Sh84]. Our main results include theories and simulations for 3D RT unstable flows in agreement with experiment. We identify sensitivities in the simulations to numerical modeling (numerical mass diffusion and numerical surface tension) and to physical modeling (transport and surface tension effects and mixed phase thermal equilibrium). The sensitivities arise from a chaotic and nonconvergent interfacial area for the surface separating the two fluids. Such macro observables as the overall growth of the RT mixing zone show sensitive behavior. This same issue is explored more extensively in 2D RM simulations, where sensitivity to physical modeling and to numerical modeling show up in the temperature and species concentration probability density functions (PDFs). As observed in [2:Pitsch], these PDFs are a necessary input for the modeling of combustion.
2. Front tracking: control of numerical diffusion
Front tracking is a numerical algorithm in which a dynamically moving discontinuity surface is represented numerically as a lower dimensional grid, moving freely through a volume grid. Software to support this concept has been improved over a number of years, and is now highly robust, and openly available for downloading. With the use of this algorithm, finite difference expressions are not computed across discontinuity interfaces. Rather, the physically based jump conditions (the weak solution to the differential equations) provide the connection between the regions on each side of the interface. Not only does this algorithm eliminate numerical mass (or thermal) diffusion, but it can be modified to allow arbitrary amounts of (physical) mass diffusion [3:LiuLiGli05], a capability that is useful when the physical mass diffusion is small but nonzero.
Rigorous testing and comparison to alternate methods for interface handling have demonstrated the quality of the front tracking method [4:DuFixGli05].
From the point of view of numerical algorithms, we use a higher order Godunov scheme (MUSCL) and from the point of view of turbulence simulations, the method is implicit LES, or ILES, i.e. no explicit subgrid models have been employed.
3. 3D RT stimulations
We demonstrate sensitivity of the RT growth rate to transport (miscible fluids) and surface tension (immiscible fluids) and to the numerical analogues of these effects. With setting of the physical modeling parameters to their experimental values and with Front Tracking to control the numerical issues, we obtain agreement with experiment for the overall RT growth rate [5:GeoGli 06,6:LiuGeoBo06] . With and g denoting gravity, the self similar bubble side penetration distance h has the form , where is the dimensionless bubble side growth rate. See Table 1. Moreover, the simulations agree with other experimentally observed quantities: the bubble radius and the bubble height fluctuations. In the self similar flow regime, these quantities also have a growth rate proportional to , leading to the definition of an and an to characterize the growth of the bubble radii and the bubble height fluctuations.
Table 1. Comparison of RT growth rates.
Fluid type / Experiment / Simulation / TheoryMiscible: / 0.070 / 0.070
Immiscible: / 0.050-0.077 / 0.067 / 0.060
Immiscible: / 0.01 / 0.01 / 0.01
Immiscible: / 0.028 / 0.034 / 0.025
Simulation results of others are summarized in [7:alpha, 8:c&c]. Briefly, these investigators find a growth rate half or less of the experimental value. Uncertainties in the initial conditions has been mentioned as a possible cause of the discrepancy. Here we duplicate the discrepancy in our own untracked simulation [9:Li93], and in this simulation we show that numerical mass diffusion is the primary cause of the discrepancy, in the sense that if its effects are compensated for, even in the untracked simulations, experimentally consistent results are obtained [10:GeoGli05]. Numerical surface smoothing in the untracked simulations also plays a role in the discrepancy.
Sensitivity to both physical modeling and numerical effects is displayed in Fig. 1. The simulated value of varies by a factor of 25% with change of dimensionless surface tension, a factor of 3 with change of numerical method (tracked vs. untracked), and over a factor of 50% within the untracked simulations.
Figure 1. Bubble growth rate vs. dimensionless surface tension. Comparison among experiment and several simulations.
Comparison of tracked to untracked simulations (which disagree with experiments and with tracked simulations in the overall growth rate) clearly indicate that numerical mass diffusion in untracked simulations is a major contributor to the discrepancy between tracked and untracked simulations. To identify numerical mass diffusion as the leading cause of this discrepancy, we determine a time dependent Atwood numberfrom the simulation (tracked and untracked) itself [10:George]. Using , we redetermine and find consistency among all simulations and between simulation and experiment. See Figure 2.
The use of a time dependent Atwood number is a post processing step that can be added to any simulation. For this purpose, we outline the main steps in its definition. In each horizontal layer and for a fixed time step, we determine extreme densities or densities of extreme values of the two fluid concentrations. These densities serve to define a time and layer dependent Atwood number . For analysis of the bubble growth rate, we average over horizontal layers in the upper half of the bubble region, which in most cases is about the top quarter of the mixing layer. This average defines . Because is defined in terms of extreme, not average values, it gives mixing information different from the local mixing rate .
Presently missing from our simulations are two possibly canceling effects: the role of unobserved noise in the initial conditions and the role of subgrid models, to introduce additional transport from unresolved scales in an ILES simulation.
4. Theory to predict growth rates
4.1 Bubble merger models
Bubble merger models were first introduced by Sharp and Wheeler [11:SW]. They were given quantitative validity [12:GS] with the assumption of a bubble interaction velocity derived from the relative positions of
neighboring bubbles. These models are used to predict the RT bubble side mixing zone growth rate as well as the bubble width, both in agreement with experiment. As we formulate them [13:chaos], they have no free parameter. These models postulate an ensemble of bubbles of different heights and (due to bubble merger) an evolving mean radius. The excellent agreement of the model with experiment and simulation is displayed in Table 1. Note that agreement comprises three statistical measures of the mixing rate: growth rates for bubble height, bubble radii and bubble height fluctuations. Specifically, the bubble height to diameter ratio is given by in predictions of the model and in agreement with experiments.
Figure 2. Above: Two simulations with ideal physics, tracked and untracked, compared to one (tracked) with physical mass diffusion. Below: The same simulations, but compared using a time dependent Atwood number.
The two key inputs to the model are a formula for the velocity of each bubble and the criteria and method of bubble merger. The method of merger is the simplest. It is by preservation of cross sectional area among the pair of merging bubbles. The criteria for merger is a function of the velocities of the two bubble candidates for merger. When one of them acquires a negative velocity, it is merged into its (generally larger) neighbor.
This brings us to the formula for the bubble velocity. It is composed of two parts. The first is the single bubble velocity, i.e. the velocity of the bubble tip in a periodic array of bubbles of the same radius. This quantity has been extensively researched, and in dimensionless terms is known as the Froude number. It is independent of the Atwood number, and is well characterized, experimentally, theoretically and through simulation. The second component to the bubble velocity, the bubble interaction velocity, has to do with non-periodic perturbations of the ensemble of bubbles occupying the light fluid edge of the growing RT mixing layer. We assume a hexagonal array for the bubbles, in which each bubble is surrounded by six neighbor bubbles. The influence of each of these neighboring bubbles on the velocity of the central bubble is treated additively. For this purpose, we pick randomly from the ensemble a candidate bubble (i.e. a bubble height) to interact with a given bubble. For a chosen pair of bubbles with a given relative height displacement, we compute an interaction velocity (which can be positive or negative, relative to the basic single bubble motion). It is positive if the central bubble is more advanced (with a larger height), and otherwise negative. To compute this interaction velocity, we regard the interacting bubble pair as a bubble-spike combination at a larger length scale. This bubble-spike combination has a velocity determined by the previously mentioned single bubble theory, but now applied at intermediate times and not at the time of terminal velocity. The bubble interaction velocity, which is perhaps the most phenomenological aspect of this model, has been studied separately through experimental data analysis and simulations of bubble interactions [14:GliLi88,15:GliLiMen90].
The above conceptual definition of the bubble merger model allows the derivation of the formula
for the bubble growth rate in terms of the radius and height fluctuation growth rates. Through solution of the self-similarity conditions (the RNG fixed point), we obtain the values for these growth rates as in Table 1, and complete the definition of the model.
Here for hexagonal geometry bubbles [16:Abarazi]. Note is taken of uncertainties in the “post terminal” oscillations observed in the terminal velocity [17,18:Lin] and extrapolation of the formulas for bubble velocity to values of less than one. Also is a parameter relating the new to the old bubble radii after an area preserving bubble merger. The minor degree of uncertainty for this ratio arises from its weak dependence on fluctuations of the bubble radius, a phenomena not included in the model; for uniform bubbles, .
A related bubble merger model [19;svartz] is based on the conservation of probability for a distribution of bubbles of different sizes. The key input, the bubble merger rate, as a function of the merging bubble radii, is supplied presumably as in [20:alon] by a potential theory analysis of bubble merger. This model does not agree with experimental values for the bubble width to bubble height ratio, a fact that has been noted [21:Dimonte]. Both models predict values for close to experiment. Our model has a richer structure, in that the merger rate depends on the height separation as well as the bubble radii. We also use the knowledge of height fluctuations to advance the mixing zone edge by the velocity of the extreme bubble, not the average bubble velocity, in contrast to [19]. These differences presumably account for the difference in the conclusions of the models and the superior agreement our model achieves in comparison to experiment.
4.2 Buoyancy drag models
Buoyancy drag models predict bubble and spike side mixing zone growth rates for both RT and RM mixing, both asymptotic and transient. They have been considered by a number of authors [22,23:Alon 24:DS,25:D,26:DS]. We introduce the notation , or to denote the position of the edge of the mixing zone at which fluid vanishes, In this notation, bubbles () correspond to . In a self similar scaling regime, we expect
.
Here sign conventions are fixed by an assumption that the light fluid is above the heavy and . The buoyancy equations are ordinary differential equations for and the above are constraints on the choice of parameters in these equations.
We have favored an approach which minimizes the number of free parameters in the buoyancy drag equations. For , we write the equation derived from force balance among inertial, buoyancy and drag forces, with the inertial term corrected for added mass effects.
Here is the index complementary to .
Even this minimal equation has four undetermined parameters, the four values of the drag coefficient for and for RT and RM cases. Our main result is to reduce this number to one, for example the RT drag coefficient, or equivalently the RT bubble growth rate coefficient determined above from the bubble merger model. This goal is achieved for most values of , and a uniform theory valid for all values of requires one additional parameter. With these choices, we predict all available RT and RM bubble and spike data for the full range of allowed values of A, and achieve full agreement with experiment and with the theoretically required A = 1 spike values [27:Cheng, 28:Cheng-a].