# On direct Arbitrary-Lagrangian-Eulerian ADER-WENO finite volume schemes for the HPR model of nonlinear hyperelasticity

###### Abstract

This paper is concerned with the numerical solution of the unified first order hyperbolic formulation of continuum mechanics proposed by Peshkov & Romenski [PeshRom2014] (HPR model), which is based on the theory of nonlinear hyperelasticity of Godunov & Romenski [GodunovRomenski72, Godunov:2003a]. Notably, the governing PDE system is symmetric hyperbolic and fully consistent with the first and the second principle of thermodynamics. The nonlinear system of governing equations of the HPR model is large and includes stiff source terms as well as non-conservative products. In this paper we solve this model for the first time on moving unstructured meshes in multiple space dimensions by employing high order accurate one-step ADER-WENO finite volume schemes in the context of cell-centered direct Arbitrary-Lagrangian-Eulerian (ALE) algorithms.

The numerical method is based on a WENO polynomial reconstruction operator on moving unstructured meshes, a fully-discrete one-step ADER scheme that is able to deal with stiff sources [DumbserEnauxToro], a nodal solver with relaxation to determine the mesh motion, and a path-conservative technique of Castro & Parés for the treatment of non-conservative products [Pares2006, Castro2006]. We present numerical results obtained by solving the HPR model with ADER-WENO-ALE schemes in the stiff relaxation limit, showing that fluids (Euler or Navier-Stokes limit), as well as purely elastic or elasto-plastic solids can be simulated in the framework of nonlinear hyperelasticity with the same system of governing PDE. The obtained results are in good agreement when compared to exact or numerical reference solutions available in the literature.

On direct Arbitrary-Lagrangian-Eulerian ADER-WENO finite volume schemes for the HPR model of nonlinear hyperelasticity

Walter Boscheri , Michael Dumbser, Raphaël Loubère

Laboratory of Applied Mathematics, Department of Civil, Environmental and

Mechanical Engineering, University of Trento, I-38123 Trento, Italy.

CNRS and Institut de Mathématiques de Toulouse (IMT)

Université Paul-Sabatier, Toulouse, France

Received …

KEY WORDS: high order direct Arbitrary-Lagrangian-Eulerian finite volume schemes, hyperbolic Peskhov & Romenski model (HPR model), nonlinear hyperelasticity, stiff source terms, non-conservative products, ADER-WENO schemes on unstructured meshes, high order of accuracy in space and time, hyperbolic conservation laws, fluid mechanics and solid mechanics, continuum mechanics

^{1}

^{1}footnotetext: Correspondence to: walter.boscheri@unitn.it (W. Boscheri), raphael.loubere@math.univ-toulouse.fr (R. Loubère), michael.dumbser@unitn.it (M. Dumbser).

## 1 Introduction

The aim of this paper is the numerical solution of the unified first order hyperbolic formulation of continuum mechanics proposed by Peshkov & Romenski [PeshRom2014], denoted as HPR model in the following, which is based on the theory of nonlinear hyperelasticity of Godunov and Romenski [GodunovRomenski72, Godunov:2003a], and which describes fluid mechanics and solid mechanics at the same time in one single system of governing partial differential equations (PDE). In the HPR model the viscous stresses are computed from the so-called distortion tensor , which is one of the primary state variables in this first order system. The appealing property of the HPR model is its ability to describe within the same mathematical framework the behavior of inviscid and viscous compressible Newtonian and non-Newtonian fluids with heat conduction, and, at the same time, the behavior of elastic and elasto-plastic solids. In this model fluids as well as solids are modeled via a stiff source term that accounts for strain relaxation in the evolution equations of the distortion tensor. In addition, heat conduction is included using a first order hyperbolic evolution equation of the thermal impulse which allows the heat flux to be retrieved. The governing system of PDEs is symmetric hyperbolic and fully consistent with the first and the second principle of thermodynamics, as detailed in [PeshRom2014, Dumbser_HPR_16]. However, this system has a large number of equations, is nonlinear and it includes stiff source terms and also non-conservative products.

Consequently, the numerical solution of such a large multi-dimensional system on moving meshes is a big challenge. For this purpose, in this work we propose to employ a high order accurate multi-dimensional ADER-WENO finite volume scheme in the context of direct Arbitrary-Lagrangian-Eulerian (ALE) algorithms. This scheme is constructed with a high order WENO polynomial reconstruction operator on unstructured meshes [DumbserKaeser06b, Dumbser2007204], a one-step space-time ADER integration [titarevtoro, titarevtoro2, Toro:2006a] that is suitably extended for dealing with stiff sources [DumbserEnauxToro, HidalgoDumbser], a nodal solver with relaxation to determine the mesh motion [MaireRezoning, Lagrange2D, Lagrange3D, LagrangeMHD], and a path-conservative integration technique for the treatment of non-conservative products, following the ideas of Castro & Parés [Pares2006, Castro2006], which have been recently extended to the moving-mesh framework in [LagrangeNC, Lagrange3D, ALE-MOOD2]. The proper treatment of boundary conditions is of paramount importance for these simulations on moving meshes. We will pay special attention to them in this work.

In this paper we intend to show that, although the HPR model may seem to be more complex and difficult to solve than other classical ones (Euler & Navier-Stokes equations, linear elasticity or nonlinear hypo-elasticity with plastic strain), the high order ADER-WENO-ALE schemes which allow for a proper treatment of non-conservative terms and stiff source terms [LagrangeNC, Lagrange3D] are an appropriate candidate for this task. Therefore, we will present numerical results obtained with ADER-WENO-ALE schemes for the HPR model in the stiff relaxation limits showing that fluids (Euler or Navier-Stokes limits) as well as pure elastic and elasto-plastic solids can be simulated. In these different situations — fluids, elastic and elasto-plastic solids — which usually require a different mathematical model for each situation, we will numerically prove that the high order accurate ADER-WENO-ALE algorithm is able to reproduce existing exact or numerical reference solutions even for very demanding test cases. These test problems involve shocks (viscous or inviscid ones), contacts and rarefactions in fluids, along with reversible or irreversible deformations in elasto-plastic solids.

The rest of this paper is organized as follows. Section 2 introduces the unified first order hyperbolic Peshkov-Romenski (HPR) model of continuum mechanics, which is numerically solved in this paper. Section 3 presents the high order accurate ADER-WENO-ALE schemes devoted to solve general hyperbolic systems of conservation laws with stiff source terms and non-conservative products. Boundary conditions are discussed in Section 4, while numerical experiments are carried out in Section 5, which also contains a detailed description of these test cases, as well as the obtained numerical results with associated comments. Note that the numerical experiments are designed so that the scheme solves two extreme limits of the HPR model, namely inviscid and viscous fluids (i.e. the compressible Euler equations for gasdynamics and the compressible Navier-Stokes equations) as well as elastic / elasto-plastic solids. Conclusions and perspectives are proposed in the last Section 6.

## 2 The HPR model: a unified first order hyperbolic approach to continuum mechanics

In this work we consider the so-called Hyperbolic Peshkov-Romenski (HPR) model [PeshRom2014], which is the first successful attempt to build a unified formulation of continuum mechanics under a first order symmetric hyperbolic form that includes classical fluid mechanics and solid mechanics just as two special limiting cases of the same formulation. We refer to the recent work of Dumbser et al. [Dumbser_HPR_16], where a detailed introduction to this model is given and where the HPR model has been solved numerically for the first time using high order accurate Eulerian ADER-WENO and ADER-DG schemes on fixed grids, and where many numerical examples have been provided. The HPR model also includes a hyperbolic formulation of heat conduction and it can be written under the form given in [Dumbser_HPR_16] as follows:

(1a) | |||

(1b) | |||

(1c) | |||

(1d) | |||

(1e) |

The solutions of the above PDE system fulfill also the additional conservation of total energy

(2) |

At this point we emphasize that the system above is an overdetermined system of PDE, hence in the numerical solution of the above model we solve the total energy conservation equation (2) and not the entropy equation (1e). Such a choice is mandatory for overdetermined systems. We use the following notation: is the mass density, is the velocity vector, is the distortion tensor, is the thermal impulse vector, is the entropy, is the pressure, is the total energy potential, is the Kronecker delta, is the symmetric viscous shear stress tensor, is the temperature, is the heat flux vector and and are positive scalar functions depending on the strain dissipation time and the thermal impulse relaxation time , respectively. The dissipative terms and on the right hand side of the evolution equations for , and are defined as and , respectively. Accordingly, the viscous stress tensor and the heat flux vector are directly related to the dissipative terms on the right hand side via and . Note that , , and denote the partial derivatives , , and ; they are the energy gradients in the state space or the thermodynamic forces. The Einstein summation convention over repeated indices is implied throughout this paper.

These equations express the mass conservation (1a), the momentum conservation (1b), the time evolution for the distortion tensor (1c), the time evolution for the thermal impulse (1d), the time evolution for the entropy (1e), and the total energy conservation (2). The PDE governing the time evolution of the thermal impulse (1d) looks similar to the momentum equation (1b), where the temperature takes the role of the pressure . Therefore we refer to this equation as the thermal momentum equation.

To close the above system, the total energy potential must be specified. This potential definition will then generate all constitutive fluxes (i.e. non advective fluxes) and source terms by means of its partial derivatives with respect to the state variables. As a consequence the energy potential specification is fundamental for the model formulation.

In order to specify , following [PeshRom2014, Dumbser_HPR_16] we note that there are three scales: the molecular scale, referred to as the microscale; the scale of the material elements, called here mesoscale; and the main flow scale, that is the macroscale. As a consequence it is assumed that the total energy is decomposed into three terms, each of them representing the energy distributed in its corresponding scale, that is:

(3) |

The specific kinetic energy per unit mass refers to the macroscale part of the total energy. The internal energy is related to the kinetic energy of the molecular motion and it is sometimes referred to as the equilibrium energy because it is the only energy which does not disappear in the thermodynamic equilibrium when meso- and macro-scopic dynamics are absent, but only molecular dynamics remains. In this paper, for we will use either the ideal gas equation of state

(4) |

or the Mie-Grüneisen equation of state

(5) |

where has the meaning of the adiabatic sound speed, and are the specific heat capacities at constant volume and at constant pressure, respectively, which are related by the ratio of specific heats . Moreover is the reference mass density and is the reference (atmospheric) pressure. For the mesoscopic, or non-equilibrium, part of the total energy we adopt a simple quadratic form

(6) |

with

(7) |

Here, is the deviator, or the trace-free part, of the tensor and is its trace, is the unit tensor and is the characteristic velocity of propagation of transverse perturbations. In the following we shall refer to it as the shear sound velocity. The characteristic velocity of heat wave propagation is related to the variable .

The fundamental frame invariance principle implies that the total energy can only depend on vectors and tensors by means of their invariants. Hence,

where and , and therefore , as well as the total energy , are only a function of invariants of and .

The algebraic source term on the right-hand side of equation (1c) describes the shear strain dissipation due to material element rearrangements, and the source term on the right-hand side of (1d) describes the relaxation of the thermal impulse due to heat exchange between material elements. Once the total energy potential is specified, all fluxes and source terms have an explicit form. Thus, for the energy given by (6), we have , hence the shear stresses are explicitly given by

(8) |

and the strain dissipation source term becomes

(9) |

where we have chosen , with the determinant of and being the strain relaxation time, also called the particle-settled-life (PSL) time in [Frenkel1955, PeshRom2014]. In other words, this time scale characterizes how long a material element is connected with its neighbor elements before rearrangement occurs. The determinant of must satisfy the constraint

(10) |

where is the density at the reference configuration, see [PeshRom2014]. Furthermore, from the energy potential the heat flux vector follows from as

(11) |

For the thermal impulse relaxation source term, we postulate that yielding

(12) |

The previous formula contains another characteristic relaxation time which is associated to heat conduction. The motivation for this particular choice of and is the connection with classical Navier-Stokes-Fourier theory in the stiff limit and , see [Dumbser_HPR_16] for details.

As shown in [PeshRom2014, Dumbser_HPR_16], the HPR model is compatible with the first and second law of thermodynamics and it constitutes a hyperbolic system of PDEs. For a detailed discussion of the hyperbolicity of nonlinear hyperelasticity, see [Ndanou2014]. For a discussion on the symmetric hyperbolic structure, see [Dumbser_HPR_16] and references therein. In other words the Cauchy problem for the system (1) is well-posed. A detailed discussion of the intrinsic nature of this model can be found in [PeshRom2014, Dumbser_HPR_16] and we refer the interested reader to these references. Further work on nonlinear hyperelasticity can be found e.g. in [GodunovRomenski72, Godunov:1995a, Godunov:2003a, Kluth2010, Pesh2010, Pesh2015, BartonRomenski2010, Barton2013, BartonRom2012]. In this paper we assume the model as given and our goal is to solve it numerically in an accurate, robust and efficient way on moving unstructured meshes using one of the most advanced high order accurate ADER-WENO direct Arbitrary-Lagrangian-Eulerian schemes that it currently available [Lagrange1D, Lagrange2D, LagrangeNC, Lagrange3D].

## 3 High order accurate direct ADER-WENO-ALE schemes for hyperbolic PDE

As already mentioned, the HPR model is a large nonlinear system of hyperbolic balance laws which contains non-conservative products and stiff source terms. To solve this system we consider the arbitrary high order accurate ADER-WENO direct Arbitrary-Lagrangian-Eulerian (ALE) finite volume schemes derived in [Lagrange2D, LagrangeNC, Lagrange3D] that we refer to as ADER-WENO-ALE in the rest of the paper. The HPR model (1) can be cast into the following general formulation which holds in multiple space dimensions :

(13) |

where is the vector of conserved variables, is the conservative nonlinear flux tensor, is the purely non-conservative part of the system written in block-matrix notation and is the vector of algebraic source terms. We furthermore introduce the abbreviation to simplify the notation in some parts of the manuscript.

In our moving-mesh framework the computational domain is discretized at any time level by a set of moving and deforming simplexes . denotes the total number of elements and the union of all elements is referred to as the mesh configuration of the domain: We assume that the computational domain continuously changes in time. Because of this fundamental assumption we adopt the mapping between the physical element to the reference element defined in the reference coordinate system . As usual, the reference element is taken to be the unit triangle in 2D or the unit tetrahedron in 3D, see [Lagrange2D, Lagrange3D].

For any finite volume scheme, data are represented by piecewise constant cell averages both in space and time. As a consequence we define at each time level within the control volume the mean value of the state vector as

(14) |

where is the volume of element . High order of accuracy in space is obtained by means of a polynomial reconstruction technique that provides piecewise high order WENO polynomials from the known cell averages (see next Section 3.1). High order of accuracy in time is further achieved by applying a local space-time discontinuous Galerkin predictor method starting from the high accurate WENO reconstruction polynomials (see Section 3.2). Both techniques are now introduced.

### 3.1 Polynomial reconstruction

#### 3.1.1 Single stencil reconstruction.

The reconstruction operator generates piecewise polynomials of degree which are computed for each element considering the so-called reconstruction stencil and its associated known cell averages. The reconstruction stencil is composed of a number of neighbor elements of , which is bigger than the smallest number

(15) |

needed to reach the nominal order of accuracy in space dimensions, according to [StencilRec1990, Olliver2002, friedrich, kaeserjcp, DumbserKaeser06b]. As suggested in [DumbserKaeser06b, Dumbser2007204], for an unstructured mesh we usually take , with representing the number of space dimensions. The stencil called is defined as , where is a local index counting the elements in the stencil and is a mapping from the local index to the global index of the element in . We rely on the orthogonal Dubiner-type basis functions [Dubiner, orth-basis, CBS-book], defined on the reference element , to explicitly write the high order accurate reconstructed polynomial as

(16) |

where the mapping from to the reference coordinate system is considered and the denote the unknown degrees of freedom, also called expansion coefficients. The procedure to determine the degrees of freedom demands the integral conservation for the reconstruction on each element belonging to stencil , that is

(17) |

The above relations (17) yield an overdetermined linear system of equations for the unknowns that can be solved using either a constrained least squares technique (LSQ), see [DumbserKaeser06b], or a more sophisticated singular value decomposition (SVD) algorithm [DumbserKaeser06b, ADER_MOOD_14].

#### 3.1.2 WENO procedure.

As stated by the Godunov theorem [godunov], linear monotone schemes are at most of order one and if the scheme is required to be high order accurate and non-oscillatory, it must be nonlinear. In this work we consider the pragmatic polynomial WENO approach that has also been adopted in [friedrich, kaeserjcp, DumbserKaeser06b, Dumbser2007204, AboiyarIske, MixedWENO2D, MixedWENO3D, LagrangeNC, Lagrange2D, Lagrange3D, LagrangeNC, LagrangeQF, LagrangeMHD, LagrangeMDRS] to supplement the linear polynomial reconstruction procedure previously described with a nonlinearity. For optimal WENO schemes, see [balsarashu, HuShuTri, shi, ZhangShu3D, WENOsubcell, Semplice2016, Cravero2015]. Seven or nine reconstruction stencils are first determined for and , respectively, and they are further used to compute the associated different polynomials for each cell of the computational domain. These stencils are supposed to cover sufficiently enough “directions” in order to “catch” local oscillatory phenomena. Next, these seven or nine polynomials are blended together using nonlinear weights to obtain the actual high order WENO polynomials . This rather classical procedure has already been described in [DumbserKaeser06b, Dumbser2007204, Lagrange2D, Lagrange3D] and in all the aforementioned references, consequently we omit the details in this paper. However we highly recommend the interested readers to consult these references.

### 3.2 Local space-time Discontinuous Galerkin predictor on moving curved meshes

The reconstructed polynomials computed at time are then evolved during one time step locally within each element , without needing any neighbor information, but still solving the original PDEs (13). As a result one obtains piecewise space-time polynomials of degree , denoted by , that allow the scheme to achieve high order of accuracy even in time. An element-local weak space-time formulation of the governing equations (13) is employed, following the approach developed in the Eulerian framework on fixed grids by Dumbser et al. in [DumbserEnauxToro, USFORCE2, HidalgoDumbser]. According to [DumbserEnauxToro, HidalgoDumbser, DumbserZanotti, Lagrange3D] we adopt the local space-time Discontinuous Galerkin predictor method due to the presence of stiff source terms in the governing equations (1). Let and be the spatial coordinate vectors defined in the physical and in the reference system, respectively, and let and be the corresponding space-time coordinate vectors. Let furthermore be a space-time basis function defined by the Lagrange interpolation polynomials passing through the space-time nodes , which are defined by the tensor product of the spatial nodes of classical conforming high order finite elements in space and the Gauss-Legendre quadrature points in time. Following [DumbserPNPM], the local solution , the fluxes , the source term and the non-conservative products , are approximated within the space-time element with

(18) |

Since the Lagrange interpolation polynomials lead to a nodal basis, we evaluate the degrees of freedom of , and from in a pointwise manner as

(19) |

with representing the gradient of at node . An isoparametric approach is adopted, where the mapping between the physical space-time coordinate vector and the reference space-time coordinate vector is represented by the same basis functions used for the discrete solution . Consequently we have and , where are the degrees of freedom of the spatial physical coordinates of the moving space-time control volume, which are unknown, while denote the known degrees of freedom of the physical time at each space-time node . The mapping in time is simply linear: , then , with denoting the current time. is the time step and it is computed under a classical Courant-Friedrichs-Levy number (CFL) stability condition of the form

(20) |

where is the insphere diameter of element and corresponds to the maximum absolute value of the eigenvalues computed from the solution in . For the HPR model (1) the sound speed is computed according to [PeshRom2014] as

(21) |

On unstructured meshes the CFL stability condition for explicit upwind schemes must satisfy the inequality .

We want the governing PDE formulation (13) to be written in the space-time reference system , hence we first define the Jacobian of the space-time transformation from the physical to the reference element and its inverse:

(22) |

Furthermore let us introduce the nabla operator in the reference space and in the physical space as:

(23) |

and two integral operators

that denote the scalar products of two functions and over the spatial reference element at time and over the space-time reference element , respectively.

The system of balance laws (13) is then reformulated in the reference coordinate system with the following compact notation

(24) |

where we have introduced the unified term by using the inverse of the associated Jacobian matrix (22) and the gradient notation (23). The numerical approximation of is computed by the same isoparametric approach (18), i.e. . Inserting this approximation and (18) into (24), then multiplying (24) with a space-time test function and further integrating the resulting equation over the space-time reference element , one obtains a weak formulation of the original governing system (13):

The term on the left hand side can be integrated by parts in time considering the initial condition of the local Cauchy problem , yielding

(25) |

that simplifies to

(26) |

with the following more compact matrix-vector notation:

(27) |

De facto equation (26) constitutes an element-local nonlinear system of algebraic equations for the unknown space-time expansion coefficients ^{†}^{†}†
This system is solved using the following iterative scheme
where denotes the iteration number. Stiff algebraic source terms are implicitly discretized, see [DumbserEnauxToro, DumbserZanotti, HidalgoDumbser]. .

Together with the solution, we have to evolve the geometry of the space-time control volume which moves in time. The motion of the nodes of element is described by the ODE system

(28) |

with denoting the local mesh velocity. Our direct Arbitrary-Lagrangian-Eulerian (ALE) method allows the mesh velocity to be chosen independently from the fluid velocity. Following the same philosophy as for the solution, the velocity inside element is also expressed in terms of the space-time basis functions as , with the notation . The local space-time DG method is used again to solve (28) for the unknown coordinate vector , according to [Lagrange2D, LagrangeNC], hence

(29) |

where is given by the mapping based on the known vertex coordinates of simplex at time . The above system is iteratively solved together with (26).

Once the above procedure is performed for all cells, an element-local predictor for the numerical solution , for the fluxes , for the non-conservative products , for the source term and also for the mesh velocity is available. This procedure is carried out locally for each cell, consequently it remains to update the mesh motion globally, by assigning a unique velocity vector to each node. To address this issue, in the next section a local nodal solver algorithm for the velocity together with an embedded rezoning technique are presented.

### 3.3 Mesh motion

The aim of any ALE scheme is to follow as closely as possible the material motion. This motion can generate highly deformed cells specifically when fluids or gases are considered. That may drastically reduce the admissible timestep, or, worse, may lead to tangled elements. In order to guarantee good resolution properties for contact waves and material interfaces together with a good geometrical mesh quality, the mesh velocity must be chosen carefully. When natural evidences emanate from the motion of the material boundary conditions, such a mesh velocity can be inferred. However in the general case, specifically for fluids and gases, we adopt a suitable Lagrangian nodal solver technique [Despres2005, Maire2009, chengshu1, chengshu2] to assign a unique velocity vector to each node accurately representing the “true” material velocity. Notice that since we are dealing with a direct ALE formulation the mesh velocity is a degree of freedom. As a consequence we could run our ALE code in a pure Eulerian regime by setting the mesh velocity to zero, or in an almost Lagrangian regime by setting the velocity to an local average of the computed Lagrangian velocities. We could also force any sort of intermediate or artificial mesh motion leading de facto to a so-called ALE motion. In this work the simple nodal solver of Cheng and Shu is used [chengshu1, chengshu2] and the rezoning strategy exposed in [MaireRezoning, Lagrange3D] is employed to locally improve the mesh quality. The final mesh configuration, i.e. the vertex coordinates at the new time level are then computed relying on the relaxation algorithm presented in [MaireRezoning].

### 3.4 Finite volume scheme

The same approach already developed in two and three space dimensions discussed in [Lagrange2D, LagrangeNC, Lagrange3D] is briefly summarized here. To begin with, the governing PDE (13) is more compactly reformulated using a space-time divergence operator :

(30) |

where the space-time flux tensor and the system matrix are given by and . For the computation of the state vector at the new time level , the balance law (30) is integrated over a four-dimensional space-time control volume , which after the application of the theorem of Gauss yields

(31) |

The non-conservative products are treated with the path-conservative approach of Castro and Parés, see [Toumi1992, Pares2006, Castro2006, Castro2008, Rhebergen2008, ADERNC, USFORCE2, OsherNC, LagrangeNC], for a non-exhaustive overview, hence leading to

(32) |

where a new term has been introduced in order to take into account the jumps of the solution on the space-time element boundaries . This term is computed by the path integral

(33) |

where the integration path in (33) is chosen according to [Pares2006, Castro2006, USFORCE2, OsherNC] to be a simple straight-line segment, i.e. , and are the conserved variables in element and its direct neighbor , respectively. Moreover denotes the outward pointing space-time unit normal vector on the varying space-time volume .

Let denote the Neumann neighborhood of simplex , which is the set of directly adjacent neighbors that share a common face with . The space-time volume is composed of space-time sub-volumes , each of them defined for each face of , and two more space-time sub-volumes, and , that represent the simplex configuration at times and , respectively (see [Lagrange3D] for details). Therefore the space-time volume involves overall a total number of space-time sub-volumes, i.e.

(34) |

Each of the space-time sub-volumes is mapped to a reference element in order to simplify the integral computation. For the configurations at the current and at the new time level, and , we use the mapping from the physical to the reference element. The space-time unit normal vectors simply read for and for , since these volumes are orthogonal to the time coordinate. For the lateral sub-volumes we adopt a linear parametrization to map the physical volume to a -dimensional space-time reference prism [Lagrange3D].

Starting from the old vertex coordinates and the new ones , that are known from the mesh motion algorithm described in Section 3.3, the lateral sub-volumes are parametrized using a set of linear basis functions that are defined on a local reference system which is oriented orthogonally w.r.t. the face of , e.g. the reference time coordinate is orthogonal to the reference space coordinates that lie on . The temporal mapping is simply given by , hence and . The lateral space-time volume is defined by six vertices of physical coordinates . The first three vectors are the nodes defining the common face at time , while the same procedure applies at the new time level . Therefore the six vectors are given by

(35) |

The parametrization for reads

(36) |

with , and and the linear basis functions given by

(37) |

The coordinate transformation is associated with a matrix that reads

(38) |

with . Let represent the unit vector aligned with the -th axis of the physical coordinate system and let denote the -th component of vector . The determinant of produces at the same time the quantity of the space-time sub-volume and the space-time normal vector , as

(39) |

where the Levi-Civita symbol has been used according to the usual definition

(40) |

and with

The final one-step direct ALE ADER-WENO finite volume scheme takes the following form:

(41) |

where in the term the Arbitrary-Lagrangian-Eulerian numerical flux function is embedded, as well as the path-conservative jump term, which allows the discontinuity of the predictor solution that occurs at the space-time boundary to be properly resolved also in the presence of non-conservative products. The volume integrals in (41) are approximated using multidimensional Gaussian quadrature rules [stroud] of suitable order of accuracy and the term is evaluated relying on a simple ALE Rusanov-type scheme [Lagrange1D, Lagrange2D, Lagrange3D] as

</ |