A FLUID-STRUCTURE INTERACTION SCHEME APPLIED TO THE ANALYSIS AND DESIGN OF A WIND TURBINE GENERATOR SYSTEM – AN ENGINEERING SOLUTION

Nowadays, looking to reduce global warming and improve sustainability energy, wind turbine generators are a real alternative for clean power generation. To withstand wind forces, aerogenerators should be analyzed properly to ensure a good operation in the system lifetime. Several approaches can be established to perform analysis and design of this kind of structures. In this paper, an engineering solution considering fluid structure interaction is presented. Wind action is modelled thru a stabilized fluid flow formulation, while the structure (wind turbine) is solved with geometrically nonlinear shell elements with only translation degrees of freedom. In both cases, the finite element method is used to found a solution. Interaction between both formulations are capable to reproduce spin of turbine blades, including self weight forces and time history analysis to obtain stresses and acting forces for several operational conditions.


INTRODUCTION
The implementation of a coupled problem can be done using two different global strategies, which are the monolithic methods and the partitioned methods.In monolithic methods, the discretized fluid-structure system is solved together with the mesh movement system in a single iteration loop, leading to a very large system of nonlinear equations to be solved simultaneously.Some advantages of this method are that it ensures stability and convergence of the whole coupled problem.On the contrary, in simultaneous solution procedures the time step has to be equal for all subsystems, which may be inefficient if different time scales are presented for the problem.An important disadvantage is the considerably high computing time effort required to solve each algebraic system and sometimes the necessity to develop new software and solution methods for the coupling method.A monolithic approach to FSI is presented by Hübner et al. [18].
In partitioned methods, the application of existing appropriate and sophisticated solvers for either structural or fluid subsystems will continue to be used.These methods enjoy great popularity due to the simplified coupling procedure in many cases.Partitioned methods are divided into weak or loose coupling algorithms and strong or implicit coupling schemes.Weak algorithms are also known as staggered or explicit schemes.The major drawbacks of partitioned methods are lack of accuracy and stability problems, which sometimes may diverge from the solution.
Partitioned methods were introduced by Park and Felippa [29].The key idea for these methods is described in Felippa et al. [13].An interesting and simple example showing a complete description for partitioned methods can be found in Valdés [37].Partitioned solutions with staggered coupling algorithms are developed by Farhat et al. [12] to be used in aeroelastic wing problems.Strong coupling of partitioned algorithms are applied to large displacements 2D structural problems coupled to viscous incompressible fluids by Wall and Ramm [40] and Wall [38].They also applied the same method for a coupled fluid structure environment with an initially flat three-dimensional shell model as given in Wall and Ramm [41].Other large displacements structural problems interacting with incompressible fluids are detailed in Mok [26], Mok and Wall [25] and Tallec and Mouro [32].FSI with large displacements applied to wind problems is developed by Rossi [30], Wüchner and Bletzinger [43], Badia [2] and Wüchner [42].
More sophisticated developments on strong partitioned methods for FSI problems can be found in Steindorf [31], Matthies and Steindorf [22], [23], [24] and Tezduyar et al. [34].A study on strong coupling partitioned methods for FSI applied to hemodynamic problems can be found in Nobile [27], Causin et al. [4], and Fernández and Moubachir [14].Strong coupling of fluid-structure interaction including free surfaces is studied in Dettmer [9] and Wall et al. [39].Recently, a new approach based on Robin transmissions conditions for fluid-structure interaction problems is given in Badia et al. [1].

FORMULATION
The governing equations for the couple incompressible fluid-structure problem consist of the momentum equations together with the continuity equation.However the fluid and the structural parts of the domain must be treated differently.Then the problem is split into the fluid test functions over the fluid domain and the solid test function operating over the structural part.
For convenience, the boundary of the coupled problem is divided into the Dirichlet boundary for the fluid Γ D f and the Dirichlet boundary for the solid Γ D s , the Neumann boundary for the fluid Γ N f and the Neumann boundary for the solid Γ N s , and a common interface boundary Γ I between the fluid and the solid.Then the boundary of the couple problem is Also the coupled problem domain is divided into a solid part Ω s and a fluid part Ω f , where Ω = Ω s ∪ Ω f .In particular, solid displacements u s are given as a function of the solid motion x s as where X ∈ Ω s 0 with Ω s 0 being the reference solid domain.The fluid mesh motion is defined as an extension inside the fluid domain Ω f of the solid displacements u s located at the interface boundary Γ I .In other words, fluid displacements are found by where different forms to make this extension are i.e. the Laplacian method or the pseudoelastic structural technique, among others.Fluid mesh displacements yield a fluid mesh ve-locity, v(u f ) = v mesh , that has to be integrated into the fluid formulation.Since the Lagrangian solid displacements move the eulerian fluid domain then our fluid formulation has to be modified to include this movement.In this work we take into account the fluid movement using the Arbitrary Lagrangian Eulerian (ALE) formulation, as given in [3], [10].Before writing the continuous formulation of the couple fluid-structure interaction problem, the subspace test functions for the fluid, with a homogeneous Dirichlet boundary condition, are defined by and the subspace test functions for the structure are expressed as Note that the test functions for the fluid vanish at the interface and inside the solid subdomain, while the solid test functions are nonzero on the interface and extend inside the fluid subdomain.In this way, the kinematic continuity for the displacement and velocity field is imposed as Dirichlet boundary conditions on the fluid by the interface, and kinetic continuity for the traction is given as Neumann boundary conditions on the structure at the interface.
The structural problem is obtained by considering the space of the solid test functions on the momentum equations, yielding which shows the interface traction obtained directly from the momentum conservation equation and not considered as an additional independent equation.A detailed mathematical explanation can be found in [37], [32].
The fluid problem is obtained by substituting into the momentum equation the corresponding test functions, which vanish in the solid part, yielding after integrating by parts In this equation, mesh movement is taking into account with the convective term c j which is given by the difference v j − v mesh .This equation must be complemented by adequate boundary conditions dictated by the fluid space test functions.Since fluid test functions vanish at the interface, then displacement and velocity fields can only be governed by imposing the kinematic compatibility given by the Dirichlet conditions The continuity equation for the fluid part remains the same that for the eulerian fluid problem.Besides for the fluid formulation, a stabilization technique has to be added, which in this work is taken with the FIC formulation of Oñate [28].New development of the FIC formulation can be found in Lynga [21].
Among the wide possibilities to update meshes in the fluid field, the most common is the Laplacian method which can be found in [3].Other methods are based on the pseudostructural system, which can be done through the elastic spring analogy, i.e. see [11] and [8], or by solving the elasticity equation as a pseudo-elastic system, i.e. see [20], [3] and [5].In this work, two different ways for the mesh movement of the fluid subdomain are implemented: the Laplacian and the pseudo-elastic method.
In this way, the mesh is considered as another system.Therefore the fluid-structure interaction problem yields a coupled three-field system: the fluid, the solid and the mesh movement.

FINITE ELEMENT DISCRETIZATION FOR SHELLS
The discretization used for the rotation-free thin shell finite element is given for the total Lagrangian formulation.A detailed description of this rotation-free element is found in Flores and Oñate [15], [16], Valdés [37], and Valdés and Oñate [36].Here only the general idea and principal results are presented.It is important to remark that in case of studying an orthotropic shell, first a principal fiber orientation is needed, as given in Valdés [37].Fig. 1 shows a general shell mesh before the fiber orientation (left) and after the fiber orientation (rigth) where red-dashed lines show vectors of principal fiber orientation e d for an initially out-ofplane membrane structure together with its unit base vectors e i .A deep study of membranes using the fiber orientation methodology can be found in Valdés et al. [35].Bending effect for this rotation-free shell triangle element is given by the displacement field of one element and all nodes of adjacent elements, as shown in the patch of Fig. 2. Membrane forces for the rotation-free element are computed taking into account only the main element of the patch.
The path description is as follows: • Element number is inside a circle.
• Nodes of the main element (M) are numbered locally as 1, 2 and 3.
• Main element sides are defined by its local node opposite to the side.
• Adjacent elements are numbered with the number associated to the common side 1, 2 and 3. • The remaining nodes of the patch are numbered locally as 4, 5, and 6 corresponding to nodes on adjacent elements (1), ( 2) and ( 3) respectively.
In this work, a local coordinate system is given by the local fiber Cartesian base system as shown in Fig. 1, and given in Valdés [37].Then the base system for each finite element is given by the unit vectors e f ib 1 , e f ib 2 and the normal e f ib 3 .The choice for this local system allow us to compute shells with orthotropic material definitions.
The virtual internal work for shells can be expressed by where the Green-Lagrange strain tensor is given by where measures membrane strains and measures bending strains.
From the virtual internal work, internal forces for the membrane part in curvilinear coordinates can be written as Using the Voigt notation to express the internal forces in curvilinear coordinates, Eq. ( 13), yields where the strain-displacement matrix B B B memb I is given by Finally, this strain-displacement matrix is changed from curvilinear to Cartesian coordinates and then using the fiber orientation methodology, it is rotated to fiber direction yielding the membrane strain-displacement matrix B memb ab .
From the virtual internal work, forces for the bending part in curvilinear coordinates can be written as Using the Voigt notation to express the internal forces in curvilinear coordinates, Eq. ( 16), yields where the variation of the bending strain tensor in Voigt notation and expressed in Cartesian coordinates can be written explicitly as where nJ α is the normal of the each side of the main element of the patch, x J ,α are the Cartesian derivatives following the principal fiber orientation and xh ,α are the contravariant base vectors.Finally the internal forces for the rotation-free shell element are expressed by Now the complete structural problems is expressed as where M is the mass matrix and ü is the acceleration vector of the nodal variables.An implicit solution of Eq. ( 20) requires the derivative of the internal forces yielding a tangent matrix.In Flores and Oñate [15] the complete tangent matrix including material nonlinearity is developed for the isotropic rotation-free shell element.

FINITE ELEMENT DISCRETIZATION FOR FLUIDS
Finite element discretization of the incompressible flow equations is presented using the Galerkin-type weak form.Oñate [28] gives a full description of the incompressible fluid formulation using the finite calculus FIC stabilization technique and the fractional step method.A complete description for the fluid problem emerged from the continuum mechanics equations including finite element discretization, stabilization and time integration schemes is given in Valdés [37].From the conservation of linear momentum equation, given by Eq. ( 6), we can find after discretization where with c h being the convective velocity given by the difference between the spatial velocity v h and the mesh velocity v h mesh .From the continuity equation, the divergence of the velocity yields after discretization which is the incompressibility condition for fluids where Equations ( 21) and ( 26) describe the incompressible Navier-Stokes equations for fluids.Since the original works of Chorin [6] and Temam [33], fractional step methods for the incompressible Navier-Stokes equations have earned widespread popularity because of the computational efficiency given by the uncoupling of the pressure from the velocity field.A detailed stability analysis of fractional step methods for incompressible flows is given in Codina [7].
The easiest form to understand the development of the fractional step method is from the expression of the incompressible Navier-Stokes equations, where Eq. ( 21) can be split and yield the equivalent incompressible flow equations where ṽn+1 is an auxiliary velocity variable and the following essential approximation has been taken: K(c n+1 )v n+1 ≈ K(c n+1 )ṽ n+1 .From Eq. ( 29), v n+1 can be expressed in terms of ṽn+1 yielding Substituting this last equation into Eq.( 26) yields Now observe that G T M −1 G represent an approximation to the Laplacian operator, as mentioned in Codina [7], given by L IJ = ρ (∇N I , ∇N J ).
Finally the incompressible fluid flow equations to be solved using the θ f -family integration scheme and using the FIC stabilization technique, developed by Oñate [28], yield where with 1 2 ≤ θ f ≤ 1 for unconditionally stable implicit schemes and τ i being the intrinsic time parameters defined in Oñate [28].

FLUID-STRUCTURE INTERACTION
With the developed equations for the fluid solver F and the solid solver S, the problem consist on finding the appropriate method to solve the fluid-structure interaction (FSI) problem.Two different global methods for FSI problems can be used: monolithic methods and partitioned methods.In this work the last technique is described.
Numerical simulation of FSI problems is not only difficult because of the problems associated with the fluid or structure solution, but also because of the coupling interface between these two fields which sometimes represents another challenge.These difficulties depend strongly on the added mass effect introduced by the fluid over the structure, as given in Causin et al. [4].When the structure density ρ s is much larger than the fluid density ρ f , the added mass effect of the fluid is not significant and the problem can be solved with staggered partitioned schemes or strong coupling partitioned techniques as in the case of aeroelasticity.However when the structure and fluid densities are of the same order, as in the case of hemodynamics, the added mass effect of the fluid over the structure is very important and the coupling algorithm to be used must be a strong coupling partitioned scheme with relaxation, or even better techniques as the exact Newton or inexact Newton method for strong coupling problems, or the recently partitioned procedures based on Robin transmissions conditions given by Badia et al. [1].
To explain these methods, assume that the fluid unknowns are grouped together in the vector x and that the nonlinear iterative fluid solver is written as where y represent the displacements of the structure, which also define the current configuration of the interface.The nonlinear iterative solid solver is given by where x defines the nodal fluid velocities and pressure including the traction fluid forces at the interface for the structure.In this work, the Block Gauss Seidel method is used which consists in iterating the fluid and solid solvers independently as shown next where p and q the number of times that the solvers F and S are repeated respectively.Also for i = 1, x 0 n+1 = x n and y 0 n+1 = y n .Once the tolerance or maximum number of iterations is reached, the time step is advanced and the process is repeated.

STRONG COUPLING WITH RELAXATION
In this work, the strong coupling block Gauss-Seidel partitioned method with relaxation has been implemented.The structural solver S is now referred to as the computational solid dynamic (CSD) solver.The fluid solver F solves the fluid equations plus the movement of the mesh, yielding in a high-cost task from the computational point of view.Therefore the fluid solver F is split into the computational mesh dynamic (CMD) solver that moves the interior nodes of the finite element mesh of the fluid subdomain, and the computational fluid dynamic (CFD) solver that uses the fractional step method together with the ALE technique incorporated in the momentum equation.All solvers of each field use an implicit scheme.
There are several methods in the literature to accelerate the solution of the problem by means of relaxation, i.e.Mok [26].In this work the Aitken method is implemented, i.e.Irons and Tuck [19], which yields excellent solutions with simple modifications to the code.
In order to compute the coupled fluid-structure interaction problem, a unified algorithmic framework for the whole procedure is presented next.Considering known all values of solid, fluid and mesh at time step t n , the new step t n+1 is found following the simple steps given ahead: 1. Advance time step: t n+1 = t n + ∆t 2. Set iteration i = 1 3. Compute interface predictor displacement with one of the following methods (a) Structural predictor: Solve the CSD problem to find ûn+1 |Γ I ,i at the interface from u n+1 |Ω s ,i with a predicted external force given by: i. Pressure: • Order 1:  The Aitken method of relaxation is based on Aitken's acceleration method for vectors.The method can be easily implemented in any code which then ensures converge of the coupled problem for adequate time step parameters.The Aitken relaxation parameter is computed with the following algorithm, as given in Wall et al. [39].
(c) Compute Aitken optimal relaxation parameter wi = 1 − μn+1 i More sophisticated and computationally expensive methods, such as the gradient method, lead to solutions as good as the Aitken method for fluid-structure interaction problems.Additional computational cost for this technique is insignificant since only vector operations over the interface nodes are performed.

Wind Turbine Generator System
This work consists of several parts and aims to find the influence of local topography on the efficiency of wind turbines.The first part is to determine the pressure exerted by the wind on the turbine.Figure 3 shows the wind pressure over the turbine blades.This pressure is transferred to the structural solver as an external force, and then the wind turbine is solved with those forces.
Figure 4 shows the wind pressure over the turbine blades with the whole aerogenerator system.Movement of blades can be plot for its norm displacement, as shown in Figure 5.A complete cicle for rotation of the turbine blades is presented in Figure 6.As can be seen from these figures, more work needs to be done to find out how efficient wind generator systems are with different terrain conditions.Next parts of this work will be presented shortly.

CONCLUSIONS
Many fluid-structure interaction problems are solved with staggered coupling techniques, which are only acceptable for problems where the added mass effect do not have influence in the structure.However this kind of problem may become unstable for long time periods of study.To avoid this problem, strong coupling partitioned schemes are advised.
However when the added mass effect plays an important role in the structure, partitioned methods with block Newton schemes are an excellent choice.Another option is the block Jacobi or block Gauss-Seidel method with relaxation techniques.These last two methods are useless if they do not include the relaxation technique.In this work, the strong coupling block Gauss-Seidel partitioned method with a relaxation technique based on Aitken's method is implemented, yielding an excellent solution to the examples presented.
Finally, a detailed algorithm for the fluid-structure interaction problem using the strong coupling block Gauss-Seidel partitioned method is presented.In there, a relaxation technique with Aitken's method is used.However, other relaxation techniques can be added with minor modifications.

5 .
mesh velocities vn+1 i (b) CFD solver: Solve fluid i. Transfer vn+1 i to the fluid solver ii.Solve the CFD problem and find v n+1 i , p n+1 i iii.Compute fluid stress tensor at interface σ f |Γ I (c) CSD solver: Solve structure i. Transfer σ f |Γ I to the solid solver and compute structure forces t f ii.Solve the CSD problem and find u n+1 i+1 (d) Relaxation phase i. Compute optimal relaxation parameter wi via Aitken method ii.Relaxation of predicted interface position with ûn+1 |Γ I ,i+1 = (1 − wi )û n+1 |Γ I ,i + wi u n+1 |Γ I ,i+1 (e) Advance iteration: i = i + 1 (f) Check convergence.If reached, go to 5, else go to 4 Check time step.If end of time not reached, go to 1, else end of calculation

Figure 3 .
Figure 3. Wind Pressure on Turbine Blades