Static Response of Double-Layered Pipes via a Perturbation Approach

: A double-layered pipe under the effect of static transverse loads is considered here. The mechanical model, taken from the literature and constituted by a nonlinear beam-like structure, is constituted by an underlying Timoshenko beam, enriched with further kinematic descriptors which account for local effects, namely, ovalization of the cross-section, warping and possible relative sliding of the layers under bending. The nonlinear equilibrium equations are addressed via a perturbation method, with the aim of obtaining a closed-form solution. The perturbation scheme, tailored for the speciﬁc load conditions, requires different scaling of the variables and proceeds up to the fourth order. For two load cases, namely, distributed and tip forces, the solution is compared to that obtained via a pure numeric approach and the ﬁnite element method.


Introduction
Beams are structural members which are massively used in several application fields. Examples come from civil engineering, e.g., bridges and buildings, from aerospace engineering, e.g., helicopters and aircrafts, from industrial engineering, e.g., turbines and automotive constructions, and many other fields. The requirement for high structural performance of such members led designers to optimize the choice of the constituting materials and their ways of use. Historically, the advent of composite materials [1,2], for instance, represented a turning point which brought a clear enhancement in the structural efficiency of the members; concurrently, it required a redefinition of the classical beam theories (see [3]), which, however, were mostly developed in linear field and for isotropic materials. In this context, in [4], an in-depth treatment of beam theory is presented, with argumentation addressed to thin-walled beams and composite material applications. In [5], the one-dimensional theory of composite thin-walled beams is consistently derived from a three-dimensional continuum in the framework of elasticity theory.
Accounting for local distortion effects in beams-for example, changes in shapes of the cross-sections, such as warping and ovalization-becomes crucial when dealing with thin-walled members. For instance, the contribution of warping in non-uniform torsion is explicitly addressed in Vlasov's theory [6], whereas Brazier's theory [7] explains softening effects in bending of pipes as a consequence of ovalization of the cross-section, also in the presence of a soft core [8]. The aforementioned aspects may be addressed in the framework of the generalized beam theory (GBT; see, e.g., [9]) as well. It represents a compelling instrument where, still dealing with one-dimensional continua, the local effects are apriori determined and described by assumed functions, amplified by parameters which are evaluated through equilibrium conditions. In this context, the choice of the assumed functions, as an outcome of the cross-section analysis, denotes a main step, as discussed in [10][11][12], where the use of free-dynamical modes of equivalent frames is proposed, and in case of curvilinear cross-sections [13].
Beam-like structures are also conveniently used as homogeneous counterparts of lattices, when the latter are made of periodic repetitions of a modular cell in one direction. This is the case, for instance, for tall buildings where a shear-shear-torsional beam model can be used for both static [14,15] and dynamic behavior analysis [16][17][18]. In these papers, the kinematic and static (or dynamic) problems relevant to the beam model (referred to as coarse model) are directly defined, while the response law for the homogeneous material, as usual, is deduced through an identification procedure from a refined model.
In [19] a beam-like model is proposed to address nonlinear statics and dynamics of a tubular beam; there, as an extension of the classical Euler-Bernoulli beam model, further kinematic descriptors besides the classical ones are introduced in order to account for local distortion, whereas the constitutive law is obtained from a consistent three-dimensional fiber model. Along the same lines, in [20], nonlinear statics of double-layered pipes are addressed using an underlying Timoshenko beam, combined to a local distortion model which draws inspiration from GBT. In particular, the local distortion variables introduced therein represent amplitudes of assumed functions, which describe changes in shape of the annular cross-section, and the obtained coupled equilibrium equations in terms of the kinematic descriptors are solved using numerical tools.
In this paper, starting from the model of a double-layered pipe developed in [20], and with the aim of obtaining a closed-form description of the static behavior, the nonlinear equilibrium equations in terms of kinematic descriptors are addressed through a specific perturbation scheme. As an essential tool to get a clear insight of the possible occurring phenomena, the closed-form solution is obtained for two case-studies, namely, for both distributed and tip loads. A systematic comparison of the outcomes to pure numerical approaches and finite element models is carried out, so as to validate the steps performed.
The paper is organized as follows: In Section 2 the beam-like model of a pipe is briefly recalled and the equilibrium equations are shown; in Section 3 a perturbation scheme is applied to get the closed form of the solution for two case-studies; in Section 4 numerical results are shown; and finally, conclusions are drawn in Section 5.

Model Description
The object of this study is a double-layered beam with a thin annular cross-section (see Figure 1). The beam length is l, the average radius of the annular cross-section is R and the two plies have thicknesses h i and h e , where the subscripts i and e stand for internal and external, respectively (h i , h e R). A thinner inter-layer is present as well, constituted by the adhesive, and its thickness is h a h i , h e . A beam-like model is used here, with the aim of roughly describing the in-plane static behavior of the aforementioned pipe, under the action of transversal forces. The model was taken from the literature (see [20] for details), and it is briefly recalled here for the sake of completeness. It is constituted by an initially straight axis-line, whose direction coincides with the unitary vectorā x , and which is spanned by the abscissa s ∈ [0, l]. Cross-sections are initially transverse to the axis-line and are connected to it at their center of gravity. With reference to Figure 2, they initially lie on the plane spanned by the unitary vectors a y andā z , where (ā x ,ā y ,ā z ) is the canonical basis. The points of the axis-line perform displacements u(s) = u(s)ā x + v(s)ā y , i.e., limited to the (ā x ,ā y )-plane, and the crosssections perform independent rotations about axisā z of amplitude ϑ(s), as in a planar Timoshenko beam.
Besides (u(s), v(s), ϑ(s)), further scalar kinematic descriptors are defined in the model, in order to account for local distortions of the beam, i.e., changes in shape for the crosssections: indicated as (a p (s), a w (s), a g (s)), the physical interpretation of the distortion variables is given in [20] and comes from the expression of the displacement of a generic point of the beam, seen as a corresponding three-dimensional body and functionally introduced for homogenization purposes. Specifically, as in the spirit of generalized beam theory (GBT-see, e.g., [9]), each of them represents the amplitude of a trial function, defined in the cross-section domain and describes an elementary distortion, namely, ovalization (a p (s)) and warping (a w (s)) of the cross-section itself, and longitudinal sliding of the plies (a g (s)) due to local opposite bending, as shown in Figure 3. The strain-displacement relationship is consistently introduced: classical strain components for Timoshenko beam are defined, namely, longitudinal strain ε 0 , shear strain γ 0 and bending curvature κ 0 , plus strain components relevant to each distortion variables, namely, α j , which are equal to the distortion variable themselves, and β j , which are their gradient (for j = p, w, g). The strain-displacement relationship follows: where prime stands for s-derivative. Equation (1) is combined to geometrical boundary conditions, evaluated at the restrained cross-sections = 0 and/ors = l.
whereȗs,vs,θs,ȃ js are imposed displacements. Time derivative (indicated with the dot) of Equation (1) provides the strain-rates: which are functional to define the internal virtual power, after introducing mechanical quantities dual to each strain-rate component: normal force N, shear force T and bending moment M, and distortion forces D j and bi-forces B j . Imposing the validity of the virtual power theorem, i.e., equating the internal power to the external one, where imposed external forces spend power in the velocity field, one gets: where p u , p v are external forces per unit length in directionā x ,ā y , respectively, c is external bending moment per unit length, q j are external distortion forces per unit length whereas, at cross-sections = 0 and/ors = l, P us , P vs are external tip forces in directionā x ,ā y , respectively, Cs is external tip bending moment, Q js are external tip distortion forces. After localization of Equation (4), the following equilibrium equations are obtained: together with the boundary conditions: It is worth noting that in both the strain-displacement equations, Equations (1) and (2), and the equilibrium equations, Equations (5) and (6), the global variables (i.e., those relevant to the Timoshenko beam model) are uncoupled to the local variables (i.e., those relevant to the cross-section distortion).
The constitutive relationship is evaluated in case of hyperelastic materials and comes from the definition of a strain potential, whose expression is identified after a homogenization procedure from a suitable three-dimensional continuum (see [20]). It turns out to describe a nonlinear elastic material (up to the third order), where coupling between global and local variables actually occurs: Elastic coefficients c n , n = 1, . . . , 31 are defined in Appendix A. It is worth noting that the coefficients c n depend on the longitudinal and transversal elastic moduli of both internal and external plies, namely, E i and E e , G i and G e , and on their thickness h i , h e ; furthermore, coefficients c 17 , c 21 depend also on the transversal elastic modulus of the adhesive G a and its thickness h a .
Equations (1), (2), (5)-(7) define the elastic problem for the homogeneous beam. Substitution of Equation (1) in (7), and then in (5), provides the following equilibrium equations in terms of displacement (where only the linear part is made explicit due to the large quantity of nonlinear terms): while the nonlinear terms, expanded up to the third order, are collected in the functions F n (u, v, ϑ, a p , a g , a w ), n = 1, . . . , 6, whose expressions are shown in Appendix B.
A clamped-free beam is analyzed, leading to the following boundary conditions: and The functions G n (u, v, ϑ, a p , a g , a w ), n = 1, . . . 6, collect the nonlinear boundary terms truncated up to the third order and are shown in Appendix B.

Perturbation Analysis
The nonlinear boundary value problem (8)-(10) is addressed via a perturbation method, so as to get an asymptotic expression for the solution. The cases of applied distributed transverse load p v and tip load P v are considered, whereas the other loads are assumed as null, i.e., Furthermore, similar but not equal geometric and mechanic features of the two plies are assumed, i.e., Young modulus E i E e , transverse elastic modulus G i G e and thickness h i h e .
The dependent variables are series-expanded after introducing the small positive parameter, ε 1, as: where the odd and even natures of the dominant part of v, ϑ and of u, a p , a w are exploited, respectively, in case of transverse load application; more specifically, the perturbation scheme is based on the idea that the loads directly trigger the global bending problem (variables v, ϑ), at the first step; then, in turn, at the second step, the bending problem triggers the axial and local distortion problems (variables u and a p , a w , a g ) by the nonlinear coupling; finally, higher order corrections to the global bending and local distortion problems, respectively, are evaluated. Since the coefficients c 5 , c 11 , c 12 , c 13 , c 14 tend to zero in case of exactly equal internal and external plies ( , then they are scaled at order , as: giving rise to vanishing linear coupling between the variables a g and v, ϑ, and suggesting for the former to start at the second order (Equation (12) 6 ). Moreover, the applied loads are considered of order , so that they appear at the the leading order, i.e., Under the aforementioned assumptions, after substituting Equations (12)- (14) in Equations (8)-(10) and letting to zero the terms at the same powers of , the following perturbation equations and boundary conditions are obtained, at order : at order 2 : c 19 a p,2 (l) + c 20 a w,2 (l) = 0 c 30 a w,2 (l) = 0 c 25 a g,2 (l) = −c 11 ϑ 1 (l) (16) at order 3 : and at order 4 : c 19 a p,4 − c 17 a p,4 + c 20 a w,4 = (c 12 − c 13 )a g,2 ϑ 1 + c 10 a p,2 ϑ 2 1 + 2c 18 ϑ 1 ϑ 3 − c 13 a g,2 ϑ 1 c 30 a w,4 − c 26 a w,4 − c 20 a p,4 = c 3 a g,2 a g,2 + c 14 a g,2 ϑ 1 + (c 28 − 2c 31 )a w,2 a w,2 + c 15 a w,2 ϑ 2 4 (0) = a g,4 (0) = a w,4 (0) = 0 c 1 u 4 (1) = −c 3 a g,2 (l) 2 + c 6 a w,2 (l)ϑ 1 (l)ϑ 1 (l) + c 5 a g,2 (l)ϑ 1 (l) − c 2 a w,2 (l) 2 + (c 4 − c 1 )ϑ 3 (l)v 1 (l) c 19 a p,4 (l) + c 20 a w,4 (l) = −c 13 a g,2 (l)ϑ 1 (l) c 30 a w,4 (l) = −c 31 a w,2 (l) 2 c 25 a g,4 (l) = −c 3 a g,2 (l)a w,2 (l) − c 12 a p,2 (l)ϑ 1 (l) − c 11 ϑ 3 (l) (18) which can be chain-solved. It is worth noting how Equation (15) represents the bending boundary value problem for a linear Timoshenko beam under the assigned load. Its solution appears at the right-hand side of system (16), which rules the longitudinal displacement and the local distortion, where the equations in a p,2 , a w,2 are linearly coupled. Then, Equation (17) provides corrections to the bending problem, as function of first order bending, longitudinal displacement and local distortion. Finally, Equation (18) provides the correction for the longitudinal displacement and local distortion variables.

Solution for the Distributed Load Case
In case of uniformly distributed vertical load, solution to problem (15) is: As when solving classical Timoshenko beam problem, Equation (19) is obtained pursuing the following steps: (1) integration of both the sides of Equation (15) 1 to evaluate c 4 (v 1 − ϑ 1 ) (one arbitrary constant is introduced); (2) substitution of that in Equation (15) 2 and double integration of both the sides of it to obtain ϑ 1 (two more arbitrary constants are introduced); (3) substitution of the obtained ϑ 1 in the integrated Equation (15) 1 and one more integration of both its sides, to evaluate v 1 (one more arbitrary constant is introduced). The four arbitrary constants are evaluated imposing the four boundary conditions (b.c.) in Equation (15) 3,4 .
To continue the chain-solving procedure, Equation (19) is substituted to the right-hand side of Equation (16). The solution for variable u 2 is directly obtained with two steps of integration of Equation (16) 1 and use of its b.c.: and for variable a g,2 , after double integration of Equation (16) 4 and use of its b.c.: The solution to the coupled problem Equation (16) 2,3 in the variables a p,2 and a w,2 is written here in compact form, due to its large dimensions. In particular, its expression is evaluated with the method of variation of the constants: a vector of state variables is defined as y 2 = (a p,2 , a w,2 , a p,2 , a w,2 ) T , where the superscript T indicates the transpose, and Equations (16) 2,3 are written in the state-space form as: y 2 (s) = Ay 2 (s) + f 2 (s) (22) where: Then, its solution turns out to be: where: with λ i eigenvalues of A, which are: which is an algebraic non-homogeneous system in the unknown α, with non-singular system matrix. Equation (29) is solved to get α, and despite the large shape, full expression of the solution (24) is easily obtained with the help of an algebraic manipulator software [21]. Solutions to Equations (17) and (18), i.e., at cubic and quartic perturbation orders, are obtained with the same procedure as above, but they are not reported here for the sake of brevity.

Solution for the Tip Load Case
Analogously to what done in the previous section, in case of tip vertical load, solution to problem (15) is: After substituting Equation (30) into the right-hand side of Equation (16), the solutions for variables u 2 and a g,2 are easily evaluated, as: and: However, to obtain expressions for a p,2 and a w,2 , still Equations (24) and (29) are called for, where the load vector is now: The solutions of the cubic and fourth order problem can be found in closed form as well, but as for the distributed load case, they are not reported here.

Numerical Results
A case-study is considered to validate the perturbation procedure. The geometrical and mechanical parameters adopted for the pipe are: length l = 2 m; average radius of the cross-section R = 0.1 m; thicknesses of the outer and inner layers h e = h i = 2.0 mm, respectively. The Young's modulus for the outer layer is E e = 3.0 × 10 8 Pa, and for the inner layer it is E i = 3.0 × 10 7 Pa; both materials are assumed to have a zero Poisson's ratio, so that related effects are neglected (anyway, the contribution of Poisson's modulus is known to produce quantitative modifications on the Brazier effect of order less than 10%, in case of common materials (see [7])). The thin adhesive layer, interposed between the inner and outer ones, has a thickness of h a = 0.1 mm and a shear modulus equal to G a = 1.5 × 10 6 Pa.
The uniformly distributed load is assumed p v = 174 N/m, whereas the tip load is The obtained solution is superimposed on a numerical one, retrieved by applying a finite difference scheme to the nonlinear boundary value problem (8)-(10), after having divided the domain into n = 100 equally-spaced nodes, and using centered differences in internal domains and backward differences in the final node, where natural boundary conditions are assigned. Moreover, a further comparison of the solution to the outcomes of a Finite Element model implemented on a commercial software [22] is given: in the FEM-based approach, the pipe is modeled as an assembly of multi-layer shell elements (quadratic serendipity type), assuming finite kinematics and making use of the equivalent single layer theory (ESL) [1,2]; details are given in Appendix D.

Distributed Load
The outcomes for the first load case are shown as functions of s/l in Figure 4, where (1) the perturbation solution is marked by the dark blue lines; (2) the numerical one obtained via the finite difference method is marked by light blue lines; and (3) that provided by the FEM model is marked by gray lines (note that variables a g and a w as given are not directly observable by the FEM model). As expected, the effect of the distributed vertical load is evident on the variables v and ϑ (Figure 4b,c), which are directly activated at the linear order. However, due to the nonlinear coupling, the longitudinal displacement u is activated as well (Figure 4a), together with the local distortion variables (Figure 4d-f). In particular, the cross-section ovalization amplitude a p attains its maximum approximately at s/l = 0.35, where the cross-section warping a w assumes zero value. Finally, the inter-layer sliding a g , whose behavior is shown in Figure 4e, shows an almost linear behavior along the most part of the domain and a boundary layer close to the clamp (at s = 0), where it suddenly goes to zero. It is remarked that the perturbation solution, with higher order corrections, well captures the response of the structure: marginal differences appear with the finite difference solution, since the curves are overlapped in the entire domain. Satisfying agreement is found in the comparison with the FEM solution, as well.
In Figure 5 the evolution of the displacement at the tip v(l)/l (in abscissa) is shown when increasing the intensity of the loadp := p v l πR(G e h e +G i h i ) (in ordinate): the softening behavior due to the Brazier effect is well caught by the finite difference solution, and the perturbation solution correctly catches most of the curve except nearby the limit point. There, the agreement is less satisfying, perhaps due to the large distance of the limit point from the origin, and hence to the necessity of higher order approximation. The FE solution shows satisfying agreement as well, even if it is not very precise nearby the limit point, due to the approaching numerical instability. The dashed line in Figure 5 indicates the load condition corresponding to the previous Figure 4.

Concentrated Load
The response of the structure to the concentrated load P v at s = l is now analyzed. The results are illustrated in Figure 6 in terms of the same quantities discussed in the previous case and the same color legend.
As observed before, the perturbation solution well approximates the one obtained via the finite difference approach and it is also in good agreement with the FEM solution. The applied load P v activates directly the variables v and ϑ (Figure 6b,c), which in turn activate the longitudinal displacement u (Figure 6a) and the local distortion variables (Figure 6d-f). In particular, the inter-layer sliding a g exhibits an almost constant evolution, except for the presence of the boundary layer at the clamp ( Figure 6e); furthermore, for the same variable, the finite difference solution reveals a low amplitude oscillation around a mean value that is instead well captured by the perturbation solution.

Conclusions
Nonlinear statics of a double-layered pipe under the action of transverse forces is addressed here through a perturbation method. The beam-like model, taken from the literature, provides nonlinear equilibrium equations in terms of kinematic descriptors, where those relevant to the classic Timoshenko beam are combined to further variables and introduced to account for local distortion. The perturbation scheme, tailored for the two loading conditions, which consist of distributed and tip forces, respectively, requires different scaling for the involved variables, in order to exploit their even or odd nature. The obtained closed-form solutions, evaluated up to the fourth perturbation order, guarantee very good agreement with those obtained via pure numeric tools, and concurrently, represent a valid tool to discuss the nature of the response, and potentially to optimize the choice of the relevant parameters.
In both the load conditions, the non-negligible contribution of the ovalization, which is directly coupled to the warping amplitude, appears as correctly determined by the asymptotic solution; the same happens for the small but significant effect of the longitudinal sliding, which shows a boundary layer close to the clamped cross-section. Moreover, the longitudinal displacement, which is triggered as a consequence of the nonlinear coupling, is non-negligible as well.
It is worth noting that in both the analyzed cases, good agreement of the outcomes provided by the perturbation method was obtained, even if the numerical values of the parameters slightly violated the hypotheses under which the scaling was chosen, being E e = 10E i . This aspect represents a typical strength of the perturbation methods.
Moreover, the perturbation method also guides one to an enriching mechanical interpretation of the phenomena, as a consequence of the scaling proposed here, namely: (1) the external loads directly trigger the bending problem; (2) bending in turn induces ovalization, warping and axial displacement, as secondary and nonlinear effects; (3) the local distortion in turn gives high-order corrections to the leading global bending; (4) high-order correction to the local distortion is finally evaluated.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix B. Details about the Governing Equations
The expressions of the functions F i with i = 1, . . . , 6, in Equation (8) are: where terms in G i , with i = 1, . . . , 6, are evaluated in s = l.

Appendix C. Eigenvectors of the (a p , a w ) Problem
The matrix V is composed by the eigenvectors ϕ i , i = 1, . . . , 4, of the matrix A, namely,

Appendix D. Finite Element Model
The tubular structure is modeled in a FE software ( [22]; see the sketch in Figure A1a), which adopts the equivalent single layer (ESL) theory. It allows one to take in to account the composite nature of structure and the actual configuration of the through-thickness stacking sequence provided that the overall thickness is sufficiently small (∑ n h n R, being n the number of plies). The structure is thus modeled as an assembly of shell elements (quadratic serendipity), and a convergence sensitivity study has been conducted (results not shown here) to find a suitable compromise between the solution accuracy and computational burden, leading to a mesh of 5040 elements that is represented in Figure A1b.
Of course, the displacement measures descending from the FE model cannot be directly compared to those of the proposed model, whose displacement parameters are referred to the beam axis-line. However, a proper comparison can be conducted after a post-processing of the FE model specific displacements. In particular, with reference to Figure A1a, the load is applied in the −y direction, and exploiting the symmetry of the problem about the (x, y)-plane: • The longitudinal displacement u is equal to the x-displacement of the points along the E,W generator lines; • The transverse displacement v is equal to the y-displacement of the points along the E,W generator lines; • The section rotation ϑ is retrieved as the difference between the x-displacement of the points along the N and S lines, divided by 2R; • The ovalization amplitude a p is retrieved as the mid-difference between the ydisplacement of the points along the S and N lines.
Unfortunately, a g and a w cannot be straightforwardly extracted from the FE model, as it is done for the other variables; however, for the purposes of the present work, the comparison in terms of u, v, ϑ and a p is sufficient to evaluate the agreement between the proposed solution and the FEM simulation.