Numerical Simulation of Wind Turbine Rotors Autorotation by Using the Modified LS-STAG Immersed Boundary Method

A software package is developed for numerical simulation of wind turbine rotors autorotation by using the modified LS-STAG level-set/cut-cell immersed boundary method. The level-set function is used for immersed boundaries description. Algorithm of level-set function construction for complex-shaped airfoils, based on Bézier curves usage, is proposed. Also, algorithm for the level-set function recalculation at any time without reconstructing the Bézier curve for each new rotor position is described. The designed second-order Butterworth low-pass filter for aerodynamic torque filtration for simulations using coarse grids is presented. To verify themodified LS-STAGmethod, the flow past autorotating Savonius rotor with two blades was simulated at Re = 1.96 ⋅105.


Introduction
The rotor is the first element in the chain of functional elements of a wind turbine.Its aerodynamic and dynamic properties, therefore, have a decisive influence on the entire system in many respects.The designer faces the problem of finding the relationship between the actual shape of the rotor, for example, the number of rotor blades or the airfoil of its blades, and its aerodynamic properties.
To simulate rotor's dynamics [1][2][3][4][5][6] and, in particular, its autorotation, there is a need to solve coupled aeroelastic problems.Such problems are complicated for numerical solution, since it is necessary to take into account interference between the flow and moving immersed body.So, the aim of this research is to develop an efficient numerical method for simulation of wind turbine rotors autorotation.
In case of massive body, when its average density is significantly higher than density of the flow, coupled aeroelastic problems can be solved using "step-by-step" weakcoupling numerical algorithm, firstly simulating flow around the body, which moves according to known parameters, and then computing the dynamics of the body under known hydrodynamic loads.Such case is considered in this research.
Immersed boundary methods [7] are suitable for numerical simulation in coupled aeroelastic problems, since they do not require a coincidence of cell edges and boundaries of the computational domain and allow solving problems when domain shape is irregular or it changes significantly in the simulation process due to aeroelastic body motion.The main advantage of these methods is that the mesh should not be reconstructed at every time step.
In the present study, the LS-STAG cut-cell immersed boundary method [8] is used for rotors autorotation simulation.This method permits solving problems on the Cartesian grid.The immersed boundary is represented with the levelset function [9].Linear systems resulting from the LS-STAG discretization of the Navier-Stokes or Reynolds-averaged Navier-Stokes equations are solved using the BiCGStab method [10] with the ILU and Multigrid [11,12] preconditioning.An original algorithm for the solver cost-coefficient estimation [13] is used for the optimal parameters of the multigrid preconditioner choice.

Governing Equations
The problem is considered in 2D unsteady case when the flow around a rigid airfoil is assumed to be viscous and incompressible.The continuity and momentum equations are as follows: Here k is the velocity vector,  is the pressure,  is the time,  is the flow density, and ] is the flow kinematic viscosity.
The boundary conditions on the outer boundaries of the computational domain are as follows: Here k ∞ is the velocity vector on the inlet boundary and n is the unit outer normal vector.The boundary conditions on the surface line of the airfoil are no-slip conditions: Here k ib is the velocity of the immersed boundary.
To simulate the rotation of wind turbine rotors, the following dynamics equation is being solved: Here  is the rotation angle of the rotor,  is the polar inertia moment of the rotor,  is the viscous friction coefficient in the axis, and  flow is the aerodynamic torque.Two-dimensional problem is considered, so  flow =   (Figure 1).

Numerical Method
The Cartesian mesh with cells Ω , = ( −1 ,   ) × ( −1 ,   ) is introduced in the rectangular computational domain.It is considered that Γ , is the face of Ω , cell and x  , = (   ,    ) is its center.Unknown components  , and V , of velocity  vector k are computed in the middle of fluid parts of the cell faces.These points are the centers of control volumes Ω  , = (   ,   +1 ) × ( −1 ,   ) and Ω V , = ( −1 ,   ) × (   ,   +1 ) with faces Γ  , and Γ V , , respectively (Figure 2).The main idea of this numerical method is described in our previous papers, for example, in [14].
Cells which the immersed boundary intersects are the socalled cut-cells (Figures 3 and 4).These cells contain the solid part together with the liquid one.The level-set function  is used for immersed boundary Γ ib description.The boundary Γ ib is represented by a line segment on the cut-cell Ω , .Locations of this segment's endpoints are defined by a linear interpolation of the variable  , = (  ,   ).
A few notes should be mentioned about construction of the level-set function.The level-set function cannot be constructed analytically for rotor with complex shape, that is, Darrieus rotor, as for circular airfoil and for simple rotor shapes, for example, for Savonius rotor (Figures 5 and 6).For this reason, it is necessary to approximate the rotor surface line with a curve, which would make it possible to simulate both smooth parts of the boundary and the sharp edges.Moreover, it is desirable that the distance from an arbitrary point to the boundary can be calculated easily.An efficient approach for the level-set function construction for an airfoil of arbitrary shape, which corresponds to the mentioned requirements, is described in [15].The airfoil's surface line approximation by using Bézier curve is taken as a basis of the developed algorithm.
In order to avoid the Bézier curve reconstruction at every time step for the airfoil surface line, the following approach can be used.At the beginning of the computation, the Bézier curve and its derivative should be built for the rotor blade International Journal of Rotating Machinery airfoil at position  = 0.Then, the level-set function   in point (  ,   ) can be computed as   (  ,   ) =  0 (x  , ỹ ) for the rotor rotated by angle   counterclockwise.Here  0 is the level-set function built for the rotor surface line at position  = 0 by using the Bézier curve and (x  , ỹ ) is the image of the point (  ,   ) after clockwise rotation on the angle   .Thus, the level-set function can be recalculated at any time without reconstructing the Bézier curve for each new rotor position.
The cell-face fraction ratios   , and  V , are introduced.They take values in interval ( According to the concept of the LS-STAG method, normal Reynolds stress components are sampled on the base mesh (similar to pressure discretization) and shear ones are sampled in the upper right corners of the base mesh cells.
The time integration of the differential algebraic system resulting from continuity and momentum equations sampling in space is performed with a semi-implicit Euler scheme [8].Predictor step leads to discrete analogues of the Helmholtz equation for velocities prediction at the next time layer.Corrector step leads to discrete analogue of the Poisson equation for pressure correction.After computation of the flow variables, the dynamics equation for the airfoil motion (4) should be solved.
The rotor angular velocity is  = α .So difference analogue of (4) can be written in the following form: The aerodynamic torque at the -th time step can be computed as follows: Here and Quad ib , are the quadratures of the shear stress which depend on the cut-cells [8].
The value of the rotor angular velocity at the next time step is computed from (6).New rotor position and the immersed boundary velocity can be defined by using this value: Such approach allows reconstructing the level-set function and all matrices and the source terms required for the next time step.When investigating some complicated physical phenomenon, preliminary qualitative estimations for the considered construction behavior are very important.They allow finding domains for grid refinement, prediction of the structure's dynamic, and estimation of the CFL number.But large fluctuations in the aerodynamic loads can occur on the coarse grid.So the values of the torque acting on the airfoil should be filtered by the low-pass filter (Figure 7).For this purpose, a second-order Butterworth low-pass filter [16] is designed.
It is necessary to explain the reasons for choosing this filter.Filters with infinite impulse response are less expensive from a computational point of view than filters with finite impulse response.The frequency response of the Butterworth filter is maximally flat (i.e., has no ripples) in the passband and rolls off towards zero in the stopband [16].Compared with a Chebyshev type I/type II filter or an elliptic filter, the Butterworth filter has a slower roll-off and thus will require a higher order to implement a particular stopband specification, but Butterworth filters have a more linear phase response in the passband than Chebyshev type I/type II and elliptic filters can achieve [16].
The transfer function of the Butterworth second-order filter is as follows: Here  ∈ C,  = √ −1,   = √ 10   /10 − 1, and   is a distortion level in the passband.As practice shows, filters of higher order can lead to the appearance of numerical instability in the filtered signal.For this function, (0) =   .Therefore the following transfer function corresponds to the normalized Butterworth second-order filter: In order to control the filter cutoff frequency, it is necessary to use the following transfer function: The sampling frequency of the torque is equal to   = 1/Δ in the numerical simulation (Δ is a time discretization step).It is necessary that the filter cut-off frequency is equal to   = 5 Hz and the suppression level on the cut-off frequency is equal to   = 3 dB.So, the following condition imposed on the required filter frequency response function | LP ( ⋅ )| (Figure 8) must be satisfied for the required filter transfer function  LP (): Solution of (13) leads to the following result: Thus, the desired filter transfer function can be written as follows: To obtain the corresponding digital filter coefficients, there is a need to use the bilinear transformation [16]: Here  ∈ R. It is known that for the second-order low-pass filter with infinite impulse response.Therefore, it can be obtained from formulae (16) and ( 17) that the designed digital filter coefficients are the following: So, the filtered value of the torque acting on the airfoil from the flow at the time  +1 = ( + 1)Δ is

Numerical Experiments
We considered autorotation of Savonius rotor with two blades (Figure 5) at Re = 1.96 ⋅ 10 5 .The radius of rotation  was measured from the axis of rotation to the outer edge of the buckets.The turbine-swept area   is equal to  2 .Reynolds number is based on rotor diameter.The following parameters were used in the simulation: (20) It should be noted that the above-proposed algorithm for level-set reconstruction can be easy to apply for rotors of other shapes, that is, Darrieus rotor.
Averaging was carried out over 16 dimensionless time units.This time is enough for the rotor to make full turn with the smallest value of average rotor angular velocity  = 0.4.At the chosen values of the parameters, tip speed ratio is equal to torque coefficient value   is obtained by averaging of nonstationary dependency   () over time, where Similarly, the power coefficient   is obtained by averaging of nonstationary dependency   () over time, where A number of computations have been performed on nonuniform grid 544 × 496.The uniform mesh block with space discretization step Δℎ = /128 was used in the vicinity of the rotor.Time discretization step was equal to Δ = 10 −5 .Computations were performed on a server based on the Intel C610 platform using the Intel Xeon E5-1620 V3 4-core processor (3.5 GHz) with Hyper-Threading support (8 logical cores).The server is equipped with 16 GB of ECC DDR4-2133 RAM and two hard drives (2 TB), united in a RAID1 disk volume.This server is running Windows Server 2012 R2 operating system.To simulate 1 dimensionless time unit, about 24 hours are required.The computed values of the power coefficient   were compared with the experimental data [17,18] presented in [18].As can be seen from Figure 9, the computational results are in good agreement with experiment.An example of simulated nonstationary dependency   () and computed value of the corresponding power coefficient   at Re = 1.96 ⋅ 10 5 and  = 0.9 is shown in Figure 10.

Conclusions
We have the following conclusions (i) A software package is developed for numerical simulation of wind turbine rotors autorotation by using the modified LS-STAG level-set/cut-cell immersed boundary method.This software package can be used for rotors of other shapes (Savonius rotor, Darrieus rotor, etc.).
(ii) Algorithms for level-set function construction and recalculation are described.
(iii) A second-order Butterworth low-pass filter is designed for aerodynamic torque filtrating at simulation on the coarse grid.
(iv) Simulation of flow past autorotating Savonius rotor with two blades is considered at Re = 1.96 ⋅ 10 5 .
(v) Computational results are in good qualitative agreement with the experimental data.

Figure 1 :
Figure 1: The airfoil of irregular shape with one rotational degree of freedom and schematic view of vortex wake behind it.

Figure 2 :
Figure 2: Staggered arrangement of the variables on the modified LS-STAG mesh.

Figure 3 :
Figure 3: Example of a cut-cell.

Figure 4 :Figure 5 :
Figure 4: Locations of the variables discretization points in case of generic cells on the modified LS-STAG mesh: (a) Cartesian fluid cell; (b) north trapezoidal cell; (c) northwest pentagonal cell; (d) northwest triangle cell.

Here 2 Figure 7 :
Figure 7: The time dependency for the torque  flow =   () (on the coarse grid 272 × 292 with time step Δ = 0.0001) and filtered time dependency  flow =   () for the torque.

Figure 8 :
Figure 8: The required filter frequency response function.