Form 836 (7/06) Development of a Sub-Scale Dynamics Model for Pressure Relaxation of Multi-Material Cells in Lagrangian Hydrodynamics

We have extended the Sub-Scale Dynamics (SSD) closure model for multi-fluid computational cells. Volume exchange between two materials is based on the interface area and a notional interface translation velocity, which is derived from a linearized Riemann solution. We have extended the model to cells with any number of materials, computing pressure-difference-driven volume and energy exchange as the algebraic sum of pairwise interactions. In multiple dimensions, we rely on interface reconstruction to provide interface areas and orientations, and centroids of material polygons. In order to prevent unphysically large or unmanageably small material volumes, we have used a flux-corrected transport (FCT) approach to limit the pressure-driven part of the volume exchange. We describe the implementation of this model in two dimensions in the FLAG hydrodynamics code. We also report on Lagrangian test calculations, comparing them with others made using a mixed-zone closure model due to Tipton, and with corresponding calculations made with only single-material cells. We find that in some cases, the SSD model more accurately predicts the state of material in mixed cells. By comparing the algebraic forms of both models, we identify similar dependencies on state and dynamical variables, and propose explanations for the apparent higher fidelity of the SSD model.


Introduction
Multimaterial ALE and Lagrange calculations may need to account for mixed cells, which contain multiple pure materials meeting at one or more interfaces within the cell.In this case a closure model is required to partition the cell volume and internal energy.Absent such a model, the variables associated with the mesh are generally insufficient to determine that partition uniquely.
We have extended the Subscale Dynamics (SSD) closure model [1,2] to multiple dimensions and enabled it to handle any number of materials in a cell.We have implemented it in one and two dimensions in the Lagrange/ALE hydrocode FLAG [3].We compare below the performance of the SSD model with a well-known closure model due to Tipton [4][5][6][7], as also implemented in FLAG.Finding considerable similarity as well as some significant differences between the two, we discuss the reasons for this in terms of the equations defining both models.

Notation
Let V be the volume of a single computational cell.Then within that cell, material i has volume V i , volume fraction f i = V i /V, mass density ρ i , pressure p i , sound speed c i and a e-mail: alanh@lanl.govcompliance (bulk modulus) B i = ρ i c 2 i .Superscripts denote instants or intervals within the timestep from t n to t n+1 ; to wit, superscript 0 denotes time t n , + is t n+1/2 , 1 is t n+1 ; and a and b refer to the intervals from t n to t n+1/2 and from t n to t n+1 , respectively.

Sub-Scale Dynamics (SSD) model
The sub-scale dynamics (SSD) model, inspired by the work of Delov and Sadchikov [8], Goncharov and Yanilkin [9] and Barlow [10], estimates the material volume changes based on interface motion, using a Riemann solution for velocity.We have extended the model to cells with any number of materials, computing the pressure-differencedriven volume and energy change of each material as the algebraic sum of pairwise interactions.In multiple dimensions, we use interface reconstruction to provide interface areas.A few of the swept volumes are depicted as rectangles in Figure 1, where S ik is the area of the interface between materials i and k.The corresponding interface velocity (directed from i to k) is so the swept volumes are Then the new material volumes at t n+1/2 (predictor step) and t n+1 (corrector) result from changes due to dilation (or compression) of the entire zone, and increments resulting from pressure differences between materials.The limiters C {a,b} ik ∈ [0, 1] multiplying the swept volumes in the last expression will be defined below.
During the corrector step, internal energies are also updated by adding work terms in terms of the Riemann pressure in which we define an averaging weight W i = 1/(ρ i c i ) and a normal velocity difference ∆u ik = (u i − u k ) • nik , where u i is the velocity at the centroid of material i, and nik is the unit normal to the interface, directed from i to k.
It is necessary to limit material volume changes; physics requires that V + i ≤ V + , and for numerical stability we desire that f + i ≥ f 0 i for some number ∈ [0, 1] (taken equal to 0.25 in our calculations).This is achieved by computing the limiters in equation ( 5) in an FCT-like way.For the predictor step, The corrector-step limiters C b ik are computed in the same way, replacing

Tipton model
A simple heuristic to motivate the Tipton closure model is that if we linearize equations of state as we can solve for the equilibrium pressure and volume fractions in the cell.In order to prevent instantaneous relaxation to pressure equilibrium, the model uses a rescaled compliance in which the second parenthesized term is motivated by artificial-viscosity-like considerations, and L is a characteristic length scale (width) of the computational cell.For this model is it useful to define the averaging weight w i = f i /D i and mean rescaled compliance and pressure Then the predictor step (as implemented in FLAG) is where α is a stability parameter ≤ 1, and the result is limited so that |∆ f a i | ≤ 0.25 f 0 i .The corrector step could be done by computing ∆ f b i as in (15), based on the updated material state p + i , D + i .However, for efficiency reasons the implementation in FLAG simply approximates ∆ f b i = 2∆ f a i .Internal energy is updated by adding the work done by a zonal-average effective pressure at t n+1/2 Note that Tipton's model uses no subcell geometrical information.
3 Test Problems

Expanding Bubble
We have tested the SSD model in the FLAG hydrocode (which also contains the Tipton model).Figure 2, colored by specific internal energy, shows 2D Lagrange calculations of an expanding gas bubble surrounded by lowerpressure gas.We made SSD and Tipton model calculations on a regular rectangular mesh with premixed cells at the bubble surface, while the pie mesh pertains to a calculation with no mixed cells.We see that the SSD model predictions are similar to those of the unmixed calculation, and possibly a little better than those of the Tipton model.

Test with a Three-material Zone
We would expect the SSD model to treat cells containing three or more materials better than the Tipton model does; in such a cell, there are so many possible geometrical arrangements of materials and interfaces that a model that lacks subcell geometrical information has little hope of doing the right thing. Figure 3 shows a physically onedimensional shock tube calculation, in which the left material is set up as two different materials with identical properties and initial states.The problem includes a threematerial cell and a large number of two-material cells.A short time into the calculation, we see that the SSD model (Figure 4) has allowed all the vertical interfaces to advance to the right by the same distance, thus appropriately maintaining planar symmetry.The Tipton model (in Figure 5) breaks that symmetry.For comparision, we have also run the same problem without multimaterial cells (Figure 6).

Two-material Shock Tube Problems
We also investigated how the two closure models compare in two-material cells.For clarity, we used Lagrangian calculations of physically one-dimensional problems to investigate this.One setup was the Sod shock tube, containing gamma-law gases with gamma values and initial conditions as given in Table 1.In another setup, we replaced the low-pressure gas with a Grüneisen-law fluid [11] with properties typical of copper; details are shown in Table 2.
(We shall refer to this problem as the "Cu" shock tube   problem.)Both shock tube problems were run in Lagrange mode in three different ways-with no mixed zones and no closure model ("clean"), and with mixed zones at the interface, with both closure models.The initial mesh for the Sod problems is shown in Figure 7.
The Sod shock tube calculations with the Tipton closure model became unstable and halted with a tangled mesh.To improve stability, model parameter α was decreased from 1 to 0.5 and finally to 0.25 before the prob-Table 1.Initial conditions for Sod shock tube problem (CGS units).
Table 2. Initial conditions for shock tube problem with Grüneisen-law fluid (CGS units).
γ-Law Gas Grüneisen-Law Fluid    lem would run to completion.In the "Cu" shock tube runs, it was not necessary to decrease α.
Time histories of pressure, density and specific energy for the Sod problem are shown in figures 8-10.For the clean calculations, the state of the zones adjoining the interface are plotted.For the mixed calculations, the state of each material in a mixed zone are shown.(In all calculations that ran to completion, the four zones on one side of the interface, or straddling it for mixed runs, had identical histories.)We see from Figure 8 that all completed calculations attained the same equilibrium pressure, but the   approach to equilibrium differed.Regarding the clean calculation as likely the most accurate, the SSD model was more accurate than the Tipton model, and the latter became less accurate as α was decreased to attain stability.Similar conclusions follow from the results for density and energy, although the material 2 equilibrium density and energy may be a little better in the Tipton model calculation.
Figures 11-13 show the results of the "Cu" shock tube problem.As before, the SSD model is closer than Tipton's model to the "clean" results.For this shock tube, the Tipton model does not get the right density and energy even in the late-time limit.We speculate that this poorer performance (compared with the Sod problem) may be due to the fact that the Grüneisen-law pressure is not proportional to density [11].In fact, Tipton's model predicts a mechanically impossible pressure in the Grüneisen material in the two-material zone.We can see this in Figure 14, which shows mesh plots of the central region of each of the three "Cu" shock tube calculations, colored by pressure.Figure 15 shows pressure as a function of time in both sections of the two-material cell and in the two single-material cells to the left and right of that cell (or, for the clean calculation, in the three cells on each side of the material interface).We see that the Grüneisen material in a two-material zone maintains for a considerable period a lower pressure than either the γ-law gas to its left in the same cell or the pure Grüneisen material in the cell on its right-even though it started with a pressure intermediate between those two neighboring pressures.This is mechanically impossible.However, the SSD model does not exhibit this anomaly; it matches the clean calculation well.

Algebraic Comparison of Closure Models
We would like to understand the differences in behavior between the Tipton and SSD models-but upon comparing model philosophies, let alone equations, it is not clear why the two models should have any similarities at all.Tipton's model is based on the pressure equilibrium state, and deliberately relaxes toward it, while the SSD model has no notion of the pressure equilibrium state.The Tipton model depends explicitly on the preexisting material volume fractions within each cell, while the SSD model does not even refer to them.On the other hand, the SSD model uses interface areas and orientations and the velocities of individual materials (by interpolating mesh velocity to the centroid of each material polygon), and estimates interface velocities.The Tipton model knows nothing of material velocities or the geometry or motion of the interfaces.How can such dissimilar models behave similarly at all?

Volume Relaxation
In order to facilitate algebraic comparison of the models, we note that the Tipton model volume change can be expressed in the form (3) if we approximate ∆V b = 2∆V a and neglect ∆V b ∆ f b i .Then the zone-dilation and pressuredriven terms are Fig. 15.Pressure versus time for material near interface in "Cu" shock tube, for material at positions indicated by white rectangles in Figure 14.Traversing a white rectangle from left to right corresponds to plot colors magenta → red → green → blue.In the Tipton plot, the position of the green curve below all the others shows that the predicted state of the Grüneisen material in the mixed zone is unphysical for many cycles.
Clearly, both models will give the same material volumes at t n+1/2 if they have the same values of ∆V b i,d and ∆V b i,p .The volume changes ∆V b i,d due to zone dilation are given by equations ( 4) and (18), so those quantities will be of the same order of magnitude, and approximately equal if the materials in the zone have roughly equal scaled compliances D i , that is, roughly equal compliances and sound speeds [see Eqn.(13)].It remains to understand the term ∆V b i,p that accounts for pressure-driven exchange of volumes between materials.
Consider a cell containing only two materials, the most common situation requiring closure modeling.We simplify the comparison by observing that the timestep must be bounded by the Courant limit, so the second parenthesized term in (13) dominates the first.If we neglect the first term, the pressure-driven volume change in the Tipton model becomes while for the SSD model we get We find that both expressions give the volume change as a product of a dimensionless number of the order of one, an area, the pressure difference divided by a combination of state variables, and the timestep.Thus, in this approximate treatment of the two-material case, the differences in pressure relaxation between the two models boil down to differences in the first three of those factors (since the fourth, ∆t, is the same for both models).The initial, dimensionless, factor is C b 12 in the SSD model; this is unity unless FCT-like limiting is required, in which case it is between 0 and 1.The corresponding factor in Tipton's model is 2/α, which is never less than 2. This is presumably one reason why we found the initial rate of volume relaxation in the Tipton model calculations to be greater than that of the SSD model, and why it increased with decreasing α (Figure 9).
The area factor is the interface area S 0 12 in the SSD model, consistent with the picture that the interface sweeps out a volume ∆V b i,p as it moves through space.However, the area factor V 0 /L 0 in the Tipton model is simply the cell volume divided by a characteristic length of the cell (L 0 = √ V 0 in two-dimensional FLAG calculations), which may be regarded as an order-of-magnitude approximation to the interface area.Lacking all information about subcell geometry, the Tipton model can do no better than this.
The third factor, depending on the states of the two materials, is remarkably similar between the two models.One difference is that the SSD model depends on the predictorstep state at t n+1/2 , while the Tipton model uses the state at t n ; this is due to the way the latter model was implemented, as described earlier.The other difference is the appearance of volume fractions in the Tipton denominator.If the volume fractions are about equal, the Tipton expression looks very much like the SSD expression.
Note that in the Tipton model, if It is partly due to the lack of this desirable feature that the SSD model requires limiters C ik .

Energy Partition
In both models, the energy partition among materials is specified by a pdV work expression.Here again the Tipton model suffers from a lack of information about the internal structure of the zone, and is forced to use a single effective pressure (17) in its work calculation (16).The SSD model uses a more plausible work calculation (6) in which the volume dilation term for each material uses that material's pressure, and the pairwise volume exchange terms each use the corresponding pairwise Riemann pressure.

Conclusions
The sub-scale dynamics model has been extended to treat mixed cells with any number of materials, and to multiple dimensions.FCT-like limiters have been introduced for robustness and stability.The model has been implemented and tested in FLAG, and found to give results comparable to and in some cases superior to those of Tipton's model.Although its derivation from interface geometry and dynamics is very different from the approach to pressure equilibrium that motivates the Tipton model, we have found close algebraic correspondence between the two models as applied to two-material cells.Their analogous form explains the observed similarity between SSD and Tipton model calculations.Nevertheless, to the extent that interface reconstruction faithfully represents subcell morphologies, and the Riemann solution, subcell dynamics, the SSD model expressions are better suited to predict the resulting state of the cell.The Tipton model relies on overall cell properties (area, pressure) as surrogates for the subcell information it lacks.
In the Tipton model, p is a weighted mean of material pressures (see Eq. 14), and we have expressed the Riemann pressure (7) in the SSD model in a similar way.One difference between models emerges from a comparison of the weight factors used in those means.The SSD weight W i depends only on material states, implying that the interface moves as if it separated two infinite media.The averaging weight w i in the Tipton model depends explicitly on volume fraction, implying that every interface can "feel" the extent of the materials it separates.As we noted above, this should enable the Tipton model to avoid "overshoots" when one volume fraction is very small.On the other hand, when each material has a volume comparable to the cell volume, the interface should not detect the far boundaries of the materials within a Courant-limited timestep.In most cases, we expect this issue to confer an advantage on the SSD model.
Clearly there is much more to be understood about these models, and closure modeling in general.The test calculations presented here indicate strengths and weaknesses to be explored further.By extending the algebraic analysis, we should be able to isolate and understand model dependence on single phenomena and independent variables.In addition, we must study the SSD model in ALE hydro, which is the most common environment requiring closure modeling.

Fig. 1 .
Fig. 1.A single mixed cell: material polygons and exchange volumes of the SSD model.

Fig. 3 .
Fig. 3. Initial mesh (central portion), two-material shock tube set up as three materials.Yellow lines are reconstructed interfaces, not mesh edges.

Fig. 4 .
Fig. 4. SSD model calculation, two-material shock tube set up as three materials.

Fig. 7 .
Fig. 7. Central section of initial mesh for Sod shock tube problems.For the "clean" calculations, the yellow rectangle encloses two unmixed zones.For the closure model calculations, it is a single two-material zone.Figures 8-10 show the state of the materials in that rectangle as functions of time.

Fig. 10 .
Fig. 10.Specific energy at interface in Sod shock tube versus time, in γ-law gases 1 and 2.

Fig. 12 .
Fig. 12. Density at interface in "Cu" shock tube versus time, in γ-law gas and Cu-like fluid.

Fig. 13 .
Fig. 13.Specific energy at interface in "Cu" shock tube versus time, in γ-law gas and Cu-like fluid.

Fig. 14 .
Fig. 14.Central portion of meshes from "Cu" shock tube calculations at t = 0.01, colored by pressure.Material interfaces are in magenta.Black boxes enclose two-material cells at the interface, to which two cells of the clean calculation correspond.Note that the Tipton model has underestimated the pressure of the Grüneisen material in the two-material cells.White rectangles refer to Figure 15.