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

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 FG shells.

Keywords:
Variable thickness shells; Nonlinear perturbation solution; Cylindrical pressure vessel; Hyperelastic FGM; Mooney-Rivlin model

1 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 anti-vibration 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. Anani and 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 Gauss-hypergeometric 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 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.

2 BASIC FORMULATIONS

2.1 Shear deformation theory

Consider a thick-walled axisymmetric cylindrical shell with variable thickness (of the outer surface) subjected to non-uniform internal pressure in the reference configuration as Figure 1. Geometry of the shell could be described in the terms of cylindrical polar coordinates $r$, $\theta$and $x$:

$r_{\text{i}}\le r\le r_{\text{o}}(x),0\le \theta \le 2\pi ,0\le x\le L$ (1)

where $r_i$ and $r_o(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)$:

$r=R(x)+z,-\frac{h(x)}{2}\le z\le \frac{h(x)}{2}\Rightarrow dr=dz,(r,\theta ,x)\Rightarrow (z,\theta ,x)$ (2)

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:

$R(x)=r_i+\frac{h(x)}{2},r_o(x)=r_i+h(x)$ (3)

Figure 1:
Geometry and loading parameters in cross section of the shell.

The general axisymmetric displacement field, in the first-order Mirsky-Hermann's theory could be expressed on the basis of radial displacement $U_z$ and axial displacement $U_x$, as follows

$U_z(z,x)=w(x)+z\psi (x),U_{\theta }=\mathrm{0,}{U}_x(z,x)=u(x)+z\phi (x)$ (4)

where $w(x)$and $u(x)$are the displacement components of the middle surface. Also, $\psi (x)$and $\phi (x)$are the rotational components used to determine the displacement field.

The matrix representation of deformation gradient tensor $\left[F\right]$ in geometrically nonlinear kinematic is (Zhu et al. 2010)

$\left[F\right]=\left[\begin{array}{ccc}1+\psi & 0& w^{\prime }+{\psi }^{\prime }z\\ 0& 1+\frac{w+\psi z}{R+z}& 0\\ \phi & 0& 1+u^{\prime }+{\phi }^{\prime }z\end{array}\right]$ (5)

where${\left(\right)}^{\prime }=\partial \left(\right)/\partial x$. Consequently, the right Cauchy-Green deformation tensor, $\left[C\right]={\left[F\right]}^T\left[F\right]$ and its principal invariants $I_{\mathrm{1,2,3}}$ in cylindrical polar coordinate basis are given by (Tanveer and Zu 2012)

$\{\begin{array}{l}\left[C\right]=\left[\begin{array}{ccc}{C}_{zz}& 0& C_{zx}\\ 0& C_{\theta \theta }& 0\\ C_{xz}& 0& C_{xx}\end{array}\right],C_{zz}={\left(1+\psi \right)}^2+{\phi }^2,C_{\theta \theta }={\left(1+\frac{w+\psi z}{R+z}\right)}^2\\ C_{xx}={\left(w^{\prime }+{\psi }^{\prime }z\right)}^2+{\left(1+u^{\prime }+{\phi }^{\prime }z\right)}^2,C_{zx}=C_{xz}=\left(1+\psi \right)\left(w^{\prime }+{\psi }^{\prime }z\right)+\phi \left(1+u^{\prime }+{\phi }^{\prime }z\right)\end{array}$ (6)

$I_1=C_{zz}+C_{\theta \theta }+C_{xx},I_2=C_{zz}{C}_{\theta \theta }+C_{\theta \theta }{C}_{xx}+C_{zz}{C}_{xx}-C_{zx}^2,I_3=C_{zz}{C}_{\theta \theta }{C}_{xx}-C_{\theta \theta }{C}_{zx}^2=J^2$ (7)

Jacobian which is known as volume ratio of deformation has the following terms (det is determinant operator):

$J=\det \left(F\right)=\left(\frac{1}{R(x)+z}\right)\left[\left(R(x)+z+w+\psi z\right)\left(\left(1+\psi \right)\left(1+u^{\prime }+{\phi }^{\prime }z\right)-\phi \left(w^{\prime }+{\psi }^{\prime }z\right)\right)\right]$ (8)

The Green-Lagrange strain tensor can be defined as $E=1/2\left(C-I\right)$ ($\left[I\right]$is the identity tensor). Considering Voigt notation of Green-Lagrange strains ${\left\{\epsilon \right\}}^T=\left\{\begin{array}{cccc}{\epsilon }_{zz}& {\epsilon }_{\theta \theta }& {\epsilon }_{xx}& {\gamma }_{zx}\end{array}\right\}$, Its components are as follows:

$\{\begin{array}{l}{\epsilon }_{zz}=\psi +\frac{{\psi }^2}{2}+\frac{{\phi }^2}{2},{\epsilon }_{xx}=u^{\prime }+\frac{{u^{\prime }}^2}{2}+\frac{{w^{\prime }}^2}{2}+\left({\phi }^{\prime }+u^{\prime }{\phi }^{\prime }+w^{\prime }{\psi }^{\prime }\right)z+\left(\frac{{{\phi }^{\prime }}^2}{2}+\frac{{{\psi }^{\prime }}^2}{2}\right)z^2\\ {\epsilon }_{\theta \theta }=\frac{w+\psi z}{R(x)+z}+\frac{{\left(w+\psi z\right)}^2}{2{\left(R(x)+z\right)}^2},{\gamma }_{zx}=2{\epsilon }_{zx}=\phi +u^{\prime }\phi +w^{\prime }+w^{\prime }\psi +\left({\psi }^{\prime }+\psi {\psi }^{\prime }+\phi {\phi }^{\prime }\right)z\end{array}$ (9)

2.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 $\left(W(I_1,I_2,I_3)\right)$. For incompressible material, the determinant of deformation gradient is equal to unity $\left(J=1\right)$. 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 $I_1,I_2$ and $J$ (Ghaemi et al. 2006):

$W(I_1,I_2,J)=W^{*}(I_1,I_2)+\tilde{W}(J)+cH(J)$ (10)

$W^{*}$ in Eq. (10) is the response of material to distortional part of the deformation. In the present study, the hyperelastic material of the shell is assumed to be isotropic and Non-homogenous with two-term Mooney-Rivlin material description in nearly incompressible condition which is a suitable hyperelastic model for rubbers. It has the following form (Batra and Bahrami 2009):

$W^{*}(I_1,I_2)=C_{1n}\left(I_1-3\right)+C_{2n}\left(I_2-3\right)$ (11)

where $C_{1n}$ and $C_{2n}$ are material constants resulting from experimental tests. Non-homogenous and isotropic FG hyperelastic materials have different properties in terms of points. The changes of properties in functionally graded rubbers are generally considered in radial and/or longitudinal directions. As the constants of Mooney-Rivlin model have relation with initial shear modulus $G_n=2\left(C_{1n}+C_{2n}\right)$, the variations of these two constants could be considered with power-law 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 and Bahrami 2009, Anani and Rahimi 2016):

$C_{1n}\left(x,z\right)=C_1{\left(\frac{R(x)+z}{r_i}\right)}^n,C_{2n}\left(x,z\right)=C_2{\left(\frac{R(x)+z}{r_i}\right)}^n$ (12)

where $C_1$ and $C_2$ 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:

$\tilde{W}(J)=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\lambda {\left\{G(J)\right\}}^2$ (13)

$G(J)$ in Eq. (13) is a penalty function which has to satisfy the conditions $G(J)=0\iff J=1$ and $\lambda >0$ is a penalty parameter which can be estimated by experimental data proportional to the material properties and is known as compressibility parameter (Silva and Bittencourt 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 $\lambda =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).

$k=\frac{P}{\Delta V/V_0}=\frac{P}{J-1}$ (14)

$V_0$ and $\Delta V$are the reference volume and volume changes through deformation, respectively. Various relations is recommended for estimating the value of incompressibility parameter; for instant $k=2G(1+\nu )/3(1-2\nu )$. The common part of similar relations is definition of the bulk modulus based on material constants ($C_{1n}$and $C_{2n}$in current model) or initial shear modulus$(G_n)$and Poisson’s ratio $(\nu )$. In nearly incompressible rubber materials, Poisson’s ratio with respect to compressibility intensity consider about $\nu =0.49-0.499$ which is constant in shell. The order of graded compressibility parameter $k_n$ (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):

$\nu =0.49\to k_n\propto \left(C_{1n}+C_{2n}\right)\times {10}^{2}to\nu =0.499\to k_n\propto \left(C_{1n}+C_{2n}\right)\times {10}^3$ (15)

Therefore, FG bulk modulus can be considered similar to Eq. (12).

$k_n\left(x,z\right)=k{\left(\frac{R(x)+z}{r_i}\right)}^n$ (16)

$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 $J\to 1\Rightarrow \tilde{W}(J)\to \mathrm{1,}{\tilde{W}}^{\prime }(J)\to \mathrm{0,}{\tilde{W}}^{″}(J)\to k$conditions (Doll and Schweizerhof 2000).

Considering zero values of displacement components in the reference configurations (initial state) with Eqs. (5), (6) and (7) lead in $\left[F\right]=\left[C\right]=\left[{\rm I}\right]$ and $\left(I_1,I_2,I_3,J\right)=\left(\mathrm{3,3,1,1}\right)$. In the third term of the right hand side of Eq. (10), constant $c$and function $H(J)$ with the condition $H(J)=\mathrm{0,}{H}^{\prime }(J)=1\iff J=1$ only guarantee the stress free reference configuration with no physical meaning. Moreover, in the general case of nearly incompressible hyperelasticity (as Mooney-Rivlin material), hydrostatic pressure $P=k(J-1)$ does not vanish even at the natural state. The first condition $\left(H(J)=0\right)$corresponds to the incompressibility constraint $J=1$ and the second condition $\left(H^{\prime }(J)=1\right)$ is necessary for giving the meaning of pressure to the constant multiplier of $H$as $c=-p_0$. Then the initial value of $P$ (i.e.$p_0$), which has no clear physical meaning, must be introduced to make the initial stress zero. In the current study, the function $H(J)$and $G(J)$are considered as (Doll and Schweizerhof 2000, Ghaemi et al. 2006):

$H(J)=\ln (J),G(J)=J-1$ (17)

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).

$W=C_{1n}\left(I_1-3\right)+C_{2n}\left(I_2-3\right)-p_{0n}\ln (J)+\frac{k_n}{2}{\left(J-1\right)}^2$ (18)

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 and Weichert 2000).

$S=2\frac{\partial W\left(I_1,I_2,J\right)}{\partial C}=2\left[C_{1n}+C_{2n}{I}_1\right]I-2C_{2n}C+\left[k_{n}J\left(J-1\right)-p_{0n}\right]C^{-1}$ (19)

Eq. (19) relates the right Cauchy-Green strain tensor $\left[C\right]$ to the second Piola-Kirchhoff stress tensor $\left[S\right]$ through constitutive relation. $\left[I\right]$ is the identity tensor. The initial stress is zero if the hydrostatic pressure vanishes at the natural state, and vice versa. Recalling the assumption of stress-free reference configuration, Eq. (19) result in $p_{0n}=2\left(C_{1n}+2C_{2n}\right)$ (Holzapfel 2000, Bijelonja et al. 2005). Thus, the multiplier $p_{0n}$ 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.

2.3 Governing equations

Considering the boundary surface $A_0$ of the body consists of two parts $A_0{}_y$and $A_0{}_{\sigma }$ where displacements $\overrightarrow{y}$ and forces ${\overrightarrow{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).

$\delta \Pi =\delta {\Pi }_{EXT}^{}-\delta {\Pi }_{INT}^{}=0={\displaystyle \underset{V_0}{\iiint }\delta \overrightarrow{y}.\left({\rho }_0{\overrightarrow{b}}_0-{\rho }_0\overrightarrow{a}\right)d{V}_0}+{\displaystyle \underset{A_{0\sigma }}{\iint }\delta \overrightarrow{y}.{\overrightarrow{t}}_{0}d{A}_0}-{\displaystyle \underset{V_0}{\iiint }P:\delta Fd{V}_0}$ (20)

where ${\rho }_0$, $V_0$ and ${\overrightarrow{b}}_0$are density, volume and body force in undeformed configuration, respectively. $\overrightarrow{a}$ is dynamic acceleration. $\left[P\right]$and $\left[F\right]$are first Piola-Kirchhoff stress and deformation gradient tensors, respectively. Kinematically admissible virtual deformation variables are understood to be the variations $\delta \overrightarrow{y}$ and $\delta F$ which are subjected to the constraints (GRAD is gradient operator)

$\delta F=GRAD\delta \overrightarrow{y}(in{V}_0)\&\delta \overrightarrow{y}=0(on{A}_{0y})$ (21)

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:\delta F=S:\delta 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 $P_i(x)$applying at the internal surface$A_i$ of cylinder:

$\delta {\Pi }_{EXT}^{}={\displaystyle \underset{A_0}{\iint }\delta y.t_{0}d{A}_0}={\displaystyle \underset{A_i}{\iint }{P}_i(x){\delta U_z|}_{r=r_i}d{A}_i},d{A}_i=r_{i}d\theta dx$ (22)

$d{A}_i$is the internal surface element. Considering displacement components from Eq. (4), we can rewrite Eq. (22):

$\delta {\Pi }_{EXT}^{}={\displaystyle \underset{0}{\overset{L}{\int }}2\pi P_i(x)}\left(R(x)-\frac{h(x)}{2}\right)\left(\delta w-\frac{h(x)}{2}\delta \psi \right)dx$ (23)

The internal virtual work in material description can be expressed from Eq. (20) and energy conjugate variables

$\delta {\Pi }_{INT}^{}={\displaystyle \underset{V_0}{\iiint }S:\delta Ed{V}_0}={\displaystyle \underset{V_0}{\iiint }{S}^{ij}\delta {\epsilon }_{ij}d{V}_0,d{V}_0=r(x)drd\theta dx=\left(R(x)+z\right)dzd\theta dx}$ (24)

$d{V}_0$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:

$\delta {\Pi }_{INT}^{}=2\pi {\displaystyle \underset{0}{\overset{L}{\int }}{\displaystyle \underset{-h(x)/2}{\overset{+h(x)/2}{\int }}\left[S^{zz}\delta {\epsilon }_{zz}+S^{\theta \theta }\delta {\epsilon }_{\theta \theta }+S^{xx}\delta {\epsilon }_{xx}+S^{zx}\delta {\gamma }_{zx}\right]R(x)\left(1+\frac{z}{R(x)}\right)dzdx}}$ (25)

The variation of strain tensor components can be calculated by the usage of variational principles:

$\{\begin{array}{l}\delta {\epsilon }_{xx}=\delta u^{\prime }+u^{\prime }\delta u^{\prime }+w^{\prime }\delta w^{\prime }+(\delta {\phi }^{\prime }+u^{\prime }\delta {\phi }^{\prime }+{\phi }^{\prime }\delta u^{\prime }+w^{\prime }\delta {\psi }^{\prime }+{\psi }^{\prime }\delta w^{\prime })z+\left({\phi }^{\prime }\delta {\phi }^{\prime }+{\psi }^{\prime }\delta {\psi }^{\prime }\right)z{}^2\\ \delta {\epsilon }_{zz}=\delta \psi +\psi \delta \psi +\phi \delta \phi ,\delta {\epsilon }_{\theta \theta }=\frac{\left(\delta w+z\delta \psi \right)}{{\left(R(x)+z\right)}^2}\left(R(x)+z+w+\psi z\right)\\ \delta {\gamma }_{zx}=\delta \phi +u^{\prime }\delta \phi +\phi \delta u^{\prime }+\delta w^{\prime }+w^{\prime }\delta \psi +\psi \delta w^{\prime }+\left(\delta {\psi }^{\prime }+\psi \delta {\psi }^{\prime }+{\psi }^{\prime }\delta \psi +\phi \delta {\phi }^{\prime }+{\phi }^{\prime }\delta \phi \right)z\end{array}$ (26)

The stress resultants are defined as follows:

${\left\{\begin{array}{ccc}{N}_z& M_z& Q_z\end{array}\right\}}^T={\displaystyle \underset{-h(x)/2}{\overset{+h(x)/2}{\int }}{S}^{zz}{\left\{\begin{array}{ccc}1& z& z^2\end{array}\right\}}^T\left(1+\frac{z}{R(x)}\right)dz}$ (27)

${\left\{\begin{array}{ccc}{N}_{\theta }& M_{\theta }& Q_{\theta }\end{array}\right\}}^T={\displaystyle \underset{-h(x)/2}{\overset{+h(x)/2}{\int }}{S}^{\theta \theta }{\left\{\begin{array}{ccc}1& z& z^2\end{array}\right\}}^{T}dz}$ (28)

${\left\{\begin{array}{ccc}{N}_{\theta }{}^{*}& M_{\theta }{}^{*}& Q_{\theta }{}^{*}\end{array}\right\}}^T={\displaystyle \underset{-h(x)/2}{\overset{+h(x)/2}{\int }}{S}^{\theta \theta }{\left\{\begin{array}{ccc}1& z& z^2\end{array}\right\}}^T\left(\frac{R(x)}{R(x)+z}\right)dz}$ (29)

${\left\{\begin{array}{ccc}{N}_x& M_x& Q_x\end{array}\right\}}^T={\displaystyle \underset{-h(x)/2}{\overset{+h(x)/2}{\int }}{S}^{xx}{\left\{\begin{array}{ccc}1& z& z^2\end{array}\right\}}^T\left(1+\frac{z}{R(x)}\right)dz}$ (30)

${\left\{\begin{array}{ccc}{N}_{zx}& M_{zx}& Q_{zx}\end{array}\right\}}^T=K_S{\displaystyle \underset{-h(x)/2}{\overset{+h(x)/2}{\int }}{S}^{zx}{\left\{\begin{array}{ccc}1& z& z^2\end{array}\right\}}^T\left(1+\frac{z}{R(x)}\right)dz}$ (31)

In the last equation, $K_S$ is shear correction factor which is applying in the stress resultant derived from shear stresses because of preventing stress overestimation. We consider $K_S=5/6$ in the present study (Ghannad et al. 2013). Substituting strain invariants from Eq. (26) into Eqs. (23) and (25), considering $\delta {\Pi }_{EXT}^{}=\delta {\Pi }_{INT}^{}$ 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:

$\frac{d}{dx}\left[R(x)\left(N_x\left(1+u^{\prime }\right)+M_x{\phi }^{\prime }+N_{zx}\phi \right)\right]=0$ (32)

$\frac{d}{dx}\left[R(x)\left(M_x\left(1+u^{\prime }\right)+Q_x{\phi }^{\prime }+M_{zx}\phi \right)\right]-R(x)\left(N_z\phi +N_{zx}\left(1+u^{\prime }\right)+M_{zx}{\phi }^{\prime }\right)=0$ (33)

$\frac{d}{dx}\left[R(x)\left(N_x{w}^{\prime }+M_x{\psi }^{\prime }+N_{zx}\left(1+\psi \right)\right)\right]-N_{\theta }-\frac{1}{R(x)}\left(N_{\theta }{}^{*}w+M_{\theta }{}^{*}\psi \right)=-P_i(x)\left(R(x)-\frac{h(x)}{2}\right)$ (34)

$\begin{array}{l}\frac{d}{dx}\left[R(x)\left(M_x{w}^{\prime }+Q_x{\psi }^{\prime }+M_{zx}\left(1+\psi \right)\right)\right]-M_{\theta }-R(x)\left(N_z\left(1+\psi \right)+N_{zx}{w}^{\prime }+M_{zx}{\psi }^{\prime }\right)\\ -\frac{1}{R(x)}\left(M_{\theta }{}^{*}w+Q_{\theta }{}^{*}\psi \right)=P_i(x)\frac{h(x)}{2}\left(R(x)-\frac{h(x)}{2}\right)\end{array}$ (35)

3 ANALYTICAL SOLUTION

3.1 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):

$\overline{x}=\frac{x}{L},\overline{h}=\frac{h}{h_0},\overline{z}=\frac{z}{h_0},{\overline{r}}_i=\frac{r_i}{h_0},\overline{R}=\frac{R}{h_0},\overline{u}=\frac{u}{h_0},\overline{\phi }=\phi ,\overline{w}=\frac{w}{h_0},\overline{\psi }=\psi ,\epsilon =\frac{h_0}{L},{\overline{P}}_i=\frac{P_i{\overline{r}}_i{}^n}{\epsilon C_1}$ (36)

Dimensionless quantities of parameters are marked with over-bar $(\overline{})$ from now. Perturbation parameter$\epsilon$ should be normally considered as the ratio of minimal geometrical dimension of the structure $(h_0)$to maximal one$(L)$in order to be small quantity. $h_0$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 ${\epsilon }^i$ coefficients. Considering each coefficient result in displacement quasi-vector$\left\{\overline{y}\left(\overline{x}\right)\right\}$. 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 $\overline{x}$should be rewrite as follows:

$\frac{d}{dx}=\frac{1}{L}\frac{d}{d\overline{x}},\frac{d^2}{d{x}^2}=\frac{1}{L^2}\frac{d^2}{d{\overline{x}}^2}$ (37)

In shear deformation theory, the differentials and integrations are with respect to $\overline{x}$and $\overline{z}$, respectively. Therefore, for simplification and abbreviation of representing equations, a dimensionless integral could be defined:

$\overline{I}\overline{I}(i,j)=\frac{II(i,j)}{h_0{}^{i+j+1}}={\displaystyle \underset{-\overline{h}(\overline{x})/2}{\overset{+\overline{h}(\overline{x})/2}{\int }}{\overline{z}}^i{\left(\overline{R}(\overline{x})+\overline{z}\right)}^{j}d\overline{z}}$ (38)

which is the function of geometrical parameters $\overline{R}(\overline{x})$,$\overline{h}(\overline{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 $(c_0)$ should be either dimensionless as other parameters. Therefore, it is replaced with its dimensionless equivalent ${\overline{c}}_0=(c_0{\overline{r}}_i{}^n)/(\epsilon C_1)$ where $\epsilon$ is bookkeeping perturbed parameter. Furthermore, as there is no $\overline{u}$ in governing equations unlike $d\overline{u}/d\overline{x}$, we define $\overline{v}=\epsilon (d\overline{u}/d\overline{x})$as new unknown intermediate parameter. Therefore, one can calculate $\overline{u}$ indirectly from intermediate parameter $\overline{v}$by relation $\overline{u}=(1/\epsilon ){\displaystyle \int \overline{v}d\overline{x}+c_7}$where $c_7$ is integral constant. The intentional constants ${\overline{c}}_0$ and $c_7$will be calculated from boundary conditions. The following parameters need to be defined based on material constants of internal layer $C_1,C_2,k$and geometrical parameter of shear correction factor $K_s$because of abbreviation in representing inner and outer equations.

$\{\begin{array}{l}{C}_2/C_1=\overline{C_2},k/C_1=\overline{k},\left(2\left(C_1+C_2\right)\right)/C_1=\overline{C_{12}},\left(4\left(C_1+C_2\right)-k\right)/C_1=\overline{C{k}_1},\left(4\left(C_1+2C_2\right)-k\right)/C_1=\overline{C{k}_2},\\ \left(2\left(C_1+C_2\right)-3k\right)/C_1=\overline{C{k}_3},\left(2\left(C_1+C_2\right)-k\right)/C_1=\overline{C{k}_4},\left(4\left(C_1+C_2\right)-5k\right)/C_1=\overline{C{k}_5},\left(4\left(C_1+C_2\right)-6k\right)/C_1=\overline{C{k}_6},\\ \left(8\left(C_1+2C_2\right)-7k\right)/C_1=\overline{C{k}_7},\left(2\left(2C_1+3C_2\right)-k\right)/C_1=\overline{C{k}_8},\left(2\left(2C_1+5C_2\right)-6k\right)/C_1=\overline{C{k}_9},\left(1-K_s\right)={\widehat{K}}_s\end{array}$ (39)

3.2 Outer expansion

The outer expansion of solution is considered a uniform series of $\epsilon$ as ${\overline{y}}_O\left(\overline{x},\epsilon \right)=\epsilon \left({\overline{y}}_{O1}\left(\overline{x}\right)+\epsilon {\overline{y}}_{O2}\left(\overline{x}\right)\right)$. The subscript “O” 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 $\epsilon$, result in the first and second order equations of outer solution. In this section${\left(\right)}^{\prime }=d\left(\right)/d\overline{x}$.

$O({\epsilon }^1):\left[A_O\right]\left\{{\overline{y}}_{O1}\right\}=\left\{F_{O1}\right\},O({\epsilon }^2):\left[A_O\right]\left\{{\overline{y}}_{O2}\right\}=\left\{F_{O2}^{}\right\},\left\{F_{O2}^{}\right\}=\left\{F_{O2}^{II}\right\}+\left\{F_{O2}^{I{I}^{\prime }}\right\}$ (40)

where $\left[A_O\right],\left\{F_{O1}\right\}$and $\left\{F_{O2}^{}\right\}$are coefficient matrices, non-homogeneity vectors of first and second order equation, respectively. $\left\{{\overline{y}}_{Oi}\right\}$ are unknown displacement vectors in “i”th-order of outer solution. $\left\{F_{O2}^{}\right\}$ consist of two vectors $\left\{F_{O2}^{II}\right\}$and $\left\{F_{O2}^{I{I}^{\prime }}\right\}$ correspond to $\overline{I}\overline{I}(\overline{x})$ and derivative of $\overline{I}\overline{I}(\overline{x})$, respectively. These vectors would be defined in the ^{Appendix A}. Other non-zero components of the mentioned matrix and vectors are

$\left\{{\overline{y}}_{O1}\right\}={\left\{{\overline{v}}_{O1},{\overline{\phi }}_{O1},{\overline{w}}_{O1},{\overline{\psi }}_{O1}\right\}}^T,\left\{{\overline{y}}_{O2}\right\}={\left\{{\overline{v}}_{O2},{\overline{\phi }}_{O2},{\overline{w}}_{O2},{\overline{\psi }}_{O2}\right\}}^T$ (41)

$\{\begin{array}{l}{\left[A_O\right]}_{11}=\overline{k}\overline{I}\overline{I}(\mathrm{0,}n+1),{\left[A_O\right]}_{22}=-K_s\overline{C_{12}}\overline{I}\overline{I}(\mathrm{0,}n+1),{\left[A_O\right]}_{13}=-{\left[A_O\right]}_{31}=-\overline{C{k}_1}\overline{I}\overline{I}(\mathrm{0,}n)\\ {\left[A_O\right]}_{14}=-{\left[A_O\right]}_{41}=-\overline{C{k}_1}\left(\overline{I}\overline{I}(\mathrm{1,}n)+\overline{I}\overline{I}(\mathrm{0,}n+1)\right),{\left[A_O\right]}_{33}=-\overline{k}\overline{I}\overline{I}(\mathrm{0,}n-1)\\ {\left[A_O\right]}_{34}={\left[A_O\right]}_{43}=\overline{C{k}_1}\overline{I}\overline{I}(\mathrm{0,}n)-\overline{k}\overline{I}\overline{I}(\mathrm{1,}n-1),{\left[A_O\right]}_{44}=\overline{C{k}_1}\overline{I}\overline{I}(\mathrm{1,}n)-\overline{k}\left(\overline{I}\overline{I}(\mathrm{0,}n+1)+\overline{I}\overline{I}(\mathrm{2,}n-1)\right)\end{array}$ (42)

${\left\{F_{O1}\right\}}_1={\overline{c}}_0,{\left\{F_{O1}\right\}}_2=\mathrm{0,}{\left\{F_{O1}\right\}}_3=-{\overline{P}}_i(\overline{x})\left(\overline{R}(\overline{x})-\frac{\overline{h}(\overline{x})}{2}\right),{\left\{F_{O1}\right\}}_4={\overline{P}}_i(\overline{x})\left(\frac{\overline{R}(\overline{x})\overline{h}(\overline{x})}{2}-\frac{{\overline{h}}^2(\overline{x})}{4}\right)$ (43)

The solutions of the algebraic equations (40) are as follows:

$\left\{{\overline{y}}_{O1}\right\}={\left[A_O\right]}^{-1}\left\{F_{O1}\right\},\left\{{\overline{y}}_{O2}\right\}={\left[A_O\right]}^{-1}\left\{F_{O2}^{}\right\}$ (44)

3.3 Inner expansion

As the outer solutions don’t satisfy the B.C., we can conclude existence of boundary layers at $\overline{x}=\mathrm{0,1}$. Therefore, fast variable $({\tilde{x}}_{\alpha })$ 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.

$\overline{x}=0\to \alpha =\mathrm{0,}{\tilde{x}}_0=\frac{\overline{x}}{\epsilon }(leftboundary),\overline{x}=1\to \alpha =\mathrm{1,}{\tilde{x}}_1=\frac{\left(\overline{x}-1\right)}{\epsilon }(rightboundary)$ (45)

Subscript “$\alpha =\mathrm{0,1}$” shows $\overline{x}=\alpha$ in the boundary. Perturbation parameter appears in new definitions of differential based on fast variables${\tilde{x}}_{\alpha }$:

$\frac{d}{d{\tilde{x}}_{\alpha }}=\epsilon \frac{d}{d\overline{x}},\frac{d^2}{d{\tilde{x}}_{\alpha }{}^2}={\epsilon }^2\frac{d^2}{d{\overline{x}}^2}$ (46)

In cylinder with variable thickness and non-uniform pressure, it is necessary to derive Taylor expansion for all the parameters of axial function $\overline{\Lambda }(\tilde{x})$as:

$\overline{\Lambda }(\tilde{x})={\overline{\Lambda }}_{\alpha }+\epsilon {\tilde{x}}_{\alpha }D{\overline{\Lambda }}_{\alpha }+\mathrm{...}$ (47)

Taylor expansion with ${\overline{\Lambda }}_{\alpha }=\overline{\Lambda }(\overline{x}=\alpha )$and $D{\overline{\Lambda }}_{\alpha }={d\overline{\Lambda }(\overline{x})/d\overline{x}|}_{\overline{x}=\alpha }$definition substitutes $\overline{x}$ by new fast variable ${\tilde{x}}_{\alpha }$ of inner expansion. Taylor expansion should be written for the following parameters:

$\{\begin{array}{l}\overline{h}(\tilde{x})={\overline{h}}_{\alpha }+\epsilon {\tilde{x}}_{\alpha }D{\overline{h}}_{\alpha },\overline{R}(\tilde{x})={\overline{R}}_{\alpha }+\epsilon {\tilde{x}}_{\alpha }D{\overline{R}}_{\alpha },\\ {\overline{P}}_i(\tilde{x})={\overline{P}}_i{}_{\alpha }+\epsilon {\tilde{x}}_{\alpha }D{\overline{P}}_i{}_{\alpha },\overline{I}{\overline{I}}_{\alpha }(\tilde{x})=\overline{I}{\overline{I}}_{\alpha }+\epsilon {\tilde{x}}_{\alpha }D\overline{I}{\overline{I}}_{\alpha }\end{array}\left(\alpha =\mathrm{0,1}\right)$ (48)

Substitution of inner expansion ${\overline{y}}_{\alpha }\left({\tilde{x}}_{\alpha },\epsilon \right)=\epsilon \left({\overline{y}}_{\alpha 1}\left({\tilde{x}}_{\alpha }\right)+\epsilon {\overline{y}}_{\alpha 2}\left({\tilde{x}}_{\alpha }\right)\right)$ in governing equations with mentioned changes in section 3.1 and considering terms with the same order of $\epsilon$ result in inner equations at boundary$\alpha$:

$\{\begin{array}{l}O({\epsilon }^1):\left[A_{\alpha 1}\right]\frac{d^2}{d{\tilde{x}}_{\alpha }{}^2}\left\{{\overline{y}}_{\alpha 1}\right\}+\left[A_{\alpha 2}\right]\frac{d}{d{\tilde{x}}_{\alpha }}\left\{{\overline{y}}_{\alpha 1}\right\}+\left[A_{\alpha 3}\right]\left\{{\overline{y}}_{\alpha 1}\right\}=\left\{F_{\alpha 1}\right\}\\ O({\epsilon }^2):\left[A_{\alpha 1}\right]\frac{d^2}{d{\tilde{x}}_{\alpha }{}^2}\left\{{\overline{y}}_{\alpha 2}\right\}+\left[A_{\alpha 2}\right]\frac{d}{d{\tilde{x}}_{\alpha }}\left\{{\overline{y}}_{\alpha 2}\right\}+\left[A_{\alpha 3}\right]\left\{{\overline{y}}_{\alpha 2}\right\}=\left\{F_{\alpha 2}^{}\right\}\\ \left\{F_{\alpha 2}^{}\right\}=\left\{F_{\alpha 2}^{I{I}_{\alpha }}\right\}+\left\{F_{\alpha 2}^{D{A}_{\alpha }}\right\}+\left\{F_{\alpha 2}^{DI{I}_{\alpha }}\right\}+\left\{F_{\alpha 2}^{P_{\alpha }}\right\}+\left\{F_{\alpha 2}^{D{P}_{\alpha }}\right\}\end{array}$ (49)

$\left[A_{\alpha 1}\right],\left[A_{\alpha 2}\right]$and $\left[A_{\alpha 3}\right]$ are coefficient matrices at the boundary $\alpha$. $\left\{F_{\alpha 1}^{}\right\}$and $\left\{F_{\alpha 2}^{}\right\}$are non-homogeneity vectors of differential equation at the first and second order equation for the boundary $\alpha$, respectively. Unknown displacement vectors in “i”th-order of inner solution at the boundary $\alpha$ are $\left\{{\overline{y}}_{\alpha i}\right\}$.$\left\{F_{\alpha 2}^{I{I}_{\alpha }}\right\}$and $\left\{F_{\alpha 2}^{DI{I}_{\alpha }}\right\}$, in$\left\{F_{\alpha 2}^{}\right\}$, are correspond with $\overline{I}{\overline{I}}_{\alpha }(\tilde{x})$ and $D\overline{I}{\overline{I}}_{\alpha }(\tilde{x})$ (in Taylor expansion), respectively. $\left\{F_{\alpha 2}^{D{A}_{\alpha }}\right\}$derive in non-homogeneity of second order equation which resulting from Tayler expansion of coefficient matrices in first order equation because of variable thickness. $\left\{F_{\alpha 2}^{P_{\alpha }}\right\}$and $\left\{F_{\alpha 2}^{D{P}_{\alpha }}\right\}$, include ${\overline{P}}_i{}_{\alpha }({\tilde{x}}_{\alpha })$ and $D{\overline{P}}_i{}_{\alpha }({\tilde{x}}_{\alpha })$(in Taylor expansion), are resulting from variable thickness and non-uniform pressure, respectively. Considering too many terms, $\left\{F_{\alpha 2}^{}\right\}$ vector would be defined in the ^{Appendix A}. Other non-zero components of the matrices and vectors defined below. In this section${\left(\right)}^{\prime }=d\left(\right)/d{\tilde{x}}_{\alpha }$ at the boundary $\alpha$.

$\left\{{\overline{y}}_{\alpha 1}\right\}={\left\{{\overline{v}}_{\alpha 1},{\overline{\phi }}_{\alpha 1},{\overline{w}}_{\alpha 1},{\overline{\psi }}_{\alpha 1}\right\}}^T,\left\{{\overline{y}}_{\alpha 2}\right\}={\left\{{\overline{v}}_{\alpha 2},{\overline{\phi }}_{\alpha 2},{\overline{w}}_{\alpha 2},{\overline{\psi }}_{\alpha 2}\right\}}^T$ (50)

$\{\begin{array}{l}{\left[A_{\alpha 1}\right]}_{22}=\overline{k}\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1),{\left[A_{\alpha 1}\right]}_{33}=K_s\overline{C_{12}}\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)\\ {\left[A_{\alpha 1}\right]}_{34}={\left[A_{\alpha 1}\right]}_{43}=K_s\overline{C_{12}}\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1),{\left[A_{\alpha 1}\right]}_{44}=K_s\overline{C_{12}}\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1)\end{array}$ (51)

$\{\begin{array}{l}{\left[A_{\alpha 2}\right]}_{12}={\left[A_{\alpha 2}\right]}_{21}=\overline{k}\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1),{\left[A_{\alpha 2}\right]}_{23}=-{\left[A_{\alpha 2}\right]}_{32}=-K_s\overline{C_{12}}\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)-\overline{C{k}_1}\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n)\\ {\left[A_{\alpha 2}\right]}_{24}=-{\left[A_{\alpha 2}\right]}_{42}=-K_s\overline{C_{12}}\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1)-\overline{C{k}_1}\left(\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n)+\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1)\right)\end{array}$ (52)

$\{\begin{array}{l}{\left[A_{\alpha 3}\right]}_{11}=\overline{k}\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1),{\left[A_{\alpha 3}\right]}_{22}=-K_s\overline{C_{12}}\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1),{\left[A_{\alpha 3}\right]}_{13}=-{\left[A_{\alpha 3}\right]}_{31}=-\overline{C{k}_1}\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n)\\ {\left[A_{\alpha 3}\right]}_{14}=-{\left[A_{\alpha 3}\right]}_{41}=-\overline{C{k}_1}\left(\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)+\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n)\right),{\left[A_{\alpha 3}\right]}_{33}=-\overline{k}\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n-1)\\ {\left[A_{\alpha 3}\right]}_{34}={\left[A_{\alpha 3}\right]}_{43}=\overline{C{k}_1}\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n)-\overline{k}\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n-1),{\left[A_{\alpha 3}\right]}_{44}=2\overline{C{k}_1}\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n)-\overline{k}\left(\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)+\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n-1)\right)\end{array}$ (53)

${\left\{F_{\alpha 1}\right\}}_1={\overline{c}}_0,{\left\{F_{\alpha 1}\right\}}_2=\mathrm{0,}{\left\{F_{\alpha 1}\right\}}_3=-{\overline{P}}_{i\alpha }\left({\overline{R}}_{\alpha }-\frac{{\overline{h}}_{\alpha }}{2}\right),{\left\{F_{\alpha 1}\right\}}_4={\overline{P}}_{i\alpha }\left(\frac{{\overline{R}}_{\alpha }{\overline{h}}_{\alpha }}{2}-\frac{{\overline{h}}_{\alpha }^2}{4}\right)$ (54)

Eqs. (49) are systems of coupled non-homogenous differential equations with constant coefficients. Each equation has general $\left({\left\{\right\}}_{gen.}\right)$and particular $\left({\left\{\right\}}_{par.}\right)$solution:

$\left\{{\overline{y}}_{\alpha 1}\right\}={\left\{{\overline{y}}_{\alpha 1}\right\}}_{gen.}+{\left\{{\overline{y}}_{\alpha 1}\right\}}_{par.},\left\{{\overline{y}}_{\alpha 2}\right\}={\left\{{\overline{y}}_{\alpha 2}\right\}}_{gen.}+{\left\{{\overline{y}}_{\alpha 2}\right\}}_{par.}$ (55)

Considering $m_{\alpha }$and $\left\{V_{\alpha }\right\}$ as eigenvalues and eigenvectors, respectively, and substituting general solution of exponential form ${\left\{{\overline{y}}_{\alpha }\right\}}_{gen.}=\left\{V_{\alpha }\right\}\exp (m_{\alpha }{\tilde{x}}_{\alpha })$in homogenous part of Eqs. (49) lead in an eigenvalue problem:

$\left(\left[A_{\alpha 1}\right]m_{\alpha }{}^2+\left[A_{\alpha 2}\right]m_{\alpha }+\left[A_{\alpha 3}\right]\right)\left\{V_{\alpha }\right\}=\left\{0\right\}$ (56)

The necessary condition for existing the solution of Eq. (56) is zero value of the coefficient determinant which is the characteristic equation of the system. Six non-zero roots of it are the eigenvalues$(m_{\alpha }{}_i)$. Substituting roots in Eq. (56) lead in corresponding eigenvectors$(V_{\alpha }{}_i)$. The eigenvalues and eigenvectors are complex conjugate. Considering Van-Dyke’s matching principle (Nayfeh 1981), the solution should be finite at ${\tilde{x}}_{\alpha }\to \infty$. Therefore, in left boundary$(\alpha =0)$ eigenvalues with positive real part and in right boundary $(\alpha =1)$ eigenvalues with negative real part are omitted. The general solution of the boundary $\alpha$ could be calculated.

${\left\{{\overline{y}}_{\alpha 1}\right\}}_{gen.}={\left\{{\overline{y}}_{\alpha 2}\right\}}_{gen.}={\left\{{\overline{y}}_{\alpha }\right\}}_{gen.}={\displaystyle \sum _{i=1}^3{c}_{\alpha }{}_i}\left\{V_{\alpha }{}_i\right\}\exp (m_{\alpha }{}_i{\tilde{x}}_{\alpha })$ (57)

where $c_{\alpha }{}_i(i=\mathrm{1,2,3})$are three constant at each boundary which could be calculated by boundary conditions. The particular solution of first order equation in Eq. (49) is simply calculated by ${\left\{{\overline{y}}_{\alpha 1}\right\}}_{par.}={\left[A_{\alpha 3}\right]}^{-1}\left\{F_{\alpha 1}\right\}$. As $\left\{F_{\alpha 2}^{}\right\}$is consist of nonlinear polynomial terms$\left({\left\{\right\}}^{pol.}\right)$ and exponential terms by the equal roots $(\exp (m_{\alpha }{}_i{\tilde{x}}_{\alpha }))$and unequal roots $(\exp (q_{\alpha }{}_j{\tilde{x}}_{\alpha }))$with characteristic equations based on $O({\epsilon }^1)$solution, the particular solution of $O({\epsilon }^2)$is calculated by undetermined coefficients method as follows:

$\{\begin{array}{l}{\left\{{\overline{y}}_{\alpha 2}\right\}}_{par.}={\left\{{\overline{y}}_{\alpha 2}\right\}}_{par.}^{pol.}+{\left\{{\overline{y}}_{\alpha 2}\right\}}_{par.}^{\exp (mi)}+{\left\{{\overline{y}}_{\alpha 2}\right\}}_{par.}^{\exp (qj)},{\left\{{\overline{y}}_{\alpha 2}\right\}}_{par.}^{pol.}=\left\{B_2\right\}{\tilde{x}}^2+\left\{B_1\right\}\tilde{x}+\left\{B_0\right\}\\ {\left\{{\overline{y}}_{\alpha 2}\right\}}_{par.}^{\exp (mi)}={\displaystyle {\sum }_i\left({\left\{B_2\right\}}_{mi}{\tilde{x}}^2+{\left\{B_1\right\}}_{mi}\tilde{x}+{\left\{B_0\right\}}_{mi}\right)\exp (m_{\alpha }{}_i{\tilde{x}}_{\alpha })}\\ {\left\{{\overline{y}}_{\alpha 2}\right\}}_{par.}^{\exp (qj)}={\displaystyle {\sum }_j\left({\left\{B_2\right\}}_{qj}{\tilde{x}}^2+{\left\{B_1\right\}}_{qj}\tilde{x}+{\left\{B_0\right\}}_{qj}\right)\exp (q_{\alpha }{}_j{\tilde{x}}_{\alpha })}\end{array}$ (58)

where subscript “i” and “j” shows the number of equal and unequal roots with characteristic equations. Substituting ${\left\{{\overline{y}}_{\alpha 2}\right\}}_{par.}$in Eq. (49) -$O({\epsilon }^2)$ lead in undetermined coefficients $\left\{B_1\right\}$,$\left\{B_1\right\}$and $\left\{B_2\right\}$.

3.4 Composite solution

In the MAE method, the composite solution $\left(\left\{{\overline{y}}_{comp.}\right\}\right)$is the summation of three calculated expansions (one outer

$\left\{{\overline{y}}_O\right\}$and two inner$\left\{{\overline{y}}_{\alpha =0}\right\},\left\{{\overline{y}}_{\alpha =1}\right\}$) minus the overlapped parts of them. Outer solution at $\overline{x}\to \mathrm{0,1}$ and inner solutions at ${\tilde{x}}_{\alpha }\to \pm \infty$are overlapped and these common parts have to be removed from composite solution. Therefore,

$\left\{{\overline{y}}_{comp.}\right\}=\left\{{\overline{y}}_O\right\}+\left\{{\overline{y}}_{\alpha =0}\right\}+\left\{{\overline{y}}_{\alpha =1}\right\}-\left\{{\overline{y}}_O^{\alpha =0}\right\}-\left\{{\overline{y}}_O^{\alpha =1}\right\}$ (59)

where $\left\{{\overline{y}}_O^{\alpha =0}\right\}$and $\left\{{\overline{y}}_O^{\alpha =1}\right\}$are common parts of inner and outer solutions at two ends of the shell which can be determined by definition of intermediate variable or Van-Dyke’s matching principle (Nayfeh 1981). Eight constants, consist of three constants in general solution of each boundary and two constants ${\overline{c}}_0$and$c_7$, should be calculated by the boundary conditions. The clamped boundary conditions in “i”th-order perturbation solution are:

$\overline{x}=\mathrm{0,1}\left({\tilde{x}}_{\mathrm{0,1}}=0\right)\to {\overline{u}}_{\alpha i},{\overline{\phi }}_{\alpha i},{\overline{w}}_{\alpha i},{\overline{\psi }}_{\alpha i}=0(i=\mathrm{1,2})$ (60)

Finally, the unknown vector$\left\{\overline{y}\right\}=\left\{{\overline{y}}_{comp.}\right\}=\left\{\overline{u},\overline{\phi },\overline{w},\overline{\psi }\right\}$ which consists of dimensionless displacement field components would be obtained in terms of $\overline{x}$and $\overline{z}$variables. Considering Eq. (4) and ${\overline{U}}_{z,x}=U_{z,x}/h_0$, the dimensionless radial and axial displacements can be calculated. Using Eqs. (5-9) would yield $\left[F\right],\left[C\right],I_{\mathrm{1,2,3}},J$and $\left[E\right]$, respectively. The hydrostatic pressure, strain energy function and second Piola-Kirchhoff stress distribution could be calculated by using Eqs. (14), (18) and (19). The relation $\left[\sigma \right]=\left(1/J\right)\left[F\right]\left[S\right]{\left[F\right]}^T$ would result in Cauchy stress components. The analytical solution has been carried out by writing the program in MAPLE 18 software.

4 RESULTS AND DISCUSSION

4.1 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 $C_{10},C_{01}$and $d$should be defined for ANSYS software. $C_{10}$and $C_{01}$corresponding to $C_{1n}$and $C_{2n}$could be calculated from Eq. (12) in each layer. $d$is material incompressibility parameter by relation with bulk modulus as $k_n.d=2$ (Montella et al. 2014, Gao et al. 2009) which could be calculated from Eq. (16) in each layer. For non-uniform 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.

4.2 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: $r_i=47mm$and $L=400mm$. 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 $C_1$ and $C_2$in Table 1 and Eq. (15) result the variation range of $k$about $1-100MPa$ in nearly incompressible limit. Important remark is that for $k<2MPa$, the eigenvalues of characteristic equation no longer have conjugate complex form and MAE solution diverge. Hence, in Table 1, we set $k=10MPa$ 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 $5kPa-13kPa$. Dimensionless Cauchy stresses and hydrostatic pressure are defined as:

$\overline{\sigma }=\frac{\sigma }{P_{i0}},\overline{P}=\frac{P}{P_{i0}}$ (61)

$P_{i0}$is considered $9kPa$which is the average of maximum and minimum pressure distributions. The variations of thickness profiles are in the range of $6mm-12mm$. Tables 2 and 3 show the characteristic of applied non-uniform pressure and thickness profiles, respectively.

Table 1:
The characteristic of material constants.

Material Constants ID	Reference	$C_1(MPa)$	$C_2(MPa)$	$k(MPa)$	Material Sample
$M{C}_1$	Selvadurai and Shi (2012)	0.242	0.142	10	Natural gum rubbers
$M{C}_2$	Montella et al. (2014)	0.251	0.263	36	TDM 600
$M{C}_3$	Sussman and Bathe (1987)	0.177	0.045	10	Rubber bushes
$M{C}_4$	Boyce and Arruda (2000)	0.16	0.015	10	Rubbers
$M{C}_5$	Batra and Bahrami (2009)	0.1858	-0.01935	10	Rubberlike materials
$M{C}_6$	Dias et al. (2014)	0.423	-0.262	10.5	Silicone elastomers
$M{C}_7$	Vossoughi and Tozeren (1998)	0.09	0.118	10	polyurethane rubbers
$M{C}_8$	Kiendl et al. (2015)	0.1848	0.0264	10	Rubbers
$M{C}_9$	Tanveer and Zu (2012)	0.552	0.138	10	Rubber seals

Table 2:
The characteristic of non-uniform internal pressure profiles.

Pressure Profile ID	Pressure Profile Relation	Load Constants
$P_{i0}$	$P_i(\overline{x})=P_{i0}=cons\tan t$	$P_{i0}=9kPa$
$P_{i1}$	$P_i(\overline{x})=P_a-\left(P_a-P_b\right)\overline{x}$	$P_a=13kPa,P_b=5kPa$
$P_{i2}$	$P_i(\overline{x})=P_a-\left(P_a-P_b\right){\overline{x}}^2$	$P_a=13kPa,P_b=5kPa$
$P_{i3}$	$P_i(\overline{x})=P_a-\left(P_a-P_b\right){\overline{x}}^3$	$P_a=13kPa,P_b=5kPa$
$P_{i4}$	$P_i(\overline{x})=P_a-4\left(P_a-P_{ab}\right)\left(\overline{x}-{\overline{x}}^2\right)$	$P_a=P_b=5kPa,P_{ab}=13kPa$
$P_{i5}$	$P_i(\overline{x})=P_a-\left(P_a-P_{ab}\right)\sin \left(\pi \overline{x}\right)$	$P_a=P_b=13kPa,P_{ab}=5kPa$

Table 3:
The characteristic of non-uniform thickness profiles.

Thickness Profile ID	Thickness Profile Relation	Geometry Constants
$h_0$	$h(\overline{x})=h_0=cons\tan t$	$h_0=6mm$
$h_1$	$h(\overline{x})=h_a-\left(h_a-h_b\right)\overline{x}$	$h_a=12mm,h_b=6mm$
$h_2$	$h(\overline{x})=h_a-\left(h_a-h_b\right){\overline{x}}^2$	$h_a=12mm,h_b=6mm$
$h_3$	$h(\overline{x})=h_a-\left(h_a-h_b\right){\overline{x}}^3$	$h_a=12mm,h_b=6mm$
$h_4$	$h(\overline{x})=h_a-4\left(h_a-h_{ab}\right)\left(\overline{x}-{\overline{x}}^2\right)$	$h_a=h_b=6mm,h_{ab}=12mm$
$h_5$	$h(\overline{x})=h_a-\left(h_a-h_{ab}\right)\sin \left(\pi \overline{x}\right)$	$h_a=h_b=12mm,h_{ab}=6mm$

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 and Gasser 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 $5kPa(40mmHg)-13kPa(100mmHg)$ which are the mean blood vessel pressures of human soft tissues. $100mmHg$is the mean of systolic/diastolic pressure and $40mmHg$may be occurred in hypotension pressure of arteries (Abdessamad et al. 2018, Humphrey and O’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 non-uniform 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).

4.3 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 define the validity range of the current solution. Material constants, applied pressures and geometry of the cylinder are considered in the range of $5<\overline{R}<20$ and $1/200<P_i/\left(C_1+C_2\right)<1/50$ (Batra and Bahrami 2009, Humphrey and O’Rourke 2015, Azar et al. 2018). The difference percentage of dimensionless radial and axial displacement resulting from the numerical and analytical solution i.e. $Diff{\overline{U}}_{z,x}(\%)=\left|\left({\overline{U}}_{z,x}{}^{MAE}-{\overline{U}}_{z,x}{}^{FEM}\right)/{\overline{U}}_{z,x}{}^{FEM}\right|\times 100$ for mentioned $\overline{R}$and $P_i/(C_1+C_2)$ are less than 8%. $\overline{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 $\overline{R}$because of intensifying nonlinear behavior of the cylinder while for small $\overline{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 $\overline{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 $\overline{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 $P_i/(C_1+C_2)$and $k$or decrease in $C_1+C_2$ 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 $\left(\overline{x}=0.5\right)$, 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.

Table 4:
Analytical and numerical values of maximum radial displacements for different material constants.

Mat. ID	${\overline{U}}_z{,}_{\max }$	$P_i/(C_1+C_2)=1/50$				$P_i/(C_1+C_2)=1/100$				$P_i/(C_1+C_2)=1/200$
		$\overline{R}=5$	$\overline{R}=10$	$\overline{R}=20$		$\overline{R}=5$	$\overline{R}=10$	$\overline{R}=20$		$\overline{R}=5$	$\overline{R}=10$	$\overline{R}=20$
$M{C}_1$	FSDT	0.0734	0.2991	1.2127		0.0363	0.1481	0.5943		0.0182	0.0736	0.2952
	FEM	0.0759	0.3047	1.2699		0.0370	0.1488	0.5970		0.0184	0.0737	0.2930
$M{C}_2$	FSDT	0.0691	0.2859	1.1754		0.0339	0.1410	0.5760		0.0169	0.0700	0.2812
	FEM	0.0736	0.2913	1.2210		0.0359	0.1407	0.5746		0.0178	0.0706	0.2780
$M{C}_3$	FSDT	0.0701	0.2917	1.1863		0.0347	0.1437	0.5810		0.0173	0.0714	0.2876
	FEM	0.0730	0.2957	1.2319		0.0358	0.1442	0.5812		0.0180	0.0717	0.2847
$M{C}_4$	FSDT	0.0685	0.2876	1.1771		0.0340	0.1421	0.5763		0.0169	0.0706	0.2861
	FEM	0.0726	0.2921	1.2234		0.0356	0.1430	0.5766		0.0179	0.0712	0.2830
$M{C}_5$	FSDT	0.0681	0.2869	1.1762		0.0338	0.1417	0.5760		0.0168	0.0704	0.2855
	FEM	0.0725	0.2930	1.2214		0.0349	0.1419	0.5742		0.0177	0.0710	0.2824
$M{C}_6$	FSDT	0.0676	0.2832	1.1371		0.0336	0.1408	0.5682		0.0167	0.0702	0.2836
	FEM	0.0722	0.2902	1.1780		0.0362	0.1421	0.5683		0.0178	0.0708	0.2812
$M{C}_7$	FSDT	0.0696	0.2902	1.1871		0.0345	0.1433	0.5810		0.0172	0.0712	0.2875
	FEM	0.0743	0.2969	1.2375		0.0360	0.1440	0.5813		0.0181	0.0716	0.2850
$M{C}_8$	FSDT	0.0697	0.2899	1.1861		0.0346	0.1425	0.5798		0.0172	0.0712	0.2875
	FEM	0.0742	0.2970	1.2360		0.0365	0.1440	0.5799		0.0180	0.0716	0.2846

Table 5:
Analytical and numerical values of maximum axial displacements for different material constants.

Mat. ID	${\overline{U}}_x{,}_{\max }$	$P_i/(C_1+C_2)=1/50$				$P_i/(C_1+C_2)=1/100$				$P_i/(C_1+C_2)=1/200$
		$\overline{R}=5$	$\overline{R}=10$	$\overline{R}=20$		$\overline{R}=5$	$\overline{R}=10$	$\overline{R}=20$		$\overline{R}=5$	$\overline{R}=10$	$\overline{R}=20$
$M{C}_1$	FSDT	0.0147	0.0453	0.1096		0.0073	0.0231	0.0623		0.0037	0.0117	0.0331
	FEM	0.0151	0.0469	0.1089		0.0076	0.0242	0.0621		0.0038	0.0125	0.0332
$M{C}_2$	FSDT	0.0138	0.0453	0.1162		0.0069	0.0227	0.0634		0.0035	0.0113	0.0332
	FEM	0.0149	0.0460	0.1080		0.0073	0.0239	0.0611		0.0037	0.0120	0.0333
$M{C}_3$	FSDT	0.0138	0.0445	0.1139		0.0069	0.0226	0.0628		0.0034	0.0114	0.0329
	FEM	0.0149	0.0452	0.1068		0.0075	0.0239	0.0611		0.0037	0.0121	0.0331
$M{C}_4$	FSDT	0.0136	0.0445	0.1148		0.0068	0.0225	0.0630		0.0034	0.0113	0.0329
	FEM	0.0148	0.0457	0.1069		0.0075	0.0237	0.0603		0.0037	0.0122	0.0332
$M{C}_5$	FSDT	0.0135	0.0441	0.1135		0.0068	0.0224	0.0626		0.0034	0.0113	0.0327
	FEM	0.0148	0.0456	0.1050		0.0073	0.0239	0.0602		0.0037	0.0121	0.0329
$M{C}_6$	FSDT	0.0134	0.0408	0.0960		0.0066	0.0216	0.0584		0.0034	0.0111	0.0316
	FEM	0.0145	0.0432	0.0973		0.0073	0.0224	0.0587		0.0037	0.0121	0.0318
$M{C}_7$	FSDT	0.0139	0.0454	0.1166		0.0069	0.0228	0.0640		0.0034	0.0114	0.0332
	FEM	0.0150	0.0466	0.1085		0.0075	0.0232	0.0613		0.0038	0.0123	0.0335
$M{C}_8$	FSDT	0.0137	0.0444	0.1134		0.0068	0.0226	0.0628		0.0034	0.0114	0.0329
	FEM	0.0149	0.0460	0.1063		0.0075	0.0238	0.0610		0.0037	0.0123	0.0328

Figure 2:
(a) Longitudinal strain distribution along axial direction at different layers ($P_i(\overline{x})=1.5P_{i0},h(\overline{x})=h_0,M{C}_5$) and (b) circumferential strain distribution along thickness for different pressure ($\overline{x}=\mathrm{0.5,}h(\overline{x})=h_0,M{C}_4$).

In order to show the displacements and stresses distribution in FG hyperelastic cylinder with variable thickness under non-uniform internal pressure, linear increase of thickness counter $\overline{x}$direction is considered proportional to linear pressure profile with positive inhomogeneity index $(n=2)$. The constants of FG Mooney-Rivlin model for rubber are considered as $M{C}_9$. Figure 3 illustrate displacement contour resulted from FE modeling in 3/4 expansion of the shell section. The distribution of dimensionless displacements resulted from FSDT and FEM are plotted in Figure 4 at different layers. Boundary conditions and radial direction of applying pressure cause ${\overline{U}}_x<{\overline{U}}_z$ 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 $\overline{z}\le 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 $0.993<J<1.010$. Dimensionless radial Cauchy stress distribution in middle layer resulted from first and second order MAE solution are depicted in Figure 6. $O({\epsilon }^1)$ solution is suitable for linear analysis while $O({\epsilon }^2)$ 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.

Figure 3:
Contour plot of (a) radial and (b) axial displacement resulted from FEM in – expansion of shell section

Figure 4:
Dimensionless (a) radial and (b) axial displacement distribution along axial direction at different layers

Figure 5:
Dimensionless (a) circumferential (b) axial (c) shear Cauchy stress and (d) hydrostatic pressure distribution along axial direction at different layers ($n=\mathrm{2,}{P}_i(\overline{x})=P_{i1},h(\overline{x})=h_1,M{C}_9$).

Figure 6:
Dimensionless radial Cauchy stress distribution at middle layer along axial direction for different MAE order.

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$M{C}_9$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\le n\le +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 non-homogeneous 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 7:
The effect of inhomogeneity index on dimensionless (a) radial displacement (b) axial displacement (c) hydrostatic pressure (along axial direction at internal layer) and (d) material properties (along radial direction) ($P_i(\overline{x})=P_{i1},h(\overline{x})=h_1$).

Table 6:
Maximum displacements and hydrostatic pressure of internal, middle and external layer for different inhomogeneity index ($P_i(\overline{x})=P_{i1},h(\overline{x})=h_1$).

	${\overline{U}}_{z}at\overline{x}=0.15$				${\overline{U}}_{x}at\overline{x}=0.65$				$\overline{P}at\overline{x}=0$
	$\overline{z}=-0.5$	$\overline{z}=0$	$\overline{z}=+0.5$		$\overline{z}=-0.5$	$\overline{z}=0$	$\overline{z}=+0.5$		$\overline{z}=-0.5$	$\overline{z}=0$	$\overline{z}=+0.5$
$n=+4$	0.0744	0.0697	0.0650		-0.0224	-0.0222	-0.0219		9.425	2.064	-5.296
$n=+2$	0.0932	0.0873	0.0813		-0.0271	-0.0267	-0.0264		11.177	1.965	-7.247
$n=0$	0.1152	0.1078	0.1004		-0.0325	-0.0320	-0.0315		13.105	1.631	-9.843
$n=-2$	0.1403	0.1312	0.1221		-0.0387	-0.0380	-0.0373		15.201	0.970	-13.262
$n=-4$	0.1685	0.1575	0.1465		-0.0457	-0.0447	-0.0438		17.459	-0.128	-17.715

4.4 Effect of pressure and geometry profiles

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 $M{C}_9$. 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 $P_{i1}$,$P_{i2}$ and $P_{i3}$ 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 $P_{i1}$and thickness $h_1$ 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 $P_{i4}$ and $P_{i5}$ 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 $P_{i4}$ and $P_{i5}$ reduce ${\overline{U}}_z$ and $\overline{P}$ around $\overline{x}=0$. Maximum pressure applied in middle and end of cylinder for $P_{i4}$ and $P_{i5}$, respectively, in addition to descending thickness cause maximum of ${\overline{U}}_z$and $\overline{P}$near $\overline{x}=1$. It can be concluded from Figure 8(a) and (c) that ${\overline{U}}_z$and $\overline{P}$ patterns along the length of shell follow the pattern of the applied pressure function.

Figure 8:
The effect of non-uniform pressure profile on dimensionless (a) radial displacement (b) axial displacement (c) hydrostatic pressure along axial direction at middle layer ($n=\mathrm{2,}h(\overline{x})=h_1,M{C}_9$) and (d) internal pressure distribution along axial direction.

The influences of thickness profile function on the distribution of displacements and hydrostatic pressure under $P_{i1}$loading are illustrated in Figures 9- 11 for heterogeneous cylinder with $n=2$ and $M{C}_9$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, ${\overline{U}}_z$and $\overline{P}$ for $h_4$ profile intensify near $\overline{x}=0$ because of lower thickness in addition to maximum pressure. It is obviously observed that constant thickness $(h_0)$have the greatest displacements and stresses under$P_{i1}$. 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, $P_{i5},h{}_5$and $P_{i4},h{}_4$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 $P_{i5},h{}_5$ and $P_{i1},h{}_1$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.

Figure 9:
The effect of thickness profile on dimensionless radial displacement ($n=\mathrm{2,}{P}_i(\overline{x})=P_{i1},M{C}_9$).

Figure 10:
The effect of thickness profile on dimensionless axial displacement ($n=\mathrm{2,}{P}_i(\overline{x})=P_{i1},M{C}_9$).

Figure 11:
The effect of thickness profile on dimensionless hydrostatic pressure ($n=\mathrm{2,}{P}_i(\overline{x})=P_{i1},M{C}_9$).

Table 7:
Numerical results for similar pressure and thickness profiles in various sections at internal layer ($n=\mathrm{2,}M{C}_9$).

	${\overline{U}}_z$				${\overline{U}}_x$				$\overline{P}$
	$\overline{x}=0.15$	$\overline{x}=0.5$	$\overline{x}=0.85$		$\overline{x}=0.05$	$\overline{x}=0.25$	$\overline{x}=0.7$		$\overline{x}=0$	$\overline{x}=0.5$	$\overline{x}=1$
$P_{i0},h{}_0$	0.1255	0.1214	0.1255		0.0200	0.0067	-0.0053		17.751	5.402	17.751
$P_{i1},h{}_1$	0.0921	0.0819	0.0699		0.0155	-0.0099	-0.0268		11.177	3.876	10.419
$P_{i2},h{}_2$	0.0924	0.0865	0.0732		0.0167	-0.0067	-0.0283		11.448	4.186	10.893
$P_{i3},h{}_3$	0.0920	0.0882	0.0761		0.0173	-0.0040	-0.0265		11.455	4.335	11.325
$P_{i4},h{}_4$	0.0809	0.0909	0.0809		0.0205	0.0170	-0.0141		11.743	4.407	11.743
$P_{i5},h{}_5$	0.0914	0.0614	0.0914		0.0054	-0.0165	0.0158		3.741	3.407	3.741

5 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<\overline{R}<20$ and $P_i/\left(C_1+C_2\right)<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 $\overline{R}$because of intensifying nonlinear behavior of the cylinder while for small $\overline{R}$, the accuracy of shear deformation theory decrease in analyzing thick cylindrical shells. An increase in $P_i/(C_1+C_2)$and $k$or decrease in $C_1+C_2$ 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.

References

• Abdessamad, M., Mohamed, H., Mohamed, A., (2018). Analytical modeling of a descending aorta containing human blood flow. Defect and Diffusion Forum 384: 117-129.
• Anani, Y., Rahimi, G.H., (2015). Stress analysis of thick pressure vessel composed of functionally graded incompressible hyperelastic materials. International Journal of Mechanical Sciences 104: 1-7.
• Anani, Y., Rahimi, G.H., (2016). Stress analysis of rotating cylindrical shell composed of functionally graded incompressible hyperelastic materials. International Journal of Mechanical Sciences 108-109: 122-128.
• Azar, D., Ohadi, D., Rachev, A., Eberth, J.F., Uline, M.J., Shazly, T., (2018). Mechanical and geometrical determinants of wall stress in abdominal aortic aneurysms: A computational study. PLoS ONE 13(2): e0192032.
• Başar, Y., Weichert, D., (2000). Nonlinear Continuum Mechanics of Solids, Springer (Berlin).
• Batra, R.C., (2011). Material tailoring and universal relations for axisymmetric deformations of functionally graded rubberlike cylinders and spheres. Mathematics and Mechanics of Solids 16(7): 729-738.
• Batra, R.C., Bahrami, A., (2009). Inflation and eversion of functionally graded nonlinear elastic incompressible cylinders. International Journal of Non-Linear Mechanics 44(3): 311-323.
• Bijelonja, I., Demirdžic, I., Muzaferija, S., (2005). A finite volume method for large strain analysis of incompressible hyperelastic materials. International Journal for Numerical methods in Engineering 64: 1594-1609.
• Bilgili, E., (2003). Controlling the stress-strain inhomogeneities in axially sheared and radially heated hollow rubber tubes via functional grading. Mechanics Research Communications 30(3): 257-266.
• Bilgili, E., (2004). Modelling mechanical behavior of continuously graded vulcanized rubbers. Plastics Rubber and Composites 33(4): 163-169.
• Boyce, M.C., Arruda, E.M., (2000). Constitutive models of rubber elasticity: a review. Rubber Chemistry and Technology 73(3): 504-523.
• Dias, V., Odenbreit, C., Hechler, O., Scholzen, F., BenZineb, T., (2014). Development of a constitutive hyperelastic material law for numerical simulations of adhesive steel-glass connections using structural silicone. International Journal of Adhesion & Adhesives 48: 194-209.
• Doghri, I., (2000). Mechanics of Deformable Solids: Linear, Nonlinear, Analytical and Computational Aspects. Springer (Berlin).
• Doll, S., Schweizerhof, K., (2000). On the development of volumetric strain energy functions. Journal of Applied Mechanics 97: 17-21.
• Eipakchi, H.R., (2010). Third-order shear deformation theory for stress analysis of a thick conical shell under pressure. Journal of Mechanics of materials and structures 5(1): 1-17.
• Eipakchi, H.R., Rahimi, G.H., Esmaeilzadeh Khadem, S., (2003). Closed form solution for displacements of thick cylinders with varying thickness subjected to non-uniform internal pressure. Structural Engineering and mechanics 16(6): 731-748.
• Gao, H., Long, Q., Graves, M., Gillard, J.H., Li, Z.Y., (2009). Carotid arterial plaque stress analysis using fluid-structure interactive simulation based on in-vivo magnetic resonance images of four patients. Journal of Biomechanics 42: 1416-1423.
• Ghadiri Rad, M.H., Shahabian, F., Hosseini, S.M., (2015). Geometrically nonlinear elastodynamic analysis of hyper-elastic neo-Hooken FG cylinder subjected to shock loading using MLPG method. Engineering Analysis with Boundary Elements 50: 83-96.
• Ghaemi, H., Behdinan, K., Spence, A., (2006). On the development of compressible pseudo-strain energy density function for elastomers Part 1. Theory and experiment. Journal of Materials Processing Technology 178: 307-316.
• Ghannad, M., Rahimi, G.H., Nejad, M.Z., (2013). Elastic analysis of pressurized thick cylindrical shells with variable thickness made of functionally graded materials. Composites: Part B 45: 388-396.
• Gharooni, H., Ghannad, M., Nejad, M.Z., (2016). Thermo-elastic analysis of clamped-clamped thick FGM cylinders by using third-order shear deformation theory. Latin American Journal of Solids and Structures 13(4): 750-774.
• Holzapfel, G.A., (2000). Nonlinear Solid Mechanics, a Continuum Approach for Engineering, Wiley (New York).
• Holzapfel, G.A., Gasser, T.C., (2007). Computational stress-deformation analysis of arterial walls including high-pressure response. International Journal of Cardiology 116: 78-85.
• Humphrey, J.D., O’Rourke, S.L., (2015). An Introduction to Biomechanics Solids and Fluids, Analysis and Design, second ed., Springer (New York).
• Jabbari, M., Nejad, M.Z., Ghannad, M. (2016). Thermo-elastic analysis of axially functionally graded rotating thick truncated conical shells with varying thickness. Composites: Part B 96: 20-34.
• Kiendl, J., Hsu, M.C., Wu, M.C.H., Reali, A., (2015). Isogeometric Kirchhoff-Love shell formulations for general hyperelastic materials. Computer Methods in Applied Mechanics and Engineering 291: 280-303.
• Kim, J.H., Avril, J.H., Duprey, A., Favre, J.P., (2012). Experimental characterization of rupture in human aortic aneurysms using full-field measurement technique. Biomechanics and Modeling in Mechanobiology 11(6): 841-854.
• Lally, C., Dolan, F., Prendergast, P.J., (2005). Cardiovascular stent design and vessel stresses: a finite element analysis. Journal of Biomechanics 38: 1574-1581.
• Łos, M.J., Panigrahi, S., Sielatycka, K., Grillon, C., (2018). Successful Biomaterial-Based Artificial Organ-Updates on Artificial Blood Vessels. In: Los, M., Hudecki, A., Wiechec, E., (Eds.), Stem Cells and Biomaterials for Regenerative Medicine (203-222), Academic Press (United States).
• Mihai, L.A., Goriely, A., (2017). How to characterize a nonlinear elastic material? A review on nonlinear constitutive parameters in isotropic finite elasticity. Journal of Royal Society A 473(2207): 20170607.
• Montella, G., Calabrese, A., Serino, G., (2014). Mechanical characterization of a Tire Derived Material: experiments, hyperelastic modeling and numerical validation. Construction and Building Materials 66: 336-347.
• Nayfeh, A.H., (1981). Introduction to Perturbation Techniques, Wiley (New York).
• Nejad, M.Z., Jabbari, M., Ghannad, M., (2015). Elastic analysis of axially functionally graded rotating thick cylinder with variable thickness under non-uniform arbitrarily pressure loading. International Journal of Engineering Science 89: 86-99.
• Selvadurai, A.P.S., Shi, M., (2012). Fluid pressure loading of a hyperelastic membrane. International Journal of Non-Linear Mechanics 47: 228-239.
• Silva, C.A.C., Bittencourt, M.L., (2008). Structural shape optimization of 3D nearly-incompressible hyperelasticity problems. Latin American Journal of Solids and Structures 5(2): 129-156.
• Simo, J.C., Taylor, R.L., (1982). Penalty function formulations for incompressible nonlinear elastostatics. Computer Methods in Applied Mechanics and Engineering 35: 107-118.
• Sussman, T., Bathe, K.J., (1987). A finite element formulation for nonlinear incompressible elastic and inelastic analysis. Computers & Structures 26(112): 357-409.
• Tanveer, M., Zu, J.W., (2012). Non-linear vibration of hyperelastic axisymmetric solids by a mixed p-type method. International Journal of Non-Linear Mechanics 47: 30-41.
• Vossoughi, J., Tozeren, A., (1998). Determination of an effective shear modulus of aorta. Russian Journal of Biomechanics 1-2: 20-36.
• Xie, J., Zhou, J., Fung, Y.C., (1995). Bending of Blood Vessel Wall: Stress-Strain Laws of the Intima-Media and Adventitial Layers. Journal of Biomechanical Engineering 117: 136-145.
• Zhu, Y., Luo, X.Y., Ogden, R.W., (2010). Nonlinear axisymmetric deformations of an elastic tube under external pressure. European Journal of Mechanics- A/Solids 29(2): 216-229.

APPENDIX A

The non-homogeneity vectors of $O\left({\epsilon }^2\right)$ equations in outer and inner expansions are as follows:

$\{\begin{array}{l}{\left\{F_{O2}^{II}\right\}}_1=\overline{C{k}_3}[\left(\overline{I}\overline{I}(\mathrm{2,}n-1)+\overline{I}\overline{I}(\mathrm{0,}n+1)\right){\overline{\psi }}_{O1}{}^2+2\overline{I}\overline{I}(\mathrm{1,}n-1){\overline{w}}_{O1}{\overline{\psi }}_{O1}+2\left(\overline{I}\overline{I}(\mathrm{1,}n)+\overline{I}\overline{I}(\mathrm{0,}n+1)\right){\overline{v}}_{O1}{\overline{\psi }}_{O1}\\ +2\overline{I}\overline{I}(\mathrm{0,}n){\overline{v}}_{O1}{\overline{w}}_{O1}+\overline{I}\overline{I}(\mathrm{0,}n-1){\overline{w}}_{O1}{}^2]+\overline{C{k}_7}\left[\overline{I}\overline{I}(\mathrm{1,}n){\overline{\psi }}_{O1}{}^2+\overline{I}\overline{I}(\mathrm{0,}n){\overline{w}}_{O1}{\overline{\psi }}_{O1}\right]\\ -\overline{k}\left[\overline{I}\overline{I}(\mathrm{1,}n+1){\overline{\phi }}_{O1}{}^{\prime }+2\overline{I}\overline{I}(\mathrm{0,}n+1){\overline{v}}_{O1}{}^2\right]-{\widehat{K}}_s\overline{C_{12}}\overline{I}\overline{I}(\mathrm{0,}n+1){\overline{\phi }}_{O1}{}^2,{\left\{F_{O2}^{I{I}^{\prime }}\right\}}_1=0\end{array}$ (A-1)

$\{\begin{array}{l}{\left\{F_{O2}^{II}\right\}}_2=\overline{C{k}_1}[\left(\overline{I}\overline{I}(\mathrm{2,}n)+\overline{I}\overline{I}(\mathrm{1,}n+1)\right){\overline{\psi }}_{O1}{}^{\prime }+\overline{I}\overline{I}(\mathrm{1,}n){\overline{w}}_{O1}{}^{\prime }-\left(\overline{I}\overline{I}(\mathrm{1,}n){\overline{\psi }}_{O1}+\overline{I}\overline{I}(\mathrm{0,}n){\overline{w}}_{O1}\right){\overline{\phi }}_{O1}\\ +{\widehat{K}}_s\overline{I}\overline{I}(\mathrm{0,}n+1){\overline{v}}_{O1}{\overline{\phi }}_{O1}]+K_s\overline{C{k}_2}\left(\overline{I}\overline{I}(\mathrm{1,}n){\overline{\psi }}_{O1}+\overline{I}\overline{I}(\mathrm{0,}n){\overline{w}}_{O1}\right){\overline{\phi }}_{O1}\\ +K_s\overline{C_{12}}\left(\overline{I}\overline{I}(\mathrm{0,}n+1){\overline{w}}_{O1}{}^{\prime }+\overline{I}\overline{I}(\mathrm{1,}n+1){\overline{\psi }}_{O1}{}^{\prime }\right)-\overline{k}\overline{I}\overline{I}(\mathrm{1,}n+1){\overline{v}}_{O1}{}^{\prime }+{\widehat{K}}_s\overline{k}\overline{I}\overline{I}(\mathrm{0,}n+1){\overline{\phi }}_{O1}{\overline{\psi }}_{O1}\\ {\left\{F_{O2}^{I{I}^{\prime }}\right\}}_2=\overline{C{k}_1}[\overline{I}{\overline{I}}_{}{}^{\prime }(\mathrm{2,}n){\overline{\psi }}_{O1}+\overline{I}{\overline{I}}_{}{}^{\prime }(\mathrm{1,}n){\overline{w}}_{O1}+\overline{I}{\overline{I}}_{}{}^{\prime }(\mathrm{1,}n+1){\overline{\psi }}_{O1}]-\overline{k}\overline{I}{\overline{I}}_{}{}^{\prime }(\mathrm{1,}n+1){\overline{v}}_{O1}\end{array}$ (A-2)

$\{\begin{array}{l}{\left\{F_{O2}^{II}\right\}}_3=-\overline{C{k}_3}[\overline{I}\overline{I}(\mathrm{0,}n)\left({\overline{v}}_{O1}{}^2+{\overline{\psi }}_{O1}{}^2\right)+2\left(\overline{I}\overline{I}(\mathrm{1,}n-1){\overline{\psi }}_{O1}+\overline{I}\overline{I}(\mathrm{0,}n-1){\overline{w}}_{O1}\right)\left({\overline{v}}_{O1}+{\overline{\psi }}_{O1}\right)]\\ -\overline{C{k}_1}\overline{I}\overline{I}(\mathrm{1,}n){\overline{\phi }}_{O1}{}^{\prime }-\overline{C{k}_7}\overline{I}\overline{I}(\mathrm{0,}n){\overline{v}}_{O1}{\overline{\psi }}_{O1}-K_s\overline{C_{12}}\overline{I}\overline{I}(\mathrm{0,}n+1){\overline{\phi }}_{O1}{}^{\prime }+2\overline{C_2}\overline{I}\overline{I}(\mathrm{0,}n){\overline{\phi }}_{O1}{}^2\\ +2\overline{k}[\overline{I}\overline{I}(\mathrm{2,}n-2){\overline{\psi }}_{O1}{}^2+2\overline{I}\overline{I}(\mathrm{1,}n-2){\overline{w}}_{O1}{\overline{\psi }}_{O1}+\overline{I}\overline{I}(\mathrm{0,}n-2){\overline{w}}_{O1}{}^2]\\ {\left\{F_{O2}^{I{I}^{\prime }}\right\}}_3=-K_s\overline{C_{12}}\overline{I}{\overline{I}}_{}{}^{\prime }(\mathrm{0,}n+1){\overline{\phi }}_{O1}\end{array}$ (A-3)

$\{\begin{array}{l}{\left\{F_{O2}^{II}\right\}}_4=-\overline{C{k}_1}\left(\overline{I}\overline{I}(\mathrm{2,}n)+\overline{I}\overline{I}(\mathrm{1,}n+1)\right){\overline{\phi }}_{O1}{}^{\prime }-K_s\overline{C_{12}}\overline{I}\overline{I}(\mathrm{1,}n+1){\overline{\phi }}_{O1}{}^{\prime }-\overline{C{k}_3}[\overline{I}\overline{I}(\mathrm{1,}n)\left(3{\overline{\psi }}_{O1}{}^2+{\overline{v}}_{O1}{}^2\right)\\ +\overline{I}\overline{I}(\mathrm{2,}n-1){\overline{\psi }}_{O1}\left(3{\overline{\psi }}_{O1}+2{\overline{v}}_{O1}\right)+\overline{I}\overline{I}(\mathrm{0,}n+1){\overline{v}}_{O1}\left(2{\overline{\psi }}_{O1}+{\overline{v}}_{O1}\right)+\overline{I}\overline{I}(\mathrm{1,}n-1){\overline{w}}_{O1}\left(4{\overline{\psi }}_{O1}+2{\overline{v}}_{O1}\right)\\ +2\overline{I}\overline{I}(\mathrm{0,}n){\overline{w}}_{O1}{\overline{\psi }}_{O1}+2\overline{I}\overline{I}(\mathrm{0,}n-1){\overline{w}}_{O1}{}^2]-\overline{C{k}_7}\left[\overline{I}\overline{I}(\mathrm{0,}n){\overline{v}}_{O1}{\overline{w}}_{O1}+2\overline{I}\overline{I}(\mathrm{1,}n){\overline{v}}_{O1}{\overline{\psi }}_{O1}\right]\\ -2\overline{C_2}\overline{I}\overline{I}(\mathrm{1,}n){\overline{\phi }}_{O1}{}^2+2\overline{k}[\overline{I}\overline{I}(\mathrm{1,}n-2){\overline{w}}_{O1}{}^2+\left(\overline{I}\overline{I}(\mathrm{3,}n-2)+\overline{I}\overline{I}(\mathrm{0,}n+1)\right){\overline{\psi }}_{O1}{}^2+2\overline{I}\overline{I}(\mathrm{2,}n-2){\overline{w}}_{O1}{\overline{\psi }}_{O1}]\\ {\left\{F_{O2}^{I{I}^{\prime }}\right\}}_4=-K_s\overline{C_{12}}\overline{I}{\overline{I}}_{}{}^{\prime }(\mathrm{1,}n+1){\overline{\phi }}_{O1}\end{array}$ (A-4)

${\left\{F_{\alpha 2}^{P_{\alpha }}\right\}}_1={\left\{F_{\alpha 2}^{P_{\alpha }}\right\}}_2=\mathrm{0,}{\left\{F_{\alpha 2}^{P_{\alpha }}\right\}}_3=-{\overline{P}}_{i\alpha }{\tilde{x}}_{\alpha }\left(D{\overline{R}}_{\alpha }-\frac{D{\overline{h}}_{\alpha }}{2}\right),{\left\{F_{\alpha 2}^{P_{\alpha }}\right\}}_4=\frac{{\overline{P}}_{i\alpha }{\tilde{x}}_{\alpha }}{2}\left(D{\overline{h}}_{\alpha }{\overline{R}}_{\alpha }+{\overline{h}}_{\alpha }D{\overline{R}}_{\alpha }-{\overline{h}}_{\alpha }D{\overline{h}}_{\alpha }\right)$ (A-5)

${\left\{F_{\alpha 2}^{D{P}_{\alpha }}\right\}}_1={\left\{F_{\alpha 2}^{D{P}_{\alpha }}\right\}}_2=\mathrm{0,}{\left\{F_{\alpha 2}^{D{P}_{\alpha }}\right\}}_3=-D{\overline{P}}_{i\alpha }{\tilde{x}}_{\alpha }\left({\overline{R}}_{\alpha }-\frac{{\overline{h}}_{\alpha }}{2}\right),{\left\{F_{\alpha 2}^{D{P}_{\alpha }}\right\}}_4=\frac{D{\overline{P}}_{i\alpha }{\tilde{x}}_{\alpha }{\overline{h}}_{\alpha }}{2}\left({\overline{R}}_{\alpha }-\frac{{\overline{h}}_{\alpha }}{2}\right)$ (A-6)

$\begin{array}{l}{\left\{F_{\alpha 2}^{I{I}_{\alpha }}\right\}}_1=\overline{C{k}_3}[2\left(\left(\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1)+\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n)\right){\overline{\psi }}_{\alpha 1}+\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n){\overline{w}}_{\alpha 1}\right){\overline{\phi }}_{\alpha 1}{}^{\prime }+\left(\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)+\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n-1)\right){\overline{\psi }}_{\alpha 1}{}^2\\ +2\left(\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n)+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)\right){\overline{v}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}+2\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n-1){\overline{w}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}+2\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n){\overline{v}}_{\alpha 1}{\overline{w}}_{\alpha 1}\\ +\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n-1){\overline{w}}_{\alpha 1}{}^2]-\overline{C_{12}}\left[K_s\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1){\overline{\phi }}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime }+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)\left(K_s{\overline{\phi }}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime }+{\widehat{K}}_s{\overline{\phi }}_{\alpha 1}{}^2\right)\right]\\ -\overline{k}{[2\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1)\left({\overline{\phi }}_{\alpha 1}{}^{\prime }\right)}^2+4\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1)\left(4{\overline{v}}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }-{\overline{\phi }}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime }\right)+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)\left(2{\overline{v}}_{\alpha 1}{}^2-{\overline{\phi }}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime }\right)]\\ +\overline{C{k}_7}\left[\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n){\overline{\psi }}_{\alpha 1}{}^2+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n){\overline{w}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}\right]\end{array}$ (A-7)

$\begin{array}{l}{\left\{F_{\alpha 2}^{D{A}_{\alpha }}\right\}}_1={\tilde{x}}_{\alpha }\overline{C{k}_1}[\left(D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n)+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)\right){\overline{\psi }}_{\alpha 1}+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n){\overline{w}}_{\alpha 1}]\\ -{\tilde{x}}_{\alpha }\overline{k}\left(D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1){\overline{\phi }}_{\alpha 1}{}^{\prime }+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1){\overline{v}}_{\alpha 1}\right)\end{array}$ (A-8)

${\left\{F_{\alpha 2}^{DI{I}_{\alpha }}\right\}}_1=0$ (A-9)

$\begin{array}{l}{\left\{F_{\alpha 2}^{I{I}_{\alpha }}\right\}}_2=\overline{C{k}_1}\left(\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n){\overline{\psi }}_{\alpha 1}+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n){\overline{w}}_{\alpha 1}+{\widehat{K}}_s\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1){\overline{v}}_{\alpha 1}\right){\overline{\phi }}_{\alpha 1}+K_s\overline{C{k}_2}[\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n){\overline{\psi }}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime }\\ +\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n)\left({\left({\overline{w}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}\right)}^{\prime }+{\overline{\psi }}_{\alpha 1}{\overline{\phi }}_{\alpha 1}\right)+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n)\left({\overline{\phi }}_{\alpha 1}+{\overline{w}}_{\alpha 1}{}^{\prime }\right){\overline{w}}_{\alpha 1}]+\overline{C{k}_5}\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1){\overline{\psi }}_{\alpha 1}{}^{\prime }{\overline{\phi }}_{\alpha 1}{}^{\prime }\\ +K_s\overline{C{k}_4}[\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)\left({\overline{\psi }}_{\alpha 1}+{\overline{v}}_{\alpha 1}\right){\overline{w}}_{\alpha 1}{}^{\prime }+\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1)\left({\overline{v}}_{\alpha 1}+{\overline{\psi }}_{\alpha 1}\right){\overline{\psi }}_{\alpha 1}{}^{\prime }]-K_s\overline{C_{12}}[\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1){\overline{\psi }}_{\alpha 1}{}^{\prime \text{}\prime }\\ +\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1){\overline{w}}_{\alpha 1}{}^{\prime \text{}\prime }]{\overline{\phi }}_{\alpha 1}+\overline{C{k}_6}[\overline{I}{\overline{I}}_{\alpha }(\mathrm{3,}n){\left({\overline{\psi }}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }\right)}^{\prime }+\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1){\overline{\psi }}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime \text{}\prime }\\ +\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1)\left({\left({\overline{\psi }}_{\alpha 1}{\overline{v}}_{\alpha 1}\right)}^{\prime }+{\overline{\psi }}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime }\right)+\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n-1){\overline{w}}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime }+\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n-1){\left({\overline{w}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}\right)}^{\prime }\\ +\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n){\left({\overline{w}}_{\alpha 1}{\overline{v}}_{\alpha 1}\right)}^{\prime }+\overline{I}{\overline{I}}_{\alpha }(\mathrm{3,}n+1){\overline{\psi }}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime }+\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n)\left({\left({\overline{v}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}\right)}^{\prime }+{\left({\overline{w}}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }\right)}^{\prime }\right)]\\ -\overline{k}[{\widehat{K}}_s\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1){\overline{\phi }}_{\alpha 1}{\overline{\psi }}_{\alpha 1}+\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1)\left(4{\left({\overline{v}}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }\right)}^{\prime }+K_s{\overline{\psi }}_{\alpha 1}{}^{\prime }{\overline{\phi }}_{\alpha 1}{}^{\prime }-{\overline{\phi }}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime \text{}\prime }\right)\\ +\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1)\left(4{\overline{v}}_{\alpha 1}{\overline{v}}_{\alpha 1}{}^{\prime }-{\overline{\phi }}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime \text{}\prime }+{\widehat{K}}_s\left({\overline{w}}_{\alpha 1}{}^{\prime }+{\overline{\phi }}_{\alpha 1}\right){\overline{\phi }}_{\alpha 1}{}^{\prime }\right)+4\overline{I}{\overline{I}}_{\alpha }(\mathrm{3,}n+1){\overline{\phi }}_{\alpha 1}{}^{\prime }{\overline{\phi }}_{\alpha 1}{}^{\prime \text{}\prime }]\\ +\overline{C{k}_7}\left[2\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n){\overline{\psi }}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime }+\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n){\left({\overline{w}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}\right)}^{\prime }\right]\end{array}$ (A-10)

$\begin{array}{l}{\left\{F_{\alpha 2}^{D{A}_{\alpha }}\right\}}_2={\tilde{x}}_{\alpha }\overline{C{k}_1}\left[\left(D\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n)+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1)\right){\overline{\psi }}_{\alpha 1}{}^{\prime }+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n){\overline{w}}_{\alpha 1}{}^{\prime }\right]+{\tilde{x}}_{\alpha }{K}_s\overline{C_{12}}[D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1){\overline{\psi }}_{\alpha 1}{}^{\prime }\\ +D\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)\left({\overline{\phi }}_{\alpha 1}+{\overline{w}}_{\alpha 1}{}^{\prime }\right)]-{\tilde{x}}_{\alpha }\overline{k}\left(D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1){\overline{v}}_{\alpha 1}{}^{\prime }+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1){\overline{\phi }}_{\alpha 1}{}^{\prime \text{}\prime }\right)\end{array}$ (A-11)

${\left\{F_{\alpha 2}^{DI{I}_{\alpha }}\right\}}_2=\overline{C{k}_1}\left[\left(D\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n)+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1)\right){\overline{\psi }}_{\alpha 1}+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n){\overline{w}}_{\alpha 1}\right]-\overline{k}\left(D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1){\overline{v}}_{\alpha 1}+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1){\overline{\phi }}_{\alpha 1}{}^{\prime }\right)$ (A-12)

$\begin{array}{l}{\left\{F_{\alpha 2}^{I{I}_{\alpha }}\right\}}_3=\overline{C{k}_1}[\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n){\overline{\psi }}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime \text{}\prime }+\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n)\left({\overline{w}}_{\alpha 1}{}^{\prime }{\overline{\psi }}_{\alpha 1}{}^{\prime }+{\overline{w}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime \text{}\prime }+{\overline{\psi }}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime \text{}\prime }\right)+2\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n){\overline{w}}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime \text{}\prime }\\ -{\widehat{K}}_s\left(\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1){\left({\overline{\psi }}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime }\right)}^{\prime }+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1){\left({\overline{w}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}\right)}^{\prime }\right)]-\overline{C{k}_7}\left(\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n){\overline{v}}_{\alpha 1}+\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n){\overline{\phi }}_{\alpha 1}{}^{\prime }\right){\overline{\psi }}_{\alpha 1}\\ +4\overline{C_2}\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n){\overline{\phi }}_{\alpha 1}{}^2-\overline{C{k}_2}[\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n)\left({\widehat{K}}_s{\overline{\phi }}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime }+K_s\left({\left({\overline{w}}_{\alpha 1}{}^{\prime }\right)}^2+{\overline{w}}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime \text{}\prime }+{\overline{w}}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }\right)\right)\\ +\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n)({\widehat{K}}_s{\overline{\phi }}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime }+K_s\left(2{\overline{w}}_{\alpha 1}{}^{\prime }{\overline{\psi }}_{\alpha 1}{}^{\prime }+{\overline{w}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime \text{}\prime }+{\overline{\psi }}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime \text{}\prime }+{\overline{\psi }}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }\right))\\ +K_s\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n)\left({\left({\overline{\psi }}_{\alpha 1}{}^{\prime }\right)}^2+{\overline{\psi }}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime \text{}\prime }\right)]-\overline{C{k}_3}[2\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n-1){\overline{\psi }}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }+2\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n-1)\left({\overline{v}}_{\alpha 1}+{\overline{\psi }}_{\alpha 1}\right){\overline{w}}_{\alpha 1}\\ +2\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n-1)\left({\overline{w}}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }+{\overline{\psi }}_{\alpha 1}{}^2+{\overline{v}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}\right)+2\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n){\overline{v}}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n)\left({\overline{v}}_{\alpha 1}{}^2+{\overline{\psi }}_{\alpha 1}{}^2\right)\\ +2\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n){\left({\overline{\phi }}_{\alpha 1}{}^{\prime }\right)}^2]-K_s\overline{C{k}_4}\left[\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1){\left({\overline{\phi }}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }\right)}^{\prime }+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1){\left({\overline{v}}_{\alpha 1}{\overline{\phi }}_{\alpha 1}+{\overline{\phi }}_{\alpha 1}{\overline{\psi }}_{\alpha 1}\right)}^{\prime }\right]\\ +\overline{C{k}_8}\left[\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n){\left({\overline{\psi }}_{\alpha 1}{}^{\prime }\right)}^2+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n){\left({\overline{w}}_{\alpha 1}{}^{\prime }\right)}^2\right]+{\widehat{K}}_s\overline{k}[\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1){\left({\overline{\phi }}_{\alpha 1}{}^{\prime }{\overline{\psi }}_{\alpha 1}{}^{\prime }\right)}^{\prime }\\ +\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1){\left({\overline{v}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime }+{\overline{\phi }}_{\alpha 1}{}^{\prime }{\overline{w}}_{\alpha 1}{}^{\prime }\right)}^{\prime }+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1){\left({\overline{v}}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime }\right)}^{\prime }]\\ +2\overline{k}\left[\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n-2){\overline{\psi }}_{\alpha 1}{}^2+2\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n-2){\overline{w}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n-2){\overline{w}}_{\alpha 1}{}^2\right]\end{array}$ (A-13)

$\begin{array}{l}{\left\{F_{\alpha 2}^{D{A}_{\alpha }}\right\}}_3=-{\tilde{x}}_{\alpha }\overline{C{k}_1}\left[D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n){\overline{\phi }}_{\alpha 1}{}^{\prime }+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n)\left({\overline{\psi }}_{\alpha 1}+{\overline{v}}_{\alpha 1}\right)\right]-{\tilde{x}}_{\alpha }{K}_s\overline{C_{12}}[D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1){\overline{\psi }}_{\alpha 1}{}^{\prime \text{}\prime }\\ +D\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)\left({\overline{\phi }}_{\alpha 1}{}^{\prime }+{\overline{w}}_{\alpha 1}{}^{\prime \text{}\prime }\right)]+{\tilde{x}}_{\alpha }\overline{k}\left[D\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n-1){\overline{w}}_{\alpha 1}+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n-1){\overline{\psi }}_{\alpha 1}\right]\end{array}$ (A-14)

${\left\{F_{\alpha 12}^{DI{I}_{\alpha }}\right\}}_3=-K_s\overline{C_{12}}\left[D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1){\overline{\psi }}_{\alpha 1}{}^{\prime }+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)\left({\overline{\phi }}_{\alpha 1}+{\overline{w}}_{\alpha 1}{}^{\prime }\right)\right]$ (A-15)

$\begin{array}{l}{\left\{F_{\alpha 2}^{I{I}_{\alpha }}\right\}}_4=-\overline{C{k}_1}[\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1){\overline{\phi }}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime }+{\widehat{K}}_s\left(\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1){\overline{w}}_{\alpha 1}{}^{\prime \text{}\prime }+\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1){\overline{\psi }}_{\alpha 1}{}^{\prime \text{}\prime }\right){\overline{\psi }}_{\alpha 1}-\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n){\overline{w}}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime \text{}\prime }\\ -\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n)\left({\left({\overline{w}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime }\right)}^{\prime }+{\overline{\psi }}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime \text{}\prime }\right)-\overline{I}{\overline{I}}_{\alpha }(\mathrm{3,}n){\overline{\psi }}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime \text{}\prime }]-\overline{C{k}_9}\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1)\left({\overline{v}}_{\alpha 1}+{\overline{\psi }}_{\alpha 1}\right){\overline{\phi }}_{\alpha 1}{}^{\prime }\\ -\overline{C{k}_2}[\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n)\left(K_s\left({\left({\overline{w}}_{\alpha 1}{}^{\prime }\right)}^2+{\overline{w}}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }+{\overline{w}}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime \text{}\prime }\right)+{\widehat{K}}_s{\overline{\phi }}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime }\right)\\ +\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n)\left(K_s\left({\left({\overline{w}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime }\right)}^{\prime }+{\overline{\psi }}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime \text{}\prime }+{\overline{\psi }}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }\right)+{\widehat{K}}_s{\overline{\phi }}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime }\right)\\ +K_s\overline{I}{\overline{I}}_{\alpha }(\mathrm{3,}n)\left({\left({\overline{\psi }}_{\alpha 1}{}^{\prime }\right)}^2+2{\overline{\psi }}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime \text{}\prime }\right)]-\overline{C{k}_3}[\overline{I}{\overline{I}}_{\alpha }(\mathrm{3,}n){\left({\overline{\phi }}_{\alpha 1}{}^{\prime }\right)}^2+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n-1){\overline{w}}_{\alpha 1}{}^2\\ +2\left(\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n-1)+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)\right){\overline{v}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}+2\left(2\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n-1)+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n)\right){\overline{w}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}\\ +3\left(\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n)+\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n-1)\right){\overline{\psi }}_{\alpha 1}{}^2+\left(\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)+\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n)\right){\overline{v}}_{\alpha 1}{}^2+2\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n){\overline{v}}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }\\ +2\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n-1){\overline{v}}_{\alpha 1}{\overline{w}}_{\alpha 1}+2\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n-1){\overline{w}}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }+2\overline{I}{\overline{I}}_{\alpha }(\mathrm{3,}n-1){\overline{\psi }}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }+\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1){\left({\overline{\phi }}_{\alpha 1}{}^{\prime }\right)}^2]\\ -\overline{C{k}_4}[\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1)\left({\widehat{K}}_s{\left({\overline{\psi }}_{\alpha 1}{}^{\prime }\right)}^2+K_s{\left({\overline{\phi }}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }\right)}^{\prime }\right)+K_s\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1)\left({\left({\overline{\phi }}_{\alpha 1}{\overline{v}}_{\alpha 1}\right)}^{\prime }+{\overline{\psi }}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }\right)]\\ +\overline{C_{12}}\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)\left({\widehat{K}}_s{\left({\overline{w}}_{\alpha 1}{}^{\prime }\right)}^2+K_s{\overline{\phi }}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime }\right)+\overline{C{k}_8}\left[\overline{I}{\overline{I}}_{\alpha }(\mathrm{3,}n){\left({\overline{\psi }}_{\alpha 1}{}^{\prime }\right)}^2+\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n){\left({\overline{w}}_{\alpha 1}{}^{\prime }\right)}^2\right]\\ -\overline{C{k}_7}[2\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n){\overline{\psi }}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }+2\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n)\left({\overline{v}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}+{\overline{w}}_{\alpha 1}{\overline{\phi }}_{\alpha 1}{}^{\prime }\right)+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n){\overline{v}}_{\alpha 1}{\overline{w}}_{\alpha 1}]\\ +{\widehat{K}}_s\overline{k}[\overline{I}{\overline{I}}_{\alpha }(\mathrm{3,}n+1){\left({\overline{\phi }}_{\alpha 1}{}^{\prime }{\overline{\psi }}_{\alpha 1}{}^{\prime }\right)}^{\prime }+\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1)\left({\left({\overline{\phi }}_{\alpha 1}{}^{\prime }{\overline{w}}_{\alpha 1}{}^{\prime }\right)}^{\prime }+{\left({\overline{v}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}{}^{\prime }\right)}^{\prime }\right)\\ +\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1)\left({\left({\overline{v}}_{\alpha 1}{\overline{w}}_{\alpha 1}{}^{\prime }\right)}^{\prime }+\left({\overline{w}}_{\alpha 1}{}^{\prime }+{\overline{\phi }}_{\alpha 1}\right){\overline{\psi }}_{\alpha 1}{}^{\prime }\right)]+2\overline{k}[\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n-2){\overline{w}}_{\alpha 1}{}^2\\ +\left(\overline{I}{\overline{I}}_{\alpha }(\mathrm{3,}n-2)+\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)\right){\overline{\psi }}_{\alpha 1}{}^2+2\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n-2){\overline{w}}_{\alpha 1}{\overline{\psi }}_{\alpha 1}]\end{array}$ (A-16)

$\begin{array}{l}{\left\{F_{\alpha 2}^{D{A}_{\alpha }}\right\}}_4=-{\tilde{x}}_{\alpha }\overline{C{k}_1}[D\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)\left({\overline{v}}_{\alpha 1}\right)+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n)\left(2{\overline{\psi }}_{\alpha 1}+{\overline{v}}_{\alpha 1}\right)+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n)\left({\overline{w}}_{\alpha 1}\right)\\ +\left(D\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n)+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1)\right){\overline{\phi }}_{\alpha 1}{}^{\prime }]-{\tilde{x}}_{\alpha }{K}_s\overline{C_{12}}\left[D\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1){\overline{\psi }}_{\alpha 1}{}^{\prime \text{}\prime }+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1)\left({\overline{\phi }}_{\alpha 1}{}^{\prime }+{\overline{w}}_{\alpha 1}{}^{\prime \text{}\prime }\right)\right]\\ +{\tilde{x}}_{\alpha }\overline{k}\left[\left(D\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n-1)+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{0,}n+1)\right){\overline{\psi }}_{\alpha 1}+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n-1){\overline{w}}_{\alpha 1}\right]\end{array}$ (A-17)

${\left\{F_{\alpha 2}^{DI{I}_{\alpha }}\right\}}_4=-K_s\overline{C_{12}}\left[D\overline{I}{\overline{I}}_{\alpha }(\mathrm{2,}n+1){\overline{\psi }}_{\alpha 1}{}^{\prime }+D\overline{I}{\overline{I}}_{\alpha }(\mathrm{1,}n+1)\left({\overline{\phi }}_{\alpha 1}+{\overline{w}}_{\alpha 1}{}^{\prime }\right)\right]$ (A-18)

Edited by

Editor:

Rogério José Marczak

Publication Dates

• Publication in this collection
  28 Oct 2019
• Date of issue
  2019

History

• Received
  11 Mar 2019
• Reviewed
  10 July 2019
• Accepted
  08 Sept 2019
• Published
  10 Sept 2019