Numerical methods for interface coupling of compressible and almost incompressible media. ^{†}^{†}thanks: This work was supported in part by NSF grant DMS1216732 (RJL, MdR) and National Council of Science and Technology of Mexico [CONACyT] (MdR).
Abstract
Many experiments in biomedical applications and other disciplines use a shock tube. These experiments often involve placing an experimental sample within a fluidfilled 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 threedimensional (2D axisymmetric) model of these equations is solved using highresolution shockcapturing 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).
.0cm
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 shocktube, 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 shocktube, 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 onedimensional 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 threedimensional (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 transmissionbased 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 onedimensional case [9] and perform a convergence analysis for the twodimensional 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 largescale 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 threedimensional equations into twodimensional axisymmetric Euler equations, which in cylindrical coordinates (,,) take the form
\hb@xt@.01(2.1) 
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 twodimensional 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
\hb@xt@.01(2.2) 
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.
Material  

Air (Ideal gas EOS)  1.4  0.0 
Plastic (polystyrene)  1.1  4.79 
Water  7.15  0.3 
3 Numerical methods
The Euler equations are a nonlinear hyperbolic system of conservation laws, so they can be efficiently solved with highresolution shockcapturing 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 onedimensional Riemann problem for a system of conservation laws like Euler equations can be stated as
\hb@xt@.01(3.1)  
The Euler equatons in this work are solved by implementing a hybrid Riemann HLLCexact type approximate solver for onedimensional Euler equations with interfaces. This solver couples an HLLC approximate Riemann solver to an exact Riemann solver for the Tammann EOS and an EulerianLagrangian 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 wellknown solution to the Euler equations for an ideal EOS [25, 45], we expect our solution will consist of two acoustic waves, the 1wave and 3wave (rarefactions or shocks), and a contact discontinuity, the 2wave 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 1wave, a shock wave or a rarefaction. The state and will be connected by a 2wave, the contact discontinuity with equal pressure and velocity but different density on both sides. The states and are connected by a 3wave, 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 (HartenLaxvan LeerContact) 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 onedimensional 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 RankineHugoniot 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,
\hb@xt@.01(3.2) 
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 onedimensional Euler equations with an ideal gas EOS. However, we want to implement the HLLC solver with the Tammann EOS across an airwater or airplastic 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 multidimensional 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 ,
\hb@xt@.01(3.3) 
For instance, assume we are running a onedimensional 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.
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 onedimensional Riemann problem for the Euler equations with the Tammann EOS. We want to solve the onedimensional Euler equations,
\hb@xt@.01(3.4) 
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 1wave and 3wave (rarefactions or shocks), and a contact discontinuity, the 2wave between them. The system will have four different solution states, separated by the three waves. In order to figure out if the 1wave and 3wave 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 1wave. In a similar manner, we also know the 3wave should be a shock or rarefaction, so we can calculate ; therefore, we define
\hb@xt@.01(3.5)  
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
\hb@xt@.01(3.6) 
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 1wave and 3wave. 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 1wave and 3wave 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]:

1rarefaction, 3rarefaction: and
, 
1shock, 3rarefaction
, 
1rarefaction, 3shock
, 
1shock, 3shock: and
,
where the index indicates if the was obtained by using the RankineHugoniot 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 RankineHugoniot equations and the Riemann invariants. We will denote the speed of the 1wave, , the 2wave, , and the 3wave .
3.2.1 RankineHugoniot conditions for shock waves
As we know the 1wave and the 3wave could each be a shock. In that case, the velocity of the wave, i.e. the shock, will be given by the RankineHugoniot conditions. We will generalize this method for the 1wave velocity and the 3wave velocity , by employing , with .
The RankineHugoniot 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],
\hb@xt@.01(3.7)  
\hb@xt@.01(3.8)  
\hb@xt@.01(3.9) 
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
\hb@xt@.01(3.10)  
\hb@xt@.01(3.11) 
As , from these two equations we can obtain the wave speeds in terms of and ,
\hb@xt@.01(3.12) 
Though, we still need to find and , so we substitute equation (LABEL:ql0) and (LABEL:qr0) into (LABEL:RHmom) to immediately obtain
\hb@xt@.01(3.13)  
\hb@xt@.01(3.14) 
where is defined to simplify future notation with and . Note that and that since and (. Solving for we obtain the equations,
\hb@xt@.01(3.15)  
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,
\hb@xt@.01(3.16) 
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 RankineHugoniot condition (LABEL:RHene).
Finding and : From equation (LABEL:RHene), we can obtain
\hb@xt@.01(3.17) 
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,
\hb@xt@.01(3.18) 
Replacing this result into equation (LABEL:eq:qk), we obtain in terms of and known variables,
\hb@xt@.01(3.19) 
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), 1shock and 3shock 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 2wave are the pressure and the normal velocity . The Riemann invariants for the 1wave and 3wave are the entropy and the quantities,
\hb@xt@.01(3.20)  
\hb@xt@.01(3.21) 
correspondingly. The speed of sound is obtained by applying equation (LABEL:eq:sos) to the Tammann EOS,
\hb@xt@.01(3.22) 
As the entropy is invariant, we can use the Tammann EOS isentropic relation to obtain the density in the middle states,
\hb@xt@.01(3.23) 
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 1rarefaction wave along the rays is then
and for a 3rarefaction 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,
\hb@xt@.01(3.24) 
with,
\hb@xt@.01(3.25) 
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 onedimensional implementations of these methods. In the next Sections, we will study twodimensional 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 fractionalstep 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 RungeKutta or TRBDF2. In a similar manner, the simplest approach to solve the two dimensional system is dimensional splitting. This is done again with a fractionalstep method to split the two dimensional problem up into a sequence of onedimensional problems alternating between solving and . For more details and different splitting algorithms see [25].
Although dimensional splitting is simple to implement, we can obtain secondorder 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.
The twodimensional 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 twodimensional 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.
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].
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 wavepropagation 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 twodimensional 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.
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 leftgoing and rightgoing “fluctuations” and . Each fluctuation , is then further split into downgoing and upgoing 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 upgoing fluctuation from the transverse splitting is based on eigenvectors of and , while the downgoing 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,
\hb@xt@.01(4.1) 
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 upgoing and downgoing fluctuations for are obtained by expanding the fluctuation in terms of these eigenvectors or waves,