Proper formulation of viscous dissipation for nonlinear waves in solids

arXiv:1303.1617v1 [cond-mat.soft] 7 Mar 2013 Proper formulation of viscous dissipation for nonlinear waves in solids Michel Destradea , Giuseppe Saccomandib , Maurizio Vianelloc a School of Mathematics, Statistics and Applied Mathematics, National University of Ireland Galway, Ireland. b Dipartimento di Ingegneria Industriale, Università degli Studi di Perugia, 06125 Perugia, Italy. c Dipartimento di Matematica, Politecnico di Milano, Piazza Leonardo da Vinci 32, I-20133 Milano, Italy. Abstract In order to model nonlinear viscous dissipative motions in solids, acoustical physicists usually add terms linear in Ė, the material time derivative of the Lagrangian strain tensor E, to the elastic stress tensor σ derived from the expansion to the third- (sometimes fourth-) order of the strain energy density E = E(tr E, tr E 2 , tr E 3 ). Here, it is shown that this practice, which has been widely used in the past three decades or so, is physically wrong for at least two reasons, and that it should be corrected. One reason is that the elastic stress tensor σ is not symmetric while Ė is symmetric, so that motions for which σ + σ T 6= 0 will give rise to elastic stresses which have no viscous pendant. Another reason is that Ė is frame-invariant, while σ is not, so that an observer transformation would alter the elastic part of the total stress differently than it would alter the dissipative part, thereby violating the fundamental principle of material frame-indifference. These problems can have serious consequences for nonlinear shear wave propagation in soft solids, as seen here with an example of a kink in almost incompressible soft solids. 1 1 Introduction and main statement Nonlinear elastic wave propagation is a subject of considerable interest for many scientific and industrial applications such as geophysical exploration, soft tissue acoustics, and the dynamics of rubbers, silicones, and gels. From a theoretical point of view the mathematics and mechanics of nonlinear wave phenomena is a classical yet still active subject of research where many outstanding problems are awaiting a definitive systematic treatment. In Physical Acoustics, the expansion of the strain energy density E to include nonlinear corrections is often attributed to Landau. Hence at the “third-order” in the Lagrangian strain tensor E, we write E = µI2 + 21 λI12 + 31 AI3 + BI1 I2 + 13 CI13 , (1) where Ik = tr E k (k = 1, 2, 3) are the strain invariants, λ and µ are the (second-order) Lamé coefficients, and A, B, C are the (third-order) Landau coefficients. From an historical point of view, it is interesting to note that this notation by Landau seems to have appeared first in the Theory of Elasticity (1986), co-written with Lifshitz, as an unnumbered equation in an exercise[1]. In 1937, Landau[2] had already provided the third-order expansion of E, that time denoting the third-order elastic constants by the letters P ′ , Q′ , R′ . That same year, Murnaghan[3] also proposed a third-order expansion using different strain invariants, which is still in use today (mostly by geophysicists). However the paternity of the third-order expansion can be traced further back in time, at least to a 1925 paper by Brillouin[4], who in fact seems to have also been the first to use the letters A, B, C for third-order elastic constants. In any event, Landau and Lifshitz denote by xk the Lagrangian coordinate components, by uk the mechanical displacement components, and they derive the equations of motion in Cartesian coordinates as follows, ∂ 2 ui ∂σik = ρ0 2 , ∂xk ∂t (2) where ρ0 is the mass density in the undeformed configuration, and the stress tensor is related to the strain through σik = ∂E . ∂ (∂ui /∂xk ) 2 (3) (It is easy to recognize that this is the “first Piola-Kirchhoff stress tensor” of Classical Continuum Mechanics.) For the reader’s convenience, we recall that the full Lagrangian strain tensor E has Cartesian components   1 ∂ui ∂uk ∂uj ∂uj Eik = , (4) + + 2 ∂xk ∂xi ∂xi ∂xk and can be rewritten as Eik = 21 (Fji Fjk − δik ) . (5) Here, Fik are the Cartesian components of F , the deformation gradient, defined as ∂yi ∂ui Fik = = δik + , (6) ∂xk ∂xk where yi (xk ) are the space coordinates of the current position of material point xk , for the deformation given by yi (xk ) = xi + ui (xk ). In another part of their book (more precisely: at the end of §34 in Chapter 5) Landau and Lifchitz model viscous dissipation by adding to the elastic stress tensor σ a “viscosity stress tensor” σ ′ , with Cartesian components ′ σik = 2η(Ėik − 13 δik Ėll ) + ζ Ėll δik = (ζ − 32 η) Ėll δik + 2η Ėik , (7) where η > 0 and ζ > 0 are the shear and bulk viscosity coefficients, respectively, and the superposed dot denotes the time derivative. Note that it is completely unambiguous from the context of this part of their book that Landau and Lifshitz are speaking here of the infinitesimal theory of viscoelasticity, and that Ė in Eq. (7) is the time derivative of the infinitesimal strain tensor (i.e. the linear part of Eq. (4)), see their unnumbered equation between (34.2) and (34.3). It might seem at first glance that adding σ ′ to σ, and taking the strain components to be those of the finite Lagrangian strain (4) instead of the infinitesimal strain, would be a first, logical step towards the inclusion of nonlinear dissipative effects. However there are at least two problems with this seemingly anodyne approach. One problem is that σ ′ is symmetric while σ is not. Then the resultant total stress tensor σ +σ ′ has a visco-elastic symmetric part: (σ +σ T )/2+σ ′ , but a purely elastic anti-symmetric part: (σ−σ T )/2, where the superscript T denotes the transpose. It is thus impossible to model antisymmetric viscous stress effects in motions and boundary conditions with this formulation. 3 The other problem is that σ and σ ′ are made objective, or frameindifferent, in two different and irreconcilable ways. Indeed, it is well known that there are transformation rules to follow in order to ensure that the directions associated with a tensor are unaltered by an observer transformation. To summarize, two observers, one associated with a frame with position x and time t and the other associated with a frame with position x∗ and time t∗ are said to be equivalent when they are connected by (see e.g. Chadwick[5]), x∗ = c(t) + Q(t)x, t∗ = t − a, (8) where Q is a proper orthogonal tensor, c is a vector, and a is a constant scalar. The transformation rule for the elastic part of the stress tensor (3) is that σ transforms into σ given by σ = Qσ, (9) in order to ensure frame-indifference, whilst for the viscous part (7), it is σ′ = σ′, (10) because E – and hence Ė – is observer-invariant (this is not so obvious from Eq. (4), but becomes more so in view of Eq. (5), as shown in Ref. [5]). Clearly, it is not possible for the composite stress tensor σ + σ ′ to comply with both requirements simultaneously and appear the same to two equivalent observers. This state of affair violates the fundamental principle of objectivity, a cornerstone of rational mechanics, e.g. see Gurtin[6]. Now, there are two ways to reconcile σ and σ ′ in order to construct a stress tensor coherent with respect to symmetry and objectivity. One way is to make σ behave like σ ′ , by replacing it with F −1 σ, the (symmetric) second Piola-Kirchhoff stress tensor; the other way is to make σ ′ behave like σ, by replacing it with σ̄ ′ = F σ ′ . Evidently, both courses of action are eventually equivalent. In weakly nonlinear elasticity theory, governing equations have been around for σ longer than for σ ′ and we henceforth take advantage of these and concentrate on the consequences of modifying σ ′ rather than σ. We thus propose the following form for the viscous stress, σ̄ ′ = F σ ′ , with components ′ ′ σ̄ik = Fij σjk
= ζ − 23 η Fik Ėll + 2ηFij Ėjk    
∂ui ∂ui 2 = ζ − 3 η δik + Ėll + 2η δij + Ėjk , ∂xk ∂xj 4 (11) and investigate the consequences of using σ̄ ′ instead of σ ′ in the equations of motion. In passing, it is interesting to show the form taken by the constitutive equation for the Cauchy stress tensor τ̄ ′ , related to σ̄ ′ by ′ j τ̄ik′ = σ̄is Fks , (12) where j = det F . For this we recall that[5] Ėik = Fli dlj Fjk , (13) where dlj is the Cartesian component of d, the Eulerian stretching tensor,   1 ∂vl ∂vj dlj = . (14) + 2 ∂yj ∂yl and that the (symmetric) left Cauchy-Green strain tensor b is defined as bjl = Fji Fli , (15) so that, in view of Eq. (13), the expression for Ėii can be written as Ėii = Fli dlj Fji = bjl dlj = tr(bd). (16) Thus, Eq. (11) yields
′ σ̄ik = ζ − 23 η Fik Ėll + 2ηFij Ėjk
= ζ − 32 η Fik bjl dlj + 2ηFij Flj dls Fsk
= ζ − 23 η Fik bjl dlj + 2ηbil dls Fsk . (17) From Eqs. (12), (13) and (16) we compute the viscous part τ̄ik′ of the Cauchy stress as ′ ′ j τ̄im = σ̄ik Fmk
= ζ − 32 η Fik Fmk blj djl + 2ηbil dls Fsk Fmk (18)
= ζ − 23 η bim bil dls + 2ηbil dls bsm , or, in absolute notation,
j τ̄ ′ = ζ − 32 η b tr(bd) + 2ηbdb. (19) It is quite clear that this is a polynomial isotropic function of the two tensor variables b and d, which can be shown to be a special case of the universal representation for such functions (see, e.g., the classical treatise by Truesdell and Noll [7, Sect. 13] and the paper by Rivlin[8]). It is reassuring to notice that expression (19) for the viscous part of Cauchy stress τ̄ ′ falls perfectly within the class of constitutive equations for finite viscoelasticity. This makes much more plausible our proposal (11) for a modified version of Eq. (7). 5 2 Wave motion When it comes to study wave propagation, attention is usually focused on special motions, for which the general equations of motion simplify greatly. One might then be led to believe that no relevant corrections should be made on the governing equations as a consequence of the substitution of Eq. (7) with the more appropriate version Eq. (11). In this section we show that this is not the case, indeed, and that some relevant and quite complicated changes should be made to the propagation equations, even for simple (bulk) motions. We discuss the propagation in a soft material of a transverse wave described by a displacement field ui (xk , t) of the form u1 = u(z, t), u2 = u3 = 0, (20) where z is the third Cartesian coordinate. This is the choice made by Zabolotskaya et al.[9], who added a viscous part in the form (7) to an elastic constitutive equation (see their Eq. 34). Straightforward computations show that     1 0 uz 0 0 uz 1 E =  0 0 0 , F = 0 1 0  , (21) 2 2 0 0 1 uz 0 uz (here, for compactness, we denote partial derivatives with subscripts, while in other formulas we use the more explicit notation, so that a comparison with some relevant references is made easier). The time derivative of E and its trace are thus given by   0 0 uzt 1 0 , tr(Ė) = uz uzt , (22) Ė =  0 0 2 uzt 0 2uz uzt and, moreover,  1 0 1 0 1 F Ė = 2 0 0  uu 1  z zt 0 = 2 uzt   0 0 uzt uz 0  0  0 0 uzt 0 2uz uzt 1  0 uzt + 2u2z uzt  0 0 0 2uz uzt 6 (23) Now, since by Eq. (11) σ̄ ′ = (ζ − 23 η) tr(Ė) F + 2ηF Ė, (24) we find that the Cartesian components of the viscous part σ̄ ′ of the PiolaKirchhoff stress tensor are given by   0 ηuzt + (ζ + 34 η)u2z uzt (ζ + 31 η)uz uzt .  0 (ζ − 32 η)uz uzt 0 (25) 4 ηuzt 0 (ζ + 3 η)uz uzt In their Section III, Zabolotskaya et al.[9] investigate the propagation of waves described by Eq. (20) in an almost incompressible soft material with elastic strain energy E = µI2 + 31 AI3 + DI22 , (26) where µ, A, D are second-, third-, and fourth-order elastic constants, respectively. Terms of order higher than the third in the strain ∂u/∂z are then neglected and the governing equations reduce to the single equation  3 ∂ 2u ∂ ∂u ∂ 2u , (27) ρ0 2 = µ 2 + γ ∂t ∂z ∂z ∂z where ρ0 is the constant density and γ = µ + 21 A + D, see Ref. [9, Eq. 12]. ′ Next, in Section V of Ref.[9], a viscous stress is added, in the form of σik as defined here by Eq. (7), and an appropriate modification for the governing equation (27) is deduced, which leads to the introduction of an additional viscous term:  3 ∂ 2u ∂ ∂u ∂ 3u ∂ 2u ρ0 2 = µ 2 + γ +η 2 . (28) ∂t ∂z ∂z ∂z ∂z ∂t} | {z | {z } elastic viscous This equation can also be found in the Refs.[10, 11, 12], for instance. Our aim here is to verify which modifications would be required for the wave equation (27) if the added viscous Piola-Kirchhoff stress were defined by Eq. (11), as we suggest, rather than by Eq. (7). In other words: what difference does it ′ make to the wave propagation equations to consider such a stress tensor σ̄ik , ′ rather than σik ? 7 In order to answer this question we first compute the divergence of σ̄ ′ , as "   2 # 3 2 ∂ u ∂u ∂ u ∂ (Divσ̄ ′ )1 = η 2 + (ζ + 34 η) , ∂z ∂t ∂z ∂z∂t ∂z (29) (Divσ̄ ′ )2 = 0,   ∂ ∂u ∂ 2 u (Divσ̄ ′ )3 = (ζ + 34 η) . ∂z ∂z ∂z∂t Here we remark that the first and last of these components lead to two non-trivial equations of motion, in contrast to the situation encountered in Ref. [9], where there was only one equation. The second equation here may however be made to disappear if we were to consider soft solids to be perfectly incompressible, and had thus to introduce an arbitrary Lagrange multiplier (for instance, see Ref. [13] for a rational inclusion of perfect incompressibility in weakly nonlinear elasticity, and Ref. [14] for a derivation of the equations of motion in compressible and incompressible materials and the possible decoupling of longitudinal from transverse waves.) Next, we readily derive the wave equation "  3   2 # 2 ∂ 2u ∂ ∂u ∂u ∂ u ∂ 3u ∂ 2u ∂ + η 2 + (ζ + 34 η) ρ0 2 = µ 2 + γ . ∂t ∂z ∂z ∂z ∂z ∂t ∂z ∂z∂t ∂z | {z } | {z } elastic viscous (30) The viscous part of this equation is remarkably more complex than the simpler term η∂ 3 u/∂z 2 ∂t shown in Eq. (28), which is Eq. (36) of Ref. [9]. Notice that the last addendum in Eq. (30) is of third order in the strain ∂u/∂z, exactly as the nonlinear elastic term ∂/∂z (∂u/∂z)3 . Thus, since Eq. (27) was obtained in Ref.[9] by neglecting terms of order higher than three, it is our opinion that, in principle, the additional viscous term in (30) should be kept and considered for a coherent discussion of wave propagation. Dimensional analysis reveals that this term would be negligible if fifth-order elastic constants were much larger than lower-order constants. Although this might be the case for some solids, it can be shown by following the steps presented in Ref.[13] and Ref. [18] that for almost incompressible solids, all elastic constants are of the same order as µ (see Ref.[19] for experimental evidence). To illustrate the potential influence of the extra term in Eq.(30), we focus on a staple of acoustic nonlinearity: the finite amplitude, traveling 8 kink solution. First we rewrite the equation of motion for the strain[15]: w ≡ ∂u/∂z, as  
∂2 ∂ 2w ∂ 2w ∂2 ∂ 3w 3 2 ∂w 4 ρ0 2 = µ 2 + γ 2 (w) + η 2 + ζ + 3 η w . (31) ∂t ∂z ∂z ∂z ∂t ∂z 2 ∂t Then we perform the following changes of variables and of function: r √ ρ0 µ µ γ ξ= z, τ = t, W = w, η η µ (32) to obtain
∂2 ∂ 2W ∂2 ∂ 3W µ ∂ 2W 3 4 = + (W ) + ζ + + η 3 ∂τ 2 ∂ξ 2 ∂ξ 2 ∂ξ 2 ∂τ γ ∂ξ 2  W 2 ∂W ∂τ  Looking for a traveling wave solution in the form √ W (ξ, τ ) = c2 − 1 ω(x), where x = (c2 − 1)(τ − ξ/c), . (33) (34) and c > 1 is the arbitrary, non-dimensional speed, we find the following equation for the amplitude ω,


′′
′′ where α ≡ (µ/γ) ζ + 34 η c2 − 1 . (35) ω ′′ = ω 3 + ω ′′′ + α ω 2 ω ′ , It can be integrated twice for a finite kink with tails such that ω(−∞) = 0, ω(∞) = 1, ω ′ (±∞) = ω ′′ (±∞) = 0, to give
(36) ω 1 − ω 2 = (1 + α ω 2 )ω ′ . Now the variables of this equation can be separated, and by centring the kink so that ω(0) = 1/2, we obtain its inverse definition,  α+1 1 4 2 2 . (37) (1 − ω ) 2ω 3 It is a simple matter to construct the ω − x curves and then to invert them to generate the x − ω curves. Explicit inversions include s  2 3 x 1 3 1+ ω(x) = √ , e (38) − e−x , 4 4 1 + 3e−2x x = − ln at α = 0, 1, respectively. At α = 0, there is no extra term in Eq.(30); as soon as α 6= 0, the difference between the incorrect formulation and the proper formulation of viscous effects is felt, as illustrated by Figure 1. 9 Figure 1: Finite amplitude travelling transverse strain kink in a viscous soft solid. The equation of motion has been fully non-dimensionalized. The parameter α is a measure of the consequence of properly incorporating viscous effects in the formulation of the equations of motion. At α = 0 (dotted curve), the viscous formulation does not allow for the principle of objectivity to be respected. As α grows (here, α = 1, 2, 3 in turn), the wave front displays a gentler slope and the wave is more attenuated. 3 Concluding remarks Viscous stresses are often introduced in the formulation of nonlinear wave motion problems in order to prevent the formation of shocks, because they confer a parabolic character to the equations of motion for continua. However, as pointed out earlier by Antman [16], the choices made historically in the literature for these stresses sometimes turn out to be physically unacceptable because their material responses are affected by rigid motions. Here we showed that simply extending the linear Kelvin-Voigt model of differential visco-elasticity from linearized to finite elasticity is not a straightforward process. In particular, one choice is to take the strain-rate effects to be described by the time-derivative of the full Lagrangian strain instead of the infinitesimal strain. With that choice must come great care in formulating a corresponding viscous stress which obeys the fundamental principle of objectivity in the same manner that the elastic stress does. We saw here that this compatibility can be achieved by pre-multiplying the linear expansion of the viscous stress in terms of u̇ by the deformation gradient F , which complicates greatly the equations of motion, as illustrated in Section 2. An alternative constitutive assumption for modeling viscous effects is to take the Cauchy stress tensor to be linear in d, the Eulerian stretching tensor defined in Eq. (14). That assumption is perfectly coherent with the fundamental principles of mechanics, including frame-invariance, and is aligned 10 with modeling of viscous effects in Fluid Mechanics and the emergence of the Navier-Stokes equations. It has also been used to model nonlinear wave propagation in solids[14, 15, 17]. It does not complicate the equations of motion excessively. Acknowledgements We are grateful to the Istituto Nazionale di Alta Matematica (Italy) for support through the PRIN 2009 Matematica e meccanica dei sistemi biologici e dei tessuti molli for the second author and a short-term Visiting Professor position for the first author. We thank the anonymous referees for their constructive comments. References [1] L.D. Landau and E.M. Lifshitz, Theory of Elasticity, 3rd ed. (Pergamon Press, New York, 1986), p.107. [2] L.D. Landau and G. Rumer, “Über Schallabsorption in festen Körpern (On supersonic absorption in solid bodies),” Sov. Phys. B3, 18–27 (1937). [3] F.D. Murnaghan, “Finite deformations of an elastic solid,” Am. J. Math. 59, 235–260 (1937). [4] L. Brillouin, “Sur les tensions de radiation (On tensions of radiation),” Ann. Phys. 4, 528-86 (1925). [5] P. Chadwick, Continuum Mechanics pages 129-135 (Dover Publications, New York, 1999). [6] M.E. Gurtin, An Introduction to Continuum Mechanics pages 139-145 (Academic Press, New York, 1981). [7] C. Truesdell and W. Noll, “The non-linear field theories of mechanics,” In S. Flügge, editor, Handbuch der Physik, volume III/3, pages 1-591 (Springer Verlag, Berlin, 1965). 11 [8] R. S. Rivlin, “Further remarks on the stress-deformation relations for isotropic materials,” J. Rational Mech. Anal. 4, 681-702 (1955). [9] E. A. Zabolotskaya, M. F. Hamilton, Y.i A. Illinski, and G. D. Meegan, “Modeling of nonlinear shear waves in soft solids,” J. Acoust. Soc. Am. 116, 2807-2813 (2004). [10] M.S. Wochner, M.F. Hamilton, Y.A. Ilinskii, and E.A. Zabolotskaya, “Cubic nonlinearity in shear wave beams with different polarizations,” J. Acoust. Soc. Am. 123, 2488-2495 (2008). [11] M. Rénier, J.-L. Gennisson, C. Barrière, D. Royer, and M. Fink, “Fourth-order shear elastic constant assessment in quasi-incompressible soft solids,” Appl. Phys. Lett. 93, 101912 (2008). [12] X. Jacob, S. Catheline, J.-L. Gennisson, C. Barrière, D. Royer, and M. Fink, “Nonlinear shear wave interaction in soft solids,” J. Acoust. Soc. Am. 122, 1917-1926 (2007). [13] M. Destrade and R.W. Ogden, “On the third- and fourth-order constants of incompressible isotropic elasticity,” J. Acoust. Soc. Am. 128, 33343343 (2010). [14] M. Destrade and G. Saccomandi, “On finite amplitude elastic waves propagating in compressible solids,” Phys. Rev. E 72 016620 (2005). [15] M. Destrade and G. Saccomandi, “Solitary and compact-like shear waves in the bulk of solids,” Phys. Rev. E. 73, 065604 (2006). [16] S.S. Antman, “Physically unacceptable viscous stresses,” Zeits. Angew. Math. Phys. (ZAMP) 49 980-988 (1998). [17] M. Destrade, P. M. Jordan, and G. Saccomandi, “Compact travelling waves in viscoelastic solids,” Europhys. Letters 87 48001 (2009). [18] R. W. Ogden, “On isotropic tensors and elastic moduli,” Proc. Cambr. Phil. Soc. 75, 427–436 (1974). [19] S. Catheline, J.-L. Gennisson, M. Fink, “Measurement of elastic nonlinearity of soft solid with transient elastography,” J. Acoust. Soc. Am. 114, 3087–3091 (2003). 12