Static and Dynamic Behavior of Circular Cylindrical Shell Made of Hyperelastic Arterial Material

Static and dynamic responses of a circular cylindrical shell made of hyperelastic arterial material are investigated. The material is modeled as a combination of Neo-Hookean and Fung hyperelastic materials. Two pressure loads are implemented: distributed radial force and deformation-dependent pressure. The static responses of the shell under these two different loads differ essentially at moderate strains, while the behavior is similar for small loads. The main difference is in the axial displacements that are much larger under distributed radial forces. Free and forced vibrations around pre-loaded configurations are analyzed. In both cases the nonlinearity of the single-mode (driven mode) response of the pre-loaded shell is quite weak but a resonant regime with co-existing driven and companion modes is found with more complicated nonlinear dynamics.


Introduction
The study of pressurized circular cylindrical shells has long history. Firstly, the problem was formulated for the case of infinitely small deformations and linear elasticity by Lamé [21], whose name the problem carries now. Literature reviews of early [24] and recent [2] works on circular cylindrical shells made of conventional materials under pressure are available. Here we focus only on the problem with material nonlinearity.
In the middle of the XX th century, some progress was made in the problem incorporating nonlinear materials [25]. In [25], axially symmetric deformations are assumed and the formula describing the internal radius of long cylinders made of semi-linear materials and subjected to an external pressure is obtained. The closed-form stress distribution within a cylinder of incompressible Mooney-Rivlin hyperelastic material is also provided.
In the recent years some experimental studies on hyperelastic circular cylindrical shells have been published [23,14,15]. Most of the modern computational works on the problem [16,1,12,30,31,32,20] employ finite elements. This allows taking into account pressure as the follower load, which is crucial in the case of soft materials, since they usually experience large deformations.
However, the capabilities of commercial finite elements software are rather limited when dealing with the nonlinear dynamics of structures. So works on dynamic behavior of hyperelastic shells usually do not follow this approach. Most early and some new studies [18,19,29,11,28] use very strong assumption of known (and simple) shape of the shell after the deformation.
The meshless approach that the authors proposed in the previous studies [10,6,9] is intended to overcome these difficulties. In the present study we extend the method for the case of a thick shell made of special type of hyperelastic material (combination of Neo-Hookean and Fung materials) that is able to reproduce the actual behavior of an arterial tissue [17]. The shell is subjected to combination of static and dynamic internal pressure loads. We also follow the works [5] and study two types of pressure loads -dead distributed radial force and deformation-dependent pressure, which are similar for small deflections but differ significantly for large strains.

Geometrical relations and description of motion
We consider the thick circular cylindrical shell depicted in Fig. 1 and made of hyperelastic arterial material. The displacements u; v; w of points L h w u v R Â Figure 1: Investigated shell and corresponding coordinate system located on the middle surface of the shell are defined in the axial x, angular Â , and radial z directions of the cylindrical coordinate system, respectively; x and Â are the rotations of the transverse normals at the middle surface about the Â and x axes, respectively. R is the radius of the middle surface, h is the thickness and L is the length of the shell in axial direction.
In the present study, a single-layer shell with radius-to-thickness ratio equal to 6 is investigated, which belongs the thick shell category for which a shear deformation theory should be used in order to obtain accurate results. The strain-displacement relations for the nonlinear higher-order shear deformation theory, as developed by Amabili and Reddy [7] are: where where " x , " Â , " xÂ , " xz , and " Âz are the components of the Green (Lagrangian) strain tensor for shells. We employ Lagrange equations to describe the behavior of the shell: d dt L P q n Á L q n D Q n ; n D 1; : : : ; N; (3) where L D T … is the Lagrange's functional, T , the kinetic energy of the shell, …, the potential deformation energy; Q n are the generalized forces, and N is the number of the generalized coordinates q n : Q n D F=q n , where F is virtual work done by external forces. A static pressure and a dynamic radial point force will be considered. The dot stands for differentiation with respect to time. The potential and kinetic energies are given by [7]: and … D • V W dV (5) where W is the strain energy density (SED), V , the volume of the shell, and , the mass-density of the shell material. The displacements and rotations are expanded into truncated series:

Hyperelastic material
The nonlinear elasticity of soft biomaterials is usually described by hyperelastic laws. In most cases, such materials are assumed to be incompressible [17,26]; the incompressibility hypothesis is used in our study. In the present study we use the combination of Neo-Hookean term to describe the isotropic part of the tissue response and the Fung term, which describes the anisotropic response of the fibers [17]: where E is the Young's modulus of the shell material; I 1 is the first invariant of the right Cauchy-Green deformation tensor C D 2E C I, and E is Green-Lagrange strain tensor: The reason for using the combination of a Neo-Hookean term and the Fung term to describe the material is related to the mechanical property of biological soft tissues, like aortic layers. In fact, these types of tissues present an almost elastic behavior for small to medium strains, which can be accurately described by a Neo-Hookean hyperelastic model. For large strains instead the material stiffness is exponentially increasing, due to the fact that fibers (e.g. elastin and collagen fibers in aortic layers) that are originally loose in the soft tissue become to stretch; those fibers are much stiffer than the surrounding soft tissue and give a sudden increase of the global stiffness that can be accurately described by an exponential law.
The first invariant is expressed as We also need the third invariant The material is assumed to be incompressible, so the components of (8) are not independent. The incompressibility condition J D 1 yields Expression (11) is further substituted in (8) and (7b).

Distributed radial force versus area-dependent actual pressure
Two pressure-like loads are compared. In both cases, we assume that the load is applied on the middle surface of the shell. The first one is a distributed radial force that does not depend on the deformations of the shell, i.e. this is a dead load. It always acts in the radial direction z and the total load acting on the shell is area-independent. The expression for virtual work done by this force is [4] F D Z L 0 Z 2 0 P wR dÂ dx; (12) where P is the magnitude of the distributed radial forces per surface area.
The second type of load is the actual pressure, which is dependent on the area of the deformed surface and the resulting force is always orthogonal to it. The expression for this type of load applied to circular cylindrical shell is [5] F D

Local Expansion of the Strain Energy Density
Equations (2) together with " z in (11) are not polynomial in strains, which essentially complicates the investigation of the shell behavior. The analysis is simplified by introducing a transformation of strain energy density (SED) (7b) in order to derive approximate governing equations in the form of ordinary differential equations with polynomial nonlinearities of order not higher than three. First, we describe the method for the case of static loading. A static configuration (for example, the undeformed configuration) is assumed to be known through the generalized coordinates q .0/ . A new configuration defined in the vicinity of the known one q D q .0/ C˛q .1/ ; (14) is to be calculated, where˛is a formal small parameter and the incremental generalized coordinate vector q .1/ is unknown. Expression (14) is substituted into (2) and the resulting strain-displacement relationships are subsequently inserted into (7a). As a result, the SED is obtained as a function of q .1/ and˛. Then, we expand the SED into a truncated Maclaurin series 0/ / C : : : (15) where W .i/ is the i -th power polynomial of the components of the vector q .1/ . The formal only small parameter˛is introduced to show that the incremental generalized coordinate vector q .1/ is small with respect to q .0/ . Expansion (15) is truncated after the term W (5), so we have the SED as a 4-th power polynomial in q .1/ , and hence, Lagrange equations (3) are equations with quadratic and cubic nonlinearities only. There is a variety of available numerical techniques to solve the nonlinear algebraic equations obtained from Lagrange equations for static problems [8]; Newton's method is preferred in the present study. However, this model is able to describe the behavior of the shell only locally, i.e. in the vicinity of a given configuration q .0/ around which the SED is expanded. Accordingly, we name this approach the local models method (LMM).
Once the algebraic equations are numerically solved, the new configuration q given by Eq. (14) is used as the initial guess q .0/ for the next iteration. By successive iterations of local models, the final deformed configuration is found. Iterations continue until the required deflection or applied force is reached. Afterwards, the nonlinear dynamics is analyzed around the previously found static configuration, denoted q .0/ , using Lagrange equations (3) with the SED in Eq. (15). It must be observed that in case of dynamics, nonlinear coupled ordinary differential equations are obtained and their numerical solution requires dedicated treatment. More details on the LMM can be found in [10]. The material parameters are similar to those of the human aorta adventitia layer as listed by Holzapfel [17]: E D 51 900 Pa, C D 471 10 6 Pa, c 11 D 37:7, c 12 D 58, and c 22 D 63:8. Note that the material properties in the longitudinal and circumferential directions here are permuted with respect to [17] and to the adventitia layer in the human aorta. The shell is simply supported on its edges and the boundary conditions are w D 0, v D 0, Â D 0, N x D 0, and M x D 0 where N x is the axial stress resultant per unit length and M x , the axial stress moment resultant per unit length, i.e.
As a simplifying assumption, in (16) we take x as the simple linear Cauchy stress [4] x D and Poisson's ratio for a incompressible material is D 1=2.
Expressions (18) are substituted into (2) and then in Eqs. (4), (5), and (7a). This yields the expressions for kinetic and potential energies as polynomial functions in the generalized coordinates and the derivation of the Lagrange equations (3) is then straightforward; they are not explicitly provided for conciseness purposes. In the numerical example, note that generalized coordinates with three-subscripts fw i;n;c .t /, w i;n;s .t /, u i;n;c .t/, u i;n;s .t/, v i;n;c .t/, v i;n;s .t/, x i;n;c .t/, x i;n;s .t/, Â i;n;c .t /, Â i;n;s .t /g are used instead of the single-subscript coordinates fq i .t/g in order to specify their physical meaning. The generalized forces Q n are obtained by differentiation with respect to generalize coordinates (i) of expression (12) in case of dead load (radial distributed forces) or (ii) of expression (13) in case of follower load (actual pressure). The constant part of Q n in the latter case coincides with the expression for Q n in the former case, but in the latter case Q n includes also linear and quadratic terms in the generalized coordinates that are not present for dead load.

Static analysis
Since both loads are axially symmetric, only 15 DOF corresponding to the axisymmetical modes (n D 0) in (18) are considered here.
The force versus central deflection curves for both radial force and actual pressure obtained with the LMM (see Section 5) is shown in Fig. 2. For small deflections, the curves almost superimpose, but the actual pressure curve becomes almost vertical for large deflections because the stiffness grows very quickly. Also, the inflection point (i.e. the point at which the static behavior suddenly changes) appears at smaller central deflection.  The three-dimensional distorted shells corresponding to points A and B in Fig. 2 are shown in Fig. 3. In case of radial force, the axial displacements are very large and of the order of magnitude of the radial displacements; they induce a significant length contraction. In case of actual pressure, axial displacements are much smaller. This can be explained by the fact that the actual pressure induced force always acts orthogonal to the deformed shell surface, giving a large axial load component which balances the shell contraction by almost the same amount in the present investigation.  Fig. 2) 6.2. Dynamic analysis Since real arteries are never perfect, perturbations are introduced in the linear stiffness matrix to break the symmetry of the system. As a result, the eigenfrequencies corresponding to cos.nÂ / and sin.nÂ / in Eq. (18) are now distinct thus neutralizing the pitchfork bifurcation giving rise to driven and companion mode participation in the nonlinear frequency response [4]. However, the two natural frequencies begin extremely close, a perturbed traveling wave in circumferential direction still emerges from perturbed driven and companion modal participations. The dynamic response in the frequency range centered on the first non-axisymmetric mode eigenfrequency is now explored. Since the shell is long, this first non-axisymmetric mode corresponds to n D 2 in Eq. (18) [22]. First, we investigate the free and forced nonlinear vibrations of the shell statically pre-loaded around the deformed configuration B in the Fig. 2. The comparison with the static solution shows that the local model is accurate for deflections up to 0:2h which cannot be considered small for thick shells.
The first eigenfrequency in this case is 189:27 rad=s and it corresponds to the axial vibrations with predominant coordinate u 1;0 , and is therefore not further investigated since vibration modes with predominant radial displacement are of interest. The second and third eigenfrequencies are 339:03 rad=s and 339:06 rad=s respectively and corresponds to the w 1;n;c and w 1;n;s predominant modes; those are the modes that we are interested in. The harmonic balance method [27] is employed to find the periodic solutions of system (3). Every generalized coordinate is expanded as a truncated series of cos.j! /, j D 0; : : : ; N h , where N h indicates the harmonic balance expansion order. A convergence analysis shows that for N h D 3 leads to an acceptable approximate solution.
The free response backbone curves of the shell pre-loaded by radial distributed forces (point B) are shown in Fig. 4. The forced vibrations are studied with the computer program AUTO [13]  dedicated to the numerical integration of coupled nonlinear ordinary differential equations. This software combines pseudo-arclength continuation and collocation method to continue the solution starting from an initial solution point. The trivial configuration with zero external dynamic load is used as a starting point. The solution is first traced with respect to the excitation amplitude when the excitation frequency is kept fixed. Once the desired excitation level is reached, the solution is continued with respect to the excitation frequency and the full frequency range of interested is spanned at a fixed excitation level. The external dynamic force P t cos.!t/ is taken in the form of a periodic radial point load, applied at x D L=2, having an harmonic component P t D 0:76 N and frequency ! in the vicinity of the natural frequency ! n of the investigated mode; the frequency ratio is defined as D !=! n . Modal damping ratio n D D 0:005 is inserted in Eq. (3). The forced vibration response for the bending coordinate w 1;n;c is shown in Fig. 4 together with the free response backbone curve. The backbone curves at moderate amplitudes are very close to the linear ones. The forced vibration response exhibits a loop.
To compare the amplitudes of bending and in-plane displacements the frequency responses for coordinates u 1;n;c and v 1;n;c are exposed in Fig. 5. Their responses are similar to the response of  1;n;c with magnitudes ten and two times smaller, respectively.
Next, we study the free and forced vibrations of the shell pre-loaded by actual pressure around the deformed configuration marked with point A in Fig. 2. The deformed shape of the shell is shown in Fig. 3[left]. The local model is accurate for deflections up to 0:2h in this case too. The lowest eigenfrequency of 159:32 rad=s corresponds to vibrations with predominant mode u 1;0 and is not investigated further. The first bending eigenfrequencies are 264:99 rad=s and 265:1 rad=s for participations w 1;n;s and w 1;n;c , respectively. Again, the harmonic balance method with N h D 3 is used to study free vibrations. The backbone curve is shown in Fig. 6 together with the forced vibrations obtained with AUTO. A periodic radial point load, applied at x D L=2, with P t D 0:375 N, As in the previous case, the responses of u 1;n;c and v 1;n;c are very similar to the response of the main bending mode w 1;n;c but their amplitudes are approximately 8 and 1.65 times smaller than the amplitude of w 1;n;c .
In order to see both driven and companion mode active, a larger dynamic excitation is necessary, as shown in Fig. 8 for P t D 0:76 N. In this case, the vibration amplitude might exceed the validity limits of expansion (15).

Conclusions
In this paper, the behavior of a circular cylindrical shell under static and dynamic pressures is studied. The shell is made of an arterial biomaterial described through a hyperelastic law that captures the principal feature of the soft tissue nonlinearity: a steep increase in stiffness after certain strain threshold.
Two types of pressure-like loads are analyzed: the dead load in the form of distributed radial forces and actual pressure which is an area-dependent follower load always normal to the deformed surface. At small deflections, the difference between these two loads is insignificant but increases with larger deflections. The shell under actual pressure exhibits much steeper increase in stiffness due to axial stretching coupled to circumferential stretching. The shell under radial forces shrinks in length essentially, which is not happening for the case of actual pressure. The reason of this is the rotation of the normal in case of actual pressure, so that the load acts also in axial direction.
The dynamics of the pre-loaded shell is studied under the hypothesis of very small imperfections to make the study more realistic. So, no bifurcations are observed due to the broken symmetry of the system. In both cases, the nonlinearity of the single-mode (driven mode) response of the pre-loaded shell is quite weak, but a resonant regime with both driven and companion modes active has been found with more complicated nonlinear dynamics.
The approach presented in this study targets the analysis of blood vessel deformations and thus can be used in case of composite thick shells too. The approach is also applicable to the cases of laminated composite shells, different boundary conditions and various types of loading, including thermal loading. In order to improve the model in the case of human aorta, the material properties and boundary conditions should be assessed more carefully; moreover, viscoelastisity and fluid-structure interaction have to be taken into account. The present model is intended for shells, so that the length to radius radio should not be excessive in order to avoid that the fundamental mode of the system becomes a beam mode (n D 1, where n is the number of circumferential waves of the fundamental mode).