<figure><img src="image_1.png" /></figure>

HAL
open science

# An energy-conserving strategy for coupling of finite element method and SPH-ALE method for transient fluid-structure interaction

Zhe Li, Julien Leduc, Alain Combescure, Francis Leboeuf

► To cite this version:

Zhe Li, Julien Leduc, Alain Combescure, Francis Leboeuf. An energy-conserving strategy for coupling of finite element method and SPH-ALE method for transient fluid-structure interaction. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717070

HAL Id: hal-01717070

https://hal.science/hal-01717070

Submitted on 25 Feb 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# An energy-conserving strategy for coupling of finite element method and SPH-ALE method for transient fluid-structure interaction

Zhe LI<fnref n="1" />*, Julien LEDUC<fnref n="2" />, Alain COMBESCURE<fnref n="3" />, Francis LEBOEUF<fnref n="4" />

**Abstract** — This article presents a coupling strategy for transient fluid-structure interaction. The proposed method is applied to a mono-dimensional test case where we study the phenomenon of propagation of shock waves across the fluid-structure interface. The solid sub-domain is discretized by Finite Element Method in Total Lagrangian Formulation with Newmark time integrator, whereas for the fluid sub-domain we use the mesh-less method SPH-ALE with 2^nd^ order Runge-Kutta scheme. At the fluid-structure interface, we impose a continuity condition for velocity to ensure that the interface energy is zero during the whole period of numerical simulation. This coupling method can thus allow us to preserve the minimum order of accuracy in time of the used time integrators for each sub-domain. A good agreement is found between the numerical result and the analytical solution in the 1D shock wave propagation test case. Finally, a multi-dimensional example is presented.

**Keywords** — fluid-structure interaction, SPH-ALE, interface energy, time integrators

## 1 Introduction

An important branch of multi-physics problems is the fluid-structure interaction. Generally the coupling strategies of the two different physical domains can be classified into two major types [6]: monolithic and partitioned procedures. In monolithic procedures, the fluid and solid equations are solved simultaneously, which is quite difficult when different solvers are used for the two sub-domains. The partitioned procedures can overcome this limit. However, there always exists a time lag between the integration of fluid and structure [5]. As a consequence, this partitioned method is typically energy increasing, hence sometimes numerically unstable [7].

As presented in [1], from the energy point of view, we can preserve the minimum level of the order of accuracy of the coupled system, as long as the interface energy is ensured to be zero during the numerical simulation. In [2] Mahjoubi et al. have proposed an energy conserving method to couple heterogeneous time integrators in structural dynamics. This strategy has been successfully used for coupling 3D problems with many sub-domains.

In the present paper, we propose a monolithic approach for coupling finite element method (solid) and a hybrid SPH-ALE method [4] (fluid), using different time integrators: Newmark scheme and 2^nd^ order Runge-Kutta scheme. By imposing the same mean value of normal velocity of the two sub-domains at the fluid-structure interface, we can ensure rigorously the zero interface energy condition, hence we can preserve the order of accuracy in time as well as the numerical stability.

The paper is organized as follows: we first present the governing equations and discretization methods for the fluid and solid sub-domains in Section 2 and Section 3. In Section 4, we describe the proposed coupling strategy which can ensure the zero interface energy condition. Then, we give two 1D test cases compared with the analytical solution in Section 5. A 3D example is also given to show the feasibility of the proposed coupling method for multi-dimensional cases. Finally, the conclusion is offered.

<footnote n="1">¹ Ecole Centrale de Lyon, Laboratoire de Mécanique des Fluides et d'Acoustique, zhe.li@ec-lyon.fr</footnote>

<footnote n="2">² ANDRITZ HYDRO, Research & Development, julien.leduc@andritz.com</footnote>

<footnote n="3">³ INSA de Lyon, Laboratoire de Mécanique des Contacts et des Structures, alain.combescure@insa-lyon.fr</footnote>

<footnote n="4">⁴ Ecole Centrale de Lyon, Laboratoire de Mécanique des Fluides et d'Acoustique, francis.leboeuf@ec-lyon.fr</footnote>

<footnote>* Auteur correspondant</footnote>
# 2 SPH-ALE method for fluid field

## 2.1 Governing equations

We consider a non-viscous, quasi-incompressible fluid on the domain $\Omega_f$, which is governed by the Euler equations in ALE integral form [8]

$$ \left\{ \begin{array}{l} \left. \frac{\partial}{\partial t} \right|_{\chi_0} \int_{\Omega_t} \rho_f \, d\Omega + \int_{\Gamma_t} \rho_f (\mathbf{v}_f - \mathbf{v}_0) \cdot \mathbf{n} \, d\Gamma = 0 \\ \\ \left. \frac{\partial}{\partial t} \right|_{\chi_0} \int_{\Omega_t} \rho_f \mathbf{v}_f \, d\Omega + \int_{\Gamma_t} \rho_f \mathbf{v}_f (\mathbf{v}_f - \mathbf{v}_0) \cdot \mathbf{n} \, d\Gamma = \int_{\Omega_t} (-\nabla p_f + \rho_f \mathbf{g}) \, d\Omega \end{array} \right. \quad (1) $$

where $\rho_f$ denotes the fluid density, $\mathbf{v}_f$ the fluid velocity vector, $p_f$ the fluid pressure. $\mathbf{g}$ is the body force, here it is the gravity vector. $\Omega_t$ represents an arbitrary time-varying volume bounded by a closed surface $\Gamma_t$, which is the referential domain in ALE setting. In this referential domain, $\chi_0$ denotes the coordinate and $\mathbf{v}_0$ is the arbitrary velocity vector. Finally, “$|\chi_0$” means holding the referential coordinate fixed.

In addition, the Tait equation is used as the equation of state for quasi-incompressible fluid [9]

$$ p_f = B \left[ \left( \frac{\rho_f}{\rho_f^{\text{ref}}} \right)^{\gamma} - 1 \right] \quad (2) $$

where $B = \rho_f^{\text{ref}}(c_f^{\text{ref}})^2/\gamma$ and $\gamma = 7$, with $\rho_f^{\text{ref}}$ being the reference density and $c_f^{\text{ref}}$ the reference speed of sound for the fluid.

## 2.2 SPH-ALE method

To discretize the governing equations (1) in space, we apply the hybrid SPH-ALE approach proposed by Vila [4], which combines the Smoothed Particle Hydrodynamics (SPH) and the Arbitrary Lagrangian-Eulerian methods (ALE). Rewriting (1) in a concise form

$$ \left. \frac{\partial}{\partial t} \right|_{\chi_0} \int_{\Omega_t} \Phi \, d\Omega + \int_{\Omega_t} \nabla \cdot (\mathbf{F}_E - \mathbf{v}_0 \otimes \Phi) \, d\Omega = \mathbf{S} \quad (3) $$

where,

$$ \Phi = \begin{cases} \rho_f, & \mathbf{F}_E = \begin{cases} \rho_f \mathbf{v}_f \\ \rho_f \mathbf{v}_f \otimes \mathbf{v}_f + p_f \mathbf{I} \end{cases} \quad \text{and} \quad \mathbf{S} = \begin{cases} 0 \\ \int_{\Omega_t} \mathbf{g} \, d\Omega \end{cases} \end{cases} \quad (4) $$

with $\mathbf{I}$ being the identity tensor, and we will note that $\mathbf{F} = \mathbf{F}_E - \mathbf{v}_0 \otimes \Phi$.

Consider that the whole fluid domain $\Omega_f$ is discretized into $N_f$ fluid “particles”. Each fluid “particle” $\Omega_i$ can be an arbitrary time-varying domain.

The semi-discrete fluid equations write

$$ \frac{d}{dt}(\omega_i \Phi_i) + \omega_i \sum_{k \in \partial D_i} \mathbf{n}_k \cdot W_{ik} (\mathbf{F}_k + \mathbf{F}_i) s_k + \omega_i \sum_{j \in D_i} \nabla_i W_{ij} \cdot (\mathbf{F}_j + \mathbf{F}_i) \omega_j = \mathbf{S}_i \quad (5) $$

where $\omega_i$ denotes the volume of $\Omega_i$, $\Phi_i$ the volume average value of $\Phi$ in each $\Omega_i$, $W_{ij}$ the smoothing kernel function, $\mathbf{n}_k$ the normal vector of support domain truncated by solid boundary (Fig. 1), $s_k$ the area of surface element “k”. $(\mathbf{F}_j + \mathbf{F}_i)$ is estimated by $2G_E(\Phi_i, \Phi_j)$ [4] with

$$ \left\{ \begin{array}{l} \mathbf{G}_E(\Phi_i, \Phi_j) = \mathbf{F}_E(\Phi_{ij}(\lambda_0^{ij})) - \mathbf{v}_0(\mathbf{x}_{ij}, t) \otimes \Phi_{ij}(\lambda_0^{ij}) \\ \Phi_{ij}(\lambda_0^{ij}) = \Phi_E(\lambda_0^{ij}, \Phi_i, \Phi_j) \\ \lambda_0^{ij} = \mathbf{v}_0(\mathbf{x}_{ij}, t) \cdot \mathbf{n}_{ij} \\ \mathbf{x}_{ij} = \frac{\mathbf{x}_i + \mathbf{x}_j}{2} \end{array} \right. \quad (6) $$

with $\Phi_E$ denoting an intermediate status obtained by solving the “moving” Riemann problem between the two fluid states: $\Phi_i$ and $\Phi_j$.
<figure><img src="image_2.png" /><figcaption>Fig. 1: The support domain truncated by solid wall.</figcaption></figure>

In [3], J.C. Marongiu has initially proposed a method for calculating fluid pressure at solid walls, which is coherent with the use of Riemann solvers inside of the fluid domain. In this method, the term $(\mathbf{F}_k + \mathbf{F}_i)$ is calculated by

$$2\mathbf{G}_{E,ik} = 2[\mathbf{F}_E(\mathbf{\Phi}_{E,ik}) - \mathbf{v}_0(\mathbf{x}_k) \otimes \mathbf{\Phi}_{E,ik}] \quad (7)$$

where $\mathbf{\Phi}_{E,ik}$ denotes the state of variables obtained by resolving a “partial Riemann problem” [10]. The fluid pressure at solid wall “$p_k$” writes

$$p_k = \sum_{i \in D_k} \omega_i 2 p_{E,ik} W_{ik} \quad (8)$$

where $p_{E,ik}$ is the intermediate fluid pressure.

# 3 Finite element method for solid field

## 3.1 Governing equations

The balance equation of linear momentum for the solid sub-domain $\Omega_s$ writes

$$\int_{\Omega_s} \left( \rho_s \frac{d^2 \mathbf{u}_s}{dt^2} - \mathbf{\rho}_s \mathbf{b} - \nabla \cdot \boldsymbol{\sigma}_s \right) d\Omega = 0 \quad (9)$$

where $\rho_s$ denotes the solid density, $\mathbf{u}_s$ the displacement vector, $\mathbf{b}$ the body force, and $\boldsymbol{\sigma}_s$ the Cauchy stress tensor. And we note in this paper

$$\left\{ \begin{array}{l} \mathbf{v}_s = \frac{d\mathbf{u}_s}{dt} \\ \mathbf{a}_s = \frac{d^2\mathbf{u}_s}{dt^2} \end{array} \right. \qquad (10)$$

with $\mathbf{v}_s$ and $\mathbf{a}_s$ being the velocity and acceleration vectors.

## 3.2 Discretization method

To discretize the governing equation in space, we apply the Total Lagrangian Finite Element formulation [11]. The semi-discrete linear momentum equations write

$$\mathbf{M}_s \mathbf{a}_s = \mathbf{f}_{ext} - \mathbf{f}_{int} \quad (11)$$

where $\mathbf{M}_s$ is the mass matrix, $\mathbf{f}_{ext}$ and $\mathbf{f}_{int}$ represent the external and internal nodal forces, respectively. For external nodal force, we can write

$$\mathbf{f}_{ext} = \mathbf{L}_p \Lambda \quad (12)$$
<table><caption>Table 1: Two different Newmark schemes used for structure.</caption><thead><tr><th>Schemes</th><th>β</th><th>γ</th></tr></thead><tbody><tr><td>Implicit Newmark</td><td>0.25</td><td>0.50</td></tr><tr><td>Explicit Newmark</td><td>0.00</td><td>0.50</td></tr></tbody></table>

where $L_p$ denotes the geometry operator, and $\Lambda = [p_1 \ p_2 \ \dots \ p_k \ \dots]^T$ with $p_k$ being the pressure applied at the surface element “k”. In linear geometry cases, i.e. the geometry of structure is considered to be constant in time, $L_p$ will not change as time evolves.

For internal nodal force, in linear case, $f_{int}$ can be simply written as

$$f_{int} = K_s u_s \qquad (13)$$

where $K_s$ is the stiffness matrix which is also constant in time.

The Newmark scheme is used as the time integrator for the structure

$$ \left\{ \begin{aligned} \mathbf{u}_s^{n+1} &= \mathbf{u}_s^n + \Delta t \mathbf{v}_s^n + \frac{\Delta t^2}{2} [(1 - 2\beta)\mathbf{a}_s^n + 2\beta \mathbf{a}_s^{n+1}] \\ \mathbf{v}_s^{n+1} &= \mathbf{v}_s^n + \Delta t [(1 - \gamma)\mathbf{a}_s^n + \gamma \mathbf{a}_s^{n+1}] \end{aligned} \right. \qquad (14) $$

with $\beta$ and $\gamma$ being the Newmark scheme's coefficients. Choosing different coefficients, we have two types of Newmark schemes as presented in Table 1.

# 4 Energy conserving coupling strategy

## 4.1 Zero interface energy condition

As presented in [1], when coupling two physical domains with two different time integrators, one can preserve the minimal order of accuracy in time for the coupled system, as long as the interface energy is ensured to be zero during the numerical simulation. For example, when one wants to couple two different time integrators which are both second order accurate in time, the coupled system can have a second order of accuracy, if neither energy injection nor energy dissipation occurs at the interface. In contrast, if the zero interface energy condition can not be ensured, the coupled system will have probably a first order of accuracy or a zero order of accuracy, or even an instable result.

The increment of interface energy over the time interval $t \in [t^n, t^{n+1}]$ is defined as

$$ \Delta W_l = \int_{t^n}^{t^{n+1}} \int_{\Gamma_l} [\mathbf{n}_s \cdot (-p_s \mathbf{I}) \cdot \mathbf{v}_s + \mathbf{n}_f \cdot (-p_f \mathbf{I}) \cdot \mathbf{v}_f] \, d\Gamma \, dt \qquad (15) $$

where $p_s$ and $p_f$ represent the pressure applied to each sub-domain at the interface, and $p_s = p_f = p_k$.

Supposing that $p_k$ and $s_k$ are piece-wise constant in time, one can write (15) as

$$ \Delta W_l^d = \sum_{k=1}^{N_k} \bar{p}_k \bar{s}_k \int_{t^n}^{t^{n+1}} \mathbf{n}_k \cdot [\mathbf{v}_s(\mathbf{x}_k) - \mathbf{v}_f(\mathbf{x}_k)] \, dt \qquad (16) $$

with $\bar{X}$ being the mean value of the variable $X$ in $[t^n, t^{n+1}]$.

Hence, if we impose that the two sub-domains (fluid and solid) have the same mean value of normal velocity at each interface element “k”

$$ \overline{\mathbf{n}_k \cdot \mathbf{v}_s(\mathbf{x}_k)} - \overline{\mathbf{n}_k \cdot \mathbf{v}_f(\mathbf{x}_k)} = 0 \qquad (17) $$

we can ensure rigorously the zero interface energy condition $\Delta W_l = 0$ over the time step $t \in [t^n, t^{n+1}]$, and thus for the whole period of numerical simulations.

## 4.2 Coupling algorithm

The overall coupling procedure is presented in Fig. 2, and the coupling algorithm is presented as follows:
<figure><img src="image_3.png" /><figcaption>Fig. 2: The coupling procedure.</figcaption></figure>

I-1: The fluid solver receives the already known interface status $U_I^n$, and then calculates the fluid status $U_f^{n+1/2}$;

I-2: The coupler uses $U_f^{n+1/2}$ and the solid status $U_s^n$ to solve the system of equations at $t^{n+1/2}$ in order to obtain the interface status $U_I^{n+1/2}$;

II-1: The fluid solver gets $U_I^{n+1/2}$ and finishes the Runge-Kutta scheme to calculate $U_f^{n+1}$;

II-2: The coupler solves the system of equations at $t^{n+1}$ with $U_s^n$ and $U_f^{n+1}$ to calculate the interface status $U_I^{n+1}$;

II-3: The solid solver receives $U_I^{n+1}$ and uses it as the imposed boundary condition to update to $U_s^{n+1}$.

# 5 Numerical results

## 5.1 Mass-spring system coupled with a column of water

Firstly, we couple a mass-spring system with a 1D water tube (Fig. 3). At the other side of this tube, we impose the movement of the solid wall: $x_B(t) = A_m [1 - \cos(\omega t)]$, with $A_m = 2.5 \times 10^{-4}$ m and $\omega = 2000$ rad/s. The length of the tube $L_f = 1$ m, which is discretized into 200 fluid particles. The mass $M_s = 0.8$ kg, the spring stiffness $K_s = 8000$ N/m. The reference parameters for the fluid are chosen to $\rho_f^{ref} = 1000$ kg/m$^3$ and $c_f^{ref} = 1500$ m/s. Initially, the whole system is at rest, and the result is shown in Fig. 4.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Linear test case – mass-spring system coupled with a column of water.</figcaption></figure>

The objective of this linear test case is to verify if the coupling strategy can preserve the order of accuracy in time. Since the used time integrators for the solid and fluid domains are both second order accurate in time, hence if one ensures the zero interface energy condition, one can obtain a second order accuracy in time for the coupled system.

To verify numerically the order of accuracy in time, we choose to apply the method used in [7], which determines the observed order of accuracy by calculating

$$p = \ln \left( \frac{\| X_{4\tau}^{num} - X_{2\tau}^{num} \|}{\| X_{2\tau}^{num} - X_{\tau}^{num} \|} \right) / \ln(2) \quad (18)$$
<figure><img src="image_5.png" /><figcaption>Fig. 4: Coupling result of the linear test case – (a) evolution in time of the mass point velocity; (b) velocity profiles in fluid domain at three different moments ($\Delta t = 10^{-6}$ s).</figcaption></figure>

<table><caption>Table 2: The observed order of accuracy in time of the coupling result.</caption><thead><tr><th>Variables (X)</th><th>$p_{\infty}$</th><th>$p_{L_2}$</th></tr></thead><tbody><tr><td>$u_s$</td><td>1.9313</td><td>2.0059</td></tr><tr><td>$\rho_f(x = 25\%L_f)$</td><td>2.1125</td><td>2.1058</td></tr><tr><td>$v_f(x = 25\%L_f)$</td><td>1.8934</td><td>1.7871</td></tr></tbody></table>

where $p$ denotes the observed order of accuracy in time, $X$ a certain variable calculated by using different time steps ($\tau$, $2\tau$ and $4\tau$) and $|| ||$ means “$L_2$-norm” or “$\infty$-norm”. Table 2 shows that the coupled system possesses a second order of accuracy in time.

## 5.2 1D propagation of shock wave

In the second test case, we replace the mass-spring system by a 1D linear beam (Fig 5). The initial length $L_s^0 = 1$ m, the initial solid density $\rho_s^0 = 2700$ kg/m$^3$, and the initial section area $A_s^0 = 1$ m$^2$. The Young’s modulus $E_s = 67.5$ GPa. The solid beam is discretized into 200 elements, and $\Delta t = 10^{-6}$ s.

<figure><img src="image_6.png" /><figcaption>Fig. 5: 1D linear beam coupled with a column of water: propagation of shock wave across the fluid-structure interface.</figcaption></figure>

The initial configuration of this test case is shown in Fig. 6. The initial discontinuity is located at the fluid-structure interface. The total period of calculation $T = 100\Delta t$. In Fig. 7 and Fig. 8, one can observe a good agreement between the numerical result and the analytical solution.

## 5.3 3D propagation of shock wave – linear structure

Finally, a 3D example is given to show the feasibility for multi-dimensional problems of the presented coupling approach. Fig. 9 shows a water filled tank which is made of five rigid walls and one deformable plate clamped at the four sides. The thickness of the plate $d = 0.01$ m, $a = 0.1$ m, the Young’s modulus $E_s = 100$ GPa, the density $\rho_s = 2700$ kg/m$^3$, and the Poisson’s ratio $\nu = 0.3$. A structural mesh is used to discretize the plate into $20(a) \times 20(a) \times 10(d)$ elements. The length of the tank $L_f = 0.2$ m, and there are $20 \times 20 \times 40$ fluid particles. Initially, the system is at rest, a discontinuity of fluid pressure is located at $x = 15\%L_f$ (Fig. 10).

Normally, the central point of the plate has the maximum amplitude of velocity in $X$. Fig. 11 shows the numerical result of this central point. The total period of simulation $T = 2 \times 10^4\Delta t = 0.004$ s. As time evolves, obvious numerical dissipation can be observed. Notice that the numerical dissipation is
<figure><img src="image_7.png" /><figcaption>Fig. 6: Initial configuration for the test case with 1D linear beam: (a) initial profile of fluid pressure ($p_f$) and solid stress ($-\sigma_s$); (b) initial profile of fluid and solid velocity.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7: Comparison with the analytical solution (Implicit Newmark scheme).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8: Comparison with the analytical solution (Explicit Newmark scheme).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9: Configuration of the 3D test case.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10: Initial fluid pressure field.</figcaption></figure>
introduced only by the fluid solver, since we ensure the zero interface energy condition, and the used
Newmark scheme does not dissipate.

<figure><img src="image_12.png" /><figcaption>Fig. 11: Evolution in time of the central point’s velocity in X.</figcaption></figure>

**6 Conclusion**

An energy conserving coupling approach is presented for transient fluid-structure interaction with dif-
ferent time integrators. Using such coupling strategy, neither numerical energy injection nor dissipation
occurs at the interface during the simulation. Two 1D test cases show that this coupling approach can
preserve the order of accuracy in time for the coupled system, and with this method one can correctly
calculate the interface status when a shock wave impacts onto the fluid-structure interface. A 3D example
is presented briefly to show the feasibility of this approach for multi-dimensional cases.

**References**

[1] A. Combescure, A. Gravouil. A numerical scheme to couple subdomains with different time-steps for predominantly linear transient analysis. *Computer Methods in Applied Mechanics and Engineering* 2002; **191**: 1129-1157.

[2] N. Mahjoubi, A. Gravouil, A. Combescure, N. Greffet. A monolithic energy conserving method to couple heterogeneous time integrators with incompatible time steps in structural dynamics. *Computer Methods in Applied Mechanics and Engineering* 2011; **200**: 1069-1086.

[3] J.C. Marongiu, F. Leboeuf, J. Caro, E. Parkinson. Free surface flows simulations in Pelton turbines using an hybrid SPH-ALE method. *Journal of Hydraulic Research* 2010; **48**: 40-49.

[4] J.P. Vila. On particle weighted methods and smooth particle hydrodynamics. *Mathematical models and methods in applied sciences* 1999.

[5] Blom FJ. A monolithic fluid-structure interaction algorithm applied to the piston problem. *Computer Methods in Applied Mechanics and Engineering* 1998; **167**: 369-391.

[6] Farhat C, Rallu A, Wang K, Belytschko T. Robust and provably second-order explicit-explicit and implicit-explicit staggered time-integrators for highly non-linear compressible fluid-structure interaction problems. *International Journal for Numerical Methods in Engineering* 2010; **84**: 73-107.

[7] Michler C, Hulshoff SJ, van Brummelen EH, de Borst R. A monolithic approach to fluid-structure interaction. *Computers & Fluids* 2004; **33**: 839-848.

[8] Donea J, Huerta A, Ponthot JP, Rodríguez-Ferran A. *Arbitrary Lagrangian-Eulerian Methods*, John Wiley & Sons, 2004.

[9] Macdonald JR. Some simple isothermal equations of state. *Reviews of Modern Physics* 1966; **38**: 669-679.

[10] Dubois F. Partial Riemann problem, boundary conditions, and gas dynamics. *Absorbing Boundaries and Layers, Domain Decomposition Methods: Applications to Large Scale Computations* 2011: 16-77.

[11] Belytschko T, Liu WK, Moran B. *Nonlinear Finite Elements for Continua and Structures*. John Wiley & Sons, 2000.