An Isotropic Material Remap Scheme for Eulerian Codes

Raymond L. Bell

Computational Shock and Multi-Physics Department

Sandia National Laboratories

Albuquerque, NM 87185

ABSTRACT

Shock Physics codes in use at many Department of Energy (DOE) and Department of Defense (DoD) laboratories can be divided into two classes; Lagrangian Codes (where the computational mesh is ‘attached’ to the materials) and Eulerian Codes (where the computational mesh is ‘fixed’ in space and the materials flow through the mesh). These two classes of codes exhibit different advantages and disadvantages. Lagrangian codes are good at keeping material interfaces well defined, but suffer when the materials undergo extreme distortion which leads to severe reductions in the time steps. Eulerian codes are better able to handle severe material distortion (since the mesh is fixed the time steps are not as severely reduced), but these codes do not keep track of material interfaces very well. In an Eulerian code the developers must design algorithms to track or reconstruct accurate interfaces between materials as the calculation progresses. However; there are classes of calculations where an interface is not desired between some materials, for instance between materials that are intimately mixed (dusty air or multiphase materials). In these cases a material interface reconstruction scheme is needed that will keep this mixture separated from other materials in the calculation, but will maintain the mixture attributes.

This paper will describe the Sandia National Laboratories Eulerian Shock Physics Code known as CTH, and the specialized isotropic material interface reconstruction scheme designed to keep mixed material groups together while keeping different groups separated during the remap step.

KEYWORDS.

Advection, Hydrocode, Shock Physics, Interface Reconstruction, Explicit Time Integration, Half Index Shifted Momenta

1. Code Description

CTH Definition

CTH (a one, two, and three dimensional Eulerian code) has been under development and in use for several decades [1]. It is the direct descendant of the one and two dimensional Eulerian code known as CSQ (developed in the 1970’s) [2]. CSQ evolved from the one dimensional Lagrangian code known as CHART-D, which stands for Combined Hydro And Radiation Transport Diffusion. To get the acronyms out of the way, CSQ is CHART-D squared (going from one dimension to two dimensions) and CTH is CSQ raised to the Three Halves (going from two dimensions to three dimensions). CTH is essentially CHART-D raised to the third (going all the way from one dimensional capability to three dimensional capability).

CTH numerically solves the partial differential equations describing the conservation of mass, momentum, and energy (see Eq’s 1, 2, and 3). It does this in a structured Eulerian mesh fixed in space and uses equations of state (EOS) to close the coupled system of equations.

Where: ρ = mass density

V = velocity vector

σ = stress tensor

Q = artificial viscosity

P = cell pressure

E = specific energy

F = applied body force

S = energy source term

CTH has been licensed to over 340 institutions. Most are DOE Laboratories, DoD Laboratories, or DoD/DOE affiliated contractors.

Computer Platforms.

Initially CTH was targeted for CRAY super computers. The code was highly ‘vectorized’. Since current super computers are massively parallel, recent modifications to the code have emphasized message passing interface (MPI) algorithms. CTH scales quite nicely on “single instruction multiple data” (SIMD) platforms. CTH is routinely run on virtually any platform available from massively parallel supercomputers, to Linux clusters, to laptops running windows. Calculations containing more than a billion cells have been run on SNL massively parallel computers.

Solution Sequence.

Like many Eulerian shock physics codes, CTH uses a two step explicit solution scheme. The first step, referred to as the Lagrangian step, advances the momentum and energy in the mesh through a single time step assuming that the mesh moves with the materials. The second step, referred to as the remap or Eulerian step, remaps the distorted cell values back onto the original Eulerian mesh.

Lagrangian Step. CTH uses a staggered grid in both temporal and spatial meshes. The velocities are defined at, and perpendicular to, the cell faces (the spatial half index locations) at times that are half a timestep behind the cell centered values (pressures, masses, energies, stress deviators, etc.). Gradients of pressure (and divergence of stress deviators) are used to advance the face centered velocities through a timestep during the Lagrangian step, so they end up a half timestep ahead of the cell centered quantities midway through the step. The velocity gradients are then used to advance the cell centered quantites (including the stress deviators) through a timestep. At the end of the Lagrangian step the cell centered quantities are advanced to time(n+1) and the face centered velocities are at time(n+1/2). This scheme is sometimes referred to as a ‘leapfrog scheme’.

Eulerian (Remap) Step. The second step is the remap step (also known as the advection step). This step uses the velocity values to determine the volume overlaps between the distorted mesh and the original Eulerian mesh. These volume overlaps (advection volumes) are then filled with material (masses, momenta, energies, etc.) from the original donor cells. One of the difficulties in this step is determining the proper materials to fill these advection volumes when the donor cell is a mixed material cell (for this step, void is considered a material). If the donor is a single material cell the answer is trivial, since only one material is available. The division between materials when the donor cell contains more than one material will be explained below. In any case, CTH uses a monotone, second order van Leer algorithm [3] to determine the values moved from the donor to the recipient cells. This is accomplished by independently performing advection in each coordinate direction (operator split). The order of these advection operations is permuted to prevent directional bias.

Half Index Shifted (HIS) Momenta. The momenta are advected using the half index shifted (HIS) momentum scheme developed by Benson [6]. This algorithm requires the use of a number of scratch arrays to keep track of the HIS velocities during the remap step. Instead of explicitly using the staggered spatial mesh, each velocity is shifted half a spatial index into each cell sharing the cell face (CTH requires two additional scratch arrays for one dimensional calculations, four for two dimensional calculations, and six for three dimensional calculations). The HIS velocities are used to define HIS momenta (defined at the cell centers). Doing this precludes having to explicitly define the staggered mesh and interpolating advection masses on this staggered mesh. These momenta are advected with the total mass moving from the donor to the recipient cells. After the advection step is completed in a direction, the face centered velocities are recovered by dividing the sum of the residual HIS momenta in each cell sharing a face by the sum of the updated masses in these two cells. The scratch arrays are then free to be used elsewhere.

After the advection step is complete, new material pressures, material temperatures, cell pressures, cell temperatures, and cell sound speeds are defined by calling the material EOS for each material in a computational cell. Multi-material calculations use the multi-material pressure (MMP) options. The cell pressures are usually defined as the volume fraction average of the material pressures in the cell (an exception is the MMP2 option, uses material compressibility as the weighting factor). Cell sound speed is defined as the volume fraction average of the sound speeds determined for each material in the cell. The cell temperature is defined as the heat capacity weighted average of the material temperatures for the materials present in the cell. The next timestep will be determined using the velocities and sound speeds from the cells in the calculation (the cell constituents may be modified below).

Database Modification Step. The final step is referred to as the ‘database modification’ step. Here the user may specify various material discards (material density, energy, temperature values that are somehow ‘bad’, resulting in this material being discarded and replaced by void on a cell by cell basis) or velocity additions (to keep the mesh properly positioned over the volume of interest, usually for penetration or fragment formation calculations). If any materials are ‘discarded’, the EOS step is repeated for the cells involved. The time and cycle numbers are incremented and this procedure (starting again with the Lagrangian step) is repeated for the next timestep.

Periodically (user controlled) restart files are written and plot files are produced, until a user specified problem termination event occurs (the number of cycles desired is reached or the problem stop time is reached).

CTH Code Family.

CTH consists of a family of codes. The CTHGEN program is used to initialize a numerical simulation mesh. CTHGEN is being deprecated because the adaptive mesh refinement (AMR) option is not compatible with this program. When running an AMR calculation, the CTHGEN function is accomplished by CTH. Flatmesh (no AMR) calculations may also be run ‘genless’.

The CTH program actually integrates the simulation through time. The CTHPLT program is used to plot spatial results (color bands or contours of the mesh at given times). HISPLT is the name of the program that is used to plot variables at given tracer points as a function of time.

CTHPLT and HISPLT are being replaced with newer, more robust capabilities known as SPYPLT and SPYHIS. With SPYPLT users may specify ‘plotting on the fly’ rather than waiting to post-process their graphics. By default these utilities produce jpeg images.

CTHED is a utility that can be used to interactively read CTH restart files to determine the values of variables on a cell by cell basis. When running CTH interactively a menu driven option similar to CTHED (known as the ‘hello’ menu) allows users to look in detail at the contents of individual cells. When ‘hello’ is entered in the window where CTH is running the code pauses and allows the user to browse around the mesh.

CTHREZ is the program that may be used to manually rezone the results of a calculation into a new mesh, using the cell variables of the old calculation to initialize a new calculation. When doing the rezone, the mesh size, dimension, and number of materials may be changed. Typically a calculation that starts out one dimensional may be rezoned to a two dimensional cylindrical calculation at a time when it is no longer one dimensional (if a shock encounters a boundary for instance). Subsequently the two dimensional cylindrical calculation may be rezoned into a three dimensional mesh when a shock encounters a non symmetrical boundary (or material interface).

Equations of State and Constitutive Models.

CTH handles up to 20 separate materials in a simulation. It may be run in a pure hydrodynamic mode or, if some of the materials in the calculation support stress, it can take stress deviators into account when calculating the energy and momentum in a simulation.

Equations of State. Many EOS models with predefined material constants are available in CTH. The primary EOS models are ideal gas, Mie Gruneisen, and SESAME tables. Also available for high explosive (HE) materials are the Jones-Wilkins-Lee (JWL) model and some new history variable based models.

A P-Lambda model and a P-Alpha model (embedded in the SESAME and Mie Gruneisen models) are available to model porous material crush-up.

Constitutive Models. Numerous constitutive models with predefined material parameters are available. These models include a simple elastic-plastic model, a geologic material model, several advanced constitutive models for metals, a concrete model, and a couple of ceramic (brittle material) models. CTH also has several specialized fracture models.

Adaptive Mesh Refinement (AMR).

Developers have added an adaptive mesh refinement (AMR) capability that may be used on both single processor and multi-processor platforms. This capability allows users to run much larger calculations involving more diverse length scales. AMR utilizes a ‘patch’ block based scheme. All AMR blocks in a calculation are logically identical (they all contain the same number of cells in the same configuration). The difference is in the cell sizes and the coordinates of the volume being discritized by various blocks. Adjacent blocks may contain cells that are identical, half the size, or twice the size compared to their neighbor block cells. A volume of space in the simulation may be coarsely resolved (if the user is not particularly interested in the results in that area) or it may be highly refined (if the area is of interest). The mesh resolution for a volume of space is determined by user defined ‘refinement indicators’. These indicators are defined by the user to do such things as resolve material interfaces, resolve shocks moving through the mesh, resolve HE materials, etc. The resolution in adjacent blocks is always adjusted to comply with the two to one rule.

Material Interface Reconstruction.

Since CTH is an Eulerian code, the developers have created some material interface reconstruction algorithms. The simplest is the Simple Line Interface Calculation (SLIC) reconstruction algorithm developed at Lawrence Livermore National Laboratory (LLNL). The most complex is the Sandia Modified Youngs’ Reconstruction Algorithm (SMYRA) developed at Sandia National Laboratories (SNL) in Albuquerque, NM [4]. SMYRA was developed from the algorithms designed by David Youngs at the Atomic Weapons Research Establishment in England [5]. All material interface reconstruction methods in CTH use the material (and void) volume fraction information that is stored for each cell. Fig 1 illustrates a three cell patch with two advection volumes and two materials. Note how the advection volumes are divided by the material interface.

2. The Isotropic Remap Option

Motivation.

The impetus for this scheme was the difficulty encountered when remapping multiphase materials. Our initial scheme involved summing the volume fractions, masses, and energies for the parent and child materials into the parent material locations (along with extra variables describing the ‘child to total’ ratios), then zeroing the child values. After the remap step, the child and parent values were recovered using the advected extra variable values. This method required several extra variables and was not easily extended to other mixtures (gas mixtures, etc.).

Interface Reconstruction.

When running a calculation, CTH will attempt to define an interface between all materials during the remap step. If a calculation includes materials that are somehow thoroughly mixed (multiphase materials or a gas mixture for instance) it is often desired that no interface be defined between these materials (see Fig 2 a). Users now have the ability to define groups of materials that should be kept together during remap. The code will then define an interface between groups of materials (all materials will be assigned to a group, some groups will contain only a single material, see Fig 2 b). This scheme will work with any interface reconstruction method used in CTH. If a group contains more than one material, the group volume fraction is determined by simply summing the volume fractions of the group members. The group volume fraction for a single material group is simply the volume fraction for that material.

Group Advection. Advection volumes for the remap step are calculated using the cell face centered velocities. These advection volumes are then sub-divided using the interface reconstruction algorithms. If using the isotropic remap option, the advection volumes are sub-divided by groups rather than individual materials. Once the advection volume for a given group is determined, the advection volume for each material member of the group is simply determined by the relative volume fractions of the group member materials in the donor cell (see Fig 2 c). This procedure has an additional option of using a second order slope on the relative volume fractions in the upstream, donor, and downstream cells to determine the actual division in the group advection volume. In either case the actual volume of each material that is to be advected is uniquely determined.