Computational modelling of a large dimension wind farm cluster using domain coupling

The accuracy of Computational Fluid Dynamics (CFD) models for Atmospheric Boundary Layer (ABL) flows relies largely on the placement of the domain boundaries and the quality of the imposed flow conditions, the inlet boundary in particular. Exploiting the parabolic nature of many ABL flows and of CFD modelled ABL flow in particular, a precursor simulation is used as source of flow data to improve the target domain's inlet flow description over the standard synthetic boundary conditions, one-directionally coupling the solutions to the two simulations. Using the approach, a case of flow over a two wind farm offshore cluster is modelled using two small coupled simulations, matching the results of a single simulation including the full cluster at a significant computational time saving, in the order of 70%. Further savings were shown to be possible by reducing the resolution of the precursor simulation, with negligible impact on the results at the target domain.


Introduction
CFD engineering tools for modelling of ABL flows are limited by the tool's ability to simulate the full extent of the flow detail and scale: the model's domain can only represent a finite portion of the terrain surrounding the area of interest (AOI) and of the atmosphere's height. This poses the problem of boundary placement in relative proximity of the AOI and estimating the flow at those boundaries. Depending on the complexity of the surroundings, the synthetic boundary conditions the modeller is forced to impose may bear little resemblance to real world flow conditions. Increasing the distance between boundaries and AOI minimizes the issue but leads to increased mesh node count and exponentially increased computation time requirements, whereas improved boundary conditions can boost model accuracy with no direct computational cost.
Typically micro-scale ABL flow modelling is done by assuming the wind is flowing steadily for an indeterminate amount of time, in a certain direction over a large (larger than that of depicted by the domain) span of terrain, providing a flow solution to very precise ideal set of atmospheric conditions. This is achieved by considering flow through a wall bounded "channel", with boundary conditions set for all surfaces forming the box-domain over the terrain boundary: typically inlet and outlet boundaries plus the side and top boundaries as zero-shear walls. While the flow is unrestricted by these boundary conditions in the inner reaches of the domain, they force a dominant flow direction and so remove some of degrees of liberty from the overall flow description.
Forced to comply to boundary conditions in a "channel" flow, ABL flows (like other high Reynolds number flows) tend to be parabolic in behaviour: flow is predominant in a main direction, in which transport is dominated by convection, diffusion and cross-stream pressure variations are negligible and no reverse flow occurs [1]. This way, in simple flow cases where these conditions are met it is acceptable to consider the stream-wise coordinate to have a oneway character and flow to be dependant on upstream conditions alone. Linearised Reynolds averaged Navier-Stokes (RaNS) equation solvers (arguably the most common form of CFD ABL flow modelling) take advantage of this very fact, allowing them to produce accurate results quickly [2].
The inflow description is thus a dominant factor in AOI flow solution, more so than any other boundary condition. The commonly used ABL inlet boundary conditions [3,4] are adequate when inflow is known to be free from obstructions (terrain, but also canopy and/or wind turbines) and heterogeneous across the whole inlet boundary, but in most cases it is easy to determine that it is not so. In those cases improving the inlet boundary description is key to increasing accuracy of the AOI results, and using other CFD models as sources for boundary conditions an accessible form of doing so. The model-chain approach [6] uses one-directional CFD model coupling to improve flow description on all boundaries: it discards the constraints set by the "channel" flow approach, allowing another CFD flow solution to be the source of flow description at the domain's boundaries, significantly improving the model accuracy at without the need for domain growth. The use of precursor simulations in LES modelling [5] is another example of overcoming the flaws of synthetic boundary conditions, by generating the realistic turbulence required to accurately represent the ABL.

Approach & methods
Unlike linearized RaNS equation solvers, CFD models solving the equations in their full elliptical form are a far more comprehensive form of solving ABL flows, since diffusion and pressure terms are not discarded so information propagation is not prioritized in any single direction. The CFD software VENTOS R /2 used in this work is of that kind, geared towards solving wind flow problems over complex terrain [7,8]. The elliptical RaNS equations are solved using the SIMPLE algorithm [1], with a two equation k-ε model for turbulence closure and the QUICK scheme for the convective term discretization. It uses the "channel" approach to describe the flow domain, with a terrain-following structured mesh, classic inlet and outlet boundary conditions, symmetry condition for the lateral boundaries and a zero-shear condition for the top boundary. Using this approach, while the "channel" approach ensures there is a dominant flow direction, solving the RaNS equations in elliptic formulation ensures any terrain included in the domain is relevant to the global solution, its presence communicated on all directions through pressure and diffusive transport.

Adjusting domain boundaries
Given the advantages of solving the RaNS equation in their elliptical form, the modeller has the temptation of adding to the terrain included in a simulation, increasing the domain. With structured meshes, increasing domain size while maintaining high resolution at the AOI will lead to a rapid increase in node count as the domain grows. By increasing the node count and the size of the linear equation system, the number of unknowns grows with it and the corresponding solver matrix becomes far more rigid, computation times growing exponentially and convergence issues frequently surfacing. So there are significant numerical advantages to keeping domain size and node count reasonable, most importantly greater calculation efficiency. To that end, the AOI solution is prioritized by being solved with a finer mesh resolution than near the domain boundaries, and the domain typically restricted to the least necessary to represent local flow.
On the other hand, as previously highlighted, it is still desirable to include as much complex flow surrounding the AOI as possible, with the drawback of increasing the node count if no loss mesh resolution at the AOI is not an option. To achieve optimal use of available computational resources, it is useful to prioritize the directions in which the domain should be expanded, since moving boundaries away from the AOI will not all yield the same effect on AOI results: • inlet -Traditional boundary layer inlet conditions impose a log-law velocity profile normal to the boundary, with zero transverse velocity components. While suitable to flow over simple topography, in complex cases these conditions can be at odds with the neighbouring boundaries and flow equations. Given that the inlet boundary conditions are fully defined, these disagreements can promote unnatural flow structures by neighbouring boundaries or further downstream into the domain. The dominance of the main flow direction ensures these errors propagate downstream and influence the flow at the AOI, so moving the inlet boundary in such cases is strongly advised. In addition to to this, any complex flow that possibly occurs upstream of the inlet boundary is discarded, moving the boundary to include it will always improve the real-world relevance of the results at the AOI. • outlet -The outlet boundary condition assumes the flow is one-directional and forces a null-gradient solution, so negative velocities and strong transverse flows break the flow assumptions at the outlet. While this type of situation should be avoided, the parabolic nature of channel ABL flows means the flow at the AOI is mostly unaffected by downstream flow behaviour, so unless the boundary condition is forcing an unnatural solution in an area of complex terrain, moving the outlet boundary away from the AOI will give little improvement on the results at the AOI. • lateral -The "channel" flow ensures terrain outside these boundaries is in line in the main flow direction with the AOI, so it should not be a direct influence on the AOI flow. Despite that fact, the lateral boundary condition forces the flow to go tangentially to them so should there be reason to think nearby surface obstacles cause the flow to majorly deviate from that orientation, meaning they are artificially constricting or expanding the flow and the effect propagating inwards towards the AOI, moving them further away from the AOI could have a positive effect on results. • top -The traditional top boundary condition caps the boundary layer, forcing zero vertical velocity and zero gradient at the surface. In very simple terrain the top boundary condition should be placed at the intended boundary layer height, but should the terrain obstacles constrict the flow through the "channel" (artificially accelerating it), it should be displaced further upwards into the atmosphere.
In simple flow cases (over simple topography, ground coverage and any other atmospheric effects) then, expanding the domain in any direction should have negligible effect on AOI results: domain can be kept to a minimum, the AOI plus a small buffer margin between it and the synthetic boundary conditions. On the other hand, usually in complex flows the obstacles and complex phenomena extend far beyond the domain boundaries, so there is no obvious point in any direction after which the flow becomes simple. In those cases the domain should be expanded as far as the available computational resources allow, reducing the AOI flow sensitivity to lateral, top and outlet boundary displacements -which usually requires relatively small changes -and by moving the inlet boundary to include as much upstream complexity as possible, since it will always be relevant to AOI results. One-directional coupling of RaNS domains Exploiting the numerical advantages of smaller linear equation systems and the parabolic nature of ABL flows, we attempt to solve a large flow case using two smaller simulations: one covering the AOI and its close surroundings, another covering as much upstream complex flow as necessary/possible. If the assumptions behind the parabolic flow approximation hold true in the region, coupling those two simulations through the downstream domain's inlet boundary should allow the results of the coupled simulations to closely match those of the equivalent larger domain. Most importantly, solving the case with the two sequentially run small simulations should take far less time than doing so with a single, much larger and numerically cumbersome, simulation covering both AOI and the upstream flow, even if the total node count is comparable or greater in the coupled case. The domain providing the upstream flow solution will be a form of producing an improved description of the flow at the inlet boundary location, which unlike the standard synthetic boundary conditions, takes into account all (or as much as possible, again limited by reason and computational resource availability) the relevant upstream obstructions. This precursor simulation can be made as large or small and its domain mesh as fine or course as necessary, always independently of the target simulation (focusing on the AOI), which then can be optimized in its domain dimensions and mesh to provide the best possible results where they are of most interest. Again this will potentially further compound the computational advantages of such a solution versus solving the whole case in a single domain, where flexibility in allocation of computational resources in space is much more limited.
Coupling the two simulations is done by taking solution data from the concluded precursor simulation (all solved variables apart from pressure, plus velocity gradients) and using it as static inlet boundary condition on the main simulation, interpolating between meshes. Since copied data is from the mean flow field (the part solved by the RaNS equations), the transfer is lossless apart from interpolation errors. To ensure the validity of some approximations and best results, some domain choice limitations apply:

Rødsand II/Nysted wind farm cluster
The modelled flow case sits in the Baltic sea, just off the Danish coast, where two closely spaced offshore WFs (wind farms) operate forming a cluster problem where for the prevailing wind directions one WF operates in the wake of the other [9]. This cluster is formed by: The two WFs are separated by 3 kms in the East-West orientation (see Figure 1) and benifit from a large fetch of undisturbed inflow from the Eastern sector. While under Eastern sector flows, Nysted will be operating in undisturbed inflow but not Rødsand II: its inflow will be partially affected by the wake from Nysted's WTs, an effect that cannot be disregarded considering the proximity of the two WFs. Simulating the operation of both combined implied modelling a cluster spanning more than 20 kms in length and with 162 operating WTs, double the problem size of modelling Rødsand II alone at 12 kms in length and 90 WTs. Using the coupling technique proposed here, a classical simulation covered the upstream flow (the operating Nysted WF) acts as a precursor simulation, providing the necessary improved inlet conditions to model the flow at the target simulation (covering the Rødsand II WF), solving the equivalent cluster flow problem using two smaller mathematical problems, with potential computational resource savings. Between WFs flow will be characterized by the slowly recovering wakes, with very small streamwise gradients and no source of pressure perturbations, meaning the parabolic assumption hold validity and making it an ideal location for the domain coupling to occur.
The VENTOS R /2 software version used here is equipped with a non-linear wake model based on Froude's Actuator Disk concept [10], capable of modelling the operation of multiple WT's iteratively with the RaNS solver [11]. The development of this wake model, along with the data behind proposed offshore flow case, came as part of the work developed under the EU project EERA-DTOC [12,13].

Simulated domains and Coupled runs
A single direction from the Eastern sector was covered, 97 • . It coincides with the alignment of the Nysted WF rows, and guarantees a totally clear inflow to Nysted and a partially waked inflow to Rødsand II. Inflow was calibrated to achieve approximately 6 m/s and 8% turbulence intensity at hub height. Four different domains (detailed in Table 1) covering part or all of the WF cluster were created, serving to simulate the case with several different runs:  • Coupled Low.Res. -as in the Coupled run, but using downgraded quality data for the inlet boundary condition by coupling the target simulation (RSII domain) to a low resolution precursor, which solves flow over Nysted using the NY LR domain. • Regular -a classical approach to modelling the Rødsand II flow without a precursor simulation, using synthetic inlet boundary conditions and a domain covering only the Rødsand II WF (RSII domain), disregarding the effect Nysted has on inflow. The main mesh characteristics -maximum resolution, domain dimensions, mesh expansion factors -were kept equal or close to equal between all runs' domains, and target and cluster domain meshes were tuned to match in the Rødsand II WF. The meshing approach used ensures there is a constant horizontal resolution area covering the totality of the WTs included in the domain at the peak resolution (visible in all three domains in Figure 2), expanding towards the domain boundaries. The maximum vertical resolution is set to a value which sufficiently describes the near wall flow, and since the flow is otherwise a boundary layer over a flat surface, the WT rotor diameter defines the scale of the relevant flow physics: maximum horizontal resolution was set to a quarter of the largest included WT's diameter, which from previous work [11] is thought to be sufficient to accurately describe the momentum effects in the WT vicinity (vertical resolution is always below this inside the rotor span). The NY LR domain differs by halving peak resolution on all directions relative to those of the NY domain. All simulations on the domains of Table 1 were run in steady-state formulation and until residuals on all solved equations dropped below 5 × 10 −4 , where the solution was considered to have converged. Other relevant solver characteristics (relaxation factor, TDMA solver configuration) were also kept the same across all simulations.
The RSII-NY/NY LR domain coupling section (the RSII inlet section) is located between WFs, just downstream of Nysted and outside the constant resolution region of both domains: this area is one of the areas where reduction in node count is possible between Cluster and Coupled runs (observable also by the individual domain node counts in Table 1), maintaining high-resolution domain coverage on the high priority areas (the WFs).

Flow solution
The wind flow description for both wind speed and turbulence intensity across the target inlet ( Figure 3) changes dramatically from synthetic boundary conditions to the full Cluster run solution, and the latter is very closely matched by the Coupled run solution. The distinct wake signatures of Nysted's 9 WT rows are clear, given the short gap between Nysted WF and RSII domain's inlet boundary and the perfect alignment between inflow and WT rows. The vertical flow profiles corroborate that observation, with the synthetic boundary conditions being, inside the rotor span, a very close match to clear inflow conditions on the precursor and cluster solutions. Some flow acceleration is observable outside the WF's wake compared to synthetic boundary condition, observable left and right of the WF signature on the hub height horizontal profiles as well as on the clear inflow vertical profiles in Figure 3, essentially caused    Figure 3: Wind speed and turbulence intensity at the target domain's inlet, across the domain at hub heigh (above) and vertical profiles at free-stream and in the wake of the middle Nysted WT row (below), when sourced from the precursor solutions, full cluster solution or synthetic boundary conditions. Consequence of such similarity between solutions at RSII's inlet section, the overall domain solutions are also very close between Cluster and Coupled runs. In Figure 4 only minute

WF power estimation
Rødsand II WF's power estimation is only very lightly affected by the division of the case into two separate simulations, Table 2 showing that the total WF power increases by only 0.2 % in the Coupled run compared with the Cluster run. This contrasts starkly with the Regular run's result which, unlike Coupled and Cluster runs, does not consider the upstream flow and thus is unable to reproduce the power loss from operating in partly waked inflow, overall power estimation differing by more than 10% to that of the Cluster run. Reducing the quality of the precursor simulation had a clear impact on the upstream WF's total power estimation, but the Rødsand II solution appeared mostly unresponsive to the change, even when the Nysted overall WF power estimation degrades significantly, from -0.2% difference to 6.9%. This goes in agreement with the lack significant change observed in both the inlet profile and inner domain solution when sourcing flow description for the target's inlet boundary from the NY or NY LR domains' solution.

Conclusions
The importance of the upstream flow conditions to a given AOI dictates that the quality of the inflow description should be a priority in any CFD calculation. The approach tested here showed that it is possible reach great improvements over the standard synthetic inlet boundary conditions available to flow modellers, producing an inflow description that reflects upstream conditions at a sensible computational cost. Using the precursor-sourced inlet boundary conditions as form of coupling, results showed that splitting the WF cluster case into two simulations had a limited effect on the flow field and as a consequence on the overall WF power estimation, both on the downstream (Rødsand II) and upstream (Nysted) WFs. There were differences found in individual WT's power estimations, but only on the upstream WF results where cluster and coupled run meshes were non-matching, which raises a known issue of the WT wake model used where WT performance estimation is sensitive to the WT placement within the mesh. More importantly, the aggregate computational cost of the coupled solutions was a fraction of the single domain's. By focusing resources on the domain covering the AOI (Rødsand II), further computational time savings were possible by reducing the quality (and thus cost) of the upstream solution, with negligible impact on the AOI results. The traditional boundary conditions will remain relevant and ultimately still be a limiting factor in CFD modelling, but reducing their impact on results should be a priority of the modeller, and being able to do so with reasonable increase in computational resources is of great use in high complexity cases. The work showed that the principles upon which domain coupling is founded are valid in offshore conditions, where coupling happens in more controlled flow with a single perturbation, in this case in the form of a WF wake. Further work would focus on studying the limits of the technique in onshore cases, where orography as well as heterogeneous surface coverage can introduce added degrees of complexity.