A correction method for large deflections of cantilever beams with a modal approach
 ^{1}DTU Wind and Energy Systems, Technical University of Denmark, Frederiksborgvej 399, 4000 Roskilde, Denmark
 ^{2}Ørsted, Nesa Allè 1, 2820 Gentofte, Denmark
 ^{1}DTU Wind and Energy Systems, Technical University of Denmark, Frederiksborgvej 399, 4000 Roskilde, Denmark
 ^{2}Ørsted, Nesa Allè 1, 2820 Gentofte, Denmark
Correspondence: Ozan Gözcü (mehgo@orsted.com)
Hide author detailsCorrespondence: Ozan Gözcü (mehgo@orsted.com)
Modalbased reducedorder models are preferred for modeling structures due to their computational efficiency in engineering problems. One of the important limitations of the classic modal approaches is that they are geometrically linear. This study proposes a fast correction method to account for geometric nonlinearities which stem from large deflections in cantilever beams. The method relies on precomputed correction terms and thus adds negligibly small extra computational efforts during the time domain response analyses. The accuracy of the method is examined on a straightbeam model and International Energy Agency (IEA) 15 MW wind turbine blade model. The results show that the proposed method increases the accuracy of modal approaches significantly in secondary deflections due to nonlinearities such as axial and torsional motions for the two studied cases.
Reducedorder models (ROMs) based on a modal approach are used in many structural engineering problems such as wind turbine blades (Hansen, 2015), aircraft wings (Bisplinghoff et al., 2013) and spacecraft (Marshall and Pellegrino, 2021) due to their computational efficiency and reasonable accuracy. These ROMs are based on the smalldeflection assumption; in other words they use constant stiffness, mass and damping matrices which are not updated by deflections unlike nonlinear models. Hence, the accuracy of modalbased ROMs reduces as deflections increase and their errors become significant for applications with large deflections such as long and flexible wind turbine blades or aircraft wings. Moreover, the error in structural response amplifies the error in aeroelastic response and load analysis due to the coupled nature of the problem.
The large deflection effects on aeroelastic stability and loads for wind turbine blades (Kallesøe, 2011; Beardsell et al., 2016; Collier and Sanz, 2016; Riziotis et al., 2008) and aircraft wings (Cesnik et al., 2014) are now well known. Although these effects can be modeled in some aeroelastic tools by using geometric nonlinear structural solvers (Larsen and Hansen, 2019; DNV, 2016; Wang et al., 2017; Bauchau, 2009) where the stiffness matrix terms are a function of deflections or in other words nonlinear, linear modalbased ROMs are still in use even for structures with large deflections such as wind turbine blades (Branlard, 2019; Jonkman and Buhl, 2005; Branlard and Geisler, 2022) due to their speed. Therefore, it is very desirable to increase accuracy with negligible extra computational demand for such tools. This is fully aligned with the main focus of the present work.
The focus of this study is cantilever beam structures and their reducedorder models (ROMs) based on modal approaches used in coupled simulations such as aeroelasticity and load simulation of wind turbines and aircraft. A new correction method for moderately large bending deflection effects is proposed for modalbased ROMs. The method has possible minimum computational cost during the coupled analysis, since it includes only some correction terms which do not require any extra iteration during the response analysis.
Kim et al. (2009)Wang et al. (2013)Jain et al. (2017)Rutzmoser et al. (2017)Wu et al. (2018)Gözcü and Dou (2020)There are many studies in the literature on geometrically nonlinear ROMs, and most of them focus on clamped–clamped beams or simply supported panels. Table 1 shows some prominent works with cantilever structures from the literature together with the methods and examples used in them. The methods shown in Table 1 are either intrusive or nonintrusive. In the former, the nonlinear stiffness matrix terms in modal space are derived analytically, while the latter assumes a cubic relation between modal amplitudes and internal forces and the coefficients are determined with static analysis. As shown in the table, most of the work uses von Kármán kinematics, which are very useful for deriving the nonlinear stiffness terms as a function of displacements, and it captures the moderately large deflections accurately when plate or shell elements are used.
Intrusive and nonintrusive methods require reduction bases where the displacements are a nonlinear function of modal amplitudes. In the given studies a quadratic relation between displacements and modal amplitudes is used with modal derivatives or expansion modes which are called “dual modes” in some studies. All of these studies, except Gözcü and Dou (2020), are limited to beam models with forces applied in a single direction or forces in two bending directions with bending deflections less than 2 % of the beam span. Moreover, the methods used in existing studies require much more computational time than linear ROMs due to the iterative solutions required by nonlinear stiffness matrices. On the other hand, our proposed method uses a linear reduction basis and linear stiffness matrix for computing deflections and corrects the linear deflection results with some quadratic vectors for capturing large bending deflection effects. It may be less accurate compared to the ROMs with nonlinear stiffness terms; however it has a very similar computational speed to the existing linear ROMs and is easy to implement into existing aeroelastic tools compared to nonlinear models.
The focus of nonlinear ROMs in the literature is generally to use the same model for response and stress analysis. On the other hand, this study aims to come with a fast and accurate ROM correction method that is suitable only for response analysis and not for stress analysis. The main contributions of this paper are

the selection of the important physics for the given applications and thus justification of particular method choice (i.e., displacementonly correction),

a proof of method suitability using a highfidelity model with prescribed external loads on both simple and complex geometries.
The paper is divided into the following sections. The kinematics of the cantilever beam problem relevant to this study are explained in Sect. 2, and the proposed method is explained in Sect. 3. Example cases are introduced and their results are given together with discussion in Sect. 4, and the conclusion is given in Sect. 5.
This section explains the kinematics of cantilever beams and geometric nonlinear effects due to large bending deflections for symmetric beams and initially curved beams such as wind turbine blades.
Although there are many research studies about geometrically nonlinear ROMs, most of them focus on clamped–clamped beams or simply supported plates. In the clamped–clamped case, the nonlinearity arises due to the length change where a lateral deflection actually alters the length of the structure which causes an additional stiffness effect such as hardening or softening depending on the sign of length change. This effect is referred to as the membrane effect or bending–extension coupling (Touzé et al., 2021). In the case of cantilever beams, a lateral deflection due to lateral forces does not result in length change (no axial strain), so there is no bending–extension coupling for cantilever beams. However, the free end of a cantilever beam is displaced in the beam span direction to keep the length constant when it bends as shown in Fig. 1. In other words, the cantilever beams with lateral loading go through large rotations which do not result in any strain.
The cantilever beam in its deflected state can display a bending–torsion coupling behavior that is not seen in the undeflected state of the straight symmetric cantilever beam. Figure 2 depicts these kinematics using a straight cantilever beam. The lateral force (F_{x}) deflects the beam to “State1”. Subsequently, the force perpendicular to F_{x} is applied to the beam which is now at “Final State”. It is observed that the beam deflects not only laterally and vertically but also in the torsion direction, even though there is no torsion load. Geometrically linear models (such as ROMs based on a modal approach) can capture bending deflections but not the torsion deflection for this problem, since their stiffness matrices are not updated by new states or deflections.
These kinematics become even more prominent for applications such as wind turbine blades, which are initially curved structures due to prebend. They already have bending–torsion couplings at their undeflected states, but the lateral deflections change the magnitude as well as the direction of this coupling. Figure 3 shows the torsion deflection due to the combination of flapwise and edgewise deflections of a wind turbine blade with prebend at its undeflected position. The edgewise and torsion motions are already coupled for a blade with prebend as shown in the figure. However this coupling first reduces and then changes sign with flapwise bending deflections. As the blade goes through sinusoidal edgewise bending deflections due to gravity and rotation, the torsion induced by edgewise–torsion coupling also alters blade loads and aeroelastic stability (Kallesøe, 2011). These coupling effects become significant for very flexible wind turbine blades. In the case geometrically linear models like classical ROMs based on modal reduction are used, the change in torsion–edgewise coupling cannot be captured correctly for large blade deflections.
Here, the kinematics of a cantilever beam due to bending deflections are explained and kinematics due to axial–bending coupling are left without mentioning since the proposed method is developed only for secondary effects due to large bending deflections of cantilever beams. Hence, the method cannot capture stiffening effects with centrifugal force (Wallrapp and Schwertassek, 1991; Branlard and Geisler, 2022) which comes from bending–extension coupling. However, centrifugal stiffening is already corrected for linear models in existing tools such as Flex (Branlard, 2019) and FAST (Jonkman and Buhl, 2005). In this study, mode shapes at an undeformed state without centrifugal stiffening effects are use, since the mentioned existing tools compute blade mode shapes at undeflected states. However, the method is not limited with undeflected state modal analysis, and it can also work with mode shapes computed at any deflected state.
The nonlinear geometric effects of a cantilever beam can be captured by nonlinear ROMs using different methods (intrusive or nonintrusive; see Table 1) and a reduction basis where linear modes are used together with quadratic vectors such as expansion modes (Hollkamp and Gordon, 2008) or modal derivatives (Idelsohn and Cardona, 1985). It is observed that effects of nonlinear stiffness terms such as stiffening are not significant in bending directions for moderately large deflections of a cantilever beam. On the other hand, the geometric nonlinear effects explained in Sect. 2 are still apparent for cantilever beams.
These kinematics together with other nonlinear effects such as centrifugal stiffening are captured very accurately by geometrically nonlinear models which come with much higher computational cost compared to linear ROMs. Our method uses linear ROMs which use linear mode shapes, stiffness, mass and damping matrices for response calculation, and quadratic vectors are added to the linear response results to capture large bending deflection effects. In the studies where intrusive methods are used, modal derivatives are generally preferred as quadratic vectors for deriving nonlinear stiffness terms in modal space. For nonintrusive methods, the terms for the cubic relation between modal amplitudes and reduced stiffness terms are obtained from some static solutions and the reduction basis includes generally expansion modes. In this study, both modal derivatives and expansion modes are investigated as correction terms for large bending effects. The formulation and calculation process of these vectors are explained below.
3.1 Modal derivatives
The modal derivatives (MDs) are quadratic vectors which include secondary effects that occur due to large deflections (geometric nonlinearities) (Idelsohn and Cardona, 1985; Wu and Tiso, 2014). The quadratic relation needs to be written between physical displacements (u) and modal amplitudes (q) for defining modal derivatives. So, they can be thought of as a second derivative of physical displacements (u) with respect to modal amplitudes (q) when the displacements are represented as a function of modal amplitudes. When Taylor series expansion of the displacements around an equilibrium state (u_{0}=u(q_{0})) is written, the first term (u(q_{0})) is the equilibrium state, the first derivative represents the linear mode shapes and the second derivative term includes the modal derivatives as shown in Eq. (1).
In this study, the equilibrium state is taken as the initial state where u_{0}=0 (undeflected state) and Eq. (1) can be written in terms of mode shapes and their derivatives at the initial state as
where Φ is linear mode shapes, u is the displacement vector and q is the modal amplitude vector. The displacement vector is a function of time (t) and position (x), whereas mode shapes are only a function of position (x) and modal amplitudes are only function of time (t). So, they can be written as u(t,x), q(t) and Φ(x); however their time and position relations are not included in the equations for simplicity. The linear mode shapes (Φ) and corresponding natural frequencies (ω) can be found by the generalized eigenvalue solution,
where K and M are the tangent stiffness and mass matrices and ω_{i} and ϕ_{i} are the ith eigenvalue and corresponding eigenvector (mode shape vector). The tangent stiffness is the Jacobian of the internal force with respect to the displacement. It is a generalization of the linear stiffness matrix in the nonlinear case. Here the tangent stiffness and mass matrices are computed at the initial state (u_{0}=0) for the purpose of solving eigenvalue problems. As seen in the later derivation in this section, the tangent stiffness at the deflected state is also required to compute the modal derivatives. The modal derivatives are calculated by taking the derivative of the eigenvalue problem shown in Eq. (3) with respect to modal amplitude q_{j}.
To determine the derivative of ϕ_{i} and ω_{i} with respect to the jth modal amplitude q_{j}, another equation is needed. This equation can be chosen as Eq. (5), where ϕ_{i} is the massnormalized mode shape vector.
When Eqs. (4) and (5) are combined, the modal derivative of the ith mode shape vector ϕ_{i} and natural frequency ω_{i} with respect to the jth modal amplitude can be determined by
Equation (6) contains all the terms that are required to compute the modal derivatives. The derivation of modal derivatives is similar to the derivation of the sensitivity of eigenmodes and eigenfrequencies with respect to a design variable in structural optimization.
For computation of modal derivatives, the terms related to inertia effects (i.e., mass matrix and its derivatives) are generally ignored, since their contribution to the modal derivatives is very limited (Rutzmoser et al., 2017). The derivative of the eigenvalue can also be assumed to be zero due to the fact the eigenfrequencies of the slender cantilever beams are not sensitive to the vibration amplitude. These assumptions lead to static modal derivatives which are symmetric, whereas modal derivatives computed from Eq. (6) are not necessarily symmetric.
For conciseness, the static modal derivatives are simply called modal derivatives (MDs) hereafter. In Eq. (7), the derivative of the stiffness matrix with respect to the jth modal amplitudes, i.e., $\frac{\partial \mathbf{K}}{\partial {q}_{j}}$, is needed to compute the modal derivatives. In this study, the derivative of the stiffness matrix (K) with respect to the jth modal amplitudes (q_{j}) is computed by central finite differences as
where K(ϕ_{j}δ_{j}) is the tangent stiffness matrix when the system displacements equal ϕ_{j}δ_{j}. The tangent stiffness matrix at the deflected state of the structure by a given amplitude of δ_{j} is computed with a geometrically nonlinear beam solver based on the corotational formulation (https://doi.org/10.5281/zenodo.7533686) in Krenk (2005). The calculation of $\frac{\partial \mathbf{K}}{\partial {q}_{j}}$ from Eq. (8) is the only step where a geometrically nonlinear solver is required. When the stiffness matrix derivatives are ready, the computation of Eq. (7) is straightforward, since the rest of the equation consists of the stiffness matrix and linear eigenvectors at the undeflected state. The tangent stiffness matrix at the undeflected state is actually the linear stiffness matrix which is used in existing modal approaches. There are M_{MD} number of modal derivatives which are computed for M number of linear mode shapes. The relation between M and M_{MD} can be written as
Figures 4 and 5 are used to give visual understanding of the modal derivatives. Figure 4 shows bending of a straight beam with an airfoil crosssection in one direction and its representation by a linear bending mode and its modal derivative by itself. The linear bending mode shape does not have any displacement in the axial (spanwise) direction, so the total (curved) length of the beam increases. The modal derivative vector of the bending mode shape with respect to its modal amplitude includes the axial displacement effect. When the linear bending mode and its modal derivative are summed, the axial displacement and bending effects are captured together.
Figure 5 shows the modal derivatives that are similar to that in Fig. 4 except that the deflections occur in two bending directions at the same time. Hence, together with the axial displacements, there is also a torsion effect due to the couplings at the deflected state. The linear modal approach is only capable of capturing the bending deflections in two lateral (x and y) directions. The axial displacements are captured by the modal derivatives of each bending mode shape with respect to its own modal amplitudes in Fig. 5, whereas the torsion effects are included in the crossmodal derivative vectors which are the derivatives of one bending mode shape with respect to the modal amplitudes of another bending mode in the other direction.
3.2 Expansion modes
Expansion modes are similar to modal derivatives, but they are computed by a nonintrusive approach (Hollkamp and Gordon, 2008). The expansion modes also contain the secondary effects such as the axial and torsion motions caused by the bending motions. These vectors are computed from the difference between the deflections computed by linear ROMs and the nonlinear deflections. The difference is written as a function of the selected order of modal amplitudes, and a quadratic (secondorder) relation is used in this study. Nonlinear deflections for static cases are computed in a geometrically nonlinear solver. In the expansion mode approach, the physical displacements are represented as
where Φ and Φ_{EM} are the matrices whose columns are the linear modes and the expansion modes, x is position, and t is time. The expansion mode amplitudes q_{EM} are quadratic functions of the bending modes and can be written in terms of the M linear mode amplitudes as
where q_{i} is the ith linear mode amplitudes. The relation between the number of linear modes (M) and the number of corresponding expansion modes (M_{EM}) is
So, there are the same number of modal derivatives and expansion modes for the same number of linear modes, since both of these vectors represent the quadratic relation between deflections and modal amplitudes.
While the above Eq. (10) defines the expansion modes, the procedure to compute the expansion modes is given below.
First, a set of static analyses are performed. The applied forces are chosen as a combination of two linear mode shapes (Φ), which can be written as
where λ_{i} and λ_{j} are the prescribed scaling factors for the ith and jth linear mode shape. The nonlinear displacements (u) are computed with a geometrically nonlinear solver for the loads f.
Second, the amplitudes of the linear modes, i.e., q, are computed from the applied loads as
where K and K_{r} are the full and reduced stiffness matrices and f and f_{r} are the full and reduced force vectors. When the applied forces (f) are known, q can be computed.
Once we have the linear modal amplitudes (q), the corresponding expansion modal amplitudes (q_{EM}) can be computed by using Eq. (11).
Then, the nonlinear displacements, the amplitudes of the linear modes and the amplitudes of the expansion modes can be collected into three matrices: $\mathbf{U}\in {\mathbb{R}}^{N\times {N}_{k}}$, $\mathbf{Q}\in {\mathbb{R}}^{M\times {N}_{k}}$ and ${\mathbf{Q}}_{\mathrm{EM}}\in {\mathbb{R}}^{{M}_{\mathrm{EM}}\times {N}_{k}}$. Here each column of the matrices corresponds to a static load case. All load case results can be written in matrix form as
Equation (15) is similar to Eq. (10) except that Eq. (15) is in matrix form where the nonlinear displacements U, the linear modes Φ and the amplitudes of the linear modes Q are known from the above steps. The only unknown in Eq. (15) is the expansion modes Φ_{EM}.
Finally, the expansion modes can be determined by using the leastsquares method:
where $\Vert \phantom{\rule{0.25em}{0ex}}\Vert $ denotes the Euclidean or L^{2} norm of the vector and it is solved for modal expansion matrix Φ_{EM}.
It is noted that the expansion modes are very similar to the modal derivatives. They are identical except the numeric differences if the sets of the ith and jth are symmetric and the force amplitudes λ_{i} and λ_{j} are small enough. The computational cost of the expansion modes is more expensive than the modal derivatives, since their calculation requires solving more static problems than the calculation of the modal derivatives. The difference in the computational cost is small when the number of linear modes (M) is small. However, this difference becomes significant when the number of the linear modes is large.
Further, it is noted that one can extend Eq. (11) to include higherorder expansion modes corresponding to higherorder polynomial terms of the amplitudes of the linear modes, for example, the expansion modes for the cubic and quadratic terms of the linear modal amplitudes.
3.3 Numerical implementation
As mentioned before, the modal derivatives or the expansion modes are not included in the reduction space. Consequently, only the linear mass, stiffness and damping matrices are used in the dynamic response analysis. This is not the case for other studies in the literature (see Table 1), which have used the nonlinear stiffness terms in the modal space.
Algorithm 1 shows the procedure to obtain a structural response with correction vectors for pure structural response analysis where the forces are functions of time and for coupled response analysis where the forces are functions of both time and displacements. In this study, modal derivative and expansion mode vectors are used as correction vectors.

Model generation step (time independent and done once)

Generate the linear ROM from the finite element model:

Compute the linear modes (ϕ_{i}) from $\left(\mathbf{K}{\mathit{\omega}}_{i}^{\mathrm{2}}\mathbf{M}\right){\mathit{\varphi}}_{i}=\mathbf{0}$

Compute the reducedorder stiffness K_{r}, mass M_{r} and damping C_{r} matrices by Galerkin projection:
K_{r}=Φ^{T}KΦ, M_{r}=Φ^{T}MΦ, C_{r}=Φ^{T}CΦ


Compute the correction vectors:
For the modal derivatives see Sect. 3.1
For the expansion modes see Sect. 3.2


Pure structural response analysis (forces are functions of time)

Compute the modal amplitudes for the linear ROM by solving the equation of motion:
${\mathbf{M}}_{\mathrm{r}}\ddot{\mathit{q}}\left(t\right)+{\mathbf{C}}_{\mathrm{r}}\dot{\mathit{q}}\left(t\right)+{\mathbf{K}}_{\mathrm{r}}\mathit{q}\left(t\right)={\mathbf{\Phi}}^{T}\mathit{f}\left(t\right)={\mathit{f}}_{\mathrm{r}}\left(t\right)$ 
Compute the displacements by using the linear modes, the correction vectors and the linear modal amplitudes:
With MDs : $\mathit{u}(x,t)=\mathbf{\Phi}\left(x\right)\cdot \mathit{q}\left(t\right)+\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial \mathbf{\Phi}}{\partial \mathit{q}}\cdot \mathit{q}\left(t\right)\right)\cdot \mathit{q}\left(t\right)$
With EMs: $\mathit{u}(x,t)=\mathbf{\Phi}\left(x\right)\mathit{q}\left(t\right)+{\mathbf{\Phi}}_{\mathrm{E}}\left(x\right){\mathit{q}}_{\mathrm{E}}\left(q\right)$
The calculation of displacements (u(x,t)) does not have to be done at each time step and can be performed as a postprocess after the time simulations.


Coupled response analysis (forces are functions of time and displacements)

Compute modal amplitudes for linear ROM by solving equation of motion:
${\mathbf{M}}_{\mathrm{r}}\ddot{\mathit{q}}\left(t\right)+{\mathbf{C}}_{\mathrm{r}}\dot{\mathit{q}}\left(t\right)+{\mathbf{K}}_{\mathrm{r}}\mathit{q}\left(t\right)={\mathbf{\Phi}}^{T}\mathit{f}(t,\mathit{u})={\mathit{f}}_{\mathrm{r}}(t,\mathit{u})$ 
Compute displacements by using linear modes, correction vectors, and the linear modal amplitudes:
With MDs: $\mathit{u}(x,t)=\mathbf{\Phi}\left(x\right)\cdot \mathit{q}\left(t\right)+\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial \mathbf{\Phi}}{\partial \mathit{q}}\cdot \mathit{q}\left(t\right)\right)\cdot \mathit{q}\left(t\right)$
With EMs: $\mathit{u}(x,t)=\mathbf{\Phi}\left(x\right)\mathit{q}\left(t\right)+{\mathbf{\Phi}}_{\mathrm{E}}\left(x\right){\mathit{q}}_{\mathrm{E}}\left(q\right)$ 
Update the loads (f(t,u)) and apply the loads

The proposed method is demonstrated in the case studies of a straight cantilever beam and the blade of the International Energy Agency (IEA) 15 MW reference wind turbine under static and dynamic loads. The displacement results are given to evaluate the effectiveness of the correction methods based on the modal derivatives and expansion modes to capture the large bending deflection effects. The reducedorder model results are compared with HAWC2 results for both static and dynamic cases. HAWC2 (Larsen and Hansen, 2019) is an aeroservohydroelastic load analysis tool for wind turbines and has been developed by DTU Wind and Energy Systems. It uses multibody formulation with Timoshenko beam and can capture geometric nonlinearities (Pavese et al., 2016). In HAWC2, the z axis is along the axial direction, whereas the x and y directions are lateral directions (see Fig. 3 for the HAWC2 coordinate system).
4.1 Straight beam
A cantilever straightbeam model is used for static and dynamic load cases. Table 2 shows the general properties of the beam, whose crosssection properties are constant along the beam length and shear coefficients are very high compared to the bending stiffness values, so it behaves like an Euler–Bernoulli beam.
Figure 6 shows x, y, z displacement and torsion motion components of the first two modeshape vectors of the beam. Since the beam is straight and has no material coupling, both bending mode shapes have only one lateraldirection displacement without any axial or torsion component. The first mode shape is in the x direction, and the second one is in the y direction.
Figure 7 shows modal derivative (MD) and expansion mode (EM) vectors for first two mode shapes. The results are very similar for MD and EM vectors. These vectors show the sensitivity of a mode shape with respect to a modal amplitude. The vector names are given so that the first number represents the mode number for the linear modeshape vector whose sensitivity is computed with respect to the modal amplitude of the mode number, which is given as the second number in the vector names. Hence, MD11, MD22, EM11 and EM22 show the sensitivity of a mode shape with respect to its own modal amplitude, whereas MD12 and EM12 illustrate the sensitivity of a mode shape with respect to another mode's modal amplitude. Since MD21 and MD12 are the same for static modal derivatives, only MD12 results are shown here. This symmetry case is also valid for expansion modes (EMs), since they are also computed for static cases. MDii and EMii vectors have only axial displacements since their mode shapes are only in one lateral direction. On the other hand, MDij and EMij have only torsion motions, which represent the coupling between two lateral directions.
The first test was carried with static loads in the x direction. The applied force vectors are determined from the stiffness matrix (K), first mode shape Φ_{1} and amplification factor (λ) as
where mode shape Φ_{1} is normalized according to its maximum value.
Expansion mode shape and modal derivative results are the same for this test case; therefore they are given together in Table 3, which shows the tip displacements for the linear ROM without corrections (“Lin.”) and with modal derivative (“MD”) and expansion mode (“EM”) corrections and HAWC2 results. Modal derivative, expansion mode and linear ROM results are the same for lateral (x) deflections, since uncoupled mode shapes result in zero values in lateraldirection components (see Fig. 7) of modal derivative and expansion mode vectors. Linear ROM error in the x direction increases as the deflections increases, and the error is less than 5 % for deflections of 20 % beam length. Moreover the linear model cannot capture any axial (z) displacement, whereas correction factors work quite well. The axial (z) deflection is a secondary effect and is captured by corrections coming from modal derivatives and expansion modes, whereas the linear model estimates zero displacement in the z direction and 100 % error for all solutions. Modal derivative and expansion mode errors are increasing by deflections and are much larger than bending (x) direction errors. The two main sources of error in the secondary (z) direction are truncation errors coming from Eq. (2) and the error in the bending (x) direction, which is even amplified for secondary effects due to the formulation.
Figure 8 shows x and z (axial) positions of the structural nodes computed from HAWC2, linear ROM, MD and EM models for λ=2 and λ=3 load cases. The quadratic corrections capture the secondary effects all along the beam length.
The second test case includes static loads in the x and y directions. The load components are determined by similar formulation given in Eq. (17). The xforce components are amplified by λ=2.5, and the yforce components are amplified by λ=1.0. Figure 9 shows x, y, z displacements and torsion deflections along the beam length. The linear model estimates x and y displacements quite accurately, and quadratic vectors do not alter results in these directions, since they do not have components in these directions (see Fig. 7). On the other hand they capture the secondary effects in the z and torsion directions due to geometric nonlinearities. The linear ROM cannot capture any axial (z) displacements and torsion deflections. On the other hand, MD and EM corrections result in accurate estimation of the secondary effect even though the loads are applied in two directions. Hence, they can also capture the couplings between different directions.
The last load case with the straightbeam model includes dynamic loads in two lateral directions. The xdirection load is constant with λ=2, whereas the ydirection load is computed as
where M is the mass matrix, g is the acceleration vector and ω is a given rotation frequency of the rotor. Note that the beam does not rotate; only the load varies with the given frequency. So, there is no axial (z) force due to centrifugal effects. For this example the gravity vector has only y acceleration with a value of 9.81 m s^{−2}, and ω is taken as 1 rad s^{−1}. The simulation time step is 0.01 s.
Figure 10 shows tip x, y and z displacements and torsion deflections for 100 s. The linear model can capture dynamics in the y direction very accurately, and its xdirection results are very close to HAWC2 results without any fluctuations after 60 s, which cannot be captured without nonlinear stiffness terms (Gözcü and Dou, 2020). However these fluctuations are very small compared to the overall displacements in the x direction. On the other hand, linear ROM results are again zero in axial and torsion directions, whereas MD and EM corrections work very well even for the initial transition period. Moreover, there is no phase difference between HAWC2 and MD/EM results, which indicates that the nonlinear stiffness effects are very small for cantilever beams with moderate deflections. The fluctuations in the x direction due to nonlinear stiffness term effects show the limitation of this method even though the fluctuations are very small.
4.2 IEA 15 MW turbine blade
The blade of the IEA 15 MW wind turbine is modeled as beams in HAWC2 and FAST for load analysis (Gaertner et al., 2020). It is a more complex structure in comparison to the beam studied in the previous section. The blade has material coupling terms and an initially curved shape due to prebend (mostly in the y direction), aerodynamic twist and being backswept (mostly in the x direction). Hence, the blade mode shapes are coupled in the x, y, z and torsion directions unlike the straightbeam example. Besides, their crosssection stiffness properties change over the span and material couplings between lateral bending directions, and torsion brings additional couplings in mode shapes. Table 4 shows the relevant blade properties.
The straightbeam example mode shapes are uncoupled so that modeshape vectors include nonzero values only in one lateral bending direction. However, the mode shapes of the turbine blade are coupled due to its geometry and crosssection material couplings. Figure 11 shows the first two mode shapes of the blade. The first mode shape of the blade is mainly in the y (flapwise) direction with 0.5 Hz, whereas the second mode shape is mainly in the x (edgewise) direction with around 0.7 Hz. Moreover, both mode shapes have components in all directions including the z (axial) and torsion directions. The second mode shape's torsion coupling is much stronger than the first mode shape's coupling, whereas it has weaker coupling in the axial direction than the first mode shape.
Figure 12 shows x, y, z and torsion components of modal derivative and expansion mode vectors for the first two mode shapes. MD and EM vectors are very similar, and they include components in all directions due to the couplings. The difference between MD and EM vectors comes from the numerical computation process. As seen in Sect. 3.2, expansion modes (EMs) are computed via the leastsquares method, which uses all displacement and modal amplitude results gathered in matrices, whereas each modal derivative (MD) is computed independently. Since the IEA 15 MW reference wind turbine blade mode shapes are coupled, there are nonzero values at offdiagonal terms in the displacement matrices used during the computation of EMs. This leads to some numerical differences between MDs and EMs for the IEA 15 MW blade case.
The blade loads include aerodynamic loads at a steady 11 m s^{−1} wind speed, which gives the highest thrust force. The steady aerodynamic loads for a symmetric rotor (no tilt, no yaw) are time independent (static), and their torsion components are not included for this example so that the torsion due to modeshape couplings and dynamic forcing is clearer for this example. On top of the aerodynamic loads, the weight of the blade is applied in the x (edgewise) direction with the formulation given in Eq. (18), where ω is 1.0 rad s^{−1}. There are 15 linear mode shapes for the blade ROM, and quadratic vectors for the first three modes are included only, since these modes have very high modal amplitudes among all modes (see Table 5).
Figure 13 shows the tip displacement results of the blade in the edgewise (x), flapwise (y), axial (z) and torsion directions for linear ROM, HAWC2 and correction models. Since most of the thrust force is in the flapwise direction, the largest displacements occur in that direction with a mean displacement around 13.4 m. The fluctuations in the y direction are mostly due to couplings between mode shapes; in other words they are mostly due to linear stiffness effects and not the nonlinear stiffness effects unlike in the case for straightbeam example (see Fig. 10, xdirection results). Since nonlinear stiffness effects are very small in x and y displacements, the results are very similar in terms of magnitude and phase for all models.
The secondary effects become clear in the z (axial) and torsion directions. The linear model has 0.84 m mean tip axial displacement for the last 50 s due to modeshape couplings (see Fig. 11); however the HAWC2 model has −0.41 m mean tip axial displacements for the same time period. Hence, the linear ROM estimates a 2.5 m larger rotor diameter compared to HAWC2. The modal derivative correction model has −0.29 m and expansion mode correction model has −0.23 m mean axial (z) displacements. They also represent the fluctuations in the z directions much more accurately compared to the linear model, which has almost no fluctuations for the last 50 s. The linear model maximum tip torsion error compared to HAWC2 for the last 50 s is about 1.36^{∘}. The maximum tip torsion errors in MD and EM models for the last 50 s are 0.63 and 0.26^{∘}. The correction modes also capture the phase of torsion motion more accurately than the linear model, which is out of phase with HAWC2 in the torsion direction. In this study, the load time series are prescribed and not updated with deflection. In other words, the results are not obtained for an aeroelastic analysis. On the other hand, the importance of torsion motion will increase for an aeroelastic analysis. The difference in mean torsion results is critical in static aeroelastic analysis where the steady operating conditions are determined. On the other hand, the difference in torsion fluctuations will affect the damageequivalent loads (Gözcuü and Verelst, 2020) and aeroelastic stability analysis (Kallesøe, 2011).
Figures 14 and 15 show the displacement results over the blade span at the two selected time steps where minimum and maximum torsion deflections are obtained in the last 50 s from HAWC2 (shown with dashed black lines). These time steps also correspond to maximum and minimum edgewise displacements. At 70.56 s, where the maximum tip torsion occurs, linear ROM results are quite accurate except z and x displacement (see Fig. 14). On the other hand, MD and EM results are much better than linear ROM results in all directions for that moment. The linear ROM has 1.29 m error in the z direction at the blade tip, whereas MD and EM errors are less than 0.17 m for the blade tip in Fig. 15. The linear ROM predicts a nonphysical elongation in the z direction because of the nonzero motion in the z direction in the linear mode shapes (particularly “Mode1”) shown in Fig. 11. In contrast, the two correction methods in this study are not limited by the linear mode shapes. Both of the correction methods provide a correct estimation of the shortening of the blade in the axial direction.
Figure 14 shows displacement and torsion results along the blade span at the maximum tip torsion moment (67.51 s in Fig. 13). The linear ROM is quite accurate in the x and y directions; however it has low accuracy in the z and torsion directions. Its total torsion error over the blade span is 33.6^{∘}, whereas MD error is 7.7^{∘} and EM error is 1.7^{∘} along the blade span compared to HAWC2 results. The zdirection results are similar to minimumtorsionmoment results in Fig. 15.
One of the findings in Figs. 10 and 13 is that the frequency of the vibration is not affected by the large deflections of the blade. Even the linear model can predict the vibration frequency or the time period of the periodic response in good agreement with HAWC2 simulation in the main deflection directions. In practice, this means that no correction is needed for the vibration frequency in the case of large deflections. This also suits the application of the proposed correction method to a variety of slender cantilever structures including wind turbine blades. Besides, this allows using correction terms to capture secondary effects due to geometric nonlinearities for moderately large deflections.
This study introduced a correction approach for capturing the geometric nonlinear effects in the reducedorder models of cantilever beam structures which go through moderately large deflections. The approach is examined based on two different types of correction vector (i.e., the modal derivatives and expansion modes), leading to similar results. Both correction vectors can be precomputed by using a geometrically nonlinear beam solver and then used in the dynamic response analysis. The advantage of this approach is its low computational cost during the coupled response analysis. It is suitable for integration into the aeroelastic analysis tools that use the modalbased reducedorder models for the cantilever beam structures such as wind turbine blades and airplane wings.
The proposed method is demonstrated by using two case studies, including a simple cantilever beam and the blade for the IEA wind 15 MW reference turbine. For the simple cantilever beam, its mode shapes are uncoupled in the sense that the displacements in the two bending modes are orthogonal to each other. For the turbine blade, its mode shapes are coupled because of its curved and twisted geometry and complex crosssectional properties. In the case study of the simple cantilever beam, the correction vectors based on the modal derivatives and expansion modes have identical results and they captured geometrically nonlinear effects accurately for the lateral deflections for up to 25 % of the beam length. In the case study of the blade of the IEA 15MW reference turbine, modal derivative and expansion mode vectors have small deviations from each other, since they are computed by different numerical processes. The blade axial displacements estimated by the linear ROM are in the opposite direction to the axial displacements in the HAWC2 results. In other words, the linear ROM predicts a nonphysical elongation of the IEA 15 MW blade span for the given load case. In contrast, the two correction vectors in this study give a good estimation of the axial motion of the blade tip. Further, the torsion displacements obtained by the linear ROM are out of phase with the HAWC2 results, whereas the torsion displacements obtained by using the two correction vectors in this study are in phase with the HAWC2 results. The torsion phase difference between the nonlinear solver HAWC2 and the linear ROM is crucial for aeroelastic stability and load analysis of wind turbine blades. The comparison of these results shows that the proposed correction method is capable of accurately estimating the nonlinear dynamic response of the cantilever beam structures with large bending deflections, and it is also suitable for the wind turbine blades with curved and twisted geometry and complex crosssectional properties.
The study can be extended with an implementation of the correction method into an aeroelastic tool that uses the reducedorder models for airplane wings or wind turbine blades. Further, the proposed correction method can also be implemented based on machine learning methods such as neural network models or adaptive kriging methods.
The IEA Wind 15 MW reference wind turbine model can be accessed from GitHub (https://doi.org/10.5281/zenodo.6330754; Barter et al., 2022). The beam formulation and corotational formulation for modeling structures and computing MDs and EMs can be accessed from the GitLab of DTU Wind and Energy Systems (https://doi.org/10.5281/zenodo.7533686; Gözcuü et al., 2022).
OG and EB conceived of the presented idea, developed theory and performed all computations. SD was involved in theoretical discussion and helped in writing formulations and idea presentation. All authors discussed the results and contributed to the final manuscript.
DTU Wind and Energy Systems develops, supports and distributes HAWC2 on commercial terms.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This study has received funding from DTU Wind and Energy Systems and Ørsted A/S.
This work was funded by the DTU and Ørsted cofinanced project called “Large Blade Deflections Modelling” (LaBDeM).
This paper was edited by Mingming Zhang and reviewed by William Collier and Emmanuel Branlard.
Barter, G., Gaertner, E., Bortolotti, P., Abbas, N. J., Rinker, J., Zahle, F., Branlard, E., Lapadron, Hall, M., Dzalkind, and Issraman: IEAWindTask37/IEA15240RWT: IEA Wind 15MW Reference Wind Turbine update, v1.1, Zenodo [data set], https://doi.org/10.5281/zenodo.6330754, 2022. a
Bauchau, O. A.: DYMORE user's manual: Formulation and finite element implementation of beam elements, University of Maryland, Maryland, https://www.nrel.gov/docs/fy06osti/38230.pdf (last access: August 2022), 2009. a
Beardsell, A., Collier, W., and Han, T.: Effect of linear and nonlinear blade modelling techniques on simulated fatigue and extreme loads using Bladed, J. Phys.: Conf. Ser., 753, 042002, https://doi.org/10.1088/17426596/753/4/042002, 2016. a
Bisplinghoff, R. L., Ashley, H., and Halfman, R. L.: Aeroelasticity, Courier Corporation, ISBN 13:9780486691893 ISBN 10:0486691896, 2013. a
Branlard, E. and Geisler, J.: A symbolic framework to obtain midfidelity models of flexible multibody systems with application to horizontalaxis wind turbines, Wind Energ. Sci., 7, 2351–2371, https://doi.org/10.5194/wes723512022, 2022. a, b
Branlard, E. S. P.: Flexible multibody dynamics using joint coordinates and the Rayleigh‐Ritz approximation: The general framework behind and beyond Flex, Wind Energy, 22, 877–893, https://doi.org/10.1002/we.2327, 2019. a, b
Cesnik, C. E., Palacios, R., and Reichenbach, E. Y.: Reexamined structural design procedures for very flexible aircraft, J. Aircraft, 51, 1580–1591, 2014. a
Collier, W. and Sanz, J. M.: Comparison of linear and nonlinear blade model predictions in Bladed to measurement data from GE 6 MW wind turbine, J. Phys.: Conf. Ser., 753, 082004, https://doi.org/10.1088/17426596/753/8/082004, 2016. a
DNV: BLADED 4.8 Theory Manual, DNV GL, 2016. a
Gaertner, E., Rinker, J., Sethuraman, L., Zahle, F., Anderson, B., Barter, G. E., Abbas, N. J., Meng, F., Bortolotti, P., Skrzypinski, W., Scott, G. N., Feil, Ro., Bredmose, H., Dykes, K., Shields, M., Allen, C., and Viselli, A.: IEA wind TCP task 37: definition of the IEA 15megawatt offshore reference wind turbine, Tech. rep., NREL – National Renewable Energy Lab., Golden, CO, USA, https://doi.org/10.2172/1603478, 2020. a
Gözcü, O. and Dou, S.: Reduced order models for wind turbine blades with large deflections, J. Phys.: Conf. Ser., 1618, 052046, https://doi.org/10.1088/17426596/1618/5/052046, 2020. a, b, c
Gözcü, O. and Verelst, D. R.: The effects of blade structural model fidelity on wind turbine load analysis and computation time, Wind Energ. Sci., 5, 503–517, https://doi.org/10.5194/wes55032020, 2020. a
Gözcü, O., Barlas, E., and Richardt, J. D.: Corotational beam formulation, Zenodo [data set], https://doi.org/10.5281/zenodo.7533686, 2022. a
Hansen, M.: Aerodynamics of wind turbines, Routledge, ISBN 9781138775077, 2015. a
Hollkamp, J. J. and Gordon, R. W.: Reducedorder models for nonlinear response prediction: Implicit condensation and expansion, J. Sound Vibrat., 318, 1139–1153, 2008. a, b
Idelsohn, S. R. and Cardona, A.: A reduction method for nonlinear structural dynamic analysis, Comput. Meth. Appl. Mech. Eng., 49, 253–279, 1985. a, b
Jain, S., Tiso, P., Rutzmoser, J. B., and Rixen, D. J.: A quadratic manifold for model order reduction of nonlinear structural dynamics, Comput. Struct., 188, 80–94, 2017. a
Jonkman, J. M. and Buhl Jr., M. L.: FAST User's Guide, NREL – National Renewable Energy Laboratory, https://www.nrel.gov/docs/fy06osti/38230.pdf (last access: August 2022), 2005. a, b
Kallesøe, B. S.: Effect of steady deflections on the aeroelastic stability of a turbine blade, Wind Energy, 14, 209–224, https://doi.org/10.1002/we.413, 2011. a, b, c
Kim, K., Khanna, V., Wang, X., and Mignolet, M.: Nonlinear reduced order modeling of flat cantilevered structures, in: Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 4–7 May 2009, Palm Springs, California, https://doi.org/10.2514/6.20092492, 2009. a
Krenk, S.: Nonlinear Modelling and Analysis of Structures and Solids, Cambridge University Press, https://doi.org/10.1017/CBO9780511812163, 2005. a
Larsen, T. and Hansen, A.: How 2 HAWC2, the user's manual, Tech. Rep. R1597 (ver. 12.8)(EN), DTU Wind Energy, http://tools.windenergy.dtu.dk/HAWC2/manual/How2HAWC2_12_9.pdf (last access: 12 January 2023), 2019. a, b
Marshall, M. and Pellegrino, S.: ReducedOrder Modeling for Flexible Spacecraft Deployment and Dynamics, in: AIAA Scitech 2021 Forum, 11–15 and 19–21 January 2021, online, p. 1385, https://doi.org/10.2514/6.20211385, 2021. a
Pavese, C., Kim, T., Wang, Q., Jonkman, J., and Sprague, M. A.: Hawc2 and beamdyn: Comparison between beam structural models for aeroservoelastic frameworks, Technical Report, NREL – National Renewable Energy Lab., https://www.osti.gov/biblio/1304572 (last access: 12 January 2023), 2016. a
Riziotis, V., Voutsinas, S., Politis, E., Chaviaropoulos, P., Hansen, A. M., Madsen, A., and Rasmussen, F.: Identification of structural nonlinearities due to large deflections on a 5 MW wind turbine blade, in: Proceedings of the EWEC, 31 March–3 April 2008, Brussels, Belgium, 8 pp., https://orbit.dtu.dk/en/publications/identificationofstructuralnonlinearitiesduetolargedeflect#:~:text=Author,BIBTEX,Harvard (last access: 12 January 2023), 2008. a
Rutzmoser, J. B., Rixen, D. J., Tiso, P., and Jain, S.: Generalization of quadratic manifolds for reduced order modeling of nonlinear structural dynamics, Comput. Struct., 192, 196–209, 2017. a, b
Touzé, C., Vizzaccaro, A., and Thomas, O.: Model order reduction methods for geometrically nonlinear structures: a review of nonlinear techniques, Nonlin. Dynam., 105, 1141–1190, https://doi.org/10.1007/s11071021066939, 2021. a
Wallrapp, O. and Schwertassek, R.: Representation of geometric stiffening in multibody system simulation, Int. J. Numer. Meth. Eng., 32, 1833–1850, https://doi.org/10.1002/nme.1620320818, 1991. a
Wang, Q., Sprague, M. A., Jonkman, J., Johnson, N., and Jonkman, B.: BeamDyn: a highfidelity wind turbine blade solver in the FAST modular framework, Wind Energy, 20, 1439–1462, https://doi.org/10.1002/we.2101, 2017. a
Wang, X., Perez, R. A., and Mignolet, M. P.: Nonlinear reduced order modeling of complex wing models, in: Proceedings of the 54th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 8–11 April 2013, Boston, Massachusetts, AIAA 20131520, https://doi.org/10.2514/6.20131520, 2013. a
Wu, L. and Tiso, P.: Modal Derivatives based Reduction Method for Finite Deflections in Floating Frame, in: Proceedings – WCCM XI: 11th World Congress on Computational Mechanics; ECCM V: 5th European Conference on Computational Mechanics; ECFD VI: 6th European Conference on Computational Fluid Dynamics, 20–25 July 2014, Barcelona, Spain, 1–11, https://congress.cimne.com/iacmeccomas2014/admin/files/filePaper/p3020.pdf (last access: 12 January 2023), 2014. a
Wu, L., Tiso, P., Tatsis, K., Chatzi, E., and van Keulen, F.: A modal derivatives enhanced Rubin substructuring method for geometrically nonlinear multibody systems, Multibody Syst. Dynam., 45, 57–85, 2018. a