Swelling Induced Twist in Hyperelastic Tubes Due to Spiral Patterned Biasing Fibers in the Cross Section

Simple fiber reinforcing patterns can serve to guide deformations in specialized ways if the material experiences expansion due to some sort of swelling phenomenon. This occurs even when the only activation is via the material swelling itself; the fibers being a passive hyperelastic material embedded in a swellable hyperelastic matrix. Using anisotropic hyperelasticity where the usual incompressibility constraint is generalized to model swelling, we consider such fiber guided deformation in the context of a circular cylinder subject to uniform swelling. The material is taken to be transversely isotropic with a fiber pattern corresponding to helical spirals in each cross section. This paper extends previous work which had examined a traction free outer radius that expanded while the inner radius was held fixed. Because of the spiral pattern, the tube in these previous studies exhibited increasing twist as the swelling proceeded. The problem considered here takes both inner and outer radius as free surfaces, thus causing the amount of radial expansion itself to be unknown. It is found that the spiral fiber pattern again induces a twist, and that this pattern also influences the nature of the radial expansion.


Introduction
Within the context of anisotropic hyperelasticity we consider the effect that special stiff fiber patterning can have on deformation as the material changes volume.This volume change could arise from a variety of sources and is here referred to as swelling.For the purposes of the analysis presented in this paper the swelling volume change is regarded as a prescribed (input) field quantity.Such a modeling approach has been shown to be useful in various actuator applications (e.g., [9,10,19]), as well as in the modeling of soft biological tissue under certain types of swelling [14-16, 27, 28].The mathematical treatment of the latter has features that are similar to the modeling of biological growth even though growth would typically manifest itself over longer time scales [11][12][13].Conversely, in other applications and physical systems it may be that the evolving stress field will have a significant effect upon the swelling mechanism.In such cases a more involved modeling approach would likely be called for (see e.g., [1,8,21,23,25]).
The present work extends [5] and [6] which addressed a situation in which nonsymmetric helical fiber patterning gave rise to an azimuthal shear deformation as the material becomes swollen.The azimuthal shear deformation is one in which cylindrical domains twist as a function of radial distance.This contrasts with torsional deformations in which cylindrical domains twist as a function of axial distance (which can also be induced by swelling [4]).The boundary value problem of azimuthal shear has been extensively studied in incompressible hyperelasticity where applied tractions generate a moment that causes a relative twisting between an inner radial location and an outer radial location [2,3,7,17,24].Subsequently, [5] and [6] show how swelling can generate the same type of deformation even when the external surfaces are free of load induced twisting moments provided that a nonsymmetric helical fiber patterning is present.The stiff fibers, which do not experience the same swelling volume change as the background matrix material, twist the system in order to elongate less than they otherwise would if there were to be no azimuthal shearing.The works [5] and [6] examined this issue for an annular cylinder in which only the outer radius was free to expand because the inner radius was bonded to a rigid core cylinder.The present work considers a geometry in which both the inner and outer radii are free to deform.This makes for a more involved boundary value problem than that considered in [5] and [6].The purpose of the present paper is to clarify this distinction and to provide the additional mathematical treatment so as to analyze this more complicated state of affairs.
The development is organized as follows.Section 2 provides the general framework for the overall swelling treatment, along with the boundary value problem formulation in a manner that allows for the consideration of rather general hyperelastic strain energy densities.As mentioned above, the continuum treatment is in terms of anisotropic hyperelasticty.This means that there are both matrix and fiber contributions to the overall strain energy density for an effectively homogenized continuum.Section 3 then provides the specialization to a strain energy density consisting of an I 1 power law for the swellable matrix contribution and a standard reinforcing function in terms of I 4 for the nonswelling fiber contribution.
The matrix part of the strain energy density allows for the consideration of a neo-Hookean type swelling model as a special subcase, and Sect. 4 develops the corresponding analysis in order to provide specific examples.As swelling proceeds there is both an inner and outer radius expansion, as well as an increasing overall twist.The amount of twist will depend upon the fiber winding angle α.There will be no twist at all for the special symmetric fiber patterning that involves either radially aligned fibers (α = 0) or circumferentially aligned fibers (α = π/2).Among other results, twist vs. swelling response curves are obtained for the intermediate helical fiber pattern with α = π/4.The analysis is extended in Sect. 5 so as to consider matrix material response that is not neo-Hookean in nature.This includes both highly softening and highly stiffening matrix material response in the large deformation regime.Whereas the neo-Hookean type matrix response led to certain analytical simplifications, the more general response models considered in this section require an enhanced computational treatment.For certain highly stiffening behaviors it is shown that there is an internal maximum in the twist vs. swelling response.Final commentary is provided in Sect.6.

Swelling and Rotational Shearing of a Circular Tube with in-Plane
Fiber Biasing We consider a hollow circular cylinder that, prior to swelling, occupies the region The cylinder is made of a swellable hyperelastic material that contains fine scale reinforcing fibers that are distributed in a helical spiral fashion on planes of constant Z.At each location these fibers are taken to all align in the same direction.In other words there is only one family of fibers.Hence, from the continuum perspective, the material is locally transversely isotropic.The fiber pattern is taken to be uniform in Z and .Its reference configuration direction can therefore be described by a unit vector where, in general, α = α(R).The special case in which α is independent of R gives an initial fiber pattern in terms of logarithmic spirals [17] as shown in Fig. 1.Swelling gives a volume change that is described in terms of a swelling field v that is the ratio of the local swollen material volume to the unswollen volume.Attention is restricted to radially symmetric swelling v = v(R).The swelling condition gives the kinematic constraint where F = ∂x/∂X is the deformation gradient.In the absence of swelling (v = 1) this reduces to the usual constraint associated with the incompressible material theory.The corresponding hyperelastic material theory with swelling corresponds to a material stored energy density function W = W (F, v) and a Cauchy stress tensor given by where I is the identity tensor and p = p(X) is the hydrostatic pressure associated with the constraint (2.2).It still remains to specify boundary conditions on the inner and outer radii R = R i and R = R o .In this work both of these surfaces will be taken to be traction free.The swelling Fig. 1 A sketch of 12 equally spaced fiber paths in the reference configuration for a geometry with R i = 1 and R o = 2.The fiber paths are given by (2.1) with α = π/4 and hence α is independent of R in this case.Consequently fiber patterning is that of a logarithmic spiral [17] field v then becomes the sole driving agent for the deformation.In particular, the tube will both expand (because of the swelling) and twist (because of the fiber biasing).
The just described setup motivates the consideration of axially symmetric azimuthal shear deformations of the form: ( The functions r(R) and g(R) need to be determined.The locations of the deformed inner and outer radius are denoted by r i := r(R i ) and r o := r(R o ).The same type of deformation was considered for the boundary value problems examined in [5] and [6].Those works considered a different boundary condition at R = R i , namely a zero displacement condition associated with a fixed inner radius.Such a condition immediately determines the function r(R).In contrast, the problem under study here involves a tube that is free at both its inner and outer radius.This makes for a somewhat more involved analytical development.
The deformation gradient F for the deformation (2.4) is given by where the prime in r and g denotes differentiation with respect to R. Consequently, detF = r r/R which is now used in the swelling condition (2.2) to yield which can be integrated to obtain where the as yet unknown (2.8) Note that once c 1 is determined for a given v then r(R) is known.In [5] and [6] there was no need to introduce the equivalent of c 1 since the fixed boundary condition at the inner radius would immediately lead to it taking on the known value R 2 i .The plane strain condition z = Z confines all of the swelling expansion to be in (r, θ ) cross sectional planes.It is useful to be aware of the possible special situation for which this in-plane swelling expansion is the same in all such planar directions.This shall be referred to as equidirectional in-plane swelling.This would occur, for example, in an isotropic material cylinder (no fibers) that was first allowed to freely expand in all directions (i.e., F = v 1/3 I), and then was squeezed between frictionless plates in the z-direction so as to restore the original length.This would make F = v 1/2 (e r ⊗ e R + e θ ⊗ e ) + e z ⊗ e Z and T = T zz e z ⊗ e z where T zz is the confining stress generated by this process.In view of (2.5) such an equidirectional expansion corresponds to the additional requirements that g = 0 and r = r/R.The second of these makes r = v 1/2 R and c 1 = vR 2 i .In general c 1 will be found to differ from vR 2  i in the problem under consideration here, meaning that the swelling will not be equidirectional in the plane.
Returning back to the case of general c 1 , it follows from (2.5) that the left Cauchy-Green tensor B = FF T is given by B = (r ) 2 e r ⊗ e r + kr (e r ⊗ e θ + e θ ⊗ e r ) + k 2 + χ −2 e θ ⊗ e θ + e z ⊗ e z , (2.9) where we have introduced χ = χ(R) and k = k(R) in the form The right Cauchy-Green tensor As in a hyperelastic treatment without swelling, the stored energy density W depends upon F only through C. Since the material composing the cylinder is both transversely isotropic and subject to a pointwise volume constraint, the stored energy density W depends on C only through the invariants I 1 , I 2 , I 4 and I 5 .Note that the local swelling condition makes It will henceforth be assumed that W is independent of both I 2 and I 5 .For the deformation (2.4),I 1 is given by (2.13) Note that the unit vector M is mapped into (2.14) Of special significance going forward is the stretch in the fiber direction √ I 4 where I 4 is the invariant (2.15) Because we consider W = Ŵ (I 1 , I 4 , v), the derivative of W with respect to F in (2.3) is accomplished by differentiating Ŵ with respect to the two invariants I 1 , I 4 with the aid of the chain rule in the usual fashion.This yields where the notation is simplified by replacing Ŵ with W . Then by virtue of the forms (2.11) and (2. ) .17) ) The Cauchy stress in the absence of body forces is required to satisfy the equilibrium equation div T = 0. (2.20) Because v = v(R), the particular equation of equilibrium (2.20) associated with the axial direction is satisfied identically provided that p = p(R).The equations of equilibrium in the r and θ directions reduce to ) The second of these, (2.22), integrates to T rθ = c 2 r 2 where c 2 is yet another constant.
The development thus far hews closely to that of the swell-twist boundary value problem studied in [5] and [6].In those papers the inner boundary R = R i was held fixed and only the outer boundary was traction free.The fixed inner boundary in [5] and [6] then made c 1 = R 2 i in (2.7) and the traction free outer boundary then made c 2 = 0 and hence T rθ = 0 for all R. Thus the equivalent of χ is fully known (there was no need to even introduce such a concept in [5] and [6]), simplifying the problem considerably by reducing it to a problem for a single unknown function that could either be taken as k or g in the above formulation (with the connection k = rg ).
Now however, we study the effect of releasing the fixed boundary constraint at R = R i and replacing it with a traction free condition in its own right.Hence the tube is subjected to the following boundary conditions and The zero shear traction part of the boundary conditions, either (2.23) 2 or (2.24) 2 , in conjunction with T rθ = c 2 r 2 gives c 2 = 0. Consequently, for all R, which was also the case in [5] and [6].Consequently the two equations (2.25) and (2.19) combine to give It is to be emphasized that it is the swelling field v = v(R) that completely drives the deformation in this case.Specifically, because the initial (reference) configuration, which is both undeformed and unswollen, is stress free, it follows that χ = 1 and g = 0 in the absence of swelling.This contrasts with the azimuthal shear investigations [2,3,7,17] in which there was no swelling, and the deformation was driven by an applied twisting moment (or equivalently an imposed twist).This made T rθ = 0 in [2,3,7,17].Of these, [3] is the closest in spirit to the boundary value problem that is now under consideration here.In [3] there is no swelling and there is also a zero normal traction condition T rr = 0 on the surfaces r i and r o .Thus the first of (2.23) and (2.24) was stipulated in [3], but not the second of each.From a practical point of view that leads to the question of what kind of mechanism acting at the boundary in [3] could apply such a boundary condition while still allowing unhindered radial deformation.Such an issue does not arise in the present investigation because the traction free boundary condition is the most natural kind that there is.
It is to be remarked that the absence of any displacement boundary condition in the context of the deformation (2.4) will allow an arbitrary rigid body rotation about the origin to be added to any solution of the boundary value problem.We therefore fix the rigid body rotation by taking It is immediate from (2.16) and (2.17) that T rr − T θθ is free of the pressure term p. Hence integration of (2.21) gives an integral expression for T rr that does not contain p.To perform the integration it is convenient to rewrite (2.21) as which follows from using (2.6) and the chain rule.Integrating with respect to R and equating the resulting expression for T rr with the expression for T rr given by (2.16) allows for the extraction of the pressure function p(R): Here c 3 is the associated integration constant.It can be found by expressing the boundary condition T rr (r i ) = 0 in terms of p(R i ) upon using (2.16).Comparing the resulting expression for p(R i ) to the evaluation of (2.28) at R = R i then shows that c 3 = 0.It remains to enforce the boundary condition T rr (r o ) = 0. Now (2.16) evaluated at r = r o in conjunction with (2.28) For our purposes this important governing equation shall be written as where (2.30) By using (2.16) and (2.17) it follows that H (c 1 ) is expressed directly in terms of derivatives of W as (2.31) We are now in a position to summarize the problem as follows: • We seek a function k(R) and a constant c 1 so as to satisfy the conditions (2.26) and (2.29).Here (2.26) is a field equation and (2.29) is a single value equation.Even so it is important to observe that both equations are coupled for the determination of k(R) and c 1 .
Note that H is a function of c 1 and a functional of k(R) even though our notation H (c 1 ) only emphasizes the former.It is also to be remarked that these equations still permit a general dependence α = α(R).Once k(R) and c 1 are known then g(R) is obtained by integrating k = r g using (2.7) and (2.27).
Although (2.26) and (2.29) are coupled equations, they do have distinct physical interpretations.Equation (2.29) can be viewed as ensuring the proper baseline pressure value that permits meeting the zero normal traction condition on both the inner and outer radius.While this causes T rr to vanish on the boundaries, in general T rr will not vanish in the interior.In contrast, (2.26) ensures that the shear stress T rθ vanishes not only on the boundaries but also in the interior.In view of this latter observation it will be useful to express (2.26) as a balance between the constituent shear stresses, namely with individual matrix and fiber contributions to the shear stress, It is worth pointing out that the intuitive expectation is that twist should only occur if fibers are present in a biased way.In particular, no fibers should mean no twist.This is indeed ensured by the formulation.In addition, no fibers should also mean that the resulting radial expansion is simply that of equidirectional in-plane swelling.To verify these claims note that if no fibers are present then ∂W ∂I 4 = 0 making T (f ) rθ = 0. Accordingly (2.32) requires that T (m)  rθ = 0, which is satisfied by taking k = 0, i.e., no twist.Then, putting k = 0 into the integrand for (2.31), which now only contains the term with ∂W ∂I 1 , we observe that the factor which corresponds to the equidirectional in-plane swelling response.

Constitutive Specifications
We now provide more specific constitutive assumptions in order to proceed with a mathematical treatment for some specific physically motivated examples.It will be assumed that swellable matrix-fiber composite has a separable energy density form: where W m and W f are energy densities associated with matrix and fiber components respectively.Because W f does not contain v, the specification (3.1) models non-swelling fibers embedded within a swellable hyperelastic matrix.
The matrix part of the strain energy density will be given by the following swelling generalization of the Knowles power law material Here μ, b and n are positive material constants; in particular μ is the shear modulus.In the absence of swelling (when v = 1) this W m reduces to the hyperelastic form as employed in [18] and [24].The above swellable form was used in our previous study [6].The special case n = 1 gives the swellable neo-Hookean material as used in a variety of modeling scenarios [26,29,30].
For n = 1 this gives the linear-in-k response T (m) rθ = 1 2 μ χ k, in which case b has dropped out.For n = 1 the parameter b plays a role in determining the range of small k for which the response remains close to the n = 1 neo-Hookean type of response.To see this consider a hypothetical scenario in which r(R) is given by the equidirectional in-plane swelling response r = v 1/2 R which would occur in the absence of fibers.In this scenario χ = v −1/2 and I 1 = 1 + 2v + k 2 whereupon (3.3) becomes Figure 2 shows the associated T (m) rθ versus k response using (3.4) for swelling value v = 2.It is apparent from both panels in this figure that the qualitative behavior is highly dependent upon the parameter n.If n > 1 then the departure from linearity is stiffening in k.If 1 2 ≤ n < 1 then the departure from linearity is softening in k.Parameter values 0 < n < 1  2 involves an eventual decreasing response ("collapse in k") which is potentially useful in modeling phenomena such as phase transformation and damage, although we do not focus on such phenomena here.The left panel of Fig. 2 is for b = 1 and shows a clear range of k where all response curves remain close to the n = 1 neo-Hookean type response.As indicated by the right panel, larger values of b give a more immediate departure from the neo-Hookean type response.
It is important to remark that the solution to the boundary value problem when fibers are present will generally have I 1 = 1 + 2v + k 2 and so will not give the equidirectional in-plane swelling response described by (3.4).Thus although (3.3) will hold in what follows, it will generally be the case that (3.4) is, at best, a useful approximation.
Turning to the fiber constitutive response, a simple model for the fiber energy function, W f (I 4 ), is given by where γ > 0 is the fiber modulus.The form (3.5) is sometimes referred as the standard fiber reinforcing model [20].By using (3.2) and (3.5) in (2.26), the field governing equation takes the form where I 1 and I 4 are given respectively by (2.13) and (2.15).The single value governing equation (2.29) takes the form 4 Response for the Case of a Neo-Hookean Type Matrix (n = 1) 2) reduces to a neo-Hookean type swelling response for the matrix material.
In particular, material parameter b completely disappears from the analysis.This relatively simple and informative case is examined in this section.The field governing equation (3.6) for n = 1 simplifies to where I 4 is given by (2.15).Because this I 4 is a quadratic polynomial in k it follows that (4.1) reduces to a cubic equation for k.The coefficients depend upon both R and c 1 because of the presence of χ in both (2.15) and (4.1).Cubic equations of this broad form often arise in the solution of boundary value problems for hyperelastic materials involving the standard reinforcing (3.5) because this causes the shear stress contribution from the fiber constituent to naturally depend upon the local value of shear in a cubic manner [22].It is found that there is always one real root and a complex conjugate pair of roots to this cubic equation.In particular, a sign analysis of the coefficients shows that there is exactly one real positive root k for all possible parameter values.This leads to an analytical expression involving cube roots for k = k(R, c 1 ) > 0. Using this in (2.31) it follows that H (c 1 ) is a simple function of c 1 , i.e., it no longer operates as a functional.Consequently, (2.29) has become a single uncoupled equation for the determination of c 1 .In particular, it is now the case that In the above expression, r, χ , k and I 4 are all dependent upon R (the integration variable) with expressions given in terms of the as yet undetermined constant c 1 .
If v = 1 then the choice c 1 = R 2 i makes χ = 1, k = 0 and I 4 = 1 which in turn gives H (c 1 ) = 0.This confirms that v = 1 gives r(R) = R.Moreover, k = 0 makes g = 0, which, in conjunction with the condition g(R i ) = 0, gives g(R) = 0 for all R. Clearly, this is an expected consistency result: no swelling gives no twist.
Table 1 Values of radial expansion parameter c 1 = r 2 i and amount of twist parameter g(R o ) as swelling proceeds for the n = 1 neo-Hookean type matrix material.Other parameter values are More generally, the governing equation H (c 1 ) = 0 must be treated numerically.It is also to be noted that if α is independent of R then it is easier to integrate (4.2).Numerical examples in what follows are restricted to this case of constant α.
Figure 3 displays four separate graphs of H (c 1 ) when α = π/4.Each graph corresponds to a different uniform swelling value v, namely: v = 0.5 (deswelling case), v = 1 (noswelling case), and v = 1.5 and v = 2 (swelling cases).The zero crossing locations provide the sought after roots c 1 to (2.29) which solve the boundary value problem for the indicated material and geometric parameters.Numerical values for these c 1 roots are given in Table 1.
The example under consideration in Fig. 3 and Table 1 takes R i = 1.This can be viewed as choosing R i as the characteristic length for the purposes of nondimensionalization.The choice R i = 1 is convenient for the interpretation of c 1 because an equidirectional in-plane swelling will then cause c 1 to be equal to v.While the previous considerations have established that R i = 1 will cause c 1 to be equal to v when v = 1, there is no reason to expect that c 1 will be equal to v when swelling takes place (i.e., no reason to expect equi-directional in-plane swelling).This is confirmed by Table 1 where c 1 is consistently larger than v when v = 1.Because c 1 = r 2 i this shows for this example that the radial expansion under swelling gives more radial expansion than that which would occur under equidirectional in-plane swelling.The converse occurs for the deswelling case of v < 1 where there is less radial contraction than that which would occur under equidirectional in-plane deswelling.Turning now to the twist, it follows from (2.10) using (2.8) and (2.27) that where k(R) = k(R, c 1 ) is the unique real root that follows explicitly from (4.1).The total amount of twist proceeding from the inner radius to the outer radius is given by g ).These values are also shown in Table 1 confirming that increasing swelling gives increasing twist.
Figure 4 provides a pictorial view of these deformations, showing both the expanded cross section (contracted for deswelling) as well as the deformed location of a typical fiber.The deformed fiber location shows the additional twisting displacement that occurs compared to how the fiber would deform in the absence of twist.As the matrix swells, the additional twisting displacement causes the fiber to displace in an opposite sense to that of Fig. 5 Plot of H (c 1 ) when v = 2 for six selected values of the fiber stiffness parameter γ when n = 1 (neo-Hookean type swellable matrix).Other parameter values continue as α = π/4, μ = 1, Table 2 Effect of γ on the swelling-twist response when v = 2, again using the n = 1 matrix model.Other parameter values continue as the original spiral helix (in this case the original helical pattern involves counter clockwise spirals as one moves outward, whereupon the swelling induced twist is clockwise).This serves to straighten the fiber, allowing it to elongate less than that which would occur if no twisting took place.The presented results are for both α and v independent of R. As shown in this relatively simple setting, the swelling response is sensitive to fiber orientation α, the cylinder wall thickness ratio R o /R i and the relative stiffness γ /μ of the nonswelling fiber constituent to the swelling matrix constituent.A systematic examination of these effects was conducted in the previous studies [5] and [6] which focused on the relatively simpler boundary value problem where the inner tube radius was attached to a rigid inner core.For the more challenging boundary value problem under consideration here, the relative stiffness effect is particularly striking.It is easy to analyze in this n = 1 case because of the previously described analytical solution k(R, c 1 ) that follows explicitly from (4.1).Taking v = 2, Fig. 5 displays the associated H (c 1 ) graphs for six different values of fiber stiffness parameter γ while keeping μ = 1.The other parameter values are as before and indicated in the figure caption.As one can expect, the curves are ordered with respect to γ , with the curves moving rightward as γ increases.For each stiffness value γ , the value of c 1 giving H (c 1 ) = 0 is shown in Table 2 along with the swelling induced twist g(R o ).As one would expect, there is no twist when there are no fibers (γ = 0), and the swelling induced twist g(R o ) increases as the fibers become relatively stiffer (i.e., as γ increases).Observe also that if γ = 0 then c 1 = 2 = v which makes the azimuthal stretch r/R coincident with the radial stretch r , thereby confirming the anticipated equidirectional in-plane swelling when no fibers are present.Conversely, the presence of fibers (γ > 0) causes the azimuthal stretch r/R to differ from the radial stretch r .

Response for the Case of a Power Law Matrix with n = 1
We now turn to consider a matrix constituent response that differs from that of the swellable neo-Hookean model.This is done by taking n = 1 in (3.2).It was shown in the previous section that the neo-Hookean type case gave certain simplifications in the analysis of (3.6).A similar simplification occurs for n = 2 and so that case is addressed first.We then move on to consider the more involved case where n is neither one nor two.

Shear Stiffening Response (n = 2)
If n = 2 then governing equation (3.6) takes the form While somewhat more involved than the n = 1 equation (4.1) due to an extra term, this extra term still results in an equation that is cubic in k.This provides a second case in which it is feasible to explicitly solve (3.6) for k = k(R, c 1 ).Other values of n (other than 1 and 2) give more involved equations that do not make it feasible to obtain similarly explicit expressions for k as a function of R and c 1 .As before, this removes the functional dependence of H on k making H a direct function of c 1 .The second governing equation H (c 1 ) = 0 now proceeds with with r, χ , k, I 4 and now also I 1 all depend upon R with expressions given in terms of the as yet undetermined constant c 1 .
As was the case for n = 1, the no swelling case v = 1 gives c 1 = 1 and also no twist.For v = 1 this equation is treated numerically in a similar fashion to the n = 1 case.For n = 2 the value of the constitutive parameter b in (3.2) plays a prominent role, where it is recalled from Fig. 2 that larger values of b cause a swifter departure from the linear response associated with n = 1.We consider an example with b = 10.
Figure 6 is the n = 2, b = 10 counterpart to the previous Fig. 3.The corresponding values of the expansion parameter c 1 and the overall twist g(R o ) are given in Table 3. Except for the baseline no-swelling case v = 1 where there is no radius change and no twist, it is again noted that c 1 > v. However now the c 1 values are closer to v than they were for the n = 1 case given in Table 1.Thus the radial expansion is somewhat closer to that of an equidirectional in-plane swelling.Also, we note that if v = 1 then |g(R o )| for n = 2 is less than the corresponding |g(R o )| for n = 1.Thus there is also less twist when n = 2.Both of these observations are consistent with the increased matrix stiffening associated for n = 2, which has the effect of lessening the effect of the fiber constituent.4 shows broad qualitative agreement.As already remarked, the fibers exert less relative influence for n = 2 than for n = 1.This effect increases with b.Conversely, as b → 0 the response of the n = 2 material approaches the n = 1 response (recall that the n = 1 response is independent of b).This effect of b is shown in Fig. 8.Note that the b = 10 curve shows that the power law matrix stiffening can eventually result in the overall twist reaching a maximum at a particular swelling value, followed by a gradual recovery of the previously accumulated twist.

A Modified Shooting Method for General n
As discussed in Sect.5.1, the governing equation (3.6) cannot be solved explicitly for k = k(R, c 1 ) except for the special cases n = 1 or n = 2.These special cases had the beneficial effect of causing the integrand in (3.7) to depend upon c 1 in an explicit, albeit complicated, manner.This allowed for numerical evaluation of the integral in (3.7) by means of standard software packages.For other values of n we find that such an approach is no longer tenable.Accordingly, an alternative more direct computational method is employed for solving the problem numerically for other values of n in the material model (3.2).  and g(R o ) for different values of n when the swelling parameter is v = 2, all taking b = 10.Other parameter values continue as  Divide the interval [R i , R o ] into N equal subintervals and let R be the length of each subinterval.Note that R = (R o − R i )/N .Suppose that all geometry parameters R i , R o , α, swelling parameter v and material parameters μ, γ , b, n are given.Suppose c 1 is assigned a value that will used as a "shooting parameter".Then r(R) is completely determined via (2.8) and hence χ(R) is also completely determined from (2.10).At each subinterval endpoint R j , j = 0, 1, . . ., N, the value k(R j ) can now be obtained from numerical solution of (3.6).This enables numerical evaluation of the integral in (3.7) by means of a Riemann sum.Governing equation (3.7) thereupon generates where G(c 1 , k(R), R) is the corresponding integrand.Search and update on the shooting parameter c 1 can now proceed until the above equation is met to within a desired tolerance.We find that it is important to take R sufficiently small.In our calculations with R i = 1 and R o = 2 we have taken R = 0.0001.This enabled an estimate of c 1 to within four decimal places which in turn led to (5.1) being satisfied to within 0.000001.Then with c 1 in hand, g (R j ) follows from (2.10) using k(R j ) and r(R j ).Numerical integration then gives g(R).
The above procedure has been tested against the more analytic solution procedure employed when n = 1 and n = 2.The two methods are in good agreement.Using this algorithm, we consider matrix softening models: n = 0.75, n = 0.5, n = 0.25, as well as the additional stiffening material n = 3.The most extreme softening model n = 0.25 when used without fibers gives the collapse in shear phenomena previously remarked upon in the context of the discussion around (3.3).We take b = 10 in order to sufficiently depart from the n = 1 response when v = 2. Table 4 then shows the associated v = 2 responses for these different n.The corresponding results are depicted pictorially in Fig. 9.
It can be seen from both Table 4 and Fig. 9 that the net twist decreases as n increases.This is again consistent with a matrix stiffening because that is akin to fibers that are relatively less stiff.This serves to inhibit the overall twist as n increases.As illustrated in Fig. 10, these trends persist for all values of the swelling parameter v.

Concluding Remarks
Consistent with previous studies, we have examined how standard fiber reinforcing (i.e., inactive/passive fibers) can guide material deformation in a specialized way as the overall material expands due to swelling.These specialized deformations serve to lessen the stretch in the fiber direction compared to that which would otherwise take place in the absence of fibers.This gives rise to localized shears whose action over finite distances generates the specialized deformation modes.In this work the specialized deformation was azimuthal shearing in a circular tube.The fiber pattern which enables this deformation during a swelling expansion is that of a helical spiral.Straightening of the fibers then "unspirals" the pattern, giving rise to the twist.The present work examines these issues in the context of a hyperelastic swelling treatment.
Previous related work [5] and [6] in a similarly patterned tube considered a fixed displacement boundary condition on the tube's inner surface while the outer surface was traction free.In that case the swelling expansion was purely outward from the fixed inner radius.This immediately determined the radial part of the displacement, leaving only the twisting portion unknown.In the present work, both the inner and outer tube radii are traction free so that the radial portion of the displacement must also be determined.Because the overall swelling that drives the problem is stipulated, the unknown radial expansion is reduced to the determination of a single constant, denoted by c 1 .The mathematical treatment for the total deformation then reduces to the determination of the constant c 1 along with the determination of the localized shearing field, here denoted by k(R).In contrast, only the equivalent of k(R) needed finding in [5] and [6].
In addition to the fiber patterning, the specific details of the resulting deformation are dependent upon the modeling choice for the matrix material part of the strain energy W m and the fiber portion of the strain energy W f .The specific choices considered here enable certain simplifications to the resulting mathematical problem (two coupled equations -an algebraic equation and a functional equation -for the determination of c 1 and k(R)).A means for solving these equations are described in this work, enabling the determination of the overall twist as a function of swelling for the spiral fiber patterns under consideration.In most cases, including the case of a neo-Hookean type matrix containing reinforcing described by the standard fiber model, the amount of twist was found to increase monotonically with the amount of swelling.However certain other cases corresponding to a matrix that was not of neo-Hookean type displayed alternative interesting behavior.In particular, a matrix that more rapidly stiffened in comparison to the fibers was found to exhibit a local maximum in the twisting response as a function of swelling.
14) for B and m it follows that T = T rr e r ⊗ e r + T rθ (e r ⊗ e θ + e θ ⊗ e r ) + T θθ e θ ⊗ e θ + T zz e z ⊗ e z , with nontrivial components

Fig. 2
Fig. 2 Matrix shear stress response T (m) rθ versus amount of shear k using the power law model (3.4) with μ = 1 and v = 2.The qualitative behavior changes with constitutive exponent n.The linear response range decreases with b as can be seen by comparing b = 1 (left panel) with b = 10 (right panel)

Fig. 4
Fig. 4 Expansion/contraction of the cylinder cross-section and twist induced displacement of a fiber path when n = 1 for selected values of the swelling parameter: v = 0.5 (top left panel), v = 1 (top right panel), v = 1.5 (bottom left panel), v = 2 (bottom right panel).The dashed curve is the hypothetical location associated with no twist (i.e., only the radial displacement).Other parameter values are α = π/4, μ = 1, γ = 1, R i = 1 and R o = 2

Fig. 8
Fig. 8 Overall twist |g(R o )| versus swelling v comparing the n = 1 and n = 2 material models for the matrix constituent.Other parameter values are α = π/4, μ = 1, γ = 1, R i = 1 and R o = 2.The n = 1 response does not depend upon b whereas the n = 2 does depend upon b.As b → 0 the n = 2 response converges to that for n = 1

Fig. 9 4 Fig. 10
Fig. 9 Radial expansion and twist displacement of a fiber path when v = 2 for the power law material (3.2) showing the effect of n while keeping b = 10.The six values of n are: n = 0.25 (top left panel), n = 0.5 (top right panel), n = 0.75 (middle left panel), n = 1 (middle right panel), n = 2 (bottom left panel), n = 3 (bottom right panel).Other parameter values continue as R i = 1, R o = 2, μ = 1, γ = 1, and α = π/4 where (R, , Z) are the polar coordinate parameters associated with this reference configuration.The values R o and R i obeying R o > R i > 0 give the initial outer and inner radii of the cylinder.Let e R , e and e Z denote orthonormal unit vectors in this polar coordinate description of the reference configuration.Polar coordinates in the deformed configuration are denoted by (r, θ, z) with corresponding orthonormal unit vectors e r , e θ , e z .Attention is restricted to plane deformations with z = Z, as might occur for example if the cylinder is constrained in the Z direction by frictionless plates at Z = ±L.