General Stiffness Matrix for a Thin-Walled, Open-Section Beam Structure

Abdelraouf M. Sami Alsheikh¹, D. W. A. Rees²

¹Department of Mechanical and Aerospace, Brunel University, London, UK.
²Department of Design, Brunel University, London, UK.

DOI: 10.4236/wjm.2021.1111015                

Abstract

This paper is to review the theory of thin-walled beam structures of the open cross-section. There is scant information on the performance of structures made from thin-walled beam elements, particularly those of open sections, where the behavior is considerably complicated by the coupling of tensile, bending and torsional loading modes. In the combined loading theory of thin-walled structures, it is useful to mention that for a thin-walled beam, the value of direct stress at a point on the cross-section depends on its position, the geometrical properties of the cross-section and the applied loading. This applies whether the thin-walled section is closed or open but this study will be directed primarily at the latter. Theoretical analyses of structures are fairly well established, considered in multi-various applications by many scientists. However, due to the present interest in lightweight structures, it is necessary to specify where the present theory lies. It does not, for example, deal with compression and the consequent failure modes under global and local buckling. Indeed, with the inclusion of strut buckling failure and any other unforeseen collapse modes, the need was perceived for further research into the subject. Presently, a survey of the published works has shown in the following: 1) The assumptions used in deriving the underlying theory of thin-walled beams are not clearly stated or easily understood; 2) The transformations of a load system from arbitrary axis to those at the relevant centre of rotation are incomplete. Thus, an incorrect stress distribution may result in; 3) Several methods are found in the recent literature for analyzing the behaviour of thin-walled open section beams under combined loading. These reveal the need appears for further study upon their torsion/flexural behaviour when referred to any arbitrary axis, a common case found in practice. This review covers the following areas: 1) Refinement to existing theory to clarify those observations made in 1 - 3 above; 2) Derivation of a general elastic stiffness matrix for combined loading; 3) Calculation of the stress distribution on the cross-section of a thin-walled beam. A general transformation matrix that accounts for a load system applied at an arbitrary point on the cross-section will be published in a future paper.

Keywords

Thin-Walled Open Sections, Shear Centre, Warping, Bi-Moment, Sectorial Area Properties

1. Introduction

Beams composed from thin plates are popular with designers, the advantages are that they are easy to produce and assemble. Moreover, their performance under different force systems gives high efficiency in terms of an optimal weight reduction to strength ratio. Several methods are found in the recent literature for analyzing the behaviour of thin-walled beams under combined loading [1] - [6]. From these, a need appears for further study upon their combined torsion/flexural behaviour. What investigated here specifically is the deformation response referred to any arbitrary axis, a common case found in practice.

Earlier work is instructive in this regard. Wagner [7] in 1936 recognized the significance of thin-walled structures, presenting the essential theory, which perhaps indicated to Vlasov [8]. The need to establish the main assumptions of elastic deformation of thin-walled beams as follows:

1) The cross-section is rigid (undeformable), so under bending moments, its middle section will move only as one body in its plane. Consequently, there are no normal and tangential stresses in the direction perpendicular to a given fixing.

2) Under twisting moments, points on an unrestrained cross-section will move at different rates in the axial direction, thus the plane of a cross-section does not remain plane, it becomes warped. This is in contrast with the Euler-Bernoulli assumption of plane sections remaining plane.

3) The application of a longitudinal force at an arbitrary position will cause warping displacement. Since this force may not be replaced by a statically equivalent longitudinal force, the beam will be subjected to a self-balancing set of longitudinal forces. This differs from the elementary theory of beam bending, which does not account for shear strain and stress due to warping.

4) The shearing deformations in the cross-section defined by the plane between the longitudinal axis (x) and the periphery (s) will be infinitesimal (see Figure 1).

The theory of Vlasov was published firstly in 1940, in Russian and translated into English, in 1961 [8]. The elements of Vlasov theory were also outlined by Zbirohowski-Ko’scia in 1967 [9] but apparently, these still remain unknown to many engineers, and scientists in the West, particularly as means to solve the load bearing capacity of thin-walled structures.

Figure 1. Warping of thin-walled open section under torsion.

Related studies by Krahula [10], Renton [11], Krajcinovic [12], Rajasekaran [13], Baigent, Hancock [14], and Al-Sheikh [15], employed matrix analysis for thin-walled beams, which depended upon the solution to the Vlasov’s differential equation. This present review is pertinent to more recent applications of Vlasov’s theory for thin-walled structures.

Chen and Hu [16] formed a general stiffness matrix for a thin-walled beam element, under combined torsion, and bending loads, based upon this beam’s linear, axial warping displacement proposed by Kollbrunner and Hajdin [17].

Dvorkin, et al. [18] presented a finite element for the analysis of a thin-walled open section beam structure. The element is based on Vlasov’s beam theory, which shows that for a linear elastic analysis, a single element is very effective for nonuniform torsion. Here no numerical integration is required, given that a single element can satisfy the equilibrium and compatibility conditions required.

The concept of a strip plate was introduced by Sorin and Bogdan, [19] to account for the effect of bending induced by torsion. The corresponding “macro” finite element matched experimental data more closely than could the element from classical theory provided.

Wen and Kuo [20] provided a general function for the torsional warping of a thin-walled open-section beam. To achieve this, the assumptions made in the Vlasov theory were combined with those assumptions made in Kirchhoff’s plate/shell theory.

Emre Erkmen [21] claimed a new technique to solve the equations of equilibrium for a thin-walled structure under combined loading. An orthogonal coordinate system was adopted to uncouple the differential equations that apply to a thin-walled beam element. Therein, lay the possibility to introduce several modeling schemes with different limitations specified in each case.

Vetyuko [22] advanced a direct approach based on the principles of Lagrangean mechanics. A reduction to small deformations led to linear equations for a stability analysis of a beam subjected to axial or transversal loading.

Wang et al. [23] established a relationship between the warping torque and the restrained shear rotation for torsion of an open thin-walled beam. An energy method was used to develop, to a first-order, a relation based on Vlasov’s theory.

To describe the transversal deformation of a thin-walled beam cross-section with its axial warping, a new formulation of a 3D beam element was presented by Khalil and Cespedes [24]. Using “enriched kinematics”, overcame the error within the assumption that the plane section remains plane under torsion.

An approach based upon a higher-order beam theory was proposed by Choi et al. [25]. This was applied to a thin-walled tubular beam with a varying quadrilateral cross-section. Thereby, they accounted for a lateral distortion induced by the Poisson’s effect and axial warping displacement.

A comprehensive study made by Drillenburg [26], verified the manual calculation of torsional stresses and deformations from the classical theory as compared to the output of a frame analysis program “SCIA Engineer”. Very impressive hand calculations, included in the report, agreed well with those of the program [26].

Finally, a deformable cross-section, based upon generalized nonlinear beam theory, was presented by Duan and Zhao [27]. This theory applies to rigorous numerical analyses of thin-walled structures. Here it was employed for elastic and elastic-plastic analyses of thin-walled members undergoing arbitrary deformations, such as large deflections, finite rotations, distortion within local buckling and out-of-plane warping.

2. Torsion Warping Theory

2.1. List of Symbols

A is the centre of twist for the cross section.

$a_y,a_z,y\text{and}z$ are the Cartesian coordinates of points A, and S, respectively.

$v,w,v_s\text{and}{w}_s$ are the displacements in the y- and z-directions of points A, and S, respectively.

$\bar{\eta },\bar{\xi }$ are the tangential and perpendicular displacements of the point S, to the contour.

u is the longitudinal displacement of an arbitrary point S, of Cartesian coordinate (y, z) and sectorial coordinate $\omega$.

u_o is the integration constant, which represents the longitudinal displacement of the origin of the peripheral coordinate, s, at x = const. Physically, u_o is the longitudinal displacement of the beam due to pure extension only, depending on the variable x.

$v^{\prime },w^{\prime }$ are rates of change of displacement of point S, (with respect to X-axis), in the y and z-directions respectively.

${\phi }^{\prime }$ is the rate of change of twist of the beam with respect to the X-axis.

${\phi }^{\prime }\omega$ is a product giving longitudinal warping displacement of point S due to the rate of change of twist.

N is the axial force. ${\bar{M}}_y$ is the total bending moment about Y-axis, ${\bar{M}}_z$ is the total bending moment about Z-axis, A is the total area of the cross section and t is the wall thickness of the thin-walled beam.

${\displaystyle {\int }_{A}y\text{d}A}=S_z$ first moment of area of the cross-section about Z-axis.

${\displaystyle {\int }_{A}z\text{d}A}=S_y$ first moment of area of the cross-section about Y-axis.

${\displaystyle {\int }_A{y}^2\text{d}A}=I_z$ second moment of area of the cross-section about Z-axis.

${\displaystyle {\int }_A{z}^2\text{d}A}=I_y$ second moment of area of the cross-section about Y-axis.

${\displaystyle {\int }_{A}yz\text{d}A}=I_{yz}$ product moment of area of the cross-section

${\displaystyle {\int }_A\omega \text{d}A}=S_{\omega }$ first moment of sectorial area.

${\displaystyle {\int }_A\omega z\text{d}A}=S_{\omega z}$ product moment of sectorial area about Y-axis.

${\displaystyle {\int }_A\omega y\text{d}A}=S_{\omega y}$ product moment of sectorial area about Z-axis.

${\displaystyle {\int }_A{\omega }^2\text{d}A}=\Gamma$ second moment of sectorial area.

2.2. Theoretical Summary

When a thin-walled beam of open-section is subjected to torsion (i.e. twisting moments) at its ends, points on the cross-section of the beam move at different rates in the axial direction (see Figure 1). Thus, the plane of a cross-section becomes warped. If one or more cross-sections are restrained against warping, the state of stress in the beam will be totally different to that analyzed by St. Venant’s theory of unrestrained torsion [28] [29]. There the thin walled beam can undergo longitudinal extension as a result of warping only. Otherwise, with restraints to free warping, longitudinal stresses are created. Vlasov [8] (pp. 4-5, p. 11) showed that there are also shear stresses produced in such a beam subjected only to a longitudinal force applied at an arbitrary point on the cross-section. These shear stresses vary according to distance away from the point at which the force is applied. To make this problem clear, consider a longitudinal compressive force P applied to the free-end of a thin walled I section cantilever. In Figure 2(a), P is positioned at the right end of the top I-beam flange (point n).

(a)(b)(c)(d)(e)

Figure 2. Cantilever of I-cross section loaded at point n.

It can be seen in Figures 2(b)-(e) how the cantilever suffers four different types of stress in the axial direction:

1) Normal stress due to the four direct forces P/4, applied at the far ends of the cross-section, (Figure 2(b)).

2) Bending stresses due to bending moment about the Z-axis (Figure 2(c)) and about the Y-axis (Figure 2(d)).

3) Stresses due to skew-symmetry of four forces of value P/4, two bending moments equal in magnitude, but opposite in sign, acting in the plane of the two flanges of the I-beam. These stresses will be known as warping stresses (Figure 2(e)).

The splitting of the force within the four diagrams provides a resultant force P as it is applied to the position n. Bending moment resultants equal Pa and Pb about the axes z and y respectively.

3. Displacements and Strains on Thin-Walled Beams of Open Section

3.1. Out of Plane and Lateral Displacements

As was stated by Vlasov [8] (pp. 7-9), the assumption of a non-deformable cross-section will lead to the fact that any out of plane displacements of the wall of a cross-section must vanish.

Also, lateral displacements in the z- and y-directions, as well as any rotations the beam may undergo, are functions of the longitudinal direction x only. In addition to this the axial displacement of a beam will be a function of both the longitudinal axis x and the profile coordinate of the cross-section s. It is also assumed that the x-axis and s remain orthogonal while the beam undergoes deformations, i.e., $\gamma =0$ (see Figure 3).

3.2. The Kinematics of Displacements

If a thin-walled beam of open section undergoes deformation, the arbitrary point B on the profile Figure 4, which is rigidly connected, connected to the instantaneous centre of twist A, results the following kinematical relations for both points displacements:

$v_b=v_A-\left(b_y-a_y\right)+\left[\left(b_z-a_z\right)\sin \phi +\left(b_y-a_y\right)\cos \phi \right]$

$w_b=w_A-\left(b_z-a_z\right)+\left[\left(b_z-a_z\right)\cos \phi +\left(b_y-a_y\right)\sin \phi \right]$

If $\phi$ is small then, $\cos \phi \approx 1$, and $\sin \phi \approx \phi$. Hence,

Figure 3. Orthogonal axes for thin shell element.

Figure 4. Displacements and rotations of Element.

$v_b=v+\left(b_z-a_z\right)\phi$ (1)

$w_b=w-\left(b_y-a_y\right)\phi$ (2)

In last two equations v, and w, are the coordinate’s values of the instantaneous centre of twist A, in y, and z directions respectively. The displacements of an arbitrary point S, in the middle surface of the cross-section when moved to point S’, Figure 5, are as follows:

$v_s=v+\left(z-a_z\right)\phi$ (3)

$w_s=w-\left(y-a_y\right)\phi$ (4)

From Figure 5, the following kinematic relations apply:

$r=\left(z-a_z\right)\sin \alpha -\left(y-a_y\right)\cos \alpha$ (5)

$r_n=\left(z-a_z\right)\cos \alpha -\left(y-a_y\right)\sin \alpha$ (6)

$\bar{\eta }=w_s\cos \alpha -v_s\sin \alpha$ (7)

$\bar{\xi }=w_s\sin \alpha -v_s\cos \alpha$ (8)

Now substituting Equations (3) and (4) into Equations (7) and, (8). We have:

$\bar{\eta }=\left[w-\left(y-a_y\right)\phi \right]\cos \alpha +\left[v+\left(z-a_z\right)\phi \right]\sin \alpha$ (9)

$\bar{\xi }=\left[w-\left(y-a_y\right)\phi \right]\sin \alpha -\left[v+\left(z-a_z\right)\phi \right]\cos \alpha$ (10)

Figure 5. Displacements and rotations of the contour of a thin–walled Beam.

Rewriting Equations (9) and (10):

$\bar{\eta }=v\sin \alpha +w\cos \alpha +\left[\left(z-a_z\right)\sin \alpha -\left(y-a_y\right)\cos \alpha \right]\phi$ (11)

$\bar{\xi }=-v\cos \alpha +w\sin \alpha -\left[\left(z-a_z\right)\cos \alpha +\left(y-a_y\right)\sin \alpha \right]\phi$ (12)

Recalling Equations (5) and (6) and substituting into Equations (11) and (12) we get:

$\bar{\eta }=v\sin \alpha +w\cos \alpha +r\phi$ (13)

$\bar{\xi }=-v\cos \alpha +w\sin \alpha -r_n\phi$ (14)

Orthogonality between the longitudinal X-axis and the peripheral co-ordinate s, can be explained by observing the displacements of line elements cut-out from the beam in Figure 6. The sum of angles ${\gamma }_1$ and ${\gamma }_2$ are the total change of

angle between both axes x and s, ${\gamma }_1=\frac{\partial u}{\partial s}$, and ${\gamma }_2=\frac{\partial \bar{\eta }}{\partial x}$. Their sum describes the engineering shear strain in radians:

$\gamma ={\gamma }_1+{\gamma }_2=\frac{\partial u}{\partial s}+\frac{\partial \bar{\eta }}{\partial x}$ (15)

Given that the total shear angle $\gamma$ is very small, and can be neglected, then

$\gamma =0$. This assumption is quite feasible physically since it means that, $\frac{\partial u}{\partial s}$, and, $\frac{\partial \bar{\eta }}{\partial x}$ are, sensibly, equal in magnitude. Their very small difference allows us to write the following equation:

$\frac{\partial u}{\partial s}=-\frac{\partial \bar{\eta }}{\partial x}$ (16)

Figure 6. Shear deformations of shell line elements.

Now apply the expression for $\bar{\eta }$ from Equation (13) and find its rate of change w r t, the variable x:

$\frac{\partial \bar{\eta }}{\partial x}=v^{\prime }\sin \alpha +w^{\prime }\cos \alpha +{\phi }^{\prime }r$ (17)

where prime designates ( $\partial /\partial x$ ), noting that from Figure 7, the following relations,

$\text{d}y=\sin \alpha \text{d}s$, $\text{d}z=\cos \alpha \text{d}s$ and $\text{d}\omega =r\text{d}s$

From Equation (16), and (17) we get:

${\displaystyle \int \text{d}u}=-{\displaystyle \int \left(v^{\prime }\text{d}y+w^{\prime }\text{d}z+{\phi }^{\prime }\text{d}\omega \right)}$ (18)

It should be obvious that, v, w, and $\phi$ are functions only of x, whereas u is a function of both, x and s. The integration constant for u, from Equation (18), can be considered as a function $f\left(x\right)$, resulting in following equation:

Here $f\left(x\right)$, the integration constant, is given the symbol, $u_o\left(x\right)$, so that

$u\left(s,x\right)=u_o\left(x\right)-v^{\prime }\left(x\right)y-w^{\prime }\left(x\right)z-{\phi }^{\prime }\left(x\right)$ (20)

Equation (20) was derived using only kinematic relations of Equations (5)-(8). It reveals the following:

- The first three terms express the Bernoulli-Navier law, in which any plane cross section before deformation, remains plane after deformation.

- The fourth term in this equation, is the relative longitudinal displacement due to warping of the beam cross-section, where, u is the longitudinal displacement of an arbitrary point S having Cartesian coordinates (y, z) and sectorial coordinate $\omega$.

Figure 7. Sectorial semi-area dω/2 for coordinate point S.

Recall that u_o is the integration constant, which represents the longitudinal displacement of the origin of the peripheral coordinate, s, at x = constant. This is equivalent to u_o as the longitudinal displacement of the beam due to pure extension only, depending on the variable x. In Equations (19) and (20) v’, w’ are rates of change of the displacements u and $\bar{\eta }$ of point S w. r. t. x in the y- and z-directions respectively. Also, ${\phi }^{\prime }$ is the rate of change of twist of the beam with respect to the X-axis. Hence the product ${\phi }^{\prime }\omega$ represents the longitudinal displacement of point S due to the rate of change of twist. Note here that the sectorial coordinate ω is a function of the peripheral coordinate s, so the magnitude of the warping displacement at x = const, is ${\phi }^{\prime }\omega$.

This shows that, having assumed that shear strain is absent, the warping varies linearly on the cross-section according to the law of sectorial area. Let us continue in this derivation to justify that assumption, but before doing so, it is very important to know and understand the law of sectorial area and its characteristic properties.

4. The Law of Sectorial Area

The Determination of Sectorial Area Expressions

Let ${\omega }_A$ and ${\omega }_D$ be sectorial areas of the profile shown in Figure 8, their poles are A and D respectively. The sectorial area is twice the area swept by the sectorial radius r from the point $S_o\left(y_o,z_o\right)$ to the point $S_o\left(y_o+\Delta y,z_o+\Delta z\right)$, where S_o is the point from which the profile co-ordinate s is measured, i.e. ${\omega }_A={\omega }_D=0$, at s = 0 (see Figure 9).

Now from Figure 9 twice the shaded area is given by the following equation:

$\text{d}{\omega }_A=\left(z_o-a_z\right)\text{d}y-\left(y_o-a_y\right)\text{d}z$ (21)

Figure 8. Sectorial coordinate of point S.

Figure 9. Geometry of sectorial coordinate of point S.

where a_y, a_z are the coordinates of pole A referred to the coordinate system OYZ, in Figure 8. Similarly, for pole D (see Figure 8):

$\text{d}{\omega }_D=\left(z_o-d_z\right)\text{d}y-\left(y_o-d_y\right)\text{d}z$ (22)

Integrating Equations (21) and (22) we get:

${\omega }_A=\left(z_o-a_z\right)y-\left(y_o-a_y\right)z+C_1$ (23)

${\omega }_D=\left(z_o-d_z\right)y-\left(y_o-d_y\right)z+C_2$ (24)

Since C₁ and C₂ are the integration constants, they are found from imposing the condition:

${\omega }_A=0$, at $y=y_o$, and $z=z_o$

$0=\left(z_o-a_z\right)y_o-\left(y_o-a_y\right)z_o+C_1\text{or}{C}_1=a_z{y}_o-a_y{z}_o$

Substituting the expression for C₁ into Equation (23), we have:

${\omega }_A=\left(z_o-a_z\right)\left(y-y_o\right)-\left(y_o-a_y\right)\left(z-z_o\right)$ (25)

which is consistent with the conditions mentioned before: ${\omega }_A=0$, at $y=y_o$, and $z=z_o$. Similarly, for ${\omega }_D$ :

${\omega }_D=\left(z_o-d_z\right)\left(y-y_o\right)-\left(y_o-d_y\right)\left(z-z_o\right)$ (26)

Subtracting Equation (26) from Equation (25) we get:

${\omega }_A={\omega }_D+\left(a_y-d_y\right)\left(z-z_o\right)-\left(a_z-d_z\right)\left(y-y_o\right)$ (27)

For the calculation of ${\omega }_A$ and ${\omega }_D$ we have A and D as two arbitrary poles. Now, certain conditions will be imposed on a pole point, which we will call the orthogonality conditions. They are:

$\begin{array}{l}{S}_{\omega }={\displaystyle {\int }_A\omega \text{d}A}=0\\ S_{\omega y}={\displaystyle {\int }_A\omega y\text{d}}A=0\\ S_{\omega z}={\displaystyle {\int }_A\omega z\text{d}A}=0\\ \Gamma ={\displaystyle {\int }_A{\omega }^2\text{d}A}\end{array}$ (28)

Let us call the pole A, the principal sectorial pole when satisfying Equations (28). Pole D is any arbitrary point, found from multiplying Equation (27) by dA and integrating:

${\displaystyle {\int }_A{\omega }_A\text{d}A}={\displaystyle {\int }_A{\omega }_D\text{d}A}+\left(a_y-d_y\right){\displaystyle {\int }_A\left(z-z_o\right)\text{d}A}-\left(a_z-d_z\right){\displaystyle {\int }_A\left(y-y_o\right)\text{d}A}$ (29)

If the origin of the coordinate system OYZ coincides with the centroid of the cross-section we have:

${\displaystyle {\int }_{A}z\text{d}A}={\displaystyle {\int }_{A}y\text{d}A}=0$

and Equation (29) reads:

$\frac{S_{\omega D}}{A}=\left(a_y-d_y\right)z_o-\left(a_z-d_z\right)y_o$ (30)

Substituting Equation (30) into (27) we have:

${\omega }_A={\omega }_D+\left(a_y-d_y\right)z-\left(a_z-d_z\right)y-\frac{S_{\omega D}}{A}$ (31)

Also, multiplying Equation (27) by ydA and zdA respectively and integrating over the whole area of the cross-section A, we have:

${\displaystyle {\int }_A{\omega }_{A}y\text{d}A}={\displaystyle {\int }_A{\omega }_{D}y\text{d}A}+\left(a_y-d_y\right){\displaystyle {\int }_A\left(z-z_o\right)y\text{d}A}-\left(a_z-d_z\right){\displaystyle {\int }_A\left(y-y_o\right)y\text{d}A}$ (a)

and:

${\displaystyle {\int }_A{\omega }_{A}z\text{d}A}={\displaystyle {\int }_A{\omega }_{D}z\text{d}A}+\left(a_y-d_y\right){\displaystyle {\int }_A\left(z-z_o\right)z\text{d}A}-\left(a_z-d_z\right){\displaystyle {\int }_A\left(y-y_o\right)z\text{d}A}$ (b)

Now imposing the integrals in Equations (28) the Equations (a), and (b), will read:

$a_y-d_y=-\left(\frac{S_{\omega Dz}{I}_z-S_{\omega Dy}{I}_{yz}}{I_y{I}_z-I_{yz}^2}\right)$ (32)

$a_z-d_z=\left(\frac{S_{Dy}{I}_y-S_{Dz}{I}_{yz}}{I_y{I}_z-I_{yz}^2}\right)$ (33)

${\Gamma }_A={\Gamma }_D+\left(a_y-d_y\right)S_{\omega Dz}-\left(a_z-d_z\right)S_{\omega Dy}-\frac{S_{\omega D}^2}{A}$ (34)

5. Derivation of the Normal Stress Formulae

Assumptions of Elastic Beams and the Stress Formulae

It was assumed that the displacements of the thin-walled beams of open section are too small in comparison with the cross-sectional dimensions [8], thus the assumptions of elastic beams are valid, hence:

$\sigma =\epsilon E$ (35)

where $\sigma$ is the normal stress, $\epsilon$ is the strain, and E is Young’s modulus. Now recall Equation (20) along with $\epsilon =\frac{\partial u}{\partial x}$. Then, from Equation (35) we get:

$\sigma =E\left({u^{\prime }}_o-v^{″}y-w^{″}z-{\phi }^{″}\omega \right)$ (36)

where: $\text{'}\equiv \frac{\partial }{\partial x}$ and $\text{'}\text{'}\equiv \frac{{\partial }^2}{\partial x^2}$

The normal stresses $\sigma ={\sigma }_x$ in Equation (36) is a function of x and s (see Figure 10). Here one can introduce the system of generalized forces and moments acting on cross-section at x = const. As follows:

$N={\displaystyle {\int }_A\sigma \text{d}A}$ (37)

${\bar{M}}_y=–{\displaystyle {\int }_A\sigma z\text{d}A}$ (38)

${\bar{M}}_z={\displaystyle {\int }_A\sigma y\text{d}A}$ (39)

Figure 10. Axial stress within a thin-walled beam cross-section.

where,

N is the axial force, ${\bar{M}}_y$ is the total bending moment about the Y-axis, ${\bar{M}}_z$ is the total bending moment about Z-axis, A is the total area of the cross section, and t is the wall thickness.

$\text{d}A=t\text{d}s$ (40)

Given that the second moment of area represents the bending stiffness for a cross section of a structural member, the final integral represents the warping stiffness of the cross section of that member.

Substituting the value of normal stress $\sigma$ from Equation (36) into Equations (37)-(39), having assumed that u_o, v, w and $\phi$ are functions of x only, this gives:

$N=E\left({u^{\prime }}_{o}A-{\phi }^{″}{S}_{\omega }\right)$ (41)

${\bar{M}}_y=E\left(v^{″}{I}_{yz}+w^{″}{I}_y+{\phi }^{″}{S}_{\omega z}\right)$ (42)

${\bar{M}}_z=-E\left(v^{″}{I}_z+w^{″}{I}_{yz}+{\phi }^{″}{S}_{\omega y}\right)$ (43)

Solving Equations (41)-(43) for ${u^{\prime }}_o$, $v^{″}$, and ^{$w^{″}$} , then substituting into Equation (37), gives the normal stress in a general form:

$\begin{array}{c}\sigma =\frac{N}{A}+\frac{{\bar{M}}_z{I}_y+{\bar{M}}_y{I}_{yz}}{I_y{I}_z-I_{yz}^2}y-\frac{{\bar{M}}_y{I}_z+{\bar{M}}_z{I}_{yz}}{I_y{I}_z-I_{yz}^2}z\\ -\left({\omega }_D–\frac{S_{\omega yD}{I}_y-S_{\omega zD}{I}_{yz}}{I_y{I}_z-I_{yz}^2}y-\frac{S_{\omega zD}{I}_z-S_{\omega yD}{I}_{yz}}{I_y{I}_z-I_{yz}^2}z–\frac{S_{\omega D}}{A}\right)E{\phi }^{″}\end{array}$ (44)

for which the sectorial coordinate $\omega$ was chosen arbitrarily. Here Equation (44) applies to the normal stress in a beam section that rotates about an arbitrary point D in which ${\omega }_D$, $S_{\omega yD}$ , $S_{\omega zD}$, $S_{\omega D}$ appear. This approach is different to that given by Vlasov [8]. Specifically, he took $\bar{M}$ to represent the total external bending moment on the cross-section at x = constant. This includes the bending moments due to the axial and transverse forces applied at positions away from centroid, in addition to the direct bending moment. Here the terms between two parentheses, in the Equation (44) equal to the sectorial area $\omega$ of the profile, calculated with respect to, the centre of rotation. Substituting from Equations (31)-(33), simplifies Equation (44):

$\sigma =\frac{N}{A}+\frac{{\bar{M}}_z{I}_y+{\bar{M}}_y{I}_{yz}}{I_y{I}_z-I_{yz}^2}y-\frac{{\bar{M}}_y{I}_z+{\bar{M}}_z{I}_{yz}}{I_y{I}_z-I_{yz}^2}z-E\omega {\phi }^{″}$ (45)

for which S_y, and S_z vanish when the origin of the Cartesian system lies at the centroid of the beam cross-section. Multiplying both sides of Equation (45) by $\omega \text{d}A$, and integrating over the whole area of the cross-section enables the application of the orthogonality condition from Equations (28). That is:

${\displaystyle {\int }_A\sigma \omega \text{d}A}=-E\Gamma {\phi }^{″}$ (46)

Let us call the integral ( ${\displaystyle {\int }_A\sigma \omega \text{d}A}$ ) the bi-moment and give it a symbol B, hence,

$B=-E\Gamma {\phi }^{″}\to -{\phi }^{″}=\frac{B}{E\Gamma }$ (47)

Substituting ${\phi }^{″}$ into Equation (45) we have the final form of general normal stress expression in a thin-walled beam of open section:

$\sigma =\frac{N}{A}+\frac{{\bar{M}}_z{I}_y+{\bar{M}}_y{I}_{yz}}{I_y{I}_z-I_{yz}^2}y-\frac{{\bar{M}}_y{I}_z+{\bar{M}}_z{I}_{yz}}{I_y{I}_z-I_{yz}^2}z+\frac{B}{\Gamma }\omega$ (48)

when the coordinate system OYZ coincides with the principal axes of the section area, then $I_{yz}=0$, and Equation (48), can be written in the following simpler form:

$\sigma =\frac{N}{A}+\frac{{\bar{M}}_z}{I_z}y-\frac{{\bar{M}}_y}{I_y}z+\frac{B}{\Gamma }\omega$ (49)

Equation (49) was in fact derived by Vlasov [8], by applying the orthogonality conditions from the beginning considering that the sectorial coordinate $\omega$ is associated with the principal sectorial pole. In contrast, the present approach starts with an arbitrary sectorial pole D and through mathematical verification leads directly to this result. Another point which can be mentioned is that the definition of the bi-moment transpires from Equations (46), and (47) whereas Vlasov had to define the bi-moment one stage before. The usefulness of this approach can be seen, by establishing the equation of total bi-moment when a force system is applied at a point away from the principal sectorial pole.

6. Stress Distribution on Thin Walled Beams of Open Section

6.1. Assumptions

The distribution of stress in an open-section is quite different from that in solid or closed section. A few assumptions are made since the thickness of the wall is very small compared to the outer dimension, Vlasov [8], Oden and Ripperger [28]. The first assumption is that the normal stresses ${\sigma }_x$ (Figure 11(a)), is uniform over the wall thickness. The second assumption is that the stresses and strains normal to the wall surface are very small and hence they are neglected, i.e.:

${\sigma }_n={\tau }_{xn}={\tau }_{sn}={\epsilon }_n={\gamma }_{xn}={\gamma }_{sn}=0$ (50)

Only one shearing stress is considered, this is τ_xs, (Figure 11(b)), which is due to two different modes of deformation, Vlasov [8]. The first mode is due to external torsional moments, transverse loads, with non-uniform axial deformation of a free-to-warp cross section, which result in a linear distribution of shearing stresses over the wall thickness (see Figure 11(c)). This will be given the symbol ${\left({\tau }_{xs}\right)}_T$. The other mode of deformations is due to the lateral shear forces in the direction of the tangent on the contour arc (see Figure 11(d)). This will result in uniform tangential stresses over the wall thickness and will be given the symbol ${\left({\tau }_{xs}\right)}_V$. Hence the total shear stress (Figure 11(b)) over the cross-section will be:

${\tau }_{xs}={\left({\tau }_{xs}\right)}_T+{\left({\tau }_{xs}\right)}_V$ (51)

Figure 11. Distribution of shear stresses on thin-walled shell.

6.2. Determination of Tangential Stress

The condition of equilibrium of an infinitesimally small shell element of the beam will be used to determine the relationship between the normal stresses and tangential stresses as seen in Figure 12:

$\partial \left(\sigma t\right)\text{d}s+\partial \left(\tau t\right)\text{d}x+P_x\text{d}x\text{d}s=0$ (52)

Dividing Equation (52) by dx, and integrating equation along the profile S, we get:

$\tau =\left(\frac{1}{t}\right)\left[S_o\left(x\right)-{\displaystyle {\int }_s{P}_x\text{d}S}-{\displaystyle {\int }_s\frac{\partial \sigma }{\partial x}t\text{d}s}\right]$ (53)

where: $P_x=P_x\left(x,s\right)$ is the projection of an external surface load per unit area in x-direction and S_o(x) is an integration constant which can be determined from the initial conditions, at S = 0. Thus,

$\tau \left(x,0\right)=\left(\frac{1}{t\left(0\right)}\right)S_o\left(x\right)\text{or}{S}_o\left(x\right)=t\left(0\right)\cdot \tau \left(x,0\right)$

where t = t(s) is the wall thickness. Recalling Equation (36) and substituting into Equation (53) we get:

$\tau =\left(\frac{1}{t}\right)\left[S_o\left(x\right)-{\displaystyle {\int }_s{P}_x\text{d}S}-E{u^{″}}_{o}A\left(s\right)+E{v}^{‴}{s}_z\left(s\right)+E{w}^{‴}{s}_y\left(s\right)+E{\phi }^{‴}{s}_{\omega }\left(s\right)\right]$ (54)

where: $\text{d}A=t\text{d}S$

For a general equation, recalling Equation (48) and substitute it into Equation (54), we get:

$\begin{array}{c}\tau =\left(\frac{1}{t}\right)[S_o\left(x\right)-{\displaystyle {\int }_s{P}_x\text{d}S}-\frac{N^{\prime }}{A}A\left(s\right)-\frac{{{\bar{M}}^{\prime }}_z{I}_y+{{\bar{M}}^{\prime }}_y{I}_{yz}}{I_y{I}_z-I_{yz}^2}{s}_z\left(s\right)\\ +\frac{{{\bar{M}}^{\prime }}_y{I}_z+{{\bar{M}}^{\prime }}_z{I}_{yz}}{I_y{I}_z-I_{yz}^2}{s}_y\left(s\right)-\frac{B^{\prime }}{\Gamma }{s}_{\omega }\left(s\right)]\end{array}$ (55)

where: $N^{\prime }=\frac{\text{d}N}{\text{d}x}$, ${\bar{M}}^{\prime }=\frac{\text{d}\bar{M}}{\text{d}x}$, $B^{\prime }=\frac{\text{d}B}{\text{d}x}$ For a useful but less general expression

Figure 12. Shell element under in-plane forces.

assume $N^{\prime }=0$, for a constant axial force along the X-axis, also if the axes y, z are the principal axes, $I_{yz}=0$, hence Equation (55) can be written as follows:

$\tau =\left(\frac{1}{t}\right)\left[S_o\left(x\right)-{\displaystyle {\int }_0^s{P}_x\text{d}S}-\frac{{{\bar{M}}^{\prime }}_z}{I_y}{S}_z\left(s\right)+\frac{{{\bar{M}}^{\prime }}_y}{I_z}{S}_y\left(s\right)-\frac{B^{\prime }}{\Gamma }{S}_{\omega }\left(s\right)\right]$ (56)

We note that in Equations (55) and (56) the values S_z(s), S_y(s), and $S_{\omega }\left(s\right)$ are the values of limited integration between s = 0, which represents the starting point on the profile, up to an arbitrary distance s on the cross-section, at which the shearing stresses $\tau$ is to be determined. Also, $S_{\omega }$ is referred to the principal sectorial pole (shear centre). The thin-walled beam of open section is subjected to a general system of transverse loads, and a distributed external torque of intensity m per unit length. Since the bi-moment B on the structure is self-equilibrating, the equations of equilibrium for the rest of the forces will not be affected (see Figure 13).

Thus,

$\begin{array}{l}\frac{\text{d}{Q}_y}{\text{d}x}=-P_y,\frac{\text{d}{Q}_z}{\text{d}x}=-P_z\\ \frac{\text{d}{\bar{M}}_z}{\text{d}x}=Q_y,\frac{\text{d}{\bar{M}}_y}{\text{d}x}=Q_z\end{array}$ (57)

Substituting Equation (57) into Equations (55) and (56) we have:

$\begin{array}{l}\tau =\left(\frac{1}{t}\right)[S_o\left(x\right)–{\displaystyle {\int }_0^s{P}_x\text{d}s}-\frac{N^{\prime }}{A}A\left(s\right)-\frac{Q_y{I}_y+Q_z{I}_{yz}}{I_y{I}_z-I_{yz}^2}{S}_z\left(s\right)\\ +\frac{Q_z{I}_z+Q_y{I}_{yz}}{I_y{I}_z-I_{yz}^2}{S}_y\left(s\right)-\frac{T_{\omega }}{\Gamma }{S}_{\omega }\left(s\right)]\end{array}$ (58)

where $T_{\omega }$ from Equation (58), is defined by the following:

$T_{\omega }=B^{\prime }=\frac{\text{d}B}{\text{d}x}=–E\Gamma {\phi }^{‴}$ (59)

Figure 13. Open section thin-walled element, under combined load.

7. Differential Equations of Equilibrium for a Thin-Walled Beam

The state of equilibrium of a beam under a set of internal stresses and external forces can be shown in Figure 14, where an element of thickness t, profile length Δs, and axial length Δx, was examined for equilibrium under the prescribed forces and stresses, which include internal normal and shearing stresses, external torsional moment and body forces. These equilibrium conditions are:

${\displaystyle \sum F_x}={\displaystyle {\int }_s\frac{\partial \left(\sigma t\right)}{\partial x}\text{d}x\text{d}s}+\left(q_L-q_K+P_x\right)\text{d}x=0$ (60)

${\displaystyle \sum F_y}={\displaystyle {\int }_s\frac{\partial \left(\tau t\right)}{\partial x}\text{d}x\text{d}s\sin \alpha }+P_y\text{d}x=0$ (61)

${\displaystyle \sum F_z}={\displaystyle {\int }_s\frac{\partial \left(\tau t\right)}{\partial x}\text{d}x\text{d}s\cos \alpha }+P_z\text{d}x=0$ (62)

${\displaystyle \sum M_D}={\displaystyle {\int }_s\frac{\partial \left(\tau t\right)}{\partial x}\text{d}x\text{d}s\left[\left(z_o-d_z\right)\sin \alpha -\left(y_o-d_y\right)\cos \alpha \right]}+{T^{\prime }}_v\text{d}x=0$ (63)

The first element of Equation (60) gives the total change in magnitude between axial forces exerted on the shell element (dx × ds), and P_x is the axial surface load applied on the whole cross-section of unit width, where q_L and q_k are the shear stresses acting on the both extreme edges of the element in direction x, i.e., q_L = q_L(x), and q_k = q_k(x).

Figure 14. External and internal forces on shell element.

The first element in Equations (61) and (62) are the projections of the differential of the shear force upon the y- and z-axes, where, P_y = P_y(x) and P_z = P_z(x) are external of transverse loads in the y- and z-directions.

Equation (63) represents the equilibrium of moments about arbitrary centre of rotation D, while the first term in this equation represents the moment of the shear stresses acting on the sides of the element (dx × ds), tangential to contour S. The second term is the differential (primed) of torsional moment T_v on the shell element of surface area (dx × ds).

We have seen from Figure 7; the following relationships apply:

$\text{d}y=\sin \alpha \text{d}s,\text{d}z=\cos \alpha \text{d}s\text{and}\text{d}{\omega }_D=\left(z_o-d_z\right)\text{d}y-\left(y_o-d_y\right)\text{d}z$

Dividing Equations (60)-(63), by dx, and taking into consideration these relationships, shows:

${\displaystyle {\int }_L\frac{\partial \left(\sigma t\right)}{\partial x}\text{d}s}+P_x+q_L+q_K=0$ (64)

${\displaystyle {\int }_L\frac{\partial \left(\tau t\right)}{\partial x}\text{d}y}+P_y=0$ (65)

${\displaystyle {\int }_L\frac{\partial \left(\tau t\right)}{\partial x}\text{d}z}+P_z=0$ (66)

${\displaystyle {\int }_L\frac{\partial \left(\tau t\right)}{\partial x}\text{d}{\omega }_D}+{T^{\prime }}_v=0$ (67)

It is assumed that the thickness t is a function of the variable s, hence (tds) in Equation (64) can be considered as dA. At the limit of integration s = L this becomes A and we get:

${\displaystyle {\int }_A\frac{\partial \sigma }{\partial x}\text{d}A}+P_x+q_L+q_K=0$ (68)

Integrating by parts the Equations (65)-(67):

${\left|\frac{\partial \left(\tau t\right)}{\partial x}y\right|}_K^L-{\displaystyle {\int }_{L}y\frac{\partial }{\partial S}\left[\frac{\partial \left(\tau t\right)}{\partial x}\right]\text{d}s}+P_y=0$ (69)

${\left|\frac{\partial \left(\tau t\right)}{\partial x}z\right|}_K^L-{\displaystyle {\int }_{L}z\frac{\partial }{\partial S}\left[\frac{\partial \left(\tau t\right)}{\partial x}\right]\text{d}s}+P_z=0$ (70)

${\left|\frac{\partial \left(\tau t\right)}{\partial x}{\omega }_D\right|}_K^L-{\displaystyle {\int }_L{\omega }_D\frac{\partial }{\partial S}\left[\frac{\partial \left(\tau t\right)}{\partial x}\right]\text{d}s}+{T^{\prime }}_v=0$ (71)

The values of the shear flow ( $q=\tau t$ ) per unit length in the longitudinal direction at the extreme points K and L of the element in Figure 14, are, s = s_K and s = s_L. These must equal to the shear stress per unit length in the lateral direction. Hence, if follows:

q_L (longitudinal) = q_L (lateral) and q_K (longitudinal) = q_K (lateral)

Also, from Equations (69)-(71), with $q^{\prime }=\frac{\text{d}q}{\text{d}x}$ we write:

${\left|\frac{\partial \left(\tau t\right)}{\text{d}x}y\right|}_K^L=\left({q^{\prime }}_L{y}_L-{q^{\prime }}_K{y}_K\right)$ (72)

${\left|\frac{\partial \left(\tau t\right)}{\text{d}x}z\right|}_K^L=\left(q^{\prime }{z}_L-{q^{\prime }}_L{z}_K\right)$ (73)

${\left|\frac{\partial \left(\tau t\right)}{\text{d}x}{\omega }_D\right|}_K^L=\left({q^{\prime }}_L{\omega }_{DL}-{q^{\prime }}_K{\omega }_{DK}\right)$ (74)

Substituting Equations (72)-(74), into Equations (69)-(71) respectively and applying the partial derivative identity: $\frac{\partial }{\text{d}s}\left[\frac{\partial \left(\tau t\right)}{\text{d}x}\right]=\frac{\partial }{\text{d}x}\left[\frac{\partial \left(\tau t\right)}{\text{d}S}\right]$, leads to:

${\displaystyle {\int }_{L}y\frac{\partial }{\partial x}\left[\frac{\partial \left(\tau t\right)}{\partial s}\right]\text{d}s}-P_y-\left({q^{\prime }}_L{y}_L-{q^{\prime }}_K{y}_K\right)=0$ (75)

${\displaystyle {\int }_{L}z\frac{\partial }{\partial x}\left[\frac{\partial \left(\tau t\right)}{\partial s}\right]\text{d}s}-P_z-\left({q^{\prime }}_L{z}_L-{q^{\prime }}_K{z}_K\right)=0$ (76)

${\displaystyle {\int }_L{\omega }_D\frac{\partial }{\partial x}\left[\frac{\partial \left(\tau t\right)}{\partial s}\right]\text{d}s}-{T^{\prime }}_v-\left({q^{\prime }}_L{\omega }_{DL}-{q^{\prime }}_K{\omega }_{DK}\right)=0$ (77)

Now recall Equations (53), and (37) and solve with Equations (68), (75), (76) and (77) to give the quartet:

$E{u^{″}}_o{\displaystyle {\int }_A\text{d}A}-E{v}^{‴}{\displaystyle {\int }_{A}y\text{d}A}-E{w}^{‴}{\displaystyle {\int }_{A}z\text{d}A}-E{\phi }^{‴}{\displaystyle {\int }_A{\omega }_D\text{d}A}+\left(q_L-q_K\right)+P_x=0$ (78)

$\begin{array}{l}{\displaystyle {\int }_L\frac{\partial P_x}{\partial x}y\text{d}s}+E{u^{‴}}_o{\displaystyle {\int }_{A}y\text{d}A}-Ev\text{'}\text{'}\text{'}\text{'}{\displaystyle {\int }_A{y}^2\text{d}A}-Ew\text{'}\text{'}\text{'}\text{'}{\displaystyle {\int }_{A}yz\text{d}A}\\ -E\phi \text{'}\text{'}\text{'}\text{'}{\displaystyle {\int }_{A}y{\omega }_D\text{d}A}+P_y+\left({q^{\prime }}_L{y}_L-{q^{\prime }}_K{y}_K\right)=0\end{array}$ (79)

$\begin{array}{l}{\displaystyle {\int }_L\frac{\partial P_x}{\partial x}z\text{d}s}+E{u^{‴}}_o{\displaystyle {\int }_{A}z\text{d}A}-Ev\text{'}\text{'}\text{'}\text{'}{\displaystyle {\int }_{A}yz\text{d}A}-Ew\text{'}\text{'}\text{'}\text{'}{\displaystyle {\int }_A{z}^2\text{d}A}\\ -E\phi \text{'}\text{'}\text{'}\text{'}{\displaystyle {\int }_{A}z{\omega }_D\text{d}A}+P_z+\left({q^{\prime }}_L{z}_L-{q^{\prime }}_K{z}_K\right)=0\end{array}$ (80)

$\begin{array}{l}{\displaystyle {\int }_L\frac{\partial P_x}{\partial x}{\omega }_D\text{d}s}+E{u^{‴}}_o{\displaystyle {\int }_A{\omega }_D\text{d}A}-Ev\text{'}\text{'}\text{'}\text{'}{\displaystyle {\int }_{A}y{\omega }_D\text{d}A}-Ew\text{'}\text{'}\text{'}\text{'}{\displaystyle {\int }_{A}z{\omega }_D\text{d}A}\\ -E\phi \text{'}\text{'}\text{'}\text{'}{\displaystyle {\int }_A{\omega }_D^2\text{d}A}+{T^{\prime }}_v+\left({q^{\prime }}_L{\omega }_{DL}-{q^{\prime }}_K{\omega }_{DK}\right)=0\end{array}$ (81)

For consistency of the last four equations, derive Equation (78), and substitute each integrand’s symbol within Equations (78)-(81) we get:

$E{u^{‴}}_{o}A-Ev\text{'}\text{'}\text{'}\text{'}{S}_z-Ew\text{'}\text{'}\text{'}\text{'}{S}_y-E\phi \text{'}\text{'}\text{'}\text{'}{S}_{\omega D}+\left({q^{\prime }}_L-{q^{\prime }}_K\right)+{P^{\prime }}_x=0$ (82)

${\displaystyle {\int }_L\frac{\partial P_x}{\partial x}y\text{d}s}+E{u^{‴}}_o{S}_z-Ev\text{'}\text{'}\text{'}\text{'}{I}_z-Ew\text{'}\text{'}\text{'}\text{'}{I}_{yz}-E\phi \text{'}\text{'}\text{'}\text{'}{S}_{Dz}+P_y+\left({q^{\prime }}_L{y}_L-{q^{\prime }}_K{y}_K\right)=0$ (83)

${\displaystyle {\int }_L\frac{\partial P_x}{\partial x}z\text{d}s}+E{u^{‴}}_o{S}_y-Ev\text{'}\text{'}\text{'}\text{'}{I}_{yz}-Ew\text{'}\text{'}\text{'}\text{'}{I}_y-E\phi \text{'}\text{'}\text{'}\text{'}{S}_{Dy}+P_z+\left({q^{\prime }}_L{z}_L-{q^{\prime }}_K{z}_K\right)=0$ (84)

$\begin{array}{l}{\displaystyle {\int }_L\frac{\partial P_x}{\partial x}{\omega }_D\text{d}s}+E{u^{‴}}_o{S}_{\omega D}-Ev\text{'}\text{'}\text{'}\text{'}{S}_{\omega Dz}-Ew\text{'}\text{'}\text{'}\text{'}{S}_{\omega Dy}-E\phi \text{'}\text{'}\text{'}\text{'}{\Gamma }_D\\ +{T^{\prime }}_v+\left({q^{\prime }}_L{\omega }_{DL}-{q^{\prime }}_K{\omega }_{DK}\right)=0\end{array}$ (85)

By assuming that, the origin of the system of axes coincides with the centroid of the cross section then,

$S_z={\displaystyle {\int }_{A}y\text{d}A}=0$ and $S_y={\displaystyle {\int }_{A}z\text{d}A}=0$

Applying this condition to Equations (82)-(85), provides the following solution:

$\begin{array}{l}\frac{S_{\omega D}}{A}\left(E\phi \text{'}\text{'}\text{'}\text{'}{S}_{\omega D}-C_1\right)+E{S}_{\omega Dy}\left(\frac{S_{\omega Dy}{I}_y-S_{\omega Dz}{I}_{yz}}{I_y{I}_z-I_{yz}^2}\right)\phi \text{'}\text{'}\text{'}\text{'}+S_{\omega Dy}\left(\frac{C_3{I}_{yz}-C_2{I}_y}{I_y{I}_z-I_{yz}^2}\right)\\ +E{S}_{\omega Dz}\left(\frac{S_{\omega Dz}{I}_z-S_{\omega Dy}{I}_{yz}}{I_y{I}_z-I_{yz}^2}\right)\phi \text{'}\text{'}\text{'}\text{'}+S_{\omega Dz}\left(\frac{C_2{I}_{yz}–C_3{I}_z}{I_y{I}_z-I_{yz}^2}\right)-E{\Gamma }_D\phi \text{'}\text{'}\text{'}\text{'}+C_4=0\end{array}$ (86)

where:

$C_1=\left({q^{\prime }}_L-{q^{\prime }}_K\right)+{P^{\prime }}_x$

$C_2=\left({q^{\prime }}_L{y}_L-{q^{\prime }}_K{y}_K\right)+P_y+{\displaystyle {\int }_L\frac{\partial P_x}{\partial x}y\text{d}s}$

$C_3=\left({q^{\prime }}_L{z}_L-{q^{\prime }}_K{z}_K\right)+P_z+{\displaystyle {\int }_L\frac{\partial P_x}{\partial x}z\text{d}s}$

$C_4=\left({q^{\prime }}_L{\omega }_{DL}-{q^{\prime }}_K{\omega }_{DK}\right)+{T^{\prime }}_v+{\displaystyle {\int }_L\frac{\partial P_x}{\partial x}{\omega }_D\text{d}s}$

Re-writing Equation (86), we get:

$\begin{array}{l}\left\{\frac{S_{\omega D}^2}{A}+S_{\omega Dy}\left(\frac{S_{\omega Dy}{I}_y-S_{\omega Dz}{I}_{yz}}{I_y{I}_z-I_{yz}^2}\right)+S_{\omega Dz}\left(\frac{S_{\omega Dz}{I}_z-S_{\omega Dy}{I}_{yz}}{I_y{I}_z-I_{yz}^2}\right)-{\Gamma }_D\right\}E\phi \text{'}\text{'}\text{'}\text{'}\\ -C_2\left(\frac{S_{\omega Dy}{I}_y-S_{\omega Dz}{I}_{yz}}{I_y{I}_z-I_{yz}^2}\right)-C_3\left(\frac{S_{\omega Dz}{I}_z-S_{\omega Dy}{I}_{yz}}{I_y{I}_z-I_{yz}^2}\right)-C_1\left(\frac{S_{\omega D}}{A}\right)+C_4=0\end{array}$ (87)

Substituting Equations (32), and (33) into Equation (87), leads to:

$\begin{array}{l}\left\{\frac{S_{\omega D}^2}{A}+\left(a_z-d_z\right)S_{\omega Dy}-\left(a_y-d_y\right)S_{\omega Dz}-{\Gamma }_D\right\}E\phi \text{'}\text{'}\text{'}\text{'}\\ -\left(a_z-d_z\right)C_2+\left(a_y-d_y\right)C_3-C_1\left(\frac{S_{\omega D}}{A}\right)+C_4=0\end{array}$ (88)

Substitute the expressions for, C₁, C₂, C₃, and C₄, into Equation (88), and applying Equation (34):

$\begin{array}{l}-E{\Gamma }_A\phi \text{'}\text{'}\text{'}\text{'}-\left(a_z-d_z\right)\left[\left({q^{\prime }}_L{y}_L-{q^{\prime }}_K{y}_K\right)+P_y+{\displaystyle {\int }_L\frac{\partial P_x}{\partial x}y\text{d}s}\right]\\ +\left(a_y+d_y\right)\left[\left({q^{\prime }}_L{z}_L-{q^{\prime }}_K{z}_K\right)+P_z+{\displaystyle {\int }_L\frac{\partial p_x}{\partial x}z\text{d}s}\right]\\ -\frac{S_{\omega D}}{A}\left({q^{\prime }}_L–{q^{\prime }}_K+{p^{\prime }}_x\right)+\left({q^{\prime }}_L{\omega }_{DL}-{q^{\prime }}_K{\omega }_{DK}\right)+{\displaystyle {\int }_L\frac{\partial P_x}{\partial x}{\omega }_D\text{d}s}+{T^{\prime }}_v=0\end{array}$ (89)

Equation (89) can be rewritten as follows,

$\begin{array}{l}-E{\Gamma }_A\phi \text{'}\text{'}\text{'}\text{'}+{q^{\prime }}_L\left[{\omega }_{DL}-\left(a_z-d_z\right)y_L+\left(a_y-d_y\right)z_L–\left(\frac{S_{\omega D}}{A}\right)\right]\\ -{q^{\prime }}_K\left[{\omega }_{DK}-\left(a_z-d_z\right)y_K+\left(a_y-d_y\right)z_K–\left(\frac{S_{\omega D}}{A}\right)\right]\\ -\left[\left(a_z-d_z\right)p_y+\left(a_y-d_y\right)p_z\right]+{p^{\prime }}_x\left[{\omega }_D-\left(a_z-d_z\right)y+\left(a_y-d_y\right)z-\left(\frac{S_{\omega D}}{A}\right)\right]\\ -{\displaystyle {\int }_L\frac{\partial P_x}{\partial x}\left[\left(a_z-d_z\right)y-\left(a_y-d_y\right)z+{\omega }_D\right]\text{d}s}+{T^{\prime }}_v=0\end{array}$ (90)

Now substituting Equation (31), into Equation (90), we get:

$\begin{array}{l}-E{\Gamma }_A\phi \text{'}\text{'}\text{'}\text{'}+\left({q^{\prime }}_L{\omega }_{AL}-{q^{\prime }}_K{\omega }_{AK}\right)+\left[\left(a_y-d_y\right)p_z-\left(a_z-d_z\right)p_y\right]+{p^{\prime }}_x{\omega }_A\\ -{\displaystyle {\int }_L\frac{\partial p_x}{\partial x}\left[\left(a_z-d_z\right)y-\left(a_y-d_y\right)z+{\omega }_D\right]\text{d}s}+{T^{\prime }}_v=0\end{array}$ (91)

By considering p_x is only a function of x, Equation (93), reads:

$-E{\Gamma }_A\phi \text{'}\text{'}\text{'}\text{'}+\left({q^{\prime }}_L{\omega }_{AL}-{q^{\prime }}_K{\omega }_{AK}\right)+\left[\left(a_y-d_y\right)p_z-\left(a_z-d_z\right)p_y\right]+{p^{\prime }}_x{\omega }_A+{T^{\prime }}_v=0$ (92)

In Equation (92), we need to define the last term ( ${T^{\prime }}_v$ ) which represents the differential torsional moment on the finite strip of width dx, from the conventional theory of pure torsion, which is mainly related to St-Venant’s hypothesis that, this moment will only result from the torsional shearing stresses, then St-Venant’s torque will be proportional to the rate of change of twisting angle w. r. t. the longitudinal dimension x, i.e. with ${\phi }^{\prime }=\text{d}\phi /\text{d}x$, this leads to the following:

$T_v=GJ{\phi }^{\prime }$ (93)

and:

${T^{\prime }}_v=GJ{\phi }^{″}$ (94)

where, G is the modulus of rigidity, and J is the polar moment of area [29]. Substituting Equation (94) into (92), to read as follows,

$E{\Gamma }_A\phi \text{'}\text{'}\text{'}\text{'}-GJ{\phi }^{″}=\left({q^{\prime }}_L{\omega }_{AL}-{q^{\prime }}_K{\omega }_{AK}\right)+\left[\left(a_y-d_y\right)p_z-\left(a_z-d_z\right)p_y\right]+{p^{\prime }}_x{\omega }_A$ (95)

Subscript letter A is used to indicate that A is the point of the principal sectorial pole.

Equation (95) is the final general form of the differential equation concerning the torsional moment effect of the thin-walled beam. For the terms p_x,p_y, p_z,q_L, andq_K refer to Figure 14. The right-hand side of this equation, reveals three different sources of torsional moments:

1) The first one is the contribution of shearing stresses at sides L and K. This term may vanish, when for thin-walled beams, both points lie at the extreme sides of the open cross-section, or, when this equation is applied to a closed cross section.

2) The second term is due to the external transverse loads p_y, and p_z acting on the beam.

3) The third term is due to the rate of change of the axial force p_x, with respect to the longitudinal direction x, exerted on the beam at an arbitrary point on the cross section, apart from the principal sectorial pole A.

Equations (78), (83), (84), and (85), re-appear as follows:

$-EA{u^{″}}_o=P_x$ (96)

$E{I}_{z}v\text{'}\text{'}\text{'}\text{'}+E{I}_{yz}w\text{'}\text{'}\text{'}\text{'}=P_y$ (97)

$E{I}_{yz}v\text{'}\text{'}\text{'}\text{'}+E{I}_{y}w\text{'}\text{'}\text{'}\text{'}=P_z$ (98)

$E{\Gamma }_A\phi \text{'}\text{'}\text{'}\text{'}-GJ{\phi }^{″}=m$ (99)

Equations (96)-(98) were already known from the conventional strength of materials theory developed earlier by Oden and Ripperger [28] and Timoshenko and Gere [29] [30].

As was mentioned above the additional differential Equation (99), was found originally by Vlasov [8]. The significance of this present derivation is that it began by considering an arbitrary pole point as a centre of rotation. Consequently, this led to the fact that the centre of rotation of an open thin-walled section is unique, identified as the principal sectorial pole. Also, found here is that, in derivation of the equation for shear stress, there appeared a bi-moment integral term B in Equations (46) and (47). B transpired through mathematical processes applied within Equations (41)-(49). It was not imposed as a term required to satisfy any of the three elements of elasticity theory, namely: equilibrium, compatibility and boundary conditions [29] [30].

8. Derivation of Torsion Warping Element Stiffness Matrix

The theory of torsional warping couples the twist mode, defined by φ, with the warping mode defined by dφ/dx (see Figure 15).

Thus, the stiffness matrix for a number of nodes i and j, may be defined by:

${\left[T_i,B_i,T_j,B_j\right]}^{\text{T}}=\left[K_w\right]{\left[{\phi }_i,{{\phi }^{\prime }}_i,{\phi }_j,{{\phi }^{\prime }}_j\right]}^{\text{T}}$ (100)

Figure 15. Nodal torsional forces and displacements of thin–walled beam element.

From Equation (99) we have:

$\phi \text{'}\text{'}\text{'}\text{'}-K^2{\phi }^{″}=\frac{m}{E\Gamma }$ (101)

where $K^2=\frac{GJ}{E\Gamma }$.

The solution to the homogeneous differential Equation (101) for m = 0, is:

$\phi =A_1\sinh \left(kx\right)+A_2\cosh \left(kx\right)+A_{3}x+A_4$ (102)

where A₁, A₂, A₃, and A₄ are arbitrary constants, dependent on the boundary conditions.

Now by differentiating $\phi$ with respect to x three times, we have:

${\phi }^{\prime }=A_{1}k\cosh \left(kx\right)+A_{2}k\sinh \left(kx\right)+A_3$ (103)

${\phi }^{″}=A_1{k}^2\sinh \left(kx\right)+A_2{k}^2\cosh \left(kx\right)$ (104)

${\phi }^{‴}=A_1{k}^3\cosh \left(kx\right)+A_2{k}^3\sinh \left(kx\right)$ (105)

Recalling Equations (47) and from (93) for the bimoment and torsional moment, Figure 15 shows:

$B=-E\Gamma {\phi }^{″}$ (106)

$T=T_v+T_w$ or,

$T=T_v+B^{\prime }$

$T=GJ{\phi }^{\prime }+E\Gamma {\phi }^{‴}$ (107)

Substituting Equations (103)-(105) into Equations (106) and (107) leads to:

$T=A_{3}E\Gamma K^2$ (108)

$B=-E\Gamma K^2+A_{3}E\Gamma K^{2}x+A_{4}E\Gamma K^2$ (109)

Using Figure 15, and Equations (106) and (107) we define the constants A₁, A₂, A₃, and A₄ in terms of nodal displacements, ${\phi }_i$, ${{\phi }^{\prime }}_i$, ${\phi }_j$, and ${{\phi }^{\prime }}_j$ as follow:

$\begin{array}{l}{A}_1=\frac{1}{K{D}_o}[-\left(K\sinh \alpha \right){\phi }_i+\left(\cosh \alpha -1-\alpha \sinh \alpha \right){{\phi }^{\prime }}_i\\ +\left(K\sinh \alpha \right){\phi }_j-\left(\cosh \alpha -1\right){{\phi }^{\prime }}_j]\end{array}$ (110)

$A_2=\frac{1}{K{D}_o}\left[K\left(\cosh \alpha -1\right){\phi }_i+n{{\phi }^{\prime }}_i+K\left(\cosh \alpha -1\right){\phi }_j+\left(\sinh \alpha -\alpha \right){{\phi }^{\prime }}_j\right]$ (111)

$A_3=\frac{1}{K{D}_o}\left[\left(K\sinh \alpha \right){\phi }_i+\left(\cosh \alpha -1\right){{\phi }^{\prime }}_i-\left(K\sinh \alpha \right){\phi }_j+\left(\cosh \alpha -1\right){{\phi }^{\prime }}_j\right]$ (112)

$\begin{array}{l}{A}_4=\frac{1}{K{D}_o}[K\left(\cosh \alpha -1-\alpha \sinh \alpha \right){\phi }_i-n{{\phi }^{\prime }}_i+K\left(\cosh \alpha -1\right){\phi }_j\\ -\left(\sinh \alpha -\alpha \right){{\phi }^{\prime }}_j]\end{array}$ (113)

where, $D_o=2\cosh \alpha -2-\alpha \sinh \alpha$, $n=\cosh \alpha -\sinh \alpha$ and $\alpha =k\mathcal{l}$. Substituting Equations (110) to (113) into Equations (108) and (109) provides the torque T and the bi-moment B:

$T=\frac{E\Gamma }{D_o}{K}^2\left[\left(K\sinh \alpha \right){\phi }_i+\left(\cosh \alpha -1\right){{\phi }^{\prime }}_i-\left(K\sinh \alpha \right){\phi }_j+\left(\cosh \alpha -1\right){{\phi }^{\prime }}_j\right]$ (114)

$\begin{array}{c}B=-E\Gamma K^2\phi +\frac{E\Gamma }{D_o}{K}^2[\left\{\left(\cosh \alpha -1\right)-\left(1-\frac{x}{\mathcal{l}}\alpha \sinh \alpha \right)\right\}{\phi }_j\\ +\left\{x\left(\cosh \alpha -1\right)-\frac{\alpha \left(\cosh \alpha -\sinh \alpha \right)}{K}\right\}{{\phi }^{\prime }}_i\\ +\left\{\left(\cosh \alpha -1\right)-xK\sinh \alpha \right\}{\phi }_j+\left\{x\left(\cosh \alpha -1\right)-\frac{\left(\sinh \alpha -\alpha \right)}{K}\right\}{{\phi }^{\prime }}_j]\end{array}$ (115)

(115)

Equations (114) and (115) represent general equations of the torsional moments and bi-moments along the beam element from $x=0$ to $x=\mathcal{l}$, the sign convention will now be used in order that a positive twisting angle ${\phi }_i$ will be associated with a positive torque (twisting moment) T_i, when the remaining three degrees of freedom: ${{\phi }^{\prime }}_i$, ${\phi }_j$ and ${{\phi }^{\prime }}_j$, each have a value of zero. Conversely, a positive twisting angle ${\phi }_j$ will be associated, conventionally, with a positive torque T_j again where ${\phi }_i$, ${\phi }_j$ and ${{\phi }^{\prime }}_j$ are zeros.

The same rule will be used for the warping mode and the bi-moment, i.e. a positive warping mode ${\phi }_i$, is associated with a positive bi-moment B_i where ${\phi }_i$, ${\phi }_j$ and ${{\phi }^{\prime }}_j$, are zeros Similarly for ${\phi }_j$ and B_j.

Applying this rule to Equations (114) and (115), bearing in mind that the other signs in the equations follow suit, they appear in the matrix form given in Table 1:

One must mention that this element stiffness matrix, is the same as was produced by Krahula [10], and Renton [11]. Also, the above matrix, represents only terms concerning the torsional moment and bi-moment. To include the stiffness terms concerning the six degrees of freedom in the conventional strength of materials theory, we may call readily on the stiffness matrix, produced by many authors,

Table 1. Torsional and bi-moment stiffness matrix.

including Prezemieniecki [31] [32] and Zienkiewicz, [33]. See Table 2, as example of the standard 12 × 12 beam element stiffness matrix [34].

The elements in rows/columns 1 and 7 of this matrix arise from an axial tension applied to each bar element node. Those in 3 and 9 arise from nodal bending and shear. The fourth and tenth rows and columns represent the torsional stiffness, (GJ/L) which refers to a torsional moment that twists the beam around its longitudinal axis. The corresponding force and displacement vectors for Table 2 are:

$f^e={\left[f_x{}_{\text{1}}{q}_y{}_{\text{1}}{q}_z{}_{\text{1}}{t}_x{}_{\text{1}}{m}_y{}_{\text{1}}{m}_z{}_{\text{1}}{f}_x{}_{\text{2}}{q}_y{}_{\text{2}}{q}_z{}_{\text{2}}{t}_x{}_{\text{2}}{m}_y{}_{\text{2}}{m}_z{}_{\text{2}}\right]}^{\text{T}}$

${\delta }^e={\left[u_1{v}_1{w}_1{\phi }_1{\theta }_{y1}{\theta }_{z1}{u}_{x2}{v}_{z2}{w}_y{}_2{\phi }_2{\theta }_{y2}{\theta }_{z2}\right]}^{\text{T}}$

where the beam bending rotations are: θ_y = dw/dx and θ_z= dv/dx. Within the element’s force vector $f^e$, for each section node 1 and 2: f_x is the axial tensile force aligned with the x-direction, (replacing N above) and t_x is an axial torque (replacing T) applied about the X-axis, see Figure 16.

In addition, q_y andq_z are transverse shear forces aligned with section co-ordinates y and z, and m_y and m_z are bending moments applied about the y- and z-axes at each of the element’s ends. The sense of actions f_x, q_y andq_z is positive for positive co-ordinates x, y and z and so too t_x, m_y and m_z are positive by the right.

Table 2. Beam element’s conventional stiffness matrix.

Figure 16. Nodal forces for a thin–walled beam element under combined loading.

In addition, q_y andq_z are transverse shear forces aligned with section co-ordinates y and z, and m_y and m_z are bending moments applied about the y- and z-axes at each of the element’s ends. The sense of actions f_x, q_y andq_z is positive for positive co-ordinates x, y and z and so too t_x, m_y and m_z are positive by the right-hand screw rule. The displacement vector ${\delta }^e$ identifies deformation with the four loadings as follows: displacement u with f_x, v and w under q_y and q_z respectively, angular twist φ under t and rotations θ_y and θ_y under m_z and m_y respectively. Both φ andθ are given in radian.

It is necessary to modify Table 2 to include the Wagner torque and bi-moment derived here. The torsional stiffness term GJ/L is only valid for thick-walled and solid beams. It is no longer valid for thin-walled beam sections where it should be replaced by Equation (114). The hyperbolic terms appear within the rows (and columns) four and eleven in Table 3. Also, Equation (115) for the bi-moment has been added to the stiffness matrix creating the seventh and fourteenth rows (and columns). The elements shown represent the resistance to the warping of cross-section from the bi-moment. The full 14 × 14 matrix, so modified for a thin-walled open section beam, is given as Table 3. The accompanying element’s force $f^e$ and displacement ${\delta }^e$ vectors for the thin-walled beam’s stiffness matrices $K^e$, ( $f^e=K^e{\delta }^e$ ) given in Table 3, are as follows:

$f^e={\left[f_x{}_{\text{1}}{q}_y{}_{\text{1}}{q}_z{}_{\text{1}}{t}_x{}_{\text{1}}{m}_y{}_{\text{1}}{m}_z{}_{\text{1}}{b}_1{f}_x{}_{\text{2}}{q}_y{}_{\text{2}}{q}_z{}_{\text{2}}{t}_x{}_{\text{2}}{m}_y{}_{\text{2}}{m}_z{}_{\text{2}}{b}_2\right]}^{\text{T}}$

${\delta }^e={\left[u_1{v}_1{w}_1{\phi }_1{\theta }_{y1}{\theta }_{z1}{{\phi }^{\prime }}_1{u}_{x2}{v}_{z2}{w}_y{}_2{\phi }_2{\theta }_{y2}{\theta }_{z2}{{\phi }^{\prime }}_2\right]}^{\text{T}}$

For consistency lower case t_x replaces T and b replaces B in Equations (114) and (115). It is pertinent to recognize how these additional terms arise. Figure 16 shows the origin of co-ordinates x, y and z at O. If O is the area centroid of the cross section and Oz is an axis of symmetry then the element will receive tension, bending shear and torsion under the six nodal forces acting at each of its

Table 3. General stiffness matrix for thin-walled beam element under combined loading (torsion, bending and shear).

nodes 1 and 2, as given in Figure 16. Here the force q_z passes through the shear centre E since this will lie on Oz but q_y does not. Relocating q_y to pass through E requires an accompanying torsion t_x about the flexural axis E, which is also the axis of the centre of twist, that point in the cross section which does not rotate. Equivalently, t_x may be returned to the axis Ox as shown given that here Oz and Oy are principal axes for the section’s second moments of area originating at the section’s centroid.

The force vector $f^e$ above also applies generally where O does not lie at the centroid. Here it is necessary to identify moments m_y and m_z about the area centroid principal axes and to refer shear forces q_z and q_y to the shear centre E with an accompanying “Wagner torsion” and bi-moment referred to axis O. In this way the force and displacement vectors may be referred to a common set of co-ordinates O, x, y and z as indicated within the elements of the matrix given in Table 3. This has been the case considered generally above. Reductions apply only when the origin of co-ordinates lies at the shear centre E in Figure 16 for then the beam will bend and shear but not twist. Taking the centre of rotation for cross-section to coincide with the shear centre is the basis of Wagner torsional theory of thin sections [1]. It is noted that some [35] have disputed this claiming that the former centre lies at a position which ensures a minimum in the store of potential energy.

9. Conclusions

The theory of thin-walled beams of open sections as reviewed suggests certain refinements and adequate definitions to accompany the theory. These were made in order to clarify the ambiguities that appeared in the original as it was. Thus, based on Vlasov’s theory, and through mathematical derivation, a general differential equation of beam element as well as its general solution was derived. The above work created the possibility to combine the axial, and transverse forces, to include the effect of extension/compression, bending and torsion on the warping of the cross-section. Here a stiffness matrix for the element has been derived. The work proved the following points:

Taking the centre of rotation for cross-section to coincide with the shear centre is the basis of Wagner’s torsional theory of thin sections. Both the angle of twist of the beam $\phi$, and the warping of its cross-section ${\phi }^{\prime }$ are functions of three terms: 1) the shear flow multiplied by the value of sectorial area at the point; 2) the torsional moment about the shear centre; and 3) the axial load multiplied by the sectorial area at the point of its application.

Further work upon that given in [15] is required to derive a general transformation matrix corresponding to the loading considered for the stiffness matrix given in Table 3.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Guo, Z.-H. (1981) A Unified Theory of Thin-Walled Structures. Journal of Structural Mechanics, 9, 179-197. https://doi.org/10.1080/03601218108907382
[2]	Landiero, F. and Savoia, M. (1990) Shear Strain Effects in Flexure and Torsion of Thin-Walled Beams with Open or Closed Cross Section. Thin Walled Structures, 10, 87-119. https://doi.org/10.1016/0263-8231(90)90058-7
[3]	Fatmi, R.E. (2007) Non-Uniform Warping Including the Effects of Torsion and Shear Forces, Part I: A General Beam Theory. The International Journal of Solids and Structures, 44, 5912-5929. https://doi.org/10.1016/j.ijsolstr.2007.02.006
[4]	Wang, Z.R., Zhao, J.C. and Gong, J.H. (2011) A New Torsion Element of Thin-Walled Beams Including Shear Deformation. Applied Mechanics and Materials, 94-96, 1642-1645. https://doi.org/10.4028/www.scientific.net/AMM.94-96.1642
[5]	Vieira, R.F., Virtuoso, F.B. and Pereria, E.B.R. (2014) A Higher Order Model for Thin-Walled Structures with Deformable Cross Sections. The International Journal of Solids and Structures, 51, 575-598. https://doi.org/10.1016/j.ijsolstr.2013.10.023
[6]	Pavazza, R., Matokovic, A. and Vukasovic, M. (2020) A Theory of Torsion of Thin-Walled Beams of Arbitrary Open-Section with Influence of Shear. Mechanics Based Design of Structures and Machines, 1-36.
[7]	Wagner, H. (1936) Torsion and Buckling of Open Sections. T. M. No. 807, US National Advisory Committee for Aeronautics.
[8]	Vlasov, V.Z. (1961) Thin-Walled Elastic Beams. Oldbourne Press, London, Chapter 1.
[9]	Zbirohowski-ko’scia, K. (1967) Thin-Walled Beams from Theory to Practice. Crosby Lockwood & Son Ltd., London.
[10]	Krahula, J.L. (1967) Analysis of Bent and Twisted Bars Using the Finite Element Method. AIAA Journal, 5, 1194-1197. https://doi.org/10.2514/3.4163
[11]	Renton, J.D. (1960) A Direct Solution of the Torsional-Flexural Buckling of Axially Loaded Thin-Walled Bars. The Institution of Structural Engineers, London.
[12]	Krajcinovic, D. (1969) A Consistent Discrete Elements Technique for Thin-Walled Assemblages. International Journal of Solids and Structures, 5, 639-662. https://doi.org/10.1016/0020-7683(69)90085-7
[13]	Rajasekaran, S. (1977) Finite Element Method for Plastic Beam-Columns. In: Chem, W.-F. and Atsuta, T., Eds., Theory of Beam-Columns, Vol. 2, Space Behaviour and Design, McGraw-Hill International Book Company, New York, 539-608.
[14]	Baigent, A.H. and Hancock, G.J. (1982) Structural Analysis of Assemblages of Thin-Walled Members. Engineering Structures, 4, 207-216. https://doi.org/10.1016/0141-0296(82)90010-4
[15]	Al-Sheikh, A. (1985) Behaviour of Thin-Walled Structures under Combined Loads. Ph.D. Thesis, Loughborough University of Technology, Loughborough.
[16]	Chen, B. and Hu, Y. (1988) The Torsional Stiffness Matrix of a Thin-Walled Beam and Its Application to Beams under Combined Loading. Computers and Structures, 28, 421-431. https://doi.org/10.1016/0045-7949(88)90080-6
[17]	Kollbrunner, C.F. and Hajdin, N. (1966) Wolbkraft-torsion dunnwandiger Stabe mit geschlossenem profil., Schweizer Stahlbauverband Mitteiungen der Technischen Komm, 32.
[18]	Dvorkin, E.N., Celentano, D., Cuitino, A. and Gioia, G. (1989) A Vlassov Beam Element. Computers & Structures, 33, 187-196. https://doi.org/10.1016/0045-7949(89)90140-5
[19]	Sorin, D.M. and Bogdan, I.E. (1999) Study of Warping Torsion of Thin-Walled Beams with Open Cross-Section Using Macro-Elements. International Journal for Numerical Methods in Engineering, 44, 853-868. https://doi.org/10.1002/(SICI)1097-0207(19990228)44:6<853::AID-NME534>3.0.CO;2-P
[20]	Wen, Y.L. and Kuo, M.H. (2003) More General Expression for the Torsional Warping of a Thin-Walled Open-Section Beam. International Journal of Mechanical Sciences, 45, 831-849. https://doi.org/10.1016/S0020-7403(03)00137-1
[21]	Erkmen, R.E. (2006) Finite Element Formulations for Thin-Walled Members. Ph.D. Thesis, Department of Civil Engineering University of Ottawa, Ontario.
[22]	Vetyuko, Y. (2008) Direct Approach to Elastic Deformations and Stability of Thin-Walled Rods of Open Profile. Springer-Verlag, Berlin.
[23]	Wang, et al. (2012) Restrained Torsion of Open Thin-Walled Beams Including Shear Deformation Effects. Journal of Zhejiang University—Science A: Applied Physics & Engineering, 13, 260-273. https://doi.org/10.1631/jzus.A1100149
[24]	Khalil, M.F. and Cespedes, X. (2014) A New Beam Element with Transversal and Warping Eigenmodes. Computers and Structures, 131, 12-33. https://doi.org/10.1016/j.compstruc.2013.10.001
[25]	Choi, I.S., Jang, G.W., Choi, S., Shin, D. and Kim, Y.Y. (2017) Higher Order Analysis of Thin-Walled Beams with Axially Varying Quadrilateral Cross Sections. Computers and Structures, 179, 127-137. https://doi.org/10.1016/j.compstruc.2016.10.025
[26]	Drillenburg, M. (2017) Torsion and Shear Stresses Due to Shear Centre Eccentricity. SCIA Engineer Delft University of Technology.
[27]	Duan, L. and Zhao, J. (2019) A Geometrically Exact Cross-Section Deformable Thin-Walled Beam Finite Element Based on Generalized Beam Theory. Computers & Structures, 218, 32-59. https://doi.org/10.1016/j.compstruc.2019.04.001
[28]	Oden, J.T. and Ripperger, E.A. (1980) Mechanics of Elastic Structures. Second Edition, McGraw-Hill Book Company, New York.
[29]	Timoshenko, S.P. (1956) Strength of Materials. Third Edition, Part II, D. Van Nostrand Co. Inc., Princeton.
[30]	Timoshenko, S.P. and Gere, J.M. (1961) Theory of Elastic Stability. Second Edition, McGraw-Hill Book Co. Inc., New York.
[31]	Przemieniecki, J.S. (1972) Matrix Analysis of Local Instability in Plates, Stiffened Panels and Columns. International Journal for Numerical method in Engineering, 5, 151-300. https://doi.org/10.1002/nme.1620050207
[32]	Przemieniecki, J.S. (1973) Finite Element Structural Analysis of Local Instability. AIAA Journal, 11, 33-39. https://doi.org/10.2514/3.50433
[33]	Zienkiewicz, O.C. (1971) The Finite Element Method in Engineering Science. McGraw-Hill, London.
[34]	Gjelsvik, A. (1981) The Theory of Thin-Walled Beams. John Wiley and Sons, Hoboken.
[35]	Bleich, F. and Bleich, H.H. (1952) Buckling Strength of Metal Structures. McGraw-Hill Book Company, Inc., New York.