Efficient three-dimensional geometrically nonlinear analysis of variable stiffness composite beams using strong Unified Formulation

https://doi.org/10.1016/j.tws.2021.107672 Received 15 December 2020; Received in revised form 8 March 2021; Accepted 10 Available online 29 March 2021 0263-8231/© 2021 The Authors. Published by Elsevier Ltd. This is an open access (http://creativecommons.org/licenses/by/4.0/). March 2021 article under the CC BY license S.O. Ojo, G. Zucco and P.M. Weaver Thin-Walled Structures 163 (2021) 107672 S E S Fig. 1. Cartesian frame of reference for (a) Laminate stack, (b) constant stiffness orientation, and (c) variable stiffness orientation. numerical method to evaluate interlaminar stresses in VS composites. Similarly, a p-version finite element (FE) model was proposed in [26,27] to investigate the natural modes of vibration, stresses and nonlinear bending deflection of VS laminates using a Reddy-type thirdorder shear deformation theory, although Groh and Weaver [28] later showed that the Reddy-type model leads to static inconsistencies at clamped boundaries. Furthermore, Soriano and Díaz [29] studied the failure process of VS composites using 3D FE analysis combined with continuum damage mechanics. Moreover, to achieve computationally efficient outcomes in the analysis of VS laminates, Demasi et al. [30] adopted equivalent single layer, layerwise and zigzag models based on FE methods. Despite the improved efficiency over 3D FE models for some simple laminate stacking sequences studied in [30], the potential of these models for analysis of complex and arbitrary lay-ups remains unclear. In addition, FE-based models often require fine meshes which may be computationally expensive to extensively exploit the merits of VS designs. To this end, high-order methods which can guarantee high levels of accuracy combined with faster solution than FE methods become desirable. The Differential Quadrature Method (DQM) has been widely applied for analysis of CS composite structures, see [31–33] for example, due to its promising characteristics of spectral accuracy and fast computation. In the context of VS applications, the Generalised Differential Quadrature (GDQ) has been used to study free vibration and static response of doubly-curved shell structures reinforced by curvilinear fibres in [34,35] while Groh and Weaver applied DQM for accurate stress analysis of VS laminates based on a third-order zig-zag theory within a Hellinger-Reissner mixed variational framework [36]. Furthermore, Luan et al. [37] employed a fifth-order zigzag theory within a HellingerReissner mixed formulation and a mixed inverse DQM approach for accurate stress analysis of VS laminates under different loading conditions. Although DQM implementation of the zigzag formulation for analysis of VS composites leads to significant computational savings over 3D FE solutions, further studies by Patni et al. [38] show that high fidelity models with enriched kinematics are required to resolve the discrepancies in the prediction of transverse stresses which affect delamination. To this end, a high-order serendipity-Lagrange based 2 FE model within Unified Formulation (UF) framework was proposed in [38], which yields accurate predictions for linear static analysis of VS laminates. To a large extent, the theory of UF has been employed for a variety of mechanics problems ranging from linear to nonlinear analyses (see [39–41]), due to the generic qualities of the kinematics adopted that allow for arbitrary expansion of displacement variables, leading to accurate prediction of local and global response of composite structures. With regards to geometrically nonlinear analysis, Pagani and Carrera [40,41] demonstrated that UF implemented within a FE framework is able to capture large displacements/rotations, in-plane deformation, localised buckling, cross-sectional warping, bending–torsion coupling and other high-order phenomena in thin and thick composite structures. Consistent with the benefits of UF for geometrically nonlinear analysis by Pagani and Carrera, Ojo and Weaver [42] proposed a DQMbased geometrically nonlinear strong UF (SUF) model to extend the capabilities of the linear strong UF model proposed in [31] for analysis of CS composites. Outcomes show significant improvement in terms of accuracy and efficiency over the FE-based model developed in [41]. In a bid to explore the potential of the newly developed geometrically nonlinear SUF, this study extends the capability of SUF to investigate the structural response of VS laminates undergoing large deflections. Considering the merits of VS designs for a variety of structural applications, it is important to understand the role of VS properties in achieving improved performance especially in nonlinear, large deflection applications. With this aim in mind, SUF is employed to investigate the effect of tow-steering on the behaviour of composite laminates under large bending and compressive loads. In addition, since efficient prediction of nonlinear 3D stresses is an important requirement to fully characterise large deflections in engineering structures, the potential of SUF for computationally efficient nonlinear analysis of composite structures is assessed comparatively with the models in the literature and ABAQUS model. The rest of the paper is arranged as follows. The basic formulation for geometrically nonlinear SUF and subsequent discretisation by DQM are presented in Sections 2 and 2.1. Derivations of the constitutive and geometric components of the tangent stiffness matrix are presented in Section 2.2 while Section 3 is dedicated to description of the arc-length routine employed for solution of the nonlinear system. In Section 4, settings for numerical test are explained while mathematical description of the constitutive relations of variable fibre architecture is also outlined. Numerical results and discussions are further presented in Section 4 before the final concluding remarks are detailed in Section 5. 2. Nonlinear strong unified formulation According to the theory of Unified Formulation (UF), the displacement components [ ux, uy, uz ] in the three dimensions, x, y and z axes of the beam-type structure of length L in the cartesian coordinate system are represented by the relation (see Fig. 1), u(x, y, z) = Fτ (x, z)uτ (y), (1) where Fτ is a function that captures the cross-sectional deformation and can be expanded to any order τ for enrichment of the beam’s kinematical description. The Serendipity Lagrange expansion (SLE) introduced in [43] idealises the beam’s cross-sectional behaviour accurately and efficiently without the need for remeshing or loss of numerical stability. Therefore, this work adopts the SLE function to represent Fτ and readers are referred to [43] for its detailed formulation. The Green– Lagrange strain components E = [ Exx, Eyy, Ezz, Eyz, Exz, Exy ] and the econd-Piola Kirchoff stress components S = [ Sxx, Syy, Szz, Syz, Sxz, Sxy ] take the form: = ( Blτ + Bnlτ ) uτ (y), (2) = ?̃?E, (3) S.O. Ojo, G. Zucco and P.M. Weaver Thin-Walled Structures 163 (2021) 107672 w v r δ w a f δ A G e Blτ = ⎡


Introduction
Compared to conventional materials, carbon fibre reinforced polymer (CFRP) composite laminates offer high-stiffness, good fatigue resistance and good strength-to-weight properties that make them desirable for advanced structural applications in many industries such as aerospace, sport and medical [1,2]. With excellent design flexibility with respect to fibre orientation, stacking sequence, through-thickness material, and geometric properties, it is feasible to tailor composite materials to meet specific requirements [2,3]. In terms of material design variables, CFRP laminates offer the potential for constant stiffness (CS) or variable stiffness designs (VS) [4,5]. While a single domain structure with uniform, optimised stacking sequence is often the goal of CS designs, a variation in stiffness properties as a function of spatial location is the purpose of VS designs [5,6]. Therefore, VS designs possess greater design space that can potentially achieve better performance compared to CS designs. The current ongoing development of fibre placement technologies such as Automated Fibre Placement (AFP) [7,8], Tailored Fibre Placement (TFP) [9,10] and Continuous Tow Shearing (CTS) [11,12] has bolstered the production of VS laminates for industrial applications. numerical method to evaluate interlaminar stresses in VS composites. Similarly, a p-version finite element (FE) model was proposed in [26,27] to investigate the natural modes of vibration, stresses and nonlinear bending deflection of VS laminates using a Reddy-type thirdorder shear deformation theory, although Groh and Weaver [28] later showed that the Reddy-type model leads to static inconsistencies at clamped boundaries. Furthermore, Soriano and Díaz [29] studied the failure process of VS composites using 3D FE analysis combined with continuum damage mechanics. Moreover, to achieve computationally efficient outcomes in the analysis of VS laminates, Demasi et al. [30] adopted equivalent single layer, layerwise and zigzag models based on FE methods. Despite the improved efficiency over 3D FE models for some simple laminate stacking sequences studied in [30], the potential of these models for analysis of complex and arbitrary lay-ups remains unclear. In addition, FE-based models often require fine meshes which may be computationally expensive to extensively exploit the merits of VS designs. To this end, high-order methods which can guarantee high levels of accuracy combined with faster solution than FE methods become desirable.
The Differential Quadrature Method (DQM) has been widely applied for analysis of CS composite structures, see [31][32][33] for example, due to its promising characteristics of spectral accuracy and fast computation. In the context of VS applications, the Generalised Differential Quadrature (GDQ) has been used to study free vibration and static response of doubly-curved shell structures reinforced by curvilinear fibres in [34,35] while Groh and Weaver applied DQM for accurate stress analysis of VS laminates based on a third-order zig-zag theory within a Hellinger-Reissner mixed variational framework [36]. Furthermore, Luan et al. [37] employed a fifth-order zigzag theory within a Hellinger-Reissner mixed formulation and a mixed inverse DQM approach for accurate stress analysis of VS laminates under different loading conditions. Although DQM implementation of the zigzag formulation for analysis of VS composites leads to significant computational savings over 3D FE solutions, further studies by Patni et al. [38] show that high fidelity models with enriched kinematics are required to resolve the discrepancies in the prediction of transverse stresses which affect delamination. To this end, a high-order serendipity-Lagrange based FE model within Unified Formulation (UF) framework was proposed in [38], which yields accurate predictions for linear static analysis of VS laminates.
To a large extent, the theory of UF has been employed for a variety of mechanics problems ranging from linear to nonlinear analyses (see [39][40][41]), due to the generic qualities of the kinematics adopted that allow for arbitrary expansion of displacement variables, leading to accurate prediction of local and global response of composite structures. With regards to geometrically nonlinear analysis, Pagani and Carrera [40,41] demonstrated that UF implemented within a FE framework is able to capture large displacements/rotations, in-plane deformation, localised buckling, cross-sectional warping, bending-torsion coupling and other high-order phenomena in thin and thick composite structures. Consistent with the benefits of UF for geometrically nonlinear analysis by Pagani and Carrera, Ojo and Weaver [42] proposed a DQMbased geometrically nonlinear strong UF (SUF) model to extend the capabilities of the linear strong UF model proposed in [31] for analysis of CS composites. Outcomes show significant improvement in terms of accuracy and efficiency over the FE-based model developed in [41]. In a bid to explore the potential of the newly developed geometrically nonlinear SUF, this study extends the capability of SUF to investigate the structural response of VS laminates undergoing large deflections. Considering the merits of VS designs for a variety of structural applications, it is important to understand the role of VS properties in achieving improved performance especially in nonlinear, large deflection applications. With this aim in mind, SUF is employed to investigate the effect of tow-steering on the behaviour of composite laminates under large bending and compressive loads. In addition, since efficient prediction of nonlinear 3D stresses is an important requirement to fully characterise large deflections in engineering structures, the potential of SUF for computationally efficient nonlinear analysis of composite structures is assessed comparatively with the models in the literature and ABAQUS model. The rest of the paper is arranged as follows.
The basic formulation for geometrically nonlinear SUF and subsequent discretisation by DQM are presented in Sections 2 and 2.1. Derivations of the constitutive and geometric components of the tangent stiffness matrix are presented in Section 2.2 while Section 3 is dedicated to description of the arc-length routine employed for solution of the nonlinear system. In Section 4, settings for numerical test are explained while mathematical description of the constitutive relations of variable fibre architecture is also outlined. Numerical results and discussions are further presented in Section 4 before the final concluding remarks are detailed in Section 5.

Nonlinear strong unified formulation
According to the theory of Unified Formulation (UF), the displacement components [ , , ] in the three dimensions, , and axes of the beam-type structure of length in the cartesian coordinate system are represented by the relation (see Fig. 1), where is a function that captures the cross-sectional deformation and can be expanded to any order for enrichment of the beam's kinematical description. The Serendipity Lagrange expansion (SLE) introduced in [43] idealises the beam's cross-sectional behaviour accurately and efficiently without the need for remeshing or loss of numerical stability. Therefore, this work adopts the SLE function to represent and readers are referred to [43]  take the form: =̃, , , where Eqs. (4a)-(4c) are given in Box I. Invoking the principle of virtual work, the equilibrium relations for static deformation can be realised as, where and represent, respectively, the virtual strain energy and virtual work of external loadings. The virtual strain energy is further explicitly expressed as, Applying partial integration-by-parts to Eq. (6) after substituting for Green-Lagrange strain components using Eqs. (2)-(4), the virtual strain energy relation becomes, In this work, it assumed that the fibre orientation varies along the axial dimension, , such that the material coefficient tensor,̃( ), is a function of . Therefore, the product rule is applied to derivatives of stresses along the -axis in Eq. (7), since stresses are function of ( ) and̃( ). In compact form, Eq. (7) is rewritten as, where represents the first-order derivative of the material stiffness tensor , and the terms containing account for VS properties which vanish in the case of CS laminates, since is constant in the axial dimension of the structure. The derivation of for SUF is discussed in Section 4. Based on Eq. (10), the secant stiffness operator and its boundary components are respectively expressed as, Explicit expressions for the components of and are given in Appendix.
Note: subscripts , , , represents linear, linear-nonlinear, nonlinear-linear, and nonlinear-nonlinear components of the secant stiffness operator while and represent the order of SLE function.
Based on Eq. (5), Eq. (10) and (11), the governing equation and boundary conditions for geometrically nonlinear static analysis of a general beam-type composite structure are given as, where ( ) consists of external load contributions from volume forces, line forces and concentrated force while accounts for external load contributions from surface forces at the edges of the beam. Equation (12) constitutes a 1D nonlinear strong UF system of equations which can be expanded to any order for solution of geometrically nonlinear problems.

Differential quadrature method ( DQM )
In line with the DQM approach described in [44], the solution of the one-dimensional system of equations, i.e., Eq. (12), is assumed to take a form in which and are approximated as, while higher derivatives of are approximated as, where ( ), (1) and (2) are weighting coefficients of the control variables, with their first and second derivatives respectively, defined by sets of base polynomials (see [44]). Based on Eq. (14), the governing Eq. (12) can be generalised as, Equation (15) constitutes the fundamental nucleus of the nonlinear system governing equations and boundary conditions which can be expanded to any arbitrary order of , = 1, … , and , = 1, … , to form the global system of equations,

Derivation of the tangent stiffness matrix
Equation (16) represents a nonlinear system of equations which require a nonlinear solution algorithm to determine the equilibrium state. Therefore, an incremental linearised scheme described in Section 3 is implemented here. To determine the incremental displacement in the linearised system of equations, the tangent stiffness matrix, , is computed which is typically derived from the linearisation of the variational formulation as follows, Using the procedure outlined in Section 2, the linearisation of constitutive relation is derived as, while the geometric component of the tangent stiffness matrix is obtained through second-order linearisation of the Green-Lagrange strain as, , whereas̃1 and̃2 are given in Box III.
Note: and * are derivative operators for displacement and stress variables, respectively.
In Eq. (19), diag ⟨̃1 ⟩ and diag ⟨̃2 ⟩ are 3 × 3 diagonal matrices, whose diagonal terms are the components of the column vectors̃1 and̃2 , respectively. Therefore, explicitly, the geometric stiffness operator for the interior and boundaries of the 1D structure is expressed as, where is a 3 × 3 identity matrix. Thẽterms in Eq. (21a) represent corresponding stress components with the material stiffness tensor replaced by its first-order derivative .
In accordance with Eqs. (17)- (21), linearisation of the virtual strain energy relation can be rewritten as, and are the fundamental nuclei of the tangent stiffness operator in the interior and boundaries of the beam, respectively. After applying DQM discretisation, the tangent stiffness operators can be generalised as and which can be expanded to any order to form the global tangent stiffness matrix, .
Box III.

Arc-length solution scheme
Solving Eq. (16) requires finding the equilibrium configuration of the structure under the action of applied forces. In this case, Eq. (16) is reconfigured, in line with the arc-length scheme described in [45], as where is the reference load, is the unknown scaling load parameter, and ( , ) is the global residual of nodal forces which is expressed as, Eqs. (23)-(24) constitute the essential elements to implement the arclength method. For the sake of brevity, readers are referred to [45] for a detailed routine of the arc-length procedure.

Numerical examples
To assess the accuracy and efficiency of the proposed SUF model, we consider geometrically nonlinear static analysis of composite laminated beams with geometric and material configurations given in Table 1. In addition, three load configurations are considered as illustrated in Fig. 2, where different CS and VS composite configurations are subject to clamped-free conditions with either uniformly distributed surface load (Q1), or compression (Q2), or shear load (Q3). Numerical validation of our SUF model against models proposed in the literature and ABAQUS 3D FE models is performed for geometrically nonlinear static analysis of some CS laminates to compare respective computational efficiencies. Furthermore, analysis of VS laminates using the SUF model is performed in which the fibre angle, (y), which varies along the axial direction of each ply is defined according to the representation proposed by Gürdal and Olmedo [19] as, where 1 and o represent, respectively, the fibre angle at the edge and midspan of the beam-type structure. Subject to Eq. (25), additional contributions to the structural stiffness matrix of VS laminates result from the first-order derivative of the material stiffness tensor, which is given at a material point, , according to the DQM representation as, It is noted that Gürdal and Olmedo's representation of (y) leads to a discontinuous fibre angle at the beam's midspan since there is a change of fibre orientation in this position. Therefore, accurate representation of the derivative of the material stiffness tensor requires separation of a Material orientation is rotated around the Z-axis by 90 • to match the reference configuration in Fig. 1 represents the 1 laminate orientation angle at the laminate edge and o is the laminate orientation angle at the midspan.

Validation of the SUF model
The test configurations for CS laminates used for numerical validation are given in Table 2, in which geometric, material and laminate nomenclatures are defined according to Table 1. In this context, three tests are considered as detailed in Table 2.

Isotropic cantilevered beam under shear load (Test 1)
Test 1 comprises a clamped-free isotropic beam with geometric, material and load configurations given in Table 1. The SUF model is implemented using a 35-node single domain DQM-based 1D beam element while one fifth-order SLE 2D element is used to describe the cross-sectional kinematics leading to a total of 2415 degrees of freedom (DOF). It is evident from Fig. 3 that SUF estimates of global axial and transverse deflections agree well with the weak Unified Formulation (WUF Patni ) model proposed in [46]. In addition, good agreement is recorded for axial and shear stresses in Fig. 4, highlighting the accuracy of the SUF model. Although, WUF Patni and the present SUF model use the same order of degrees of freedom for this configuration, it is instructive to note that the WUF Patni model requires 4209 DOF to achieve convergence. Therefore, the SUF model is able to save almost 50% of the computational effort compared to the WUF Patni model. More importantly, when compared with the results of the 3D model obtained by ANSYS software [46], which uses 139,623 DOF, the SUF model saves up to two orders of degrees of freedom to achieve convergence, highlighting the efficiency of the proposed model.

Constant stiffness composite under bending load (Test 2)
Test 2 consists of a two-layer cross-ply laminate with structural and load configurations described according to Tables 1 and 2, and Fig. 2b, respectively. The laminate which is subject to a clamped-free boundary condition and a constant transverse load per unit area, p 0 , undergoes large bending deflections. Results obtained from SUF are benchmarked  against Pagani and Carrera's layerwise model (herein referred to as Lmodels) proposed in [41] and ABAQUS 3D FE models. According to Fig. 5, SUF estimates of transverse deflection throughout the loading regimes compare well with results obtained from L-and ABAQUS 3D FE models. Comparatively, SUF achieved convergence of the solution for the nonlinear problem with a third-order 2D SLE element for each layer coupled with a 17-node DQM beam element leading to 1020 DOF, while it was reported in [41] that the L-model requires a total of 5124 DOF to accurately capture the global and local response of this test configuration. Furthermore, the ABAQUS 3D model is completed with 43,200 C3D20R reduced integration 3D solid elements which gives 578,397 DOF in total. Therefore, SUF demonstrates improved efficiency over L-and ABAQUS models for this study. It should be noted that ABAQUS results were obtained using a somewhat coarse mesh due to limitation of computational capacity, hence the minor discrepancy between results obtained from SUF and L-model and estimates from the ABAQUS model.

Postbuckling response of CS composite beam under compression (Test 3)
In this setting (see Table 2), the cross-ply laminate is subject to compression load P under a clamped-free condition (see Fig. 2c) in which the postbuckling behaviour of the composite structure is investigated. Global deflection of the laminate is then validated against both an Lmodel and an ABAQUS 3D FE model. Instructively, to the best of the authors' knowledge, this is the first study to present ABAQUS validation of this laminate configuration. According to Fig. 6, the potential of SUF to accurately predict the global response of the structure from the linear regime through the deep postbuckling regime is evident, as the result satisfactorily agrees with estimates from both L-and ABAQUS 3D FE models. The ABAQUS model was undertaken using 43,200 C3D20R reduced integration solid elements constituting 578,397 DOF while SUF achieved convergence with 2520 DOF (consisting of a fifth-order 2D element per layer and a 21-node DQM-based 1D beam element). In terms of accuracy, SUF and the L-model are better than the ABAQUS model which is attributed to the limitation of computational capacity under current hardware limitations to sufficiently refine the ABAQUS model. Nonetheless, all three models show good general agreement throughout the loading regimes. However, SUF shows superior efficiency to ABAQUS models since far fewer degrees of freedom are required to attain convergence.

Analysis of variable stiffness ( vs) composites
VS composites can be considered by tow steering of individual plies in a two-layer cross-ply laminate. Two laminate configurations were selected comprising: (a) 0 • bottom layer steered at 10 • at the midspan and 90 • topmost layer steered at 80 • at the midspan and (b) 0 • bottom layer steered at 10 • at the midspan and 90 • topmost layer steered at 70 • at the midspan. Essentially, the properties of these VS laminates differ by steering the 90 • layer at different angles at the midspan while the steering of the 0 • layer is fixed at an orientation angle of 10 • at the midspan. These configurations allow the effect of steering of  individual layers to be investigated for enhanced performance. The inplane variation of material stiffness, modelled as a continuous function within a ply along the axial dimension leads to extra terms to those from CS considerations in the SUF's structural stiffness matrix, which depend on the first-order derivatives of material stiffnesses denoted as̃. Figs. 7 and 8 show variations of the first-order derivatives of axial and in-plane shear stiffnesses of individual layers along the fibre direction, . With these variations, it is possible to explore the potential of SUF for capturing different coupling effects as well as other nonclassical effects in the VS laminates. In general, the converged results of all VS laminates reported in this section are realised with a fifthorder SLE element for each layer and single domain DQM-based beam element with 27-35 nodes, leading to a total number of 3240-4200 DOF. Test configurations for the analysis of VS laminates are given in Table 3.

ABAQUS modelling of VS composites
Variable stiffness (VS) composites cannot be modelled directly in ABAQUS. Therefore, a Python script [47] was used to define the equivalent problem in ABAQUS using the Riks algorithm with threedimensional finite elements. In particular, the beam is discretised by using a C3D20R element, a 20-node quadratic brick with reduced integration. For each element, three integration points through the thickness is specified. For all ABAQUS results presented in this section, each beam is discretised with 43,200 3D solid elements, which corresponds to a total of 578,397 DOF, where each cubic element has side length of 50 mm. In this manner, volumetric locking effects typical of three-dimensional FE are avoided [48].
For the beam with load Q1 configuration, a surface load is applied without a following rotation. Whereas, for both load Q2 and Q3 configurations, the load is introduced using a multi-point constraint (MPC) approach as shown in Fig. 9. More precisely, all nodes lying on the vertical surface of the beam are connected to a control node where the concentrated force (in the Y direction for the load Q2 configuration or in the Z direction for the load Q3 configuration, respectively) is applied. For all geometries under consideration, boundary conditions = = = 0 are applied to the entire surface of the beam lying on the XZ plane at Y = L (see Fig. 9). To the best of our knowledge, for both load Q2 and Q3 configurations, no work has been presented in the literature showing comparisons between equilibrium curves obtained with UF and those obtained with ABAQUS using the Riks algorithm with three-dimensional FE analysis.

VS Composite beam under bending load (Test 4)
Fig. 10 reveals good agreement of global deflections for different VS laminates between ABAQUS 3D FE and SUF models throughout the loading regimes. These results demonstrate the potential of SUF for accurate prediction of the behaviour of VS structures under large displacements and rotations. In addition, SUF demonstrates superior computational efficiency compared to ABAQUS since SUF uses 3240 DOF (consisting of one fifth-order SLE 2D element per layer and a 27node DQM 1D beam element) to converge against 578,397 DOF for ABAQUS model.
Under a transversely distributed load, the VS laminates L2 and L3 show negligible deviation in their global and local behaviours, as illustrated in Figs. 10 and 11. Considering that the variable stiffness property in laminates L2 and L3 affects the axial and shear stiffnesses of the beam according to Figs. 7 and 8, it can be deduced that there is a minimal contribution of variable fibre orientation to the deformation of the beam under bending load. Moreover, compared with the crossply CS laminate, there is no substantial change in the global bending behaviour. However, Fig. 11 shows that there is a glaring deviation between the local behaviour of VS and CS structures in this context, which is evident in the transverse shear stress distribution of the 0 • towsteered layer. This deviation is primarily attributed to discrepancies in the local bending behaviour of VS and CS structures. Since a 0 • towsteered layer is much stiffer than a 90 • tow-steered layer, the local  bending behaviour of the 3D structure is dominated by the response of the 0 • layer. Fig. 12 shows satisfactory agreement between ABAQUS 3D FE and SUF models for this test configuration. The global and local responses of the cantilevered VS laminates, i.e, L2 and L3, under a shear load are similar to those under bending in terms of the role of the 0 • tow-steered layer in the overall structural response. Notably, global and local deformations are controlled by bending. The global response of the VS laminates is dominated by the contributions of the 0 • tow-steered layer, whereas the responses of the CS and VS structures differ only by local bending. The effects of in-plane or out-of-plane shear couplings are negligible in the structural response of the two VS configurations due to steering of the 90 • layers (see Fig. 13).

Postbuckling response of VS composite beam (Test 6)
Due to the discrete representation of the continuous tow-steered fibre, ABAQUS requires a fine mesh to accurately capture the response of VS laminates considered herein. With 43,200 C3D20R elements, the ABAQUS model requires at least 578,397 DOF to converge to estimates from SUF which uses up to 4200 DOF to achieve convergence (see Fig. 14).
Comparing the structural responses of VS and CS laminates in Figs. 15 and 16, it can be observed that tow-steering shifts the structural response of the CS laminates from a stiff to a less stiff configuration at different loading regimes. In line with Figs. 7 and 8, tow-steering affects the in-plane and out-of-plane shear properties of the composite beam. Due to the variation of these properties along the axial dimension of the composite beam, different sections of the beam may experience different modes of deformation under large compression loads. Consequently, the VS laminates may exhibit different variable stiffness structural configurations under different loading phases, a phenomenon that is attributed to fluctuations of the contributions of the bendingtorsion and compression-bending couplings to the behaviour of VS laminates. In the less stiff configuration after the critical buckling load, the VS structures exhibit greater compliance to the moderate compression load than the CS laminate due to higher contribution of compression-bending coupling to the structural response. This phenomenon leads to higher stresses in the VS laminates compared to CS laminates in this loading regime, as shown in Fig. 15. As the compression of the beam becomes large, the VS laminates experience large rotations and displacements, and bending-torsion effects become the dominant mechanism of deformation. In this regard, according to Fig. 14, the CS laminate experiences higher bending deformation compared to VS laminates leading to a stiffer response of the VS structures. This observation is corroborated by Fig. 16, where VS laminates exhibit lower levels of transverse shear stress than CS laminates. Therefore, tow-steering can be used to achieve variable structural response under different loading regimes.

Conclusion
A geometrically nonlinear strong Unified Formulation (SUF) has been developed for the three-dimensional large deflection analysis of constant stiffness (CS) and variable stiffness (VS) composite structures. For accurate representation of the continuous tow-steered fibre trajectory by the SUF model, the material domain is separated into two parts to enable the fibre angle discontinuity at the beam's midspan to be captured. On the other hand, ABAQUS 3D FE modelling of VS structures is accomplished by a discrete representation of the fibre orientation, thus necessitating a fine mesh for convergence. The proposed geometrically nonlinear SUF is shown to accurately capture the behaviour of VS and CS laminates undergoing nonlinear large deflections due to different loading conditions. In addition, SUF exhibits superior efficiency over weak-formulation-based models in the literature as well as from our own ABAQUS 3D FE solutions. The SUF model requires a maximum of 2520 DOF and 4200 DOF, respectively, to capture the 3D behaviour of the CS and VS laminates, thereby outperforming the computational efficiencies, for similar levels of accuracy, of the weak-formulation SLE-UF model [46], weak-formulation layerwise-UF model [41] and ABAQUS 3D FE model by 44%-99%. Based on the outcomes of the large deflection analysis, it is shown that, under different loading conditions, VS laminates of similar structural configurations exhibit different behaviours under different loading conditions in terms of the roles of individual tow-steered plies. Under bending and shear loading conditions, the structural responses of VS and CS structures differ by local stress redistribution. However, under compression, VS properties may facilitate variable structural configurations under different loading regimes, highlighting the potential of tow-steering to tailor and enhance properties of composite structures.
Although not affecting our conclusions, it is noted that due to the limitation of computational capacity, optimum refinement of the ABAQUS 3D FE models could not be accomplished leading to small discrepancies between SUF and ABAQUS estimates, with the former providing slightly greater accuracy.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.