1. Introduction
Wind energy is a substantial paradigm of human technology in which atmospheric air flows are leveraged for environmentally neutral energy generation. In the effort to reduce wind energy costs, the wind energy sector is directed towards increasing the size and flexibility of modern wind turbine (WT) rotors. To that end, rotor blades are becoming larger and more slender, leading to high complexity regarding the analysis of their loading and the description of the flow field in their wake. In connection to the above need for highly flexible designs, rotor analysis is today directed towards high fidelity numerical tools, applying elaborated, physically motivated models for the accurate prediction of the response of the rotor blades. This requires the use of highly accurate computational models both for the elasto–dynamic and the aerodynamic analysis of the blades, which, combined together, comprise the field of aeroelasticity.
In terms of aerodynamic analysis, computational fluid dynamics (CFD) methodologies are considered to be the most common option among several numerical approaches in terms of accuracy. The advantage of employing CFD in rotor analysis is related to its capability to account in maximum detail for viscous and turbulence effects. These are dominant in wind farm simulations, where rotor wake evolution has an impact on the performance of downwind turbines. Furthermore, viscous effects are important when studying the rotor wake interaction with the boundary layer developed on surrounding bodies, such as the WT tower or even the ground. Such interactions affect the rotor performance and the blades’ loading.
Despite being accurate and reliable, the increased computational cost of conventional Eulerian CFD methodologies renders them prohibitive for aeroelastic design simulations. As a remedy, research is today directed towards the development of novel CFD methodologies that share the same level of accuracy as the conventional methodologies under reduced computational requirements. In this direction, domain decomposition seems to be an attractive technique provided that the best performing formulation in terms of cost and accuracy is employed in each sub-domain. The proposed hybrid CFD methodology combines a standard finite volume Eulerian CFD approach close to the solid boundaries with a Lagrangian CFD approach in vorticity formulation for the rest of the domain [
1,
2,
3].
The flow–field may be discretised in two different ways, thus classifying CFD methodologies into Eulerian and Lagrangian methodologies. On the one hand, Eulerian approaches discretise the entire computational domain in “stationary” nodes through which the fluid moves and on whom the flow properties are recorded. On the other hand, in Lagrangian methods, the fluid is discretised in “numerical particles”. In this way, the flow-field is described in a material approach by tracking the motion and following the flow properties of the fluid particles.
Eulerian methodologies [
4,
5,
6] have gained much popularity and reliability, thanks to the accurate consideration of solid-wall boundaries, as no-penetration and no-slip conditions are both satisfied in maximum detail. However, this is not the case regarding the far-field boundaries. Truncation of the flow at a finite distance, where the boundary conditions typically approximate those at infinity, is a numerical approximation that leads to errors and is considered as significant drawback in external aerodynamics applications. Furthermore, domain truncation is usually combined with gradual grid coarsening, which increases numerical diffusion and adds errors as well. This may be catastrophic in wind farm applications, where a detailed description of the wake of the leading WTs is crucial for the correct estimation of the total power production. Local grid refinement seems like an obvious remedy, but it comes with a substantial increase in computational cost. Moreover, the consideration of the interaction between independently moving bodies is usually implemented through the application of special methodologies (e.g., overset, chimera or sliding grids), which penalize computational cost as well.
The alternative would be to solve the Navier–Stokes equations in the Lagrangian formulation, as particle methods do. They are grid–free, self–adaptive and have (in theory) zero numerical diffusion. The intuitive option would be to define particles that carry mass, momentum and energy, as the smooth particle hydrodynamic (SPH) methods do [
7,
8,
9]. Alternatively, vorticity [
10,
11] and dilatation [
12,
13] may be used as the primary flow variables leading to vortex methodologies (VM). Due to their robustness in pressure variations, VM are quite popular in external aerodynamics applications [
14,
15,
16]. Consequently, they are considered suitable for WT applications [
17,
18,
19,
20,
21] as well. A major challenge for particle methods is the treatment of solid-wall boundaries [
22,
23]. Another drawback is the fact that the computational cost rises proportionally to
(
N is the number of particles), thus rendering long–lasting simulations computationally prohibitive. To overcome this problem, methodologies such as the particle-mesh (PM) [
24] may be adopted in order to enhance performance.
To sum up, Lagrangian methodologies are more effective in the far-field region, whereas Eulerian methodologies are considered to be ideal for the region close to the solid-wall boundaries. It is therefore reasonable to combine the two methodologies through a domain decomposition technique in order to enhance accuracy and reduce computational cost. In general, the sub-domains may either overlap or not. Strong viscous-inviscid interaction models [
25,
26] and RANS-vortex coupled models [
27,
28] are examples of completely overlapping hybrid methodologies, while the method presented in [
29] is an example of limited overlapping over a buffer area.
In the current project, an in-house hybrid solver, named HoPFlow, is presented and assessed in WT applications. A preliminary version of the solver and results concerning low subsonic turbulent flows around airfoils are presented in [
1]. Implementation details have been presented analytically in [
3], where compressibility effects have been also taken into account. Recently, the hybrid solver has been used in
, WT applications where an axial flow test case of the New MEXICO experimental campaign [
30,
31] was simulated. Preliminary results have been presented in [
32], focusing on aerodynamic load prediction and near-wake flow field computations. In this study, the numerical parameters of the presented solver that affect load predictions (e.g., blade surface mesh, time-step value, PM discretisation length) are analyzed in detail. The main motivation of the present work is to validate the flow-field characteristics close to solid boundaries using measured datasets and thus assess the coupling procedure applied in the hybrid solver. As indicated in the discussion above, it is noted that the value of the present method lies in applications in which there are multiple, mutually interacting regions of interest within the computational space [
33] (e.g., multiple-rotor configurations). The Lagrangian domain serves as the coupling domain between these regions of interest in the same way that overset grids function in standard Eulerian CFD simulations. Lagrangian methods are advantageous compared to overset approaches, as they minimize numerical diffusion and therefore preserve wake structures and flow disturbances.
The test case is the run no. 266, which is a 14.7 m/s axial flow case at 425 rpm of the New MEXICO experimental campaign [
30,
31] that corresponds to a tip speed ratio of
= 6.81 and a pitch angle of
nose down. The New MEXICO WT model rotor is a 4.5 m diameter 3–bladed rotor. It consists of three different airfoil shapes at the root (DU91-W2-250), mid-span (RISØ A1-21) and tip (NACA 64418) region of the blade according to
Table 1. The twist and chord distribution of the blade is shown in
Table 2. Turbulent transition is triggered with trip tapes of 5 mm width and 0.2 mm thickness placed at 10% chord of both pressure and suction sides of the blade. A thorough comparison of blade loads is performed against experimental measurements and computational results from the Eulerian counterpart (MaPFlow) of the hybrid solver. Blade loads are depicted as integrated rotor loads, span–wise distribution of aerodynamic loads and pressure distribution at specific span–wise positions. More results from numerical investigations concerning test cases from the MEXICO and the New MEXICO experimental campaigns can be found in [
34,
35,
36,
37,
38,
39,
40].
2. Methodology
The rationale behind the application of the hybrid CFD solver HoPFlow is to combine an Eulerian approach close to solid-wall boundaries with a Lagrangian one for the rest of the domain. In this way, both the solid-wall and the exact far-field boundary conditions are satisfied in maximum accuracy within the Eulerian and Lagrangian framework, respectively. The Lagrangian particles are distributed over the whole computational domain, overlapping with the Eulerian computational cells close to solid boundaries (see
Figure 1). The Eulerian part of the hybrid solver solves the compressible Navier–Stokes equations under a cell-centered finite-volume approach in a confined region
around solid-wall boundaries
. The Lagrangian part solves the compressible flow equations as well, in their material form, based on particle representation of the essential flow quantities, e.g., mass, pressure, dilatation and vorticity [
12]. Lagrangian particles are distributed over the entire computational space and, hence, inside the Eulerian domain as well. The coupling procedure consists of two parts: (i) the Lagrangian particles solution is interpolated to the ghost cells of the Eulerian domain to define the correct boundary conditions on its far-field
(see
Figure 2); (ii) the Eulerian solution is used in order to correct the solution of the Lagrangian particles that lie within the Eulerian domain (see
Figure 3). The corrected Lagrangian particles information is then used in order to update the whole Lagrangian field and thus ensure that it will be a smooth extension of the Eulerian one. More details concerning the coupling procedure are given in
Section 2.3.
2.1. The Eulerian Solver
MaPFlow [
41] is an in-house typical Eulerian CFD solver which solves the compressible unsteady Reynolds averaged Navier–Stokes (URANS) equations under a cell-centered finite volume spatial discretization scheme. MaPFlow can handle both structured and unstructured grids; it is parallelized under the MPI protocol, and the grid partitioning is performed using the METIS library [
42]. The convective fluxes are evaluated by solving the preconditioned local Riemann problem between the neighboring cells of each face, using the Roe’s approximate Riemann solver [
43] with the Venkatakrishnan limiter [
44]. The viscous fluxes are discretized using a central 2nd order scheme. For the reconstruction of variables at the interface, a piecewise linear interpolation scheme is used. The evaluation of the spatial gradients of the primitive variables is done using the Green–Gauss formula, with a centered scheme approximation. Multiple options are available for turbulence modeling, such as the one-equation model of Spalart–Allmaras [
45] or the two-equation model
of Menter [
46]. Regarding laminar to turbulent flow transition modeling, the
model of Menter [
47] is used. A delayed detached eddy simulation (DDES) approach is also implemented in MaPFlow, following the suggestions of [
48]. Unsteady simulations are performed through an implicit second-order backwards difference scheme [
49], along with a dual time-stepping technique [
50] in order to facilitate convergence. Finally, the implicit operator inversion is accomplished with the use of the Gauss–Seidel iterative method alongside the reverse Cuthill–Mckee reordering scheme.
2.2. The Lagrangian Solver
In a Lagrangian formulation (material coordinates), the flow-field is represented by following the evolution of a number of particles along their trajectories. In that sense, particles act as flow marker points that are assigned with volume
and carry mass
, dilatation
, vorticity
and pressure
, regarded as the volume integrals of the continuous flow quantities density
, dilatation
, vorticity
and pressure
p respectively. In material coordinates, the flow equations take the form:
where
denotes the material time derivative, and
indicates evaluation at the position of particle
p.
denotes the divergence of the viscous stress tensor, and
is the kinematic viscosity, which here is assumed constant.
The flow equations are supplemented with the Helmholtz’s decomposition theorem (
2), which states that every velocity field
can be expressed as the sum of a rot-free potential part
and a div-free vortical one
, alongside a constant velocity component representing the undisturbed velocity field at infinity
. The potential part is defined through a scalar potential
and is associated with the compressibility effects expressed by the dilatation of the flow
, whereas the vortical part is defined through a vector potential (stream-function)
which is associated with the free vorticity of the flow
. Consequently, the scalar and vector potential satisfy the Poisson Equation (
3).
By using Green’s theorem, the velocity field
can be expressed in integral form:
where
is the Green’s function for the Laplace operator,
and
. Computational cost is dominated by the convolution integral in (
4). For
N particles, the associated cost is proportional to
, which can easily explode as
N becomes large and the intended duration of the simulation is long. To make things even worse, when boundaries are present, the surface convolutions in (
4) must be also evaluated. In order to reduce computational cost, the PM technique is employed, and the Poisson Equation (
3) is solved for the scalar potential
and the stream function
. In such a manner, computational cost is minimized from
to
. Computational performance is also enhanced by using Cartesian grids in order to discretize the Lagrangian domain, thus enabling the employment of fast Poisson solvers [
51]. Particularly, in HoPFlow, the James–Lackner algorithm is used [
52].
The PM framework is also employed in order to evaluate the right-hand side (RHS) of (
1). The Lagrangian particles solution is interpolated to the PM nodes, and the desired differentiations are easily computed through finite difference schemes. Consequently, the RHS terms are first evaluated on the PM nodes, and then they are interpolated back to the particles positions. Afterwards, time marching is performed through a standard 4th order Runge–Kutta explicit scheme. In every sub-step of Runge–Kutta, intermediate convection steps are carried out, requiring intermediate evaluations of velocity, which are also conducted with the usage of the PM technique.
At the end of every time-step, remeshing is applied in order to recover full coverage of the computational domain and ensure a regular distribution of the numerical particles. Remeshing is a standard procedure in order to prevent excessive concentration or spreading of particles and, in this way, preserve the consistency and accuracy of the numerical solution.
For a given number of particles , the sub-steps taken in the Lagrangian time-step can be listed as follows:
- Step 1:
Project on the PM grid and obtain ;
- Step 2:
Solve and obtain ;
- Step 3:
Calculate on the PM grid the RHS terms of (
1), e.g.,
;
- Step 4:
Interpolate all grid-based data at the particle positions ;
- Step 5:
Update all particle properties (integrate (
1) in time);
- Step 6:
Re-mesh if needed.
2.3. The Hybrid Solver
The hybrid solver, HoPFlow, couples the two distinct, previously described, Eulerian and Lagrangian solvers. The Eulerian solution is used only in the regions close to the solid boundaries, whereas the Lagrangian one needs to be valid on the whole computational space, overlapping with the Eulerian one (see
Figure 1), in order to satisfy the true far-field boundary conditions.
As mentioned before, the Eulerian and Lagrangian solutions are coupled in two ways. First of all, the Lagrangian part shall undertake to provide the proper flow conditions on the outer boundaries of the Eulerian domain
. In order to do so, the Lagrangian solution is interpolated from the PM nodes (or, in a more generic approach, from the Lagrangian particle positions) to the ghost cells of the Eulerian grid (see
Figure 2). The fluxes at the Eulerian boundary
may now be evaluated from the Riemann invariants by taking into account the correct flow information on both sides of
. Using the correct boundary conditions, MaPFlow shall be capable of describing the flow-field close to the wall boundaries in the detail that is provided by the Eulerian framework.
The closure of the coupling is achieved by correcting the flow information on the PM nodes (more generally on the Lagrangian particles) that lie within the Eulerian domain
and then by updating the whole Lagrangian field. In order to do so, the Eulerian solution is transformed into particles that carry mass, dilatation, vorticity, pressure and volume
. The Eulerian particles need to be densely populated and regularly placed within the Eulerian computational cells, so that full coverage of the PM nodes is ensured. The flow quantities of the Eulerian particles are interpolated from the cell-centered values based on a purely geometric approach using iso-parametric finite element approximations. The presence of solid boundaries is taken into account as surface (singular) particles that carry dilatation
and vorticity
, but no pressure, volume and mass. These particles only affect the solution of the Poisson Equation (
3) as contribution in
and
in the surface convolution related to boundary terms (
4) and must not be convected during time marching.
The steps followed for the coupling between the Lagrangian and the Eulerian solver are depicted in the flow chart displayed in
Figure 4. The corresponding implementation details have been thoroughly analyzed in [
3].
4. Conclusions
In this paper, the hybrid Lagrangian–Eulerian solver HoPFlow has been presented, along with details about its distinct solvers and the coupling between them. The Eulerian solver, MaPFlow, solves the compressible Navier–Stokes equations under a cell-centered finite-volume discretization scheme, while the Lagrangian solver uses numerical particles that carry mass, pressure, dilatation and vorticity as flow markers in order to represent the flow-field by following their trajectories. The calculation of the velocity field is performed by utilizing the decomposition theorem introduced by Helmholtz. Computational performance is enhanced by utilizing the particle mesh (PM) methodology in order to solve the Poisson equations for the scalar potential and the stream function .
The hybrid solver is tested in 3-D unsteady simulations concerning the axial flow around the wind turbine model rotor used in the New MEXICO experimental campaign. Simulation results are presented as integrated rotor loads, span-wise distribution of aerodynamic loads and gauge pressure distribution at various span-wise positions along the rotor blades. Comparison is made against experimental measured data and computational results produced with the Eulerian solver.
Under a PM discretisation with 5 nodes per chord length (in this analysis the characteristic chord length is assumed to be the on at of the blade), the hybrid solver predicts the loading close to the blade tip with great accuracy. This is regarded to be very important, as the aerodynamic loads close to the blade tip tend to dominate the whole rotor performance and the dynamic excitation of the blades in aeroelastic simulations. Deviations from the experimentally measured values are observed close to the blade root, which become more pronounced in the middle region of the blade. However, very good agreement is shown with the corresponding Eulerian solver results. These discrepancies are attributed to the absence of the spinner geometry in both hybrid and Eulerian solver simulations. Furthermore, it needs to be highlighted that the Lagrangian formulation used for the wake description, alongside with employment of vorticity as the primary flow quantity of the particles, reduces numerical diffusion significantly compared to the Eulerian solver. For this reason, the wake is resolved much more efficiently in the hybrid solver simulations, thus resulting in greater wake-induced velocities. Consequently the hybrid solver ends up with a better estimation of the aerodynamic loads along the whole blade span that is closer to the experimental values.
Summarizing the above, the results presented in this work indicate that the accuracy of the boundary layer solution (near-body flow-field) of the hybrid solver is comparable with that produced by a standard Eulerian CFD code. This confirms that the coupling method that determines the boundary conditions for the confined Eulerian grid is adequate and consistent. Nevertheless, the cost of the hybrid approach is overwelming for a single rotor simulation. A remedy for moderating the computational cost is the application of multi-level/telescopic PM grids [
15], finer in the vicinity of the body and coarser in the far-domain. The benefit for paying this increased computational cost (using fine PM grids in the entire computational domain) lies in unsteady applications where interaction phenomena (e.g. rotor-rotor interactions) are prominent or in cases where the accurate characterization of the far-wake dynamics is required. It is also attractive for aeroelastic analyses in which a standard Eulerian methodology relies on the use of overset grids and the corresponding coupling between different sub-domains, which penalizes computational cost.