48 KiB
PAPER • OPEN ACCESS
Section force calculation in flexible substructures modelled by wind turbine design tool Bladed
To cite this article: Erik Nim et al 2024 J. Phys.: Conf. Ser. 2767 052043
View the article online for updates and enhancements.
You may also like
-Correction of periodic displacement nonlinearities by two-wavelength
interferometry
Angus Bridges, Andrew Yacoot, Thomas Kissinger et al.
Odd Random Phase Electrochemical
Impedance Spectroscopy to Study the Corrosion Behavior of Hot Dip Zn and ZnAlloy Coated Steel Wires in Sodium
Chloride Solution
Gopal Ji, Lucía Fernández Macía, Bart Allaert et al.
-New evidence and impact of electron
transport non-linearities based on new perturbative inter-modulation analysis M. van Berkel, T. Kobayashi, H. Igami et al.
Section force calculation in flexible substructures modelled by wind turbine design tool Bladed
Erik \mathbf{Nim}^{1} , John Roadnight2, Hassan Moharram2, William Collier2 and Chr. Sigurd L. Jensen
1 DNV Denmark A/S, Ecopark, Bautavej 1A, 8210 Aarhus V, Denmark
2 DNV Services UK Limited, One Linear Park, Bristol, BS2 0PS, United Kingdom
E-mail: erik.nim@dnv.com
Abstract. Aeroelastic simulation tools are widely used to predict the coupled response for onshore, offshore and floating wind turbines. A key output of the analysis is the section forces in the flexible substructures, such as blades, towers, and floating foundations. Individual substructures are typically modelled as a single linear flexible body, which is a part of a multibody system used for modelling the complete wind turbine. However, in some floating foundation concepts, in particular those with pre-tensioned cables, geometric non-linearities in the substructure can have a significant effect. With the conventional method for calculating section forces, such non-linearities are not captured accurately. In this paper, a new equilibriumbased method for calculation of section forces is presented, as part of the wind turbine design tool Bladed. This method relies on equilibrium between the section forces and the applied loads. The section forces are determined in terms of Lagrange multipliers, which represent the internal constraint forces in the structural elements. This method is valid even when significant geometric non-linearities are present, with reduced order models, and for statically indeterminate substructures. A case study is presented demonstrating the self-consistency of the new method for the OC4 semi-submersible floating platform, modified to include pre-tensioned cables. A verification against ANSYS is also presented for some elementary cases.
1. Introduction and motivation
Aeroelastic simulation tools are widely used to predict the coupled response for onshore, offshore and floating wind turbines. In wind turbine design tool Bladed [1], the structural model is based on a multibody system approach, where the complete wind turbine structure is modelled as an assembly of rigid and flexible bodies representing typical turbine substructures such as blades and tubular towers, as well as complex offshore support structures such as jackets and floating foundations.
A key functionality for a flexible body is the calculation of interior section forces (also known as internal forces). In the conventional method for calculating section forces, the calculation is carried out in the undeflected configuration. This is a good assumption in many cases but can lead to inaccuracies in bending and torsion moments when geometrical non-linearities become more significant. The inaccuracy is observed most clearly when calculating section forces in floating foundations, which can have significant non-linearities introduced by certain features, e.g., pre-tensioned cables. The present work represents a novel contribution to address these inaccuracies in an efficient and consistent way.
This paper proposes a new equilibrium-based method for the calculation of section forces in flexible bodies. The theoretical approach is described in section 2. The implementation in Bladed is verified against ANSYS in section 3, and a case study is presented for a floating substructure with pre-tensioned cables in section 4.
2. Theory summary
In this section, a new method for calculating section forces in flexible bodies is described. The method is based on an equilibrium approach often used in ad-hoc manual analysis of simple problems, but it has also been used as a basis for more general methods applicable to both statically determinate and indeterminate space frames [2][3].
2.1. Modelling of flexible bodies
In Bladed, the flexible bodies are modelled as structural frames using a linear finite element method [4][5] combined with the multibody system theory described by Shabana [6]. A structural frame is modelled by conventional finite element technique using two-node Timoshenko beam elements with 12 degrees of freedom (DOFs) along with a few additional elements for modelling specific structural parts, such as rigid bodies and bar elements used for modelling cables. An explicit geometric stiffness model is introduced to account for the second-order effects of the applied loads [4], but this approach has some limitations and is generally not sufficiently accurate. A traditional Craig-Bampton type modal truncation method is introduced to reduce the total number of degrees of freedom of the flexible body and to filter out unwanted high-frequency modes [7][8].
It is well known that the introduction of a modal approach for system reduction in multibody systems introduces a problem for representing the section forces sufficiently accurately, and it has therefore been proposed to calculate the section forces in a post-processing step, e.g., using finite element codes [9][10]. In Bladed, the section forces in the flexible bodies are calculated in a separate step at the modal truncated state resulting from dynamic or static analyses. For dynamic analyses, this post-processing step must be performed after every time step, which means that computational performance is critical.
2.2. Conventional calculation of section forces with the finite element method
With the displacement-based finite element method, the section forces are calculated at the end nodes of the beam elements using the element stiffness matrix. Possible geometric constraints, specified by the user to constrain some displacement patterns such as the axial extension of tubular towers, are modelled by a set of linear constraint equations with an associated set of Lagrange multipliers.
The conventional linear method for calculating section forces does not include any second-order effects of the loads applied to the flexible body due to fact that the loads are applied to the reference state and does not account for any deflections. This introduces an inaccuracy in the calculated bending and torsion moments, which can introduce significant errors in some load components, e.g., the torsion moment along a highly flexible blade. A geometric stiffness matrix is added to the material stiffness matrix to include some second-order effects [4].
2.3. A new equilibrium-based method for section force calculation
The most simple and accurate method for the calculation of section forces in a flexible body at a known deformed state is derived by equilibrating the section forces with the applied loads. This method allows the calculation to be carried out using only the loading, the actual deflection state and the geometry as well as the element connectivity, without using any stiffness (or damping) properties. However, the method only applies to statically determinate flexible bodies such as blades and tubular towers. For indeterminate flexible bodies including most support structures, the section forces depend on the stiffness properties, which consequently must be included in the calculation.
The new equilibrium-based method employs the multibody system approach used by Bladed for modelling of the complete wind turbine structure [1][11]. With this approach, a flexible body is modelled as an assembly of rigid and flexible elements that are interconnected at a number of userdefined nodes. A subset of these nodes is selected as global nodes used for connecting the flexible body to the neighbouring structure, whereas the remaining interior nodes are only used internally.
In the following, it is assumed that the deformation of a flexible body element is described by a set of generalised strains, which for a two-node beam element without constraints is conveniently collected in a vector \pmb{\varepsilon}^{\mathrm{e}} . In case that any deformation mode, such as torsion or axial extension, is constrained, the number of generalised strains is reduced accordingly. The absolute position and orientation of a node is described by a position vector and a rotation matrix to enable description of finite rotations. It is important to note that these quantities are dependent variables that are determined by the generalised strains. Although it is not possible to describe the absolute position and orientation of the element nodes by a vector because finite rotations are not vector quantities, a variation, i.e., a small virtual change, of these quantities can be described by a common vector \delta\mathbf{x}_{\mathrm{r}}^{\mathrm{e}} , which for a two-node beam element is a 12-by-1 vector. The relation between a variation \delta\mathbf{x}_{\mathrm{r}}^{\mathrm{e}} of the nodal position and orientation and a variation \delta\pmb{\varepsilon}^{\mathrm{e}} of the generalised strains defines the geometric constraints for the element. These constraint relations are written in a linear form
\mathbf{C}_{\mathrm{r}}^{\mathrm{e}}\,\delta\mathbf{x}_{\mathrm{r}}^{\mathrm{e}}+\mathbf{C}_{\mathrm{\varepsilon}}^{\mathrm{e}}\,\delta\mathbf{\varepsilon}^{\mathrm{e}}=\mathbf{o}
where \mathbf{C}_{\mathrm{r}}^{\mathrm{e}} and \mathbf{C}_{\varepsilon}^{\mathrm{e}} are geometric constraint matrices that generally depend on the deformation of the element. The symbol 𝐨 introduced in (1) represents a null vector with matching dimension. In general, the constraint matrices must be derived for all types of elements based on the geometry and the internal kinematics of the element. The derivation of the expressions for the constraint matrices for the used beam element was done using the co-rotational approach based on six natural deformation modes [12].
The derivation of the dynamic equilibrium equations for a flexible body element is based on the principle of virtual work similar to what is done for the multibody system modelling the complete structure [1][11]. It turns out that the element contribution to the virtual work can be written in the simplified form
\delta W^{\mathrm{e}}=\left[\begin{array}{l}{\delta\mathbf{x}_{\mathrm{r}}^{\mathrm{e}}}\\ {\delta\mathbf{\varepsilon}^{\mathrm{e}}}\end{array}\right]^{T}\left(\left[\begin{array}{l}{\mathbf{C}_{\mathrm{r}}^{\mathrm{e}^{T}}}\\ {\mathbf{C}_{\varepsilon}^{\mathrm{e}T}}\end{array}\right]\lambda^{\mathrm{e}}-\left[\begin{array}{c}{\mathbf{p}_{\mathrm{r}}^{\mathrm{e}}+\mathbf{f}_{0}^{\mathrm{e}}}\\ {\mathbf{p}_{\varepsilon}^{\mathrm{e}}-\mathbf{K}_{\varepsilon\varepsilon}^{\mathrm{e}}\,\varepsilon^{\mathrm{e}}}\end{array}\right]\right)
where the vectors {\bf p}_{\mathrm{r}}^{\mathrm{e}} and {\bf p}_{\varepsilon}^{\mathrm{e}} represent the applied element loads including inertial loads conjugate to respectively nodal motion and generalised strain, while the stiffness matrix {\bf K}_{\varepsilon\varepsilon}^{\mathrm{e}} represent the elastic properties on the assumption of linear elastic material. In addition, the vector \lambda^{\mathrm{e}} represents the yet unknown Lagrange multipliers associated with the constraint matrices, while the vector \mathbf{f}_{0}^{\mathrm{e}} contains the element section forces, which originate from the connection to neighbouring elements. The element equilibrium equations corresponding to the virtual work expression (2) is found by requiring \delta W^{\mathrm{e}}=0 for any combination of the variations \delta\mathbf{x}_{\mathrm{r}}^{\mathrm{e}} and \delta\pmb{\varepsilon}^{\mathrm{e}} :
\left[\begin{array}{c}{\mathbf{C}_{\mathrm{r}}^{\mathrm{e}^{T}}}\\ {\mathbf{C}_{\mathrm{\varepsilon}}^{\mathrm{e}T}}\end{array}\right]\boldsymbol{\lambda}^{\mathrm{e}}=\left[\begin{array}{c}{\mathbf{p}_{\mathrm{r}}^{\mathrm{e}}+\mathbf{f}_{0}^{\mathrm{e}}}\\ {\mathbf{p}_{\mathrm{\varepsilon}}^{\mathrm{e}}-\mathbf{K}_{\mathrm{\varepsilon}\mathrm{\varepsilon}}^{\mathrm{e}}\varepsilon^{\mathrm{e}}}\end{array}\right]
The resulting system for the complete flexible body is derived with the conventional finite element technique by assembling the element equations. To this end, the variation of the nodal positions and orientations for the complete flexible body are collected in a common 6N. -by-1 vector \delta\mathbf{x}_{\mathrm{r}} , where N denotes the total number of nodes of the flexible body. In addition, all generalised strain components are collected in a common N_{\varepsilon} -by-1 strain vector \pmb{\varepsilon} , where N_{\varepsilon} denotes the total number of generalised strains, while the Lagrange multipliers are collected in a common N_{\mathsf{c}} -by-1 vector \lambda , where N_{\mathsf{c}} denotes the total number of geometric constraints.
The geometric constraints for the complete flexible body are found by assembling the element constraints (1):
\mathbf{C}_{\mathrm{r}}\,\delta\mathbf{x}_{\mathrm{r}}+\mathbf{C}_{\varepsilon}\,\delta\mathbf{\varepsilon}=\mathbf{o}
where \mathbf{C}_{\mathrm{r}} is the N_{\mathsf{c}} -by- .6N constraint matrix for nodal displacement and rotation, while \mathbf{C}_{\varepsilon} is the N_{\mathsf{c}} -by$N_{\varepsilon}$ constraint matrix for the generalised strains.
The virtual work for the complete flexible body simply becomes \delta W=\textstyle\sum_{\mathrm{e}}\delta W^{\mathrm{e}} , where \delta W^{\mathbf{e}} are the element virtual work contribution (2). It is straightforward to show that the resulting virtual work contribution can be written in the form
\delta W=\left[\begin{array}{c}{\delta\mathbf{x}_{\mathrm{r}}}\\ {\delta\mathbf{\varepsilon}}\end{array}\right]^{T}\left(\left[\begin{array}{c}{\mathbf{C}_{\mathrm{r}}^{\mathbf{\Gamma}}}\\ {\mathbf{C}_{\varepsilon}^{\mathbf{\Gamma}}\mathbf{\Gamma}}\end{array}\right]\lambda\mathbf{\varepsilon}-\left[\begin{array}{c}{\mathbf{p}_{\mathrm{r}}}\\ {\mathbf{p}_{\varepsilon}-\mathbf{K}_{\varepsilon\varepsilon}\mathbf{\varepsilon}}\end{array}\right]\right)
where the flexible body load vectors \mathbf{p_{r}} and \mathbf{p_{\varepsilon}} , stiffness matrix \mathbf{K}_{\varepsilon\varepsilon} , and Lagrange multipliers \lambda are found by assembling the corresponding element quantities. The flexible body equations corresponding to the virtual work contribution (5) becomes
\left[\begin{array}{c}{\mathbf{C}_{\mathrm{r}}^{\mathrm{\Delta}T}}\\ {\mathbf{C}_{\mathrm{\varepsilon}}^{\mathrm{\Delta}T}}\end{array}\right]\boldsymbol{\lambda}=\left[\begin{array}{c}{\mathbf{p}_{\mathrm{r}}}\\ {\mathbf{p}_{\mathrm{\varepsilon}}-\mathbf{K}_{\mathrm{\varepsilon}\varepsilon}\;\varepsilon}\end{array}\right]
2.3.1. Statically determinate flexible bodies
For statically determinate flexible bodies, the number of constraints equals the number of nodal DOFs, i.e. 6 times the number of nodes. In this case, the Lagrange multipliers \lambda for the flexible body are calculated from the linear system \mathbf{C}_{\mathrm{r}}^{\textit{T}}\hat{\mathbf{\alpha}}=\mathbf{p}_{\mathrm{r}} appearing from (6) with the assumption that the constraint matrix \mathbf{C}_{\mathrm{r}} is invertible. The extracted element Lagrange multipliers \lambda^{\mathrm{e}} are then used for calculating the resulting section forces at the element end nodes by the linear system {\bf f}_{0}^{\mathrm{e}}={\bf C}_{\mathrm{r}}^{\mathrm{e}^{T}}{\bf\lambda}^{\mathrm{e}}-{\bf p}_{\mathrm{r}}^{\mathrm{e}} derived from (3). With this approach, the section forces are determined in terms of the Lagrange multipliers, which represent the internal constraint forces in the elements.
It is important to note that the flexible body constraint matrix \mathbf{C_{\mathrm{r}}} in (4) and the element constraint matrix \mathbf{C}_{\mathrm{r}}^{\mathrm{e}} in (1) must be evaluated for the modal truncated deflection state used during time integration or the steady state analysis of the complete system. In particular, the modal truncation can exclude all modes in which case the flexible body appears to be rigid.
2.3.2. Statically indeterminate flexible bodies
For a statically indeterminate flexible body, the number of constraints is larger than the number of nodal coordinates, i.e., N_{\mathrm{c}}>6N , which means that it is not possible to calculate the section forces from equilibrium alone, as is the case for a statically determinate flexible body. The degree of static indeterminacy, in the following denoted as D_{s} , is commonly used for quantifying the number of redundant forces, which in the present context means that D_{\mathrm{s}}=N_{\mathrm{c}}-6N . It is clear that {{D}_{s}}=0 for a statically determinate flexible body, while D_{s}>0 in the statically indeterminate case.
The approach for calculating section forces in statically indeterminate flexible bodies is based on subdividing the N_{\mathsf{c}} geometric constraints for a flexible body into a set of 6N determinate constraints and the remaining D_{\mathrm{s}} indeterminate constraints. In the following, the Lagrange multipliers associated with the determinate constraints are collected in the 6N-by-1 vector \lambda_{\mathrm{r}} , while the Lagrange multipliers associated with the remaining indeterminate constraints are collected in the D_{\mathrm{s}} -by-1 vector \lambda_{\mathrm{e}} . The subdivision of the Lagrange multipliers and hence the constraint relations is conveniently described in terms of a permutation matrix [\mathbf{P}_{\mathrm{r}}\;\mathbf{P}_{\mathrm{e}}] in the form
\lambda=\left[\mathrm{\bf~P}_{\mathrm{r}}\;\mathrm{\bf~P}_{\mathrm{e}}\right]\left[\lambda_{\mathrm{r}}\right]=\mathrm{\bf~P}_{\mathrm{r}}\,\lambda_{\mathrm{r}}+\mathrm{\bf~P}_{\mathrm{e}}\,\lambda_{\mathrm{e}}
where \mathbf{P}_{\mathrm{r}} is a N_{\mathsf{c}} -by- .6N partial permutation matrix for identifying the determinate constraints while \mathbf{P_{e}} is a N_{\mathsf{c}} -by- \cdot D_{s} partial permutation matrix for identifying the indeterminate constraints. The corresponding subdivision of the constraint relations (4) results in the system
\begin{array}{r}{\mathbf{C}_{\mathrm{rr}}\,\delta\mathbf{x}_{\mathrm{r}}+\mathbf{C}_{\mathrm{r}\varepsilon}\,\delta\mathbf{\varepsilon}=\mathbf{o}}\\ {\mathbf{C}_{\mathrm{er}}\,\delta\mathbf{x}_{\mathrm{r}}+\mathbf{C}_{\mathrm{e}\varepsilon}\,\delta\mathbf{\varepsilon}=\mathbf{o}}\end{array}
where \mathbf{C}_{\mathrm{rr}}=\mathbf{P}_{\mathrm{r}}^{\,\,T}\mathbf{C}_{\mathrm{r}} and {\bf C}_{\mathrm{r}\mathrm{s}}={\bf P}_{\mathrm{r}}^{~T}{\bf C}_{\mathrm{s}} , while \mathbf{C}_{\mathrm{er}}=\mathbf{P}_{\mathrm{e}}^{\ \ T}\mathbf{C}_{\mathrm{r}} and {\bf C}_{\bf e\mathrm{s}}={\bf P}_{\bf e}^{\mathrm{~\tiny~T~}}{\bf C}_{\bf s} . With this subdivision, the virtual work expression (5) takes the form
\delta W=\left[\begin{array}{c}{\delta\mathbf{x_{r}}}\\ {\delta\varepsilon}\end{array}\right]^{T}\left(\left[\begin{array}{c c}{\mathbf{C_{rr}}^{T}}&{\mathbf{C_{er}}^{T}}\\ {\mathbf{C_{r\varepsilon}}^{T}}&{\mathbf{C_{e\varepsilon}}^{T}}\end{array}\right]\left[\begin{array}{c}{\lambda_{\mathrm{r}}}\\ {\lambda_{\mathrm{e}}}\end{array}\right]-\left[\begin{array}{c}{\mathbf{p_{r}}}\\ {\mathbf{p_{\varepsilon}}-\mathbf{K_{\varepsilon\varepsilon}}\,\varepsilon}\end{array}\right]\right)
In principle, the subdivision of the geometric constraints into determinate and indeterminate partitions corresponds to the introduction of a set of virtual structural cuts that makes an indeterminate flexible body statically determinate and enables the calculation of section forces by an equilibrium approach. It is assumed that the subdivision of the constraints is carried out in such a way that the square matrix \mathbf{C}_{\mathrm{rr}} is invertible, which ensures that the constraint equation (8) can be solved with respect to the variation \delta\mathbf{x}_{\mathrm{r}} of the nodal position and orientation. With this assumption, it follows that the subdivided constraint relations (8) and (9) can be written as
\delta\mathbf{x}_{\mathrm{r}}=\mathbf{T}_{\mathrm{r}\mathrm{{s}}}\delta\mathbf{\varepsilon}
\bigl(\mathbf{C}_{\mathrm{e}\varepsilon}+\mathbf{C}_{\mathrm{er}}\,\mathbf{T}_{\mathrm{r}\varepsilon}\bigr)\delta\mathbf{\varepsilon}=\mathbf{o}
where \mathbf{T}_{\mathbf{r}\mathbf{s}} is a transformation matrix that is calculated from the linear system \mathbf{C}_{\mathrm{rr}}\,\mathbf{T}_{\mathrm{r\varepsilon}}=-\mathbf{C}_{\mathrm{r}\mathrm{\varepsilon}} .
In order to eliminate the variation \delta\mathbf{x}_{\mathrm{r}} of the nodal position and orientation in the virtual work expression (10), the transformation (11), together with the trivial identity \delta\pmb{\varepsilon}=\delta\pmb{\varepsilon} , are written in the common linear form
\left[\begin{array}{l}{\delta\mathbf{x}_{\mathrm{r}}}\\ {\delta\mathbf{\varepsilon}}\end{array}\right]=\left[\begin{array}{l}{\mathbf{T}_{\mathrm{r}\varepsilon}}\\ {\mathbf{I}}\end{array}\right]\delta\mathbf{\varepsilon}
where the symbol 𝐈 represents the identity matrix with matching dimensions. With this transformation, the virtual work expression (10) can be written in terms of the generalised strain variations \delta\pmb{\varepsilon} as
\delta W=\delta\varepsilon^{T}\left(\left(\mathbf{C_{e\varepsilon}}+\mathbf{C_{er}}\,\mathbf{T_{r\varepsilon}}\right)^{T}\lambda_{\mathrm{e}}-\left(\mathbf{p_{\varepsilon}}+\mathbf{T_{r\varepsilon}}^{T}\mathbf{p_{r}}-\mathbf{K_{\varepsilon\varepsilon}}\varepsilon\right)\right)
The resulting system of equilibrium equations corresponding to the virtual work expression (14) is found by requiring \delta W=0 for any combination of the generalised strain variations \delta\pmb{\varepsilon} . It is convenient to write the resulting system together with the constraint relation (12) in terms of a linear system
\left[\begin{array}{c c}{\mathbf{K}_{\mathrm{{e}}\varepsilon}^{*}}&{\mathbf{C}_{\mathrm{e}\varepsilon}^{*\ {T}}}\\ {\mathbf{C}_{\mathrm{e}\varepsilon}^{*}}&{\mathbf{0}}\end{array}\right]\left[\begin{array}{l}{\pmb{\varepsilon}^{\prime}}\\ {\lambda_{\mathrm{e}}}\end{array}\right]=\left[\begin{array}{l}{\mathbf{p}_{\mathrm{{e}}}^{*}}\\ {\mathbf{o}}\end{array}\right]
where {\bf C}_{\mathrm{e}\varepsilon}^{\ast}={\bf C}_{\mathrm{e}\varepsilon}+{\bf C}_{\mathrm{er}}{\bf T}_{\mathrm{r}\varepsilon} and \mathbf{p}_{\mathrm{{z}}}^{\ast}=\mathbf{p}_{\mathrm{{z}}}+\mathbf{T}_{\mathrm{{r}{\boldsymbol{\mathrm{\varepsilon}}}}}^{\;\;\;T}\mathbf{p}_{\mathrm{{r}}} , while \mathbf{K}_{\varepsilon\varepsilon}^{*}=\mathbf{K}_{\varepsilon\varepsilon} for notational consistency. This linear system enables the calculation of the indeterminate Lagrange multipliers \lambda_{\mathrm{e}} as well as the generalised strains \pmb\varepsilon^{\prime} corresponding to the applied loads {\bf p}_{\varepsilon}^{*} . Here the symbol \pmb\varepsilon^{\prime} is introduced to emphasise that the generalised strains calculated from the linear system (15) define a refined displacement field rather than the modal truncated displacement field used in the static or the dynamic analysis of the complete wind turbine. In particular, this approach allows an accurate calculation of the section forces for simulations with a low number of mode shapes (possibly none) for describing elastic deformations of the flexible body component.
In principle, the linear system (15) can be solved directly using Gaussian elimination, but it turns out that the augmented stiffness matrix (i.e., the coefficient matrix) on the left-hand side is ill-conditioned. A numerically stable solution method was derived based on the constraint elimination method [5], where the generalised strains \pmb\varepsilon^{\prime} are subdivided into a set of strains to be eliminated and a set of remaining independent strains. With this approach, the linear system (15) is transformed into a well-conditioned system, where the only unknowns are the independent generalised strains.
The determinate Lagrange multipliers \lambda_{\mathrm{r}} corresponding to the vector \lambda_{\mathrm{e}} of indeterminate Lagrange multipliers is calculated from the linear system {{\bf{C}}_{\mathrm{{rr}}}}^{T}{\lambda_{\mathrm{{r}}}}={{\bf{p}}_{\mathrm{{r}}}}-{{\bf{C}}_{\mathrm{{er}}}}^{T}{\lambda_{\mathrm{{e}}}} that is found from the virtual work expression (10) by requiring \delta W=0 for any combination of the variations \delta\mathbf{x}_{\mathrm{r}} . The complete vector \lambda of Lagrange multipliers can then be assembled by (7), which means that the section forces can be calculated similarly to what is done in the statically determinate case.
3. Verification against ANSYS
Figure 1 illustrates the general form of the elementary test models that are used for comparison between the section forces computed by Bladed and ANSYS. The models are statically indeterminate structural frames comprising of 4 cylindrical beam members interconnected at 4 vertices.
Figure 1. General definition of elementary verification model. Note that S1 to S5 locations for section force extraction are coincident with member ends.
As indicated in Figure 1, two different model configurations (denoted as 1 and 2) have been considered. In Model 2, the right and left vertices (V2 and V4) are joined by the addition of a bar element with no rotational degrees of freedom at each end.
All members are assigned tubular sectional properties with the outer diameter of 4\;\mathrm{m} and 20\;\mathrm{mm} thickness, with steel material with elastic modulus E=200 GPa, Poisson’s ratio \mathbf{0}=0.3 , shear modulus G=76.9\;\mathrm{GPa} .
Three different load cases are solved in ANSYS and Bladed, as listed in Table 1, representing different combinations of concentrated point loads and distributed line loads applied to the models.
The Bladed load cases are individually solved as static non-linear analyses based on the Bladed multibody system in a separate application. All load cases in ANSYS are individually solved as static non-linear analyses, with loads applied over a single load step.
Bladed is solved with the geometric non-linearity option of ‘Support Structure Geometric Stiffness Model’ set as ‘Axial loads only’ [1]. Geometric non-linearity is switched on in ANSYS using the ‘NLGEOM’ option [13].
Nodal deflections are extracted at each of the vertices V1 to V4. Section forces and moments are extracted from the member section locations labelled S1 to S5 in Figure 1. These section locations are located immediately adjacent to the vertices.
Table 1. Definition of load cases. All loads are given in the global frame
<html>| Model | Load case | Point loads on V2 (kN) (Right vertex) | Point loads on V3 (kN) (Top vertex) | Line loads on Mbr2 (kN/m) (Upper right member) | Line loads on Mbr3 (kN/m) (Upper left member) | ||||||||
| Px | Py | Pz | Px | Py | Pz | Lx | Ly | Lz | Lx | Ly | Lz | ||
| 1 | 1A | 50,000 | 5000 | 500 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |
| 1 | 1B | 0 | 0 | 0 | 100,000 | 10,000 | 0 | 1000 | 1 | 1 | -1000 | 0 -1 | -1 |
| 2 | 2A | 0 | 0 | 0 | 0 | 0 | 0 | -1000 | 0 | -1 | -1000 | 0 | 1 |
3.1. Beam element discretisation
Each frame member of the verification models is defined in Bladed using a single two-node beam element. The outer frame members in the ANSYS models are meshed with three-node BEAM189 elements [13] with two element divisions per frame member. For Model 2, there is a horizontal bar element between the left and right vertices. This bar is modelled using a LINK180 element in ANSYS [13]. This element has three translational degrees of freedom at each end node, and no rotational stiffness, so it does not transfer any moment loading.
The choice of two BEAM189 element divisions per frame member in ANSYS is made following a mesh sensitivity study, which showed small differences between force and moment outputs for 1 and 2 element divisions, but negligible difference between results for 2 and 10 element divisions per frame member.
3.2. Verification results
Relative differences between the Bladed and ANSYS computed section forces, moments, displacements and rotations are shown in Table 2 for load case 1A, Table 3 for load case 1B and Table 4 for load case 2A. The difference in the peak section force values is less than 0.3\% for all three load cases and differences in the peak section moments between 0.0\% and 1.4\% . Where more significant percentage errors are shown for individual force and moment components, this is considered acceptable where the absolute values of these components are orders of magnitude smaller than the peak values.
Differences in the Bladed and ANSYS peak displacements and rotations for load case 1A agree to within 0.7\% and 0.6\% , respectively. For vertices where the absolute value of displacement and rotation are smaller than the peak, but not insignificant, differences are up to 5.6\% and 7.0\% , respectively. For load case 1B, differences at the peak displacement and rotations are 4.5\% and 2.9\% , respectively. Finally for load case 2A, the difference for the peak displacement is 0.7\% and 1.7\% for peak rotation.
Some discrepancies in predictions of nodal displacements and rotation are expected due to differences in the formulations for geometric non-linearity for the two software. Most importantly for the present study, however, the new method for calculation of section forces and moments in Bladed is shown to give very good agreement with ANSYS.
Table 2. Main results for load case 1A given in terms of displacements/rotations and section forces/moments at four different sections, i.e., S1 at the support, S2 at the right vertex, S3 at the left vertex, and S4 at the top. All results are given in the global frame.
<html>| Displacements s(m) | Rotations (deg) | Force components (kN) | Moment components (kNm) | ||||||||||
| Tool | Ux | Uy | Uz | x | y | z | Fx | Fy | Fz | Mx | My | Mz | |
| S1 | ANSYS | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -55,261 | -14,272 | -498.2 | -11,461 | 10,495 | 762,603 |
| Bladed | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -55,247 | -14,264 | -497.8 | -11,471 | 10,573 | 766,338 | |
| Rel. diff. | - | - | -0.0% | -0.1% | -0.1% | +0.1% | +0.7% | +0.5% | |||||
| S2 | ANSYS | 2.52 | -2.72 | 0.09 | 0.13 | -0.10 | -5.50 | -5,261 | -9,272 | 1.7 | 894 | -912 | -281,053 |
| Bladed | 2.66 | -2.60 | 0.09 | 0.12 | -0.11 | -5.53 | -5,247 | -9,264 | 2.2 | 940 | -932 | -281,404 | |
| Rel.diff. | +5.6% | +4.2% | +0.8% | +7% | -7% | -0.6% | -0.3% | -0.1% | +32% | +5% | +2% | +0.1% | |
| S3 | ANSYS | -2.00 | -1.86 | 0.00 | 0.06 | -0.05 | -5.51 | -5,261 | -9,273 | 1.7 | 1,835 | -1,395 | -11,852 |
| Bladed | -1.93 | -1.93 | 0.00 | 0.06 | -0.06 | -5.51 | -5,247 | -9,264 | 2.2 | 1,860 | -1,390 | -12,612 | |
| Rel.diff. | -3.6% | +4.1% | +0.4% | +21% | +0.0% | -0.3% | -0.1% | +25% | +1.3% | -0.4% | +6.4% | ||
| S4 | ANSYS | 0.50 | -4.60 | 0.07 | 0.10 | -0.07 | -3.83 | -5,261 | -9,272 | 1.8 | 1,126 | -962 | 146,111 |
| Bladed | 0.69 | -4.57 | 0.08 | 0.09 | -0.08 | -3.85 | -5,247 | -9,264 | 2.2 | 1,080 | -910 | 146,020 | |
| Rel.diff. | +39% | -0.7% | +15% | -3.3% | +17% | +0.3% | -0.3% | -0.1% | +22% | -4.1% | -5.4% | -0.1% | |
Table 3. Main results for load case 1B given in terms of displacements/rotations and section forces/moments at four different sections, i.e., S1 at the support, S2 at the right vertex, S3 at the left vertex, and S4 at the top. All results are given in the global frame.
<html>| Displacements (m) | Rotations (deg) | Force components (kN) | Moment components (kNm) | ||||||||||
| Tool | Ux | Uy | Uz | x | y | z | Fx | Fy | Fz | Mx | My | Mz | |
| ANSYS | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -68,766 | -4,766 | -20 | -585 | 378 | 922,425 | |
| S1 | Bladed Rel. diff. | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -68,945 | -4,827 | -19 | -587 | 384 | 926,881 |
| - | - | +0.3% | +1.3% | -3.6% | +0.4% | +1.4% | +0.5% | ||||||
| S2 | ANSYS | 2.20 | -2.32 | 0.004 | 0.009 | -0.003 | -1.08 | 68,776 | 4,769 | 20 | 61 | -25 | 827,599 |
| Bladed | 2.29 | -2.23 | 0.003 | 0.006 | -0.002 | -1.07 | 68,945 | 4,827 | 19 | 77 | 10 | 831,745 | |
| Rel.diff. | +4.5% | +3.9% | +22% | +29% | +16% | +0.7% | +0.2% | +1.2% | -3.6% | +27% | -140% | +0.5% | |
| ANSYS | 0.92 | 0.94 | -0.005 | 0.009 | 0.004 | -1.10 | 31,238 | 5,233 | -20 | 33 | 13 | -580,803 | |
| S3 | Bladed | 0.94 | 0.92 | -0.004 | 0.007 | 0.004 | -1.08 | 31,055 | 5,173 | -19 | 55 | -10 | -577,186 |
| Rel.diff. | +1.8% | -2.1% | -20% | -26% | -14% | -2.0% | -0.6% | -1.1% | -3.6% | +70% | -179% | -0.6% | |
| ANSYS | 3.12 | -1.39 | -0.001 | 0.012 | 0.001 | -2.00 | 73,666 | 5,279 | 23 | 12 | -150 | 698,553 | |
| S4 | Bladed | 3.24 | -1.33 | -0.001 | 0.010 | 0.001 | -2.06 | 73,481 | 5,216 | 23 | 10 | -110 | 703,899 |
| Rel.diff.+3.7% | -4.4% | -30% | -22% | -20% | +2.9% | -0.3% | -1.2% | +3.1% | -17% | -26% | +0.8% | ||
Table 4. Main results for load case 2A given in terms of displacements/rotations and section forces/moments at five different sections, i.e., S1 at the support, S2 at the right vertex, S3 at the left vertex, S4 at the top and S5 in the bar element at the right vertex. All results are given in the global frame.
<html>| Displacements (m) | Rotations (deg) | Force components (kN) | Moment components (kNm) | ||||||||||
| Tool | Ux | Uy | Uz | x | @y | z | Fx | Fy | Fz | Mx | My | Mz | |
| ANSYS | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 42,426 | 38,066 | 12.23 | 635.2 | -471.4 | -55,708 | |
| S1 | Bladed Rel. diff. | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 42,426 -0.0% | 38,108 | 12.08 | 636.6 | -471.5 | -55,287 -0.8% |
| +0.1% | -1.3% | +0.2% | +0.0% | ||||||||||
| S2 | ANSYS Bladed | -0.105 -0.105 | 0.036 | -0.005 | -0.010 | 0.005 | -0.335 | 42,427 | -22,279 | 22.63 | 73.5 | 109.7 | 80,646 79,775 |
| Rel.diff. | -0.2% | 0.036 +0.2% | -0.005 -1.0% | -0.010 | 0.005 | -0.329 | 42,426 -0.0% | -22,214 -0.3% | 22.43 | 77.3 | 108.3 -1.2% | -1.1% | |
| -1.9% | -0.4% | +1.7% | -0.9% | +5% | |||||||||
| S3 | ANSYS | -0.105 | -0.036 | 0.005 | -0.010 | -0.005 | 0.335 | -42,426 | 38,065 | 12 | -74 | 110 | 80,643 |
| Bladed | -0.105 | -0.036 | 0.005 | -0.010 | -0.005 | 0.329 | -42,426 | 38,108 | 12 | -77 | 108 | 79,775 | |
| Rel. diff. | -0.2% | -0.2% | +1.3% | +1.9% | +0.4% | -1.7% | +0.0% | +0.1% | -1.2% | +4.4% | -1.7% | -1.1% | |
| S4 | ANSYS | -0.179 | 0.000 | 0.000 | -0.014 | 0.000 | 0.000 | 0 | -22,278 | -20.13 | 0 | 54.6 | 108,353 |
| Bladed | -0.178 | 0.000 | 0.000 | -0.014 | 0.000 | 0.000 | 0 | -22,214 | -19.99 | 0 | 49.8 | 106,802 | |
| Rel. diff. | -0.7% | +0.4% | -0.3% | -0.7% | -8.7% | -1.4% | |||||||
| S5 | ANSYS | -0.105 | 0.036 | -0.005 | -0.010 | 0.005 | -0.335 | 0 | 60,344 | -10.23 | 0 | 0 | 0 |
| Bladed | -0.105 | 0.036 | -0.005 | -0.010 | 0.005 | -0.329 | 0 | 60,322 | -10.35 | 0 | 0 | 0 | |
| Rel. diff. | -0.2% | +0.2% | -1.0% | -1.0% | -0.4% | +1.7% | -0.0% | +1.2% | |||||
4. Case study: floating wind turbine support structure
The Offshore Code Comparison Collaboration (OC4) project aimed to improve the understanding and development of offshore wind energy technologies [14]. As part of the OC4 project, one of the designs examined was the DeepCwind floating semi-submersible system. This design is analysed at full scale to support the National Renewable Energy Laboratory's (NREL) offshore 5 MW baseline wind turbine [15]. The turbine serves as a representation of a large-scale, multi-megawatt structure. The DeepCwind floating semi-submersible system was selected for this case study as it is suitable for representing a statically indeterminate floating structure with multiple members. The hydrodynamic loads calculation is based on potential theory with viscous drag modelled by Morison’s equation to model the hydrodynamic damping [1].
In this example, the floating platform model has been modified to simulate a challenging scenario for calculating section forces, requiring the use of the new equilibrium-based method for section force calculation described in this paper. The unmodified OC4 DeepCwind validation model is shown in Figure 2(a). For this case study, the traditional cantilevering tower has been replaced with a slender tower supported by three pre-tensioned cables, as shown in Figure 2(b). The modified tower has a diameter of 3\textrm{m} and walls that are 19\,\mathrm{mm} thick. The cables are modelled using a massless material and are pre-tensioned to 13\;\mathrm{MN} in each cable. Additional point masses are added to the tower nodes to maintain the total mass and the mass centre of the original design. The tower height interface sections remain the same as in the original OC4 model configurations detailed in [14]. It is important to note that this configuration has not been validated through DeepCwind scale model testing or numerically in the OC4 project. Instead, it is presented in this case study as a conceptual illustration.
Figure 2: (a) OC4 DeepCwind floater design (b) modified design with cable support (c) modal reference node position at tower base.
In Bladed, the floating platform and tower are considered as a single flexible body, which can undergo rigid body motions and linear structural deformations. For the purpose of calculation of Craig Bampton mode shapes [7][8], it is necessary to specify a node to be treated as fixed in space in order to exclude rigid body modes. This so-called 'modal reference node' is chosen to be located at the tower base, as shown in Figure 2(c). If the applied loading used for calculating the section forces is not in equilibrium, a non-zero reaction force will be generated at the modal reference node. This, in turn, leads to a discontinuity in the member forces at the sections on either side of the node. The purpose of this example is to demonstrate this discontinuity and how the new equilibrium-based method can serve as a better basis for calculating section forces compared to the conventional method.
A simulation was carried out with the turbine idling at a wind speed of 10\;\mathrm{m/s} and no waves. The section force output location and coordinate system at the modal reference node are shown in Figure 2(c). Table 5 presents the bending moment about z on the element ends on either side of this node. It is observed that the conventional method shows a 9.7\% discontinuity in bending moment in the structure, whereas the equilibrium-based method shows no discontinuity. The latter is the expected result as no loads are applied at the modal reference node.
Table 5. Bending moment just above and below the modal reference node (MRN) of the floater.
<html>| Bending moment | |||
| Method | Location BelowMRN | about z (MNm) 0.964 | Rel. diff. (%) |
| Conventional | AboveMRN | 1.057 | 9.7% |
| Equilibrium | BelowMRN | 1.052 | 0.0% |
| Above MRN | 1.052 | ||
5. Conclusion
A new equilibrium-based approach for calculating internal section forces in flexible bodies has been developed and presented.
The method was used in the wind turbine design tool Bladed as an integrated part of the multibody system modelling the complete wind turbine structure. A separate method for calculating the section forces is mandatory for this application, where the deformation of the flexible bodies is described by a truncated modal expansion with a relatively low number of mode shapes (or none). The method is particularly relevant for statically indeterminate structures with significant geometrical non-linear effects, such as floating wind turbine support structure with pre-tensioned cables, where conventional methods can yield inaccuracies of the calculated section forces. In particular, this improves the confidence of practical wind turbine structural analyses in case a substructure is subject to large point loads, which may originate from pre-tensioned cables.
The new method has been verified by comparison against ANSYS, and the significance for section force calculation demonstrated through a case study based on a modified version of the DeepCwind floating platform.
References
[1] DNV AS 2023 Bladed 4.15 theory documentaion (online) https://dnvgldocs.azureedge.net/Bladed%204.15%20Theory%20Manual/index.html
[2] Pellegrino P and Calladine C R 1986 Matrix analysis of statically and kinematically indeterminate frameworks Int. J. Solids Structures 22 pp 409-28
[3] Pellegrino P 1992 Structural computations with the singular value decomposition of the equilibrium matrix Int. J. Solids Structures 30 pp 3025-35
[4] Przemieniecki J S 1968 Theory of Matrix Structural Analysis (New York: McGraw-Hill)
[5] Cook R D, Malkus D S and Plesha M E Conceps and Application of Finite Element Analysis, 3^{r d} edition (New York: John Wiley & Sons)
[6] Shabana A A 2014 Dynamics of Multibody Systems \mathcal{I}^{t h} edition (Cambridge: Cambridge University Press)
[7] Craig Jr R R and Bampton M C C 1968 Coupling of substructures for dynamic analysis AIAA Journal 6 pp 1313-19
[8] Craig Jr R R 2000 Coupling of substructures for dynamic analysis: an overview AIAA/ASME/ ASCE/AHS/ASC 41st structures, structural dynamics and materials conference (Reston, VA: American Institute of Aeronautics and Astronautics)
[9] Schwertassek R, Wallrapp O and Shabana A A 1999 Flexible Multibody Simulation and choice of shape functions Nonlinear Dynamics 20 pp 361-80
[10] Schwertassek R, Dombrowski, S V and Wallrapp O 1999 Modal representation of stress in flexible multibody simulation Nonlinear Dynamics 20 pp 381-99
[11] Meijaard J P 2022 An extended modelling technique with generalized strains for flexible multibody systems Multibody System Dynamics 57 pp 133-55
[12] Krenk S 2009 Non-linear Modeling and Analysis of Solids and Structures (Cambridge: Cambridge University Press)
[13] ANSYS, Inc. 2013 ANSYS Mechanical APDL Element Reference and Theory Reference manuals, Release 14.5.7 (Canonsburg, PA: SAS IP)
[14] Robertson A, et al 2014 Definition of the Semisubmersible Floating System for Phase II of OC4, NREL/TP-5000-60601 (Golden, CO: National Renewable Energy Laboratory)
[15] Jonkman J, Butterfield S, Musial W and Scott G 2009 Definition of a 5-MW Reference Wind Turbine for Offshore System Development, NREL/TP-500-38060 (Golden, CO: National Renewable Energy Laboratory)