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

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 unfo-reseen 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 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.


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

List of Symbols
A is the centre of twist for the cross section. u is the longitudinal displacement of an arbitrary point S, of Cartesian coordinate (y, z) and sectorial coordinate ω . 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 w ′ ′ are rates of change of displacement of point S, (with respect to X-axis), in the y and z-directions respectively. ϕ′ is the rate of change of twist of the beam with respect to the X-axis.

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

Out of Plane and Lateral Displacements
As was stated by Vlasov [8] (pp. 7-9), the assumption of a non-deformable crosssection 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 zand 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., 0 γ = (see Figure 3).

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: If ϕ is small then, cos 1 ϕ ≈ , and sinϕ ϕ ≈ . Hence,   Rewriting Equations (9) and (10): Recalling Equations (5) and (6) and substituting into Equations (11) and (12) we get: 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 1 γ and 2 γ are the total change of angle between both axes x and s, 1 . Their sum describes the engineering shear strain in radians: Given that the total shear angle γ is very small, and can be neglected, then 0 γ = . This assumption is quite feasible physically since it means that,  Now apply the expression for η from Equation (13) and find its rate of change w r t, the variable x: where prime designates ( x ∂ ∂ ), noting that from Figure It should be obvious that, v, w, and ϕ 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 x , resulting in following equation: Here ( ) f x , the integration constant, is given the symbol, 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 ω .  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 η of point S w. r. t. x in the yand z-directions respectively. Also, ϕ′ is the rate of change of twist of the beam with respect to the X-axis. Hence the product ϕ ω ′ 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 ϕ ω ′ .
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. ω ω = = , at s = 0 (see Figure 9). Now from Figure 9 twice the shaded area is given by the following equation:   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):

The Law of Sectorial Area The Determination of Sectorial Area Expressions
Integrating Equations (21) and (22) we get: Since C 1 and C 2 are the integration constants, they are found from imposing the condition: which is consistent with the conditions mentioned before: 0 Similarly, for D ω : Subtracting Equation (26) from Equation (25) we get: For the calculation of A ω and 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: 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) If the origin of the coordinate system OYZ coincides with the centroid of the cross-section we have: and Equation (29) reads: Substituting Equation (30) into (27) we have: Also, multiplying Equation (27) by ydA and zdA respectively and integrating over the whole area of the cross-section A, we have: and: Now imposing the integrals in Equations (28)

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: where σ is the normal stress, ε is the strain, and E is Young's modulus. Now recall Equation (20) The normal stresses (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: is the total bending moment about Z-axis, A is the total area of the cross section, and t is the wall thickness.
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 σ from Equation (36) for which the sectorial coordinate ω was chosen arbitrarily. Here Equation (44) applies to the normal stress in a beam section that rotates about an arbitrary S ω appear. This approach is different to that given by Vlasov [8]. Specifically, he took 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) (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.

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 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.: 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 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 ( ) xs V τ . Hence the total shear stress (Figure 11

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: Dividing Equation (52) by dx, and integrating equation along the profile S, we get: where: ( ) x ′ = For a useful but less general expression We note that in Equations (55)  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 τ is to be determined. Also, S ω 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).
where T ω from Equation (58), is defined by the following:

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: 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). World Journal of Mechanics The first element in Equations (61) and (62) are the projections of the differential of the shear force upon the yand z-axes, where, P y = P y (x) and P z = P z (x) are external of transverse loads in the yand 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: 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: The values of the shear flow ( q 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) Substituting Equations (72)-(74), into Equations (69)-(71) respectively and applying the partial derivative identity: Now recall Equations (53), and (37) and solve with Equations (68), (75), (76) and (77) to give the quartet: For consistency of the last four equations, derive Equation (78), and substitute each integrand's symbol within Equations (78)-(81) we get: By assuming that, the origin of the system of axes coincides with the centroid of the cross section then, where:  I I I  I I I   S I S I  S I S I  S  C  C  C  C  A  I I I  I I Substitute the expressions for, C 1 , C 2 , C 3 , and C 4 , into Equation (88), and applying Equation (34): Equation (89) can be rewritten as follows, Now substituting Equation (31), into Equation (90), we get: By considering p x is only a function of x, Equation (93), reads: In Equation (92), we need to define the last term ( v T ′ ) 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.
where, G is the modulus of rigidity, and J is the polar moment of area [29]. Subs- 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 , and q 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.
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].

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:  From Equation (99) we have: The solution to the homogeneous differential Equation (101) for m = 0, is: where A 1 , A 2 , A 3 , and A 4 are arbitrary constants, dependent on the boundary conditions. Now by differentiating ϕ with respect to x three times, we have: Recalling Equations (47) and from (93) for the bimoment and torsional moment, Figure 15 shows: Substituting Equations (103)-(105) into Equations (106) and (107) leads to: Using Figure 15, and Equations (106) and (107) we define the constants A 1 , A 2 , A 3 , and A 4 in terms of nodal displacements, i ϕ , i ϕ′ , j ϕ , and j ϕ′ as follow: Equations (114) and (115) represent general equations of the torsional moments and bi-moments along the beam element from 0 x = to x =  , the sign convention will now be used in order that a positive twisting angle i ϕ will be associated with a positive torque (twisting moment) T i , when the remaining three degrees of freedom: i ϕ′ , j ϕ and j ϕ′ , each have a value of zero. Conversely, a positive twisting angle j ϕ will be associated, conventionally, with a positive torque T j again where i ϕ , j ϕ and j ϕ′ are zeros.
The same rule will be used for the warping mode and the bi-moment, i.e. a positive warping mode i ϕ , is associated with a positive bi-moment B i where i ϕ , j ϕ and j ϕ′ , are zeros Similarly for 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,   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: where the beam bending rotations are: θ y = dw/dx and θ z = dv/dx. Within the element's force vector e f , 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 and q 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 yand z-axes at each of the element's ends. The sense of actions f x , q y and q 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 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 e f 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.

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 ϕ , and the warping of its cross-section ϕ′ 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.