Table of Contents Author Guidelines Submit a Manuscript
International Journal of Rotating Machinery
Volume 2017, Article ID 6418108, 7 pages
Research Article

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

Applied Mathematics Department, Bauman Moscow State Technical University (BMSTU), 5 2nd Baumanskaya, Moscow 105005, Russia

Correspondence should be addressed to Ilia K. Marchevsky; ur.liam@yksvehcramaili

Received 7 April 2017; Revised 21 July 2017; Accepted 14 August 2017; Published 11 October 2017

Academic Editor: Lin Zhu

Copyright © 2017 Ilia K. Marchevsky and Valeria V. Puzikova. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


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 [16] 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).

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 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].

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 (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 .

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 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.

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 . 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.

Figure 7: The time dependency for the torque (on the coarse grid with time step ) and filtered time dependency 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: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 :

Figure 8: The required filter frequency response 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.

Figure 9: Comparison of computed power coefficients for a two-bucket Savonius rotor at with experimental data (Sheldahl et al. [17] at and Shankar [18] at ).
Figure 10: Time evolution of the power coefficient at and . The red line corresponds to the average value .

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.


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).


  1. P.-A. Krogstad and J. A. Lund, “An experimental and numerical study of the performance of a model turbine,” Wind Energy, vol. 15, no. 3, pp. 443–457, 2012. View at Publisher · View at Google Scholar · View at Scopus
  2. N. N. Sorensen, J. A. Michelsen, and S. Schreck, “Navier-stokes predictions of the NREL phase VI rotor in the NASA Ames 80-by-120 wind tunnel,” in Proceedings of the 2002 ASME Wind Energy Symposium, pp. 94–105, Reno, Nev, USA, January 2002. View at Scopus
  3. K. Rogowski and R. Maroński, “CFD computation of the savonius rotor,” Journal of Theoretical and Applied Mechanics, vol. 53, no. 1, pp. 37–45, 2015. View at Publisher · View at Google Scholar · View at Scopus
  4. R. Gupta and K. K. Sharma, “Flow physics of a combined Darrieus-Savonius rotor using computational fluid dynamics (CFD),” Int. Research J. Eng. Science, Technology and Innovation, vol. 1, pp. 1–13, 2012. View at Google Scholar
  5. L. Oggiano, “CFD Simulations on the NTNU Wind Turbine Rotor and Comparison with Experiments,” in Proceedings of the Renewable Energy Research Conference, RERC 2014, pp. 111–116, June 2014. View at Publisher · View at Google Scholar · View at Scopus
  6. M. H. Mohamed, A. M. Ali, and A. A. Hafiz, “CFD analysis for H-rotor Darrieus turbine as a low speed wind energy converter,” Engineering Science and Technology, an International Journal, vol. 18, no. 1, pp. 1–13, 2015. View at Publisher · View at Google Scholar · View at Scopus
  7. R. Mittal and G. Iaccarino, “Immersed boundary methods,” in Annual review of fluid mechanics. Vol. 37, vol. 37 of Annu. Rev. Fluid Mech., pp. 239–261, Annual Reviews, Palo Alto, CA, 2005. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  8. Y. Cheny and O. Botella, “The LS-STAG method: A new immersed boundary/level-set method for the computation of incompressible viscous flows in complex moving geometries with good conservation properties,” Journal of Computational Physics, vol. 229, no. 4, pp. 1043–1076, 2010. View at Publisher · View at Google Scholar · View at MathSciNet
  9. S. Osher and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, vol. 153, Springer, New York, NY, USA, 2003. View at Publisher · View at Google Scholar · View at MathSciNet
  10. H. A. van der Vorst, “BI-CGSTAB: a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems,” SIAM Journal on Scientific and Statistical Computing, vol. 13, no. 2, pp. 631–644, 1992. View at Publisher · View at Google Scholar · View at MathSciNet
  11. P. Wesseling, An Introduction to Multigrid Methods, John Wiley & Sons, Chichester, UK, 1992. View at MathSciNet
  12. J. van Kan, C. Vuik, and P. Wesseling, “Fast pressure calculation for 2D and 3D Time dependent incompressible flow,” Numerical Linear Algebra with Applications, vol. 7, no. 6, pp. 429–447, 2000. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  13. I. Marchevsky and V. Puzikova, “OpenFOAM iterative methods efficiency analysis for linear systems solving,” Proceedings of the Institute for System Programming of RAS, vol. 24, pp. 71–86, 2013. View at Publisher · View at Google Scholar
  14. V. V. Puzikova and I. K. Marchevsky, “Application of the LS-STAG immersed boundary method for numerical simulation in coupled aeroelastic problems,” in Proceedings of the Joint 11th World Congress on Computational Mechanics, WCCM 2014, the 5th European Conference on Computational Mechanics, ECCM 2014 and the 6th European Conference on Computational Fluid Dynamics, ECFD 2014, pp. 1995–2006, esp, July 2014. View at Scopus
  15. V. V. Puzikova, “Construction of Level-Set Function for an Airfoil of Arbitrary Topology when Modelling a Flow past It Using the LS-STAG Method,” in Herald of the Bauman Moscow State Technical University. Natural Sciences, pp. 163–173, 2012 (Russian). View at Google Scholar
  16. R. G. Lyons, Understanding Digital Signal Processing, Pearson Education, 2010.
  17. R. E. Sheldahl, B. F. Blackwell, and L. V. Feltz, “Wind tunnel performance data for two- and three-bucket savonius rotors,” J Energy, vol. 2, no. 3, pp. 160–164, 1978. View at Publisher · View at Google Scholar · View at Scopus
  18. P. N. Shankar, National Aeronautical Laboratory, AE-TM-3-76, National Aeronautical Laboratory, Bangalore, India, 1976.