Fully consistent CFD methods for incompressible flow computations

Nowadays collocated grid based CFD methods are one of the most efficient tools for computations of the flows past wind turbines. To ensure the robustness of the methods they require special attention to the well-known problem of pressure-velocity coupling. Many commercial codes to ensure the pressure-velocity coupling on collocated grids use the so-called momentum interpolation method of Rhie and Chow [1]. As known, the method and some of its widely spread modifications result in solutions, which are dependent of time step at convergence. In this paper the magnitude of the dependence is shown to contribute about 0.5% into the total error in a typical turbulent flow computation. Nevertheless if coarse grids are used, the standard interpolation methods result in much higher non-consistent behavior. To overcome the problem, a recently developed interpolation method, which is independent of time step, is used. It is shown that in comparison to other time step independent method, the method may enhance the convergence rate of the SIMPLEC algorithm up to 25 %. The method is verified using turbulent flow computations around a NACA 64618 airfoil and the roll-up of a shear layer, which may appear in wind turbine wake.


Introduction
During the last decade a large effort has gone into developing CFD tools for prediction of wind turbine aerodynamics. In 1999, the flow over the NREL Phase II was computed using overset grid method by Duque [2]. Later, accurate predictions of the NREL Phase IV rotor were performed by Sørensen [3] and Johansen [4] and the stalled rotor prediction was also presented by Duque in [5]. To decrease the computational costs Xu and Sankar proposed in [6] an approach where small zones surrounding each blade were solved, whereas the rest of the domain was treated using a significantly less expensive full potential solver. Rather than modeling the entire rotor Pape in [7] modeled a single blade omitting the tower and nacelle. To predict the rotor-tower interactions, computations of the fully resolved rotors were performed by Zahle using incompressible overset grid method in [8]. Compressible sliding grid method was used by Gomez-Iradi in [9] for computations of the NREL Phase IV rotor. CFD modeling of laminar-turbulent transition for the wind turbine rotors was presented by Sørensen in [10]. In addition to the computations on structured grids computations using unstructured grid solvers were presented by Sezer-Udol in [11] and Potsdam in [12]. Wake simulations behind the wind turbines using CFD methods have been also an intensive field of study. Recent comparison of CFD wake simulations with the MEXICO experimental data can be found e.g. in Bechmann [13]. To decrease the computational costs, actuator disk and actuator line models have been applied for the wake modeling as can be seen in [14][15][16].
Nowadays, most of the incompressible CFD tools for wind turbine computations are based on the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) to link the velocity and pressure. To achieve a high robustness of the SIMPLE algorithm there is a common practice to employ it on collocated grids. One of the complexities of the collocated grid-based algorithms is the well known problem of the pressure-velocity decoupling. To overcome the problem, vast majority of the methods used today in commercial codes are based on the so-called momentum interpolation methods, initially proposed by Rhie and Chow [1].
Since the last two decades various modifications of the momentum interpolations have been presented in order to ensure an accurate solution on collocated grids. Originally the Rhie-Chow interpolation was developed for steady flow computations and is known to possess a dependence of velocity underrelaxation parameter at convergence. The problem with the dependence was solved in [17,18], nevertheless it was later shown in [19,20] that if the Rhie-Chow interpolation is used for unsteady flow computations, pressure wiggles appear for small time steps. An interpolation method for unsteady flow computations free from the pressure wiggles was later proposed independently by Choi [21] and Shen et al. [19]. Note that the method of Shen et al possesses the same properties as the method of Choi, but contrary to the Choi's method, it is based on second order scheme in time. As shown in [22] both methods of Shen et al and Choi possess a weak dependence of time step and relaxation parameter at convergence. To overcome the difficulty, several methods, which are independent of time step and relaxation parameter, were proposed in [20,[23][24][25]. Nowadays in spite of the existence of the time step independent methods, the standard methods of Choi and Shen et al are widely used as can be seen in example in [26][27][28][29].
On collocated grids both the solution accuracy and the convergence rate of the SIMPLE-like algorithms strictly depend on the choice of the interpolation method. Originally the SIMPLElike algorithms, such as SIMPLE [30], SIMPLEC [31] and PISO [32], were developed for the staggered grids, where mass flux interpolation is not necessary. In most of the computational codes used in engineering application the choice of the SIMPLE-like algorithm and the choice of the mass flux interpolation method have been done independently [33,34]. In the current paper, it is shown that using such the approach convergence rate is not always optimal and in order to achieve a high efficiency of the SIMPLE-like algorithms, the mass flux interpolation should be used consistently with the algorithms.
By taking an example of the SIMPLEC algorithm we will demonstrate that the usage of interpolation methods, which are fully compatible with the SIMPLEC algorithm, results in a convergence rate up to 25 % higher than the rate of the standard SIMPLEC algorithm. Theoretical justification of this fact is given in [35], whereas here by using a typical turbulent flow computations we will show that the proper choice of interpolation method may also increase the accuracy of the SIMPLEC algorithm.
Standard interpolation methods on collocated grids, such as the methods of Choi [21] and Shen et al [19], are known to result in solution, dependent of time step at convergence. In this work the magnitude of the dependence is estimated and is shown to be negligible for a typical turbulent flow computation. Nevertheless, when coarse grids are used the inconsistency of the solution becomes non-negligible. To overcome the problem the recently developed Modified Momentum Interpolation method (MMI) of Kolmogorov et al [35], will be used. The efficiency of the standard and the MMI interpolation methods will be tested in the roll-up case of a shear layer vortex, which appears in wind turbine wake. For the standard interpolation methods of Shen [19,36], the magnitude of the solution dependence of time step will be measured in application to the turbulent flow field around a NACA 64618 airfoil. For the test case the MMI method is proved to result in the solutions, which are independent of time step at convergence.

Discretization of the Navier-Stokes equations
The discretization of the incompressible Navier-Stokes equations for unsteady flow computations on collocated grids is presented below. First, the momentum equations are discretized. Second, the discrete continuity equation is presented with one of the recently developed interpolation methods of Kolmogorov et al in [35].

Momentum equations
Using the second order backward difference scheme in time and grouping two momentum equations in the x-and the y-directions together, the system of the underrelaxed momentum equations on collocated grids can be expressed for a control volume p in the following form: where vector υ p denotes the vector of the velocity flow field (u p , v p ) T , whereas vector S p indicates the explicitly treated source terms (S x p , S y p ) T . Terms A p , A nb are the diagonal and non-diagonal matrix terms accounting for the discrete convective and diffusion terms, the term A V p is the coefficient of the time derivative and equal to ρdVp τ , where τ is the time step and ρdV p is the control volume mass. The term A p accounts for the diagonal term of the matrix of the momentum equations, which are underrelaxed using the spatial term A p similarly to [19,24,36] as below: where α is the velocity underrelaxation parameter.
The superscripts m and n are the subiteration-and time step-counters respectively, such that the solution at time step n + 1 is obtained at convergence. In order to compute the flow field at the subiteration step m + 1, the coefficients A p , A nb and A p are taken from the former subiteration step m. In the notations of the coefficients the superscript counter m is dropped in the following section.

Mass flux interpolation methods
To ensure flow field continuity the momentum equations in Eq. (1) have to be solved together with the continuity equation: which essentially represents the fact that the sum of mass fluxes f m+1 k through the control volume faces k equals to zero.
On collocated grids the cell face fluxes are not available. To identify the fluxes one of the widely spread methods is the momentum interpolation method originally proposed by Rhie and Chow [1]. Nowadays, for unsteady flow computations there is a common practice [26][27][28][29] to employ the modified Rhie-Chow interpolations, proposed by Choi or Shen et al in [19,21]. It is known that, the methods are weakly dependent of time step at convergence, but the magnitude of the dependence in real life applications is not known yet.
Nowadays, there exist several interpolation methods on collocated grids, resulting in solutions, independent of time step at convergence [20,[22][23][24][25]. Most of the methods are aimed for  [35]. From one hand the MMI method results in the solution, which is independent of time step at convergence. From the other hand, as will be later seen in the paper, it may increase the convergence rate of the SIMPLEC algorithm up to 25% in comparison to other interpolation method existing in literature.
According to the MMI method the mass flux at some cell face k is defined as below: where [ ] k denotes linear interpolation from cell centers to cell face k, the term χ k and the vector h p m+1 are defined respectively as: where the new parameters γ and β are the constant parameters described below. In Eq.  The MMI method above results in solution independent of time step and relaxation parameter, as will be seen in the result section. The method has two forms, namely the first form with γ = 0, β = 0 and the second form with γ = 1, β > 0. The two forms of the method possess different properties as described below: 1) The MMI method with γ = 0 and β = 0 is similar to an existing method of Pascau in [25], but contrary to Pascau's method, it is based on momentum equations, which are discretized using second order backward difference scheme in time. As shown in [35], the MMI method in this form is fully compatible with SIMPLE algorithm.
2) The MMI method with γ = 1 and β > 0 is another interpolation method, independent of time step and relaxation parameter at convergence. According to [35], for the robust performance of the method, the parameter β has to be equal 0.04. Contrary to other time step independent methods, the MMI method in this form is fully compatible with SIMPLEC algorithm. As will be later seen in the result section, if the MMI method is applied with SIMPLEC algorithm, it becomes advantageous in the convergence speed over other time step independent interpolation method.

Results
Two test cases are computed, namely an idealized case of a shear layer [37], which appears in wind turbine wake, and the turbulent flow around a NACA 64618 arifoil. The steady and unsteady state solutions are achieved when the residuals, computed in L 1 norm, are reduced by a factor of 10 8 and 10 6 , respectively. For the tests the velocity underrelaxation parameter α = 0.8 is used.

Roll-up of a shear layer vortex
An idealized case of a shear layer, which appears in wind turbine wake, is computed at Re = 100 using SIMPLEC algorithm and the MMI interpolation method in Eq. (4). To achieve the high efficiency of the SIMPLEC algorithm, the MMI method is used with γ = 1 and β = 0.04 (see Section 2.2). To verify the efficiency of the MMI method, the method is compared with the interpolation method of Pascau [25].
The flow is initialized in the domain of a square of [0 π] with periodic boundary conditions as following: , y > π δ = π/15, = 0.05   An example of the flow field at a dimensionless time t = 8 is shown in Fig. 1. The solutions of the two methods are compared at the dimensionless time t = 4 and the errors are measured in comparison to the reference solution of the MMI method obtained on the fine grid of 512 2 cells and with a small time step of 0.005. The computations are performed on a sequence of the successively coarsened grids. For each coarse grid, the error is computed by interpolating the reference solution on the coarse grid and subtracting the solution from that grid. As seen from Fig. 2, the second order of the spatial accuracy is obtained and the error tolerances of both methods are nearly identical. The corresponding work load of the methods, measured in CPU seconds, is plotted in Fig. 3. The efficiencies of the methods can be compared if the work loads are compared for same accuracy levels. Therefore, the work load for each of the method is found first as a function of tolerance error. Then the ratio of the work loads of the MMI method, W ork M , and Pascau's method, W ork P , is computed and plotted in Fig. 4.
As seen in Fig. 4, the work load of the MMI method is lower then the work load of Pascau's method for nearly all accuracy levels. For the highest accuracy level, the efficiency of the MMI method is 25% higher in comparison to Pascau's method. As seen in Fig. 4 there is a narrow accuracy region where the MMI method is less efficient. This is explained by the fact that the optimal β parameter of the MMI method depends on grid resolution. In spite of the fact that the optimal β is not known in advance, the results of the test case and the deliberate analysis of the optimal β presented in [35] show that for general use β = 0.04 can be employed to ensure higher efficiency of the MMI method.

Turbulent flow around a NACA 64618 airfoil
The turbulent flow around NACA 64618 airfoil at zero angle of attack is computed at Re = 1.6 · 10 6 using SIMPLEC algorithm. The grid of O-type is used with boundaries placed in a distance of 20 chords from the airfoil. Shen et al [19] MMI, γ = 0, β = 0 (a) Shen's method [19] and the MMI method.   [19,36] in the turbulent flow computations around a NACA 64618 airfoil on grids with 64 x 32 and 128 x 64 cells. Two standard interpolation methods of Shen et al in [19] and [36] are used as the representatives of interpolation methods, resulting to the solutions, dependent of time step at convergence. To compare with the standard methods, the MMI method in Eq. (4) is used.
The k − ω SST turbulence model is used on two relatively coarse grids with 64 x 32 and 128 x 64 cells. For the two grids the maximum y + at points one cell away from the airfoil equals to 2.6 and 0.4, respectively. To measure the time step dependence of the standard methods, the lift coefficient is compared against an experimental value in [38], which is equal to 0.44. As seen from Figs. 5(a) and 5(b), contrary to the MMI method, the solutions of the standard methods of Shen, are dependent of time step at convergence. For the standard methods, the change of the error due to the time step dependence is about 1-2% on the grid with y + = 0.4, whereas for the grid with y + = 2.6 the error variations may achieve up to 5.5%, as seen from Fig. 5(b).
It should be noted, that a typical grid set-up for the turbulent flow computations is based on a grid with 128 cells in the normal direction. Results of computations on a grid with 256 x 128 cells using the standard and the MMI methods are shown in Fig. 6. It is seen from the figure that on such fine grid, the time step dependence of the error contributes about 0.5% into the total error.
The solution, obtained using the MMI method with γ = 0, β = 0 on the grid with 256 x 128 cells is seen to be about 4% different from the experimental value. The solution is less accurate than the solution obtained on the coarser grid using the MMI method with γ = 1, β = 0.04. This is explained by the fact that for the steady state problems, the MMI method with γ = 1 and β = 0.04 may result in superconvergence, as was reported in [35].

Conclusions
It is concluded, that for a typical turbulent flow computation the standard momentum interpolation methods still can be used. For the methods the solution dependence of the time step is negligible. However, in optimization tasks, where coarse grids are used, the standard methods may result in non consistent solution behavior. The inconsistency may become crucial in 3D computations, where employing the fine grids becomes computationally demanding. For the turbulent flow computations, the MMI method of Kolmogorov et al [35] was proved to result in solutions, independent of time step at convergence. The results of the unsteady flow computations have also shown, that when SIMPLEC algorithm is used, the MMI method, results in up to 25% higher convergence rate, than an existing method. In general, the MMI method can be used on both coarse and fine grids to ensure both high convergence rate and high accuracy of the SIMPLE-like algorithms.