Reﬁned Beam Theory for Geometrically Nonlinear Pre-Twisted Structures

: This paper proposes a novel fully nonlinear reﬁned beam element for pre-twisted structures undergoing large deformation and ﬁnite untwisting. The present model is constructed in the twisted basis to account for the effects of geometrical nonlinearity and initial twist. Cross-sectional deformation is allowed by introducing Lagrange polynomials in the framework of a Carrera uniﬁed formulation. The principle of virtual work is applied to obtain the Green–Lagrange strain tensor and second Piola–Kirchhoff stress tensor. In the nonlinear governing formulation, expressions are given for secant and tangent matrices with linear, nonlinear, and geometrically stiffening contributions. The developed beam model could detect the coupled axial, torsional, and ﬂexure deformations, as well as the local deformations around the point of application of the force. The maximum difference between the present deformation results and those of shell/solid ﬁnite element simulations is 6%. Compared to traditional beam theories and ﬁnite element models, the proposed method signiﬁcantly reduces the computational complexity and cost by implementing constant beam elements in the twisted basis. theory [32], Saint-Venant beam theory [33], and shell ﬁnite element results. For the present method, equal or more than 10 B2 elements are required to obtain results consistent with those obtained by the shell model. The maximum difference of the torsional rotational angle between the 1 L9 + 10 B2 model and ABAQUS model is 1%, indicating that the present method is more accurate than traditional beam theory. With a total end load of 338 kN, the untwisted angle is approximately 0.72rad (41.25 ◦ ). Therefore, the deformed structure is nearly straight after untwisting.


Introduction
The increasing application of initially twisted structures has attracted numerous researchers due to their structural strengthening capabilities [1,2] and highly coupled mechanical properties [3]. However, these twisted structures may be subjected to large loads, such as centrifugal forces, actuator forces to achieve shape morphing, and deploying forces to achieve foldable designs. In-depth investigations usually lack computationally friendly tools to capture their complex nonlinear behaviors-i.e., large bending-bending-torsional deformation, cross-sectional in-plane and out-plane deformation. The development of refined beam theories may provide a proper choice to consider nonlinear and coupled effects with reasonable computational cost. The refinement of the cross-sectional displacement field could advance the beam theories in two respects: (i) the distribution of generalized displacements over the cross section is approximated by two-dimensional polynomials in the cross-sectional domain, and (ii) the choice of expansion functions allows the precise mapping of the cross section.
In the pioneering work of Giavotto et al. [4], Saint-Venant warping was used to refine the cross-sectional field description in the three-dimensional beam formulation. This technique was later extended to the case of large displacement formulation [5] and curved and twisted beam formulation [6]. Exact geometrical theories further separated the beam model to a displacement field of the reference axis and a rotational field of the cross-section to incorporate torsional warping [7]. By characterizing the reference axis with arbitrary cross section differently, it was possible to treat large deformation analysis of space-curved beams with eccentric cross section [8], finite stretches and rotations [9], co-rotational beam elements with finite 3D rotations [10], continuous translation and rotation [11], and thinwalled beams formulated using the special Euclidean group [12]. To obtain a generalized Figure 1 shows the geometrical description of a pre-twisted structure defined in four different coordinate systems. They are the Cartesian orthonormal basis (e x , e y , e z ) located at one end of the pre-twisted structure, the twisted orthonormal basis (e ξ , e y , e η ) translating and rotating along the longitudinal direction to trace the structural twist, and the covariant basis (g 1 , g 2 , g 3 ) and contravariant basis (g 1 , g 2 , g 3 ) describing the deformation gradient in the twisted basis. Figure 1. Geometrical description of a pre-twisted structure in the reference configuration with the Cartesian basis (ex, ey, ez), the twisted basis (eξ, ey, eη), the covariant basis (g1, g2, g3), and the contravariant basis (g 1 , g 2 , g 3 ).
The twisted base vectors are related to the Cartesian coordinate system by the transformation matrix T as follows: T e e e (1) where the transformation matrix is expressed as To derive the differential operators in the twisted basis, the local covariant basis is introduced as follows: , , / , / , / y g g g X X X (4) From Equations (5)- (7), the relationship between the covariant basis and the twisted basis can be written as  (5) and in matrix form as ξ η φη φξ where the transformation matrix is expressed as T =   cos(φy + θ) 0 sin(φy + θ) 0 1 0 − sin(φy + θ) 0 cos(φy + θ)   (2) with φ being the twist angle per unit length and θ as the phase angle in the y = 0 plane. For the sake of simplification, the notations c = cos(φy + θ) and s = sin(φy + θ) are introduced in the following analysis. The position vector of any point of the structure in the twisted basis is defined as X(x, y, z) = ξe ξ (y) + ye y + ηe η (y) To derive the differential operators in the twisted basis, the local covariant basis is introduced as follows: (g 1 , g 2 , g 3 ) = (∂X/∂ξ, ∂X/∂y, ∂X/∂η) From Equations (5)- (7), the relationship between the covariant basis and the twisted basis can be written as g 1 = e ξ (y), g 2 = −φηe ξ (y) + e y + φξe η (y), g 3 =e η (y) (5) and in matrix form as The underlined covariant basis (g 1 , g 2 , g 3 ) is a local basis and varies with the material point considered. This is an obvious consequence of the curvedness of the coordinate system. This is the result of Equation (4), where the vectors of the covariant basis at a point are tangent to the three coordinate curves passing through it.
It is worth noting that in the covariant basis, the components of the vectors and the tensors are defined locally. For instance, the covariant metric tensor, defined by g mn = g m · g n , is given by By introducing g m · g n = δ n m , ∀m, n = 1, 2, 3., the contravariant basis (g 1 , g 2 , g 3 ) is defined. The contravariant metric tensor components relate the covariant counterparts by the following expression: Using Equations (7) and (8), the contravariant components of the metric tensor is derived as follows: Then, we can define the local contravariant basis by the following: The components of the base vectors are specifically written as g 1 = e ξ (y) + φηe y , g 2 = e y , g 3 = −φξe y + e η (y) (11) and in the matrix form as The Christoffel symbols of the second kind Γ k ij , i, j, k = 1, 2, 3, associated with the covariant and contravariant bases by Γ k ij = g i,j ·g k , are given by

Nonlinear Green Lagrange Strain
The deformation components of an arbitrary point in the twisted basis can be described as For the displacement vector u ξ , u y , u η of a specific point, the covariant and contravariant components are specified by the following expressions: The 3D Green-Lagrange strain tensor is given by half the increment in the 3D metric tensor. When written in the covariant basis of the curvilinear coordinate system, the strain-displacement relationship is given by where Substituting Equations (13) and (15) into Equation (17), the first-order covariant derivatives of the covariant and contravariant displacement components are obtained as follows: From Equations (16) and (18), the Green-Lagrange strain vector can be written in the contravariant basis as follows: where the expressions of b cov,l and b cov,nl are given in Equations (A1) and (A2), respectively.
To derive the components of the strain vectors in the twisted basis, the relationship between the deformation gradients in the twisted basis and the covariant basis is established as follows: where Greek subscripts α and β denote components ξ, y, and η in the twisted basis. Accordingly, Therefore, the Green-Lagrange strain tensor can be expressed in the twisted basis as where the expressions of b l and b nl are given in Equations (A3) and (A4), respectively.
Here, the coefficients of the differential operators used in the twisted basis (i.e., ∂ξ, ∂η, Λ) are independent of the axial coordinate. For initially twisted helical structures, the Aerospace 2022, 9, 360 6 of 21 translational invariance of the differential operator coefficients is essential to assume that the beam mechanical properties do not vary along the axial direction in the twisted basis [28].

Constitutive Law
A constitutive law of linearly elastic material establishes a linear relationship between the second Piola-Kirchhoff stress tensor σ and the Green-Lagrange strain tensor ε by Since these tensors are second order, the elasticity tensor is given in the covariant basis by the following fourth-order tensor: In the twisted basis, the elasticity tensor can be defined as where the Greek subscripts α, β, δ, and γ denote the components ξ, y, and η in the twisted basis. By introducing the Kronecker δ, the tensor components in the twisted basis are obtained as Finally, the elasticity coefficients that relate the stress vector to the strain vector in the twisted basis can be given by The above expression coincides with the one obtained in the Cartesian basis, as the twisted basis is orthonormal.

Carrera Unified Formulation
To include geometrical nonlinearity of thin structures, the displacement field is enriched with the higher order terms in the framework of the Carrera unified formulation (CUF). Herein, the displacement vector u at any point is approximated by the following expression: where q sj is the generalized displacement vector, N j represents the shape functions along the axial direction, and F s denotes the expansion functions of the beam cross section, with M u as the order of F s . Here, we use Einstein's summation convention-i.e., summation signs are not written for all indices that repeat once as a subscript and once as a superscript in an expression. Additionally, subscripts and superscripts s and j are called dummy indices, since they can be replaced by any other letters-i.e., τ and i. Several different types of polynomials could be used to construct the expansion functions describing the cross-sectional geometry, including Taylor expansions (TE), Lagrange expansions (LE), and Hierarchical Legendre expansions (HLE). The TE expands the polynomial globally over the cross section from the beam axis, while the LE beam theories discretize the cross section physical surface into a number of local expansion sub-domains, whose polynomial degree depends on the type of Lagrange expansion employed. Threenode linear L3, four-node bilinear L4, nine-node quadratic L9, and sixteen-node cubic L16 polynomials have been developed in the framework of CUF. HLE models combine the main features of the previous beam theories, including the hierarchy of the high-order terms of TE and the geometric discretization of the beam section surface of LE. For the sake of brevity, their expressions are not included here, but they can be found in the literature [29].
In the framework of CUF, the strain vector in the twisted basis is calculated as follows: where the expressions of B

Equations of Motion
Substituting the stress and strain vectors, the variation of the elastic potential energy δU is given by According to the definition, the secant stiffness matrix K ijτs s relates to the strain energy by where K ijτs l represents the linear component, K ijτs lnl and K ijτs nll are the first-order nonlinear contributions, and K ijτs n ln l contains the second-order nonlinearities. They are expressed as Since the secant stiffness matrix is generally not symmetrical, the numerical results are merely with low orders of convergence. A symmetrical tangent stiffness matrix can be derived by linearizing the virtual variation of the strain energy as follows: From Equation (29), the linearization of the virtual variation of the strain is calculated as where B * nl_1 and B * nl_2 are given in Equations (A9) and (A10), respectively. Then, the first contribution of the virtual variation of the strain energy δ(δU) is derived as where the expression of K ijτs σ is given in Equation (A11). The second term of the virtual variation of the strain energy δ(δU) is defined as Combining Equations (34), (36) and (37), the tangent stiffness matrix is given by where K ijτs l , K ijτs T1 , and K ijτs σ are the linear stiffness matrix, the nonlinear stiffness matrix, and the geometrically stiffening stiffness matrix, respectively. Unlike Ks, all the contributions of the tangent stiffness matrix are symmetrical.
When the pre-twist contributions are not included, the secant and tangent stiffness matrices have the same expression as those given by [30]. If the pre-twist terms are included, the strong coupling of nonlinear axial and cross-sectional deformations is reflected.
As for the external work, the contributions from volume forces g, surface forces p, line forces q, and concentrated force P at point Q are considered. The variation of the external work is written as From Equation (32), the nonlinear governing equation is formulated as follows: where the generalized force vector F τi is expressed as

Arc-Length Method
Using Crisfield's arc-length method, Equation (41) is rewritten as where R τi represents the generalized out-of-balance force vector, and λ is the scalar load parameter. Along the response curve of the nonlinear deformation, the incremental expression of the governing equation becomes At any incremental step of the converged solution, the iteration direction is orthogonal to the tangent of the response curve. By fixing the length of the incremental step, the constraint equation with a constant arc length is formulated as where ∆S is the finite increment of the arc length. The generalized displacement increment at iteration m from step p to p + 1 can be expressed as q Then, the quadratic constraint equation at iteration m from step p to p + 1 is given by The appropriate generalized displacement increment solution can be evaluated by imposing the least angle between q p,m sj and q p,m+1 sj .

Results
The nonlinear deformation simulation of three kinds of pre-twisted structures are presented in this section. Figures 2-4 show the simplified geometries of pre0twisted structures with rectangular, arc profile, and airfoil profile cross sections, respectively. All those structures are assumed to be uniformly pre-twisted and the dimensions are listed in Table 1.

Results
The nonlinear deformation simulation of three kinds of pre-twisted structures are presented in this section. Figures 2-4 show the simplified geometries of pre0twisted structures with rectangular, arc profile, and airfoil profile cross sections, respectively. All those structures are assumed to be uniformly pre-twisted and the dimensions are listed in Table  1.
The rectangular cross section is characterized by the width b1 and the height h1, and the pre-twist center P1 coincides with the geometrical center C1. The physical properties of the arc profile cross sections are taken from widely investigated twisted cylindrical panels [31]. These dimensions give shell configurations with the same width b2 and thickness h2, and the pre-twist center P2 is located at the center of the middle arc chord. The arc profile cross-sectional structures are distinguished from each other by the arc angle θ2 and the pre-twist angle k2. As for the airfoil profile section, the same NACA-8405 geometry is selected in accordance with the authors' previous paper [2]. The 4 digits define the maximum camber, its position, and the maximum thickness, and the pre-twist center P3 is located at the geometrical center C3 of the airfoil section. The material of these structures is assumed to be aluminum alloy, with Young's modulus of 70GPa and Poisson's ratio of 0.3.

Pretwisted Cantilever with Rectangular Cross Section
A thin-walled, rectangular, cross-sectional cantilever with 45° pre-twist is first considered. A uniformly distributed longitudinal load (338 kN in total) is applied to the freeend cross section, resulting in the untwisting of the structure.
In our model, depicted in Figure 5b, one L9 Lagrange element is implemented to describe the cross-sectional displacement field. In the longitudinal direction, 5, 10, and 20 two-node (B2) beam elements are applied for discretization to study the convergence rate. The corresponding degree of freedoms (DOFs) are 162, 297, and 567. The obtained nonlinear deformation results are compared with those obtained by shell models from

Pretwisted Cantilever with Rectangular Cross Section
A thin-walled, rectangular, cross-sectional cantilever with 45° pre-twist is first considered. A uniformly distributed longitudinal load (338 kN in total) is applied to the freeend cross section, resulting in the untwisting of the structure.
In our model, depicted in Figure 5b, one L9 Lagrange element is implemented to describe the cross-sectional displacement field. In the longitudinal direction, 5, 10, and 20 two-node (B2) beam elements are applied for discretization to study the convergence rate. The corresponding degree of freedoms (DOFs) are 162, 297, and 567. The obtained nonlinear deformation results are compared with those obtained by shell models from  The rectangular cross section is characterized by the width b 1 and the height h 1 , and the pre-twist center P 1 coincides with the geometrical center C 1 . The physical properties of the arc profile cross sections are taken from widely investigated twisted cylindrical panels [31]. These dimensions give shell configurations with the same width b 2 and thickness h 2 , and the pre-twist center P 2 is located at the center of the middle arc chord. The arc profile cross-sectional structures are distinguished from each other by the arc angle θ 2 and the pretwist angle k 2 . As for the airfoil profile section, the same NACA-8405 geometry is selected in accordance with the authors' previous paper [2]. The 4 digits define the maximum camber, its position, and the maximum thickness, and the pre-twist center P 3 is located at the geometrical center C 3 of the airfoil section. The material of these structures is assumed to be aluminum alloy, with Young's modulus of 70GPa and Poisson's ratio of 0.3.

Pretwisted Cantilever with Rectangular Cross Section
A thin-walled, rectangular, cross-sectional cantilever with 45 • pre-twist is first considered. A uniformly distributed longitudinal load (338 kN in total) is applied to the free-end cross section, resulting in the untwisting of the structure.
In our model, depicted in Figure 5b, one L9 Lagrange element is implemented to describe the cross-sectional displacement field. In the longitudinal direction, 5, 10, and 20 two-node (B2) beam elements are applied for discretization to study the convergence rate. The corresponding degree of freedoms (DOFs) are 162, 297, and 567. The obtained nonlinear deformation results are compared with those obtained by shell models from ABAQUS, a popular commercial software. To achieve stable convergence, the shell model consists of 988 quadrilateral elements shown in Figure 5c, with a total DOF of 6468.
x FOR PEER REVIEW 12 of 22 ABAQUS, a popular commercial software. To achieve stable convergence, the shell model consists of 988 quadrilateral elements shown in Figure 5c, with a total DOF of 6468.  Figure 6 shows the 3D structural configurations obtained from 10 B2 + 1 L9 CUF model and ABAQUS shell model, respectively. Results show that the nonlinear axial-torsional coupling of the cantilever is accurately captured. The whole structure undergoes a large untwist deformation under the applied longitudinal load. This interesting phenomenon is attributed to the coupling of torsion and extension through pre-twist. Actually, the axial stress component becomes inclined in the twisted basis and generates a net torque that twists the structure in the direction opposite to its initial torsion. The rotational center of the cross section coincides with the pre-twist center of the structure. Moreover, the free-end section keeps the rectangular profile, which reflects the mechanical property of the thin-walled structure.  Figure 6 shows the 3D structural configurations obtained from 10 B2 + 1 L9 CUF model and ABAQUS shell model, respectively. Results show that the nonlinear axial-torsional coupling of the cantilever is accurately captured. The whole structure undergoes a large untwist deformation under the applied longitudinal load. This interesting phenomenon is attributed to the coupling of torsion and extension through pre-twist. Actually, the axial stress component becomes inclined in the twisted basis and generates a net torque that twists the structure in the direction opposite to its initial torsion. The rotational center of the cross section coincides with the pre-twist center of the structure. Moreover, the free-end section keeps the rectangular profile, which reflects the mechanical property of the thin-walled structure.
sional coupling of the cantilever is accurately captured. The whole structure undergoes a large untwist deformation under the applied longitudinal load. This interesting phenomenon is attributed to the coupling of torsion and extension through pre-twist. Actually, the axial stress component becomes inclined in the twisted basis and generates a net torque that twists the structure in the direction opposite to its initial torsion. The rotational center of the cross section coincides with the pre-twist center of the structure. Moreover, the free-end section keeps the rectangular profile, which reflects the mechanical property of the thin-walled structure.   Figure 7 plots the rotation of the beam free-end section obtained by the present theory, Euler-Bernoulli beam theory [32], Saint-Venant beam theory [33], and shell finite element results. For the present method, equal or more than 10 B2 elements are required to obtain results consistent with those obtained by the shell model. The maximum difference of the torsional rotational angle between the 1 L9 + 10 B2 model and ABAQUS model is 1%, indicating that the present method is more accurate than traditional beam theory. With a total end load of 338 kN, the untwisted angle is approximately 0.72rad (41.25 • ). Therefore, the deformed structure is nearly straight after untwisting.
Aerospace 2022, 9, x FOR PEER REVIEW 13 of 22 Figure 7 plots the rotation of the beam free-end section obtained by the present theory, Euler-Bernoulli beam theory [32], Saint-Venant beam theory [33], and shell finite element results. For the present method, equal or more than 10 B2 elements are required to obtain results consistent with those obtained by the shell model. The maximum difference of the torsional rotational angle between the 1 L9 + 10 B2 model and ABAQUS model is 1%, indicating that the present method is more accurate than traditional beam theory. With a total end load of 338 kN, the untwisted angle is approximately 0.72rad (41.25°). Therefore, the deformed structure is nearly straight after untwisting.

Pretwisted Cantilever with arc Profile Cross Section
Pre-twisted structures with different cross-sectional arc angles and pre-twist angles are studied. First, a uniformly distributed axial load (347 kN in total) is applied to the freeend cross section.
As shown in Figure 8b, one L9 element is implemented to construct the physical boundary and the second-order displacement field within the cross section. Since L9 set polynomials can be seen as parabolic expansions plus cubic and quartic terms, the use of these polynomials allows for an accurate description of the second-order cross-sectional field. In addition, 10 B2 elements are exploited along the longitudinal direction.
The present method formulates the governing equation using axially invariant ele-

Pretwisted Cantilever with arc Profile Cross Section
Pre-twisted structures with different cross-sectional arc angles and pre-twist angles are studied. First, a uniformly distributed axial load (347 kN in total) is applied to the free-end cross section.
As shown in Figure 8b, one L9 element is implemented to construct the physical boundary and the second-order displacement field within the cross section. Since L9 set polynomials can be seen as parabolic expansions plus cubic and quartic terms, the use of these polynomials allows for an accurate description of the second-order cross-sectional field. In addition, 10 B2 elements are exploited along the longitudinal direction.  Figure 9 provides the 3D structural configurations of a 30° pre-twisted cantilever with a 90° arc-shaped cross section. Comparison results are given for 10 B2 + 1 L9 CUF model and ABAQUS shell model with uniform end loads. It is shown that the present model can accurately capture the nonlinear axial-torsional coupling of the cantilever, with the torsional center consistent with the pre-twist center. The deformed structure is almost straight after untwisting and the free-end section maintains the arc profile.  The present method formulates the governing equation using axially invariant elements with a separated structural pre-twist. Thus, we can use the same elements to discretize pre-twisted structures with any pre-twist angle, as long as their cross-sectional geometry remains the same. In other words, for pre-twisted structures with different pre-twist angles, we only need to mesh the structures once. Therefore, the complexity and cost of constructing the governing equation are greatly reduced.
The present nonlinear simulation results are compared with those obtained from commercial shell models. As shown in Figure 8c, the shell model in ABAQUS uses the same discretization of 1008 quadrilateral elements. This discretization with 6468 DOFs is found to be the proper choice for achieving stable convergence. Figure 9 provides the 3D structural configurations of a 30 • pre-twisted cantilever with a 90 • arc-shaped cross section. Comparison results are given for 10 B2 + 1 L9 CUF model and ABAQUS shell model with uniform end loads. It is shown that the present model can accurately capture the nonlinear axial-torsional coupling of the cantilever, with the torsional center consistent with the pre-twist center. The deformed structure is almost straight after untwisting and the free-end section maintains the arc profile.
Then, the effect of geometric differences on the nonlinear deformation results is investigated. Figures 10 and 11 show the end section torsional rotational angles of 30 • and 60 • pre-twisted shells with different arc angles, respectively. The rotation angles calculated by the present theory are consistent with those obtained by ABAQUS. For cantilevers with 30 • , 60 • , and 90 • arc profile sections, the maximum differences are 3.2%, 6%, and 2.5%, respectively. Figure 9 provides the 3D structural configurations of a 30° pre-twisted cantilever with a 90° arc-shaped cross section. Comparison results are given for 10 B2 + 1 L9 CUF model and ABAQUS shell model with uniform end loads. It is shown that the present model can accurately capture the nonlinear axial-torsional coupling of the cantilever, with the torsional center consistent with the pre-twist center. The deformed structure is almost straight after untwisting and the free-end section maintains the arc profile. Then, the effect of geometric differences on the nonlinear deformation results is investigated. Figures 10 and 11 show the end section torsional rotational angles of 30° and 60° pre-twisted shells with different arc angles, respectively. The rotation angles calculated by the present theory are consistent with those obtained by ABAQUS. For cantilevers with 30°, 60°, and 90° arc profile sections, the maximum differences are 3.2%, 6%, and 2.5%, respectively. As shown in Figure 10, the torsional angle shows a maximum reduction of 23% as the arc angle increases from 30° to 90°. Correspondingly, this reflects a 30% increase in the torsional stiffness, while in Figure 11, the maximum reduction in the torsional angle is 15%, and the torsional stiffness increases by 17%. In all, a pre-twisted structure with larger cross-sectional curvature has higher stiffness and undergoes less untwisting. Therefore, structural stiffening can be achieved by increasing the structural curvature. However, this stiffening decreases when the initial twist increases.   60° pre-twisted shells with different arc angles, respectively. The rotation angles calculated by the present theory are consistent with those obtained by ABAQUS. For cantilevers with 30°, 60°, and 90° arc profile sections, the maximum differences are 3.2%, 6%, and 2.5%, respectively. As shown in Figure 10, the torsional angle shows a maximum reduction of 23% as the arc angle increases from 30° to 90°. Correspondingly, this reflects a 30% increase in the torsional stiffness, while in Figure 11, the maximum reduction in the torsional angle is 15%, and the torsional stiffness increases by 17%. In all, a pre-twisted structure with larger cross-sectional curvature has higher stiffness and undergoes less untwisting. Therefore, structural stiffening can be achieved by increasing the structural curvature. However, this stiffening decreases when the initial twist increases.   As shown in Figure 10, the torsional angle shows a maximum reduction of 23% as the arc angle increases from 30 • to 90 • . Correspondingly, this reflects a 30% increase in the torsional stiffness, while in Figure 11, the maximum reduction in the torsional angle is 15%, and the torsional stiffness increases by 17%. In all, a pre-twisted structure with larger cross-sectional curvature has higher stiffness and undergoes less untwisting. Therefore, structural stiffening can be achieved by increasing the structural curvature. However, this stiffening decreases when the initial twist increases.
Finally, the large deformation analysis of a pre-twisted cantilever with central and edge point loads is performed. Here, 4 L9 cross-sectional elements and 10 B2 axial elements are used to solve the problem. As shown in Figure 12, the present model can accurately capture the local deformations around the point of application of the force. In Figure 13, the present model can accurately capture coupled axial, torsional, and flexure deformations, as well as the local deformations around the point of application of the force. Finally, the large deformation analysis of a pre-twisted cantilever with central and edge point loads is performed. Here, 4 L9 cross-sectional elements and 10 B2 axial elements are used to solve the problem. As shown in Figure 12, the present model can accurately capture the local deformations around the point of application of the force. In Figure  13, the present model can accurately capture coupled axial, torsional, and flexure deformations, as well as the local deformations around the point of application of the force.

Pretwisted Cantilever with Airfoil Profile Cross Section
A pre-twisted airfoil profile cantilever with NACA 8405 dimensions is investigated. A uniformly distributed longitudinal load (636 kN in total) is applied to the free-end cross section.
Ten B2 and 8 L9 elements are found to be proper discretization choices to simulate the nonlinear deformation. Figure 14b plots the L9 element distribution in the cross section. In the 0-5% section of the camber line, the curvature of the thickness distribution changes relatively rapidly. Therefore, 3 L9 elements are implemented to construct the cross-sectional displacement field. As for the remaining 95% section, it is appropriate to select 5 L9 elements to describe the smooth camber line and thickness distribution. In the ABAQUS solid model, 5900 hexahedral elements are needed to achieve convergence. Each element contains 8 3-DOF nodes and the total number of DOFs is 26,469, as shown in Figure 14c. Finally, the large deformation analysis of a pre-twisted cantilever with central and edge point loads is performed. Here, 4 L9 cross-sectional elements and 10 B2 axial elements are used to solve the problem. As shown in Figure 12, the present model can accurately capture the local deformations around the point of application of the force. In Figure  13, the present model can accurately capture coupled axial, torsional, and flexure deformations, as well as the local deformations around the point of application of the force.

Pretwisted Cantilever with Airfoil Profile Cross Section
A pre-twisted airfoil profile cantilever with NACA 8405 dimensions is investigated. A uniformly distributed longitudinal load (636 kN in total) is applied to the free-end cross section.
Ten B2 and 8 L9 elements are found to be proper discretization choices to simulate the nonlinear deformation. Figure 14b plots the L9 element distribution in the cross section. In the 0-5% section of the camber line, the curvature of the thickness distribution changes relatively rapidly. Therefore, 3 L9 elements are implemented to construct the cross-sectional displacement field. As for the remaining 95% section, it is appropriate to select 5 L9 elements to describe the smooth camber line and thickness distribution. In the ABAQUS solid model, 5900 hexahedral elements are needed to achieve convergence. Each element contains 8 3-DOF nodes and the total number of DOFs is 26,469, as shown in Figure 14c.

Pretwisted Cantilever with Airfoil Profile Cross Section
A pre-twisted airfoil profile cantilever with NACA 8405 dimensions is investigated. A uniformly distributed longitudinal load (636 kN in total) is applied to the free-end cross section.
Ten B2 and 8 L9 elements are found to be proper discretization choices to simulate the nonlinear deformation. Figure 14b plots the L9 element distribution in the cross section. In the 0-5% section of the camber line, the curvature of the thickness distribution changes relatively rapidly. Therefore, 3 L9 elements are implemented to construct the cross-sectional displacement field. As for the remaining 95% section, it is appropriate to select 5 L9 elements to describe the smooth camber line and thickness distribution. In the ABAQUS solid model, 5900 hexahedral elements are needed to achieve convergence. Each element contains 8 3-DOF nodes and the total number of DOFs is 26,469, as shown in Figure 14c. Comparison results are given for the torsional rotation of the end section and the 3D deformed configuration of the whole structure. The results presented in Figures 15 and 16 show good agreement. The maximum difference in the torsional rotational angle is about 4% with the same applied load. Therefore, the present model shows good accuracy in capturing the nonlinear axial-torsional coupling. The rotational center of the cross section is consistent with the pre-twist center of the structure.  Comparison results are given for the torsional rotation of the end section and the 3D deformed configuration of the whole structure. The results presented in Figures 15 and 16 show good agreement. The maximum difference in the torsional rotational angle is about 4% with the same applied load. Therefore, the present model shows good accuracy in capturing the nonlinear axial-torsional coupling. The rotational center of the cross section is consistent with the pre-twist center of the structure. Comparison results are given for the torsional rotation of the end section and the 3D deformed configuration of the whole structure. The results presented in Figures 15 and 16 show good agreement. The maximum difference in the torsional rotational angle is about 4% with the same applied load. Therefore, the present model shows good accuracy in capturing the nonlinear axial-torsional coupling. The rotational center of the cross section is consistent with the pre-twist center of the structure.

Discussion
This paper presents a new refined beam formulation for geometrically nonlinear twisted structures to accurately predict their large deformations and finite rotati Strain-displacement and stress-strain relationships are constructed in the twisted bas account for the effects of geometrical nonlinearity and initial twist. The secant and tan stiffness matrices of the beam element with geometrical nonlinearities are derived unified way owing to the scalable characteristics of CUF. Three kinds of pre-twisted st tures with different cross sections are investigated to study their large deflection sponses. The specific conclusions are summarized as follows: 1. L9 elements have advanced local mapping capabilities to describe pre-twisted st tures with different cross-sectional geometries. One L9 element is suitable for structing the cross-sectional displacement field of pre-twisted structures with profile sections, while only eight L9 elements can describe the sharp curva change and the displacement field within the cross sections of airfoil profile twisted structures. The maximum difference between the present deformation res and those from commercial simulations is 6%. 2. The stiffness of a pre-twisted structure can be enhanced by increasing its crosstional curvature. For pre-twisted structures with arc profile cross sections, the st tural stiffness can be increased by up to 30% as the arc angle increases from 30 90°. However, this enhancement is reduced for structures with larger pretwist. 3. The untwisted center of the pre-twisted structure coincides with its pre-twisted ter. Moreover, the free-end cross section of a thin-walled, pre-twisted structure k its original profile.
Compared with traditional beam theories and finite element models, the propo method significantly reduces computational complexity and cost by constructing ax invariant elements with separated structural pre-twist. The same discretization ca used to solve large deformation analysis of pre-twisted structures with different pre-t angles. In addition, the developed refined beam model can describe shell-like prope of pre-twisted, thin-walled structures, including coupled axial, torsional, and flexure formations, as well as the local deformations around the point of application of the fo

Discussion
This paper presents a new refined beam formulation for geometrically nonlinear pre-twisted structures to accurately predict their large deformations and finite rotations. Strain-displacement and stress-strain relationships are constructed in the twisted basis to account for the effects of geometrical nonlinearity and initial twist. The secant and tangent stiffness matrices of the beam element with geometrical nonlinearities are derived in a unified way owing to the scalable characteristics of CUF. Three kinds of pre-twisted structures with different cross sections are investigated to study their large deflection responses. The specific conclusions are summarized as follows: 1.
L9 elements have advanced local mapping capabilities to describe pre-twisted structures with different cross-sectional geometries. One L9 element is suitable for constructing the cross-sectional displacement field of pre-twisted structures with arc profile sections, while only eight L9 elements can describe the sharp curvature change and the displacement field within the cross sections of airfoil profile pre-twisted structures. The maximum difference between the present deformation results and those from commercial simulations is 6%.

2.
The stiffness of a pre-twisted structure can be enhanced by increasing its crosssectional curvature. For pre-twisted structures with arc profile cross sections, the structural stiffness can be increased by up to 30% as the arc angle increases from 30 • to 90 • . However, this enhancement is reduced for structures with larger pretwist.

3.
The untwisted center of the pre-twisted structure coincides with its pre-twisted center. Moreover, the free-end cross section of a thin-walled, pre-twisted structure keeps its original profile.
Compared with traditional beam theories and finite element models, the proposed method significantly reduces computational complexity and cost by constructing axially invariant elements with separated structural pre-twist. The same discretization can be used to solve large deformation analysis of pre-twisted structures with different pre-twist angles. In addition, the developed refined beam model can describe shell-like properties of pre-twisted, thin-walled structures, including coupled axial, torsional, and flexure deformations, as well as the local deformations around the point of application of the force.
Author Contributions: Y.H.: Conceptualization, methodology, software, validation, investigation, resources, data curation, writing-original draft preparation, visualization; Y.Z.: Conceptualization, software, formal analysis, investigation, resources, data curation, writing-review and editing, visualization, supervision, project administration, funding acquisition; H.L.: Software, validation, investigation, resources, data curation, writing-review and editing, visualization. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
Linear operator matrix of the Green-Lagrange strain vector in the contravariant basis Nonlinear operator matrix of the Green-Lagrange strain vector in the contravariant basis b cov,nl = 1 2 Linear operator matrix of the Green-Lagrange strain vector in the twisted basis Nonlinear operator matrix of the Green-Lagrange strain vector in the twisted basis where, Λ = ∂ y + φη∂ ξ − φξ∂ η .
In the framework of CUF, linear operator matrix of the strain vector in the twisted basis In the framework of CUF, nonlinear operator matrices of the strain vector in the twisted basis u ξ,ξ F s,ξ N j u y,ξ F s,ξ N j u η,ξ F s,ξ N j φ 2 u ξ F s N j + Λu ξ Λ F s N j Λu y Λ F s N j φ 2 u η F s N j + Λu η Λ F s N j u ξ,η F s,η N j u y,η F s,η N j u η,η F s,η N j u ξ,ξ F s,η N j + u ξ,η F s,ξ N j u y,ξ F s,η N j + u y,η F s,ξ N j u η,ξ F s,η N j + u η,η F s,ξ N j Λu ξ F s,η N j + u ξ,η Λ F s N j Λu y F s,η N j + u y,η Λ F s N j Λu η F s,η N j + u η,η Λ F s N j u ξ,ξ Λ F s N j + Λu ξ F s,ξ N j u y,ξ Λ F s N j + Λu y F s,ξ N j u η,ξ Λ F s N j + Λu η F s,ξ N j with Λ F s N j = F s N j,y − φξF s,η N j + φηF s,ξ N j (A8) Operator matrices of the linearization for the virtual variation of the strain: F τ,ξ F s,ξ N i N j F τ,ξ F s,ξ N i N j F τ,ξ F s,ξ N i N j φ 2 + Λ 2 F τ N i F s N j Λ 2 F τ N i F s N j φ 2 + Λ 2 F τ N i F s N j F τ,η F s,η N i N j F τ,η F s,η N i N j F τ,η F s,η N i N j F τ,ξ F s,η N i N j + F τ,η F s,ξ N i N j F τ,ξ F s,η N i N j + F τ,η F s,ξ N i N j F τ,ξ F s,η N i N j + F τ,η F s,ξ N i N j Λ(F τ N i )F s,η N j + F τ,η N i Λ F s N j Λ(F τ N i )F s,η N j + F τ,η N i Λ F s N j Λ(F τ N i )F s,η N j + F τ,η N i Λ F s N j F τ,ξ N i Λ F s N j + Λ(F τ N i )F s,ξ N j F τ,ξ N i Λ F s N j + Λ(F τ N i )F s,ξ N j F τ,ξ N i Λ F s N j + Λ(F τ N i )F s,ξ N j