Plane-strain waves in nonlinear elastic solids with softening

TL;DR

This paper investigates the propagation of nonlinear elastic waves in damaged, softening media like concrete and rocks, using a combined theoretical and numerical approach to understand wave behavior and solitary wave generation.

Contribution

It introduces a novel constitutive model that integrates Murnaghan hyperelasticity with slow dynamics to analyze nonlinear wave propagation in damaged solids.

Findings

Softening affects wave propagation characteristics.

Numerical simulations show formation of solitary waves.

Material softening influences wave speed and shape.

Abstract

Propagation of elastic waves in damaged media (concrete, rocks) is studied theoretically and numerically. Such materials exhibit a nonlinear behavior, with long-time softening and recovery processes (slow dynamics). A constitutive model combining Murnaghan hyperelasticity with the slow dynamics is considered, where the softening is represented by the evolution of a scalar variable. The equations of motion in the Lagrangian framework are detailed. These equations are rewritten as a nonlinear hyperbolic system of balance laws, which is solved numerically using a finite-volume method with flux limiters. Numerical examples illustrate specific features of nonlinear elastic waves, as well as the effect of the material's softening. In particular, the generation of solitary waves in a periodic layered medium is illustrated numerically.

Tables

Table 1. Table 1: Physical parameters of concrete.

${\rho }_0$ (kg m$-3$)	$\lambda$ (GPa)	$\mu$ (GPa)	$𝔩$ (GPa)	$𝔪$ (GPa)	$\gamma$ (J m $-3$)	${\tau }_1$ (J m$-3$ s)
$2400$	$12.7$	$17.5$	$-3007$	$-2283$	$4.0\times {10}^{-2}$	$2.0\times {10}^{-6}$

Equations

$$\mathbf{F}=\mathbf{grad}{\mathbf{x}}_t=\mathbf{I}+\mathbf{grad}\mathbf{u},$$

$$\dot{\mathbf{F}}=\mathbf{grad}\mathbf{v},$$

$${\rho }_0\dot{\mathbf{v}}=\mathbf{div}\mathbf{P}+{\mathbf{f}}^v,\text{where}\mathbf{P}=\det (\mathbf{F})\mathbf{\sigma }\cdot {\mathbf{F}}^{-\top }$$

$$\mathbf{E}=\frac{1}{2}\left(\mathbf{grad}\mathbf{u}+{\mathbf{grad}}^{\top }\mathbf{u}+{\mathbf{grad}}^{\top }\mathbf{u}\cdot \mathbf{grad}\mathbf{u}\right).$$

$${\rho }_{0}e=\left(1-g\right)W(\mathbf{E})+\Phi (g),$$

$$W(\mathbf{E})=\frac{\lambda +2\mu }{2}{E_{\text{I}}}^2-2\mu E_{\text{II}}+\frac{\mathfrak{l}+2\mathfrak{m}}{3}{E_{\text{I}}}^3-2\mathfrak{m}{E}_{\text{I}}{E}_{\text{II}}+\mathfrak{n}{E}_{\text{III}},$$

$$\Phi (g)=-\frac{1}{2}\gamma \ln (1-g^2),$$

$$\mathbf{P}=\left(1-g\right)\mathbf{F}\cdot \frac{\partial W}{\partial \mathbf{E}},\text{and}{\tau }_1\dot{g}=W(\mathbf{E})-{\Phi }^{\prime }(g),$$

$$\frac{\partial W}{\partial \mathbf{E}}={\alpha }_0\mathbf{I}+{\alpha }_1\mathbf{E}+{\alpha }_2{\mathbf{E}}^2,$$

$${\alpha }_0$$

$${\alpha }_1$$

$${\alpha }_2$$

$$\left(\begin{array}{c}{u}_{i,j}\end{array}\right)=\left(\begin{array}{ccc}{u}_{1,1}& u_{1,2}& 0\\ u_{2,1}& u_{2,2}& 0\\ 0& 0& 0\end{array}\right).$$

$$\frac{\partial W}{\partial E_{ij}}={\tilde{\alpha }}_0{\delta }_{ij}+{\tilde{\alpha }}_1{E}_{ij}$$

$${\tilde{\alpha }}_0$$

$${\tilde{\alpha }}_1$$

$$P_{ij}=(1-g)\left({\delta }_{im}+u_{i,m}\right)\left({\tilde{\alpha }}_0{\delta }_{mj}+{\tilde{\alpha }}_1{E}_{mj}\right)$$

$${\partial }_t\mathbf{q}+{\partial }_x\mathbf{f}(\mathbf{q})+{\partial }_y\mathbf{g}(\mathbf{q})=\mathbf{r}(\mathbf{q})+\mathbf{s},$$

$$\mathbf{f}(\mathbf{q})=-\left(\begin{array}{c}{v}_1\\ 0\\ v_2\\ 0\\ P_{11}/{\rho }_0\\ P_{21}/{\rho }_0\\ 0\end{array}\right),\mathbf{g}(\mathbf{q})=-\left(\begin{array}{c}0\\ v_1\\ 0\\ v_2\\ P_{12}/{\rho }_0\\ P_{22}/{\rho }_0\\ 0\end{array}\right),\mathbf{r}(\mathbf{q})=\frac{1}{{\tau }_1}\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\\ W-{\Phi }^{\prime }(g)\end{array}\right),\mathbf{s}=\frac{1}{{\rho }_0}\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ {\mathbf{f}}^v\cdot {\mathbf{e}}_1\\ {\mathbf{f}}^v\cdot {\mathbf{e}}_2\\ 0\end{array}\right).$$

$${\mathbf{f}}^{\prime }(\mathbf{q})=-\left(\begin{array}{ccccccc}& & & & 1& 0& 0\\ & & & & 0& 0& 0\\ & & & & 0& 1& 0\\ & & & & 0& 0& 0\\ Q_{1111}& Q_{1112}& Q_{1121}& Q_{1122}& 0& 0& G_{11}\\ Q_{2111}& Q_{2112}& Q_{2121}& Q_{2122}& 0& 0& G_{21}\\ 0& 0& 0& 0& 0& 0& 0\end{array}\right),$$

$$(\Delta c/c{)}_P$$

$$(\Delta c/c{)}_S$$

$${\mathbf{q}}_{i,j}^{n+1}={\mathbf{q}}_{i,j}^n+\Delta {\mathbf{q}}_{\text{FV}}^n+\Delta t\left(\mathbf{r}({\mathbf{q}}_{i,j}^n)+{\mathbf{s}}_{i,j}^n\right),$$

$${\mathcal{H}}_x{\mathbf{q}}_{i,j}^n$$

$${\mathcal{H}}_y{\mathbf{q}}_{i,j}^n$$

$$\Delta {\mathbf{q}}_{\text{FV}}^n=\frac{1}{2}\left({\mathcal{H}}_x{\mathcal{H}}_y+{\mathcal{H}}_y{\mathcal{H}}_x-2\right){\mathbf{q}}_{i,j}^n,$$

$$\text{Co}=\underset{\begin{array}{c}0⩽i⩽N_x\\ 0⩽j⩽N_y\end{array}}{\max }\max \left\{{\varrho }_{{\mathbf{f}}^{\prime }}({\mathbf{q}}_{i,j}^n)\frac{\Delta t}{\Delta x},{\varrho }_{{\mathbf{g}}^{\prime }}({\mathbf{q}}_{i,j}^n)\frac{\Delta t}{\Delta y}\right\}⩽1,$$

$${\mathbf{A}}_{i+1/2,j}={\mathbf{f}}^{\prime }\left(\frac{1}{2}({\mathbf{q}}_{i,j}^n+{\mathbf{q}}_{i+1,j}^n)\right)$$

$${\mathbf{q}}_{i+1,j}^n-{\mathbf{q}}_{i,j}^n=\sum _{k=1}^7{\alpha }_{i+1/2,j}^k{\mathbf{p}}_{i+1/2,j}^k=\sum _{k=1}^7{\mathcal{W}}_{i+1/2,j}^k,$$

$${\mathbf{f}}_{i+1/2,j}^L$$

$${\mathbf{f}}_{i+1/2,j}^H$$

$$+\frac{1}{2}{c}_S\left(1-\frac{\Delta t}{\Delta x}{c}_S\right)\left(\varphi ({\theta }_{i+1/2,j}^3){\mathcal{W}}_{i+1/2,j}^3+\varphi ({\theta }_{i+1/2,j}^4){\mathcal{W}}_{i+1/2,j}^4\right).$$

$${\theta }_{i+1/2,j}^{1,3}=\frac{{\mathcal{W}}_{i+3/2,j}^{1,3}\cdot {\mathcal{W}}_{i+1/2,j}^{1,3}}{{\mathcal{W}}_{i+1/2,j}^{1,3}\cdot {\mathcal{W}}_{i+1/2,j}^{1,3}},{\theta }_{i+1/2,j}^{2,4}=\frac{{\mathcal{W}}_{i-1/2,j}^{2,4}\cdot {\mathcal{W}}_{i+1/2,j}^{2,4}}{{\mathcal{W}}_{i+1/2,j}^{2,4}\cdot {\mathcal{W}}_{i+1/2,j}^{2,4}},$$

$$\mathbf{q}(x,y,0)=\left\{\begin{aligned} &V\,\big{(}0,0,0,0,\phantom{-\!}\sin\varphi,-\!\cos\varphi,0\big{)}^{\top}&&\text{if}\;x_{\varphi}<0\,,\\ &V\,\big{(}0,0,0,0,-\!\sin\varphi,\phantom{-\!}\cos\varphi,0\big{)}^{\top}&&\text{if}\;x_{\varphi}>0\,.\end{aligned}\right.$$

Full text

††journal: Wave Motion

1 Introduction

Geomaterials such as rocks and concrete exhibit nonlinear features at small strains $\varepsilon\sim 10^{-6}$. In particular, longitudinal resonance experiments reveal the generation of higher-order harmonics. Besides this phenomenon, the material exhibits also a softening with increasing strain amplitudes [1, 2]. After the forcing is stopped, the material recovers gradually its initial stiffness. This softening/recovery process is not instantaneous. The transient regime, named “slow dynamics” in the corresponding literature [3, 4, 5, 6], has a characteristic time much larger than the period of the forcing.

To describe nonlinear elastic behavior, the finite-strain theory is a self-consistent framework. Various constitutive laws of isotropic hyperelastic material express the stress as a function of a strain tensor, and are compatible with Hooke’s law in the infinitesimal strain limit [7, 8, 9]. Among them, the Murnaghan law [10] is frequently used to describe acoustic nonlinearity in rocks and concrete, and values of the parameters can be found in the literature [11, 12, 13]. However, hyperelasticity does not account for the long-time relaxation of elastic constants.

Several models from the literature describe the softening of the material by a dependence of the elastic constants on a scalar variable $g$ [14, 15, 16]. In the present study, the internal-variable model [16] is used, where the stress issued from Murnaghan’s law is multiplied by $1-g$. Hence, $g\equiv 0$ corresponds to the hyperelastic case, and the material softens as $g$ increases. The evolution of $g$ is governed by a first-order differential equation compatible with the principles of thermodynamics. Both softening processes $\dot{g}\geqslant 0$ and recovery processes $\dot{g}\leqslant 0$ are possible [16], accordingly to the experimental observations.

In recent works [17], the corresponding one-dimensional equations of motion have been solved numerically using a finite-volume method with flux limiters [18]. The numerical method has been validated with respect to reference solutions, and qualitative agreement with experimental observations has been obtained. The present article describes a similar method in the two-dimensional plane-strain case. The numerical method is well-suited to the computation of nonlinear waves in the Lagrangian framework, and it can be used for various hyperelastic material models (cf. the related study [19] and references therein).

The article is organized as follows. In Section 2, the equations of motion are detailed, which are rewritten as a nonlinear hyperbolic system of balance laws in two space variables. The numerical method is presented in Section 3. Section 4 shows 2D numerical results illustrating the nonlinear wave propagation in elastic solids with softening. The case of a homogeneous medium is considered, as well as the case of a periodic layered medium with solitary wave solutions.

2 Governing equations

2.1 Lagrangian hyperelasticity with softening

We consider an homogeneous continuum. A particle initially located at some position $\bm{x}_{0}$ of the reference configuration moves to a position $\bm{x}_{t}$ of the deformed configuration. The deformation gradient is a second-order tensor defined by (see e.g. [7, 8, 20, 9])

[TABLE]

where $\bm{u}=\bm{x}_{t}-\bm{x}_{0}$ denotes the displacement field and $\mathbf{grad}$ is the gradient with respect to the material coordinates $\bm{x}_{0}$ (Lagrangian gradient). In the reference configuration, the deformation gradient (1) is equal to the metric tensor $\bm{I}$. Here, the Euclidean space is described by an orthonormal basis $(\bm{e}_{1},\bm{e}_{2},\bm{e}_{3})$ and a Cartesian coordinate system $(O,x,y,z)$, so that the matrix of the coordinates of $\bm{I}$ is the identity matrix.

The Lagrangian representation of motion is used. Hence, the material time derivative $\dot{\bm{F}}=\partial_{t}\bm{F}$ of the deformation gradient satisfies

[TABLE]

where $\bm{v}=\dot{\bm{u}}$ is the velocity field. The conservation of mass implies ${\rho_{0}}/{\rho}=\det\bm{F}$, where $\rho$ denotes the mass density in the deformed configuration, and $\rho_{0}$ denotes the mass density in the reference configuration. Self-gravitation and heat conduction are neglected, so that the motion is driven by the conservation of momentum

[TABLE]

is the first Piola–Kirchhoff stress tensor, and $\mathbf{div}$ denotes the divergence with respect to the material coordinates $\bm{x}_{0}$. The Cauchy stress tensor $\bm{\sigma}=(\det\bm{F})^{-1}\bm{P}\cdot\bm{F}^{\top}\!=\bm{\sigma}^{\top}\!$ is detailed later on through a specification of $\bm{P}$. The term $\bm{f}^{v}$ is an external volume force applied to the material.

In hyperelasticity, the only variables of state are the specific entropy and a strain tensor. Here, the Green–Lagrange strain tensor $\bm{E}=\frac{1}{2}(\bm{F}^{\top}\!\cdot\bm{F}-\bm{I})$ is used, i.e.

[TABLE]

An internal variable $g\in[0,1[$ accounting for the softening of the material is added to the previous list of variables of state. We define the internal energy by unit of reference volume as [16]

[TABLE]

where $e$ is the specific internal energy. In (5), the potential $W(\bm{E})$ is the strain energy function of Murnaghan’s law [10]

[TABLE]

which is expressed as a function of the strain invariants ${E_{\text{I}}}=\text{tr}\,\bm{E}$, {E_{\text{II}}}=\frac{1}{2}\big{(}(\text{tr}\,\bm{E})^{2}-\text{tr}(\bm{E}^{2})\big{)}, and ${E_{\text{III}}}=\det\bm{E}$. The constants $\lambda$, $\mu$ are the Lamé parameters and the constants $\mathfrak{l}$, $\mathfrak{m}$, $\mathfrak{n}$ are the Murnaghan coefficients (third-order elastic constants). In the case where the Murnaghan coefficients in (6) equal zero, the strain energy of the Saint Venant–Kirchhoff model is recovered. The potential $\Phi(g)$ represents a storage energy. Basic requirements are the convexity of $\Phi$ and $\Phi^{\prime}(0)=0$, where $\Phi^{\prime}$ is the derivative of $\Phi$ (see [16]). Moreover, an asymptote at $g=1$ is introduced to avoid the destruction of the material. A suitable expression is

[TABLE]

where $\gamma>0$ is an energy per unit volume. The energy (7) is quadratic $\Phi(g)\simeq\frac{1}{2}\gamma g^{2}$ in the limit $g\to 0$.

Under the previous assumptions, the simplest set of constitutive equations with softening which is thermodynamically admissible reads [16]

[TABLE]

where $\tau_{1}>0$ in J m$-3$ s is a material parameter. The first equation in (8) is the mechanical constitutive law, which reduces to the case of hyperelasticity if $g\equiv 0$. According to the second equation in (8) which governs the evolution of the variable $g$, the classical theory of nonlinear elastodynamics is recovered if $\tau_{1}\to{+\infty}$. Similarly to [16, 17], we rewrite the tensor derivative $\partial W/\partial\bm{E}$ of the strain energy using the invariants’ tensor derivatives:

[TABLE]

where

[TABLE]

In the next subsection, we detail the case of plane strain. The case of uniaxial strain is addressed in [17].

2.2 The plane-strain assumption

The displacement field $\bm{u}$ is independent of $z$, and its component $u_{3}$ along $\bm{e}_{3}$ is zero. In the basis of unit tensors $(\bm{e}_{i}\otimes\bm{e}_{j})_{1\leqslant i,j\leqslant 3}$, the matrix of coordinates of the displacement gradient $\mathbf{grad}\,\bm{u}$ is therefore

[TABLE]

Using the Einstein notation with indices in $\{1,2\}$, the coordinates of the Green–Lagrange tensor (4) write E_{ij}=\frac{1}{2}\big{(}u_{i,j}+u_{j,i}+u_{p,i}u_{p,j}\big{)}. Its invariants in (6) are $E_{\text{I}}=E_{nn}$, E_{\text{II}}=\frac{1}{2}\big{(}{E_{\text{I}}}^{2}-E_{ij}E_{ij}\big{)}=\epsilon_{ij}E_{1i}E_{2j} and $E_{\text{III}}=0$, where $\epsilon_{ij}$ is the Levi-Civita symbol of $\mathbb{R}^{2}$. The Cayley–Hamilton theorem applied to the restriction of $\bm{E}$ to $\mathbb{R}^{2}\times\mathbb{R}^{2}$ reads $E_{im}E_{mj}-E_{\text{I}}E_{ij}+E_{\text{II}}\delta_{ij}=0$, where $\delta_{ij}$ is the Kronecker delta. Hence, the expression (9) of $\partial W/\partial\bm{E}$ becomes

[TABLE]

in the basis of unit tensors $(\bm{e}_{i}\otimes\bm{e}_{j})_{1\leqslant i,j\leqslant 2}$, where

[TABLE]

The components $P_{ij}$ of the first Piola–Kirchhoff stress tensor $\bm{P}$ in (8) are therefore

[TABLE]

under the plane strain assumption, which does not depend upon the third Murnaghan coefficient $\mathfrak{n}$.

When the geometric nonlinearities are negligible, the Green–Lagrange strain tensor (4) is linearized with respect to $\mathbf{grad}\,\bm{u}$, i.e. \bm{E}\simeq\frac{1}{2}\big{(}\mathbf{grad}\,\bm{u}+\mathbf{grad}^{\top}\bm{u}\big{)}=\bm{\varepsilon} reduces to the infinitesimal strain tensor. The coordinates $E_{ij}$ of $\bm{E}$ are replaced by the coordinates \varepsilon_{ij}=\frac{1}{2}\big{(}u_{i,j}+u_{j,i}\big{)} of $\bm{\varepsilon}$. Moreover, the first Piola–Kirchhoff stress tensor $\bm{P}$ is linearized with respect to $\mathbf{grad}\,\bm{u}$ too, i.e. $\bm{F}\cdot\partial W/\partial\bm{E}\simeq\partial W/\partial\bm{\varepsilon}$. Hence, the equation P_{ij}=(1-g)\,\big{(}\tilde{\alpha}_{0}\delta_{ij}+\tilde{\alpha}_{1}\varepsilon_{ij}\big{)} replaces (14). Under this assumption, linear elastodynamics is recovered if $g\equiv 0$ (i.e., $\tau_{1}\to{+\infty}$ in (8)), and if the Murnaghan coefficients $\mathfrak{l}$, $\mathfrak{m}$ in (13) are zero (see A for details).

2.3 System of balance laws

Under the plane-strain assumption, the equations of motion (2)-(3)-(8) are rewritten as a two-dimensional system of balance laws

[TABLE]

where $\mathbf{q}=(u_{1,1},u_{1,2},u_{2,1},u_{2,2},v_{1},v_{2},g)^{\top}\!$ is the vector of unknowns. The expressions of the flux functions, the relaxation function and the source term are

[TABLE]

In (16), the Piola–Kirchhoff stress components $(P_{ij})_{1\leqslant i,j\leqslant 2}$ depend on $(u_{i,j})_{1\leqslant i,j\leqslant 2}$ and $g$ according to (14). The strain energy $W$ depends on $(u_{i,j})_{1\leqslant i,j\leqslant 2}$ according to (6).

The Jacobian matrix of the flux component $\mathbf{f}$ along the $x$-axis is

[TABLE]

where only three strips are displayed (everywhere else, the coefficients in the matrix are zero). The expression of the coefficients $Q_{ijk\ell}$ in (17) defined by $\rho_{0}Q_{ijk\ell}=\partial P_{ij}/\partial u_{k,\ell}$ is detailed in the A, as well as the expression of the coefficients $G_{ij}$ defined by $\rho_{0}G_{ij}=\partial P_{ij}/\partial g$. A similar Jacobian matrix $\mathbf{g}^{\prime}(\mathbf{q})$ is obtained for the flux component $\mathbf{g}$ along the $y$-axis. These matrices are diagonalized in the A. The spectrum of both matrices has the form $\{-c_{P},c_{P},-c_{S},c_{S},0,0,0\}$. In the case of Murnaghan hyperelasticity, the eigenvalues $c_{P}$, $c_{S}$ can be complex [17], so that the system (15)-(16) is not unconditionally hyperbolic (see e.g. [21] for discussions on hyperbolicity in hyperelasticity). Here, we restrict ourselves to configurations where the eigenvalues $c_{P}>c_{S}>0$ are real. Thus, $c_{P}$ and $c_{S}$ correspond to the velocities of compression waves and shear waves, respectively.

Plane waves

We assume furthermore that the displacement field is invariant along a direction, say $\bm{e}_{2}$, so that $\bm{u}$ does not depend on $y$. In this case, the vector of unknown reduces to $\mathbf{q}=(u_{1,1},u_{2,1},v_{1},v_{2},g)^{\top}\!$ — the second and fourth rows of (15)-(16) are zero — and the flux $\mathbf{g}$ along $y$ is zero. The Jacobian matrix $\mathbf{f}^{\prime}(\mathbf{q})$ is obtained from (17) by removing the second and fourth rows, as well as the second and fourth columns. Doing so, two zero eigenvalues are removed from the spectrum, which reduces to $\{-c_{P},c_{P},-c_{S},c_{S},0\}$.

We consider the case of Murnaghan material $g\equiv 0$, with the parameters in Table 1. The latter, found in [13], have been measured on concrete. Let us introduce the relative variation $\Delta c/c$ of the sound velocities $c_{P}$ and $c_{S}$ with respect to the case of Hooke’s law, where $c_{P}=\sqrt{(\lambda+2\mu)/\rho_{0}}\approx 4458$ m/s and $c_{S}=\sqrt{\mu/\rho_{0}}\approx 2700$ m/s. Fig. 1 displays the evolution of $\Delta c/c$ with respect to the compression strain $u_{1,1}$, when the shear strain $u_{2,1}$ is set to zero. One observes that the variations of $c_{P}$ with respect to $u_{1,1}$ are much larger than the variations of $c_{S}$. This is confirmed by the Taylor series approximations

[TABLE]

represented as dotted lines in Fig. 1 (the magnitude of $u_{1,1}$ in the figure has been chosen for graphical reasons). These approximations show also that the shear strain $u_{2,1}$ has much less influence than the compression strain $u_{1,1}$ on the variations of the sound velocities.

3 Numerical resolution

3.1 Numerical strategy

In the examples presented later on, the physical domain is assumed unbounded. We consider a finite numerical domain $[0,L_{x}]\times[0,L_{y}]$. It is discretized using a regular grid in space with mesh size $\Delta x$ in the $x$ direction, and $\Delta y$ in the $y$ direction. The coordinates of the nodes are $(x_{i},y_{j})=(i\,\Delta x,j\,\Delta y)$, where $0\leqslant i\leqslant N_{x}$ and $0\leqslant j\leqslant N_{y}$. The total number of nodes is $(N_{x}+1)\times(N_{y}+1)$, where $N_{x}=L_{x}/\Delta x$ and $N_{y}=L_{y}/\Delta y$. A variable time step $\Delta t=t_{n+1}-t_{n}$ is introduced. Therefore, $\mathbf{q}(x_{i},y_{j},t_{n})$ denotes the solution to (15) at the grid node $(i,j)$ and at the $n$th time step. Numerical approximations of the solution are denoted by $\mathbf{q}_{i,j}^{n}\simeq\mathbf{q}(x_{i},y_{j},t_{n})$.

The non-homogeneous system of balance laws (15) is integrated explicitly in time:

[TABLE]

where the approximation $\mathbf{s}_{i,j}^{n}$ of the source term $\mathbf{s}$ is specified later on. The increment $\Delta\mathbf{q}_{\text{FV}}^{n}$ is deduced from the integration of $\partial_{t}\mathbf{q}+\partial_{x}\mathbf{f}(\mathbf{q})+\partial_{y}\mathbf{g}(\mathbf{q})=\mathbf{0}$ over one time step. Usually, one has $\Delta\mathbf{q}_{\text{FV}}^{n}=(\mathcal{H}_{x}+\mathcal{H}_{y}-2)\,\mathbf{q}_{i,j}^{n}$, where the discrete operators

[TABLE]

involve the fluxes $\mathbf{f}_{i+1/2,j}^{n}$, $\mathbf{g}_{i+1/2,j}^{n}$ of a 2D finite-volume scheme [18]. Here, a second-order symmetric dimensional splitting [22] is used instead. That is to say, $\mathcal{H}_{x}$ and $\mathcal{H}_{y}$ correspond to the integration of $\partial_{t}\mathbf{q}+\partial_{x}\mathbf{f}(\mathbf{q})=\mathbf{0}$ and $\partial_{t}\mathbf{q}+\partial_{y}\mathbf{g}(\mathbf{q})=\mathbf{0}$ over one time step, so that (20) involves the fluxes $\mathbf{f}_{i+1/2,j}^{n}$, $\mathbf{g}_{i+1/2,j}^{n}$ of a 1D finite-volume scheme. The increment $\Delta\mathbf{q}_{\text{FV}}^{n}$ is computed according to

[TABLE]

where $\mathcal{H}_{x}\mathcal{H}_{y}$ denotes the composition of the operators $\mathcal{H}_{x}$ and $\mathcal{H}_{y}$.

The numerical fluxes in (20)-(21) are computed according to the flux-limiter method [18, 17] described in the next subsection. This finite-volume scheme is well-suited for nonlinear wave propagation and second-order accurate. The operators $\mathcal{H}_{x}$ and $\mathcal{H}_{y}$ are stable under the Courant–Friedrichs–Lewy (CFL) condition

[TABLE]

where Co is the maximum Courant number in the $x$ and $y$ directions. The spectral radius $\varrho_{\mathbf{f}^{\prime}}(\mathbf{q})$ of $\mathbf{f}^{\prime}(\mathbf{q})$ corresponds to $c_{P}$ (expression detailed in the A), ditto the spectral radius $\varrho_{\mathbf{g}^{\prime}}(\mathbf{q})$ of $\mathbf{g}^{\prime}(\mathbf{q})$. The stability of the scheme (19) is also restricted by the spectral radius of the Jacobian matrix $\mathbf{r}^{\prime}(\mathbf{q})$. As in 1D [17], the stability limits imply that the scheme (19) is stable under the classical CFL condition (22). Hence, given a spatial discretization and a Courant number $\text{Co}\leqslant 1$, the value of the time step $\Delta t$ is imposed by (22).

3.2 Flux limiter

We describe now the flux-limiter scheme [18, 17]. Since the computation of the numerical fluxes in the $x$ and $y$ directions is similar, only the numerical flux $\mathbf{f}_{i+1/2,j}^{n}$ in the $x$ direction is detailed here. To do so, we introduce the Jacobian matrix

[TABLE]

at the arithmetic mean of the grid node values in the $x$ direction. The jump of the numerical solution $\mathbf{q}_{i+1,j}^{n}-\mathbf{q}_{i,j}^{n}$ along $x$ is decomposed in the basis of right eigenvectors $\{\mathbf{p}^{k}_{i+1/2,j},k=1,\dots,7\}$ of $\mathbf{A}_{i+1/2,j}$,

[TABLE]

which correspond to the eigenvalues $\{-c_{P},c_{P},-c_{S},c_{S},0,0,0\}$ (cf. detailed expressions in the A).

The numerical flux in (20) is the sum of a first-order flux and a second-order limited correction, $\mathbf{f}_{i+1/2,j}^{n}=\mathbf{f}_{i+1/2,j}^{L}+\mathbf{f}_{i+1/2,j}^{H}$, where

[TABLE]

The coefficients $\theta_{i+1/2,j}^{k}$ where $k=1,\dots,4$ express the upwind variation of the jump (24) in the $k$th characteristic field,

[TABLE]

and $\phi$ denotes the minmod limiter function $\phi(\theta)=\max\{0,\min\{1,\theta\}\}$. As such, the weights $\phi(\theta_{i+1/2,j}^{k})$ are designed to avoid spurious oscillations in the numerical solution. Since the eigenvalues indexed by $k=5,\dots,7$ in the decomposition of the jump (24) are zero, the corresponding terms $\bm{\mathcal{W}}^{k}_{i+1/2,j}$ do not appear in the numerical flux (25).

To carry out one iteration in time (19)-(21) at some grid node $(i,j)$, the numerical values of $\mathbf{q}$ at the grid nodes $(i-2,\dots,i+2)\times(j-2,\dots,j+2)$ are required (25). Therefore, two columns and two rows of “ghost cells” are added on the left, the right, the top, and the bottom of the numerical domain. If not specified differently, a zero-order extrapolation of the numerical values is used to update the ghost cell values at each step of (21). This procedure is detailed in Section 21.8 of [18], and is used here to simulate outflow boundary conditions (i.e., an infinite physical domain).

4 Numerical experiments

In the following numerical examples, the Courant number (22) is set to $\text{Co}=0.9$. If not specified otherwise, the physical parameters are given in Table 1. The parameters $\gamma$, $\tau_{1}$ have been chosen so as to obtain significant effects of the softening at the scale of the simulation. The numerical domain is defined by $L_{x}=L_{y}=0.4~{}\text{m}$, and is discretized using $N_{x}=N_{y}=800$ points in each direction.

4.1 Murnaghan hyperelasticity

The first example focuses on nonlinear elastodynamics, i.e., no softening occurs in the material. In (15)-(16), the source term and the relaxation function are removed ($\mathbf{s}=\mathbf{0}$, $\tau_{1}\to{+\infty}$). We consider a Riemann problem with initial data $\mathbf{q}(x,y,0)$, where the material is initially undeformed and opposite transverse velocities with amplitude $V$ are applied:

[TABLE]

The variable $x_{\varphi}=(x-x_{s})\cos\varphi+(y-y_{s})\sin\varphi$ is the $x$-abscissa of a new coordinate system, corresponding to a rotation by an angle $\varphi$ and a translation by $(x_{s},y_{s})$ of the original one. Here, the origin is set at $(x_{s},y_{s})=(L_{x},L_{y})/2$, the rotation angle is $\varphi=15^{\circ}\!$, and the velocity amplitude is $V=0.1$ m/s. To reduce discretization artifacts due to the oblique discontinuity, the average value of (27) over the cell $[x_{i-1/2},x_{i+1/2}]\times[y_{j-1/2},y_{j+1/2}]$ is initially set at the grid node $(i,j)$.

Figures 2 and 3 illustrate the coupling between plane shear waves and plane compression waves in hyperelasticity [9, 23], contrary to linear elasticity where both types of waves are decoupled. Fig. 2 displays a map of $W^{1/8}$ at $t=0.015$ ms, where $W$ is the strain energy (6) obtained numerically with the above method. Fig. 3 display the evolution of the rotated longitudinal velocity $(v_{1})_{\varphi}=v_{1}\cos\varphi+v_{2}\sin\varphi$ and the rotated transverse velocity $(v_{2})_{\varphi}=-v_{1}\sin\varphi+v_{2}\cos\varphi$, along the solid line displayed in Fig. 2. In the case of Hooke’s law of linear elasticity, the solution to the initial-value problem (27) writes

[TABLE]

with $c_{S}=\sqrt{\mu/\rho_{0}}\,$, and only shear waves propagate (solid line in Fig. 3). In the hyperelastic case, Fig. 3b shows that shear waves are generated from the initial data (27), but faster compression waves are also generated. This amplitude-dependent nonlinear effect is better observed at large amplitudes. The small-amplitude perturbations in the numerical solution are caused by the discretization of the oblique discontinuity (27).

4.2 Softening

In this second example, we consider the full system (15)-(16). The material is initially undeformed and at rest, $\mathbf{q}(x,y,0)=\mathbf{0}$, and a volume force $\bm{f}^{v}$ is used for the forcing $\mathbf{s}$. The volume force is an acoustic point source along $x$ with expression $\bm{f}^{v}=A^{v}\sin(2\pi f_{c}t)\,\delta(x-x_{s})\,\delta(y-y_{s})\,\bm{e}_{1}$, where $\delta$ is the Dirac delta function, $A^{v}$ is the amplitude, and $f_{c}$ is the characteristic frequency. Usually, the increment $\mathbf{s}_{i,j}^{n}$ in (19) is obtained by averaging the source term $\mathbf{s}|_{t=t_{n}}$ of (15)-(16) over the cell $[x_{i-1/2},x_{i+1/2}]\times[y_{j-1/2},y_{j+1/2}]$. Here, we approximate the Dirac deltas by a truncated Gaussian function to avoid strain concentration at the source. Thus,

[TABLE]

where $d_{ij}$ is the distance between $(x_{i},y_{j})$ and $(x_{s},y_{s})$. Denoted by the indicator function $\mathbf{1}_{d_{ij}\leqslant R}$, the function’s support is a disk with radius $R=c_{P}/(7.5\,f_{c})$, where $c_{P}=\sqrt{(\lambda+2\mu)/\rho_{0}}$ is the speed of linear compression waves. The width parameter of the Gaussian function is chosen such that $\sigma_{c}=R/2$. The point load has amplitude $A^{v}=0.5$ kN/m, frequency $f_{c}=100$ kHz, and it is located at the nearest grid node of the domain’s center: $(x_{s},y_{s})\simeq(L_{x},L_{y})/2$. The source (29) is switched on at $t=0$, and switched off at $t=0.04$ ms. Two receivers R1-R2 record the numerical solution during the simulation. R1 is located at $(x_{r},y_{r})=(0.2,0.22)$ m, and R2 is located at $(x_{r},y_{r})=(0.2,0.27)$ m.

Figure 4 illustrates the effect of the softening on the wave propagation. Fig. 4a displays a map of the strain energy $\left(1-g\right)W$ at the time $t=0.04$ ms, which shows the propagation of cylindrical waves. Denoted by a bullet point and by a square, the receivers R1-R2 are located in a region of the plane where mainly shear waves propagate. Figs. 4b-4c show the effect of the softening at the position of the receivers for several forcing amplitudes $A^{v}$. One observes that $g$ increases while the wave passes by the receiver, and that it relaxes towards zero afterwards. This softening/recovery process is all the more important as the forcing amplitude is large. The characteristic time of the slow dynamics $\tau_{1}/\gamma\approx 0.05$ ms corresponds to the characteristic time of the recovery [16, 17]. In Figs. 4b-4c, one observes the distortion of the velocity signal during propagation, and its delay due to the increase of $g$. The recorded signals are similar to experimental ones obtained in a longitudinal configuration [24].

4.3 Periodic layered medium

Similarly to [25, 26, 27], a periodic layered medium is considered. The basis vector $\bm{e}_{1}$ which orientates the $x$-axis is normal to the interfaces. The layers have the same thickness $d=1$ cm and perfect contact is assumed between neighbor layers (continuity of the displacement and of the stress). As illustrated in Fig. 5, the linear material properties $\Pi\in\{\rho_{0},\lambda,\mu\}$ vary in space according to

[TABLE]

where $\bar{\Pi}\in\{\bar{\rho}_{0},\bar{\lambda},\bar{\mu}\}$ denotes the reference parameters given in Table 1. The nonlinear material properties $\mathfrak{l}$, $\mathfrak{m}$, $\gamma$, $\tau_{1}$ do not vary in space. A velocity pulse with frequency $f_{c}=10$ kHz and amplitude $V=1$ m/s is imposed at the boundary $x=0$:

[TABLE]

Since the configuration is invariant along $y$, the problem is one-dimensional. The wavelength of compression waves $c_{P}^{\text{eff}}/f_{c}\approx 39~{}\text{cm}\gg d$ is deduced from the effective sound speed

[TABLE]

in the linear layered medium.

The numerical domain defined by $L_{x}=1$ m is discretized using $N_{x}=3000$ points, i.e. each layer is represented by 30 points. The boundary condition (31) is implemented accordingly to Sec. 7.3.4 of [18] up to $t=1/f_{c}$. For $t>1/f_{c}$, periodic boundary conditions are applied (see Sec. 7.1 of [18]). A modification of the numerical method described in Sec. 3 is introduced to account for spatially-varying coefficients (B).

Figure 6 illustrates the effect of spatially-varying coefficients on the numerical solution up to $t=2.5$ ms. Figs. 6a-6b display seismograms obtained by unwrapping the velocity signals $v_{1}(0,y,t)$ recorded at the abscissa $x=0$. The “offset” corresponds to the effective propagation distance from the abscissa $x=0$ where the source (31) is imposed. Fig. 6a is obtained with the Murnaghan parameters $\mathfrak{l}$, $\mathfrak{m}$ in Table 1 while the softening variable $g$ is equal to zero (no softening: $\tau_{1}\to{+\infty}$). Wavefront steepening is observed, leading to a series of oscillations. Fig. 6b shows the same output when the softening variable $g$ evolves according to (8), with the parameters $\gamma=10^{5}$ $\text{J}\,\text{m}^{-3}$ and $\tau_{1}=5.0$ $\text{J}\,\text{m}^{-3}\,\text{s}$. The parameters $\gamma$, $\tau_{1}$ are chosen such that the variable $g$ reaches values of 3% during the simulation, while the characteristic time $\tau_{1}/\gamma=5\times 10^{-5}$ s is the same as in Table 1. With respect to the case without softening (Fig. 6a), a longer propagation distance is needed until the wave separates into a series of oscillations, and the amplitudes of oscillations are modified.

The speed of the largest oscillations over the range of propagation distances $[6,8]$ m is deduced from Figs. 6a-6b and reported in Fig. 6c, where the abscissa is the amplitude of each oscillation at the propagation distance $7$ m. This procedure is repeated with a smaller input amplitude $V=0.9$ m/s. A nearly linear evolution of the propagation speed with respect to the amplitude is observed, as is the case for solitary waves. At small amplitudes, the propagation speed is close to the effective speed of sound $c_{P}^{\text{eff}}$ (32). Similar observations are reported in [27] where a different constitutive law is used. Numerically, Fig. 6c shows that the softening does not suppress the solitary waves, but it modifies the relationship between their amplitude and their propagation speed.

5 Conclusion

Within the Lagrangian finite-strain theory, the constitutive model used in this study expresses the stress as a function of the Green–Lagrange strain tensor and a softening variable $g$. Also, an evolution equation for $g$ is provided. The system of partial differential equations so-obtained writes as a nonlinear hyperbolic system of balance laws, so that finite-volume methods can be applied. If the softening is neglected, then the material follows Murnaghan’s law, where a shear wave excitation induces the propagation of smaller-amplitude compression waves. Otherwise, the propagation of a perturbation is responsible for the softening of the material, which recovers gradually its initial stiffness after excitation has stopped. Several numerical examples are considered along the paper. The latter have been chosen to serve modeling purposes rather than experimental ones, related to the observation of hysteresis and long-time relaxation in geophysics and non-destructive testing. In a periodic layered medium, solitary waves are numerically observed, as known in the case of nonlinear elastodynamics.

Now, let us mention possible future works. In the Lagrangian plane-strain case, the numerical method presented here can be used for various hyperelastic constitutive models, such as neo-Hookean, Mooney–Rivlin, Ogden, etc. [7, 8]. For such models, it would be interesting to solve analytically the symmetric shear impact problem in Fig. 2 for validation purposes. Variations of the evolution equation of $g$ in (8) can be considered as well. The influence of the softening on the smoothness of solutions could be investigated numerically and theoretically. To go further towards realistic configurations, viscoelastic attenuation should be accounted for, e.g. in a similar fashion to [28]. Concerning the nonlinear layered medium, it would be interesting to derive the corresponding nonlinear dispersive wave equations by homogenization, as done in [29, 27, 30] and related works. For the development of higher-order methods such as ENO or WENO schemes [31], the eigendecomposition of the fluxes provided in the A is a useful result.

Acknowledgments

This work was supported by the interdisciplinary mission of CNRS (INFINITI). The project leading to this publication has received funding from Excellence Initiative of Aix-Marseille University - A*MIDEX, a French “Investissements d’Avenir” programme. It has been carried out in the framework of the Labex MEC. Opération réalisée avec le concours du Programme d’investissements d’avenir du Gouvernement français dont la gestion a été confiée à l’Andra.

Appendix A Jacobian matrices of the flux

A.1 Expression of the coefficients

We use the Einstein summation convention with indices in $\{1,2\}$. To encompass both cases with and without geometric nonlinearity in a single equation, we introduce a parameter $\Theta\in\{0,1\}$ such that $\Theta=1$ corresponds to finite strain and $\Theta=0$ corresponds to infinitesimal strain. Hence, the coordinates of the strain tensor (4) are written $E_{ij}=\varepsilon_{ij}+\frac{1}{2}\Theta\,u_{p,i}u_{p,j}$, where $\varepsilon_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i})$ are the coordinates of the infinitesimal strain tensor. Moreover, the components of the Piola–Kirchhoff stress tensor (14) are written P_{ij}=(1-g)\left(\delta_{im}+\Theta\,u_{i,m}\right)\big{(}\tilde{\alpha}_{0}\delta_{mj}+\tilde{\alpha}_{1}E_{mj}\big{)}, i.e.

[TABLE]

where $\tilde{\alpha}_{0}$, $\tilde{\alpha}_{1}$ are given in (13). The coefficients (17) of the Jacobian matrices $\mathbf{f}^{\prime}(\mathbf{q})$ and $\mathbf{g}^{\prime}(\mathbf{q})$ satisfy $\rho_{0}Q_{ijk\ell}=\partial P_{ij}/\partial u_{k,\ell}$ and $\rho_{0}G_{ij}=\partial P_{ij}/\partial g$. In the present case of Murnaghan material with softening (33), one has

[TABLE]

where

[TABLE]

The case of Hookean solids is recovered if $g\equiv 0$, geometric nonlinearity is neglected ($\Theta=0$), and the Murnaghan coefficients $\mathfrak{l}$, $\mathfrak{m}$ are zero. In this case, Eq. (34) gives $\rho_{0}Q_{ijk\ell}=\lambda\delta_{ij}\delta_{k\ell}+\mu(\delta_{ik}\delta_{j\ell}+\delta_{jk}\delta_{i\ell})$.

A.2 Eigendecomposition

We provide an eigendecomposition of the Jacobian matrices $\mathbf{f}^{\prime}(\mathbf{q})$ and $\mathbf{g}^{\prime}(\mathbf{q})$ of the fluxes. The hyperelastic case without softening is recovered by removing the last row and the last column of each matrix in the following paragraphs.

Flux along the $x$-axis

The Jacobian matrix (17) of $\mathbf{f}$ at the linear average (23) is diagonalized. Let us write $\mathbf{A}_{i+1/2,j}=\mathbf{P}\mathbf{\Lambda}\mathbf{P}^{-1}$ where $\mathbf{P}$ is an invertible real matrix, and $\mathbf{\Lambda}$ is a diagonal real matrix. The matrix of eigenvalues $\mathbf{\Lambda}=\text{diag}(-c_{P},c_{P},-c_{S},c_{S},0,0,0)$ satisfies

[TABLE]

where the plus sign gives the expression of $c_{P}$ (compressional waves), and the minus sign gives the expression of $c_{S}$ (shear waves). The first four right eigenvectors $\mathbf{p}^{k}_{i+1/2,j}$ of $\mathbf{A}_{i+1/2,j}$ used in (24)-(25) are the first four columns of $\mathbf{P}$, where

[TABLE]

with the coefficients

[TABLE]

The matrix $\mathbf{P}$ is invertible provided that its determinant is nonzero, i.e. $Q_{1121}Q_{2111}\neq Q_{1111}Q_{2121}$ and $Q_{1121}Q_{2111}\neq-\frac{1}{4}(Q_{1111}-Q_{2121})^{2}$. Let us consider each equality case:

if $Q_{1121}Q_{2111}=Q_{1111}Q_{2121}$, then the eigenvalues of $\mathbf{f}^{\prime}(\mathbf{q})$ satisfy $c_{S}=0$. Therefore, the reduced system of conservation laws for plane waves propagating along $x$ is not strictly hyperbolic (eigenvalues $\{{-c_{P},c_{P},-c_{S},c_{S},0}\}$);

• 2.

if $Q_{1121}Q_{2111}=-\frac{1}{4}(Q_{1111}-Q_{2121})^{2}$, then the eigenvalues of $\mathbf{f}^{\prime}(\mathbf{q})$ satisfy $c_{P}=c_{S}$, which is impossible for the same reason.

Therefore, the previous eigendecomposition is valid over the domain of strict hyperbolicity. The first four left eigenvectors $\mathbf{l}_{i+1/2,j}^{k}$ of $\mathbf{A}_{i+1/2,j}$ are the first four rows of $\mathbf{P}^{-1}$, where

[TABLE]

with the coefficients

[TABLE]

The coefficients $\alpha^{k}_{i+1/2,j}$ in (24) are equal to the scalar products $\alpha^{k}_{i+1/2,j}=\mathbf{l}^{k}_{i+1/2,j}\cdot\left(\mathbf{q}_{i+1,j}^{n}-\mathbf{q}_{i,j}^{n}\right)$.

Flux along the $y$-axis

Similarly to (23), we introduce the Jacobian matrix $\mathbf{g}^{\prime}(\mathbf{q})$ of $\mathbf{g}$ at the linear average $\mathbf{B}_{i,j+1/2}=\mathbf{g}^{\prime}\left(\frac{1}{2}(\mathbf{q}_{i,j}^{n}+\mathbf{q}_{i,j+1}^{n})\right)$, and provide an eigendecomposition $\mathbf{B}_{i,j+1/2}=\mathbf{P}\mathbf{\Lambda}\mathbf{P}^{-1}$. The matrix of eigenvalues $\mathbf{\Lambda}=\text{diag}(-c_{P},c_{P},-c_{S},c_{S},0,0,0)$ satisfies

[TABLE]

where the plus and minus signs give the expressions of $c_{P}$ and $c_{S}$, respectively. With similar notations as (36), we have

[TABLE]

with the coefficients

[TABLE]

A similar analysis shows that $\mathbf{P}$ is invertible over the domain of strict hyperbolicity:

[TABLE]

with the coefficients

[TABLE]

Appendix B Spatially-varying coefficients

The problem described in Fig. 5 is one-dimensional. Thus, the variable $y$ and the index $j$ do not appear in the present description of the numerical method, which is a modification of (25) to account for spatially-varying coefficients [25, 26, 27] (cf. Chap. 9 of [18]). The conservation laws (15) are rewritten in terms of ${\bf q}=(u_{1,1},\rho_{0}v_{1},g)^{\top}\!$. The modified method amounts to choosing the eigenvalues $\{{-c_{P}},c_{P},0\}$ and corresponding eigenvectors ${\bf p}_{i+1/2,j}^{k}$, $k\in\{1,3\}$ of ${\bf A}_{i+1/2}$ in a downwind fashion. Hence, $c_{P}$ is computed at the abscissa $x_{i}$ for $k=1$ and at the abscissa $x_{i+1}$ for $k=2$. Since both left-going and right-going waves do not have equal absolute speeds at the cell interface $x_{i+1/2}$, writing the low-order part of the algorithm (25) as a flux difference is no longer possible: the “wave-propagation form” is used instead. To guarantee the continuity of the velocity and of the stress at the material interfaces, a modified expression of the coefficients $\alpha^{k}_{i+1/2}$ defining the vectors $\bm{\mathcal{W}}^{k}_{i+1/2}$ is used. In particular, the sum of these vectors does not equal the jump ${\bf q}_{i+1}-{\bf q}_{i}$ anymore, as was the case in a homogeneous medium (24). To avoid numerical instability, a transmission-based limiter is implemented, which modifies the expression (26) of the coefficients $\theta^{k}_{i+1/2}$.