Numerical Model for a Geometrically Nonlinear Analysis of Beams with Composite Cross-Sections

Abstract

This paper presents a beam model for a geometrically nonlinear stability analysis of the composite beam-type structures. Each wall of the cross-section can be modeled with a different material. The nonlinear incremental procedure is based on an updated Lagrangian formulation where in each increment, the equilibrium equations are derived from the virtual work principle. The beam model accounts for the restrained warping and large rotation effects by including the nonlinear displacement field of the composite cross-section. First-order shear deformation theories for torsion and bending are included in the model through Timoshenko’s bending theory and a modified Vlasov’s torsion theory. The shear deformation coupling effects are included in the model using the six shear correction factors. The accuracy and reliability of the proposed numerical model are verified through a comparison of the shear-rigid and shear-deformable beam models in buckling problems. The obtained results indicated the importance of including the shear deformation effects at shorter beams and columns in which the difference that occurs is more than 10 percent.

1. Introduction

Due to their appealing characteristics, composites have occupied a considerable place in the design process of load-bearing structures. Weight, load-bearing capacity, functionality, construction cost, energy efficiency, and resistance to the chemical processes are some of the characteristics now able to be optimized by the use of different composite materials [1,2,3], including the currently widely applied FRP composites [4,5,6,7].

Furthermore, thin-walled structures provide additional possibilities for optimizing structures in some of the above-mentioned elements. Although there is a wider area for the optimization of such structures, the response of such a structure to the effects of external loading is much more complex. It should be emphasized that there is a tendency of such a structure to lose the stability of the deformation form and the appearance of the buckling [8,9]. Therefore, to produce optimal design solutions for such structures, an exact determination of the limit state of the stability of the deformation forms, or of the buckling strength, is of the utmost interest [10,11,12]. We note that the analytical solutions of such complex structures are limited to only simple cases and primarily deal with just one beam/column [13,14,15], and the development of an accurate numerical model is imposed as a necessity [11,12].

Although for slender beam-type structural members, shear deformations do not play a critical role in the determination of the buckling limit state and, respectively, the Euler-Bernoulli and Vlasov theories for bending and torsion could be applied [16,17], a combination of composite materials and various external loading conditions can provide more pronounced shear deformations [18,19,20]. Several geometric nonlinear analyses of composite beam structures with the influence of shear deformations are presented in [10,11,12]. It should be emphasized here that these numerical models include higher numbers of degrees of freedom over the cross-section area, resulting in accurate models capable of dealing with local stability issues, but they are very computationally demanding and, thus, are hardly applicable for frame-type structures. In order to simplify the numerical models and include the shear deformation effects applicable for frame structures, some researchers have applied more conventional beam models [21,22,23,24,25]. For asymmetric thin-walled cross-sections at which the principal bending and shear axes do not coincide, bending–torsion coupling effects arise and are more significant in stability strength [21,22,23,24,25,26].

In [27], composite frames were considered where a geometrically nonlinear beam model and shear deformation effects were presented. Cross-section walls were modeled by cross-ply laminates, and the material was same in each wall of the cross-section. As an extension to that work, an improved beam formulation will be presented herein which takes into account a change in material along the walls of a thin-walled cross-section. The material is assumed to be constant along the thickness of the cross-section wall, but it can be different for each wall. Since the coupling between the normal and shear stresses is not considered, the current beam model is not applicable for angle-ply composites. This will be a topic in our future research activities.

The presented beam model also includes the following assumptions: the cross-section is not deformable in its own plane and the beam member is prismatic and straight, while the external loads are conservative and static, the material obeys Hooke’s law, translations and rotations are allowed to be large but strains are treated as small, and the shear coupling effects are introduced through corresponding shear correction factors [21,22,25,27,28].

The element geometric stiffness is obtained using the nonlinear displacement field of a cross-section. The numerical procedure is formulated in the framework of an updated Lagrangian (UL) formulation [17,29]. To describe large rotations, the nonlinear displacement field includes the second-order displacement terms. To preserve the equilibrium conditions at the frame joint, the internal bending and torsion moments are described as semi-tangential [30,31]. The shear locking [32] does not occur in the beam model due to the interdependent shape functions for the deflections and their slopes and for the twist rotations and the warping parameters. The generalized displacement control method is used in the incremental-iteration algorithm [31]. To update the nodal orientations, special transformations are used due to the semi-tangential description of the internal moments [33,34].

2. Kinematics of the Beam Element

A straight beam element with an arbitrary open thin-walled composite cross-section is presented in Figure 1. The coordinate x- and y-axes are the principal inertial axes of the cross-section, while the longitudinal z-axis passes through the centroid O of each cross-section. The shear center S of the cross-section is defined by the coordinates x_s and y_s. All the cross-section properties are material-weighted.

Displacement increments are described by the rigid-body translations of the cross-sections w_O, u_S, and v_S, and by the rigid-body rotations of the cross-sections φ_z, φ_x, and φ_y. θ is the warping parameter. The displacement field of the cross-section is defined as follows:

$$u=u_0+\tilde{\mathsf{\phi }}\hspace{0.17em}{r}_0-\tilde{\mathsf{\omega }}\hspace{0.17em}s+\frac{1}{2}{\tilde{\mathsf{\phi }}}^2\hspace{0.17em}{r}_0-\frac{1}{2}\hspace{0.17em}\tilde{\mathsf{\phi }}\hspace{0.17em}\tilde{\mathsf{\omega }}\hspace{0.17em}s,$$

(1)

where

$$u_0=\left\{\begin{array}{c}{w}_{\mathrm{O}}+\Omega \hspace{0.17em}\theta \\ u_{\mathrm{S}}\\ v_{\mathrm{S}}\end{array}\right\},\hspace{1em}{r}_0=\left\{\begin{array}{c}0\\ x\\ y\end{array}\right\},\hspace{1em}s=\left\{\begin{array}{c}0\\ x_{\mathrm{S}}\\ y_{\mathrm{S}}\end{array}\right\},\hspace{1em}\tilde{\mathsf{\phi }}=\left[\begin{array}{ccc}0& -{\phi }_{\mathrm{y}}& {\phi }_{\mathrm{x}}\\ {\phi }_{\mathrm{y}}& 0& -{\phi }_{\mathrm{z}}\\ -{\phi }_{\mathrm{x}}& {\phi }_{\mathrm{z}}& 0\end{array}\right],\hspace{1em}\tilde{\mathsf{\omega }}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& -{\phi }_{\mathrm{z}}\\ 0& {\phi }_{\mathrm{z}}& 0\end{array}\right].$$

(2)

As can be seen, the displacement field contains nonlinear terms arising from the large rotation effects [35]. The associated Green–Lagrange strain tensor can be written as:

$$E\approx e+\mathsf{\eta }+\tilde{e}$$

(3)

where

$$\begin{array}{l}e=\frac{1}{2}\left\{\nabla \left(u_0+\tilde{\mathsf{\phi }}\hspace{0.17em}{r}_0-\tilde{\mathsf{\omega }}\hspace{0.17em}s\right)+{\left[\nabla \left(u_0+\tilde{\mathsf{\phi }}\hspace{0.17em}{r}_0-\tilde{\mathsf{\omega }}\hspace{0.17em}s\right)\right]}^{\mathrm{T}}\right\},\\ \mathsf{\eta }=\frac{1}{2}{\left[\nabla \left(u_0+\tilde{\mathsf{\phi }}\hspace{0.17em}{r}_0-\tilde{\mathsf{\omega }}\hspace{0.17em}s\right)\right]}^{\mathrm{T}}\nabla \left(u_0+\tilde{\mathsf{\phi }}\hspace{0.17em}{r}_0-\tilde{\mathsf{\omega }}\hspace{0.17em}s\right),\\ \tilde{e}=\frac{1}{2}\left\{\nabla \left(\frac{1}{2}{\tilde{\mathsf{\phi }}}^2\hspace{0.17em}{r}_0-\frac{1}{2}\hspace{0.17em}\tilde{\mathsf{\phi }}\hspace{0.17em}\tilde{\mathsf{\omega }}\hspace{0.17em}s\right)+{\left[\nabla \left(\frac{1}{2}{\tilde{\mathsf{\phi }}}^2\hspace{0.17em}{r}_0-\frac{1}{2}\hspace{0.17em}\tilde{\mathsf{\phi }}\hspace{0.17em}\tilde{\mathsf{\omega }}\hspace{0.17em}s\right)\right]}^{\mathrm{T}}\right\}\end{array}$$

(4)

It should be noted that the strain increments ε_x, ε_y, and γ_xy in Equation (4) are equal to zero due to the applied beam theories. Since the shear deformation effects are considered, it is valid for the linearized shear strain increments [25]:

$$e_{zy}^{\mathrm{SD}}=\frac{\mathrm{d}{v}_{\mathrm{S}}}{\mathrm{d}z}+{\phi }_x,\hspace{1em}{e}_{zx}^{\mathrm{SD}}=\frac{\mathrm{d}{u}_{\mathrm{S}}}{\mathrm{d}z}-{\phi }_y,\hspace{1em}{\theta }^{\mathrm{SD}}=\frac{\mathrm{d}{\phi }_z}{\mathrm{d}z}+\theta ,$$

(5)

where the right superscript ‘SD’ denotes the quantities arising from the shear deformability.

3. Cross-Section Properties

In Figure 1b, an arbitrary open cross-section composed of a number of thin walls is shown. Each wall is composed of a single orthotropic ply, $\widehat{x}$ and $\widehat{y}$ denote the reference axes, and $\tilde{x}$ and $\tilde{y}$ are the centroid axes while x and y are the principal ones. Furthermore, s and n are the contour and thickness coordinates with their origins in O_ω while R1 is the total number of walls. The incremental constitutive equations valid for a cross-section wall can be written as:

$$\begin{array}{l}\left\{\begin{array}{c}{\sigma }_{zi}\\ {\tau }_{zsi}\end{array}\right\}=\left[\begin{array}{cc}{\overline{Q}}_{11i}& 0\\ 0& {\overline{Q}}_{66i}\end{array}\right]\left\{\begin{array}{c}{e}_z\\ e_{zs}\end{array}\right\},\\ e_{zs}=e_{zs}^{\mathrm{SV}}+e_{zs}^{\mathrm{SD}},\hspace{1em}{e}_{zs}^{\mathrm{SV}}=2\hspace{0.17em}\theta \hspace{0.17em}n,\hspace{1em}{e}_{zs}^{\mathrm{SD}}=e_{zx}^{\mathrm{SD}}\frac{\mathrm{d}x}{\mathrm{d}s}+e_{zy}^{\mathrm{SD}}\frac{\mathrm{d}y}{\mathrm{d}s}+{\theta }^{\mathrm{SD}}\frac{\mathrm{d}\Omega }{\mathrm{d}s}\end{array}$$

(6)

In Equation (6), ${\overline{Q}}_{11i}$ and ${\overline{Q}}_{66i}$ denote the transformed reduced stiffness coefficients [36,37], $e_z$ is the linearized normal strain increment, and $e_{zs}^{\mathrm{SV}}$ is the linearized shear strain increment due the St. Venant torsion and $e_{zs}^{\mathrm{SD}}$ is the linearized shear strain increment occurring at bending with shear and warping torsion, respectively [28,38], while Ω denotes the principal warping function. The reference moduli can be determined as follows:

$$\begin{array}{l}{\overline{Q}}_{11\mathrm{R}}\hspace{0.17em}=\frac{1}{A}{\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}{\overline{Q}}_{11i}\hspace{0.17em}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}=\frac{1}{A}{\displaystyle \sum _{i=1}^{\mathrm{R}1}\left({\overline{Q}}_{11i}\hspace{0.17em}{t}_i\hspace{0.17em}{l}_i\right)}\\ {\overline{Q}}_{66\mathrm{R}}\hspace{0.17em}=\frac{1}{A}{\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}{\overline{Q}}_{66i}\hspace{0.17em}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}=\frac{1}{A}{\displaystyle \sum _{i=1}^{\mathrm{R}1}\left({\overline{Q}}_{66i}\hspace{0.17em}{t}_i\hspace{0.17em}{l}_i\right)}\end{array}$$

(7)

In Figure 2, the local axes (s, n) of each wall should be transformed to the coordinate system of reference axes ($\widehat{x}$, $\widehat{y}$) of the cross-section as:

$$\begin{array}{l}\widehat{x}={\widehat{x}}_{\mathrm{R}}-s\sin {\widehat{\alpha }}_i+n\cos {\widehat{\alpha }}_i\\ \widehat{y}={\widehat{y}}_{\mathrm{R}}+s\cos {\widehat{\alpha }}_i+n\sin {\widehat{\alpha }}_i\end{array}$$

(8)

In the figure, A is an arbitrary point while R denotes the referent point placed on the centerline on the arbitrary position along the wall (i). Applying Equation (8), the modulus-weighted cross-section area and the modulus-weighted first moments of the area can be obtained as follows:

$$A^{\ast }={\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\displaystyle \sum _{i=1}^{\mathrm{R}1}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}{t}_i\hspace{0.17em}{l}_i},$$

(9)

$${\widehat{S}}_x^{\ast }={\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}\widehat{y}\hspace{0.17em}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\displaystyle \sum _{i=1}^{\mathrm{R}1}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}\left({\widehat{y}}_{\mathrm{R}}\hspace{0.17em}{t}_i\hspace{0.17em}{l}_i+\frac{1}{2}{t}_i\hspace{0.17em}{l}_i^2\cos {\widehat{\alpha }}_i\right)},\mathrm{and}$$

(10)

$${\widehat{S}}_y^{\ast }={\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}\widehat{x}\hspace{0.17em}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\displaystyle \sum _{i=1}^{\mathrm{R}1}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}\left({\widehat{x}}_{\mathrm{R}}\hspace{0.17em}{t}_i\hspace{0.17em}{l}_i-\frac{1}{2}\hspace{0.17em}{t}_i\hspace{0.17em}{l}_i^2\sin {\widehat{\alpha }}_i\right)}.$$

(11)

The coordinates of the centroid with reference to the reference axes of the cross-section can be obtained as:

$${\widehat{x}}_{\mathrm{O}}=\frac{{\widehat{S}}_y^{\ast }}{A^{\ast }},\hspace{1em}{\widehat{y}}_{\mathrm{O}}=\frac{{\widehat{S}}_x^{\ast }}{A^{\ast }}.$$

(12)

After the position of the cross-section centroid is determined, the cross-section axes are repositioned to the centroid. All the above cross-sectional properties need to be recalculated for the new position of the axes. Analogously, the modulus-weighted second moments of area can be derived as:

$${\tilde{I}}_x^{\ast }={\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}{\tilde{y}}^2\hspace{0.17em}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\displaystyle \sum _{i=1}^{\mathrm{R}1}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}{t}_i\hspace{0.17em}{l}_i\left[{\left({\tilde{y}}_{\mathrm{R}}+\frac{1}{2}{l}_i\cos {\alpha }_i\right)}^2+\frac{1}{12}{l}_i^2{\cos }^2{\alpha }_i+\frac{1}{12}{t}_i^2{\sin }^2{\alpha }_i\right]},$$

(13)

$${\tilde{I}}_y^{\ast }={\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}{\tilde{x}}^2\hspace{0.17em}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\displaystyle \sum _{i=1}^{\mathrm{R}1}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}{t}_i\hspace{0.17em}{l}_i\left[{\left({\tilde{x}}_{\mathrm{R}}-\frac{1}{2}{l}_i\sin {\alpha }_i\right)}^2+\frac{1}{12}{l}_i^2{\sin }^2{\alpha }_i+\frac{1}{12}{t}_i^2{\cos }^2{\alpha }_i\right]},\mathrm{and}$$

(14)

$${\tilde{I}}_{xy}^{\ast }={\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}\tilde{x}\tilde{y}\hspace{0.17em}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\displaystyle \sum _{i=1}^{\mathrm{R}1}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}{t}_i\hspace{0.17em}{l}_i\left[{\tilde{x}}_{\mathrm{R}}{\tilde{y}}_{\mathrm{R}}+\frac{1}{2}{l}_i\left({\tilde{x}}_{\mathrm{R}}\cos {\alpha }_i-{\tilde{y}}_{\mathrm{R}}\sin {\alpha }_i\right)+\frac{1}{6}\sin 2{\alpha }_i\left(\frac{1}{4}{t}_i^2-l_i^2\right)\right]}.$$

(15)

After that, the angles of the principal axes can be obtained from the well-known formula:

$$\tan 2\chi =\frac{2{\tilde{I}}_{xy}^{\ast }}{{\tilde{I}}_y^{\ast }-{\tilde{I}}_x^{\ast }}.$$

(16)

Establishing the angle from Equation (12), the cross-section properties should be recalculated for the principal axes of the cross-section.

The warping function is defined for the pole set at the centroid of the cross-section and the associated origin is set according to Figure 1b. From Figure 2b, the difference in contour warping function can be calculated as:

$$\Delta \hspace{0.17em}\overline{\omega }={\overline{h}}_{\mathrm{R}}{l}_i,\hspace{1em}{\overline{h}}_{\mathrm{R}}=x_{\mathrm{R}}\cos {\alpha }_i+y_{\mathrm{R}}\sin {\alpha }_i,$$

(17)

while, from Figure 2c, the difference in thickness warping function is:

$$\Delta \hspace{0.17em}\overline{\overline{\omega }}=-\left({\overline{\overline{h}}}_{\hspace{0.17em}\mathrm{R}}+l_i\right)n,\hspace{1em}{\overline{\overline{h}}}_{\mathrm{R}}=y_{\mathrm{R}}\cos {\alpha }_i-x_{\mathrm{R}}\sin {\alpha }_i.$$

(18)

By combining Equations (13) and (14), the expression for the warping function of an arbitrary point A is derived as:

$$\omega ={\omega }_{\hspace{0.17em}\mathrm{R}}+\Delta \hspace{0.17em}\overline{\omega }+\Delta \hspace{0.17em}\overline{\overline{\omega }}={\omega }_{\hspace{0.17em}\mathrm{R}}+{\overline{h}}_{\hspace{0.17em}\mathrm{R}}s-{\overline{\overline{h}}}_{\mathrm{R}}n-s\hspace{0.17em}n.$$

(19)

The warping function from Equation (15) is used to derive, respectively, the modulus-weighted first and second sectorial moments as:

$$S_{\omega }^{\ast }={\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}\omega \hspace{0.17em}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\displaystyle \sum _{i=1}^{\mathrm{R}1}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}{t}_i\hspace{0.17em}{l}_i\left({\omega }_{\mathrm{R}}+\frac{1}{2}{l}_i\hspace{0.17em}{\overline{h}}_{\mathrm{R}}\right)},$$

(20)

$$\begin{array}{ll}{I}_{x\omega }^{\ast }& ={\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}x\omega \hspace{0.17em}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\displaystyle \sum _{i=1}^{\mathrm{R}1}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}{t}_i\hspace{0.17em}{l}_i[\left({\omega }_{\mathrm{R}}+\frac{1}{2}{\overline{h}}_{\mathrm{R}}\hspace{0.17em}{l}_i\right)\left(x_{\mathrm{R}}-\frac{1}{2}{l}_i\sin {\alpha }_i\right)}\\ & -\frac{1}{12}{\overline{h}}_{\mathrm{R}}\hspace{0.17em}{l}_i^2\sin {\alpha }_i-\frac{1}{12}{t}_i^2\cos {\alpha }_i\left({\overline{\overline{h}}}_{\mathrm{R}}+\frac{1}{2}{l}_i\right)]\end{array}, \mathrm{and}$$

(21)

$$\begin{array}{ll}{I}_{y\omega }^{\ast }& ={\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}y\omega \hspace{0.17em}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\displaystyle \sum _{i=1}^{\mathrm{R}1}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\hspace{0.17em}\mathrm{R}}}{t}_i\hspace{0.17em}{l}_i[\left({\omega }_{\mathrm{R}}+\frac{1}{2}{\overline{h}}_{\mathrm{R}}\hspace{0.17em}{l}_i\right)\left(y_{\mathrm{R}}+\frac{1}{2}{l}_i\cos {\alpha }_i\right)}\\ & +\frac{1}{12}{\overline{h}}_{\mathrm{R}}\hspace{0.17em}{l}_i^2\cos {\alpha }_i-\frac{1}{12}{t}_i^2\sin {\alpha }_i\left({\overline{\overline{h}}}_{\mathrm{R}}+\frac{1}{2}{l}_i\right)]\end{array},$$

(22)

and afterwards, the shear center position coordinates with reference to the principal axes can be calculated as:

$$x_{\mathrm{S}}=\frac{I_{y\omega }^{\ast }}{I_x^{\ast }},\hspace{1em}{y}_{\mathrm{S}}=-\frac{I_{x\omega }^{\ast }}{I_y^{\ast }}.$$

(23)

Knowing the position of the shear center, the principal warping function is derived by the following transformation [26]:

$$\begin{array}{l}\Omega =\omega +y_{\mathrm{S}}\hspace{0.17em}x-x_S\hspace{0.17em}y-S_{\omega }^{\ast }/A^{\ast }={\Omega }_{\hspace{0.17em}\mathrm{R}}+r\hspace{0.17em}s-q\hspace{0.17em}n-s\hspace{0.17em}n,\\ {\Omega }_{\hspace{0.17em}\mathrm{R}}={\omega }_{\hspace{0.17em}\mathrm{R}}+y_{\mathrm{S}}\hspace{0.17em}{x}_{\mathrm{R}}-x_{\mathrm{S}}\hspace{0.17em}{y}_{\mathrm{R}}-S_{\omega }^{\ast }/A^{\ast }\\ r=\left(x_{\mathrm{R}}-x_{\mathrm{S}}\right)\cos {\alpha }_i+\left(y_{\mathrm{R}}-y_{\mathrm{S}}\right)\sin {\alpha }_i,\hspace{1em}\\ q=\left(y_{\mathrm{R}}-y_{\mathrm{S}}\right)\cos {\alpha }_i-\left(x_{\mathrm{R}}-x_{\mathrm{S}}\right)\sin {\alpha }_i,\end{array}$$

(24)

and the modulus-weighted warping constant can be determined as:

$$I_{\omega }^{\ast }={\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}{\Omega }^2\hspace{0.17em}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\displaystyle \sum _{i=1}^{\mathrm{R}1}\frac{{\overline{Q}}_{11i}}{{\overline{Q}}_{11\mathrm{R}}}{t}_i\hspace{0.17em}{l}_i\left[{\left({\Omega }_{\hspace{0.17em}\mathrm{R}}+\frac{1}{2}r\hspace{0.17em}{l}_i\right)}^2+\frac{1}{12}{l}_i^2\left(r^2+\frac{1}{3}{t}_i^2\right)+\frac{1}{12}q\hspace{0.17em}{t}_i^2\left(q+l_i\right)\right]}.$$

(25)

4. Governing Equations

Applying the conventional engineering theories for torsion and bending, the stress resultant increments, shown in Figure 1a, can be calculated as follows:

$$\begin{array}{l}{F}_z={\displaystyle \underset{A}{\int }{\sigma }_z\hspace{0.17em}\mathrm{d}A}\hspace{0.17em},\hspace{1em}{F}_x={\displaystyle \underset{A}{\int }{\tau }_{zx}\hspace{0.17em}\mathrm{d}A}\hspace{0.17em},\hspace{1em}{F}_y={\displaystyle \underset{A}{\int }{\tau }_{zy}\hspace{0.17em}\mathrm{d}A},\hspace{1em}{M}_x={\displaystyle \underset{A}{\int }{\sigma }_z\hspace{0.17em}y\hspace{0.17em}\mathrm{d}A},\\ M_y=-{\displaystyle \underset{A}{\int }{\sigma }_z\hspace{0.17em}x\hspace{0.17em}\mathrm{d}A},\hspace{1em}{M}_{\omega }={\displaystyle \underset{A}{\int }{\sigma }_z\hspace{0.17em}\omega \hspace{0.17em}\mathrm{d}A},\hspace{1em}{M}_z=T_{\mathrm{SV}}+T_{\omega },\\ \hspace{0.17em}{T}_{\omega }={\displaystyle \underset{A}{\int }\left({\tau }_{zx}\frac{\partial \omega }{\partial y}+{\tau }_{zy}\frac{\partial \omega }{\partial x}\right)\hspace{0.17em}\mathrm{d}A},\hspace{1em}\overline{K}={\displaystyle \underset{A}{\int }{\sigma }_z\hspace{0.17em}\left[{\left(x-x_{\mathrm{S}}\right)}^2+{\left(y-y_{\mathrm{S}}\right)}^2\right]\hspace{0.17em}\mathrm{d}A},\\ T_{\mathrm{SV}}={\displaystyle \underset{A}{\int }\left[{\tau }_{zy}\left(\hspace{0.17em}x-x_{\mathrm{S}}-\frac{\partial \omega }{\partial y}\right)+{\tau }_{zx}\left(\hspace{0.17em}y-y_{\mathrm{S}}+\frac{\partial \omega }{\partial x}\right)\right]\hspace{0.17em}\mathrm{d}A}\end{array}$$

(26)

where F_z is the axial force, F_x and F_y are the shear forces, T_SV is the St. Venant torsion moment, T_ω is the warping torsion moment, M_x and M_y are the bending moments, M_ω is the bimoment, and $\overline{K}$ is the Wagner coefficient [17,38]. For a thin-walled composite cross-section, the St. Venant torsion moment can be calculates as [38]:

$$T_{\mathrm{SV}}=2\theta {\displaystyle \int {\displaystyle \underset{-t/2}{\overset{t/2}{\int }}{\overline{Q}}_{66i}\left(\frac{t_i^2}{4}-n^2\right)\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\overline{Q}}_{66\mathrm{R}}\hspace{0.17em}\theta I_t^{\ast },$$

(27)

where $I_t^{\ast }$ is the modulus-weighted torsional moment of inertia, and it is defined as:

$$I_t^{\ast }=2{\displaystyle \int {\displaystyle \underset{-t/2}{\overset{t/2}{\int }}\frac{{\overline{Q}}_{66i}}{{\overline{Q}}_{66\mathrm{R}}}\left(\frac{t_i^2}{4}-n^2\right)\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\displaystyle \sum _{i=1}^{\mathrm{R}1}\frac{{\overline{Q}}_{66i}}{{\overline{Q}}_{66\mathrm{R}}}\frac{1}{3}{l}_i\hspace{0.17em}{t}_i^3}.$$

(28)

From Equations (4), (6), and (20), expressions for the axial force, the bending moments, and the bimoment for a thin-walled composite cross-section can be derived as follows:

$$F_z={\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}{\sigma }_z\hspace{0.17em}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\overline{Q}}_{11\mathrm{R}}\left(A^{\ast }\frac{\mathrm{d}{w}_{\mathrm{O}}}{\mathrm{d}z}+S_x^{\ast }\frac{\mathrm{d}{\phi }_x}{\mathrm{d}z}-S_y^{\ast }\frac{\mathrm{d}{\phi }_y}{\mathrm{d}z}+S_{\omega }^{\ast }\frac{\mathrm{d}\theta }{\mathrm{d}z}\right),$$

(29)

$$M_x={\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}{\sigma }_{z}y\hspace{0.17em}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\overline{Q}}_{11\mathrm{R}}\left(S_x^{\ast }\frac{\mathrm{d}{w}_{\mathrm{O}}}{\mathrm{d}z}+I_x^{\ast }\frac{\mathrm{d}{\phi }_x}{\mathrm{d}z}-I_{xy}^{\ast }\frac{\mathrm{d}{\phi }_y}{\mathrm{d}z}+I_{y\omega }^{\ast }\frac{\mathrm{d}\theta }{\mathrm{d}z}\right),$$

(30)

$$M_y=-{\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}{\sigma }_{z}x\hspace{0.17em}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\overline{Q}}_{11\mathrm{R}}\left(-S_y^{\ast }\frac{\mathrm{d}{w}_{\mathrm{O}}}{\mathrm{d}z}-I_{xy}^{\ast }\frac{\mathrm{d}{\phi }_x}{\mathrm{d}z}+I_y^{\ast }\frac{\mathrm{d}{\phi }_y}{\mathrm{d}z}-I_{x\omega }^{\ast }\frac{\mathrm{d}\theta }{\mathrm{d}z}\right),\mathrm{and}$$

(31)

$$M_{\omega }={\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}{\sigma }_z\omega \hspace{0.17em}\hspace{0.17em}\mathrm{d}n\hspace{0.17em}\mathrm{d}s}}={\overline{Q}}_{11\mathrm{R}}\left(S_{\omega }^{\ast }\frac{\mathrm{d}{w}_{\mathrm{O}}}{\mathrm{d}z}+I_{y\omega }^{\ast }\frac{\mathrm{d}{\phi }_x}{\mathrm{d}z}-I_{x\omega }^{\ast }\frac{\mathrm{d}{\phi }_y}{\mathrm{d}z}+I_{\omega }^{\ast }\frac{\mathrm{d}\theta }{\mathrm{d}z}\right),$$

(32)

where $\mathrm{d}{w}_{\mathrm{O}}/\mathrm{d}z$, $\mathrm{d}{\phi }_x/\mathrm{d}z$, $\mathrm{d}{\phi }_y/\mathrm{d}z$, and $\mathrm{d}\theta /\mathrm{d}z$ are the generalized deformation increments.

Since the following is valid tor the principal cross-sectional axes and the principal warping function:

$$S_x^{\ast }=S_y^{\ast }=I_{xy}^{\ast }=S_{\omega }^{\ast }=I_{x\omega }^{\ast }=I_{y\omega }^{\ast }=0,$$

(33)

the generalized deformation increments from Equations (29)–(32) can be rewritten as:

$$\frac{\mathrm{d}{w}_{\mathrm{O}}}{\mathrm{d}z}=\frac{F_z}{{\overline{Q}}_{11\mathrm{R}}{A}^{\ast }},\hspace{1em}\frac{\mathrm{d}{\phi }_x}{\mathrm{d}z}=\frac{M_x}{{\overline{Q}}_{11\mathrm{R}}{I}_x^{\ast }},\hspace{1em}\frac{\mathrm{d}{\phi }_y}{\mathrm{d}z}=\frac{M_y}{{\overline{Q}}_{11\mathrm{R}}{I}_y^{\ast }},\hspace{1em}\frac{\mathrm{d}\theta }{\mathrm{d}z}=\frac{M_{\omega }}{{\overline{Q}}_{11\mathrm{R}}{I}_{\omega }^{\ast }},$$

(34)

while, for the normal stress increment, one can obtain:

$${\sigma }_{z\mathrm{R}}={\overline{Q}}_{11\mathrm{R}}\hspace{0.17em}{e}_z=\frac{F_z}{A^{\ast }}+\frac{M_x}{I_x^{\ast }}y-\frac{M_y}{I_y^{\ast }}x+\frac{M_{\omega }}{I_{\omega }^{\ast }}\omega .$$

(35)

Since the shear stress increment is uniformly distributed over a wall thickness, it can be expressed as:

$${\tau }_{zs\mathrm{R}}^{\mathrm{SD}}={\overline{Q}}_{66\mathrm{R}}\hspace{0.17em}{e}_{zs}^{\mathrm{SD}}=\frac{F_y\hspace{0.17em}{S}_{x1}^{\ast }}{I_x^{\ast }{t}_i}+\frac{F_x\hspace{0.17em}{S}_{y1}^{\ast }}{I_y^{\ast }{t}_i}+\frac{T_{\omega }{S}_{\omega 1}^{\ast }}{I_{\omega }^{\ast }{t}_i},$$

(36)

where $S_{x1}^{\ast }$, $S_{y1}^{\ast }$, and $S_{\omega 1}^{\ast }$ are, respectively, the modulus-weighted first moments of the area and the first sectorial moment for the cross-section cut off.

The virtual incremental strain energy due to shear deformability is derived using Equations (6) and (27) as follows:

$$\mathsf{\delta }{\mathcal{U}}_{\hspace{0.17em}\mathrm{E}}^{\mathrm{SD}}={\displaystyle \underset{0}{\overset{l}{\int }}{\displaystyle \underset{s}{\int }{\displaystyle \underset{-t/2}{\overset{t/2}{\int }}{\tau }_{zsi}^{\mathrm{SD}}\hspace{0.17em}\mathsf{\delta }{e}_{zs}^{\mathrm{SD}}\hspace{0.17em}\mathrm{d}n\hspace{0.17em}\mathrm{d}s\hspace{0.17em}\mathrm{d}z}}}={\displaystyle \underset{0}{\overset{l}{\int }}{\displaystyle \underset{s}{\int }{\overline{Q}}_{66\mathrm{R}}\frac{{\overline{Q}}_{66i}}{{\overline{Q}}_{66\mathrm{R}}}{t}_i\hspace{0.17em}{e}_{zs}^{\mathrm{SD}}\hspace{0.17em}\mathsf{\delta }{e}_{zs}^{\mathrm{SD}}\mathrm{d}s\hspace{0.17em}\mathrm{d}z}}={\displaystyle \underset{0}{\overset{l}{\int }}{\displaystyle \underset{s}{\int }{t}_i\hspace{0.17em}{\tau }_{zs\mathrm{R}}^{\mathrm{SD}}\frac{\mathsf{\delta }{\tau }_{zs\hspace{0.17em}\mathrm{R}}^{\mathrm{SD}}}{{\overline{Q}}_{66\mathrm{R}}}\frac{{\overline{Q}}_{66i}}{{\overline{Q}}_{66\mathrm{R}}}\hspace{0.17em}\mathrm{d}s\hspace{0.17em}\mathrm{d}z}},$$

(37)

i.e.,

$$\begin{array}{ll}\mathsf{\delta }{\mathcal{U}}_{\hspace{0.17em}\mathrm{E}}^{\mathrm{SD}}& =\frac{1}{{\overline{Q}}_{66\mathrm{R}}}{\displaystyle \underset{0}{\overset{l}{\int }}\{\frac{k_x}{A^{\ast }}{F}_x\hspace{0.17em}\mathsf{\delta }{F}_x+\frac{k_y}{A^{\ast }}{F}_y\hspace{0.17em}\mathsf{\delta }{F}_y+\frac{k_{\omega }}{I_t^{\ast }}{T}_{\omega }\hspace{0.17em}\mathsf{\delta }{T}_{\omega }+\frac{k_{xy}}{A^{\ast }}\left(F_x\hspace{0.17em}\mathsf{\delta }{F}_y+F_y\hspace{0.17em}\mathsf{\delta }{F}_x\right)}\\ & +\frac{k_{x\omega }}{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}\left(F_x\hspace{0.17em}\mathsf{\delta }{T}_{\omega }+T_{\omega }\hspace{0.17em}\mathsf{\delta }{F}_x\right)+\frac{k_{y\omega }}{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}\left(F_y\hspace{0.17em}\mathsf{\delta }{T}_{\omega }+T_{\omega }\hspace{0.17em}\mathsf{\delta }{F}_y\right)\}\mathrm{d}z\end{array}.$$

(38)

In Equation (29), $k_x$, $k_y$, $k_{\omega }$, $k_{xy}$, $k_{x\omega }$, and $k_{y\omega }$ denote the shear correction factors calculated as:

$$\begin{array}{l}\frac{k_x}{A^{\ast }}=\frac{1}{{\left(I_y^{\ast }\right)}^2}{\displaystyle \underset{s}{\int }\frac{{\left(S_{y1}^{\ast }\right)}^2}{t_i}\mathrm{d}s}\hspace{0.17em},\hspace{1em}\frac{k_y}{A^{\ast }}=\frac{1}{{\left(I_x^{\ast }\right)}^2}{\displaystyle \underset{s}{\int }\frac{{\left(S_{x1}^{\ast }\right)}^2}{t_i}\mathrm{d}s},\hspace{0.17em}\hspace{1em}\frac{k_{\omega }}{I_t^{\ast }}=\frac{1}{{\left(I_{\omega }^{\ast }\right)}^2}{\displaystyle \underset{s}{\int }\frac{{\left(S_{\omega 1}^{\ast }\right)}^2}{t_i}\mathrm{d}s}\hspace{0.17em}\\ \frac{k_{xy}}{A^{\ast }}=\frac{1}{I_x^{\ast }\hspace{0.17em}{I}_y^{\ast }}{\displaystyle \underset{s}{\int }\frac{S_{x1}^{\ast }\hspace{0.17em}{S}_{y1}^{\ast }}{t_i}\mathrm{d}s}\hspace{0.17em},\hspace{1em}\frac{k_{x\omega }}{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}=\frac{1}{I_y^{\ast }\hspace{0.17em}{I}_{\omega }^{\ast }}{\displaystyle \underset{s}{\int }\frac{S_{y1}^{\ast }\hspace{0.17em}{S}_{\omega 1}^{\ast }}{t_i}\mathrm{d}s}\hspace{0.17em},\hspace{1em}\frac{k_{y\omega }}{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}=\frac{1}{I_x^{\ast }\hspace{0.17em}{I}_{\omega }^{\ast }}{\displaystyle \underset{s}{\int }\frac{S_{x1}^{\ast }\hspace{0.17em}{S}_{\omega 1}^{\ast }}{t_i}\mathrm{d}s}\end{array}$$

(39)

By using Equations (5) and (6) in Equation (20), the shear forces and the warping torsion moment increments can be expressed in terms of the linearized shear strain increments:

$$\left\{\begin{array}{c}{F}_{\mathrm{x}}\\ F_{\mathrm{y}}\\ T_{\mathsf{\omega }}\end{array}\right\}={\overline{Q}}_{66\mathrm{R}}\left[\begin{array}{ccc}\frac{A^{\ast }}{K_{\mathrm{x}}}& \frac{A^{\ast }}{K_{\mathrm{xy}}}& \frac{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}{K_{\mathrm{x}\omega }}\\ \frac{A^{\ast }}{K_{\mathrm{xy}}}& \frac{A^{\ast }}{K_{\mathrm{y}}}& \frac{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}{K_{\mathrm{y}\omega }}\\ \frac{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}{K_{\mathrm{x}\omega }}& \frac{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}{K_{\mathrm{y}\omega }}& \frac{I_t^{\ast }}{K_{\omega }}\end{array}\right]\left\{\begin{array}{c}{e}_{zx}^{\mathrm{SD}}\\ e_{zy}^{\mathrm{SD}}\\ {\theta }^{\mathrm{SD}}\end{array}\right\},$$

(40)

where

$$\left[\begin{array}{ccc}\frac{A^{\ast }}{K_x}& \frac{A^{\ast }}{K_{xy}}& \frac{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}{K_{x\omega }}\\ \frac{A^{\ast }}{K_{xy}}& \frac{A^{\ast }}{K_y}& \frac{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}{K_{y\omega }}\\ \frac{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}{K_{x\omega }}& \frac{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}{K_{y\omega }}& \frac{I_t^{\ast }}{K_{\omega }}\end{array}\right]={\left[\begin{array}{ccc}\frac{k_x}{A^{\ast }}& \frac{k_{xy}}{A^{\ast }}& \frac{k_{x\omega }}{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}\\ \frac{k_{xy}}{A^{\ast }}& \frac{k_y}{A^{\ast }}& \frac{k_{y\omega }}{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}\\ \frac{k_{x\omega }}{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}& \frac{k_{y\omega }}{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}& \frac{k_{\omega }}{I_t^{\ast }}\end{array}\right]}^{-1}.$$

(41)

In Equation (32), k_ij and K_ij are dimensionless constants, while the explicit solutions of the integrals from Equation (30) are provided in Appendix A.

5. Finite Element Formulation

A thin-walled beam element, as shown in Figure 1, consists of two nodes, and each node has seven degrees-of-freedom. The corresponding nodal displacement and force vectors are:

$$u^e=\left\{\begin{array}{c}{u}_{\mathrm{A}}^e\\ u_{\mathrm{B}}^e\end{array}\right\};\hspace{1em}{\left(u_i^e\right)}^{\mathrm{T}}=\left\{\begin{array}{ccccccc}{w}_i& u_i& v_i& {\phi }_{z\hspace{0.17em}i}& {\phi }_{x\hspace{0.17em}i}& {\phi }_{y\hspace{0.17em}i}& {\theta }_i\end{array}\right\},\hspace{1em}i=\mathrm{A}, \mathrm{B}\mathrm{and}$$

(42)

$$f^e=\left\{\begin{array}{c}{f}_{\mathrm{A}}^e\\ f_{\mathrm{B}}^e\end{array}\right\};\hspace{1em}{\left(f_i^e\right)}^{\mathrm{T}}=\left\{\begin{array}{ccccccc}{F}_{z\hspace{0.17em}i}& F_{x\hspace{0.17em}i}& F_{y\hspace{0.17em}i}& M_{z\hspace{0.17em}i}& M_{x\hspace{0.17em}i}& M_{y\hspace{0.17em}i}& M_{\omega \hspace{0.17em}i}\end{array}\right\},\hspace{1em}i=\mathrm{A}, \mathrm{B},$$

(43)

where the right superscript e denotes the e-th finite element. The virtual work principle for a beam element deforming between the last calculated configuration C1 and the current unknown configuration C2 can be expressed as follows:

$$\mathsf{\delta }{\mathcal{U}}_{\hspace{0.17em}\mathrm{E}}+\mathsf{\delta }{\mathcal{U}}_{\hspace{0.17em}\mathrm{G}}={\mathsf{\delta }}^2\mathcal{W}-{\mathsf{\delta }}^1\mathcal{W},$$

(44)

where, respectively, the virtual incremental elastic strain and geometric energy are determined as:

$$\mathsf{\delta }{\mathcal{U}}_{\mathrm{E}}\approx {\displaystyle \underset{1_V}{\int }{\mathsf{\delta }}_1{e}^{\mathrm{T}}:\hspace{0.17em}{}_{1}C:{}_{1}e^1\mathrm{d}V} \mathrm{and}$$

(45)

$$\mathsf{\delta }{\mathcal{U}}_{\mathrm{G}}={\displaystyle \underset{1_V}{\int }{}_1{}^{1}S:{\mathsf{\delta }}_1{\mathsf{\eta }}^1\hspace{0.17em}\mathrm{d}V}+{\displaystyle \underset{1_V}{\int }{}_1{}^{1}S:{\mathsf{\delta }}_1{\tilde{e}}^1\hspace{0.17em}\mathrm{d}V}-{\displaystyle \underset{A_{\mathrm{t}}}{\int }{{}_1{}^{1}t}_{\mathrm{n}}\cdot {\mathsf{\delta }}_1{\tilde{\mathbf{u}}}^1\mathrm{d}{A}_{\mathrm{t}}}.$$

(46)

In Equation (35), the quantities on the right-hand side represent the virtual incremental work performed by the nodal forces, i.e.,

$${\mathsf{\delta }}^2\mathcal{W}={\displaystyle \underset{1_{A_{\mathrm{t}}}}{\int }{}_1{}^{2}t{}_{\mathrm{n}}\cdot {\mathsf{\delta }}_1{{\mathbf{u}}_{\mathrm{ldf}}}^1\mathrm{d}{A}_{\hspace{0.17em}\mathrm{t}}}={\left({\mathsf{\delta }}_1{u}^e\right)}^{\mathrm{T}}{}_1{}^{2}f{}^e\mathrm{and}$$

(47)

$${\mathsf{\delta }}^1\mathcal{W}={\displaystyle \underset{1_{A_{\mathrm{t}}}}{\int }{}_1{}^{1}t{}_{\mathrm{n}}\cdot {\mathsf{\delta }}_1{{\mathbf{u}}_{\mathrm{ldf}}}^1\mathrm{d}{A}_{\hspace{0.17em}\mathrm{t}}}={\left({\mathsf{\delta }}_1{u}^e\right)}^{\mathrm{T}}{}_1{}^{1}f{}^e.$$

(48)

In the above equations, S is the second Piola–Kirchhoff stress tensor, t_n represents the surface tractions, and C is the incremental constitutive tensor. According to the applied UL incremental description, all the quantities occurring in Equations (36)–(39) are expressed with reference to the configuration C1. The index notation adopted in the presented work is the same as that in [31], i.e., the left superscript refers to the configuration in which a quantity occurs while the left subscript refers to the configuration in which the quantity is measured.

Furthermore, the virtual incremental elastic strain energy can be decomposed as:

$$\mathsf{\delta }{\mathcal{U}}_{\hspace{0.17em}\mathrm{E}}=\mathsf{\delta }{\mathcal{U}}_{\hspace{0.17em}\mathrm{E}}^{\mathrm{EB}}+\mathsf{\delta }{\mathcal{U}}_{\hspace{0.17em}\mathrm{E}}^{\mathrm{SD}},$$

(49)

where the first term on the right-hand side denotes the incremental energy obtained by the Euler–Bernoulli and St. Venant beam theories for bending and torsion, respectively, i.e.,

$$\mathsf{\delta }{\mathcal{U}}_{\hspace{0.17em}\mathrm{E}}^{\mathrm{EB}}={\displaystyle \underset{0}{\overset{l}{\int }}{\overline{Q}}_{11\mathrm{R}}\left(A^{\ast }\frac{\mathrm{d}{w}_O}{\mathrm{d}z}\mathsf{\delta }\frac{\mathrm{d}{w}_O}{\mathrm{d}z}+I_x^{\ast }\frac{\mathrm{d}{\phi }_x}{\mathrm{d}z}\mathsf{\delta }\frac{\mathrm{d}{\phi }_x}{\mathrm{d}z}+I_y^{\ast }\frac{\mathrm{d}{\phi }_y}{\mathrm{d}z}\mathsf{\delta }\frac{\mathrm{d}{\phi }_y}{\mathrm{d}z}+I_{\omega }^{\ast }\frac{\mathrm{d}\theta }{\mathrm{d}z}\mathsf{\delta }\frac{\mathrm{d}\theta }{\mathrm{d}z}+I_t^{\ast }\theta \hspace{0.17em}\mathsf{\delta }\theta \right)\mathrm{d}z},$$

(50)

while the second energy term is due to the shear deformability which, according to Equations (29) and (31), can be defined as:

$$\begin{array}{ll}\mathsf{\delta }{\mathcal{U}}_{\hspace{0.17em}\mathrm{E}}^{\mathrm{SD}}& ={\displaystyle \underset{0}{\overset{l}{\int }}{\overline{Q}}_{66\mathrm{R}}[\frac{A^{\ast }}{K_x}{e}_{zx}^{\mathrm{SD}}\hspace{0.17em}\mathsf{\delta }{e}_{zx}^{\mathrm{SD}}+\frac{A^{\ast }}{K_y}{e}_{zy}^{\mathrm{SD}}\hspace{0.17em}\mathsf{\delta }{e}_{zy}^{\mathrm{SD}}+\frac{I_t^{\ast }}{K_{\omega }}{\theta }^{\mathrm{SD}}\hspace{0.17em}\mathsf{\delta }{\theta }^{\mathrm{SD}}+\frac{A^{\ast }}{K_{xy}}\left(e_{zx}^{\mathrm{SD}}\hspace{0.17em}\mathsf{\delta }{e}_{zy}^{\mathrm{SD}}+e_{zy}^{\mathrm{SD}}\hspace{0.17em}\mathsf{\delta }{e}_{zx}^{\mathrm{SD}}\right)+}\\ & +\frac{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}{K_{x\omega }}\left(e_{zx}^{\mathrm{SD}}\hspace{0.17em}\mathsf{\delta }{\theta }^{\mathrm{SD}}+{\theta }^{\mathrm{SD}}\hspace{0.17em}\mathsf{\delta }{e}_{zx}^{\mathrm{SD}}\right)+\frac{\sqrt{A^{\ast }\hspace{0.17em}{I}_t^{\ast }}}{K_{y\omega }}\left(e_{zy}^{\mathrm{SD}}\hspace{0.17em}\mathsf{\delta }{\theta }^{\mathrm{SD}}+{\theta }^{\mathrm{SD}}\hspace{0.17em}\mathsf{\delta }{e}_{zy}^{\mathrm{SD}}\right)]\mathrm{d}z\end{array}$$

(51)

Applying a linear interpolation for the axial displacement and a cubic interpolation for the deflections and twist rotations, respectively, the rigid-body displacements of the cross-section can be defined as:

$$\left\{\begin{array}{c}{w}_{\mathrm{O}}\\ u_{\mathrm{S}}\\ v_{\mathrm{S}}\\ {\phi }_z\\ {\phi }_x\\ {\phi }_y\\ \theta \end{array}\right\}=\left[\begin{array}{cccc}{N}_{w1}& 0& N_{w2}& 0\\ N_{u1}& N_{u2}& N_{u3}& N_{u4}\\ N_{v1}& N_{v2}& N_{v3}& N_{v4}\\ N_{{\phi }_{z}1}& N_{{\phi }_{z}2}& N_{{\phi }_{z}3}& N_{{\phi }_{z}4}\\ N_{{\phi }_{x}1}& N_{{\phi }_{x}2}& N_{{\phi }_{x}3}& N_{{\phi }_{x}4}\\ N_{{\phi }_{y}1}& N_{{\phi }_{y}2}& N_{{\phi }_{y}3}& N_{{\phi }_{y}4}\\ N_{\theta 1}& N_{\theta 2}& N_{\theta 3}& N_{\theta 4}\end{array}\right]\left[\begin{array}{ccccccc}{w}_{\mathrm{A}}& u_{\mathrm{A}}& v_{\mathrm{A}}& {\phi }_{z\mathrm{A}}& v_{\mathrm{A}}& u_{\mathrm{A}}& {\phi }_{z\mathrm{A}}\\ 1& {\phi }_{y\mathrm{A}}& {\phi }_{x\mathrm{A}}& {\theta }_{\mathrm{A}}& {\phi }_{x\mathrm{A}}& {\phi }_{y\mathrm{A}}& {\theta }_{\mathrm{A}}\\ w_{\mathrm{B}}& u_{\mathrm{B}}& v_{\mathrm{B}}& {\phi }_{z\mathrm{B}}& v_{\mathrm{B}}& u_{\mathrm{B}}& {\phi }_{z\mathrm{B}}\\ 1& {\phi }_{y\mathrm{B}}& {\phi }_{x\mathrm{B}}& {\theta }_{\mathrm{B}}& {\phi }_{x\mathrm{B}}& {\phi }_{y\mathrm{B}}& {\theta }_{\mathrm{B}}\end{array}\right],$$

(52)

where the shape-functions occurring in the first matrix on the right-hand side of Equation (43) are provided in Appendix B. Substituting Equations (36)–(43) into Equation (35), the incremental equilibrium equations for the beam element in C2 can be rewritten as:

$${}_1{}^{2}f{}^e={}_1{}^{1}f{}^e+\left({}^{1}k_{\mathrm{E}}^e+{}^{1}k_{\mathrm{G}}^e\right){}_{1}u{}^e,$$

(53)

where $k_{\mathrm{E}}^e$ is the incremental elastic stiffness matrix of the beam element while $k_{\mathrm{G}}^e$ is the corresponding incremental geometric stiffness matrix. To perform the force recovery, the so-called conventional approach [39] is adopted in this work, i.e.,

$${}_2{}^{2}f{}^e={}_1{}^{2}T{}^e\hspace{0.17em}\left[{}_1{}^{1}f{}^e+\left({}^{1}k_{\mathrm{E}}^e+{}^{1}k_{\mathrm{G}}^e\right){}_{1}u{}^e\right],$$

(54)

where ${}_1{}^{2}T{}^e$ denotes the incremental transformation matrix reported in [18]. After performing the assembly procedure, the incremental equilibrium of a beam-type structure can be obtained as:

$$\left({}^{1}K_{\mathrm{E}}+{}^{1}K_{\mathrm{G}}\right)\hspace{0.17em}{}_{1}U={}^{2}P-{}^{1}P,\hspace{1em}{}^{1}K_{\mathrm{E}}={\displaystyle \sum _e{\left({}^{1}t^e\right)}^{\mathrm{T}}\hspace{0.17em}{}^{1}k_{\mathrm{E}}^e\hspace{0.17em}{}^{1}t^e},\hspace{1em}{}^{1}K_{\mathrm{G}}={\displaystyle \sum _e{\left({}^{1}t{}^e\right)}^{\mathrm{T}}\hspace{0.17em}{}^{1}k_{\mathrm{G}}^e\hspace{0.17em}{}^{1}t^e},$$

(55)

where ¹K_E is the incremental elastic stiffness matrix, ¹K_G is the incremental geometric stiffness matrix, U is the incremental displacement vector, ¹P is the external load vector applied on the structure at C₁, ²P is the external load vector applied on the structure at C₂, and ¹t^e is the transformation matrix for the e-th beam element from the local coordinate system (¹z, ¹x, and ¹y) in C₁ to the global coordinate system (Z, X, and Y).

6. Numerical Examples

A computer program, named THINWALL v.17 [25], is developed on the basis of the numerical procedures presented in this paper. The nonlinear stability analysis is performed in the form of incremental-iterative solver based on the generalized control method [31]. To initiate the buckling, a small perturbation load should be introduced. Nodal coordinates and orientations of the cross-section axes of each beam element are updated at the end of each load step [25]. To stop the iteration, energy criteria is adopted with the tolerance value of 10⁻⁶.

Three beam models are considered in the analysis in order to emphasize the significance of the shear effects on the stability behavior of the structure under consideration. In the first beam model, the shear deformability effects are ignored entirely, and in the presented results, it is labelled with ‘SR’. The second beam model includes shear deformation coupling effects, using the procedures presented in this paper, and it is labelled with ‘SD’. The last beam model ignores the shear deformation coupling effects, and in the provided diagrams and tables, it is denoted by ‘SD1’. The structural material used in all the examples is graphite-epoxy (AS4/3501), with E₁ = 144 GPa, E₂ = 9.65 GPa, G₁₂ = 4.14 GPa, and v₁₂ = 0.3, as reported in [11]. The general cross-section shape used in the examples is shown in Figure 3. Each wall of the cross-section is of the same thickness t and is composed of single AS4/3501 ply. The ply orientations are provided in Figure 3.

6.1. Example 1

Figure 4a shows a cantilever column under an axial force F. A mono-symmetric cross-section of 6 × 10 × 6 × 1 cm (b₁ × h × b₂ × t) is adopted with the properties provided in

Table 1. Cross-section properties of Examples 1 and 3.

[0°/90°/0°]	$A^{\ast }$ (cm2)	$I_x^{\ast }$ (cm4)	$I_y^{\ast }$ (cm4)	$I_t^{\ast }$ (cm4)	$I_{\omega }^{\ast }$ (cm4)	xc (cm)	xg (cm)	xs (cm)
	22	532.343	72.525	7.333	1718.819	3.159	4.364	5.767
	Kx	Ky	Kω	Kyω	${\overline{Q}}_{11\mathrm{R}}$ (GPa)	${\overline{Q}}_{66\mathrm{R}}$ (GPa)
	2.1504	2.3969	0.0217	−0.4999	82.932	4.14

. A force F is applied at the geometric centroid of the cross-section. Since the neutral axis passes through the material-weighted centroid, the force acting through the geometric centroid behaves as the off-axis load. The lowest buckling load pertains to the torsional-flexural (TF) mode, as shown in Figure 4c. The second buckling mode is the flexural (F) mode, as shown in Figure 4b. Both modes are analyzed using five different mesh configurations, each consisting of one, two, four, eight, and sixteen beam elements. The results obtained are shown in Figure 5a and are compared with the NX Nastran’s laminate shell results.

As one can see, the shear locking does not occur in any of the shear deformable models and the SD model’s results approach the shell solution. The results obtained by the SD1 model underestimate the column buckling strength in TF mode, leaving the construction on the safe side if shear deformation effects are ignored. Since the beam model assumes a rigid cross-section and is unable to model the cross-section distortion, some differences between the beam model and the shell model are expected for higher loads and shorter beams where the cross-section distortion affects the critical load, as can be observed in the F mode results. Furthermore, the buckling strength reductions due to the shear deformability are 8% and 12% in comparison with the shear rigid model for the TF and F modes, respectively. The buckling strength increase due to the shear coupling effects is 6% for the TF mode, and there is no influence of the shear couplings for the F mode. In Figure 5b, the variations in the value of buckling load against the column slenderness ratio are shown. It can be seen that the differences between the results obtained by the SR and SD beam models are decreasing as the slenderness ratio increases.

6.2. Example 2

In this example, a cantilever beam from Figure 6a is considered. The lateral force F is acting at the free end in the Y-direction through the shear center of the asymmetric cross-section 4 × 10 × 2 × 1 cm. The cross-section properties are shown in

Table 2. Example 2’s cross-section properties.

[0°/90°/0°]	$A^{\ast }$ (cm2)	$I_x^{\ast }$ (cm4)	$I_y^{\ast }$ (cm4)	$I_t^{\ast }$ (cm4)	$I_{\omega }^{\ast }$ (cm4)	xs (cm)	ys (cm)	α (°)
	16	342.486	17.701	5.333	186.781	2.782	−1.769	6.404
	${\overline{Q}}_{11\mathrm{R}}$ (GPa)	${\overline{Q}}_{66\mathrm{R}}$ (GPa)	Kx	Ky	Kω	Kxy	Kxω	Kyω
	60.031	4.14	3.2748	1.6373	0.0358	−36.7315	−1.2571	−0.7571

.

Four different mesh configurations, each consisting of two, four, eight, and sixteen beam elements, are used to analyze the lateral-torsional buckling mode in Figure 6b. Since the cross-section is asymmetric, a difference in the buckling load for the positive and negative direction of the force is expected [38]. To initiate the buckling, a small perturbation force of intensity 0.001 F, acting through the shear center in the X-direction, is applied at the free end of the beam.

In Figure 7, the load-deflection curves are presented for the mesh configuration of the sixteen beam elements and for both of the force directions. The convergence study is shown in Figure 8. The critical forces obtained by NX Nastran’s laminate shell model are F_cr = 5.75 kN for the force acting in the +Y direction and F_cr = 13.873 kN for the force acting in the −Y direction. From the figures, it can be seen that all the mesh models recognize the buckling load. Significant differences appear in the post-buckling range above F > 3F_cr and F > 2.5F_cr, respectively. As can be seen, the reductions in the buckling strength due to the shear deformability are 2% and 8% for the positive and negative force directions, respectively.

6.3. Example 3

In this example, an L-frame from Figure 9a is considered. the frame is loaded with the force F acting through the geometric centroid of the cross-section at the free end. A mono-symmetric cross-section of 6 × 10 × 6 × 1 cm is adopted with the properties shown in Table 1. Warping is fully restrained at the frame corner. Four mesh configurations, each consisting of one, two, four, and eight beam elements per frame member, are used for the analysis of the lateral-torsional buckling mode, as shown in Figure 9b.

To initiate the buckling, a small perturbation force of intensity 0.001 F, acting in the X-direction, is introduced. The obtained results are shown in Figure 10. The critical force obtained by the NX Nastran’s laminate shell model is F_cr = 1.81 kN. From the figures, it can be observed the significant overestimation of the SR buckling strength differs from the SD model by 12%. In addition, it can be seen that all the mesh models recognize the buckling state very well.

7. Conclusions

In this paper, a refined shear deformable beam model for the nonlinear stability analysis of composite frame structures has been presented. The materials can vary over the walls of the cross-section, where unique expressions have been derived for the properties of such a composite cross-section. Shear deformations and shear deformation couplings have been included by the use of the six shear correction factors. The UL procedure has been used to describe the nonlinear behavior of such frame structures. The equilibrium equations for the two-node beam element have been derived on the basis of the nonlinear displacement field of the cross-section. In the incremental-iterative procedure, a generalized control method has been used where the nodal coordinates and orientations of the cross-section axes of each beam element are updated at the end of each load step. In the force recovery phase of the incremental-iterative solver, conventional approach based on the concept of the semi-tangential rotations have been adopted. Three benchmark examples have been presented and the results have been explained in detail. As expected, shear deformations can play a critical role in the design of composite frame structures, especially for structures with low slenderness ratios. For very short structures, cross-section distortion should be modeled in order to acquire accurate values for the critical force. In pure flexural and torsional buckling modes, there is no influence of the coupling effects on the shear deformations. In addition, from the presented results, it can be concluded that the presented shear deformable beam model experiences no shear locking.

The topic of our future activities will be to extend the proposed nonlinear finite element algorithm to large displacement problems of angle-ply laminates and functionally graded materials, with special emphasis on the buckling problems.