Swelling-induced twisting and shearing in fiber composites: the effect of the base matrix mechanical response

If a helical network of fibers is embedded in a swellable matrix, and if the fibers themselves resist swelling, then a change in the amount of swelling agent will cause a corresponding twisting motion in the material. This effect has recently been analyzed in highly deformable soft material tubes using the theory of hyperelasticity, suitably modified to incorporate the swelling effect. Those studies examined the effect of spiral angle and fiber-to-matrix inherent stiffness in the context of a ground state matrix material that exhibited classical neo-Hookean behavior in the absence of swelling. While such a ground state material is nonlinear in general, its shear response is linear. As we describe here, it is this shear response that governs the matrix contribution to the twist-swelling interaction. Because gels, elastomers, and even biological tissue can exhibit complex ground state behavior in shear—behavior that may depart significantly from a linear response—we then examine the effect of alternative ground state behaviors on the twist-swelling interaction. The range of behaviors considered includes materials that harden in shear, materials that soften in shear, materials that have an ultimate shear stress bound, and materials that collapse in shear. Matrix materials that either soften or collapse in shear are found to amplify the twisting effect.


Introduction
By embedding biassing fibers in a matrix material that is highly absorbant, it is possible to generate specialized deformation modes as the material swells. This is because the fibers naturally guide the deformation in what is essentially a prearranged way. In a previous paper [3], we analyzed this effect in tubes containing spiral patterns of fibers that lie in the plane of the tube's cross-section. Swelling of the matrix then cause an overall tube twist to take place even though no external forces are applied. Effectively, the internal swelling acts as a body force to induce the twist.
In conjunction with other modes of deformation that can be induced by swelling, such as bending [8] and torsion [1], this suggests that a wide variety of deformation modes could be engendered by intelligent pattern design for fiber biasing Hasan [2,4]. While fiber pattern biassing is inherently a design freedom that could be exploited, other design freedoms should not be overlooked. One of these is the qualitative mechanical behavior of the swellable base material, i.e., the matrix in which the fibers are embedded.
Here it is to be noted that previous studies often focus on matrix behavior that has an essential neo-Hookean type character and which reduces to neo-Hookean behavior in the absence of swelling. While such a ground state material is nonlinear in general, its shear response is linear. As we describe here, it is this shear response that governs the matrix contribution to the twist-swelling interaction. Because gels, elastomers, and even biological tissue can exhibit complex ground state behavior in shear -behavior that may depart significantly from a linear response -it is our purpose here to examine the effect of alternative ground state behaviors on the twist-swelling interaction. The range of behaviors considered includes materials that harden in shear, materials that soften in shear, materials that have an ultimate shear stress bound, and materials that collapse in shear.
To this end, we present the overall continuum mechanics formulation in Section 2 where the focus is on the relevant constitutive treatment for the mechanical behavior due to swelling. The overall twisting boundary value problem of interest is formulated in Section 3. It is in fact a swelling version of what is known as the boundary value problem for azimuthal shear, although for our purposes we often refer to it as twisting of a tube due to swelling. Complexities associated with this boundary value problem, while surmountable numerically as we show in later sections, then lead us to examine a simpler problem, namely swelling-induced shearing of a layer. The main advantage is that the layer problem is one of homogeneous deformation. This simpler problem has many useful analogies to the swelling-induced tube twisting problem. This analogy, and certain basic results, is described in Section 4. Up to this point in the paper, the matrix material constitutive relation has concentrated on the neo-Hookean type response. A far more general constitutive model for the matrix, which is a swellable version of the Knowles I 1 power law material [6], is introduced in Section 5. The layer problem for homogeneous deformation under swelling is examined thoroughly for this class of materials. We construct curves showing how the layer shears as a function of swelling. This provides certain important insights into how best to approach the tube twisting problem for the case of the power-law matrix material, and a thorough study of this problem is conducted in Section 6. Now curves of layer twisting as a function of swelling are obtained. Both the layer curves and the twisting curves have certain remarkable features. These are discussed and explained in the concluding Section 7.

Hyperelastic framework for treating swellable matrices and nonswelling fibers
Our examination is based on swellable hyperelastic materials. We consider materials with a single family of fibers at each material point, where the direction can vary from point to point (as in the spiral fiber patterns considered in [3]).
The swellable hyperelastic materials are treated in terms of the swelling volume constraint det F = v where F is the gradient of the deformation mapping from reference locations X to deformed locations x and v is the local swelling value. If v is equal to one everywhere, then this describes an incompressible material. Locations where v > 1 are swollen locations. For a variety of physical processes in soft materials, swelling is determined by chemical or electro-chemical influences that largely decouple from mechanical considerations. In such a case, v can be regarded as a specified field. Moreover, because det F = v is a pointwise constraint, the Cauchy stress will contain a constraint reaction pressure −pI where I is the identity tensor and p is the scalar hydrostatic pressure. The value of p is determined directly from the equations of equilibrium in conjunction with the boundary conditions. In particular, as is true for constraint reactions in general, p is not subject to a constitutive prescription.
Once transient effects associated with possible swelling agent diffusion have abated, it is presumed that the equilibrium mechanical response is elastically dominated. The elastic energy density associated with this response is W = W (C, v, X) where C = F T F is the right Cauchy-Green deformation tensor. The Cauchy stress is then given by In conjunction with the equilibrium equation div T = 0, (2.2) this generates partial differential equations that govern the deformation and the pressure field p = p(X). These are to be solved subject to appropriate boundary conditions. Because we are considering a single direction of fibers at each material point, even though that direction may vary within the material due to fiber patterning, the dependence of W on C is only through the invariants I 1 , I 2 , I 3 , I 4 , and I 5 . The local swelling condition eliminates I 3 because I 3 = v 2 . The fiber direction itself is specified by a unit vector M in the unswollen reference configuration. The deformed image of this vector is m = FM. The vector m, which is generally not a unit vector, gives the fiber direction in the deformed configuration, and its length gives the local fiber stretch.
In this paper, we restrict attention to materials in which the dependence of W upon I 1 , . . . I 5 takes the specific form where w m and w f are energy densities associated with matrix and fiber parts respectively. This models fibers that swell very little within a matrix that is able to swell appreciably. In Eq. 2.3, the invariants I 1 and I 4 are given by I 1 = I : C and I 4 = M · CM = m · m. This makes where w m := ∂w m /∂I 1 and w f := ∂w f /∂I 4 , and B is the left Cauchy Green deformation tensor B = FF T . An alternative situation, fibers that swell in a matrix that does not appreciably swell, can be described by a converse to Eq. 2.3, W = w m (I 1 ) + w f (I 4 , v). For example, this might describe cellulose fibers within a rubbery matrix. However, here our focus is always on Eq. 2.3.
Standard material models in the form of trusted mathematical expressions for w m and w f can then be used for mathematical simulation and design modeling. Among these are the forms that we used in our previous work [3] where we considered the well-known neo-Hookean form, suitably modified so as to treat swelling, for the matrix constituent, namely Here μ > 0 and q are material constants. Specifically μ is the matrix shear modulus for the nonswollen base material. The parameter q allows for possible modulus variation as a function of swelling. We henceforth disregard that possible effect and so take q = 0 in all equations and discussion going forward.
In [3], we also focussed upon where γ ≥ 0 is the fiber modulus. The form (2.6) is often referred to as the standard fiber reinforcing model. Whatever form for w f is chosen, because √ I 4 is the fiber stretch, it is reasonable to require that w f (1) = 0 and w f (I 4 ) ≥ 0. These restrictions cause the reference length of the fiber to also be the natural length of the fiber (the length at which it bears no load). On the other hand, there is no similar requirement on the first derivative of w m as is evidenced by the fact that the neo-Hookean form used in Eq. 2.5 has a constant nonzero first derivative. When Eqs. 2.5 and 2.6 are used together, the dimensionless ratio γ /μ is the essential measure of the relative stiffness of the fiber constituent with respect to the matrix constituent.
In this paper, we retain the form used in Eq. 2.6 for w f but consider alternative forms for w m . Attention is restricted to the plane problem. This means that all the displacement and all of the swelling are in an (X, Y )-plane. Such a situation would arise for example if the domain is a uniform cylinder with a general cross-section Ω in the (X, Y )-plane. In the axial direction, say Z, the cylinder is then confined to be between frictionless plates, say at Z = ± 1 2 L. The fiber direction M is similarly presumed to always lie in the (X, Y )-plane in a manner that is independent of Z. The material and swelling can then vary in the (X, Y ) cross-section while being uniform in Z. This readily realizable state of affairs applies to a potential wide variety of situations and in particular describes the azimuthal shear swelling deformation examined in [3].

Swelling and azimuthal shear
We recall the azimuthal shear boundary value problem for tube swelling considered in [3]. The reference domain is an annulus which is described in polar coordinates Ω = {(R, Θ) : A ≤ R ≤ B, 0 ≤ Θ < 2π }. Here A and B are positive numbers representing the inner and outer radii of the cylinder in the reference configuration. The fiber orientation is also described using polar coordinates M = cos α e R + sin α e Θ , (3.1) where azimuthal symmetry permits α = α(R). The special case where α is independent of R describes a fiber pattern in the form of a geometric spiral [5]. To illustrate, Fig. 1 shows six equally spaced fiber paths in the unswollen reference configuration when α is independent of R. Swelling deforms the geometry. Allowing v = v(R) preserves the azimuthal symmetry. The inner surface of the cylinder is presumed to be fixed as would occur, for example, if it were bonded to a rigid cylinder. The outer surface R = B is presumed to be able to freely expand and so traction free boundary conditions are considered on this surface. The just described set up motivates the consideration of the azimuthal shear deformation: where (r, θ ) are polar coordinate parameters associated with the deformed configuration. The functions r(R) and g(R) need to be determined. No deformation takes place in the third direction, which is formally described by the plane deformation condition z = Z. The boundary conditions for the in-plane deformation on the inner surface are that Boundary conditions due to the traction free outer surface will be clarified later. Clearly, swelling (when v > 1) away from the inner surface will make r(R) > R for locations R > A. The fully 3-D deformation gradient tensor F for the deformation (3.2) is F = r e r ⊗ e R + rg e θ ⊗ e R + r R e θ ⊗ e Θ + e z ⊗ e Z (3.4) where the prime in r and g denotes differentiation with respect to R. Here e r , e θ , e z are unit basis vectors in cylindrical polar coordinates of the deformed configuration. Because det F = r r/R, the local swelling condition det F = v gives r r/R = v which can be integrated to yield Thus, r(R) is fully determined by the swelling in a manner that is independent of the fiber reinforcing. At this juncture, it is convenient to introduce λ = λ(R) and k = k(R) via Here λ is the azimuthal stretch, which is known because of Eq. 3.5. In contrast, k is still to be determined. Observe also that the r appearing in Eq. 3.4 can be eliminated by writing r = v/λ. Direct calculation now gives C, B, m, I 1 , and I 4 in terms of λ, k, and v. In particular, The Cauchy stress tensor (2.4) thus takes the form T = T rr e r ⊗ e r + T θθ e θ ⊗ e θ + T zz e z ⊗ e z + T rθ (e r ⊗ e θ + e θ ⊗ e r ). This is consistent with the frictionless confining plate boundary condition on the fixed Z-surfaces. The in-plane stress components are The equilibrium (2.2) in the z-direction is automatically satisfied by taking p = p(R). In the r and θ directions, the equilibrium equations give ∂T rr ∂r

14)
∂T rθ ∂r These are to be solved subject to the free surface condition on R = B, which is simply The first of these serves as a boundary condition for Eq. 3.14 which in turn serves to determine the pressure field p(R).
Because we do not here choose to examine the normal stresses in detail, we need not consider those equations in what follows. Our interest thus focuses on Eq. 3.15, which integrates to T rθ = c 1 r 2 where c 1 is the integration constant.
In view of the second of Eq. 3.16, it then follows that c 1 = 0 and hence that for all R. Consequently, the Eqs. 3.17 and 3.13 together give where I 1 is given by Eq. 3.9 and I 4 is given by Eq. 3.10. In particular, both I 1 and I 4 are quadratic in the unknown k = k(R). On the presumption that Eq. 3.18 is solvable at each location R then, because k(R) = rg (R), one can integrate to obtain the azimuthal shear field g(R) on A ≤ R ≤ B. This in turn gives the overall twist at the outer periphery as Note that if there is no swelling at all (i.e., if v = 1 for all R) then Eq. 3.5 gives r = R making λ = 1 and hence I 4 = k 2 cos 2 α + 2k cos α sin α + 1. In this case, k = 0 makes I 4 = 1 which causes w f (I 4 ) = 0. Thus, k = 0 solves Eq. 3.5 if v = 1 for all R. In other words, if there is no swelling, then there is no shearing. Equation 3.18 allows for v = v(R) and α = α(R). It is also useful to note that k = 0 solves Eq. 3.18 if either α = 0 or α = π/2. Thus, locations R where the fibers are oriented purely radially or purely azimuthally are not locally sheared.
However, if 0 < α < π/2, then the possible nonlinear nature of w m as a function I 1 , and w f as a function of I 4 , can cause Eq. 3.18 to become a highly nonlinear equation for k = k(R). For the case where w m and w f are given by Eqs. 2.5 and 2.6, certain analytical simplifications result, and this was exploited in the methodology of [3].

Local analogue: combined biaxial stretch and simple shear with swelling (a simple layer)
The deformation (3.2) described in the previous section is locally one of plane biaxial stretch accompanied by simple shear. We can write such a deformation in a local (X, Y ) coordinate system as Here λ x , λ y and K are each fixed constants. The local (X, Y )-system identifies X with Θ and Y with R. Because M as given by Eq. 3.1 makes angle α with respect to the R direction, it follows that the locally defined version of M takes the form The deformation gradient tensor F for the deformation (4.1) is given by Because det F = λ x λ y , the local swelling condition det F = v gives λ y = v/λ x and we use this to eliminate λ y in what follows. Using this, we may calculate C, B, and m, the latter two of which become and (4.6) Similar correspondences hold for I 1 , I 4 , and the stress quantities. This motivates the consideration of plane-strain shearing of a semi-infinite layer that is also subject to swelling. The reference cross-section in the 2-D plane of interest is now Ω = {(X, Y ) : −∞ < X < ∞, 0 ≤ Y ≤ H } and the fibers are oriented as in Eq. 4.2. As in the azimuthal shear problem, the third Z-direction is constrained, thus rendering the problem planar in the same way. Symmetry with respect to X requires that both v and α must be independent of X.
Boundary conditions must be specified on Y = 0 and Y = H . To put Y = 0 in correspondence with the inner circle R = A in the azimuthal shear problem, we stipulate x = X on Y = 0. For the deformation (4.1), this causes λ x = 1. The upper surface Y = H is the analogue of the outer circle. However, here we temporarily allow not only expansion in this direction, but also the possibility of a shear traction τ . It is to be noted that dependence of either v or α upon Y is consistent with this scenario. This is analogous to the fact that the azimuthal shear development of the previous section allowed for the possibility of either v or α to depend upon R. For the present layer problem, if either v or α depend upon Y , then λ y = v/λ x = v would generally vary as well. All of these conditions would then generally make it necessary to let K = K(Y ) in order to satisfy the equations of equilibrium.
For our purposes here, it is sufficient to restrict attention to v and α that are uniform throughout the layer. Then the homogeneous deformation where K is constant satisfies the equilibrium equations provided that p is also constant. Both the amount of swelling v and the applied shear traction τ enter into the determination of K. The deformation itself is thus specialized from Eq. 4.1 to The unit vector M in the fiber direction is deformed into The value of K is found by requiring that the shear stress T xy (which is a function of K and v even though it is constant within layer) must match the loading value τ . This condition is and We have thus arrived at the local analogue problem that is summarized in Fig. 2 as a plane strain swelling deformation of a layer. Both the swelling v and the shear traction τ drive the deformation and so determine K. If there is no shearing (i.e., if K = 0), then swelling (v > 1) expands the layer and elongates the fibers. Shearing (K = 0) can counteract this elongation by rotating the fibers. The amount of shear K is determined by the combined action of τ (which acts to shear the matrix) and the fiber rotational tendency (to relieve the fiber elongation). The solution to Eq. 4.9 determines K on the basis of the trade-offs associated with those two effects.
If τ = 0 then v is the sole driving agent for the amount of shear K. Increasing v in the absence of shearing then elongates the fiber. The direction of shearing that reduces this elongation is one that rotates the fibers toward the vertical. Then K > 0 gives clockwise rotation and K < 0 gives counterclockwise rotation. In order to work with K > 0 (clockwise rotation), we shall take −π/2 < α < 0, and this is the situation shown in Fig. 2. For sufficiently small swelling, this rotation can completely eliminate the fiber elongation, however that is energetically costly with respect to the matrix deformation. The optimum value of K as determined by Eq. 4.9 is one that minimizes the combined energetic cost of matrix shearing and fiber elongation.
For the layer deformation (4.7), the effective shearinduced rotation of the fibers causes them to make an angle of α * with respect to the vertical (y-direction), where tan α * = m x m y = sin α + K cos α v cos α , (4.12) and we note that the combination K = 0 and v = 1 gives α * = α. In other words, as one would expect, absence of both shearing and swelling preserves the fiber orientation. Also, even without shearing (K = 0), a large swelling can cause α * to formally approach zero because of the v in the denominator of Eq. 4.12.
We now focus on the mechanics of how swelling (v > 1) induces a nonzero amount of shear K without any applied shear traction (τ = 0).
As a first example consider the case where w m and w f are given by Eqs. 2.5 and 2.6. These stipulations, along with τ = 0, cause the Eq. 4.9 for K as a function of v to become μK + 2γ K cos 2 α + sin α cos α (I 4 − 1) = 0 (4.13) or, in view of Eq. 4.11, with A 3 = 2γ cos 4 α, A 2 = 6γ cos 3 α sin α, If v = 1 this reduces to μK +γ 2K cos 2 α(K cos α + sin α)(K cos α + 2 sin α) = 0 (4.15) and the solution is simply K = 0. More generally for v > 1, this cubic also has one real root K. Beginning with K = 0 when v = 1 the response curve of K vs. v has initial slope where the sign is based upon taking −π/2 ≤ α ≤ 0. The curves are then monotonically increasing. For very large v, the O(v 2 ) terms in Eq. 4.14 become dominant and one Fig. 3 Response curves of K vs. v as given by the roots of Eq. 4.14. Increase in swelling v generates an increased shearing K although there is an upper bounding limit associated with the fully vertical deformed fiber orientation. Panels display different scales for the horizontal axis. Here α = −π/6 and curves are shown for different fiber-to-matrix stiffness ratios γ /μ. The dashed line shows the asymptote as v → ∞ which is given by K = − tan α = 0.577, a value that is common for all of the curves. The approach to this asymptote is more rapid as γ /μ increases obtains K → −tan α as v → ∞. Additional large v analysis of Eq. 4.14 further yields K = −tan α+ μ sec 4 α tan α 2γ as v → ∞. As the swelling becomes large, not only do the fibers approach a vertical orientation, but they do so in a manner such that m x = m · e x → 0. Figure 3 displays a number of such response curves, all for α = −π/6, showing the effect of the stiffness ratio γ /μ. Each individual response curve approaches the dashed line asymptote K = − tan α as v → ∞. This approach is more rapid for the graphs with larger γ /μ. As γ /μ → ∞ the associated limiting response is the inextensible fiber limit. In this limit, the response curve is found as the root of the equation I 4 = 1. This causes the cubic equation to degenerate into a quadratic, and K actually attains the value − tan α when v = sec α (not as v → ∞). In this inextensible fiber limit, the fibers have become aligned with the y-direction when v = sec α. For the inextensible fiber case, there are no solutions if v > sec α because then the two conditions of detF = v and I 4 = 1 are incompatible.
The qualitative nature of the curves shown in Fig. 3 is similar to the response curves shown in Fig. 7 of [3] for the azimuthal shear problem (the twist-swelling problem). This is the anticipated result since a motivation for considering the layer problem is as a simple analogue to the twistswelling problem. Here it is to be noted that the present X-stretch condition λ x = 1 is the analogue of a condition that would take λ = r/R = 1 in the twist-swelling problem. For the twist-swelling boundary value problem, the condition λ = 1 holds on the inner radius R = A, however λ > 1 for A < R ≤ B. Thus, the layer problem is most representative of the twist-swelling problem near the inner surface.

Power law matrix material and the associated swelling response for the layer
We now examine the effect of the base matrix mechanical behavior on the swelling-shear interaction. In this section, we do this for the layer problem that was just presented in the previous Section 4. Then in the next section, we examine the full azimuthal shear boundary value problem based on an analysis of Eq. 3.18. For these purposes, we continue to consider w f in the form Eq. 2.6. For w m we take where μ > 0, b > 0, n > 0 and q are material constants. This is a swellable version of the Knowles power law incompressible hyperelastic model presented in [6]. In particular, μ is the infinitesimal shear modulus for the unswollen material. For n = 1 and v = 1, this model reduces to the conventional neo-Hookean material. For n = 1 and general v, it reduces to Eq. 2.5. As was the case with Eq. 2.5, we henceforth take q = 0 in all discussion going forward. The incompressible (v = 1) version of Eq. 5.1 allows for the consideration of a wide variety of different material qualitative behavior by simply varying the power law parameter n. To see this, consider the layer with no fibers and no swelling. Then the relation between shear stress τ and the amount of shear K is given by Eq. 4.9 upon taking w f = 0 (no fibers) and v = 1 (no swelling) with w m given by Eq. 5.1. This results in thus retrieving the original result as given in Eq. 4.3 of [6]. The response behavior given by Eq. 5.2 is indicated in Fig. 4. For n = 1, it gives the well-known linear shear stress response associated with the neo-Hookean form. For n > 1, the result is hardening (the curve steepens). For n < 1, it is softening. For 1/2 < n < 1, this softening still results in τ → ∞ as K → ∞. For n = 1/2, the shear stress becomes bounded by a finite ultimate stress. For 0 < n < 1/2, the now more exaggerated softening response leads to a collapse in shear wherein the stress response rises to a maximum and then subsequently tends to zero. where I 1 and I 4 continue to be given respectively by Eqs. 4.10 and 4.11. For n = 1, the Eq. 5.3 reduces to Eq. 4.13. As before, v = 1 implies K = 0 and we again seek to determine the response curves for K as a function of v. The general case for n = 1 leads to a highly nonlinear equation for K as a function of v. There is however one value of n = 1 that permits similar analytical progress and that is the case n = 2. For n = 2, the Eq. 5.3 is again cubic in K and so of the general form Eq. 4.14, but now with different coefficients that are given by bμ + 2γ cos 4 α, A 2 = 6γ cos 3 α sin α, The initial slope of the response curve is again given by Eq. 4.16 but now making use of the above values for A 1 and A 0 . The cubic equation also permits a large v analysis for the n = 2 response and one finds that Thus, as v → ∞, the n = 2 response curve approaches the asymptote K = K (∞,n=2) . With one limiting exception, this asymptotic value is strictly less than the value − tan α that was obtained for the n = 1 neo-Hookean case (viz. Eq. 4.17). The limiting exception is the inextensible fiber limit γ /μ → ∞ in which case K (∞,n=2) → − tan α.
Using Eqs. 5.4 in 4.8 gives m ∼ sin α − 8γ cos 4 α sin α bμ + 8γ cos 4 α e x + v cos α e y , as v → ∞. (5.5) As the swelling becomes large, the fibers again approach a vertical orientation, but now m · e x = m x no longer tends to zero. Figure 5 shows the response curves for both n = 1, n = 2, as well as for several other values of n that do not permit the direct analytical treatment. As in Fig. 3, we take α = −π/6. The horizontal asymptote for n = 2 is evident. Both the n = 2 curve and the n = 3 curve have an internal maximum.

Azimuthal shear response with the power law matrix material
We now return to the full azimuthal shear problem, where now we use the power law material model (5.1) for the matrix response. The reference fiber angle α and the swelling v are both taken to be uniform (independent of R). The latter in conjunction with Eq. 3.5 gives r(R) in the form (6.1) The governing equation for g(R) follows by substituting from Eqs. 5.1 and 2.6 into Eq. 3.18 yielding where I 1 and I 4 are given by Eqs. 3.9 and 3.10 and so also bring in additional terms involving k = rg (R). Thus, Eq. 6.2 is a first-order ODE for g(R) subject to the condition g(A) = 0. When there is neither swelling nor fibers, the corresponding boundary value problem for a tube composed of such a material and subject to an externally imposed twist was studied in [7]. The additional complications due to swelling and the presence of fibers motivated the local analysis of the previous section. It is thus noted that Eq. 6.2 as an ODE for g (R) mirrors the layer Eq. 5.3-which was a scalar equation for the amount of shear K. In particular, the special cases n = 1 and n = 2 give a cubic structure to Eq. 6.2 which facilitates the integration in those two cases.

Linear response in shear (the neo-Hookean type response, n = 1)
If n = 1, then Eq. 5.1 reduces to the standard neo-Hookean material in the absence of swelling. In particular, material parameter b completely disappears from the constitutive form. When swelling is present, we exactly recover the problem studied previously in [3] and the reader is referred to that work for complete details. For our present purposes, representative results are shown in Fig. 6 for a sequence of deformations corresponding to increasing v. If there was no shearing, then each location would only displace radially, giving the hypothetical fiber locations as shown by the dashed lines. However, the localized shearing, which causes the fibers to elongate less, results in the actual fiber locations as given by the solid lines. Thus there is an overall twist.
For a given B/A and γ /μ, the outer surface net twist φ o = |g(B)| is dependent on the fiber winding angle α.  Table 1 gives values of the maxima. The dashed line shows the upper bound as given by Eq. 6.3 There is no twist for either radially oriented fibers (α = 0) or circumferentially oriented fibers (α = π/2). Figure 7 shows how twist φ o = |g(B)| varies with α on the interval 0 < α < π/2. For each v, the value of φ o is maximized for some intermediate value of α. This maximizing value of α varies with v as shown in the accompanying table (Table 1).
The theoretical upper bound value for the twist φ o occurs if the deformed fibers become perfectly radial in orientation. This upper bound value for the outer surface twist g(B) associated with deformed fibers that are straight and radial (dfsr) is (in radians) This bounding curve is also shown in Fig. 7. As v increases, the various curves in Fig. 7 approach this upper bound value, albeit in a highly nonuniform fashion in view of "pinning condition" φ o ≡ 0 at α = π/2. In other words, circumferential fibers (α = π/2) give no shearing tendency upon swelling. However, highly inclined fibers (α = π/2 − with > 0 very small but nonzero) become formally unbounded in length as → 0. Straightening these long fibers thus gives a formally unbounded twist φ o in the limit as → 0. Table 1 The location of maximum outer twist φ o vs. fiber winding angle α (both in degrees) of the five twist-swelling response curves in Fig. 7

Stiffening in shear (n = 2)
In the layer analysis of Section 5, the case n = 2 was the other instance of the Knowles power law matrix material that gave analytical simplification. Similarly, the case n = 2 gives certain analytical simplification to the governing Eq. 6.2 because it reduces to where B 3 = bμ 4 + γ cos 4 α, B 2 = γ (2 + λ) cos 3 α sin α, Other values of n (other than 1 and 2) when used in Eq. 6.2 give highly nonlinear equations for k (with embedded radical expressions) that require numerical root finding procedures (as considered shortly). However, this n = 2 case permits an explicit solution for k and hence g . Numerical integration starting from g(A) = 0 then provides g(R). Figures 8, 9 and 10 show a sequence of swelling-induced twist deformations corresponding to increasing v. Each figure is analogous to Fig. 6, the difference being that now n = 2 instead of n = 1. For the n = 1 case, the parameter b dropped out. The parameter b does not drop out for the case n = 2, and this accounts for the three separate figures showing the effect of b. Taking n = 2 and b = 1 (Fig. 8) to compare against n = 1 (Fig. 6) shows relatively similar  Table 2 gives values of the maxima. The dashed line shows the upper bound as given by Eq. 6.3 twist behavior, with the n = 2 results lagging somewhat. Increasing b for n = 2 decreases the twist (Fig. 9) while decreasing b increases the twist (Fig. 10). Curves of overall twist as a function of fiber winding angle at fixed v, analogous to the curves in Fig. 7 for n = 1, can also be constructed for the case n = 2. One such set is shown in Fig. 11. These curves are the counterparts to the curves in Fig. 7, but now for n = 2 and b = 1. In particular, we observe how the location of the maxima does not shift as much with v as it did for the n = 1 curves. The theoretical upper bound value for φ o as v increases continues to be given by Eq. 6.3. However, now, in contrast to the n = 1 case, we find that the curves remain well below this bounding value as v increases ( Table 2).
The effect of fiber-to-matrix stiffness is shown in Fig. 12. Each curve of twist φ o as a function of swelling v is for a different value of γ /μ, all other parameters being common to all of the curves. As v increases, each curve approaches a separate asymptote, whose ultimate twist value increases as the relative fiber stiffness increases. For example, the ultimate twist value associated with γ /μ = 1 is given by

Other values of n in Eq. 5.1
When n is different than 1 or 2, the governing Eq. 6.2 is not readily solved for k by algebraic means. Thus a numerical root finding procedure is utilized prior to the numerical integration to obtain g(R). In particular, as g(A) = 0, a recursive tangent line approximations can be used to approximate g(R) for A ≤ R ≤ B. For this purpose, the interval [A, B] is divided into some number of equal subintervals where the size of each interval is dR. Figures 13 and 14 show the swelling-induced twist deformation for the softening materials with n = 0.75 and n = 0.5, the latter of which is the case with the ultimate shear stress limit. The unswollen (starting) v = 1 configuration is therefore the same as that shown in previous Figs. 6, 8, 9, and 10. A comparison between all of the different presented cases for n is found in Fig. 15 which clearly shows how increasing n decreases the amount of twist. A verification of the numerical solution procedure that generated all of these figures is presented in Fig. 16. That figure shows the convergence of the numerical solution to that obtained by the analytical (algebraic) root finding procedure for the algebraically solvable case n = 2.
By using the above solution procedure, it is then possible to obtain graphs of twist at the outer surface φ o versus swelling v for the full range of n values. In this regard, the graphed information is like that presented in Fig. 12 but now for different exponents n rather than different stiffness ratios. Figure 17 is representative of such graphs, here using the values n = 0.25, 0.5, 0.75, 1, 2, and 3. The graphs  for n ≤ 1 all approach a common asymptote, and this asymptote is given by φ df sr o from Eq. 6.3. In contrast, the cases n = 2 and n = 3 approach asymptotes that are strictly below that theoretical upper bounding value.

Discussion and conclusions
Whereas previous studies have concentrated on the effect of fiber stiffness and patterning on eliciting complex deformation as a material swells, the present work focuses specific attention on the effect of the base matrix response, i.e., the mechanical behavior of the material in which the fibers are embedded. We show that the base matrix material response also has a significant effect on both simple layer shearing and circular tube twisting as the material swells.
For circular tube twisting, we consider the relative rotation of the tube cross-section when the tube wall contains spirally patterned fibers that proceed from the inner to the outer radius. Because the local deformation is one of simple shear superposed on biaxial deformation, the tube twist findings of Section 6 are anticipated on the basis of certain key understandings of the matrix-fiberswelling mechanics obtained in Sections 4 and 5 for the more straight forward problem of a simple material layer that swells. For the simpler layer problem, the fibers proceed through the layer cross-section along straight lines. Both horizontally aligned fibers and vertically aligned fibers provide no shearing tendency as the layer swells. However, inclined fibers will cause shearing. In all cases throughout this study, the fibers were modeled on the basis of w f as given by Eq. 2.6 Section 4 restricted attention to the matrix material model w m given by Eq. 2.5 with q = 0 and provided the most intuitive key understanding, namely that the generation of shear under simple swelling is due to the presence of fibers, and the effect is magnified as the fibers become relatively stiffer with respect to the matrix. In this case, the ratio γ /μ provided a complete characterization of this relative stiffness. This led to the shear-swelling behavior shown in Fig. 3 with curves for different γ /μ approaching a common asymptote as v becomes large. The rate of approach increased with γ /μ and the common asymptotic value has magnitude tan α. This asymptote is associated with a shearing that places the fibers in the minimal length position, i.e., along the Y axis. Small values of α mean the originally unswollen state involves nearly vertical fibers, and so there is minimal total shearing capacity. Values of α near ±π/2 mean the originally unswollen state involves nearly horizontal fibers and so there is a large shearing capacity and in fact the asymptote's magnitude | tan α| becomes unboundedly large because the fibers become unboundedly long as they proceed from Y = 0 to Y = H in the limit as α → ±π/2. However, if α = ±π/2 then the fibers are purely horizontal, and there is no shearing tendency because the fibers no longer traverse the layer thickness.
All of these results are mirrored in the tube problem for swelling-induced twist where φ o is the analogue of K and the curves of φ o vs. v are the analogue of the curves of K vs. v for the layer problem. The maximum twisting capacity is now given by φ df sr o from Eq. 6.3, which not only involves dependence on α but also on the tube aspect ratio B/A. This is reflective of the fact that the deformed fiber orientation under swelling varies with radius, and every location must have its fibers brought into radial alignment before the twisting capacity is formally exhausted. As in the case of the layer, if w f and w m are given by Eqs. 2.6 and 2.5 then the asymptote φ df sr o is independent of γ /μ, although it is approached more rapidly as γ /μ becomes large. The asymptotic value is formally zero (no twisting capacity) if the fibers in the unswollen state are already radial (α = 0), whereas the asymptotic value tends to infinity as α → ±π/2 because now the unswollen state involves unboundedly large spirals of fiber paths that are nearly circular. However, if α = ±π/2 then the unswollen spirals degenerate to closed perfect circles, and there is no longer any twist under pure swelling, just as the purely horizontal fibers caused the layer problem to lose its shearing tendency.
The above summary results are for w m given by Eq. 2.5 which means that the matrix has a straightforward linear relation in simple shear between shear stress and the amount of shearing. Polymers, gels, and the ground substance in biological tissue may exhibit far more complex behavior in simple shear, specifically some hardening or softening is not uncommon. Furthermore, under softening, the shear traction may or may not have an ultimate value. Finally, softening may be so severe as to elicit a collapse in shear phenomenon as might occur in phase transition or damage scenarioseither permanent or temporary (reversible/healable).
To investigate these effects, w m was replaced with Eq. 5.1 where now the power-law exponent n allows a simple distinction between all of these cases. In addition, the special value n = 1 recovers the original model (2.5). The results of Section 5 show these significant effects in the context of the layer analysis. Here it is immediate that matrix hardening (n > 1), at least as it impacts upon the relative fiber-to-matrix stiffness, would correspond to fiber softening. The converse also holds. Thus, increasing n could be expected to have qualitatively similar effects as decreasing γ /μ, and vice versa. This is indeed the case; however, the layer analysis shows that the effect of varying n is in fact more profound. Here the analytical result (5.4) for n = 2 showed that the asymptote itself, not just the rate of approach to the asymptote, was dependent on γ /μ. Moreover, the full shearing capacity -meaning a situation where the asymptotic value of the K vs. v curves for the layer has magnitude| tan α| -is only achieved as γ /μ becomes unboundedly large. For finite γ /μ, the asymptote height remains below that associated with exhausting this capacity.
This effect is mirrored in the swelling tube findings of Section 6, where Fig. 12 clearly shows how the asymptotes of the various φ o vs. v curves are below the upper bound provided by φ df sr o . For the n = 2 tube twist problem, because it was the solution to a boundary value problem, we were unable to determine this asymptotic value by pure analysis and so found it numerically. This is in contrast to the layer problem, where the relative simplicity of the homogeneous deformation allowed us to extract the asymptotic value K (∞,n=2) < − tan α for finite values of γ /μ.
For values of n different from 1 and 2, we appealed to more extended computational procedures for both the layer and the tube problem. Because n = 1 already gave asymptotes that fully exhausted either the shearing capacity or the twisting capacity, a similar full use of this deformation capacity was found for any n ≤ 1 because by softening the matrix it relatively stiffens the fibers. This can only increase the asymptotes beyond their n = 1 heights. However, since the n = 1 asymptotes already involved full use of the shearing/twisting capacity, such full usage is maintained for n < 1. Thus all the softening cases made use of the full capacity. For this reason, the n ≤ 1 curves in both Figs. 5 and 17 show common upper bound (full usage capacity) asymptotes for the n ≤ 1 cases. Finally, the n = 2 and n = 3 curves in Figs. 5 and 17 also show similar commonality. Now however the n = 2 and n = 3 asymptotes are below that associated with full usage of either layer shearing capacity (Fig. 5) or tube twisting capacity (Fig. 17) as the swelling becomes large.