Structural analysis of composite tubes using a meshless analytical dimensional reduction method

A polynomial‐based model in conjunction with dimensional reduction method are presented to perform cross‐sectional analysis and to determine strain distribution in composite tubes under bending loading. For beam structures with tubular cross‐section, the Variational Asymptotic Method (VAM) has been employed to decompose a three‐dimensional (3D) elasticity problem into a two‐dimensional cross‐sectional analysis and a one‐dimensional analysis along the length. This greatly reduces computational time as compared to 3D Finite Element Method (FEM). There also exists publically available Variational Asymptotic Beam Sectional Analysis (VABS), a FEM cross‐sectional analysis tool based on VAM. For VABS, FEM mesh for the beam section needs to be generated. For the case of composite tubes with many layers, the mesh generation consumes efforts and is unnecessary. We introduce a new meshless dimensional reduction method for the analysis of composite tubes. This method utilizes Pascal polynomials in polar coordinates to model the warping functions. Using this, one can obtain stiffness constants and 3D strains of the composite tube. This method is straightforward, meshless, with similar computation time as VABS, and is much more efficient than conventional 3D FEM. The accuracy of the proposed method is examined by comparing the obtained results with 3D FE (ANSYS), VABS, literature, and experiment.

should be modeled using one element. For a structure containing 50 layers or more, with a cross-section diameter of about 50 mm and a length of about 1200 mm, the number of elements can go up to several hundred thousand elements. In one case study (shown in Section 5.3.1 in this paper), using ANSYS (commercial software) 3D FE requires 3192 k elements with a computer run time of 55 mi, while similar results can be obtained within 15 s using a more efficient method (also presented later in this paper).
For structures where there exists one (or two) dimensions (such as thickness or radius of the cross-section) that are much smaller than the other dimension (such length), one approach to reduce the computation effort in deformation analysis is the dimensional reduction method. Many previous researchers have used this approach.
Giavotto et al. 2 presented a numerical dimensional reduction method based on FE for cross-sectional discretization of an arbitrary cross-section of a beam. They extracted a 6 × 6 stiffness matrix of the cross-section with considering the shear effects. They developed a computer code HANBA2. Borri et al 3 extended the application of the method by Giavotto et al. to the case of naturally curved and twisted beams. The most significant work in the analysis of composite beams based on FE was developed by Hodges et al. [3][4][5] The Variational Asymptotic Method (VAM) is employed to minimize the energy functional of the beam. Hence, there is no ad hoc assumption in the 3D beam analysis. The approach includes all the nonclassical effects and splits the original 3D geometrically nonlinear problem of a beam into a one-dimensional (1D) global beam analysis and a two-dimensional (2D) cross-sectional analysis. Generally, based on first and second approximations of the strain energy, a certain beam model is obtained including classical beam model which is represented by a 4 × 4 stiffness matrix without effects of transverse shears, and a generalized Timoshenko beam model which leads to a 6 × 6 stiffness matrix which takes into account transverse shear effect. 6,7 Berdichevsky et al., 8 applied VAM to cross-sectional analysis of thin-walled anisotropic beams which led to a closed-form solution for the 4 × 4 stiffness matrix. Hodges et al. 9 used VAM to analyze of initially curved and twisted composite beams. The FEM is used to solve the nonlinear equations for nonlinear static deformation and linearized free vibration about the static state of the deformation. Popescu and Hodges 10 presented a FE cross-sectional analysis of beams which accounts for shear effects using VAM. Both the classical as well as Timoshenko beam models results for different cases of cross-section geometry have been presented.
A FE-based computer code called Variational Asymptotic Beam Sectional Analysis (VABS) was developed for 2D analysis based on VAM. 11 VABS extracts cross-sectional stiffness constants needed for 1D analysis of a FEA software. It also carries out recovery solution to provide stress and strain contour distribution in the cross-section. 4,7 It uses 2D FE for the cross-sectional analysis. The computational time is decreased considerably compared to 3D FE analysis. 11 This can be considered as the main advantage of the VABS over 3D FE analysis.
In the case of dimensional reduction method for tubular cross-sections using VAM, several works have been carried out. Based on the Timoshenko refinement of the VAM, Popescu and Hodges 10 compared the obtained stiffness constants of a composite Circumferentially Uniform Stiffness (CUS) beam using FE for cross-sectional discretization with those obtained Rehfield et al. 12 in which they used nonclassical beam theory. Ghafari and Rezaeepazhand 13 applied dimensional reduction method to obtain a 4 × 4 classical beam model stiffness matrix of an isotropic circular cross-section using Nonuniform Rational B-spline for cross-sectional discretization. They also applied the method to composite box beams. Jiang and Yu 14 developed nonlinear cross-sectional analysis of a hyperelastic material under finite deformation using VABS. They studied structural behavior of a composite pipe consisting of two layers, and having cross section in the form of a box beam with rounded ends.
Previous work as seen from the literature for the stress analysis of composite tubular beams has focused on beams with fairly complex geometries sch as composite box beams. For such structure, the use of FE to mesh the complex geometry of the cross-section is essential. For circular sections, such as that of the cross tube used for helicopters or similar structures, the use of FE to mesh the cross-section may not be necessary. The use of a simpler method may facilitate simpler input, and better insight into the understanding of the behavior of the structure. This is the objective of the current paper. That is to use the procedure of the VAM so that dimensional reduction is exploited to reduce the computational effort, but not using the FE meshing of the VABS. Rather, the warping function for the deformation of the section is represented by a polynomial (Pascal polynomial) in polar coordinates. Since no meshing is required, it is termed a "meshless" method.
In order to introduce the Pascal polynomial method, due to the complexity of the procedure, some aspects of the VAM is represented here for clarity and completeness.

VAM PROCEDURE
The VAM consists of many steps. Details of the steps are presented below.

2.1
Step 1: Identification of small parameters and formulation for 3D strains In step 1, three aspects are considered. One is the identification of the small parameters essential for the application of VAM. The small parameters are the maximum of the axial strain max and the thickness-to-radius ratio (h/r o ) and the radius-to-length ratio (r o /L). 15 It is assumed that the h/r o ≪ 1 and r o /L ≪ 1. So that once can conclude that h/L ≪ 1 and the ratio between the small dimension and the large dimension (h/L) of the structure, where h can be the thickness or diameter of the section, and L is the length. The other aspect is the formulation to obtain relations between the 3D strains in terms of warping functions and the 1D strains, using vector mechanics. According to kinematics of the beam deformation presented by Yu et al. 4 based on the decomposition of the rotation tensor, the 3D strain field for a beam without initial twist and curvature can be written the following form.
where 1, 2, and 3 refer to the x 1 , x 2 , x 3 coordinate axes as shown in Figure 2, respectively. This 3D strain can be considered as follows. where where, [I] is the identity matrix and [O] is the null matrix. The 1D strains of the reference line can be expressed as the following where is strain vector of the reference line, which consists of the extensional strain 11 , the twist about the axis x 1 ( 1 ), bending curvature about the axis x 2 ( 2 ), and bending curvature about axis x 3 ( 3 ). w is the warping function displacements of the cross-section of the beam, which can be written in terms of its components as where w 1 represents the out-of-plane warping displacement at points on the cross-section in direction x 1 , w 2 , and w 3 are in-plane warping displacements of an arbitrary point in directions x 2 and x 3 , respectively. The calculated warping functions in terms of the 1D strains of the reference line of the beam results in the transformation of the 3D model into a 1D beam model which maintains the 3D effects. Using the 3D strain field above, the dimensional reduction process can be performed. w ′ the partial derivative of w with respect to x 1 as shown in Figure 1.

2.2
Step 2: Strain energy expression In step 2, 3D strain energy expressions are formulated, first in terms of the 3D strains, and then in terms of the 1D strain and curvature of the reference line. The strain energy of the beam U can be expressed as 6 F I G U R E 1 Beam and its coordinates where Ω is the cross-sectional area.
[D] is the matrix of material stiffness which carries the information of the material property and fiber orientation and layer orientation.

Step 3: Assumption for the warping field of the cross-section
In Equation (6), the unknowns are the warping functions w (x 1 , x 2 , x 3 ), and the strains of the reference line . The VABS approach to provide solution is to assume that the warping function is a product two entities, one is some function of x 2 and x 3 , and the other is function of only on x 1 . The warping field is discretized as where [W(x 2 , x 3 )] is the matrix of finite element shape functions, and V (x 1 ) as a column vector of the nodal values of the warping displacement over the cross-section in VABS method.

2.4
Step 4: Asymptotically correct refined theory According to the VAM, we now consider the perturbation in the warping field as follows. 16 where w 0 has the order of max and w 1 has the order of max h/L. We can write and, Therefore, the strain field in Equation (2) can be expressed as In the strain field of Equation (11), the first and second terms have the order of max , the third and fourth terms have the order of max h/L and the last term has the order of max (h/L) 2 . By considering Equations (10,11) and Equation (6), the strain energy per unit length can be expressed as follows.

Step 5: First and second approximation of the strain energy
The expression of the energy in Equation (12) has terms of different orders. By dropping terms of smaller magnitudes, different levels of solution can be obtained.

First approximation
The first approximation in the warping field in the VAM comes from keeping the terms up to order 2 max in the strain energy of Equation (12). Therefore, the order-zero of the linear strain density is as follows.
The strain energy in Equation (14) is minimized with respect to the unknown warping V 0 using the following relation.
This yields the classical beam model in which the effects of transverse shear are neglected. Substituting Equation (14) into Equation (15) leads to Therefore, the unknown warping coefficients according to the order-zero of the strain energy is obtained as It is difficult to invert the matrix [D hh ] in Equation (17). To find the inverse of this matrix, a mathematical Pseudo method 17 is used. Using the obtained column matrix of Equation (17) and the order-zero of the strain energy Equation (14), the first approximation of strain energy and stiffness constants of the classical beam model can be achieved as follows.
where S is a 4 × 4 stiffness matrix that expresses the beam cross-sectional stiffness coefficients. For a cross-section with area Ω, second moment of inertia I, polar moment of inertia J, elastic modulus E and shear modulus G, the stiffness constant S 11 or (EΩ) is extensional stiffness, S 22 or (GJ) is torsional stiffness, S 33 or (EI) x 2 and S 44 or (EI) x 3 are cross-sectional bending stiffness about x 2 and x 3 coordinates, respectively. Moreover, S 12 is the extension-twist coupling, S 13 and S 14 are the extension-bending coupling stiffness. In the case of isotropic material, the cross-sectional stiffness matrix can be in the form expressed in (A1).

Second approximation
The first approximation which leads to the classical beam model may not be accurate sufficiently for analysis of composite tubular beams with various lay-ups through the section. The second approximation of the strain energy with considering the transverse shear effect (Timoshenko refinement) may be accurate to predict the behavior of composite tubes. The strain energy up to order O( 2 max h 2 ∕L 2 ) can be written as It is noted that in the equation above the terms 2V T 0 [D hh ]V 1 and 2V T 1 [D h ] cancel out each other. 16 Also, in Equation (19), the term V T 0 [D hh ]V 0 is not mentioned in the strain energy obtained by Hodges. 3 However, this term is too low order. Integrating by parts of Equation (19) leads to Equation (20) as follows.
The energy functional Equation (20) can be minimized to result in the Euler-Lagrange equation with the same procedure as Equation (15) so that Using the Equations (21, 17), V 1 can be calculated as The strain energy functional Equation (19) can be obtained in terms of , ′ , and ′′ by substituting Equation (17 and 22) to into Equation (19). We have The use of Equation (23) includes derivatives of the 1D strains which is complicated. Therefore, in order to eliminate the derivatives of the 1D strains, the equilibrium equations are used to make a relationship between the strains and their derivatives. To do so, the transformation to the generalized Timoshenko theory is introduced in which the strain energy is converted to the generalized Timoshenko strain energy to take into account the shear strain effects. This procedure is carried out by deriving kinematic relations between asymptotically obtained 1D strains of Equation (4) and the generalized Timoshenko strains. The details of the procedure are expressed in References 3,7. It is noted that in the generalized Timoshenko beam model, the cross-section plane is not normal to the reference line in contrast to the classical beam model. The generalized Timoshenko strain energy is presented as ] T contains the Timoshenko beam model generalized strains.
The relation between and is where The 1D linear equilibrium equations are employed as ] T are force and moment resultants. F 1 is the sectional axial force, The relation between cross-sectional stress resultants and the 1D generalized strains can be written as Using the Equations (25, 26, 28, and 30), the derivatives of the 1D strains can be calculated in terms of and s . The details of the procedure are explained in Hodges. 3 The unknown stiffness matrices [X], [Y ], and [G] can be expressed as follows.
By calculating the stiffness matrices of Equation (31), one can rearrange the achieved stiffness constants to the form of strain energy of Equation (32) as follows.
In the strain energy of Equation (32) the effect of shear is considered where = [ 11 2 12 2 13 ] T . Moreover, the matrix S contains the Timoshenko cross-sectional stiffness constants. S 11 is the extensional stiffness, S 22 and S 33 are shear stiffness along x 2 and x 3 , respectively and S 44 is the torsional stiffness, S 55 (EI) x 2 and S 66 (EI) x 3 are bending stiffness about x 2 and x 3 , respectively. The off-diagonal elements are coupling stiffness constants, that is, S 14 is the extension-twist coupling and so on.

One-dimensional finite element beam analysis
The strain energy Equation (32) which is applied for 1D beam analysis, includes 3D deformation effects of the cross-section. Therefore, the 1D static analysis of the tubular beam can be performed using this strain energy equation. The 1D strain energy of beam element is as follows ] T and L e is the element length. We can calculate the stiffness of an element if we express the strains in Equation (33) in terms of degrees of freedom of beam element. For a Timoshenko beam model, we can use an element with three nodes which can prevent shear locking. 18 The element and its coordinate are shown in Figure 2. The element has 6 degrees of freedom at each node including three displacements (u, v, and w) and three cross-sectional rotations ( x 1 , x 2 , x 3 ) in the x 1 , x 2 , and x 3 , respectively. Hence, 18 degrees of freedom for the element is considered. In this model, the following polynomials for the approximation of the degree of freedom in the length is used. 18 For a Timoshenko beam model, the 1D strains can be written as It is noted that the x 2 and x 3 are the bending rotations of the beam section with respect to x 2 and x 3 , respectively, while 2 12 and 2 13 are rotation of beam section about x 3 and x 2 for transverse shear deformation, respectively. Using strain-displacement Equation (35) and Equation (34), the rotations can be expressed as The column matrix of degrees of freedom of the beam element will be as follows The unknown coefficients matrix is as the following The unknowns ( 1 − 18 ) in the Equation (34) are expressed in terms of the nodal displacement vector by substituting u, v, w, x 1 , x 2 , and x 3 at all three nodes of the beam element. we have The matrix [T] has an order of 18 × 18 includes the element nodal coordinates. The 1D strain vector can be expressed in terms of using Equation (32) as where matrix [N] has order of 6 × 18 which is a function of x 1 . Finally, the strain vector can be expressed in terms of nodal displacement vector using Equation (40) as The strain energy relation of the element can be written as The stiffness matrix of the beam element is The nodal force vector caused by external loads and moments can be expressed as where f is the vector of nodal force and f i

Step 8: Recovery analysis
The 3D strain distribution through the thickness of the tube can be obtained using Equations (2 and 26), and the use of obtained from Equation (A2),. 3

Novelty of the present work
In Equation (7), the warping functions are assumed to be a product of functions W(x 2 , x 3 ) and V(x 1 ). There are sections of complex geometries such as beams with C section as shown in Figure 3. To handle these complex geometries, Yu et al., 4 developed cross-sectional analysis of a composite box beams using VABS. The software performs FEA for cross-sectional discretization. It provides the equivalent stiffness for the section to be used in the subsequent 1D analysis (step 7) and recovery of the strains (step 8).
For composite beams of circular sections consisting of many layers such as that shown in Figure 4(A), the procedure to discretize the section using VABS is complex. The cross-section modeling procedure includes: (a) Defining two laminate sections as upper and lower laminate sections (the upside of the horizontal solid line is the upper laminate section and The modeling process of the circular cross-section in the VABS is time-consuming, particularly for parametric study. This is also not necessary. The analysis of composite circular cross-section can be more straightforward and faster using polynomial based as shown below.

Dimensional reduction using Pascal polynomial method
Note that in Equation (6), the first equality expresses the strain energy in terms of the 3D strains, and the second equality expresses the strain energy in terms of the warping functions w, its derivative with respect to the axis x 1 (w ′ ), and the 1D strains of the reference line. For the analysis of composite tube, a cylindrical coordinate system is introduced. Figure 5(A) shows the beam Cartesian coordinates (x 1 , x 2 , and x 3 ), material coordinates (e 1 , e 2 , and e 3 ) where e 1 is along the fiber direction, and the beam reference line. The position of an arbitrary point in the cross-section is shown by a radial position r and a circumferential angle in the polar coordinate system. represents the fiber orientation angle in its plane (winding angle). Figure 5(B) represents the cross-section of a single-layered tube with the inner and outer radius of r in and r out , respectively. The angle denotes the layer orientation to axis x 2 .
The matrix of stiffness coefficients [C] within this local coordinate system is given by Equation (A3). A relation between [C] and [D] is required. This is done through two coordinate transformations. In the first transformation, the axis e 1 is rotated by an angle α to coincide with x 1 . The transformation matrix [T ] is given as Equation (A4). At the end of this rotation, the axes e 2 , e 3 become e ′ 2 , e ′ 3 and these may not coincide with the axes x 2 , x 3 . The angle between the axis e ′ 2 and x 2 which takes into account the orientation of the layer is used. A rotation of the axis e ′ 2 by angle to coincide with axis x 2 is done through a rotation matrix [T ] as shown in Equation (A5). As such one has The relation between the circumferential angle and layer orientation can be introduced by a change of variable as = − 2 . In the discretization of warping field of Equation (7) used by VABS, the function [W(x 2 , x 3 )] consists of FE shape functions, and the vector V is a column matrix of the nodal values of the warping displacement over the cross-section. 3 In the proposed Pascal polynomial approach, w k (x 1 , x 2 , x 3 ) is proposed to be a simple polynomial as where k = 1, 2, 3. w k can be expressed in the vector form as The cross-sectional discretization in the polynomial method would be as follows.
where the vector V k consists of unknown Pascal coefficients a i which is a function of axial coordinate and should be determined. The Pascal polynomials as the warping functions are used for the whole section and the strain energy will be minimized with respect to a vector of coefficients V k (x 1 ). It should be pointed out that in VABS which applies FE, F I G U R E 5 (A) Coordinate systems of a composite tubular beam, (B) The layer orientation of (β) in the circumference of a one-layer tube with inner and outer radius r in and r out the warping at the common nodes at the interface between layers must set be to the same value, while in the Pascal polynomial method, the warping functions will be considered for the whole domain so that the problem of continuity of warping at the interfaces will be ensured. Using the polar coordinates, the transformation can be as (x 2 = r cos , x 3 = r sin ). For polynomial with order m (i.e., m = 4), we have W k = [1 r cos( ) r sin( ) r 2 cos 2 ( ) r 2 cos( )sin( ) r 2 sin 2 ( ) r 3 cos 3 ( ) r 3 cos( )sin 2 ( ) r 3 cos 2 ( ) sin( ) r 3 sin 3 ( ) r 4 cos 4 ( ) r 4 cos( )sin 3 ( ) r 4 cos 3 ( ) sin( ) r 4 cos 2 ( )sin 2 ( ) r 4 sin 4 ( )], and, The order of the polynomial can determine the accuracy of solution. A large number of layers will require higher polynomial order to assure accuracy. From Equation (4), w = [w 1 w 2 w 3 ] T , the matrix form for the warping function is where Ω i is the area of the ith layer. The use of VAM and VABS can greatly reduce the computational time and efforts for beams or tubes with complex cross-section. 4 For tubular composites with a large number of layers, the use of Pascal polynomial simplifies the input data and provides better insight into the problem. This is shown in the case studies below.

VALIDATION OF THE RESULTS AND DISCUSSION
In order to validate the present method of solution, a program in Maple 15 has been developed named AMTSM (Analytical Meshless Tubular Sectional Method). To do so, various tubular sections with different lay-ups are utilized as case studies. The parametric study to investigate the effect of polynomial order and lay-up sequence on the accuracy of the method is presented.

Aluminum tube
An aluminum 6061 T6 19 isotropic tubular section with the elastic modulus E = 69.7 GPa, Poisson ratio = 0.33, an outer diameter D out = 88.9 mm and the thickness h = 3.175 mm is studied. The obtained stiffness constants reported in Table 1

Composite tube with lay-up [20/ -70/20/(-70) 2 /20], 10
As a second case study, a thin CUS composite tubular section is studied and the stiffness constants are compared with those in Reference 10 which used an earlier version of VABS, Reference 12, and VABS 3.8. The obtained results are tabulated in Table 2 12 There are errors in the VABS (Timoshenko) results obtained by Popescu and Hodges 10 with respect to both present (Timoshenko) and Rehfield et al. 12 This error was modified in the higher versions of VABS proposed by Yu et al. 4 To this end, we analyzed this tube using VABS 3.8 with 7588 elements and the obtained results are reported. It can be seen that there is a large difference in the obtained stiffness constants between Timoshenko solution and Classical Solution, which can be due to the higher-order approximation of energy in Timoshenko and considering the shear effect. There is a sign difference in the calculated S 14 and S 36 of Rehfield et al. 12    The tube has inner and outer diameters D in = 56 mm and D out = 98 mm, respectively. The tube cross-section is meshed using VABS 3.8 with 9098 elements as shown in Figure 7. The obtained stiffness constants are reported in Table 3. There is a large difference in the obtained flexural stiffness (S 55 and S 66 ) between Cl. and Timosheko models. In addition, the obtained stiffness constants are compared with those of the 3D elasticity solution of Lekhnitskii for multilayered composite tubes proposed by Jolicoeur and Cardou 20

Effect of polynomial order (m), and lay-up sequence, on the accuracy of the stiffness constants of the cross-section
The selected order of Pascal polynomial in the warping field, and lay-up sequence may affect the accuracy of not only the stiffness constants of the tube but also on the 3D strains of the tube. To this end, we examined composite tubes with material properties E 1 = 140 GPa, E 2 = E 3 = 10 GPa, v 12 = v 13 = 0.31, v 23 = 0.33 and G 12 = G 13 = G 23 = 5.53 GPa. The outer and inner radiuses of the tubes are r out = 59 mm and r out = 39 mm, respectively. The thickness of each layer is assumed to be 2 mm. Note that all the following reported stiffness constants are in SI.

Unidirectional composite tube [0 10 ]
The effect of polynomial order on the accuracy of stiffness constants of a unidirectional tube is investigated in Table 4.
Values are compared with those obtained from VABS 3.8. Except for S 22 of (m = 2), the other stiffness constants are close to which obtained by VABS 3.8. The obtained S 22 and S 33 of orders (m = 5 and 7) are more accurate compared to other orders. It is noted that both the obtained Timoshenko and Classical are the same.

Composite tubes with lay-up [90 5 /0 5 ]
The effect of polynomial order (m) on the accuracy of the obtained flexural stiffness of cross-ply tube [90 5 /0 5 ], is demonstrated in Table 5. It is observed that the obtained Classical and Timoshenko flexural stiffness of the tube are similar when polynomial order increases. Increasing the polynomial order from 4 to 7 does not have a significant effect on the accuracy of the predicted flexural stiffness for tubes. Also, both the Timoshenko and Classical common stiffness constants are the same. The obtained results converge for [90 5 /0 5 ] with order (m = 7) of the polynomial.

Composite tubes with lay-up [±60 5 ]
The influence of polynomial order on the nonzero stiffness constants of the tube [±60 5 ] is investigated, as shown in Table 6. It is observed that for order (m = 2) of the polynomial, except for S 14 , there is no significant difference between the obtained (Cl.) results and those obtained by VABS 3.8. By increasing the order of the polynomial, the values of S 14 become close to VABS 3.8. The two stiffness constants S 22 and S 33 of the order (m = 2) have a large difference compared with VABS 3.8. As the polynomial order increases, the difference decreases noticeably. A very good agreement between the results of the present approach with (m = 7) compared to VABS 3.8 is observed. In addition, it can be observed that the Classical and Timoshenko stiffness constants for this type of lay-up sequence [± 5 ] are the same.

Composite tubes with lay-up [45 5 /15 5 ]
To further investigate the effect of lay-up sequence on the accuracy of the present approach, we evaluate flexural stiffness of angle-ply [45 5 /15 5 ] tube. Table 7 shows the effect of polynomial order on the stiffness of tubes [45 5 /15 5 ]. Also, it is seen that increasing Pascal polynomial order, increases the accuracy of stiffness prediction. There is a large difference between the obtained flexural stiffness (S 55 and S 66 ) of the Classical and Timoshenko for both the methods. Which means that combination of angle-ply layers causes difference between the Classical and Timoshenko results. This difference between Classical and Timoshenko is due to the effect of shear strains which are considered in Timoshenko beam model. It is observed that the other common stiffness constants of the Classical and Timoshenko are the same. The percentage difference of the m = 2 values compared to VABS 3.8 is large. By increasing the polynomial order to (m = 5 or 7), the percentage difference decreases. The maximum of the percentage difference between the present results and VABS 3.8 is 7.690 for S 22, whereas the minimum of the percentage difference is 0.485 for flexural stiffness (S 66 ).

Distribution of 3D strains of [±60 5 ] tube
The distribution of longitudinal strain Γ 11 , Γ 22 (normal strain along x 2 ), Γ 33 (normal strain along x 3 ) and out-of-plane shear strain Γ 12 through the thickness at = 2 of the tube [±60 5 ] with the material properties and dimensions of the Section 5.2 under M 2 = 1 KN.m are demonstrated in Figure 13. The order (m = 7) of the polynomial is used for the analysis. An excellent correlation between the present method and VABS 3.8 can be observed. It is noted that the tube is meshed with 97,800 elements in ANSYS 3D and the cross-section is meshed with 3740 2D elements in VABS 3.8.

Validation of present 1D FE analysis
The deflection of the tubular composite beam calculated using the present polynomial (Section 4) and 1D FE (Section 3.    In the ANSYS modeling, beam 189 element type was selected which has three nodes and 3 degrees of freedom at each node. The boundary conditions are considered as (u = v = w = 0, x 1 = x 3 = 0) at two ends of the beam. The distribution of transverse displacement in the length of the tube is demonstrated in Figure 14, with the experimental point at mid length of the tube. The convergence study shows that the obtained displacement converges with 40 number of 1D elements. It can be seen that a good agreement is obtained.

Parametric study
In Figures 15-17, we investigated the influence of the fiber angle α on the change of stiffness constants S 11 , S 14 , S 36 , S 33 , S 55 , and S 66 of the composite tubes with lay-up [± 5 ]. The dimensions and material properties are similar to Section 5.2. It is seen from Figure 15 that as the fiber angle increases, the extensional stiffness (S 11 ) decreases. The change of S 11 is considerable until the fiber angle of = 45 • . The change of extensional-torsional (S 14 ) and torsional-bending (S 36 ) stiffness versus the fiber angle are depicted in Figure 16. At the fiber angle = 30 • , the maximum of the S 14 Figure 17 show that the torsional stiffness (S 33 ) is maximum at = 45 • and the bending stiffness (S 55 and S 66 ) decreases by increasing the fiber angle.

CONCLUSION
The conclusion can be made along two different aspects: The accuracy of the simple Pascal polynomial method, and the effect of the lay-up sequence on the need to use refined theory, as follows • It was shown that the simple Pascal polynomial method can provide results that are as accurate as the VABS method and more computationally efficient than 3D FEA. This has been demonstrated for the determination of the stiffness constants, the through-section strains, and the deflection of the beam. This has also been demonstrated for tubes of different lay-up sequences, from simple lay-up sequence such as cross-plies to complicated lay-up sequences with different fiber angle orientations. For example, Section 5.3.1 shows the present cross-sectional analysis can substantially reduce the computational cost compared to 3D FEM. As such, the present method would be beneficial for structural optimization • For tubes made of cross-ply layers, that is, [0 m /90 n ], there is no significant difference between the stiffness constants obtained using Classical model and Timoshenko model. This is due to the absence of shear coupling terms.
• In order to analyze cross-ply tubes, it is not necessary to use higher degree of the polynomial.
• For tubes with the lay-up sequence [± 5 ], there is no difference between the obtained flexural stiffness of Classical and Timoshenko models.
• For tubes with lay-up sequence made of layers of different fiber orientation such as [45/15], there is a significant difference in the flexural stiffness of Timoshenko and Classical models obtained by both the Pascal polynomial method and VABS 3.8.
• To achieve higher precision of 3D strains for tubes with combination of angle ply layers, the order of applied polynomial needs to be increased, (m ≥ 5).
• The extensional stiffness (S 11 ) decreases by increasing the fiber angle of a composite tube.
• The magnitude of extensional-torsional (S 14 ) and torsional-bending (S 36 ) of a composite tube are maximum at = 30 • of the fiber angle.
• The torsional stiffness (S 33 ) of a composite tube is maximum at = 45 • of the fiber angle.
• As the fiber angle increases, the bending stiffness (S 55 and S 66 ) of a composite tube decreases.
• The present 1D FE formulation agrees well with 1D FE of ANSYS.