HYPERELASTIC SWELLING OF SPHERES AND CYLINDERS AND ITS GENERALIZATION TO ELASTIC INTERNALLY BALANCED MATERIALS By Vahid Zamani A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering – Doctor of Philosophy 2018 ABSTRACT HYPERELASTIC SWELLING OF SPHERES AND CYLINDERS AND ITS GENERALIZATION TO ELASTIC INTERNALLY BALANCED MATERIALS By Vahid Zamani Swelling as a notion of free volume change is typically due to some added mass procedures. We use modified constitutive laws that incorporate swelling into a continuum mechanics treatment. By incorporating local volume change (swelling) as a parametric constraint into the conventional theory of hyperelasticity it is possible to model a variety of swelling effects. We consider these effects in the study of certain boundary value problems for spherical and cylindrical finite deformations. In addition to the traditional hyperelastic model, we also employ a relatively new type of constitutive treatment, termed internal balance. The theory of internally balanced materials employs energy minimization to obtain an additional balance principle to treat more complex behaviors. This is useful when conventional elastic behavior is modified by substructural reconfiguration. Hence, we also formulate our problems in the context of the internally balanced material theory for the case of cylindrical deformation where the results are compared to that of the conventional hyperelastic model. For thick spherical shells, the incompressible hyperelastic Mooney-Rivlin constitutive model allows for response to pressure-inflation that could either be globally stable (a monotonic pressure-radius graph) or could instead involve instability jumps of various kinds as pressurization proceeds. The latter occurs when the pressure-radius graph is not monotonic, allowing for a snap-through bifurcation that gives a sudden burst of inflation. Internal swelling of the material that makes up the shell wall will generally change the response. Not only does it alter the quantitative pressure-inflation relation but it can also change the qualitative stability response, allowing burst phenomena for certain ranges of swelling and preventing burst phenomena for other ranges of swelling. These issues are examined both for the case of uniform swelling for the case of a spatially varying swelling field. For cylindrical deformations, we examine the finite strain swelling of a soft solid plug within a rigid tube of circular cross section. The eventual channel wall contact as the swelling proceeds generates a confinement pressure that increases as the plug expands. We consider plug geometries that incorporate an internal channel as well as a simpler case of a solid plug. For the case of a plug with a channel, the wall contact now gives a deformation in which swelling combines axial lengthening with internal channel narrowing. Of particular interest is the closing behavior as the swelling proceeds and we treat the problem using asymptotic expansions. Finally, the same problem is examined in the context of the internal balance constitutive theory. To my great parents, Seyfollah and Behjat, who were my first teachers, and to my wonderful wife Mahdieh and my adorable sisters, and dedicated to our little angel who is on her way coming to this world. iv ACKNOWLEDGEMENTS I would like to express my deepest respects to my graduate adviser Dr. Thomas Pence who patiently mentored me throughout the years, who showed me the paradigm of professional language and strength of simplicity. I always figure him in my mind as an exemplar of trust and honesty that reinvigorated my hope for a better face of the world. I also should highly appreciate my other committee members and MSU engineering faculties and staff who were indeed helpful in this part of my education. My dear love Mahdieh educated me and planted the finest of hope in me. During the hardest of the years of my education her true love made me a better person for our lives. Without her love this achievement would not have come true. My parents and my sisters, who are my true origin and I could not live a lovely life without them, whom I gave them pain by being far far away from my home country during the years, are of my best wishes for the rest of my life. v PREFACE "What we observe is not nature in itself, but nature exposed to our method of questioning." W. Heisenberg vi TABLE OF CONTENTS LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . . . . . . . . . . . . . . . . . ix KEY TO ABBREVIATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv CHAPTER I Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Hyperelastic swelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2 Isotropic material behavior subject to swelling . . . . . . . . . . . . . . . . . . . 12 CHAPTER II Burst Instability in Uniform Swelling of Hyperelastic Spherical Shells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.1 Kinematics for radial inflation of a spherical shell . . . . . . . . . . . . . . . . 18 2.2 The role of swelling in the pressure-inflation relation . . . . . . . . . . . . . 24 2.2.1 Uniform expansion occurs for homogeneous swelling in the absence of pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.2.2 Qualitative behavior in the absence of swelling . . . . . . . . . . . . 26 2.2.3 Quantitative determination of the inflation behavior . . . . . . . . 28 2.3 The swellable Mooney-Rivlin material . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.3.1 Inflation behavior prior to swelling . . . . . . . . . . . . . . . . . . . . . 35 2.3.2 Inflation graph sequences for increasing swelling . . . . . . . . . . 38 2.3.3 Swelling dependent material stiffness parameters . . . . . . . . . . 40 2.3.4 Effect of the constitutive exponents m and n on the transitional swelling value va/c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.3.5 Numerical illustration with swelling-dependent material stiffness parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 2.4 Swelling induced burst . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 CHAPTER III Non-Uniform Swelling Field of Hyperelastic Spherical Shells . . 56 3.1 A family of swelling fields with the same added mass . . . . . . . . . . . . . 58 3.2 Kinematics of the spherical deformation with non-uniform swelling field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.3 Behavior type independent of distribution for fixed added mass . . . . . 64 3.4 Inflation behavior instability due to mass redistribution . . . . . . . . . . . . 70 CHAPTER IV Channel Confinement Swelling of Hyperelastic Plugs and Tubes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.2 The homogeneous deformation of laterally confined swelling . . . . . . . 75 4.3 The confinement boundary value problem for an annular plug . . . . .. . 79 4.4 Annulus contact for the neo-Hookean type swelling model . . . . . . . . . 84 4.4.1 Existence and uniqueness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.4.2 Numerical demonstrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.4.3 Asymptotic behavior for large swelling values . . . . . . . . . . . . 95 vii CHAPTER V Channel Confinement Swelling of Internally Balanced Material Plugs and Tubes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.1 Internally balanced elastic materials . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.1.1 Retrieving the hyperelastic theory . . . . . . . . . . . . . . . . . . . . . 105 5.1.2 Isotropic materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.2 Homogeneous deformation for the internally balanced material model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.3 Cylindrical deformation in internally balanced material model . . . . 123 5.3.1 Internal balance requirement . . . . . . . . . . . . . . . . . . . . . . . . . . 124 5.3.2 Stress and equilibrium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 5.3.3 Explicit solution of the special case v=β2/3 . . . . . . . . . . . . . . 132 CHAPTER VI Concluding Remarks, Speculations and Broader Connections . . . . . . . . . . . . . . . . . . . .. . . . . . . . 136 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . 141 viii LIST OF FIGURES Figure 2.1 Inflation graphs showing three qualitatively different types of be- havior (a)-(c) in the absence of swelling. These particular graphs cor- respond to W given by (1.16), all with thickness ratio ξ = 0.5. The differences are due to the values of d1 and d2. Here: d1 = 4d2 (top); d1 = 9d2 (middle); and d2 = 0 (bottom). . . . . . . . . . . . . . . . . . Figure 2.2 Graphs for G(η, v) corresponding to the three inflation curves in Figure 2.1. The function G(η, v) is computed on the basis of (2.32) using (1.17) and taking v = 1. This is equivalent to using (1.16) and ultimately gives the expression (2.43) that we examine in more detail later. . . . Figure 2.3 G graphs from Fig. 2.2 (solid) along with the corresponding H graphs (dashed) for thickness ratio ξ = 0.5. Each point on a G graph is shifted to the left to give a corresponding point on the H graph. This shift is small if ξ is close to one (a thin shell). Here, because ξ is not very close to one, the nonuniform shift distorts the curves, however the basic monotonicity properties do not change. . . . . . . . . . . . . . . Figure 2.4 Qualitative behavior of the inflation graph for the Mooney-Rivlin model W = d1(I1 − 3) + d2(I2 − 3) as a function of material parameter α = d1/(d1 + d2) and thickness ratio ξ = Ri/Ro. The curve ξ = ξa/c provides a transition between type (c) and type (a)-behavior. . . . . . . Figure 2.5 Inflation graphs for the Mooney-Rivlin-type model (1.17) using (2.45) with α = 0.85, and thickness ratio ξ = Ri/Ro = 0.3. All the inflation graphs exhibit type (a) behavior. . . . . . . . . . . . . . . . . Figure 2.6 Inflation graphs for the same material as in Fig. 2.5 (i.e., (1.17) and (2.45) with α = 0.85), but now the thickness ratio ξ = 0.7. This corresponds to a relatively thinner walled structure. All the inflation graphs now exhibit type (c) behavior. . . . . . . . . . . . . . . . . . . . Figure 2.7 Transitional swelling value va/c versus ξ = Ri/Ro for the Mooney- Rivlin-type model (1.17) with parameters µ, α, m and n in (2.49). The transitional swelling value va/c is independent of µ and is dependent on n and m only via the difference m− n. These plots are for m− n = 2/3. For a given α-curve the inflation graph exhibits type (c) behavior if (ξ, v) is in the region above the curve and type (a) behavior if (ξ, v) is in the region below the curve. . . . . . . . . . . . . . . . . . . . . . . . . . . 27 30 32 37 38 39 45 ix Figure 2.8 Transitional values for swelling va/c versus ξ = Ri/Ro for Mooney- Rivlin-type model (1.17) with parameters µ, α, m and n in (2.49). The transitional swelling value va/c is dependent on α as shown but is in- dependent of µ. The curves are dependent on m and n only via the difference m − n. This figure is for m − n = −2/3. For a given α-curve the inflation graph exhibits type (c) behavior when (ξ, v) is in the region that is below the curve and type (a) behavior in the region that is above the curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 2.9 Inflation graphs for the Mooney-Rivlin-type model (1.17) using (2.49) with α = 0.85, m = 2/3, n = 0, and thickness ratio ξ = 0.3. The inflation graphs exhibit the type (a) behavior for 1 ≤ v < 1.54 and type (c) behavior for v > 1.54. . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 2.10 Inflation graphs for the Mooney-Rivlin-type model (1.17) using (2.49) with α = 0.85, m = 2/3, n = 0, and thickness ratio ξ = 0.7. The inflation graphs exhibit type (c) behavior for all v ≥ 1. . . . . . . . . . Figure 2.11 Inflation burst caused by increasing v at fixed ∆P = 0.258µ for the inflation graphs from Fig. 2.6. Prior to swelling the pressurization ∆P = 0.258µ has given a mild radial increase (from si = ri/Ri = 1 to si = ri/Ri = 1.14 on the v = 1 curve). Now increasing v at this fixed ∆P gives a continuous increase of si with v (dashed red line) until encountering the inflation graph for v = 2 where there is a local maxi- mum. Further increase of v requires a jump across to the other increasing branch of the v = 2 curve (solid red segment). This corresponds to an inflation burst with radial increase from si = 2.37 to si = 4.32. . . . . Figure 2.12 Inflation burst showing si = ri/Ri vs. v at ∆P = 0.258µ caused by an increase in the swelling parameter v. Locations denoted by • provide correlation with the inflation graphs depicted in Figure 2.11. . . Figure 2.13 Inflation burst caused by increasing v at fixed ∆P = 1.16µ for the inflation graphs from Fig. 2.9. Initially, the radius increases contin- uously, first with v = 1 as ∆P increases from zero to 1.16µ and then at this fixed ∆P as v increases to v = 2 (dashed red segment). At v = 2 there is a jump from the first increasing branch to the second increasing branch after which a continuous increase is again the case. . . . . . . . Figure 3.1 The family of swelling distributions (3.3) with constants (3.12) for ξ = 0.5 that are parameterized with respect to vi such that v(R) ≥ 1. All distributions have the same amount of overall added mass associated with vuni = 1.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 49 50 51 52 55 62 x Figure 3.2 Type (c) behavior in response to pressure-inflation for the six different swelling fields in figure 3.1; using fixed properties (ξ, α) = (0.5, 0.86) in the Mooney-Rivlin-type material model (1.17) and (2.46). Figure 3.3 Type (a) behavior in response to pressure-inflation for the six different swelling fields in figure 3.1; using fixed properties (ξ, α) = (0.5, 0.83) in the Mooney-Rivlin-type material model (1.17) and (2.46). Figure 3.4 Qualitative behavior of the inflation graph for the Mooney-Rivlin model with uniform swelling fields, previously shown in Figure 2.4. Here 18 points are chosen along the transition curve (red dots) to be used in (3.29) for numerical integration. The inflation graphs with the circled parameter choices are shown in Figs. 3.2, 3.3, 3.5 and 3.6. . . . . . . . 67 68 69 Figure 3.5 Inflation behavior for the two limits of swelling distribution (vi,min, vi,max) in figure 3.1; using fixed parameters chosen from Figure 3.4 with (ξ, α) = (0.1, 0.94) (on the left) showing type (a) behavior and (ξ, α) = (0.1, 0.96) (on the right) showing type (c) behavior. The Mooney-Rivlin-type ma- . . . . . . . . . . . . . . . . terial model is based on (1.17) and (2.46). 70 Figure 3.6 Inflation behavior for the two limits of swelling distribution (vi,min, vi,max) in figure 3.1; using fixed parameters chosen from Figure 3.4 with (ξ, α) = (0.9, 0.82) (on the left) showing type (a) behavior and (ξ, α) = (0.9, 0.83) (on the right) showing type (c) behavior. The Mooney-Rivlin-type ma- terial model is based on (1.17) and (2.46). . . . . . . . . . . . . . . . . Figure 3.7 Inflation burst due to redistribution of the fixed added mass in type (c) behaviors; using fixed properties (ξ, α) = (0.5, 0.86) in the Mooney-Rivlin-type material model (1.17). . . . . . . . . . . . . . . . Figure 4.1 Solid cylinder with the radius Ro and the confinement pipe with . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . radius Rc > Ro. Figure 4.2 Representation of the plug swelling within a rigid pipe. For v < λ3 lat the plug has yet to make contact with the pipe wall. Contact first occurs when v = R3 lat all further swelling is directed into Z-direction extension. This generates the pressure Plat between the plug and the pipe wall. . . . . . . . . . . . . . . . . . . . . lat. For v > λ3 c/R3 o = λ3 Figure 4.3 Confinement Pressure Plat as a function of swelling v in a solid cylinder taking the outer radius Ro = 1 (so that λlat = Rc). The graphs are for the pipe radii Rc = 1.51/3 on the left and Rc = 21/3 on the right. . . . . . . . . . . . Here Plat is normalized by the material modulus µ. xi 71 72 76 76 78 Figure 4.4 Axial elongation λz as a function of swelling v for the expanding annular plug with inner radius Ri = 1/2 and outer radius Ro = 1 (so that λlat = Rc and ζ = 1/2). Wall contact occurs when v = λ3 lat. The graphs are for two separate cases of outer pipe radius: Rc = 1.145 and Rc = 1.260, which are chosen so as to give contact v values of 1.5 and 2, respectively. The slope of the curves immediately after contact are given by (4.31). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 Figure 4.5 Deformed inner radius ri as a function of swelling v for the wall contact cases from Fig. 4.4. Prior to wall contact the swelling v < λ3 lat and the homogeneous deformation causes ri to increase. After contact, when v > λ3 lat, the inner radius is monotonically decreasing with swelling. 91 Figure 4.6 Confinement pressure Plat as a function of swelling v for the wall contact cases from Figs. 4.4 and 4.5. Here Plat is normalized by the . . material modulus µo. The slope at first contact is given by (4.32). 92 Figure 4.7 Confinement pressure Plat for a solid plug (Ri = 0) and for tubes (annulus plugs with Ri = 1/10, 1/4, 1/2, 3/4) as a function of swelling v. In all cases Ro = 1 (so that λlat = Rc). The confining radius is Rc = 1.51/3 so as to be consistent with one of the cases shown previously in Figs. 4.4, 4.5 and 4.6. Plat is again normalized by the material modulus µo. The zero channel radius limit (ζ = 0) recovers the neo-Hookean type curves in Fig. 4.3. Thinner walls (larger ζ) give less contact pressure for the same value of swelling. . . . . . . . . . . . . . . . . . . . . . . . . . Figure 4.8 Asymptote of the deformed inner radius ri as a function of swelling v in a hollow tube taking the inner radius Ri = 1/2 and the outer radius Ro = 1 (so that λlat = Rc and ζ = 1/2) with Rc = 1.51/3. The graphs are for the pure numerical solution (red curve in the middle) and the asymptotic expansion (4.38) with only the leading term (green curve on the bottom) and with the additional first correction term (blue curve on the top). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 99 Figure 5.1 Confinement pressure Plat = Π(λlat, v) as a function of swelling v taking ˆα = 1, α∗ = 2 and Ro = 1 (so that λlat = Rc). The three graphs on the left are for Rc = 1.51/3. The three graphs on the right are for Rc = 21/3. For each Rc the graph of Plat is given for three different values of α. Here Plat is normalized by α + (ˆαα∗)/(ˆα + α∗). Because of (5.56) this accounts for the common slope at the rightmost departure (Rc = 21/3) being less than the common slope at the leftmost departure (Rc = 1.51/3). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 xii Figure 5.2 Three solutions of ρ(s; v, β, λz) to the requirement (5.67) with one representative set of v = 2 and λz = 2. Two of the non-trivial solutions are obtained for different ration β = 2 and β = 1/2. If β is taken as β = v2/3 = 22/3 then the obvious solution is ρ(s) ≡ 0. . . . . . . . . . . 128 Figure 5.3 Axial elongation λz on a hollow tube as a function of swelling v taking α = 0 ˆα = 1, α∗ = 2 and Ri = 1/2, Ro = 1 (so that λlat = Rc). The graphs are the solution to (5.89) and are for two separate cases of outer pipe radius: Rc = 1.51/3 and Rc = 21/3, which are chosen so as to give contact v values of 1.5 and 2, respectively. . . . . . . . . . . . . . 131 Figure 5.4 Inner radius ri on a hollow tube as a function of swelling v taking α = 0 ˆα = 1, α∗ = 2 and Ri = 1/2, Ro = 1 (so that λlat = Rc). The graphs are the solution to (4.9) and are for two separate cases of outer pipe radius: Rc = 1.51/3 and Rc = 21/3, which are chosen so as to give contact v values of 1.5 and 2, respectively. . . . . . . . . . . . . . . . . 132 Figure 5.5 Confinement pressure Plat on a hollow tube as a function of swelling v taking α = 0 and β = 1/2 (ˆα = 1, α∗ = 2) and Ri = 1/2, Ro = 1 (so that λlat = Rc) for longer graphs with Rc = 1.51/3 and Rc = 21/3. The graphs of β = 4, 10 and ∞ are also plotted for Rc = 1.51/3. Note that the graph of β = ∞ retrieves the hyperelastic behavior shown in Figure 4.6 based on (5.45). Here Plat is normalized by ˆαα∗/(ˆα + α∗). 133 xiii KEY TO ABBREVIATIONS α: Mooney-Rivilin material parameter ˆα, α∗: internal balance material parameters β: internal balance material parameter ratio F: deformation gradient tensor ξ, ζ: thickness ratio parameter ˆξ, ξ∗: internal balance deformation variables I: identity tensor λi: principal stretches λlat: lateral elongation λz: axial elongation µ: elastic modulus p: Lagrange multiplier Plat: lateral confinement pressure q: internal balance Lagrange multiplier Ri: inner radius Ro: outer radius Rc: confinement radius ρ: internal balance variable s: radial deformation ratio (r/R) T: Cauchy stress tensor Φ(i): strain energy densities v: swelling amount xiv CHAPTER I Introduction Swelling can be viewed as a general phenomena that represents free volume change, typically due to mass addition resulting from some diffusive or transport mechanism. It is seen that many biological tissues and cells exhibit this type of volume and shape change as a result of biological growth, hydration and mass exchange. In many cases osmotic pressure is the causal agent that drives water and other mass transport across bio-membranes such as those surrounding red blood cells and inter-cellular vesicles and lipids Graf et al. (1995); Vinod Kumar and Demeke (2011); Li et al. (2013). For our purposes, swelling is regarded as a general process that encompasses free-volume change at the microscopic level. This would typically be due to mass addition but other fine scale processes of a mechanical or chemical nature can also be regarded as contributing to volume change. Polymers, elastomers and hydrogels naturally swell when exposed to liquid or when subject to high humidity Treloar (1975); Stuart et al. (2010); Drozdov (2013). Biological tissues and cells exhibit volume and shape change under similar processes of hydration and mass exchange Van der Sman (2015), but also more generally as a result of biological growth Goriely et al. (2010); Sadik et al. (2016). The scientific literature on swelling is vast, and can be approached from a vari- ety of perspectives: material science, physical chemistry, continuum mechanics, etc. We focus on continuum mechanics in which case the literature is still extensive (and 1 overlapping): mixture theory, biphasic theory, transport theory, poroelasticity, large strain at the outset, small strain on top of a base swollen state, etc, as described for example by a variety of related approaches (see e.g., Bowen (1980); Wineman and Rajagopal (1992); Ateshian (2007); Hong et al. (2008); Markert et al. (2008); Duda et al. (2010); Chester and Anand (2010); Pence (2012); Drozdov et al. (2013); Sel- vadurai and Suvorov (2016)). The treatment given in this thesis is of a generalized hyperelastic nature. By this we mean it specifically does not seek to model the mi- crostructural mechanisms associated with the swelling process, although, as discussed for example in Baek and Pence (2011), with suitable modification it could be related to the broader frameworks mentioned above . The mechanical consequences of these swelling-involved processes are significant. Swelling can lead to qualitative changes in the material’s mechanical properties such as the deformability, stability of the overall structure (possibly triggering various bifurcation phenomena associated with localization, buckling and other forms of non- uniqueness), maximum stretch and rupture of membranes and also the capacitance of cells. These are mainly because of the fact that material properties and structure of the solid may be altered as a result of the mass exchanges, leading to associated re- configurations in stresses and deformations and in some cases instabilities. Hence, the significant impacts on the mechanical properties of the materials require attentions in the study of the interaction between swelling and other mechanical effects. This has been a major issue in many recent studies. This includes the study of Li et al. (2013) where an impulsive-like forcing is employed to quantify the yield strain thresh- old of red blood cell’s membranes before rupture. Vinod Kumar and Demeke (2011) uses spherically symmetric deformation with neo-Hookean model to analyze infected red blood cell’s membrane stress-stretch behavior. Nagel et al. (2009) investigates the effects of cellular stiffening and changes in material properties on the damage 2 evolution in deep tissue pressure-induced injury. The effect of osmotic swelling on the conductance and capacitance of hepatocytes’ membrane is studied in Graf et al. (1995). Van der Sman (2015) represents a hyperelastic neo-Hookean model for hy- dration of cellular tissue under pressurization of the internal cavity where an elastic shell undergoes inhomogeneous deformation. Other related studies include those of Gibbons and Klug (2008) and Evans et al. (2003). The interaction of mechanical responses and absorption of a swelling agent arises in a variety of contexts that go beyond biology. In particular, it associates with many important phenomena in porous absorbent or polymeric solids, elastomers and hydro- gels exposed to any surrounding liquid. Among those, hydrogels play an increasingly important role in a wide range of applications, such as in tissue engineering and smart optical systems Stuart et al. (2010). With regards to this kind of interaction, many continuum theories and methods have been developed over the years. In par- ticular Deng and Pence (2010); Baek and Pence (2011); Ben Amar and Ciarletta (2010); Duda et al. (2010); Chester and Anand (2010); Duda et al. (2011); McMa- hon et al. (2010); Goriely et al. (2010); Van der Sman (2015) treat the emerging problems and capture new behaviors in mechanical responses. Theoretical analysis of residual stresses in growing elastic bodies have shown that swelling-induced stresses can initiate mechanical instabilities such as elastic cavitation Pence and Tsai (2006); McMahon et al. (2010). The possible role of mechanical stress in the opening of cavities in elastic cylindrical model of aerenchyma tissue has also been considered in Goriely et al. (2010). As an another recent example Baek and Pence (2011) em- ployed a variational method to obtain governing equations of large deformations of elastomeric gels. Their model was employed to treat both saturated and unsaturated gels in equilibrium and subject to loading at the gel interface. Similar treatments on this kind of gels is found in Deng and Pence (2010). Also swelling instability of 3 surface-attached gels has been studied in Ben Amar and Ciarletta (2010) as a model of soft tissue growth. They have proposed a theoretical framework for the mechanical treatment of the growth of soft materials under geometrical constraints. In the rest of this chapter we review the formulation of swelling in a hyperelas- tic framework and its implementation in isotropic materials. In this context we also introduce a somewhat broader framework (see Eq. (1.5)) and some new material models (see Eq. (1.12)) to incorporate swelling which involves a straightforward gen- eralization of the conventional hyperelastic theory (Pence and Tsai (2005b); Tsai et al. (2004)). In this regard, chapters II to V constitute the new contributions of this research. In Chapter II we focus attention on a class of swellable hyperelastic materials in order to examine how the constitutive theory affects the spherically symmetric expansion of a pressurized hollow sphere. The inflation response of a hyperelastic sphere in the absence of swelling is a classical problem in finite deformation contin- uum mechanics. As is well known, the resulting pressure-expansion response is not always given by a monotonically increasing graph. The sphere may be thick or thin. The thin wall limit corresponds to a hyperelastic swellable membrane. In the absence of swelling the class of materials corresponds to a classical Mooney-Rivlin material. The inflation response for a hollow sphere composed of the classical Mooney-Rivlin material has been extensively studied. We especially focus on a detailed characteri- zation by Carroll (1987) that provides conditions that determine when the inflation response is monotone versus when it is not. We use this characterization to describe the non-swelling response of the sphere problem under consideration here. We then generalize the analysis so as to determine how swelling affects the outcome. In par- ticular, we exhibit swelling induced transitions between monotone and non-monotone 4 inflation curves. When this happens, various inflation jump events can be triggered. We provide a systematic framework for understanding and predicting these transi- tions, and we discuss the ramifications of these transitions in terms of a snap-through type bifurcation (a swelling induced burst). In Chapter III, we again consider a spherical inflation but now under a nonuni- form swelling distribution. A family of swelling distributions is defined such that each member of the family has the same added swelling mass. We then examine the stability behavior within a family with common added mass but different spatial dis- tributions of that mass. Numerical analysis of several examples support a hypothesis that the qualitative behavior is independent of the distribution so long as the amount of added mass remains fixed. Numerical analysis of several examples that support this hypothesis are presented. Finally, an inflation instability is demonstrated within a certain range of material parameters. In Chapter IV we continue to employ a conventional hyperelastic model that con- siders the impact of swelling for the different geometry of cylindrical deformation. In this regard we obtain the responses of both a solid plug, and of a plug with an internal channel, both of which are confined by a rigid outer wall. Sufficiently large swelling leads to deformation with wall contact in which case the confinement pressure is de- termined. For the case of a plug with a channel, asymptotic analysis is employed to investigate the channel closure. It is also the purpose of this thesis to explore the theory of internally balanced elas- tic materials for describing complex deformational response in solids that are subject to swelling. Hence, in Chapter V , we study the boundary value problem of cylindrical deformation that is laid out in Chapter IV but in the context of internally balanced 5 material theory. As such, we review the internal balance theory and the incorporation of swelling in order to formulate the problem and investigate the aspects of the theory that accounts for swelling. In order to plot a picture of the structure of the formulations in this study we should note that the ground framework of this study is in the context of hypere- lasticity. The equations are obtained by minimizing the energy functional involving the elastic energy density and thus our framework is based on the equilibrium defor- mations that are governed by the stress equations of equilibrium. These equations associated with unknown parameters and specific boundary conditions can sometimes be solved analytically and in other cases are solved with the help of numerical treat- ments although we mainly seek to explore analytical solutions and explanations rather than studying numerical procedures. 6 1.1 Hyperelastic swelling As is standard in finite deformation continuum mechanics, let X be a generic position vector in a reference configuration ΩX that is regarded as the state of an unloaded body prior to swelling. The deformation under consideration maps the coordinates of the reference configuration to the coordinates x in the current deformed configuration denoted by Ωx. The gradient of the map x = χ(X) is the tensor F = ∂x/∂X. The loading is described in the standard way in terms of boundary tractions and body forces. In addition the material expands and contracts, and this is described in terms of change in its natural free volume. This volume change is referred to as swelling in this study, and can vary from point to point in the form of a swelling field v = v(X) with v > 1 for swelling expansion and 0 < v < 1 for deswelling. Throughout this thesis, v is treated as a prescribed quantity. In other words it can be viewed as a control parameter. It follows that the appropriate volume constraint on the deformation is detF = v, (with v > 0). (1.1) A theory that uses (1.1) provides a generalization of the conventional theory for in- compressible materials. In the context of hyperelasticity such a framework is demon- strated in Pence and Tsai (2005b, 2006); Tsai et al. (2004); Zamani and Pence (2017) for isotropic material behavior and in Demirkoparan and Pence (2007a,b, 2015a,b, 2017) for anisotropic materials. This framework has been utilized to guide the design of actuators Fang et al. (2011) and to analyze the behavior of biological soft tissue Gou and Pence (2016). The generality of the condition (1.1) permits the phenomenological modeling of a 7 variety of physical swelling processes. More specificity as regards the precise swelling mechanism will generally lead to additional conditions on the physical model ex- pressed in terms of additional constraints and equations. For example, an elastomer in a solvent bath that swells due to solvent uptake might take v in (1.1) in the form v = 1 + ϑmol cref where cref is the concentration of solvent molecules within the elas- tomer network (per unit reference volume) and ϑmol is the molecular volume of an individual solvent molecule. In this case our treatment can apply, but subject to v ≥ 1. Additional equations beyond those that we present here will then enter for the purpose of determining cref on the basis of poroelastic diffusion. In addition, such a description presumes a completely dry (fully dessicated) reference configuration for the elastomeric constituent (which would in general differ from a natural (stress-free) configuration). Whether or not such a reference configuration description is useful may then be highly problem dependent. For additional discussion on this and other issues with respect to such more specific solute/solvent systems we defer to expert sources such as Drozdov and Christiansen (2013). Complete detail on the general hyperelastic modeling treatment employed here is provided and specialized in Section 1.2 to the case of isotropic solid materials that swell. Generalizations of conventional hyperelastic models, such as the neo-Hookean model, so as to incorporate swelling are given. The hyperelastic treatment of swelling is based on an elastic energy density W as a function of both F and v. The hyperelastic energy density is frame-invariant and this requires W to depend on F only through the right Cauchy-Green deformation tensor C = FT F. Thus W = W (C, v). In the absence of body forces the equilibrium equation is div T = 0 where T is the Cauchy stress tensor, T = 2 v F ∂W ∂C FT − pI. (1.2) 8 Here p is the hydrostatic pressure associated with the volume constraint (1.1). It is to be emphasized that p in (1.2) arises as a Lagrange multiplier associated with the constraint (1.1). In particular, it is not simply related to the type of pore pressure variable that is present in say hyperelastic mixture theory. In using (1.2) it is important to realize that W gives the stored energy density with respect to a reference frame that is unswollen. Alternatively, one could consider a reference configuration in which the material is uniformly expanded. Such a uniform expansion would have deformation gradient F∗ = v1/3I. Thus the mapping from the unswollen reference configuration to the current configuration that passes through the uniformly swollen reference configuration is given by F = ˆFF∗, F∗ = v1/3I. (1.3) (1.4) with This makes ˆF = v−1/3F and det ˆF = 1. Analysis can now proceed with respect to the deformation gradient ˆF measured from the uniformly swollen reference configuration, and so one can similarly define ˆC = ˆFT ˆF and hence ˆC = v−2/3C. Let ˆW ( ˆC, v) be the stored energy density per unit volume with respect to this uniformly swollen reference configuration. The overall stored energy calculated by integrating ˆW using the swollen reference frame must be equal to the overall stored energy calculated by integrating W using the unswollen reference frame. This means that W = v ˆW . The material is incompressible with respect to the uniformly swollen reference configuration, and so the Cauchy stress 9 calculated using ˆW is T = 2 ˆF ˆFT − ˆpI ∂ ˆW ∂ ˆC (1.5) where ˆp is the hydrostatic pressure associated with the constraint det ˆF = 1. This constraint can now be interpreted as the standard incompressibility constraint that the deformation (as measured from the uniformly swollen reference configuration) is isochoric. Consistency of these different procedures demands that (1.5) and (1.2) yield the same result for T. This is immediately verified by examining the derivative expression on the right side of (1.2) using F = v1/3 ˆF, C = v2/3 ˆC, W (C, v) = v ˆW ( ˆC, v) and invoking the chain rule: 2 v F ∂W ∂C FT = 2 v v1/3 ˆF v ˆW ( ˆC, v) (cid:16) (cid:16) (cid:17) ∂ (cid:124) ∂ ˆC (cid:16) (cid:17)T v1/3 ˆF = 2 ˆF ∂ ˆW ∂ ˆC ˆFT . (1.6) (cid:17) ∂ ˆC ∂C(cid:124)(cid:123)(cid:122)(cid:125) v−2/3I(cid:2)I (cid:125) (cid:123)(cid:122) v1/3 ∂ ˆW ∂ ˆC From this1 it follows that taking ˆp = p causes (1.5) and (1.2) to yield the same result for T. As the reader is probably aware, a multiplicative decomposition of the form F = ˆFF∗ from (1.3) is commonly used in continuum mechanics to describe a va- riety of physical processes. Having invoked it here to describe swelling by taking the special case F∗ = v1/3I and requiring det ˆF = 1 retrieves the well known formulation attributed to Flory (see e.g., the discussion in Chester and Anand (2010)). In addi- tion, one can now contemplate other descriptions for swelling that invoke F = ˆFF∗ with det ˆF = 1 but which weakens the condition F∗ = v1/3I to something less restric- tive, such as detF∗ = v. This could offer a route to the consideration of swelling processes in which a uniform expansion in all directions is not energetically preferred, such as might be the case if there are microstructural processes that would serve to 1 I (cid:2) I in (1.6) is the fourth order identity tensor, i.e., (I (cid:2) I)ijkl = δikδjl . 10 favor a non equi-axial expansion. Such considerations are formulated in and explained by the internal balance theory in chapter V. 11 1.2 Isotropic material behavior subject to swelling The equivalent equations (1.2) and (1.5) apply to both isotropic and anisotropic behavior. The latter includes fiber reinforced materials, in which case especially designed patterns of fiber stiffening can elicit specialized deformation modes as the material swells. This study focuses solely on isotropic behavior, in which case the dependence of W on C can only be through the invariants of C, I1 = C : I = tr C, I2 = 1 2 (cid:0)(C : I)(C : I) − C : C(cid:1) = 1 (cid:0)(tr C)2 − tr C2(cid:1). 2 (1.7) The third invariant I3 = det C is equal to v2 by virtue of (1.1) and thus W = W (I1, I2, v). This causes the Cauchy stress tensor in (1.2) to take the form (cid:18) ∂W ∂I1 (cid:19) T = 2 v + I1 ∂W ∂I2 B − 2 v ∂W ∂I2 B2 − pI (1.8) where B is the left Cauchy-Green deformation tensor B = FFT . The tensor B is symmetric positive definite and so has positive eigenvalues, say λ2 3. Then λ1 ≥ 0, λ2 ≥ 0, λ3 ≥ 0 are the principal stretches as measured from the unswollen 2 and λ2 1, λ2 reference configuration. The tensors B and C have the same eigenvalues and it there- fore follows from (1.7) that I1 = λ2 1 + λ2 2 + λ2 3, I2 = λ2 1λ2 2 + λ2 2λ2 3 + λ2 3λ2 1 and λ1λ2λ3 = v. For a general isotropic hyperelastic material subject to the swelling constraint (1.1), the stored energy density can equivalently be taken as a function of these stretches. In this case we shall write W = ¯W (λ1, λ2, λ3, v). The Cauchy stress tensor T is symmetric and its eigenvalues are the principle stresses, denoted by T1, T2 and T3. In the absence of swelling a standard requirement on any isotropic hyperelastic constitutive model is the well known Baker-Ericksen 12 inequality (Ti − Tj)(λi − λj) > 0, (i (cid:54)= j, λi (cid:54)= λj, no sum). (1.9) This is a requirement that the maximum (minimum) principle stress direction cor- relates with the maximum (minimum) stretch direction. For a swellable isotropic material the same logic continues to apply and so we presume that (1.9) holds for the materials under consideration here. (cid:18) ∂ ˆW ∂ ˆI1 T = 2 (cid:19) The corresponding result that follows from (1.5) is that ˆW = ˆW ( ˆI1, ˆI2, v) with + ˆI1 ∂ ˆW ∂ ˆI2 ˆB − 2 ∂ ˆW ∂ ˆI2 ˆB2 − pI (1.10) and ˆB = ˆF ˆFT = v−2/3B, ˆI1 = v−2/3I1, ˆI2 = v−4/3I2. Note in particular, that there is no obvious concept here of needing to place a ˆ over the T since the Cauchy traction is uniquely defined irrespective of reference configuration. As will be discussed in what follows, the use of the unswollen reference configura- tion is often preferred when dealing with complex geometry, loads or swelling fields, which argues for a development based on (1.8). On the other hand, use of (1.10) corresponds more directly to the classical incompressible theory (because det ˆF = 1), with its familiar hyperelastic models (corresponding now to ˆW ). Such familiar models can be pulled back to the unswollen reference configuration using W = v ˆW . Thus for a familiar form ˆW = ˆW ( ˆI1, ˆI2, v) it becomes natural to consider W (I1, I2, v) = v ˆW (v−2/3I1, v−4/3I2, v), (1.11) where also any material parameters in the familiar ˆW are allowed to become de- pendent upon v. For example, the conventional Mooney-Rivlin form ˆW ( ˆI1, ˆI2) = a1( ˆI1 − 3) + a2( ˆI2 − 3) with non-negative a1 and a2 becomes, under this transforma- 13 tion, W (I1, I2, v) = v a1(v) (cid:18) I1 (cid:19) v2/3 − 3 (cid:18) I2 (cid:19) v4/3 − 3 + v a2(v) (1.12) where the notation a1(v) and a2(v) is because these parameters can now be depen- dent upon v. Constitutive laws similar to (1.12) have previously been used to study a variety of boundary value problems involving mass addition and volumetric change. This includes the studies of Pence and Tsai on swelling induced cavity formation in tubes Pence and Tsai (2005a) and spheres Pence and Tsai (2006). Amar and Goriely Ben Amar and Goriely (2005) make use of a constant material parameter version of (1.12) in the context of more generalized processes of growth to investigate insta- bilities in the inflation response of spherical shells when the shell wall experiences anisotropic growth. The study Demirkoparan and Pence (2008) analyzes swelling in- duced twist in fiber reinforced composites for materials where the matrix constituent swells (and is described by (1.12)) but the fibrous constituent does not swell and so admits to an alternative constitutive law. In the current and next chapter we will ultimately make use of (1.12) to study the interaction of swelling and a mechanical response when the material parameters d2 and d2 are also swelling dependent. Specif- ically, this model is considered in the next chapter and in Zamani and Pence (2017) for the case a1(v) = c1vm0 and a2(v) = c2vn0 with fixed values c1 ≥ 0, c2 ≥ 0, m0 and n0. The mathematical equivalent to the constitutive law (1.12) can be expressed in the form W (I1, I2, v) = 1 2 µoαvm (cid:18) I1 (cid:19) v2/3 − 3 (cid:18) I2 (cid:19) v4/3 − 3 , µo(1 − α)vn + 1 2 with µo > 0 and 0 ≤ α ≤ 1. Here, m and n are the material parameters and identified empirically. In this way µo becomes the infinitesimal shear modulus in the unswollen 14 reference configuration. The focus in the next chapter will be on the influence of m and n in various instability phenomena in balloons and pressurized spherical shells as such a material swells. The pure I1 special case is found by taking α = 1, yielding W (I1, v) = 1 2µov(m−2/3)(I1 − 3v2/3), T = µov(m−5/3)B − pI. (1.14) The no-swelling case (v = 1) then retrieves a neo-Hookean treatment. For later purposes of numerical demonstration we shall consider m = 2/3 in (1.14), giving W (I1, v) = 1 2µo(I1 − 3v2/3), T = B − pI. µo v (1.15) The particular constitutive model that will be used for specific examples in next chapter is motivated by the well known Mooney-Rivlin model WMR(I1, I2) = d1(I1 − 3) + d2(I2 − 3), (1.16) in the classical incompressible theory where the positive constants d1 and d2 are em- pirically determined material parameters. The generalization of the Mooney-Rivlin model for swelling is to keep the basic form (1.16) while now letting d1 and d2 de- pend upon v. For this purpose we shall in what follows consider examples using the constitutive model (1.12) W (I1, I2, v) = d1 (cid:19) (cid:18) I1 v2/3 − 3 (cid:18) I2 (cid:19) v4/3 − 3 , d1 = va1(v), d2 = va2(v), + d2 (1.17) with d1 ≥ 0, d2 ≥ 0 and d1 + d2 > 0. Recall that the reason for the scaling I1/v2/3 and I2/v4/3 in (1.17) is that I1/v2/3 = I2/v4/3 = 3 for an equiaxial free expansion F = v1/3I. This enables certain algebraic simplifications. When (1.17) holds the 15 Cauchy stress tensor (1.2) becomes (cid:19) T = 2 v5/3 + I1 d2 v7/3 B − 2 (cid:18) d1 (cid:18) d1 Calculating Ti and Tj from (1.18) one then obtains Ti = 2 v5/3 + (λ2 j + λ2 k) d2 v7/3 d2 v7/3 B2 − pI, (cid:19) i − p, λ2 B = FFT . (1.18) (i (cid:54)= j (cid:54)= k (cid:54)= i). (1.19) Using this result it follows for the material model (1.17) that (cid:18) d1 (cid:19) (Ti − Tj)(λi − λj) = 2 v5/3 + λ2 k d2 v7/3 (λi − λj)2(λi + λj) (1.20) Thus the Baker-Ericksen type condition (1.9) is automatically satisfied when W is given by (1.17) because d1 ≥ 0 and d2 ≥ 0. 16 CHAPTER II Burst Instability in Uniform Swelling of Hyperelastic Spherical Shells In chapters II and III we consider spherical symmetric deformation in a general Mooney-Rivlin type material in which swelling is taken into account. The types of Mooney-Rivlin materials that are distinguished with respect to the response to applied internal and external pressure (inflation) is central to this study. In the clas- sical hyperelastic setting (no swelling) a basic classification of materials with regard to their qualitative pressure-inflation behavior is given in Carroll (1987). Pence and Tsai (2006) generalize Carroll’s procedure so as to account for swelling and show that when swelling is uniform the basic qualitative response is preserved. I briefly review those results and then make use of the literature as a base to generalize the classifi- cation for the material types that were given by Carroll. Consequently, in this thesis, the more general type of materials are studied for which even homogeneous swelling field can cause instability to the type of behavior that they exhibit in response to the pressure-swelling interaction. Then in chapter III the effects of a non-homogeneous swelling field are investigated where it is shown that the precise way in which an over- all fixed amount of swelling is distributed spatially can alter the quantitative response while preserving key features of the qualitative response. 17 2.1 Kinematics for radial inflation of a spherical shell Using the framework just described in Section 1.1 we consider a finite thickness spher- ical shell with inner radius Ri > 0 and outer radius Ro > Ri prior to any loading or any swelling. Attention is restricted to radially symmetric swelling v = v(R). The loading is taken to consist of applied pressures Pi and Po on the inner and outer boundaries. These symmetric conditions motivate the consideration of the symmetric deformation for radial inflation r = r(R), θ = Θ, φ = Φ (2.1) on Ri ≤ R ≤ Ro, 0 ≤ Θ < 2π, 0 ≤ Φ ≤ π where the radial inflation function r(R) is to be determined. Thus (2.1) is a map from reference spherical coordinates (R, Θ, Φ) to deformed spherical coordinates (r, θ, φ). Let {eR, eΘ, eΦ} represent the unit basis vectors in the reference configuration and let {er, eθ, eφ} represent the unit basis vectors in the deformed configuration. It follows from (2.1) that the deformation gradient is given by F = λr(er ⊗ eR) + λθeθ ⊗ eΘ + λφeφ ⊗ eΦ (2.2) with λr = r (cid:48) and λθ = λφ = r/R in which prime (cid:48) denotes the differentiation with (cid:48) respect to R (r = dr/dR). Hence C = FT F = λ2 r(eR ⊗ eR) + λ2 θ(eΘ ⊗ eΘ + eΦ ⊗ eΦ), B = FFT = λ2 r(er ⊗ er) + λ2 θ(eθ ⊗ eθ + eφ ⊗ eφ), (2.3) (2.4) and we observe that λr and λθ = λφ are the principal stretches. in which prime (cid:48) denotes the differentiation with respect to R (r (cid:48) = dr/dR). The 18 swelling constraint (1.1) is v = (cid:48) r2r R2 , Integrating (2.5) from the inner radius Ri to a generic radial value R gives R(cid:90) r3 = r3 i + 3 v(ξ)ξ2dξ (2.5) (2.6) Ri where ri = r(Ri). More generally (2.6) provides the map r = r(R) in terms of the single parameter ri which still needs to be determined. The Cauchy stress tensor takes the form with T = Trr(er ⊗ er) + Tθθ(eθ ⊗ eθ + eφ ⊗ eφ) Trr = Tθθ = 2 v 2 v ∂W ∂I1 ∂W ∂I1 λ2 r + λ2 θ + 4 v 2 v ∂W ∂I2 ∂W ∂I2 θ − p, λ2 rλ2 (λ2 r + λ2 θ)λ2 θ − p. (2.7) (2.8) The equilibrium equations divT = 0 gives that p = p(R) along with the requirement dTrr dr + 2 r (Trr − Tθθ) = 0. (2.9) The specified pressures Pi and Po at the inner and outer surfaces yield the boundary conditions where ro = r(Ro). Trr (cid:12)(cid:12)ri (cid:12)(cid:12)ro = −Pi, Trr = −Po, (2.10) Using (2.5) and (2.8) the equilibrium equation (2.9) provides an ordinary differen- tial equation for p(R) which can be integrated. The two boundary conditions (2.10) determine the integration constant that emerges from this integration as well as the 19 scalar ri. Once ri is so determined the whole kinematics r = r(R) then follows from (2.6) and consequently the deformation gradient tensor F is fully known. This consti- tutes the most obvious solution procedure for determining the output response r(R) as a function of the input swelling field v(R) and the input pressure values Pi and Po. There is however a shortcut that gets to the same result by making direct use of the stored energy density in terms of the principle stretches ¯W (λr, λθ, λφ, v). It is based on a straight forward adaptation of a well known procedure from hyperelasticity when no swelling is present. This earlier procedure corresponds to the special case v ≡ 1 in the present treatment. Because W (I1, I2, v) = ¯W (λr, λθ, λφ, v) with λθ = λφ one may employ the chain rule for the differentiation in (2.8) in the form ∂W ∂I1 = ∂ ¯W ∂λr ∂λr ∂I1 + 2 ∂ ¯W ∂λθ ∂λθ ∂I1 , with ∂λr ∂I1 = r + λ2 λ2 θ r − λ2 θ) 2λr(λ2 , ∂λθ ∂I1 = λθ θ − λ2 r) 2(λ2 . (2.11) (2.12) A similar differentiation applies with respect to I2 and on this basis one confirms that Trr = λr v ∂ ¯W ∂λr − p, Tθθ = Tφφ = λθ v ∂ ¯W ∂λθ − p. Introduce the variable s = r/R, (2.13) (2.14) which, in view of (2.2), is the biaxial stretch λθ = λφ. Also let si = ri/Ri and so = ro/Ro and note that (2.6) gives (cid:18) Ri (cid:19)3 Ro(cid:90) Ri s3 i + 3 R3 o v(R)R2dR. (2.15) s3 o = Ro In this regard, the following derivation will be useful for the formulations in the next 20 section ds dr = ds/dR dr/dR = Rr (cid:48) − r R2r (cid:48) = v − s3 Rv , (2.16) where the last step from the above chain of equations is due to (2.5). The relation (2.16) together with (2.5) gives (cid:48) λr = r = v/s2, λθ = λφ = r/R = s, (2.17) and hence ¯W (λr, λθ, λφ, v) = ¯W (v/s2, s, s, v). Now define w(s, v) = ¯W (v/s2, s, s, v) (2.18) and calculate the derivative ∂w(s, v) ∂s = ∂ ¯W ∂λr ∂λr ∂s + 2 ∂ ¯W ∂λθ ∂λθ ∂s = −2v s3 ∂ ¯W ∂λr + 2 ∂ ¯W ∂λθ . (2.19) The combination of (2.17), (2.19) and (2.13) gives Trr − Tθθ = −(s/2v)∂w(s, v)/∂s so that the equilibrium equation (2.9) can be expressed as dTrr dr = s rv ∂w(s, v) ∂s . so(cid:90) (2.20) (2.21) This is now integrated with the aid of (2.10) and (2.16) to yield ∆P ≡ Pi − Po = 1 v − s3 ∂w(s, v) ∂s ds. si Because so is determined by si from (2.15) it follows that ∆P from (2.21) is indeed a function of si. In the absence of swelling, meaning that v = 1 identically, (2.21) retrieves a standard expression from conventional incompressible, isotropic hypere- lasticity (see (7.18) of Green and Shield (1950) and (5.3.21) of Ogden (1997)). The 21 above swelling generalization is equivalent to that given in (22) of Pence and Tsai (2006). In general the swelling field v could depend on position within the shell wall, i.e., v = v(R). This would then require v to be treated as a function of s for the purpose of the integration in (2.21), say v(R) = ˆv(s). Such a treatment will be developed later in chapter III but not do so in this chapter. In this chapter we restrict atten- tion to homogeneous swelling in the shell wall. This means that v is constant as a function of R, however such a v could vary with time in a quasi-static fashion. Thus in this chapter we now restrict attention to homogeneous swelling where v is a time- dependent parameter. However, note that throughout this work, ∆P and v serve as independent control variables in both cases of uniform and varying swelling fields. It now follows from (2.6) that r =(cid:0)r3 i )(cid:1)1/3 . i + v(R3 − R3 For the case of material model (1.17), one obtains that (cid:18) v2/s4 + 2s2 v2/3 (cid:19) − 3 d1 + (cid:18) s4 + 2v2/s2 v4/3 w(s, v) = (2.22) (cid:19) − 3 d2 (2.23) and this in turn puts (2.21) in the form so(cid:90) ∆P = 4 v − s3 (cid:20)(cid:18) s v2/3 − v4/3 s5 (cid:19) d1 + (cid:18) s3 v4/3 − v2/3 s3 (cid:19) (cid:21) d2 ds. (2.24) si Note that (2.24) continues to allow for the possibility of d1 = d1(v) and d2 = d2(v). Equation (2.24) in conjunction with i )(cid:1)1/3 (2.25) (cid:0)R3 so = 1 Ro i s3 i + v(R3 o − R3 22 provides the general relation between: amount of swelling v, applied pressure ∆P , and inner radius expansion ri = si Ri for a material with swelling dependent stored energy density (1.17). Although the integration associated with (2.24) could certainly be performed, the resulting lengthy analytical expression does not provide much in- sight. Instead, we obtain results by generalizing ideas put forward by Carroll in the context of the incompressible theory (in which there is no swelling concept). 23 2.2 The role of swelling in the pressure-inflation relation In the absence of swelling, the inflation of a pressurized spherical shell is a classi- cal problem that has been widely studied within the theory of incompressible finite hyperelasticity, i.e., with the constraint detF = 1. As first discussed in detail by Green and Shield (1950), a radially symmetric spherical inflation is possible in every isotropic homogeneous incompressible hyperelastic material. Ultimately, it is given by the v = 1 version of (2.21) and (2.23). This permits the construction of an inflation graph, which is a plot of ∆P as a function of ri. A basic discussion on different qual- itative forms for the inflation graph is given by Carroll (1987). This in turn allows one to identify different material classes. As we now show, these concepts readily generalize so as to provide similarly useful organizing concepts when swelling takes place. 2.2.1 Uniform expansion occurs for homogeneous swelling in the absence of pressure While the discussion in Carroll (1987) made no reference to the swelling concept, the concept is easily introduced into the treatment. Namely, there is now a separate inflation graph for each value of v. As v is increased continuously, it generates a family of inflation graphs in a continuous fashion. We consider the basic features of this family of graphs for a material obeying the Baker-Ericksen type condition (1.9). Using (2.17) this condition, henceforth referred to as the B-E condition, becomes (Trr − Tθθ)(λr − λθ) = s3 − v 2v s ∂w ∂s = −(v − s3)2 (cid:125) (cid:123)(cid:122) (cid:124) 2vs ≤ 0 1 v − s3 (cid:124) (cid:123)(cid:122) ∂w ∂s (cid:125) integrand of (2.21) . (2.26) Thus the B-E condition (1.9) gives that the integrand in (2.21) is negative at all locations where λr (cid:54)= λθ. What happens if λr = λθ? Isolated locations where λr and λθ coincide have no effect on the overall integral. This leaves a case in which λr and 24 λθ coincide on some interval in s. In that case r (cid:48) = r/R = v1/3 on that interval. Then, because s = r/R = v1/3 which is a single value, the interval is in fact just a single point. We conclude that (1.9) ensures that the integrand of (2.21) is negative except at possible isolated locations that do not effect the evaluation of the integral. It is useful to remark upon the special case for which λr = λθ for all R within the shell wall. This corresponds to r (cid:48) = r/R = v1/3 for all R, i.e., r3 = vR3 throughout the spherical shell. This represents a uniform expansion. For a uniform expansion it follows that s = r/R = v1/3 for all R so that so = si = v1/3. Consequently the limits of integration in (2.21) are identical and thus ∆P = 0. This gives the baseline result: • In the absence of pressurization (∆P = 0), homogeneous swelling causes the sphere to undergo uniform expansion (r3 = vR3). This gives λθ = λr and Trr = Tθθ for all R. Such a homogeneous swelling expansion is represented by the point (si, ∆P ) = (v1/3, 0) on the inflation graph. Suppose now that r3 > vR3, i.e., an amount of inflation that exceeds uniform expansion. This will be the case if ri > v1/3Ri. Then s3 = r3/R3 > v and it follows from (2.16) that ds/dr < 0 and hence so < si. Thus the integral in (2.21), because the integrand is negative, gives ∆P > 0. Furthermore λθ = r/R > v1/3, λr = r(cid:48) = vR2/r2 < v1/3 and hence Tθθ > Trr for all R. Consequently one obtains another useful result, the first part of which is also intuitive: • A positive pressurization ∆P > 0 causes the inflation to exceed that of the uniform expansion due to the swelling alone. Then λr < v1/3 < λθ and Trr < Tθθ for all R. In a similar fashion, it follows that ∆P < 0 gives an inflation that is less than that of a swelling induced uniform expansion. In such a case one might also expect various 25 wrinkling type instabilities that break the spherical symmetry. For this reason we restrict attention to ∆P ≥ 0. For the same reason we also avoid the consideration of de-swelling (0 < v < 1). 2.2.2 Qualitative behavior in the absence of swelling For a specified value of v the form of the pressure-inflation graph is determined by the stored energy density W using (2.21). The resulting relation between ∆P and si = ri/Ri is dependent on the shell thickness. This shell thickness will be characterized by the thickness ratio parameter ξ def= Ri/Ro, 0 < ξ < 1. (2.27) The thin shell limit is then ξ → 1. The other limit of ξ → 0 can be viewed either as a microvoid in a finite body or as a spherical hole in an infinite body. If v = 1, i.e. no swelling at all, then we are in the domain of conventional incom- pressible hyperelasticity and the problem reduces to one that has been extensively studied. In this conventional incompressible hyperelastic context, Carroll (1987) iden- tifies three different types of behavior which he names type (a), type (b) and type (c). These three types are diagrammed in Figure 2.1 and are described as follows: • type (a) behavior: ∆P increases monotonically with increasing ri; • type (b) behavior: ∆P increases to a maximum value and then decreases to a nonnegative asymptotic value; • type (c) behavior: ∆P first increases to a local maximum, then decreases to a local positive minimum before again monotonically increasing. In the hyperelastic theory without swelling, certain stored energy forms W always give an inflation graph with type (a) behavior. The Mooney-Rivlin material (1.16) specialized to d1 = 0 and d2 > 0 is such a material. 26 Other stored energy forms always give an inflation graph with type (b) behavior; the neo-Hookean material, meaning a Mooney-Rivlin material (1.16) with d2 = 0 and d1 > 0, is such a material. Finally there are certain stored energy forms that give type (a) behavior if ξ is relatively small but give type (c) behavior if ξ is relatively large (close to one). The Mooney-Rivlin material (1.16) with d1 = 10d2 is such a material. For these materials there is a transitional value of the thickness ratio Ri/Ro = ξ, say ξa/c, such that ξ < ξa/c implies type (a) behavior and ξ > ξa/c implies type (c) behavior. Alterna- tively stated, these materials have a type (c) inflation graph for thin shell geometries but have a type (a) inflation graph for thick shell geometries. While the above clas- sification framework was established by Carroll (1987) for the incompressible theory (v ≡ 1) we now use it to describe the swelling materials under consideration. Figure 2.1 Inflation graphs showing three qualitatively different types of behavior (a)-(c) in the absence of swelling. These particular graphs correspond to W given by (1.16), all with thickness ratio ξ = 0.5. The differences are due to the values of d1 and d2. Here: d1 = 4d2 (top); d1 = 9d2 (middle); and d2 = 0 (bottom). 27 type (c) type (a)type (b) local minimumlocal maximumsi=1; DP=0246810-0.50.00.51.0si=riRiDPΜ 2.2.3 Quantitative determination of the inflation behavior type For any fixed value v > 1 the inflation graph will continue to display one of the various behaviors shown in Figure 2.1. However changes in v could cause a transition from one behavior type to another. For this reason it is useful to obtain a more quantitative characterization of the conditions that distinguish the different graph behaviors. The presence of either a local maximum or a local minimum in the inflation graph is dependent on whether the derivative d dsi (∆P ) vanishes for some value of si. It follows from (2.21) that this derivative is given by (cid:18) 1 s2 o (cid:12)(cid:12)(cid:12)(cid:12)so ∂w ∂s (cid:19) , (cid:12)(cid:12)(cid:12)(cid:12)si − 1 s2 i ∂w ∂s d dsi (∆P ) = s2 i v − s3 i where use has been made of (2.6) and the connections o = s3 s3 i ξ3 + v(1 − ξ3), dso dsi = ξ3 s2 i s2 o . (2.28) (2.29) Equation (2.28) shows that the inflation graph will have a zero slope location only if the following condition is met: d dsi (∆P ) = 0 ⇔ 1 s2 o ∂w ∂s (cid:12)(cid:12)(cid:12)(cid:12)so (cid:12)(cid:12)(cid:12)(cid:12)si . = 1 s2 i ∂w ∂s (2.30) To make use of this condition let η be the similarity variable v/s3. Because we restrict attention to ∆P ≥ 0 it then follows that r ≥ v1/3R and hence 0 < η = v/s3 ≤ 1. Next define the auxiliary function G(η, v) def= 1 2v1/3η2/3 ∂w(s, v) ∂s 28 (cid:12)(cid:12)(cid:12)(cid:12)s=(η/v) . −1/3 (2.31) (2.32) Using (2.32) it follows that condition (2.30) is equivalently expressed as d dsi (∆P ) = 0 ⇔ G(ηi, v) = G(ηo, v), (2.33) with ηi = v/s3 i , ηo = v/s3 o. (2.34) Conditions (2.13), (2.19) and (2.32) enable a physical interpretation for the function G in terms of swelling v, stretch s = λθ, and stresses Tθθ, Trr, namely (cid:12)(cid:12)(cid:12)η=vR3/r3 G(η, v) (Tθθ − Trr). = v2 λ3 θ (2.35) Now for any fixed value v, the development in Section (2.2.2) showed that a uniform expansion takes place if ∆P = 0. This means that s = r/R = v1/3 for all R so that, in particular, so = si = v1/3 and hence ηi = ηo = 1. Because all of the principle stretches are then coincident one also obtains from (2.13) that Tθθ = Trr. Consequently (2.35) indicates that G(1, v) = 0. (2.36) Conversely, ∆P > 0 gives si > so > v1/3 which in turn implies s3 i > s3 o > v and hence 0 < v/s3 i < v/s3 o < 1. It follows that the first arguments of G in (2.33) are ordered 0 < ηi < ηo < 1 when ∆P > 0. (2.37) Also in this case the B-E inequality (1.9) gives Tθθ > 0 > Trr at each location R of the nonuniform spherical expansion. Hence (2.35) yields G(η, v) > 0 for 0 < η < 1. (2.38) 29 Figure 2.2 Graphs for G(η, v) corresponding to the three inflation curves in Figure 2.1. The function G(η, v) is computed on the basis of (2.32) using (1.17) and taking v = 1. This is equivalent to using (1.16) and ultimately gives the expression (2.43) that we examine in more detail later. Figure 2.2 shows graphs for G(η, v) corresponding to each of the three inflation curves displayed in Figure 2.1 computed on the basis of (2.32) taking v = 1. The three graphs are ordered from top to bottom in the same way as the graphs in Figure 2.1. Note that each of the three graphs in Fig. 2.2 obey both of the conditions (2.36) and (2.38). The two top curves tend to ∞ as η → 0 while the bottom curve goes to zero at η = 0. The top curve is monotone decreasing, the middle curve is decreasing- increasing-decreasing, and the bottom curve is increasing-decreasing. As we discuss next, these behaviors correlate with the type (a), type (c) and type (b) behaviors exhibited in Figure 2.1. To make direct contact with condition (2.33) we define the function that maps ηi to ηo for the specific thickness ratio ξ under consideration. In view of (2.29) this function is ˆηo(ηi, ξ) def= ηi ξ3 + ηi(1 − ξ3) . (2.39) 30 d2=19Ξ=0.5d1=1d2=0d2=140.00.20.40.60.81.00.00.51.01.52.02.53.03.5ΗG The function ˆηo is now used to define the composite function H def= {G ◦ ˆηo}, i.e., H(η, v, ξ) def= G(cid:0)ˆηo(η, ξ), v(cid:1). (2.40) We now have two functions: G(η, v) defined in (2.32) and H(η, v, ξ) defined in (2.40). In general these are different functions of their first argument η. An exception occurs when ξ = 1. This is because (2.39) gives ˆηo(ηi, 1) = ηi which in turn provides G(ηi, v) = H(ηi, v, 1). The stationary value characterization (2.33) is now expressed as d dsi (∆P ) = 0 ⇔ G(ηi, v) = H(ηi, v, ξ). (2.41) We seek to determine under what circumstances, namely for what values (ηi, v, ξ), the condition (2.41) is met. For this reason we now, for the rest of this section, use ηi for the first argument of both G and H. The previous result (2.38) establishes that G(ηi, v) is a strictly positive function of ηi on 0 < ηi < 1. The function H(ηi, v, ξ) is similarly strictly positive on 0 < ηi < 1. Also (2.36) gives that H(1, v, ξ) = G(1, v) = 0. Fix the value v and consider the graphs of G(ηi, v) and H(ηi, v, ξ) as a function of ηi on the interval 0 < ηi ≤ 1 for different values of ξ. Because G(1, v) = H(1, v, ξ) = 0 these graphs meet at the endpoint ηi = 1. However because of inequality (2.37): The graph of H is shifted to the left of the graph of G on the interval 0 < ηi < 1. The amount of this shift is nonuniform in ηi and is dependent upon ξ. In the thin shell limit ξ → 1 this shift becomes vanishingly small. Figure 2.3 shows such a leftward shift for each of the three G graphs from Figure 2.2. In particular, each of the Fig. 2.2 graphs is redisplayed as solid curve. The left shifted graphs are displayed as dashed curves of the same color. We take ξ = 0.5 because this gives the 31 thickness ratio associated with the curves from Figure 2.1. Because ξ = 0.5 is not close to one (i.e., the shell is thick) the amount of leftward shift is large and this causes the H curves to become distorted relative to the original G curves. However, what is not changed for each same color pair is the monotonicity properties: decreasing for the blue pair, decreasing-increasing-decreasing for the red pair, increasing-decreasing for the green pair. In other words the monotonicity properties of H as a function of ηi do not vary with ξ and so can be regarded as inherited from the original function G. Figure 2.3 G graphs from Fig. 2.2 (solid) along with the corresponding H graphs (dashed) for thickness ratio ξ = 0.5. Each point on a G graph is shifted to the left to give a corresponding point on the H graph. This shift is small if ξ is close to one (a thin shell). Here, because ξ is not very close to one, the nonuniform shift distorts the curves, however the basic monotonicity properties do not change. Condition (2.41) holds if and only if the graph of G intersects the graph of H somewhere on the interval 0 < ηi < 1. If such an intersection occurs, then the associated value of ηi locates either a maximum or a minimum in the corresponding inflation graph. We now consider the consequences of this observation for each of the 32 0.00.20.40.60.81.00.00.51.01.52.02.53.03.5ΗiG,H three different forms of G shown in Figure 2.2: The first and simplest G graph form is one that is monotonically decreasing as a function of ηi. In this case the graph of G cannot intersect its shifted image H. Consequently, • forms for W that generate a monotonically decreasing G as a function of ηi always give a type (a) inflation graph. This case is represented by the pair of blue curves in Figure 2.3. The solid and dashed blue curves do not intersect, and consequently the corresponding inflation graph in Figure 2.1 is monotone. The second graph form for G is one in which it is decreasing-increasing-decreasing as a function of ηi as represented by the red curve of Figure 2.2. It then follows that a small shift to the left of this graph will result in two intersections of the original graph with its shifted image. The associated inflation graph will then have a local maximum followed by a local minimum, in other words type (c) behavior. In this graphical construction the amount of shift increases as ξ decreases from ξ = 1, i.e., as the shell gets thicker. This is represented by the pair of red curves in Figure 2.3. There are two points of intersection, and these correspond to the local maximum and local minimum of the middle curve in Figure 2.1. Eventually however, the amount of shift will be sufficient to cause the shifted graph to clear itself of any intersection with the original graph. The specific shift associated with just losing this intersection involves the two intersection points com- ing together at a single special intersection point where the two graphs now have a common tangent. At this special shift, not only does (2.41) hold, but also ∂ ∂ηi G(ηi, v) = ∂ ∂ηi H(ηi, v, ξ). (2.42) This now becomes an equation for the value of ξ associated with a transition from 33 type (c) behavior to type (a) behavior. Consequently: • forms for W that generate a G function that is decreasing-increasing-decreasing on 0 < ηi < 1 give type (c) inflation behavior for thin shells and type (a) inflation behavior for thick shells. The special transition value ξ = ξa/c is found by simultaneous solution of (2.41) and (2.42) for ηi and ξ. The third and final graph form for G is one that obeys G = 0 at ηi = 0 and which then increases with ηi before decreasing back to zero at ηi = 1. In this case all left shifted curves for H will have exactly one intersection with the original G curve. Hence there will be a type (b) inflation graph for all values of ξ. Consequently: • forms for W that generate a G function that is increasing-decreasing on 0 < ηi < 1 always give a type (b) inflation graph. It is important to realize that the above inflation graph characterization in terms of G has focused on the effect of ηi and ξ irrespective of the value of v. In other words the homogeneous swelling value v has been regarded as a fixed parameter in all of the above discussion. However, for a given stored energy density W (I1, I2, v) the conclusions based on the above G graph treatment for one value of v could differ from the conclusions obtained for a different value of v. For example, the graph of G could be monotonically decreasing with ηi for values of v at and near v = 1, but could be decreasing-increasing-decreasing for relatively larger values of v. In such a case, if the shell is sufficiently thin, such a W would lead to a type (a) inflation graph when the material is unswollen but would give a type (c) inflation graph when the material is swollen. It is this issue of swelling induced changes in qualitative behavior to which we now turn our attention. 34 2.3 The swellable Mooney-Rivlin material We illustrate the effect of a changing amount of swelling using the material consti- tutive law (1.17). Thus for v = 1 this retrieves the familiar Mooney-Rivlin form (1.16). 2.3.1 Inflation behavior prior to swelling Setting v = 1 in G(η, v) gives a function that shall be denoted by g(η), i.e., g(η)def=G(η, 1). This removes v from consideration and effectively reduces the analysis procedure to that described by Carroll (1987) for conventional incompressible isotropic hyperelas- tic materials (i.e., no volume change). For the conventional Mooney-Rivlin material (1.16) this function is given by g(η) = 2d1(η1/3 − η7/3) + 2d2(η−1/3 − η5/3). (2.43) Note that g(1) = 0. Also g(η) → ∞ as η → 0 provided that d2 > 0. However if d2 = 0 then g(0) = 0. Indeed the cases used in generating the graphs in Figures 2.1 - 2.3 all corresponded to this specific example. By considering the equivalent of the derivative of g, Carroll (1987) shows how the monotonicity of this g is dependent upon the parameter ratio d2/d1. In particular, the following critical value (d2/d1)|cr def= max 0<η<1 (cid:20)η−2/3 − 7η4/3 5η2/3 + η−4/3 (cid:21) ≈ 0.215 (2.44) has special significance. Carroll shows that if d2/d1 is greater than this critical value then g(η) is monotone decreasing on 0 < η ≤ 1, but if d2/d1 is less than this critical value then g(η) is decreasing-increasing-decreasing. For our purposes it is convenient to examine the resulting consequences after 35 expressing d1 and d2 in the form which in turn gives d1 = 1 2αµ, d2 = 1 2(1 − α)µ 1 − α α , d2 d1 = α = d1 d1 + d2 . (2.45) (2.46) The reason for introducing (2.45) and (2.46) is that it makes µ > 0 the shear modulus. Indeed using (2.45) in (1.16) gives an alternative standard way of writing the Mooney- Rivlin energy form. The parameter α is in the interval 0 ≤ α ≤ 1. The special ratio of d2/d1 given in (2.44) corresponds to the critical value (cid:12)(cid:12)(cid:12)(cid:12)cr (cid:12)(cid:12)(cid:12)(cid:12)cr αcr = d1 d1 + d2 = 1 1 + (d2/d1) ≈ 0.823. (2.47) The derivative dg/dη = g (η) that is computed from (2.43) has two roots if α is in the range αcr < α < 1 and has no roots if α is in the range 0 ≤ α < αcr. This identifies (cid:48) the behavior of the inflation graph as follows: • If 0 ≤ α < αcr then the function g is monotonically decreasing with η. The inflation graph has no stationary value and so gives type (a) behavior for all ξ. • If αcr < α < 1 then the function g is decreasing-increasing-decreasing. The behavior is either type (a) or type (c) depending on whether ξ is greater or less than a transitional value ξa/c = ξa/c(α). If ξ > ξa/c then the behavior is type (c). If ξ < ξa/c then the behavior is type (a). • If α = 1 then d2 = 0 and g(0) = 0. This is the neo-Hookean special case and the function g has only one stationary value. The behavior is then type (b) for all ξ. As described in the discussion following (2.42), the transition value of ξ when αcr < 36 α < 1 can be obtained by simultaneous solution of (2.41) and (2.42). This gives a value ξa/c = ξa/c(α) for each value of α in the range αcr < α < 1 when v = 1. The curve ξ = ξa/c(α) is plotted in Figure 2.4. Any ordered pair (ξ,α) that is above the curve ξ = ξa/c corresponds to a structure (characterized by ξ) composed of a material (characterized by α) that gives an inflation graph having type (c) behavior. Conversely, ordered pairs (ξ,α) below the curve ξ = ξa/c correspond to a structure- material combination with a type (a) inflation graph. All of this follows directly from Carroll’s work in Carroll (1987). Figure 2.4 Qualitative behavior of the inflation graph for the Mooney-Rivlin model W = d1(I1 − 3) + d2(I2 − 3) as a function of material parameter α = d1/(d1 + d2) and thickness ratio ξ = Ri/Ro. The curve ξ = ξa/c provides a transition between type (c) and type (a)-behavior. 37 Αcr=0.823DPririDPtype (c) behaviortype (a) behaviorriDPtype (a) behaviorΞ=Ξac0.00.20.40.60.81.00.750.800.850.900.951.00ΞΑ 2.3.2 Inflation graph sequences for increasing swelling When swelling is present we consider the generalization of (1.16) that is given by (1.17). Then g as given by (2.43) generalizes to G(η, v) = 2(η1/3 − η7/3)d1(v) + 2(η−1/3 − η5/3)d2(v). (2.48) The direct correspondence between (2.48) and (2.43) is due to the scalings I1/v2/3 and I2/v4/3 in (1.17). This allows the analysis of G in (2.48) to proceed in a similar fashion to the previous analysis of (2.43). The main difference is that now we must account for the possible dependence of the ratio d2/d1 upon v. Figure 2.5 Inflation graphs for the Mooney-Rivlin-type model (1.17) using (2.45) with α = 0.85, and thickness ratio ξ = Ri/Ro = 0.3. All the inflation graphs exhibit type (a) behavior. We begin by considering the case where d1 and d2 are independent of v. According to (2.46) this is equivalent to α being independent of v. In this case homogeneous swelling has no effect on the type of inflation graph. Such a result is consistent with 38 v=2v=5v=3v=1v=1.4Ξ=0.3Α=0.852468100.00.51.01.52.02.53.0siDPΜ remarks given in Pence and Tsai (2006). As a first example, consider the material parameter α = 0.85. Then, because α = 0.85 > αcr, there is a transition value of ξa/c which, according to Figure 2.4, is given by ξa/c = 0.47. We now separately consider ξ = 0.3 < ξa/c (a relatively thick walled structure) and ξ = 0.7 > ξa/c (a relatively thin walled structure). The pair (ξ, α) = (0.3, 0.85) is in the type (a) behavior region of Figure 2.4, and so the v = 1 inflation graph for ξ = 0.3 displays type (a) behavior. This inflation graph is shown in Figure 2.5 along with the inflation graphs for an increasing sequence of v values. Because d2/d1 is independent of v all of the inflation graphs are monotonic. Turning to the pair (ξ, α) = (0.7, 0.85), which is in the type (c) behavior region of Figure 2.4, it follows that the v = 1 inflation graph for ξ = 0.7 displays type (c) behavior. Figure 2.6 shows this inflation graph along with those for a similarly increasing sequence of swelling values v. Because d2/d1 is again independent of v all of these graphs exhibit type (c) behavior. Figure 2.6 Inflation graphs for the same material as in Fig. 2.5 (i.e., (1.17) and (2.45) with α = 0.85), but now the thickness ratio ξ = 0.7. This corresponds to a relatively thinner walled structure. All the inflation graphs now exhibit type (c) behavior. 39 v=3v=2v=5v=1v=1.4Ξ=0.7Α=0.852468100.00.20.40.60.81.0siDPΜ 2.3.3 Swelling dependent material stiffness parameters More generally the parameters d1 and d2 may be swelling dependent, i.e., d1 = d1(v) and d2 = d2(v). It then follows that the ratio d1(v)/d2(v) changes with the amount of swelling. This can lead to the inflation graph behavior changing its type as v increases. To demonstrate consider materials for which the Mooney-Rivlin parameters d1(v) and d2(v) in (1.17) have the form d1 = 1 2µα vm and d2 = 1 2µ(1 − α) vn, (2.49) where µ > 0 and 0 ≤ α ≤ 1 are fixed material constants. This is consistent with (2.45) as can be seen by taking v = 1. Equation (2.49) introduces the additional exponent parameters m and n. The choice m = 0 and n = 0 then formally retrieves the case that was just examined in Section 2.3.2 with both d1 and d2 independent of v. For m (cid:54)= 0 and n (cid:54)= 0 the material behavior remains dependent on the ratio 1 − α α d2 d1 = vn−m. (2.50) Thus if m = n then the ratio d2/d1 is independent of v and the inflation graph behavior does not change with v. However if m (cid:54)= n then the behavior of the function G, which is now determined by d2/d1, depends on the amount of swelling v. The critical value of (d2/d1)cr = 0.215 from (2.44) continues to distinguish be- tween monotonic and non-monotonic graphs G. In this regard, for any fixed ma- terial parameter α, one may solve (2.50) for the special swelling value v that is associated with (d2/d1)cr. Define this special value of v as vA↔C. Making the re- placements d2/d1 → (d2/d1)cr and v → vA↔C in (2.50) and solving for vA↔C yields 40 vA↔C = vA↔C(m − n, α) with vA↔C(m − n, α) def= (cid:18) α 1 − α (cid:19) 1 n−m (d2/d1)cr = (0.215) n−m(cid:0) α 1 1 − α (cid:1) 1 n−m . (2.51) Now working through the various possibilities it follows that: if m < n then the graph of G is and: if m > n then the graph of G is   monotone whenever v > vA↔C, non-monotone whenever v < vA↔C, monotone whenever v < vA↔C, non-monotone whenever v > vA↔C. (2.52) (2.53) For the case of a non-monotone G graph, as discussed in Sections 2.2.2 and 2.3.1, there is a special value ξa/c of the thickness ratio ξ that gives the transition between type (a) and type (c) behavior. It is obtained by solving simultaneously the two equations (2.41) and (2.42). For given α, m and n, this value is a function of the swelling amount, hence we can write ξa/c = ξa/c(v). Such a function is directly useful if one seeks to determine the effect of a fixed amount of swelling as applied to a range of different structures, each with a different shell thickness. However the more practical problem involves a fixed structure that is subject to a changing amount of swelling. This motivates an inverting of the relation ξ = ξa/c(v) 41 to obtain v = va/c(ξ). The value of swelling va/c = va/c(ξ) demarcates the transition between type (a) behavior and type (c) behavior for the given shell geometry. From the material perspective, va/c(ξ) will depend on α, m and n. In fact, like vA↔C as given in (2.51), the dependence of va/c upon m and n will be in terms of m − n, i.e., va/c = va/c(ξ, m − n, α). However, unlike vA↔C which is independent of ξ and given by the simple form (2.51), the function va/c is dependent upon ξ and not given by a similarly simple expression. In fact the connection between these two is that va/c(ξ, m − n, α) = vA↔C(m − n, α). (2.54) (cid:12)(cid:12)(cid:12)(cid:12)ξ=1 The qualitative form of the curves va/c as a function of ξ depends on whether m > n or m < n. This is thoroughly discussed in the next section 2.3.4. In particular, it explains why (2.54) holds, and how this leads to the conclusions (2.52) and (2.53). This allows a detailed accounting for how the inflation graph varies with v beginning from the originally unswollen value v = 1 and then predicting if and when the inflation graph changes its behavior type as v increases. 2.3.4 Effect of the constitutive exponents m and n on the transitional swelling value va/c The inflation graph of ∆P vs. si for the Mooney-Rivlin swelling model that combines (1.17) with (2.49) displays either type (a) or type (c) behavior depending on the thickness ratio ξ = Ri/Ro and the swelling value v. If for fixed ξ it is possible that v alone can cause such a transition, then this transition happens when v = va/c. The transition value va/c is sensitive to the constitutive parameters α, m and n in (2.49), however it is not sensitive to µ. In (2.54) it is stated that the connection between the function va/c and the function 42 vA↔C is va/c(ξ, m − n, α) (cid:12)(cid:12)(cid:12)(cid:12)ξ=1 = vA↔C(m − n, α). (2.54) This can be understood as follows: type (c) behavior is associated with a graph for G that is not monotonic. Any transition from a monotonic graph to a nonmonotonic graph for G must take place at a value of v for which the graph develops an inflection point with zero slope. The condition for this determines vA↔C. On the other hand for a finite thickness shell the condition that determines va/c is the simultaneous solution of (2.41) and (2.42). The conditions (2.41) and (2.42) depend on the thickness ratio ξ because this dictates the amount that the graph of G shifts to the left in order to generate the H graph. This shift becomes vanishingly small in the thin shell limit ξ → 1. In order for the match condition (2.41) to hold under a vanishingly small shift it is required that any such location is one at which the graph of G has zero slope. Similarly, for the matching slope condition (2.42) to hold under a vanishingly small shift requires a zero curvature location. A location with both zero slope and zero curvature is the defining condition for vA↔C. Consequently, vA↔C is the same as va/c in the thin shell limit ξ = 1. For a finite thickness shell (ξ < 1) the values of vA↔C and va/c will no longer be the same. Here it is useful to recall the diagram in Figure 2.4 which, for v = 1, served to determine the specific thickness ratio ξ associated with the (a) to (c) behavior transition for values of α that were in the special range permitting both behaviors. When swelling is present any such transition is sensitive to both ξ and v. It is then useful to construct curves of va/c as a function of the structural parameter ξ for fixed material parameters α, m, n. Given a particular shell geometry constructed of a specific material, one can then locate the appropriate point on such a va/c curve for the purpose of determining the transitional swelling value. The form of these curves are qualitatively different depending on whether m > n or m < n. We now describe in more detail these two separate cases. 43 2.3.4.1 Dependence of va/c on ξ for m > n If m > n then the Mooney-Rivlin swelling model (1.17) with (2.49) gives va/c > vA↔C. The sphere exhibits type (a) behavior for v < va/c and type (c) behavior for v > va/c. For fixed constitutive parameters α, m, n the difference |va/c − vA↔C| decreases for relatively thinner shells (i.e., as ξ = Ri/Ro increases). In particular, va/c → vA↔C as ξ → 1. These features are apparent in Figure 2.7 which plots the dependence of va/c upon ξ for exponent choices m = 2/3 and n = 0. Because va/c depends on m and n only through their difference, the Fig. 2.7 plots apply more generally to any m and n values obeying m − n = 2/3. The different curves correspond to different values of α. Each curve is monotonically decreasing from infinity (as ξ → 0) to the value of vA↔C at ξ = 1. Curves for values of α > 0.823 are everywhere above the line va/c = 1. This is because vA↔C > 1 when α > 0.823. In contrast, because α < 0.823 makes vA↔C < 1 it follows that the curves for α < 0.823 cut the line v = 1. Because we limit attention to v ≥ 1 the portions that continue into v < 1 are shown as dashed. The α value of 0.823 is the value of αcr that was first introduced in (2.47) in the context of the standard neo-Hookean model. By virtue of (2.54) it also serves to make vA↔C(m− n, αcr) = 1 because of the direct way in which the standard incompressible model ((1.16) with (2.45)) was generalized to the swelling model ((1.17) with (2.49)). The curves shown in Figure 2.7 correspond to m − n = 2/3. Curves with similar qualitative behavior are obtained provided that m > n. In particular, the α value of 0.823 is always associated with vA↔C = 1. Spherical shells with m > n and α > 0.823 have va/c > 1. They exhibit type (a) behavior for 1 ≤ v < va/c and exhibit type (c) behavior for for v > va/c. 44 Figure 2.7 Transitional swelling value va/c versus ξ = Ri/Ro for the Mooney-Rivlin- type model (1.17) with parameters µ, α, m and n in (2.49). The transitional swelling value va/c is independent of µ and is dependent on n and m only via the difference m − n. These plots are for m − n = 2/3. For a given α-curve the inflation graph exhibits type (c) behavior if (ξ, v) is in the region above the curve and type (a) behavior if (ξ, v) is in the region below the curve. 45 Α=0.80Α=0.823Α=0.85Α=0.90Α=0.95vA«C=1.26Α=0.75vA«C=1.94Α=0.7vA«C=2.820.00.20.40.60.81.0012345Ξvac 2.3.4.2 Dependence of va/c on ξ for m < n One may similarly construct curves of va/c vs. ξ for the case in which m < n. Now the curves are increasing with ξ instead of decreasing. Each curve continues to approach the value vA↔C as ξ → 1, however now they increase from the value zero at ξ = 0. The other major difference is in the significance of these curves. Namely, the spherical shells now have a type (a) behavior in the region above the curves (v > va/c) and have a type (c) behavior in the region below the curves (v < va/c). Figure 2.8 Transitional values for swelling va/c versus ξ = Ri/Ro for Mooney-Rivlin- type model (1.17) with parameters µ, α, m and n in (2.49). The transitional swelling value va/c is dependent on α as shown but is independent of µ. The curves are dependent on m and n only via the difference m− n. This figure is for m− n = −2/3. For a given α-curve the inflation graph exhibits type (c) behavior when (ξ, v) is in the region that is below the curve and type (a) behavior in the region that is above the curve. Such curves are displayed in Figure 2.8 which plots the dependence of va/c upon ξ for exponent choices m = 0 and n = 2/3. More generally this figure also applies to any m and n values obeying m − n = −2/3. Qualitatively similar curves hold for 46 Α=0.95Α=0.9Α=0.85Α=0.98Α=0.823Α=0.4Α=0.6Α=0.8Α=0; vac=0vA«C=2.68vA«C=1.340.00.20.40.60.81.00.00.51.01.52.02.53.0Ξvac any m and n obeying m < n. The value α = 0.823 continues to retain its special significance because of its continued association with the condition vA↔C = 1. It is now the case that curves for α < 0.823 are always confined to the region va/c < 1. Because the condition v < 1 is not being considered, the α < 0.823 curves are shown as dashed over their entire length. Thus if α < 0.823 (and m < n) then any spherical shell has type (a) behavior for v ≥ 1. Conversely curves for α > 0.823 cut the line v = 1. The portions of these curves that are below the value v = 1 are again shown as dashed. If α > 0.823 and the shell is sufficiently thick then it has a type (a) inflation graph for all v ≥ 1; this is because va/c < 1. However if the shell is sufficiently thin then va/c > 1; this means that it has a type (c) inflation graph for 1 ≤ v < va/c and a type (a) inflation graph for v > va/c. Consequently in such a case a quasi-static increase in v from the unswollen state v = 1 will generate a transition from type (c) to type (a) behavior as v passes through the special value va/c. 2.3.5 Numerical illustration with swelling-dependent material stiffness parameters The inflation behavior of the sphere with swelling dependent material stiffness can be illustrated by considering the same α and ξ values associated with Figures 2.5 and 2.6 but now allowing for m (cid:54)= 0 and n (cid:54)= 0. For this purpose we first consider the case m < n that is obtained by taking m = 0 and n = 2/3. In particular, consider two subcases corresponding respectively to a thick shell (ξ = 0.3) and to a relatively thinner shell (ξ = 0.7). Thus the two subcases correspond to (ξ, α, m, n) = (0.3, 0.85, 0, 2/3) and to (ξ, α, m, n) = (0.7, 0.85, 0, 2/3). The v = 1 curve for the first subcase is identical to the v = 1 type (a) curve from Figure 2.5. Similarly, the v = 1 curve in the second subcase matches the v = 1 type (c) curve from Figure 2.6. However the curves for v > 1 will no longer match the curves shown in these respective figures. One finds 47 for the first subcase, that with ξ = 0.3, the type (a) inflation graph that is present for v = 1 persists for all increasing v. This aspect mirrors the situation in Figure 2.5 even though the individual curves for v > 1 are different. In the second subcase of ξ = 0.7 one finds that the inflation graph is originally type (c) for v = 1 but it eventually transitions to type (a) behavior as v increases. This transition occurs at v = 1.27, a result that can be predicted on the basis of the procedure for determining va/c that is described in the appendix. The case m > n can be handled similarly. For this purpose consider m = 2/3 and n = 0, again for the respective thick and thin shell values ξ = 0.3 and ξ = 0.7. Once again the v = 1 curves match the v = 1 curves from Figures 2.5 and 2.6 respectively. Once again the v > 1 curves do not match the curves in these two figures. In fact, Figures 2.9 and 2.10 show a v > 1 curve sequence for these two respective cases. For the subcase of ξ = 0.7 one finds that the original type (c) behavior at v = 1 will persist as v increases (Fig. 2.10). In contrast, for the case of ξ = 0.3 one finds that the original type (a) behavior at v = 1 will transition to type (c) behavior as v increases (Fig. 2.9). This transition occurs at v = 1.54, where, again, such a result can be understood in detail on the basis of the treatment given in the appendix. 48 Figure 2.9 Inflation graphs for the Mooney-Rivlin-type model (1.17) using (2.49) with α = 0.85, m = 2/3, n = 0, and thickness ratio ξ = 0.3. The inflation graphs exhibit the type (a) behavior for 1 ≤ v < 1.54 and type (c) behavior for v > 1.54. 2.4 Swelling induced burst Each of the previous Figures 2.5 - 2.10 shows a sequence of inflation graphs for a given shell thickness ratio ξ composed of a given model material (α, m, n). Such a figure can be used to gauge how the sphere expands as a function of increasing swelling v. If the pressurization is fixed during the swelling, then a quasi-static increase in v corresponds to moving between different curves on the same figure along the horizontal line determined by the stipulated ∆P . For continuously increasing v the associated increase in si will also be continuous so long as all of the curves in the sequence are monotonically increasing. However if the curves are not all monotone increasing then there is clearly the possibility of a discontinuity in si. For example, consider again Figure 2.6. The inflation graphs for all v are non- monotone (type (c)) and the si interval of graphical decrease varies with v. Figure 2.11 49 Figure 2.10 Inflation graphs for the Mooney-Rivlin-type model (1.17) using (2.49) with α = 0.85, m = 2/3, n = 0, and thickness ratio ξ = 0.7. The inflation graphs exhibit type (c) behavior for all v ≥ 1. identifies the specific pressurization ∆P that corresponds to the local maximum for v = 2. Its value is ∆P = 0.258µ. Starting with an unswollen and unpressurized sphere (v = 1, ∆P = 0) consider first an increase in pressure from ∆P = 0 to ∆P = 0.258µ while the sphere remains unswollen. The inflation response corresponds to climbing the v = 1 curve to ∆P = 0.258µ with a relatively small increase in si from si = 1 to si = 1.14. Now holding this pressurization fixed let v increase. Then one may proceed in se- quence through all of the curves from the original v = 1 curve to the curve for v = 2. During this sequence there is a corresponding continuous increase in si. However, increasing v beyond v = 2 cannot proceed with a continuous increase in si because the local maximum signals the onset of an interval in si corresponding to v < 2. This interval proceeds from si = 2.37 to si = 4.32. While this interval precludes a continu- ous increase in si as v increases through v = 2 it does permit a discontinuous increase 50 Figure 2.11 Inflation burst caused by increasing v at fixed ∆P = 0.258µ for the inflation graphs from Fig. 2.6. Prior to swelling the pressurization ∆P = 0.258µ has given a mild radial increase (from si = ri/Ri = 1 to si = ri/Ri = 1.14 on the v = 1 curve). Now increasing v at this fixed ∆P gives a continuous increase of si with v (dashed red line) until encountering the inflation graph for v = 2 where there is a local maximum. Further increase of v requires a jump across to the other increasing branch of the v = 2 curve (solid red segment). This corresponds to an inflation burst with radial increase from si = 2.37 to si = 4.32. from si = 2.37 to si = 4.32 at v = 2. After such a jump in si it is then again possible for si to increase continuously because the inflation graphs again become ordered so as permit si to increase with v. Figure 2.12 shows directly the corresponding radial increase with swelling (si as a function of v). The jump in si corresponds to a “burst of inflation” of limited extent (it concludes at si = 4.32). Such a burst, which can also be described as a snap-through, is due to the presence of a local maximum in the sequence of inflation graphs. This gives multi-valued choices for si when an inflation graph exhibits two increasing branches separated by a decreasing branch. 51 v=2v=5v=3v=2.5v=2.3v=1v=1.42468100.00.20.40.60.81.0siDPΜ Figure 2.12 Inflation burst showing si = ri/Ri vs. v at ∆P = 0.258µ caused by an increase in the swelling parameter v. Locations denoted by • provide correlation with the inflation graphs depicted in Figure 2.11. Under such circumstances some kind of burst is inherent in the mechanical de- scription. However, the description is potentially ambiguous as regards the value of v at which the jump occurs. For example, we have just described a jump from si = 2.37 to si = 4.32 when both ∆P = 0.258µ and v = 2. However, for ∆P = 0.258µ the inflation graphs become multi-valued in si for values of v < 2 and so the question arises, “why not jump before v = 2?”. In other words, while v = 2 is the maximum value of v that permits one to avoid a discontinuity, there is always the possibility of executing an earlier jump. Such issues have been extensively studied in conventional hyperelasticity M¨uller and Strehlow (2004) (i.e., no swelling). Then for a single type (c) inflation graph an increase in ∆P eventually provokes a jump to the second increasing branch for the simple reason that the first increasing branch has a maximum permissible ∆P 52 Inflation burst at v=2Dsi=1.95DP=0.258Μ1.01.52.02.53.012345678vsi value. This jump could occur at the local maximum or it could occur before the local maximum. Viewing such jumps as a type of phase transition it can be shown that an energy minimal quasi-static process of ∆P increase predicts that the transition occurs prior to attaining the maximum. Specifically it occurs at the value of ∆P associated with the “Maxwell line” construction Ericksen (1975). On the other hand, a transition that occurs at the local maximum upon loading (and at the local minimum upon unloading) is consistent with a notion that the prevailing phase can, under carefully controlled conditions, be preserved even though distantly related states of deformation may now lower the system free energy. In other words, if the system is not subject to large disturbances then jumps will occur at extrema of the inflation graph because it is only then that the inevitable small disturbances provoke a jump to a more energetically favorable configuration. Such considerations continue to apply to the notion of swelling induced burst that we have been describing. In particular, the sequence of inflation graphs depicted in Figure 2.6 leads to a situation where, at fixed ∆P , a continuous increase in v will give some kind of abrupt change in inflation. Whether this occurs at the local maximum of an inflation graph or whether it occurs prior to such a local maximum is then to be answered on the basis of a more refined treatment. This includes energetic stability analysis such as that described in Ericksen (1975) as well as the consideration of less symmetric deformations (such as those with new modes of localized deformation Kyri- akides and Chang (1990)). More generally, one can employ a broader thermodynamic framework that allows for supplemental physical considerations (e.g., an additional kinetic relation), as well as additional theoretical considerations from the outset (e.g., inertial dynamics, finer scale physics, a statistical physics treatment of fluctuations). Finally, it is worth remarking that the notion of pressure control is itself likely to be an idealization, and that other forms of control, such as one based on controlling a set mass of sealed in gas Alexander (1971), can lead to different predictions on how 53 transitions occur between different points on an inflation graph. The inflation burst illustrated in Figure 2.11 was based on the inflation graphs for the case that was presented in Figure 2.6. In that case all of the inflation graphs for v ≥ 1 involved type (c) behavior. Thus one could possibly argue that the swelling induced burst could have been anticipated on the basis of the original unswollen v = 1 inflation graph. However, in general it would be premature to draw conclusions on either the presence or absence of swelling induced burst just on the basis of the v = 1 inflation graph. For example, the unswollen v = 1 inflation graph in Figure 2.9 exhibits type (a) behavior. Thus if v = 1 then a continuous increase in pressure will result in a continuous expansion and so by itself provides no indication of a burst possibility. However swelling induced inflation burst can still occur. This is shown in Figure 2.13 for the example of Figure 2.9. Starting on the v = 1 inflation graph with ∆P = 1.16µ we consider a subsequent increase in v. The value ∆P = 1.16µ is chosen for this discussion because it gives the local maximum on the v = 2 inflation graph (other values could similarly be considered). Holding ∆P at this fixed value, an inflation burst is triggered at v = 2 in a manner similar to that previously depicted in the example of Figures 2.11 and 2.12. In that previous example the inflation graph behavior was type (c) for all values of ∆P prior to the burst. In the present example, the swelling induced burst involves inflation graphs that transition from “benign” type (a) graphs to “burstible” type (c) graphs as the swelling proceeds. A converse phenomena is also possible if the v = 1 unswollen graph is type (c) and which then transitions to type (a) as the swelling proceeds. This was the case for (ξ, α, p, q) = (0.7, 0.85, 0, 2/3) that was discussed in Section 2.3.3 right after equation (2.54). In such a case it is found that certain loading sequences which alternate pres- surization with strategically placed episodes of swelling and deswelling enable burst avoidance. This contrasts with the inevitability of burst if all of the pressurization 54 Figure 2.13 Inflation burst caused by increasing v at fixed ∆P = 1.16µ for the inflation graphs from Fig. 2.9. Initially, the radius increases continuously, first with v = 1 as ∆P increases from zero to 1.16µ and then at this fixed ∆P as v increases to v = 2 (dashed red segment). At v = 2 there is a jump from the first increasing branch to the second increasing branch after which a continuous increase is again the case. takes place at fixed v. 55 v=3v=2v=1v=1.4246810120.00.51.01.52.02.53.0siDPΜ CHAPTER III Non-Uniform Swelling Field of Hyperelastic Spherical Shells As it was mentioned in chapter II, the swelling field v could in general depend on position within the shell wall, i.e., v = v(R). Such a treatment is developed in this Chapter. In order to make a distinction between the varying (non-uniform) and uni- form swelling fields we use the notation vuni for the uniform swelling field in this Chapter. It was shown that for the Mooney-Rivlin-type model (1.17) while the uni- form swelling field is present and the material parameters are independent of swelling, only the constant ratio κ = d2/d1 and the thickness ratio ξ = Ri/Ro identify the clas- sification of the inflation behavior (a), (b) or (c) and the material types A, B or C. This is for example presented in Figure 3.29 wherein it is shown that different amounts of uniform swelling field value vuni has no effect on the material and behavior type. Conversely, the uniform swelling field value vuni > 1 can change the behavior of the inflation response only when the material parameters d1 and d2 are dependent on swelling amount vuni (in the form of different power laws). The question then arises as to how to determine the type of the behavior when the swelling field v is no longer distributed uniformly and varying with R. In this chapter we study the inflation behavior of such fields when the material parameters of the model (1.17) are independent of local volume change v(R). In such an event where v = v(R), it is required to invert the relation s = r(R)/R and write R = ˆR(s) and hence v = ˆv(s) for the purpose of the integration in (2.21). Additionally, in order to characterize 56 specific non-uniform swelling fields we introduce a new variable which is the overall added volume or equivalently absorbed mass due to the swelling process. This will allow us to describe family of swelling fields with the same added mass. This in turn enables us to study as how such nonuniform swelling fields will affect the inflation behavior. 57 3.1 A family of swelling fields with the same added mass To investigate the characteristics of varying swelling fields we introduce the overall added swelling volume ∆V Ro(cid:90) (v(R) − 1)R2dR, Ri ∆V = 4π (3.1) which is the total amount of volume change due to the mass absorption during the swelling. We limit our attention to consider only swelling fields that are motivated by observation of steady state distributions of liquids in porous media. For such swelling fields, the distribution of v(R) is governed by Laplace’s equation that in the reference configuration is ∇2v(R) = 1 R2 ∂ ∂R ∂R (cid:0)R2 ∂v(R) (cid:1) ≡ 0. (3.2) (3.3) (3.4) The general distribution of the swelling field that satisfies (3.2) has the form v(R) = A R + B in which A and B are constants. Let vi ≡ v(Ri) = A Ri + B, vo ≡ v(Ro) = A Ro + B where vi is the amount of swelling at R = Ri and vo is the amount of swelling at R = Ro. It then follows that A = vo − vi 1/Ro − 1/Ri , B = voRo − viRi Ro − Ri . Entering (3.1) with (3.3) one obtains ∆V R3 i = 2π Ri (ξ−2 − 1)A + 4π 3 (ξ−3 − 1)(B − 1). 58 (3.5) (3.6) Note that for any fixed ∆V and fixed geometry (Ri and ξ) in (3.6), it is seen a linear relation between the parameters A and B. Now imagine any two different mass distributions represented by (Ap,Bp) and (Aq,Bq) but with the same ∆V . The linear relation (3.6) requires that After some simplification Ap − Aq Bp − Bq = − 4π(ξ−3 − 1)/3 2π(ξ−2 − 1)/Ri . Ap − Aq Bp − Bq = −2(ξ−3 − 1) 3(ξ−2 − 1) Ri. Now define the radial location Rm as Rm := 2(ξ−3 − 1) 3(ξ−2 − 1) Ri = 2(ξ3 − 1) 3(ξ2 − 1) Ro (3.7) (3.8) (3.9) where it can be easily shown that since 0 < ξ < 1 then Ri < Rm < Ro. It then follows from rearranging (3.8) with (3.9) that Ap Rm + Bp = Aq Rm + Bq (3.10) This is in fact the amount of swelling at the radius Rm where according to (3.3) is v(Rm) for both distributions. This means that both distribution have the same amount of swelling at the same radial location Rm. Since the choice of distributions was otherwise arbitrary it is concluded that all swelling distributions v(R) given by (3.3) that absorb the same fixed added mass ∆V have the same amount of swelling at location Rm. In order to characterize the distributions (3.3) in terms of vi and vo instead of A 59 and B note that plugging (3.5) into (3.6) and solving for vo gives ξ(ξ − 1)(2ξ + 1) 2 − ξ 2(1 − ξ3) 2 − ξ . vi + + (3.11) (cid:19) (cid:18) ξ3 2 − ξ vo = 3∆V 2πR3 i It follows that (cid:0) (cid:18) 3∆V (cid:0) 2πR3 i 3∆V 2πR3 i A = B = (cid:1) + (cid:1) + vi (cid:0) (cid:0) 3ξ − 3ξ3 (cid:1) + vi ξ3 − 3ξ + 2 2ξ3 − 2 −ξ3 + 3ξ − 2 (cid:1) + 2(1 − ξ3) ξ3 − 3ξ + 2 ξ3 −ξ3 + 3ξ − 2 ξ3 ξ3 − 3ξ + 2 2(1 − ξ3) −ξ3 + 3ξ − 2 (cid:19) Ri (3.12) Thus for a given ∆V > 0, Ri, and Ro = Riξ−1 we may view (3.3) with A and B given by (3.12) as a family of swelling fields for the same added mass but different distributions with tuning parameter vi. As vi changes, the same overall added mass ∆V is distributed through the spherical shell in different ways. We restrict considerations to swelling fields (3.3) such that v(R) ≥ 1 at all loca- tions. This will be the case if both vi ≥ 1 and vo ≥ 1. Setting vo = 1 in either (3.6) or (3.11) and solving for vi we obtain (cid:12)(cid:12)vo=1 = 1 + vi 3∆V 2πR3 i (cid:18) ξ3(1 − ξ) 2ξ4 − 3ξ3 + ξ (cid:19) ≡ vmax i (3.13) Thus (3.3) with A and B given by (3.12) is parameterized by vi on the interval 1 ≤ vi ≤ vmax . i If vo = vi then A = 0 and we retrieve a uniform distribution of the kind studied in Chapter II. For a given ∆V the uniform distribution is associated with the vi value that is found by substituting vi = vo = vuni in (3.11) and solving for the special value 60 vuni. This gives (cid:0) 3∆V 4πR3 i ξ3 1 − ξ3 (cid:1). vuni = 1 + (3.14) If vi = vuni as given by (3.14) then all of the results from Chapter II apply for the given added mass ∆V . If 1 ≤ vi ≤ vuni then the added mass is more concentrated near the outer surface. If vuni ≤ vi ≤ vmax then the added mass is more concentrated i at the inner surface. Conversely, (3.14) shows that the added mass associated with this uniform distribution is ∆V R3 i = 4 3 π(vuni − 1) 1 − ξ3 ξ3 , (3.15) For any fixed amount of ∆V , the amount of vo is identified from (3.11) in terms of vi and this defines a family of swelling fields v(R) from (3.3) that can be param- eterized only by the amount vi. One example that will also be used for numerical showcase in the next section is shown in Figure 3.1. The figure shows a family of swelling distributions where they all have the same overall added mass, say 30 per- cent of the original volume and thus vuni = 1.3. In the following section the inflation behavior of such distribution families is studied. 61 Figure 3.1 The family of swelling distributions (3.3) with constants (3.12) for ξ = 0.5 that are parameterized with respect to vi such that v(R) ≥ 1. All distributions have the same amount of overall added mass associated with vuni = 1.3. 3.2 Kinematics of the spherical deformation with non-uniform swelling field The symmetric spherical deformation that is described in (2.1) is assumed here as a response to radial inflation that is also subject to the non-uniform swelling field (3.3). In fact, the deformation map (2.6) is now considered with the radial dependent v = v(R). It then follows that r3 = r3 i + 3 (cid:19) + B R(cid:90) (cid:18)A ζ Ri and this simplifies to ζ 2dζ, (3.16) r3 = r3 i + 3 2 A(R2 − R2 i ) + B(R3 − R3 i ) (3.17) 62 vi,min=1vuni=1.3vi,max=2.05vo,max=1.41vo,min=11.01.21.41.61.82.01.01.21.41.61.82.02.2RRivR wit A and B given by (3.12). 63 3.3 Behavior type independent of distribution for fixed added mass The inflation behavior for varying swelling fields is considered in this section. This makes the use of (2.21) in the form so(cid:90) ∆P = 1 ˆv(s) − s3 ∂w(s, v) ∂s ds. (3.18) si This relation can also be verified by an alternative proof with energy argument to confirm that this relation indeed holds even in the case of non-uniform swelling. The deformation field (2.1) in the absence of body forces is determined by the minimization of the total potential energy E = Estore − Eload, (3.19) where Estore is the stored energy with respect to reference configuration that includes the effect of elastic deformations and is expressed as Ro(cid:90) Estore = 4π ¯W R2dR (3.20) Ri where ¯W is the local strain energy density as defined before (2.11). The work func- tional associated with the prescribed surface tractions is denoted by Eload. Hence the work due to the internal pressure loading is given by Eload = 4 3π(r3 i − R3 i )Pi − 4 o − R3 o)Po 3π(r3 (3.21) The requirement of stationary total potential energy E with respect to the only un- 64 known ri provides (Estore − Eload) = 0. ∂ ∂ri Thus it follows that because the radial position R is independent of ri then Estore = 4π ∂ ∂ri and from (3.21) Ro(cid:90) Ri ∂ ¯W ∂ri R2dR Eload = 4πr2 i (Pi − Po). ∂ ∂ri Equating the two relations (3.23) and (5.61) one obtains ∆P ≡ Pi − Po = 1 r2 i ∂ ¯W ∂ri R2dR. Ro(cid:90) Ri By recalling (2.16) and using the connections from dR Rs2 = ds v − s3 , ∂s ∂ri = r2 i Rr2 , (3.22) (3.23) (3.24) (3.25) (3.26) we convert the relation (3.25) to (3.18). In the event that the relation s = r(R)/R is not explicitly invertible, that is R = ˆR(s) is not available, it is required to use (3.18) in the reference configuration and perform the integration with respect to R. For convenience we use this integration in this Chapter. Moreover, instead of converting (3.25) into (3.18) we use the connections (2.5) and (3.26) to rewrite (3.25) with the more convenient derivations of ¯W with respect to s that uses the definition (2.18). It 65 follows that Ro(cid:90) Ri ∆P = (cid:12)(cid:12)(cid:12)(cid:12)s=r/R dR R r2 ∂w(s, v) ∂s (3.27) Note that from (2.6) we know the map r = r(R; ri). This provides the integrand in (3.27) in terms of R and the parameter ri. This in turn makes (3.27) be indeed a relation for the applied pressure and inner deformed radius. In order to employ the inflation relation (3.27) we use the Mooney-Rivlin-type model (1.17) with (2.46). It follows that ∂w(s, v) ∂s = 2αµ (cid:18) r Rv2/3 − R5v4/3 r5 (cid:19) (cid:18) r3 R3v4/3 − R3v2/3 r3 (cid:19) + 2(1 − α)µ and thus the inflation relation (3.27) becomes (cid:12)(cid:12)(cid:12)(cid:12)s=r/R (cid:18) Ro(cid:90) α(cid:0) R r2 (3.28) (cid:1)(cid:19) dR (3.29) ∆P/µ = 2 Ri r Rv(R)2/3 − R5v(R)4/3 r5 r3 R3v(R)4/3 − R3v(R)2/3 r3 (cid:1) + (1 − α)(cid:0) where r = r(R) is from (3.17), v = v(R) is given by the swelling field (3.3) and A and B are according to (3.12). This provides the relation between the applied pres- sure ∆P and the inner radius ri. The closed form solution to this integration is not available, however the numerical integration is performed for special cases that are presented in the following. The varying swelling fields for numerical integration (3.29) will be selected such that they all have one fixed overall added mass ∆V . The choice of swelling fields is the same as was shown in Figure 3.1. As for the first set of numerical result, the other parameters in the integration are set to (ξ, α)=(0.5, 0.86). This is chosen such that, 66 Figure 3.2 Type (c) behavior in response to pressure-inflation for the six different swelling fields in figure 3.1; using fixed properties (ξ, α) = (0.5, 0.86) in the Mooney- Rivlin-type material model (1.17) and (2.46). Figure 2.4. Also note that this is associated with the added mass of (∆V /R3 with respect to vuni = 1.3, type (c) behavior is expected from the results plotted in i ) (cid:39) 8.8 given by (3.15). This set is considered in Figure 3.2 in which the solution to (3.29) for the material model (1.17) with (2.46) is plotted. As it is presented in this figure, all of the family of different distributions with fixed added mass showing type (c) behavior consistent with the type (c) behavior that is identified for the uniform swelling of this family with vuni = 1.3. For the same family of swelling distribution represented in Figure 3.1 we now consider the second set of parameters of (ξ, α) = (0.5, 0.83) such that type (a) behavior is expected for the uniform swelling vuni = 1.3 according to Figure 2.4. As it is seen in Figure 3.3 the solution shows type (a) behavior for the whole distribution family. The two above numerical examples show that the distribution family of the fixed 67 vi,min=1vi,max=2.0512345670.00.20.40.60.81.01.2siDPΜ Figure 3.3 Type (a) behavior in response to pressure-inflation for the six different swelling fields in figure 3.1; using fixed properties (ξ, α) = (0.5, 0.83) in the Mooney- Rivlin-type material model (1.17) and (2.46). added volume with the representative uniform swelling vuni = 1.3 followed the same behavior type as the uniform swelling field depicted, depending on the parameters (ξ, α). The latter case was proved in Chapter II to be identified from Figure 2.4. Our hypothesis is that the same type of behavior occurs for all (A, B)-pairs that forms a family of yielding same added volume ∆V and the type of behavior of the family is identified based on the behavior of vuni with (ξ, α). To test this hypothesis we consider 18 additional cases chosen so as to be close to transition curve shown in Figure 3.4. In each case we calculated the response for vi = 1 and vi = vi,max. The numerical results of the inflation curves with the two circled dots in the middle of this figure have been already shown in Figs. 3.2 and 3.3. The numerical results of the inflation curves are showcasing in Figs. 3.5 and 3.6 with four other circled dots from Figure 3.4 as additional examples. For any parameter-point that was chosen above the transition curve both of the representative distributions (vi = 1 and vi = vi,max) 68 vi,min=1vi,max=2.0512345670.00.20.40.60.81.01.2siDPΜ Figure 3.4 Qualitative behavior of the inflation graph for the Mooney-Rivlin model with uniform swelling fields, previously shown in Figure 2.4. Here 18 points are chosen along the transition curve (red dots) to be used in (3.29) for numerical integration. The inflation graphs with the circled parameter choices are shown in Figs. 3.2, 3.3, 3.5 and 3.6. show type (c) behavior. In contrast, for any parameters chosen below the transition curve both distributions (vi = 1 and vi = vi,max) show type (a) behavior. It is seen that at least these numerical cases support our hypothesis. 69 DPririDPtype (c) behaviortype (a) behaviorriDPtype (a) behaviorΑcr=0.823Ξtransient0.00.20.40.60.81.00.750.800.850.900.951.00ΞΑ Figure 3.5 Inflation behavior for the two limits of swelling distribution (vi,min, vi,max) in figure 3.1; using fixed parameters chosen from Figure 3.4 with (ξ, α) = (0.1, 0.94) (on the left) showing type (a) behavior and (ξ, α) = (0.1, 0.96) (on the right) showing type (c) behavior. The Mooney-Rivlin-type material model is based on (1.17) and (2.46). 3.4 Inflation behavior instability due to mass redistribution The inflation response (3.29) is numerically obtained for the swelling distributions family (3.3) and A and B given by (3.12) with fixed added mass corresponding to vuni = 1.3 (Figure 3.1), in which the family is parameterized by inner amount of swelling according to 1 ≤ vi ≤ vmax = 2.05. For the material model (1.17) with i constant parameters, the inflation graphs of this solution to (3.29) show behavior of type (c) for the set ξ = 0.5 and α = 0.86 and this is graphed in Figure 3.2. Hence the family shows type (c) behavior for the range from vi = 1, where the added mass is more concentrated on the outer layer of shell, to vi = vmax i = 2.05 where the added mass is more absorbed toward the inner layer of the shell. Thus for all of these dis- tributions the inflation graph shows type (c) behavior independent of the swelling field. Here, the possibility of an inflation instability in inner radius can be captured of the same kind that was illustrated in Figure 2.11 for the case of uniform swelling. However, in the current case the redistribution of the swelling field can cause the burst. For the case of uniform swelling with constant material parameters the quasi- static increase in swelling amount led to the movement of the inflation curves of type (c) behavior such that for a fixed applied pressure it provided the jump in the inner 70 Ξ=0.1Α=0.94vi=1vi,max=2.0551015202501234siDPΜΞ=0.1Α=0.96vi=1vi,max=2.05510152025300.51.01.52.0siDPΜ Figure 3.6 Inflation behavior for the two limits of swelling distribution (vi,min, vi,max) in figure 3.1; using fixed parameters chosen from Figure 3.4 with (ξ, α) = (0.9, 0.82) (on the left) showing type (a) behavior and (ξ, α) = (0.9, 0.83) (on the right) showing type (c) behavior. The Mooney-Rivlin-type material model is based on (1.17) and (2.46). radius. Similarly, for the case of a nonuniform swelling fields that generate a family of type (c) behavior it is the redistribution of the field that provides the movements of the inflation curves and this can causes the inflation jump. The procedure of the inflation instability is represented numerically where one aspect of Figure 3.2 is plotted in Figure 3.7. It is seen that the maximum pres- sure associated with the branch vi = 1 is ∆P/µ = 0.74 and the maximum pressure reached on the curve for vi = 2.05 is ∆P/µ = 0.67. Remaining on the branch vi = 1, the increase of pressure from ∆P = 0 to slightly higher pressures cause the radius to increase continuously. Now the pressure continues to increase and then remains fixed at some value between 0.67 < ∆P/µ < 0.74 and here we choose ∆P/µ = 0.71 for showing the example of inner radius jump. At this fixed pressure if the total added mass starts to redistribute such that it is more absorbed towards the inner layer it means that inner amount of swelling increases from vi = 1 to some amount vi > 1. This increase in vi moves the corresponding inflation graph to the right and down at the fixed pressure. This in turns moves the solution to si = ri/Ri closer to the top of the corresponding curve until it eventually reaches to the top of the 71 Ξ=0.9Α=0.82vi,max=2.05vi=1123450.000.050.100.150.20siDPΜΞ=0.9Α=0.83vi,max=2.05vi=1123450.000.050.100.150.20siDPΜ inflation graph where vi = 1.25 which corresponds to the maximum value for the pressure ∆P = 0.71. At this point spherical shell experiences a sudden expansion burst and a rapid jump in the inner radius. It follows that, at a fixed pressure, the redistribution of a fixed added mass can cause a shell expansion instability that is similar to that which was previously displayed in Figure 2.11. The important point is that now the overall amount of swelling does not change, just the way it is distributed. Figure 3.7 Inflation burst due to redistribution of the fixed added mass in type (c) be- haviors; using fixed properties (ξ, α) = (0.5, 0.86) in the Mooney-Rivlin-type material model (1.17). 72 vi,min=1vi,max=2.05vi=1.25DP/Μ=0.74DP/Μ=0.67DP/Μ=0.71vi,max=2.05vi,min=10.81.01.21.41.61.82.02468viri12345670.00.20.40.60.81.01.2siDPΜ CHAPTER IV Channel Confinement Swelling of Hyperelastic Plugs and Tubes 4.1 Introduction The swelling of highly deformable gels and other soft rubbery materials within a con- fined space will exert forces on the walls of their container. Conversely, the confining walls will act on the soft solid, causing it to distort as it continues to expand. Under- standing these effects are important for the design of soft material actuation devices. Similarly, certain biological processes, either of long term growth, or of short term swelling, give rise to tissue distortion as organs and vessels impinge. In this chapter we model this type of process in a simple setting where the contain- ment is due to an open ended rigid cylinder with circular cross-section – an open tube – within which an initially unswollen soft solid expands. As one might anticipate, this leads to a formulation in cylindrical coordinates (unlike the spherical coordinate considered in chapters II and III). Free swelling, which is a simple homogeneous deformation, takes place until the solid makes contact with the confining tube wall. After contact, we seek to determine the shape of the expanding solid as it continues to swell while being confined. At the same time we also seek to determine the traction forces that this expanding solid exerts on the confining tube. 73 The expansion of a solid cylindrical plug of original radius Ro within a rigid cylindrical tube of radius Rc > Ro is studied in Section 4.2 where pressure-swelling response graphs are obtained for a conventional hyperelastic model. In Section 4.3 we consider the effect of a plug that has an internal channel, i.e., an annular plug. Now the deformation after wall contact is no longer one of homogeneous deforma- tion. The associated boundary value problem is formulated, and in Section 4.4 this boundary value problem is solved for the case of a neo-Hookean type swelling model. Wall contact now gives a deformation in which swelling combines axial lengthen- ing with internal channel narrowing. Of particular interest is the closing behavior of the internal channel as the swelling proceeds. Treating the associated boundary value problem provides asymptotic expressions for the channel radius closing and the contact pressure in the large swelling regime. 74 4.2 The homogeneous deformation of laterally confined swelling We first examine the case in which the swelling material is a solid cylindrical plug prior to swelling. Subsequent expansion takes place freely inside a rigid pipe until it makes contact with the pipe wall (Figure 4.1). After so plugging the pipe all further expansion is necessarily along the pipe axis, whereupon the associated wall contact pressure is of particular interest. With respect to a common origin, take circular cylindrical coordinates (R, Θ, Z) in the reference configuration and (r, θ, z) in the current configuration. Let {eR, eΘ, eZ} and {er, eθ, ez} denote unit basis vectors in the reference and deformed configurations, respectively. Consider a circular disk of the material with original (i.e. unswollen) radius Ro and original length L when v = 1. Swelling is then described by an increase in v. In this chapter we restrict attention to spatially uniform swelling, so that v is independent of X. Free swelling F = v1/3I can then be expressed in the above radial coordinates as r = v1/3R, θ = Θ, z = v1/3Z. Moreover, I1 = 3v2/3, I2 = 3v4/3, B = v2/3I so that (1.8) gives where the notation |v1/3I denotes evaluation at the above values for F = v1/3I. By taking p = 2v−1/3 ∂W , the stress tensor T = 0 and so all surfaces ∂I1 +4v1/3 ∂W ∂I2 are traction free. We now suppose that an original unswollen disk is placed in a rigid pipe with an inner channel radius Rc > Ro. Then free swelling may proceed so long as v < R3 c/R3 o. However for v > R3 c/R3 o the pipe provides a lateral confinement and all subsequent volume change must be accommodated by lengthening of the disk in the Z-direction because the disk plugs the pipe as shown in Figures 4.1 and 4.2. 75 T = 2v−1/3 ∂W ∂I1 + 4v1/3 ∂W ∂I2 (cid:19) I − p (cid:12)(cid:12)(cid:12)(cid:12)v1/3I (cid:18) (cid:12)(cid:12)(cid:12)(cid:12)v1/3I (cid:12)(cid:12)(cid:12)(cid:12)v1/3I (cid:12)(cid:12)(cid:12)(cid:12)v1/3I Figure 4.1 Solid cylinder with the radius Ro and the confinement pipe with radius Rc > Ro. Figure 4.2 Representation of the plug swelling within a rigid pipe. For v < λ3 lat the plug has yet to make contact with the pipe wall. Contact first occurs when v = R3 lat all further swelling is directed into Z-direction extension. This generates the pressure Plat between the plug and the pipe wall. lat. For v > λ3 c/R3 o = λ3 If the contact between the plug and the pipe is frictionless then for v > R3 c/R3 o it follows that T = −Plat(er ⊗ er + eθ ⊗ eθ), F = λlat(er ⊗ eR + eθ ⊗ eΘ) + ez ⊗ eZ, v λ2 lat (4.1) where we have defined the fixed constant λlat ≡ Rc/Ro. (4.2) This λlat is the free swelling stretch that just causes all-around contact of the swollen specimen with the pipe’s inner wall. This fixed constant (4.2) will appear often in what follows. 76 Consequently, the condition v > R3 c/R3 o will be written as v > λ3 lat when the contact condition is met. The Plat appearing in (4.1)1 is the lateral confining pressure. In contrast to the known and fixed λlat, the lateral confining pressure Plat is as yet unknown. Consistency with free swelling ensures that Plat = 0 when v = λ3 lat (free- swelling that has just made contact with the wall). Our interest is in determining the dependence of Plat upon v for v > λ3 lat (see the Fig. 4.2 schematic). It is useful to note that the problem can also be described with respect to fixed Cartesian coordinates with an orthonormal basis {e1, e2, e3}. Letting the e3 direc- tion coincide with the Z-direction of the cylindrical geometry it follows that (4.1) is equivalent to T = −Plat(e1 ⊗ e1 + e2 ⊗ e2), F = λlat(e1 ⊗ e1 + e2 ⊗ e2) + e3 ⊗ e3. v λ2 lat This form for F makes B = λ2 lat + v2/λ4 lat ≡ ˜I1(v) and I2 = 2v2/λ2 2λ2 T11(e1 ⊗ e1 + e2 ⊗ e2) + T33e3 ⊗ e3 with lat(e1 ⊗ e1 + e2 ⊗ e2) + v2/λ4 late3 ⊗ e3 so that I1 = lat ≡ ˜I2(v). Thus (1.8) gives T = lat + λ4 (cid:12)(cid:12)(cid:12)(cid:12) ˜I1, ˜I2 (cid:12)(cid:12)(cid:12)(cid:12) ˜I1, ˜I2 ∂W ∂I1 ∂W ∂I1 T11 = 2 T33 = 2 λ2 lat v v λ4 lat + 2 + 4 ∂W ∂I2 ∂W ∂I2 (cid:12)(cid:12)(cid:12)(cid:12) ˜I1, ˜I2 (cid:12)(cid:12)(cid:12)(cid:12) ˜I1, ˜I2 (cid:18) λ4 lat v + v λ2 lat (cid:19) − p, − p. v λ2 lat In polar coordinates Trr = Tθθ = T11 and Tzz = T33. Invoking −Plat = T11 and T33 = 0 to eliminate p gives Plat = 2 (cid:32) (cid:12)(cid:12)(cid:12)(cid:12) ˜I1, ˜I2 ∂W ∂I1 (cid:12)(cid:12)(cid:12)(cid:12) ˜I1, ˜I2 (cid:33)(cid:18) v λ4 lat (cid:19) . − λ2 lat v + λ2 lat ∂W ∂I2 This is the relation of the lateral confining pressure exerted on the solid cylinder due 77 Figure 4.3 Confinement Pressure Plat as a function of swelling v in a solid cylinder taking the outer radius Ro = 1 (so that λlat = Rc). The graphs are for the pipe radii Rc = 1.51/3 on the left and Rc = 21/3 on the right. Here Plat is normalized by the material modulus µ. to the all-around contact after the cylinder has plugged the pipe (v > λ3 o). lat = R3 c/R3 Turning to the specific material model introduced in Section 1.2 this now gives the result of using (1.15) directly, namely, Plat = µo (cid:18) v λ4 lat (cid:19) . − λ2 lat v (4.3) We note that for v = λ3 lat the lateral pressure (4.3) gives Plat = 0 which is consistent with the fact that the cylinder has no lateral pressure when it has just made contact with the surrounding pipe. Graphs of the confining pressure Plat (normalized by µo) as a function of the swelling v for this material model is shown in Figure 4.3. For this figure we take Ro = 1 and the two different values Rc = 1.51/3 and Rc = 21/3 so as to initiate contact either at v = 1.5 or v = 2 (equivalently, we consider λlat = 1.51/3 and 21/3). 78 1.01.52.02.50.00.20.40.60.81.01.2SwellingamountvNormalizedConfinementPressure 4.3 The confinement boundary value problem for an annular plug The problem for laterally confined swelling as just considered in Section 4.2 is a prob- lem for homogeneous deformation. Thus all of the various scalar and tensor quantities: F, T, p, while varying with v in a parametric fashion, were always independent of spatial position x. As such, the equilibrium equation div T = 0 was automatically satisfied. The situation changes if we replace the originally unswollen solid disk of radius Ro with an annular disk of outer radius Ro and inner radius Ri > 0. This originally unswollen annular disk is again placed in a channel of radius Rc > Ro and allowed to swell. Under these circumstances there will once again be a period of free swelling, specifically as v increases from its initial value v = 1 to the value v = R3 c/R3 o = λ3 lat during which the deformation gradient is F = v1/3I and this free swelling with the homogeneous deformation generates no stresses and hence T = 0. Note that this free swelling period of deformation is again expressed as r = v1/3R, θ = Θ and z = v1/3Z. The outer radius attains the value r(Ro) = Rc when v = λ3 lat and it is again at this point that the solid cylinder makes contact with the pipe wall. Further swelling can now be accommodated not only by axial lengthening but also by continued change in the inner radius of the disk. The resulting inhomogeneous deformation will no longer render div T = 0 in a trivial fashion. It is to this issue that we now turn our attention. We consider the following somewhat more general symmetric deformation r = r(R), θ = Θ, z = λzZ, (4.4) for Ri ≤ R ≤ Ro, 0 ≤ Θ < 2π, − 1 2L ≤ Z ≤ 1 2L. With respect to this 79 deformation define ri = r(Ri) and ro = r(Ro). The deformation gradient is given by F = r(cid:48)(er ⊗ eR) + r R (eθ ⊗ eΘ) + λz(ez ⊗ eZ) (4.5) in which prime ((cid:48)) denotes the derivative with respect to R, i.e., r(cid:48) = dr/dR. Both C and B are diagonal with respect to their polar basis sets. The principal stretch are clearly λr = r(cid:48), λθ = r/R and λz. Because the swelling is prescribed in the amount v it follows from (1.1) that v = rr (cid:48) λz R , and because the swelling v is spatially constant, this can be integrated to give r2 = r2 i + (R2 − R2 i ). v λz (4.6) (4.7) Once contact is made with the rigid pipe wall the outer radius of the plug stays fixed at the radius Rc. Thus once v ≥ λ3 lat it follows that ro = r(Ro) = λlatRo = Rc. (4.8) The inner radius is then determined from (4.7) and (4.8) in terms of λz as (cid:18) (cid:19)1/2 ri = c − v R2 λz R2 o + v λz R2 i . (4.9) The condition ri ≥ 0 gives the lower bound λz ≥ R2 o − R2 R2 c i v. (4.10) The symmetry of the deformation (4.4) makes the stress tensor take the form 80 T = Trrer ⊗ er + Tθθeθ ⊗ eθ + Tzzez ⊗ ez. The contact pressure at the rigid channel wall R = Rc will still be denoted by Plat, giving Trr|r=Rc = −Plat where again Plat needs to be determined. Unlike the spatially constant swelling characterized by (4.1) when no channel is present, it will now be the case that all three quantities Trr, Tθθ, Tzz are functions of r, or, equivalently, of R. It will also generally be the case that Trr is no longer equal to Tθθ. The stress equations of equilibrium then give that the hydrostatic pressure p is a function of r (or R), and that Trr and Tθθ relate to each other via ∂Trr ∂r + 1 r (Trr − Tθθ) = 0. (4.11) The boundary conditions for this problem are (4.8) at the contact surface in conjunc- tion with the traction free conditions on the remaining free surfaces. On the inner radius this gives Trr = 0. However there is the possibility that the inner channel completely closes in which case this surface collapses to a line segment and ceases to provide a boundary condition. Thus the original reference inner radius at R = Ri is subject to the condition Trr|ri = 0 if ri > 0. (4.12) The remaining boundaries are the annular caps at z = − 1 2λzL and z = 1 2λzL. The traction free condition associated with such a surface then formally requires the van- ishing of Tzz for all R obeying Ri < R < Ro (equivalently all r obeying ri < r < Rc). However because Tzz is no longer constant, such a pointwise condition cannot be met. This reflects the fact that z = λzZ in (4.4) is too simplistic of an assumption on the detailed nature of the inhomogeneous deformation. We invoke the usual remedy of St. Venant’s principle in which (4.4) approximates the deformation away from the caps (in a long cylinder approximation). Consequently, we abandon the point-wise condition Tzz = 0 and replace it with a condition that the resultant axial force must 81 vanish i.e., Rc(cid:90) ri TzzdA = 0 =⇒ 2π Rc(cid:90) ri Tzzrdr = 0. (4.13) We proceed to investigate this boundary value problem for the case in which W = W (I4, v) as is the case for the material model (1.15). We then have from (1.8), (4.5) and (4.6) that the polar coordinate stress components are 1 (cid:12)(cid:12)(cid:12)I lat (cid:12)(cid:12)(cid:12)I lat (cid:12)(cid:12)(cid:12)I lat 1 1 Trr = Tθθ = Tzz = 2 v ∂W ∂I1 2 v 2 v ∂W ∂I1 ∂W ∂I1 r2 v2R2 λ2 zr2 − p, R2 − p, z − p. λ2 (4.14) (4.15) (4.16) The notation is to indicate that I1 is evaluated at the value I lat 1 = I lat 1 (R, ri, λz, v) with 1 (R, ri, λz, v) = (r(cid:48))2 + I lat r2 R2 + λ2 z = v2R4 + (vR2 − vR2 λzR2(vR2 − vR2 i )2 i + λzr2 i + λzr2 i ) + λ2 z where use has been made of (4.6) and (4.7) in arriving at the final expression. We now integrate the equation of equilibrium (4.11) in the form r(cid:90) (cid:124) ∂Trr ∂r (cid:123)(cid:122) dr (cid:125) ri Trr(r)−Trr(ri) r(cid:90) ri = − (Trr − Tθθ) dr. 1 r Making use of (4.12), (4.14), (4.15) then gives Trr(r) = − 2 v ∂W ∂I1 (r(cid:48))2 − r2 R2 (cid:19) dr r . (4.17) r(cid:90) ri (cid:18) (cid:12)(cid:12)(cid:12)I lat 1 82 The integration can be performed with the aid of (4.6), giving Trr as a function of the parameters λz, ri and v. In particular, evaluating this expression at r = ro = Rc gives Plat, which at this stage is a function of λz, ri and v. To render everything as a function of v it is necessary to determine ri and λz as functions of v. One equation for this purpose is given by (4.9), however another equation is also needed. This additional equation is provided by (4.13). For this purpose, Tzz follows from (4.16) where p = p(r) follows from (4.14) using the explicit form for Trr that has already been obtained from (4.17). In the next section we demonstrate these procedures for the case of material model (1.15). 83 Examining each integral Ja and Jb separately using (4.6) and (4.7) we have that 4.4 Annulus contact for the neo-Hookean type swelling model The derivative ∂W /∂I1 = 1 2µo for the material model (1.15) so that (4.17) becomes Trr(R) = −µo v (r(cid:48))2 dr r + µo v dr r µo v (−Ja + Jb) . (4.18) r(cid:90) (cid:124) ri (cid:123)(cid:122) Ja (cid:125) r(cid:90) ri r(cid:90) (cid:124) ri r2 R2 (cid:123)(cid:122) Jb = (cid:125) (cid:18) λz v r(cid:90) v2 λ2 z R2 r3 dr = v2 λ2 z (r2 − r2 i ) + R2 i Ja = (cid:19)2 dr r(cid:90) (cid:18) vR r(cid:90) rλz ri v λz v λz v λz = = = + (cid:21) (cid:21) dr r (cid:20) r (cid:20) r ri ri ri ln ln = (cid:18) r v2 λ2 z + + v2 2λ2 z v2 2λ2 z (cid:18) (cid:18) ri dr r3 (cid:19) r(cid:90) (cid:19)(cid:18) 1 (cid:19) ri r2 i − r2 i λz v + R2 i − r2 i λz v − R2 r2 + + R2 i R2 i r2 i (cid:19) − 1 r2 (cid:19) dr r3 (4.19) (4.20) and r(cid:90) ri 1 R2 Jb = (cid:122)(cid:125)(cid:124)(cid:123) (vR/λz) dR r dr = R(cid:90) Ri v λz (cid:20) R (cid:21) Ri . dR R = v λz ln Combining the previous results (4.18), (4.19), (4.20) now yields (cid:20) ri R (cid:21) r Ri (cid:19) (cid:18) R2 r2 − R2 i r2 i Trr(R) = µo λz ln + µov 2λ2 z . (4.21) The hydrostatic pressure p now follows from (4.14) and (4.21), whereupon the other 84 stress components from (4.15) and (4.16) become (cid:19) (cid:18) R2 (cid:19) (cid:18) R2 − µov λ2 z r2 (cid:18) r2 (cid:19) R2 z − µov λ2 λ2 z + µo λz ln r2 + (cid:20) ri R (cid:21) (cid:20) ri R r Ri ln µo λz + r Ri (cid:21) + µov 2λ2 z (cid:18) R2 (cid:19) (cid:18) R2 r2 − R2 r2 − R2 µov 2λ2 z i r2 i i r2 i . (cid:19) , (4.22) Tθθ(R) = Tzz(R) = µo v µo v It remains to determine ri and λz in terms of v. The axial stress expression (4.22) enables one to express the zero axial load con- dition (4.13) in the form (cid:18) 1 z − v λ2 2λ2 z R2 i r2 i v (cid:19) Rc(cid:90) (cid:124) (cid:123)(cid:122) (cid:125) rdr ri Jc Here we have defined Rc(cid:90) (cid:124) Rc(cid:90) ri ri − v 2λ2 z Jc = (cid:19) (cid:18) R2 (cid:123)(cid:122) r2 Jd ri Rc(cid:90) (cid:124) (cid:1) , rdr − 1 λz (cid:125) (cid:0)R2 rdr = 1 2 c − r2 i (cid:18) ln (cid:20) rRi (cid:123)(cid:122) riR Je (cid:21)(cid:19) rdr (cid:125) = 0. (4.23) and Rc(cid:90) ri λz v λz 2v Jd = = = Rc(cid:90) (cid:18)λz v ri = R2 dr r (cid:18) (r2 − r2 i ) + R2 i (cid:19) (cid:20)Rc (cid:21) (cid:19) dr r Jc + (cid:0)R2 r2 i (cid:18) i − λz R2 v (cid:1) + 1 2 ln ri i − λz R2 v c − r2 i (cid:19) ln (cid:21) (cid:20)R2 c r2 i , r2 i 85 (cid:18) ln (cid:20) R Ri Rc(cid:90) ri (cid:21)(cid:19) (vR/λz) dR (cid:122)(cid:125)(cid:124)(cid:123) r dr R dR (cid:26) (cid:20)R2 c r2 i (cid:0)R2 (cid:21) −1 4 − v λz o − R2 i R2 o ln (cid:20)Ro (cid:21)(cid:27) Ri R2 o ln 1 2 (cid:21)(cid:17) , (cid:1) + (cid:20)R2 o R2 i r dr = r dr − ri ln riR (cid:18) (cid:18) (cid:20) rRi (cid:21)(cid:19) Rc(cid:90) (cid:20) r (cid:21)(cid:19) Rc(cid:90) (cid:26) (cid:0)R2 (cid:16) − R2 c − r2 (cid:124) −1 4 c + r2 1 4 ln ri i ri and Je = = = = so that ri (cid:20) r (cid:21)(cid:19) (cid:20) R (cid:21)(cid:19) (cid:21)(cid:27) Ri − v λz +R2 c ln ri Ri ln ln (cid:18) Rc(cid:90) (cid:18) Ro(cid:90) (cid:20) Rc ri (cid:125) o − R2 i ) (cid:20) R2 c ln R2 c ln c r2 i r dr − v λz (cid:1) + 1 2 R2 (R2 (cid:123)(cid:122) R2 c i + v λz Je = 1 4 (cid:21) (cid:21) (cid:20)R2 o R2 i . − v 4λz R2 o ln The axial force balance (4.23) thus contains terms either with or without ln[·] func- tions. Using the expressions for Jc and Jd in the right side of (4.23) one finds that the terms that do not contain logarithm factors sum to (cid:18) λ2 z v 1 2 ( R2 (cid:124) (cid:123)(cid:122) (cid:125) c − r2 o−R2 i i )/λz v(R2 ) (cid:19) − v 2λ2 z R2 i r2 i − 1 2λz = v 4λ2 z o − R2 i ) (R2 (cid:18)2λ3 z v − v λz R2 i r2 i (cid:19) − 1 . Using the expressions for Jd and Je in the right side of (4.23) one finds that the terms that contain the logarithm factors sum to (cid:18) 1 4λz − v λz (cid:124) R2 i + r2 (cid:123)(cid:122) (cid:125) i − R2 c −vR2 o/λz (cid:19) ln (cid:21) (cid:20) R2 c r2 i (cid:21) (cid:20)R2 o R2 i + v 4λ2 z R2 o ln (cid:20) r2 i R2 o cR2 R2 i (cid:21) . = v 4λ2 z R2 o ln Combining these results and dividing by vR2 o/4λ2 z puts the zero axial load condition 86 into the form (cid:21) (cid:20) r2 i R2 o cR2 R2 i (cid:18) + ln 1 − R2 i R2 o (cid:19)(cid:18)2λ3 z v − vR2 i λzr2 i (cid:19) − 1 = 0. We now eliminate ri in this equation by substituting from (4.9). The various radial distances can be normalized out by introducing the initial thickness ratio ζ ≡ Ri/Ro < 1 and recalling that λlat = Rc/Ro. The result of this long string of calculations that began with (4.23) is: (cid:18) (cid:19) (cid:20) 1 ζ 2 − v ln (cid:0)1 − ζ 2 (cid:1)(cid:21) λzλ2 lat ζ 2 − (1− ζ 2) 1 − 2λ3 v z + vζ 2 lat + v(ζ 2 − 1) λzλ2 = 0. (4.24) Equation (4.24) is a single equation for the determination of λz = λz(v) at any swelling value v ≥ λ3 lat associated with wall contact. In terms of these normalizations the restriction (4.10) on λz now writes itself as λz ≥ (1 − ζ 2)v/λ2 lat. (4.25) Once λz is known then ri follows from (4.9). Finally, with the two primary unknowns ri and λz so determined, the lateral confinement pressure Plat follows from the evalu- ation of (4.21) at R = Ro. Specifically, Plat = Plat(v, ζ, λlat) is found to take the form (cid:18) (cid:21) + µov 2λ2 z latζ 2 λzλ2 lat − v(1 − ζ 2) λzλ2 (cid:20) Plat = µo 2λz ln using λz = λz(v). (cid:19) − 1 λ2 lat , (4.26) λzζ 2 lat − v(1 − ζ 2) λzλ2 The initial contact value v = λ3 lat causes (4.24) to be satisfied with λz = λlat for all ζ. These values then make Plat = 0 thus confirming consistency with all of the conditions associated with first wall contact. 87 4.4.1 Existence and uniqueness The natural question that arises is whether (4.24) has a unique solution for λz obeying (4.25) for all v > λ3 lat. For this purpose we define G(λz; v, ζ, λlat) ≡ ln (cid:1)(cid:21) , (cid:0)1 − ζ 2 ζ 2 (cid:20) 1 ζ 2 − v (cid:18) 1 − 2λ3 v z + λzλ2 lat (cid:19) . (4.27) (4.28) H(λz; v, ζ, λlat) ≡ (1 − ζ 2) vζ 2 lat − v(1 − ζ 2) λzλ2 Now the solution to (4.24) for λz is equivalent to the solution of G = H. Note that at first wall contact, where v = λ3 lat and λz = λlat, one finds that both G and H vanish, i.e., G(λlat; λ3 lat, ζ, λlat) = 0, H(λlat; λ3 lat, ζ, λlat) = 0, confirming that G = H at the first wall contact. The following two lemmas will be used to establish the existence of a solution. Lemma 1. If λz ↓ (1 − ζ 2) v λ2 lat then G < H. Proof. The notation ↓ means that the approach is from above. We do not simply evaluate λz at (1 − ζ 2)v/λ2 lat because λz ↓ (1 − ζ 2)v/λ2 lat ⇒ G → −∞, and H → ∞, which will be sufficient to establish the result. To see that indeed G → −∞ take λz = (1 − ζ 2 + δ1)v/λ2 lat whereupon G = ln[δ1/(ζ 2(1 − ζ 2 + δ1))] → −∞ as δ1 ↓ 0. To see that indeed H → ∞ take λz = δ2 + (1 − ζ 2)v/λ2 lat whereupon (cid:18) v2ζ 2λ4 H = (1−ζ 2) lat + δ2vλ6 lat − 2δ2(v − vζ 2 + δ2λ2 lat)3 → (1−ζ 2) (cid:19) (cid:18) v2ζ 2λ4 lat δ2vλ6 lat (cid:19) → ∞ δ2vλ6 lat 88 as δ2 ↓ 0. (cid:3) Lemma 2. If λz → ∞ then G > H. Proof. It follows that lim λz→∞ G(λz; v, ζ, λlat) = ln (cid:20) 1 (cid:21) ζ 2 > 0, since 0 < ζ < 1. Also λz→∞ H(λz; v, ζ, λlat) → (1 − ζ 2) lim thus confirming the result. (cid:3) (cid:19) (cid:18) −2λ3 v z → −∞, Thus by continuity the graphs of G and H must intersect on λz > (1 − ζ 2)v/λ2 lat thereby ensuring existence. Furthermore, if it is also the case that G is increasing and H is decreasing everywhere on this range, then the intersection will be unique. This is ensured by the next lemma. Lemma 3. For fixed v > 0, 0 < ζ < 1, λlat > 0 it follows that dG dλz > 0, dH dλz < 0 for all λz > (1 − ζ 2) v λ2 lat . Proof. The derivative of (4.27) is dG dλz = v (1 − ζ 2) latλz − v (1 − ζ 2)) λz (λ2 . (4.29) On the right hand side of this equation the numerator v(1−ζ 2) is positive. The factor in the denominator (λ2 latλz − v(1 − ζ 2)) is positive because λz > (1 − ζ 2)v/λ2 lat. Thus this derivative is positive and the first condition is confirmed. 89 Figure 4.4 Axial elongation λz as a function of swelling v for the expanding annular plug with inner radius Ri = 1/2 and outer radius Ro = 1 (so that λlat = Rc and ζ = 1/2). Wall contact occurs when v = λ3 lat. The graphs are for two separate cases of outer pipe radius: Rc = 1.145 and Rc = 1.260, which are chosen so as to give contact v values of 1.5 and 2, respectively. The slope of the curves immediately after contact are given by (4.31). The derivative of (4.28) is dH dλz = −(cid:2)1 − ζ 2(cid:3)(cid:32)(cid:20)6λz v (cid:34) (cid:21) 2 + vλ2 latζ 2 latλz − v (1 − ζ 2))2 (λ2 (cid:35)(cid:33) , (4.30) where all of the terms in brackets [·] are positive. Hence this derivative is negative and the lemma is proved. (cid:3) 4.4.2 Numerical demonstrations As an example, consider the axial stretch λz as determined by (4.24) for a case where Ri = 1/2 and Ro = 1 (making ζ = 0.5). Figure 4.4 plots λz as a function of v for the two different pipe confinement radii considered previously in Section 4.2: 90 Figure 4.5 Deformed inner radius ri as a function of swelling v for the wall contact cases from Fig. 4.4. Prior to wall contact the swelling v < λ3 lat and the homogeneous deformation causes ri to increase. After contact, when v > λ3 lat, the inner radius is monotonically decreasing with swelling. Rc = 1.51/3 = 1.145 and Rc = 21/3 = 1.260. Consequently, λlat = Rc/Ro = Rc. Wall contact takes place when v = λ3 lat = R3 c and hence either v = 1.5 or v = 2. These aspects all mirror the situation considered previously with respect to Figure 4.3, although of course there was no channel in that case (formally Ri = 0) whereas now we are taking Ri = 0.5. Prior to wall contact (v < λ3 lat) the homogeneous free-swelling deformation gives λz = v1/3 and this relation is the first part of the graphs in Figure 4.4. At wall con- tact the axial elongation noticibly increases. By implicitly differentiating (4.24) with respect to v and evaluating the result at v = λ3 lat one obtains the slope immediately after first contact dλz dv (cid:12)(cid:12)(cid:12)f c = (cid:18) 1 + ζ 2 (cid:19) 1 1 + 3ζ 2 λ2 lat > 1 3λ2 lat , (4.31) where the value 1 3λ−2 contact with the wall. lat is the slope of the free-swelling part of the curve just before 91 Rc=1.513 = -1+Ζ2Ri1+3Ζ2Λlat2>-0.164slopewall contactsRc=213123450.00.10.20.30.40.50.60.7SwellingvInnerRadiusri Figure 4.6 Confinement pressure Plat as a function of swelling v for the wall contact cases from Figs. 4.4 and 4.5. Here Plat is normalized by the material modulus µo. The slope at first contact is given by (4.32). 92 Prior to wall contact the inner radius ri is increasing with v according to the simple homogeneous deformation rule ri = v1/3Ri. At wall contact this increase is reversed so that after wall contact ri is decreasing with v. These decreasing values are found by using the previously obtained λz = λz(v) in the kinematic condition (4.9). Hence the channel achieves a maximum of ri given by λlatRi = RcRi/Ro at first wall contact when v = (Rc/Ro)3. This is shown in Figure 4.5 for the same cases considered in Fig. 4.4. Differentiating (4.9) and using (4.31) gives the decreasing slope immediately after first contact (cid:12)(cid:12)(cid:12)f c dri dv (cid:18) 1 − ζ 2 (cid:19) Ri 1 + 3ζ 2 λ2 lat < 0, = − thus showing the abrupt turn-around in the radial deformation. As was the case for the plug without a channel (Section 4.2) the surface R = Ro is traction free until wall contact. After becoming confined by the wall at Rc the contact pressure is given by (4.26). The monotonic increase of this Plat with v is shown in Figure 4.6 for the same cases as considered in Figs. 4.4 and 4.5. Differentiating (4.26) and using (4.31) it is found that the slope of these curves just after contact is given by (cid:12)(cid:12)(cid:12)f c = (cid:18) 1 − ζ 2 (cid:19) 2µ 1 + 3ζ 2 λ4 lat dPlat dv > 0. (4.32) This can be compared to the slope of Plat at first contact for the plug without a channel that is obtained by differentiation of (4.3) and evaluation at v = λ3 lat yielding lat. As one might have anticipated, this matches (4.32) in the ζ → 0 limit corresponding to zero channel radius. More generally, the initial slope (4.32) of the 2µo/λ4 contact pressure is decreasing with ζ. This is consistent with the expectation that the contact force generated by the expanding annulus would decrease as its wall becomes thinner. One would similarly anticipate that the pressure curves would exhibit this 93 Figure 4.7 Confinement pressure Plat for a solid plug (Ri = 0) and for tubes (annulus plugs with Ri = 1/10, 1/4, 1/2, 3/4) as a function of swelling v. In all cases Ro = 1 (so that λlat = Rc). The confining radius is Rc = 1.51/3 so as to be consistent with one of the cases shown previously in Figs. 4.4, 4.5 and 4.6. Plat is again normalized by the material modulus µo. The zero channel radius limit (ζ = 0) recovers the neo-Hookean type curves in Fig. 4.3. Thinner walls (larger ζ) give less contact pressure for the same value of swelling. 94 decreasing contact pressure with wall thickness not only at first contact, but for all values of swelling. This is confirmed in Figure 4.7 where the contact pressure as a function of v is shown for a variety of thicknesses ratios ζ. 4.4.3 Asymptotic behavior for large swelling values The examples of Section 4.4.2 raise an obvious question as to the channel closing behavior in the large v limit. Specifically, when viewing a graph like that given in Fig. 4.5 one can ask which of the three logical alternatives occurs: (a) the channel closes at some finite v, (b) the channel never fully closes but approaches a zero radius as v tends to infinity, or (c) the channel radius approaches a finite positive asymptote as v tends to infinity? To answer this question we first perform a v → ∞ asymptotic analysis of the condition G = H in order to extract the large v behavior of λz. To this end we first observe from (4.25) that the bound v/λz ≤ λ2 lat/(1 − ζ 2) gives that λz → ∞ as v → ∞. If λz is order Cvm in this limit then m ≥ 1. Furthermore if m > 1 then C > 0 and if m = 1 then C ≥ (1 − ζ 2)/λ2 lat. Consider a case in which λz ∼ Cvm with m > 1. It would then follow from (4.27) and (4.28) that (cid:20) 1 (cid:21) ζ 2 G = ln + o(1) and H = −(1 − ζ 2) 2λ3 z v + O(1) as v → ∞. (4.33) This case is now discarded because (4.33) is inconsistent with meeting the condition G = H in the large v limit. It is thus concluded that λz = Cv + o(v) with C ≥ (1 − ζ 2)/λ2 lat as v → ∞. 95 If the inequality is strict, i.e., if C > (1 − ζ 2)/λ2 lat, this then makes (cid:20) 1 − 1 − ζ 2 Cλ2 lat (cid:21) G = ln + o(1) + O(1), and H = (1 − ζ 2) 1 − 2C 3v2 + o(v2) + (cid:18) (cid:19) + O(1) , ζ 2 Cλ2 lat − (1 − ζ 2) so that once again it is not possible to meet the condition G = H in the large v limit. Therefore the remaining possibility is that λz ∼ Cv with C = (1−ζ 2)/λ2 lat whereupon we write λz = Cv (1 + 1(v)) with C = (1 − ζ 2)/λ2 lat and 1(v) = o(1) as v → ∞. (4.34) It then follows that (cid:20) 1 (cid:21) ζ 2 (cid:20) 1 (cid:21) 1 + 1 + ln G = ln = ln[1] + O(1) and H = (1 − ζ 2) = (1 − ζ 2) (cid:18) (cid:18) 1 − 2 v C 3v3(1 + 31 + O(2 1)) + 1 − 2C 3v2 + O(1v2) + (cid:19) vζ 2 v(1 − ζ 2)(1 + 1) − v(1 − ζ 2) 1 1 ζ 2 1 − ζ 2 (cid:19) = −2(1 − ζ 2)C 3v2 + ζ 2 1 + O(1) + o(v2). Hence the condition G = H gives ln[1] + O(1) = −2(1 − ζ 2)C 3v2 + ζ 2 1 + O(1) + o(v2), as v → ∞. (4.35) There are three main terms in (4.35) and at this stage it is not obvious which two 96 are dominant, and hence in balance with each other, as v → ∞. We seek the leading order behavior in the function 1(v) so as to provide the appropriate balance between whichever two terms are in fact dominant. There are three possible ways to balance the expression (4.35): II ln[1] = −2(1 − ζ 2)C 3v2 + ζ 2 1 . I III (4.36) Balance possibility I gives 1 ∼ e−2(1−ζ2)C3v2, so that the two balanced terms are O(v2). However, this makes the third (unbalanced) term ζ 2/1 ∼ ζ 2e2(1−ζ2)C3v2 which dominates the O(v2) terms as v → ∞. Hence balance I is inconsistent and so removed from further consideration. Balance possibility II involves two terms that are of fundamentally different orders and hence cannot balance each other. This balance possibility is therefore also discarded. Balance possibility III yields 1 ∼ ζ 2 2(1 − ζ 2)C 3 1 v2 , v → ∞ as which is consistent with the requirement that 1(v) = o(1) as v → ∞. The remaining unbalanced term in (4.36) is O(ln[v]) and hence satisfies the requirement that it is dominated by the O(v2) terms that are in balance for this possibility. Hence III provides a consistent balance, and it is the only consistent balance from among the three possibilities diagrammed in (4.36). Thus the 1(v) appearing in (4.34) can now be written as 1(v) = ζ 2 2(1 − ζ 2)C 3 1 v2 (1 + 2(v)) with 2(v) = o(1) as v → ∞ 97 where C continues to be by given (1 − ζ 2)/λ2 lat. Continuing the analysis now using (cid:19) ζ 2 2(1 − ζ 2)C 3 1 v2 (1 + 2) (cid:18) λz = Cv 1 + it is found that G becomes simply G = −2 ln[v] + O(1). On the other hand, the determination of the leading behavior for H requires a se- quence of tedious manipulations that results in H = (1 − ζ 2) − 3(1 − ζ 2)ζ 2 latC 3 − 2C 3(1 − ζ 2)2v2 + O(2) + O(v−2) + O(2 λ6 2v2). Enforcing the condition G = H now leads to the conclusion that 2 ∼ λ6 lat (1 − ζ 2)4 1 v2 ln[v]. Consolidating all of these results provides the following asymptotic expression for λz as a function of v, λz = 1 − ζ 2 λ2 lat (cid:124) (cid:123)(cid:122) (cid:125) C v + λ4 latζ 2 2(1 − ζ 2)3 (cid:18) 1 (cid:19) v + λ10 latζ 2 2(1 − ζ 2)7 (cid:18) 1 v3 ln[v] (cid:19) + o (cid:19) (cid:18) 1 v3 ln[v] . (4.37) Employing this expansion of λz in the lateral pressure expression (4.26) it is found that Plat = µC 2v + O(ln[v]/v) ∼ µ(1 − ζ 2)2v/λ4 lat. Employing (4.37) in (4.9) now yields the corresponding expansion for ri by using (cid:19)1/2 (cid:18) 1 − Cv λz (cid:19)1/2 (cid:18) 1 − 1 1 + 1 = Roλlat ri = Roλlat = Roλlat 1/2 1 (1 + 1)−1/2, 98 Figure 4.8 Asymptote of the deformed inner radius ri as a function of swelling v in a hollow tube taking the inner radius Ri = 1/2 and the outer radius Ro = 1 (so that λlat = Rc and ζ = 1/2) with Rc = 1.51/3. The graphs are for the pure numerical solution (red curve in the middle) and the asymptotic expansion (4.38) with only the leading term (green curve on the bottom) and with the additional first correction term (blue curve on the top). 99 Wall contactatv=Λlat351015200.00.10.20.30.40.50.6SwellingvInnerRadiusri with 1 again defined as in (4.34). It follows that √ 2 2 ri Ro = λ4 lat ζ (1 − ζ 2)2 (cid:18) 1 (cid:19) v √ 2 4 + λ10 lat ζ (1 − ζ 2)6 (cid:18) 1 v3 ln[v] (cid:19) + o (cid:18) 1 (cid:19) v3 ln[v] , (4.38) which shows for this material model how the channel closes in the large v asymptotic limit. This is depicted in Figure 4.8 where the asymptotic result is compared with the pure numerical treatment. 100 CHAPTER V Channel Confinement Swelling of Internally Balanced Material Plugs and Tubes 5.1 Internally balanced elastic materials Chapters II to IV considered swelling in the purely hyperelastic context. In this chapter we broaden the material constitutive theory so as to place swelling in the framework of internally balanced materials. The theory of internally balanced elas- tic materials as described in Demirkoparan et al. (2014); Demirkoparan and Pence (2015a); Hadoush et al. (2017, 2014, 2015) makes use of the general multiplication decomposition of the deformation gradient F = ˆFF∗ that was also used here in (1.3). Such a decomposition, which was originally invoked in the setting of finite deforma- tion plasticity due to dislocation slip Kr¨oner (1959); Lee (1969), is now commonly used in other mechanical descriptions of finite deformations, including thermoelastic deformation Vujoˇsevi´c and Lubarda (2002), crystallographic transformation plasticity Havner (1992); Forest et al. (2014) and biological growth Chen and Hoger (2000). In such treatments the first factor in (1.3) typically describes the action of the physical effect under consideration and the second factor ˆF corresponds to some sort of elastic recovery or accommodation. This allows ˆF to be described via energy minimization, whereas F∗ need not in general have a corresponding energetic framework. Indeed it may be the case that the phenomena described by F∗ is inherently dissipative. In con- trast, the theory of internally balance elastic materials as studied in Demirkoparan 101 et al. (2014); Demirkoparan and Pence (2015a); Hadoush et al. (2017, 2014, 2015) specifically attributes an energy minimal variational structure to both factors in the decomposition (1.3). The internally balanced material theory as described in Hadoush et al. (2017, 2014, 2015) involved no material constraint whereas the internally balanced material theory as described in Demirkoparan et al. (2014); Demirkoparan and Pence (2015a) involved the material constraint of incompressibility. It is our purpose here to consider the effect of the constraint (1.1) in the internal balance theory. Recall the discussion in page 10 that indicated the possibility of considering functions F∗ in (1.3) that were not prescribed to be equal to v1/3I. For our purpose we instead specify that the individual factors in (1.3) are to be subject to the individual constraints detˆF = 1, detF∗ = v. (5.1) The special case of (5.1) with v = 1 then recovers the treatment as given in Demirkoparan et al. (2014); Demirkoparan and Pence (2015a). One motivation for the particular constraint choice (5.1) is that it makes the de- composition F = ˆFF∗ a generalization of the simpler decomposition F = v1/3 ˆF with det ˆF = 1. This simpler decomposition is a special case of F = ˆFF∗ as can be seen by taking F∗ = v1/3I. The simpler decomposition F = v1/3 ˆF with det ˆF = 1 has served as a useful description of material swelling in a more conventional hyperelas- tic framework Tsai et al. (2004); Pence and Tsai (2005b); Demirkoparan and Pence (2007a); Fang et al. (2011). In contrast the framework presented here does not require that F∗ = v1/3I; hence it offers the possibility of modeling more complicated swelling phenomena. 102 Mirroring conventional hyperelasticity let W denote the mechanical energy storage density per unit volume in the reference configuration. This W is to depend on both ˆF and F∗ and thus any such W is described in terms of alternative and equivalent constitutive specifications W = Φ(0)(F, F∗, v) = Φ(1)(F, ˆF, v) = Φ(2)( ˆF, F∗, v), (5.2) with F, ˆF and F∗ related via (1.3). A specification for any one of the energy density functions Φ(0), Φ(1) or Φ(2) permits determination of the other two by simply requiring consistency with (1.3). In what follows we shall omit v in the argument list of these energy density functions in order to provide a simpler notation. Continuum mechanics principles such as material frame indifference and material symmetry place specific restrictions on the way in which these energy densities may depend upon F, F∗, and ˆF. The principle of material frame indifference specifically gives that W = Φ(3)(C, C∗), with C = FT F, C∗ = (F∗)T F∗, ˆC = ˆFT ˆF. (5.3) (5.4) Although ˆC does not appear in (5.3) it is introduced in (5.4) because of subsequent developments. By virtue of (5.1) it follows that detC = detC∗ = v2 and det ˆC = 1. It is to be remarked that C and C∗ do not determine ˆC, and thus one may not conclude that W can be written as a function of C and ˆC. The function Φ(3) can always be expressed as a symmetric function of the tensors C and C∗, and it is presumed that this is the case in going forward Demirkoparan et al. (2014). 103 (cid:90) (cid:90) (cid:90) As in conventional hyperelasticity, equilibrium deformations minimize the total stored energy minus the work of any body forces and prescribed surface tractions. The total stored energy is computed by the volume integral of W taken with respect to the reference configuration. In terms of Φ(0) this energy integral is E(χ, F∗) = Φ(0)(∇χ, F∗)dVR. (5.5) B The work functional W(χ) is a volume integral over the region that is subject to body forces in conjunction with a surface integral over the portion of the boundary that is subject to prescribed surface tractions. The minimization is both with respect to deformations χ and with respect to the factor F∗. Incorporating the constraints (5.1) by means of Lagrange multiplier fields p = p(X) and q = q(X) the functional that is to remain stationary with respect to variations of χ, F∗, p and q can be taken as I(χ, F∗) = E(χ, F∗) − W(χ) − p(det∇χ − v)dVR − q(det F∗ − v)dVR. (5.6) B B For simplicity in going forward, body forces will be neglected and each boundary point is presumed to be either a free surface or else subject to a prescribed displacement. This makes the work functional W(χ) = 0. Consideration of body forces and consid- eration of more general boundary conditions is easily accomplished by specifying an appropriate W(χ) Demirkoparan et al. (2014). Minimization of I with respect to χ leads to divT = 0, T = 2 v F ∂Φ(3) ∂C FT − pI, (5.7) where div is the divergence operator with respect to the deformed configuration de- scribed by locations x. Minimization with respect to F∗ provides the new aspect of the treatment. This 104 leads to the internal balance requirement in the form of the tensor equation ∂Φ(0) ∂F∗ − q(detF∗)F∗−T = 0. By using (5.3), the internal balance (5.8) is manipulated into the form 2 ∂Φ(3) ∂C∗ − qvC∗−1 = 0. (5.8) (5.9) The partial derivatives with respect to C and C∗ in (5.7)2 and (5.9) involve holding the other tensor fixed (the standard meaning). For the case of no swelling, meaning v = 1, the above equations all reduce to those as given previously in Demirkoparan et al. (2014); Demirkoparan and Pence (2015a). 5.1.1 Retrieving the hyperelastic theory The purely hyperelastic theory of swelling with W = W (C, v) and T given by (1.2) can be recovered form the above framework in which (1.1) is imposed as a constraint. This theory makes no explicit use of the decomposition (1.3), although it can be connected to the present framework by imposing the additional requirement F∗ = v1/3I. In such a case the tensor C∗ in (5.4) becomes a prescribed quantity. As such, the notion of minimizing with respect to either F∗ or C∗ becomes moot so that (5.9) is no longer generated. Alternatively, taking C∗ = v2/3I and temporarily allowing variation with respect to v we have that ∂Φ(3) ∂C∗ = ∂Φ(3) ∂v dv dC∗ (5.10) where, in view of (5.1) and (5.4), dv dC∗ = d dC∗ (det C∗)1/2 = 1 2 (det C∗)1/2C∗−1 = vC∗−1, 1 2 (5.11) 105 so that (5.9) becomes a requirement that v ∂Φ(3) ∂v C∗−1 − qvC∗−1 = 0, which is satisfied identically provided that q takes on the pleasing form q = ∂Φ(3) ∂v . (5.12) (5.13) Thus by a variety of arguments it is confirmed that (5.9) plays no role if F∗ is stip- ulated in the form F∗ = v1/3I. This leaves (5.7) with Φ(3) = Φ(3)(C) which retrieves the theory of hyperelasticity as explored in a variety of papers such as Tsai et al. (2004); Pence and Tsai (2005b); Demirkoparan and Pence (2007b); Fang et al. (2011); Demirkoparan and Pence (2015b); Gou and Pence (2016). 5.1.2 Isotropic materials The framework as summarized in (1.1), (1.3), (5.1), (5.3), (5.7) and (5.9) is a straight- forward generalization of a theory of incompressible internally balanced elastic ma- terials, as described in Demirkoparan and Pence (2015a), to a theory of internally balanced elastic materials subject to prescribed volumetric change. In Demirkoparan et al. (2014) and Demirkoparan and Pence (2015a) the incompressible theory is con- sidered in detail where it is shown that the resulting formulation describes a material that is isotropic if W = Φ(4)(I1, I2, I∗ 1 , I∗ 2 , ˆI1, ˆI2). (5.14) The same logic continues to apply in the present setting (the change from det F∗ = 1 to det F∗ = v in (5.1) has no effect on the formal argument) and thus (5.14) provides constitutive description for an isotropic material in the internally balanced treatment when swelling is present. Two questions immediately arise. The first is whether (5.14) 106 is consistent with (5.3). It is consistent because ˆI1 = C∗−1 : C, and (cid:0)(C∗−1 : C)(C∗−1 : C) − (C∗−1C) : (CC∗−1)(cid:1), ˆI2 = 1 2 (5.15) (5.16) which allows ˆI1 and ˆI2 to be calculated directly from C and C∗. The second question is to what extent (5.14) is consistent with the usual isotropy stipulation in terms of a material symmetry group that consists of all proper rotations. This issue is also discussed in Demirkoparan et al. (2014) where it is shown that (5.1) is consistent with material isotropy in the standard sense. However because of the decomposition (1.3) one may introduce secondary stipulations on the individual parts, including a notion of interstitial symmetry. As discussed in Demirkoparan et al. (2014) the form (5.14) is not generally consistent with the additional notion of isotropic interstitial symmetry. However, the form (5.14) does become consistent with the notion of isotropic intersti- tial symmetry if the argument list of Φ(4) in (5.14) is restricted so as to include only 1 , I∗ (I∗ For the purpose of using (5.7)2 and (5.9) for materials obeying (5.14) it is noted 2 , ˆI1, ˆI2). that ∂Φ(4) ∂I2 (cid:0)I1I − C(cid:1) + (cid:0) ˆI1C∗−1 − C∗−1CC∗−1(cid:1) ∂Φ(4) ∂ ˆI1 ∂Φ(3) ∂C = ∂Φ(4) ∂I1 I + + ∂Φ(4) ∂ ˆI2 C∗−1 (5.17) 107 and ∂Φ(3) ∂C∗ = ∂Φ(4) ∂I∗ 1 I + ∂Φ(4) ∂I∗ 1 I − C∗(cid:1) − ∂Φ(4) (cid:0)I∗ C∗−1CC∗−1 ∂ ˆI1 2 (cid:0) ˆI1C∗−1CC∗−1 − C∗−1CC∗−1CC∗−1(cid:1). (5.18) Alternative expressions may be obtained through use of the connection C∗−1(CC∗−1)n = F∗−1 ˆCnF∗−T . The Cauchy stress tensor (5.7)2 thus becomes −∂Φ(4) ∂ ˆI2 (cid:0) ∂Φ(4) (cid:0) ∂Φ(4) ∂I1 ∂ ˆI1 + 2 v T = 2 v + I1 ∂Φ(4) ∂I2 ∂Φ(4) ∂I2 B2 (cid:1)B − 2 (cid:1) ˆB − 2 v + ˆI1 ∂Φ(4) ∂ ˆI2 ˆB2 − pI ∂Φ(4) ∂ ˆI2 (5.19) v with B = FFT , B∗ = F∗(F∗)T , ˆB = ˆF ˆFT . (5.20) Although B∗ does not appear in (5.19) it is introduced in (5.20) because of subsequent developments. Also (5.1) gives detB = detB∗ = v2 and det ˆB = 1. The internal balance (5.9) is cast into a convenient form by the following sequence of operations: premultiplication by F, postmultiplication by FT , substitution from (5.18), and use of the Cayley Hamilton theorem to express ˆB3 in terms of lower powers. The result is −2I∗ 1 ∂Φ(4) ∂I∗ 2 FC∗FT + 2(cid:0) ∂Φ(4) ∂I∗ 1 (cid:1)B + I∗ 1 ∂Φ(4) ∂I∗ 2 −¯qv ˆB − 2 ∂Φ(4) ∂ ˆI1 ˆB2 + 2 ∂Φ(4) ∂ ˆI2 I = 0, (5.21) 108 where ¯q is a redefined multiplier in place of q (formally ¯qv = qv + 2 ˆI2∂Φ(4)/∂ ˆI2). A natural configuration for an unconstrained material is one for which the Cauchy stress vanishes, i.e., T = 0. For a constrained material, a natural configuration is one for which the Cauchy stress is restricted to take the form of the constraint stress. In conventional isotropic incompressible hyperelasticity, the reference configuration is automatically a natural configuration, meaning that the stress tensor takes the form of a hydrostatic pressure. An analogous result holds for the internally balanced material. To see this for a given v we consider the state of uniform (equi-axial) swelling F = v1/3I. The Cauchy stress (5.19) is then (cid:16) (cid:16)∂Φ(4) ∂ ˆI1 + 2 v T = 2v−1/3 ∂Φ(4) ∂I1 + 4v1/3 ∂Φ(4) ∂I2 − p (cid:17)(cid:12)(cid:12)(cid:12)F=v1/3I + ˆI1 ∂Φ(4) ∂ ˆI2 ˆB − 2 v ∂Φ(4) ∂ ˆI2 (cid:12)(cid:12)(cid:12)F=v1/3I ˆB2. (5.22) (cid:17)(cid:12)(cid:12)(cid:12)F=v1/3I I Here the subscript F = v1/3I indicates I1 = 3v2/3 and I2 = 3v4/3 with the four re- maining arguments ( ˆI1, ˆI2, I∗ 2 ) as yet undetermined. In addition F = v1/3I casts 1 , I∗ the internal balance (5.21) into the form −2v2/3I∗ ∂Φ(4) ∂I∗ (cid:12)(cid:12)(cid:12)F=v1/3I +2(cid:0)v2/3 ∂Φ(4) 2 1 ∂I∗ 1 (cid:12)(cid:12)(cid:12)F=v1/3I C∗ − 2 ∂Φ(4) ∂ ˆI1 + v2/3I∗ 1 ∂Φ(4) ∂I∗ 2 + ∂Φ(4) ∂ ˆI2 ˆB2 − ¯qv ˆB (cid:1)(cid:12)(cid:12)(cid:12)F=v1/3I I = 0. (5.23) This internal balance (5.23) in conjunction with ˆFF∗ = v1/3I, det ˆF = 1, and detF∗ = v is now a set of equations for ˆF, F∗ and ¯q. A solution to this set is given by ˆF = I, F∗ = v1/3I and (cid:16)∂Φ(4) ∂ ˆI1 ¯q = − 2 v − ∂Φ(4) ∂ ˆI2 −v2/3 ∂Φ(4) ∂I∗ 1 +3(v2−v4/3) ∂Φ(4) ∂I∗ 2 109 (cid:17)(cid:12)(cid:12)(cid:12)(3v2/3,3v4/3,3,3,3v2/3,3v4/3) . (5.24) The subscript (3v2/3, 3v4/3, 3, 3, 3v2/3, 3v4/3) in (5.24) denotes the evaluation (I1, I2, ˆI1, ˆI2, I∗ (3v2/3, 3v4/3, 3, 3, 3v2/3, 3v4/3). The particular value of ¯q as given above plays no for- 1 , I∗ 2 ) = mal role other than to show that it allows the satisfaction of the internal balance (5.24). Then entering (5.22) with ˆB = I and the full evaluation (I1, I2, ˆI1, ˆI2, I∗ (3v2/3, 3v4/3, 3, 3, 3v2/3, 3v4/3) one obtains 1 , I∗ 2 ) = T = −¯pI (5.25) where ¯p is given by a sum of terms where, with the exception of the constraint term containing −p, each term involves the evaluation of a derivative of Φ(4) with respect to one of the six arguments (I1, I2, ˆI1, ˆI2, I∗ 2 ). Since the constraint pressure 1 , I∗ p in this expression is still arbitrary it follows that ¯p in (5.25) is similarly arbitrary. Taking ¯p = 0 one obtains from (5.25) that free swelling T = 0 gives an all-around equiaxial volume change F = v1/3I. Furthermore since the volume v is determined by constraint, it follows more generally from (5.25) that equiaxial volume change is consistent with an arbitrary hydrostatic pressure ¯p. All of these results exactly mirror that of conventional isotropic hyperelasticity with a volume constraint. In particular, as in conventional isotropic incompressible hyperelasticity, the reference configuration is automatically a natural configuration. For the internally balanced theory the hyperelastic energy expression (1.15) sug- gests consideration of the form Φ(4)(I1, I2, I∗ 1 , I∗ 2 , ˆI1, ˆI2) = (I1 − 3v2/3) + α 2 α∗ 2 1 − 3v2/3) + (I∗ ( ˆI1 − 3) ˆα 2 (5.26) where the moduli parameters obey α ≥ 0, α∗ ≥ 0, ˆα ≥ 0. Just as the modulus µ in the hyperelastic model (1.15) was allowed to depend on v, the moduli in the internally balanced material (5.26) are also allowed to depend upon v. For the internally balanced material with stored energy density (5.26), the Cauchy 110 stress follows from (5.19) as T = α v B + ˆB − pI. ˆα v The internal balance condition (5.21) reduces to α∗B − ˆα ˆB2 − qv ˆB = 0, (5.27) (5.28) where, recalling that ¯qv = qv + 2 ˆI2∂Φ(4)/∂ ˆI2, we return to the use of q for notational ease. Equivalent forms for (5.28) are and α∗B∗ − ˆα ˆC − qvI = 0, α∗C∗2 − ˆαC − qvC∗ = 0. (5.29) (5.30) If ˆα = 0 and α∗ > 0 it follows from (5.27) that T = (α/v)B − pI so that if v = 1 then the model is equivalent to the usual neo-Hookean material model in hyperelasticity with modulus α. The internal balance (5.28) then formally gives ˆB = v−2/3B and q = α∗v−1/3. If α∗ = 0 and ˆα > 0 it follows from (5.28) that ˆB = I and q = −ˆα/v. This gives T = (α/v)B + (ˆα/v − p)I and ˆα/v can be absorbed into p. Hence if v = 1 the model with α∗ = 0 and ˆα > 0 is also equivalent to a neo-Hookean material with modulus α. Additional understanding of the range of qualitative behaviors follows from con- sideration of the three different two-term dominant balances that may occur between the various terms of the internal balance (5.28). For this purpose we rewrite (5.28) as 2 ˆα ˆB + qv ˆB = α∗ B, (5.31) 111 and inquire into the determination of both ˆB (obeying det ˆB = 1) and q for specified ˆα, α∗, v and B (obeying detB = v2). Consider first a dominant balance between the second and third terms of (5.31), that is 2 ˆα ˆB + qv ˆB = α∗ B. (5.32) Such a balance becomes dominant as α∗ >> ˆα. Specifically consider α∗ → ∞ at finite ˆα. Introduce β = ˆα/α∗, ρ(X) = q(X)v(X)/α∗. (5.33) Here we have written q and v as possible functions of X to emphasize that possibility in the formulation and solution of boundary value problems. We could also have written ˆα and α∗ as a possible function of X so as to treat the case of nonuniform and composite materials. In going forward, we shall cease to indicate this possible dependence upon X. Using the notation (5.33) the internal balance equation (5.31) is recast into the form 2 β ˆB + ρ ˆB = B. (5.34) In the limit β → 0 the solution pair (ρ, ˆB) such that (5.34) is satisfied with det ˆB = 1 is ρ = v2/3 and ˆB = v−2/3B. On the presumption that both ρ and ˆB can be expanded in powers of β it is found from the O(β) analysis of (5.34) (using also det ˆB = 1) that the O(β) corrections are respectively: −(v−2/3/3)I1 and −v−2B2 +(v−2/3)I1B. Hence returning to original variables it is concluded that ˆB = v−2/3B + ˆα α∗ (cid:16)1 3 v−2I1B − v−2B2(cid:17) (cid:16)(cid:16) ˆα α∗ (cid:17)2(cid:17) , + O 112 q = α∗v−1/3 − 1 3 ˆαv−5/3I1 + ˆα O (cid:17) , (cid:16) ˆα α∗ (5.35) as α∗ → ∞ at finite ˆα. Next consider a dominant balance between the first and second terms, that is 2 ˆα ˆB + qv ˆB = α∗ B. (5.36) Such a balance becomes dominant as ˆα >> α∗. Specifically consider ˆα → ∞ at finite α∗. Introduce γ = α∗/ˆα, κ = qv/ˆα, to rewrite equations (5.31) in the form 2 ˆB + κ ˆB = γB. (5.37) (5.38) In the limit γ → 0 the solution pair (κ, ˆB) such that (5.38) is satisfied with det ˆB = 1 is κ = −1 and ˆB = I. On the presumption that both κ and ˆB can be expanded in powers of γ it is found from the O(γ) analysis of (5.38) (using also det ˆB = 1) that the O(γ) corrections are respectively: (1/3)I1 and B − (1/3)I1I. Hence returning to original variables it is concluded that ˆB = I + α∗ ˆα (cid:16) (cid:17) I1I + O (cid:16)(cid:16) α∗ ˆα (cid:17)2(cid:17) , B − 1 3 q = − ˆα v + α∗ v 1 3 I1 + α∗ O (cid:16) α∗ (cid:17) ˆα , (5.39) 113 as ˆα → ∞ at finite α∗. Finally, consider a dominant balance between the first and third terms, that is 2 ˆα ˆB + qv ˆB = α∗ B. (5.40) definite tensors have a unique symmetric positive-definite square root (the symbol This balance gives ˆB ∼ (cid:112)(α∗/ˆα)B where it is recalled that symmetric positive- √· denotes this unique tensor square-root). The scalar q drops out of the balance, and this loss has as a consequence that the condition det ˆB = 1 is met if and only if α∗/ˆα = v−2/3. More formally, this balance becomes dominant as α∗ → ˆαv−2/3 and it is an exact balance for α∗ = ˆαv−2/3. Specifically, it follows that ˆB = v−1/3 √ B, q = 0, (5.41) when α∗ = ˆαv−2/3. The three cases examined above show how the theory is sensitive to the ratio β = ˆα/α∗. With respect to q the above results summarize as  q → −∞, = 0, +∞, if ˆα/α∗ → 0, if ˆα/α∗ = v2/3, if ˆα/α∗ → ∞, (5.42) in a manner that is independent of B. For values of β = ˆα/α∗ that do not correspond to those given in (5.42) it would generally be the case that q is also dependent upon B. Turning to the Cauchy stress, the result (5.42) for q in conjunction with (5.27) gives 114  T → (αv−1 + ˆαv−5/3)B − pI, αv−1B + ˆαv−4/3 (αv−1 + α∗v−1)B − ¯pI, √ B − pI, if α∗ → ∞ at finite ˆα, if ˆα = α∗v2/3, if ˆα → ∞ at finite α∗. (5.43) For the case of ˆα → ∞ at finite α∗ the scalar ¯p = p− ˆαv−1 +α∗v−1I1/3. Since however p is arbitrary the extra term in ¯p can be absorbed into p. It is to be emphasized that the internal balance theory considered here need not be tied to the swelling phenomena. In other words the isochoric material case (formally the specialization of all previous results to v = 1) is a viable material description in its own right. Indeed the previous works Demirkoparan et al. (2014); Demirkoparan and Pence (2015a) considered only this case. Thus it is worth remarking that the results (5.43) for the isochoric specialization become  T → (α + ˆα)B − pI, √ B − pI, αB + ˆα (α + α∗)B − pI, if α∗ → ∞ at finite ˆα, if ˆα = α∗, if ˆα → ∞ at finite α∗. (5.44) Comparing (5.43) with (1.15)2 shows that both extreme cases of ˆα/α∗ → 0 and ˆα/α∗ → ∞ give neo-Hookean behavior in the limit. Specifically this correspondence with the neo-Hookean hyperelastic theory with swelling requires the identifications α + ˆα v2/3 α + α∗ µ ↔ if α∗ → ∞ at finite ˆα, if ˆα → ∞ at finite α∗. (5.45) In contrast, as shown by the case ˆα = α∗, the general case gives behavior that does 115 not coincide with the hyperelastic neo-Hookean theory. 116 5.2 Homogeneous deformation for the internally balanced material model Equation (4.3) provided an explicit relation for Plat in terms of v for the hyperelastic swelling theory. We now inquire into the extent to which a similar relation can be obtained for the internally balanced material theory. To this end it is necessary to decompose F in (4.3) according to (1.3). By virtue of symmetry it is presumed that ˆF and F∗ can be written in terms of positive parameters ˆξ and ξ∗ in the forms F∗ = ξ∗(e1⊗e1 +e2⊗e2)+ v ξ∗2 e3⊗e3, ˆF = ˆξ(e1⊗e1 +e2⊗e2)+ e3⊗e3. (5.46) 1 ˆξ2 These forms are consistent with the constraints (5.1). They are also consistent with the decomposition (1.3) provided that ˆξξ∗ = λlat. (5.47) Also ˆξ and ξ∗ may vary with v even though λlat is fixed. Using (4.3)1 two nontrivial scalar equations then emerge from the diagonal entries of the stress tensor (5.27): −Plat = α v λ2 lat + ˆξ2 − p, ˆα v 0 = αv λ4 lat + ˆα v ˆξ4 − p. (5.48) Two nontrivial scalar equations also emerge from the diagonal entries of the internal balance (5.31): α∗λ2 lat − ˆα ˆξ4 − qv ˆξ2 = 0, α∗v2 λ4 lat − ˆα ˆξ8 − qv ˆξ4 = 0. (5.49) For given λlat and v > λ3 lat equations (5.47)-(5.49) provide five scalar relations among the five unknown quantities: Plat, ˆξ, ξ∗, p and q. As shown below, this allows solutions for Plat, ˆξ, ξ∗, p and q in terms of v and λlat. In particular this then gives Plat as a 117 function of v and λlat such that this Plat represents the lateral confining pressure for the case v > λ3 lat. In order to verify the above claim, eliminate q between the two equations of (5.49) so as to obtain α∗v2 ˆξ8 + ˆαλ4 lat ˆξ6 − α∗λ6 lat ˆξ2 − ˆαλ4 lat = 0. (5.50) For a given pair (λlat, v) with λlat > 0 and v > 0 this polynomial equation in ˆξ has a unique positive root by virtue of Descarte’s rule of signs. This defines a function ˆΞ(λlat, v) such that ˆξ = ˆΞ(λlat, v) solves (5.50). This function also depends upon the ratio ˆα/α∗ (we do not show this in the argument list for ˆΞ). Recall also that it is permissible to allow either ˆα or α∗ to depend upon v provided that these moduli remain positive for all v. In such a case the dependence of ˆΞ upon v must also take this moduli dependence into account. Such a dependence has no effect upon the gen- erality of the results given here. Eliminating p between the two equations of (5.48) gives Plat = α − λ2 lat v + ˆα v (cid:17) − ˆξ2(cid:17) (cid:16) 1 ˆξ4 . (5.51) Using (5.47) and (5.50) it follows that an alternative form for this result is (cid:16) v λ4 lat (cid:16) v λ4 lat (cid:17) + α∗(cid:16) 1 ξ∗4 − ξ∗2(cid:17) − λ2 lat v Plat = α . (5.52) This function (5.51) with ˆξ = ˆΞ(λlat, v) defines a function (cid:16) v λ4 lat (cid:17) (cid:16) Π(λlat, v) = α − λ2 lat v + ˆα v 1 ˆΞ4(λlat, v) − ˆΞ2(λlat, v) (5.53) (cid:17) such that Plat = Π(λlat, v) gives the lateral confining pressure. The result is gen- eral so as to accommodate any dependence of α, ˆα or α∗ upon v provided that the 118 requirements α ≥ 0, ˆα ≥ 0 or α∗ ≥ 0 always hold. Because lateral confinement occurs for v > λ3 lat, it follows that v > λ3 lat defines the physically meaningful range for the argument pair (λlat, v). The value λlat = v1/3 corresponds to the free swelling threshold at which the plug has swollen by an amount that is just sufficient to allow all-around contact with the inner channel wall of the confining rigid pipe. This motivates the observation that if λlat = v1/3 then the posi- tive root ˆξ to (5.50) is ˆξ = 1. Thus ˆΞ(v1/3, v) = 1. Use of this result in (5.53) gives Π(v1/3, v) = 0. This is the expected result, namely Plat = 0 at initial material contact with the wall, i.e., when v = λ3 lat. The general dependence of confining pressure Plat on swelling v when v exceeds the contact value λ3 lat then follows from Plat = Π(λlat, v) using (5.50) and (5.53). The detailed dependence of Π(λlat, v) upon v is sensitive to how the moduli α, ˆα and α∗ vary with v. Figure 5.1 then shows Plat vs. v as determined from solution to (5.50) and (5.53) using chosen parameter values. The initial slope of the curves in Figure 5.1 can be determined analytically. For this purpose let us suppose that the moduli α, ˆα and α∗ are all independent of v. Now (5.53) can be differentiated with respect to v to get ∂Π(λlat, v) ∂v = α (cid:16) 1 λ4 lat (cid:16) (cid:17) − ˆα v2 − λ2 lat v2 1 ˆΞ4(λlat, v) (cid:17) − ˆΞ2(λlat, v) (cid:16) − 2ˆα v 2 ˆΞ5(λlat, v) + ˆΞ(λlat, v) (cid:17)∂ ˆΞ(λlat, v) ∂v . (5.54) As pointed out in the previous paragraph, ˆΞ(v1/3, v) = 1 when λlat = v1/3. Therefore ˆΞ(λlat, λ3 lat) = 1. Moreover (5.50) can be differentiated with respect to v and then evaluated at (λlat, λ3 lat) to get 119 Figure 5.1 Confinement pressure Plat = Π(λlat, v) as a function of swelling v taking ˆα = 1, α∗ = 2 and Ro = 1 (so that λlat = Rc). The three graphs on the left are for Rc = 1.51/3. The three graphs on the right are for Rc = 21/3. For each Rc the graph of Plat is given for three different values of α. Here Plat is normalized by α + (ˆαα∗)/(ˆα + α∗). Because of (5.56) this accounts for the common slope at the rightmost departure (Rc = 21/3) being less than the common slope at the leftmost departure (Rc = 1.51/3). 120 11.522.500.10.20.30.40.50.6SwellingNormalized Confinement Pressure α=0α=0.5α=1 −α∗ lat + ˆαλlat) 3(α∗λ3 . = lat) ∂ ˆΞ ∂v (cid:12)(cid:12)(cid:12)(λlat,λ3 (cid:12)(cid:12)(cid:12)(λlat,λ3 lat) ∂Π ∂v Consequently use of (5.55) in (5.54) and algebraic manipulations give (cid:16) = 2 λ4 lat α + ˆαα∗ lat + ˆα α∗λ2 (cid:17) . (5.55) (5.56) This in turn describes the slope in Figure 5.1 when the pressure Plat departs from zero due to first wall contact. In order to recover the hyperelastic result, recall that (5.45) indicates how the hyperelastic neo-Hookean swelling theory follows from two special limits of the inter- nally balanced material theory. Thus it must be the case that (4.3) will follow from (5.50) and (5.51) in these two special limits. We now verify that this is indeed the case. Consider first the limit α∗ → ∞ at finite ˆα. It then follows that the dominant balance in (5.50) gives ˆξ → λlat v1/3 as α∗ → ∞, so that (5.47) gives ξ∗ → v1/3. Thus (5.51) gives Plat → (α + ˆα v2/3 ) (cid:16) v λ4 lat (cid:17) , − λ2 lat v (5.57) (5.58) which matches (4.3) under the identification connection from (5.45). Second consider the limit ˆα → ∞ at finite α∗. It then follows that the dominant balance in (5.50) gives ˆξ → 1 as ˆα → ∞ (5.59) 121 so that (5.47) gives ξ∗ → λlat. Because (5.51) would generate the indeterminate form of infinity multiplied by zero we use the alternative form (5.52) to obtain Plat → (α + α∗)(cid:0) v (cid:1) − λ2 lat v λ4 lat (5.60) which matches (4.3) under the correspondence from (5.45). 122 5.3 Cylindrical deformation in internally balanced material model In this section we use the internally balanced material model for the same boundary value problem of the confined hollow tube that was considered in section 4.3. In this regard the cylindrical deformation (4.4) that is subject to the volume change (4.6) is again considered. Recall that the tube is still confined such that the (4.8), (4.9) and (4.10) still hold. The relation (4.9) again provides that ri = ri(λz) and this makes the axial force requirement (4.13) be a single relation to obtain the ultimate unknown λz. Here in this section in the internally balanced material model that is given by (5.26) for simplicity it is assumed that α = 0. The left Cauchy-Green deformation tensor B is given by B = FFT =  (cid:48)2 r 0 0 0 (r/R)2 0 0 0 λ2 z  = r (cid:48)2(er⊗er)+(r/R)2(eθ⊗eθ)+λ2 z(ez⊗ez). (5.61) Moreover, notice that it follows from (5.28) that an eigenvector of B is also an eigen- vector of ˆB and this argument is shown in Demirkoparan et al. (2014). Hence it is reasonable to seek solutions in which ˆB is symmetric positive-definite in the form  ˆB = ˆBrr 0 0 0 ˆBθθ 0 0 0 ˆBzz  = ˆBrr(er ⊗ er) + ˆBθθ(eθ ⊗ eθ) + ˆBzz(ez ⊗ ez), (5.62) where its components satisfy the internal balance equations (5.31). It then follows 123 from (5.62) and (5.27) that the stress tensor T is of the form  = Trr(er ⊗ er) + Tθθ(eθ ⊗ eθ) + Tzz(ez ⊗ ez) 0 0  T = where Trr 0 Tθθ 0 0 0 Tzz (5.63) (5.64) Trr = Tθθ = Tzz = ˆBrr − p, ˆBθθ − p, ˆBzz − p. ˆα v ˆα v ˆα v This will reduce the equilibrium equation into (4.11) that is associated with the boundary condition (4.12). The boundary value problem for the deformation field (4.4) with given swelling amount v, the geometric configuration Ri, Ro and Rc and the stress field (5.64) is then formulated to seek the seven unknowns as follows: axial elongation λz, the outer confinement pressure Plat, scalar functions p and q and a symmetric positive-definite tensor function ˆB (with the three unknown component functions of ˆBrr, ˆBθθ and ˆBzz) such that the following equations are satisfied: the equilibrium equation (4.11), the zero resultant axial force , the requirement (5.1) and three scalar equations of Internal Balance that come from diagonal entries in (5.28). Once λz is determined then ri is obtained directly from (4.9). 5.3.1 Internal balance requirement We consider an approach in which q and ˆBrr, ˆBθθ and ˆBzz are regarded as the primary and λz, Plat and the hydrostatic pressure function p are regarded as the secondary unknowns. Notice that instead of q, we seek for ρ = qv/ˆα. Also because we assume finite and positive values for ˆα and α∗ we choose to pick the notations in (5.33). The 124 primary unknowns are then formulated in terms of the secondary unknowns. Once the secondary unknowns are solved, the primary unknowns are consequently obtained. To formulate the four primary unknowns we use the three equations of Internal Balance (5.31) and the requirement of (5.1). Thus using the assumptions (5.61) and (5.62) in the internal balance (5.31), it gives three quadratic and decoupled equations relating ˆBrr, ˆBθθ and ˆBzz to ρ as follows (cid:48)2/β = 0, rr + ρ ˆBrr − r ˆB2 θθ + ρ ˆBθθ − r2/(βR2) = 0, ˆB2 zz + ρ ˆBzz − λ2 ˆB2 z/β = 0. (5.65) These equations can be solved for positive ˆBrr, ˆBθθ and ˆBzz. Making one additional replacement of r (cid:48) = Rv/(rλz) from (4.6) and introducing s = r/R these solutions are , (5.66) (cid:19) (cid:18) (cid:18) (cid:18) − ρ +(cid:112)ρ2 + 4v2/(βs2λ2 (cid:19) − ρ +(cid:112)ρ2 + 4s2/β (cid:19) − ρ +(cid:112)ρ2 + 4λ2 . , z/β z) ˆBrr = ˆBθθ = ˆBzz = 1 2 1 2 1 2 Now we need one more equation in order to solve for the four primary unknowns. To this end, det ˆB = 1 requires that ˆBrr ˆBθθ ˆBzz = 1 which because of (5.66) becomes − ρ +(cid:112)ρ2 + 4v2/(βs2λ2 z) (cid:19)(cid:18) − ρ +(cid:112)ρ2 + 4s2/β (cid:19)(cid:18) − ρ +(cid:112)ρ2 + 4λ2 z/β (cid:19) (cid:18) 1 8 = 1. (5.67) This brings additional relation between ρ, s, v, λz and β and shows that ρ is indepen- dent of θ and z. The immediate question then arises if there is any unique ρ which satisfies (5.67). To investigate that let’s define h(ρ) = 1 8 h1(ρ)h2(ρ)h3(ρ), (5.68) 125 in which h1(ρ) = −ρ +(cid:112)ρ2 + 4R2v2/(βr2λ2 h2(ρ) = −ρ +(cid:112)ρ2 + 4r2/(βR2) > 0, h3(ρ) = −ρ +(cid:112)ρ2 + 4λ2 z/β > 0, z) > 0, (5.69) and seek the solution ρ to h(ρ) = 1 (5.70) which is equivalent to (5.67). To this end observe that 8 dh dρ (cid:48) (cid:48) (cid:48) 2h3 + h1h2h 1h2h3 + h1h = h 3 with 1 = (−1 + (cid:48) h 2 = (−1 + (cid:48) h 3 = (−1 + (cid:48) h ρ (cid:112)ρ2 + 4R2v2/(βr2λ2 (cid:112)ρ2 + 4r2/(βR2) ρ(cid:112)ρ2 + 4λ2 ) < 0. ρ ) < 0, z) ) < 0, z/β (5.71) (5.72) This makes the channel argument dh/dρ < 0 for all real ρ. Note also from (5.67) that ρ→−∞h(ρ) = ∞, lim lim ρ→∞h(ρ) = 0. (5.73) It is concluded that the equation h(ρ) = 1 has a unique solution ρ which is a function of s = r/R and also the parameters v, λz, β and Ri. We recall that all of these parameters are given with the exception of λz. It follows that we can write ρ = ρ(s; v, β, λz). (5.74) To determine the sign of ρ notice that from (5.68) we have h(ρ)|ρ=0 = v/β3/2. Hence, for the solution of h(ρ) = 1, it is seen that v/β3/2 < 1 gives ρ < 0 and v/β3/2 > 1 126 gives ρ > 0. In other words the sign of ρ which solves (5.67) depends on the value of v/β3/2 = v(α∗/ˆα)3/2. Hence, in any given problem ρ will either be restricted to vary through only positive values or through only negative values. If this ρ is taken such that ρ = ρ(s; v, β, λz) solves (5.67), then from (5.66), the components of ˆB are found to be functions of s. With the help of the map r = r(R) from (4.4), the function ρ and all the components of ˆB are obtained and found to be functions of s and the still unknown parameters λz and ri(λz). Depending on which aspect is to be emphasized we either write ρ = ρ(s) or ρ = ρ(s, λz). The analytical solution to (5.67) is not available. Two numerical solutions of ρ = ρ(s; v, β, λz) are graphed in the Figure 5.2. It is seen that for β < v2/3 the solution ρ(s) is positive and if β > v2/3 the solution ρ(s) is negative. Moreover, the solutions are finite over the channel range of s, it goes to zero when s → 0 and asymptotes to zero when s → ∞. In order to identify the special value of the parameter s = sext that makes the extremum of ρext = ρ(sext) we can take the derivatives of the two sides of (5.67) with respect to s and use dρ/ds(cid:12)(cid:12)s=sext (cid:113) (cid:113) (cid:18) 1 s4 ext = 0. This yields (cid:19) (cid:18) extλ2 z) − ρ + (cid:113) ρ2 ext + 4s2 ext/β ρ2 ext + 4s2 (cid:19) One obvious root of this equation for sext is sext = (cid:112)v/λz that is shown in Figure ρ2 ext + 4v2/(βs2 ρ2 ext + 4v2/(βs2 − ρext + ext/β = extλ2 z) (5.75) (cid:113) λ2 z v2 5.2. 127 Figure 5.2 Three solutions of ρ(s; v, β, λz) to the requirement (5.67) with one repre- sentative set of v = 2 and λz = 2. Two of the non-trivial solutions are obtained for different ration β = 2 and β = 1/2. If β is taken as β = v2/3 = 22/3 then the obvious solution is ρ(s) ≡ 0. 5.3.2 Stress and equilibrium To find the secondary unknowns λz, Plat and p, we first show that p is a function of R alone. We notice that since the components of ˆB are independent of θ and z, regarding (5.64) and equilibrium equations from (4.11) the components of T and consequently p, are also independent of θ and z. Hence it follows that p = p(R) and by presuming that the map r = r(R) is invertible we can write the function p as a function of s = r/R and hence p = p(s). The function p(s) will then be solved from the first equilibrium equation (4.11) and satisfying the boundary condition (4.12) where T is given by (5.27). The parameter ri(λz) is obtained via (4.9) and recall that the unknown parameter λz is the solution to the requirement (4.13). Thus the procedure is as follows: after determining ρ(s, λz) the stress components to be used 128 in (4.11) are from (5.27). If we rewrite the stress components from (5.64) it follows Trr = Tθθ = Tzz = ˆBrr − p(s), ˆBθθ − p(s), ˆBzz − p(s). ˆα v ˆα v ˆα v Thus the components of ˆB are used from (5.66) into (5.76) and it gives (cid:0) − ρ(s, λz) +(cid:112)ρ(s, λz)2 + 4v2/(βs2λ2 (cid:0) − ρ(s, λz) +(cid:112)ρ(s, λz)2 + 4s2/β(cid:1) − p(s), (cid:0) − ρ(s, λz) +(cid:112)ρ(s, λz)2 + 4λ2 z/β(cid:1) − p(s). z)(cid:1) − p(s), Trr = Tθθ = Tzz = ˆα 2v ˆα 2v ˆα 2v (5.76) (5.77) (5.78) (5.79) (5.80) (5.81) Hence (4.11) becomes the differential equation (cid:18)(cid:112)ρ(s, λz)2 + 4v2/(βs2λ2 z) −(cid:112)ρ(s, λz)2 + 4s2/β (cid:19) = 0, (5.82) ∂ ∂r Trr + ˆα 2rv and notice that with the help of the connection ∂ ∂r Trr = 1 (cid:48) r d dR Trr, into (5.82) it follows that R(cid:90) R(cid:90) (cid:18)(cid:112)ρ(s, λz)2 + 4s2/β −(cid:112)ρ(s, λz)2 + 4v2/(βs2λ2 z) dTrr = ˆα 2rv (5.83) (cid:19) r (cid:48) dR. (5.84) Ri Ri The left hand side of the relation (5.85) is reduced with the help of the boundary condition (4.12) and in the right hand side we can use (4.6) to replace r(cid:48) and perform 129 the integration to obtain Trr(s) = − ˆα 2 s(cid:90) si φ(ρ(s, λz), v, λz, β) ds, (5.85) 1 vs − λzs3 φ(ρ(s, λz), v, λz, β) = in which it is defined that si = ri/Ri and also (cid:18)(cid:112)ρ(s, λz)2 + 4v2/(βs2λ2 z) −(cid:112)ρ(s, λz)2 + 4s2/β bound (cid:112)v/λz) or increasing (with upper bound (cid:112)v/λz) for Ri ≤ R ≤ Ro. Since concluded that si < so. Consequently it follows that si < so <(cid:112)v/λz. This guaran- the boundary value problem is defined such that Pi = 0 and Po = Plat > 0 then it is It can be confirmed that the parameter s = r/R is either decreasing (with lower (cid:19) . (5.86) tees a non-zero denominator in the function φ in (5.86) and it allows the integration (5.90) to be determined. Once the radial component of the Cauchy stress tensor is obtained via (5.85), the hydrostatic pressure p(s) follows from (5.79) and other stress components are also given by (5.80) and (5.81). We have (cid:18)(cid:112)ρ(s, λz)2 + 4s2/β −(cid:112)ρ(s, λz)2 + 4v2/(βs2λ2 s(cid:90) z) φ(ρ(s, λz), v, λz, β) ds, (cid:19) (cid:19) (5.87) (5.88) Tθθ(s) = Tzz(s) = ˆα 2v − ˆα 2 si ˆα 2v − ˆα 2 z/β −(cid:112)ρ(s, λz)2 + 4v2/(βs2λ2 z) (cid:18)(cid:112)ρ(s, λz)2 + 4λ2 s(cid:90) φ(ρ(s, λz), v, λz, β) ds. si 130 Figure 5.3 Axial elongation λz on a hollow tube as a function of swelling v taking α = 0 ˆα = 1, α∗ = 2 and Ri = 1/2, Ro = 1 (so that λlat = Rc). The graphs are the solution to (5.89) and are for two separate cases of outer pipe radius: Rc = 1.51/3 and Rc = 21/3, which are chosen so as to give contact v values of 1.5 and 2, respectively. Moreover, λz is the only unknown in all above relations because the parameter ro is prescribed as ro = Rc for v > λ3 lat. To solve for λz recall that the requirements of volume conversation (4.9) and also the zero axial resultant load (4.13) with the axial component of Cauchy stress (5.85) are used. The latter is λlat(cid:90) ri/Ri (cid:12)(cid:12)(cid:12)(cid:12)s= r R 2π Tzz(s) rdr ≡ 0, (5.89) in which the definition so = ro/Ro = Rc/Ro = λlat is again used and recall that ri(λz) is given by (4.9). Two different numerical solutions of (5.89) are graphed in Figure 5.3. As it is seen the solution is consistent with conventional M-R model that was obtained and shown in Figure 4.4. This leads to obtaining the inner radius from (4.9) that is shown in Figure 5.4. Once λz and ri(λz) are determined the lateral confinement pressure follows from Trr|λlat= −Plat which is equivalent to 131 Figure 5.4 Inner radius ri on a hollow tube as a function of swelling v taking α = 0 ˆα = 1, α∗ = 2 and Ri = 1/2, Ro = 1 (so that λlat = Rc). The graphs are the solution to (4.9) and are for two separate cases of outer pipe radius: Rc = 1.51/3 and Rc = 21/3, which are chosen so as to give contact v values of 1.5 and 2, respectively. λlat(cid:90) si Plat = ˆα 2 φ(ρ(s, λz), v, λz, β) ds. (5.90) The values of lateral pressure are also plotted in Figure 5.5. This plot also shows the variation of lateral pressure with respect to changing the material parameter ratio β = ˆα/α∗. One special curve is for β → ∞ that is obtained with very large values of ˆα at a finite α∗. Recall from (5.45) that in this limit the finite material parameter α∗ retrieves the elastic modulus µ. This figure shows how the lateral pressure also recovers the hyperelastic behavior that was shown in Figure 4.7. 5.3.3 Explicit solution of the special case v = β3/2 We recall the special case of the relation (5.66) when v = β3/2 for which the solution ρ(s) is identical to zero for all si ≤ s ≤ so and this corresponds to the second special case introduced in (5.42). Note that in this case the axial stress component 132 Figure 5.5 Confinement pressure Plat on a hollow tube as a function of swelling v taking α = 0 and β = 1/2 (ˆα = 1, α∗ = 2) and Ri = 1/2, Ro = 1 (so that λlat = Rc) for longer graphs with Rc = 1.51/3 and Rc = 21/3. The graphs of β = 4, 10 and ∞ are also plotted for Rc = 1.51/3. Note that the graph of β = ∞ retrieves the hyperelastic behavior shown in Figure 4.6 based on (5.45). Here Plat is normalized by ˆαα∗/(ˆα+α∗). given by (5.88) simplifies to Tzz(s) = ˆα β1/2 (cid:18) λz β3/2 − 1 λzsi (cid:19) , (5.91) which is independent of the normalized radius s. The requirement of the zero axial load invokes that this radial stress is zero for all pointwise radial locations. Note that this result is unlike the M-R behavior that was described in Chapter IV and also in contrast with the internally balanced material model with swelling amounts v (cid:54)= β3/2 where the pointwise zero axial stress could not be captured. This pecu- liar behavior shows that with the internally balanced material, the special value of swelling makes the axial stress field vanishes while the assumption z = λzZ still holds. 133 The condition (5.91) provides the relation λz = β3/4 s1/2 i , (5.92) that in conjunction with (4.9) obtains λz. It follows that λz is the solution to the single equation that is a forth order polynomial z − β3/2 λ4 λ2 lat (1 − ζ 2)λ3 z − β3ζ 2 λ2 lat = 0. (5.93) Since there is only one sign change to the coefficients of the polynomial it follows that there exists only one root to (5.93) which is the solution to λz. Note that once λz is obtained from (5.93) then the radial stress from (5.85) simplifies to φ(ρ(s, λz), v, λz, β) ds φ(0, v, λz, β) ds (cid:18) 2v sλzβ1/2 − 2s β1/2 (cid:19) ds (5.94) 1 vs − λzs3 Trr(s) = − ˆα 2 = − ˆα 2 = − ˆα 2 si s(cid:90) s(cid:90) s(cid:90) s(cid:90) si si s(cid:90) si ds vs2 − λzs4 ˆα β1/2 si ˆα = = ds v − λzs2 − ˆαv (cid:19) (cid:18)1 λzβ1/2 − λ2 z β3/2 . λzβ1/2 s It is clear that with (5.92) this radial stress vanishes at s = si and it is decreasing (compression) toward the outer radius at which so = λlat and the lateral pressure is given by Plat = ˆα λzβ1/2 (cid:18) λ2 β3/2 − 1 z λlat (cid:19) . (5.95) 134 If this pressure is normalized by ˆαα∗/(ˆα + α∗) and is renamed to ¯Plat it follows that (cid:18) λ2 β3/2 − 1 z ¯Plat = β + 1 λzβ1/2 (cid:19) λlat . (5.96) In addition, the circumferential stress form (5.87) reduces to Tθθ(s) = ˆα β2 (s − λz) . (5.97) Since si ≤ s ≤ so ≤ λlat ≤ λz this hoop stress is also compressive. Moreover note that in a case where this special case coincides with the first wall contact we have s = λlat, lat and λz = λlat for all si ≤ s ≤ so and hence all three stress components vanish. v = λ3 135 CHAPTER VI Concluding Remarks, Speculations and Broader Connections In this thesis we have shown how swelling can be treated in the framework of hyper- elastic theory as a parametric volume change. We used the multiplicative decompo- sition of the deformation gradient to distinguish between swelling and other elastic accommodations within different processes. This allows for introducing the incorpo- ration of swelling in the constitutive relations both in the conventional hyperelastic and in the internal balance theory. In the context of conventional hyperelastic consti- tutive model we considered spherical and cylindrical geometries and in the former case allowed the material parameters to depend on the swelling amount. In the context of the internally balance material theory we considered the cylindrical deformation and showed how swelling can lead to different response than that of the conventional hyperelastic theory. The possibility of modeling a variety of swelling effects has immediate conse- quences in several areas of application, many of which are motivated by biomechani- cal modeling. Others have direct bearing on possible mechanical devices that achieve their control by soft material action. The combined effect of external loading and internal swelling can give rise to complicated states of deformation. Even in the simple setting of a spherical shell subject to combinations of simple pressure with uniform through-thickness (homogeneous) swelling, instabilities can arise that might 136 not otherwise be present if either pressurization or swelling was acting by itself. In this thesis we have explored the role that swelling can have on eliciting qualitative changes in the pressure-expansion inflation response of spherical deformation. Unlike in the context of conventional hyperelastic theory where the swelling is imposed as prescribed deformation, the internally balanced theory makes use of a deformation gradient decomposition to obtain the portion of deformation that is responsible for swelling. The additional requirement that obtains the swelling deformation is a nat- ural consequence of energy minimization with respect to the new decomposition. We explored aspects of these two constitutive theories for cylindrical deformation of con- fined plug and tube. In the first and second chapters considerations have been limited to the deforma- tions of the spherical symmetry. Generalizing methods of analysis pioneered by Car- roll in the context of incompressible hyperelasticity we have examined a rather straight forward constitutive model, one which is motivated by the well known Mooney-Rivlin model in the incompressible theory, so as to incorporate swelling dependent stiffness parameters. We have shown how certain dependencies preserve the overall qualita- tive nature of the inflation process independent of the amount of swelling, whereas other dependencies do not. In the latter case we have provided general rules, il- lustrated with examples, showing how certain constitutive forms cause a monotonic (benign) inflation response in the absence of swelling to become nonmonotonic (burst- inducing) as the swelling proceeds. Alternative constitutive forms have the opposite effect, burst-inducing inflation response in the absence of swelling can be mitigated into benign inflation response as the swelling proceeds. It is also important to note that the development of surface roughness due to swelling has been observed in solid hydrogel spheres. This may then give way to an aspherical and facetted surface mor- phology as the swelling proceeds Bertrand et al. (2016). These are typically confined 137 to distinct ranges of overall swelling and may be connected to nonuniform states of internal hydration (v = v(R) in the notation of the third chapter), especially if the identified swelling ranges differ on the basis of whether the overall fluid content is increasing or decreasing. It is interesting to speculate in the possible broader connections of this work to physical and biological phenomena. The overall considerations of the present study, as well as possible future studies that bring to bear the techniques in the above referenced works, give rise to the prospect that swelling, when viewed as a control variable, could be manipulated so as to tune the inflation response of spherical shells and membranes. This includes the possibility of both triggering and avoiding in- stances of inflation burst. On this basis, one may even speculate to what extent such processes of regulation might be present in biological systems. For example, colonies of soft celled creatures are capable of rapidly undergoing complex shape changes. This includes the green alga volvox in the shape of a spherical shell. At a crucial point in their embryonic development, volvox essentially turn themselves inside-out in a process that is conjectured to be triggered by cell shape change at a specific latitude on the shell H¨ohn et al. (2015); Haas and Goldstein (2015). An intriguing issue in this context is the extent to which the global conditions of the type examined here might possibly abet the resulting snap-through process. In the final two chapters we have examined the finite strain swelling of cylindrical plugs (with original outer radius Ro) within a rigid cylindrical tube of radius Rc > Ro, all with circular cross-section. For a simple plug with no internal cavities or voids the wall contact deformation involves only a simple axial lengthening. When the plug is not simply connected, in our case because of an internal channel of original radius Ri, the deformation after contact involves a combination of axial lengthening and channel closing. A boundary value problem must be solved in order 138 to determine the relative importance of these two ways to accommodate the swelling volume increase. For a conventional swelling model based on a generalized neo- Hookean response we show that this boundary value problem has a unique solution. Asymptotic analysis of the problem shows how the channel radius closing depends upon v as the swelling becomes large. Figure 4.8 shows the various trends for the plug’s internal channel radius as the swelling proceeds. Prior to wall contact the internal channel radius grows linearly with v causing ri to reach its maximum size when the channel makes contact with the rigid tube wall. Additional swelling then constricts the internal channel in a manner that can be addressed analytically both for the immediate post-contact values of v and then again for sufficiently large v. While these trends are intuitive for the problem formulated here they can also provide useful insight for various hypothetical elastic systems involving tubes with multiple layers. One reason for pointing out these connections is because certain biological tube organs have such a multi-layer structure, indeed, for the purpose of mechanics mod- eling, arteries are often regarded as three layered tubes Holzapfel et al. (2000). So is the windpipe (the trachea) Gou and Pence (2016). Moreover, it is typically the inner layer that is the most subject to distress and disease, thus exposing it to various possibilities for inflammation, edema or other phenomena that at least have some aspect of swelling. The analogy with our above hypothetical tube is however a highly imperfect one – it would typically be the swelling layer that is the least stiff and not the central layer. In any event, the type of formalism employed here might at least provide some insight into the behavior of such biological systems when subject to inflammatory swelling. Other soft tissue cylindrical structures in the body that are subject to swelling, although not in the sense of inflammation but rather in the sense of normal function, include the cervix, especially as it undergoes remodeling during pregnancy Myers et al. (2015). As indicated by the above references, all of these 139 systems are the subject of current biomechanics based modeling. Finally brain swelling (hydrocephalus), while being confined to more of a spherical (not tubular) enclosure, is also the object of recent mechanics based modeling Wilkie et al. (2011). In this case the three layer system identification is again interesting, if we identify Ri < R < Ro with the swelling brain, and R > Rc with the surrounding skull. Detailed finite strain treatments relating to this issue are now being developed. The approach presented in this thesis could aid in the formulation of more precise mathematical models of this important human health concern. 140 BIBLIOGRAPHY 141 BIBLIOGRAPHY Alexander, H. (1971), Tensile instability of initially spherical balloons, International Journal of Engineering Science, 9 (1), 151–160. Ateshian, G. A. (2007), On the theory of reactive mixtures for modeling biological growth, Biomechanics and Modeling in Mechanobiology, 6, 423–445. Baek, S., and T. J. Pence (2011), Inhomogeneous deformation of elastomer gels in equilibrium under saturated and unsaturated conditions, Journal of the Mechanics and Physics of Solids, 59 (3), 561–582. Baek, S., and A. R. Srinivasa (2004), Diffusion of a fluid through an elastic solid un- dergoing large deformation, International Journal of Non-linear Mechanics, 39 (2), 201–218. Ben Amar, M., and P. Ciarletta (2010), Swelling instability of surface-attached gels as a model of soft tissue growth under geometric constraints, Journal of the Mechanics and Physics of Solids, 58 (7), 935–954. Ben Amar, M., and A. Goriely (2005), Growth and instability in elastic tissues, Journal of the Mechanics and Physics of Solids, 53 (10), 2284–2319. Bertrand, T., J. Peixinho, S. Mukhopadhyay, and C. W. MacMinn (2016), Dynamics of swelling and drying in a spherical gel, Physical Review Applied, 6 (6), 064,010. Bowen, R. M. (1980), Incompressible porous medium models by the use of the theory of mixtures, International Journal of Engineering Science, 18, 1129–1148. Carroll, M. M. (1987), Pressure maximum behavior in inflation of incompressible elastic hollow spheres and cylinders, Quarterly of applied mathematics, 45 (1), 141– 154. Chen, Y.-C., and A. Hoger (2000), Constitutive functions of elastic materials in finite growth and deformation, Journal of Elasticity, 59 (1), 175–193. Chester, S. A., and L. Anand (2010), A coupled theory of fluid permeation and large deformations for elastomeric materials, Journal of the Mechanics and Physics of Solids, 58 (11), 1879–1906. Demirkoparan, H., and T. J. Pence (2007a), The effect of fiber recruitment on the swelling of a pressurized anisotropic non-linearly elastic tube, International Journal of Non-Linear Mechanics, 42 (2), 258–270. Demirkoparan, H., and T. J. Pence (2007b), Swelling of an internally pressurized nonlinearly elastic tube with fiber reinforcing, International journal of solids and structures, 44 (11), 4009–4029. 142 Demirkoparan, H., and T. J. Pence (2008), Torsional swelling of a hyperelastic tube with helically wound reinforcement, Journal of Elasticity, 92 (1), 61–90. Demirkoparan, H., and T. J. Pence (2015a), Finite stretching and shearing of an internally balanced elastic solid, Journal of Elasticity, 121 (1), 1–23. Demirkoparan, H., and T. J. Pence (2015b), Magic angles for fiber reinforcement in rubber-elastic tubes subject to pressure and swelling, International Journal of Non-Linear Mechanics, 68, 87–95. Demirkoparan, H., and T. J. Pence (2017), Swelling–twist interaction in fiber- reinforced hyperelastic materials: the example of azimuthal shear, Journal of En- gineering Mathematics, pp. 1–22. Demirkoparan, H., T. J. Pence, and H. Tsai (2014), Hyperelastic internal balance by multiplicative decomposition of the deformation gradient., Archive for Rational Mechanics & Analysis, 214 (3). Deng, H., and T. J. Pence (2010), Equilibrium states of mechanically loaded saturated and unsaturated polymer gels, Journal of Elasticity, 99 (1), 39–73. Drozdov, A., et al. (2013), Constitutive equations in finite elasticity of swollen elas- tomers, International Journal of Solids and Structures, 50 (9), 1494–1504. Drozdov, A. D. (2013), Finite elasticity of nanocomposite hydrogels, Composite In- terfaces, 20 (9), 673–692. Drozdov, A. D., and J. d. Christiansen (2013), Constitutive equations in finite elas- ticity of swollen elastomers, International Journal of Solids and Structures, 50, 1494–1504. Duda, F. P., A. C. Souza, and E. Fried (2010), A theory for species migration in a finitely strained solid with application to polymer network swelling, Journal of the Mechanics and Physics of Solids, 58 (4), 515–529. Duda, F. P., A. C. Souza, and E. Fried (2011), Solvent uptake and cavitation, Journal of the Mechanics and Physics of Solids, 59 (11), 2341–2354. Ericksen, J. L. (1975), Equilibrium of bars, Journal of Elasticity, 5 (3), 191–201. Evans, E., V. Heinrich, F. Ludwig, and W. Rawicz (2003), Dynamic tension spec- troscopy and strength of biomembranes, Biophysical Journal, 85 (4), 2342–2350. Fang, Y., T. J. Pence, and X. Tan (2011), Fiber-directed conjugated-polymer torsional actuator: nonlinear elasticity modeling and experimental validation, IEEE/ASME Transactions on Mechatronics, 16 (4), 656–664. Forest, S., K. Ammar, B. Appolaire, N. Cordero, and A. Gaubert (2014), Micro- morphic approach to crystal plasticity and phase transformation, in Plasticity and Beyond, pp. 131–198, Springer. 143 Gibbons, M. M., and W. S. Klug (2008), Influence of nonuniform geometry on nanoin- dentation of viral capsids, Biophysical Journal, 95 (8), 3640–3649. Goriely, A., D. E. Moulton, and R. Vandiver (2010), Elastic cavitation, tube hol- lowing, and differential growth in plants and biological tissues, EPL (Europhysics Letters), 91 (1), 18,001. Gou, K., and T. J. Pence (2016), Hyperelastic modeling of swelling in fibrous soft tissue with application to tracheal angioedema, Journal of mathematical biology, 72 (1-2), 499–526. Graf, J., M. Rupnik, G. Zupancic, and R. Zorec (1995), Osmotic swelling of hepato- cytes increases membrane conductance but not membrane capacitance., Biophysical Journal, 68 (4), 1359. Green, A. E., and R. T. Shield (1950), Finite elastic deformation of incompressible isotropic bodies, in Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, vol. 202, pp. 407–419, The Royal Society. Haas, P. A., and R. E. Goldstein (2015), Elasticity and glocality: initiation of embryonic inversion in volvox, Journal of The Royal Society Interface, 12 (112), 20150,671. Hadoush, A., H. Demirkoparan, and T. J. Pence (2014), Updated lagrange formula- tion of internally balanced compressible hyperelastic materials, in Qatar Foundation Annual Research Conference, 1, p. HBPP0271. Hadoush, A., H. Demirkoparan, and T. J. Pence (2015), A constitutive model for an internally balanced compressible elastic material, Mathematics and Mechanics of Solids, p. 1081286515594657. Hadoush, A., H. Demirkoparan, and T. J. Pence (2017), Finite element analysis of internally balanced elastic materials, Computer Methods in Applied Mechanics and Engineering. Havner, K. S. (1992), Finite plastic deformation of crystalline solids, Cambridge University Press. H¨ohn, S., A. R. Honerkamp-Smith, P. A. Haas, P. K. Trong, and R. E. Goldstein (2015), Dynamics of a volvox embryo turning itself inside out, Physical Review Letters, 114 (17), 178,101. Holzapfel, G. A., T. C. Gasser, and R. W. Ogden (2000), A new constitutive frame- work for arterial wall mechanics and a comparative study of material models, Jour- nal of Elasticity, 61, 1–48. Hong, W., X. Zhao, J. Zhou, and Z. Suo (2008), A theory of coupled diffusion and large deformation in polymer gels, Journal of the Mechanics and Physics of Solids, 56, 1779–1793. 144 Kr¨oner, E. (1959), Allgemeine kontinuumstheorie der versetzungen und eigenspan- nungen, Archive for Rational Mechanics and Analysis, 4 (1), 273–334. Kyriakides, S., and Y.-C. Chang (1990), On the inflation of a long elastic tube in the presence of axial load, International Journal of Solids and Structures, 26 (9-10), 975–991. Lee, E. H. (1969), Elastic-plastic deformation at finite strains, ASME. Li, F., C. U. Chan, and C. D. Ohl (2013), Yield strength of human erythrocyte membranes to impulsive stretching, Biophysical Journal, 105 (4), 872–879. Markert, B., B. Monastyrskyy, and W. Ehlers (2008), Fluid penetration effects in porous media contact, Continuum Mechanics and Thermodynamics, 20, 303=315. McMahon, J., A. Goriely, and M. Tabor (2010), Spontaneous cavitation in growing elastic membranes, Mathematics and Mechanics of Solids, 15 (1), 57–77. M¨uller, I., and P. Strehlow (2004), Rubber and rubber balloons: paradigms of thermo- dynamics, vol. 637, Springer Science & Business Media. Myers, K. M., H. Feltovich, E. Mazza, J. Vink, M. Bajka, R. Wapner, T. Hall, and M. House (2015), The mechanical role of the cervix in pregnancy, Journal of Biomechanics, 48, 1511–1523. Nagel, T., S. Loerakker, and C. W. J. Oomens (2009), A theoretical model to study the effects of cellular stiffening on the damage evolution in deep tissue injury, Computer Methods in Biomechanics and Biomedical Engineering, 12 (5), 585–597. Ogden, R. W. (1997), Non-linear elastic deformations, Courier Corporation. Pence, T. J. (2012), On the formulation of boundary value problems with the incom- pressible consituents constraint in finite deformation poroelasticity, Mathematical Methods in the Applied Sciences, 35, 1756–1783. Pence, T. J., and H. Tsai (2005a), Swelling-induced microchannel formation in non- linear elasticity, IMA journal of applied mathematics, 70 (1), 173–189. Pence, T. J., and H. Tsai (2005b), On the cavitation of a swollen compressible sphere in finite elasticity, International Journal of Non-Linear Mechanics, 40 (2), 307–321. Pence, T. J., and H. Tsai (2006), Swelling-induced cavitation of elastic spheres, Math- ematics and Mechanics of Solids, 11 (5), 527–551. Sadik, S., A. Angoshtari, A. Goriely, and A. Yavari (2016), A geometric theory of nonlinear morphoelastic shells, Journal of Nonlinear Science, 26 (4), 929–978. Selvadurai, A. P. S., and A. P. Suvorov (2016), Coupled hydro-mechanical effects in a poro-hyperelastic material, Journal of the Mechanics and Physics of Solids, 91, 311–333. 145 Stuart, M. A. C., et al. (2010), Emerging applications of stimuli-responsive polymer materials, Nature Materials, 9 (2), 101–113. Treloar, L. R. G. (1975), The physics of rubber elasticity, Oxford University Press, USA. Tsai, H., T. J. Pence, and E. Kirkinis (2004), Swelling induced finite strain flexure in a rectangular block of an isotropic elastic material, Journal of Elasticity, 75 (1), 69–89. Van der Sman, R. G. M. (2015), Hyperelastic models for hydration of cellular tissue, Soft Matter, 11 (38), 7579–7591. Vinod Kumar, K., and F. Demeke (2011), Analysis of mechanical behavior of red blood cell membrane with malaria infection, World Journal of Mechanics, 2011. Vujoˇsevi´c, L., and V. A. Lubarda (2002), Finite-strain thermoelasticity based on multiplicative decomposition of deformation gradient, Theoretical and Applied Me- chanics, (28-29), 379–399. Wilkie, K. P., C. S. Drapaca, and S. Sivaloganathan (2011), A nonlinear viscoelas- tic fractional derivative model of infant hydrocephalus, Applied Mathematics and Computation, 217, 8693–8704. Wineman, A. S., and K. R. Rajagopal (1992), Shear induced redistribution of fluid within a uniformly swollen nonlinear elastic cylinder, International Journal of En- gineering Science, 30, 1583–1595. Yavari, A., and A. Goriely (2013), Nonlinear elastic inclusions in isotropic solids, in Proc. R. Soc. A, vol. 469, p. 20130415, The Royal Society. Zamani, V., and T. Pence (2017), Swelling, inflation, and a swelling-burst instability in hyperelastic spherical shells, International Journal of Solids and Structures, 125, 134–149. 146