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

Abstract

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.

1. 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 $\mathsf{\Delta }={\left(u_1,\hspace{0.33em}{\delta }_1,\hspace{0.33em}{\theta }_1,\hspace{0.33em}{u}_2,\hspace{0.33em}{\delta }_2,\hspace{0.33em}{\theta }_2\right)}^{\mathrm{T}}$ and the element-end forces ${\mathit{F}}_{ext}={\left(F_1,\hspace{0.33em}{V}_1,\hspace{0.33em}{M}_1,\hspace{0.33em}{F}_2,\hspace{0.33em}{V}_2,\hspace{0.33em}{M}_2\right)}^{\mathrm{T}}$ of a beam-column element (where u_i, d_i, and q_i are element-end axial displacement, translational displacement, and rotation angle, respectively; F_i, V_i, and M_i are element-end axial load, shear force, and bending moment, respectively, as shown in Figure 1) as

$$\left[{\mathit{K}}_e\right]\Delta ={\mathit{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 beam-column element can be formulated as Equations (2) and (3), respectively.

$$\left[{\mathit{K}}_{\mathbf{e},\mathbf{1}\mathbf{st}\mathbf{-}\mathbf{ord}}\right]=\left(\begin{array}{cccccc}\frac{EA}{L}& & & -\frac{EA}{L}& & \\ & 12\frac{EI}{L^3}& 6\frac{EI}{L^2}& & -12\frac{EI}{L^3}& 6\frac{EI}{L^2}\\ & & 4\frac{EI}{L}& & -6\frac{EI}{L^2}& 2\frac{EI}{L}\\ & & & \frac{EA}{L}& & \\ & & & & 12\frac{EI}{L^3}& -6\frac{EI}{L^2}\\ \mathrm{sym}.& & & & & 4\frac{EI}{L}\end{array}\right)$$

(2)

$$\left[{\mathit{K}}_{\mathbf{e}\mathbf{,}\mathbf{2}\mathbf{nd}\mathbf{-}\mathbf{ord}}\right]=\left(\begin{array}{cccccc}\frac{EA}{L}& & & -\frac{EA}{L}& & \\ & T_{\mathrm{c}}(\lambda )\frac{EI}{L^3}& Q_{\mathrm{c}}(\lambda )\frac{EI}{L^2}& & -T_{\mathrm{c}}(\lambda )\frac{EI}{L^3}& Q_{\mathrm{c}}(\lambda )\frac{EI}{L^2}\\ & & S_{\mathrm{c}}(\lambda )\frac{EI}{L}& & -Q_{\mathrm{c}}(\lambda )\frac{EI}{L^2}& C_{\mathrm{c}}(\lambda )\frac{EI}{L}\\ & & & \frac{EA}{L}& & \\ & & & & T_{\mathrm{c}}(\lambda )\frac{EI}{L^3}& -Q_{\mathrm{c}}(\lambda )\frac{EI}{L^2}\\ \mathrm{sym}.& & & & & S_{\mathrm{c}}(\lambda )\frac{EI}{L}\end{array}\right).$$

(3)

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).

$$\begin{array}{l}{T}_{\mathrm{c}}(\lambda )={\lambda }^3\sin \lambda /{\varphi }_{\mathrm{c}}(\lambda )\hfill \\ Q_{\mathrm{c}}(\lambda )={\lambda }^2\left(1-\cos \lambda \right)/{\varphi }_{\mathrm{c}}(\lambda )\hfill \\ S_{\mathrm{c}}(\lambda )=\lambda \left(\sin \lambda -\lambda \cos \lambda \right)/{\varphi }_{\mathrm{c}}(\lambda )\hfill \\ C_{\mathrm{c}}(\lambda )=\lambda \left(\lambda -\sin \lambda \right)/{\varphi }_{\mathrm{c}}(\lambda )\hfill \\ {\varphi }_{\mathrm{c}}(\lambda )=2-2\cos \lambda -\lambda \sin \lambda \hfill \end{array}$$

(4)

$$\lambda =\sqrt{\frac{P_{\mathrm{c}}{L}^2}{EI}}=\pi \sqrt{\frac{P_{\mathrm{c}}}{P_{\mathrm{E}}}}$$

(5)

Using the element stiffness matrix, many different types of analyses can be conducted, including the traditional analysis for the element deformations and internal forces (e.g., [21,22]), as well as elastic buckling and second-order stability analyses (e.g., [15,16,17,18]).

The matrix stiffness method can also be used for three-dimensional analysis of beam-columns. For second-order analysis, Ekhande [23] presented an exact element stiffness matrix associated with the 12 element-end deformations/rotation angles of 3D beam-columns as Equation (6). The matrix has been used to conduct stability analysis of 3D structures [19].

$$\left(\begin{array}{cccccccccccc}\frac{EA}{L}& & & & & & -\frac{EA}{L}& & & & & \\ & T_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^3}& & & & Q_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^2}& & -T_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^3}& & & & Q_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^2}\\ & & T_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^3}& & -Q_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^2}& & & & -T_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^3}& & -Q_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^2}& \\ & & & \frac{GJ}{L}& & & & & & -\frac{GJ}{L}& & \\ & & & & S_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L}& & & & Q_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^2}& & C_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L}& \\ & & & & & S_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L}& & -Q_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^2}& & & & C_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L}\\ & & & & & & \frac{EA}{L}& & & & & \\ & & & & & & & T_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^3}& & & & -Q_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^2}\\ & & & & & & & & T_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^3}& & Q_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^2}& \\ & & & & & & & & & \frac{GJ}{L}& & \\ & \mathrm{sym}.& & & & & & & & & S_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L}& \\ & & & & & & & & & & & S_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L}\end{array}\right)\left(\begin{array}{c}{u}_{x1}\\ u_{y1}\\ u_{z1}\\ {\phi }_1\\ {\theta }_{y1}\\ {\theta }_{z1}\\ u_{x2}\\ u_{y2}\\ u_{z2}\\ {\phi }_2\\ {\theta }_{y2}\\ {\theta }_{z2}\end{array}\right)=\left(\begin{array}{c}{F}_{x1}\\ F_{y1}\\ F_{z1}\\ T_1\\ M_{y1}\\ M_{z1}\\ F_{x2}\\ F_{y2}\\ F_{z2}\\ T_2\\ M_{y2}\\ M_{z2}\end{array}\right)$$

(6)

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; ${\phi }_i$ denotes the element-end angle of twist; ${\theta }_{yi}$ and ${\theta }_{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 element-end 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.

2. 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].

$$B_{\omega }=-E{I}_{\omega \omega }\frac{d^2\phi }{d{x}^2}$$

(7)

$$T_s=GJ\frac{d\phi }{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

$$T_{\mathrm{w}}=T-T_s=T-GJ\frac{d\phi }{dx}$$

(9)

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

$$\frac{dT}{dx}=-\tau$$

(10)

where τ denotes the distributed torque.

2.

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

$$T_{\mathrm{w}}dx=\left(T-GJ\frac{d\phi }{dx}\right)dx=Tdx-GJd\phi$$

(11)

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 σ_ndAρ_sdφ 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).

$$\underset{A}{\int }\left({\sigma }_{\mathrm{n}}dA\cdot {\rho }_{S}d\phi \cdot {\rho }_S\right)={\sigma }_{\mathrm{n}}d\phi {\displaystyle \underset{A}{\int }{\rho }_S{}^{2}dA}=\frac{P}{A}{I}_{\mathrm{p},S}d\phi$$

(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).

$$Tdx-GJd\phi +\frac{P}{A}{I}_{\mathrm{p},S}d\phi =d{B}_{\omega }\Rightarrow T=\frac{d{B}_{\omega }}{dx}-\left(P{I}_{\mathrm{p},S}/A-GJ\right)\frac{d\phi }{dx}$$

(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).

$$E{I}_{\omega \omega }\frac{d^4\phi }{d{x}^4}+\left(P{I}_{\mathrm{p},S}/A-GJ\right)\frac{d^2\phi }{d{x}^2}=\tau$$

(14)

$$\begin{array}{cc}\frac{d^4\phi }{d{x}^4}+{\frac{{\lambda }_{\mathrm{c}}}{L^2}}^2\frac{d^2\phi }{d{x}^2}=\frac{\tau }{E{I}_{\omega \omega }}& \mathrm{for} P{I}_{\mathrm{p},S}/A\ge GJ\\ \frac{d^4\phi }{d{x}^4}-{\frac{{\lambda }_{\mathrm{t}}}{L^2}}^2\frac{d^2\phi }{d{x}^2}=\frac{\tau }{E{I}_{\omega \omega }}& \mathrm{for} P{I}_{\mathrm{p},S}/A<GJ\end{array}$$

(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.

$$\begin{array}{cc}{\lambda }_{\mathrm{c}}=\sqrt{\frac{P{I}_{\mathrm{p},S}/A-GJ}{E{I}_{\omega \omega }/L^2}}& \mathrm{for} P{I}_{\mathrm{p},S}/A\ge GJ\\ {\lambda }_{\mathrm{t}}=\sqrt{\frac{GJ-P{I}_{\mathrm{p},S}/A}{E{I}_{\omega \omega }/L^2}}& \mathrm{for} P{I}_{\mathrm{p},S}/A<GJ\end{array}$$

(16)

For null distributed torque (τ = 0), the general solution for the differential equation of equilibrium is given by

$$\phi =\left\{\begin{array}{cc}{Q}_1\cos \left({\lambda }_{\mathrm{c}}x/L\right)+Q_2\sin \left({\lambda }_{\mathrm{c}}x/L\right)+Q_{3}x+Q_4& \mathrm{for} P{I}_{\mathrm{p},S}/A\ge GJ\\ Q_{1\mathrm{t}}\cosh \left({\lambda }_{\mathrm{t}}x/L\right)+Q_{2\mathrm{t}}\sinh \left({\lambda }_{\mathrm{t}}x/L\right)+Q_{3\mathrm{t}}x+Q_{4\mathrm{t}}& \mathrm{for} P{I}_{\mathrm{p},S}/A<GJ\end{array}\right.$$

(17)

where Q₁, Q₂, Q₃, and Q₄ (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)].

Based on Equations (7) and (13), the element stress resultants (bimoment and torque) can be derived as Equations (18) and (19), respectively.

$$\frac{B_{\omega }}{E{I}_{\omega \omega }}=-\frac{d^2\phi }{d{x}^2}=Q_1{\left(\frac{{\lambda }_{\mathrm{c}}}{L}\right)}^2\cos \frac{{\lambda }_{\mathrm{c}}x}{L}+Q_2{\left(\frac{{\lambda }_{\mathrm{c}}}{L}\right)}^2\sin \frac{{\lambda }_{\mathrm{c}}x}{L}$$

(18)

$$\frac{T}{E{I}_{\omega \omega }}=\frac{d\left(B_{\omega }/E{I}_{\omega \omega }\right)}{dx}-{\left(\frac{{\lambda }_{\mathrm{c}}}{L}\right)}^2\frac{d\phi }{dx}=-{\left(\frac{{\lambda }_{\mathrm{c}}}{L}\right)}^2{Q}_3$$

(19)

Analysis associated with the element deformations and stress resultants of torsion members can then be conducted using Equations (17)–(19).

3. 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.

3.1. 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).

$$\left(\begin{array}{c}{\phi }_1\\ {\phi }_2\\ {\omega }_1\\ {\omega }_2\end{array}\right)=\left(\begin{array}{c}\phi (0)\\ \phi (L)\\ \omega (0)\\ \omega (L)\end{array}\right)=\left(\begin{array}{cccc}1& 0& 0& 1\\ \cos {\lambda }_{\mathrm{c}}& \sin {\lambda }_{\mathrm{c}}& L& 1\\ 0& {\lambda }_{\mathrm{c}}/L& 1& 0\\ -{\lambda }_{\mathrm{c}}\sin {\lambda }_{\mathrm{c}}/L& {\lambda }_{\mathrm{c}}\cos {\lambda }_{\mathrm{c}}/L& 1& 0\end{array}\right)\left(\begin{array}{c}{Q}_1\\ Q_2\\ Q_3\\ Q_4\end{array}\right)$$

(20)

The element-end stress resultants (torque and bimoment) are formulated in a matrix form as Equation (21) based on Equations (18) and (19).

$$\left(\begin{array}{c}{T}_1\\ T_2\\ B_{\omega 1}\\ B_{\omega 2}\end{array}\right)=\left(\begin{array}{c}-T\left(0\right)\\ T\left(L\right)\\ B_{\omega }\left(0\right)\\ -B_{\omega }\left(L\right)\end{array}\right)=E{I}_{\omega \omega }\left(\begin{array}{cccc}0& 0& {\left({\lambda }_{\mathrm{c}}/L\right)}^2& 0\\ 0& 0& -{\left({\lambda }_{\mathrm{c}}/L\right)}^2& 0\\ {\left({\lambda }_{\mathrm{c}}/L\right)}^2& 0& 0& 0\\ -{\left({\lambda }_{\mathrm{c}}/L\right)}^2\cos {\lambda }_{\mathrm{c}}& -{\left({\lambda }_{\mathrm{c}}/L\right)}^2\sin {\lambda }_{\mathrm{c}}& 0& 0\end{array}\right)\left(\begin{array}{c}{Q}_1\\ Q_2\\ Q_3\\ Q_4\end{array}\right)$$

(21)

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).

$$\begin{array}{l}\left(\begin{array}{c}{T}_1\\ T_2\\ B_{\omega 1}\\ B_{\omega 2}\end{array}\right)=E{I}_{\omega \omega }\left(\begin{array}{cccc}0& 0& {\left({\lambda }_{\mathrm{c}}/L\right)}^2& 0\\ 0& 0& -{\left({\lambda }_{\mathrm{c}}/L\right)}^2& 0\\ {\left({\lambda }_{\mathrm{c}}/L\right)}^2& 0& 0& 0\\ -{\left({\lambda }_{\mathrm{c}}/L\right)}^2\cos {\lambda }_{\mathrm{c}}& -{\left({\lambda }_{\mathrm{c}}/L\right)}^2\sin {\lambda }_{\mathrm{c}}& 0& 0\end{array}\right)\left(\begin{array}{cccc}\frac{1-\cos {\lambda }_{\mathrm{c}}}{{\varphi }_c}& \frac{\sin {\lambda }_{\mathrm{c}}-{\lambda }_{\mathrm{c}}\cos {\lambda }_{\mathrm{c}}}{k{\varphi }_c}& -\frac{1-\cos {\lambda }_{\mathrm{c}}}{{\varphi }_c}& \frac{{\lambda }_{\mathrm{c}}-\sin {\lambda }_{\mathrm{c}}}{k{\varphi }_c}\\ \frac{k\sin {\lambda }_{\mathrm{c}}}{{\varphi }_c}& \frac{1-\cos {\lambda }_{\mathrm{c}}}{{\varphi }_c}& -\frac{k\sin {\lambda }_{\mathrm{c}}}{{\varphi }_c}& \frac{1-\cos {\lambda }_{\mathrm{c}}}{{\varphi }_c}\\ -\frac{\sin {\lambda }_{\mathrm{c}}}{{\varphi }_c}& \frac{1-\cos {\lambda }_{\mathrm{c}}-{\lambda }_{\mathrm{c}}\sin {\lambda }_{\mathrm{c}}}{k{\varphi }_c}& \frac{\sin {\lambda }_{\mathrm{c}}}{{\varphi }_c}& -\frac{1-\cos {\lambda }_{\mathrm{c}}}{k{\varphi }_c}\\ \frac{1-\cos {\lambda }_{\mathrm{c}}-{\lambda }_{\mathrm{c}}\sin {\lambda }_{\mathrm{c}}}{{\varphi }_c}& -\frac{\sin {\lambda }_{\mathrm{c}}-{\lambda }_{\mathrm{c}}\cos {\lambda }_{\mathrm{c}}}{k{\varphi }_c}& \frac{1-\cos {\lambda }_{\mathrm{c}}}{{\varphi }_c}& -\frac{{\lambda }_{\mathrm{c}}-\sin {\lambda }_{\mathrm{c}}}{k{\varphi }_c}\end{array}\right)\left(\begin{array}{c}{\phi }_1\\ {\phi }_2\\ {\omega }_1\\ {\omega }_2\end{array}\right)\\ =\left[{\mathit{K}}_{\mathbf{e}}\right]\left(\begin{array}{c}{\phi }_1\\ {\phi }_2\\ {\omega }_1\\ {\omega }_2\end{array}\right)\end{array}$$

(22)

where [K_e] is the element stiffness matrix for torsion and warping, which can be simplified as

$$\begin{gathered} \left[{\mathit{K}}_{\mathbf{e}}\right]=\left(\begin{array}{cccc}{T}_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}& -T_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}& Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^2}& Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^2}\\ & T_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}& -Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^2}& -Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^2}\\ & & S_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L}& C_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L}\\ \mathrm{sym}.& & & S_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L}\end{array}\right) \\ \begin{array}{l}{T}_c\left({\lambda }_c\right)={\lambda }_c{}^3\sin {\lambda }_c/{\varphi }_c\left({\lambda }_c\right)\\ Q_c\left({\lambda }_c\right)={\lambda }_c{}^2\left(1-\cos {\lambda }_c\right)/{\varphi }_c\left({\lambda }_c\right)\\ S_c\left({\lambda }_c\right)={\lambda }_c\left(\sin {\lambda }_c-{\lambda }_c\cos {\lambda }_c\right)/{\varphi }_c\left({\lambda }_c\right)\\ C_c\left({\lambda }_c\right)={\lambda }_c\left({\lambda }_c-\sin {\lambda }_c\right)/{\varphi }_c\left({\lambda }_c\right)\\ {\varphi }_c\left({\lambda }_c\right)=2-2\cos {\lambda }_c-{\lambda }_c\sin {\lambda }_c\end{array} \end{gathered}$$

(23)

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.

$$\begin{array}{l}{T}_t\left({\lambda }_t\right)={\lambda }_t{}^3\sinh {\lambda }_t/{\varphi }_t\left({\lambda }_t\right)\\ Q_t\left({\lambda }_t\right)={\lambda }_t{}^2\left(\cosh {\lambda }_t-1\right)/{\varphi }_t\left({\lambda }_t\right)\\ S_t\left({\lambda }_t\right)={\lambda }_t\left({\lambda }_t\cosh {\lambda }_t-\sinh {\lambda }_t\right)/{\varphi }_t\left({\lambda }_t\right)\\ C_t\left({\lambda }_t\right)={\lambda }_t\left(\sinh {\lambda }_t-{\lambda }_t\right)/{\varphi }_t\left({\lambda }_t\right)\\ {\varphi }_t\left({\lambda }_t\right)=2-2\cosh {\lambda }_t+{\lambda }_t\sinh {\lambda }_t\end{array}$$

(24)

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 ${\lambda }_{\mathrm{c}}=\sqrt{\left(P{I}_{\mathrm{p},S}/A-GJ\right)/\left(E{I}_{\omega \omega }/L^2\right)}$, 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).

The linear approximations of the transcendental element stiffness functions T_c, Q_c, S_c, and C_c are given as Equation (25) based on the Taylor series of the functions at λ_c = 0. The approximations are plotted as the dashed lines in 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 −π²EI_ωω/L² and π²EI_ωω/L².

$$\begin{array}{l}{T}_{\mathrm{c}}({\lambda }_{\mathrm{c}})\approx 12-\frac{6}{5}{\lambda }_{\mathrm{c}}{}^2=12+\frac{6}{5}\frac{GJ-P{I}_{\mathrm{p},S}/A}{E{I}_{\omega \omega }/L^2}\\ Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})\approx 6-\frac{1}{10}{\lambda }_{\mathrm{c}}{}^2=6+\frac{1}{10}\frac{GJ-P{I}_{\mathrm{p},S}/A}{E{I}_{\omega \omega }/L^2}\\ S_{\mathrm{c}}({\lambda }_{\mathrm{c}})\approx 4-\frac{2}{15}{\lambda }_{\mathrm{c}}{}^2=4+\frac{2}{15}\frac{GJ-P{I}_{\mathrm{p},S}/A}{E{I}_{\omega \omega }/L^2}\\ C_{\mathrm{c}}({\lambda }_{\mathrm{c}})\approx 2+\frac{1}{30}{\lambda }_{\mathrm{c}}{}^2=2-\frac{1}{30}\frac{GJ-P{I}_{\mathrm{p},S}/A}{E{I}_{\omega \omega }/L^2}\end{array}$$

(25)

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].

$$\left[{\mathit{K}}_{\mathbf{e}}\right]=\frac{E{I}_{\omega \omega }}{L^3}\left(\begin{array}{cccc}12& -12& 6L& 6L\\ & 12& -6L& -6L\\ & & 4L^2& 2L^2\\ \mathrm{sym}.& & & 4L^2\end{array}\right)+\frac{GJ-P{I}_{\mathrm{p},S}/A}{L}\left(\begin{array}{cccc}\frac{6}{5}& -\frac{6}{5}& \frac{L}{10}& \frac{L}{10}\\ & \frac{6}{5}& -\frac{L}{10}& -\frac{L}{10}\\ & & \frac{2L^2}{15}& -\frac{L^2}{30}\\ \mathrm{sym}.& & & \frac{2L^2}{15}\end{array}\right)$$

(26)

3.2. 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 (${\omega }_1={\omega }_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.

$$\left[{\mathit{K}}_{\mathbf{e},\mathbf{tor}}\right]=T_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]$$

(27)

where the torsional stiffness k_tor for this restrained-restrained warping condition can be expressed as

$$k_{\mathrm{tor}}=T_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}$$

(28)

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.

$$\begin{array}{l}\left(\begin{array}{cc}{T}_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]& Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^2}\left[\begin{array}{cc}1& 1\\ -1& -1\end{array}\right]\\ \mathrm{sym}.& \frac{E{I}_{\omega \omega }}{L}\left[\begin{array}{cc}{S}_{\mathrm{c}}({\lambda }_{\mathrm{c}})& C_{\mathrm{c}}({\lambda }_{\mathrm{c}})\\ C_{\mathrm{c}}({\lambda }_{\mathrm{c}})& S_{\mathrm{c}}({\lambda }_{\mathrm{c}})\end{array}\right]\end{array}\right)\left(\begin{array}{c}\left(\begin{array}{l}{\phi }_1\\ {\phi }_2\end{array}\right)\\ \left(\begin{array}{l}{\omega }_1\\ {\omega }_2\end{array}\right)\end{array}\right)=\left(\begin{array}{c}\left(\begin{array}{l}{T}_1\\ T_2\end{array}\right)\\ \left(\begin{array}{l}0\\ 0\end{array}\right)\end{array}\right)\\ \Rightarrow \left\{T_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]-{\left[Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})\right]}^2\frac{E{I}_{\omega \omega }}{L^3}\left[\begin{array}{cc}1& 1\\ -1& -1\end{array}\right]{\left[\begin{array}{cc}{S}_{\mathrm{c}}({\lambda }_{\mathrm{c}})& C_{\mathrm{c}}({\lambda }_{\mathrm{c}})\\ C_{\mathrm{c}}({\lambda }_{\mathrm{c}})& S_{\mathrm{c}}({\lambda }_{\mathrm{c}})\end{array}\right]}^{-1}{\left[\begin{array}{cc}1& 1\\ -1& -1\end{array}\right]}^{\mathrm{T}}\right\}\left(\begin{array}{l}{\phi }_1\\ {\phi }_2\end{array}\right)=\left(\begin{array}{l}{T}_1\\ T_2\end{array}\right)\\ \Rightarrow \left(T_{\mathrm{c}}({\lambda }_{\mathrm{c}})-\frac{2{\left[Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})\right]}^2}{S_{\mathrm{c}}({\lambda }_{\mathrm{c}})+C_{\mathrm{c}}({\lambda }_{\mathrm{c}})}\right)\frac{E{I}_{\omega \omega }}{L^3}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\left(\begin{array}{l}{\phi }_1\\ {\phi }_2\end{array}\right)=\left(\begin{array}{l}{T}_1\\ T_2\end{array}\right)\end{array}$$

(29)

By using the element stiffness functions in Equation (23), the torsional stiffness k_tor for this free-free warping condition is derived as

$$k_{\mathrm{tor}}=\left(T_{\mathrm{c}}({\lambda }_{\mathrm{c}})-\frac{2{\left[Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})\right]}^2}{S_{\mathrm{c}}({\lambda }_{\mathrm{c}})+C_{\mathrm{c}}({\lambda }_{\mathrm{c}})}\right)\frac{E{I}_{\omega \omega }}{L^3}=-{\lambda }_{\mathrm{c}}{}^2\frac{E{I}_{\omega \omega }}{L^3}=\frac{GJ}{L}-\frac{P{I}_{\mathrm{p},S}}{AL}$$

(30)

For members with restrained warping at one end (${\omega }_1=0$) and free warping at the other end (B_ω2 = 0), the torsion stiffness matrix can be derived as follows.

$$\begin{array}{l}\left(\begin{array}{ccc}{T}_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}& -T_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}& Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^2}\\ & T_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}& -Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^2}\\ \mathrm{sym}.& & S_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L}\end{array}\right)\left(\begin{array}{c}{\phi }_1\\ {\phi }_2\\ {\omega }_2\end{array}\right)=\left(\begin{array}{c}{T}_1\\ T_2\\ B_{\omega 2}=0\end{array}\right)\\ \Rightarrow \left\{T_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]-{\left[Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})\right]}^2\frac{E{I}_{\omega \omega }}{L^3}\left[\begin{array}{c}1\\ -1\end{array}\right]{\left[S_{\mathrm{c}}({\lambda }_{\mathrm{c}})\right]}^{-1}{\left[\begin{array}{c}1\\ -1\end{array}\right]}^{\mathrm{T}}\right\}\left(\begin{array}{l}{\phi }_1\\ {\phi }_2\end{array}\right)=\left(\begin{array}{l}{T}_1\\ T_2\end{array}\right)\\ \Rightarrow \left(T_{\mathrm{c}}({\lambda }_{\mathrm{c}})-\frac{{\left[Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})\right]}^2}{S_{\mathrm{c}}({\lambda }_{\mathrm{c}})}\right)\frac{E{I}_{\omega \omega }}{L^3}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\left(\begin{array}{l}{\phi }_1\\ {\phi }_2\end{array}\right)=\left(\begin{array}{l}{T}_1\\ T_2\end{array}\right)\end{array}$$

(31)

By using the element stiffness functions in Equation (23), the torsional stiffness k_tor for this restrained-free warping condition is derived as

$$k_{\mathrm{tor}}=T_{\mathrm{F}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}=\left(T_{\mathrm{c}}({\lambda }_{\mathrm{c}})-\frac{{\left[Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})\right]}^2}{S_{\mathrm{c}}({\lambda }_{\mathrm{c}})}\right)\frac{E{I}_{\omega \omega }}{L^3}=\frac{-{\lambda }_{\mathrm{c}}{}^2}{1-\tan {\lambda }_{\mathrm{c}}/{\lambda }_{\mathrm{c}}}\frac{E{I}_{\omega \omega }}{L^3}$$

(32)

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.

$$T_{\mathrm{F}}({\lambda }_{\mathrm{c}})\approx 3+\frac{6}{5}\frac{GJ-P{I}_{\mathrm{p},S}/A}{E{I}_{\omega \omega }/L^2}$$

(33)

In summary, the torsional stiffnesses of members with three typical element-end warping conditions can be formulated as

$$k_{\mathrm{tor}}=\left\{\begin{array}{cc}{T}_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}\approx \frac{6}{5}\left(\frac{GJ}{L}-\frac{P{I}_{\mathrm{p},S}}{AL}\right)+12\frac{E{I}_{\omega \omega }}{L^3}& restrained-restrained\hspace{0.33em}warping\\ -{\lambda }_{\mathrm{c}}{}^2\frac{E{I}_{\omega \omega }}{L^3}=\frac{GJ}{L}-\frac{P{I}_{\mathrm{p},S}}{AL}& free-free\hspace{0.33em}warping\\ T_{\mathrm{F}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}\approx \frac{6}{5}\left(\frac{GJ}{L}-\frac{P{I}_{\mathrm{p},S}}{AL}\right)+3\frac{E{I}_{\omega \omega }}{L^3}& restrained-free\hspace{0.33em}warping\end{array}\right.$$

(34)

where the linear approximations are based on Section 3.1. Equation (34) shows that 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. 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.

4. 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.

4.1. 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).

$$\det \left[{\mathit{K}}_{\mathbf{s}}({\lambda }_{\mathrm{c},\mathrm{cr}})\right]=0$$

(35)

$$\sqrt{\frac{P_{\mathrm{cr}}{I}_{\mathrm{p},S}/A-GJ}{E{I}_{\omega \omega }/L^2}}={\lambda }_{\mathrm{c},\mathrm{cr}}\Rightarrow P_{\mathrm{cr}}=\left(GJ+{\lambda }_{\mathrm{c},\mathrm{cr}}{}^{2}E{I}_{\omega \omega }/L^2\right)A/I_{\mathrm{p},S}$$

(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. Analysis of torsion-warping-axial stability problems with typical element-end boundary conditions.

Case.	Boundary Condition at Node 1	Boundary Condition at Node 2	Condition for Buckling	Dimensionless Factor at Buckling State (λc,cr)	Buckling Axial Force (Pcr)
1	φ1 = 0	φ2 = 0	Sc−Cc = 0	λc,cr = π	$P_{\mathrm{cr}}=\left(GJ+{\pi }^{2}E{I}_{\omega \omega }/L^2\right)A/I_{\mathrm{p},S}$
2	${\phi }_1=0, {\omega }_1$= 0	/	TcSc−Qc2 = 0	λc,cr = π/2	$P_{\mathrm{cr}}=\left(GJ+{\pi }^{2}E{I}_{\omega \omega }/4L^2\right)A/I_{\mathrm{p},S}$
3	${\phi }_1=0, {\omega }_1$= 0	${\omega }_2$=0	Tc = 0	λc,cr = π	$P_{\mathrm{cr}}=\left(GJ+{\pi }^{2}E{I}_{\omega \omega }/L^2\right)A/I_{\mathrm{p},S}$
4	${\phi }_1=0, {\omega }_1$= 0	φ2 = 0	Sc = 0	λc,cr ≈ π/0.7	$P_{\mathrm{cr}}\approx \left(GJ+2.05{\pi }^{2}E{I}_{\omega \omega }/L^2\right)A/I_{\mathrm{p},S}$
5	${\phi }_1=0, {\omega }_1$= 0	${\phi }_2=0, {\omega }_2$= 0	φc = 0	λc,cr = 2π	$P_{\mathrm{cr}}=\left(GJ+4{\pi }^{2}E{I}_{\omega \omega }/L^2\right)A/I_{\mathrm{p},S}$

.

For a torsion restrained column (φ₁ = 0, φ₂ = 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

$$\det \left[\begin{array}{cc}{S}_{\mathrm{c}}({\lambda }_{\mathrm{c}})& C_{\mathrm{c}}({\lambda }_{\mathrm{c}})\\ C_{\mathrm{c}}({\lambda }_{\mathrm{c}})& S_{\mathrm{c}}({\lambda }_{\mathrm{c}})\end{array}\right]=0\iff S_{\mathrm{c}}({\lambda }_{\mathrm{c}})-C_{\mathrm{c}}({\lambda }_{\mathrm{c}})=0\Rightarrow {\lambda }_{\mathrm{c},\mathrm{cr}}=\pi$$

(37)

For a column with torsion and warping restraints only at one node (φ₁ = 0, ${\omega }_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

$$\det \left[\begin{array}{cc}{T}_{\mathrm{c}}({\lambda }_{\mathrm{c}})& -Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})\\ -Q_{\mathrm{c}}({\lambda }_{\mathrm{c}})& S_{\mathrm{c}}({\lambda }_{\mathrm{c}})\end{array}\right]=0\iff T_{\mathrm{c}}({\lambda }_{\mathrm{c}})S_{\mathrm{c}}({\lambda }_{\mathrm{c}})-Q_{\mathrm{c}}{({\lambda }_{\mathrm{c}})}^2=0\Rightarrow {\lambda }_{\mathrm{c},\mathrm{cr}}=\pi /2$$

(38)

For a column with torsion and warping restraints at one node (φ₁ = 0, ${\omega }_1$ = 0) and a warping restraint at the other node (${\omega }_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

$$T_{\mathrm{c}}({\lambda }_{\mathrm{c}})=0\Rightarrow {\lambda }_{\mathrm{c},\mathrm{cr}}=\pi$$

(39)

For a column with torsion and warping restraints at one node (φ₁ = 0, ${\omega }_1$= 0) and a torsion restraint at the other node (φ₂ = 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

$$S_{\mathrm{c}}({\lambda }_{\mathrm{c}})=0\Rightarrow {\lambda }_{\mathrm{c},\mathrm{cr}}\approx \pi /0.7$$

(40)

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

$${\Phi }_{\mathrm{c}}({\lambda }_{\mathrm{c}})=0\Rightarrow {\lambda }_{\mathrm{c},\mathrm{cr}}=2\pi$$

(41)

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.

4.2. 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 φ₂, ${\omega }_1$, ${\omega }_2$, and ${\omega }_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 (φ₁, φ₂, ${\omega }_1$, and ${\omega }_2$) (or (φ₂, φ₃, ${\omega }_2$, and ${\omega }_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, φ₂, ${\omega }_1$, and ${\omega }_2$) and element II (row/column 1, 3, and 4 of [K_e] associated with the unconstrained element end deformations, φ₂, ${\omega }_2$, and ${\omega }_3$), the structural stiffness matrix [K_s] is formulated as Equation (43).

$$\left[{\mathit{K}}_{\mathbf{s}}\right]\left(\begin{array}{c}{\phi }_2\\ {\omega }_1\\ {\omega }_2\\ {\omega }_3\end{array}\right)=\left(\begin{array}{c}{T}_{\mathrm{mid}}\\ 0\\ 0\\ 0\end{array}\right)$$

(42)

$$\left.\begin{array}{l}\left[{\mathit{K}}_{\mathbf{e},\mathbf{I}}\right]\left\{2,3,4\right\}\to \left[{\mathit{K}}_{\mathbf{s}}\right]\left\{1,2,3\right\}\\ \left[{\mathit{K}}_{\mathbf{e},\mathbf{I}\mathbf{I}}\right]\left\{1,3,4\right\}\to \left[{\mathit{K}}_{\mathbf{s}}\right]\left\{1,3,4\right\}\end{array}\right\}\Rightarrow \left[{\mathit{K}}_{\mathbf{s}}\right]=\left(\begin{array}{cccc}2T_{\mathrm{t}}({\lambda }_{\mathrm{t}})\frac{E{I}_{\omega \omega }}{L^3}& -Q_{\mathrm{t}}({\lambda }_{\mathrm{t}})\frac{E{I}_{\omega \omega }}{L^2}& & Q_{\mathrm{t}}({\lambda }_{\mathrm{t}})\frac{E{I}_{\omega \omega }}{L^2}\\ & S_{\mathrm{t}}({\lambda }_{\mathrm{t}})\frac{E{I}_{\omega \omega }}{L}& C_{\mathrm{t}}({\lambda }_{\mathrm{t}})\frac{E{I}_{\omega \omega }}{L}& \\ & & 2S_{\mathrm{t}}({\lambda }_{\mathrm{t}})\frac{E{I}_{\omega \omega }}{L}& C_{\mathrm{t}}({\lambda }_{\mathrm{t}})\frac{E{I}_{\omega \omega }}{L}\\ \mathrm{sym}.& & & S_{\mathrm{t}}({\lambda }_{\mathrm{t}})\frac{E{I}_{\omega \omega }}{L}\end{array}\right)$$

(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

$${\phi }_2=\frac{\left(1-\tanh {\lambda }_{\mathrm{t}}/{\lambda }_{\mathrm{t}}\right)T_{\mathrm{mid}}}{2{\lambda }_{\mathrm{t}}{}^{2}E{I}_{\omega \omega }/L^3}$$

(44)

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 φ₂, 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).

$${\displaystyle {\int }_0^L{T}_{\mathrm{s}}dx}={\displaystyle {\int }_0^L\left(GJ\frac{d\phi }{dx}\right)dx}=GJ{\phi }_2$$

(45)

$${\displaystyle {\int }_0^L{T}_{\mathrm{w}}dx}={\displaystyle {\int }_0^L\left(T-T_{\mathrm{s}}\right)dx}=B_{\omega ,\mathrm{mid}}$$

(46)

$${\displaystyle {\int }_0^{L}Tdx}={\displaystyle {\int }_0^L{T}_{\mathrm{mid}}/2dx}=T_{\mathrm{mid}}L/2$$

(47)

Then, the ratios of Equation (45) to Equation (47) and Equation (46) to Equation (47) are formulated as Equations (48) and (49), respectively. Equation (48) represents the ratio of the accumulation of St. Venant torque along the element length to the accumulation of total torque, and Equation (49) represents the ratio of the accumulation of warping restraint torque to the accumulation of total torque. The two ratios are plotted in relationships with the factor for St. Venant torsion. As shown in Figure 7, with the increase in the factor for St. Venant torsion (i.e., increase in GJ/EI_ωω), the ratio of the accumulation of warping restraint torque decreases, while the ratio of the accumulation of St. Venant torque increases. Therefore, for thin-walled cross-sections with relatively large St. Venant torsion rigidity GJ (e.g., closed cross-sections), the St. Venant torque could be dominant in the total torque; in contrast, for cross-sections with relatively small St. Venant torsion rigidity GJ (e.g., open cross-sections), the warping restraint torque could be dominant in the total torque.

$$\frac{{\displaystyle {\int }_0^L{T}_{\mathrm{s}}dx}}{{\displaystyle {\int }_0^{L}Tdx}}=\frac{GJ{\phi }_2}{T_{\mathrm{mid}}L/2}=\frac{GJ\frac{\left(1-\tanh {\lambda }_{\mathrm{t}}/{\lambda }_{\mathrm{t}}\right)T_{\mathrm{mid}}}{2{\lambda }_{\mathrm{t}}{}^{2}E{I}_{\omega \omega }/L^3}}{T_{\mathrm{mid}}L/2}=1-\tanh {\lambda }_{\mathrm{t}}/{\lambda }_{\mathrm{t}}$$

(48)

$$\frac{{\displaystyle {\int }_0^L{T}_{\mathrm{w}}dx}}{{\displaystyle {\int }_0^{L}Tdx}}=\frac{B_{\omega ,\mathrm{mid}}}{T_{\mathrm{mid}}L/2}=1-\frac{GJ{\phi }_2}{T_{\mathrm{mid}}L/2}=\tanh {\lambda }_{\mathrm{t}}/{\lambda }_{\mathrm{t}}$$

(49)

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.

4.3. 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.

$$\mathbf{\Delta }={\left(u_{x1},\hspace{0.33em}{u}_{y1},\hspace{0.33em}{u}_{z1},\hspace{0.33em}{\phi }_1,\hspace{0.33em}{\theta }_{y1},\hspace{0.33em}{\theta }_{z1},\hspace{0.33em}{\omega }_1,\hspace{0.33em}\hspace{0.33em}{u}_{x2},\hspace{0.33em}{u}_{y2},\hspace{0.33em}{u}_{z2},\hspace{0.33em}{\phi }_2,\hspace{0.33em}{\theta }_{y2},\hspace{0.33em}{\theta }_{z2},{\omega }_2\right)}^{\mathrm{T}}$$

(50)

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.

$${\mathit{F}}_{\mathrm{ext}}={\left(F_{x1},\hspace{0.33em}{F}_{y1},\hspace{0.33em}{F}_{z1},\hspace{0.33em}{T}_1,\hspace{0.33em}{M}_{y1},\hspace{0.33em}{M}_{z1},\hspace{0.33em}{B}_{\omega 1},\hspace{0.33em}{F}_{x2},\hspace{0.33em}{F}_{y2},\hspace{0.33em}{F}_{z2},\hspace{0.33em}{T}_2,\hspace{0.33em}{M}_{y2},\hspace{0.33em}{M}_{z2},\hspace{0.33em}{B}_{\omega 2}\right)}^{\mathrm{T}}$$

(51)

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:

$$\left[\begin{array}{ccccccccccccc}\frac{EA}{L}& & & & & & & -\frac{EA}{L}& & & & & \\ & T_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^3}& & & & Q_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^2}& & & -T_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^3}& & & & Q_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^2}\\ & & T_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^3}& & -Q_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^2}& & & & & -T_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^3}& & -Q_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^2}& \\ & & & T_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}& & & Q_{\mathrm{c}}\left({\lambda }_{\mathrm{c}}\right)\frac{E{I}_{\omega \omega }}{L^3}& & & & -T_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}& & & Q_{\mathrm{c}}\left({\lambda }_{\mathrm{c}}\right)\frac{E{I}_{\omega \omega }}{L^3}\\ & & & & S_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L}& & & & & Q_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^2}& & C_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L}& \\ & & & & & S_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L}& & & -Q_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^2}& & & & C_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L}\\ & & & & & & S_{\mathrm{c}}\left({\lambda }_{\mathrm{c}}\right)\frac{E{I}_{\omega \omega }}{L^3}& & & & -Q_{\mathrm{c}}\left({\lambda }_{\mathrm{c}}\right)\frac{E{I}_{\omega \omega }}{L^3}& & & C_{\mathrm{c}}\left({\lambda }_{\mathrm{c}}\right)\frac{E{I}_{\omega \omega }}{L^3}\\ & & & & & & & \frac{EA}{L}& & & & & \\ & & & & & & & & T_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^3}& & & & -Q_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^2}\\ & & & & & & & & & T_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^3}& & Q_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^2}& \\ & & & & & & & & & & T_{\mathrm{c}}({\lambda }_{\mathrm{c}})\frac{E{I}_{\omega \omega }}{L^3}& & & -Q_{\mathrm{c}}\left({\lambda }_{\mathrm{c}}\right)\frac{E{I}_{\omega \omega }}{L^3}\\ & & & & & & & & & & & S_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L}& \\ & \mathrm{sym}.& & & & & & & & & & & S_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L}\\ & & & & & & & & & & & & & S_{\mathrm{c}}\left({\lambda }_{\mathrm{c}}\right)\frac{E{I}_{\omega \omega }}{L^3}\end{array}\right]$$

(52)

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.

$$\left(\begin{array}{cccccccccccc}\frac{EA}{L}& & & & & & -\frac{EA}{L}& & & & & \\ & T_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^3}& & & & Q_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^2}& & -T_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^3}& & & & Q_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^2}\\ & & T_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^3}& & -Q_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^2}& & & & -T_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^3}& & -Q_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^2}& \\ & & & k_{tor}& & & & & & -k_{tor}& & \\ & & & & S_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L}& & & & Q_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^2}& & C_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L}& \\ & & & & & S_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L}& & -Q_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^2}& & & & C_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L}\\ & & & & & & \frac{EA}{L}& & & & & \\ & & & & & & & T_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^3}& & & & -Q_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L^2}\\ & & & & & & & & T_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^3}& & Q_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L^2}& \\ & & & & & & & & & k_{tor}& & \\ & \mathrm{sym}.& & & & & & & & & S_{\mathrm{c}}\left({\lambda }_y\right)\frac{E{I}_{yy}}{L}& \\ & & & & & & & & & & & S_{\mathrm{c}}\left({\lambda }_z\right)\frac{E{I}_{zz}}{L}\end{array}\right)\left(\begin{array}{c}{u}_{x1}\\ u_{y1}\\ u_{z1}\\ {\phi }_1\\ {\theta }_{y1}\\ {\theta }_{z1}\\ u_{x2}\\ u_{y2}\\ u_{z2}\\ {\phi }_2\\ {\theta }_{y2}\\ {\theta }_{z2}\end{array}\right)=\left(\begin{array}{c}{F}_{x1}\\ F_{y1}\\ F_{z1}\\ T_1\\ M_{y1}\\ M_{z1}\\ F_{x2}\\ F_{y2}\\ F_{z2}\\ T_2\\ M_{y2}\\ M_{z2}\end{array}\right)$$

(53)

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.

5. 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:

• 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).
• 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.
• 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.