A hyperbolic framework for shear sound beams in nonlinear solids

In soft elastic solids, directional shear waves are in general governed by coupled nonlinear KZK-type equations for the two transverse velocity components, when both quadratic nonlinearity and cubic nonlinearity are taken into account. Here we consider spatially two-dimensional wave fields. We propose a change of variables to transform the equations into a quasi-linear first-order system of partial differential equations. Its numerical resolution is then tackled by using a path-conservative MUSCL-Osher finite volume scheme, which is well-suited to the computation of shock waves. We validate the method against analytical solutions (Green's function, plane waves). The results highlight the generation of odd harmonics and of second-order harmonics in a Gaussian shear-wave beam.


Introduction
Understanding the propagation and diffraction of sound emitted by a directional source has been an important concern of the nonlinear acoustics community since the late 1960s [1]. For this purpose, the celebrated Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation was derived from the equations of fluid dynamics by introducing an appropriate scaling. Valid in the paraxial region of a directive acoustic source (e.g., a transducer), this nonlinear parabolic equation describes how sound beams spread with increasing propagation distance, incorporating harmonic generation and attenuation due to nonlinearity and dissipation effects. The same approach was used for the equations of nonlinear Lagrangian elastodynamics, leading to similar partial differential equations (see the review by Norris [2]).
In soft incompressible solids, the experimental observation of nonlinear shear waves has been reported in the literature [3,4]. Along with these observations, the generation of mainly odd harmonics and shocks has also been reported. To explain these features, Zabolotskaya et al. [5] showed that plane shear waves with a single transverse displacement component are governed by a Burgers-like equation when cubic nonlinearity is taken into account and quadratic nonlinearity is ignored. In that case, the shear waves are linearly polarized and the motion is purely antiplane.
However, incompressible solids with cubic nonlinearity only, and no quadratic nonlinearity, are modelled by a very special constitutive law [6], not representative of real-world materials. This limitation is resolved by considering transverse shear waves with an arbitrary polarization. Then, two coupled KZK-type equations are obtained [7]. If the corresponding wave fields have variations in the transverse direction (that is, if the plane-wave assumption is relaxed), then the governing equations of motion include both quadratic and cubic nonlinearity. Related works show that the second harmonic can be generated in this configuration [8], which is in agreement with more recent measurements [9].
No analytical solution is known for this system of coupled nonlinear partial differential equations, and most of the above-mentioned studies rely on quasi-analytical approaches. Hence, Zabolotskaya et al. [5] estimate the generation of harmonics by using a space-dependent harmonic expansion, and by performing harmonic balance. Wochner et al. [7] use a similar harmonic expansion along with Green's function expansions. Finally, Destrade et al. [8] implement a perturbation method based on a small amplitude parameter.
While analytical results are of great interest, computational approaches may be more versatile. The numerical resolution of KZK-type equations was addressed by Hamilton et al. [10], by matching near-field and far-field Fourier series expansions (Bergen code). Other frequency-domain approaches [11,12] were then followed by time-domain methods, in particular by making use of operator splitting [13]. Pinton and Trahey [14] combined operator splitting with shock-capturing Godunov-type methods to provide accurate shock-wave solutions. To the present authors' knowledge, no method has yet been successful in solving the quadratic-cubic nonlinear system [7] describing directional shear-wave motion in soft elastic solids.
In this article, we consider spatially two-dimensional wave fields. After a brief presentation of the governing equations, we introduce a change of unknowns and of dependent variables that transforms the system at hand into a quasi-linear system of first-order partial differential equations (Section 2). In particular, the physical time variable is used instead of the retarded time. Since the differential system so-obtained is non-conservative, particular care is required when computing shock-wave solutions. Indeed, a naive upwind scheme would lead to inaccurate wave speeds [15]. In this study, we implement a finite volume method based on the path-conservative Osher Riemann solver and on MUSCL reconstruction [16,17] (Section 3). Although it does not involve operator splitting, the scheme accounts naturally for nonlinearity, coupled motion and beam diffraction. We validate the method by using dedicated analytical solutions which are summarized in the C (Green's function, plane waves). Numerical simulations of directional wave beams illustrate the generation of odd and second-order harmonics (Section 4), as predicted by Destrade et al. [8].
One benefit of the first-order formulation introduced in this study is the potential to use advanced computational methods, including high-order adaptive schemes based on ADER or WENO approaches (see Refs. [18,19] and references therein). Moreover, such differential systems of hydrodynamic type have been studied extensively, and dedicated integrability criteria are known [20]. Prospective applications encompass the study of traumatic brain injury [9], as well as related medical imaging techniques.

Governing equations
We introduce the deformation gradient tensor F = ∂x/∂X , where x represents the position of a particle in the deformed configuration, and X represents its position in the undeformed configuration [2,21,22]. The Lagrangian specification of motion is used throughout the present document, so that spatial differential operators are always computed with respect to X . The components X = (X , Y , Z ) of the position are expressed with respect to an orthonormal basis (e 1 , e 2 , e 3 ) of the Euclidean space, and a Cartesian coordinate system is chosen. Introducing the displacement field u = x − X = (u 1 , u 2 , u 3 ), we write the deformation gradient as F = I + grad u, where I is the identity tensor. Consequently, we also have where v = ∂ t u is the particle velocity.
In this paper we consider incompressible hyperelastic materials, for which the constraint of no volume dilatation is prescribed at all times, so that the mass density ρ is constant. The deformation is also governed by the equation of motion [2,21,22] where f is the density of body force per unit volume. The dependence of the first Piola-Kirchhoff stress tensor P with F is specified by the constitutive law. For incompressible solids, the constitutive law may be expressed as P = −pF − + ∂W /∂F , where p is a Lagrange multiplier due to incompressibility and W is the strain energy density. For instance, the strain energy of homogeneous and isotropic incompressible solids may be expanded as [7,8] in terms of the invariants I k = tr E k of the Green-Lagrange strain tensor This finite-strain tensor is linked to the right Cauchy-Green deformation tensor C = F F through the relation E = 1 2 (C − I ). The material parameters of (4) are the initial shear modulus µ (second Lamé coefficient) and the higherorder elastic constants A, D (Landau constants of quadratic and cubic nonlinear elasticity, respectively). The constitutive law may be rewritten as P = F S, where S = −pC −1 + ∂W /∂E is the second Piola-Kirchhoff stress tensor. Note that when computing the tensor derivative ∂W /∂E , we must keep in mind that the incompressibility constraint introduces a dependence of one invariant I k with respect to the two others (see [23] and A). In the Appendix, we derive the expression of the Cauchy stress tensor σ as in terms of B = F F , the left Cauchy-Green deformation tensor, from which the expression of P = σF − is deduced. The coefficients in Eq. (6) are detailed in the Appendix (Eq. (32)). Using the Cayley-Hamilton theorem, the previous stress-strain relationship may be written in terms of I , B and B −1 up to a redefinition of the arbitrary Lagrange multiplier [8]. However, we keep the present form to avoid the computation of inverse matrices when working out the equations of motion later on.

Scaling the equations of motion
Z X Y Figure 1: Typical directional wave fields of this study; Sketch of equal-phase surfaces issued from a line source located on the Z -axis.
Similarly to Wochner et al. [7], we introduce the following scaling where is a small dimensionless parameter and c = µ/ρ is the shear wave speed of linear elasticity. No scaling is assumed for the Z coordinate, as the field variables are assumed invariant with respect to Z (see Fig. 1). Here, the fields u, p, f are functions of the coordinates X , t , while the new variables U ,p,f depend onX ,t . Note that the present paraxial approximation for directive sources differs significantly from the geometric acoustics/optics approximation (ray theory) [24], even though formal similarities may be found. Following the transformation rules (7), we rewrite the variables F , v and the equations of motion (1)-(3) as in terms of the new coordinates, where the Jacobian matrices have the following components: In terms of the displacement field U , we have and ρ 2 U 1,tt = 2 P 11,1 + P 12,2 − 1 c P 11,t + 3f where the components of P are deduced from the scaled components of the deformation gradient tensor F . Here, partial differentiation is specified using subscript notation (after the commas). Integer subscripts1-3 denote partial differentiation with respect to the components of the position vectorX , while the subscriptt denotes differentiation in time. At leading (quadratic) order in , the incompressibility constraint (2) amounts to the substitution U 1,t = cU 2,2 in the above equations. By transforming back to the original spatial coordinates X and by making appropriate substitutions, the same equations as Eqs. (10)-(11) of [7] are obtained at cubic order as but with additional spatial symmetries due to invariance along the Z -coordinate. Here, the displacement u depends on the coordinates X ,t , which is standard but slightly abusive notation compared to the initial definitions (7). At leading (quadratic) order, the equation giving the Lagrange multiplier in unbounded domain reads p = ρβ 2 u 2,t 2 + u 3,t 2 . In Eq. (11), we introduced the quadratic terms coefficient β 2 and the cubic terms coefficient β 3 [7,8] For later use, we introduced a parameter α ∈ {0, 1} in Eq. (11) that gives the possibility to discard the diffraction term. Note that if the configuration is invariant along the transverse Y -axis, then the diffraction term vanishes as well as the quadratic term [4, 5].

First-order recast
Let us introduce the displacement's partial derivatives v i = u i ,t , θ i = u i ,1 , ε = u 2,2 and ϑ = u 3,2 with i = 2, 3. Thus, the system (11) is rewritten as By making use of the equality of mixed partial derivatives, four kinematic relationships between the strains u i , j and the velocities u i ,t are derived. The equations form a first-order PDE system aq ,1 +bq ,2 +c(q)q ,t = s in terms of the vector of unknowns q = (v 2 , θ 2 , ε, v 3 , θ 3 , ϑ) . The matrix c(q) = c L +c NL (q) is decomposed as the sum of a constant part c L and of a non-constant part c NL (q), which vanishes if the parameters of nonlinearity β 2 , β 3 are zero. The corresponding matrices are detailed in the B. Note that the matrix c NL (q) does not depend on θ 2 , θ 3 , and that it vanishes if q → 0. In Eqs. (11)-(13), the unknowns u, q are functions of X andt (explicit dependence has been dropped for sake of conciseness). Now, the transformation from the retarted timet to the real time t is carried out. Thus, we introduce the displacement field u(X ,t ) = η(X , t ). Using differentiation rules, one shows that q satisifes v i = η i ,t and θ i = η i ,1 + v i /c, while the differential definitions of ε, ϑ are the same whether u or η is used. Next, we introduce the vector The transformation matrix T is defined in such a way that the variables q = T p are modified according to The first-order PDE system reads a p ,1 +b p ,2 +c (q)p ,t = s in terms of the original time variable t , where the prime denotes right-multiplication by T and where c (q) = c L + 1 c a + c NL (q). For later use, we compute the determinant ∆ of the matrix c (q), and the expression is obtained. In particular, we see that this expression is nonzero even if β 2 and β 3 are both equal to zero, or if (v 2 , ε, v 3 , ϑ) is sufficiently close to zero. As long as the matrix c (q) is not singular, we can left-multiply our system by c (q) −1 to rewrite the equations of motion as a quasi-linear first-order system of balance laws, We give the expressions of the above matrices and vectors in B.
Hyperbolicity Let us assume that β 2 and β 3 both equal zero, so that the system (11) becomes linear. In this case, the displacement fields u 2 and u 3 decouple. Moreover, they satisfy the same linear partial differential equation u ,1t = c 2 (α 2 u ,22 + f /µ). According to the theory of second-order differential equations in three independent variables (see e.g. [25], Sec. III.3.1), this equation is hyperbolic. If the real time variable t is used instead of the retarded timet , then the previous equation may be rewritten as (u ,t + cu ,1 ) ,t − 1 2 α 2 c 2 u ,22 = g where g = 1 2 f /ρ. We then recognize a modified wave equation in (Y , t ) coordinates, where the classical term u ,t t has been replaced by the time derivative of a transport term along X . Analytical solutions to this equation are detailed in C. Now we study the system of balance laws (15) where the coefficients β 2 , β 3 are set to zero. The spectral properties of the matrix A show that the characteristic speed along the X -axis is +c. Similarly, the spectral properties of B yield the characteristic speeds ±αc/ 2 along the Y -axis. Moreover, the hyperbolicity property can be deduced from the spectrum of the linear combination n 1 A + n 2 B where n = (n 1 , n 2 ) is a unit vector [15,26]. This property generalises to the quasi-linear system (15) with arbitrary β 2 , β 3 , as long as (v 2 , ε, v 3 , ϑ) stays in the vicinity of the origin. More detailed conditions can be derived analytically in the particular case β 2 = 0, β 3 > 0, where the non-singularity of the determinant (14) implies v 2 2 + v 3 2 < c 2 /β 3 . If this condition is satisfied, then the eigenvalues of n 1 A(p) + n 2 B(p) are real for any unit vector n in both cases α ∈ {0, 1}, so that hyperbolicity is ensured.
We note that the equations of motion have a structure similar to Lagrangian elastodynamics. Indeed, the full system of nonlinear elasticity reads as a conservative first-order system [27], where the conserved variables p are the displacement gradients and the velocities. Nevertheless, the number of state variables is smaller than in the case of two-dimensional elastodynamics [28] due to assumed deformation and symmetries. An additional remarkable feature is the one-way nature of the wave propagation, which results from the scaling procedure. In the present directional wave beam model, a non-conservative system is naturally obtained. Specific theoretical and numerical difficulties arise with such quasi-linear systems of partial differential equations, due to the presence of the non-conservative products A(p) p ,1 and B(p) p ,2 . A dedicated numerical method is presented in the next section.

Numerical resolution
In the examples presented later on, the physical domain is assumed unbounded. We consider a finite numerical domain for X = (X , Y ) in [−1, 1] × [−1, 1]/ 2. It is discretized using a regular grid in space with mesh size ∆x in the X -direction, and ∆y in the Y -direction. The coordinates of the nodes are (x i , y j ) = (i ∆x, j ∆y), where 0 i N x and 0 j N y . The total number of nodes is (N x + 1) × (N y + 1), where N x = 2/∆x and N y = 2/∆y denote the number of cells in each direction. A variable time step ∆t = t n+1 − t n is introduced. Therefore, p(x i , y j , t n ) denotes the solution to (15) at the grid node (i , j ) and at the nth time step. Numerical approximations of the solution are denoted by p n i , j p(x i , y j , t n ).

Finite volume method
The system of balance laws (15) is integrated explicitly according to the following updating formula: where the approximation S n i , j of the source term S(p) is specified later on. This formula is used for the interior cells, while pseudo-absorbing boundary conditions are implemented at the boundaries of the numerical domain (Sec. 21.8.5 p. 488 of [15]).
As described in the next paragraph, the jumps Df ± i ∓1/2, j , Dg ± i , j ∓1/2 and the corrections δf i , j , δg i , j in (16) are computed according to a path-conservative finite volume method with slope limiters [16][17][18]. The method avoids spurious oscillations and is nearly second-order accurate in space and time on smooth solutions. Moreover, it does not suffer the severe limitations of the "naive" non-conservative upwind method regarding non-smooth solutions (see [15], p. 238). Note that due to the nonlinearity of the system, shocks might form even if loadings are smooth.
The numerical flux differences are computed by applying the MUSCL-Hancock procedure componentwise [17], in combination with a path-conservative Osher scheme [16]. The method consists of the following steps: 1. We construct the left ('−') and right ('+') linearly extrapolated values at the cell interfaces (x i +1/2 , y j ) and (x i , y j +1/2 ) as follows (Sec. 14.4.2 and 16.5 of [17]): The coefficients of the matrices A(p) and B(p) are provided in B. The vectors are limited differences of p along the X -and Y -directions -in other words, they are limited slopes multiplied by the mesh size. Here, the monotonized central-difference limiter defined by MC(a, b) = 1 2 (sgn a + sgn b) min 2 |a|, 2 |b|, 1 2 |a + b| (19) is applied componentwise, where sgn denotes the sign function [15].
2. The evaluation of the jumps is performed according to a path-conservative Osher scheme [16] Df ± i + 1 2 , j ) , where a linear integration path is used. A spectral decomposition of the matrix is obtained numerically (as well as a spectral decomposition of B(p)), e.g. by using the eigen function of the Julia Language [29]. Here, R(p), Λ(p) denote the matrices of right eigenvectors of A(p) and corresponding eigenvalues, respectively. The matrices A ± (p) = R(p) Λ ± (p) R(p) −1 (22) and B ± (p) are obtained by taking the positive part '+' or the negative part '−' of the eigenvalues to ensure correct upwinding. The corrections are designed to ensure consistency with conservative methods [18]. Note that these important terms are not included in Ref. [16]. Similarly to (20), a linear integration path was used. The integrals (20)-(23) are computed numerically by using the three-point Gauss-Legendre quadrature rule [16], which proved sufficient to reach second-order accuracy.
Fromm-type methods are recovered by replacing the MC function (19) with the linear average (a, b) → 1 2 (a + b), while the projections (a, b) → a or b yield either Beam-Warming or Lax-Wendroff-type methods [17]. The first-order pathconservative Osher scheme is recovered by replacing the limiter function (19) with the zero function. We note that the non-conservative upwind scheme is then recovered if we choose to evaluate the integrals of Eq. (20) by using downwind-biased Riemann sums.
Empirically, the method is observed to be stable under the Courant-Friedrichs-Lewy (CFL) condition where Co is the maximum Courant number in the X and Y directions. The spectral radii A = max |Λ| and B are deduced from the spectral decomposition (21) of the system's matrices. The stability of the scheme (16) is also restricted by the spectral radius of the Jacobian matrix S (p). For the examples presented hereinafter, a comparison of the stability limits implies that the scheme (16) is stable under the classical CFL condition (24). Hence, given a spatial discretization and a Courant number Co, the value of the time step ∆t is updated at each iteration according to Eq. (24).

Validation
The previous method is applied to a set of test cases with analytical solutions are detailed in C. The first two tests correspond to the two-dimensional linear case where α = 1, β 2 = 0, β 3 = 0, and ρ and µ are taken in Table 1. The third test is performed in a one-dimensional nonlinear case where α = 0, β 2 = 0 (cubic nonlinearity only), and the other parameters are specified in Table 1. The values in Table 1 are representative of pig brain matter [8,30]. Unless stated otherwise, the Courant number is Co = 0.45.
Linear initial-value problem The source term S n i , j is zero. The initial data p 0 i , j with wavelength λ = 0.2 m is obtained by computing the cell averages of in m/s. Thus, the initial data consists of a smooth sinusoidal bump and a rectangular bump, which propagate along the direction of angle ϕ = π/4 (see C). In theory, the initial data is translated diagonally with constant speed 1+ 3 2 2 c. This is illustrated in Fig. 2, which displays the numerical solution. It was obtained by iterating the time-stepping formula (16) up to t ≈ 0.1 s. Fig. 2a displays the numerical solution obtained with N x = N y = 200. This figure shows that both parts of the wave are captured well. Fig. 2b shows error measurements in L 2 -norm performed along the line Y = 0 for −0.05 < X < 0.2 (smooth bump), where N x = N y varies from 25 to 800. The experimental order of accuracy is evidenced by the slope of the error curve. With the present method, second-order accuracy is obtained.
Linear non-homogeneous problem The initial data p 0 i , j is zero. The only non-zero component of the body force f is f 2 = 2ρg (N/m 3 ). The corresponding non-zero component of the source term S is of the form S 1 = g (see B). Here, we consider a sinusoidal point source g (X , t ) = aδ(X )δ(Y )s(t ) with s(t ) = sin(ωt ) for times t > 0. Computing the cell averages of S, we have Indeed, the cell average of Dirac deltas δ(X )δ(Y ) produces Kronecker symbols δ i ,i 0 δ j , j 0 divided by the cell's surface area ∆x∆y. Here, the source is localised at the origin, i.e. x i 0 = 0 and y j 0 = 0 when the numbers of cells N x , N y are even integers. The source has amplitude a = 0.2 m 3 /s 2 and angular frequency ω = 20π rad/s. Figure 3 represents the solution obtained numerically for N x = N y = 400, where the time-stepping formula (16) was iterated up to t ≈ 0.4 s. As time increases, a directional wave beam propagates along the X -direction.
Comparisons between the numerical solution and the analytical solution of Eq. (44) in Figs. 3a-3b show that the numerical method produces consistent results. Due to diffraction, the long-time velocity amplitude decreases as X −1/2 with the distance of propagation X . The slight amplitude mismatch between numerical and analytical computations is due to the numerical diffusion of the MUSCL scheme. The small phase mismatch is caused by the explicit integration of the source. All these numerical artifacts vanish as the mesh is refined.

Nonlinear initial-value problem
This configuration is spatially one-dimensional, and the source term S n i , j is zero. The initial data p 0 i , j = p • (x i , y j ) with wavelength λ = 0.4 m reads in m/s. The time-stepping formula (16) is applied up to t ≈ 0.2 s, with a Courant number Co = 0.95. This value, larger than 0.5, does not induce numerical instability due to the one-dimensional nature of the problem. The initial data (27) consists of a rectangular bump, which is invariant along Y (see analytical developments in C).
Similarly to Burgers' equation with rectangular data [31], the solution is made of a rarefaction and a shock, which interact after a certain amount of time. This is illustrated in Fig. 4, where the solution obtained with N x = 200 is displayed at three different times. The number of cells in the Y -direction includes twice the schemes stencil, i.e. N y = 5. The position of the discontinuity in Fig. 4c is obtained quasi-analytically, by numerical integration of the Rankine-Hugoniot condition. We note that the method captures both the rarefaction and the shock wave, and that the latter is well-located.

Harmonic generation in Gaussian beams
The main goal of this section is to investigate harmonic generation numerically by solving boundary-value problems. This configuration is closely related to other studies in the literature, and it provides a natural way to control velocity amplitudes. Similarly to Destrade et al. [8], we consider two-dimensional motions with Gaussian sound beams.
The computational domain is reduced to X = (X , Y ) in [0, 0.6] × [−0.6, 0.6]/ 2. Therefore, the mesh size is now deduced from N x = 0.6/∆x and N y = 1.2/(∆y 2). The initial data p 0 i , j is zero, and the boundary data is specified at the domain's left boundary X = 0. Numerically, the boundary condition is imposed by implementing an incoming wave condition (Sec. 7.3.2 of [15]). The domain's top, bottom and right boundaries have the same absorbing properties as in the previous section. Receivers are placed every 0.05 m along the line Y = 0 to record the signal in time.

Linear case
Here the nonlinearity coefficients β 2 , β 3 are set to zero. The boundary data is that of a pure anti-plane shear beam, with only non-zero component v 2 (0, Y , t ) = ah(Y )s(t ), where the spatial evolution is a Gaussian function h(Y ) = exp(−( ω c Y ) 2 ). In practice, the function h(Y ) is truncated at the distance 3 c ω from the origin, where the Gaussian has sufficiently vanished. We take a causal sinusoidal signal s(t ) = sin(ωt ) with amplitude a = 0.7 m/s and angular frequency ω = 20π rad/s. The solution is obtained numerically for N x = 250 and N y = 500, where the time-stepping formula (16) was iterated up to t ≈ 0.6 s. With the present grid, we have 63 points per wavelength at the fundamental frequency, and 10 points per wavelength at the sixth harmonic frequency. A snapshot of the final numerical solution is shown in Figure 5a, and a video of the simulation is provided in the supplementary material. To estimate whether a given harmonic amplitude is significant or not, we measure the harmonic amplitudes along the beam axis. Figure 5b displays the evolution of the harmonic amplitudes with the propagation distance. The sine and cosine Fourier coefficients a n , b n at the angular frequency ω were computed over the last period of signal by numerical integration (trapezoidal rule). Then, the harmonic amplitudes H n = (a n , b n ) 2 were divided by the theoretical harmonic amplitude H 1 of the first harmonic, see analytical solution (50) in C. At each receiver, a small amount of undesired harmonics -mainly odd ones -is spuriously generated by the numerical procedure. Therefore, in the nonlinear cases below, only harmonic amplitudes larger than those in Fig. 5b will be considered to be physically significant.

Cubic nonlinearity only
The nonlinearity coefficient β 2 is set to zero, while β 3 is taken from Table 1. The other parameters are the same as in the linear case. Without the quadratic nonlinearity coefficient β 2 , the system (11) governing displacement components decouples and v 3 remains equal to zero. Figure 6 illustrates the generation of odd harmonics with increasing propagation distances. In the farfield, numerical results show that harmonic generation slows down as waves propagate, due to the combined effects of nonlinearity and wave diffraction (diminution of wave amplitudes). In the absence of diffraction (α = 0), a shock would have formed at the distance X s ≈ 0.08 m [5] (see also C). By making amplitudes decrease as X −1/2 , diffraction prevents wave breaking, and the solution keeps smooth during the simulation.

Quadratic and cubic nonlinearity
Both nonlinearity coefficients β 2 , β 3 are taken from Table 1. In the present configuration, the full system is solicited, and the results are included in Fig. 6. We note that the picture is very similar to the purely cubic case, up to the fact that the second harmonic is slightly more present. This observation confirms that such shear waves produce mainly odd harmonics, and it is consistent with the analysis of Destrade et al. [8] when the initial data is a pure anti-plane shear beam.

Conclusion
We introduced a numerical method that solves the coupled nonlinear partial differential equations governing 2D directional shear waves in elastic solids. We proposed a change of variables leading to a hyperbolic system of firstorder partial differential equations, where the time variable is the physical time t . The resulting system is in quasilinear form, which suggests that specific numerical methods can be implemented. We showed that a MUSCL-Osher path-conservative Godunov-type method provides a second-order shock-capturing algorithm. The slope-limiting procedure prevents spurious oscillations from being produced around discontinuities. Numerical examples illustrate how the algorithm can be used to study the propagation of nonlinear shear-wave beams. The method has great potential. It could be used to investigate how various polarizations of v 2 and v 3 lead to different harmonic generation features. As inferred by Destrade et al. [8], the second harmonic may be generated more substantially if the velocity field v 3 was not zero at the boundary. According to their calculations, by enhancing the coupling between both components of the velocity field through the quadratic term, the second harmonic may reach magnitudes similar to the fifth harmonic.
The method could also be extended to other KZK-type equations which include dissipation [6,14]. It could be extended to nonlinear viscoelastic and anisotropic materials [32][33][34], and to slightly compressible materials. Moreover, the modelling of directional wave beams in pre-stressed solids -a.k.a. materials submitted to conditioning [35] -is an open problem [36].
The implementation of such a method in the full three-dimensional case remains a challenge. The development of an efficient code using parallelization, higher-order methods and adaptive mesh refinement would be useful for applications [18,19].  (30), where the index i increases -respectively, the index j decreasesfrom the left (i = 0) to the right ( j = 0).
[40] P. D. Lax, Hyperbolic systems of conservation laws and the mathematical theory of shock waves, in: CBMS-NSF

A Using the principal invariants
Let us express the principal invariants I C , II C , III C of the right Cauchy-Green tensor C = F F in terms of the invariants I k = tr E k : or inversely, The incompressibility constraint (2) imposes III C ≡ 1, which implies that one invariant I k depends on the two others. Substituting the expressions (29) in the strain energy function (4) leads to a fourth-order Rivlin series [8] which coefficients are given in Table 2. Using the chain rule along with the tensor derivatives of the principal invariants ∂I C /∂C = I and ∂II C /∂C = I C I −C , the constitutive law is written as with the coefficients The above constitutive law is rewritten in terms of the Cauchy stress tensor σ = F SF in Eq. (6).

B System matrices
The matrices of the first-order system in retarded timet are specified below: The source term has components s = 1 2 c/µ (0, 0, f 2 , 0, 0, f 3 ) . The matrices A(p) = [A i j ], B(p) = [B i j ] and the vector S(p) = [S i ] of the quasi-linear system of balance laws (15) in real time t are deduced from the above arrays. Non-zero coefficients are detailed below: Initial value problems The initial value problem p(X , 0) = p • (X ) can be solved analytically in terms of the Riemann invariantsw 1 andw 2 of (51), which satisfy the scalar transport equationsw i ,t +w i ,1 /σ i = 0. Solutions to smooth initial value problems can be expressed in implicit form up to the breaking time by applying the method of characteristics. In particular, if the invariantw 1 is constant in space and time, then both velocity components of (51) are transported at constant speed 1/σ 2 -that is to say, waves propagate linearly at the same speed. Similarly, if we assume thatw 2 is a constant, then the velocity components of (51) are advected non-linearly at the speed 1/σ 1 withw 1 = (1 +w 2 2 )v 2 2 (or equivalently,w 1 = (1 +w −2 2 )v 2 3 ). The system (51) decouples, and it can be rewritten in conservation form as and where artanh is the inverse hyperbolic tangent function. Thus, we note that the linear advection equation is recovered at small amplitudes. Solving the Riemann problem of (52) for shock and rarefaction waves requires particular care, since the artanh function is neither convex nor concave (see e.g. Ref. [41] and references therein). To avoid complications, the example considered in this document involves data located on the same side of the inflection point.

Boundary value problems
The same method applies for the boundary value problem p(0, Y , t ) = p • (t ). The Riemann invariantsw 1 andw 2 of (51) satisfyw i ,1 + σ iwi ,t = 0. Solutions to smooth boundary value problems can be expressed in implicit form up to the shock distance. If the only nonzero component at the boundary is v 2 (0, Y , t ) = a sin(ωt ), then the shock distance reads X s = c 3 /(β 3 ωa 2 ) [5,40].