An efficient modelling approach based on a rigorous cross-sectional analysis for analysing box girder bridge superstructures

Abstract

A novel analysis technique is introduced for efficient modelling of box girder bridge decks. The general three-dimensional equations used to accurately define the deformation of these complex beam-like slender structures are decoupled into a two-dimensional cross-sectional problem and a one-dimensional beam problem through decomposition of the three-dimensional strain field. The two-dimensional cross-sectional problem is solved by a two-dimensional finite element analysis considering in-plane as well as out-of-plane warping displacements of the beam section. This gives an accurate constitutive relationship of the one-dimensional beam problem without making any major assumptions as is often done in usual beam theories. The one-dimensional beam problem is solved by a one-dimensional beam finite analysis and the results obtained are used to recover three-dimensional stress, strain and displacement fields accurately. Numerical examples of box girder bridge deck systems having thin-walled sections are solved by the proposed approach to show its performance.

Keywords

box girder bridge deck finite element analysis warping

Introduction

Beam-like slender structures have long been used in civil engineering. A common application of such structures is for bridges with box girder deck systems which are basically thin-walled beams having closed or a combination of closed and open sections. The behaviour of these thin-walled structures under an arbitrary loading scenario is quite complex, primarily due to the cross-sectional warping (out-of-plane warping) and the distortion of the sections (in-plane warping). In recent years, accurate calculation of warping displacements has been the subject of many research studies since the variation of warping displacements over a cross-section does not follow a standard pattern and it changes from one case to another depending on the cross-sectional profile and the loading pattern.

The establishment and development of finite element methods (FEMs) have made it possible to more accurately predict the behaviour of these structures, but this approach has some disadvantages. For example, it is very time-consuming and requires significant computational resources (i.e. computational and storage memory), especially for elaborate geometries. In addition, this much effort is excessive in the initial design stages. Moreover, all comprehensive FEM software packages are very expensive compared to simpler beam-modelling tools. The alternative is to model these structures as a beam which makes the analysis highly efficient computationally, but a model based on a classical beam theory such as Euler–Bernoulli hypothesis cannot predict the behaviour of these structures well because it assumes that a beam cross-section will remain plane and perpendicular to beam axis after deformation and so cannot capture the local deformations in the form of warping and the effect of shear deformation. The shear deformation effects are incorporated in Timoshenko’s beam theory, but this theory also cannot capture the warping deformation. Saint Venant tried to incorporate the effect of torsional deformation which is absent in the above theories, but eventually he could only consider the effect of uniform (unconstrained) warping (out-of-plane) without any shear deformation. A better representation of the deformation of thin-walled beams is made in Vlasov’s (1961) theory which addresses the problem of constrained warping (out-of-plane), but not the in-plane warping (distortion) of the beam section. A number of researchers (e.g. Bauld and Tzeng, 1984; Chandra et al., 1990; Lee, 2005; Sheikh and Thomsen, 2008) have used Vlasov’s theory to make some advancements, but these approaches need prior knowledge of the pattern of warping displacements which can change from one case to another. Subsequently, the effect of distortion in thin-walled box girders has been studied (e.g. Kermani and Waldron, 1993; Razaqpur and Li, 1991, 1994), but the techniques they developed for modelling the distortion are based on ad hoc assumptions and cannot be used in general cases.

In this context, the concept of beam sectional analysis proposed by Giavotto et al. (1983) and later extended and employed by Borri and Merlini (1986), Ghiringhelli and Mantegazza (1994) and few others (e.g. Blasques, 2014) seems to be the most attractive. In this approach, the complete three-dimensional (3D) elasticity problem defining the actual behaviour of the beam-like structures is decomposed to a two-dimensional (2D) beam section analysis and a one-dimensional (1D) beam analysis without the need for any ad hoc assumptions. The method only assumes that the cross-sectional dimensions are small compared to the beam length which is true for slender beam structures. The 2D beam sectional analysis is carried out using a 2D finite element (FE) discretisation where the effects of in-plane warping as well as out-of-plane warping are considered. The 2D FE analysis generates the exact constitutive matrix (or stiffness matrix) of the beam cross-section which ensures a proper coupling between the different modes of deformation. This cross-section stiffness matrix is then used in the 1D beam analysis which can be carried out using a standard 1D beam FE model. The stress resultants obtained in the 1D beam analysis together with the results of the 2D cross-sectional analysis are used in combination to obtain the warping displacements and finally recover the 3D stress and displacement fields of the beam. The computational efficiency of this approach is remarkable in terms of the prediction of the 3D response of these structures.

A similar approach employed by Hodges et al. (1992), called the variation asymptotic method (VAM), is based on a rigorous mathematical foundation. The VAM was introduced by Berdichevskii (1979) who applied this method to shell structures which helped to reduce the 3D problem into a 2D problem consistently utilising the shell thickness as the small parameter which is a key concept of this method. Cesnik and Hodges (1997) extended this method for beam analysis and developed variational asymptotic beam section (VABS) analysis which is appealing due to its mathematical consistency, but the method is significantly complex with respect to its mathematical treatment. On the other hand, the approach proposed by Giavotto et al. (1983) is less complex but with similar capabilities as that of VABS.

It is interesting to note that these methods so far have been only applied in the analysis of aerospace-related structures such as helicopter, wind turbine blades and so on. To the knowledge of the authors, no one has attempted to take advantage of such methods for the analysis of box girder bridge decks and similar structures.

In this article, an attempt has been made to develop Giavotto’s technique for the analysis of box girder bridge decks and similar structures having closed or a combination of open and closed sections. The 2D cross-sectional problem is solved using eight-node quadratic isoparametric elements, whereas three-node isoparametric linear elements are used to solve the 1D beam problem. A computer code was developed in FORTRAN to implement the different steps associated with the 2D sectional analysis, the 1D beam analysis as well as the recovery of the beam 3D response. In order to test the performance of the proposed analysis technique, a number of numerical examples including solid and thin-walled box girders with different cross-sectional configurations have been solved of which some are reported here. Detailed 3D FE analyses of these beams have been also carried out using the commercially available FE code, ABAQUS (which has been recognised to be capable of giving accurate predictions for box girder response), and the results obtained are used to validate the results produced by the proposed technique.

Mathematical formulation

The formulation is based on the assumption that the beam-like structure has a prismatic slender geometry and does not have any abrupt variation of the cross-sectional geometry and material properties. Similarly, the loading should be such that it will only produce a small or gradual variation of deformations and strains along the beam length. The formulation is based on small deformation theory considering elastic material behaviour.

2D sectional analysis

Figure 1 shows a typical beam-like structure where the Cartesian coordinate axes (x-y-z) are used as the reference system. The displacement vector at an arbitrary point of a beam section ${u}$ can be divided into two components as

{\begin{matrix} u_{x} \\ u_{y} \\ u_{z} \end{matrix}} = {\begin{matrix} v_{x} \\ v_{y} \\ v_{z} \end{matrix}} + {\begin{matrix} w_{x} \\ w_{y} \\ w_{z} \end{matrix}} or {u} = {v} + {w}

(1)

where ${v}$ is the displacement vector due to translations and rotations of the section without having any distortions and ${w}$ is due to the generalised 3D warping of the beam section. The first component of the displacement vector ${v}$ can be expressed in terms of translations ${U} = [\begin{matrix} U_{x} & U_{y} & U_{z} \end{matrix}]^{T}$ and rotations ${θ} = [\begin{matrix} θ_{x} & θ_{y} & θ_{z} \end{matrix}]^{T}$ of the section as

{v} = [\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix} \begin{matrix} 0 & 0 & - y \\ 0 & 0 & x \\ y & - x & 0 \end{matrix}] {\begin{matrix} {U} \\ {θ} \end{matrix}} = [H] {Δ}

(2)

Figure 1.

A typical beam-like structure with reference coordinate system.

For the warping displacements, the beam section is discretised with 2D FEs which help to express these displacements as

{w} = [N] {d}

(3)

where ${d}$ contains the nodal parameters of $w_{x}$ , $w_{y}$ and $w_{z}$ and $[N]$ is their interpolation function (Bathe, 1982).

The 3D stress–strain relationship is written as

{σ} = [D] {ε}

(4)

where the stress and strain vectors are arranged as ${σ} = [\begin{matrix} σ_{x} & σ_{y} & τ_{xy} \end{matrix} \begin{matrix} τ_{xz} & τ_{yz} & σ_{z} \end{matrix}]^{T}$ and ${ε} = [\begin{matrix} ε_{x} & ε_{y} & γ_{xy} \end{matrix} \begin{matrix} γ_{xz} & γ_{yz} & ε_{z} \end{matrix}]^{T}$ , respectively, and $[D]$ is the corresponding constitutive matrix.

The stress resultants ( ${P} = [\begin{matrix} {F}^{T} & {M}^{T} \end{matrix}]^{T}$ ) at a section of the beam consist of three forces ( ${F} = [\begin{matrix} F_{x} & F_{y} & F_{z} \end{matrix}]^{T}$ ) and three moments ( ${M} = [\begin{matrix} M_{x} & M_{y} & M_{z} \end{matrix}]^{T}$ ) which are produced by the last three components of the stress vector ${τ} = [\begin{matrix} τ_{xz} & τ_{yz} & σ_{z} \end{matrix}]^{T}$ treated as tractions acting on the beam section having an area of A. The stress resultants can be obtained from

{P} = \int [H]^{T} {τ} dA

(5)

The strain vector is now expressed in terms of derivatives of the displacements where the derivatives with respect to x and y (sectional coordinates) are separated from the derivatives with respect to z (axial coordinate) as

\begin{matrix} {ε} = & {\begin{matrix} ε_{x} \\ ε_{y} \\ γ_{xy} \\ γ_{xz} \\ γ_{yz} \\ ε_{z} \end{matrix}} = [\begin{matrix} \partial / \partial x & 0 & 0 \\ 0 & \partial / \partial y & 0 \\ \partial / \partial y & \partial / \partial x & 0 \\ 0 & 0 & \partial / \partial x \\ 0 & 0 & \partial / \partial y \\ 0 & 0 & 0 \end{matrix}] {\begin{matrix} u_{x} \\ u_{y} \\ u_{z} \end{matrix}} \\ + [\begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}] \frac{\partial}{\partial z} {\begin{matrix} u_{x} \\ u_{y} \\ u_{z} \end{matrix}} = [\partial_{xy}] {u} + [S] {u}_{, z} \end{matrix}

(6)

Equations (1)–(3) are substituted in the above equation which leads to

\begin{matrix} {ε} = & [\partial_{xy}] [H] {Δ} + [S] [H] {Δ}_{, z} \\ + [\partial_{xy}] [N] {d} + [S] [N] {d}_{, z} \end{matrix}

(7)

As $[\partial_{xy}] [H]$ in the above equation can be written as [S][H][T] where [T] is a 6×6 square matrix having all elements zero except T(1, 5) = −1 and T(2, 4) = 1, the strain vector can be rewritten as

\begin{matrix} {ε} = & [S] [H] [T] {Δ} + [S] [H] {Δ}_{, z} \\ + [B] {d} + [S] [N] {d}_{, z} \end{matrix}

(8)

The first two terms of the above equation are rearranged to express the strain vector in its final form as

{ε} = [S] [H] {Ψ} + [B] {d} + [S] [N] {d}_{, z}

(9)

where ${Ψ}$ is the sectional strain vector consists of two shear strains, one axial strain, two curvatures and one torsional strain. The sectional strain vector can simply be expressed as

{Ψ} = [T] {Δ} + {Δ}_{, z}

(10)

Equilibrium equations

Now the virtual work principle ( $δ W_{i} = δ W_{e}$ ) can be applied to get the equilibrium equations of the beam section. For this purpose, a small element of the beam in the form of a small slice having a length of dz is taken. Assuming no variation of stresses and strains within the small length dz (along z), the internal virtual work can be written as

δ W_{i} = dz \int δ {ε}^{T} {σ} dA

(11)

The external virtual work will be caused by the stresses or tractions ${τ}$ acting on the two cross-sections of the slice and $δ W_{e}$ can be written as

\begin{matrix} δ W_{e} & = {[\int δ {u}^{T} {τ} dA + \frac{\partial}{\partial z} (\int δ {u}^{T} {τ} dA) dz]}_{at (z + dz)} \\ + {[- \int δ {u}^{T} {τ} dA]}_{at z} \\ = \frac{\partial}{\partial z} (\int δ {u}^{T} {τ} dA) dz \\ = dz (\int δ {u}_{, z}^{T} {τ} dA + \int δ {u}^{T} {τ}_{, z} dA) \end{matrix}

(12)

So the final form of the virtual work can be written as

\int δ {u}_{, z}^{T} {τ} dA + \int δ {u}^{T} {τ}_{, z} dA = \int δ {ε}^{T} {σ} dA

(13)

After substitution of equations (1)–(5) and equations (9) and (10) in the above equation, it can be rearranged as

\begin{matrix} {\begin{matrix} δ {d}_{, z} \\ δ {d} \\ δ {Ψ} \end{matrix}}^{T} [\begin{matrix} [A] & [C] & [E] \\ {[C]}^{T} & [L] & [M] \\ {[E]}^{T} & {[M]}^{T} & [R] \end{matrix}] {\begin{matrix} {d}_{, z} \\ {d} \\ {Ψ} \end{matrix}} \\ = {\begin{matrix} δ {d}_{, z} \\ δ {d} \\ δ {Ψ} \end{matrix}}^{T} {\begin{matrix} {Q} \\ {Q}_{, z} \\ {P} \end{matrix}} + δ {Δ}^{T} ({P}_{, z} - [T]^{T} {P}) \end{matrix}

(14)

where $[A] = \int {[N]}^{T} {[S]}^{T} [D] [S] [N] dA$ , $[C] = \int {[N]}^{T} {[S]}^{T} [D] [B] dA$ , $[E] = \int {[N]}^{T} {[S]}^{T} [D] [S] [H] dA$ , $[L] = \int {[B]}^{T} [D] [B] dA$ , $[M] = \int {[B]}^{T} [D] [S] [H] dA$ , $[R] = \int {[H]}^{T} {[S]}^{T} [D] [S] [H] dA$ , ${Q} = \int {[N]}^{T} {τ} dA$ and ${Q}_{, z} = \int {[N]}^{T} {τ}_{, z} dA$ . It should be noted that the above integrations are carried out for the individual elements and the resulting matrices and vectors for all the elements modelling the beam section are assembled together to form equation (14). For simplicity, the notations used for all quantities at the element level are the same as those at the section level.

For any arbitrary value of $δ {d}$ , $δ {d}_{, z}$ , $δ {Ψ}$ and $δ {Δ}$ , equations (14) is valid if

[A] {d}_{, z} + [C] {d} + [E] {Ψ} = {Q}

(15)

[C]^{T} {d}_{, z} + [L] {d} + [M] {Ψ} = {Q}_{, z}

(16)

[E]^{T} {d}_{, z} + [M]^{T} {d} + [R] {Ψ} = {P}

(17)

{P}_{, z} = [T]^{T} {P}

(18)

Equations (15) and (16) can be replaced with a single equation by differentiation of equation (15) with respect to z and this will lead to the final form of the equilibrium equations for the beam section as

\begin{matrix} [A] {d}_{, zz} + ([C] - [C]^{T}) {d}_{, z} \\ + [E] {Ψ}_{, z} - [L] {d} - [M] {Ψ} = 0 \end{matrix}

(19)

[E]^{T} {d}_{, z} + [M]^{T} {d} + [R] {Ψ} = {P}

(20)

Solution of the equilibrium equations

Equations (19) and (20) have two solutions where the complementary solution corresponds to ${P} = 0$ (i.e. self-balanced condition), which lead to an Eigen problem. The Eigen values can be used to estimate the dissipation length of end effects which are usually small (in the order of the largest dimension of the beam section). From this physical perspective, Giavotto et al. (1983) defined this as the extremity solution. For the particular solution, which is defined as the central solution by Giavotto et al. (1983), equations (19) and (20) along with their derivatives with respect to z can be rewritten as

\begin{matrix} [L] {d} & + [M] {Ψ} = [A] {d}_{, zz} \\ + ([C] - [C]^{T}) {d}_{, z} + [E] {Ψ}_{, z} \end{matrix}

(21)

[M]^{T} {d} + [R] {Ψ} = - [E]^{T} {d}_{, z} + {P}

(22)

\begin{matrix} [L] {d}_{, z} & + [M] {Ψ}_{, z} = [A] {d}_{, zzz} \\ + ([C] - [C]^{T}) {d}_{, zz} + [E] {Ψ}_{, zz} \end{matrix}

(23)

[M]^{T} {d}_{, z} + [R] {Ψ}_{, z} = - [E]^{T} {d}_{, zz} + {P}_{, z}

(24)

A further derivative of equations (23) and (24) will eliminate ${P}_{, zz}$ according to the equilibrium equations of the beam which state that the resulting cross-sectional forces vary linearly along the beam. This leads to the solutions of ${d}$ and ${Ψ}$ in the form of linear combinations of polynomials in z which may have a highest order of n. If equations (23) and (24) are differentiated further up to (n − 1) times, the right-hand side of the last two equations will be zero since $(\partial^{(n + 1)} {d}) / (\partial z^{(n + 1)}) = 0$ and $(\partial^{(n + 1)} {Ψ}) / (\partial z^{(n + 1)}) = 0$ as the highest order of ${d}$ and ${Ψ}$ is n. So the solution of the left-hand side quantities $(\partial^{(n)} {d}) / (\partial z^{(n)})$ and $(\partial^{(n)} {Ψ}) / (\partial z^{(n)})$ of these two equations will be zero which will again make the right-hand sides of the previous two equations zero. If this process is continued successively, the second and subsequent higher-order derivatives of ${d}$ and ${Ψ}$ will be found to be zero. With this information, equations (21) and (24) can be simplified as

[L] {d} + [M] {Ψ} = ([C] - [C]^{T}) {d}_{, z} + [E] {Ψ}_{, z}

(25)

[M]^{T} {d} + [R] {Ψ} = - [E]^{T} {d}_{, z} + {P}

(26)

[L] {d}_{, z} + [M] {Ψ}_{, z} = 0

(27)

[M]^{T} {d}_{, z} + [R] {Ψ}_{, z} = {P}_{, z}

(28)

For a given value of the stress resultant vector ( ${P}$ ), ${P}_{, z}$ can simply be calculated using equation (18) and with this known ${P}_{, z}$ , equations (27) and (28) can be solved simultaneously to evaluate ${d}_{, z}$ and ${Ψ}_{, z}$ . These will make the right-hand side of equations (25) and (26) known which can finally be solved to get ${d}$ and ${Ψ}$ .

Constraints

The above equations cannot be solved without some constraints being imposed as the displacement formulation defined in equation (1) is six times indeterminate. This is due to the use of ${v}$ (equation (2)) in addition to ${w}$ (equation (3)) to define the displacements of a point within the beam section. Actually, ${w}$ could be used to define the total displacement vector ${u}$ through the 2D FE approximation, but ${w}$ is used only for the warping displacements and ${v}$ is used in addition to define ${u}$ . In the present case, the number of constraints is six as ${v}$ is defined in terms of ${Δ}$ (equation (2)) which has six components. The constraints can be defined on the basis that the warping displacements should not contribute towards any translation or rotation of the beam section and can be expressed as

\begin{matrix} \sum_{i = 1}^{n} d_{3 i - 2} = 0, \sum_{i = 1}^{n} d_{3 i - 1} = 0, \sum_{i = 1}^{n} d_{3 i} = 0, \\ \sum_{i = 1}^{n} - z_{i} d_{3 i - 1} + y_{i} d_{3 i} = 0, \sum_{i = 1}^{n} z_{i} d_{3 i - 2} - x_{i} d_{3 i} = 0, \\ \sum_{i = 1}^{n} - y_{i} d_{3 i - 2} + x_{i} d_{3 i - 1} = 0 \end{matrix}

(29)

where $d_{3 i - 2}$ , $d_{3 i - 1}$ and $d_{3 i}$ are nodal parameters (ith node) corresponding to $w_{x}$ , $w_{y}$ and $w_{z}$ , respectively, $(x_{i}, y_{i}, z_{i})$ is the nodal coordinate vector and n is the total number of nodes used to model the beam section. These six constraints for the nodal displacement parameters ${d}$ are also applicable to ${d}_{, z}$ and all these constraints can be expressed as

\begin{matrix} [c] {d} = {0}, [c] {d}_{, z} = {0} and \\ [c] = [\begin{matrix} [I_{3}] & [I_{3}] & \dots & [I_{3}] & \dots & [I_{3}] \\ [h_{1}] & [h_{1}] & \dots & [h_{i}] & \dots & [h_{n}] \end{matrix}] \end{matrix}

(30)

where

[h_{i}] = [\begin{matrix} 0 & 0 & y_{i} \\ 0 & 0 & - x_{i} \\ - y_{i} & x_{i} & 0 \end{matrix}]

and $[I_{3}]$ is the identity matrix of order 3.

The Lagrangian multiplier technique is used to impose the above constraints (equation (30)) with the equilibrium equations (equations (25)–(28)) as

\begin{matrix} [\begin{matrix} [L] & [M] & {[c]}^{T} \\ {[M]}^{T} & [R] & [0] \\ [c] & [0] & [0] \end{matrix}] {\begin{matrix} {d} \\ {ψ} \\ {λ_{1}} \end{matrix}} \\ = [\begin{matrix} ([C] - {[C]}^{T}) & [E] \\ - {[E]}^{T} & [0] \\ [0] & [0] \end{matrix}] {\begin{matrix} {d}_{, z} \\ {ψ}_{, z} \end{matrix}} + {\begin{matrix} {0} \\ {P} \\ {0} \end{matrix}} \end{matrix}

(31)

[\begin{matrix} [L] & [M] & {[c]}^{T} \\ {[M]}^{T} & [R] & [0] \\ [c] & [0] & [0] \end{matrix}] {\begin{matrix} {d}_{, z} \\ {ψ}_{, z} \\ {λ_{2}} \end{matrix}} = {\begin{matrix} {0} \\ {P}_{, z} \\ {0} \end{matrix}}

(32)

where ${λ_{1}}$ and ${λ_{2}}$ contain Lagrange multipliers.

Cross-section stiffness matrix

Equations (31) and (32) have solutions which are linear and homogeneous functions of ${P}$ and can be written as

\begin{matrix} {d} = [X_{1}] {P}, {d}_{, z} = [X_{2}] {P}, {ψ} = [Y_{1}] {P}, \\ {ψ}_{, z} = [Y_{2}] {P}, {λ_{1}} = [Λ_{1}] {P}, {λ_{2}} = [Λ_{2}] {P} \end{matrix}

(33)

where $[X_{1}]$ , $[X_{2}]$ , $[Y_{1}]$ , $[Y_{2}]$ , $[Λ_{1}]$ and $[Λ_{2}]$ are linear operators consisting of some constants. With the help of equations (18) and (33), equations (31) and (32) can be written as

\begin{matrix} [\begin{matrix} [L] & [M] & {[c]}^{T} \\ {[M]}^{T} & [R] & [0] \\ [c] & [0] & [0] \end{matrix}] {\begin{matrix} [X_{1}] \\ [Y_{1}] \\ [Λ_{1}] \end{matrix}} \\ = [\begin{matrix} ([C] - {[C]}^{T}) & [E] \\ - {[E]}^{T} & [0] \\ [0] & [0] \end{matrix}] {\begin{matrix} [X_{2}] \\ [Y_{2}] \end{matrix}} + {\begin{matrix} [0] \\ [I] \\ [0] \end{matrix}} \end{matrix}

(34)

[\begin{matrix} [L] & [M] & {[c]}^{T} \\ {[M]}^{T} & [R] & [0] \\ [c] & [0] & [0] \end{matrix}] {\begin{matrix} [X_{2}] \\ [Y_{2}] \\ [Λ_{2}] \end{matrix}} = {\begin{matrix} [0] \\ {[T]}^{T} \\ [0] \end{matrix}}

(35)

This is the final form of the governing equations which are solved to calculate the linear operators and these are utilised to calculate the sectional stiffness/constitutive matrix.

Using equations (11) and (12), the virtual work expression can now be written as

\frac{\partial}{\partial z} (\int δ {u}^{T} {τ} dA) = \int δ {ε}^{T} {σ} dA

(36)

The left-hand side of the above equation represents the external virtual work and may be restated in terms of the stress resultant vector ${P}$ and sectional flexibility matrix $[F_{s}]$ as

d {P}^{T} [F_{s}] {P} = \int δ {ε}^{T} {σ} dA

(37)

Considering equations (4), (7) and (33), the above equation can be expressed as

\begin{matrix} d {P}^{T} [F_{s}] {P} \\ = \int δ {P}^{T} ({[Y_{1}]}^{T} {[H]}^{T} {[S]}^{T} + {[X_{1}]}^{T} {[B]}^{T} + {[X_{2}]}^{T} {[N]}^{T} {[S]}^{T}) [D] \\ ([S] [H] [Y_{1}] + [B] [X_{1}] + [S] [N] [X_{2}]) {P} dA \end{matrix}

(38)

The integrations in the above equation are carried out similar to those in equation (14) and it can be expressed as

\begin{matrix} δ {P}^{T} [F_{s}] {P} = δ {P}^{T} \\ {[\begin{matrix} [X_{1}] \\ [X_{2}] \\ [Y_{1}] \end{matrix}]}^{T} [\begin{matrix} [L] & [C] & [M] \\ {[C]}^{T} & [A] & [E] \\ {[M]}^{T} & {[E]}^{T} & [R] \end{matrix}] [\begin{matrix} [X_{1}] \\ [X_{2}] \\ [Y_{1}] \end{matrix}] {P} \end{matrix}

(39)

so that the stiffness matrix of the beam section can be written as

\begin{matrix} [K_{s}] = [F_{s}]^{- 1} = \\ {({[\begin{matrix} [X_{1}] \\ [X_{2}] \\ [Y_{1}] \end{matrix}]}^{T} [\begin{matrix} [L] & [C] & [M] \\ {[C]}^{T} & [A] & [E] \\ {[M]}^{T} & {[E]}^{T} & [R] \end{matrix}] [\begin{matrix} [X_{1}] \\ [X_{2}] \\ [Y_{1}] \end{matrix}])}^{- 1} \end{matrix}

(40)

1D beam analysis and 3D stress recovery

The 1D beam analysis is carried out using beam FEs where a C⁰ continuous displacement-based formulation is adopted for the derivation of these elements considering bi-axial bending, bi-axial shear, torsion and axial deformations. Each element has three nodes which allow a quadratic variation of the six displacement components ${Δ}$ as defined in equation (2). The sectional stiffness matrix as obtained in equation (40) is used as the constitutive matrix for the derivation of the stiffness matrix of the beam elements where ${Ψ}$ (equation (10)) is used as the generalised strain vector of the beam. In order to avoid the shear locking problem and give better stress prediction without any oscillation, the field consistent technique (Prathap and Babu, 1986) has been used. In the 1D FE analysis, the nodal displacements are first obtained which are utilised to evaluate the generalised strain vector ${Ψ}$ and the stress resultant vector ${P} (= [K_{s}] {Ψ})$ at any desirable section of the beam. With the stress resultant vector ${P}$ , the 3D stresses can easily be recovered using equations (4), (9) and (33).

Numerical examples

A computer program was developed in FORTRAN to implement the proposed method and solve numerical examples of beams having different cross-sections (mostly thin-walled box girders commonly found in bridge decks) to demonstrate the performance of the method. The method has the capability of accommodating anisotropic and inhomogeneous materials, but the material properties used in all these examples are isotropic and homogeneous for simplicity. More specifically, the material is steel with an elastic modulus E = 200 GPa and Poisson ratio ν = 0.3. According to this technique, the cross-section analysis was first carried out to evaluate the sectional stiffness matrix which was then used in the 1D beam FE analyses, and the 3D displacement and stress field were finally recovered in all the cases. For the validation of the present results, these structures were also analysed with the ABAQUS software, where solid elements were used to create a detailed 3D FE model of these beams. The minimum number of elements (or degree of freedom (DOF)) necessary for the convergence of the results is reported for each example and load case.

The following sections give examples of the significant reductions in computational effort that can be achieved in terms of the reduced DOFs using the proposed method of analysis over those associated with detailed 3D FE modelling.

A cantilever beam having a solid square section

The beam has a length of 20 meters (m) while the breadth and depth of its cross-section is 1.0 m. In order to implement the sectional analysis, a 2D FE analysis using eight-node quadratic elements is carried out. The results obtained for the section stiffness parameters corresponding to three different mesh sizes are presented in Table 1, which shows a good and monotonic convergence of all the stiffness parameters. For the validation of the present results, the stiffness parameters used in a conventional beam analysis based on the classical beam theory are presented in Table 1. The torsional constant (J) in the classical beam theory is obtained from the analytical solutions based on membrane analogy (Timoshenko and Goodier, 1970). The convergence of the present results is achieved by 196 elements (i.e. 14×14 element mesh in Table 1) and they agreed well with the classical beam theory with a maximum difference of 0.2%.

Table 1.

Section stiffness parameters of the solid section (1×1 m).

Stiffness component	K₁₁ = K₂₂ (N)	K₃₃ (N)	K₄₄ = K₅₅ (N m²)	K₆₆ (N m²)
Present^a (5×5)	6.5988×10¹⁰	200.00×10⁹	16.805×10⁹	11.204×10⁹
Present^a (10×10)	6.4378×10¹⁰	200.00×10⁹	16.726×10⁹	10.9298×10⁹
Present^a (14×14)	6.4158×10¹⁰	200.00×10⁹	16.699×10⁹	10.836×10⁹
Classical beam theory	(κGA)^b 6.4102×10¹⁰	(EA) 200.00×10⁹	(EI) 16.667×10⁹	(GJ) 10.815×10⁹
Diff. (%)	+0.087	+0.0	+0.2	+0.18

Mesh size (number of element in x-direction×number of element in y-direction).

(κ) shear correction factor (Gere and Timoshenko, 1990).

With the converged values of the sectional stiffness matrix corresponding to a mesh of 14×14 elements (No of DOF = 1935), the 1D beam analysis was carried out for two load cases: (1) a shear force $F_{y} = - 1.0 kN$ applied at the free end and (2) a torsional moment $M_{z} = - 1.0 kN m$ applied at the free end. The 1D FE analysis was carried out using 15 three-node beam elements. The displacement components at the free end of beam are presented in Tables 2 and 3 and the variations of stress components at the free end and mid-span sections are plotted in Figures 2 to 6.

Table 2.

Displacement components (mm) at the free end for example 1 (Load case 1).

Displacement component	$u_{y}$	$u_{z}$ (x = 0.5 m, y = 0.5 m)	$u_{z}$ (x = 0, y = 0.5 m)
Present	1.5983×10⁻¹	5.9826×10⁻³	5.9831×10⁻³
Classical beam theory	1.6031×10⁻¹	5.9928×10⁻³	5.9928×10⁻³
ABAQUS	1.5975×10⁻¹	5.9875×10⁻³	5.9883×10⁻³

Table 3.

The angle of twist at the free end for example 1 (Load case 2).

Component	Classical beam theory	Present	ABAQUS
Angle of twist	1.8492×10⁻⁶	1.8439×10⁻⁶	1.8464×10⁻⁶

Figure 2.

Normal stress (σ_zz kN/m²) at the mid-span section (Load case 1): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 3.

Transverse shear stress (σ_yz kN/m²) at the mid-span section (Load case 1): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 4.

Transverse shear stress (σ_xz kN/m²) at the mid-span section (Load case 1): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 5.

Transverse shear stress (σ_yz kN/m²) at the free end section (Load case 2): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 6.

Transverse shear stress (σ_xz kN/m²) at the free end section (Load case 2): (a) proposed method and (b) 3D FE model (ABAQUS).

For the validation of the present results, the beam was analysed with ABAQUS using 20-node quadratic hexahedral solid elements (C3D20) with a mesh size of 14×14×140 (No of DOF = 377,775) with an aspect ratio of 1×1×2 for an element. For Load case 1, the shear force was applied as a uniform surface traction at the free end while the twisting moment for Load case 2 was applied at the centroid of the free end using a feature in ABAQUS known as multiple point constraint (MPC). Timoshenko’s beam theory using the shear correction factor of $κ = 10 (1 + ν) / (12 + 11 ν)$ was used to evaluate the transverse displacement $u_{y}$ of the beam at its free end for Load case 1. Similarly, the axial displacement $u_{z}$ at the top fibre of the free end is evaluated using a 2D analytical solution of Timoshenko and Goodier (1970). The angle of twist at the free end for Load case 2 is also calculated using an analytical solution based on membrane analogy (Timoshenko and Goodier, 1970). The results obtained from these analytical solutions as well as 3D FE analyses using ABAQUS are presented along with those obtained by the proposed method similarly in Tables 2 and 3 and Figures 2 to 6. Even with significantly reduced DOF, the results were found in all the cases to be in excellent agreement with the 3D FE and analytical solutions. As can be seen from Figures 3 and 4, there is a variation of transverse shear stresses which cannot be captured using the classical beam theory. It is worth mentioning that in all examples stress results were obtained for each node and reported as they are and no stress smoothing was conducted as for commercial packages like ABAQUS.

Thin-walled single cell box beam

A thin-walled closed section box beam as shown in Figure 7 was analysed for three load cases where the first two load cases are identical to those of the previous example and the third load case consists of a force $F_{y} = - 1.0 kN$ acting at the middle of one of the web plates (right side) as shown in Figure 7. The 2D cross-sectional analysis was carried out with 196 eight-node quadratic elements having a size of 0.02×0.02 m (DOF = 2946) while the 1D analysis is carried out using 15 beam elements (DOF = 186) which ensured the convergence of results. The beam was also analysed with ABAQUS using 11,700 numbers of solid elements (C3D20) where each element has an aspect ratio of 1.25:1:1 and a total number of DOF of 148,980. Note that the proposed method has roughly 50 times fewer DOFs than the ABAQUS model. For Load case 1, the vertical force $F_{y} = - 1.0 kN$ (total) was applied as a uniform surface traction on two web plates while the twisting moment $M_{z} = - 1.0 kN m$ in Load case 2 was represented in the form of two equal and opposite vertical forces which applied as uniform traction on two web plates at the free end of the box beam. The vertical force $F_{y} = - 1.0 kN$ in Load case 3 was also applied as a uniform traction, but it was only applied on the right side web plate at the free end (Figure 7). The displacement components at different locations of the beam obtained from the proposed method and the 3D FE analysis using ABAQUS are presented in Tables 4 and 5 for the three load cases. Table 4 shows that there is a variation of the longitudinal displacement (u_z) along the flange which is reflected by the displacement values at two critical points, points A and B (Figure 7). This variation is caused by the warping of the thin-walled section which cannot be obtained via the conventional beam models. The results obtained by the present technique have a maximum difference of 2.1% when compared with the 3D FE results. Table 5 compares the performance of the method in the vicinity of the load and away from the load for the eccentrically loaded beam (Load case 3). The maximum difference between the results obtained for both methods is 1.8% at the free end while this difference is only 0.28% at the mid-span.

Figure 7.

Thin-walled closed section box beam (dimensions are in m).

Table 4.

Displacement components at the free end of the beam (Load cases: 1 and 2).

	Load case 1			Load case 2
Component	$u_{z} (mm)$ at A (Figure 7)	$u_{z} (mm)$ at B (Figure 7)	$u_{y} (mm)$ at A and B	$θ_{z}$
Present	3.9739×10⁻²	3.9686×10⁻²	1.0667	6.8156×10⁻⁶
ABAQUS	4.0611×10⁻²	3.9741×10⁻²	1.0695	6.6522×10⁻⁶

Table 5.

Deflection $u_{y} (mm)$ at free end and mid-span section of the beam (Load case 3).

Location (Figure 7)	Mid-span section		Free end section
Location (Figure 7)	A	B	A	B
Present	3.3649×10⁻¹	3.3497×10⁻¹	1.0701	1.0667
ABAQUS	3.3556×10⁻¹	3.3555×10⁻¹	1.0911	1.0695

Similarly, the variations of stress components obtained from the proposed method and the 3D FE analysis at different sections of the beam are plotted in Figures 8 to 15. In these cases, it can be seen that the agreement is very good.

Figure 8.

Normal stress (σ_zz kN/m²) at the mid-span section of the beam (Figure 7) (Load case 1): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 9.

Transverse shear stress (σ_yz kN/m²) at the mid-span section of the beam (Figure 7) (Load case 1): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 10.

Transverse shear stress (σ_xz kN/m²) at the mid-span section of the beam (Figure 7) (Load case 1): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 11.

Transverse shear stress (σ_xz kN/m²) at a section 5 m away from the free end of the beam (Figure 7) (Load case 2): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 12.

Transverse shear stress (σ_yz kN/m²) at a section 5 m away from the free end of the beam (Figure 7) (Load case 2): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 13.

Normal stress (σ_zz kN/m²) at the mid-span section of the beam (Figure 7) (Load case 3): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 14.

Transverse shear stress (σ_yz kN/m²) at the mid-span section of the beam (Figure 7) (Load case 3): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 15.

Transverse shear stress (σ_xz kN/m²) at the mid-span section of the beam (Figure 7) (Load case 3): (a) proposed method and (b) 3D FE model (ABAQUS).

Thin-walled two-cell box beam

A thin-walled two-cell box girder as shown in Figure 16 is solved for two load cases as those of the previous example.

Figure 16.

Thin-walled two-cell rectangular box beam (dimensions are in m).

For the 2D cross-sectional analysis, a single 2D element was used along the thickness direction of all flange and web plates while the number of elements used along the other direction of these plates was 60 and 25 for each flange and web, respectively. The 1D beam analysis was carried out using 15 beam elements. This led to a total DOF (DOFs) of 3012 and 186 for the 2D sectional analysis and 1D beam analysis, respectively. On the other hand, the number of DOF required for the detailed 3D FE model of the beam was 351,051 (more than 100 times the DOFs for the proposed method). For the 3D FE analysis (ABAQUS), the shear force ( $F_{y} = - 1.0 kN$ ) was applied as a uniform traction on the three web plates at the free end of the beam, whereas in the other load case, the twisting moment ( $M_{z} = - 1.0 kN m$ ) was applied at the centroid of the free end with the use of MPC (ABAQUS). The displacements at key points of the beam’s free end obtained by the proposed approach and the 3D FE analysis are presented in Tables 6 and 7 for the two load cases. The results are in very good agreement with a maximum difference of only 0.3%.

Table 6.

Displacement components (mm) at the free end (Load case 1).

Method of analysis	$u_{z} (mm)$ at A (Figure 16)	$u_{z} (mm)$ at B (Figure 16)	$u_{y} (mm)$
Present	1.4057×10⁻²	1.4034×10⁻²	3.6136×10⁻¹
ABAQUS	1.4104×10⁻²	1.4027×10⁻²	3.6033×10⁻¹

Table 7.

Displacement components at the free end (Load case 2).

Method of analysis	u_x (mm) at A (Figure 16)	u_y (mm) at A (Figure 16)
Present	1.6059×10⁻³	3.7625×10⁻³
ABAQUS	1.6071×10⁻³	3.7652×10⁻³

Figures 17 to 21 illustrate the distribution of normal and shear stresses at two representative sections (mid-span section and a section 5 m away from the free end) of the beam for the two load cases, respectively. The stress distributions on the cross-section for both load cases predicted by the present method agreed very well with the 3D FE results. In Figure 18, both of the methods show the localised change of the shear stress distribution at the flange–web junctions.

Figure 17.

Normal stress (σ_zz kN/m²) at the mid-span section of the beam (Figure 16) (Load case 1): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 18.

Transverse shear stress (σ_yz kN/m²) at the mid-span section of the beam (Figure 16) (Load case 1): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 19.

Transverse shear stress (σ_xz kN/m²) at the mid-span section of the beam (Figure 16) (Load case 1): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 20.

Transverse shear stress (σ_xz kN/m²) at a section 5 m away from the free end of the beam (Figure 16) (Load case 2): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 21.

Transverse shear stress (σ_yz kN/m²) at a section 5 m away from the free end of the beam (Figure 16) (Load case 2): (a) proposed method and (b) 3D FE model (ABAQUS).

This beam was also analysed under a uniformly distributed load of 1 kN/m along the beam length, applied at the centre line of the top flange over the entire beam axis. The vertical deflection at the free end of the beam obtained by the proposed approach is 2.73 mm which is in close agreement with the 2.71 mm calculated by the 3D FE model (0.39% difference).

Thin-walled girder having combined open and closed sections

Bridge decks having combined closed and open sections are quite common and such a structure as shown in Figure 22 was considered in this section.

Figure 22.

Bridge girder with a combination of closed and open sections (dimensions are in m).

This problem was also solved under two load cases. For convergence of results, a mesh of 202 eight-node quadratic FEs with a total of 3027 DOFs was required for the 2D cross-sectional analysis. 1D analysis was performed using 15 beam elements and a total of 186 DOFs. In contrast, the DOF required for the 3D FE model was 435,351 (almost 150 times more). The same two load cases were applied as for the previous 3D FE model example; the shear force ( $F_{y} = - 1.0 kN$ ) is applied as uniform traction on the three web plates at the free end of the beam for the first load case and the twisting moment ( $M_{z} = - 1.0 kN m$ ) is applied using MPC in ABAQUS.

The displacements at points on the beam’s free end obtained by both analysis techniques are reported in Tables 8 and 9. Table 8 shows that the present method captures the variation of the longitudinal displacement (u_z) (due to distortion as well as out-of-plane warping) along the flanges. These results also agree well with those obtained by the 3D FE analysis (maximum difference = 0.67%). Table 9 shows the performance of the present approach in estimating the horizontal and vertical displacements (u_x and u_y) of the thin-walled beam section due to the twisting moment (maximum difference = 2.7%). Since these displacement variations are due to the in-plane and out-of-plane warping of the thin-walled sections, they cannot be captured by conventional beam models. However, the present analysis technique, which is consistent with the 3D deformation of the structure, can predict the structural response of the beam in a degree of detail which is not feasible with a conventional beam analysis method. At the same time, the degree of computational efficiency enjoyed by the present approach cannot be achieved in the 3D FE analysis.

Table 8.

Displacement components (mm) at the free end (Load case 1).

Method of analysis	$u_{z}$ at A (Figure 22)	$u_{z}$ at B (Figure 22)	$u_{z}$ at C (Figure 22)	$u_{z}$ at D (Figure 22)	$u_{z}$ at E (Figure 22)	$u_{y}$ at C (Figure 22)
Present	8.6317×10⁻³	8.7078×10⁻³	8.7280×10⁻³	−1.3326×10⁻²	−1.3294×10⁻²	−2.8419×10⁻¹
ABAQUS	8.5746×10⁻³	8.6933×10⁻³	8.7284×10⁻³	−1.3363×10⁻²	−1.3292×10⁻²	−2.8344×10⁻¹

Table 9.

Displacement component (mm) at the free end (Load case 2).

Method of analysis	$u_{x}$ at A (Figure 22)	$u_{x}$ at D (Figure 22)	$u_{y}$ at A (Figure 22)	$u_{y}$ at D (Figure 22)
Present	1.5211×10⁻³	−1.6879×10⁻³	6.8154×10⁻³	3.7592×10⁻³
ABAQUS	1.4799×10⁻³	−1.7071×10⁻³	6.7688×10⁻³	3.7334×10⁻³

The stress distributions at mid-span section of the girder obtained by both methods are plotted in Figures 23 to 27 for the two load cases. These results again emphasise the ability of the present analysis technique to predict the 3D response of these beams accurately with a significant reduction in computational effort.

Figure 23.

Normal stress (σ_zz kN/m²) at the mid-span section of the beam (Figure 22) (Load case 1): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 24.

Transverse shear stress (σ_yz kN/m²) at the mid-span section of the beam (Figure 22) (Load case 1): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 25.

Transverse shear stress (σ_xz kN/m²) at the mid-span section of the beam (Figure 22) (Load case 1): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 26.

Transverse shear stress (σ_xz kN/m²) at a section 5 m away from the free end of the beam (Figure 22) (Load case 2): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 27.

Transverse shear stress (σ_yz kN/m²) at a section 5 m away from the free end of the beam (Figure 22) (Load case 2): (a) proposed method and (b) 3D FE model (ABAQUS).

Figure 26 shows that the shear stress in the flange plates in the open part of the section is close to zero. This is expected as closed section boxes will absorb the maximum torsional moment due to significantly higher torsional rigidity of these boxes.

Conclusion

In this article, an innovative modelling technique is presented for analysis of thin-walled box girder bridge deck structures which is consistent with 3D elasticity solution. The technique is suitable for analysis of any slender beam-like structure. The method decomposes the general 3D problem defining the behaviour of beam-like slender structures into a 2D cross-sectional analysis and a 1D beam analysis, which are carried out using a 2D plane FE model and a 1D beam FE model, respectively. This decomposition makes the presented approach computationally more efficient than a conventional 3D FE analysis. For the examples considered in this article, fewer DOFs (e.g. 150 times) were required to give similarly accurate predictions. Importantly, the accuracy of the solution is also not affected as the method does not adopt any major assumptions which are quite common in most of the beam theories. From 2D analysis, a cross-sectional stiffness matrix is obtained which can account for all of the geometrical couplings in multi-cell bridge girders and material couplings in anisotropic beams. Apart from roto-translational displacements, the 2D sectional analysis can obtain 3D warping displacements consisting of out-of-plane warping as well as in-plane warping (i.e. sectional distortion). Once the cross-sectional stiffness matrix is evaluated using the 2D FE analysis, it can be repeatedly used for the 1D analysis of the beam for different loading, boundary and other conditions. This method can be specifically used in the preliminary design of these structures when the engineer needs to find the contribution of different actions to the overall response of the structure.

The proposed modelling technique was used to solve numerical examples of thin-walled steel box girders which were also analysed via detailed 3D FEM. An example of a girder having a solid rectangular section was also solved at the beginning in order to demonstrate its general applicability. The results were shown to have a very good correlation in all situations. Considering the level of accuracy and efficiency required for the analysis and design of bridge superstructures, the presented modelling approach seems to have very good potential in its application for static and dynamic analyses of a wide variety of bridge configurations (e.g. straight, curved, composite, etc.).

Footnotes

Acknowledgements

The authors would like to acknowledge that the research presented in this paper is supported by the Australian Government through Endeavour Postgraduate Award bestowed to the lead author of this paper. The authors would also like to acknowledge that the technical discussions with Dr. Jose Pedro Blasques from National Laboratory for Sustainable Energy of Technical University of Denmark at different stages of the present research have really helped to complete the work successfully. The authors sincerely thank Dr. Blasques for his unconditional help.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship and/or publication of this article.

References

Bathe

(1982) Finite Element Procedures in Engineering Analysis. Englewood Cliffs, NJ: Prentice Hall.

Bauld

Tzeng

(1984) A Vlasov theory for fiber-reinforced beams with thin-walled open cross sections. International Journal of Solids and Structures 20(3): 277–297.

Berdichevskii

(1979) Variational-asymptotic method of constructing a theory of shells: PMM vol. 43, no.4, 1979, pp. 664–687. Journal of Applied Mathematics and Mechanics 43(4): 711–736.

Blasques

(2014) Multi-material topology optimization of laminated composite beams with eigenfrequency constraints. Composite Structures 111: 45–55.

Borri

Merlini

(1986) A large displacement formulation for anisotropic beam analysis. Meccanica 21(1): 30–37.

Cesnik

CES

Hodges

(1997) VABS: a new concept for composite rotor blade cross-sectional modeling. Journal of the American Helicopter Society 42(1): 27–38.

Chandra

Stemple

Chopra

(1990) Thin-walled composite beams under bending, torsion, and extensional load. AIAA Journal 27(7): 619–626.

Gere

Timoshenko

(1990) Mechanics of Materials PWS. Amsterdam: KENT Publishing Company, Elsevier Science BV.

Ghiringhelli

Mantegazza

(1994) Linear, straight and untwisted anisotropic beam section properties from solid finite elements. Composites Engineering 4(12): 1225–1239.

10.

Giavotto

Borri

Mantegazza

. (1983) Anisotropic beam theory and applications. Computers & Structures 16(1–4): 403–413.

11.

Hodges

Atilgan

Cesnik

CES

. (1992) On a simplified strain energy function for geometrically nonlinear behaviour of anisotropic beams. Composites Engineering 2(5–7): 513–526.

12.

Kermani

Waldron

(1993) Analysis of continuous box girder bridges including the effects of distortion. Computers & Structures 47(3): 427–440.

13.

Lee

(2005) Flexural analysis of thin-walled composite beams using shear-deformable beam theory. Composite Structures 70(2): 212–222.

14.

Prathap

Babu

(1986) Field-consistent strain interpolations for the quadratic shear flexible beam element. International Journal for Numerical Methods in Engineering 23(11): 1973–1984.

15.

Razaqpur

(1991) Thin-walled multicell box-girder finite element. Journal of Structural Engineering 117(10): 2953–2971.

16.

Razaqpur

(1994) Refined analysis of curved thin-walled multicell box girders. Computers & Structures 53(1): 131–142.

17.

Sheikh

Thomsen

(2008) An efficient beam element for the analysis of laminated composite beams of thin-walled open and closed cross sections. Composites Science and Technology 68: 2273–2281.

18.

Timoshenko

Goodier

(1970) Theory of Elasticity. New York: McGraw-Hill, p. 51, vol. 19.

19.

Vlasov

(1961) Thin-Walled Elastic Beams. Washington, DC: Office of Technical Services, U.S. Department of Commerce.