A unified hyperbolic formulation for viscous fluids and elastoplastic solids

TL;DR

This paper introduces a unified hyperbolic PDE framework for modeling both viscous fluids and elastoplastic solids, incorporating finite particle length scales and particle rearrangements for more accurate continuum descriptions.

Contribution

It presents a novel hyperbolic formulation that unifies fluid and solid mechanics by explicitly modeling particle deformability and rearrangements, unlike classical models.

Findings

Successfully models fluid and solid flows within a single system

Captures particle rearrangement processes through strain dissipation time

Numerical examples validate the approach's reliability

Abstract

We discuss a unified flow theory which in a single system of hyperbolic partial differential equations (PDEs) can describe the two main branches of continuum mechanics, fluid dynamics, and solid dynamics. The fundamental difference from the classical continuum models, such as the Navier-Stokes for example, is that the finite length scale of the continuum particles is not ignored but kept in the model in order to semi-explicitly describe the essence of any flows, that is the process of continuum particles rearrangements. To allow the continuum particle rearrangements, we admit the deformability of particle which is described by the distortion field. The ability of media to flow is characterized by the strain dissipation time which is a characteristic time necessary for a continuum particle to rearrange with one of its neighboring particles. It is shown that the continuum particle length scale is intimately connected with the dissipation time. The governing equations are represented by a system of first order hyperbolic PDEs with source terms modeling the dissipation due to particle rearrangements. Numerical examples justifying the reliability of the proposed approach are demonstrated.

Equations

$$\dot{X}=-\frac{1}{\lambda }\left(X-X^{\mathrm{NSF}}\right),$$

$$\ell \sim \tau c_s.$$

$$(\mathbf{m},\mathbf{A},\rho ,\sigma ),$$

$$\mathcal{E}=\mathcal{E}(\mathbf{m},\mathbf{A},\rho ,\sigma )$$

$$&\displaystyle\frac{\partial m_{i}}{\partial t}+\frac{\partial\left(m_{i}v_{k}+[m_{l}\mathcal{E}_{m_{l}}+\rho\mathcal{E}_{\rho}-\mathcal{E}]\delta_{ik}+A_{li}\mathcal{E}_{A_{lk}}\right)}{\partial x_{k}}=0,\\ \displaystyle\frac{\partial A_{ik}}{\partial t}+\frac{\partial(A_{il}v_{l})}{\partial x_{k}}+v_{j}\left(\frac{\partial A_{ik}}{\partial x_{j}}-\frac{\partial A_{ij}}{\partial x_{k}}\right)=-\frac{\mathcal{E}_{A_{ik}}}{\theta},\\ \displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial(\rho v_{k})}{\partial x_{k}}=0,\\ \displaystyle\frac{\partial\sigma}{\partial t}+\frac{\partial(\sigma v_{k})}{\partial x_{k}}=\frac{1}{\mathcal{E}_{\sigma}\theta}\mathcal{E}_{A_{ij}}\mathcal{E}_{A_{ij}}\geq 0,$$

$$\frac{\partial \mathcal{E}}{\partial t}+\frac{\partial }{\partial x_k}\left(\mathcal{E}{v}_k+v_n\left([m_l{\mathcal{E}}_{m_l}+\rho {\mathcal{E}}_{\rho }-\mathcal{E}]{\delta }_{nk}+A_{ln}{\mathcal{E}}_{A_{lk}}\right)\right)=0.$$

$${\Sigma }_{ik}=-[m_l{\mathcal{E}}_{m_l}+\rho {\mathcal{E}}_{\rho }-\mathcal{E}]{\delta }_{ik}-A_{li}{\mathcal{E}}_{A_{lk}}$$

$$E=E^1(\rho ,s,\mathbf{A})+\frac{1}{2}{v}_i{v}_i$$

$$E=E^{\mathrm{mi}}(\rho ,s)+E^{\mathrm{me}}(\rho ,s,\mathbf{A})+E^{\mathrm{ma}}(\mathbf{v}).$$

$$E^{\mathrm{me}}=\frac{c_s^2}{4}{G}_{ij}^{\mathrm{TF}}{G}_{ij}^{\mathrm{TF}},$$

$${\sigma }_{ik}=\rho c_s^2{G}_{il}{G}_{lk}^{\mathrm{TF}}.$$

$$\frac{{\mathcal{E}}_{\mathbf{A}}}{\theta }=\frac{3}{\tau }|\mathbf{A}{|}^{\frac{5}{3}}\mathbf{A}{\mathbf{G}}^{\mathrm{TF}},$$

$${\lambda }^{\ast }=\frac{V({\omega }^{\ast })}{{\omega }^{\ast }}\approx \frac{c_{\infty }}{{\omega }^{\ast }}\sim \frac{c_s}{2\pi /\tau }\sim \tau c_s.$$

Full text

1 Introduction

This paper contains an extended abstract of the talk given at the XVI International Conference on Hyperbolic Problems Theory, Numerics, Applications (HYP2016), Aachen (Germany), August 1-5, 2016. The talk was dedicated to the unified hyperbolic formulation of fluid and solid dynamics recently proposed in HPR2016 ; DPRZ2016 . In particular, the emphasis was done on the discussion of the physical model underlying the mathematical formulation. To emphasize how important such a physical interpretation of the mathematical model is, we recall that the equations which constitute the model were proposed many years ago, back to 1970th, by Godunov and Romenski in GodRom1972 ; God1978 for modeling of large elastoplastic deformations in metals, and the equations were used until recently only in the solid dynamics context by several authors, e.g.  Rom1989 ; Resnyansky1995 ; Resnyansky2002 ; GavrFavr2008 ; BartonRom2010 ; Pesh2010 ; favrie2009solid ; Resnyansky2011constitutive ; Barton2013 ; Ndanou2014 ; Pesh2015 ; Boscheri2016 to cite just a few. Moreover, similar equations and even an idea to apply them to modeling of fluids were proposed by Besseling in Besseling1968 , but unfortunately it has never been appreciated in the fluid dynamics context nor by Besseling itself neither by the others. Perhaps, one of the reason for that the hyperbolic Godunov-Romenski equations was not even thought to be used in the fluid dynamics context is an exceptional role of the parabolic Navier-Stokes-Fourier (NSF) equations in the fluid dynamics. For example, it is believed that any mathematical model aiming to describe viscous flows has to literally coincide with the NSF equations in the diffusion regime. This should be understood as that the second order parabolic terms should appear explicitly in the PDEs and they are a fundamental hallmark of the diffusion in the mathematical description. For instance, the well-known first order hyperbolic extension of the NSF equations, the Maxwell-Cattaneo equations

[TABLE]

relax to the NSF equations as the relaxation parameter $\lambda\rightarrow 0$. Here, $X$ is a dissipative quantity in the Maxwell-Cattaneo approach, while $X^{\rm NSF}$ is the value of $X$ obtained in the framework of the NSF theory, the upper dot denotes a time derivative. This, in particular, results in that some characteristic speeds of the Maxwell-Cattaneo equations unphysically tend to infinity as $\lambda\rightarrow 0$. From the other side, as it is shown in HPR2016 ; DPRZ2016 , there are no physical reasons saying that the diffusion processes should be exclusively modeled by the second order parabolic equations, and a radically different first order hyperbolic description which is not based on the steady-state laws such as Newton’s law of viscosity or Fourier law of heat conduction is possible.

Perhaps, the right question in this context is that, after more than one hundred year history of successful use of the NSF theory, do we need at all another transport theory different from the classical parabolic approach? From a practical point of view, the answer is not clear yet, but from a physical viewpoint the answer is obviously positive. Indeed, the heart of the NSF equations, the Newton’s law of viscosity and Fourier’s law of heat conduction, are the phenomenological laws, and thus should be substituted by more physically meaningful laws. We thus would like to emphasize an important role of the physical model in that it helped us to dare to propose an alternative physically-based description of viscous dissipation. Eventually, it is necessary to note that our hyperbolic unified approach is now well established after an extensive comparison with the NSF theory in DPRZ2016 . Moreover, the model was recently extended in HPR2016elmag in order to include the interaction of matter with the electromagnetic field where we also provided an extensive comparison of the extended model with the ideal MHD and parabolic viscous resistive MHD equations.

2 Physical model

Despite we oppose our model to the classical continuum models such as the NSF equations, we underline that the proposed approach entirely relies on the conventional postulates of continuum mechanics and thermodynamics. The main difference though is that we do not assume some simplifications which are implied in the classical theories. Namely, the key difference is that the continuum particles are not treated as scaleless mathematical points but are considered as the finite volumes of a small but finite scale $\ell$. Recall that the notion of the continuum particle is central in any continuum theory. This notion relies on the longtime observations suggesting that for the macroscopic description of the dynamics of matter (gas, liquid or solid) the very detailed information about the molecular motion is irrelevant but the dynamics of ensembles of molecules instead should be considered as a dynamics of new entities of a particle nature. Of course, such particles have not to exist forever but only during a finite time which shall be considered later as an important measure of fluidity. Thus, in our approach, the continuum is represented by a system of finite scale particles (finite volumes) covering the space occupied by the media without gaps, see Fig. 1.

Once one admits or, rather to say, does not ignore that the continuum particles have a finite scale $\ell$, the description of the flow becomes straightforward because the essence of any flow phenomena in any system of finite scale particles is the process of rearrangements of these particles111From the other hand, in the context of the scaleless particles of classical continuum mechanics it is impossible to define the rearrangements because the notion of a neighboring continuum particle becomes indefinite, and thus what remains is not to describe the flow itself but rather to mimic some indirect flow indicators, such as stress–strain-rate relations, etc. Such a mimic strategy is of course admissible in the engineering problems but it is unable to give a meaningful explanation to the physical phenomena.. Thus, the central task of our approach is to find a mathematical framework to express the dynamics of the continuum particles. Further, as depicted in Fig. 1, there is no free volume between the continuum particles, and hence to allow the particle rearrangements, we have to admit the deformability of the particles, otherwise (if they would be rigid volumes) they cannot rearrange and the flow is impossible. Thus, in our approximation, the continuum particles are the structureless (homogeneous) ”soft” deformable particles. The ability of the particles to deform can be characterized by the ability to transmit the transversal perturbations, which in turn can be characterized by the shear sound speed, denoted here by $c_{s}$. As for the measure of the deformation of particles, we use the distortion field $\bm{A}=[A_{ij}]$ which maps a particle from the current deformed state to the undeformed state, see Fig. 2.

Furthermore, the ability of the particles to rearrange, or to change their neighbors, can be characterized by a time $\tau$ which is the characteristic time necessary for a given particle to rearrange with one of its neighboring particles. Because we keep the finite scale of the continuum particles in the physical model, it is then obvious that such particles can not rearrange instantaneously because of the causality principle, and hence the time $\tau$ is also finite. The time $\tau$ is a continuum interpretation of the seminal idea of the so-called particle settled life time of Frenkel Frenkel1955 , who applied it to describe the ability of liquids to flow, see also brazhkin2012two ; bolmatov2013thermodynamic ; bolmatov2015revealing and references therein.

Another nontrivial, and probably the most important, consequence of the finitennes of the particle length scale is that because the particles cannot rearrange instantaneously, there is a relative motion between the neighbouring particles, see Fig. 3. Such a relative motion assumes the existence of a slip between neighbouring particles. In turn, the transversal perturbations that carry the information about the deformation of the continuum particles cannot propagate across such slip planes without a loss of information. This results in that the distortion field is incompatible222The incompatibility condition for $\bm{A}$ is $\bm{B}:={\rm curl(\bm{A})}\neq 0$, where $\bm{B}$ is a so-called Burgers tensor which is interpreted as the number density of the slips (defects) between continuum particles. The term ${\rm curl}(\bm{A})$ also emerges in the time evolution for $\bm{A}$., or not integrable God1978 ; GodRom2003 ; HPR2016 . Such a loss of information is represented by a dissipation term in the time evolution for the distortion field which “dissipates” shear deformation stored in $\bm{A}$. This term is proportional to $1/\tau$, and thus time $\tau$ is also refereed to as the characteristic strain dissipation time in our papers HPR2016 ; DPRZ2016 ; HPR2016elmag .

{svgraybox}

At this point, it is necessary to emphasize that the representation of the continuum by a system of finite volumes is what actually unifies all the three states of mater, gaseous, liquid and solid, because now, the problem of the continuum particle dynamics (finite volumes) is essentially a geometrical problem (deformation problem). Such a geometrical reformulation is insensible to the content of the continuum particles.

It is also clear from the above discussion that the main approximation of our physical model is the treatment of the continuum particles as structureless homogeneous elastic volumes. However, as it is shown in HPR2016 ; DPRZ2016 such an approach is a very precise approximation as long as the characteristic wave length $\lambda$ of the mechanical perturbations is larger than the particle length scale $\ell$. Moreover, as it is proven in DPRZ2016 , the knowledge of only the continuum particle dynamics is sufficient to build a unified flow theory for gases and liquids which incorporates the Newtonian behavior of viscous fluids as a particular case. On the contrary, the exceptionally different molecular dynamics of gases and liquids suggests that not the molecular dynamics is responsible for the mathematical form of the transport laws (identical in both cases) but a dynamics at a larger scale, which we believe is the scale of the continuum particles. Thus, the knowledge of the length scale $\ell$ is extremely important for understanding of the limits of applicability of our physical model, and as will be shown later with the dispersion analysis, the continuum particle length scale is of the order of $\tau\,c_{s}$, i.e.

[TABLE]

If one needs to deal with a problem solution to which strongly depends on the dynamics at a scale for which $\lambda\sim\ell$ or even $\lambda<\ell$ then it is necessary to enlarge the model by providing a more accurate description of the perturbation propagation inside of the continuum particles.

3 Mathematical model

The governing PDEs are formulated for the following volume average quantities

[TABLE]

where $\bm{m}=[m_{i}]=\rho\bm{v}$ is the momentum density, $\rho$ is the mass density, $\bm{v}=[v_{i}]$ is the velocity vector, $\bm{A}=[A_{ij}]$ is the distortion field, and $\sigma=\rho s$ is the entropy density, while $s$ is the specific entropy. Also an exceptional role is played by the total energy density potential

[TABLE]

which plays the role of a generating potential as discussed in details in DPRZ2016 .

The system of governing PDEs can be written as

[TABLE]

while the energy conservation law reads as

[TABLE]

As in all our previous papers HPR2016 ; DPRZ2016 ; HPR2016elmag , the notations such as $\mathcal{E}_{\rho}$, $\mathcal{E}_{m_{i}}$, $\mathcal{E}_{A_{ij}}$, $\mathcal{E}_{\sigma}$ are used to denote the partial derivatives $\partial\mathcal{E}/\partial\rho$, $\partial\mathcal{E}/\partial m_{i}$, etc. Thus, to specify all the terms in the equations, that is to close the system, one needs to specify the total energy $\mathcal{E}$. Also, $\theta\sim\tau$ is a relaxation parameter which will be specified later. These two scalar functions, $\mathcal{E}$ and $\theta$, are the only degrees of freedom in the model formulation. For example, the non-dissipative part of the PDEs, i.e. all the differential terms which are collected on the left-hand side, is complete in the sense that no differential terms can be added or removed and the only possibility to modify something is to change the potential $\mathcal{E}$. The dissipative part of the equations is the only algebraic source terms on the right-hand side which depend on the specification of the energy and the dissipation parameter $\theta$.

The non-advective momentum flux

[TABLE]

is the stress tensor. Its form is completely defined by the structure of the time evolution equations while its further specification depends solely on the choice of the energy $\mathcal{E}$. Here, the scalar $p=m_{l}\mathcal{E}_{m_{l}}+\rho\mathcal{E}_{\rho}-\mathcal{E}$ can be refereed to as the pressure which coincides with the classical hydrodynamic pressure for equilibrium flows. Indeed, if one introduces the specific total energy density $E$ as $\mathcal{E}=\rho E$, for which the following decomposition is usually assumed

[TABLE]

then $p=m_{l}\mathcal{E}_{m_{l}}+\rho\mathcal{E}_{\rho}-\mathcal{E}=\rho^{2}E^{1}_{\rho}$, exactly as in our previous paper HPR2016 . The last term in (7) represents the viscous stresses or elastic stresses in the case of solid dynamics.

For the further specification of the total energy potential $\mathcal{E}$, we note that there are three scales involved in the physical model formulation described in Introduction. Namely, the molecular scale, called here microscale; the scale of the continuum particles, called here mesoscale; and the flow scale, or observable macroscale. We thus assume that $E$ is a sum of three terms each of which represents the amount of energy stored on the corresponding scale

[TABLE]

The terms $E^{\rm mi}$ and $E^{\rm ma}$ are conventional. They are the kinetic energy $E^{\rm ma}(\bm{v})=\frac{1}{2}v_{i}v_{i}$, which represents the part of the total energy stored in the macroscale, and an internal energy $E^{\rm mi}(\rho,s)$ represents the kinetic energy of the molecular motion. In DPRZ2016 ; HPR2016elmag , we used an ideal gas equation or stiffened gas equation of state for $E^{\rm mi}$ to model gases or liquids and solids, respectively

For the mesoscopic, or non-equilibrium, part of the total energy, we shall use a quadratic form

[TABLE]

where $G_{ij}^{\rm TF}=G_{ij}-G_{ii}/3$ is the deviator of the tensor $G_{ij}=A_{ki}A_{kj}$, $c_{s}$ is the characteristic velocity of propagation of transversal perturbations, we call it here shear sound velocity. In general, $c_{s}$ is a function of $\rho$ and $s$.

With such a specification of the term $E^{\rm me}$, the explicit form of the viscous/elastic stress (the last term in (7)) which we denote by $\sigma_{ik}$ is

[TABLE]

The mesoscopic energy $E^{\rm me}$ also defines DPRZ2016 the dissipation terms as

[TABLE]

where we use $\theta=\tau c_{s}^{2}/3|\bm{A}|^{\frac{5}{3}}$ for $\theta$ and $|\bm{A}|$ to denote the determinant of $\bm{A}$. In general, $\tau$ is a function of the state variables $\tau=\tau(\rho,s,\bm{A})$ while for Newtonian fluids it can be taken to be constant as shown in DPRZ2016 through a formal asymptotic analysis. In particular, the dependence of $\tau$ on $\bm{A}$ defines the non-Newtonian properties of fluids or controls the transition from elastic to plastic regime in solids God1978 ; GodRom2003 ; BartonRom2010 ; Pesh2010 , see also numerical examples in the following Section 4.

4 Numerical results

In this section, we demonstrate that the proposed model can be applied to modeling of nonequilibirum effects in gases as well as to modeling of viscous fluid flows and elastoplastic deformation in metals.

4.1 Non-equilibrium sound wave propagation in a viscous gas

We first study the propagation of plane acoustic waves of an angular frequency $\omega$ in a viscous gas. As it is well known the presence of the dissipative process gives rise to the phenomena called dispersion when the wave phase speed $V$ depends on the frequency of the wave, $V=V(\omega)$. This dependency is defined by the dispersion relation for a given model. The dispersion relation for the proposed hyperbolic model can be found in HPR2016 in Section 2.2.2. The phase velocity $V(\omega)$ and the attenuation factor for equations (5), (9) and (10) are presented in Fig. 4 for Helium.

As can be seen in Fig.4 (left), at a frequency $\omega^{*}=2\pi/\tau$ the dispersion almost disappears and the phase velocity $V(\omega)$ tends to a constant value $c_{\infty}=\sqrt{c_{0}^{2}+4c_{s}^{2}/3}$ (see also HPR2016 ) called a high frequency limit for the sound speed. This experimentally defined value can be used to estimate the shear sound speed $c_{s}$, and subsequently to estimate the dissipation time $\tau$ from the relation $\eta=\frac{\rho}{6}c_{s}^{2}\tau$ for the shear viscosity $\eta$.

The dispersion disappearance of $V(\omega)$ is fully conditioned by the physical model underlying the mathematical formulation. Indeed, because the continuum particles have the finite scale $\ell$, the behavior of $V(\omega)$ should change when the wave length $\lambda$ becomes comparable with the particle size, $\lambda=\lambda^{*}\sim\ell$. We thus can use this fact to estimate the particle length scale $\ell$ as $\ell\sim\tau c_{s}$. Indeed,

[TABLE]

Thus, for the experimental data presented in Fig.4, the continuum particle length scale can be estimated as $\ell\approx 4.6\cdot 10^{-7}$ m. For this, we took $c_{\infty}=2\,052$ m/s, shear viscosity $\eta=2\cdot 10^{-5}$ Pa$\cdot$s and mass density $\rho=0.16$ kg/m3 which gives us $c_{s}=\sqrt{(c_{\infty}^{2}-c_{0}^{2})3/4}\approx 1\,543$ m/s and $\tau=6\eta/(\rho c_{s}^{2})\approx 3\cdot 10^{-10}$ s.

Eventually, we note that there is a certain discrepancy in the attenuation factor visible in the Fig.4(right) if compared with the experimental data, while the Navier-Stokes-Fourier model (see Chapter 11 of EIT2010 for the dispersion relation for the NSF equations) shows an excellent agreement. First of all, one should note that, at high values of $\omega$, there may be a contribution to the absorption arising from diffusion in the piezoelectric receiver, as pointed out by Woods and Troughton Woods1980 (see also discussion in Chapter 11 of book EIT2010 ), so that the experimental result for the absorption factor should be considered as an upper limit to the actual value. Secondly, so far, we ignore such important processes as heat transfer and volume relaxation which of course should increase the dissipation. This will be studied elsewhere.

4.2 Viscous fluids and elastoplastic solids

In order to demonstrate the ability of the proposed unified hyperbolic model to deal with radically different behaviors of matter such as flows of viscous gases and elastoplastic deformation in solids, we consider two 2D problems. These examples merely serve to demonstrate the diversity of regimes allowed to be captured by the model while the numerical schemes we use in this paper are not very accurate such as those used in our recent papers DPRZ2016 ; HPR2016elmag ; Boscheri2016 where much more accurate results were obtained with the use of advanced high order ADER (arbitrary high order derivatives Toro:2006a ) Discontinuous Galerkin and Finite-Volume schemes, moving mesh and adaptive mesh refinement techniques. An extensive comparison against the parabolic theories like the Navier-Stokes-Fourier equations and resistive MHD model is also provided in DPRZ2016 ; HPR2016elmag .

The key parameter controlling the transition between the fluid-like and solid-like behavior is the dissipation time $\tau$. As discussed in the introduction section and in HPR2016 ; DPRZ2016 , for the elastic solids, the continuum particles do not rearrange and hence time $\tau$ is infinite, while for viscous fluids $0<\tau<\infty$. For elastoplastic solids, time $\tau$ depends on the yield strength and rapidly but continuously changes from an infinite value (in fact from a sufficiently large value) to a finite value in the plastic regime, in which the continuum particles do rearrange.

In the first example, a gravity driven Rayleigh-Taylor instability in a viscous gas confined in a rectangular domain with no-slip boundary conditions is simulated. The domain is a box $(x,y)\in[0,1/3]\times[0,1]$ which were discretized with a Cartesian mesh consisting of $200\times 600$ cells, gravitational field is directed vertically downward and has a magnitude $g=0.1$, the initial conditions are: $\bm{v}=0$, the density is taken 2 if $y>0.5+0.01\cos(6\pi x)$ and 1 otherwise, the ideal gas equation of state is used for the internal energy $E^{\rm mi}$ (see (9)) with the ratio of specific heats $\gamma=1.4$ the pressure is set to $1/\gamma$ everywhere, the shear sound speed $c_{s}$ was set to $c_{s}=c_{0}=1$. For the whole domain, we set the shear viscosity to $10^{-5}$ Pa$\cdot$s and the dissipation time to $\tau=6\eta/c_{s}^{2}\approx 6\cdot 10^{-5}$ s. Fig. 5 depicts several time instants of the simulation. The fluid-like motion (formation of the vortexes) is clearly identified. We also note that the gas sticks to the walls as no-slip boundary conditions is used. For this simulation, we use CLAWPACK software Clawpack designed specifically for hyperbolic PDEs. The numerical fluxes are obtained via the solution of an approximate Riemann problem which was solved using the eigenvalue decomposition of the Jacobi matrix for the fluxes. The second order wave propagation algorithm of CLAWPACK with the “minmod” limiter was used.

In the second numerical example, we consider an oblique high velocity collision of two solid plates presented in Fig. 6. This example is motivated by the explosion welding process Godunov1970welding . The initial angle between the plates is 13 degrees, the lower plate is at rest while the upper plate has the velocity 1500 m/s normal to the bottom face. It is assumed that there is no slip on the contact interface between the bodies. The upper plate has dimensions $1\times 0.1$ cm and is discretized with a $600\times 60$ Cartesian mesh while the lower plate has dimensions $1\times 0.3$ and is discretized with a $600\times 180$ Cartesian mesh. The both plates have the same material parameters which were taken as follows. The mass density is 7.9 kg/m3, the longitudinal sound speed is $6\,700$ m/s, the shear sound speed is $3\,150$ m/s, the dissipation time was taken as $\tau=\tau_{0}(\sigma_{0}/\sigma)^{n}$ where $\tau_{0}=0.1$ s, $\sigma_{Y}$ is a parameter which controls the transition from elastic to plastic regime and was set to $0.0056$ GPa, while $\sigma$ is the second norm of the deviator of the stress tensor. The power law index $n$ was set to $n=10$. For such parameters, the effective yield strength appears to be $\approx 0.04$ GPa. The sound speeds were taken as $6\,700$ m/s and $3\,150$ m/s for longitudinal and transversal sound speed respectively. The equation of state is given in Pesh2010 . This numerical example was performed several years ago, and model (5) was written in the Lagrangian coordinates which are well suited for the solid dynamics problems, see details in Pesh2010 . However, recently, the Eulerian equations (5) were also implemented in an ALE (arbitrary Langrangian Eulerian) code Boscheri2016 which also opens new possibilities for more efficient simulation of elastoplastic solids experiencing large deformations. Two independent Cartesian meshes were used in this simulation and the equations were solved with a standard first-order Godunov scheme with an acoustic Riemann solver, see Pesh2010 . The contact boundary requires a specific treatment. For this purpose, the contact cells have to be detected and the numerical flux on the contact interface is obtained with the same Riemann solver that used for the internal cells.

Acknowledgements.

I.P. have received funding from ANR-11-LABX-0040-CIMI within the program ANR-11-IDEX-0002-02 and a partial support from the Russian Foundation for Basic Research (grant number 16-31-00146). E.R. acknowledges a partial support by the Program N15 of the Presidium of RAS, project 121 and the Russian Foundation for Basic Research (grant number 16-29-15131). M.D. have received funding from the European Union’s Horizon 2020 Research and Innovation Programme under the project ExaHyPE, grant agreement number no. 671698 (call FETHPC-1-2014).