Development of a cross-sectional finite element for the analysis of thin-walled composite beams like wind turbine blades

Abstract

A method for structural analysis of thin-walled composite beams like wind turbine blades is presented. This method is based on the Nonhomogeneous Anisotropic Beam Section Analysis (NABSA) which consists in discretizing the beam cross section using finite elements. The proposed implementation uses 3-node line cross-sectional finite elements with nodes having rotational degrees of freedom to describe the cross-sectional warping displacements. Solutions obtained using this approach were verified against the corresponding analytical or numerical solutions. Agreement was very good to excellent for the computation of cross-sectional properties and distribution of stresses, strains and warping displacements for a broad range of possible composite beam behaviors including geometric and material couplings, open sections, multicell sections, and arbitrary laminates. For thin-walled layered structures, the proposed method provides models with fewer degrees of freedom than equivalent models based on a two-dimensional discretization of cross sections using triangular or quadrilateral elements such as conventional NABSA or VABS which suggests that computation time could be reduced.

Introduction

Wind turbine blades, as well as airplane wings and helicopter rotor blades, are mostly composite beams. Unlike classical metallic structures, these composite beams have cross-sectional non-homogeneity and offer the possibility of enabling coupling between all deformation modes that may arise from fiber orientation. This latter phenomenon can be used as a load alleviation system, for example when a bend-twist coupling exists in the behavior of a wind turbine blade (Fedorov and Berggreen, 2014). As the structures mentioned earlier are often modeled as beams, especially for aeroelastic analyses, their cross-sectional properties need to be computed. In order to consider all the possibilities associated with composite beams, models more elaborate than the classical strength of material methods applied to homogeneous cross sections are needed.

To develop such models, the Variational Asymptotic Beam Section Analysis available in commercial VABS software is a reference tool (Cesnik and Hodges, 1997; Hodges and Yu, 2007; Yu et al., 2002a, 2002b, 2012).The beam cross section is discretized using two-dimensional finite elements and a 6 × 6 cross-sectional stiffness matrix is calculated. Shear effects are therefore taken into account so the beam can be considered a Timoshenko beam. Giavotto et al. (1983), based on a different methodology, have also developed a tool to compute a 6 × 6 cross-sectional stiffness matrix for a Timoshenko beam using again a two-dimensional finite element discretization of the cross section. This framework is sometimes called Nonhomogeneous Anisotropic Beam Section Analysis (NABSA). This work has been updated in software called BECAS that also includes optimization capabilities (Blasques, 2011; Blasques et al., 2016; Blasques, 2012; Blasques and Stolpe, 2012).

However, as wind turbine blades, airplane wings and helicopter rotor blades are often thin-walled structures, it is possible to use methods based on the theory of strength of materials. Of course, these methods have to be adjusted to take into account the composite nature of these structures. These adjusted methods allow the computation of the cross-sectional stiffness matrix terms by evaluating integrals over the beam cross section based on various deformation hypotheses.

The first distinction that can be made among models based on strength of materials is the way transverse shear effects are taken into account. Formulations based on Timoshenko beams include such effects (Fernandes da Silva et al., 2011; Librescu and Song, 2006; Massa and Barbero, 1998; Pluzsik and Kollar, 2002; Saravanos et al., 2006; Saravia, 2014; Saravia et al., 2012; Sheikh and Thomsen, 2008; Zhang et al., 2012) whereas those based on Euler beams ignore them (Cárdenas et al., 2012; Cardoso et al., 2009; Gökhan Günay and Timarci, 2017; Kollar and Pluzsik, 2002; Victorazzo and De Jesus, 2016; Wang and Zhang, 2014; Wang et al., 2014; Zhang and Wang, 2014).

In both of these cases, different ways to model the torsional behavior of the beam are used. Some authors build their model based on uniform torsion, that is, Bredt theory for thin-walled beams (Fernandes da Silva et al., 2011; Kollar and Pluzsik, 2002; Saravanos et al., 2006; Saravia, 2014; Saravia et al., 2012; Victorazzo and De Jesus, 2016; Wang et al., 2014; Zhang et al., 2012). This means that the beam has a constant section loaded by torques at the ends only and that support effects are neglected. However, this theory is still valid for slightly variable cross sections and slightly variable distributed torques as well as for sections away from the effects of supports. Other models are based on non-uniform torsion, that is, Vlasov torsion, where variable cross section, variable torque distribution or support effects cause restraint to the out-of-plane warping of the cross section (Cárdenas et al., 2012; Cardoso et al., 2009; Gökhan Günay and Timarci, 2017; Librescu and Song, 2006; Massa and Barbero, 1998; Pluzsik and Kollar, 2002; Sheikh and Thomsen, 2008; Wang and Zhang, 2014; Zhang and Wang, 2014). According to Volovoi and Hodges (2000) and Yu et al. (2005), the Vlasov correction for restrained warping has been found unnecessary for closed thin-walled beam as this effect has minimal impact on that kind of structure.

Among the cited models, some are able to consider general composite plate stacking sequences along with all coupling implied (Cárdenas et al., 2012; Gökhan Günay and Timarci, 2017; Kollar and Pluzsik, 2002; Librescu and Song, 2006; Pluzsik and Kollar, 2002; Saravanos et al., 2006; Saravia, 2014; Saravia et al., 2012; Sheikh and Thomsen, 2008; Victorazzo and De Jesus, 2016; Wang and Zhang, 2014; Zhang and Wang, 2014). This allows modeling of the bend-twist coupling of the beam. With some exceptions (Cárdenas et al., 2012; Cardoso et al., 2009; Fernandes da Silva et al., 2011; Wang et al., 2014; Zhang et al., 2012), these models take into account wall bending rigidity.

All the cited models can be used to evaluate closed single-section beams. Some are also able to model open section beams (Cardoso et al., 2009; Kollar and Pluzsik, 2002; Librescu and Song, 2006; Massa and Barbero, 1998; Pluzsik and Kollar, 2002; Sheikh and Thomsen, 2008; Victorazzo and De Jesus, 2016; Zhang et al., 2012) and others, closed multicell section beams (Cárdenas et al., 2012; Fernandes da Silva et al., 2011; Pluzsik and Kollar, 2002; Victorazzo and De Jesus, 2016; Wang and Zhang, 2014; Wang et al., 2014; Zhang et al., 2012).

Coming back to the more complex VABS and NABSA methods, the modeling of thin-walled cross section with these tools (using line finite elements instead of 2D finite elements) is also interesting. Thin-walled implementations of the VABS method have already been proposed. Some are based on classical Euler beam model with (Yu et al., 2005) or without (Volovoi and Hodges, 2000, 2002) the Vlasov correction. Others are based on Timoshenko beam model (Ferede and Abdalla, 2014; Gupta and Hodges, 2017).

The contribution of this paper is a thin-walled implementation of the theory behind NABSA. This leads to a $6\times 6$ Timoshenko cross-sectional stiffness matrix. Instead of discretizing the cross section with 2D elements with three displacement degrees of freedom per node, the section is discretized with line elements containing three translational degrees of freedom and three rotational degrees of freedom per node. This reduced number of degrees of freedom allows a reduction of computation time (not quantified herein), which is valuable in structural optimization. It also allows, at the post-processing stage, the computation of the plate generalized stress (by integrating the stress distribution over the wall thickness) which gives valuable information for the analysis of the beam under the applied load. This thin-walled formulation can also be pre-processed more easily since the direction of laminate stacking is already determined by the cross-sectional mesh.

In the following sections, we review the literature relevant to classical NABSA theory and then describe the proposed finite-element implementation of our method, which we verify by comparison with analytical and numerical (using 3D shell finite elements) solutions.

Review of the theory underlying NABSA

This section presents the NABSA theory as developed by Giavotto et al. (1983) and extended for optimization by Blasques and Stolpe (2012), Blasques (2012, 2011) and Blasques et al. (2016). The nomenclature adopted in this article is similar to that used by the latter group. In this theory, it is assumed that the beam is straight, has a constant cross section and is loaded at its ends only.

Coordinate systems

Figure 1 shows the different coordinate systems used: a global coordinate system $\mathrm{xyz}$ , an element coordinate system $\mathrm{rst}$ and a material coordinate system $123$ . The $x$ - and $y$ -axes are in the cross-sectional plane and the $z$ -axis is the beam reference axis. The $r$ -axis is parallel to $z$ and the $s$ - and $t$ -axes are respectively parallel and normal to the beam wall. Finally, the 3-axis of the material coordinate system is parallel to the $t$ -axis and the 1- and 2-axes are in the plane of the beam wall and constitute the material’s principal direction of orthotropy (1-axis being typically aligned with the direction of the fibers).

Figure 1.

The coordinate systems used in this study.

Stress and strain

The stress and strain components are used in a vectorial form defined as

$$\begin{array}{c}\mathbf{\sigma }={\left[\begin{array}{cccccc}{\sigma }_x\hfill & {\sigma }_y\hfill & {\tau }_{\mathrm{xy}}\hfill & {\tau }_{\mathrm{xz}}\hfill & {\tau }_{\mathrm{yz}}\hfill & {\sigma }_z\hfill \end{array}\right]}^T\hfill \\ \mathbf{ϵ}={\left[\begin{array}{cccccc}{ϵ}_x\hfill & {ϵ}_y\hfill & {\gamma }_{\mathrm{xy}}\hfill & {\gamma }_{\mathrm{xz}}\hfill & {\gamma }_{\mathrm{yz}}\hfill & {ϵ}_z\hfill \end{array}\right]}^T\hfill \end{array}$$

(1)

and the stress-strain relationship is

$$\mathbf{\sigma }=\mathbf{Q}\mathbf{ϵ}$$

(2)

where $\mathbf{Q}$ is the material constitutive matrix. The second half of the stress vector contains the out-of-plane components of stress and is designated as

$$\mathbf{p}={\left[\begin{array}{ccc}{\tau }_{\mathrm{xz}}\hfill & {\tau }_{\mathrm{yz}}\hfill & {\sigma }_z\hfill \end{array}\right]}^T$$

(3)

Internal loads

The beam internal load vector is defined as

$$\mathbf{V}={\left[\begin{array}{cccccc}{V}_x\hfill & V_y\hfill & N\hfill & M_x\hfill & M_y\hfill & M_t\hfill \end{array}\right]}^T$$

(4)

where $V_x$ and $V_y$ are the shear forces, $N$ is the axial force, $M_x$ and $M_y$ are the bending moments and $M_t$ is the torsional moment. The internal loads are the resultants of the out-of-plane components of stress on the blade cross section:

$$\mathbf{V}={\int }_A{\mathbf{Z}}^T\mathbf{p}\mathrm{d}A$$

(5)

where

$$\mathbf{Z}=\left[\begin{array}{cccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & -y\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill & 0\hfill & x\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & y\hfill & -x\hfill & 0\hfill \end{array}\right]$$

(6)

Generalized strain

Defining the beam reference axis displacements ${\chi }_i$ and rotations ${\phi }_i$ as

$$\mathbf{r}={\left[\begin{array}{cccccc}{\chi }_x\hfill & {\chi }_y\hfill & {\chi }_z\hfill & {\phi }_x\hfill & {\phi }_y\hfill & {\phi }_z\hfill \end{array}\right]}^T$$

(7)

and defining the Timoshenko beam generalized strains as

$$\mathbf{\kappa }={\left[\begin{array}{cccccc}{\chi \prime }_x-{\phi }_y\hfill & {\chi \prime }_y+{\phi }_x\hfill & {\chi \prime }_z\hfill & {\phi \prime }_x\hfill & {\phi \prime }_y\hfill & {\phi \prime }_z\hfill \end{array}\right]}^T$$

(8)

It can be shown that

$$\mathbf{\kappa }=\mathbf{Tr}+\mathbf{r}\mathbf{\prime }$$

(9)

where

$$\mathbf{T}=\left[\begin{array}{cccccc}0\hfill & 0\hfill & 0\hfill & 0\hfill & -1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]$$

(10)

and the prime symbol (’) denotes a first derivative with respect to $z$ .

Kinematics of the cross section

The displacement of an arbitrary point in the cross section $\mathbf{s}={\left[\begin{array}{ccc}{s}_x\hfill & s_y\hfill & s_z\hfill \end{array}\right]}^T$ is separated into two terms:

$$\mathbf{s}=\mathbf{Zr}+\mathbf{g}$$

(11)

The term $\mathbf{Zr}$ represents the displacement due to the cross-sectional rigid body motion and the term $\mathbf{g}={\left[\begin{array}{ccc}{g}_x\hfill & g_y\hfill & g_z\hfill \end{array}\right]}^T$ , the displacement due to cross section warping.

Strain-displacement relationships

By separation of the derivative with respect to $x$ and $y$ (in the cross-sectional plane) from the derivative with respect to $z$ (along the beam longitudinal axis), the strain-displacement relationship becomes

$$\mathbf{ϵ}=\mathbf{Bs}+\mathbf{Ss}\mathbf{\prime }$$

(12)

where

$$\mathbf{B}=\left[\begin{array}{ccc}\partial /\partial x\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \partial /\partial y\hfill & 0\hfill \\ \multicolumn{1}{c}{\partial /\partial y}& \partial /\partial x\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \partial /\partial x\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \partial /\partial y\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right],\mathbf{S}=\left[\begin{array}{ccc}0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right]$$

(13)

Using equations (11) and (9), given that $\mathbf{BZ}=\mathbf{SZT}$ , the strain vector becomes

$$\mathbf{ϵ}=\mathbf{SZ}\mathbf{\kappa }+\mathbf{Bg}+\mathbf{Sg}\mathbf{\prime }$$

(14)

Finally, the warping function $\mathbf{g}$ is expressed using a finite element discretization of the cross section so that

$$\mathbf{g}=\mathbf{Nu}$$

(15)

where $\mathbf{N}$ is a matrix containing the shape functions and $\mathbf{u}$ is a vector containing the warping displacements at nodes. The strain vector then becomes

$$\mathbf{ϵ}=\mathbf{SZ}\mathbf{\kappa }+\mathbf{BNu}+\mathbf{SNu}\mathbf{\prime }$$

(16)

Principle of virtual work for a beam slice and solution

Based on the virtual work (per unit length of the beam) of a beam slice

$$\underset{\delta W_i}{\underbrace{{\int }_A\delta {\mathbf{ϵ}}^T\mathbf{\sigma }\mathrm{d}A}}=\underset{\delta W_e}{\underbrace{{\int }_A(\delta {\mathbf{s}}^T\mathbf{p})\prime \mathrm{d}A}}$$

(17)

Giavotto et al. (1983), Blasques and Stolpe (2012), Blasques (2012, 2011) and Blasques et al. (2016) show that the following system is obtained for the center of the beam (far from the boundary conditions):

$$\begin{array}{c}\left[\begin{array}{ccc}\mathbf{E}\hfill & \mathbf{R}\hfill & \mathbf{D}\hfill \\ \multicolumn{1}{c}{{\mathbf{R}}^T}& \mathbf{A}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{{\mathbf{D}}^T}& \mathbf{0}\hfill & \mathbf{0}\hfill \end{array}\right]\left[\begin{array}{c}\mathbf{u}\prime \hfill \\ \multicolumn{1}{c}{\mathbf{\kappa }\prime }\\ \multicolumn{1}{c}{{\mathbf{\lambda }}_1}\end{array}\right]=\left[\begin{array}{c}\mathbf{0}\hfill \\ \multicolumn{1}{c}{{\mathbf{T}}^T\mathbf{V}}\\ \multicolumn{1}{c}{\mathbf{0}}\end{array}\right]\hfill \\ \left[\begin{array}{ccc}\mathbf{E}\hfill & \mathbf{R}\hfill & \mathbf{D}\hfill \\ \multicolumn{1}{c}{{\mathbf{R}}^T}& \mathbf{A}\hfill & \mathbf{0}\hfill \\ \multicolumn{1}{c}{{\mathbf{D}}^T}& \mathbf{0}\hfill & \mathbf{0}\hfill \end{array}\right]\left[\begin{array}{c}\mathbf{u}\prime \hfill \\ \multicolumn{1}{c}{\mathbf{\kappa }\prime }\\ \multicolumn{1}{c}{{\mathbf{\lambda }}_1}\end{array}\right]=\left[\begin{array}{c}(\mathbf{C}-{\mathbf{C}}^T)\mathbf{u}\prime +\mathbf{L}\mathbf{\kappa }\prime \hfill \\ \multicolumn{1}{c}{-{\mathbf{L}}^T\mathbf{u}\prime +\mathbf{V}}\\ \multicolumn{1}{c}{\mathbf{0}}\end{array}\right]\hfill \end{array}$$

(18)

Given the internal load vector $\mathbf{V}=[V_x{V}_{y}N{M}_x{M}_y{M}_t{]}^T$ , this system can be solved to obtain the beam deformation vector $\mathbf{\kappa }=[{\gamma }_{\mathrm{zx}}^0{\gamma }_{\mathrm{zy}}^0{ϵ}_z^0{\kappa }_x{\kappa }_y{\kappa }_z{]}^T$ and the nodal warping displacement vector $\mathbf{u}$ .

The $\mathbf{D}$ matrix adds constraints to rigid body motion of the cross section (already described by the $\mathbf{\kappa }$ vector) using the Lagrange multiplier method. Here, the retained solution differs from the one used by Giavotto et al. (1983), Blasques and Stolpe (2012), Blasques (2012, 2011) and Blasques et al. (2016). Based on the fact that the linear relations to impose between the components of the warping vector are of the same type as the relations generated by an RBE3 element (Delmas, 2012; Reese et al., 2011), the proposed method is to use this kind of interpolation element. This has the advantage of allowing the easy management of the rotation degrees of freedom of the warping vector. The method consists of imposing that the displacements and rotations of a fictitious reference node (located at the origin of the global coordinate system) are zero. Only the displacement degrees of freedom of the connected nodes (all the mesh nodes) are used. An exception arises when all the model nodes are aligned along the $x$ -axis or the $y$ -axis. In this particular case (see verification case 5 below, for instance), the rotational degree of freedom of the connected nodes around the axis on which they are aligned is used. These constraints result in relations of the following form:

$${\mathbf{D}}^T\mathbf{u}=0,{\mathbf{D}}^T\mathbf{u}\mathbf{\prime }=0$$

(19)

The other matrices needed for equation (18) are defined as

$$\begin{array}{c}\underset{6\times 6}{\mathbf{A}}=\sum _e\int {\mathbf{Z}}^T{\mathbf{S}}^T\mathbf{QSZ}\mathrm{d}A,\underset{n\times 6}{\mathbf{R}}=\sum _e\int {\mathbf{N}}_e^T{\mathbf{B}}^T\mathbf{QSZ}\mathrm{d}A,\hfill \\ \underset{n\times n}{\mathbf{E}}=\sum _e\int {\mathbf{N}}_e^T{\mathbf{B}}^T\mathbf{QB}{\mathbf{N}}_e\mathrm{d}A,\underset{n\times n}{\mathbf{C}}=\sum _e\int {\mathbf{N}}_e^T{\mathbf{S}}^T\mathbf{QB}{\mathbf{N}}_e\mathrm{d}A,\hfill \\ \underset{n\times 6}{\mathbf{L}}=\sum _e\int {\mathbf{N}}_e^T{\mathbf{S}}^T\mathbf{QSZ}\mathrm{d}A\hfill \end{array}$$

(20)

where $n$ is the number of degrees of freedom in the model (six per node). The sums are performed over each element $e$ .

Cross-sectional stiffness matrix

To get the beam cross-sectional stiffness matrix ${\mathbf{K}}_s$ and compliance matrix ${\mathbf{F}}_s$ so that $\mathbf{V}={\mathbf{K}}_s\kappa$ , $\mathbf{\kappa }={\mathbf{F}}_s\mathbf{V}$ , and ${\mathbf{K}}_s={\mathbf{F}}_s^{-1}$ , we can solve equation (18) by replacing the $\mathbf{V}$ vector by a 6 × 6 matrix $\hat{\mathbf{V}}$ , in which each column contains a linearly independent load case. This consists in solving six different load cases at the same time. The result for cross-sectional generalized strains $\hat{\mathbf{\kappa }}$ is a matrix with six columns, each one containing the results of a particular load case. Then, $\hat{\mathbf{V}}={\mathbf{K}}_s\hat{\mathbf{\kappa }}$ and

$${\mathbf{K}}_s=\hat{\mathbf{V}}{\hat{\mathbf{\kappa }}}^{-1}$$

(21)

The simplest way to proceed is to choose $\hat{\mathbf{V}}={\mathbf{I}}_6$ , that is, the six different load cases contain each a unit load component and all others are set to zero.

Note that this method of computing the beam cross-sectional stiffness matrix differs from that used by Giavotto et al. (1983), Blasques and Stolpe (2012), Blasques (2012, 2011), and Blasques et al. (2016),which was based on the principle of virtual work.

Implementation of the line finite element

In order to be able to perform a cross-sectional analysis, the matrices of equation (20) have to be evaluated. This section presents the procedure for the computation of a general $n$ -node line element based on a pure displacement formulation similar to the one used for shell finite elements (Ahmad et al., 1970).

Shape functions and element coordinate system at nodes

Figure 2 shows an example of real and reference elements for a 3-node line element. The $\mathrm{rst}$ element coordinate system is defined as follows. $\hat{\mathbf{s}}$ and $\hat{\mathbf{t}}$ are unit vectors, respectively in the $s$ and $t$ directions. They are respectively parallel and perpendicular to the element reference axis. We define the components of vectors $\hat{\mathbf{s}}$ and $\hat{\mathbf{t}}$ at node $i$ as as ${\widehat{\mathbf{s}}}^i$ and ${\hat{\mathbf{t}}}^i$ .

Figure 2.

Real and reference elements (example for $n_n=3$ ).

Coordinate of a point in the element

The $x$ -, $y$ -, and $z$ -coordinates of a point in an element can be computed as a function of the $s$ - and $t$ -coordinates as (shape functions $N_i$ are functions of $s$ ) (Bathe, 2006; Cook et al., 2002):

$$\begin{array}{c}x=\sum _{i=1}^{n_n}{N}_i{x}_i+\frac{\mathrm{ta}}{2}\sum _{i=1}^{n_n}{N}_i{\hat{t}}_x^i\hfill \\ y=\sum _{i=1}^{n_n}{N}_i{y}_i+\frac{\mathrm{ta}}{2}\sum _{i=1}^{n_n}{N}_i{\hat{t}}_y^i\hfill \\ z=\sum _{i=1}^{n_n}{N}_i{z}_i+\frac{\mathrm{ta}}{2}\sum _{i=1}^{n_n}{N}_i{\hat{t}}_z^i\hfill \end{array}$$

(22)

where $n_n$ is the number of nodes in the element. $x_i$ , $y_i$ , and $z_i$ are respectively the $x$ -, $y$ -, and $z$ -coordinates of the $i$ th node. $a$ $=a(s)$ is the element thickness (along the $t$ -direction) and $z_i$ is the $z$ -coordinate of the $i$ th node.

Jacobian matrix

The Jacobian matrix $\mathbf{J}$ is defined as:

$$\left[\begin{array}{c}\partial /\partial s\hfill \\ \multicolumn{1}{c}{\partial /\partial t}\end{array}\right]=\underset{\mathbf{J}}{\underbrace{\left[\begin{array}{cc}\partial x/\partial s\hfill & \partial y/\partial s\hfill \\ \multicolumn{1}{c}{\partial x/\partial t}& \partial y/\partial t\hfill \end{array}\right]}}\left[\begin{array}{c}\partial /\partial x\hfill \\ \multicolumn{1}{c}{\partial /\partial y}\end{array}\right]$$

(23)

Using equation (22), this matrix can be expressed as follows:

$$\mathbf{J}=\sum _{i=1}^{n_n}\left[\begin{array}{cc}{N}_{i,s}{x}_i+\frac{\mathrm{ta}}{2}{N}_{i,s}{\hat{t}}_x^i\hfill & N_{i,s}{y}_i+\frac{\mathrm{ta}}{2}{N}_{i,s}{\hat{t}}_y^i\hfill \\ \multicolumn{1}{c}{\frac{a}{2}{N}_i{\hat{t}}_x^i}& \frac{a}{2}{N}_i{\hat{t}}_y^i\hfill \end{array}\right]$$

(24)

Displacement field, shape function and derivative of shape function matrices

The displacement field can be evaluated by subtracting the final and initial positions of each node using equation (22) (Bathe, 2006).The warping displacement field is then:

$$\mathbf{g}={\mathbf{N}}_e{\mathbf{u}}_e$$

(25)

where the shape function matrix is

$${\mathbf{N}}_e=\left[\begin{array}{cccccccc}\hfill & N_i\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & -\frac{\mathrm{ta}}{2}{N}_i{\hat{t}}_y^i\hfill & \hfill \\ \multicolumn{1}{c}{\dots }& 0\hfill & N_i\hfill & 0\hfill & 0\hfill & 0\hfill & \frac{\mathrm{ta}}{2}{N}_i{\hat{t}}_x^i\hfill & \dots \hfill \\ \hfill & 0\hfill & 0\hfill & N_i\hfill & \frac{\mathrm{ta}}{2}{N}_i{\hat{t}}_y^i\hfill & -\frac{\mathrm{ta}}{2}{N}_i{\hat{t}}_x^i\hfill & 0\hfill & \hfill \end{array}\right]$$

(26)

The nodal displacement vector is

$${\mathbf{u}}_e={\left[\begin{array}{cccccccc}\dots \hfill & u_x^i\hfill & u_y^i\hfill & u_z^i\hfill & {\theta }_x^i\hfill & {\theta }_y^i\hfill & {\theta }_z^i\hfill & \dots \hfill \end{array}\right]}^T$$

(27)

$i$ ranges from 1 to $n_n$ , the number of nodes in the element.

Differentiating the shape function matrix using equation (23) yields:

$${\mathbf{B}}_e{\mathbf{N}}_e=\left[\begin{array}{cccccccc}\hfill & N_{i,x}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & \frac{a}{2}{\hat{t}}_x^i{(t{N}_i)}_{,y}\hfill & \hfill \\ \hfill & 0\hfill & N_{i,y}\hfill & 0\hfill & 0\hfill & 0\hfill & \frac{a}{2}{\hat{t}}_x^i{(t{N}_i)}_{,y}\hfill & \hfill \\ \multicolumn{1}{c}{\dots }& N_{i,y}\hfill & N_{i,x}\hfill & 0\hfill & 0\hfill & 0\hfill & -\frac{a}{2}{\hat{t}}_y^i{(t{N}_i)}_{,y}+\frac{a}{2}{\hat{t}}_x^i{(t{N}_i)}_{,x}\hfill & \dots \hfill \\ \hfill & 0\hfill & 0\hfill & N_{i,x}\hfill & \frac{a}{2}{\hat{t}}_y^i{(t{N}_i)}_{,x}\hfill & -\frac{a}{2}{\hat{t}}_x^i{(t{N}_i)}_{,x}\hfill & 0\hfill & \hfill \\ \hfill & 0\hfill & 0\hfill & N_{i,y}\hfill & \frac{a}{2}{\hat{t}}_y^i{(t{N}_i)}_{,y}\hfill & -\frac{a}{2}{\hat{t}}_x^i{(t{N}_i)}_{,y}\hfill & 0\hfill & \hfill \\ \hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & \hfill \end{array}\right]$$

(28)

where

$$\begin{array}{cc}{N}_{i,x}=j_{11}{N}_{i,s},\hfill & {(t{N}_i)}_{,x}=j_{11}t{N}_{i,s}+j_{12}{N}_i,\hfill \\ \multicolumn{1}{c}{N_{i,y}=j_{21}{N}_{i,s},}& {(t{N}_i)}_{,y}=j_{21}t{N}_{i,s}+j_{22}{N}_i\hfill \end{array}$$

(29)

and $i$ again ranges from 1 to $n_n$ . $j_{11}$ , $j_{12}$ , $j_{21}$ , and $j_{22}$ are the components of the $\mathbf{j}$ matrix, which is the inverse of the Jacobian matrix, that is, $\mathbf{j}={\mathbf{J}}^{-1}$ .

Material properties

Since the stress–strain relationship of the material must be expressed in the global coordinate system ( $\mathrm{xyz}$ ) whereas the elastic properties are known in the material coordinate system ( $123$ ), a transformation of the material constitutive law is required. The stress-strain relationship of a specially orthotropic lamina is

$${\mathbf{Q}}_{123}=\left[\begin{array}{cccccc}\frac{E_1}{1-{\nu }_{12}{\nu }_{21}}\hfill & \frac{{\nu }_{12}{E}_2}{1-{\nu }_{12}{\nu }_{21}}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \hfill & \frac{E_2}{1-{\nu }_{12}{\nu }_{21}}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \hfill & \hfill & G_{12}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \hfill & \hfill & 0\hfill & \frac{5}{6}{G}_{13}\hfill & 0\hfill & 0\hfill \\ \hfill & \mathrm{sym}.\hfill & \hfill & \hfill & \frac{5}{6}{G}_{23}\hfill & 0\hfill \\ \hfill & \hfill & \hfill & \hfill & 0\hfill \end{array}\right]$$

(30)

where ${\nu }_{21}=\frac{E_2}{E_1}{\nu }_{12}$ . The 1- and 2-axes are the principal direction of the material properties of the lamina, where the 1-axis is typically the direction of the fibers and the 3-axis is normal to the plate plane. The $5/6$ factor in front of the shear modulus in planes 13 and 23 is a correcting factor for the transverse shear stress. This value is correct for homogeneous plates and may be adjusted for composite plates (Berthelot, 2012), which was not done here. According to the plane stress hypothesis, there is no rigidity in the thickness direction. Note that a generally orthotropic stress-strain relationship (with no rigidity in the 3-direction) could be used as well.

The transformation of the constitutive matrix to the global coordinate system is achieved in two steps. The first step consists of a transformation in the element coordinate system:

$${\mathbf{Q}}_{\mathrm{rst}}={\mathbf{T}}_{\beta }^T{\mathbf{Q}}_{123}{\mathbf{T}}_{\beta }$$

(31)

and the second step, in the global coordinate system:

$$\mathbf{Q}={\mathbf{T}}_{\alpha }^T{\mathbf{Q}}_{\mathrm{rst}}{\mathbf{T}}_{\alpha }$$

(32)

for which the transformation matrices ${\mathbf{T}}_{\alpha }$ and ${\mathbf{T}}_{\beta }$ depend respectively on the angles $\alpha$ and $\beta$ defined in Figure 1.

Numerical integration

The integration of equation (20) over the area of the elements is computed numerically using Gauss formula. A full integration is used in the $s$ -direction and two integration points per layer are used in the $t$ -direction.

Fictitious rigidity for drilling degrees of freedom

Since no stiffness is associated with the node drilling degree of freedom (rotation about the $\hat{\mathbf{t}}$ -axis, normal to the beam wall), singular or ill-conditioned matrices can arise when two adjacent elements are parallel. To avoid this problem, small stiffness terms are added to the $\mathbf{E}$ matrix (in equation (18)) using a procedure similar to the one presented by Cook et al. (2002).

Post-processing

For a given load case, once the nodal displacements $\mathbf{u}$ , their derivatives $\mathbf{u}\prime$ and the beam generalized strains $\mathbf{\kappa }$ are known, the strains and stresses in each element can be computed using respectively equations (16) and (2).

The stresses and strains can be computed in the element and material coordinate systems as:

$$\begin{array}{c}{\mathbf{\sigma }}_{\mathrm{rst}}={\mathbf{T}}_{\alpha }^{-T}\mathbf{\sigma },{\mathbf{\sigma }}_{123}={\mathbf{T}}_{\beta }^{-T}{\mathbf{\sigma }}_{\mathrm{rst}},\hfill \\ {\mathbf{ϵ}}_{\mathrm{rst}}={\mathbf{T}}_{\alpha }\mathbf{ϵ},{\mathbf{ϵ}}_{123}={\mathbf{T}}_{\beta }{\mathbf{ϵ}}_{\mathrm{rst}}\hfill \end{array}$$

(33)

The shell internal loads can then be computed as integrals of the stress components over the element thickness:

$$N_m={\int }_{-a/2}^{a/2}{\sigma }_m\mathrm{d}\hat{t},Q_n={\int }_{-a/2}^{a/2}{\sigma }_n\mathrm{d}\hat{t},M_m={\int }_{-a/2}^{a/2}\hat{t}{\sigma }_m\mathrm{d}\hat{t}$$

(34)

for $m=(r,s,\mathrm{rs})$ and $n=(\mathrm{rt},\mathrm{st})$ . These integrations are evaluated numerically using two integration points per layer.

Implemented elements

For the verification presented in the next section, a three-node element ( $n_n=3$ ) was implemented using the Python programming language. This element has $n_s=3$ integration points along the element length (full integration) and $n_t=2$ integration points along each layer in the thickness direction. Preliminary results not presented here have shown that the three-node element yields better performance than its two-node counterpart by allowing better representation of actual geometry of a curved cross section wall. Its convergence rate is thus much faster.

Verification of the implemented element

This section proposes different verification cases intended to verify the performance of the proposed method for the analysis of different types of nonhomogeneous anisotropic thin-walled beams. The pre- and post-processing of the models is done with the Gmsh software (Geuzaine and Remacle, 2009). Note that all the gradient results presented (stresses, strains and shell internal loads) are computed at the center of the element, that is, at $s=0$ .

Concerning computation time, the proposed model being coded using Python, all matrices assembly and derived result calculations are done using scripts in an interpreted language (slow) and solution of the equation system is done using a compiled function (fast). It is therefore impossible to get realistic computation time data to compare to the 2D cross-sectional discretization. However, we can suppose that the 1D proposed model will result in smaller computation times for cross section using laminated wall with several layers. 2D models need one element per layer through the thickness, while the proposed 1D model needs only one node through the thickness, resulting in a model with fewer degrees of freedom.

The results presented in the next sections are shown dimensionless but use the International System of Unit: lengths, forces and stresses are respectively in m, N and Pa.

Preliminary verification cases

Some preliminary verification tests on circular and rectangular thin-walled homogeneous sections have shown that the proposed model is able to correctly manage these kinds of problems. All components of the cross-sectional stiffness matrices are the same as the corresponding analytical solutions. For the computation of the analytical transverse shear behavior, the Cowper (1966) shear correction factor is used because it is defined for static analysis such as the one performed here. Stress and strain values returned by the proposed model are also correct.

A convergence analysis has been performed on a homogeneous thin-walled circular cross section of radius $R=1$ and thickness $t=0.01$ made of an isotropic material with $E=E_1=E_2=207\times {10}^9$ , $G=G_{12}=G_{13}=G_{23}=79.3\times {10}^9$ , and $\nu ={\nu }_{12}=E/(2G)-1=0.3052$ . This cross section is centered at the origin of the global coordinate system. This kind of geometry has the advantage of allowing the computation of analytical results to compare against numerical solutions. When looking at the convergence rate for the different terms of the cross-sectional stiffness matrix, we get values of 3.77 for transverse shear stiffness ( $K_{s11}$ and $K_{s22}$ ), 3.86 for axial stiffness ( $K_{s33}$ ), 3.98 for bending stiffness ( $K_{s44}$ and $K_{s55}$ ) and 3.94 for torsional stiffness ( $K_{s66}$ ). Figure 3a shows the convergence of the torsional cross-sectional stiffness $K_{s66}$ against element length where element lengths of 1.571, 0.7854, 0.3927, 0.1963, and 0.09817 correspond respectively to 4, 8, 16, 32, and 64 elements along the circumference. When looking at the stresses and shell internal loads, they show a lower convergence rate, as expected. For this case, Figure 3b shows the convergence of the maximum axial stress ( ${\sigma }_z$ ) value for a unit bending moment against element length.

Figure 3.

Convergence rate of: (a) torsion cross-sectional stiffness $K_{s66}$ and (b) maximum axial stress ${\sigma }_z$ under a unit bending moment. $h$ is the element length and $E$ is the error of the numerical solution defined as $E=\mathrm{abs}(u-u_{\mathrm{an}})/u_{\mathrm{an}}$ where $u$ is the numerical solution and $u_{\mathrm{an}}$ is the analytical solution.

Based on this same thin-walled circular cross section, it is also interesting to study how the proposed model performs when the shell thickness is increased. Figure 4a compares the transverse shear cross-sectional stiffness with respect to element thickness for numerical and analytical solutions. Results show that the numerical results match the analytical result even for a solid cross section ( $t/R=2$ ). Other cross-sectional stiffness terms as well as stress results associated with axial, bending, and torsional loads show similar trends. When looking at maximum shear flow under a unit shear force (Figure 4b), more discrepancy is observed as the thickness increases. However, the error stays under 10% for thickness to radius ratios up to 0.5, which shows that the proposed model is able to analyze relatively thick structures. Those numerical results were obtained using 200 elements and convergence of these models was verified. It’s worth noting that the analytical results presented here are those of the real section geometry, not the thin-walled approximation.

Figure 4.

Comparison of the proposed model solution with analytical solution for increasing thickness to the radius ratio. (a) transverse shear cross-sectional stiffness $K_{s11}$ and (b) maximum shear flow $N_{\mathrm{rs}}$ under a unit shear force. $t$ is the wall thickness and $R$ is the circular cross section mean radius ( $t/R=2$ means a solid circular section with a radius of 2).

Verification case 1

For verification case 1, a thin-walled circular cross section of radius $R=1$ and thickness $t=0.01$ is studied. This cross section is centered at the origin of the global coordinate system and is made with one layer of an orthotropic material (fiber reinforced polymer) with fibers rotated by an angle of 45° relative to the beam longitudinal axis as illustrated in Figure 5. The material properties are $E_1=39.0\times {10}^9$ , $E_2=14.5\times {10}^9$ , ${\nu }_{12}=0.290$ , and $G_{12}=4.24\times {10}^9$ .

Figure 5.

Verification case 1: (a) model geometry and (b) mesh.

The cross-sectional compliance matrix can be obtained analytically by computing the generalized strain vector for six unit internal force load cases. For that analysis, a shear correction factor of 0.5 is used. The resulting compliance matrix is as follows:

$${\mathbf{F}}_s=\left[\begin{array}{cccccc}3485\hfill & 0\hfill & 0\hfill & 689.5\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 3485\hfill & 0\hfill & 0\hfill & 689.5\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1256\hfill & 0\hfill & 0\hfill & -344.9\hfill \\ \multicolumn{1}{c}{689.5}& 0\hfill & 0\hfill & 2511\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 689.5\hfill & 0\hfill & 0\hfill & 2511\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & -344.9\hfill & 0\hfill & 0\hfill & 1743\hfill \end{array}\right]\times {10}^{-12}$$

(35)

The solution obtained from the proposed model using 100 quadratic elements (mesh shown in Figure 5b) is

$${\mathbf{F}}_s=\left[\begin{array}{cccccc}2864\hfill & 0\hfill & 0\hfill & 689.5\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 2864\hfill & 0\hfill & 0\hfill & 689.5\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1256\hfill & 0\hfill & 0\hfill & -344.8\hfill \\ \multicolumn{1}{c}{689.5}& 0\hfill & 0\hfill & 2511\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 689.5\hfill & 0\hfill & 0\hfill & 2511\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & -344.8\hfill & 0\hfill & 0\hfill & 1742\hfill \end{array}\right]\times {10}^{-12}$$

(36)

These solutions are very similar, with differences of less than 0.1% except for the transverse shear compliances ( $F_{s11}$ and $F_{s22}$ ), which differ by about 18%, due probably to the incorrect evaluation of the shear correction factor in the analytical solution which is hard to evaluate due to the presence of off-axis fibers. This shows that the proposed model is able to manage correctly the transverse shear-bending and the extension-torsion couplings due to the presence of off-axis fibers.

Verification case 2

In this verification load case, we study an open thin-walled section. The cross-sectional shape is the same as verification case 1 excepted that the section is open at the point of coordinate $x=R$ and $y=0$ . This section is now made of an isotropic material with $E=E_1=E_2=207\times {10}^9$ , $G=G_{12}=G_{13}=G_{23}=79.3\times {10}^9$ , and $\nu ={\nu }_{12}=E/(2G)-1=0.3052$ .

Using the proposed model with 100 quadratic elements, we obtain the following cross-sectional stiffness matrix:

$${\mathbf{K}}_s=\left[\begin{array}{cccccc}2820\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0.8472\hfill & 0\hfill & 0\hfill & 0\hfill & -1.694\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 13.01\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 6.503\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 6.503\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -1.694\hfill & 0\hfill & 0\hfill & 0\hfill & 3.389\hfill \end{array}\right]\times {10}^9$$

(37)

As expected, we can observe a coupling between the shear force in the $y$ -direction and the torsion ( $K_{s26}$ ). If we compute the position of the shear center, we get the same values as the analytical solution for an open thin-walled beam: $x_s=-2$ and $y_s=0$ . If we transfer the cross-sectional stiffness matrix to an origin located at the shear center, the stiffness matrix becomes:

$${\mathbf{K}}_s=\left[\begin{array}{cccccc}2820\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 847.2\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 13010\hfill & 0\hfill & -26010\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 6503\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & -26010\hfill & 0\hfill & 58520\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0.1658\hfill \end{array}\right]\times {10}^6$$

(38)

and there is now a coupling between the extension and the bending about $y$ -axis ( $K_{s35}$ ), as expected.

As we can see, the torsional stiffness $K_{s66}=165.8\times {10}^3$ is very close to the analytical solution of $\mathrm{GJ}=166.1\times {10}^3$ (0.18% difference). It is also possible to compute analytically the warping displacement at the cross section opening. This analytical solution predicts that the upper and lower parts have respectively warping displacement of $u_z=\pm 18.91\times {10}^{-6}$ . As we can see in Figure 6, the proposed model solution is really close to this value (0.16% difference). If we look at the torsional flow, the analytical solution yields $M_{\mathrm{rs}}=-79.58\times {10}^{-3}$ which is exactly what the numerical model returns.

Figure 6.

Verification case 2. Warping displacement $u_z$ for $M_t=1$ .

All of the differences observed here could be reduced by increasing the number of finite elements in the mesh. This verification case shows that the proposed model correctly manages open sections and geometric shear-torsion coupling.

Verification case 3

This verification case consists of a thin-walled multicell cross section. The geometry is represented in Figure 7a and the material is the same as verification case 2. The analytical solution to this problem shows that the torsional stiffness of this cross section is $\mathrm{GJ}=10.86\times {10}^9$ . When a unit torsional moment is applied, the shear flow in the left cell is $N_{\mathrm{rs}1}=0.08088$ and the shear flow in the right cell is $N_{\mathrm{rs}2}=0.09324$ . The shear flow in the central wall is $N_{\mathrm{rs}2}-N_{\mathrm{rs}1}=0.01236$ .

Figure 7.

Verification case 3: (a) geometry and (b) shear flow $N_{\mathrm{rs}}$ for $M_t=1$ .

With 124 quadratic elements, the proposed model locates the shear center at $x=0.3788$ and $y=0$ . Transferring the cross-sectional stiffness matrix to this point yields $K_{s66}=10.86\times {10}^9$ , which is the same as computed by the analytical solution ( $\mathrm{GJ}$ ). Figure 7b shows the shear flow distribution calculated by the proposed model. With values of 0.08089, 0.09324 and 0.01235, respectively for the left cell, right cell and central wall, the error relative to the analytical solution is limited to less than 0.1%. This shows that the proposed model correctly models multicell cross sections.

Verification case 4

In the previous verification cases, the mesh line always represented the middle of the shell thickness. However, the reference geometry is often located on one of the shell surfaces. Wind turbine blades and airplane wings are examples where the reference geometry is the exterior surface. In such cases it may be more useful to be able to build the model in order for the mesh to represent the bottom of the plate, that is, the thickness is built only on one side of the mesh. This can be achieved by integrating in the thickness direction from $t=0$ to $t=2$ instead of integrating from $t=-1$ to $t=1$ .

By modeling verification case 2 this way, the mesh being located on the outer surface, the predicted shear center location is $x_c=-1.995$ and $y_c=0$ and the cross-sectional stiffness matrix, when transferred to a coordinate system centered at the shear center becomes

$${\mathbf{K}}_s=\left[\begin{array}{cccccc}2824\hfill & 0.0\hfill & 0.0\hfill & 0.0\hfill & 0.0\hfill & 00.0\hfill \\ \multicolumn{1}{c}{0.0}& 975.1\hfill & 0.0\hfill & 0.0\hfill & 0.0\hfill & 0.00\hfill \\ \multicolumn{1}{c}{0.0}& 0.0\hfill & 13010\hfill & 0.0\hfill & -25950\hfill & 00.0\hfill \\ \multicolumn{1}{c}{0.0}& 0.0\hfill & 0.0\hfill & 6503\hfill & 0.0\hfill & 00.0\hfill \\ \multicolumn{1}{c}{0.0}& 0.0\hfill & -25950\hfill & 0.0\hfill & 58270\hfill & 0.00\hfill \\ \multicolumn{1}{c}{0.0}& 0.0\hfill & 0.0\hfill & 0.0\hfill & 00.0\hfill & 0.1659\hfill \end{array}\right]\times {10}^6$$

(39)

The only considerable error is for transverse shear stiffness in the $y$ -direction ( $K_{s22}$ , $+15.1$ %). All other results are exactly the same as those obtained for verification case 2.

Verification case 5

With this verification case, we want to see if the model reacts correctly to a nonsymmetrical and unbalanced layup. Considering the difficulty of obtaining an analytical solution for that kind of layup, the problem studied here consists of a simple plate of unit width in the global $\mathrm{xz}$ -plane. The beam cross section is then a line lying along the global $x$ -axis and 12 quadratic elements are used as mesh. The layup is made of four layers of the same material as verification case 1 with orientation $[0/15/-30/90]$ for a total thickness of 0.01 (each layer having a thickness of 0.0025). Classical lamination theory allows us to see that the $\mathbf{A}$ , $\mathbf{B}$ and $\mathbf{D}$ matrices (not to be confused with the matrices of equations (13), (19), and (20)) of this laminate are fully populated indicating that couplings exist between each of the deformation modes. Subjecting this beam section to a unit axial force is equivalent to subjecting the plate to a membrane force $N_{\mathrm{rr}}=1$ . The stress distribution in this plate should then vary only along its thickness. Table 1 shows the stress distribution along the plate thickness expressed in the layer coordinate system as computed using the classical lamination theory.

Table 1.

Stress distribution in the plate for verification case 5.

Location	${\sigma }_1$	${\sigma }_2$	${\tau }_{12}$	location	${\sigma }_1$	${\sigma }_2$	${\tau }_{12}$
layer 1, bot.	63.38	−12.68	−7.535	layer 3, bot.	102.6	16.18	20.49
layer 1, mid.	87.13	−9.590	−5.069	layer 3, mid.	112.3	23.40	23.74
layer 1, top	110.9	−6.504	−2.604	layer 3, top	122.1	30.62	27.00
layer 2, bot.	95.34	−1.908	−11.08	layer 4, bot.	−22.69	73.47	−7.259
layer 2, mid.	123.0	0.02439	−10.11	layer 4, mid.	−18.69	82.39	−9.724
layer 2, top	150.6	1.957	−9.144	layer 4, top	−14.68	91.32	−12.19

The proposed model gives the same results up to four significant digits. These results are also obtained when using the option to place the nodes at the bottom of the laminate as explained in verification case 4.

Verification case 6

As a final verification case, a nonsymmetrical multicell cross section like that of wind turbine blades is analyzed. As represented in Figure 8a, the cross section uses unbalanced laminates to induce bend-twist coupling (Fedorov et al., 2009). All laminates are centered on the mesh line. So, the nodes of the airfoil contour (laminates 1–4) are shifted toward the interior in order to get the correct geometry. All layers are made of the same material as in verification case 1 with a thickness of 0.001. Using the proposed methodology, when this cross section is meshed using 200 nodes and 101 quadratic elements (see Figure 8b), the resulting cross-sectional stiffness matrix is

$${\mathbf{K}}_s=\left[\begin{array}{cccccc}121.8\hfill & -18.06\hfill & 44.09\hfill & 1.355\hfill & 8.296\hfill & 1.076\hfill \\ \multicolumn{1}{c}{-18.06}& 40.314\hfill & -8.490\hfill & -0.9221\hfill & -1.641\hfill & 1.535\hfill \\ \multicolumn{1}{c}{44.09}& -8.490\hfill & 480.8\hfill & 12.36\hfill & 74.18\hfill & 0.3842\hfill \\ \multicolumn{1}{c}{1.355}& -0.9221\hfill & 12.36\hfill & 6.545\hfill & 8.888\hfill & -0.7557\hfill \\ \multicolumn{1}{c}{8.296}& -1.641\hfill & 74.18\hfill & 8.888\hfill & 47.69\hfill & 0.1570\hfill \\ \multicolumn{1}{c}{1.076}& 1.535\hfill & 0.3842\hfill & -0.7557\hfill & 0.1570\hfill & 4.294\hfill \end{array}\right]\times {10}^6$$

(40)

Figure 8.

Verification case 6. The airfoil is a DU 97-W-300 (Timmer and van Rooij, 2003). Both shear webs are perpendicular to the airfoil chordline and are located at 0.15 on each side of the origin of the coordinate system. For the airfoil contour (laminate 1–4), the bottom of the laminate is toward the airfoil exterior and the top of the laminates are toward the airfoil interior. For the shear webs, the bottom of the laminates are toward the leading edge and the top of the laminates are toward the trailing edge: (a) geometry and (b) mesh.

These results are compared with those of a 3D shell finite element model of a beam of length 50.0. For that model that uses the Altair OptiStruct solver (see Figure 9), the beam is meshed with 299,200 nodes and 101,000 8-node shell elements (101 in the cross section and 1000 along beam length). RBE2 elements are used at both ends to apply loads and supports. The cross-sectional stiffness matrix resulting from this 3D shell finite element model is as follows:

$${\mathbf{K}}_s=\left[\begin{array}{cccccc}124.8\hfill & -18.25\hfill & 45.20\hfill & 1.383\hfill & 8.506\hfill & 1.123\hfill \\ \multicolumn{1}{c}{-18.25}& 42.28\hfill & -8.669\hfill & -0.9712\hfill & -1.681\hfill & 1.631\hfill \\ \multicolumn{1}{c}{45.20}& -8.669\hfill & 481.3\hfill & 12.38\hfill & 74.30\hfill & 0.3971\hfill \\ \multicolumn{1}{c}{1.383}& -0.9712\hfill & 12.38\hfill & 6.546\hfill & 8.895\hfill & -0.7576\hfill \\ \multicolumn{1}{c}{8.506}& -1.681\hfill & 74.30\hfill & 8.895\hfill & 47.72\hfill & 0.1589\hfill \\ \multicolumn{1}{c}{1.123}& 1.631\hfill & 0.3971\hfill & -0.7576\hfill & 0.1589\hfill & 4.299\hfill \end{array}\right]\times {10}^6$$

(41)

Figure 9.

3D shell finite element model of verification case 6: (a) full model and (b) close-up.

These results are obtained by computing the displacements and rotations of groups of nodes at different cross section and differentiate these distributions with respect to the longitudinal axis. This methodology is presented in Forcier and Joncas (2022).

When comparing both models, the differences are (absolute value of the difference normalized by the average of both values):

$$\left[\begin{array}{cccccc}2.3& 0.3& 2.1& 2.3& 2.5& 4.0\\ 0.3& 4.8& 1.2& 4.9& 2.2& 5.6\\ 2.1& 1.2& 0.1& 0.1& 0.1& 3.1\\ 2.3& 4.9& 0.1& 0.0& 0.0& 0.2\\ 2.5& 2.2& 0.1& 0.0& 0.0& 1.4\\ 4.0& 5.6& 3.1& 0.2& 1.4& 0.1\end{array}\right]\%$$

(42)

For the extension-bending behavior ( $K_{s33}$ , $K_{s44}$ , $K_{s55}$ , $K_{s34}$ , $K_{s35}$ , and $K_{s45}$ ), the differences are well under 1%. The differences are more significant for terms associated with transverse shear and torsion but are, however, limited to 6%, which is considered satisfactory.

Figure 10 shows the warping displacement results under a torsional load for both of these models. The warping displacement fields are very similar on both models and the extreme values are within a 5% difference range.

Figure 10.

Verification case 6: Warping displacement $u_z$ for $M_t=1$ : (a) proposed cross-sectional finite element model and (b) 3D shell finite element model (a slice of elements at beam mid-length is presented).

Finally, it is interesting to compare these results with those obtained when using offset nodes on the airfoil contour (nodes on the outside surface and thickness built toward the interior). When using the proposed model, the differences with the solution obtained with mid-surface nodes (equation (40)) are

$$\left[\begin{array}{cccccc}2.9& 1.1& 2.9& 0.9& 3.5& 39.9\\ 1.1& 8.6& 1.7& 18.4& 4.7& 53.5\\ 2.9& 1.7& 0.4& 1.4& 1.2& 38.5\\ 0.9& 18.4& 1.4& 0.4& 1.4& 10.3\\ 3.5& 4.7& 1.2& 1.4& 1.4& 30.5\\ 39.9& 53.5& 38.5& 10.3& 30.5& 10.0\end{array}\right]\%$$

(43)

The difference for extension and bending is small, being limited to 1.4%. When looking at transverse shear and torsion, the differences are higher but limited to 10%. Higher error values are experienced for coupling terms of lower importance ( $K_{s16}$ , $K_{s24}$ , $K_{s26}$ , $K_{s36}$ , and $K_{s56}$ ).

When comparing the difference between the mid-surface node solution (Eq. 41) and the offset node solution of the 3D shell finite element model, the extension and bending behavior are still correctly modeled, but the error is now 36% on the torsional stiffness ( $K_{s66}$ ) and in the same range for transverse shear stiffness ( $K_{s12}$ , $K_{s22}$ , and $K_{s12}$ ). This result is not surprising since the difficulty of obtaining good results for torsion from a shell finite element model using offset shells has been reported by Laird et al. (2005), Branner et al. (2007), Fedorov and Berggreen (2014), and Tavares et al. (2022).

It is interesting to note that the proposed cross-sectional finite element model is less affected by this problem than 3D shell finite elements. Furthermore, as suggested by Tavares et al. (2022), increasing the drilling degree of freedom stiffness reduces the error. In this particular case, the difference between the offset node solution and the mid-surface node solution resulting from the proposed cross-sectional finite element model presented in equation (43) can be reduced to maximum values of 1.4% for extension and bending terms ( $K_{s33}$ , $K_{s44}$ , $K_{s55}$ , $K_{s34}$ , $K_{s35}$ , and $K_{s45}$ ), 1.0% for transverse shear terms ( $K_{s11}$ , $K_{s22}$ , and $K_{s12}$ ), 1.2% for torsion term ( $K_{s66}$ ) and 10.6% for all other coupling terms. These results are satisfactory, but more research needs to be done to optimize the drilling degree of freedom stiffness values of the proposed cross-sectional finite element.

Conclusion

In this paper, a thin-walled implementation of the Nonhomogeneous Anisotropic Beam Section Analysis (NABSA) method has been proposed using a finite element formulation similar to the pure displacement formulation used for shell elements. After a recall of the NABSA method, its thin-walled implementation was presented for a general $n$ -node line element.

A 3-node element implementation was then verified against different cases for which an analytical or numerical solution was obtained. These verification cases covered a broad range of possible composite beam behaviors: geometric coupling, material couplings due to off-axis fiber, open sections, multicell sections, nodes located on the shell bottom surface, nonsymmetrical, and unbalanced laminates. The results obtained from the proposed model showed very good agreement with the analytical and numerical solutions.

The proposed method allows the computation of the beam cross-sectional stiffness matrix as do the classical NABSA or VABS methods does, but facilitates model construction by providing a natural through-the-thickness direction to define the material coordinate system. It also facilitates the computation of the stresses and strains in different laminate layers and of shell internal loads (such as shear flow, for instance). Finally, it provides models that are smaller and compute faster than do classical NABSA or VABS methods using two-dimensional discretization of cross sections with triangular or quadrangular elements.