Wrinkles in the opening angle method

We investigate the stability of the deformation modeled by the opening angle method, often used to give a measure of residual stresses in arteries and other biological soft tubular structures. Specifically, we study the influence of stiffness contrast, dimensions and inner pressure on the onset of wrinkles when an open sector of a soft tube, coated with a stiffer film, is bent into a full cylinder. The tube and its coating are made of isotropic, incompressible, hyperelastic materials. We provide a full analytical exposition of the governing equations and the associated boundary value problem for the large deformation and for the superimposed small-amplitude wrinkles. For illustration, we solve them numerically with a robust algorithm in the case of Mooney-Rivlin materials. We confront the results to experimental data that we collected for soft silicone sectors. We study the influence of axial stretch and inner pressure on the stability of closed-up coated tubes with material parameters comparable with those of soft biological tubes such as arteries and veins, although we do not account for anisotropy. We find that the large deformation described in the opening angle method does not always exist, as it can become unstable for certain combinations of dimensions and material parameters.


Introduction
One of the most effective ways to demonstrate the existence of residual stresses in biological structures is to isolate a cylindrical shape and cut it axially. Invariably it will open up, revealing that the cylinder was under a large circumferential stress, see (c) Figure 1: Cutting biological cylindrical structures radially reveals that they were under circumferential residual stresses. Left: slice of an Irish Ash tree; Middle: a green chilli pepper; Right: Equatorial slice of rat heart (taken from [10]).
In turn, one of the most successful advances of non-linear elasticity is the modeling of this stress through the so-called opening angle method. By measuring how much a tube opens up into a sector, one can reconstitute a backward scenario whereby the structure was initially an open circular sector, subsequently bent into a complete tube by the action of what can now be identified as a residual stress, see Fig.2. Hence the opening angle gives a measure of the level of residual stress for an assumed model of material behavior.
Of course many questions remain open at the end of the process and here we address the following: Is the bending deformation always possible, or is it limited by loss of stability with respect to small-amplitude static wrinkles? Moreover, can the instability be overcome by pressurization of the reconstituted tube? These issues are most relevant to Finite Element simulations of residually-stressed tubes, where buckling should be avoided as much as possible.
Here we first formulate in Section 2 the equations governing the large deformation of a coated circular sector into an intact tube, which is possibly subjected to an internal hydrostatic pressure and a uniform axial stretch. We then specialize the analysis to the case when the coating and the substrate are made of different Mooney-Rivlin materials, because the stress components can then be computed analytically. We pay particular attention to writing the boundary conditions properly (hydrostatic pressure on inner face, perfect contact at the interface, traction-free on outer face).
In Section 3 we present the algorithm implemented to solve the incremental problem of static wrinkles superimposed onto large bending, axial stretch, and pressurizing. It relies on the Stroh formulation and the Surface Impedance Matrix method, and is robust and unaffected by numerical stiffness.
Finally, Section 4 presents experimental and numerical results: first our own, achieved by gluing a silicone coating on a urethane substrate; and second those coming from the literature on soft biological tubes, although of course those cannot be accurately modeled as isotropic. In our experiments, we show that no wrinkles form when a sector of opening angle 120 • is closed, while wrinkles form before a sector of opening angle 240 • is closed. Applying the aforementioned algorithm, we show numerically that the critical opening angle at which wrinkles form is 209 • and that four wrinkles should appear along the circumference, which is consistent with the experimental results. Applying the algorithm for dimensions and material parameters comparable (with the limitation that anisotropy is not accounted for) to those of a rabbit artery, we show that, in the absence of internal pressure, wrinkles form for an opening angle of 320 • , but that these wrinkles can be eliminated by applying an internal pressure or can be delayed by the presence of an axial stretch. These results are in line with intuition and experiments made on biological tubes.

The opening angle method
Consider the sector of a soft cylindrical tube with geometry delimited in the cylindrical coordinate system {R, Θ, Z} (and orthonormal basis {E R , E Θ , E Z }) in its natural state B 0 by the region where A, C are the radii of the inner and outer faces of the sector, respectively, L is its height, and α 0 ∈ (0, 2π) is the opening angle. The stressfree circular sector consists of a stiff thin layer placed at the inner side (A ≤ R ≤ B), glued onto a thicker and softer layer located in the outer region B ≤ R ≤ C, where B is the radius of the interface between two layers, as shown on Fig.2a. From now on, the superscripts (c) and (s) refer to the coating and the substrate, respectively.
E↵ect of coating on the wrinkling in the problem a cylindrical sector deforming into an intact tube

Sigaeva Taisiya
May 15, 2017 1 Problem formulation Let us consider a bi-layered sector of cylindrical tube made of an incompressible hyperelastic material with respect to a cylindrical coordinate system {R, ⇥, Z} in its natural state B 0 . Sector's geometry is defined by the following region: where A and C are the radii of the inner and outer faces, respectively, of the sector, H is length, ↵ 0 2 (0, 2⇡) is the opening angle. The circular sector is formed by two layers: a sti↵ relatively thin layer serving as coating that is placed at the compressive side of the system (A  R  B) and a second larger and softer layer in the region B  R  C, where B is the radius of the interface between two layers, as shown on Fig.1a. Should we define n sectors using R i as in 1 Figure 2: The opening angle method: (a) An initially stress-free coated sector is subject to axial stretch and bent into (b) a residually-stressed full tube. It can also be subject to (c) an internal pressure. But is that large deformation stable?
The sector is deformed into an intact (circular cylindrical) tube with respect to a cylindrical coordinate system {r, θ, z} (with orthonormal basis {e r , e θ , e z }) by the following mapping [4] where is a measure of the opening angle and λ z ≥ 1 is the uniform axial stretch. We denote this configuration by B r and refer to it as the residually-stressed configuration. The geometry of the tube is now where a = r(A), b = r(B), c = r(C) and = λ z L is the current tube length as shown on Fig.2.
The associated deformation gradient F is The incompressibility condition, det F = 1, and one of the geometric requirements, e.g. r(A) = a, impose Taking into account the diagonal form of deformation gradient, we introduce the principal stretches such that λ 1 λ 2 λ 3 = 1 to satisfy incompressibility. We take both coating and substrate to be made of isotropic hyperelastic materials with strain energy densities W (c) , W (s) , respectively, so that the Cauchy stress σ is diagonal in the e i ⊗ e j basis, with components Here l = c, s and q (l) are the Lagrange multipliers arising from the incompressibility condition.
In the absence of body forces the only non-trivial equation of equilibrium is ∂σ rr ∂r For the boundary conditions, we assume that the inner (coated) face of the tube at r = a is under internal pressure P , that there is perfect bonding between the two layers at the interface r = b, and that the outer face at r = c is free of traction: By introducing the following quantities [4], we may rewrite the principal stretches in terms of x as λ 1 = 1/ √ kλ z x, λ 2 = kx/λ z so that the energy density for fixed λ 3 = λ z may be seen as a function of x only: integrating equilibrium equations (9) for each layer, and using boundary conditions (10), we find that the inflating pressure P is We can also determine the stress components throughout the wall, as For a given geometry of an undeformed coated sector in B 0 , the following quantities are prescribed, Then the physics of the stretched and pressurized closed-up cylinder in B r are prescribed by the given strain energy densities W (l) for coating and substrate, the given axial stretch λ z and the given inner pressure P . The new geometry is entirely determined by solving the system of three equations for the three unknowns x a , x b , x c composed by Eq.(13) and the two relations Then the state of stress is entirely determined by Eqs. (14). For illustration, in this paper we model the substrate and coating using the Mooney-Rivlin energy density; it reads where C (l) 2 > 0 are material constants and C = F T F is the right Cauchy-Green deformation tensor. This model is quite general because it recovers, at the same level of approximation [3], the most general model of isotropic, incompressible, third-order weakly non-linear elasticity, where E = 2C + I is the Green-Lagrange strain tensor, µ is the Lamé coefficient of linear elasticity, and A is the Landau coefficient of third-order elasticity (The connections between the constants are µ = C 1 + C 2 , A = −4C 1 − 8C 2 .) For the Mooney-Rivlin material (17), we have which provides explicit expressions for the stress components in Eq.(14). Hence and so on for the other components.
For an example, assume that the coating is Γ times stiffer than the substrate, in the sense that C where Γ ≥ 1 is the stiffness contrast factor. Then we consider how the stresses are distributed along the radial axis for different stiffness factors Γ. We take the case where there is no inner pressure (P = 0) and the opening angle is 139 • . In the undeformed geometry we take A = 13mm, B = 14.5mm, C = 18mm. Fig.3 illustrates the distribution of stresses along the thickness of the wall of closedup cylinders, for a uniform material (Γ = 1.0), and for two-layered solids with moderately (Γ = 3.0) and significantly (Γ = 7.0) stiffer coatings compared to substrates. We clearly observe the jump in the circumferential stresses at the interface between coating and substrate, as expected.

Wrinkling of a coated sector
Here we study the stability of a coated sector closed into a pressurized cylinder. We signal the onset of instability by the existence of small-amplitude wrinkles, solutions to the incremental equations of equilibrium. From experimental observations, we know that they should be varying sinusoidally along the circumference of the tube, with amplitude decay from the inner face to the outer face. The analysis for the existence of such wrinkles can be put together from the results of the previous section and those of Destrade et al. [4] and we omit the details to save space.
In short, the wrinkles exist when the following boundary value problem is solved for z (l) = z (l) (x), (l = s, c), the 2 × 2 Hermitian surface impedance matrix [2]. (ii) Numerical integration of the differential Riccati matrix equation (iv) Numerical integration of the differential Riccati matrix equation (21) in the coating (l = c), from x b to x a ; and (v) Target condition: In Eq.(21), † denotes the Hermitian transpose and the Stroh sub-matrices G i have components [4], where the superscript "(l)" is understood, n denotes the wrinkling mode (number of wrinkles in the circumference), and Here, in general, and in particular for the Mooney-Rivlin model, Finally, the derivation of the target condition (22) is detailed in the appendix.

Experimental & numerical results
Here we implement the stability analysis described in the previous section for two cases: polymers and biological tissues. The algorithm is illustrated in Fig.7(a). Essentially, we implement the steps (i)-(iv) and iterate over α 0 until the target condition (v) is reached. We denote by α cr = α 0 the critical opening angle at which wrinkles form when the sector is closed into an intact tube, i.e., the value of α 0 when the target condition is reached.

Results for polymers
For our first experiment, we used artificial materials, namely relatively stiff silicone (red) for the coating, urethane (black) and very soft silicone (white) for the substrate. We subjected each material to a tensile test using a MTS electromechanical material characterization machine. We then determined the Mooney-Rivlin constants by curve-fitting over a useable range of data, and found that C  Fig.4(a). We then glued a 1.6mm thick red silicone layer onto a 26.9mm thick black urethane sector (B = 23.93mm, C = 50.83mm) and produced two coated sectors, one with opening angle 120 • , the other with opening angle 240 • , see Fig.4(b) and (c). We produced similar sectors using white urethane as the substrate, see Fig.4(d).
We found that, for the black urethane substrate, no wrinkles formed when the former sector was closed ( Fig.4(b)), while for the latter sector six wrinkles formed shortly before the sector became intact (Fig.4(c)). Thus we would expect the critical opening angle at which wrinkles form when the sector becomes intact to be somewhere between 120 • and 240 • . To check this assertion, we performed the stability analysis described in the previous section for the same dimensions and material parameters as in the experiments. We found that the critical opening angle was 209 • with corresponding mode number n = 4, which supports our previous hypothesis.

Results for soft tissues
Here we perform the stability analysis using the dimensions and material parameters which are of the same order of magnitude as those of a rabbit carotid artery, as collected by Holzapfel et al. [8].
The artery consists of three layers: the intima, the media and the adventitia. However, the intima is very thin and not very stiff (at least in healthy young individuals), and so we can use our two-layer model with the dimensions [8] B − A = 0.26mm, C − B = 0.12mm, A = 1.43mm, along with an axial stretch λ z = 1.695.
For the material parameters, Holzapfel et al. [8] used an anisotropic model. Here we have only considered isotropic models, and so we set to zero Holzapfel et al.'s anisotropic parameters to make a (somewhat arbitrary) connection with their measurements. Moreover, Holzapfel et al. [8] did not consider a dependence of W on the second invariant of strain tr(C −1 ), so here we take C of the shear modulus for the artery's elastin matrix. We perform the stability analysis over a physiological pressure range [5] of 0-170 mmHg. We plot the results in Fig.7(c) for the non-dimensional measure of pressureP = P/C (s) 1 .
Then the physiological pressure range corresponds to 0 ≤P ≤ 75.5. First we plot the curves giving the critical opening angle α cr against the pressureP for increasing values of the mode number n = 2, 3, 4, . . .. Each curve is a bifurcation plot: at a given pressureP , a tube with opening angle larger than α cr will buckle when it is bent into an intact closed tube; in order not to buckle, a sector must have an opening angle which is less than the smallest critical angle from all curves. Here we find that all curves for mode numbers n ≥ 5 are all below those for n = 2, 3, 4 and are virtually indistinguishable one from another, see Fig.5. Hence our analysis does not allow us to determine the mode number precisely here, in contrast to the scenario of Section 4.1. Figure 5: (a) Plots of the critical opening angle α cr for several mode numbers n versus the non-dimensional pressureP using the material parameters and dimensions comparable with those of a rabbit artery [8], when it is subject to an axial stretch λ z = 1.695. (b) Plots of the critical circumferential stretch λ a on the inner face of the intact tube at buckling versus the pressureP .
From the plots we see that when there is no internal pressure (P = 0), only sectors with an opening angle greater than α cr 320 • will buckle when closed into an intact tube. This value is significantly above the recorded opening angle for the rabbit artery [8], which was 160 • . Hence we would expect (provided the crudeness of our modelling arteries here is overlooked) that the rabbit artery is smooth when it is not subject to internal pressure.
We also observe that as the internal pressure increases, the critical opening angle increases, with asymptotic behaviour α cr → 360 • asP → ∞. Hence buckling can be eliminated by applying an internal pressure, which is in line with our intuition and with, for example, experiments on a rat's pulmonary artery [6], see Fig.7(b).
For comparison, we also plot the curves obtained in the case of no axial stretch, λ z = 1, see Fig.6. We find that the axial stretch makes the sector more stable with respect to bending into an intact tube (the values of α cr are higher when λ z > 1 than when λ z = 1). To complete the picture, we also provide the plots of the variations of the critical circumferential stretch λ a (contraction stretch on the inner face of the intact tube). Figure 6: (a) Plots of the critical opening angle α cr for several mode numbers n versus the non-dimensional pressureP using the material parameters and dimensions comparable with those of a rabbit artery [8], when it is not subject to an axial stretch (λ z = 1.0). (b) Plots of the critical circumferential stretch λ a on the inner face of the intact tube at buckling versus the pressureP .

Discussion
Often it is assumed that a stable deformation of a sector into an intact tube exists. These "opening angle" deformations are then used to estimate the residual stresses in the material [7]. Here we have shown that, depending on the material properties and dimensions, wrinkling may occur before the sector becomes intact, which would be followed by further buckling and creases when the sector is closed. Our results have important implications for finite element reconstructions of the opening angle method. First, a stiffer coating will lead to instabilities in finite element simulations, earlier than for a homogeneous sector [4,7]. Second, if the wrinkles occur, then our analysis is a first step towards providing meaningful precursors to creases (see Fig.1 and Fig.7(b)).
We also showed that wrinkles can be eliminated by applying an internal pressure, as has been confirmed in experiments.
Our method could also be applied to other tissues such as the esophagus, which is often modeled as a two-layered structured, and in which wrinkles and creases have been observed [12]. However, it is important to consider the limitations of our model. For example, in the iliac artery of an 81 year old human, buckling of the intima in the zero-pressure state leading to delamination has been observed [9]. As has been noted, the intima, one of the three layers of the artery, becomes thicker and stiffer with age. Evidently, there are residuals stresses present leading to buckling, but clearly a three-layer model would be necessary to investigate such an occurrence. Furthermore, each layer of the artery is highly anisotropic due to the presence of collagen fibers [8], and so a more realistic model would reflect this fact.
at r = a, where for the second equality we used the incremental incompressibility equation, div u = ∂u ∂r + 1 r u + ∂v ∂θ = (rU + U + inV ) e inθ r = 0.
On the other hand, the traction is related to the displacement by the surface impedance matrix [2]: In particular, at the r = a interface, we have by (31) from which the target condition (22) follows (see Balbi and Ciarletta [1] for an early, but not entirely correct, derivation of the target condition).