NREL is a national laboratory of the U.S. Department of Energy

Office of Energy Efficiency & Renewable Energy

Operated by the Alliance for Sustainable Energy, LLC

This report is available at no cost from the National Renewable Energy

Laboratory (NREL) at www.nrel.gov/publications.

Contract No. DE-AC36-08GO28308

BeamDyn: A High-Fidelity Wind 

Turbine Blade Solver in the 

FAST Modular Framework

Preprint

Q. Wang, M. Sprague, and J. Jonkman

National Renewable Energy Laboratory

N. Johnson

Colorado School of Mines

To be presented at SciTech 2015

Kissimmee, Florida

January 5‒9, 2015

Conference Paper

NREL/CP-5000-63165

January 2015
NOTICE

The submitted manuscript has been offered by an employee of the Alliance for Sustainable Energy, LLC 

(Alliance), a contractor of the US Government under Contract No. DE-AC36-08GO28308. Accordingly, the US 

Government and Alliance retain a nonexclusive royalty-free license to publish or reproduce the published form of 

this contribution, or allow others to do so, for US Government purposes.

This report was prepared as an account of work sponsored by an agency of the United States government. 

Neither the United States government nor any agency thereof, nor any of their employees, makes any warranty, 

express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of 

any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately 

owned rights. Reference herein to any specific commercial product, process, or service by trade name, 

trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, 

or favoring by the United States government or any agency thereof. The views and opinions of authors 

expressed herein do not necessarily state or reflect those of the United States government or any agency thereof.

This report is available at no cost from the National Renewable Energy

Laboratory (NREL) at www.nrel.gov/publications.

Available electronically at http://www.osti.gov/scitech

Available for a processing fee to U.S. Department of Energy

and its contractors, in paper, from:

U.S. Department of Energy

Office of Scientific and Technical Information

P.O. Box 62

Oak Ridge, TN 37831-0062

phone: 865.576.8401

fax: 865.576.5728

email: mailto:reports@adonis.osti.gov

Available for sale to the public, in paper, from:

U.S. Department of Commerce

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22161

phone: 800.553.6847

fax: 703.605.6900

email: orders@ntis.fedworld.gov

online ordering: http://www.ntis.gov/help/ordermethods.aspx

Cover Photos: (left to right) photo by Pat Corkery, NREL 16416, photo from SunEdison, NREL 17423, photo by Pat Corkery, NREL 

16560, photo by Dennis Schroeder, NREL 17613, photo by Dean Armstrong, NREL 17436, photo by Pat Corkery, NREL 17721.

NREL prints on paper that contains recycled content.
BeamDyn: A High-Fidelity Wind Turbine Blade Solver in the 

FAST Modular Framework 

Qi Wang∗1, Nick Johnson†2, Michael A. Sprague‡1, and Jason Jonkman§1 

1National Renewable Energy Laboratory, Golden, CO 80401 

2Colorado School of Mines, Golden, CO 80401 

BeamDyn, a Legendre-spectral-finite-element implementation of geometrically exact beam theory (GEBT), 

was developed to meet the design challenges associated with highly flexible composite wind turbine blades. In 

this paper, the governing equations of GEBT are reformulated into a nonlinear state-space form to support its 

coupling within the modular framework of the FAST wind turbine computer-aided engineering (CAE) tool. 

Different time integration schemes (implicit and explicit) were implemented and examined for wind turbine 

analysis. Numerical examples are presented to demonstrate the capability of this new beam solver. An example 

analysis of a realistic wind turbine blade, the CX-100, is also presented as validation. 

I. Introduction 

In recent years, wind power installations in the United States have exceeded 60 gigawatts in generation capacity, 

and have become an increasingly important part of the overall energy portfolio. Simultaneously, the size of wind 

turbines has also increased in the quest for economies of scale. To ensure the performance and reliability of wind 

turbines, it is crucial to use computer-aided engineering (CAE) tools that are capable of analyzing wind turbine blades 

accurately and efficiently. Although modern computers enable three-dimensional (3D) analysis of a fully resolved 

blade, such analyses are too expensive for iterative and probabilistic design. However, composite wind turbine blades 

are well represented as nonlinear beam models, which can capture the deformation response under realistic operating 

conditions with high fidelity, and in a small fraction of the time required by a fully resolved 3D simulation. 

Beam models are widely used to analyze structures in which one dimension is much larger than the other two. 

Many engineering structures such as bridges, joists, and helicopter rotor blades may be modeled as beams. Similarly, 

beam models are well suited for analyzing wind turbine blades, towers, and shafts. Most wind turbine blades are 

constructed of composite materials, which are more complicated than isotropic beams to analyze because of elastic 

coupling effects, high flexibility, and initial twist/curvatures. The geometrically exact beam theory (GEBT), first proposed by Reissner 1, is a beam-deformation model useful in efficient analysis of highly flexible composite structures. 

For instance, GEBT has demonstrated efficacy in helicopter rotor analysis 2. Simo3 and Simo and Vu-Quoc 4 extended 

Reissner’s work to include 3D dynamic problems. Jeleni ́c and Crisfield5 derived a finite-element (FE) method that 

interpolates the rotation field, thereby preserving the geometric exactness of this theory. Notably, Ibrahimbegovi ́c and 

his colleagues implemented this theory for static 6 and dynamic 7 analysis. Readers are referred to Hodges 8, in which 

comprehensive derivations and discussions on nonlinear composite-beam theories can be found. Recently, a mixed FE 

formulation of GEBT along with the numerical implementation was presented by Yu and Blair 9. 

FAST is a CAE tool developed by the National Renewable Energy Laboratory (NREL) for analyzing both landbased and offshore wind turbines under realistic operating conditions. The current beam model in FAST is not capable 

of predictive analysis of highly flexible, composite wind turbine blades. Recently, FAST has been reformulated under a 

new modularized framework that provides a rigorous means by which various mathematical systems are implemented 

in distinct modules. These modules are interconnected to solve for the globally coupled dynamic responses of wind 

turbines and wind plants 10,11 . 

In this paper, a 3D displacement-based implementation of geometrically exact beam theory using Legendre spectral 

finite elements is presented. The theory is reformulated in a nonlinear state-space form for the purpose of integrating 

with the FAST framework, thereby introducing an optional high-fidelity beam model as an alternative to the current 

beam model. Several first-order form time integrators, including the explicit fourth-order Runge-Kutta (RK4) and implicit second-order Adams-Moulton (AM2), are examined. This work builds on previous efforts that demonstrated the 

implementation of GEBT and spatial discretization using Legendre spectral finite elements (LSFEs) 12–15 in analyzing 

∗Research Engineer, National Wind Technology Center, AIAA Member. Email: Qi.Wang2@nrel.gov 

†Graduate Research Assistant, Department of Mechanical Engineering. 

‡Senior Research Scientist, Computational Science Center. 

§Senior Engineer, National Wind Technology Center, AIAA Professional Member. 

1 of 17 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
′ h ̇ − F = f (1) 

′ ′ g ̇ +  ̃uh − M x +  ̃u F = m (2)  ̇ + ( ̃ ′ )

T

0 

h u ̇ = M (3) g ω

{ F } {ǫ}

= C (4) M κ

composite wind turbine blades. The paper is organized as follows. First, the theoretical foundation of the geometrically exact beam theory along with the reformulation of the governing equations into a state-space form is introduced. 

Coupling to the FAST framework is then discussed. Finally, numerical examples are provided to verify and validate 

the accuracy and efficiency of the present model for composite wind turbine blades. 

II. Geometrically Exact Beam Theory 

This section briefly reviews the geometrically exact beam theory. Further details on the content of this section can 

be found in many other papers and textbooks 8,9,16 . Figure 1 shows a beam in its initial undeformed and deformed 

states. A reference frame bi is introduced along the beam axis for the undeformed state and a frame Bi is introduced 

along each point of the deformed beam axis. The curvilinear coordinate x1 defines the intrinsic parameterization of 

the reference line. In this paper, matrix notation is used to denote vectorial or vectorial-like quantities. For example, 

B1 

B2 

B3 

Deformed State 

Undeformed State 

r

R

R 

s

r u 

x1 

b1 

b2 

b3 

Rˆ 

rˆ 

Figure 1: A beam deformation schematic. 

an underline denotes a vector u, a bar denotes unit vector n ̄, and a double underline denotes a tensor Δ. Note that 

sometimes the underlines only denote the dimension of the corresponding matrix. The governing equations of motion 

for geometrically exact beam theory can be written as 16 

where h and g are the linear and angular momenta resolved in the inertial coordinate system, respectively; F and M 

are the beam’s sectional force and moment resultants, respectively; u is the one-dimensional (1D) displacement of 

a point on the reference line; x0 is the position vector of a point along the beam’s reference line; and f and m are 

the distributed force and moment applied to the beam structure. The notation (•)′ indicates a derivative with respect 

to beam axis x1 and (•

 ̇) indicates a derivative with respect to time. The tilde operator (e•) defines a skew-symmetric 

tensor corresponding to the given vector. In the literature, it is also termed as “cross-product matrix”. For example, 

 0 −n3 n2  

n e = n3 0 −n1   −n2 n1 0 

The constitutive equations relate the velocities to the momenta and the 1D strain measures to the sectional resultants 

as 

{ } { }

2 of 17 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
where M and C are the 6×6 sectional mass and stiffness matrices, respectively (note that they are not really tensors); ǫ 

and κ are the 1D strains and curvatures, respectively; and, ω is the angular velocity vector that is defined by the rotation 

tensor R as ω = axial(R ̇ RT ). The axial vector a associated with a second-order tensor A is denoted a = axial(A) 

and its components are defined as 

a1

 A32 − A23 



1     a = axial(A) = a2 = A13 − A31 (5) 

 2

a3

 A21 − A12  

The 1D strain measures are defined as 

′ ′ {ǫ

} {

x0 + u − (R R ) ̄ı1

}

= 0 (6) κ k 

where k = axial[(RR0)′ (RR0)T ] is the sectional curvature vector resolved in the inertial basis and  ̄ı1 is the unit 

vector along x1 direction in the inertial basis. Note that these three sets of equations, including equations of motion 

Eq. (1) and (2), constitutive equations Eq. (3) and (4), and kinematical equations Eq. (6), provide a full mathematical 

description of elasticity problems. 

For a displacement-based finite-element implementation, there are six degrees of freedom at each node: three 

T displacement components and three rotation components. Here, q denotes the elemental displacement array as q = T [

u pT ]

where u is the displacement and p is the rotation-parameter vector. The acceleration array can thus be 

T T defined as a = [

u ̈ ω ̇ T ]

. For nonlinear finite-element analysis, the discretized forms of displacement, velocity, and 

acceleration are written as 

T T q(x1) = N ˆq q = [

u p T ] (7) 

T T v(x1) = N vˆ v = [

u ̇ ωT ] (8) 

T T a(x1) = N aˆ a = [

u ̈ ω ̇ T ] (9) 

where N is the shape function matrix and (ˆ·) denotes a column matrix of nodal values. 

III. State-Space Formulation 

NREL has put considerable effort into improving the modularity of FAST. The modules in the new framework can 

be coupled in one of two ways in the time domain: loose and tight. In the loose coupling scheme, data are exchanged 

between the modules at each coupling step, but each module tracks its own states and integrates its own equations 

with its own solver. In a tightly coupled time-integration scheme, each module sets up its own equations, but the states 

are tracked and integrated by a solver common to all of the modules. To enable the most flexibility, it is useful to 

create modules (such as BeamDyn) so that they can support both loose and tight coupling. More details on the FAST 

modularization framework and an initial assessment of the numerical stability, numerical accuracy, and computational 

performance of various coupling schemes can be found in10,15,17 . 

To accommodate the tight coupling scheme in the FAST modular framework, the governing equations (1) and (2) 

needed to be reformulated into a state-space form. First, these equations are recast in compact form as 

Fext FI − (F C ) ′ + FD = (10) 

where FI , F C and FD, and Fext are the inertial force term, two elastic force terms, and externally applied force term, 

respectively; their definitions are 

{

h ̇

} 

0 0

 {

h

}

FI = +  ̇ (11) g ̇ u ̃ 0 g

{F }

F C = (12) M

{ 0 }

FD = (13) ′ (x  ̃ + u  ̃′ )T

0 F

{f}

Fext = (14) m

3 of 17 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
a =  ̇v = u  ̈

ω ̇ (24) 

A  ̇xˆ(t) = f(ˆx(t), t) (25) 

 

 

The state-space form can be obtained as 

with initial condition 

where 

xˆ(0) = ˆx0 (26) 

A(ˆx(t)) = D 0 

0 M 

D(ˆx(t)) = f l

0 

NT I 3 0 

0 H 

f(ˆx(t), t) = {I l

0 NT v dx1 

F(ˆx(t), t) } 

xˆ0 = {ˆq0 

vˆ0

} 

N dx1 

(27) 

(28) 

(29) 

(30) 

{ } 

Along with the constitutive equations (3) and (4), the inertial force FI can be written explicitly as 

 

mu ̈ + (ω ̃ ̇ + ω ̃ω ̃)mη 

FI = mη ̃u ̈ + ̺ω ̇ + ω ̺̃ω

 

mI mη ̃T {

u ̈

} 

0 mω ̃η ̃T {

u ̇

}

= + (15) mη ̃ ̺ ω ̇ 0 ω ̺̃ ω

≡ Ma + Gv 

where m is the mass density per unit span; η is the center of mass location; ̺ is the moment of inertia; and I is the 

identity matrix. The definitions of the acceleration vector a and velocity vector v can be found in Eq. (9) and (8), 

respectively. Using the newly introduced matrices, the compact form of the equations of motion can be rewritten as 

M a + f(q, v, t) = 0 (16) 

where 

− Fext f(q, v,t) = FF − F C′ + FD (17) 

FF = Gv 

 

0 mω ̃η ̃T {u ̇

}

= (18) 0 ω ̺̃ ω

A weighted residual formulation is used to enforce the dynamic equilibrium conditions in Eq. (16) 

f l

− Fext NT (Ma + FF − FC′ + FD )dx1 = 0 (19) 

0 

The above equation can be recast as 

Maˆ = F(q, v,t) (20) 

where 

f l

M = NTM N dx1 (21) 

0 

f l

+ Fext F(q, v,t) = NT (−FF + F C′ − FD )dx1 (22) 

0 

To derive the state-space form of the governing equations, the state variable in first-order form x(t) is introduced 

as 

{

q(t)

}

x(t) ≡ (23) v(t)

Note that the second component of x(t) is not q ̇ but v because the angular velocity ω cannot be calculated as the time 

derivative of the rotation parameter p. The angular velocity is related to the rotation parameter by the tangent matrix 

as ω = H(p) ̇p16. Substituting the discretized quantities in Eqs. (7) to (9) into Eq. (23) and using the relation 

4 of 17 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
{ǫ ̇

}

f = μ C (31) d κ ̇

{ Fd

}

F C = (32) d Md

{ 0 }

FD = (33) d ′ ( ̃x +  ̃u ′ )T

0 Fd

Δt Δt A k+1(ˆxk+1 − xˆk − xˆ

 ̇

k) = f(ˆxk+1, tk+1) (34) 2 2 

 

K C ˆ ˆ {

Δˆq

} {Fˆ 1

}

1 1 = (35) Kˆ Cˆ Δˆv Fˆ 2 2 2 

5 of 17 

Note that the state-space form, in Eq. (25) and (26), can be solved with any number of ordinary differential equations 

(ODE) integrators for first-order-in-time systems. A viscous damping term is also implemented to account for the 

structural damping. The damping force is defined as 

where μ is a user-provided damping-coefficient diagonal matrix. The damping force can be recast in two separate 

parts, like F C and FD in the elastic force, as 

Readers are refered to Ref.[16] for more details on the damping force and its linearization. 

IV. Implementation of BeamDyn 

The state-space form of GEBT was numerically implemented into a beam solver, called BeamDyn. BeamDyn has 

two sources of input: the user-prescribed parameters and per time step data from the FAST driver code. The input 

parameters from the users are: 

1. Key points that define the natural geometry of the beam; these key points are interpolated by cubic splines 

2. Sectional constants for a 1D beam model 

3. Finite-element mesh information (order of element) 

4. Selection of static 14 or dynamic analysis 

5. Time integrators and associated parameters: RK4, AM2, or Generalized-α (for loosing coupling only). 

It is intended to couple separate instances of BeamDyn to FAST for each blade. Structurally, the coupling will 

involve passing motions of the blade root for each blade from FAST’s ElastoDyn structural module to BeamDyn; 

including position, velocity, and acceleration (both translation and rotation); and passing root loads for each blade from 

BeamDyn to ElastoDyn (both forces and moments). Aerodynamically, the coupling will involve passing motions of 

the nodes distributed along each blade from BeamDyn to FAST’s AeroDyn aerodynamics module; including position, 

velocity, and acceleration (both translation and rotation); and passing aerodynamic loads distributed along the blade 

from AeroDyn to BeamDyn. FAST will also send BeamDyn the gravity vector and initial hub motions at initialization 

to enable an initialization of BeamDyn’s states. While not needed to enable the coupling of BeamDyn to FAST, 

BeamDyn will also calculate the internal loads along the beam, calculated by subtracting the inertial loads from the 

externally applied loads. 

The implementation of GEBT with RK4 is straightforward; however, for the AM2 scheme, a linearization is 

required because of its implicit nature. By applying the trapezoidal rule, the state-space form in Eq. (25) can be recast 

as 

where the subscript denotes the discretized time step in which the value is evaluated, and Δt is the time step size. A 

linearization is needed to solve this nonlinear equation. The linearized equation is as follows 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
 

 

 

 

 

 

 

6 of 17 

The coefficient matrices are written as 

f l

ˆ NT K = K N dx1

(36) 1 1 0 

f l

ˆ NT C

= C N dx1 (37) 1 1 0 

f

l Δt Δt Δt ˆ NT K = (A − A + A + A + Q)N dx1 2 2 3 1 5 0

2 2 2 

f l Δt Δt Δt ′ ′T ′ + NT P N + N S N + NT O N dx1

(38) 

0 2 2 2 

f l Δt ˆ NT C

= (M + A ) N dx1 (39) 2

2 4 0 

where the newly introduced matrices are 

I0 

K =

3 Δ (40) t 1 0 B(p, p) + H − B(p, p ) − B(p, p ̇ ) k 2 k

Δt I 0 C

= − 3 (41) 1 2 0 I 3 

0 mω ̃ ̇ kη ̃T

A= (42) 1 0 mu ̈ 

 ̃kη ̃ + ρ ̃ ̇ ρω ̇ ωk − k

0 mω ̃η ̃T

A= (43) 2 0 mu ̃ ̇ η ̃ + ρω ̃ − ρ ω 

0 mω ̃kη ̃T

A= (44) 3 0 mu ̃ ̇ kη ̃ + ρω ̃k − ρ ωk

wT

ηT 0 ω ̃mη + ω ̃m ̃ A

= (45) 4 0 ρω

T + ω ̃ρ 

0 ω ̃ω ̃mη ̃T

A= (46) 5 0 ω ̃(ρω ̃ − ρ ω 

The O, P, Q, and S can be found in previous work by Wang et al 14. It is also noted that the B matrix is a linearization 

of the tangent matrix H defined as 

∂(H(p))a 

B(p, a) =(47) ∂p 

One major difficulty in implementing geometrically nonlinear theory is the description of rotations in 3D space. As 

discussed in many articles 16,18,19, there will be a singular point for any type of vectorial parameterization. To address 

this singular point, we implemented an algorithm proposed by Crisfield and Jeleni ́c 18, in which the relative rotations 

are interpolated instead of the total rotations within one element. Although the magnitude of the total rotation may be 

large, it is reasonable to assume that the relative rotation in one finite element is smaller than π. In the time domain, 

all the nodal rotations need to be rescaled at the same time so the highly nonlinear governing equations in Eq. (1) and 

(2) are linearized about the correct point. 

In summary, the newly developed BeamDyn has the following features: 

1. It is based on the state-space form of geometrically exact beam theory, thus, it can be easily integrated into 

FAST as a structural module with both loose-coupling and tight-coupling schemes using time integrators for 

first-order form ordinary differential equations (ODEs) 

2. The spectral finite-element method is used to discretize the space domain; exponential convergence rate can be 

expected for a smooth solution 

3. It is compatible with the FAST modularization network. Variables are categorized into input, output, states, and 

parameters. 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
7 of 17 

V. BeamDyn Verification and Validation 

In this section, several numerical cases to verify and validate BeamDyn against numerical solutions and experimental data are presented. As discussed in the introduction, modern wind turbine blade features include initial 

twists/curvatures, elastic coupling, and high flexibility. BeamDyn’s capabilities in addressing these features are examined. 

A. Initially Twisted/Curved Beams 

An initially twisted beam was examined first. A straight beam (k2 = k3 = 0) with an initial twist (k1 �= 0) is shown 

in Figure 2. The beam is linearly twisted in the positive θ1 direction from 0 degrees at the root to 90 degrees at the 

tip. Table 1 shows the material properties for A36 steel, the beam geometry, and the force applied at the free tip along 

x3 direction. The height and base values reported in the table are the height and base of the rectangular cross section. 

The beam was discretized using a seventh-order LSFE to obtain a converged results. It is widely known that the wellrefined 3D finite-element solution can be considered exact; here we use the 3D results as benchmark solutions. The 

results for the twisted beam are shown in Table 2 and compared to the baseline results obtained from extremely refined 

3D ANSYS analysis using SOLID186 elements. We see that the tip results are quite large (35% of total length) and 

that the agreement between the solid-element and beam-element solutions is very good. These results also serve to 

validate the beam-modeling approach. 

x3 

x2 x2 

x3 x3 

x2 

x1 

x3 

Figure 2: Sketch of an initially twisted beam. 

Table 1: Properties of twisted beam 

Property Value 

Elastic Modulus 200 GPa 

Shear Modulus 79.3 GPa 

Height 0.5 m 

Base 0.25 m 

Length 10 m 

Force 4000 kN 

Table 2: Comparison of tip displacements of an initially twisted beam 

u1 (m) u2 (m) u3 (m) 

BeamDyn -1.132727 -1.715123 -3.578671 

ANSYS -1.134192 -1.714467 -3.584232 

Percent Error 0.129% 0.038% 0.155% 

Next, an initially curved beam was examined. It is clear that the initial curvature plays a major role in the distribution of the elastic forces within the beam. As such, it is very important to ensure that BeamDyn is capable of modeling 

this effect properly. A widely used benchmark problem for a curved beam is the case proposed by Bathe 20, which was 

analyzed for verification. Figure 3 shows the configuration of the cantilevered curved beam being analyzed. The beam 

lies in the x1, x2 plane, the positive x1 direction and the negative x2 direction. A force of 600 pounds is applied in 

the positive x3 direction. The beam is defined by the 45-degree arc with 100-inch radius centered at 100 inches in the 

negative x2 direction. The geometry of the cross section for the curved beam is square, and the material properties can 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
be found in Ref.[20]. The beam was discretized by a fifth-order LSFE. The results of this static analysis are shown in 

Table 3 and are compared to the results published in Bathe 20 . 

Figure 3: Sketch of an initially curved beam. 

Table 3: Comparison of tip displacements of an initially curved beam 

u1 (inches) u2 (inches) u3 (inches) 

BeamDyn (one LSFE) -23.7 13.5 53.4 

Bathe-Bolourchi 20 -23.5 13.4 53.4 

It can be seen from these results that the simulations from BeamDyn for a initially curved beam match quite well 

with the published results. The results of these cases suggest that BeamDyn is capable of modeling beams with initial 

twist and curvatures. 

B. Static Analysis of CX-100 Blade 

The primary intended use of BeamDyn is analyzing anisotropic wind turbine blades, therefore, the CX-100 was chosen 

as a validation case. The CX-100 is a 9-m blade designed by Sandia National Laboratories 21. It was chosen for this 

analysis because it is well characterized, with a wealth of publicly available data regarding its construction and material 

properties. 

The cross-sectional properties as derived from VABS for this beam were provided by D.J. Luscher of Los Alamos 

National Laboratory. These properties were provided at 40 points along the beam. For example, a typical stiffness 

matrix corresponding to 2.2 m along the span of the blade is given by 



193, 000 −75.4 12.2 −75.2 −1970 −3500  

C = 103 × 













 

−75.4 

12.2 

−75.2 

−1, 970 

19, 500 

4, 760 

62.6 

67.3 

4, 760 

7, 210 

−450 

17.0 

62.6 

−450 

518 

1.66 

67.3 

17.0 

1.66 

2, 280 

11.3 

2.68 

−1.11 

−879 













 −3, 500 11.6 2.68 −1.11 −875 4, 240 

where the units associated with stiffness values are Cij (N), Ci,j+3 (N · m), and Ci+3,j+3 (N · m2) for i, j = 1, 2, 3. 

Figure 4 shows the various material layups and the geometry of the CX-100 blade. Each color represents a section 

with unique material properties. Figure 5 shows the normalized bending stiffness along the length of the blade where 

a sharp gradient in the beam axial direction can be observed. The graph shows that the bending stiffness jumped to 

about 10% of its root value in under 10% of the blade length. 

Figure 6 shows the test configuration for the static test performed at the National Wind Technology Center (NWTC) 

in Boulder, Colorado. The whiffle-tree configuration of the test apparatus applied the load at 3.00 m, 5.81 m, and 7.26 

m from the root of the blade to achieve a maximum root moment of 128.6 kN m. The loads and positions are given in 

Table 4 below. 

The out-of-plane displacements, u3, at each of the load points were tracked for the experiment and are given in 

Table 5. The BeamDyn simulation was completed using four seventh-order LSFEs and the results are given in Table 5. 

8 of 17 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
Figure 4: Material layup and geometry of the CX-100 wind turbine blade 21 . 

Figure 5: Normalized bending stiffness as a function of normalized blade length. 

Figure 6: Test configuration for the static pull test conducted at the NWTC21 . 

Table 4: Positions and applied loads during the CX-100 static loads testing at the NWTC 

Saddle # Radial Position (m) Applied Load (kN) 

1 3.00 16.9 

2 5.81 5.47 

3 7.26 5.59 

9 of 17 

American Institute of Aeronautics and Astronautics 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
1.4 

1.2 

1 

0.8 

0.6 

0.4 

0.2 

0 

LOCATION ALONG BLADE AXIS (m) 

FLAPWISE DISPLACEMENTS (m) 

Experimental Data 

BeamDyn 

NLBeam 

0 1 2 3 4 5 6 7 8 9 

10 of 17 

Table 5: Experimental and BeamDyn simulation results for the CX-100 static test 

u3 at saddle #1 (m) u3 at saddle #2 (m) u3 at saddle #3 (m) 

Experimental 0.083530 0.381996 0.632460 

BeamDyn 0.072056 0.381074 0.698850 

Percent Error 13.74% 0.24% 10.5% 

The displacements are plotted in Figure 7 and are overall in good agreement. One reason to explain the discrepancies 

between BeamDyn results and experimental data is that the coupling effects between the degrees of freedom was 

difficult to measure in the experiments. The results obtained by NLBeam, a nonlinear beam solver developed by Los 

Alamos National Laboratory, also can be found in Figure 722. While the NLBeam solution was understood to use the 

same sectional data as the BeamDyn solution, it is currently unknown why the NLBeam solution differs. The focus of 

this case was to validate BeamDyn, not verify BeamDyn against NLBeam, and no attempt has been made to find the 

source of the differences. 

Figure 7: Displacement u3 along the length of the blade for experimental data and BeamDyn simulation. 

Next, a convergence study of the tip displacements was completed for the CX-100 blade in BeamDyn. Figure 8 

shows the error as a function of the number of nodes (5 nodes represents 2 nodes per each of the 4 elements, 9 nodes 

represents 3 nodes per each of the 4 elements, etc.). The percentage error was calculated against the experimental data 

in Table 5. It can be seen that the convergence rate is not exponential as desired. This is likely due to the fact that as 

more nodes are added, different material-property stations are used; those data are not smoothly distributed in space. 

For the next simulation, the blade was meshed such that the element boundaries coincide with the locations where the 

sectional properties are defined. As previously mentioned, the cross-sectional properties for the CX-100 blade were 

given at 40 locations along the length of the blade. To have an element coincide with each sectional property, 39 LSFEs 

were needed. Figure 9 shows the results of this simulation. Each circle on the plot indicates an additional order of the 

LSFE, with the maximum being six. These results demonstrated that exponential convergence was achieved with this 

simulation, albeit with many elements. It can therefore be stated that for composite beams with sharp gradients and 

erratic data in the cross-sectional stiffness matrix, the spectral convergence is compromised, unless element boundaries 

coincident with the discontinuities. It should be noted here that although the convergence rate suffers as a result of 

sharp gradients and erratic data in the cross-sectional stiffness matrix as other finite element 22, the simulations still 

return reasonable results in comparison to the experiment data, suggesting that a lack of exponential convergence does 

not compromise the utility of BeamDyn. 

C. Dynamics of a Cantilevered Beam 

For verification of BeamDyn’s dynamic analysis capability using the AM2 integrator, a benchmark problem for damping effects in the Dymore release package is used as an example. This example involves a dynamic analysis of a 

uniform cantilevered beam under an impulsive excitation, which is shown in Figure 10. The excitation was applied at 

the free tip along both the x2 and x3 directions simultaneously. The simulation range was 0 to 1 second. The sectional 

stiffness constants can be found in Table 6 and the units were the same as those used in the previous example. This 

beam is 2.4 meters long. The inertia properties were: 1.61 kg/m for the unit mass density, 8.60×10−4 and 1.10×10−4 

kg m2 for the mass moment of inertia about x2 and x3 direction, respectively. Two cases were analyzed: one without 

damping effects and the other with damping coefficients of 8.0 × 10−4 s for all the six degrees of freedom. The time 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
Figure 8: Percent error in u3 as a function of the number of nodes. 

Figure 9: Percent error in u3 compared to a highly refined solution in BeamDyn as a function of the number of nodes 

for 39 first- to sixth-order LSFEs, coincident with sectional properties. 

11 of 17 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
12 of 17 

steps for these two cases are 1 × 10−4 and 5 × 10−4 seconds, respectively, given that the damping will help the convergence and stability during simulation. The beam was discretized by one fifth-order element in BeamDyn analysis 

for a converged result. The time histories of all the displacement and rotation components are plotted in Figure 11. 

The damping effects can be clearly observed from these results. Furthermore, the accuracy of the present results of 

damped case were examined by the root-mean-square (RMS) errors, which aggregates the magnitudes of the errors in 

predictions for various times into a single measure of predictive performance. The error was calculated using 

max

 Ln [uk − ub(tk)]

2 

k=0 3 εRMS = max Ln (48) [ub(tk)]

2 

k=0 

where ub(t) is the benchmark solution given by Dymore using eight third-order elements and 5×10−4s time increment. 

The RMS error for three displacement components were 1.80 × 10−3 , 1.32 × 10−3, and 2.70 × 10−3, respectively. 

TIME (s) 

Figure 10: Impulsive excitation. 

0.0E+00 

5.0E+02 

1.0E+03 

1.5E+03 

2.0E+03 

2.5E+03 

3.0E+03 

3.5E+03 

4.0E+03

 0 0.02 0.04 0.06 0.08 0.1 

FORCE (N) 

Table 6: Stiffness constants of the cantilevered beam in Example C 

C11 C22 C33 C44 C55 C66 

4.35E+07 1.40E+07 2.81E+06 2.81E+04 2.33E+04 2.99E+05 

D. Analysis of NREL 5-MW Blade 

The last example is to examine the efficiency and accuracy of the different time integrators that have been implemented 

in BeamDyn. The NREL 5-MW blade was analyzed in a cantilever condition, while a white-noise force is applied at 

the free tip along the flap direction. White noise was used to mimic the magnitude and frequency content of applied 

aerodynamic loads under realistic turbulent conditions. The damping coefficient is set to 10−3 s. Figure 12 shows 

the time history of the applied force and its power spectral density (PSD). The flap direction responses are plotted in 

Figure 13. Four cases were conducted in BeamDyn: the first three cases were using the AM2 time integrator with time 

step size: 2.5 × 10−2 s, 5 × 10−3 s, and 5 × 10−4 s, respectively, while the fourth case was using RK4 time integrator 

with time step sizes 5 × 10−6 s. It can be observed that for an implicit AM2 time step beyond 0.005 s, the solution 

is nearly identical to the fully resolved explicit RK4 solution. For an AM2 time step of 0.025 s, the solution remains 

stable and tracks the other solutions, but error grows at higher frequencies. The spikes at 0.7 Hz and 2 Hz correspond 

to the first and second blade flapwise natural frequencies, respectively. The spike above 5 Hz–above the frequency 

range of excitation–is brought about by nonlinear effects. 

The convergence rate against time step size for the AM2 integrator is also examined. The result from the RK4 

simulation is taken as the benchmark solution. The normalized RMS error defined in Eq. (48) is plotted in Figure 14, 

in which second-order convergence as a function of time increment size can be observed. It is also pointed out that a 

full Newton-Raphson algorithm has been implemented in BeamDyn. The total and average numbers of linear system 

solves, which are related to the computational time required for the solution, are plotted in Figure 15. It can be 

observed that while the total number of linear system solves reduces with increasing time step, the average number of 

solves per step increases. This means that the computation time cannot be cut in two by doubling the time step. 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
-0.2 

-0.15 

-0.1 

-0.05

 0

 0.05

-0.5

 0

 0.5

 1

-0.6 

-0.4 

-0.2

 0

 0.2

 0.4

 0.6

DISPLACEMENT (m) ROTATION PARAMETER DISPLACEMENT (m) 

0.08 Undamped u1 Damped u1 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Undamped u2 Damped u2 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

0.07

 0.06

 0.05

 0.04

 0.03

 0.02

 0.01

 0

ROTATION PARAMETER ROTATION PARAMETER DISPLACEMENT (m) 

-0.25 

-0.01

-0.3 -0.02 

TIME (s) TIME (s) 

(a) u1 (b) u2

0.02

-1 

Undamped u3 Damped u3 

-0.05 

Undamped p1 Damped p1 

0.01

 0

-0.01

-0.02 

-0.03 

-0.04 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

TIME (s) TIME (s) 

(c) u3 (d) p1

0.8 0.04

-0.8 

Undamped p2 Damped p2 Undamped p3 Damped p3 

-0.005

 0 0.1 0.2 0.3 0.4 0.5 0.6

TIME (s) 0.7 0.8 0.9 1 0 0.1 0.2

-0.01 

0.3 0.4 0.5 0.6 0.7 0.8

TIME (s) 

(e) p2 (f) p3

0.9 1 

0.035

 0.03

 0.025

 0.02

 0.015

 0.01

 0.005

 0

Figure 11: Tip displacement and rotation histories of a cantilever beam under impulsive excitation. 

13 of 17 

American Institute of Aeronautics and Astronautics 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
14 of 17 

American Institute of Aeronautics and Astronautics 

x 105 

6 

4 

2 

0 

−2 

−4 

−6 

−8 0 10 20 30 40 50 60 70 80 90 100 

FORCE (N) 

TIME (s) 

(a) Time history of applied force 

10−1 100 101 102 103 0 

1 

2 

3 

4 

5 

6 

7 x 109 

PSD OF APPLIED FORCE (N2 /Hz) 

FREQUENCY (Hz) 

(b) PSD of applied force 

Figure 12: Applied white noise force at the tip of a cantilevered NREL 5-MW blade. 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
0 10 20 30 40 50 60 70 80 90 100 −20 

−15 

−10 

−5 

0 

5 

10 

15 

FLAP DISPLACEMENT (m) 

TIME (s) 

Case 1, AM2 ∆t = 2.5E−2 

Case 2, AM2 ∆t = 5E−3 

Case 3, AM2 ∆t = 5E−4 

Case 4, RK4 ∆t = 5E−6 

(a) Time history of flap displacement 

4 

10−8 

10−6 

10−4 

10−2 

100 

102 

10

PSD OF FLAP DISPLACEMENT (m2 /Hz)

Case 1, AM2 ∆t = 2.5E−2 

Case 2, AM2 ∆t = 5E−3 

Case 3, AM2 ∆t = 5E−4 

Case 4, RK4 ∆t = 5E−6 

−1 0 1 2 3 10 10 10 10 10

FREQUENCY (Hz) 

(b) PSD of flap displacement 

Figure 13: Flap direction responses of a cantilevered NREL 5-MW blade under white noise excitation. 

     

) εRMS(uflap 10­4

10­3

10­2

10­1

100

O(∆t

2 )

10­4 10­3 10­2 10­1

TIME STEP SIZE (s)

Figure 14: Normalized RMS error of flapwise displacement histories as a function of time step size for AM2 time 

integrator. The dashed line shows ideal second-order convergence. 

15 of 17 

American Institute of Aeronautics and Astronautics 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
   

  

   

  

   

  

       

  

   

  

   

  

   

10

8

109

108

107 6

TOTAL NUMBER OF LINEAR SYSTEM SOLVES 106

105

3.78E+04

1.79E+05

1.56E+06

8.00E+08

AVERAGE NUMBER OF LINEAR SYSTEM SOLVES

4

2

9.46

8.97

7.78

4.00

10 0 4

Case1 Case2 Case3 Case4	 Case1 Case2 Case3 Case4

AM2 ∆t =2.5E−2 AM2 ∆t =5E−3 AM2 ∆t =5E−4 RK4 ∆t =5E−6 AM2 ∆t =2.5E−2 AM2 ∆t =5E−3 AM2 ∆t =5E−4 RK4 ∆t =5E−6

(a) Total number of linear system solves (b) Average number of linear system solves per step 

Figure 15: Total and average numbers of linear system solves in the simulations. 

VI. Summary 

This paper reported on the development of BeamDyn, a beam solver in the FAST modularization framework. 

Geometrically exact beam theory has been reformulated into the first-order state-space form so that first-order time 

integrators, as required by FAST framework for tight-coupling analysis, can be applied. Numerical examples are 

provided to verify and validate the capability of BeamDyn in analyzing initially curved/twisted beams. A validation 

example is also presented where the numerical results are compared with experimental data. In all the cases, good 

agreement was observed. The convergence rate for analysis of realistic wind turbine blade is discussed. The features 

of BeamDyn can be summarized as follows: 

1. Based on geometrically exact beam theory, BeamDyn is capable of dealing with geometric nonlinear beam 

problems with arbitrary magnitude of displacements and rotations for both static and dynamic analyses 

2. Along with a preprocessor like PreComp or VABS, BeamDyn takes full elastic coupling effects into account 

3. The governing equations are reformulated into state-space form, thus, making it amendable into FAST for tightcoupling analysis 

4. The space is discretized by spectral finite elements, which is a p-version finite element, so that exponential 

convergence rate can be expected for smooth solutions 

5. Different time integrators have been implemented in BeamDyn; users will have options based on their needs 

6. BeamDyn is implemented following the programming requirements (data structures and interfaces) of the FAST 

modularization framework. 

Future work includes the proposal of a practical approach to smooth the geometry and sectional constants to help 

achieve the expected exponential convergence rate. In addition, this module will be coupled to FAST and verification 

and validation will be conducted on a full wind turbine system. 

Acknowledgments 

This work was supported by the U.S. Department of Energy under Contract No. DE-AC36-08GO28308 with the 

National Renewable Energy Laboratory. Funding for the work was provided by the DOE Office of Energy Efficiency 

and Renewable Energy, Wind and Water Power Technologies Office. Support was also provided through a laboratory 

directed research and development grant High-Fidelity Computational Modeling of Wind-Turbine Structural Dynamics. The authors would also like to acknowledge D.J. Luscher for CX-100 sectional data. 

References 

1Reissner, E., “On one-dimensional large-displacement finite-strain beam theory,” Studies in Applied Mathematics LII, 1973, pp. 87–95. 2Hodges, D. H., Saberi, H., and Ormiston, R. A., “Development of Nonlinear Beam Elements for Rotorcraft Comprehensive Analyses,” 

Journal of the American Helicopter Society, Vol. 52, 2007, pp. 36–48. 

16 of 17 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.
3Simo, J. C., “A finite strain beam formulation. The three-dimensional dynamic problem. Part I,” Computer Methods in Applied Mechanics 

and Engineering, Vol. 49, 1985, pp. 55–70. 4Simo, J. C. and Vu-Quoc, L., “A three-dimensional finite-strain rod model. Part II,” Computer Methods in Applied Mechanics and Engineering, Vol. 58, 1986, pp. 79–116. 5Jeleni ́c, G. and Crisfield, M. A., “Geometrically exact 3D beam theory: implementation of a strain-invariant finite element for statics and 

dynamics,” Computer Methods in Applied Mechanics and Engineering, Vol. 171, 1999, pp. 141–171. 6Ibrahimbegovi ́c, A., “On finite element implementation of geometrically nonlinear Reissner’s beam theory: three-dimensional curved beam 

elements,” Computer Methods in Applied Mechanics and Engineering, Vol. 122, 1995, pp. 11–26. 7Ibrahimbegovi ́c, A. and Mikdad, M. A., “Finite rotations in dynamics of beams and implicit time-stepping schemes,” International Journal 

for Numerical Methods in Engineering, Vol. 41, 1998, pp. 781–814. 8Hodges, D. H., Nonlinear Composite Beam Theory, AIAA, 2006. 9Yu, W. and Blair, M., “GEBT: A general-purpose nonlinear analysis tool for composite beams,” Composite Structures, Vol. 94, 2012, 

pp. 2677–2689. 10Jonkman, J. M., “The new modularization framework for the FAST wind turbine CAE tool,” Proceedings of the 51st AIAA Aerospace 

Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Grapevine, Texas, January 2013. 11Jonkman, J. and Jonkman, B., “FAST v8,” https://nwtc.nrel.gov/FAST8, October 2013, [Online; accessed 29-OCTOBER-2014]. 12Wang, Q., Yu, W., and Sprague, M. A., “Gemoetrically nonlinear analysis of composite beams using Wiener-Milenkovi ́c parameters,” 

Proceedings of the 54th Structures, Structural Dynamics, and Materials Conference, Boston, Massachusetts, April 2013. 13Wang, Q. and Sprague, M. A., “A Legendre spectral finite element implementation of geometrically exact beam theory,” Proceedings of the 

54th Structures, Structural Dynamics, and Materials Conference, Boston, Massachusetts, April 2013. 14Wang, Q., Sprague, M. A., Jonkman, J., and Johnson, N., “Nonlinear Legendre spectral finite elements for wind turbine blade dynamics,” 

Proceedings of the 32nd ASME Wind Energy Symposium, National Harbor, Maryland, January 2014. 15Sprague, M. A., Jonkman, J., and Jonkman, B., “FAST modular wind turbine CAE tool: non matching spatial and temporal meshes,” 

Proceedings of the 32nd ASME Wind Energy Symposium, National Harbor, Maryland, January 2014. 16Bauchau, O. A., Flexible Multibody Dynamics, Springer, 2010. 17Gasmi, A., Sprague, M., and Jonkman, J., “Numerical stability and accuracy of temporally coupled multi physics modules in wind-turbine 

CAE tools,” Proceedings of the 51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Grapevine, 

Texas, January 2013. 18Crisfield, M. and Jeleni ́c, G., “Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element 

implementation,” Proceedings of the Royal Society, London: Mathematical, Physical and Engineering Sciences, Vol. 455, 1999, pp. 1125–1147. 19Bottasso, C. and Borri, M., “Integrating finite rotations,” Computer Methods in Applied Mechanics and Engineering, Vol. 164, 1998, 

pp. 307–337. 20Bathe, K. J. and Bolourchi, S., “Large displacement analysis of three-dimensional beam structures,” International Journal for Numerical 

Methods in Engineering, Vol. 14, 1979, pp. 961–986. 21Paquette, J., Lairdl, D., Griffith, D., and Rip, L., “Modeling and testing of 9m research blades,” 44th AIAA Aerospace Sciences Meeting, 

Vol. 19, 2006, pp. 14569–14581. 22Fleming, I. and Luscher, D. J., “A model for the structural dynamic response of the CX-100 wind turbine blade,” Wind Energy, Vol. 17, 

2013, pp. 877–900. 

17 of 17 

American Institute of Aeronautics and Astronautics 

This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.