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

Authors: Ilia K. Marchevsky; Valeria V. Puzikova
Journal: International Journal of Rotating Machinery
 (2017)

---

## Abstract

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 the modified LS-STAG method, the flow past autorotating Savonius rotor with two blades was simulated atRe=1.96·105.


---

## Body

## 1. 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–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” weak-coupling 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 level-set 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.


## 2. 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:(1)∇·v=0,∂v∂t+v·∇v=1ρ∇p+νΔv.Here v is the velocity vector, p is the pressure, t 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:(2)vinlet=v∞,∂v∂noutlet=0,∂p∂ninlet&outlet=0.Here v∞ 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:(3)vairfoil=vib,∂p∂nairfoil=0.Here vib is the velocity of the immersed boundary.To simulate the rotation of wind turbine rotors, the following dynamics equation is being solved:(4)Iα¨+kα˙=Mflow.Here α is the rotation angle of the rotor, I is the polar inertia moment of the rotor, k is the viscous friction coefficient in the axis, and Mflow is the aerodynamic torque. Two-dimensional problem is considered, so Mflow=Mz (Figure 1).Figure 1
The airfoil of irregular shape with one rotational degree of freedom and schematic view of vortex wake behind it.


## 3. Numerical Method

The Cartesian mesh with cellsΩi,j=xi-1,xi×yj-1,yj is introduced in the rectangular computational domain. It is considered that Γi,j is the face of Ωi,j cell and xi,jc=xic,yjc is its center. Unknown components ui,j and vi,j of velocity vector v are computed in the middle of fluid parts of the cell faces. These points are the centers of control volumes Ωi,ju=xic,xi+1c×yj-1,yj and Ωi,jv=xi-1,xi×yjc,yj+1c with faces Γi,ju and Γi,jv, respectively (Figure 2). The main idea of this numerical method is described in our previous papers, for example, in [14].Figure 2
Staggered arrangement of the variables on the modified LS-STAG mesh.Cells which the immersed boundary intersects are the so-called cut-cells (Figures3 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 Ωi,j. Locations of this segment’s endpoints are defined by a linear interpolation of the variable φi,j=φxi,yj.Figure 3
Example of a cut-cell.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.
(a)
(b)
(c)
(d)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 (Figures5 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.Figure 5
Schematic of the two-bucket Savonius rotor.Figure 6
Flow past autorotating Savonius rotor with two blades.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 airfoil at positionα=0. Then, the level-set function φn in point xi,yj can be computed as φnxi,yj=φ0xi~,yj~ for the rotor rotated by angle αn counterclockwise. Here φ0 is the level-set function built for the rotor surface line at position α=0 by using the Bézier curve and xi~,yj~ is the image of the point xi,yj after clockwise rotation on the angle αn. 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ϑi,ju and ϑi,jv are introduced. They take values in interval 0,1 and represent the fluid parts of the east and north faces of Γi,j, respectively. One-dimensional linear interpolation of φxi,y on the segment yj-1,yj and φx,yj on the segment xi-1,xi is used for the cell-face fraction ratios computing:(5)ϑi,ju=min⁡φi,j-1,φi,jmin⁡φi,j-1,φi,j-max⁡φi,j-1,φi,j,ϑi,jv=min⁡φi-1,j,φi,jmin⁡φi-1,j,φi,j-max⁡φi-1,j,φi,j.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:(6)Iωn+1-ωnΔt+kωn=Mzn.The aerodynamic torque at the n-th time step can be computed as follows:(7)Mzn=∑Cut-cellsΩi,jibxiC-XCFyahi,jn-yjC-YCFxahi,jn.Here xiC,yjC are coordinates of the center of Ωi,j cell, XC,YC are coordinates of the point around which the airfoil rotates, and Fxahi,jn and Fyahi,jn are the drag and lift forces acting on the solid part of the cut-cell Ωi,j at the n-th time step:(8)Fxahi,j=nxΔSi,jibpi,j-ν∂u∂xi,j-νQuadi,jib∂u∂ye→y·n→,Fyahi,j=-νQuadi,jib∂v∂xe→x·n→+nyΔSi,jibpi,j-ν∂v∂yi,j.Here nxΔSi,jib=ϑi-1,ju-ϑi,juΔyj, nyΔSi,jib=ϑi,j-1v-ϑi,jvΔxi, Δyj=yj-yj-1, Δxi=xi-xi-1, and Quadi,jib 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:(9)αn+1=αn+ωnΔt,vi,jib,n+1=-yjC-YCyjC-YC·ωnxiC-XC,xiC-XCxiC-XC·ωnyjC-YC. 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 (Figure7). For this purpose, a second-order Butterworth low-pass filter [16] is designed.Figure 7
The time dependency for the torqueMflow=Mzt (on the coarse grid 272×292 with time step Δt=0.0001) and filtered time dependency M¯flow=M¯zt for the torque.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:(10)Hs=1s-i-1/2εps+i+1/2εp=εp1+2εps+εps2.Here s∈C, i=-1, εp=10Rp/10-1, and Rp 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, H0=εp. Therefore the following transfer function corresponds to the normalized Butterworth second-order filter:(11)Hnorms=11+2εps+εps2.In order to control the filter cutoff frequency, it is necessary to use the following transfer function:(12)HLPs=HnormsΩx.The sampling frequency of the torque is equal tofd=1/Δt in the numerical simulation (Δt is a time discretization step). It is necessary that the filter cut-off frequency is equal to fs=5 Hz and the suppression level on the cut-off frequency is equal to Rs=3 dB. So, the following condition imposed on the required filter frequency response function HLPi·ω (Figure 8) must be satisfied for the required filter transfer function HLPs:(13)HLPi·Ωs=γ=10-Rs/20,Ωs=tan⁡πfsfd=tan⁡πfsΔt.Figure 8
The required filter frequency response function.Solution of (13) leads to the following result:(14)11+εp/Ωx4Ωs4=γ⟺εpΩx=1-γ24γΩs=ξ.Thus, the desired filter transfer function can be written as follows:(15)HLPs=11+2ξs+ξs2.To obtain the corresponding digital filter coefficients, there is a need to use the bilinear transformation [16]:(16)HLP1-z-11+z-1=1+2z-1+z-21+2ξ+ξ2+21-ξ2z-1+1-2ξ+ξ2z-2.Here z∈R. It is known that(17)HLP1-z-11+z-1=b0+b1z-1+b2z-21+a1z-1+a2z-2for 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:(18)b0=b2=11+2ξ+ξ2,b1=21+2ξ+ξ2,a1=21-ξ21+2ξ+ξ2,a2=1-2ξ+ξ21+2ξ+ξ2.So, the filtered value of the torque acting on the airfoil from the flow at the time tn+1=n+1Δt is(19)M¯zn+1=b0·Mzn+1+b1·Mzn+b2·Mzn-1-a1·M¯zn+a2·M¯zn-1.


## 4. Numerical Experiments

We considered autorotation of Savonius rotor with two blades (Figure5) at Re=1.96·105. The radius of rotation R was measured from the axis of rotation to the outer edge of the buckets. The turbine-swept area As is equal to πR2. Reynolds number is based on rotor diameter. The following parameters were used in the simulation: (20)V∞=1.0,R=D2=1.0,D~=0.1D,h=0.05R,k=0.0,I=10.0,α0=α0=π18.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.The dimensionless average rotor angular velocity is in the following range:(21)ω¯=0.4,…,1.6. 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(22)X∞=Rω¯V∞=ω¯;torque coefficient value C¯Q is obtained by averaging of nonstationary dependency CQt over time, where(23)CQt=2M¯ztρV∞2RAs=2M¯ztπ.Similarly, the power coefficient C¯P is obtained by averaging of nonstationary dependency CPt over time, where(24)CPt=2M¯ztω¯ρV∞3As=2M¯ztω¯π.A number of computations have been performed on nonuniform grid544×496. The uniform mesh block with space discretization step Δh=D/128 was used in the vicinity of the rotor. Time discretization step was equal to Δt=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 coefficientC¯P 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 CPt and computed value of the corresponding power coefficient C¯P at Re=1.96·105 and ω¯=0.9 is shown in Figure 10.Figure 9
Comparison of computed power coefficientsC¯P for a two-bucket Savonius rotor at Re=1.96⋅105 with experimental data (Sheldahl et al. [17] at Re=3.9⋅105 and Shankar [18] at Re=1.96⋅105).Figure 10
Time evolution of the power coefficientCPt at Re=1.96⋅105 and ω¯=0.9. The red line corresponds to the average value C¯P.


## 5. 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 atRe=1.96·105.(v)
Computational results are in good qualitative agreement with the experimental data.



---

*Source: 6418108-2017-10-11.xml*