Exact Solutions for Torsion and Warping of Axial-Loaded Beam-Columns Based on Matrix Stiffness Method

The typically-used element torsional stiffness GJ/L (where G is the shear modulus, J the St. Venant torsion constant, and L the element length) may severely underestimate the torsional stiffness of thin-walled nanostructural members, due to neglecting element warping deformations. In order to investigate the exact element torsional stiffness considering warping deformations, this paper presents a matrix stiffness method for the torsion and warping analysis of beam-columns. The equilibrium analysis of an axial-loaded torsion member is conducted, and the torsion-warping problem is solved based on a general solution of the established governing differential equation for the angle of twist. A dimensionless factor is defined to consider the effect of axial force and St. Venant torsion. The exact element stiffness matrix governing the relationship between the element-end torsion/warping deformations (angle and rate of twist) and the corresponding stress resultants (torque and bimoment) is derived based on a matrix formulation. Based on the matrix stiffness method, the exact element torsional stiffness considering the interaction of torsion and warping is derived for three typical element-end warping conditions. Then, the exact element second-order stiffness matrix of three-dimensional beam-columns is further assembled. Some classical torsion-warping problems are analyzed to demonstrate the established matrix stiffness method.


Introduction
In recent decades, research on nanomaterials (i.e., materials with internal structure of nanoscale dimension) has made great progress and is widely used in scientific research and industrial production. Some formations of oxides, metals, ceramics, and other substances have been discovered. These nanomaterials obey the fundamental laws of the classical physics governing the macroworld [1]. Nanomaterials such as graphene may provide many enhanced properties including high strength, stiffness, and light weight [2]. Researchers also utilize nanomaterials in structures to improve the mechanical behavior and other performance [3].
Although a number of computer programs such as Abaqus, OpenSees, Ansys, and MSC.Marc are readily available for the structural analysis of thin-walled nanostructural members, approaches to obtain the exact solutions [4][5][6][7][8][9][10] in closed form are helpful in many situations. The matrix stiffness method (MSM) has been found to be a suitable and systematic method for such purposes [11][12][13][14][15][16][17][18][19][20]. The basic idea of the matrix stiffness method is to establish the equilibrium relationship between the element-end displacements ∆ = (u 1 , δ 1 , θ 1 , u 2 , δ 2 , θ 2 ) T and the element-end forces F ext = (F 1 , V 1 , M 1 , F 2 , V 2 , M 2 ) T of a beam-column element (where u i , d i , and q i are element-end axial displacement, translational [K e ]∆ = F ext (1) where [K e ] is the element stiffness matrix for flexural-axial problems. For first-order (e.g., [21,22]) and second-order (e.g., [10,[13][14][15][16]19]) analysis (i.e., without and with considering the geometric nonlinearity), the element stiffness matrix of a beamcolumn element can be formulated as Equations (2)  [K e,2nd-ord ] = where E denotes the elastic modulus; A denotes the cross-sectional area; I denotes the cross-sectional moment of inertia; L denotes the element length; and T c , Q c , S c , and C c are coefficients in the element stiffness matrix as functions of factor λ, formulated as Equation (4); λ (in a flexural-axial problem) denotes a dimensionless factor for the axial compression force P c , defined as Equation (5).
The matrix stiffness method can also be used for three-dimensional analysis of beamcolumns. For second-order analysis, Ekhande [23] presented an exact element stiffness matrix associated with the 12 element-end deformations/rotation angles of 3D beamcolumns as Equation (6). The matrix has been used to conduct stability analysis of 3D structures [19].
Sc λy where G denotes the shear modulus; J denotes the St. Venant torsion constant; I yy and I zz denote the cross-sectional moments of inertia about y-axis and z-axis, respectively; λ y and λ z denote axial force factors associated with the y-axis bending and z-axis bending, respectively; u yi and u zi denote element-end translational displacements in y-direction and z-direction, respectively; ϕ i denotes the element-end angle of twist; θ yi and θ zi denote element-end rotation angles in the X-Z plane and X-Y plane, respectively; F yi and F zi denote element-end loads in the y-direction and z-direction, respectively; T i denotes the elementend torsional moment; and M yi and M zi denote element-end bending moments about the y-axis and z-axis, respectively, as shown in Figure 2. However, these researchers neglected the interaction between torsion and warping and the axial force effect in torsional analysis. They directly used GJ/L as torsional stiffness in engineering practices, which may severely underestimate the torsional stiffness of thin-walled nanostructural members without considering warping deformations. In order to present the exact element stiffness matrix for the second-order analysis of three-dimensional beam-columns considering torsion and warping, this paper investigates the interaction of the axial force, torque, and bimoment of torsion members. The exact element torsional stiffnesses are derived and shown to be significantly larger than the typically-used GJ/L. The exact element stiffness matrix associated with the element-end torsion and warping deformations can be obtained, and the exact three-dimensional element stiffness matrix can be further assembled. Some application examples of the exact element stiffness matrix for torsion and warping are also presented.
In addition, the temperature is non-negligible, which can also influence the bifurcation buckling of thin-walled structures, and the thermoelastic analysis of structures has already been presented in many papers. However, this paper is focused on solving the mechanical buckling of thin-walled members. Therefore, we list some papers [24,25] on thermoelastic analysis as a reference instead of presenting further research.

Equilibrium Analysis
To develop the matrix stiffness method, an equilibrium analysis of an axial-loaded torsion member (especially for members with a thin-walled cross-section) is conducted, as shown in Figure 3. Besides the usual assumptions of the linear theory of elasticity, the following assumptions [5,6,10,[26][27][28] of classical theory for members with a thin-walled cross-section are employed in the analysis:

•
The global cross-sectional deformation assumption, i.e., the cross-section is assumed to be perfectly rigid in its own plane while free to warp out of its plane; • The classical Kirchhoff-Love's thin plate bending assumption, i.e., straight lines normal to the mid-surface of the thin-walled plates remain straight, inextensible, and normal to the mid-surface after deformation; • Each thin-walled plate is assumed to have null mid-surface membrane shear strains (Vlasov's hypothesis) and null transverse extensions.
Based on these assumptions, a compatibility equation relating the angle of twist ϕ and the cross-sectional bimoment B ω is formulated as Equation (7). The St. Venant torque T s is formulated as Equation (8) [5,6,10,26].
T s = GJ dϕ dx (8) where I ωω denotes the warping constant. The sign conventions of B ω and T are defined in Figure 3b. Since the total cross-sectional torque T is consisted of the T s and the warping restraint torque T w , the T w can then be derived as the total cross-sectional torque T subtracted by the St. Venant torque 1.

Equilibrium of torque for element short segment
For an equilibrium analysis, a short segment of the element with length dx is analyzed in Figure 3b. The equilibrium of torque gives where τ denotes the distributed torque.

Equilibrium of bimoment for element short segment
The equilibrium of bimoment is also analyzed. The warping restraint torques [Equation (9)] at the two cross-sections with distance x of the element short segment are in the opposite direction, and they combine to a bimoment increment formulated as In addition, due to the torsion of the short segment (cross-section elevation view in Figure 3c), the uniformly distributed axial stress σ n on the top section is inclined to produce a shear stress σ n dAρ s dϕ in the horizontal plane. Therefore, the Wagner effect can be considered by taking moment of the shear stress about the shear center S, derived as Equation (12).
where σ n = P/A denotes the assumed uniformly-distributed axial stress from the axial force; ρ S denotes the distance of a point in the cross-section to the shear center S; and I p,S denotes the polar moment of inertia about the shear center S. Therefore, considering these two effects and an increment of the cross-sectional bimoment along the element length, the equilibrium of bimoment is formulated as Equation (13).
Then, combining Equations (10) and (13) and considering Equation (7), the governing differential equation of equilibrium for this torsion-warping problem can be established as Equation (14), which can be derived to a dimensionless form as Equation (15).
where λ c/t denotes a dimensionless factor for the effect of axial force and St. Venant torsion, defined as Equation (16); the subscripts "c" and "t" are associated with conditions of PI p,S /A ≥ GJ and PI p,S /A < GJ, which can be defined as generalized "axial compression" and "axial tensile" situations, respectively.
For null distributed torque (τ = 0), the general solution for the differential equation of equilibrium is given by where Q 1 , Q 2 , Q 3 , and Q 4 (Q 1t , Q 2t , Q 3t , Q 4t ) are deformation combination factors defining the possible deformation curve. In the following, the derivations will focus on the generalized "axial compression" situation. It is noted that the derivations and the results for the generalized "axial tensile" situation is very similar to the "axial compression" situation, and the main difference is the use of hyperbolic trigonometric functions instead of trigonometric functions [as shown in Equation (17)].
Analysis associated with the element deformations and stress resultants of torsion members can then be conducted using Equations (17)- (19).

Matrix Stiffness Method for Torsion and Warping
An element stiffness matrix, showing the relationship between the element-end deformations (angle and rate of twist) and the corresponding stress resultants (torque and bimoment), is formulated based on the previous section. Then, a matrix stiffness method for torsion and warping is established for the analysis of torsion members, including the approximate element torsion-warping stiffness matrix for simpler applications and the torsional stiffness analysis for three typical element-end warping conditions.

Element Stiffness Matrix for Torsion and Warping
Equations (17)- (19) can be used in the analysis associated with the element deformations and stress resultants of torsion members. In a matrix stiffness method, attentions are focused on the element-end torsion/warping deformations and the corresponding stress resultants, formulated in matrix forms for simplified and systematic analysis.
The element-end deformations (angle and rate of twist) are formulated in a matrix form as Equation (20) based on Equation (17).
The element-end stress resultants (torque and bimoment) are formulated in a matrix form as Equation (21) based on Equations (18) and (19).
The relationship between the element-end deformations (angle and rate of twist) and the corresponding stress resultants (torque and bimoment) is then formulated as Equation (22) by combining Equations (20) and (21).
where [K e ] is the element stiffness matrix for torsion and warping, which can be simplified as where the expressions for the element stiffness coefficients T c , Q c , S c , and C c are the same as that in the element stiffness matrix for a flexural-axial problem, which are formulated in Equation (4).
Equation (23) gives the element stiffness matrix based on the axial (compression) force factor λ c (Equation (16)) in the case of PI p,S /A ≥ GJ.
In the case of PI p,S /A < GJ, the element stiffness matrix can also be derived from the general solution in Equation (17) and is then formulated in the same form as Equation (23), with a change of the subscript "c" to "t". However, the element stiffness coefficients T t , Q t , S t , and C t are expressed in a form using hyperbolic trigonometric functions based on the axial (tension) force factor λ t .
The element stiffness functions T c , Q c , S c , and C c correspond to the element-end stress resultant (torque or bimoment) for a unit element-end deformation (angle or rate of twist). These coefficients are transcendental functions of the factor λ c for St. Venant torsion and axial force. By noting that λ c = PI p,S /A − GJ /(EI ωω /L 2 ), the influences of GJ−PI p,S /A on these element stiffness functions are plotted in Figure 4. As shown in Figure 4, T c , Q c , and S c increase (while C c decreases) with the increase in GJ−PI p,S /A (which corresponds to either an increasing GJ or a decreasing axial force P).  Figure 4, which shows that these linear approximations have considerably small differences from the exact transcendental element stiffness functions when GJ−PI p,S /A is in a range between −π 2 EI ωω /L 2 and π 2 EI ωω /L 2 .
Based on Equation (25), the element torsion-warping stiffness matrix [K e ] can be approximated as Equation (26) for simpler applications, which shows the linear influences of the warping constant I ωω , the St. Venant torsion constant J, and the axial force P on [K e ]. This approximated element stiffness matrix can also be derived by using an energy approach and assuming a cubic deformation shape function [11].

Torsional Stiffnesses for Three Typical Element-End Warping Conditions
In engineering practices, the element torsional stiffness has drawn most of the attentions, but the warping stiffness and the influence of element-end warping restraint on the torsional stiffness have not been sufficiently investigated.
The element torsional stiffness is usually expressed using the St. Venant torsion constant as GJ/L. However, because of the interaction between torsion and warping, the element-end torsion stiffness may vary for different warping conditions. Therefore, an element-end torsional stiffness matrix considering different warping conditions is required.
By using the element torsion-warping stiffness matrix [K e ], we can derive the torsional stiffnesses of members with typical element-end warping conditions.
For members with restrained warping at the ends (ω 1 = ω 2 = 0), the rows and columns in [K e ] related to the rate of twist and bimoment can be deleted to obtain a torsion stiffness matrix.
where the torsional stiffness k tor for this restrained-restrained warping condition can be expressed as For members with no restraint of warping at the ends (B ω 1 = B ω 2 =0), the torsion stiffness matrix that relates the element-end twisting angles and torques can be derived as follows.
By using the element stiffness functions in Equation (23), the torsional stiffness k tor for this free-free warping condition is derived as For members with restrained warping at one end (ω 1 = 0) and free warping at the other end (B ω 2 = 0), the torsion stiffness matrix can be derived as follows.
By using the element stiffness functions in Equation (23), the torsional stiffness k tor for this restrained-free warping condition is derived as where T F is a stiffness function associated with this restrained-free-warping condition, and its linear approximation is derived as Equation (33). The stiffness function T F and its approximation are also plotted in Figure 4.
In summary, the torsional stiffnesses of members with three typical element-end warping conditions can be formulated as where the linear approximations are based on Section 3.1. Equation (34) shows that the commonlyused expression GJ/L for the torsional stiffness is only valid for members with the free-warping condition and negligible axial force effect. Figure 5 compares the torsional stiffnesses for these three typical warping conditions. The effect of St. Venant torsion and axial force is varied along the horizontal axis. As shown in Figure 5, the torsional stiffness associated with the restrained-restrained warping condition may be significantly larger than the commonly-used value GJ/L. Therefore, in structural analysis, the torsional stiffness value should be carefully selected based on the element-end warping conditions.

Applications of Torsion-Warping Stiffness Matrix
In this section, some classical torsion-warping problems will be analyzed as application examples to demonstrate the established matrix stiffness method.

Elastic Buckling Analysis for Typical Torsion-Warping-Axial Stability Problems
The elastic bifurcation buckling analysis can be directly conducted using the assembled structural stiffness matrix [K s ]. The eigenproblem for an elastic buckling analysis tries to solve a nontrivial solution for the global load deformation equation [K s ]∆ = 0. Therefore, an elastic buckling analysis can be conducted by setting the determinant of the assembled structural stiffness matrix [K s ] to zero (Equation (35)). The factor λ c,cr at the buckling state can be solved, corresponding to a buckling axial force P cr formulated as Equation (36).
For torsion-warping-axial stability problems with typical element-end boundary conditions, the elastic bifurcation buckling analysis can be conveniently conducted using the matrix stiffness method, as listed in Table 1.
For a torsion restrained column (ϕ 1 = 0, ϕ 2 = 0) with free of warping at the two element ends (case 1), the unrestrained element-end deformations are the warping deformations. Therefore, the buckling state corresponds to the condition that the determinant of the submatrix for warping (third and fourth row/column of the element stiffness matrix [K e ]) equals zero, and the buckling condition in this case can be solved as For a column with torsion and warping restraints only at one node (ϕ 1 = 0, ω 1 = 0) (case 2), the unrestrained element-end deformations are the torsion and warping deformations at the other node. Therefore, the buckling state corresponds to the condition that the determinant of the stiffness matrix at the unrestrained node (second and fourth row/column of the element stiffness matrix [K e ]) equals zero, and the buckling condition in this case can be solved as For a column with torsion and warping restraints at one node (ϕ 1 = 0, ω 1 = 0) and a warping restraint at the other node (ω 2 = 0) (case 3), the only unrestrained element-end deformation is the torsion deformation at the other node. The relationship between this unrestrained torsion deformation and the corresponding element-end torque is defined by T c (λ c ). Therefore, the buckling state corresponds to the condition that T c (λ c ) = 0. Based on Figure 4, T c decreases to 0 as λ c increases to π. Therefore, the buckling condition in this case can be solved as For a column with torsion and warping restraints at one node (ϕ 1 = 0, ω 1 = 0) and a torsion restraint at the other node (ϕ 2 = 0) (case 4), the only unrestrained element-end deformation is the warping deformation at the other node. The relationship between this unrestrained warping deformation and the corresponding element-end bimoment can be defined by S c (λ c ). Therefore, the buckling state is corresponding to the condition that S c (λ c ) = 0. Based on Figure 4, S c decreases to 0 as λ c increases to 1.43π (π/0.7). Therefore, the buckling condition in this case can be solved as In addition, for a column with torsion and warping restraints at both nodes (case 5), the buckling state is corresponding to the condition that the denominator Φ c in T c , Q c , S c , and C c equals zero [15]. It can be solved that Φ c decreases to 0 as λ c increases to 2π. Therefore, the buckling condition in this case can be solved as It is noted that this section gives the same results as that from classical analyses of these torsion-warping-axial buckling problems [4,5,7,8], but the matrix analysis procedure is considerably simplified and is more systematic.

Analysis of Torsion and Warping of a Torsion Member with a Midspan Torque
This example analyzes a classical torsion-warping problem [4,5], the torsion and warping of a torsion member with a midspan torque, as shown in Figure 6. In nanostructures, torsion members may also be used. Therefore, the matrix stiffness method could be relevant to structural nanomechanics and suitable for nanostructures. In this analysis, the different torque components of the total torque (including the St. Venant torque and the warping restraint torque) are discussed. The governing equation is formulated as Equation (42), where [K s ] is the structural stiffness matrix associated with the unconstrained deformations ϕ 2 , ω 1 , ω 2 , and ω 3 . The torsion member is considered to consist of two elements with length L (elements I and II in Figure 6). For the two elements, the relationships between their element end deformations (ϕ 1 , ϕ 2 , ω 1 , and ω 2 ) (or (ϕ 2 , ϕ 3 , ω 2 , and ω 3 )) and the corresponding element end stress resultants (T I1 , T I2 , B ω I1 , and B ω I2 ) (or (T II1 , T II2 , B ω II1 , and B ω II2 )) are both governed by the element stiffness matrix [K e ] [Equation (23)]. By combining the stability stiffness matrices of element I (row/column 2 to 4 of [K e ] associated with the unconstrained element end deformations, ϕ 2 , ω 1 , and ω 2 ) and element II (row/column 1, 3, and 4 of [K e ] associated with the unconstrained element end deformations, ϕ 2 , ω 2 , and ω 3 ), the structural stiffness matrix [K s ] is formulated as Equation (43).
Then, the element-end torsion/warping deformations can be solved based on the relationship in Equation (42). The angle of twist at midspan is For this problem particularly, the element analysis can be conducted for the element I or II. The integration of different torque components (St. Venant torque and warping restraint torque) is analyzed [4,5]. The integration of the St. Venant torque T s from the element end to the midspan can be derived (based on Equation (8)) as the product of the St. Venant torsion rigidity GJ and the angle of twist at midspan ϕ 2 , as shown in Equation (45). The integration of warping restraint torque T w from the element end to the midspan gives the bimoment at midspan B ω ,mid [based on Equation (13)], as shown in Equation (46). The integration of cross-sectional total torque (T = T mid ) from the element end to the midspan can be formulated as Equation (47).   It is noted that this section gives the same results as that from the analysis in Trahair et al. [4] and Chen [5], but the matrix analysis procedure is considerably simplified and is more systematic.

Second-Order Stiffness Matrix of 3D Beam Considering Exact Torsional Stiffness
Based on the torsion-warping stiffness matrix, this section establishes the element stiffness matrix of axial-loaded three-dimensional beam-columns with a symmetric cross-section.
The element-end displacement vector of a three-dimensional beam-column with a symmetric cross-section is considered in Equation (50) to include the displacements, the torsional and rotational angles, as well as the rates of twist, as defined in Figure 8.
where (x, y, z) denote the local coordinate systems for the three-dimensional element. The associated element-end stress resultant vector is considered in Equation (51) to include the forces, bending moments and torques, as well as the bimoments.
The element stiffness matrix of three-dimensional beam-columns relates the element-end displacement vector in Equation (50) to the element-end stress resultant vector in Equation (51), and it is, therefore, a 14-degree-of-freedom stiffness matrix. McGuire et al. [11] noted that the difference between the analyses of planar system and three-dimensional system is essentially quantitative. By considering (1) the stiffness matrix [Equation (6)] associated with the member bending in both the x-y and x-z planes and (2) the torsion-warping stiffness matrix (Equation (23))established in this paper, the element stiffness matrix of three-dimensional beam-columns can be obtained: In view of Equation (52), the stiffnesses associated with the y-axis bending, the z-axis bending, and the torsion and warping are uncoupled. Therefore, row/column 2, 6, 8, 12 of Equation (52) represents the z-axis bending of the element, row/column 3, 5, 9, 11 of Equation (52) represents the y-axis bending, and row/column 4, 10, 13, 14 of Equation (52) represents the torsion and warping.
For beam-columns with the three typical element-end warping conditions discussed in Section 3.2, the 14-degree-of-freedom element stiffness matrix of three-dimensional beam-columns can be reduced to a 12-degree-of-freedom element stiffness matrix. This 12-degree-of-freedom element stiffness matrix is more typically used in a systematic analysis of three-dimensional frame systems [11,15,23] because it can be easily transformed to the element global stiffness matrix by using a transformation matrix relating the local and global coordinate systems.
In view of Equation (53), the torsional stiffness k tor in the 4 and 10 rows/columns should be determined using Section 3.2 based on the element-end warping conditions (instead of directly using the value GJ/L).
This 3D element stiffness matrix can be readily used to solve the exact solutions of the bifurcation buckling problem of 3D frames as well as the out-of-plane buckling of funicular arches.

Conclusions
This paper presented a matrix stiffness method for the analysis of torsion and warping that is particularly important in beam-columns with torsional deformations. The main works and conclusions are summarized as follows: 1.
Equilibrium analysis of an axial-loaded torsion member was conducted based on the equilibrium conditions of torque and bimoment of an element short segment, and a governing differential equation of equilibrium for the angle of twist along the member was established. The solution of the governing differential equation can be used to analyze the element deformations (angle and rate of twist) and stress resultants (torque and bimoment).

2.
The exact element stiffness matrix of the torsion member was formulated, showing the relationship between the element-end torsion/warping deformations and the corresponding stress resultants. A dimensionless factor for the effect of St. Venant torsion and axial force was defined. The element stiffness matrix for torsion and warping was readily used for the bifurcation buckling and second-order analysis of axial-loaded torsion members. The exact element second-order stiffness matrix of three-dimensional beam-columns was further assembled.

3.
Based on the element torsion-warping stiffness matrix, the exact element torsional stiffnesses considering the interaction of torsion and warping were derived for three typical element-end warping conditions. The commonly-used expression GJ/L for the torsional stiffness is only valid for members with the free-warping condition and negligible axial force effect. For members with restrained warping at the ends, the torsional stiffness can be significantly larger. For a notable axial force effect, the torsional stiffness may be reduced. Therefore, in the analysis of thin-walled nanostructural structures, the torsional stiffness value should be carefully selected based on the element-end warping conditions and the axial force level. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data presented in this article is available on request from corresponding author.

Acknowledgments:
The authors express their sincere appreciation to the reviewers of this paper for their constructive comments and suggestions.

Conflicts of Interest:
The authors declare no conflict of interest. Bending moments about the y-axis and z-axis, respectively P Axial compressive force P cr Buckling axial force Q 1 , Q 2 , Q 3 , Q 4 (Q 1t , Q 2t , Q 3t , Q 4t ) Deformation combination factors defining the possible deformation curve T c/t , Q c/t , S c/t , C c/t , T F , ϕ c/t Coefficients in the element stiffness matrix as functions of λ c/t T mid Applied midspan torque T s , T w St. Venant torque and warping restraint torque, respectively T(x), B ω (x)

Nomenclature
Torque and bimoment at the cross-section considered, respectively T 1/2 , B ω1 /2 Torques and bimoments at the two element ends, respectively u y1/2 , u z1/2 Translational displacements in y-direction and z-direction, respectively ∆, F ext Element-end deformation vector and corresponding stress resultant vector θ y1/2 , θ z1/2 Rotation angles in the X-Z plane and X-Y plane, respectively λ c/t Factor for the effect of axial force and St. Venant torsion λ c,cr Factor for axial force and St. Venant torsion at the buckling state λ y , λ z Factors for axial force associated with the y-axis bending and z-axis bending, respectively ρ S Distance of a point in the cross-section to the shear center σ n Assumed uniformly-distributed axial stress from the axial force τ Distributed torque ϕ,ω Angle and rate of twist, respectively ϕ 1/2 , ω 1/2 Angles and rates of twist at the two element ends, respectively