CFD simulation of tidal current farm by using AL model ＊

Cheng Liu, Changhong Hu

Research Institute for Applied Mechanics, Kyushu University, Fukuoka, Japan

Abstract: In this study, an efficient numerical method for predicting the wake interference of multiple turbines is presented. The actuator line (AL) model instead of the geometry-resolved method is adopted to represent the rotor. The large-eddy simulation (LES) is performed to predict wakes of multiple turbines operated in turbulent flows. An efficient immersed boundary (IB) method with moving least square reconstruction (MLS) is developed to model the nacelle and support structure of the tidal turbine. A simple wall function based on the MLS-IB method and boundary-layer equations is employed to compute the instantaneous wall shear stress. Laminar flow simulations of unsteady flows past a cylinder illustrate the accuracy of the wall function IB method. Finally, the proposed method is extended to study turbulent flow past tandem tidal rotors, in which the wake profile behind rotors is analyzed. The results are found to be in reasonable agreement with published data.

Key words: Tidal current farm, immersed boundary method, wall function model, large eddy simulation (LES), actuator line model

Introduction

Tidal current power is one of the most potential resources for future electricity generation, corresponding experimental investigations are attracting increasing interest. Mycek et al.[1]performed experiments for two tandem horizontal-axis tidal turbines (HATTs)in the IFREMER circulation flume tank, in which two ambient turbulence intensity (TI) rates (3% and 15%)are measured. Stallard et al.[2]carried out experiments to study the wake induced by staggered and aligned configurations, the wall effect was also under investigation.

In contrast, numerical approaches offer researchers a more promising way to estimate the hydrodynamic performance of multiple tidal turbines with significantly lower cost. The Reynolds-averaged Navier-Stokes (RANS) turbulence models are often used for predicting the mean performance and near-wake structure of a HATT. However, RANS model cannot provide an insight into turbulent velocity fluctuations, thus the effect of turbulence on tidal turbine loading cannot be investigated. On the other hand, although the computational expense of large eddy simulation (LES) is higher than RANS,blind test[3]shows that the velocity recovery and the ambient turbulence of a rotating turbine can be well-predicted by LES model rather than RANS.Turbulence models are often implemented based on the curvilinear or unstructured meshes for better capability of resolving boundary layers. This brings about three disadvantages, first, the accuracy of numerical scheme may be reduced due to the absent of mesh orthogonality, second, the convergence speed may be slowed down because of the stretched mesh,finally, mesh generation is cumbersome, thus adjusting the grid to adapt to numerous case-studies of a tidal farm is too time-consuming.

Recently, an alternative approach (immersed boundary method[4]) by a combination of Cartesian grid layout and local field reconstruction has been developed. By applying the IB method, the complex geometries are well-resolved while the merits of Cartesian grids are preserved. The signed-distance function is often used to represent the solid boundary implicitly[5-6]. In the direct-forcing method, the flow variables are often reconstructed through polynomial interpolation. In this way, the Dirichlet or Neumann boundary conditions to the solid surface can be imposed properly.

Although the computational cost for LES is much lower than that of directly numerical simulation (DNS)for homogeneous turbulence simulations, for boundary layers, however, requires the near-wall mesh resolution comparable to that for DNS, thus hindering the LES application to moderate Reynolds numbers.One of the solution for above problems is utilizing a wall function to impose appropriate boundary conditions for the LES[7-9]. The formulation of wall functions is often deduced from turbulent boundarylayer equations and its simplified formulations.

One objective of this research is to study the applicability of wall models for the wake interference problems of multiple turbines. The AL model with corrections to volume force calculation is introduced to represent the rotors. The MLS-IB method[10-11]considering the wall functions is proposed to model the hub and tower effect. The combination of AL model and IB method is highly efficient for configuration studies of multiple turbines. Both RANS (two equations k- ω model) and LES (Smagorinsky with wall functions) are performed for comparison.Numerical tests show that the predicted velocity deficit and turbulence intensity (TI) confirm the experimental results. In most of the cases, LES behaves better than that of RANS.

1. Mathmatical model

1.1 Governing equations

The filtered Navier-Stokes equations for moderately high Reynolds number flows are given by

In this study, governing equations (Eq. (1)) are discretized by collocated finite volume method. An implicit Euler scheme is used for time-stepping. The QUICK scheme is adopted for the flux reconstruction of advection phase, other terms use the central difference method. The pressure-implicit with splitting of operators (PISO) algorithm is used for decoupling the nonlinear pressure-velocity system into an implicit predictor and multiple explicit corrector steps. The Krylov iterative method with multi-grid preconditioner is applied for solving the linear system of pressure.

1.2 Immersed boundary method

To impose the boundary conditions sharply around the surface, a MLS-IB method[10-11]is applied to reconstruct the velocity at the forcing points adjacent to the rigid boundary. The classification of fluid points, solid points and forcing points are shown in Fig. 1(a). Noted that the sighed distance field is initialized around the solid surface regions before the start of simulation, thus the normal, distance, in/out of a solid surface can be determined quickly.

Fig. 1(a) (Color online) Classification of the points in the computational domain

Fig. 1(b) (Color online) Definition of the MLS supporting domain

For the present direct forcing method, the scalar quantity ϕ ( x) at the forcing point is determined by the moving least square approach with a supporting domain, see Fig. 1(b)

where pT( x) is the orthogonal basis function vector.The accuracy of the interpolation depends on order of the basis function. C ( x) is the coefficient vector,which can be determined by minimizing the following weighted residual

where xjstands for the fluid point and intersecting point fall into the supporting domain. n is the number of stencils used for reconstruction. W ( x - xj) is the weighed function that depends on the distance between x, xj. By solving Eq. (2), the coefficients vector C ( x) can be determined. Finally, the scalar quantity ϕ at the forcing point x can be calculated through Eq. (1). Details are described in Refs. [10-11].

As pointed out by Liu and Hu[11], it is better to transform the coordinate origin ( x) of stencils to the forcing point to maintain the numerical stability.Besides, as shown in Fig. 1(b), if other forcing points fall into the support domain, the corresponding intersection point will be used, instead of the forcing point. For the present study, a 3-D bilinear subspace assumption is adopted, p =[ x2, y2, z2, xyz , xy , yz , xz,x , y ,z,1] . Therefore, a 11×11 linear system is derived,which can be solved by Gaussian elimination method.

1.3 Wall function

One of the drawbacks that has limited the use of the IB method is due to the low-order velocity reconstructions. This makes the traditional IB method cannot be straightforwardly extended to high-Reynolds number problems. Actually, for a body-fitted grid, it is easy to cluster grid nodes in the wall normal direction to fulfil the near-wall resolution requirements. However, in terms of the Cartesian grid layout, the refinement should be performed in all three directions.Therefore, the computational load of the IB method based LES becomes extremely heavy (the number of points per cubic L:while, even by implementing the adaptive mesh refinement strategy.

An alternative approach to avoid the shortcoming is to build an approximated wall boundary conditions with an appropriate LES wall model. In the present work, a two-layer wall modelling method[9]is applied.As indicated in Fig. 2, Eq. (1) is only solved on the fluid points; for the forcing points adjacent to the wall,a boundary layer equation for the tangential velocity components is solved

where η indicates the wall normal direction, ε, ς are other two tangential directions. In this study, a simplified equilibrium stress balance model, recommended by Tessicini et al.[9]is adopted. The wall normal velocity uηis set to zero on the wall. The eddy viscositytν is computed from a mixing-length eddy viscosity model with near-wall damping

where+η represents the distance to the wall in wall units. Here k =0.4, a =19 is set for above equation.To solve boundary layer equation (Eq. (4)) numerically, we use the LES velocities at the border of the wall-layer and apply the no-slip conditions on the body surface+(η =0). It is noted that the computation oftν in Eq. (5) requires+η , which is determined by uτ. While uτis computed from Eq. (4)with undeterminedtν. As a result, Eqs. (4)-(5) are suggested to be solved simultaneously by an iterative procedure[9].

Fig. 2 (Color online) Sketch of the wall boundary reconstruction by introducing auxiliary points

1.4 Actuator line model

The AL approach is a combination of classical BEM theory and Navier-Stokes equation[12]. It is an efficient approach that can provide major wake flow characteristics of a rotating turbine. In AL approach[13],the rotor is defined in Cartesian coordinate frame with three actuator lines to represent rotating blades. A series of points along each blade’s axis are defined and each point is the center of an actuator element.The position of these points rotate according to the angular speed of the turbine at each time step. The critical point in implementing AL model is to add the force terms at the right side of the Navier-Stokes equation (Eq. (1)) to represent the rotor effect.

Fig. 3 (Color online) Comparison of mean streamwise velocities

Fig. 4 (Color online) Comparison of mean velocity fluctuations

2. Numerical result

2.1 Flow past cylinder ( Re=3 900)

For traditional gravity based tidal turbines, the rotors are supported by nacelle and tower. The diameters of nacelle and supporting structures are designed to be strong enough to resist the tidal flow.Therefore, the effect of nacelle and tower cannot be neglected especially in the near wake region. To validate the IB method and wall models of present solver, the simulation of flow past a circular cylinder at Reynolds number Re =3 900 is performed. Both the RANS model and LES are applied for the unsteady cylinder wake predictions. Solutions are compared with the reported experimental[14]and numerical results[15].

The dimensions of computational domain are 10D, 30D and 50D (D is the diameter of the cylinder) in three directions. A simple refining criterion is applied by which the finer resolution can be derived by partitioning a standard cell into 8 sub-cells. The recirculation lengthand the negative velocitypredicted by present IB method based LES agree well with the experimentaland body-fitted LES simulations

The mean streamwise velocity profiles at the cross-sectional plane are shown in Fig. 3, numerical results of seven locations in the wake are recorded.For the near wake region ( x/ D= 0.58, 1.06), the differences between the experimental and numerical are negligible. The over predicted velocity deficit can be observed for far field prediction. The streamwise and crossflow velocity fluctuations at three locations in the wake of a circular cylinder are presented in Fig.4. The LES case gives the closest predictions compared with the experimental data.

2.2 IFREMER rotor test

In order to validate the newly developed AL model, the experimental data of the HATT model test performed by Mycek et al.[1]is used for comparison.The experiment was carried out in the IFREMER flume tank, France, in which the 1/30 scaled prototype tidal turbine with different inlet TIs are studied. In the experiment, the upstream conditions of U∞=0.8 m/s, I∞= 15% are considered. The optimized tip speed ratio (TSR) is 3.67 for all of the IFREMER rotor test. Other parameters, such as the dimensions of domain, initial conditions and boundary conditions are found inRef. [1].

Fig. 5(a) (Color online) Fully developed turbulent flows in a precursor simulation

Fig. 5(b) (Color online) Velocity field of IFREMER’s double rotors interference

For the present LES, to construct inlet conditions that are statistically similar to the real oceanic turbulence, a fully developed turbulence flow field should be computed in advance, as seen in Fig. 5(a). Then a sequence of 2-D sampling planes (y-z plane for this case) with the full information of velocity are saved and used as the inlet boundary conditions in the nextstage LES simulation. Although the whole procedure is time-consuming and cumbersome, the ambient turbulence with the similar energy spectra and Reynold’s-stress profiles is expected as in the model experiment[19]. Instantaneous turbulent velocity contour of the final LES simulation is given in Fig. 5(b).

Since the rotor is placed in the turbulent flow generated by the upstream rotor, the actual ambient TI around the downstream turbine is larger than single turbine case. Figure 6 show the results of downstream turbine under high TI( I∞= 15%) cases in which the wake profiles of 7 transects (2 D-10 D behind the rotor plane) are plotted. It is observed the TI profiles predicted by the k- ω model are not consistent with experiment at x/ D≤5 (Figs. 6(a)-6(c)). The LES model seems to resolve the TI profiles better around the rotor hub (x/ D = 2) but still gives over-predict velocity recovery when x/ D=2 . For x/ D =10, the differences of the wake profiles between two simulations are not obvious (Fig. 6(d). Generally, the LES model shows better performance in predicting TI and velocity recovery than the RANS model.

3. Conclusion

In this paper, we present our numerical developments on prediction of the wake induced by the tidal turbine rotors. This study is inspired by our previous study on the IB method[5-6,11], where the MLS interpolation is introduced to approximate the wall boundary conditions. To alleviate the computational cost, the tidal rotors are represented by AL model rather than resolving the blade geometry directly. The IB based wall function is validated by numerical tests of flow past fixed cylinder. Besides, the AL model is verified by wake interference of two tidal rotors in turbulent upstream. The present numerical predictions are found to be in good agreement with published experimental and other numerical results. In our future work we will focus on extending the method to simulate large tidal turbine farms interacting with ocean boundary turbulence[20], by taking into consideration of the terrain and flow conditions in a complete tidal period.

Fig. 6 (Color online) TI profiles and velocity deficits along the rotor center for the IFREMER’s double rotors case TI=15% and 4D spacing