New nonlinear solution of nearly incompressible hyperelastic FGM cylindrical shells with arbitrary variable thickness and non-uniform pressure based on perturbation theory

author https://doi.org/10.1590/1679-78255622 Abstract In this paper, nonlinear analysis of thick cylindrical shells with arbitrary variable thickness made of hyperelastic FGM with radially variation of material properties in nearly incompressible state under non-uniform pressure loading is presented. Thickness and pressure of the shell vary in axial direction by linear and/or nonlinear functions. The governing equilibrium equations are derived based on shear deformation theory (SDT). The Mooney-Rivlin type material is considered which is a suitable hyperelastic model for rubbers. Boundary Layer Method of the perturbation theory which is known as Matched Asymptotic Expansion (MAE) is used for solving the governing equations. A new ingenious solution and formulation have been defined during current study to simplify and abbreviate the representation of inner and outer equations components in MAE. In order to validate the results of the current analytical solution, a numerical modeling based on Finite Element Method (FEM) have been investigated. Afterwards, for different rubber case studies, the effect of material constants, inhomogeneity index, geometry and pressure profiles on displacements, stresses and hydrostatic pressure distributions resulting from MAE and FEM solution have been illustrated. This approach enables insight into the nature of the deformation and stress distribution across the wall of rubber vessels and offers the potential for investigating the mechanical functionality of blood vessels such as arteries in physiological pressure range. The results prove the effectiveness of SDT and MAE combination to derive and solve the governing equations of nonlinear problems such as nearly incompressible hyperelastic of variable thickness and non-uniform pressure have been presented in special representation separately. The effect of materials constants, inhomogeneity index, geometry and pressure parameters on displacements, stresses and hydrostatic pressure distributions resulting from MAE solution have been investigated for some case studies and the results have been compared with a FEM modeling in ANSYS software. We present the equations that provide the general continuum description of the deformation and the hyperelastic stress response of the material. Current study aims to illustrate the performance of the potentials and their reliability for the prediction of the state of deformation and stress in hyperelastic FG vessels from rubbers to arteries. more uniform distribution of displacements and stresses in heterogeneous cylinder. It can be concluded that pressure profiles increment is more effective on the response of shell respect to thickness profiles variation. Furthermore, changes of concave thickness profile to convex one lead in descending maximum displacement, stresses and hydrostatic pressure. It can be concluded that radial displacement and hydrostatic pressure patterns follow the pattern of the applied pressure function along the length of shell. The behavior of hyperelastic FG vessels under non-uniform pressure distribution shows that similar profile of variable thickness and non-uniform applied pressure result in minor displacement and stress quantities and uniform distributions which could be a suitable criterion in designing thickness profile of pressurized vessels. Applying maximum pressure and consequently maximum thickness near boundaries of shell are suitable profiles in designing hyperelastic FG shells. Authors believe that current method along with studies mentioned in the literature could direct further research toward the design, optimization, and manufacture of graded rubber-like materials.


INTRODUCTION
Hyperelastic materials are quite common in many engineering applications. These materials are incompressible or almost incompressible and undergo large strains when subjected to loads. In the last decades, many constitutive models are developed for hyperelastic materials, which can be used in computational model according to the application. The Mooney-Rivlin model of hyperelastic materials can simulate most of the mechanical behaviour of the rubber materials. The model provides a good description of the mechanical properties of rubber materials when deformation is less than 150%. Rubber products are used in different industrial applications; such as rubber hose to carry fluids, rubber antivibration mountings, cylindrical pneumatic floating rubber fenders for boats and so on. Furthermore, rubber seals for sealing connectors are used to very easily seal on the internal or the external diameters of test parts which have smooth cylindrical connections. Rubber cylindrical sleeves have been used for many years successfully for label printing and have been proven of value for the established printing processes. Packer rubber produces a larger contact pressure and forms a seal between the rubber and the casing, resulting in sealing of the annular gap and isolation of the production layer. A comprehensive survey on the finite element methods of incompressible or almost incompressible hyperelastic materials can be found in many papers (Boyce and Arruda 2000). As an important research, Sussman and Bathe (1987) introduce a displacement-pressure finite element formulation for the geometrically and materially nonlinear analysis of compressible and almost incompressible solids. Simo and Taylor (1982) analyzed incompressible nonlinear elastic solids by a penalty function approach. Bijelonja et al. (2005) presented development of a displacement-pressure based finite volume formulation for modelling of large strain problems involving incompressible hyperelastic materials with a Mooney-Rivlin model. Some useful researches focus on developing strain energy function that can describe the mechanical behavior of rubber-like materials and incompressibility characteristic. For instant, Doll and Schweizerhof (2000) developed the volumetric part of the strain energy function and investigated new volumetric functions. Montella et al. (2014) presented the mechanical behavior of a Tire Derived Material in details numerically and experimentally. The problem of the finite axisymmetric deformation of a thick-walled circular cylindrical elastic tube subject to pressure is formulated for an incompressible isotropic neo-Hookean material by Zhu et al. (2010) and solved numerically by finite element library Libmesh. Tanveer and Zu (2012) presented finite amplitude transient vibration analysis of nearly in compressible hyperelastic axisymmetric solids by a mixed p-type method and solved the equations by the Newmark's method along with the Newton-Raphson iterative technique for Mooney-Rivlin material description. Kiendl et al. (2015) presented formulations for compressible and incompressible hyperelastic thin shells with plane stress condition based on energy methods and used continuous iso-geometric discretization to build the numerical solution. The common problems with mention numerical methods, as Poisson's ratio approaches 0.5, are the ill conditioning of stiffness matrix, the locking phenomena and effect of applying numerical techniques on resulted displacements and stresses.
Functionally graded materials (FGMs) are special composites, in which material properties are varied continuously and smoothly through certain direction. Thus, the discontinuities between the layers which occur in layered composites and cause stress concentration are not observed in FGMs. Graded rubber like materials attracted the attention of researchers for modeling these materials behavior under mechanical and geometrical boundary conditions. For instance, effects of material inhomogeneities on stress distributions through the thickness of circular cylinders made of rubber like materials in mechanical and thermal load was studied by Bilgili (2003). In another study, Bilgili (2004) investigated plane strain deformations of a circular cylinder made of heterogeneous neo-Hookean material with circumferential displacements prescribed on the inner and the outer surfaces. Batra and Bahrami (2009) considered cylindrical pressure vessel made of FG rubber like material under internal pressure. To discover stress components of the pressure vessel, they assumed axisymmetric radial deformations of a circular cylinder composed of FG Mooney-Rivlin material with the material parameters varying continuously through the radial direction either by a power law relation. Rahimi (2015, 2016) studied behavior of spherical and cylindrical shell made of FG rubbers by neo-Hookean model. They assumed radial variation of material properties by power law function and used classical theory (PET) and Gausshypergeometric function to derive and solve equations, respectively. Geometrically nonlinear dynamic behavior of FG thick hollow cylinder under axisymmetric mechanical shock loading is investigated using meshless local Petrov-Galerkin (MLPG) method by Ghadiri Rad et al. (2015). The FG cylinder is assumed to be made of large deformable neo-Hookean materials such as carbon-based polymers.
In optimizing a shell with respect to weight or stress distribution, one method is to use shells with varying thickness or materials properties. The literature that addresses the stresses of thick cylindrical shells with variable thickness is quite limited. Eipakchi (2010) calculated stresses and displacements of linear elastic conical shell with varying thickness under non-uniform internal pressure analytically, using shear deformation theory (SDT). Ghannad et al. (2013) presented a closed-form analytical solution for thick FGM cylindrical shells with variable thickness subjected to constant internal pressure based on the first-order shear deformation theory (FSDT) and solved the governing equations by the usage of perturbation theory. Gharooni et al. (2016) investigated thermo-elastic analysis in pressurized thick FGM cylinders with varying properties of power function based on higher-order shear deformation theory. The innovative formulations for higher-order approximation with FG function of materials properties have been presented in this research. Jabbari et al. (2016) investigated thermo-elastic analysis of rotating truncated conical shells with varying thickness made of functionally graded materials (FGMs) subjected to thermo-mechanical loading. The system of partial differential equations is semi-analytically solved by using multi-layered method (MLM). Nejad et al. (2015) presented semi-analytical solution for elastic analysis of axially functionally graded rotating thick cylindrical shells with variable thickness under non-uniform pressure by the usage of SDT and MLM.
Investigating aortic aneurysm as pressurized hyperelastic blood vessels enable scientists to evaluate the relative sensitivity of displacement and stress to geometrical and mechanical properties of the aneurysmal tissue. Furthermore, in pathologic conditions, arteries are even under more shear deformation compared to healthy vessels (Azar et al. 2018).
In clinical interventions, such as balloon angioplasty significant wall shearing may take place. Simulation of arteries under blood pressure to yield displacement and stress analysis of blood vessels could result in useful information on the behavior of the arterial tissues under shear deformation (Vossoughi and Tozeren 1998). Mihai and Goriely (2017) investigated the physical responses of nonlinear elastic materials that are generally described by parameters which are scalar functions of the deformation. They established relations between various hyperelastic material model parameters which are used to quantify nonlinear elastic responses in several hyperelastic models as rubber to soft tissues.
Although numerous studies have been carried out on nearly incompressible hyperelastic shells, no study has been carried out to date on non-uniformly pressurized cylinder with nonlinear variable thickness made of hyperelastic FGMs. In the current study, nonlinear quasi-static analysis of thick cylindrical pressure vessels with arbitrary variable thickness made of Mooney-Rivlin model of hyperelastic FGM with radially variation of material properties in nearly incompressible state under non-uniform pressure loading is presented. The variation of the thickness and pressure profiles of the vessel are considered in axial direction by linear and/or nonlinear functions. In order to improve the approximation and take into account the effect of shear stresses and strains, the general method of derivation and nonlinear analysis has been presented by using first-order shear deformation theory. Matched Asymptotic Expansion (MAE) of the perturbation theory is used for solving the governing system of nonlinear coupled differential equations with variable coefficients. A new ingenious formulation and parameters have been defined during current study to simplify and abbreviate the representation of inner and outer equations components in MAE. In addition, the terms of variable thickness and nonuniform pressure have been presented in special representation separately. The effect of materials constants, inhomogeneity index, geometry and pressure parameters on displacements, stresses and hydrostatic pressure distributions resulting from MAE solution have been investigated for some case studies and the results have been compared with a FEM modeling in ANSYS software. We present the equations that provide the general continuum description of the deformation and the hyperelastic stress response of the material. Current study aims to illustrate the performance of the potentials and their reliability for the prediction of the state of deformation and stress in hyperelastic FG vessels from rubbers to arteries.

Shear deformation theory
Consider a thick-walled axisymmetric cylindrical shell with variable thickness (of the outer surface) subjected to nonuniform internal pressure in the reference configuration as Figure 1. Geometry of the shell could be described in the terms of cylindrical polar coordinates r , θ and x : where i r and ( ) o r x , respectively, are the inner and outer radius and L is the length of the shell. The parameter r is the radius of every layer of cylinder in the reference configuration which can be replaced in terms of radius of mid-plane ( ) R x and distance of every layer with respect to mid-plane ( ) z : where ( ) h x is the thickness of the cylinder which is varying along axial direction. The following relations can be written for the geometry components of the shell: New nonlinear solution of nearly incompressible hyperelastic FGM cylindrical shells with arbitrary variable thickness and non-uniform pressure based on perturbation theory The general axisymmetric displacement field, in the first-order Mirsky-Hermann's theory could be expressed on the basis of radial displacement z U and axial displacement x U , as follows where ( , , Jacobian which is known as volume ratio of deformation has the following terms (det is determinant operator): The Green-Lagrange strain tensor can be defined as  2  2  2  2  2  2  2  zz  xx   2   zx  zx  2   u  w  u  u  w  z  z  2  2  2  2  2  2   w  z  w  z  2

Hyperelastic FGM
Hyperelastics (as rubber-like materials) are kind of materials in which the stresses are only dependent on the initial and the final configurations but independent of the load path. These materials are characterized by the existence of stored energy function which depends on the right Cauchy-Green deformation tensor through the strain invariants . For incompressible material, the determinant of deformation gradient is equal to unity ( ) = J 1 . In the present study, the extension of incompressible materials to nearly incompressible materials is considered; means that the incompressibility constraint is replaced with a penalty function correspond to the constraint. Therefore, the following strain energy function is developed based on the invariants , 1 2 I I and J (Ghaemi et al. 2006): G 2 C C , the variations of these two constants could be considered with powerlaw distribution continuously and smoothly in the radial direction. Although the variations of the FGM layers are the function of radial direction in the current research, because of outer radius variation along axial direction, the properties of outer points are determined based on thickness profile. These constants have the following form (Batra andBahrami 2009, Anani andRahimi 2016): where 1 C and 2 C are material constants at the internal surface and n is the inhomogeneity index determined empirically. In the second term of the right hand side of Eq. (10), we can write: (13) is a penalty function which has to satisfy the conditions ( ) = ⇔ = G J 0 J 1 and > 0 λ is a penalty parameter which can be estimated by experimental data proportional to the material properties and is known as compressibility parameter (Silva andBittencourt 2008, Ghaemi et al. 2006). Therefore, the assumption of almost incompressible rubber is accomplished by dropping the restriction = J 1 and including a hydrostatic work term in the strain energy function. Considering the compressibility parameter as = k λ , where k is an additional material constant representing the bulk modulus, only scales the penalty functions but does not change their shapes. In this context k can be interpreted as a penalty parameter that enforces incompressibility if large values are chosen. In this case, k is the ratio of the volumetric stress, known as hydrostatic pressure ( ) P , to the volumetric strain (Mihai and Goriely 2017 ν which is constant in shell. The order of graded compressibility parameter n k (bulk modulus) can be estimated based on compressibility intensity in FGM rubber as (Ghaemi et al. 2006, Tanveer and Zu 2012, Kiendl et al. 2015: Therefore, FG bulk modulus can be considered similar to Eq. (12).
k is the bulk modulus at the internal radius. For the volumetric part, there are many forms proposed by researchers, all of which are functions of the bulk modulus and the Jacobian of deformation. Generally, in the limiting state, the volumetric part has to satisfy (Doll and Schweizerhof 2000). Considering zero values of displacement components in the reference configurations (initial state) with Eqs. (5), (6) and (7) Finally, the strain energy per unit undeformed volume of FG hyperelastic material for two-term Mooney-Rivlin model in nearly incompressible condition is expressed as the following coupled form (Holzapfel 2000).
Consequently, constitutive equation of coupled Mooney-Rivlin model in material description and nearly incompressible, isotropic and non-homogenous conditions would result (Holzapfel 2000, Başar andWeichert 2000).  (Holzapfel 2000, Bijelonja et al. 2005. Thus, the multiplier 0n p in the case of ( ) = H J ln(J ) denotes the pressure measured in the initial volume. The components in the right hand side of Eq. (19), other than identity tensor and material constants, can be written in the displacement field components. Therefore, the relation between the second Piola-Kirchhoff stress tensor and displacement components could be derived.

Governing equations
Considering the boundary surface 0 A of the body consists of two parts y 0 A and σ 0 A where displacements y  and forces t  0 are prescribed, respectively, based on the principle of virtual work, the variation of strain energy of the elastic body is equal to the variation of external work due to loading (Kim et al. 2012, Başar and Weichert 2000, Doghri 2000. ( ) where 0 ρ , 0 V and b  0 are density, volume and body force in undeformed configuration, respectively. a  is dynamic acceleration. [ ] P and [ ] F are first Piola-Kirchhoff stress and deformation gradient tensors, respectively. Kinematically admissible virtual deformation variables are understood to be the variations y  δ and F δ which are subjected to the constraints (GRAD is gradient operator) Therefore, the principle of virtual work (Eq. (20)) is a weak formulation of the equations of motion as well as the dynamic boundary conditions. Equality of energy conjugate variables preserves its validity i.e. : : = P δ F S δ E . In the static equilibrium and absence of body forces, Eq. (20) results in the variation of external work consists of non-uniform internal pressure ( ) i P x applying at the internal surface i A of cylinder: i dA is the internal surface element. Considering displacement components from Eq. (4), we can rewrite Eq. (22): The internal virtual work in material description can be expressed from Eq. (20) and energy conjugate variables dV is the volume element. Considering Voigt notation from Eq. (9), the variation of strain energy of cylinder with variable thickness can be derived based on non-zero physical components of second Piola-Kirchhoff stress: The variation of strain tensor components can be calculated by the usage of variational principles: The stress resultants are defined as follows: In the last equation, S K is shear correction factor which is applying in the stress resultant derived from shear stresses because of preventing stress overestimation. We consider / = S K 5 6 in the present study (Ghannad et al. 2013). Substituting strain invariants from Eq. (26) into Eqs. (23) and (25) and carrying out the integration by parts, the equilibrium equations of nonlinear hyperelastic cylindrical shell with variable thickness under non-uniform internal pressure are obtained:

Perturbation theory
In this article, Boundary Layer Method of the perturbation theory which is known as Matched Asymptotic Expansion is used for solving the governing equations. The advantages of this method are fast convergence, closed form solution and compatibility with physics of shell. MAE can explain the behavior of the shell successfully even near the boundaries. The governing equations (32)-(35) for cylinder with variable thickness is a system of four nonlinear coupled differential equations with variable coefficients. Preliminary definitions, simplifications and preparations are necessary for using MAE. At first, it is necessary to convert the equations into dimensionless form for making use of the characteristic scales. The following dimensionless parameters are defined (Eipakchi 2010, Nayfeh 1981: Dimensionless quantities of parameters are marked with over-bar ( ) from now. Perturbation parameter ε should Latin American Journal of Solids and Structures, 2019, 16(8), e229 9/28 be normally considered as the ratio of minimal geometrical dimension of the structure ( ) 0 h to maximal one ( ) L in order to be small quantity. 0 h is known as the characteristic thickness which is commonly consider the smallest thickness of shell. The main idea of perturbation theory is that perturbation parameter is so small that coefficients of different power of it don't have the same order which is lead in equality for i ε coefficients. Considering each coefficient result in displacement quasi-vector ( ) { } y x . Existence of two boundary layer lead in two region of solution near boundaries (inner expansions) and a solution away from boundaries (outer expansion) (Nayfeh 1981).
In the dimensionless forms, first and second order differential based on x should be rewrite as follows: In shear deformation theory, the differentials and integrations are with respect to x and z , respectively. Therefore, for simplification and abbreviation of representing equations, a dimensionless integral could be defined: which is the function of geometrical parameters ( ) R x , ( ) h x and inhomogeneity index n . In order to solve the set of governing differential equations in next steps, the inverse of singular coefficient matrices are needed. Two changes should be applied in equations to replace these matrices with non-singular ones. We take integrate the first equation in the set of Eqs. (32)-(35). The constant of integral ) ( 0 c should be either dimensionless as other parameters. Therefore, it is replaced with its dimensionless equivalent where 7 c is integral constant. The intentional constants 0 c and 7 c will be calculated from boundary conditions. The following parameters need to be defined based on material constants of internal layer , ,

Outer expansion
The outer expansion of solution is considered a uniform series of ε as ( ) stands for outer solution. Furthermore, the subscript "1" and "2" shows first and second order expansion, respectively. Substituting the expansion in governing equations and considering the terms with the same order of ε , result in the first and second order equations of outer solution. In this section ( ) ( ) correspond to ( ) II x and derivative of ( ) II x , respectively. These vectors would be defined in the Appendix A. Other non-zero components of the mentioned matrix and vectors are The solutions of the algebraic equations (40) are as follows:

Inner expansion
As the outer solutions don't satisfy the B.C., we can conclude existence of boundary layers at x α should be considered as a new variable for regions around boundaries. Considering fast variables make it possible to measure the great variation of mechanical response around boundaries.
Subscript " x α in the boundary. Perturbation parameter appears in new definitions of differential based on fast variables  x α : In cylinder with variable thickness and non-uniform pressure, it is necessary to derive Taylor expansion for all the parameters of axial function ( ) Λ  x as: Taylor expansion with expansion. Taylor expansion should be written for the following parameters: in governing equations with mentioned changes in section 3.1 and considering terms with the same order of ε result in inner equations at boundary α :    3  3  s 12  3  3  1  11  22  13  31   3  3  1  3  14  41  33   3  3  1  3  1  34 43 44 k II 0 n 1 K C II 0 n 1 Ck II 0 n Ck II 0 n 1 II 1 n k II 0 n 1 Ck II 0 n k II 1 n 1 2Ck II 1 n k II 0 n 1 II 2 n 1 The necessary condition for existing the solution of Eq. (56) where subscript "i" and "j" shows the number of equal and unequal roots with characteristic equations. Substituting { } y 2 par.

Composite solution
In

FE modeling
In order to validate presented analytical solution and compare the results for pressurized thick cylinder with variable thickness made of nearly compressible FG hyperelastic material, a numerical solution based on Finite Element Method is investigated. The ANSYS 16 package is used in the static analysis of thick hollow cylinder with variable thickness under non-uniform internal pressure. The PLANE183 element in the axisymmetric mode, which is an element with eight nodes and two translational degrees of freedom in the axial and radial directions per each node, has been used to model the mechanical part of the analysis. It also has mixed formulation capability for simulating deformations of nearly incompressible hyperelastic material. The cylinder with variable thickness consists of some coherent homogeneous layers which properties at the contact location of the layers are the average of left and right limit of two layer boundaries. In order to model FG hyperelastic cylinder, an innovative application for multilayering the thickness in the axial direction has been performed. Homogenous layers which have identical thickness and step-variable properties have been formed by this method. In order to consider Mooney-Rivlin elastic model in nearly incompressible condition in each homogenous layer, three constants involving , 10 01 C C and d should be defined for ANSYS software. 10 C and 01 C corresponding to 1n C and New nonlinear solution of nearly incompressible hyperelastic FGM cylindrical shells with arbitrary variable thickness and non-uniform pressure based on perturbation theory  (Montella et al. 2014, Gao et al. 2009) which could be calculated from Eq. (16) in each layer. For nonuniform internal pressure, the pressure functions have been defined and applied to the internal layer nodes. Clamped boundary conditions have been exerted by preventing the nodes around the two ends of the cylinder from movement. In the next sections, the numerical results (FEM) and analytical results (MAE) have been investigated for different case studies.

Case studies
In order to illustrate the effects of material constant, gradient index, pressure loading distribution and thickness profile type on the mechanical behavior of hyperelastic pressure vessel, thick cylinders with various non-uniform thickness and internal pressure profiles made of different materials constants and inhomogeneity index have been considered. In these example problems, unless noted otherwise, we take constant geometry parameters as: . Clamped boundary conditions are applying in two ends of the shells. The materials are assumed functionally graded hyperelastic with tow-term Mooney-Rivlin model in nearly incompressible conditions. During the computation of numerical results, different material contestants according to various references are considered which have been presented in Table 1. The material constants are considered based on two state: 1) mentioned directly in these articles or 2) computed by authors of current paper from different test results (mentioned in these references). The constants of FG Mooney-Rivlin model are considered for material properties at the internal surface. Considering the values of 1 C and 2 C in Table 1 and Eq. (15) result the variation range of k about − 1 100 MPa in nearly incompressible limit. Important remark is that for < k 2 MPa , the eigenvalues of characteristic equation no longer have conjugate complex form and MAE solution diverge. Hence, in Table 1, we set = k 10 MPa for cases with no incompressibility parameter presentation (Dias et al. 2014, Kiendl et al. 2015. Internal pressure distribution and thickness profile varies non-uniformly along axial direction of shell. Various pressure profiles are applied to the cylindrical shell in the range of kPa kPa − 5 13 . Dimensionless Cauchy stresses and hydrostatic pressure are defined as: . Tables 2 and 3 show the characteristic of applied non-uniform pressure and thickness profiles, respectively.

Pressure Profile ID Pressure Profile Relation Load Constants
P P 13 P 5 Table 3: The characteristic of non-uniform thickness profiles.

Thickness Profile ID Thickness Profile Relation Geometry Constants
As current research studies the manner of pressurized rubber vessels in dimensionless state, the results of FSDT and MAE solution may be suitable for investigating some case studies of blood vessels. In particular, we use an analysis to examine the inflation of a cylindrical tube at various internal pressure profiles and to compute the evolution of the inner radius (critical layer) with the internal pressure. Although the material models of blood vessels such as arteries have commonly exponential form to model the stiffening of the soft tissues, simple model such as neo-Hookean and/or Mooney-Rivlin one is considered in some research for isotropic parts of blood vessels (Holzapfel andGasser 2007, Lally et al. 2005). Hence, case studies of variable thickness and pressure are selected similar to that of common elastic arteries (as unite layer) in current research in order to cover pressure vessels to blood vessels. The thick and thin part of the vessels covers the average thickness of the layers in common elastic arteries. The pressure profiles vary in the range of kPa( mmHg) kPa( mmHg) is the mean of systolic/diastolic pressure and mmHg 40 may be occurred in hypotension pressure of arteries (Abdessamad et al. 2018, Humphrey andO'Rourke, 2015). Considering the applicability of the rubber elasticity theory to aortic soft tissues as one layer or multilayer vessel with variable material properties along thickness, the behavior of blood vessels under nonuniform pressure distribution has been investigated from current research. Furthermore, current study will present helpful results for estimating the wall degeneration of arteries within the aneurysm wall that affects the thickness profile of the tissue, which can be mostly analyzed as variable thickness blood vessels (Xie et al. 1995, Vossoughi and Tozeren 1998, Holzapfel and Gasser 2007.

Effect of material constants and inhomogeneity index
In order to investigate the material constants effect on the approximation of the current solution and behaviour of shell, maximum values of radial and axial displacements for different material constants resulted from FSDT and FEM is presented in Tables 4 and 5. As the maximum values of displacements and maximum difference of analytical and numerical analysis in pressure vessels occur at the internal layer, the results of Tables 4 and 5  P C C are less than 8%. R values of more than 20 (thickness limit for thick cylinder) and less than 4 (very thick shells) may lead in decrease in solution accuracy. It is observed that the accuracy of MAE descend for great values of R because of intensifying nonlinear behavior of the cylinder while for small R , the accuracy of shear deformation theory decrease in analyzing thick cylindrical shells. The main reason of increasing the difference between results of shear deformation theory and other theories (finite element or plane elasticity) for small R are the effect of domination of thickness values to displacement ones which make linear (or even more order) distribution of displacements along the thickness of cylinder not to be matched to real state; so shear deformation approximation cylinders lead in less accuracy for very thick (Eipakchi et al. 2003). Ascending R lead in higher deformations and lower effect of clamped conditions on deformations and the peak points of displacements occurs far away from boundaries. An increase in i ( ) + 1 2 P C C and k or decrease in + 1 2 C C ascend the nonlinearity and descend the accuracy.
The longitudinal and circumferential components of Green-Lagrange strain are depicted in Figure 2 for homogeneous cylinder. According to Figure 2(a), the longitudinal strain resulted from numerical and analytical solution show good agreement along axial direction of different layers. It can be seen that maximum longitudinal strain occurs around the boundaries of external layer. However, the position of extremum strains could be different for low range of loading. At the middle of the cylinder away from boundaries, no significant difference of longitudinal strain is observed along the thickness. As the maximum values of circumferential strain occur at the middle of the cylinder ( ) . = x 0 5 , distribution of this strain component is shown along thickness of cylinder for different loading in Figure 2(b). Ascending the internal pressure causes decrease in accuracy of analytical strain especially around the loading layer with higher strain (internal layer). FSDT have acceptable accuracy for indirect calculation of strain components from displacements.   U U in the current case study. Because of greater radial displacements around left boundary, it can be concluded that pressure profile is more effective than thickness variation on shell displacement. According to Figure 4(b), the axial displacement at points away from the boundaries is nearly independent from radius. Layers close to maximum pressure and clamped conditions are in axial tension; therefore, axial compression is dominants at the shell except in ≤ z 0 around left boundary. Figure 5 show dimensionless Cauchy stresses and hydrostatic pressure at different layers along axial direction. Hydrostatic pressure can be considered as average of principal stresses. Considering Figure 5 confirm this fact; so hydrostatic pressure can be a suitable equivalent parameter that show shell state from the view point of stresses. Circumferential and axial stresses, similar to hydrostatic pressure, have positive values in nearly all points of the shell except around boundaries at the outer layer away from loading. The reason is that the elements are in tensile state, but clamped conditions near boundaries at the layer away from loading cause resistance against displacement which lead in compressive stresses. In this state, inner layer of the shell in contact with pressure load have higher displacements and stresses than others. It is obviously observed that the circumferential stress is the largest component of the stress at points away the boundaries while at the points near the boundaries, the axial stress is the largest one. Existence of shear stress near boundaries reveals the advantage of shear deformation theory respect to theories that neglect shear stress effect ( Figure 5(c)). Getting away from boundaries, non-uniform peaks of displacements and stresses are observed at the points where shear stresses tend to zero values. Difference between MAE and FEM results increase at the points of internal and external layers away from boundaries. Although FSDT is suitable for displacement analyzing rather than stress one, the results of MAE are more realistic around boundaries respect to FE solution. Considering Eq. (14) and hydrostatic pressure distribution lead in . is suitable for nonlinear problems. Figure 6 shows that second order solution improves results accuracy respect to first order one, especially in the current research that kinematics and constitutive relations are highly nonlinear.    In order to investigate the effect of inhomogeneity index on the response of shell, the distribution of displacements and hydrostatic pressure at internal layer (critical layer) are plotted in Figure 7(a-c). The linear thickness and pressure profile are assumed for shell with 9 MC material constants. Figure 7(d) shows the distribution of dimensionless material properties (normalized to internal layer properties) with respect to the radius variation in a heterogeneous cylinder for integer values of n which vary in the range of − ≤ ≤ + 4 n 4 . The variations of material properties with power-law distribution are continuously and smoothly in the radial direction. The extremum values of properties at outer layer points could be determined through intersection of vertical line plotted from radius of the point and graph of arbitrary n . Positive values of gradient index increase strength of material under mechanical loading toward outer layer of shell, while the reverse holds true for negative values of n . Therefore, variation of inhomogeneity index from negative to the positive causes displacements, hydrostatic pressure and consequently stresses of cylinder to be reduced. Greater values of n intensify improvement or reduction in response of FG shell respect to homogenous one. Table 6 presents the results of maximum displacements and hydrostatic pressure in three layers of shell for different inhomogeneity index. Linear decrease in radial displacement and smooth reduction in axial displacement can be observed from internal layer to the external one for different values of n . Positive values of n cause more uniform hydrostatic pressure distribution of the layers and less maximum values of hydrostatic pressure compared with negative ones. Therefore, It could be concluded that internal layer (in contact with loading) is still critical one and positive values of n are more appropriate from the viewpoint of less values and more uniform distribution of displacements and stresses in heterogeneous cylinder.
However, we are not concerned here with possible manufacture processes of FG materials, as well as experimental tests. Authors believe that, when the FG elastomers start to be widely employed in industry or in engineering applications, our formulation is a reliable numerical tool to predict their mechanical behavior (in terms of accuracy). But it is important to note that method presented here will be useful to material scientists in designing new materials, stress analysts, and designers in two states. One can use similar solution procedure to calculate displacements and stresses for FG material models with the given constants functions applied instead of Eq. (12) distribution. Furthermore, one can control the through-the-thickness distribution of displacements and stresses as objective parameters by tailoring the through-the-thickness variation of the material constants by trial and error to achieve appropriate distribution of FGM constants. In the material tailoring problem, one has found through-the thickness variation of material constants to achieve a desired variation of stress components, frequency of free vibrations, deformation or an objective function to be optimized (Batra 2011). This method is going to be extended in FG elastomers and biological tissues. Bilgili (2004) also suggest that the presence of material non-homogeneity in test specimens might be reason for the conflicting experimental results in the technical literature regarding the nature of the rubber-elastic response functions. Hence, he developed comprehensive experimental and theoretical program to characterize the response functions of nonhomogeneous rubber components and introduced a design code based on especial (power variation) material model which can explain the essential physics-chemistry behind the intended functionality. The design code yields essential information about the grading which in turn can be used as input into the design of a fabrication process. Thus, our method along with mentioned studies could direct further research toward the design, optimization, and manufacture of graded rubber-like materials.     Figure 8 shows the effect of pressure profile on the distribution of displacements and hydrostatic pressure in FG cylinder with linear variable thickness. The material is considered with = n 2 and 9 MC . Distribution of non-uniform internal pressure functions along axial direction are depicted in Figure 8(d). Investigating the response of cylinder with linear thickness profile under i1 P , i 2 P and i 3 P in Figure 8(a) and (b) show nearly the same maximum displacements; however, variations of radial displacement are proportional to pressure profiles distribution at length of the shell. It is observed that linear variation of pressure i1 P and thickness 1 h of the shell in same direction counteract each other effect. This counterbalance is weakly true for nonlinear pressure profiles with similar range of applied pressure. Figure 8(c) expresses that hydrostatic pressure under these three pressure profiles increase by intensifying nonlinearity of pressure distribution. Hydrostatic pressure of middle layer has its maximum value near right boundary unexpectedly, because the effect of descending thickness is dominant to pressure profile increments. For internal (critical) layer reverse hold true (as previous results proved), i.e. the maximum hydrostatic pressure occurs around left boundary with higher pressure. Pressure profiles of i 4 P and i 5 P which have higher non-uniformity respect to other profiles lead in larger radial displacement and hydrostatic pressure values, but less axial displacement. Fast pressure variation of i 4 P and i 5 P reduce z U and P around = x 0 . Maximum pressure applied in middle and end of cylinder for i 4 P and i 5 P , respectively, in addition to descending thickness cause maximum of z U and P near = x 1 . It can be concluded from Figure 8(a) and (c) that z U and P patterns along the length of shell follow the pattern of the applied pressure function. MC materials properties. Considering Figure 9(a) and (c) prove that changes of concave thickness profile to convex one cause reduction in radial displacement and hydrostatic pressure. However, z U and P for 4 h profile intensify near = x 0 because of lower thickness in addition to maximum pressure. It is obviously observed that constant thickness ( ) 0 h have the greatest displacements and stresses under i1 P . No considerable variations in axial displacements are observed between non-uniform thickness profiles. Comparison of Figures 8-11 reveal that pressure profiles increment is more effective on the response of shell respect to thickness profiles variation; i.e. descending pressure causes more reduction in displacements and stresses respect to increasing the thickness. Table 7 represents the numerical results for similar distribution of pressure and thickness in different sections of internal layer. The axial sections are selected based on extremum points of displacements and hydrostatic pressure distribution. It is observed that similar thickness and pressure patterns lead in uniform distributions of displacements, stresses and hydrostatic pressure along length of the shell. According to Table 7 profiles have the least and the most hydrostatic pressure, respectively. In fact, whatever maximum pressure and consequently maximum thickness are applied near clamed boundaries, more counterbalance of thickness and pressure emerges. It can be concluded that i , 5 5 P h and i , 1 1 P h are suitable profiles in designing current hyperelastic FG shell. The rupture modes in the aortic specimens are characterized by oblique tears in the circumferential direction, indicating that the failure of the aneurismal aortic tissue is mainly governed by the axial stress. The failure stress in the axial direction is much higher in the adventitia layer compared to that in the media layer (Kim et al. 2012). This means that the failure in the aneurismal aortic tissue may initiate in the media layer; i.e. inner surface of arteries are critical one. It is considered that the current methodology could be improved to assess the aortic aneurysm rupture risk based on maximal diameter or stress by modeling blood vessels of patients having an aneurysm. In recent years, researchers have developed artificial blood vessels made from special elastomer material and the usage of artificial vascular prostheses in vascular graft (Łos et al. 2018). Over time, these artificial blood vessels are replaced by endogenous material. Some parts of mentioned prostheses don't have complicated geometries and can be models as regular shells with acceptable tolerances and imperfections. Authors believe that current method could have the potential of helping researchers in the future to analyze and obtain useful information about (a) more realistic hyperelastic material models of blood vessels (artificial or natural, isotropic or anisotropic, homogenous or non-homogeneous); (b) especial variation in internal and /or external profiles of blood vessels (as variable thickness conical shell) resulted from atherosclerotic plaque, aortic aneurysm, aging deformation and so on; (c) maybe future non-homogeneous prosthesis with position dependent functionality.

CONCLUSIONS
In current research, the heterogonous hyperelastic hollow cylinders with variable thickness under non-uniform internal pressure and clamped boundary conditions have been analyzed by FSDT. Two-term Mooney-Rivlin type material in nearly incompressible condition is considered which is a suitable hyperelastic model for rubbers. The material properties are graded along the radial direction according to a power law function. Matched Asymptotic Expansion of the perturbation theory is used for solving the governing equations analytically. The advantages of this method are fast convergence, closed form solution and compatibility with physics of shell. A new ingenious formulation and parameters have been defined during current study to simplify and abbreviate the representation of inner and outer equations components in MAE. In addition, the terms of variable thickness and non-uniform pressure have been presented in separate representation. The results prove the effectiveness of FSDT and MAE combination to derive and solve the governing equations of nonlinear problems such as nearly incompressible hyperelastic shells. This approach enables insight into the nature of the deformation and stress distribution across the wall of rubber vessels and offers the potential for investigation of the mechanical functionality of arteries in physiological pressure range. The points of internal layer close to maximum pressure are critical elements in nearly incompressible hyperelastic FG cylinder with variable thickness under non-uniform internal pressure. The acceptable range of the current analysis for the geometry, loading and materials properties is about < < 4 R 20 and ( ) i .

+ <
1 2 P C C 0 015 by considering difference percentage of deformations resulted from current analytical solution and FEM less than 10%. The accuracy of MAE descend for great values of R because of intensifying nonlinear behavior of the cylinder while for small R , the accuracy of shear deformation theory decrease in analyzing thick cylindrical shells. An increase in i ( ) + 1 2 P C C and k or decrease in + 1 2 C C ascend the nonlinearity and difference percentage of numerical and analytical solution. Variation of inhomogeneity index from negative to the positive values causes reduction in displacements, hydrostatic pressure and consequently stresses of cylinder. Therefore, it could be concluded that positive values of gradient index are more appropriate from the viewpoint of less values and more uniform distribution of displacements and stresses in heterogeneous cylinder. It can be concluded that pressure profiles increment is more effective on the response of shell respect to thickness profiles variation. Furthermore, changes of concave thickness profile to convex one lead in descending maximum displacement, stresses and hydrostatic pressure. It can be concluded that radial displacement and hydrostatic pressure patterns follow the pattern of the applied pressure function along the length of shell. The behavior of hyperelastic FG vessels under non-uniform pressure distribution shows that similar profile of variable thickness and non-uniform applied pressure result in minor displacement and stress quantities and uniform distributions which could be a suitable criterion in designing thickness profile of pressurized vessels. Applying maximum pressure and consequently maximum thickness near boundaries of shell are suitable profiles in designing hyperelastic FG shells. Authors believe that current method along with studies mentioned in the literature could direct further research toward the design, optimization, and manufacture of graded rubber-like materials.