Numerical methods for interface coupling of compressible and almost incompressible media. This work was supported in part by NSF grant DMS-1216732 (RJL, MdR) and National Council of Science and Technology of Mexico [CONACyT] (MdR).

Numerical methods for interface coupling of compressible and almost incompressible media. thanks: This work was supported in part by NSF grant DMS-1216732 (RJL, MdR) and National Council of Science and Technology of Mexico [CONACyT] (MdR).

M. J. Del Razo Department of Applied Mathematics, University of Washington, Seattle, WA 98195-3925 (,    R. J. LeVeque 22footnotemark: 2

Many experiments in biomedical applications and other disciplines use a shock tube. These experiments often involve placing an experimental sample within a fluid-filled container, which is then placed inside the shock tube. The shock tube produces an initial shock that propagates through gas before hitting the container with the sample. In order to gain insight into the shock dynamics that is hard to obtain by experimental means, computational simulations of the shock wave passing from gas into a thin elastic solid and into a nearly incompressible fluid are developed. It is shown that if the solid interface is very thin, it can be neglected, simplifying the model. The model uses Euler equations for compressible fluids coupled with a Tammann equation of state (EOS) to model both compressible gas and almost incompressible materials. A three-dimensional (2D axisymmetric) model of these equations is solved using high-resolution shock-capturing methods, with newly developed Riemann solvers and limiters. The methods are extended to work on a mapped grid to allow more complicated interface geometry, and they are adapted to work with adaptive mesh refinement (AMR) for higher resolution and faster computations. The Clawpack software is used to implement the method. These methods were initially inspired by shock tube experiments to study the injury mechanisms of traumatic brain injury (TBI).


Key words. Euler equations, Tammann equation of state, compressible and almost incompressible fluid interfaces, finite volume methods, mapped grids, shock tube, traumatic brain injury

AMS subject classifications. 65Nxx, 65Mxx, 65Zxx, 76Nxx, 35Qxx

1 Introduction

A recent collaboration with experimentalists studying traumatic brain injury (TBI) at the Seattle Veterans Administration (VA) Hospital brought to our attention the need for very specific numerical methods [9]. Many experiments performed by the TBI community, as well as in other biomedical disciplines, employ a shock-tube, where they introduce samples to be studied after being exposed to a shock wave. These samples can vary from transwells filled with aqueous solution and cell culture to live mice [6, 9, 15, 16, 17, 20, 21, 37]. Within the shock-tube, the shock wave travels through highly compressible gas before hitting the sample, typically a nearly incompressible material with a fixed location in space. The physical effects of the shock wave hitting the sample are not usually evident from experimental data nor easy to obtain through experimental techniques.

The methods presented in this paper were motivated by this application, although they may be useful in other contexts as well. In order to successfully model the shock wave/sample interaction, we develop numerical methods that can couple the shock wave dynamics in compressible gas with almost incompressible materials, like plastic, water or even bone and brain. Some of these methods have already been employed in our recent collaboration [9]. In the present paper, we give a detailed explanation of the numerical methods and their implementation; we extend them to more complicated interface geometries, enhance their stability in highly refined grids, improve their resolution and efficiency using adaptive mesh refinement (AMR), and further study their convergence. This work will refer to [9] in the sections where it is relevant. Although our simulations can only model idealized scenarios, they can help provide detailed insight into the behavior of the shock wave interaction with interfaces. For instance, in our previous work [8, 9], we obtain the dynamics of a shock wave impacting an interface that models a specific TBI experiment. It also strongly suggested cavitation as a possible damage mechanism, an issue that has been a subject of extense study among the TBI community [9, 15, 27, 30, 31, 32, 33, 35, 50].

Although there is an extensive body of work on computational fluid dynamics with interfaces that is relevant and might be applicable to this type of problems, such as [4, 12, 17, 28, 34, 38, 39, 46, 48] among others, the novel methods presented here are tailored to specifically model a set of experiments performed with a shock tube.

The methods presented here are based on finite volume methods for hyperbolic problems in their wave propagation form [25] and implemented into Clawpack 5.2.2 [5]. The key ingredient in these methods is the Riemann solver, which must be specifically designed to deal with highly nonlinear waves interacting with interfaces between materials having very different properties.

In Section LABEL:sec:TBInumimp, we present the Euler Equations coupled with the Tammann Equation of State (EOS), which can be used to model the various materials involved. In Section LABEL:sec:nummethods, we develop the one-dimensional numerical method in detail, also discussing its implementation in Clawpack [5]. The exact solution in this section is also relevant for the verification studies in Section LABEL:sec:Verif. In Section LABEL:sec:2dnummethods, we extend the numerical methods to apply them to the three-dimensional (2D axisymmetric) model and on mapped grids, allowing more complicated interface geometries. Furthermore, the code was designed to employ Clawpack’s adaptive mesh refinement (AMR) [2] to improve efficiency and resolution. We also discuss the inclusion of transmission-based limiters to reduce numerical oscillations in heavily refined grids produced at the interface corner. In Section LABEL:sec:Verif, we summarize a verification study for the one-dimensional case [9] and perform a convergence analysis for the two-dimensional method. We also show that modifying the original minmod limiter can further reduce the numerical oscillations. The last Section discusses and summarizes some of the results and utility of these methods. Appendix LABEL:sec:TBI:compexp explores the question of whether a thin plastic interface between gas and liquid can be ignored altogether in numerical studies of shock tube experiments [8]. Analysis based on the nonlinear case using the numerical methods from Section LABEL:sec:nummethods suggests it can be ignored. As verification, an analysis based on exact solutions to the linear acoustics equations confirms the result.

The numerical methods and implementation details are explained in this paper; the code is available in GitHub with a BSD license [7].

2 The model

We use the nonlinear compressible Euler equations for compressible inviscid flow, which allow accurate modeling of shock wave formation and propagation. These equations model the conservation of mass, momentum, and energy and provide a direct connection to temperature, which may be important for some biomedical experiments. In this type of experiment, we are not concerned with large-scale movement of the fluid, so viscosity does not play an important role; therefore, employing the inviscid equations is appropriate. In order to model different materials, we use different parameters in the equations of state (EOS) for each material, so we can model the different materials with the same equations. The equations are solved using the methods explained in Section LABEL:sec:nummethods.

An additional advantage of experiments performed in a shock tube is that they often exhibit cylindrical symmetry along the axis that goes through the center of the shock tube. This simplifies the three-dimensional equations into two-dimensional axisymmetric Euler equations, which in cylindrical coordinates (,,) take the form


where is the density; and denote the velocities in the radial and axial direction, and respectively; is the total energy and is the pressure. These equations have the same form as the two-dimensional Euler equations with the addition of geometrical source terms (the right hand side), and are discussed further in Section LABEL:sec:2dnummethods.

2.1 Tammann equations of state

The system of equations (LABEL:eq:Eulercyl) is closed with the addition of an EOS. It is usually given as a relation between pressure, density and specific internal energy, i.e. . The most well known EOS is the one for an ideal gas , where is the ratio of heat capacities. While this EOS is very good for describing the behavior of most gases, it is not appropriate for modeling nearly incompressible materials like water or some elastic solids.

Several alternatives exist; in this work, we will use the stiffened gas EOS, also known as the Tammann EOS. This equation of state is very useful to model a wide range of fluids even in the presence of strong shock waves [11]. The Tammann EOS is given by


where and can be determined experimentally for different materials. The internal energy is related to the total energy by . The Tammann EOS and the ideal gas EOS are the same except for the extra term , where . For fluids with (atmospheric pressure), the relative change in density, when changing the pressure, is very small. Consequently, the Tammann EOS is a good approximation for nearly incompressible fluids and can also be used to model acoustic waves in some elastic solids, like plastic. For sufficiently weak shocks the Tammann EOS can be further simplified to the Tait EOS, see [11], but for greater generality we use the Tammann EOS. Table LABEL:tab:param shows the Tammann EOS parameters for the materials used in the simulations here presented.

Air (Ideal gas EOS) 1.4 0.0
Plastic (polystyrene) 1.1 4.79
Water 7.15 0.3
Table 2.1: Parameters for the Tammann EOS to model the different materials. Note the values for plastic and water are in GPa and are several orders of magnitude above the atmospheric pressure. The parameters for air and water were taken from [11]. As polystyrene is a solid, was chosen to be close to 1, and was adjusted to yield the right speed of sound in polystyrene [29].

3 Numerical methods

The Euler equations are a nonlinear hyperbolic system of conservation laws, so they can be efficiently solved with high-resolution shock-capturing finite volume methods (FVM). This is done by using the wave propagation algorithms described in [25] and implemented in Clawpack [5]. The fundamental problem to solve at each cell interface of our computation is the well known Riemann problem. A general one-dimensional Riemann problem for a system of conservation laws like Euler equations can be stated as


The Euler equatons in this work are solved by implementing a hybrid Riemann HLLC-exact type approximate solver for one-dimensional Euler equations with interfaces. This solver couples an HLLC approximate Riemann solver to an exact Riemann solver for the Tammann EOS and an Eulerian-Lagrangian description coupling at the interface. As the interfaces are represented by contact discontinuities, the HLLC solver is ideal to deal accurately with interface problems. Furthermore, the exact solver will serve as a reference solution to verify the numerical method.

From the well-known solution to the Euler equations for an ideal EOS [25, 45], we expect our solution will consist of two acoustic waves, the 1-wave and 3-wave (rarefactions or shocks), and a contact discontinuity, the 2-wave between them. The -wave refers to the wave corresponding to the -characteristic field (see [25]). This will separate our system in four states, . The left state will be connected to the state by a 1-wave, a shock wave or a rarefaction. The state and will be connected by a 2-wave, the contact discontinuity with equal pressure and velocity but different density on both sides. The states and are connected by a 3-wave, which is a shock wave or a rarefaction.

The method can be extended to two dimensions by using dimensional splitting or transverse solvers. The geometrical source terms can be resolved using a splitting method [25, 26]. In the next paragraphs, we give an overview of the modified HLLC Riemann solver and the exact Riemann solver for the Tammann EOS with discontinuous parameters.

3.1 A modified HLLC solver

The HLLC (Harten-Lax-van Leer-Contact) solver is an approximate Riemann solver for Eq. (LABEL:eq:rphll). The main idea of the HLLC solver is, given the left and right going wave speeds and by some algorithm or approximation, assume a wave configuration of three waves separating four constant states. The Riemann solution to the one-dimensional Euler equations consists of three waves, two acoustic waves with a contact discontinuity in between. The approximate solution for this method will be of the form

where is the approximate wave speed of the contact discontinuity. Assuming we can obtain and , we only need to find , and to solve the problem. These quantities can be obtained by integrating over a box in the plane using the Rankine-Hugoniot conditions and assuming constant pressure and normal velocity across the contact discontinuity, see [45]. The desired states and contact discontinuity speed are given by

where , with are the left or right density and speed in the Euler equations [45].

In order to calculate the wave speeds and , we will need to calculate the sound speed. This is where we require the EOS. A simple estimate is the one given by Davis [45] as

where is the normal velocity and is the sound speed on each side, . Note that the easiest way to calculate the speed of sound is using the EOS . It is usually given in the form,


where is the entropy, and the first derivative is taken along the isentropic curve. A possible improvement is to employ Roe averages in wave speed estimates where and are the Roe averages of the normal velocity and speed of sound respectively [10]. These Roe averages can be calculated using different configurations. Some might be more accurate when dealing with interfaces, as pointed out in [19].

The HLLC solver just discussed works well for the one-dimensional Euler equations with an ideal gas EOS. However, we want to implement the HLLC solver with the Tammann EOS across an air-water or air-plastic interface. The difference between the parameters for different materials in the Tammann EOS are of several orders of magnitude as shown in Table LABEL:tab:param. This generates instabilities in the HLLC solver, more so in the multi-dimensional setting. The instability is generated because we model the interfaces as being fixed in space; however, there is always a displacement of the contact discontinuity, i.e. the interface, even when the material is almost incompressible. The displacement is very small indeed, but it is big enough to render our numerical method unusable. In order to solve this issue, we model each material in Eulerian coordinates using the usual HLLC solver; if any of the cells is next to the interface, we modify our original HLLC or exact solver to work in Lagrangian coordinates, where the interface is actually fixed with respect to the reference frame. This is done by displacing the frame of reference by ,


For instance, assume we are running a one-dimensional simulation of the Euler equations, with a fixed interface modeled by a jump in the parameters of the EOS. The interface is aligned to the edge between cells and , the transformed Riemann solver will be as shown in Figure LABEL:fig:hllc. This will ensure the contact discontinuity velocity is zero and consequently, the interface is modeled as fixed. The wave contributions will be the correct ones since we are just modifying the wave velocity and not the solution ’s. There is, of course, an error made at the interface when coupling the two descriptions; however, as the displacements of the interface are very small due to very low compressibility, this error is small, and it doesn’t cause instabilities as before.

Figure 3.1: Transformation for the HLLC Riemann solver between grid cells and from Eulerian coordinates to Lagrangian coordinates. The transformation can be employed for other Riemann solvers too.

In order to provide better accuracy along the interface, we will also implement an exact Riemann solver for the Tammann EOS. The HLLC solver will be used to model each of the materials in Eulerian coordinates, and the exact solver will be used to solve the Riemann problems at the interface. The transformation to Lagrangian coordinates for the exact solver is equivalent to the one in (LABEL:eultolag).

3.2 Exact Riemann solver for Tammann EOS with a jump in the parameters

The Riemann problem (LABEL:eq:rphll) sometimes can also be solved exactly; the form of the solution will depend on the equations and the EOS being used. An exact solver for the Euler equations coupled with the Tammann EOS for constant parameters was given by Ivings & Toro [22]. In the next paragraphs, we obtain the exact Riemann solver for the Euler equations coupled with the Tammann EOS with different constant parameters on the left and right states. This can be extended numerically to general varying parameters, by averaging them on each cell and using this Riemann solver to provide the solution. The solver is based on the one provided in [22]; however, it extends it to include a jump in the Tammann EOS parameters between the left and right state.

We consider the one-dimensional Riemann problem for the Euler equations with the Tammann EOS. We want to solve the one-dimensional Euler equations,


where is density, velocity, the internal energy and the pressure and the subcripts denote partial derivatives with respect and . The Tammann EOS is given by where the specific internal energy and determines which coefficients to use for the EOS. The initial conditions are given by the left and right constant states and . Note that the state of the system can also be written in terms of the primitive variables by using the equation of state.

As we mentioned before, the solution of the Euler equations will consist of the 1-wave and 3-wave (rarefactions or shocks), and a contact discontinuity, the 2-wave between them. The system will have four different solution states, separated by the three waves. In order to figure out if the 1-wave and 3-wave are rarefactions or shocks, we will need to create a function of the middle state pressure that ensures the velocity across the contact discontinuity is consistent. As we know the velocity on the left state should be connected by a rarefaction or shock to , we can calculate with the jump of the velocity across the 1-wave. In a similar manner, we also know the 3-wave should be a shock or rarefaction, so we can calculate ; therefore, we define


where will change form depending if it’s a shock or a rarefaction (signs were chosen for notation consistency). As we expect these two equations yield the same contact discontinuity velocity , then


This nonlinear equation for will yield the pressure that provides consistency between the type of waves (rarefactions or shocks), their speeds and the contact discontinuity velocity . As we mentioned before, the shape of will depend on whether the states are connected by a shock wave or rarefaction. Once the has been found, the contact discontinuity velocity can be found from (LABEL:CDvel). The only remaining quantity to calculate from the primitive variables is the density. Furthermore, we also need the speeds of the 1-wave and 3-wave. Once we write the explicit equations for our system, it will be clear how to obtain these quantities.

Before writing the equations explicitly, we should first note that having a rarefaction or shock in the 1-wave and 3-wave will depend on the pressure . How can we know which one, can be answered by simple physical intuition. If the pressure is higher on the side toward which the wave is propagating, it will yield a rarefaction. If the pressure is lower, it will be a shock. In the Euler equations, this yields four possible cases for the value of equation (LABEL:eq:phi), just as in the solution using the ideal gas EOS [22, 25]:

  • 1-rarefaction, 3-rarefaction: and

  • 1-shock, 3-rarefaction

  • 1-rarefaction, 3-shock

  • 1-shock, 3-shock: and

where the index indicates if the was obtained by using the Rankine-Hugoniot equations to connect states by shocks or the Riemann invariants to connect them by rarefactions respectively.

In the next paragraphs, we derive the functions for all the four cases with and . We show how to obtain the density and the missing wave speeds. In order to do so, we employ the Rankine-Hugoniot equations and the Riemann invariants. We will denote the speed of the 1-wave, , the 2-wave, , and the 3-wave .

3.2.1 Rankine-Hugoniot conditions for shock waves

As we know the 1-wave and the 3-wave could each be a shock. In that case, the velocity of the wave, i.e. the shock, will be given by the Rankine-Hugoniot conditions. We will generalize this method for the 1-wave velocity and the 3-wave velocity , by employing , with .

The Rankine-Hugoniot conditions are in general given by , where is the vector state variable, the vector state flux and the shock velocity. For the Euler equations this can be easily rewritten as [22],


where , , and the specific enthalpy is given by with the specific internal energy that relates to the internal energy of our original variables by . We will use these relations to find the of equation (LABEL:eq:phi) and the wave speeds .

Finding and and and : We can start by defining the mass fluxes for as


As , from these two equations we can obtain the wave speeds in terms of and ,


Though, we still need to find and , so we substitute equation (LABEL:ql0) and (LABEL:qr0) into (LABEL:RHmom) to immediately obtain


where is defined to simplify future notation with and . Note that and that since and (. Solving for we obtain the equations,


Comparing to equations (LABEL:CDvel), we notice . We also notice we have almost obtained the functions we are looking for, though we still need to find and in terms of known variables.

Finding and : From equations (LABEL:ql0) and (LABEL:ql1), we know that

Solving for , substituting the solution into (LABEL:ql0) and subtituting the for , we obtain a new equation that we can solve for that yields,


with , since we repeated the same process for and obtained exactly the same equation. However, we still don’t know , for this we will need our third Rankine-Hugoniot condition (LABEL:RHene).

Finding and : From equation (LABEL:RHene), we can obtain


where the sign above is used for and the one below for , and we used equations (LABEL:ql0) and (LABEL:qr0) for the second line and (LABEL:eq:qk) for the third line. We can now substitute the specific enthalpy in equation (LABEL:enthdiff) to obtain,

As the interface is the contact discontinuity, the jump in the parameters is only across the contact discontinuity, so . Now we can solve for the unknown density,


Replacing this result into equation (LABEL:eq:qk), we obtain in terms of and known variables,


With equations (LABEL:RHphis) and (LABEL:qkfinal), we can calculate the nonlinear functions of in terms of known variables. The functions allow us to construct equation (LABEL:eq:phi) and solve it using a Newton method or other root finder in order to obtain the value of . Equations (LABEL:RHphis) will then yield the contact discontinuity speed in terms of . Further on, we can calculate and from (LABEL:qkfinal), and we can substitute in (LABEL:RHwspeeds) to obtain the corresponding wave speeds. However, this will only solve the 4th case of equation (LABEL:eq:phi), 1-shock and 3-shock solution. If any of our waves happens to be a rarefaction, we will also need to calculate the functions. This will be obtained using the Riemann invariants.

3.2.2 Riemann invariants for rarefaction waves

Riemann invariants are variables that remain constant through simple waves such as rarefactions. The Riemann invariants across the 2-wave are the pressure and the normal velocity . The Riemann invariants for the 1-wave and 3-wave are the entropy and the quantities,


correspondingly. The speed of sound is obtained by applying equation (LABEL:eq:sos) to the Tammann EOS,


As the entropy is invariant, we can use the Tammann EOS isentropic relation to obtain the density in the middle states,


Solving (LABEL:RI1) and (LABEL:RI2) for and using equations (LABEL:sndsp) and (LABEL:isentrop), we immediately obtain

As , when we compare to equations (LABEL:CDvel) we obtain the functions. The rarefaction head velocities will be given by and ; the tail velocities will be and . For numerical purposes, a simple approximate velocity is provided for and as the average between the head and tail velocity.

In order to compute the complete structure of the rarefaction wave [22, 25], we can use the Riemann invariants from Eqs. LABEL:RI1 and LABEL:RI2, along with Eq. LABEL:sndsp and the isentropic relation from Eq. LABEL:isentrop. The solution for the 1-rarefaction wave along the rays is then

and for a 3-rarefaction wave along the rays is,

Now that we know the functions for the rarefactions, we can construct the function function from (LABEL:eq:phi) for any of the 4 possible scenarios. The value of will be found by numerically finding the roots of . Note which case to employ to calculate might change in each iteration of the root finder. Once is found, , , , and can be found using the relations we just derived depending if it’s a shock or a rarefaction. As we know the three wave speeds , and and the primitive variables on all the 4 states for all the possible cases, we have the solved the Riemann problem.

3.3 Implementation into Clawpack

These methods are implemented into the Clawpack 5.2.2 software [5]. This software employs Godunov’s method [14] with high order corrections and limiters to better handle discontinuities[25]. In order to implement these methods into Clawpack, we need to write Godunov’s method in the wave propagation form. Consider a state vector , a one dimensional conservation law is given by . We partition the space in cells with index and consider the cell average at time to be . Then the Godunov method is given by,




where and are the left and right going fluctuations of the edge of cell respectively, with indicating only those values of with sign , is the number of waves, is the velocity of the characteristic of the Riemann problem at edge , the wave corresponds to the jump across that characteristic and is the limited version of the wave, see [25] for more details.

The numerical solution requires solving a Riemann problem on each cell edge of our partition in order to obtain the fluctuations. The Riemann solutions presented previously have provided the characteristic velocities , and we can calculate the waves by calculating the jump of across the characteristic. This information is calculated for each cell edge and fed into Clawpack, where the method from Eq. LABEL:eq:Godunov is implemented. Appendix LABEL:sec:TBI:compexp contains one-dimensional implementations of these methods. In the next Sections, we will study two-dimensional implementations.

4 Two dimensional axisymmetric model

The three dimensional Euler Equations with cylindrical symmetry can be solved as two dimensional axisymmetric Euler Equations with additional source terms, see Eqs. LABEL:eq:Eulercyl and Figure LABEL:fig:2drevolve. The conservation law for takes the form . In two dimensions, the numerical cell average is calculated as , where is the cell . The source terms can be solved using a fractional-step method [25] by alternating between and . The latter is an ordinary differential equation, which has an exact solution in the case of Equations LABEL:eq:Eulercyl, as shown in [9]. More complex source terms might require implementing another time stepping method like Runge-Kutta or TR-BDF2. In a similar manner, the simplest approach to solve the two dimensional system is dimensional splitting. This is done again with a fractional-step method to split the two dimensional problem up into a sequence of one-dimensional problems alternating between solving and . For more details and different splitting algorithms see [25].

Although dimensional splitting is simple to implement, we can obtain second-order accuracy and less numerical smearing simultaneously by using transverse propagation algorithms from [24]. This will require splitting the normal wave fluctuations at edge into transverse wave fluctuations and . If the normal direction is , then the normal fluctuations are calculated with the flux and the transverse ones with the flux . Our specific model will require a very special kind of transverse solvers, which have been implemented in [9]; a generalized version of these solvers will be explored in detail later in this paper.

Figure 4.1: The axisymmetric model is obtained by revolving the 2D computational grid. The inner square corresponds to the air-water interface. The inside part is filled with water and the outside part is filled with air. All the outer boundaries are modeled with non-reflecting boundary conditions. The interface location was chosen following the source of this figure [9].

The two-dimensional axisymmetric model of Eqs. LABEL:eq:Eulercyl employing a Tammann equation of state with interfaces and transverse solvers were implemented in a traumatic brain injury application in [9]. This work showed how the geometry of the interface can be very relevant and even produce cavitation effects. The set up in [9] and in this work is essentially the one shown in Figure LABEL:fig:2drevolve. A cylindrical plastic container filled with water is placed inside a shock tube. The cylindrical outer boundary corresponds to a cylindrical cross section of the shock tube. The results shown in the Appendix LABEL:sec:TBI:compexp and in [8] show the plastic interface can be neglected in the two-dimensional model.

In this work, the model implemented in [9] is extended to work with AMR capabilities in Clawpack [5, 2]. The AMR implementation requires interpolating the value from coarser grid cells into the finer ones. However, when this interpolation is done across the interface, it will cause instabilities due to the big jump in the EOS parameters across the interface. In order to address this issue, we had to make sure that when a refinement patch intersects the interface, the interpolation for the finest grids is performed only using grid cells corresponding to the same material. For instance, if we need to refine a water grid cell, which is adjacent to the air interface, we will only use the values of adjacent cells corresponding to other water grid cells to obtain the interpolated values in the refined cells. It should be noted that the interface is always aligned to the cell edges, so there are no grid cells that contain two materials. This is also true for the mapped grid case studied below. Figure LABEL:fig:cartesian shows the pressure contours for six different time points for a shock wave traveling in air and hitting a water interface fixed in space, as illustrated in Figure LABEL:fig:2drevolve. The grid is plotted on top showing AMR in action with 4 levels of refinement. The first, second and third coarser grid levels are shown explicitly. The level four refinement is plotted as patches that indicate the highest refinement. Additionally, the code allows us to add gauges to observe the pressure as a function of time at any given point.

Figure 4.2: Axisymmetric simulation pressure contour plots at six different times points , using four levels of AMR. The parameters employed to model water and air for the Tammann EOS are the ones in Table LABEL:tab:param. The pressure amplitude is given along the color bar in KPa. The interface separating air and water is marked as a thick black line, and considering the axis of symmetry is the axis, it models a cylindrical water interface immersed in air. The shock wave travels from the left to right. The first, second and third AMR grid refinement levels are plotted explicitly while the fourth level just shows the refinement patches for clarity. The pressure contours are only shown in the highest refinement level.

In addition to the implementation of these methods in [8, 9], we now show an extension of the algorithms for a mapped grid with adaptive mesh refinement (AMR).

4.1 Two dimensional model in a mapped grid

These algorithms can also be used on a mapped grid where the quadrilateral grid cells are not necessarily rectangular. We will first consider how to implement the normal Riemann solver in the mapped grid. This will require a mapping from a Cartesian grid to a quadrilateral grid, which will tell us the normal at each cell edge where we are solving the Riemann solver as well as the scaling of the edges and the scaling of the areas of the cells. The mapped normal Riemann solver can be done using the same solver as in the Cartesian case by following these steps:

  • Define a mapping;

  • Use the normal at each mapped cell edge to rotate the velocities from the computational domain into normal and transverse components in the physical domain;

  • Solve the Riemann problem as usual with the rotated velocities and calculate the waves;

  • Rotate the waves back into the computational domain;

  • Use the cell edge and area scaling to modify the algorithm in Eq. LABEL:eq:Godunov, see [25].

Figure 4.3: Computational and physical mapped grid of a circular shell inclusion based on the mapping in [3]. The mapping provides two possible circular interfaces, so considering the model is axisymmetric along the axis, it can be used to model a spherical interface or a spherical thick shell interface. The locations of two possible interfaces are shown as thick continuous and dashed lines in both domains.

The mapping of Figure LABEL:fig:grids is based on the mappings of [3]. Consider a computational point on a rectangular grid such that and . The vertical line segment from to will be mapped to a circular arc with radius that intersect the identity diagonals at and . The center of such a circular arc is then given by , and the point in the computational grid is mapped to the physical grid point by

In the mapping of Figure LABEL:fig:grids we have indicated two interfaces: the inner one at radius (1cm) and the outer one at radius (1.5cm) from the origin. The size of the square domain in the computational grid where the mapping is applied is given by a third parameter (4cm), with . The square domain is centered at the origin and the length of each side is . In order to determine the mapping, we need to choose and in the three regions defined by the two interfaces. One option that works well, as shown in Figure LABEL:fig:grids, is given by

Note this is only for the eastern sector of the computational grid, where and ; the other sections are analogous [3].

Some of the quadrilateral cells in the physical domain are nearly triangular, with two adjacent edges nearly colinear. In spite of this, the wave-propagation algorithm with transverse solvers described below works quite robustly in general as discussed further in [3]. However, when there is also a large jump in material parameters at the interface and the grids are adaptively refined there can be some stability issues as discussed further below.

Once the mapping is defined, we proceed by rotating the normal and transverse momentum components and of the Euler equations in the computational grid by using the normal at the current edge of the mapped grid, ,

where and now point in the normal and transverse direction in the physical domain (mapped grid). Using these quantities, we solve the normal Riemann solver as usual to obtain the speeds and waves and , and we rotate the waves back to the computational domain,

Finally, we scale the speeds by the edge scaling to obtain and employ the capacity function (cell area scaling) into a modified version of the algorithm in Eq. LABEL:eq:Godunov found on [25].

The transverse solvers will also be applied on the mapped grid, but this requires more careful consideration because of our treatment of interfaces with huge jumps in the Tammann EOS parameters. This will be explained in detail in the next subsection. In this work, we implemented the 2D axisymmetric model into the mapped grid of Figure LABEL:fig:grids. Although the mapping is two-dimensional and shows half circular inclusions interfaces, the axisymmetry along the axis convert these interfaces into spherical shells. This mapped grid was selected because it could be used to model a skull in computational TBI experiments. Note the mapping allows an inner interface that could even be used to model the thickness of a skull. The code is set up so arbitrary mappings with other interface geometries can be implemented.

In Figure LABEL:fig:maphighP, we show a sample simulation of the pressure contours for the mapped grid at six different points in time. It only employs one interface along the outer circular inclusion shown in the grid of Figure LABEL:fig:grids. Once again the outer part of the circular inclusion is modeled as air and the inner material as water using the same set of parameters as Figure LABEL:fig:cartesian. This figure also shows AMR in action with 4 levels of refinement, and it is also possible to add gauges to observe the pressure as a function of time at any given point in the grid. AMR does not need many additional considerations in terms of the mapped grid since it works on the computational domain, which is still Cartesian. However, it is worth mentioning that the region around the interface is refined to the highest level from the beginning of the simulation. This is to avoid instabilities caused by employing AMR along an interface with huge jumps in the parameters while using a mapped grid with almost triangular grid cells (see Figure LABEL:fig:grids). If any of the conditions is relaxed, i.e. we use a smaller jump in the parameters or use a less severe mapped grid as in Figure LABEL:fig:cartesian, this initial refinement along the interface is no longer required to avoid instabilities.

Figure 4.4: Pressure contour plots of axisymmetric simulation on a mapped grid with a circular inclusion at six different times points , using four levels of AMR. The plot is analogous to that of Figure LABEL:fig:cartesian; however, in this figure the interface separating air and water is circular, which models a spherical water interface. Also note, the region around the interface is refined from the beginning to avoid instabilities when using AMR around corners in the mapped grid.

4.2 Transverse Riemann solver in a mapped grid

A transverse solver for a Cartesian grid was implemented in [9]. In this section, we show the extension of this transverse Riemann solver for a mapped grid. This solver takes the results of a normal Riemann solver and splits it into components moving in the transverse direction. As mentioned in [9], a special transverse solver needs to be developed due to instabilities at the interface. This is based on the solver for acoustics in a heterogeneous media that is described in Section 21.5 of [25].

We recall the basic idea of a transverse solver for a constant coefficient linear hyperbolic system of equations , the jump in normal flux between adjacent cells, , is split via the normal Riemann solver into left-going and right-going “fluctuations” and . Each fluctuation , is then further split into down-going and up-going components and , based on the matrices and .

In the case of variable coefficients or nonlinear problems, the general notation and is used for these two vectors. For variable coefficient acoustics, as described in [25], the up-going fluctuation from the transverse splitting is based on eigenvectors of and , while the down-going fluctuation is based on eigenvectors of and .

At the interface with an almost incompressible liquid, it is difficult to figure out an accurate and stable implementation of the transverse Riemann problem. This is because Euler equations, with a big jump in the parameters at the interface, are extremely sensitive to instabilities. Our first approach was to expand the normal wave as a function of linearized eigenvectors corresponding to the transverse grid cells [25] of the Euler equations. However, this approach resulted in instabilities at the interface. In order to work around this issue, we will follow the same approach as [9] and derive an approximate transverse Riemann solver based on acoustic equations, which capture the acoustic waves while avoiding instabilities.

In this interface, we will mostly be concerned with the two acoustic waves. In order to derive it, let be the transverse unitary normal vector and linearize the acoustic equations around and , with and the velocity in the and direction respectively [25]. In terms of the density and momentum,


where the derivative is taken in the normal direction , is the sound speed and can be understood as a lower dimensional approximation of the transverse Jacobian for the Euler equations. Note we assumed , which is equivalent to move into a Lagrangian frame of reference.

As we might have different materials and sound speeds in the cell above or below, we calculate the eigenvectors and evaluate them according to their location. The matrix of eigenvectors is

where the sound speeds and are the eigenvalues corresponding to the first two column eigenvectors, and . The eigenvalue for the third one is 0. The subindex and refer to cells and when computing and to cells and when computing .

The up-going and down-going fluctuations for are obtained by expanding the fluctuation in terms of these eigenvectors or waves,