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 Mixing Hyunkyung 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. 2 author name et al. Correspondence/Reprint request: Prof. James Glimm, Department of Applied Mathematics and Statistics, Stony Brook University, Stony Brook NY 11794-3600, USA. E-mail : email@example.com. 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 short title 3 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 A 2 1 / 2 1 and g denoting gravity, the self similar bubble side penetration distance h has the form h Agt , 2 where b 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 t 2 , leading to the definition of an r and an h to characterize the growth of the bubble radii and the bubble height fluctuations. Table 1. Comparison of RT growth rates. Fluid type Experiment Simulation Theory Miscible: b 0.070 0.070 Immiscible: b 0.050-0.077 0.067 0.060 Immiscible: r 0.01 0.01 0.01 Immiscible: h 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 b 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 b 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. 4 author name et al. 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 number A(t ) from the simulation (tracked and untracked) itself [10:George]. Using A(t ) , we redetermine b 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 A( z, t ) . For analysis of the bubble growth rate, we average A( z, t ) 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 A(t ) . Because A(t ) 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 short title 5 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 b / 2 r 3.0 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. 6 author name et al. 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 1 1 1 b cb r 2 1/ h 2 2k 2 for the bubble growth rate b in terms of the radius r and height fluctuation h 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. short title 7 Here cb vterminal / Agr 0.532 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 A 1 formulas for bubble velocity to values of A less than one. Also k 0.42 0.43 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, k 2 1 0.414 . 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 b 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 . 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 Z i (t ) , i 1, 2 or i b, s to denote the position of the edge of the mixing zone at which fluid i vanishes, In this notation, bubbles ( b ) correspond to i 1 . In a self similar scaling regime, we expect (1)i Z i (t ) i Agt 2 (RT); ( 1)i Z i (t ) i Agt i (RM) . Here sign conventions are fixed by an assumption that the light fluid is above the heavy and g 0 . The buoyancy equations are ordinary differential equations for Z i (t ) 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 Vi Z i dVi / dt , we write 8 author name et al. the equation derived from force balance among inertial, buoyancy and drag forces, with the inertial term corrected for added mass effects. dVi C V2 (1)i Agt i i ' i dt 1 2 Here i i 3 is the index complementary to i . ' Even this minimal equation has four undetermined parameters, the four values of the drag coefficient C i for i 1, 2 and for RT and RM cases. Our main result is to reduce this number to one, for example the RT drag coefficient C1 , or equivalently the RT bubble growth rate coefficient b determined above from the bubble merger model. This goal is achieved for most values of A , and a uniform theory valid for all values of A 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]. There are two steps to our reduction in the number of parameters. At the outset, we note that an elementary analysis of the large time asymptotes for the buoyancy drag equation link the RT drag coefficients and the RT growth rates. 1 i 2 1 Ci (1 (1)i A) The first step is a postulate that the RT values of the drag coefficients C i and the growth rate coefficients i apply to the RM data also. This postulate is justified by data analysis, and appears to have excellent agreement with the data. An asymptotic analysis of the RM solutions of the buoyancy drag equation leads to the formula 1 Ci (1 (1)i A) 1 i 1 A (1 (1)i A) This formula is a prediction of RM mixing rates in terms of RT mixing rates, as it gives the value of the observed quantity i in terms of the parameters C i , previously specified by the RT grow rates. The second main step is to relate the bubble and spike data, thereby eliminating one of the two remaining parameters. Here we introduce another hypothesis: that the center of mass Z COM (t ) of the mixing zone is stationary. short title 9 This statement is frame dependent and applies in the frame of the experimental container. The buoyancy equation itself is also frame dependent, and makes the same frame assumption. This hypothesis is satisfactory for A not too close to 1 , but it is not exactly correct, especially for A 1 . Thus, we introduce a refined assumption, that the COM satisfies self similar scaling, Z COM (t ) COM Agt . 2 A further analysis shows that s / b satisfies a quadratic equation whose coefficients depend on A and COM / s . Taking the limit A 1, we can evaluate s ( A 1) 0.5 (free fall) in the solution of this quadratic equation, and obtain 1 b 7 COM ( A 1) / s ( A 1) 6 s ( A 1) 60 b 0.05 . With the values at A 1 fixed, we further postulate for 7 COM A 60 where a fit to the LEM data [27:DS] requires 10 7 Figure 3. s / b as a function of A . Comparison of experiments [25:DS] and theory [27:cheng-a]. 10 author name et al. We compare the experimental data for s / b and as a function of A in Figures 3 and 4. In Figure 5, we compare our theory for COM to inferred data from the LES experiments [25:DS]. Figure 4. vs. A ; theory and experimental data [25:DS]. Figure 5. COM vs. A ; theory and experimental data. Note COM 0 for A 0.8 . short title 11 5. Averaged equations and mixture models 5.1 Introduction and discussion of mix models Averaged equations suppress many solution details and allow simplified solutions to problems through averaging over the details of the fine scale structure. A discussion of many proposed models and closures, with a listing of the number of adjustable parameters and other key features was given in [29:Cheng]. We propose mix models based on a complete first order multiphase closure. This means that all primitive variables, averaged over each phase separately, are retained in the averaged flow equations. Each phase has an independent pressure and complete thermodynamical variables. Mixed phase thermodynamics is not part of this model. The model has the pleasant feature of being hyperbolically stable, meaning that independently of the viscosity assumed in the equation, the time propagation has real characteristics and is stable. In contrast, a closure with fewer independent variables is called a reduced model. A commonly studied reduced model with equal phase pressures has complex characteristics and is stable only with assumption of a sufficiently large viscosity. The instability has been known for many years. Detailed studies, mathematical and numerical, of the instability have been conducted [30:Keyfitz], and show phase separation (of mixed phase flow regimes into separate unmixed regions of pure phase) as the instability mode. However, not all single pressure reduced models are hyperbolically unstable. We mention the Scannapieco-Cheng model , which is hyperbolically stable and has been successfully applied to the modeling of ICF experimental data. 5.2 Formulation Let X i denote the characteristic function of phase i , so that with averages given by angle brackets, i X i is the volume average of phase i . Density k X k / k and pressure pk are defined as volume averages while velocity vk X k vz / k , energy, internal energy and entropy, Ek , ek , S k are defined as mass weighted quantities. Multiplying the primitive (Euler equations of compressible fluids) equations by X i leads to the equations: k v* k 0 t z ( k k ) ( k k vk ) 0 t z 12 author name et al. ( k k vk ) ( k k vk vk ) ( k pk ) p * k k k g t z z z ( k k Ek ) ( k k vk Ek ) ( k pk vk ) ( pv)* k k k vk g t z z z Here we have introduced expressions v * , p * and ( pv)* still to be defined, and we have omitted the within-phase Reynolds and related heat flux tensors, as these appear to be small in the data we study. Extensions to higher dimensions (assuming a preferred normal direction ( z ) to the mixing layer) have been studied. Mathematical analysis [31,32:Jin] of the unclosed equations suggests a functional form for the closure terms as convex combinations of the corresponding pure phase quantities q* 1q q2 2 q q v, p ( pv)* p *( v v ) v *( p2 2pv p1 ) ( 1pv p2v2 2pv p2v1 ) 1 pv 2 pv 2 1 1 pv where k kq , q v, p, pv k d kq k ' These equations satisfy boundary conditions at the edge of the mixing zone, and are hyperbolic. Modifications for inclusion of transport in the primitive (unaveraged) equations have been derived. Conservation of mass, momentum and total energy follows from the form of the equations. Entropy is not conserved, even for flows which are isentropic before averaging. The reason for this is that averaging itself is not isentropic, as it modifies the entropy of mixing. But there is an entropy inequality for flows which are isentropic before averaging [32:Jin]. The identity d k d k ' 1 is required for the identity k k ' 1 , and q q q q p pv so the closure has three free parameters. Of these, two ( d k , d k ) play no role for RT and RM data, and so we set these parameters to unity. There remains one free parameter. The equations are missing one condition at each edge of the mixing zone. A count of the number of characteristics shows that at those locations, there is a missing characteristic, associated with the incoming sound wave of the fluid missing in the single phase region. Thus we couple the model to the buoyancy drag equations for the mixing zone edges, to give an extra condition there. A mathematical analysis for the RT case v shows the connection of d k to the edge motions, with short title 13 Vk ' d kv Vk v in the incompressible limit. For RM flows, d k is not important in the model and is conveniently set to one. The equations, in the 1D RT incompressible case, admit a closed form solution [34,35,36:saltz] The closed form solution was extended to weakly compressible flows in a singular perturbation expansion [37:Jin]. The extension of these ideas to multiple phases was given in [38:Cheng]. 5.3 Numerical studies of mix models We solve numerically the equations for the complete first order closure [39:Jin]. In this study, we compare the solutions to the theory derived above, and we solve a sample higher dimensional flow problem. Using the RT and RM two fluid simulation data sets discussed here, we compare [40:Bo] closure model residual errors between our closure and that of Abgrall and Saurel. In this comparison, our errors are about two times smaller than theirs. The average relative error for the three closure terms for our model is 12% (closure residuals compared to RT simulation data) and 9% (RM), while for the Abrall-Saurel [41:Saurel] these same relative errors are 30% (RT) and 13% (RM). See Table 2. The Abrall-Saurel closure model errors are computed assuming no relaxation terms. If their relaxation terms are included their errors are larger. See Figure 6. Figure 6. Plot of v * (left) and ( pv)* (right) times d / dz as an exact expression computed from the 3D RT FronTier simulation data and closed expressions to define v * and ( pv)* approximately. Our own closure is compared to two interpretations of the Abgrall-Saurel closure, with and without their relaxation terms. 14 author name et al. Figure 7. Density plot for circular RM simulation (FronTier), shown at late time. Table 2. L1 norm relative errors for closure terms as compared to simulation data. v* p* (pv)* Average RT data comparison 18% 0% 18% 12% RM data comparison 7% 0% 20% 9% 6. RM code comparisons Code comparison of RAGE to FronTier was conducted for a circular Richtmyer-Meshkov problem in 2D. This model problem consists of an incoming circular shock wave crossing a perturbed circular interface, reflecting from the origin and reshocking the now strongly perturbed interface and mixing zone. The inner and outer materials are described by a stiffened gamma law gas with parameters chosen to model either lucite or air as the interior fluid and tin as the exterior fluid. For the purpose of this comparison, all transport coefficients have been set to zero. The problem is defined in detail in [42:Yu,43:masser] Convergence under mesh refinement for the FronTier code was analyzed in [42:Yu]. The solution develops new structures at each new length scale under mesh refinement, and so convergence was considered in the sense of local spatial averages. For most variables studied, a small interval of angular averaging was sufficient to allow convergence. Density contours showing the highly chaotic mixing zone at late time are presented in Figure 7. The comparison study shows approximate agreement for such macro solution features as the locations of the lead shock waves and mixing zone edges. But micro solution features such as the temperature and concentration PDFs are not comparable, and systematic mesh refinement does not appear to remove these discrepancies. It appears that ILES simulations (e.g those with short title 15 zero or under resolved regularizing transport and no explicit subgrid model) are sensitive to numerical transport coefficients and models for mixed cell thermodynamics [43:Masser07]. The differences remain even after some level of physical mass diffusion is allowed in the FronTier simulation. To illustrate the problems created by mixed cell thermodynamics, we present in Table 3 the width, in mesh units, for the mass and temperature discontinuities after reshock, for an unperturbed (radially symmetric 1D) simulation of the same physics. Use of an EOS to model an air-tin rather than a lucite-tin problem leads to significantly larger diffusion widths . Table 3. Numerical mass and thermal diffusion layer thickness in mesh units after reshock for a 1D radially symmetric RAGE simulation. Lucite-Tin mass diffusion layer 5r thermal diffusion layer 63r In Figure 8, we show a side by side comparison of the RAGE (left) and FronTier (right) simulations at a time shortly after reshock, contrasting the temperature fields at late time. Figure 8. RAGE (left) and FronTier (right). Temperature fields compared at a time following reshock for a tin- lucite simulation. We show in Fig. 9 the RAGE light fluid concentration PDF for the ideal tin- lucite simulation. The comparable FronTier simulation has a PDF delta function showing zero fluid mixture, in accordance to the physical modeling of the simulation. More comparable are the FronTier simulations of Sect. 7 with finite Reynolds and Schmidt numbers. The RAGE PDF is roughly intermediate between our Sc 1and Sc 0.1tracked simulations, a fact consistent with 16 author name et al. the mass diffusion width of the RAGE calculation, suggesting a numerical Schmidt number of perhaps 0.3. Figure 9. RAGE PDF for light fluid mass concentration, taken along a mid line of the mixing zone. 7. RM sensitivity to physical modeling 7.1 Introduction Under mesh refinememt, the interface (50% concentration isosurface) between the two fluids fails to converge, and its length is approximately 1 proportional to x after reshock. We show in Figure 10 the ratio of the interface length to mixing zone area in mesh units, i.e. with this divergent 1 factor of x scaled out. The flow is highly chaotic. On this basis, we postulate a simple model for simulation error, Error C1 x Interface Length C1C2 for the unregularized simulation, i.e. for the simulation with zero transport coefficients. Figure 10. The interface, after reshock, occupies about 20-30% of the available grid blocks, uniformly under mesh refinement, for thunregularized simulation. short title 17 The C i are order 1, apparently mesh independent quantities. This formula raises questions for the meaning of the simulation, especially for interface sensitive quantities, such as the degree of atomic scale fluid mixing. It also raises the possibility that for quantities depending on physical transport, such as Reynolds and Schmidt numbers, there is a possibility of apparent convergence to a nonunique solution selected in a code dependent manner by numerical Reynolds and Schmidt numbers. 7.2 The simulation study To investigate grid convergence and sensitivity to physical transport properties, we studied the same 2D RM problem with non-zero transport coefficients. To determine the dependence on these coefficients, we allowed them to vary over wide ranges, as needed to illustrate a range of LES and DNS simulations, and to separate numerical from physical transport [44:Lim]. See Figure 11 for parameters describing the simulations. For the separation of numerical from physical Schmidt number, the FronTier code is especially convenient, as it can simulate any Schmidt number through the controlled interface diffusion algorithm, without the expense of additional mesh refinement. Figure 11. The Reynolds number and Kolmogorov scale K mesh K / x vs. mesh Reynolds number. 7.3 The DNS-LES transition We divide the simulations into those in an LES regime and those in a DNS regime. The division between the two is not sharp, but is characterized roughly by the condition Remesh 1 or equivalently Kmesh K / x 1 , where K is the Kolmogorov length scale. We find signs of interface convergence in the DNS regime but not in the LES regime. To show these 18 author name et al. trends, we plot the mesh length to area ratio, at a fixed time t 90 , after reshock and at the beginning of the chaotic regime. In Figure 12, we see that convergence of the interface (as measured by a decrease in the mesh interface length to area ratio), depends primarily on the mesh Reynolds number or equivalently on the mesh Kolmogorov length, and is otherwise largely independent of mesh and Schmidt number. To understand numerical convergence for the DNS regime, we replot the same data ( Sc 0.1 only) in Figure 13, interface fraction (length/area) in physical units vs. Re . By this measure, the interface has converged in the DNS regime. Higher Sc numbers give similar behavior, but delayed (slower) convergence. 7.4 Physical and numerical mass diffusion Starting with a mesh Reynolds number of about 1, it appears that the interface stops growing as the mesh is further refined. Imagining that we are dealing with a mesh converged interface, we wish to understand the degree of mass diffusion to expect in the simulation. Intuitively, this is measured by the ratio of two length scales. Figure 12. The mesh interface length to area ratio (see Fig. 11) at t 90 , i.e. at the beginning of the chaotic regime, for simulations conducted with a variety of grid levels, mesh Reynolds numbers and Schmidt numbers. short title 19 Figure 13. Replot of data from Figure 11, using physical, not mesh units. The first is a correlation length, for the bubble size distribution, and the second is a diffusion length, for the diffusion distance achieved in the time which has elapsed. To define a correlation length, we model the probability a of minimum distance to exit a given phase, given a random starting point, by a Poisson distribution. Thus the PDF of distance to exit is 1exp( x / ) and here d(Probability to exit)/dx . Also corr is the correlation 1 length for this exit distance. is determined by a fit to simulation data. In Fig. 14, we show corr as a function of Re mesh , to indicate the degree of mesh convergence achieved for corr and thus for the ratio diff / corr . Figure 14. Plot of corr vs. Remesh for Sc 10,1.0,0.1 . The diffusion distance is conveniently defined as diff (t t0 ) where is coefficient of mass diffusion and t0 is the time for creation of an element of interface length. Since the interface length grows continuously with time, there is no precise way to set t0 , but in view of the strong increase in interface length shortly after reshock, we take t0 to be the time of reshock. A ratio diff / corr 1 , which is satisfied for all LES cases ( Sc 0.1 is the smallest value we consider). For DNS cases, with Sc 10 , are not diffusive, while the opposite case, DNS with Sc 1,0.1 , is diffusive. The diffusion length scale diff has converged on all meshes considered here. 7.5 Concentration PDF’s A principal finding is that concentration and temperature PDFs show dependence on the details of physical values for transport coefficients. The 20 author name et al. goal of a predictive simulation, of course, is to compute independently of the mesh, but without the prohibitive expense of DNS. According to [2:Pitsch], subgrid models will compensate for the missing scales, as far as viscosity is concerned, while the atomic level mass diffusion is fit to a convenient distribution (the beta-distribution) with coefficients determined from experimental data. For the present purposes, this prescription, especially the dependence on experimental data, is outside of the objectives of the computation. Figure15. Light fluid concentration PDF for DNS simulations and Sc 10,1.0,0.1 . According to , the subgrid model is defined in terms of filtered variables with the filter length larger than the grid spacing. We approach this latter objective indirectly, replacing the filter and its length scale larger than the mesh with large (artificial) coefficient in a physical viscosity term, to turn an LES into a (numerical) DNS. According to Figs. 12 and 13, this will cut off the formation of small interface structures under continued mesh refinement. We hope to obtain concentration PDFs by using the correct Schmidt number for the simulation. The general shape of the concentration PDFs appear to be short title 21 largely determined by the ratio diff / corr . To force this ratio to its converged value, we want to preserve the physical Schmidt number for the simulation, rather than the physical value mass diffusion directly. 7.6 Convergence of macro solution features The convergence of principal solution waves (shock and mixing zone edges) is adversely affected by the high level of viscous damping in the DNS simulations performed at moderate grid resolution. To understand this degradation of solution accuracy, we plot the relative position error vs. Re mesh in Figure 16. The relative error for each principal wave is defined as computed - exact / exact . Here is a time L1 norm; and exact is interpreted to be the result of the finest grid computed with zero transport terms. The total wave error results from averaging over the shock, mixing zone center and mixing zone width. Figure 16. Plot of relative wave error vs. Re mesh It is evident from comparison of Figures 13 and 16 that we can achieve convergence with either the micro or the macro observables, but we have not yet achieved both in the same calculation. We can predict the level of mesh needed to follow this simulation strategy and achieve convergence for both in the same simulation. It appears that this can be accomplished with hardware currently existing. But we observe that this plan involves subgrid modification of the momentum equation only with a numerically chosen viscosity, while the mass diffusivity is chosen to give a physical Schmidt number. For a more practical LES simulation, we regard both the viscosity and mass diffusion coefficients as defining a subgrid term, for the momentum and concentration equations respectively. For this strategy, we determine the ratio diff / corr from a simulation that has converged in its micro variables, as in Figure 13, and force the ratio to have the same value in a macro-converged LES simulation by adjustment (as a mass diffusion subgrid model) of the coefficient of mass diffusivity. We note that this proposed LES algorithm depends on a unique capability of front tracking: to simulate arbitrary (even very small) levels of mass diffusion. For many codes, the combined physical and subgrid levels of mass diffusion (especially for large Schmidt numer) are less than the inherent numerical mass diffusion, and so the algorithm would not be feasible, or would require additional levels of mesh refinement. 8. Discussion, comments and conclusions The good news is that many macroscopic flow observables for RM mixing appear to be insensitive to details of numerical and physical modeling. The bad news is that some important RM observables, including temperature and 22 author name et al. concentration PDFs do appear to be sensitive to numerical and physical modeling details. For RT mixing, it appears that many flow observables, including the much studied overall growth rate, are sensitive to numerical and physical modeling details. Sensitivity is generally traced to numerical mass diffusion, exacerbated in the RM case by a temperature equilibrium hypothesis for mixed states. Both numerical mass diffusion and mixed state thermal equilibrium are common features of untracked simulations. Because these mixing flows are chaotic, and have a divergent interfacial area, apparently 1 proportional to x in the absence of physical regularization, errors associated with the interface, of necessity unavoidable especially for untracked simulations, run the risk of nonconvergence. Where these simulations are regularized by grids and numerical code features, there is a risk of an apparent convergence to code dependent values. In the language of turbulent modeling, this feature of LES has been systematized through the use of filtered equations, and LES that converge to the solution of the filtered, but not the original equations. Solutions are assessed for convergence in both macro and micro variables. While we have not achieved acceptable convergence for both types of variables simultaneously, we have developed a methodology which allows prediction of mesh levels, numerical algorithms and physical or subgrid modeling procedures which in combination will achieve this goal. With subgrid modification of both the momentum and concentration equations, these algorithms are expected to yield convergence results for both macro and micro variables with feasible grids. To obtain practical LES simulations for quantities sensitive to numerical and physical modeling, the classical ideas of turbulence modeling propose filtered equations with a filter length larger than the mesh length. These equations allow convergence under mesh refinement, but to a solution of the filtered, not the original equations. In the spirit of LES, convergence of the filtered equations to the true ones is achieved (for statistically observable quantities only) by decreasing the filter length (and the mesh in proportion). We propose here to deviate from this program only in the use of an artificial value for physical viscosity in place of the filter. The purpose of this step is to control the growth of interface length under (further) mesh refinement. An alternate route would be the use of surface tension, with a numerically chosen (artificial) non-physical coefficient. With the interface thus stabilized, we then propose that the mass diffusion should be set to preserve physical values for corr / diff . It appears that the use of physical values for the Schmidt number will accomplish this goal. In any case, it is necessary to replace the step of calibration from experimental data. This proposal, of course, requires exploration beyond the information presented here. short title 23 The agreement of front tracking RT simulations with experiment for several growth rate measures (bubble height, width, height fluctuations) does not close the RT simulation issues. There remain initial condition modeling uncertainties and possible effects of missing subgrid models. Since these effects will move the growth rate in opposite directions, the possibility of canceling errors is not yet excluded. Acknowledgment This work is supported in part by the U.S. Department of Energy and the Army Research Organization. References 1. Sharp, D.H. 1984, Physica D 12, 3. 2. Pitsch, H. 2006, Ann. Rev. Fluid Mech. 38, 453. 3. Liu, X. F., Li, Y. H., Glimm, J. and Li, X. L. 2007, J. Comp. Phys. In Press. 4. Du, J., Fix, B., Glimm, J., Jia, X., Li, X., Li, Y. and Wu, L. 2006, J. Comp. Phys. 213, 613. 5. George, E., Glimm, J, Li, X. L., Li, Y. H. and Liu, X. F. 2006, Phys. Rev. E. 73, 016304. 6. Liu, X. F., George, E. Bo, W. and Glimm, J. 2006, Phys. Rev. E 73, 056301. 7. Dimonte, G,, Youngs, D. L., Dimits, A., Weber, S., Marinak, M., Wunsch, S., Garsi, C., Robinson, A., Nadrews, M., Ramaprabhu, P., Calder, A., Fryxell, B., Bielle, J., Dursi, L., MacNiece, P., Olson, K., Ricker, P., Rosner, R., Timmes, F., Tubo, H., Young, Y.-N. and Zingale, M. 2004, Phys. Fluids 16, 1668. 8. Cabot, W. and Cook, A. 2006, Nature Physics 2, 562. 9. Li, X. L. 1993, Phys. Fluids 5, 1904. 10. George, E. and Glimm, J. 2005, Phys. Fluids 17, 054101. 11. Sharp, D.H. and Wheeler, J 1961, Institute of Defense Analysis technical report, AIP Document No. PAPSFLDA-31447-50. 12. Glimm, J. and Sharp, D.H. 1990, Phys. Rev. Lett. 64, 2137. 13. Cheng, B., Glimm, J. and Sharp, D. H. 2002, Chaos, 12, 267. 14. Glimm, J and Li, X.L. 1998, Phys. Fluids 31, 2077. 15. Glimm, J., Li, X.L., Menikoff, R., Sharp, D. H. and Zhang, Q. 1990, Phys. Fluids A, 2, 2046. 16. Abarzhi, S.I. 1998, Phys. Rev.Lett. 81, 337. 17. Lin, A. 2000, State University of New York at Stony Brook Thesis. 18. Glimm, J. Li, X.L., Lin, A. D. 2002, Acta Mathematicae Applicatae Sinica 18, 1. 19. Oron, D., Sadot, O., Srebo, Y., Rikanati, A., Yedvab, Y., Alon, U. Erez, L. Erez, G. Ben-Dor, G., Levin, L.A., Ofer, D. and Shvarts, D. 1999, Laser Part. Beams 17, 465. 20. Alon, U, Hecht, J. Ofer, D and Shvarts, D. 1995, Phys. Rev. Lett. 74, 535.. 21. Dimonte, G. 2000, Phys. Plasmas, 7, 2255. 22. Alon, U., Hecht, J., Mukamel, D. and D. Shvarts 1994, Phys. Rev. Lett. 72, 2867. 23. Alon, U., Hecht, J., Ofer, D. and D. Shvarts 1995, Phys. Rev. Lett. 74, 534. 24. Dimonte, G. and Schneider, M. 1996, Phys. Rev. E 54, 3740. 25. Dimonte, G. 1999, Phys. Plasmas, 6, 2009. 24 author name et al. 26. Dimonte, G. and Schneider, M. 2000, Phys. Fluids 12, 304. 27. Cheng, B., Glimm, J. and Sharp, D. H. 2000, Phys. Lett. A 268, 366. 28. Cheng, B., Glimm, J. and Sharp, D. H. 2002, Phys. Rev. E, 66, 036312. 29. Cheng, B. 2007, In: Proceedings of the 10th International Workshop on the Physics of Compressible Turbulent Mixing, Ed M. Legrand and M. Vandenboomgaerde. Commisariat a l’Energie Atomique, Bruyerer-le-Chatel, France. 30. Keyfitz, B. 1991, SIAM J. Math. Anal. 22, 1284. 31. Scannapieco, A. and Cheng, B. 2002, Phys Lett. A 299, 49. 32. Jin, H. Glimm, J. and Sharp, H. 2006, Phys. Lett. A 353, 469. 33. Jin, H. Glimm, J. and Sharp, H. 2006, Phys. Lett. A 360, 114. 34. Glimm, J., Saltz, D. and Sharp, D. H. 1998, In: Nonlinear Partial Differential Equations, Ed. Chen, G. Q., Li, Y. and Zhu, X. World Scientific, Singapore, 124. 35. Glimm, J., Saltz, D. and Sharp, D. H. 1998, Phys. Rev. Lett. 80, 712, 36. Glimm, J., Saltz, D. and Sharp, D. H. 1999, J. Fluid Mech. 378, 119. 37. Glimm, J. and Jin, H. 2001, Bol. Soc. Bras. Mat. 32, 213. 38. Cheng, B., Glimm, J. and Sharp, D. H. 2005, Phys. Fluids 17, 087102. 39. Glimm, J. Jin, H. Laforest, M., Tangerman, F. and Zhang, Y. 2003, Multiscale Model. Simul. 1 458. 40. Bo, W., Jin, H., Kim, D., Liu, X., Lee, H., Pestieau, N., Yu, Y. Glimm, J. and Grove, J. 2007, Stony Brook University Technical Report. 41. Saurel, A. and Abgrall, R. 1999, J. Comput. Phys. 150, 425. 42. Yu, Y, Zhao, M., Lee, T., Pestieau, N., Bo, W., Glimm, J. and Grove, J.W. 2006, J. Comp. Phys. 217, 200. 43. Masser, T. 2007, Stony Brook University Ph.D. Thesis. 44. Lim, H., Yu, Y., Jin, H., Kim, D., Lee, H. Glimm, J. and D. H. Sharp 2007, Stony Brook University Technical Report.