Creep, recovery, and waves in a nonlinear fiber-reinforced viscoelastic solid
aa r X i v : . [ c ond - m a t . s o f t ] D ec Creep, recovery, and waves in a nonlinearfiber-reinforced viscoelastic solid
M. Destrade, G. Saccomandi2007
Abstract
We present a constitutive model capturing some of the experimen-tally observed features of soft biological tissues: nonlinear viscoelas-ticity, nonlinear elastic anisotropy, and nonlinear viscous anisotropy.For this model we derive the equation governing rectilinear shear mo-tion in the plane of the fiber reinforcement; it is a nonlinear partialdifferential equation for the shear strain. Specializing the equationto the quasistatic processes of creep and recovery, we find that usual(exponential-like) time growth and decay exist in general, but that forcertain ranges of values for the material parameters and for the anglebetween the shearing direction and the fiber direction, some anoma-lous behaviors emerge. These include persistence of a non-zero strainin the recovery experiment, strain growth in recovery, strain decay increep, disappearance of the solution after a finite time, and similarodd comportments. For the full dynamical equation of motion, wefind kink (traveling wave) solutions which cannot reach their assignedasymptotic limit.
Many biological, composite, and synthetic materials must be modeled asfiber-reinforced nonlinearly elastic solids. Hence, the anisotropy due to thepresence of collagen fibers in many biological materials has been studied ex-tensively within the constitutive context of fiber-reinforced materials by sev-eral authors (see for example Humphrey (2002) and the references therein.)In nonlinear elasticity, the macroscopic response of an anisotropic materialis given in terms of a strain-energy function, which itself depends on a set ofindependent deformation invariants. This formulation captures a great vari-ety of phenomena related to the behavior of fiber-reinforced materials such1s inter alia , the examination of fiber instabilities, using loss of ellipticity(see Merodio and Ogden 2002, 2003, and the references therein).Generally speaking, a reinforcement is added to a given material withthe aim of avoiding a possible failure under operating conditions. Thereforeit is important to develop a detailed study showing how to introduce rein-forcements in a material in order to control the possible development of aboundary layer structure. Our goal here is to provide a first step in this di-rection. We make several simplifications and ad hoc assumptions. First, welimit ourselves to the consideration of only one fiber direction and second, weconsider a one-dimensional motion in the bulk of an infinite body. Here themotion is linearly polarized in a direction normal to the plane containing thedirection of propagation and the direction of the fibers. We acknowledge thatmore complex anisotropies, geometries, and couplings arise in biomechanicalapplications. For instance, the mechanics of the aorta involves two families ofparallel fibers, triaxial motions, and blood flow / arterial wall coupling. How-ever we argue that some major characteristics of biological soft tissues areencompassed in the choices of transverse isotropy, of infinite extent, and of amotion governed by an ordinary differential equation. Indeed the anisotropydue to the presence of one family of parallel fibers complicates the governingequations to an extent which is only marginally less than that due to thepresence of two families of parallel fibers. Also, soft biological tissues arenearly incompressible and a (compressive) longitudinal wave is difficult toobserve; it thus make sense to focus on transverse shear motions, which areuseful in imaging technologies. Our third assumption is that the elastic strainenergy is the sum of an isotropic part and an anisotropic part (called a rein-forcing model ), in order to model an isotropic base material augmented by auniaxial reinforcement in the fiber direction . Albeit strong, this constitutiveassumption is now common and used by many authors (e.g. Triantafyllidisand Abeyaratne 1983, Qiu and Pence 1997, Merodio and Ogden 2002). Fi-nally, we assume that the solid is viscoelastic and here we assume not onlyNewtonian viscosity (proportional to the stretching tensor) but also fiber-oriented (anisotropic) viscosity. That latter assumption is strong but can beremoved from our calculations by taking a constant to be zero. We believethat it might be useful in modeling the well-documented physiological effectof stretching training in sport medicine, which is that it affects the viscosityof tendon structures but not their elasticity (Taylor et al., 1990; Kubo et al.,2002).We divided the article into the following sections. Section 2 presents theconstitutive model and the derivation of the equation governing the rectilin-ear shear motion. As expected, this equation is nonlinear in the shear strain:it is a second-order partial differential equation, with cubic nonlinearity. To2nitiate its resolution, we first look at the quasistatic experiment of recoveryin Section 3. Then we have a first-order ordinary differential equation, andwe find that it can lead to unusual behaviors when certain conditions (stronganisotropy, large angle between the shearing direction and the fibers) aremet. The same is true of the case of creep, treated in Section 4. Basically,it turns out that the nonlinearity introduces ranges of material parametersand angles for which an expected behavior – say strain growth in creep –can be turned on its head – and lead to strain decay with time in creep,say. In the course of the investigation we develop synthetic tools of analy-sis which highlight the boundaries of these ranges. They also guide us forthe resolution of the full dynamical equation of motion, which we tackle inSection 5 for traveling wave solutions. Again the solution may behave in anunexpected way, provided the anisotropy is strong enough and the fibers arein compression. Finally, Section 6 recaps the results and puts them into awider context.
We describe the motion of a body by a relation x = x ( X , t ), where x denotesthe current coordinates of a point occupied by the particle of coordinates X in the reference configuration at the time t .We introduce F = ∂ x /∂ X , the deformation gradient, and C = F T F ,the right Cauchy-Green strain tensor. We focus on incompressible materials for which all admissible deformations must be isochoric, or equivalently, forwhich the relation det F = 1 must hold at all times.The body is reinforced with one family of parallel fibers. Our first as-sumption is that the unit vector a , giving the fiber direction in the referenceconfiguration, is independent of X . The stretch along the fiber direction is √ a · Ca = √ a · a where a = Fa .We may now introduce the elastic part of our constitutive model. Weconsider the so-called standard reinforcing model , which is a quite simplegeneralization to anisotropy of the neo-Hookean model (Triantafyllidis andAbeyaratne, 1983; Qiu and Pence, 1997). For the standard reinforcing model,the strain-energy density is given by W = µ (cid:2) ( I −
3) + γ ( I − (cid:3) , where I = tr C , I = a · Ca = a · a . (2.1)Here µ > γ > elastic anisotropy parameter , and the invariant I measures the squared stretch in the fiber direction. Mechanical tests showthat the neo-Hookean strain energy function µ ( I − / γ ( I − is adequate to describe a reinforced material which penalizesdeformation in the fiber direction (Merodio and Ogden, 2003).The spatial velocity gradient L ( X , t ) associated with a motion is definedas L = grad v , where v = ∂ x /∂t is the velocity, and the stretching tensor D is defined as D = ( L + L T ). For incompressible materials, tr D = 0at all times. Newtonian viscous fluids possess a constitutive term in theform 2 ν D , where ν is a constant. For our special solid, we modulate theNewtonian viscosity with an anisotropic term, by replacing ν with ν [1 + γ ( I − γ > viscous anisotropy parameter . We showin the course of the paper that this simple choice of anisotropic viscositycaptures the essential characteristics of attenuation in soft biological fibroustissues. According to Baldwin et al. (2006), ultrasonic measurements offreshly excised myocardium show that “the attenuation coefficient was foundto increase as a function of frequency in an approximately linear mannerand to increase monotonically as a function of angle of insonification from aminimum perpendicular to a maximum parallel relative to the direction ofthe myofibers.”We are now ready to give the complete Cauchy stress tensor of our vis-coelastic, transversally isotropic material as T = − p I + µ [ B + γ ( I − a ⊗ a ] + 2 ν [1 + γ ( I − D , (2.2)where the p is the yet indeterminate Lagrange multiplier introduced by theincompressibility constraint, and B = FF T is the left Cauchy-Green tensor. We take a fixed, orthonormal triad of vectors ( i , j , k ), and call X , Y , Z thereference coordinates; hence X = X i + Y j + Z k . The triad is such that theunit vector in the fiber direction lies in the XY plane; hence, a = cos θ i + sin θ j , (2.3)(say) where θ ∈ [0 , π ] is the angle between the X -axis and the fibers.We then consider the rectilinear shearing motion , x = X + u ( Y, t ) , y = Y, z = Z, (2.4)4here the anti-plane displacement u is real and finite. Then the componentsof the gradient of deformation F and of its inverse are given by F = U
00 1 00 0 1 , F − = − U
00 1 00 0 1 , (2.5)where U = ∂u/∂Y is the amount of shear . The left and right Cauchy-Greentensors are thus B = U + 1 U U , C = U U U + 1 00 0 1 , (2.6)respectively, from which the expressions of the invariants I and I follow, I = 3 + U , I = 1 + U sin 2 θ + U sin θ. (2.7)Figure 1a shows the variations of I with θ for several values of U between0 and 1. When I > I < I with U for several values of θ ; when 0 < θ < π/
2, the fibers are always in extensionand when π − tan − (2) = 2 . < θ < π , they are always in compressionfor 0 ≤ U ≤
1. We refer to the paper by Qiu and Pence (1997) for similarfigures and closely related discussions.In the deformed configuration, we find that a = (cos θ + U sin θ ) i + sin θ j .The remaining tensors required to compute the Cauchy stress tensor (2.2)are a ⊗ a = (cos θ + U sin θ ) (cos θ + U sin θ ) sin θ θ + U sin θ ) sin θ sin θ
00 0 0 , (2.8)and D = 12 U t U t , (2.9)so that the non-zero components of T are T = − p + µ and T = − p + µ (1 + U ) + µγ ( I − θ + U sin θ ) ,T = − p + µ + µγ ( I −
1) sin θ,T = µU + µγ ( I − θ + U sin θ ) sin θ + ν [1 + γ ( I − U t . (2.10)5 I U = U = 0.75 U = 0.5 U = 0.25 q (a) I U p /29 p /16 p - arctan(2) p /4 (b) Figure 1: Variations of the squared stretch in the fiber direction: (a) withthe angle and (b) with the shear. When I >
1, the fibers are in extension;when I <
1, they are in compression.Now the equations of motion div T = ρ x tt reduce to the two scalar equa-tions: − p x + T ,y = ρu tt and − p y + T ,y = ρu tt . Differentiating the formerwith respect to y and the latter with respect to x , and eliminating p xy , wearrive at a single governing equation for the rectilinear shear motion, ρU tt = µU yy + µγ sin θ [ U (2 cos θ + U sin θ )(cos θ + U sin θ )] yy + νU tyy + νγ sin θ [ U U t (2 cos θ + U sin θ )] yy . (2.11)Using the scalings ˜ t = µt/ν and ˜ y = y/L (where L is a characteristic lengthto be specified later on a case-by-case basis), we write this equation in di-mensionless form as εU ˜ t ˜ t = U ˜ y ˜ y + γ sin θ [ U (2 cos θ + U sin θ )(cos θ + U sin θ )] ˜ y ˜ y + U ˜ t ˜ y ˜ y + γ sin θ [ U U t (2 cos θ + U sin θ )] ˜ y ˜ y , (2.12)where ε = ρµL /ν . This is the main equation of our study. For conveniencewe drop the tildes in the remainder of the paper. We also introduce thefollowing functions, f ( γ , U, θ ) = 1 + γ sin θ (2 cos θ + U sin θ )(cos θ + U sin θ ) ,g ( γ , U, θ ) = 1 + γ U sin θ (2 cos θ + U sin θ ) , (2.13)so that (2.12) is now εU tt = [ U f ( γ , U, θ ) + U t g ( γ , U, θ )] yy . (2.14)6 Nonlinear anisotropic recovery
Our first investigation is placed in the quasistatic approximation, where westudy the influence of elastic anisotropy and viscous anisotropy on the classicexperiment of viscous recovery . We imagine that the material is sheared andthat at t = 0 the shear stress is removed: T (0) = 0. Here the characteristiclength L is the displacement at t = 0 from which the material will relax tothe unstressed state.In the quasistatic case, we neglect the inertia term of (2.12) and maythus integrate it twice to give the following first-order ordinary differentialequation, U f ( γ , U, θ ) + U t g ( γ , U, θ ) = 0 . (3.1)Here we take the constants of integration to be zero, according to the contextof recovery, as explained above. We then solve the equation as Z g ( γ , U, θ ) U f ( γ , U, θ ) d U = − t + const. , (3.2)where the constant is computed so that U (0) = 1.When θ = 0, the fibers are not active with respect to the deformation andwe recover the classical result of isotropic viscoelastic recovery: U ( t ) = e − t .When θ = π/
2, the anisotropic effects are at their strongest. In that casethe integral above has a compact expression and we find: U (cid:20) γ U γ (cid:21) “ γ γ − ” = e − t . (3.3)We now take γ = 0 (no anisotropic viscosity) and γ = 1 , ,
100 (recallthat the fibers are inextensible in the limit γ → ∞ ). Figure 2a shows thatas the anisotropic effect becomes more pronounced, the recovery is quicker;in other words, the influence of elasticity becomes stronger as γ increases.Then we fix γ at 5 for instance, and look at the role played by the anisotropicviscosity, by taking in turn γ /γ = 0 .
5, 1 .
5, 2.5. We find in Figure 2b that,as expected, the viscous recovery is slower as γ increases.When θ = 0, θ = π/
2, other behaviors arise, which call for a detailedanalysis. In particular, the exponential, or near-exponential, decay towardzero as t → ∞ is not necessarily ensured, especially when the anisotropiceffects are strong and the fibers are oriented at a large angle θ > π/
2. Clearly, U t = 0 when f = 0, according to (3.2). Also, U t < f and g are ofthe same sign and U t > f and g are of opposite signs. These twofunctions are quadratic in U . If they have no real roots in U , then they are7 e -t γ = 1 γ = γ = 100 γ = 0 U(t) t (a) γ / γ = 2.5 γ / γ = 1.5 e γ / γ = 0.5 γ = 5 U(t) t (b)
Figure 2: Time recovery function when θ = π/
2: (a) γ = 0 and γ = 1, 5,100; (b) γ = 5 and γ = 0 .
5, 1 .
5, 2 .
5. The recovery function for an isotropicsolid is also plotted (dotted curve).both of the positive sign, and U t < < θ < π/ U might be an increasing function of t . This happens for f and for g when π/ < θ < π and γ ≥ θ sin θ , γ ≥ θ , (3.4)respectively. In Figure 3, the region C corresponds to the first inequality,where the delimiting curve has a vertical asymptote at θ = π/
2, a verticalasymptote at θ = π , and a minimum at θ = π/ γ = 16; we recall thatQiu and Pence (1997) showed that when γ >
16, “simple shear at certainfiber orientations involves negative shear stress in the shearing direction forcertain positive shears.” The regions B and C correspond to the secondinequality, where the delimiting curve has a vertical asymptote at θ = π/ γ = 1. First, we take both γ and γ in region A. This is the simplest case because f and g are then both positive and so U t is always negative (damped recovery).We took several representative examples in this region (say θ = π/ γ = 20, γ = 1) and checked, through integration and implicit plotting, that thegraphs are indeed of the same nature as those in Figure 2.8 γ , γ θ BCA
Figure 3: Recovery: regions where the sign of U t may change. Second, we take γ in region A, by fixing it at γ = 1, say. In that region, g > U t is the opposite of the sign of f . Thenwe take γ = 20, which is above the minimum of region C. In Figure 4, weplotted the locus for the values of U as functions of θ such that f (20 , U, θ ) = 0.Outside the resulting oval shape, f >
0, and inside, f <
0. We also plottedthe line U = 1, which intersects the oval at θ min = 2 .
136 and θ max = 2 . U (0) = 1.When θ > θ max , U ( t ) starts at 1, and decreases because f > U t <
0; as U decreases toward 0, U t tends to zero according to (3.2) , buttakes an infinite time to do so, according to (3.2) ; hence U = 0 is a horizontalasymptote and the recovery is “classical”, see plot (i) in Figure 4, traced at θ = 2 . U clearly changes sign as t increases, in contrast withe − t , traced in dotted line.)When θ min < θ < θ max , U ( t ) starts at 1 and then grows until it hits theupper side of the oval, taking an infinite time to do so; then this upper boundgives a horizontal asymptotic value, above the initial value , see plot (ii) inFigure 4, traced at θ = 2 .
2. 9 θ maxmin θθ m U > t (i) (ii)(iii)(iv) U U UU tttt U Figure 4: Types of time recovery functions for γ = 20, γ = 1 (strong elasticanisotropy). The amount of shear U starts at 1 for t = 0. Outside the ovalshape, U t < U decreases as in ( i ) and ( iv ): decay toward zero; andin ( iii ): decay toward a value >
0. Inside the oval shape, U t > U increases as in ( ii ): growth toward a value >
1. The recovery function e − t for an isotropic solid is also plotted (dotted curves).When θ m < θ < θ min , where θ m = 2 .
124 is the angle at which the ovalplot has a vertical tangent, U ( t ) starts at 1 and then decreases until it hitsthe upper side of the oval, below 1 but above 0; then this lower bound givesa horizontal asymptotic value, above zero , see plot (iii) in Figure 4, traced at θ = 2 . θ < θ m , U ( t ) starts at 1 and then decreases until zero; thenthis lower bound gives zero as a horizontal asymptotic value, see plot (iv) inFigure 4, traced at θ = 2 .
05. Notice that the second derivative changes signsthree times as t increases. 10 .3 Strong viscous anisotropy Third, we take γ outside the C region, by fixing it at γ = 1, say. In thatregion, f > U t is the opposite of the sign of g . Then we allow γ to be in region B, and thus allow g (and U t ) to changesign with increasing θ , by taking γ = 3 . U as functions of θ such that g (3 , U, θ ) = 0 and obtained the thick-line shape. Outside the shape, g >
0, and inside, g <
0. We also plottedthe horizontal line U = 1, which intersects the shape at θ min = 2 .
356 and θ max = 2 . θ = θ m = 2 . θ < θ m or θ > θ max , U ( t ) starts at 1 and decreases until zero;as U →
0, the denominator in the integral tends to zero, indicating that ittakes an infinite time to do so; hence, zero is a horizontal asymptote in thesecases. To draw Figure 5(i), we took θ = 3 . θ = 2 .
1; both graphs show a somewhat classical decay with time.However, when θ min < θ < θ max , U ( t ) starts at 1 and then grows because U t > U hits the upper face of theshape, where g = 0; then by (3.1), either U f = 0, or U t → ∞ . Clearly,the first possibility is excluded because U = 0 when it is larger than 1, and f = 0 when γ is outside the C region. It follows that U grows and hits theupper face of the shape with a vertical asymptote after a finite time (andthen stops because it cannot increase further since U t < U t = 0 on the shape, and it cannot decreasesince U t > U ( t ),traced at θ = 2 . θ m < θ < θ min , U ( t ) decays from 1 until it hits the shapefrom above after a finite time; see Figure 5(iii), traced at θ = 2 .
2. Noticehow quickly the final value is reached, compared to the isotropic exponentialrecovery.
In the case where both γ and γ are in the region C, any combinationand overlaps of the thick curves presented in Figures 4 and 5 may arise. Thetools presented in the two previous subsections are easily transposed to thosepossibilities. A special situation arises when the locus of f = 0 intersects thelocus of g = 0; then, the numerator and the denominator in (3.2) may havea common factor so that the integrand simplifies and a regular behavior mayappear. This situation is however too special to warrant further investigation,and we do not pursue this line of enquiry.11 (i)(ii)(iii)(iv) U t U ttUU t U q m q min max q U > t Figure 5: Types of time recovery functions for γ = 1, γ = 2 (strong viscousanisotropy). The amount of shear U starts at 1 for t = 0. Outside the thick-line shape, U t < U decreases toward zero as in ( i ) and ( iv ). Insidethe thick-line shape, U t > U increases rapidly as in ( ii ), until it ceasesto exist. There is also an angular region θ m < θ < θ min where U decreasesrapidly until it ceases to exist, see ( iii ). The recovery function e − t for anisotropic solid is also plotted (dotted curves).12 Nonlinear anisotropic creep
Our second investigation is again placed in the quasistatic approximation,where we now study the influence of elastic anisotropy and viscous anisotropyon the classic experiment of viscous creep . As the resulting analysis is similarto that conducted for recovery, we simply outline the main results.We imagine that the material is sheared and that the shear stress ismaintained: T ( ∞ ) = 0. Here the characteristic length L is an asymptoticvalue of the displacement. We neglect the inertial term of (2.12) and integrateit twice to give the following ordinary differential equation, U f ( γ , U, θ ) + U t g ( γ , U, θ ) = const. , (4.1)where we took the constant of the first integration to be zero, and the con-stant of the second integration to correspond to the applied (constant) shearstress, as is usual in the creep problem. More specifically, this constant istaken so that U ( ∞ ) = 1, and so is equal to f ( γ , , θ ); it follows that theequation above can be written as h ( γ , U, θ )( U −
1) + g ( γ , U, θ ) U t = 0 , (4.2)where h is defined by h ( γ , U, θ ) = [ U f ( γ , U, θ ) − f ( γ , , θ )] / ( U − γ sin θ [1 + cos θ + ( U + 1) sin θ ( U sin θ + 3 cos θ )] . (4.3)We then solve the equation as Z g ( γ , U, θ )( U − h ( γ , U, θ ) d U = − t + const. , (4.4)where the constant is computed so that U (0) = 0. Hence the equationsgoverning creep are almost identical to those governing recovery, with thedifference that f is now replaced by h .Here we are mostly concerned with the question of how, if at all, a state ofshear can be reached such that once removed, the unusual recovery behaviorsof the previous section emerge. Thus we concentrate on strong anisotropiceffects, with emphasis on strong elastic anisotropy (where the new function h is involved). We traced the regions where g and h , and thus U t , may changesigns and found that the resulting graph is similar to that of Figure 3, withthe main difference that the minimum of region C is now located at θ = 3 π/ γ = 4. Thus unusual behavior in creep may occur at much lower levelsof elastic anisotropy than in recovery (where the minimum is at γ = 16).We recall that Qiu and Pence (1997) show that when γ >
4, “simple shearat certain fiber orientations involves a nonmonotonic relation between theshear stress in the shearing direction and the amount of shear.”13 .1 Strong elastic anisotropy
We begin with the case where h plays a major role, that is when γ is greaterthan 4. For the purpose of direct comparison with the recovery problem, wetake γ = 20 and γ = 1, as in Section 3.2. Figure 6 displays the curve where h (20 , U, θ ) = 0. Outside the thick-line curve, U t >
0, and inside, U t <
0. Thecurve intersects the line U = 0 twice, at θ min = 2 .
136 and at θ max = 2 . f = 0 intersects U = 1 in Section 3.2 (seethin-line shape), because by (4.3), h (20 , , θ min ) = f (20 , , θ min ) = 0, andsimilarly h (20 , , θ max ) = f (20 , , θ max ) = 0. We also display the verticallines θ = θ M = 2 . h = 0 intersects U = 1, and θ = θ m = 2 . h = 0 has a vertical tangent. Recall that for creep, U (0) = 0.When θ > θ M , U ( t ) starts at 0 and grows toward 1; then U t tends tozero according to (4.2), but takes an infinite time to do so; hence U = 1 isa horizontal asymptote and the creep is “classical”, see plot (i) in Figure 6,traced at θ = 2 . − e − t ) is shown in dotted line.)When θ max < θ < θ M or when θ m < θ < θ min , U ( t ) starts at 0 and thengrows until it hits the oval shape, taking an infinite time to do so; then thisupper bound gives a horizontal asymptotic value, below 1 , see plot (ii) inFigure 6, traced at θ = 2 .
5, and plot (iv), traced at θ = 2 . θ min < θ < θ max , U ( t ) starts at 0 inside the oval shape and thusit decreases until it hits the lower side of the shape, taking an infinite timeto do so; then this lower bound gives a horizontal asymptotic value, below 0 ,see plot (iii) in Figure 6, traced at θ = 2 . θ < θ m , U ( t ) can again grow toward 1, see plot (v) in Figure6, traced at θ = 2 .
0. Notice however that the concavity of the curve changesas t increases. Here we remark that the function governing the strength of the viscousanisotropy, namely g , is the same for creep as it is for recovery. Thus, theregion where U t might change sign because of strong viscous anistropy isthe region B of Figure 2. Also, the locus of points where g = 0 is typicallydisplayed by the thick-line shape of Figure 5 and because g ( γ , , θ ) = 1 > U = 0. It follows that there isonly one situation where viscous anisotropy leads to anomalous creep, when θ min < θ < θ max ; then, U ( t ) starts at zero, and grows toward the thick-lineshape, which it reaches after a finite time with a vertical asymptote.14 U < t θ M θ m θ min θ max –0.008–0.006–0.004–0.0020 0.1 0.2 0.3 0.400.10.20.30.40.50.60.70.8 0.2 0.4 0.6 0.8 100.20.40.60.81 0.2 0.4 0.6 0.8 1 1.2 (i)(ii)(iii)(iv)(v)U U t U t U tU tU t
Figure 6: Types of time creep functions for γ = 20, γ = 1 (strong elasticanisotropy). The amount of shear U starts at 0 for t = 0. Outside the thick-line shape, U t > U increases as in (i) and (v): growth toward 1, and asin (ii) and (iv): growth toward a value below 1. Inside the thick-line shape, U decreases as in (iii): decay toward a negative value. The creep function1 − e − t for an isotropic solid is also plotted (dotted curves).15 .3 Pre-stretch and nonlinear anisotropic creep Here we show how anomalous creep can be avoided (amplified) by stretch-ing (compressing) the solid prior to the shear. Hence, instead of (2.4), weconsider the motion x = λ − X + λu ( Y, t ) , y = λY, z = λ − Z, (4.5)The following decomposition of the associated deformation gradient showsthat the solid is stretched by a ratio λ in the Y direction, F = F F , where F = U
00 1 00 0 1 , F = λ − λ
00 0 λ − . (4.6)(Note that F F = F F .) The kinematic quantities of Section 2.2 aremodified accordingly. In particular, I = λ − cos θ + λ sin θ + U λ sin 2 θ + U λ sin θ. (4.7)The end result is that the differential equation governing creep is changedfrom (4.2) to h λ ( γ , U, θ )( U −
1) + g λ ( γ , U, θ ) U t = 0 , (4.8)where h λ and g λ are defined by h λ ( γ , U, θ ) = λ (cid:8) γ sin θ [2 λ sin θ + 3 λ − cos θ − U + 1) sin θ ( U sin θ + 3 cos θ )] } ,g λ ( γ , U, θ ) = 1 + γ ( λ − cos θ + λ sin θ − U λ sin 2 θ + U λ sin θ ) . (4.9)Figure 7 shows the loci of h λ = 0 in the case of a strong elastic anisotropy( γ = 30, γ = 1), for several values of λ . The figure clearly shows that thepre-stretch λ can be used to control the shape of these curves: if the solidis put in compression first, and sheared for creep next, then the region ofpotential anomalous creep is increased; if it is put under tension, then thearea of the region rapidly decreases and eventually dissapears altogether. So far we have looked at how the presence of elastic and viscous fibers affectssome quasi-static processes. Typically, creep and recovery connect one state16 l = 2l = 2.7l = 0.8 l = 1 U q Figure 7: Effect of pre-stretch on time creep functions for γ = 30, γ = 1(strong elastic anisotropy). The amount of shear U starts at 0 for t = 0.Inside the thick-line shapes, U decreases; this situation may arise when thesolid is not pre-stretched ( λ = 1) or when it is compressed ( λ = 0 . U t > U increases; this situation may arise whenthe solid is in extension ( λ = 2, λ = 2 . dynamically ,by looking for traveling wave (kink) solutions.The mathematical theory of one-dimensional transverse traveling wavesin isotropic viscoelastic materials with a Kelvin-Voigt type of constitutiveequation is well grounded, see for example Nishihara (1995) for a clear andcomplete mathematical approach, or Jordan and Puri (2005) for a specificand explicit example. A traveling wave is a solution to the equations ofmotion in the form U ( Y, t ) = U ( ξ ) , ξ = Y − ct, (5.1)where c is the constant speed; also, U is such thatlim ξ →−∞ U ( ξ ) = U L , lim ξ →∞ U ( ξ ) = U R , (5.2)where U L and U R are distinct constants. In what follows, we focus on the casewhere U L = 0, U R = 1. This case is general up to a rigid translation. Herewe take the displacement corresponding to U R as the characteristic length L .Substituting (5.1) into (3.2), we obtain εc U ′′ = ( U f − cU ′ g ) ′′ , (5.3)and then by integration, cU ′ g = (cid:0) f − εc (cid:1) U + const. (5.4)By the requirement U L = 0, the constant must be zero. By the requirement U R = 1, we have f ( γ , , θ ) = εc . (5.5)This equation prompts three remarks.First, we must ensure that f ( γ , , θ ) >
0. Recall that according to (2.13), f ( γ , , θ ) = 1 + γ sin θ (2 cos θ + sin θ )(cos θ + sin θ ) , (5.6)and so, ∂f ( γ , , θ ) /∂θ = γ sin θ (4 cos θ + 9 sin θ cos θ − θ ) . (5.7)In Figure 8a we plotted the variations of [ f ( γ , , θ ) − /γ with θ , as wellas those of its derivative with respect to θ (scaled to 1 / θ , has an absolute minimum and anabsolute maximum. The minimum is at ˆ θ , say, such that tan ˆ θ is that root18 [ f ( g , 1, q) - [1/(8 g )]¶ f /¶q q ~ q Ù (a) U qq (b) Figure 8: (a) Variations with θ of [ f ( γ , , θ ) − /γ and of its derivative,showing an absolute minimum at ˆ θ = 2 . − / tan θ − θ , crossing the abscissa line at θ = 1 . x − x = 0 corresponding to π/ < ˆ θ < π ; numerically,ˆ θ = 2 . f ( γ , , ˆ θ ) = 0 for γ , we find that f ( γ , , θ ) > < γ < ˆ γ = 18 . γ > ˆ γ , there appears a rangefor θ where f ( γ , , θ ) > γ and a given θ , the wavetravels with speed c = ± p f ( γ , , θ ) /ε. (5.8)This is of course expressed in the dimensionless variables of length /L andtime × µ/ν . Turning back, if required, to physical variables, we would findthat the wave travels with the dimensional speed p µf ( γ , , θ ) /ν .The second remark is that according to (5.5) and (5.6), the wave (whenit exists) travels with maximum speed at the angle ˜ θ , say, such that tan ˜ θ isthat root of the cubic 4 + 9 x − x = 0 corresponding to 0 < ˜ θ < π/
2; numer-ically, ˜ θ = 1 . µ , γ , and γ are. This observation indicates the way for an acoustic determi-nation of the fiber orientation: if an experimental measurement of the shearwave speed can be made in every direction of a fiber-reinforced viscoelasticnonlinear material, then the fibers are at an angle ˆ θ from the direction of theslowest wave and at an angle ˜ θ from the direction of the fastest wave. We re-call that for waves in an isotropic deformed neo-Hookean material, Ericksen(1953) found that the fastest waves propagate along the direction of greatest19nitial stretch.The third remark is that when (5.5) holds, then f ( γ , U, θ ) − εc = γ U ( U −
1) sin θ [ U sin θ + 3 cos θ + sin θ ] . (5.9)Then the separation of variables, followed by integration of the first orderdifferential equation (5.4), leads to Z g ( γ , U, θ ) U ( U −
1) sin θ [ U sin θ + 3 cos θ + sin θ ] d U = γ c ξ + const. , (5.10)where the constant of integration is arbitrary; without loss of generality, wetake it to be such that U (0) = 1 / U ′ changes sign and it might not be possibleto find a solution satisfying the requirements (5.2)). We may take care ofthe numerator’s sign by considering elastic anisotropy only ( γ = 0) anddiscarding viscous anisotropy ( γ = 0); then g = 1. For the denominatorhowever, we note that U sin θ + 3 cos θ + sin θ can change sign for certainranges of U and θ . Figure 8b shows the curve U = − / tan θ −
1; on its leftside, the denominator is positive, on its right side, it is negative. Accordingly,the wave connects 0 to 1, see in Figure 9a, or is unable to do so, see Figure 9b.In that latter case, the wave front grows toward an asymptotic value whichis less than 1; a second solution exist (dotted curve) with 1 as an asymptoticvalue, but in the ξ → −∞ direction.As a final remark, we note that when γ is large enough to allow for thepossibility that g = 0 (strong viscous anisotropy), then a “singular barrier”arises, see Pettet et al. (2000). In the course of this investigation on nonlinear anisotropic creep, recovery,and waves for fiber-reinforced nonlinear elastic materials, we unearthed somecomplex mechanical responses. For some range of the constitutive parametersand for some angle ranges of the fiber arrangement, we saw that unusual andpossibly aberrant behaviors can emerge.From a mathematical point of view, we gave a detailed explanation ofthe reasons for these behaviors, by linking them to the singularities of thedetermining equations for the amount of shear.From the mechanical point of view, we pointed out that non-standardbehaviors always occur when the angle between the fiber family and the20 U (g / c )x (a) U (g / c )x (b) Figure 9: Traveling wave solution for anisotropic elasticity ( γ = 0); (a) at θ = 1 .