Accurate Stress Analysis of Variable Angle Tow Shells by High-Order Equivalent-Single-Layer and Layer-Wise Finite Element Models

New concepts of lightweight components are conceived nowadays thanks to the advances in the manufacture of composite structures. For instance, mature technologies such as Automatic Fibre Placement (AFP) are employed in the fabrication of structural parts where fibres are steered along curvilinear paths, namely variable angle tow (VAT), which can enhance the mechanical performance and alleviate the structural weight. This is of utmost importance in the aerospace field, where weight savings are one of the main goals. For that reason, shell structures are commonly found in the aerospace industry because of their capabilities of supporting external loadings. Straight-fibre composite shell structures have been studied in recent decades and, now, spatially varying composite shells are attracting the attention of manufacturers. This work analyses the mechanical behaviour of VAT composite shells subjected to different external loadings and boundary conditions. The Carrera Unified Formulation (CUF) is employed to obtain the different structural models in a systematic and hierarchic manner. The outcomes of such numerical models are discussed and compared with commercial software Abaqus.


Introduction
Nowadays, an increasing number of engineering applications, such as civil, marine, automotive and aerospace, are employing thin-walled structures-as an example, bridges and oil rings, ships, chassis and aircraft, respectively. Specifically, shells consist of curved lightweight constructions, which became very widespread in structural engineering thanks to their high performance when supporting external loads. Such outstanding mechanical properties stem from the coupling between the membrane and flexural behaviour, induced by the curvature. The geometric characteristics of shell structures, including the initial curvatures, have a direct influence on the stiffness properties [1]. Therefore, a proper design of such structures is crucial to perform accurate stress predictions under different loading and boundary conditions. The popularity of shell models is thanks to their lower computation cost compared to three-dimensional (3D) models.
As stated by Kapania [2], composite shell structures are playing a crucial role in several branches of engineering. Particularly in aerospace, composite materials offer an attractive possibility to more traditional types of construction due to their corrosion resistance, high strength-to-weight ratio, ease of formability, excellent fatigue resistance and tailoring ability. In fact, the stiffness, strength and flexibility characteristics can be controlled in different directions allowing one to obtain a significant increase in performance. Furthermore, the possibility of applying fibres that are not constrained along straight paths but can follow curvilinear paths leads to numerous advantages from the point of view of structural efficiency, increasing the design space and improving the tailoring process [3][4][5]. These new advanced composite structures, called variable angle tow (VAT), have attracted a lot of interest because they allow the manufacturing of variable stiffness composite laminates (VSCL) without discontinuity in the material while maximising the ratio of stiffness to mass.
The VAT methodology is not new, but it has recently received renewed interest due to significant technological improvements in the automatic manufacturing process. Essentially, the advanced Automated Fibre Placement (AFP) technique [6] and the Continuous Tow Shearing (CTS) process [7] allow the fibre orientation angle of a layer to vary with respect to one or more spatial directions. The literature concerning the VAT theories is vast. As an example, several progressive damages and failure analyses on classical and VAT composite panels were performed by Lopes et al. [8], showing the potentialities of the variable stiffness structures to redirect the load fluxes to the stiffer edge area to improve the structural performance. Curvilinear fibres were adopted by Hyer and Lee [9] to change the stress concentration around a hole. Stodieck et al. [10] investigated the possibility of improving the aeroelastic tailoring of a rectangular unswept laminate wing using the VAT methodology with the comparison with the classical one. Setoodeh et al. [11] adopted the VAT technique to maximise the buckling load of composite panels.
Over the years, researchers and scientists have developed several efficient shell theories. For example, the studies of Poisson, [12], Love [13], Mindlin [14], Kirchhoff [15], Reissner [16] and Cauchy [17] represent the classical formulation of the shell models. Typically, these classical theories are adopted in the commercial codes. However, the applicability of classical theories is limited to a narrow range of applications, for example, when dealing with the thin-walled structure and without local effects. On the contrary, more accurate shell formulations are needed when transverse stresses analyses are required. Recently, several higher-order 2D formulations were formulated to improve the accuracy of classical theories. For example, Reddy [18] provided a simple high-order theory for laminated composite two-dimensional (2D) structures. A shear model considering a parabolic distribution of transverse shear deformations in the thickness direction was presented by Reddy and Liu [19]. An assessment of the relevance of displacement variables in refined theories for isotropic and multilayered shells using an axiomatic/asymptotic technique was carried out by Mashat et al. [20]. Carrera [21,22] developed several refined shell theories in the framework of the Carrera Unified Formulation (CUF). Cinefra and Carrera [23] performed linear analyses of composite cylindrical structures using finite shell elements with different through-the-thickness kinematics. A useful review of methods and guidelines for the choice of the shell model was reported by Petrolo and Carrera [24]. A detailed review of the theories is not within the scope of this work. Readers are referred to [25,26] for other significant works on refined shell formulations.
There exists a limited number of works regarding the stress evaluation of VAT shells. To the best of the authors' knowledge, only a couple of works have been published about this topic, specifically, the works by Tornabene et al. [27] and Sarvestani et al. [28]. The former [27] used a similar approach to the one proposed in this paper, whilst the latter [28] adopted a semi-analytical methodology to perform hygro-thermo-mechanical analysis on thin to relatively-thick VAT composite panels. To circumvent that research absence, this paper aims to present additional stress benchmarks for future comparisons. In this manner, the main objective of this manuscript is to accurately predict the 3D linear stress state of tow-steered composite shells and provide stress benchmarks.
The proposed methodology relies on the CUF [29]. In CUF, the accuracy of the model can be fine-tuned straightforwardly since the order of the structural model is treated as input of the analysis. CUF has been adopted to obtain 2D theories, which have been then extended to the analysis of multilayered, composite plates and shells, and later for beam models [30][31][32]. In recent years, CUF has also been employed in the analysis of VAT composites. For instance, Demasi et al. [33] showed some numerical assessments on the stress distribution of VAT plates employing Equivalent Single Layer (ESL) and Layerwise (LW) theories. The vibrations and buckling of variable stiffness plates were studied by Vescovini and Dozio [34] by coupling the CUF and Ritz method. Viglietti et al. [35] introduced 1D elements for the free vibration study of VAT laminates. Lately, Pagani and Sanchez-Majano [36,37] and Sanchez-Majano et al. [38] analysed the influence of manufacturing defects on the mechanical performance of VAT composites.
This manuscript is subdivided as follows: (i) first, a description of the 2D CUF modelling approach for shells is made in Section 2, including a description of VAT over curved domains; then, Section 3 presents the numerical results obtained with CUF shell models. Finally, conclusions are drawn in Section 4.

Preliminary Considerations
In this work, the structures are modelled using refined shell models. A shell is a 2D structural element where the thickness is negligible compared to the other dimensions. Typically, this geometry is described employing an orthogonal curvilinear reference system (α, β, z), as reported in Figure 1, in which α and β indicate the in-plane surface and z the thickness direction. For the sake of brevity, the complete description of the shell formulation in the CUF domain is not within the scope of this article. The reader is referred to [32,39]. The transposed displacement, strain and stress vectors for each layer k are written as follows: The displacement-strain relations are written as: in which D represents the differential operator containing the geometrical relations between strains and displacements. This operator reads as: where ∂ α = ∂(·)/∂ α , ∂ β = ∂(·)/∂ β , ∂ z = ∂(·)/∂ z , and H α , H β are defined as: where R k α and R k β denote the radii of the middle surface of the k th layer, and A k and B k indicate the Lamé parameters. Using the constitutive equations, stresses are evaluated as: where C k is the material elastic matrix [40,41]. Since in this work VAT structures are investigated, the fibres have a general orientation function of the space coordinates, i.e., θ(α, β). Thus, we write: in which: where T represents the rotation matrix [42]. The matrix C k changes pointwise in VAT composite structures. In a VAT structure, the fibre can change continuously along a curvilinear path in each ply. In this way, the laminate has a different stiffness value at each position. A linear fibre angle variation over the lamina is employed in this work, see Figure 1, and the fibre orientation, using the notation of Gürdal [43], is formulated as follows: where the fibre path exhibits a rotation of an angle Φ with respect to a certain reference direction. The fibre orientation angle at this point is T 0 and varies along a direction α oriented by angle Φ from the original coordinate axis α. The fibre orientation reaches the value T 1 at a characteristic distance d from the reference point. By considering this rotation angle, the fibre orientation path θ(α, β) is expressed as θ(α ), in which α = αcosΦ + βsinΦ. The parameter d is equal to a/2 or b/2 when Φ = 0°or Φ = 90°, where a and b are the width and length of the 2D structure, respectively. For a better understanding of the fibre path variation along the in-plane, two of the fibre paths considered in the upcoming assessments are represented in Figure 2. For clarity, unlike commercial codes where the lamination angle is considered constant over the entire element, in the presented methodology, the material coefficients of the VAT structure are evaluated in specific Gauss points [35]. The use of the Gauss point integration technique guarantees a more accurate and efficient analysis of composite VAT structures, since the variability in the material stiffness coefficients C ijkl is accounted in several points for the same FE. Moreover, in the present model, the number of Gauss points per element can be increased independently of the FE type, and thus, it influences the number of degrees of freedom.

Kinematic Assumption, Governing Equation and FE Approximation
In this article, VAT composite shells are modelled employing refined 2D CUF models. In the CUF framework, the refinement of the theory is assumed as an input of the analysis, so low-to higher-order models can be built with ease and in a unified manner (i.e., no adhoc formulations are needed to obtain any model). The 3D displacement field is formulated as an arbitrary through-the-thickness expansion of the in-plane variables.
in which F τ represents a set of thickness expansion functions, u τ is the generalised displacement vector depending on the in-plane coordinates α and β, M indicates the order of expansion in the thickness direction and the repeated index τ denotes summation. The choice of F τ is arbitrary and determines the class of the considered 2D CUF shell model. For the sake of brevity, the reader is referred to [29] for a detailed description of the shell theories within the CUF domain. The use of laminated composite structures leads to important challenges in the design process. One of the most important assessments is undoubtedly the correct prediction of the stress distribution within the structure. In the literature, ESL and LW theories are typically employed when dealing with composite materials. In this work, both the ESL based on Taylor Expansion (TE) and LW adopting the Lagrange Expansion (LE) are considered. The acronym LDN, used in the following figures, denotes the LE of order N assumed in the z direction. The differences in the assembly procedure for ESL and LW and the behaviour of the primary variables along the thickness of the shell are reported in  In detail, in ESL the stiffness matrix is evaluated with the homogenisation technique of the properties of each layer by summing the contributions of each layer. The result is a multilayer configuration modelled as a single layer having a set of variables assumed for the entire cross-section. On the contrary, LW theories divide and expand the displacement field within each material layer. By doing so, homogenisation is carried out at the interface layer. ESL theories exhibit accurate results of the global response (fundamental vibration frequency, transverse deflection), but they are often inaccurate for the 3D stress distributions compared to the LW methodology.
Independently of the selected shell model kinematics, the finite element method (FEM) is used to approximate the in-plane generalised displacement vector employing the Lagrange shape functions N i (α, β).
in which q τi denotes the unknown nodal variables, N n indicates the number of nodes per element and i stands for summation. For completeness, the classical 2D nine-node quadratic finite element (FE) (Q9) is considered in the following analyses for the shape function in the α, β plane. The principle of virtual displacements (PVD) is used to derive the expression for the stiffness matrix. PVD states that: in which δL int represents the virtual variation of the internal strain energy (12) and δL ext is the virtual work of external loading where p sj is the vector of the applied point load components (3 × 1) and f sj is a surface force. By using Equations (10) and (12), the constitutive law (Equation (6)) and the geometrical relations which results in the following equation for the stiffness matrix: where k ijτs is the 3 × 3 fundamental nucleus (FN), see the CUF book [29] for its derivation. The mathematical expression for FN results in: Conversely, from other CUF-based works, Equation (15) cannot be split into separate integrals where the FE solution and CUF expansion are independently evaluated, and, therefore, a 3D integration is needed. Finally, the global assembled stiffness matrix K is obtained by looping through the indices i, j, τ, s.

Numerical Results
The stress analyses of VAT shell structures subjected to different loadings are discussed in this section. In particular, flat and curved panels with different material properties, lamination schemes, curvature ratios and boundary conditions were investigated. For this purpose, both the layerwise theory and equivalent single layer approaches are adopted and compared, showing the need to adopt layerwise when evaluating transverse shear stresses. The material properties considered in the following studies are reported in Table 1. These materials were selected from existing literature concerning shell-like structures [32,44].

Simply Supported VAT Flat Panel
The first numerical assessment consists of a laminated squared flat panel composed of two layers. This flat panel is simply supported on its four edges, and the stacking sequence is: θ = [0 < 90, 45 >, 0 < 0, 45 >]. A graphite/epoxy composite, whose elastic properties are available in Table 1, is employed in this structure. A graphical representation of the flat panel subjected to a uniform pressure (P z = 10 kPa) exerted on the top surface, including its dimensions, is illustrated in Figure 4. First, a convergence study on the in-plane finite element mesh is carried out. The reference results were obtained using commercial software Abaqus [45], where a mesh employing 80 × 80 × 16 solid C3D8R elements was utilised. Figure 5 shows the transverse normal (σ zz ) and (σ yz ) shear stresses for different 2D plate models, from 64 Q9 FEs up to 196 Q9 FEs, whereas 2LD3 elements were placed through the thickness direction. Table 2 gathers the six stress components evaluated at point Q (x = 0.25 m, y = 0.25 m) and z = 0.02 m for the different mesh approximations. These results suggest that a 10 × 10 Q9 mesh approximation provides an accurate evaluation of the stresses. However, by looking at Figure 5b, the lower part of the stress distribution does not match properly with the reference results. A better prediction is provided in this case by the 14 × 14 mesh approximation. Then, it is appreciated in both Figure 5a and Table 2 that the present approach does not predict the transverse shear stress distribution of the commercial software. These differences are mentioned below.   To perform an accurate stress prediction, different expansion functions are used for the thickness discretisation. Both LE and TE functions are considered in this study. The through-the-thickness stress distribution is shown in Figure 6. The corresponding stress values are enlisted in Table 3 for different 2D theories. As appreciated, the LW approach provides the more accurate stress results when 2LD3 theories are employed. The outcomes suggest that the ESL model is sufficient to evaluate the in-plane normal and shear stresses, whilst their accuracy diminishes when predicting the transverse shear stresses. Concerning the latter, TE 6, TE 7 and TE 10 provide a similar distribution compared to 2LD3 but present equal or higher DOF. However, oscillations appear close to the interface for those two theories. This means that even considering up to tenth-order terms, an accurate evaluation of these stress components is not available. Additionally, there exist differences between the Abaqus model and the 2LD3 relative to shear transverse stresses, especially for the σ yz component. These are due to the different formulations that both models present: the dedicated Abaqus model uses a solid 3D linear model (C3D8R element), whereas a LW formulation is achieved with the present approach. Indeed, Abaqus does not respect the stress-free condition at the plate's bottom and top as shown in Figure 6e.

Clamped VAT Curved Panel
The second numerical assessment considers a curved VAT panel. The structure has one meter length, an opening angle ϕ = 0.2 rad, internal radius R α = 1.25 m and thickness h = 0.05 m. The panel is clamped on its longitudinal edges and an uniform pressure p z = 10 kPa is exerted on the upper surface, as depicted in Figure 7. The stacking sequence is chosen to be θ = [0 < 0, 50 >, 90 < 0, 75 >, 45 < 0, 15 >] s . The material properties are gathered in Table 1. First, a mesh convergence analysis is performed by varying the number of in-plane finite elements, while a 6LD2 expansion theory is employed in the thickness direction. Tables 4 and 5 show the convergence of the transverse displacement and the stresses computed at points T and V, respectively. Figure 8 shows the stress distribution at point V using the present approach, whilst the convergence of the Abaqus model is available in Figure 9. Small differences between FE meshes are appreciated when accounting for the transverse displacement and the in-plane stresses σ αα and σ αβ . However, these differences are more evident for the transverse stresses σ zz and σ βz . From these results, it is inferred that a 20 × 10Q9 + 6LD2 mesh approximation thoroughly predicts the stress distribution.  In the following, stress distributions provided by high-order LE and TE expansion theories are compared with commercial software Abaqus. For such purpose, the 20 × 10Q9 FE mesh approximation is employed. The Abaqus solid model utilises 80 × 40 × 18 C3D20R quadratic elements. Table 6 and Figure 10 provide the stress components computed at point V of the structure. From the former, it is inferred that TE theories are unable to predict the shear stress σ βz , while for the remaining stress components, these expansions are in agreement with 6LD2 and 6LD3. The latter demonstrates that 6LD1 is not capable of retrieving response for any stress component. It is also shown that TE models provide an accurate evaluation of the in-plane terms while presenting difficulties when predicting the transverse components, especially the shear stresses σ αz and σ βz . The Abaqus prediction of the transverse shear stresses are not reported in Figure 10c-e since it is not feasible to obtain accurate results for such stress components. Table 6. Comparison of the in-plane and out-of-plane stresses provided by different expansion theories for the [0 < 0, 50 >, 90 < 0, 75 >, 45 < 0, 15 >] s clamped VAT shell. Stresses are computed at point V and z = 0 using the 20 × 10 Q9 mesh approximation.

Hinged VAT Shell
The third numerical assessment consists of a hinged VAT shell structure. This curved panel comprises three composite laminae where the fibres are steered following the stacking sequence θ = [90 < 30, 0 >, 0 < 30, 0 >, 90 < 30, 0 >]. A graphical representation of the VAT shell is available in Figure 11, and the material properties are enlisted in Table 1.
Specifically, a convergence analysis is first conducted to find the mesh that provides the most accurate results; then, a comparison between ESL and LW models is carried out. Finally, the effect of the fibre orientation parameters T 0 and T 1 on the stress distribution is addressed.

Mesh Convergence Analysis
A convergence study concerning the number of in-plane finite elements is conducted. Likewise, the through-the-thickness expansion order is addressed. Tables 7 and 8 show the convergence results for both the transverse displacement and the stress distribution at point R and z = 6.35 mm, and point S, respectively. Then, in Figure 12, it is appreciated that a 32 × 32 Q9 mesh approximation is necessary to retrieve an accurate evaluation of both the transverse displacement and the stress distributions. The effect of the different expansion functions on the latter is also accounted for in the right panel of Figure 12, especially in the out-of-plane stress variation. Indeed, it is demonstrated that a 6LD2 expansion is required to predict accurately such distributions.

ESL and LW Theories for the Analysis of Laminated VAT Shells
Then, a comparison between TE and LE expansion functions is carried out. Table 9 enlists the stress values at point R and z = 0 mm for the different expansion theories. Subsequently, Figure 13 shows the stress variation provided by these expansions. It is inferred that TE 2, TE 3 and TE 4 theories are able to predict the same in-plane stress distribution that high-order LD theories provide with a fraction of DOF. However, differences appear in the evaluation of the transverse stress components. For instance, TE 1 is not enough to catch the transverse normal stress distribution. Then, as appreciated in Figure 13d,e, TE 1 and TE 2 do not predict the trend of the transverse shear stresses. Moreover, TE 3 and TE 4 present difficulties when retrieving σ αz , as depicted in Figure 13d. Table 9. Comparison of the in-plane and out-of-plane stresses of the hinged VAT shell provided by the different expansion theories. Stresses are computed at point S and z = 0 mm using the 32 × 32 Q9 mesh.

Effect of the Fibre Orientations on the Stress Distribution
First, the effect of T 0 on the stress distribution, evaluated at point R, is addressed. Precisely, T 0 takes the following values T 0 = [10, 30,50,70] • . For these analyses, the 32 × 32 Q9 mesh approximation coupled with the 6LD3 expansion theory in the thickness direction is employed. Figure 14 illustrates the influence of T 0 on the stress distribution. On the one hand, it is appreciated that the magnitude of the transverse stresses is affected by T 0 . Nevertheless, a sign variation of the stress component, due to T 0 , is not foreseen. On the other hand, concerning the in-plane stresses, T 0 varies its magnitude and the stress sign. That is, in σ ββ , the second layer is subjected to compressive stresses when T 0 = [10, 30,50] • and to tension stresses when T 0 = 70 • . To study how the variation of T 1 influences the stresses over the VAT shell, the initial [90 < 30, 0 >, 0 < 30, 0 >, 90 < 30, 0 >] stacking sequence becomes [90 < 30, T 1 >, 0 < 30, T 1 >, 90 < 30, T 1 >], where T 1 = [0, 10, 50, 70] • . The stress components are evaluated at point W, where the effect of T 1 is appreciated. Figure 15 illustrates the through-the-thickness stress distribution. It is observed that the transverse stresses are not greatly affected by the fibre angle variation. On the contrary, the in-plane components show evident changes, especially the normal stresses where differences are more obvious in the middle layer.

Concluding Remarks
This paper has presented high-order finite elements to accurately evaluate linear static stresses over variable angle tow (VAT) shell structures. The Carrera Unified Formulation (CUF) constitutes the framework in which such numerical models are embedded because of its capabilities for achieving low fidelity to high-order structural models in a systematic manner. No approximations have been made on the shell geometry nor on the strains. Several geometries, boundary and loading conditions, spatially varying fibre paths and materials have been analysed and compared to commercial software Abaqus. The results suggest that: • Equivalent Single Layer (ESL) and Layerwise (LW) are in good agreement with commercial software formulations for the analysis of flat panels; • In the case of curved panels, Abaqus and ESL provide similar results for the in-plane stress prediction when compared to LW models. Nevertheless, obvious differences arise in the evaluation of the transverse stresses. Such differences can be mitigated by employing high-order TE; • The selection of the fibre path parameters (φ, T 0 , T 1 ) can be fine-tuned to conveniently alter the stress distribution at certain points of interest of the structure.
Future works will concern the linear and nonlinear behaviour of sandwich-like VAT shells and applications for structural failures, such as postbuckling.