Initial stress symmetry and its applications in elasticity

Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 rspa.royalsocietypublishing.org Initial stress symmetry and its applications in elasticity A. L. Gower1 , P. Ciarletta2,3 and M. Destrade1,4 Research Cite this article: Gower AL, Ciarletta P, Destrade M. 2015 Initial stress symmetry and its applications in elasticity. Proc. R. Soc. A 471: 20150448. http://dx.doi.org/10.1098/rspa.2015.0448 Received: 30 June 2015 Accepted: 1 October 2015 Subject Areas: biomechanics, mechanics, applied mathematics Keywords: residual stress, initial stress, biomechanics, elasticity, constitutive equations Author for correspondence: A. L. Gower e-mail: arturgower@gmail.com 1 School of Mathematics, Statistics and Applied Mathematics, National University of Ireland Galway, University Road, Galway, Republic of Ireland 2 CNRS and Institut Jean le Rond d’Alembert, UMR 7190, Université Paris 6, 4 place Jussieu case 162, Paris 75005, France 3 MOX - Politecnico di Milano and Fondazione CEN, piazza Leonardo da Vinci 32, Milano 20133, Italy 4 School of Mechanical and Materials Engineering, University College Dublin, Belield, Dublin 4, Republic of Ireland An initial stress within a solid can arise to support external loads or from processes such as thermal expansion in inert matter or growth and remodelling in living materials. For this reason, it is useful to develop a mechanical framework of initially stressed solids irrespective of how this stress formed. An ideal way to do this is to write the free energy density Ψ in terms of initial stress τ and the elastic deformation gradient F, so we write Ψ = Ψ (F, τ ). In this paper, we present a new constitutive condition for initially stressed materials, which we call the initial stress symmetry (ISS). We focus on two consequences of this condition. First, we examine how ISS restricts the possible choices of free energy densities Ψ = Ψ (F, τ ) and present two examples of Ψ that satisfy the ISS. Second, we show that the initial stress can be derived from the Cauchy stress and the elastic deformation gradient. To illustrate, we take an example from biomechanics and calculate the optimal Cauchy stress within an artery subjected to internal pressure. We then use ISS to derive the optimal target residual stress for the material to achieve after remodelling, which links nicely with the notion of homeostasis. 1. Introduction Electronic supplementary material is available at http://dx.doi.org/10.1098/rspa.2015.0448 or via http://rspa.royalsocietypublishing.org. When all loads are removed, a body can still hold a significant amount of internal stress, called the residual stress. In manufacturing, residual stress has long been 2015 The Author(s) Published by the Royal Society. All rights reserved. Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 2 ................................................... rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 noted to be detrimental to, or enhance, the performance of a material. For biological tissues, residual stress is used to self-regulate stress and strain, and ultimately preserve ideal mechanical conditions for the tissue [1,2]. In geophysics, due to gravity and active tectonics, the Earth has developed extremely high initial stresses within its crust, which greatly influence the propagation of seismic waves. Here, we use the term initial stress to broadly mean the internal stress of some reference configuration, irrespective of how the stress was formed or the boundary conditions. In this sense, residual stress is a type of initial stress. The initial stress felt by any region of a material is due to the push and pull of the surrounding regions. If any region were to be cut out from the material, the stress on its newly formed boundary would be zero, thus reducing the potential energy in the bulk. Based on this concept, Hoger developed constitutive laws for residually stressed materials (e.g. [3–5]). Hoger showed that by taking this idea to its limit and cutting the material into possibly an infinite number of disconnected regions, the material may be relieved of all of its internal stress in a configuration called the virtual stress-free state. From that configuration, a hyperelastic energy can be defined as a function of the strain from the virtual state to the current configuration. Though the use of this virtual state is technically sound, it leads to challenging calculations even for simple deformations and rarely yields analytic results, unless great simplifications are assumed. Moreover, the experimental identification of the virtual state requires cutting the material, which is not always suitable, especially for living organisms. However, using a virtual stress-free reference is routinely seen as the only viable option; quoting Chuong & Fung [6]: ‘To characterize the arterial wall or any other biological soft tissue, we need a stressfree state’. Conversely, we believe that by developing tools to work directly with initially stressed reference configurations, without the need for a stress-free state, will be very useful, especially in biomechanics. An ideal way to account for the initial stress would be to have a free energy density function Ψ = Ψ (F, τ ) written explicitly in terms of the deformation gradient F and the initial stress τ , without any a priori restrictions. For the development of constitutive laws, there is no need to distinguish between residual stresses and initial stresses, a view which is shared with Merodio et al. [7]. The initial stress τ could then be determined from elastic wave speeds [8–10] or by solving the linear equations of momentum balance. Having the free energy in the form Ψ (F, τ ) has been used to investigate how residual stress affects elastic deformations [7,11] and for a mass-growth framework [12]. Wang et al. [13] found that including one residual stress invariant in the free energy density was a simple way to model the effects of residual stress in the myocardium. Shams et al. [14,15] made progress towards a general framework for initially stressed solids. However, a major obstacle still remains: how to write Ψ (F, τ ) explicitly in terms of the combined invariants of F and τ ? (Even when the source of anisotropy is only due to the initial stress, the free energy still depends on 10 independent invariants [14].) Johnson & Hoger [5] developed representations for the Cauchy stress response σ = ς̂ (F, τ ) in terms of F and τ by first assuming a stress-free virtual state and then inverting numerically the residual stress–strain equation. However, this approach often requires solving numerical nonlinear implicit equations. In [5,16] it was applied to a material with virtual state composed by a Mooney–Rivlin strain energy. In this work, we introduce a new constitutive requirement on Ψ (F, τ ) called the initial stress symmetry (ISS). ISS restricts the constitutive form of the stress σ = ς̂(F, τ ) by providing a constitutive equation for the initial stress τ = ς̂ (F −1 , σ ) that must hold for every F and τ . To our knowledge, this symmetry has never been discussed before. The ISS also helps answer an important question: how much can the Cauchy stress be altered by adjusting the initial stress? From a modelling perspective, initial stress has been used to make the material more or less compliant [16], to control the Poynting effect [7] or to maintain an ideal internal stress [1]. Given a Ψ (F, τ ) that satisfies ISS, then τ = ς̂(F −1 , σ ) suggests that, for any choice of σ , there will exist an initial stress τ that supports σ . Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 A common approach to model the effect of residual stress is to consider a virtual stress-free configuration B̃, from which the material is deformed and ‘glued’ together to produce a residually ◦ stressed equilibrium state B. We will use this approach to model an initially stressed material. An elastic stored energy density Ψ can then be defined as a function of the deformation gradient F̃ from B̃ to the current configuration B so that Ψ = Ψ (F̃). See figure 1 for a diagram of all the configurations. In this paper, we write the free energy density Ψ as a function of the initial stress τ and of the ◦ ◦ deformation gradient F : B → B, so that Ψ = Ψ (F, τ ). Note that B is not necessarily an unloaded configuration. We argue that it is natural to consider that initial stress contributes to the potential energy stored by a material. One extreme example is the Wapa tree, which has been know to burst open once cut, possibly causing injury, due to its immense level of residual stress [18] (figure 2). We assume that F is a purely elastic deformation, but we do not make any assumptions about the origins of the initial stress τ , except that τ affects the stored energy density Ψ . Assuming that the body is incompressible, i.e. J ≡ det F = 1 at all times, the Cauchy stress tensor σ reads [19]: σ =F ∂Ψ (F, τ ) − pI, ∂F (2.1) where p is the Lagrange multiplier associated with the constraint of incompressibility and I is the identity. If the material is compressible, then p is replaced by −2I3 ∂Ψ/∂I3 , which completely expresses the dependency of Ψ on I3 ; for clarity see equation (2.6). Note that we have and will ◦ omit the possible dependence of Ψ on the position X ∈ B for the sake of simplicity. For the body ◦ in the configuration B, we have F = I and σ = τ ; thus, we require that τ= ◦ ∂Ψ (I, τ ) − pI, ∂F (2.2) ◦ where p is the value of p when F = I. We call the above equation the initial stress compatibility. The presence of initial stress generally leads to an anisotropic response of the material in ◦ reference to B. Here, we assume no other source of intrinsic anisotropy so that Ψ can be written as a function of all the independent invariants generated by τ and C = F T F, the right Cauchy–Green ................................................... 2. Initially stressed elastic materials 3 rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 The basic equations for an elastic material subject to initial stress are summarized in §2, and then the ISS is presented. ISS is satisfied automatically if the initial stress is due to an elastic deformation of a stress-free configuration, which we demonstrate in §2b. For this reason, we develop in §3 an example for Ψ (F, τ ) that satisfies ISS, by deforming an incompressible neo-Hookean material from a stress-free virtual state. In §4, we express ISS, in all generality, as nine scalar equations for an incompressible material, written in terms of Ψ , two undetermined scalars p and pτ , and the invariants of F and τ . This form of ISS makes it easier to select appropriate representations for Ψ (F, τ ). For example, we show, with a minor adjustment to the equations, how to use the scalar equations of ISS to deduce an example of Ψ (F, τ ) for a compressible material. It is commonly thought that arteries attempt to maintain a homogeneous stress gradient within their walls [17]. In §5, we calculate this optimal Cauchy stress for a simplified arterial wall model, and then show how by using ISS we can calculate the residual stress that exactly supports this optimal Cauchy stress. We finally compare these results against those of the commonly used opening-angle method [6]. Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 F° 4 ~ ................................................... rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 B° B t s F ~ F B ◦ Figure 1. The coniguration B is the current coniguration with internal stress σ , while the coniguration B is a reference coniguration with internal stress τ . The virtual stress-free state B̃ is a collection of conigurations where the body is stress-free. (Online version in colour.) Figure 2. An intact Wapa tree is under extremely high levels of residual stresses that help to maintain its ideal mechanical conditions. When the Wapa tree is cut, the release of residual stress can cause the tree to burst apart. (Image source: http://www.fao.org.) deformation tensor. Following Shams et al. [15], we take the following complete set of 10 independent invariants: and I1 = tr C, I2 = 12 [I12 − tr(C2 )], I3 = det C, (2.3) Iτ1 = tr τ , Iτ2 = 21 [Iτ21 − tr(τ 2 )], Iτ3 = det τ , (2.4) J1 = tr(τ C), J2 = tr(τ C2 ), J3 = tr(τ 2 C), J4 = tr(τ 2 C2 ). (2.5) Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 The Cauchy stress equation (2.1) can then be written as [15] 5 2 (2.6) where B = FF T is the left Cauchy–Green deformation tensor, and ΨI1 , ΨI2 , ΨJ1 , ΨJ2 , ΨJ3 , ΨJ4 are the partial derivatives of Ψ with respect to I1 , I2 , J1 , J2 , J3 , J4 , respectively. There are no partial derivatives of Ψ with respect to Iτ1 , Iτ2 and Iτ3 appearing in (2.6) because τ does not depend on F. The initial stress compatibility equation (2.2) becomes 2 ◦ ∂Ψ ∂Ψ +4 − p = 0, ∂I1 ∂I2 2 ∂Ψ ∂Ψ +4 =1 ∂J1 ∂J2 and ∂Ψ ∂Ψ +2 = 0. ∂J3 ∂J4 (2.7) Another important physical restriction that can be imposed is the strong-ellipticity condition, which is satisfied when the fourth-order tensor of instantaneous moduli A0piqj = J−1 Fpα Fqβ ∂ 2Ψ ∂Fiα ∂Fjβ satisfies A0piqj np nq mi mj > 0 for every n, m ∈ R3 , (2.8) for compressible materials, while for incompressible materials the above need only hold when n · m = 0. Imposing strong ellipticity implies that plane waves may propagate in every direction with a real valued speed [20], and other physically expected behaviour [21]. For a representation of A0piqj in terms of the invariants of τ and F, see [15]. The issue we address now is how to write Ψ explicitly in terms of the invariants (2.3)– (2.5)? We advocate three criteria. First, the free energy density Ψ should satisfy the initial stress compatibility (2.7). Second, it should satisfy strong ellipticity (2.8) for all deformations in which the material is expected to be stable [22], and we call the third criterion the ISS requirement. (a) Initial stress symmetry For convenience, let the response function ς̂ be denoted by ς̂(F, τ , p) := F ∂Ψ (F, τ ) − pI, ∂F (2.9) for every F and τ within the elastic domain of the material. The scalar p is arbitrary if the material is incompressible, and is p = −2I3 ΨI3 if the material is compressible. The ISS states that ς̂ has no preferred reference configuration. Referring to figure 1, if we take ◦ B as the reference configuration, then the Cauchy stress becomes σ = ς̂ (F, τ , p). (2.10) ◦ However, we can also take B as the reference configuration and B as the current configuration and therefore express the initial stress as τ = ς̂(F −1 , σ , pτ ), (2.11) for some scalar pτ . For a compressible material pτ = −2I3 ΨI3 , where F is replaced with F −1 τ . In a more precise form, ISS can be stated as σ = ς̂ (F, τ , p) = τT and τ = ς̂ (F −1 , σ , pτ ), (2.12) for every F and τ , such that τ and det F = 1 for incompressible materials, and p and pτ are, respectively, given by p = p̂(σ , F, τ ) and pτ = p̂(τ , F −1 , σ ) for some scalar function p̂. The boundary conditions can determine p through its dependence on σ , and analogously for pτ . The ISS agrees with the initial stress compatibility, i.e. the condition (2.2), when we have ς̂ (I, τ , p) = τ for every τ . ................................................... + 2ΨJ1 Fτ F T + 2ΨJ2 F(τ C + Cτ )F T + 2ΨJ3 Fτ 2 F T + 2ΨJ4 F(τ 2 C + Cτ 2 )F T , rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 σ = 2ΨI1 B + 2ΨI2 (I1 B − B ) − pI Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 Referring to figure 1, the virtual stress-free configuration guarantees that, for any τ and F, there ◦ exists an F such that σ = ϑ̂(B̃, p) ◦ and ◦ ◦ τ = ϑ̂(B, p), (2.13) ◦ where B = F −1 B̃F −T ; p and p are scalar fields where p depends on the boundary conditions on B ◦ ◦ ◦ and B̃, while p depends on the boundary conditions on B and B; and ϑ̂(FF T , p) := F ∂Ψ ∂(F T F) F T − pI. Note that by setting τ = 0 in equation (2.6), we see that the right side depends only on FF T and p. ◦ Johnson & Hoger [5] demonstrated that it is possible to determine B from any given τ , together ◦ with appropriate boundary conditions on B, in order to reach useful constitutive equations such as ◦ ϑ̂(F BF T , p) = ς̂(F, τ , p) ◦ for any F and B. (2.14) ◦ As equation (2.13)2 is valid for any F and B, we can substitute F for F −1 ◦ and B for B̃ (2.15) ◦ ◦ and swap the boundary conditions on B and B (these substitutions effectively swap B for B), so that equation (2.13)2 becomes σ = ϑ̂(B̃, p), where we have assumed that the above together with the boundary conditions on B has σ as the unique solution. Analogously, equation (2.13)1 ◦ ◦ becomes τ = ϑ̂(B, p). This means that, when we make the substitutions (2.15) in equation (2.14), we should also swap τ and σ , so that ◦ ◦ ◦ ϑ̂(B, p) = ς̂ (F −1 , σ , p). (2.16) The left-hand side is simply τ and the right-hand side is τ given by ISS (2.12). Finally, the above ◦ condition is identically true, since this result is valid for every F and B. Moreover, we assumed ◦ that B can be determined from τ , so equation (2.16) holds for every F and τ . ................................................... (b) Stress-free coniguration implies initial stress symmetry 6 rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 Another way to view the ISS requirement is to assume that τ and σ are due to the elastic deformation of a virtual stress-free state, which we demonstrate in §2b. Hoger emphasized many times that using the virtual stress-free state does not restrict how the residual stress was formed. The same can be said about ISS, as only F corresponds to an elastic deformation. Moreover, the ISS is not restricted to elasticity and should hold, under suitable modifications, for other constitutive equations such as those encountered in viscoelasticity. One practical outcome from the ISS is to restrict the possible constitutive choices for ς̂. For instance, in §4, we show that the choice Ψ = 12 µ(I1 − 3) + 21 (J1 − Iτ1 ) proposed in [7] does not satisfy ISS. A second practical feature arising from ISS is that we can write the initial stress as a function of the Cauchy stress (2.11). This proves useful when σ is known a priori, as will be discussed in §5, where we determine σ from the principle of homeostasis and then use equation (2.11) to derive τ . The alternative of choosing the initial stress τ first can be far more complicated. Before developing further the implications of ISS in §4, we introduce in the next section an example of stored energy Ψ (F, τ ) that satisfies ISS, stress compatibility (2.2) and strong ellipticity (2.8). Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 3. Initially stressed neo-Hookean material 7 µ (tr C̃ − 3), 2 (3.1) where µ > 0 is the constant shear modulus and C̃ = F̃ T F̃. The Cauchy stress (2.1) then becomes σ = µF̃ F̃ T − pI. (3.2) ◦ We can rewrite tr C̃ by substituting F̃ = F F and using the properties of the trace ◦ ◦ ◦ ◦ ◦ tr C̃ = tr(F̃ T F̃) = tr(F T F T F F) = tr(F F T F T F) = tr(BC), ◦ ◦ ◦ where B = F F T and C = F T F, which leads to Ψ= ◦ µ [tr(BC) − 3]. 2 (3.3) ◦ ◦ We can write B in terms of τ by evaluating equation (3.2) at F = I, σ = τ , p = p and rearranging as follows: ◦ ◦ µB = τ + pI. (3.4) ◦ ◦ ◦ To write p as a function of τ , we require that the deformation F be isochoric, so that det(µB) = ◦ ◦ det(τ + pI) = µ3 , which results in a cubic equation for p: ◦3 ◦2 ◦ p + p Iτ1 + pIτ2 + Iτ3 − µ3 = 0. (3.5) ◦ The real roots for p are  ⎧  1 T1 3/2 ⎪ ⎪ , T2 ≤ T1 , T + − I ⎪ τ 3 1 ⎪ 3 T3 ⎪ ⎪ ⎪  ⎪  ⎨1 ◦ T1 3/2 c1 T3 + c∗1 p= − Iτ1 , −T1 ≤ T2 , ⎪ T3 ⎪3 ⎪ ⎪  ⎪ ⎪1  T1 ⎪ 3/2 3/2 ∗ ⎪ ⎩ c1 T3 + c1 − Iτ1 , −T1 ≤ T2 ≤ T1 , 3 T3 where T1 = Iτ21 − 3Iτ2 ,  T3 = ◦ 3 T2 = Iτ31 − 29 Iτ1 Iτ2 + T22 − T13 − T2 − µ3 ), √ 3 1 i, and c1 = − + 2 2 27 2 (Iτ3 (3.6) (3.7) (3.8) while if T1 = 0, then p = µ − Iτ1 /3. In terms of the eigenvalues τ1 , τ2 and τ3 of τ , we can write ⎫ 3 1 ⎬ T2 = − 27 2 µ + 2 [(τ1 − τ2 ) − (τ3 − τ1 )][(τ2 − τ3 ) − (τ1 − τ2 )][(τ3 − τ1 ) − (τ2 − τ3 )] (3.9) ⎭ and T1 = 21 (τ1 − τ2 )2 + 21 (τ1 − τ3 )2 + 21 (τ2 − τ3 )2 , ................................................... Ψ= rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 Here, we derive a simple constitutive equation for an initially stressed body by assuming that both τ and σ arise from deforming an incompressible neo-Hookean material. The result will be an explicit representation of Ψ in terms of F and τ , given by equation (3.11) below, that automatically satisfies ISS, the stress compatibility (2.7) and strong ellipticity (2.8). Both ISS and stress compatibility hold because τ arises due to an elastic deformation from a stress-free state (see §2b), and strong ellipticity is satisfied because the material is a neo-Hookean solid [23]. Referring to figure 1, as the material is stress-free in B̃, we have Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 3/2 so that, when T1 = 0, we have that τ1 = τ2 = τ3 . From the above, we see that T1 ≥ 0 and T1 µ3 27/2| for any1 else if T2 ≥ 0 then T2 < |T2 + ◦ µ3 27/2| ≤ T1 . Therefore T2 < T1 , and so the condition for the first case for p in equation (3.6) is ◦ always satisfied. We discard the second and third case for p because we expect Ψ , and therefore ◦ ◦ p, to be continuous for every τ ∈ R3 . When τ1 = τ2 = τ3 , the only viable solution for p is (3.6)1 . If τ ◦ moves into a region where (3.6)2 or (3.6)3 becomes real, then p cannot change from (3.6)1 to (3.6)2 or (3.6)3 because it can be shown that (3.6)1 does not equal (3.6)2 or (3.6)3 for any τ ∈ R3 . ◦ To represent tr(BC) in terms of the invariants of τ and C, we multiply each side of equation (3.4) on the right with C and take the trace to get ◦ ◦ µtr(BC) = tr(τ C) + ptr C, (3.10) which we use to write the free energy density equation (3.3) as ◦ Ψ = 12 (pI1 + J1 − 3µ), (3.11) ◦ ◦ with p given by equation (3.6)1 . Note that in the absence of initial stress τ = 0, p = µ by (3.5), J1 = I1 , and then Ψ reduces to the classical neo-Hookean model, as expected. Equation (3.11) represents the general extension of the neo-Hookean strain energy function to an initially stressed material, resulting in a function of only five of the nine independent invariants of C and τ . Substituting (3.11) in the Cauchy stress equation (2.1), we arrive at the constitutive relation ◦ σ = pB + Fτ F T − pI. (3.12) ◦ Taking B as the reference configuration, B as the current configuration and leaving B̃ as the stress-free configuration (figure 1), the initial stress becomes σ , the Cauchy stress becomes τ and equation (3.12) becomes ◦ τ = pτ C−1 + F −1 σ F −T − pτ I, (3.13) ◦ where pτ is given by replacing τ for σ in equation (3.6)1 , and pτ is an undetermined scalar. Substituting (3.12) in (3.13), we find the connection ◦ ◦ (pτ − p)C−1 = (pτ − p)I. (3.14) As this equation must hold for every C, we conclude that ◦ ◦ pτ = p and pτ = p. (3.15) Note that, as is expected of a neo-Hookean material, the above equations do not determine p in terms of τ or F. (a) Plane strain The free energy density (3.11) is simplified when the initial strain has only planar components. ◦ ◦ ◦ Accordingly, let us assume B13 = B23 = 0 and B33 = 1, which substituted in equation (3.4) results ◦ ◦ ◦ in τ13 = τ23 = 0 and τ33 = µ − p. We let BP , CP , BP , IP , σ P and τ P be B, C, B, I, σ and τ restricted to the (x1 , x2 ) plane, respectively. From equation (3.4), we get ◦ ◦ µBP = τ P + pIP . 1 (3.16) See http://math.stackexchange.com/questions/1255883/prove-that-a2b2c2-geq-2a-b2b-c2a-c21-3/1255907#1255907. 8 ................................................... 3/2 ≥ |T2 + rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 3/2 τ1 , τ2 , τ3 ∈ R. So if T2 < 0, then clearly 3/2 T2 < T 1 , Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 ◦ To obtain p in terms of τ P , we take the determinant of each side of equation (3.16), giving ◦2 ◦ so that µ2 = p + ptr τ P + det τ P , ◦ (3.17) ◦ where we used det(µB) = µ2 . We solve the above for p to get ◦ p± = − 21 trτ ± 1 2 −4 det τ + (trτ )2 + 4µ2 , (3.18) where −4 det τ + (trτ )2 + 4µ2 = (τ1 − τ2 )2 + 4µ2 in terms of the eigenvalues of τ . The solution ◦ ◦ ◦ ◦ p = p− should be discarded because p− = −µ when τ P = 0, while equation (3.16) shows that p ◦ ◦ should be positive when τ P = 0. So for τ P = 0, the only viable solution is p = p+ . Further, as we ◦ ◦ ◦ expect p to be continuous in τ P for every τ P ∈ R2 , and p− = p+ for every τ P , we should completely ◦ p− . discard the solution For simplicity, we assume C32 = C31 = 0; then equation (3.3) reads ◦ Ψ = 21 µ(tr(BP CP ) − 2) + 12 µ(C33 − 1), (3.19) ◦ and substituting BP from equation (3.16), we get ◦ Ψ = 21 tr(τ P CP ) + 21 ptr CP − µ + 12 µ(C33 − 1), ◦ (3.20) ◦ where p = p+ given by equation (3.18). The Cauchy stress equation (2.1) becomes ◦ σ P = pBP − pIP + F P τ P F TP , ◦ (3.21) ◦ and σ33 = p − p + τ33 = µ − p, where we have used that τ33 = µ − p. For F = I, we have σ = τ and ◦ p = p. 4. The initial stress symmetry equations In this section, we show how to express the ISS condition (2.12) as a set of scalar equations that relates the free energy density Ψ ; p and pτ from equation (2.12); and the invariants of τ and C. Let us first consider a simple example, namely assuming Ψ = I1 J1 /2. One can check that ◦ this strain energy satisfies the compatibility equations (2.2) with p = tr τ , but does it satisfy the ISS (2.12)? We can use the stress in the form (2.1) to write σ = ς̂ (F, τ , p) = tr(τ C)B − pI + tr(C)Fτ F T , (4.1) and then from ISS we have that τ = ς̂(F −1 , σ , pτ ) = tr(σ B−1 )C−1 − pτ I + tr(B−1 )F −1 σ F −T . (4.2) To check that both equations (4.1) and (4.2) hold for every τ and F, we substitute σ into equation (4.2), which after some rearranging becomes (1 − I1 I2 )τ = (3J1 − 2pI2 + I1 Iτ 1 )C−1 + (I2 J1 − pτ )I, (4.3) where we have used tr(B−1 ) = I2 , which can be shown by applying the Cayley–Hamilton theorem to B with I3 = 1 (due to incompressibility). In order to satisfy equation (4.3) for every F and τ , the coefficients of τ , C−1 and I must all be identically zero (this is shown rigorously in the electronic supplementary material). For the coefficient of τ to be zero, the identity I1 I2 = 1 must hold for every F, which is obviously impossible: so we conclude that Ψ = I1 J1 /2 does not correctly furnish the Cauchy stress for every reference configuration. ................................................... ◦ rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 ◦ det(µBP ) = det(τ P + pIP ), 9 Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 In the electronic supplementary material, we reduce the ISS to nine scalar equations, following a procedure similar to the one above. They are compactly written as (4.4) and {1} {2} {3} {4} ΨJn (Q{n} ΨJσ1 + Q{n} ΨJσ2 + Q{n} ΨJσ3 + Q{n} ΨJσ4 ) = 0, (4.5) where n sums over 0, 1, 2, 3, 4, ΨJ0 := 1, ΨIk = ΨIk (F, τ ), ΨJm = ΨJm (F, τ ), (4.6) ΨJσm := ΨJm (F −1 , σ ), (4.7) and ΨIσk := ΨIk (F −1 , σ ), for k ∈ {1, 2} and m ∈ {1, 2, 3, 4} with σ given by (2.6), ⎤ 2ΨIσ2 ⎢−2Ψ σ ⎥ ⎢ I1 ⎥ ⎥ ⎢ ⎢ pτ ⎥ b=⎢ ⎥, ⎢ 0 ⎥ ⎥ ⎢ ⎣ 0 ⎦ 1 ⎡ ⎤ ΨI2 ⎥ ⎢ p/2 ⎥ ⎢ ⎥ ⎢ −ΨI1 ⎥ ⎢ {1} P =4⎢ ⎥, ⎥ ⎢ −2ΨJ2 ⎥ ⎢ ⎣ ⎦ 2I2 ΨJ2 −ΨJ1 − 2I1 ΨJ2 ⎡ ⎡ p ⎤ ⎢I p − 2Ψ − 2I Ψ ⎥ 1 I2 ⎥ I1 ⎢2 ⎥ ⎢ 2ΨI2 ⎥ ⎢ {2} P =4⎢ ⎥ ⎥ ⎢ 0 ⎥ ⎢ ⎦ ⎣ −2ΨJ1 −4ΨJ2 (4.8) and {1} {1} {2} {2} Q{1} = Q{2} = Q{1} = Q{2} = 0, (4.9) while all other matrices are given explicitly in the electronic supplementary material. Equations (4.4) and (4.5) represent six and three scalar equations, respectively, for the unknowns Ψ , p and pτ . We now look at special cases of Ψ which simplify the ISS equations and lead to more practical representations for the free energy density. (a) Free energy independent of J3 and J4 When Ψ is independent of J3 and J4 , equation (4.5) is identically zero and only the first three terms of equation (4.4) are non-zero. The fourth scalar equation of (4.4) becomes −8ΨJ2 ΨJσ1 = 0, which is true, for every F and τ , only when Ψ does not depend on J2 or does not depend on J1 . We investigate the first case in more detail below. (i) Assuming Ψ independent of J2 , J3 and J4 If Ψ does not depend on J2 , J3 and J4 , then equation (4.4) reduces to 4ΨJ1 ΨJσ1 = 1, ΨIσ2 ΨJσ1 ΨI +  2 = 0, ΨJ1 pΨJσ1 = ΨIσ1 and pτ ΨJ1 = ΨI1 , (4.10) where we have used equation (4.10)1 to derive the other three equations. We can check if these equations are consistent with the compatibility equations (2.7) by letting F = I, σ = τ and pτ = ◦ ◦ p = p. This results in ΨJ1 = ± 12 , 2ΨI1 = ±p and ΨI2 (1 ± 1) = 0, which only satisfies equations (2.7) if ◦ ΨJ1 = 12 , 2ΨI1 = p and ΨI2 = 0 for F = I. Let us apply equations (4.10) to the following example: Ψ = 21 µ(I1 − 3) + 21 (J1 − Iτ1 ), (4.11) ................................................... {4} rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 {3} b + P{1} ΨJσ1 + P{2} ΨJσ2 + ΨJn (P{n} ΨJσ3 + P{n} ΨJσ4 ) = 0 10 Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 σ = µF T F − pI + Fτ F T . (4.12) Let us consider a cube in the reference configuration, subject to the initial stress τij = τi δij , aligned with the axes of the cube. For the current configuration, we impose two clamped conditions that fix F11 = λ1 and F33 = λ3 , then F22 = λ2 = (λ1 λ3 )−1 due to incompressibility, with all the other components of F being zero. With this imposed deformation, the Cauchy stress is diagonal with components σij = δij σi . We can also prescribe the stress σ2 , as this will not alter λ1 and λ3 because they are kept fixed. We choose σ2 = τ2 , to ensure compatibility with the reference configuration when λ1 = λ3 = 1, which results in τ2 = µλ22 − p + λ22 τ2 . Clearly, it is not possible for p to satisfy the above and be fixed at p = µ so that the ISS is also satisfied. This means that the free energy density equation (4.11) either results in non-physical behaviour or does not satisfy the ISS, which implies that the constitutive relation (4.12) does not hold for every reference configuration. The choice of Ψ should satisfy equations (4.10) without restricting the deformation F, initial stress τ or p, otherwise the material is likely to exhibit non-physical behaviour. In this respect, the free energy densities W = 21 µ(I1 − 3) + 21 µ̄(J1 − tr τ )2 + 12 (J1 − tr τ )2 , equation (5.40) in [15], and W = 12 µ(I1 − 3) + 14 (J2 − trτ )2 , equation (78) in [7], where µ and µ̄ are constants, are also likely to exhibit non-physical behaviour. Below we present a free energy density Ψ for a compressible material, inspired by the neo-Hookean example in §3, that satisfies ISS without any unphysical restrictions. For a compressible material, we substitute p with −2I3 ΨI3 and pτ with −2I3−1 ΨIσ3 in the ISS equations (4.10). We will now establish under what conditions does the free energy density −1/3 Ψ = 12 g(τ , I3 )I1 I3 −1/3 + 21 f (τ , I3 ) + 12 J1 I3 (4.13) satisfy ISS (4.10) and the stress compatibility (2.7), where f and g are arbitrary scalar functions. The corresponding Cauchy stress (2.6) is −1/3 σ = g(τ , I3 )I3 −1/3 B + 2I3 ΨI3 I + I3 Fτ F T , (4.14) [g(τ , I3 )I1 + J1 ]. (4.15) with −1/3 ΨI3 = 12 fI3 (τ , I3 ) + 21 gI3 (τ , I3 )I1 I3 −4/3 − 61 I3 −1/3 1/3 The ISS equation (4.10)1 is satisfied as ΨJ1 ΨJσ1 = ( 21 )I3 ( 21 )I3 = 41 and stress compatibility (2.7)2 is satisfied as 2ΨJ1 = 1 when F = I. By substituting p for −2I3 ΨI3 , the remaining compatibility equation (2.7)1 becomes g(τ , I3 ) = −2ΨI3 evaluated at F = I. (4.16) Equation (4.16) is satisfied if we set fI3 (τ , 1) = trτ /3 − 3gI3 (τ , 1); then equation (4.16) and equation (4.14) together imply that σ = τ when F = I. The ISS equations (4.10)3 and (4.10)4 now read, respectively, − 2I3 ΨI3 = g(σ , I3−1 ) and − 2I3−1 ΨIσ3 = g(τ , I3 ), (4.17) ................................................... Example 4.1. For the free energy density expressed in equation (4.11), the Cauchy stress from equation (2.6) becomes 11 rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 used by Merodio et al. [7]. In this case, equations (4.10)1 and (4.10)2 are satisfied while equations (4.10)3 and (4.10)4 become p = pτ = µ. However, restricting p to be constant, when it is normally determined from boundary and/or initial conditions, would result in physically unexpected behaviour. This is best seen with an example. Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 −1/3 Fτ F T + g(τ , I3 )I3 B) = I3−1 det B det(τ + g(τ , I3 )I), and by using incompressibility det B = 1 we find that −1/2 1/2 I3 det(σ − 2I3 ΨI3 I) = I3 det(τ + g(τ , I3 )I). (4.18) If we choose g(τ , I3 ) such that −1/2 I3 det(τ + g(τ , I3 )I) = k, (4.19) where k is a constant, then the right-hand side of equation (4.18) is also equal to k, thus σ must satisfy 1/2 I3 det(σ − 2I3 ΨI3 I) = k. (4.20) By swapping τ for σ and F for F −1 , equation (4.19) becomes 1/2 I3 det(σ + g(σ , I3−1 )I) = k. (4.21) For g(τ , I3 ) to satisfy equation (4.19), there are three possible solutions given by replacing µ3 with 1/2 ◦ kI3 and g(τ , I3 ) with p in (3.6)1 , (3.6)2 and (3.6)3 , which we, respectively, denote by g1 (τ , I3 ), g2 (τ , I3 ) and g3 (τ , I3 ). We highlighted in §3 that only g(τ , I3 ) = g1 (τ , I3 ) is a physically viable solution. So if we choose g(τ , I3 ) = g1 (τ , I3 ), then for equation (4.17)1 to be satisfied, we need g1 (σ , I3−1 ) = −2I3 ΨI3 for every F and τ . (4.22) We will show that the above is a consequence of equations (4.20), (4.21) and stress compatibility (4.16). Due to equation (4.16), we know that, initially for F = I, equation (4.22) is true. We can also conclude from equations (4.20) and (4.21) that −2I3 ΨI3 will equal either g1 (σ , I3−1 ), g2 (σ , I3−1 ) or g3 (σ , I3−1 ) for every F. Since −2I3 ΨI3 should be a continuous function of F and τ , it cannot change from g1 (σ , I3−1 ) to another solution gk (σ , I3−1 ) because it can be shown that gk (σ , I3−1 ) = g1 (σ , I3−1 ) for every σ and I3 . For this reason, if g1 (τ , 1) = −2ΨI3 for F = I, then g1 (σ , I3−1 ) = −2I3 ΨI3 for every F. In conclusion, the model (4.13) satisfies ISS as long as g(τ , I3 ) = g1 (τ , I3 ), where g1 (τ , I3 ) is given 1/2 by replacing µ3 with kI3 in (3.6)1 . 5. Homogeneous stress gradient in a hollow cylinder It is now well acknowledged that residual stresses in living materials are vital to maintain ideal mechanical conditions. When an external load changes, growth and remodelling alter the residual stress to best adapt to the new load. An excellent example is how arteries remodel in response to the internal pressure. The residual stress in arteries is thought to protect the arterial wall against strain concentration [24] or stress concentration [1,17]. Here, we adopt the most accepted hypothesis—that the residual stress in the artery acts to minimize the stress gradient [1,25]. The reasoning behind this hypothesis is that if a given tissue grows in response to stress, then homeostasis is only possible if the tissue is under similar stress conditions throughout. So by minimizing the stress gradient, we are selecting the most homogeneous stress possible. We will also consider a simplified artery with only one layer and no shear stress applied to the interior of the artery. ................................................... −1/3 det(σ − 2I3 ΨI3 I) = det(I3 12 rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 for every F and τ , where fI3 and gI3 are, respectively, the partial derivatives of f and g with respect to I3 . Note that the above two equations are equivalent as equation (4.17)1 becomes equation (4.17)2 when we substitute F for F −1 and τ for σ . For this reason, it is sufficient to satisfy equation (4.17)1 . If F = I, then equation (4.17)1 is satisfied as it reduces to equation (4.16). The rest of this section will focus on showing that equation (4.17)1 is satisfied for every F and τ . By rearranging equation (4.14) and taking the determinant on either side, we get Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 To describe the arterial wall, we use cylindrical coordinates (r, θ, z) with the components of σ written in terms of unit basis vectors. We assume that the arterial wall retains its cylindrical symmetry when the internal pressure is removed and that there is no shear stress at the inner wall. These assumptions imply that, when a pressure is applied in the cylinder, the resulting deformation is an inflation. Due to radial symmetry, the Cauchy stress σ is independent of θ and z, so that in the absence of body forces the equilibrium equations reduce to ∂ 2 (r σθr ) = 0, ∂r and ∂ (rσzr ) = 0 ∂r ∂ (rσrr ) − σθθ = 0, ∂r (5.1) (5.2) with the boundary conditions σθr = 0, σzr = 0, σrr = −P, for r = a, (5.3) and σθr = 0, σzr = 0, σrr = 0, for r = b, (5.4) where a and b are the inner and outer radius of the loaded artery, respectively. The equilibrium equations (5.1) together with the boundary conditions lead to σθr = 0 and σzr = 0 for all r. As neither σzθ nor σzz appear in the equations of equilibrium, they are only restricted by the constitutive choice and boundary conditions on the cross-section of the artery. We use this degree of freedom to assume there is no axial torsion σzθ = 0. Note that, for the constitutive choice equation (2.6), τZZ can be chosen so that σzz is constant, while for the simpler choice of plane strain equation (3.21) (with σzz = σ33 ), σzz is determined by p. Finally, the remaining equation of equilibrium is equation (5.2). (b) Minimal stress gradient ′ , where the Here, we develop a method to minimize the stress gradient fields σrr′ and σθθ prime denotes differentiation with respect to r. We make no assumptions about any reference configuration nor do we make any constitutive choice. We only make use of the assumptions from the section above which result in both σrr and σθθ being independent of the coordinates θ and z; σθr = σzr = σzθ = 0; and, for simplicity, we will not consider σzz . Once σθθ and σrr are determined, we use ISS equation (2.11) to write the residual stresses τΘΘ and τRR in terms of σθθ , σrr and the deformation gradient in §5c. Using the equilibrium equation (5.2), we write σθθ in terms of σrr as σθθ = σrr + rσrr′ , ′ so that σθθ = 2σrr′ + rσrr′′ . (5.5) As aortas and veins are of many different sizes, for the sake of generality, let us introduce the following dimensionless variables: ξ= r−a b−a and ϕ(ξ ) = σrr (r) , P (5.6) ................................................... (a) Plane strain cylinder 13 rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 Our strategy is to first choose an optimal Cauchy stress field σ , and then to derive the residual stress from the ISS equation τ = ς̂ (F −1 , σ , pτ ). So first we calculate the Cauchy stress with near homogeneous circumferential and radial components in §5b, then we find the residual stress that supports this optimal Cauchy stress in §5c. In §5b, we also show that the optimal Cauchy stress has a simple asymptotic formula. We finally compare the results with the corresponding ones obtained using the opening angle method [6]. Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 which result in ′ σθθ = σrr′′ = ϕ ′′ 14 P (b − a)2 (5.7) P (2ϕ ′ + (ξ + α)ϕ ′′ ). b−a (5.8) Our aim is to minimize the stress gradient density 1 b P2 1 ′ 2 [(σrr′ )2 + (σθθ ) ] dr = 2 α 2 [(ϕ ′ )2 + (2ϕ ′ + (ξ + α)ϕ ′′ )2 ] dξ , b−a a a 0 with the constraint ϕ(1) = 0 and ϕ(0) = −1, giving 1 0 (5.9) ϕ ′ dξ = 1, (5.10) where we have introduced the quantity α ≡ a/(b − a). Using calculus of variations [26], we find that ϕ ′ must satisfy the Euler–Lagrange equation d ∂f ∂f − = Λ for ξ ∈ (0, 1) ∂ϕ ′ dξ ∂ϕ ′′ and ∂f = 0 for ξ = 0, 1, ∂ϕ ′′ (5.11) where f = (ϕ ′ )2 + (2ϕ ′ + (ξ + α)ϕ ′′ )2 , (5.12) and Λ is a Lagrange multiplier due to the constraint (5.10). The solution to (5.11) is given by α Λ Λ √ ϕ = − ( 13 − 3) 6 6 ′ + 2Λ √ 3( 13 − 3) α √ 13/2+1/2 √ α 13 √ 13/2+1/2 √ (α − (1 + α) 13 − (1 + α) √ √ 13/2−1/2 − (1 + α) 13/2−1(1/2) √ √ α 13 − (1 + α) 13  + ξ) √ 13/2+1/2 α(1 + α) α+ξ √13/2+1/2 We determine Λ from the constraint (5.10) to be √ √  √ 1 √ α 13/2+1/2 − (1 + α) 13/2+1/2 √ √ 1 − ( 13 − 3) (α + ξ ) 13/2+1/2 Λ=6 0 α 13 − (1 + α) 13 +√ 4 13 − 3 α √ √ 13/2−1/2 − (1 + α) 13/2−1(1/2) √ √ α 13 − (1 + α) 13  α(1 + α) α+ξ √13/2+1/2 dξ −1 . (5.13) . (5.14) Realistic values for α = a/(b − a) can be obtained from in vivo measurements of the lumen radius (inner radius) a and total artery thickness b − a. For the descending thoracic aorta of healthy men around 51 years old, the mean lumen radius is 20 mm [27] and the mean thickness is 1.4 mm [28]. These estimates give α −1 ≈ 0.07, making it appropriate to expand ϕ ′ and Λ as a power series in α −1 , which gives   Λ 1 1 + (1 − 2ξ )α −1 + (9ξ 2 − 6ξ − 1)α −2 + O(α −3 ) (5.15) ϕ′ = 2 3 and Λ = 2 + 23 α −2 + O(α −3 ), (5.16) ϕ ′ = 1 + (1 − 2ξ )α −1 + ξ (3ξ − 2)α −2 + O(α −3 ). (5.17) or simply We then substitute ϕ ′ in equations (5.5), (5.6) and (5.7) to obtain σrr = −P(1 − ξ )(1 − ξ α −1 + ξ 2 α −2 ) + O(α −3 ) (5.18) σθθ = Pα(1 + α −3 ξ 2 (4ξ − 3)) + O(α −3 ). (5.19) and ................................................... and P , b−a rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 σrr′ = ϕ ′ Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 (a) (b) 20 s rr /P – x sqq /P –0.980 16 14 –0.985 –0.990 –0.995 12 0 0.2 0.4 0.6 0.8 a –1 1.0 0.2 0.4 x 0.6 0.8 1.0 x Figure 3. The plots of the dimensionless Cauchy stress against the dimensionless radius ξ , where the colour shades from blue to red as α −1 (the wall thickness divided by the inner radius) goes from 0.05 to 0.1. (Online version in colour.) (a) (b) 14.34 P/µ 14.30 14.28 s rr /P – x –0.975 14.32 sqq /P –0.970 –0.980 –0.985 –0.990 –0.995 14.26 0 0.2 0.4 0.6 x 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 x Figure 4. Graphs of the dimensionless Cauchy stress against the dimensionless radius ξ for the opening angle method (dashed curves) and the optimal stress equations (5.18) and (5.19) (red solid curves), for which we used parameters for the descending thoracic aorta: a = 20.0 mm, b = 21.4 mm, α −1 = 0.07 and λz = 1. As the dashed curves shade from yellow to green (arrow), P/µ goes through the values 0.05, 0.20 and 0.35 with optimal opening angle φ = 65.9◦ , 204.2◦ and 261.7◦ ; stress-free reference inner radius 20.3, 25.2 and 31.1 mm; and stress-free reference outer radius 22.0, 27.8 and 34.4 mm, respectively. (Online version in colour.) We plot the above Cauchy stress for α −1 from 0.05 to 0.10, which is within the physiological range, in figure 3. A graph for σrr /P would show all the curves bunched along the same straight line, so instead we have depicted the curves for σrr /P − ξ . Figure 3 reveals that, as the relative thickness of the arterial wall increases (shading from blue towards red), the circumferential stress σθθ decreases while the radial stress σrr increases and becomes less homogeneous. A well-established method to quantify the residual stresses within arteries is the opening angle method [6]. For a neo-Hookean material (3.1), we will use the opening angle method to minimize the circumferential stress component  b  dσθθ 2 1 dr, (5.20) b−a a dr in terms of the opening angle φ, restricted to the boundary conditions (5.3) and (5.4). The only two parameters with a unit of time and mass are P and µ, so by rewriting P = P0 µ the stress will not depend on µ. The results for P0 = 0.05, P0 = 0.20 and P0 = 0.35 are compared with the optimal stress equations (5.18) and (5.19) with parameters for the descending thoracic aorta in figure 4. We see that σθθ for the opening angle method converges to the optimal stress σθθ as P/µ tends to zero, while the plots for σrr all overlap. Note however that σzz for the opening angle method is not necessarily homogeneous in r. ................................................... 18 15 –0.975 rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 a –1 Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 Here, we use finite elasticity to connect the current state with the unloaded state. Since the cylindrical symmetry of the artery is maintained when the internal pressure is removed, we describe the unloaded state with the cylindrical coordinates (R, Θ, Z). Then the deformation gradient for the unit basis vectors er , eθ , ez , ER , EΘ and EZ becomes F= r ∂r er ⊗ ER + eθ ⊗ EΘ + λz ez ⊗ EZ , ∂R R (5.21) where we let λz be a constant. Assuming the material is incompressible, we get det F = 1, or ∂r r λz = 1. ∂R R (5.22) Let the reference configuration be a hollow cylinder with inner and outer radius A and B, respectively, so that R ∂r 2 and = λ−1 , (5.23) r = (R2 − A2 )λ−1 z +a ∂R r z which we will use to replace ∂r/∂R wherever it appears. We will assume that the residual stress τ is homogeneous in Θ and Z, as the deformation, the Cauchy stress from equations (5.18) and (5.19), and the boundary conditions (5.25) and (5.26) are all homogeneous in Θ and Z. The equilibrium equations now reduce to τΘΘ = (RτRR )′ for A < R < B, (5.24) with the zero-traction boundary conditions τRR = 0 for R = A (5.25) τRR = 0 for R = B. (5.26) and After making a constitutive choice for Ψ , the residual stresses τRR and τΘΘ will be completely determined from F, σ and pτ due to ISS (2.11). (i) Ψ independent of I2 , J2 , J3 and J4 To simplify, we assume Ψ independent of I2 , J2 , J3 and J4 , so the Cauchy stress (2.6) becomes σ = 2ΨI1 FF T − pI + 2ΨJ1 Fτ F T , and, from ISS (2.12), we can swap τ , F and p with σ , F −1 (5.27) and pτ , respectively, to get τ = 2ΨIσ1 F −1 F −T − pτ I + 2ΨJσ1 F −1 σ F −T . (5.28) Substituting F from equation (5.21) into the above equation, we arrive at r2 (Ψ σ + ΨJσ1 σrr ) − pτ R2 I1 (5.29) R2 σ (Ψ + ΨJσ1 σθθ ) − pτ . r2 I1 (5.30) τRR = 2λ2z and τΘΘ = 2 We remark that τZZ has not been derived since we have not taken into consideration σzz in §5b. ................................................... (c) Residual stress 16 rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 We can now invoke finite elasticity to derive the residual stress τ in the unloaded state necessary to sustain the optimal homogeneous Cauchy stress, and then compare our results with the residual stress predicted by the opening angle method. To do so, we let Ψ be a function of both stress and strain and use the ISS, in the form of equation (2.11), to determine τ as a function of σ and F. Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 17 20 0.2 P/m 0.3 0.4 –10 –20 Figure 5. The loaded artery wall is illustrated by the red horizontal lines, while the unloaded artery wall is illustrated by the blue curves. The y-axis gives the position of the loaded and unloaded walls measured from the centre of the artery. The stress in the loaded geometry is taken to be the optimal Cauchy stress (5.18) and (5.19), and we assume a residually stressed neo-Hookean material (3.20). (Online version in colour.) To compare unloaded arteries of different sizes, we write r and R in terms of the dimensionless radius ξ , r = (b − a)ξ + a and R = A2 + λz [(a + (b − a)ξ )2 − a2 ], (5.31) which we use to rewrite the ordinary differential equation (5.24) in the form R   dR −1 dτRR dτRR dξ dτRR + τRR − τΘΘ = R + τRR − τΘΘ = R + τRR − τΘΘ = 0, dR dξ dR dξ dξ (5.32) which implies that dτRR 1 dR + (τRR − τΘΘ ) = 0. dξ R dξ (5.33) Integrating both sides in ξ and using the unloaded boundary condition (5.25) and (5.26), we reach τRR = R 0 τΘΘ − τRR dR dξ = 0 and R dξ 1 0 τΘΘ − τRR dR dξ = 0. R dξ (5.34) Equation (5.34)2 is independent of pτ and can be used to solve for one of the parameters P, a, b, A, λz or µ (note that B depends on the other parameters because equation (5.23) evaluated at R = B gives B2 = A2 + λz (b2 − a2 )). The only two parameters with a unit of time or mass are P and µ, so by rewriting P = P0 µ equations (5.34) will not depend on µ. To illustrate, we adopt the neo-Hookean material model (3.20), which results in ΨJσ1 = 21 and ◦ ◦ ◦ 2ΨIσ1 = p(σ ), where p is given by substituting τ for σ restricted to (r, θ) in p+ of equation (3.18). To calculate the residual stress that supports the optimal Cauchy stress, we substitute σθθ and σrr from equations (5.18) and (5.19) into τΘΘ and τRR in equations (5.29) and (5.30). We then use equation (5.34)2 to determine the unloaded geometry from the loaded geometry. To illustrate, we take the parameters for the descending thoracic aorta, as used in the previous section, a = 20 mm, b = 1.4 mm and λz = 1. We can then determine the unloaded inner radius A for different values of P0 , which is shown in figure 5. Surprisingly, the unloaded geometry given by the opening angle method is approximately the same as shown in figure 5. We can also investigate the residual stress in the unloaded state. Once A is determined from equation (5.34)2 , we can use equation (5.34)1 to determine τRR , and then τΘΘ from equation (5.24). To compare with the opening angle method, we choose P/µ = 0.05, 0.20 and 0.35, which all give a reasonable unloaded geometry (figure 5). The results are shown in figure 6. ................................................... 0.1 rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 10 Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 (b) sqq /P 0.2 –0.5 –1.0 –1.5 0.4 0.6 x 0.8 1.0 –0.02 –0.04 –0.06 –0.08 –0.10 –0.12 –0.14 0.4 P/m 0.6 x 0.8 1.0 Figure 6. The dimensionless residual stress against the dimensionless radius ξ . The dashed curves are from the opening angle method and the solid curves are from equation (5.24) together with the optimal stress equations (5.18) and (5.19). As P/µ goes through the values 0.05, 0.20 and 0.35 (arrows), the dashed curves shade from yellow to green and the solid curves shade from blue to red. The unloaded geometry for both methods is approximately the same as shown in igure 5. The other parameters are a = 20.0 mm, b = 21.4 mm, α −1 = 0.07 and λz = 1. (Online version in colour.) (a) P/m = 0.20 5 (b) 4 0.004 DY/P Y/P P/m = 0.20 0.006 3 2 P/m = 0.05 0.2 0.4 0.6 x 0.8 P/m = 0.05 0.002 0.2 –0.002 0.4 0.6 x 0.8 1.0 1.0 Figure 7. (a) Comparison of the dimensionless strain energy density Ψ/P versus the dimensionless radius ξ resulting from the opening angle method, yellow (green) dashed curves, with the model (3.20) whose residual stress maintains the optimal stress equations (5.18) and (5.19), the red (blue) solid curves. The parameters used are given in igure 6. (b) Ψ/P from the opening angle method minus Ψ/P from the model (3.20), which we denote as Ψ/P. (Online version in colour.) Biological tissues are very energy efficient. It is a common line of reasoning that biological systems adapt so as to minimize their potential energy. So it is possible that biological materials remodel their residual stress to lower their potential energy. Figure 7a shows the free energy density Ψ for the opening angle method and for the model (3.20) with the optimal stress equations (5.18) and (5.19). In all cases, Ψ is approximately a straight line that decreases as ξ increases away from the pressured boundary ξ = 0 (r = a). Figure 7b shows the difference between Ψ/P from the opening angle method minus Ψ/P from the model (3.20), which we denote as Ψ/P. The integral of Ψ/P over ξ ∈ [0, 1] is positive if the opening angle method has on average a larger free energy density than the model (3.20), with the optimal stress equations (5.18) and (5.19). We found that as P/µ increased, so did the total free energy in the cylinder cross-section 1 −4 for P/µ = 0.05, 0 Ψ/P(b − a)(ξ (b − a) + a) dξ ; for instance, this integral evaluates to −1.6410 2.3710−3 for P/µ = 0.20, and 4.2510−3 for P/µ = 0.35. Essentially, we can see from figures 4 and 7 that the optimal stress equations (5.18) and (5.19) with the neo-Hookean material (3.20) produce a more homogeneous stress and a lower free energy than the opening angle method, though the two methods gave very similar results. The advantage in using the ISS (2.12) for a homeostasis hypothesis for the Cauchy stress σ is twofold. First, we have the freedom to choose σ so as to satisfy the homeostasis hypothesis, which with ISS becomes a separate step from choosing a constitutive equation. Second, it can be relatively straightforward to use ISS (2.12) to quantify the residual stress needed to support the chosen σ . ................................................... P/m 18 0.2 rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 1.5 1.0 0.5 s rr /P – x (a) Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 6. Conclusion Authors’ contributions. All the ideas and the mathematics were developed by all three authors. P.C. and M.D. edited and revised the manuscript and the electronic supplementary material. A.L.G. wrote the first drafts, generated the figures and drafted the electronic supplementary material. All authors gave final approval for publication. Competing interests. We have no competing interests. Funding. Partial funding by the Irish Research Council, the Royal Society, the Hardiman Scholarship Programme at the National University of Ireland Galway and Regione Lombardia, through the ‘Start-up Packages and PhD Program’ project, is gratefully acknowledged. References 1. Fung Y. 1991 What are the residual stresses doing in our blood vessels? Ann. Biomed. Eng. 19, 237–249. (doi:10.1007/BF02584301) 2. Holzapfel GA, Ogden RW. 2003 Biomechanics of soft tissue in cardiovascular systems. CISM Courses and Lectures, vol. 441. Vienna, Austria: Springer Science & Business Media. 3. Hoger A. 1986 On the determination of residual stress in an elastic body. J. Elast. 16, 303–324. (doi:10.1007/BF00040818) 4. Hoger A. 1993 The elasticity tensors of a residually stressed material. J. Elast. 31, 219–237. (doi:10.1007/BF00044971) 5. Johnson BE, Hoger A. 1995 The use of a virtual configuration in formulating constitutive equations for residually stressed elastic materials. J. Elast. 41, 177–215. (doi:10.1007/ BF00041874) 6. Chuong C, Fung Y. 1986 Residual stress in arteries. In Frontiers in biomechanics (eds GW Schmid-Schönbein, SL-Y Woo, BW Zweifach), pp. 117–129. New York, NY: Springer. 7. Merodio J, Ogden RW, Rodríguez J. 2013 The influence of residual stress on finite deformation elastic response. Int. J. Nonlinear Mech. 56, 43–49. (doi:10.1016/j.ijnonlinmec.2013.02.010) ................................................... Data accessibility. The material supporting this article has been uploaded as part of the electronic supplementary material. rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 In order to quantify the initial stress within a solid using non-destructive experimental techniques, it is advantageous to write the free energy density in terms of the deformation gradient F and the initial stress τ as Ψ = Ψ (F, τ ). In this article, we presented a new constitutive condition, the initial stress symmetry (ISS), that aids in proposing suitable constitutive relations for Ψ = Ψ (F, τ ) by providing nine scalar equations (equations (4.4) and (4.5)). One immediate result is that guessing a functional dependence for Ψ = Ψ (F, τ ) is not a trivial task. In fact all choices for Ψ (F, τ ) used in the literature so far do not satisfy ISS. Conversely, using ISS, we proposed two simple choices for Ψ = Ψ (F, τ ), one incompressible equation (3.11) and one compressible equation (4.13), which can be viewed as initially stressed neo-Hookean materials. One consequence of ISS is that the initial stress can be derived from the Cauchy stress. Furthermore, ISS suggests that it is possible to first choose the Cauchy stress, and second to make a constitutive choice and then use equation (2.11) to determine the corresponding initial stress. We illustrated this method in §5 for a simplified arterial wall model, where we chose the ideal Cauchy stress by using a minimal stress gradient hypothesis as the homeostatic condition, and then we calculated the optimal residual stresses of the unloaded remodelled artery. It is reasonable that other tissues with smooth muscle cells may also remodel to minimize their stress gradient. This means that the approach we developed for the arterial wall could be used to help test this hypothesis for other tissues, such as those found in the heart. Since initial stresses are widespread in both inert and living matter, the proposed constitutive restriction can be used in many applications towards a non-destructive quantification of initial tensions in solids. Future research includes determining a wider class of material behaviours, e.g. taking into account natural material anisotropy, and linking the initial stress to elastic wave speeds. 19 Downloaded from http://rspa.royalsocietypublishing.org/ on November 13, 2015 20 ................................................... rspa.royalsocietypublishing.org Proc. R. Soc. A 471: 20150448 8. Shams M, Ogden RW. 2014 On Rayleigh-type surface waves in an initially stressed incompressible elastic solid. IMA J. Appl. Math. 79, 360–376. (doi:10.1093/imamat/hxs070) 9. Destrade M, Ogden RW. 2012 On stress-dependent elastic moduli and wave speeds. IMA J. Appl. Math. 78, 965–997. (doi:10.1093/imamat/hxs003) 10. Man CS, Lu W. 1987 Towards an acoustoelastic theory for measurement of residual stress. J. Elast. 17, 159–182. (doi:10.1007/BF00043022) 11. Merodio J, Ogden RW. 2015 Extension, inflation and torsion of a residually stressed circular cylindrical tube. Contin. Mech. Thermodyn. (doi:10.1007/s00161-015-0411-z ) 12. Soldatos KP. 2013 Modelling framework for mass-growth. Mech. Res. Commun. 50, 50–57. (doi:10.1016/j.mechrescom.2013.03.005) 13. Wang H, Luo X, Gao H, Ogden R, Griffith B, Berry C, Wang TJ. 2014 A modified Holzapfel– Ogden law for a residually stressed finite strain model of the human left ventricle in diastole. Biomech. Model. Mechanobiol. 13, 99–113. (doi:10.1007/s10237-013-0488-x) 14. Shams M. 2010 Wave propagation in residually-stressed materials. PhD thesis, University of Glasgow, Glasgow, UK. 15. Shams M, Destrade M, Ogden RW. 2011 Initial stresses in elastic solids: constitutive laws and acoustoelasticity. Wave Motion 48, 552–567. (doi:10.1016/j.wavemoti.2011.04.004) 16. Johnson BE, Hoger A. 1998 The use of strain energy to quantify the effect of residual stress on mechanical behavior. Math. Mech. Solids 3, 447–470. (doi:10.1177/108128659800300405) 17. Taber LA, Eggers DW. 1996 Theoretical study of stress-modulated growth in the aorta. J. Theor. Biol. 180, 343–357. (doi:10.1006/jtbi.1996.0107) 18. Détienne P, Thiel J. 1988 Monographie des wapas de guyane franćaise. Bois Forêts Trop. 216, 43–68. 19. Guillou AK, Ogden RW. 2006 Growth in soft biological tissue and residual stress development. In Mechanics of biological tissue (eds GA Holzapfel, RW Ogden), pp. 47–62. Berlin, Germany: Springer. 20. Truesdell C. 1966 Existence of longitudinal waves. J. Acoust. Soc. Am. 40, 729–730. (doi:10.1121/1.1910144) 21. Walton JR, Wilber JP. 2003 Sufficient conditions for strong ellipticity for a class of anisotropic materials. Int. J. Nonlinear Mech. 38, 441–455. (doi:10.1016/S0020-7462(01)00066-X) 22. Merodio J, Ogden RW. 2003 Instabilities and loss of ellipticity in fiber-reinforced compressible non-linearly elastic solids under plane deformation. Int. J. Solids Struct. 40, 4707–4727. (doi:10.1016/S0020-7683(03)00309-3) 23. Ogden RW. 1997 Non-linear elastic deformations. Mineola, NY: Courier Dover Publications. 24. Destrade M, Liu Y, Murphy JG, Kassab GS. 2012 Uniform transmural strain in prestressed arteries occurs at physiological pressure. J. Theor. Biol. 303, 93–97. (doi:10.1016/j.jtbi. 2012.03.010) 25. Cardamone L, Valentin A, Eberth J, Humphrey J. 2009 Origin of axial prestretch and residual stress in arteries. Biomech. Model. Mechanobiol. 8, 431–446. (doi:10.1007/s10237-008-0146-x) 26. Gregory J, Lin C. 1992 Constrained optimization in the calculus of variations and optimal control theory. London, UK: Springer. 27. Stefanadis C et al. 1995 Pressure-diameter relation of the human aorta: a new method of determination by the application of a special ultrasonic dimension catheter. Circulation 92, 2210–2219. (doi:10.1161/01.CIR.92.8.2210) 28. Mensel B, Quadrat A, Schneider T, Kühn JP, Dörr M, Völzke H, Lieb W, Hegenscheid K, Lorbeer R. 2014 MRI-based determination of reference values of thoracic aortic wall thickness in a general population. Eur. Radiol. 24, 2038–2044. (doi:10.1007/s00330-014-3188-8)