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 at .
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:Here 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 is the velocity vector on the inlet boundary and is the unit outer normal vector. The boundary conditions on the surface line of the airfoil are no-slip conditions:Here 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 is the aerodynamic torque. Two-dimensional problem is considered, so (Figure 1).
3. Numerical Method
The Cartesian mesh with cells is introduced in the rectangular computational domain. It is considered that is the face of cell and is its center. Unknown components and of velocity vector are computed in the middle of fluid parts of the cell faces. These points are the centers of control volumes and with faces and , 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 so-called 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 description. The boundary 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)
(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 (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 airfoil at position . Then, the level-set function in point can be computed as for the rotor rotated by angle counterclockwise. Here is the level-set function built for the rotor surface line at position by using the Bézier curve and 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 are introduced. They take values in interval and represent the fluid parts of the east and north faces of , respectively. One-dimensional linear interpolation of on the segment and on the segment is used for the cell-face fraction ratios computing:
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 are coordinates of the center of cell, are coordinates of the point around which the airfoil rotates, and and are the drag and lift forces acting on the solid part of the cut-cell at the -th time step:Here , , , , and 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 , , , 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, . 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 in the numerical simulation ( is a time discretization step). It is necessary that the filter cut-off frequency is equal to Hz and the suppression level on the cut-off frequency is equal to dB. So, the following condition imposed on the required filter frequency response function (Figure 8) must be satisfied for the required filter transfer function :
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 . It is known thatfor 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 is
4. Numerical Experiments
We considered autorotation of Savonius rotor with two blades (Figure 5) at . 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 . Reynolds number is based on rotor diameter. The following parameters were used in the simulation:
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: 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 . At the chosen values of the parameters, tip speed ratio is equal totorque coefficient value is obtained by averaging of nonstationary dependency over time, whereSimilarly, the power coefficient is obtained by averaging of nonstationary dependency over time, where
A number of computations have been performed on nonuniform grid . The uniform mesh block with space discretization step was used in the vicinity of the rotor. Time discretization step was equal to . 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 and is shown in Figure 10.
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 at .(v)Computational results are in good qualitative agreement with the experimental data.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research is partially supported by Russian Ministry of Education and Science (Project no. 9.2422.2017/PP) and Russian Federation President Grant for Young Russian Ph.D. Scientists (Project no. MK-7431.2016.8).