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.

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).

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 x of the current configuration. The deformation gradient is the second-order tensor defined as F=x/X. Introducing the displacement field u=xX, and the identity tensor I=[δij], whose components are represented by the Kronecker delta function, the relationship F=I+u is obtained, where denotes the gradient operator with respect to the material coordinates X=(x,y,z).

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

J=detF1.
(1)

Thus, the mass density, ρ, is constant in time. It follows also that J̇=JtrL0, where the dot denotes the material time derivative, /t. The notation trL stands for the trace Lii of the Eulerian velocity gradient tensor L=[Lij], where summation over repeated indices was assumed (Einstein notation). The components Lij=ḞikFkj1 of L=ḞF1 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(FTFI), where FT=[Fji] denotes the transpose of F. Indeed, the finite strain tensor,

E=12(u+Tu+Tuu),
(2)

is often a preferred choice in physical acoustics; see, for instance, Zabolotskaya et al. (2004). Also, we introduce its rate, Ė=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, ρv̇=P, where v=ẋ 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 [P]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.

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/ρ 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,

I1=tr(E),I2=tr(E2),I3=tr(E3).
(3)

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,

D=tr((SSe)Ė)=tr(SvĖ)0,
(4)

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

Se=pC1+WeE
(5)

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 C1=F1FT). 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{Ė}, which is a consistent choice to enforce frame-indifference (Ball, 2002; Destrade et al., 2013). We introduce a dissipation potential, Wv(E,Ė), such that

Sv=WvĖ
(6)

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

I4=tr(Ė),I5=tr(Ė2),I6=tr(Ė3),I7=tr(ĖE),I8=tr(ĖE2),I9=tr(Ė2E),I10=tr(Ė2E2).
(7)

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

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

I1=I243I3I12+2I1I223I13,
(8)

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 Ik; see the appendix in Destrade et al. (2010a). Using the differential version of the incompressibility constraint, the invariants (3)–(7) of E, Ė satisfy the particular relationship

12I4=I7+2I1I72I8(I1I2+I12)I4,
(9)

deduced from the identity trD=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(Ė) is still linear with respect to the components of the strain-rate tensor, Ė. However, Eq. (9) shows that I4 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.

In weakly nonlinear elasticity, the strain energy density function is sought in the form of a polynomial of the invariants, Ik, 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

We=μI2+13AI3+DI22,
(10)

where μ0 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, Ė, and a zeroth-order polynomial of E. This assumption amounts to selecting Wv of second order in (E,Ė) and ignoring the terms proportional to Ė that produce elastic stresses. Due to the relationships of Eqs. (8) and (9), we, therefore, keep

Wv=ηI5,
(11)

where η0 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 ζ=23η 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(Ik,) yields, respectively, the following elastic [Eq. (5)] and viscous stress contributions [Eq. (6)]:

Se=pC1+2(μ+2DI2)E+AE2,
(12a)
Sv=2ηĖ.
(12b)

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 Ė (Maugin, 1999). A sufficient condition for the restriction D0 to be always satisfied is that the viscosity, η, is non-negative.

In complement, we consider the alternative dissipation potential Wv=ηtr(D2), which is related to the mechanics of Newtonian fluids. This quantity is equal to ηtr(Π̇Ė), where we have introduced the rate of a Piola-type strain tensor, Π=C1E (see the  Appendix). The viscous stress (6) takes the Newtonian form

Sv=2ηΠ̇,
(13)

which is more commonly written as Tv=2ηD in terms of the Cauchy stress tensor, T=FSFT. 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 u, which are assumed to be arbitrarily small. The same linearisation procedure is also applied to P=Pe+Pv with P=FS. Using the previous constitutive laws, we end up with the linearised expressions PepI+2με and Pv2ηε̇ of the elastic and viscous stresses, respectively, where ε=12(u+Tu) is the infinitesimal strain tensor.

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

F=[δij+ui,j]=[10γ010001],
(14)

where γ=u/z is the shear strain. The velocity field takes the form v=[v,0,0]T, where v=u/t is the shear velocity.

In the equation of motion, ρv̇=P, the relevant first Piola–Kirchhoff stress component, P13, is deduced from the elastic part, P13e=μγ+Γγ3, where Γ=μ+A/2+D is a parameter of nonlinearity and from the viscous part. Keeping only terms up to order γ3, the viscous stress in Eq. (12b) yields P13v=η(1+2γ2)γ̇, whereas the Newtonian stress (13) yields P13v=ηγ̇. Thus, to combine both models into one expression, we will use P13v=η(1+δγ2)γ̇, where δ{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

1c22ut2=2uz2+23βz(uz)3+τz[(1+δ(uz)2)2uzt],
(15)

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

c=μρ,β=32Γμ,τ=ημ.
(16)

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

1c22γt2=2γz2+23β2z2γ3+τ2z2[(1+δγ2)γt].
(17)

Here, we have obtained the same wave equations as those that were derived in Destrade et al. (2013) for the particular bulk viscosity ζ=23η. 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 c2utt=uzz+τuzzt. Dispersion analysis shows that the complex velocity of a linear harmonic wave equals c1+iωτ, 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 c2 m/s in gels (Jacob et al., 2007), whereas β10 and τ0.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.

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̃=ϵ2z,t̃=tz/c,u=ϵũ}, where ϵ is a small parameter and ũ=ũ(z̃,t̃). Furthermore, we assume that τ is of order ϵ2. 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

ϵ3c2ũz̃t̃=ϵ3βc2(ũt̃)22ũt̃2+ϵτ23ũt̃3.
(18)

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

cvz+(1βv2/c2)vt=τ22vt2,
(19)

for the velocity v=u/t.

Up to the choice of time variable used here (i.e., the physical time, t, instead of the retarded time, t̃), 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), δτ(γ2γ̇)/z, 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, τγ̇/z. 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̃=ϵ2t,z̃=zct,u=ϵũ}, where ϵ is a small parameter. Proceeding in a fashion similar to that above, we end up with the nonlinear transport equation,

γt+c(1+βγ2)γz=τc222γz2,
(20)

where γ=u/z 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.

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, γt+cγz=(τ/2)c2γzz. It follows that the complex velocity of a harmonic wave reads (c/2)(1+1+2iωτ). 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(ωτ)—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.

FIG. 1.

(Color online) Infinitesimal strain limit. The evolution of the phase velocity with the dimensionless frequency ωτ is depicted.

FIG. 1.

(Color online) Infinitesimal strain limit. The evolution of the phase velocity with the dimensionless frequency ωτ is depicted.

Close modal

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,

g(z¯,t¯)=23βγ(z,t),t¯=t/τ,z¯=z/(cτ),
(21)

in Eq. (17), such that

2gt¯2=2gz¯2+2z¯2g3+2z¯2[(1+3δ2βg2)gt¯].
(22)

Next, we seek travelling wave solutions of the form g=ν21G(ξ), where ξ=(ν21)(t¯z¯/ν) involves the dimensionless wave velocity ν1. 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,

G=G3+(1+αG2)ddξG,
(23)

where α=32δ(ν21)/β 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=±1 by following a smooth transition that depends on the parameter α. Solutions read (Destrade et al., 2013)

ξ=ln[12G(43(1G2))(1+α)/2]
(24)

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

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

gt¯+z¯(g+12g3)=122gz¯2.
(25)

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=ν21G(χ) with χ=(ν21)(ϑt¯z¯) in Eq. (25), where ϑ=1+12(ν21) is the new dimensionless velocity (Fig. 2). Thus, we arrive at the differential equation

G=G3+ddχG,
(26)

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

G=11+3e2χ,
(27)

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

FIG. 2.

Scaled velocity, ϑ=1+12(ν21), for the “slow time” reduced model in terms of the scaled velocity ν for the full wave equation.

FIG. 2.

Scaled velocity, ϑ=1+12(ν21), for the “slow time” reduced model in terms of the scaled velocity ν for the full wave equation.

Close modal

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¯,t¯)=2β/3v(z,t)/c in Eq. (19) to get

rz¯+t¯(r12r3)=122rt¯2.
(28)

Next, we introduce r=ν21V(ψ), where ψ=(ν21)(t¯z¯/κ) involves the dimensionless velocity, κ, which is defined by the relationship κ1=112(ν21). This way, we obtain the same differential equation, V=V3+(dV/dψ), 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).

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=ϑ/ν1 on the scaled velocity as a function of ν. To ensure that the latter remains less than 5% (respectively, 1%), we obtain the requirement ν1.3 (respectively, ν1.1), which is marked by dotted lines in Fig. 2.

Now, let us observe that for a unit kink covering the range 0G1, the corresponding shear strains satisfy

0γβγmaxβ,γmax=α/δ,
(29)

where α=32δ(ν21)/β 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 γ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 γβ1.0 (respectively, γβ0.56). Note that the parameter of nonlinearity can take such values as β10 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 ν1 (or ϑ1), in terms of the maximum strain amplitude for the full wave equation with δ = 2 and its one-way approximation. According to the expression of γmax, we have the relationship ϑ1=13(γmaxβ)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.

FIG. 3.

(Color online) For the full wave equation with δ = 2 (solid line) and slow time reduced model (dashed line), we represent the evolution of the relative velocity, ν1 (respectively, ϑ1), of travelling waves in terms of the parameter γmaxβ, which involves the strain amplitude and coefficient of nonlinearity. The axes have a logarithmic scale.

FIG. 3.

(Color online) For the full wave equation with δ = 2 (solid line) and slow time reduced model (dashed line), we represent the evolution of the relative velocity, ν1 (respectively, ϑ1), of travelling waves in terms of the parameter γmaxβ, which involves the strain amplitude and coefficient of nonlinearity. The axes have a logarithmic scale.

Close modal

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).

FIG. 4.

(Color online) Steady waveforms are deduced from Eqs. (24)–(27) for increasing values of the parameter 0α3 (arrow), with the evolution of the scaled shear strain, G, in terms of the related dimensionless coordinate, ξ or χ.

FIG. 4.

(Color online) Steady waveforms are deduced from Eqs. (24)–(27) for increasing values of the parameter 0α3 (arrow), with the evolution of the scaled shear strain, G, in terms of the related dimensionless coordinate, ξ or χ.

Close modal

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=RQ(γ), where Q(γ)=c0γ1+2βg2dg, depends explicitly on the strain γ. The scalar R is an arbitrary constant, for instance, R0 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

γt+c1+2βγ2γz=0,
(30)

where we have used the equality of mixed partials, v/z=γ/t. Obviously, the lossless one-way wave equation (20) with τ = 0 is an approximation of Eq. (30) for 2βγ21.

Let us analyse this requirement in a more quantitative manner. To ensure that the relative error on the advection velocity E=(1+a)/1+2a1 for a=βγ2 remains less than 5% (respectively, 1%), we obtain the requirement a0.44 (respectively, a0.16). Application of the square root leads to the restriction γβ0.66 (respectively, γβ0.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=Q(γ) with respect to time produces

cvz+(1+2βγ2)1/2vt=0,
(31)

where the strain γ=Q1(v) can be expressed formally as a function of the velocity despite no analytical expression of the inverse function Q1 of Q is known in the present case. If |v| is small, then we can use the approximation γv/c of the strain, which follows from the asymptotic equivalence of Qcγ at small strains. Next, the (·)1/2-factor in Eq. (31) can be approximated by the polynomial expression 1βγ2, as long as 2βγ21. 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βv2/c21 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 βv2/c2 and of βγ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 γmax=K/2β, where K is a parameter. In other words, γ(z,0) equals γmax(12|z|/L) for |z|L/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=2z/L and t=2ct/L, the obtained strain solution reads

γ(z,t)γmax={1z±t1+K2((1z)2t2)1K2t2,or(1z)2t22(1z)ifK2t2=1,
(32)

for |zt|1, and is zero elsewhere, where the expressions with the top sign (respectively, the bottom sign) are defined over the domain zt1+K20 (respectively, 0). This solution is only valid for small times tt, beyond which it becomes multivalued. The critical time t=(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

γ(z,t)γmax=11+2K2t(zt1)±K2t,
(33)

for |zt|1, and is zero elsewhere, where the expression with the top sign (respectively, the bottom sign) is defined over the domain, zt(1+12K2)0 (respectively, 0). This solution is valid for small times, t2/K2, at which a shock wave is formed.

Figure 5 displays the solutions (32) and (33) at increasing times t, 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=γmax2β 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|zt±1|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.

FIG. 5.

(Color online) Inviscid case. The evolution of the waveforms from initial triangular strain with increasing dimensionless times t is depicted.

FIG. 5.

(Color online) Inviscid case. The evolution of the waveforms from initial triangular strain with increasing dimensionless times t is depicted.

Close modal

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,

[[v]]=σ[[γ]],[[P13e]]=ρσ[[v]],
(34)

where [[·]]=(·)+(·) denotes the jump of a quantity across the discontinuity. Thus, ρσ2[[γ]]=[[P13e]]. For a right-going discontinuity σ0 propagating into an undeformed domain, (γ)+=0, we find

σ=c1+23β(γ)2,
(35)

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, tγ=zf(γ), for a suitable definition of the flux f. In this case, the Rankine–Hugoniot condition for a discontinuity travelling at speed σ=[[f]]/[[γ]] yields

σ=c(1+13β(γ)2).
(36)

Obviously, Eq. (36) is an approximation of Eq. (35) for 23β(γ)21. By setting a=13β(γ)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.

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 τ̂24γ/z2t2 in the right-hand side of the wave equation (17); see Destrade and Saccomandi (2008) for related theories. With τ̂2 of order ϵ2, this assumption leads to the addition of 12τ̂23v/t3 on the right-hand side of Eq. (19) and to the addition of 12c3τ̂23γ/z3 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).

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.

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

C3IC2+IICIIII=0,
(A1)

where I,II,andIII 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 C1Ė on the right side and substitution of C=I+2E lead to

C1Ė=(III+1)Ė+(42I)Ė+4E2Ė.
(A2)

Computation of the trace entails Eq. (9), where we have used the incompressibility property, trD=tr(C1Ė)=0, the definition of the invariants [Eqs. (3)–(7)], and the relationship between I,II, and the invariants Ik 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 Π=C1E, which satisfies

Π̇=C1ĖC1=F1DFT,Π̇Ė=(C1Ė)2=F1D2F.
(A3)

These identities lead to expression (13) of the viscous stress and dissipation potential, Wv=ηtr(Π̇Ė). Squaring Eq. (A2) and computing the trace yields a lengthy expression for Wv (hint: write the Cayley–Hamilton identity for E), which satisfies WvηI7 at leading order in E, Ė. However, the Newtonian theory (13) and viscous stress (12b) might differ at higher order.

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 /z are replaced by (1/c)/t, a relationship that is only valid for travelling waves of speed c.

1.
Antman
,
S. S.
, and
Malek-Madani
,
R.
(
1988
). “
Travelling waves in nonlinearly viscoelastic media and shock structure in elastic media
,”
Quart. Appl. Math.
46
(
1
),
77
93
.
2.
Ball
,
J. M.
(
2002
). “
Some open problems in elasticity
,” in
Geometry, Mechanics, and Dynamics
, edited by
P.
Newton
,
P.
Holmes
, and
A.
Weinstein
(
Springer
,
New York
), pp.
3
59
.
3.
Berjamin
,
H.
, and
Chockalingam
,
S.
(
2022
). “
Shear shock formation in incompressible viscoelastic solids
,”
Wave Motion
110
,
102899
.
4.
Berjamin
,
H.
,
Destrade
,
M.
, and
Parnell
,
W. J.
(
2021
). “
On the thermodynamic consistency of quasi-linear viscoelastic models for soft solids
,”
Mech. Res. Commun.
111
,
103648
.
5.
Berjamin
,
H.
,
Lombard
,
B.
,
Chiavassa
,
G.
, and
Favrie
,
N.
(
2017
). “
Analytical solution to 1D nonlinear elastodynamics with general constitutive laws
,”
Wave Motion
74
,
35
55
.
6.
Carcione
,
J. M.
(
2015
).
Wave Fields in Real Media
, 3rd ed. (
Elsevier Science
,
Amsterdam, The Netherlands
).
7.
Catheline
,
S.
,
Gennisson
,
J.-L.
,
Tanter
,
M.
, and
Fink
,
M.
(
2003
). “
Observation of shock transverse waves in elastic media
,”
Phys. Rev. Lett.
91
(
16
),
164301
.
8.
Catheline
,
S.
,
Gennisson
,
J.-L.
,
Tanter
,
M.
, and
Fink
,
M.
(
2005
). “
Erratum: Observation of shock transverse waves in elastic media [Phys. Rev. Lett. 91, 164301 (2003)]
,”
Phys. Rev. Lett.
95
(
13
),
139902
.
9.
Cormack
,
J. M.
, and
Hamilton
,
M. F.
(
2018
). “
Plane nonlinear shear waves in relaxing media
,”
J. Acoust. Soc. Am.
143
(
2
),
1035
1048
.
10.
Destrade
,
M.
,
Gilchrist
,
M. D.
, and
Murphy
,
J. G.
(
2010a
). “
Onset of nonlinearity in the elastic bending of blocks
,”
J. Appl. Mech.
77
(
6
),
061015
.
11.
Destrade
,
M.
,
Gilchrist
,
M. D.
, and
Ogden
,
R. W.
(
2010b
). “
Third- and fourth-order elasticities of biological soft tissues
,”
J. Acoust. Soc. Am.
127
(
4
),
2103
2106
.
12.
Destrade
,
M.
, and
Saccomandi
,
G.
(
2008
). “
A note about waves in dissipative and dispersive solids
,” in
Waves and Stability in Continuous Media
, edited by
N.
Manganaro
,
R.
Monaco
, and
S.
Rionero
(
World Scientific
,
Singapore
), pp.
210
217
.
13.
Destrade
,
M.
,
Saccomandi
,
G.
, and
Vianello
,
M.
(
2013
). “
Proper formulation of viscous dissipation for nonlinear waves in solids
,”
J. Acoust. Soc. Am.
133
(
3
),
1255
1259
.
14.
Hamilton
,
M. F.
, and
Blackstock
,
D. T.
, eds. (
1998
).
Nonlinear Acoustics
(
Academic Press
,
San Diego
).
15.
Holzapfel
,
G. A.
(
2000
).
Nonlinear Solid Mechanics: A Continuum Approach for Engineering
(
Wiley
,
Chichester
).
16.
Jacob
,
X.
,
Catheline
,
S.
,
Gennisson
,
J.-L.
,
Barrière
,
C.
,
Royer
,
D.
, and
Fink
,
M.
(
2007
). “
Nonlinear shear wave interaction in soft solids
,”
J. Acoust. Soc. Am.
122
(
4
),
1917
1926
.
17.
John
,
F.
(
1976
). “
Delayed singularity formation in solution of nonlinear wave equations in higher dimensions
,”
Commun. Pure Appl. Math.
29
(
6
),
649
682
.
18.
Jordan
,
P. M.
, and
Puri
,
A.
(
2005
). “
A note on traveling wave solutions for a class of nonlinear viscoelastic media
,”
Phys. Lett. A
335
(
2-3
),
150
156
.
19.
MacCamy
,
R. C.
(
1970
). “
Existence, uniqueness and stability of solutions of the equation
utt=(/x)(σ(ux)+λ(ux)utt),”
Indiana Univ. Math. J.
20
(
3
),
231
238
.
20.
Maugin
,
G. A.
(
1999
).
The Thermomechanics of Nonlinear Irreversible Behaviors
(
World Scientific
,
Singapore
).
21.
Naugolnykh
,
K.
, and
Ostrovsky
,
L.
(
1998
).
Nonlinear Wave Processes in Acoustics
(
Cambridge University Press
,
Cambridge, UK
).
22.
Nazarov
,
V. E.
,
Kiyashko
,
S. B.
, and
Radostin
,
A. V.
(
2017
). “
Stationary waves in a bimodular rod of finite radius
,”
Wave Motion
75
,
72
76
.
23.
Pioletti
,
D. P.
, and
Rakotomanana
,
L. R.
(
2000
). “
Non-linear viscoelastic laws for soft biological tissues
,”
Eur. J. Mech. A-Solids
19
(
5
),
749
759
.
24.
Pucci
,
E.
, and
Saccomandi
,
G.
(
2012
). “
On the nonlinear theory of viscoelasticity of differential type
,”
Math. Mech. Solids
17
(
6
),
624
630
.
25.
Pucci
,
E.
,
Saccomandi
,
G.
, and
Vergori
,
L.
(
2019
). “
Linearly polarized waves of finite amplitude in pre-strained elastic materials
,”
Proc. R. Soc. A
475
(
2226
),
20180891
.
26.
Radostin
,
A.
,
Nazarov
,
V.
, and
Kiyashko
,
S.
(
2013
). “
Propagation of nonlinear acoustic waves in bimodular media with linear dissipation
,”
Wave Motion
50
(
2
),
191
196
.
27.
Rozanova-Pierrat
,
A.
(
2009
). “
On the derivation of the Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation and validation of the KZK-approximation for viscous and non-viscous thermo-elastic media
,”
Commun. Math. Sci.
7
(
3
),
679
718
.
28.
Saccomandi
,
G.
, and
Vianello
,
M. S.
(
2021
). “
Shear waves in a nonlinear relaxing media: A three-dimensional perspective
,”
J. Acoust. Soc. Am.
149
(
3
),
1589
1595
.
29.
Whitham
,
G. B.
(
1999
).
Linear and Nonlinear Waves
(
Wiley
,
New York
).
30.
Wochner
,
M. S.
,
Hamilton
,
M. F.
,
Ilinskii
,
Y. A.
, and
Zabolotskaya
,
E. A.
(
2008
). “
Cubic nonlinearity in shear wave beams with different polarizations
,”
J. Acoust. Soc. Am.
123
(
5
),
2488
2495
.
31.
Zabolotskaya
,
E. A.
,
Hamilton
,
M. F.
,
Ilinskii
,
Y. A.
, and
Meegan
,
G. D.
(
2004
). “
Modeling of nonlinear shear waves in soft solids
,”
J. Acoust. Soc. Am.
116
(
5
),
2807
2813
.