Exactly isochoric deformations of soft solids
eepl draft
Exactly isochoric deformations of soft solids
John S Biggins , Z Wei ,L Mahadevan , Cavendish Laboratory, JJ Thomson Ave, University of Cambridge, Cambridge, CB3 0HE, United Kingdom School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA Department of Physics, Harvard University, Cambridge, MA 02138, USA [email protected]
PACS – Static Elasticity
PACS – Elasticity mechanical properties of solids
PACS – Elasticity theory in biological physics
Abstract –Many materials of contemporary interest, such as gels, biological tissues and elas-tomers, are easily deformed but essentially incompressible. Traditional linear theory of elasticityimplements incompressibility only to first order and thus permits some volume changes, which be-come problematically large even at very small strains. Using a mixed coordinate transformationoriginally due to Gauss, we enforce the constraint of isochoric deformations exactly to develop alinear theory with perfect volume conservation that remains valid until strains become geometri-cally large. We demonstrate the utility of this approach by calculating the response of an infinitesoft isochoric solid to a point force that leads to a nonlinear generalization of the Kelvin solution.Our approach naturally generalizes to a range of problems involving deformations of soft solidsand interfaces in 2 dimensional and axisymmetric geometries, which we exemplify by determiningthe solution to a distributed load that mimics muscular contraction within the bulk of a soft solid.
A basic question in the theory of elasticity, particularlywithin the small-strain limit of linear elasticity, is the re-sponse function of an infinite solid to a point force, theGreen’s function of the solid known as the Kelvin solution[1]. Knowledge of this function allows one to calculatethe response of the solid to an arbitrary force distribu-tion as the weighted sum of point force responses. In thepresence of constraints such as incompressibility that arecommonly used to approximate the mechanical behaviorof many soft solids such as gels, elastomers and biologicaltissues, we can ask how this solution changes. In this let-ter we answer this question by constructing a theory thatis exactly isochoric and use it to construct the isochoricanalog of the Kelvin point force solution that only im-plements the constraint of incompressibility to first order.Although we focus on a point force in an infinite body forease of exposition, our approach is easily generalizable toany 2-D or axisymmetric situation.Soft solids often have bulk moduli several orders of mag-nitudes higher than their shear moduli. If a material issubject to a displacement field u , the local shape changeis described by the deformation-gradient F = I + ∇ u , andthe local volume changes by a factor of Det( F ). The sim-plest energy function for describing such a material is the compressible neo-Hookean model which, in 2-D, is E nl = 12 µ (cid:32) Tr (cid:0) F · F T (cid:1) Det( F ) − (cid:33) + 12 B (Det( F ) − (1)where µ and B are the shear and bulk moduli. In equilib-rium, balancing the gradients of the shear and bulk termswith respect to u implies that µ ∇ u ∼ B (Det( F ) − µ/B (cid:28) E = µ (Tr (cid:0) F · F T (cid:1) − , Det( F ) = 1 . (2)The constraint Det( F ) = 1 is then naturally conjugate toa pressure field P that serves as a Lagrange multiplier.Then, in the presence of an external force distribution f the total effective energy density is given by˜ E = µ Tr (cid:0) F · F T (cid:1) + P (Det( F ) − − f · u . (3)where the last two terms correspond to the work doneby these volumetric forces. Minimizing the total energyp-1 a r X i v : . [ c ond - m a t . s o f t ] J u l ohn S Biggins1, Z Wei2,L Mahadevan2,3 pq xz pq xza) b) q q p(x ,q) p(x ,q) x x z(x,q ) z(x,q ) (a) Representations of Deformations Fig. 1: A two dimensional elastic body labeled by the coordinates ( p, q ) is deformed into the target state with coordinates ( x, z ).Top diagram shows lines of constant p and q in both states, the traditional representation of an elastic deformation. Lowerdiagram shows lines of constant x and q in both states, a representation of the deformation that makes volume preservationeasy to implement. A patch of material defined by x < x < x and q < q < q is drawn in both states. (cid:82) ˜ E d x with respect to u and P yields the Euler-Lagrangeequations ∇ · ( µF − P F − T ) = − f Det( F ) = 1 , (4)where we identify µF − P F − T as the first Piola-Kirchhoff(PK1) stress-tensor. The nonlinear constraint of volumepreservation makes the equations rarely tractable analyti-cally, unless various symmetries are imposed. However, inthe small strain limit when ∇ u and ∇ P are small, we canlinearize both equations, so that in the limit |∇ u | (cid:28)
1, weget µ ∇ u − ∇ P = − f , ∇ · u = 0 , (5)which are mathematically equivalent to Stokes equationsfor creeping viscous flow, modulo the interpretation of u as the velocity field, and µ as the dynamic viscosity. In atwo-dimensional Cartesian coordinate system p − q , setting u = ˆ q ∂ p α − ˆ p ∂ q α , where α ( p, q ) is the stream functionautomatically guarantees linearized incompressibility ∇ · u = 0. Then Eqn. (5) reads as µ (cid:18) − α qqq − α ppq α qqp + α ppp (cid:19) − (cid:18) P p P q (cid:19) = − f , (6)where subscripts denotes differentiation, i.e. ( · ) x = ∂ ( · ) /∂x etc. In a region with f = 0 eliminating P revealsthat α is biharmonic so that α qqqq + 2 α ppqq + α qqqq = 0.For a point force f = f ˆq at the origin, these equationswere solved by Kelvin [1] and yield: α = − f p log (cid:0) q + p (cid:1) πµ , P = P + f q πµ ( q + p ) . (7)Such traditional incompressible elastic solutions satisfylinear incompressibility, ∇ · u = 0, not Det( F ) = 1. Thusthe resultant fields cannot be substituted into eq. (2), solinear solutions cannot be used as a starting point for evenweakly non-linear calculations. Indeed, if one substitutesthe linear result back into the full almost incompressible energy, eqn. (1), whose behavior we are ultimately inter-ested in, we find that the dilatational energy term is B (Det( F ) − = B ( ∇ · u + Det( ∇ u )) . (8)If u is the usual incompressible linear solutionthe quadratic term vanishes, but the quartic term, B (Det( ∇ u )) ∼ B ( ∇ u ) , does not. This erroneous con-tribution to the energy will become problematic when itis commensurate in size with the leading order quadraticshear term µ Tr (cid:0) ∇ u · ∇ u T (cid:1) ∼ µ ( ∇ u ) , which will occurwhen the strain approaches ∇ u ∼ (cid:112) µ/B . For an almostincompressible material ( B/µ (cid:29)
1) the solution producesunacceptably large volume changes at very small strains,and, in the limit of true incompressibility, it is only validfor asymptotically small strains. In contrast, if B ∼ µ ,the linear regime works well until ∇ u ∼
1, at which pointnon-linear geometric effects guarantee that any linear the-ory must break down. This will be an issue if we wish touse the traditional linear solution as a seed or a boundarycondition for a full finite element calculation utilizing eqn.(1), or if we wish to use it to use it as a test function forthe constrained energy, eqn. (2).To resolve these shortcomings we first implementvolume-conservation exactly by using a coordinate trans-formation originally deployed by Gauss but only recentlyintroduced into elasticity [2,3]. Linearizing the result leadsto a linear theory of elasticity with perfect volume conser-vation that yields qualitatively different results than thetheory that preserves volumes only to linear order. Con-sider a 2-D elastic reference state labeled with the Carte-sian coordinates ( p, q ) and a target state labeled with thecoordinates ( x, z ) = ( p, q )+ u , as shown in fig. 1(a)a. Con-ventionally the deformation between these states is de-scribed by x ( p, q ) and z ( p, q ), so the deformation gradientis F = (cid:32) ∂x∂p (cid:12)(cid:12) q ∂x∂q (cid:12)(cid:12) p∂z∂p (cid:12)(cid:12) q ∂z∂q (cid:12)(cid:12) p (cid:33) . (9)p-2xactly isochoric deformations of soft solidsHowever, if we describe the deformation via a mixed coor-dinate system z ( x, q ) and p ( x, q ), as sketched in fig. 1(a)b,we can implement the volume constraint exactly [2]. Tounderstand this transformation, consider a patch of ma-terial specified by q < q < q and x < x < x shownin fig. 1(a)b. In both the target and the reference statethis patch has two parallel straight sides. Incompressibil-ity requires the area A of the patch to be the same in bothstates: A = (cid:90) q q [ p ( x , q ) − p ( x , q )]d q = (cid:90) x x [ z ( x, q ) − z ( x, q )]d x. (10)If our patch is small in both dimensions, so that q = q + δq and x = x + δx , we can evaluate these integrals tofirst order as δxδq ∂p∂x (cid:12)(cid:12) q and δxδq ∂z∂q (cid:12)(cid:12) x . Equality then yields ∂p∂x (cid:12)(cid:12) q = ∂z∂q (cid:12)(cid:12) x , a condition that is satisfied if we introducea new field ψ ( x, q ) such that p ( x, q ) = ∂ψ∂q (cid:12)(cid:12)(cid:12)(cid:12) x and z ( x, q ) = ∂ψ∂x (cid:12)(cid:12)(cid:12)(cid:12) q . (11)The deformation gradient F ( x, q ) in terms of the new func-tions p ( x, q ) and z ( x, q ) is given by F = (cid:16) ∂p∂x (cid:12)(cid:12) q (cid:17) − − ∂p∂q (cid:12)(cid:12) x (cid:16) ∂p∂x (cid:12)(cid:12) q (cid:17) − ∂z∂x (cid:12)(cid:12) q (cid:16) ∂p∂x (cid:12)(cid:12) q (cid:17) − ∂z∂q (cid:12)(cid:12) x − ∂p∂q (cid:12)(cid:12) x ∂z∂x (cid:12)(cid:12) q (cid:16) ∂p∂x (cid:12)(cid:12) q (cid:17) − . (12)Henceforth we regard all functions as functions of x and q . Rewriting F ( x, q ) in terms of ψ yields F = (cid:32) ψ xq − ψ qq ψ xq ψ xx ψ xq ψ xq − ψ qq ψ xx ψ xq (cid:33) , (13)which automatically satisfies Det( F ) = 1. We now substi-tute this result into eqn. (4), at which point it is natural toregard P and f also as functions of x and q . Starting fromthe undeformed state ψ = xq , in the small strain limit, wewrite ψ ( x, q ) = xq + α ( x, q ) (14)where α (cid:28)
1, and further we assume that |∇ P | (cid:28) F and F − T are then F = (cid:18) − α xq − α qq α xx α xq (cid:19) , (15) F − T = (cid:18) α xq − α xx α qq − α xq (cid:19) . (16)We now expand eqn. (4) to leading order in both α and ∇ P to establish the governing equations of incompressiblelinear elasticity. To first order the partial derivative iden-tities ∂∂x (cid:12)(cid:12) z = ∂∂x (cid:12)(cid:12) q and ∂∂z (cid:12)(cid:12) x = ∂∂q (cid:12)(cid:12) x hold, and hence theexpansion gives µ (cid:18) − α qqq − α xxq α xqq + α xxx (cid:19) − (cid:18) P x P q (cid:19) = − f . (17) |∂ u | qq B = μ B = μ B = μ B = μ -3 -2 -1 ε Fig. 3: Fractional energy discrepancy (cid:15) between the full non-linear finite element energy densities, E nl ( F nl ) (see eqn. (1)),and the approximate energy densities derived from the linearsolutions, E nl ( F ), as a function of strain magnitude | ∂ q u q | .The blue lines use the Kelvin “incompressible” linear solution,fig. 2(b), while the red line uses the exactly incompressiblesolution, fig. 2(a). The data is taken along the material line p =0, q <
0. The Kelvin solution only agrees with the full solutionwhen strains are small compared to (cid:112) µ/B . The curve for theexactly incompressible solution is plotted using B = 100000 µ ,but all modulus ratios considered here give essentially the samecurve, showing the solution is works well for all materials withhigh bulk modulus untill strains approach one. In a region with f = 0, α is again biharmonic [2]: α qqqq +2 α xxqq + α xxxx = 0. These equations are identical to thosegoverning the traditional linear elastic equations, eq. (6),with the important identification p → x . Thus the trulyincompressible response of an infinite 2-D medium to apoint force of f = f ˆq applied at the origin ( x = q = 0) is,by comparison with eqn. (7), α = − f x log (cid:0) q + x (cid:1) πµ P = P + f q πµ ( q + x ) . (18)Similarly, any problem solved within the traditionalstreamline formalism can be imported into the strictlyincompressible formalism with the identification p → x .However, one should not conclude the two are equivalent,as we will now detail. Comparing the strictly incompress-ible and traditional solutions in fig. 2(a-c), we see thatboth produce unphysical infinite displacement and self in-tersection near the point force and are nearly identical farfrom the point of application of the force. However, thereis a significant intermediate region where the Kelvin solu-tion produces visible area changes, while the new solutiondoes not.To assess the relative merits of the two solutions, we usethe commercial finite element package ABAQUS to com-pute the fully non-linear response of a 2-d neo Hookeanmaterial to a force f concentrated on a small rigid linewith a range of bulk moduli B (cid:29) µ , obeying eqn. (1).The elastic full-space is scale free but in terms of an arbi-trary unit-length l , we assume that the force has magni-tude 2 µl and the line has radius a = 0 . l . The dimen-sionless number f / ( µa ) = 80 is thus large, meaning theforce is point-like, being highly concentrated on a smallp-3ohn S Biggins1, Z Wei2,L Mahadevan2,3 xzf (q) (a) Exact Volume Preservation (b) Conventional Kelvin Solution pq F (c) Volume changes in Kelvin so-lution (d) Fully non-linear solution Fig. 2: Deformation caused by a point force of magnitude 2 µl on a 2-D square grid with spacing 0 . l , where l is an arbitrarylength. (a) Exactly area preserving solution given by eqn. (7). (b) Conventional “incompressible” Kelvin solution showingmarked area changes near the point of application of the force. (c) Reference state colored by the area change caused by theKelvin solution. (d) Numerical fully non-linear deformation of a material following eqn. (1) with B = 100 µ . The force is appliedto a rigid disk of radius 0 . l . disk. The rigid disk is embedded in the center of a mate-rial matrix of size 10000 l × l . Due to the left-rightsymmetry, only the right half of the domain is simulatedwith symmetric boundary conditions on the left bound-ary and the other three boundaries being clamped, whilea no-slip boundary condition is applied between the rigidline and the solid. We use a mesh of linear quadrilateralelements with a total of 1641330 degrees of freedom thathas 100 evenly spaced elements along the half-length ofthe rigid line and that coarsens geometrically with dis-tance from it. In Fig. 2(d) we depict one example of sucha deformation, with B/µ = 100, and see that it is almostexactly area-preserving and lacks the infinite displacementand self-intersection of the two linear solutions. We nowcompare the deformation gradient due to the linear so-lution, F associated with the linearized incompressibilitycondition and the new exactly isochoric case, with thefull non-linear solution, F nl , at equivalent material points,by computing the fractional discrepancy between the en-ergy (cid:15) = [ E nl ( F nl ) − E nl ( F )] /E nl ( F nl ). Fig. 3 plots thisfractional discrepancy (cid:15) as a function of the strain ∂ q u q for material points beneath the point of application ofthe force and for materials with a range of ratios B/µ .We see that the traditional Kelvin solution only agreeswith the full non-linear solution when the strain is small,i.e. |∇ u | (cid:28) (cid:112) µ/B , so as the materials become moreincompressible, the traditional solution remains accurateonly for smaller and smaller strains. In contrast, the ex-actly area preserving solution only breaks down at when ∂ q u q ∼
1, at which point geometric non-linearities domi-nate so that no linear theory can work.A similar approach can be implemented in 3-D axisym-metric situations. A deformation from a reference statedescribed with the polar-coordinates ( ρ, θ, q ) to a targetstate described by ( r, φ, z ), would usually be describedby the functions r ( ρ, q ) and z ( ρ, q ), with axisymmetry requiring that φ = θ . We can implement traditional“incompressibility” ( ∇ · u = 0) by introducing a Stokesstream-function β ( r, q ) such that r = ρ − (1 /ρ ) β q and z = q + (1 /r ) β ρ . If instead we describe the deformationvia ρ ( r, q ) and z ( r, q ), we can enforce Det( F ) = 1 by in-troducing the scalar field χ ( r, q ) such that z ( r, q ) = χ r r and ρ ( r, q ) = (cid:112) χ q . (19)In this case, if there is no deformation χ = r q , so, tolinearize, we write χ = r q + β ( r, q ), where β is small.Expanding to first order in β we see that, with the identi-fication r → ρ , the displacements and deformations are al-gebraically identical to those with the conventional Stokesstream-line function.Our point force responses contain regions of diver-gent strain that are beyond any linear theory, so weclose with an example where a distributed load leadsto a finite response. Inspired by active-soft-matter sys-tems such as magnetic-elastomers and muscular-tissues,we consider a two-dimensional elastic full-space contain-ing two horizontal rigid plates, initially at q = ± a andwith width 2 w , that pull towards each other with a to-tal force f . The contribution to α from each plate isthen given by integrating the point force solution alongthe plate weighted by a local force density f ( x ), to get α = (cid:82) w − w f ( x ) x log (cid:0) a + x (cid:1) / (8 πµ )d x , where we mustchoose f ( x ) so that it integrates to f and has the cor-rect profile to keep the plates flat. For a single plate, thiswould entail f ( q ) ( x ) = 1 / (cid:112) − x [4]. For two plateswe were unable to find a closed from solution for f ( q ) ( x )but achieved almost perfect rigidity by taking the form f ( q ) ( x ) = ( A + Bx + Cx ) / (cid:112) − x and fitting A , B and C . We compare this fully area preserving solution to thetraditional-elastic one in fig. 4(a)-4(b), and observe thatthe traditional solution suffers substantial erroneous arealoss between the plates. Although neither solution solvesp-4xactly isochoric deformations of soft solids xz (a) Exact volume preservation F (b) Conventional linear elasticity Fig. 4: Deformation caused by two parallel plates, with unit length and separation, embedded in a two dimensional elasticmedium being pulled together. (a) is calculated using the exactly area preserving techniques outlined in this paper while (b)is calculated using traditional linear elasticity. The colour shows the (erroneous) fractional volume change in the traditionalsolution. the full non-linear problem, our deformation field offers anadmissible solution for consideration in a non-linear the-ory — its fully non-linear energy can be evaluated andit can serve as a starting point in an attempt, either nu-merically or analytically, to minimize the fully non-linearenergy. A detailed account that builds on our 2d analysisof exactly isochoric deformations of soft solids document-ing all the major exactly incompressible bulk and surfaceGreen’s functions for 2-D and axisymmetric elastic full andhalf spaces with and without large pre-strains is availablehere [5].A natural objection at this point might be that it is in-consistent to linearize our theory and satisfy incompress-ibility perfectly, and that the new perspective thus addsnothing to the traditional stream-line formalism. How-ever, we see considerable value in perfect volume preser-vation. Firstly, as shown above, it extends the range of va-lidity of the linear theory from strictly infinitesimal strainsto strains approaching unity. Secondly, it produces fieldsthat can be substituted into the full non-linear energy asa starting point for building energy-estimates, numericalminimization or rigorous bounds. It thus provides a bridgeconnecting linear and non-linear elasticity in incompress-ible systems. Thirdly, it will be a useful approximatemethod in contexts where the full non-linear problem isintractable and incompressibility is critical to understand-ing the key physics. Finally, our approach may be usefulin developing numerical elasticity and elastodynamic inte-grators for incompressible systems. Our approach is analo-gous to the use of symplectic integrators that have becomecommonplace in classical dynamics. These also introducea different representation of a problem that guarantees theexact preservation of key conserved quantities (often en-ergy or momentum), and their use in numerics leads tostabler, more robust and faster simulations. In the elas- ticity of compressible materials, numerically one is oftenstymied by having to take an asymptotically small time-step because the speed of sound diverges. A scheme basedon our method would allow the implementation of strictincompressibility with a finite time-step. It is by thuslinking the linear and non-linear theories of incompress-ible soft solids that we anticipate the theory outlined inhere will find application.
REFERENCES[1] W. Thomson (Lord Kelvin). Displacement due to a pointload in an indefinitely extended solid.
Mathematical andPhysical Papers (London) , 1:97, 1848.[2] M. Ben Amar and P. Ciarletta. Swelling instability ofsurface-attached gels as a model of soft tissue growth un-der geometric constraints.
J Mech Phys Solids , 58(7):935–954, 2010.[3] M.M. Carroll. On generating isochoric finite deforma-tions.
Journal of Elasticity , 88(1):1–4, 2007.[4] K.L. Johnson and K.L. Johnson.
Contact mechanics .Cambridge Univ Pr, 1987.[5] J.S. Biggins Z. Wei and L. Mahadevan.
2D and axisym-metric exactly incompressible elastic Green’s functions .arXiv:1407.1405 [cond-mat.soft] http://arxiv.org/abs/1407.1405http://arxiv.org/abs/1407.1405