Simple strain-rate viscoelasticity models of isotropic soft solid are introduced. The constitutive equations account for finite strain, incompressibility, material frame-indifference, nonlinear elasticity, and viscous dissipation. A nonlinear viscous wave equation for the shear strain is obtained exactly and corresponding one-way Burgers-type equations are derived by making standard approximations. Analysis of the travelling wave solutions shows that these partial differential equations produce distinct solutions, and deviations are exacerbated when wave amplitudes are not arbitrarily small. In the elastic limit, the one-way approximate wave equation can be linked to simple wave theory and shock wave theory, thus, allowing direct error measurements.

## I. INTRODUCTION

In nonlinear acoustics, the Burgers equation is often viewed as the simplest model equation that includes nonlinear wave propagation and diffusion effects (Whitham, 1999). This partial differential equation in space and time can be derived directly from the one-dimensional Navier–Stokes equation by dropping the pressure term or as a special case of the Westervelt equation. In addition to Burgers' equation, other one-way wave equations have been derived to describe wave propagation in fluids and solids at large amplitudes (Hamilton and Blackstock, 1998; Naugolnykh and Ostrovsky, 1998). Based on an appropriate scaling of the wave amplitude, such approximate partial differential equations describe unidirectional wave motion for slowly varying wave profiles of moderate amplitude.

One-way approximate wave equations have found applications in various areas of nonlinear acoustics. For instance, works by Radostin *et al.* (2013) and Nazarov *et al.* (2017) describe compression wave propagation in solids with bimodular elastic behaviour. Another example is the Zabolotskaya equation, which describes unidirectional plane shear wave propagation in soft solids such as gels and brain tissue (Zabolotskaya *et al.*, 2004); see also Cormack and Hamilton (2018). In these latter cases, the underlying three-dimensional constitutive theories were revisited by Destrade *et al.* (2013) as well as Saccomandi and Vianello (2021) to enforce objectivity (i.e., invariance by change of observer), leading to slight modifications of the equations of motion.

For these partial differential equations, not many analytical solutions are known. Nevertheless, it is sometimes possible to derive exact stationary wave solutions that keep an invariant wave profile throughout the motion, which occurs at a suitable constant speed. Those permanent waveforms result from the interaction between nonlinearity and dispersion (here of dissipative nature), a common feature that they share with solitary waves.

One might wonder whether it is preferable to seek closed-form travelling wave solutions by using directly the full equations of motion or their one-way approximation. As a matter of fact, both approaches have been considered separately in the above literature. The present study aims to provide evidence to advocate for a derivation of travelling waves based on the complete equations of motion, thus supporting a remark by Jordan and Puri (2005) in relation with the study by Catheline *et al.* (2003)—this remark led to the publication of an erratum that briefly discusses the validity of a particular one-way wave equation (Catheline *et al.*, 2005).^{1}

For this purpose, we consider the case of shear wave propagation in soft viscoelastic solids of strain-rate type. We derive the simplest isotropic constitutive theories that account for finite strain, incompressibility, material frame-indifference, and viscous dissipation (Sec. II). Then, this framework is applied to simple shear deformations, also known as transverse plane waves (Sec. III), including the reduction to a one-way model described by a Burgers-type equation with cubic nonlinearity. Finally, we investigate the travelling wave solutions deduced from the full equations of motion as well as those from the reduced wave equations (Sec. IV). The results show non-negligible discrepancies introduced by the reduction to unidirectional motion as soon as wave amplitudes are no longer infinitesimal. These comparisons are reconsidered in the lossless elastic limit where connections between the one-way model and other theories are established (Sec. V).

## II. STRAIN-RATE MODEL

### A. Basic equations

In what follows, we present the basic equations of Lagrangian dynamics for incompressible solids (Holzapfel, 2000). We consider a homogeneous and isotropic solid continuum on which no external body force is applied. Its motion in the Euclidean space is described by using an orthonormal Cartesian coordinate system, $(O,x,y,z)$. Thus, a particle initially located at some position ** X** of the reference configuration moves to a position

*of the current configuration. The deformation gradient is the second-order tensor defined as $F=\u2202x/\u2202X$. Introducing the displacement field $u=x\u2212X$, and the identity tensor $I=[\delta ij]$, whose components are represented by the Kronecker delta function, the relationship $F=I+\u2207u$ is obtained, where $\u2207$ denotes the gradient operator with respect to the material coordinates $X=(x,y,z)$.*

**x**For *incompressible solids*, volume does not change during deformation, which is expressed by the constraint

Thus, the mass density, *ρ*, is constant in time. It follows also that $J\u0307=J\u2009tr\u2009L\u22610$, where the dot denotes the material time derivative, $\u2202/\u2202t$. The notation $tr\u2009L$ stands for the trace *L _{ii}* of the Eulerian velocity gradient tensor $L=[Lij]$, where summation over repeated indices was assumed (Einstein notation). The components $Lij=F\u0307ikFkj\u22121$ of $L=F\u0307F\u22121$ are obtained by matrix multiplication.

Various strain tensors are defined as functions of $F$. Here, constitutive laws are expressed in terms of the Green–Lagrange strain tensor, $E=12(FTF\u2212I)$, where $FT=[Fji]$ denotes the transpose of $F$. Indeed, the finite strain tensor,

is often a preferred choice in physical acoustics; see, for instance, Zabolotskaya *et al.* (2004). Also, we introduce its rate, $E\u0307=FTDF$, obtained by differentiation with respect to time, where $D=12(L+LT)$ is the strain-rate tensor. We note that ** D** is trace-free due to incompressibility [Eq. (1)].

The motion is governed by the conservation of linear momentum equation, $\rho v\u0307=\u2207\u22c5P$, where $v=x\u0307$ is the velocity field and *ρ* is the mass density. The equation of motion involves the Lagrangian divergence of the first Piola–Kirchhoff stress tensor, $P=FS$, where $S=ST$ is the second Piola–Kirchhoff stress tensor. Those stress tensors are specified later by the provision of a constitutive law.

The present definitions are consistent with notations and conventions used in the monograph by Holzapfel (2000). In particular, the divergence of the tensor $P$ reads $[\u2207\u22c5P]i=Pij,j$ componentwise, where indices after the comma denote spatial differentiation. In some other texts, a transposed definition of the divergence is used. Then, the equation of motion involves the material divergence of the nominal stress tensor, $PT$, instead of $P$.

### B. Generalities

In this study, we consider isotropic deformable solids whose constitutive behaviour is described by the state variables $S={s,E}$, where *s* is the specific entropy. The set of variables, $S$, is coherent with the postulate of frame-indifference of the internal energy (Holzapfel, 2000). In fact, a change of observer specified by a superimposed rigid-body motion leaves $S$ invariant, as well as the internal energy, *U*. Note that the internal energy does not depend on rates of strain.

The internal energy per unit volume, *U*, is a function of state to be specified. The thermodynamic temperature is defined as the conjugate variable of *s* in the partial Legendre transform of $U/\rho $ with respect to *s* (Berjamin *et al.*, 2021). However, the explicit dependence of *U* with respect to *s* is usually omitted in the definition of a strain energy density function, $We$, such that $U=We(E)$. The strain energy, $We$, is regarded as a scalar-valued *isotropic function* of its arguments. Thus, its dependence with respect to ** E** can be reduced to a dependence with respect to three scalar invariants,

They can be used directly or other physically meaningful scalar quantities might be defined from them.

The first and second principles of thermodynamics yield the Clausius–Duhem inequality,

where $D$ is the dissipation, $S=Se+Sv$ is the total second Piola–Kirchhoff stress,

denotes the elastic part, and $Sv$ is a viscous contribution to be specified subsequently. The scalar *p* is an arbitrary Lagrange multiplier for the incompressibility constraint [Eq. (1)]; see Sec. 6.3 of Holzapfel (2000), and $C=I+2E$ is the right Cauchy–Green strain tensor, $FTF$ (i.e., its inverse is given by $C\u22121=F\u22121F\u2212T$). Therefore, according to Eq. (4), no dissipation occurs in the elastic case $S=Se$, where the viscous stress tensor $Sv$ is equal to zero.

According to the dissipation inequality (4), the viscous stress $Sv$ is a function of state and evolution variables, e.g., the set $S\u222a{E\u0307}$, which is a consistent choice to enforce frame-indifference (Ball, 2002; Destrade *et al.*, 2013). We introduce a *dissipation potential*, $Wv(E,E\u0307)$, such that

defines the viscous stress (Maugin, 1999). In general, the dissipation potential is described by additional invariants (Pioletti and Rakotomanana, 2000),

In the present study, we consider Newtonian-type viscosity models whose dissipation potential is as simple as possible.

### C. Consequences of incompressibility

First, let us investigate the consequences of the incompressibility constraint [Eq. (1)]. As noted in Jacob *et al.* (2007), the invariants (3) of ** E** are linked through

by virtue of incompressibility. This identity follows from the expression of the principal invariants of the unimodular tensor $C=I+2E$ in terms of the invariants *I _{k}*; see the appendix in Destrade

*et al.*(2010a). Using the differential version of the incompressibility constraint, the invariants (3)–(7) of

**, $E\u0307$ satisfy the particular relationship**

*E*deduced from the identity $tr\u2009D=0$; see the Appendix.

The relationship (8) means that the invariant $I1=tr(E)$ is no longer linear with respect to the components of the strain tensor $E$; instead, Eq. (8) shows that it has terms of polynomial order two and three with respect to the strain. Furthermore, due to the relationship (9), the invariant $I4=tr(E\u0307)$ is still linear with respect to the components of the strain-rate tensor, $E\u0307$. However, Eq. (9) shows that *I*_{4} is no longer invariant on the strain tensor $E$; instead, it has terms of polynomial order one, two, and three with respect to the Green–Lagrange strain.

### D. Constitutive assumptions

In weakly nonlinear elasticity, the strain energy density function is sought in the form of a polynomial of the invariants, *I _{k}*, with constant coefficients. Similar to Zabolotskaya

*et al.*(2004), we assume that the internal energy,

*U*, has a fourth-order polynomial expression with respect to the components of the strain tensor, $E$, of the form

where $\mu \u22650$ is the shear modulus (in Pa), and the coefficients *A* and *D* are higher-order elastic constants. Note in passing that this approach can be extended to the modelling of transversely isotropic elastic solids; see discussions in Destrade *et al.* (2010b).

Now, let us propose an expression for the dissipation potential. To end up with viscosity models similar to those by Destrade *et al.* (2013), we assume that the dissipation potential is a second-order polynomial expansion of the strain-rate tensor, $E\u0307$, and a zeroth-order polynomial of $E$. This assumption amounts to selecting $Wv$ of second order in $(E,E\u0307)$ and ignoring the terms proportional to $E\u0307$ that produce elastic stresses. Due to the relationships of Eqs. (8) and (9), we, therefore, keep

where $\eta \u22650$ is the shear viscosity (in Pa s). In Eq. (11), the absence of bulk viscosity, *ζ*, is due to the assumption on polynomial orders for the viscous part and the incompressibility property (9). Setting the bulk viscosity $\zeta =23\eta $ in Destrade *et al.* (2013) yields the same expressions as above.

Computation of the tensor derivatives of the potentials in Eqs. (10) and (11) by means of the chain rule for $W\u2022(Ik,\u2026)$ yields, respectively, the following elastic [Eq. (5)] and viscous stress contributions [Eq. (6)]:

Thermodynamic consistency [Eq. (4)] is ensured provided that the dissipation $D=2Wv$ is non-negative. In fact, the present dissipation potential, $Wv$, is a homogeneous function of degree two with respect to $E\u0307$ (Maugin, 1999). A sufficient condition for the restriction $D\u22650$ to be always satisfied is that the viscosity, *η*, is non-negative.

In complement, we consider the alternative dissipation potential $Wv=\eta \u2009tr(D2)$, which is related to the mechanics of Newtonian fluids. This quantity is equal to $\eta \u2009tr(\Pi \u0307E\u0307)$, where we have introduced the rate of a Piola-type strain tensor, $\Pi =C\u22121E$ (see the Appendix). The viscous stress (6) takes the Newtonian form

which is more commonly written as $Tv=2\eta D$ in terms of the Cauchy stress tensor, $T\u2022=FS\u2022FT$. This expression is also perfectly satisfactory from the point of view of material frame-indifference (Destrade *et al.*, 2013). By using similar arguments as above, thermodynamic consistency [Eq. (4)] is ensured if the viscosity, *η*, is non-negative. For further details, the interested reader is referred to Antman and Malek-Madani (1988).

In the infinitesimal strain limit, the equations of motion are linearised with respect to the components of ** u** and $\u2207u$, which are assumed to be arbitrarily small. The same linearisation procedure is also applied to $P=Pe+Pv$ with $P\u2022=FS\u2022$. Using the previous constitutive laws, we end up with the linearised expressions $Pe\u2243\u2212pI+2\mu \epsilon $ and $Pv\u22432\eta \epsilon \u0307$ of the elastic and viscous stresses, respectively, where $\epsilon =12(\u2207u+\u2207Tu)$ is the infinitesimal strain tensor.

## III. PLANE SHEAR WAVES

### A. Nonlinear viscous wave equation

Similar to Destrade *et al.* (2013), we consider simple shear deformations described by the displacement field $u=[u,0,0]T$, where $u=u(z,t)$ denotes the particle displacement along the *x*-direction. Thus, the deformation gradient tensor reads

where $\gamma =\u2202u/\u2202z$ is the shear strain. The velocity field takes the form $v=[v,0,0]T$, where $v=\u2202u/\u2202t$ is the shear velocity.

In the equation of motion, $\rho v\u0307=\u2207\u22c5P$, the relevant first Piola–Kirchhoff stress component, *P*_{13}, is deduced from the elastic part, $P13e=\mu \gamma +\Gamma \gamma 3$, where $\Gamma =\mu +A/2+D$ is a parameter of nonlinearity and from the viscous part. Keeping only terms up to order $\gamma 3$, the viscous stress in Eq. (12b) yields $P13v=\eta (1+2\gamma 2)\gamma \u0307$, whereas the Newtonian stress (13) yields $P13v=\eta \gamma \u0307$. Thus, to combine both models into one expression, we will use $P13v=\eta (1+\delta \gamma 2)\gamma \u0307$, where $\delta \u2208{0,2}$ is a parameter.

Finally, on division by the shear modulus *μ*, the *x*-component of the equation of motion produces the nonlinear wave equation

which describes transverse wave propagation along the *z*-direction, where we have introduced the notations

Spatial differentiation of Eq. (15) allows us to write a similar wave equation for the strain:

Here, we have obtained the same wave equations as those that were derived in Destrade *et al.* (2013) for the particular bulk viscosity $\zeta =23\eta $. Note in passing the presence of a nonlinear viscous term when *δ* = 2.

According to the wave equation (15), the propagation of infinitesimal shear waves is governed by the linear Kelvin–Voigt wave equation $c\u22122utt=uzz+\tau uzzt$. Dispersion analysis shows that the complex velocity of a linear harmonic wave equals $c1+i\omega \tau $, where *ω* is the angular frequency and $i$ is the imaginary unit. Thus, in the low-frequency range, such a wave propagates at the shear wave speed *c*, in the absence of nonlinearity. Typically, this sound velocity equals $c\u22482$ m/s in gels (Jacob *et al.*, 2007), whereas $\beta \u224810$ and $\tau \u22480.12$ ms at a loading frequency of 100 Hz.

For smooth initial- and boundary-value problems, the existence of global smooth solutions to Eq. (15) with *δ* = 0 is established. However, it is worth pointing out that the picture is radically different in the case *δ* = 2 for which such a solution might not exist for arbitrarily large times (MacCamy, 1970). This property is particularly inconvenient for the study of long-time quasi-static solutions, thus indicating that Eq. (13) is the more appropriate form of the viscous stress. Pucci and Saccomandi (2012) illustrate the blow-up of smooth solutions in finite time for *δ* = 2 and *β* = 0, and further links to the mathematical literature are provided therein.

### B. Slow scale approximations

Similar to Zabolotskaya *et al.* (2004) and Pucci *et al.* (2019), we proceed now to a reduction of the wave equation (15) for one-way wave propagation with slowly varying profile. We present two approximations based either on a slow space variable or a slow time variable.

#### 1. Slow space

Let us follow the scaling procedure in Zabolotskaya *et al.* (2004). For this purpose, we introduce the *scaling* defined by the change of variables ${z\u0303=\u03f52z,t\u0303=t\u2212z/c,u=\u03f5u\u0303}$, where *ϵ* is a small parameter and $u\u0303=u\u0303(z\u0303,t\u0303)$. Furthermore, we assume that *τ* is of order $\u03f52$. Note that this set of assumptions corresponds to a slowly varying profile in space.

This ansatz is then substituted into the equation of motion (15). At leading (cubic) order in *ϵ*, the motion of soft viscous solids is governed by the scalar equation

Transforming back to the initial displacement *u* and physical coordinates (*z*,*t*) leads to a reduced wave equation,

for the velocity $v=\u2202u/\u2202t$.

Up to the choice of time variable used here (i.e., the physical time, *t*, instead of the retarded time, $t\u0303$), the partial differential equation (19) is identical to the cubic Burgers-type equation in Zabolotskaya *et al.* (2004). However, the extra nonlinear viscosity term of (17), $\delta \tau \u2009\u2202(\gamma 2\gamma \u0307)/\u2202z$, is not apparent in the one-way approximation [Eq. (19)]. In fact, this additional term is lost in the rescaling procedure given that it is of higher order in *ϵ* than the leading-order viscous term, $\tau \u2009\u2202\gamma \u0307/\u2202z$. In the end, setting *δ* = 2 in the wave equation (15) does not induce any modification of the transport equation (19).

#### 2. Slow time

For later comparisons, let us derive a similar Burgers-type equation governing the evolution of the strain instead of the velocity by following Pucci *et al.* (2019). To do so, we introduce the slow time scaling based on the change of variables ${t\u0303=\u03f52t,z\u0303=z\u2212ct,u=\u03f5u\u0303}$, where *ϵ* is a small parameter. Proceeding in a fashion similar to that above, we end up with the nonlinear transport equation,

where $\gamma =\u2202u/\u2202z$ is the shear strain. Here too, after keeping leading-order terms, we have transformed back to the initial physical coordinates, (*z*,*t*). Therefore, the above partial differential equation may be viewed as a one-way approximation of the wave equation (17). Their travelling wave solutions are compared in Sec. IV; the inviscid limit, *τ* = 0, is discussed in Sec. V.

## IV. TRAVELLING WAVE SOLUTIONS

Let us begin with a brief remark on wave dispersion (Carcione, 2015). According to the slow time approximation [Eq. (20)], the propagation of infinitesimal shear waves is governed by the advection-diffusion equation, $\gamma t+c\gamma z=(\tau /2)c2\gamma zz$. It follows that the complex velocity of a harmonic wave reads $(c/2)(1+1+2i\omega \tau )$. Comparison with the viscous wave equation of Sec. III A shows that the one-way approximation is accurate in the low-frequency range only, where it provides a correct estimation of the wave velocity up to order $O(\omega \tau )$—the same analysis for Eq. (19) leads to a similar conclusion. This property is illustrated in Fig. 1, where the evolution of the phase velocities divided by *c* is displayed. With this observation in mind, we focus now on the study of nonlinear solutions.

### A. Nonlinear viscous wave equation

Let us seek travelling wave solutions to the wave equation (17), i.e., specific smooth waveforms that propagate at a constant velocity with a steady profile. In a fashion similar to that in Destrade *et al.* (2013), we introduce the following rescaled dimensionless variables and coordinates,

in Eq. (17), such that

Next, we seek travelling wave solutions of the form $g=\nu 2\u22121\u2009G(\xi )$, where $\xi =(\nu 2\u22121)(t\xaf\u2212z\xaf/\nu )$ involves the dimensionless wave velocity $\nu \u22651$. Injecting this ansatz in the above partial differential equation and integrating twice with respect to *ξ* with vanishing integration constants yields a nonlinear differential equation for the strain,

where $\alpha =32\delta (\nu 2\u22121)/\beta $ is a parameter.

From Eq. (23), one observes that travelling wave solutions to the wave equation (17) should connect the equilibrium strains *G* = 0 and $G=\xb11$ by following a smooth transition that depends on the parameter *α*. Solutions read (Destrade *et al.*, 2013)

in implicit form, where we have enforced $G(0)=12$ without loss of generality. Illustrations are provided later on.

### B. Slow time approximation

In a similar way, let us now seek travelling wave solutions to the reduced wave equation (20). Thus, we first perform the substitutions (21) to get

To obtain wave solutions that correspond to the same strain values at infinity as in Sec. IV A, we introduce a slightly different scaling. Indeed, let us inject the ansatz, $g=\nu 2\u22121\u2009G(\chi )$ with $\chi =(\nu 2\u22121)(\u03d1t\xaf\u2212z\xaf)$ in Eq. (25), where $\u03d1=1+12(\nu 2\u22121)$ is the new dimensionless velocity (Fig. 2). Thus, we arrive at the differential equation

of which the strain values zero and one are steady states. Enforcing the initial value $G=12$ at *χ* = 0 gives

which does not involve any extra parameters. One observes that this expression corresponds to the case *α* = 0 in Eqs. (23) and (24).

*Remark*. One might proceed in a similar fashion with the Burgers-type equation (19) corresponding to the *slow space* approximation. Similar to Eq. (21), we perform the substitutions $r(z\xaf,t\xaf)=2\beta /3\u2009v(z,t)/c$ in Eq. (19) to get

Next, we introduce $r=\nu 2\u22121\u2009V(\psi )$, where $\psi =(\nu 2\u22121)(t\xaf\u2212z\xaf/\kappa )$ involves the dimensionless velocity, *κ*, which is defined by the relationship $\kappa \u22121=1\u221212(\nu 2\u22121)$. This way, we obtain the same differential equation, $V=V3+(dV/d\psi )$, for the dimensionless velocity, *V*, as obtained previously for the strain [Eq. (26)]. Therefore, within the scope of this study, the slow time and slow space approximations lead to related travelling wave solutions that describe the evolution of distinct kinematic variables (strain and velocity, respectively).

### C. Comparison

Let us compare the solutions (24)–(27) obtained for the full wave equation (17) and the one-way model (20). First, one observes that these travelling waves of the same amplitude do not propagate at the same speed independently on the value of the parameter *δ* (Fig. 2). Indeed, given the expression of ϑ, we can express the relative error $E=\u03d1/\nu \u22121$ on the scaled velocity as a function of *ν*. To ensure that the latter remains less than 5% (respectively, 1%), we obtain the requirement $\nu \u22641.3$ (respectively, $\nu \u22641.1$), which is marked by dotted lines in Fig. 2.

Now, let us observe that for a unit kink covering the range $0\u2264G\u22641$, the corresponding shear strains satisfy

where $\alpha =32\delta (\nu 2\u22121)/\beta $ was introduced earlier; see Eq. (21). In other words, in the case *δ* = 2, the coefficient *α* in the differential equation (23) is related to the maximum strain $\gamma max$ of travelling waves, and these bounds are valid for both wave equations at hand due to application of the same rescaling procedure [Eq. (21)]. Thus, restrictions of the wave speed *ν* can be expressed in terms of the strain. To ensure that the velocity error, $E$, remains less than 5% (respectively, 1%), we therefore require $\gamma \beta \u22641.0$ (respectively, $\gamma \beta \u22640.56$). Note that the parameter of nonlinearity can take such values as $\beta \u224810$ for gels (Jacob *et al.*, 2007). Hence, the slow scale approximation has a very restricted validity for soft viscoelastic materials with strain-dependent shear viscosity (*δ* = 2).

This property is further illustrated in Fig. 3, where we have represented the evolution of the relative velocity $\nu \u22121$ (or $\u03d1\u22121$), in terms of the maximum strain amplitude for the full wave equation with *δ* = 2 and its one-way approximation. According to the expression of $\gamma max$, we have the relationship $\u03d1\u22121=13(\gamma max\beta )2$ in the case of the one-way approximate model, which produces lines of slope two in log-log coordinates (dashed lines in Fig. 3). However, for the full wave equation with *δ* = 2, this relationship between the wave speed *ν* and strain amplitude is not satisfied. Differences between the one-way model and the full wave equation become visible at large strains.

In Fig. 4, we display the evolution of the waveforms [Eqs. (24)–(27)] in terms of the scaled coordinates, *ξ*,*χ*. In the case of the full wave equation (23), the parameter *α* takes the values ${0,1.2,3}$. It appears that the waveforms so-obtained follow a drastically different evolution when parameters are modified. In particular, the wavefront deduced from the full wave equation (solid lines) does not necessarily exhibit the same invariance and symmetry properties as the wavefront deduced from the one-way model (dashed line).

## V. INVISCID LIMIT

### A. Simple waves

In the lossless case, exact one-way wave equations can be derived by using the method of Riemann invariants; see, for instance, the introductory example by John (1976). Such particular wave solutions, called *simple waves*, keep one Riemann invariant constant. In other words, the particle velocity $v=R\u2212Q(\gamma )$, where $Q(\gamma )=c\u222b0\gamma 1+2\beta g2\u2009dg$, depends explicitly on the strain *γ*. The scalar *R* is an arbitrary constant, for instance, $R\u22610$ in some specific boundary-value problems (Berjamin and Chockalingam, 2022), which will be assumed satisfied from now on. Spatial differentiation of the velocity then produces

where we have used the equality of mixed partials, $\u2202v/\u2202z=\u2202\gamma /\u2202t$. Obviously, the lossless one-way wave equation (20) with *τ* = 0 is an approximation of Eq. (30) for $2\beta \gamma 2\u226a1$.

Let us analyse this requirement in a more quantitative manner. To ensure that the relative error on the advection velocity $E=(1+a)/1+2a\u22121$ for $a=\beta \gamma 2$ remains less than 5% (respectively, 1%), we obtain the requirement $a\u22640.44$ (respectively, $a\u22640.16$). Application of the square root leads to the restriction $\gamma \beta \u22640.66$ (respectively, $\gamma \beta \u22640.40$), which is slightly more constraining than that in the case of viscoelastic travelling waves (Sec. IV C).

Along a simple wave, computation of the partial derivative of the velocity $v=\u2212Q(\gamma )$ with respect to time produces

where the strain $\gamma =Q\u22121(\u2212v)$ can be expressed formally as a function of the velocity despite no analytical expression of the inverse function $Q\u22121$ of *Q* is known in the present case. If $|v|$ is small, then we can use the approximation $\gamma \u2243\u2212v/c$ of the strain, which follows from the asymptotic equivalence of $Q\u223cc\gamma $ at small strains. Next, the $(\xb7)\u22121/2$-factor in Eq. (31) can be approximated by the polynomial expression $1\u2212\beta \gamma 2$, as long as $2\beta \gamma 2\u226a1$. This way, we have shown that the one-way wave equation (19) is an approximation of Eq. (31), which is obtained for *R* = 0 and $2\beta v2/c2\u226a1$ in the elastic limit. This observation is consistent with the discussions in Catheline *et al.* (2005). In summary, the lossless “slow space” and “slow time” reductions [Eqs. (19) and (20)] with *τ* = 0 are approximate governing equations for simple waves with small values of $\beta v2/c2$ and of $\beta \gamma 2$, respectively.

The general solution of the scalar wave equations (20) and (30) can be obtained by using the method of characteristics. For instance, let us consider an initial-value problem and assume that the initial strain is a triangular bump of width *L* and amplitude $\gamma max=K/2\beta $, where *K* is a parameter. In other words, $\gamma (z,0)$ equals $\gamma max(1\u22122|z|/L)$ for $|z|\u2264L/2$, and it equals zero elsewhere.

Let us first solve this problem for the simple wave governed by Eq. (30). Introducing the dimensionless coordinates $z\u2323=2z/L$ and $t\u2323=2ct/L$, the obtained strain solution reads

for $|z\u2323\u2212t\u2323|\u22641$, and is zero elsewhere, where the expressions with the top sign (respectively, the bottom sign) are defined over the domain $z\u2323\u2212t\u23231+K2\u22650$ (respectively, $\u22640$). This solution is only valid for small times $t\u2323\u2264t\u2323\u22c6$, beyond which it becomes multivalued. The critical time $t\u2323\u22c6=(1+1+K2)/K2$ is the shock formation time.

For the slow time approximation [Eq. (20)] with *τ* = 0, the solution of the initial-value problem reads

for $|z\u2323\u2212t\u2323|\u22641$, and is zero elsewhere, where the expression with the top sign (respectively, the bottom sign) is defined over the domain, $z\u2323\u2212t\u2323(1+12K2)\u22650$ (respectively, $\u22640$). This solution is valid for small times, $t\u2323\u22642/K2$, at which a shock wave is formed.

Figure 5 displays the solutions (32) and (33) at increasing times $t\u2323$, where we have set *K* = 1. Comparison of the waveforms suggests that the error increases pointwise with time. This observation can be expressed in a more quantitative fashion. In fact, a Taylor series expansion of the difference Eq. (32) − Eq. (33) with respect to the strain amplitude parameter $K=\gamma max2\beta $ shows that the slow time solution (33) approximates the simple wave (32) correctly up to a term of order $O(K4)$, pointwise. The leading-order coefficient in this expansion is proportional to $t\u2323\u2009|z\u2323\u2212t\u2323\xb11|4$, which shows that the error increases with time, and as one moves away from the domain where the solution is equal to zero (cf. Fig. 5). Possibly, computation of the global error by spatial summation of the pointwise errors would be insightful.

### B. Shock waves

The process leading to simple wave solutions requires some smoothness of the waveform, an assumption that will be relaxed in this subsection. Indeed, based on the equations of motion, the speed *σ* of a discontinuity must satisfy the Rankine–Hugoniot jump conditions,

where $[[\xb7]]=(\xb7)+\u2212(\xb7)\u2212$ denotes the jump of a quantity across the discontinuity. Thus, $\rho \sigma 2[[\gamma ]]=[[P13e]]$. For a right-going discontinuity $\sigma \u22650$ propagating into an undeformed domain, $(\gamma )+=0$, we find

where Eq. (16) and the expression of $P13e$ in Sec. III were used. Application of Liu's entropy criterion shows that such a shock wave is admissible (Berjamin *et al.*, 2017).

Let us now consider the one-way approximate wave equation (20) with *τ* = 0, which can be rewritten in conservation form, $\u2202t\gamma =\u2212\u2202zf(\gamma )$, for a suitable definition of the flux *f*. In this case, the Rankine–Hugoniot condition for a discontinuity travelling at speed $\sigma =[[f]]/[[\gamma ]]$ yields

Obviously, Eq. (36) is an approximation of Eq. (35) for $23\beta (\gamma \u2212)2\u226a1$. By setting $a=13\beta (\gamma \u2212)2$, the expression of the relative wave velocity error $E$ from Sec. V A still applies. Furthermore, the study of initial-value problems with discontinuous solutions might be conducted in a similar way as was done in Sec. IV.

## VI. CONCLUSION

For a specific strain-rate viscoelasticity theory of soft solids, we have shown that one-way approximate wave propagation models can produce significantly different travelling wave solutions than the full equations of motion as soon as the wave amplitude is not infinitesimal. Similar observations are reported in the literature in relation with shear shock formation (Berjamin and Chockalingam, 2022). In the elastic limit, we have examined the validity of one-way approximations in relation with simple wave theory and shock wave theory, thus leading to dedicated criteria involving small velocity and strain amplitudes. We conclude that these approximations should be used with care given their limited accuracy, in general. Nevertheless, they might remain useful for the interpretation of experimental results in which their validity is not always severely penalised (Catheline *et al.*, 2003, 2005).

Due to the assumption of slowly varying wave profiles, Burgers-type model equations have a limited ability to account for dispersion effects, which play an important role in viscoelasticity (Fig. 1). In addition to viscous stresses, other dispersive effects might be impacted by the rescaling procedure. For instance, one might insert the dispersive term $\tau \u03022\u2009\u22024\gamma /\u2202z2\u2202t2$ in the right-hand side of the wave equation (17); see Destrade and Saccomandi (2008) for related theories. With $\tau \u03022$ of order $\u03f52$, this assumption leads to the addition of $12\tau \u03022\u2009\u22023v/\u2202t3$ on the right-hand side of Eq. (19) and to the addition of $\u221212c3\tau \u03022\u2009\u22023\gamma /\u2202z3$ on the right-hand side of Eq. (20). Consistently, dispersion analysis shows that the slow scale equations (19) and (20) with dispersion are valid at low frequency only.

Within the present modelling framework, more general constitutive theories could be considered, e.g., to account for material anisotropy (Destrade *et al.*, 2010b). Moreover, one could consider a generalised model that accounts for stress relaxation (Saccomandi and Vianello, 2021). Last, a similar discussion could be proposed for the partial differential equations governing the diffraction of unidirectional wave beams (Wochner *et al.*, 2008); see also the case of the Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation (Rozanova-Pierrat, 2009).

## ACKNOWLEDGMENTS

The author thanks Michel Destrade (Galway, Ireland) for support and the anonymous reviewers for their valuable suggestions. This project has received funding from the European Union's Horizon 2020 research and innovation programme under Grant No. TBI-WAVES—H2020-MSCA-IF-2020 Project No. 101023950.

### APPENDIX: CONSEQUENCES OF INCOMPRESSIBILITY

This appendix is devoted to the derivation of Eq. (9). We start with the Cayley–Hamilton identity for the right Cauchy–Green tensor $C=FTF$, which reads

where $I,\u2009\u2009II,\u2009\u2009and\u2009\u2009III$ are the principal invariants of ** C**. In the case of volume-preserving motions [Eq. (1)], the tensor $C$ is unimodular, i.e., we have $III=1$. Next, multiplication of Eq. (A1) by $C\u22121E\u0307$ on the right side and substitution of $C=I+2E$ lead to

Computation of the trace entails Eq. (9), where we have used the incompressibility property, $tr\u2009D=tr(C\u22121E\u0307)=0$, the definition of the invariants [Eqs. (3)–(7)], and the relationship between $I,\u2009II$, and the invariants *I _{k}* used here (Destrade

*et al.*, 2010a).

Now, let us briefly discuss the derivation of the Newtonian viscous stress [Eq. (13)]. For this purpose, we introduce the strain tensor $\Pi =C\u22121E$, which satisfies

These identities lead to expression (13) of the viscous stress and dissipation potential, $Wv=\eta \u2009tr(\Pi \u0307E\u0307)$. Squaring Eq. (A2) and computing the trace yields a lengthy expression for $Wv$ (hint: write the Cayley–Hamilton identity for ** E**), which satisfies $Wv\u2243\eta I7$ at leading order in

**, $E\u0307$. However, the Newtonian theory (13) and viscous stress (12b) might differ at higher order.**

*E*^{1}

The wave equation proposed by Catheline *et al.* (2003) and analysed by Jordan and Puri (2005) cannot be obtained rigorously from the equations of motion unless spatial derivatives $\u2202/\u2202z$ are replaced by $\u2212(1/c)\u2202/\u2202t$, a relationship that is only valid for travelling waves of speed *c*.