the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
TOSCA – an opensource, finitevolume, largeeddy simulation (LES) environment for wind farm flows
Arjun Ajay
Dries Allaerts
Joshua Brinkerhoff
The growing number and growing size of wind energy projects coupled with the rapid growth in highperformance computing technology are driving researchers toward conducting largescale simulations of the flow field surrounding entire wind farms. This requires highly parallelefficient tools, given the large number of degrees of freedom involved in such simulations, and yields valuable insights into farmscale physical phenomena, such as gravity wave interaction with the wind farm and farm–farm wake interactions. In the current study, we introduce the opensource, finitevolume, largeeddy simulation (LES) code TOSCA (Toolbox fOr Stratified Convective Atmospheres) and demonstrate its capabilities by simulating the flow around a finitesize wind farm immersed in a shallow, conventionally neutral boundary layer (CNBL), ultimately assessing gravitywaveinduced blockage effects. Turbulent inflow conditions are generated using a new hybrid offline–concurrentprecursor method. Velocity is forced with a novel pressure controller that allows us to prescribe a desired average hubheight wind speed while avoiding inertial oscillations above the atmospheric boundary layer (ABL) caused by the Coriolis force, a known problem in wind farm LES studies. Moreover, to eliminate the dependency of the potentialtemperature profile evolution on the code architecture observed in previous studies, we introduce a method that allows us to maintain the mean potentialtemperature profile constant throughout the precursor simulation. Furthermore, we highlight that different codes do not predict the same velocity inside the boundary layer under geostrophic forcing owing to their intrinsically different numerical dissipation. The proposed methodology allows us to reduce such spread by ensuring that inflow conditions produced from different codes feature the same hub wind and thermal stratification, regardless of the adopted precursor run time. Finally, validation of actuator line and disk models, CNBL evolution, and velocity profiles inside a periodic wind farm is also presented to assess TOSCA’s ability to model largescale wind farm flows accurately and with high parallel efficiency.
Please read the corrigendum first before continuing.

Notice on corrigendum
The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.

Article
(7346 KB)

The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.
 Article
(7346 KB)  Fulltext XML
 Corrigendum
 BibTeX
 EndNote
In 2018, Ørsted, a leading company in developing, constructing and operating offshore and onshore wind farms, concluded a project aimed at understanding the limits of models and processes used for wind energy forecasts. The investigation pointed out that blockage and wake effects are currently neglected and underestimated, respectively, when performing wind power predictions (Ørsted, 2019). Blockage, also referred to as turbine or farm induction (Bleeg et al., 2018), is defined as the wind slowdown approaching the wind farm. On the other hand, wake losses are characterized by a power production deficit by waked turbines and are claimed to be underestimated both inside and especially between neighboring sites (Pedersen et al., 2022). While wind farm losses arising from individual turbine wakes have been the subject of extensive research, farm–farm wake effects have gained importance only recently (Lundquist et al., 2019; Ahsbahs et al., 2020; Schneemann et al., 2020). Specifically, as more plants are constructed in the proximity of preexisting ones, the evolution of neighboring farm wakes is an increasingly important aspect to account for and model (Nygaard et al., 2020).
Turbinelevel induction has been researched for many years (Troldborg and Meyer Forsting, 2017; Gribben and Hawkes, 2019; Branlard and Gaunaa, 2014; Branlard et al., 2020), and extensions to wind turbine clusters have been attempted using a linear superposition of individual effects (Branlard and Meyer Forsting, 2020; Segalini, 2021). However, recent studies suggest that this could underestimate – if not totally misrepresent – windfarmlevel blockage, which is heavily influenced by the mutual interaction between the wind farm and the densitystratified atmospheric boundary layer (ABL) (Smith, 2010; Wu and PortéAgel, 2017; Allaerts and Meyers, 2017, 2018, 2019; Centurelli et al., 2021). In fact, the flow deceleration in the wind farm displaces the capping inversion layer, and interfacial waves are formed. Subsequently, their energy is transported vertically and horizontally by atmospheric internal gravity waves. This mechanism triggers pressure disturbances inside the boundary layer, altering the velocity field around the wind farm.
In industry, annual energy captures are made using lowcost but fast, often analytical, reducedorder wake models (Jensen, 1983; Ainslie, 1988; Larsen, 1988; Bastankhah and PortéAgel, 2014; Niayifar and PortéAgel, 2016), aimed at capturing the gross aerodynamic processes within the farm. While they have been used effectively for hundreds of wind energy projects, the majority of these models currently struggle in accurately reproducing wind farm blockage and farm–farm wake interactions (Nygaard et al., 2022). This is classified as an industrywide issue, as overpredicting annual energy production can have a negative impact on all companies' financial estimates.
Reducedorder models need to be thoroughly validated, but comprehensive observation datasets are difficult to obtain. Numerical analyses, in particular largeeddy simulations (LESs), are able to provide such data, together with valuable insight into the physical processes. LES resolves the largest and most energetic turbulent eddies, while the smallest ones are modeled. Nevertheless, LES of large wind farms in the atmospheric boundary layer (ABL) is extremely challenging, given the breadth of scales involved, spanning from resolved turbulence eddies of a few meters to gravity waves characterized by wavelengths of several kilometers. To successfully tackle these problems, LES solvers must possess good parallel efficiency and an optimized code input–output (I/O), allowing them to operate on modern highperformance computing (HPC) architectures. Additionally, many numerical aspects have to be carefully treated, such as the coupling with mesoscale models (Haupt et al., 2023), handling wave reflections produced by the domain boundaries or the generation of a suitable timeresolved turbulent inflow (Lanzilao and Meyers, 2022a). The last two tasks can be achieved simultaneously using the concurrentprecursor method, where a simulation without wind turbines (precursor) is advanced in sync with the wind farm simulation (successor). The latter features a fringe region, where body forces are used to damp gravity wave reflections and to restore the desired turbulent inflow. At each time step, such body forces are calculated based on the concurrentprecursor instantaneous fields, leading to the precursor and successor solutions matching at the fringe region exit. More details on precursor techniques are given in Sect. 2.4, where our new hybrid method is also described.
Several LES codes have been developed by the research community so far (see Breton et al., 2017, for a review), among which only a few can effectively tackle the abovementioned application. The KU Leuven code SPWind, for example, has been successfully used for finite wind farm simulations capturing gravity wave effects (Lanzilao and Meyers, 2022b) but unfortunately is not opensource. Conversely, opensource tools, such as the Parallelized LargeEddy Simulation Model (PALM; Maronga et al., 2015), developed by the Institute of Meteorology and Climatology at Leibniz Universität Hannover (Germany), or SOWFA (the Simulator fOr Wind Farm Applications), maintained by the National Renewable Energy Laboratory (NREL), do not implement the concurrentprecursor method, making it difficult to properly simulate gravity wave effects at the same time as avoiding inlet–outlet reflections. In addition, although SOWFA has been used in several research studies in the last decade (Churchfield et al., 2012b, a; Fleming et al., 2014; Johlas et al., 2021, to name a few), it is based on OpenFOAM 6 (OpenCFD, 2018), a generalpurpose set of libraries that are not specifically designed to efficiently run at scale. SOWFA's greatest limitation is its massive generation of output files. In particular, a directory containing all simulation fields is generated for each processor at run time. When dealing with thousands of processors, the number of produced files can easily saturate, within a few checkpoint iterations, the maximum file count in many HPC architectures. To address these shortcomings, NREL has started the ExaWind project (Min et al., 2022), but the latter is not yet at a final production stage, and it does not feature the concurrentprecursor technique or a method to model complex terrains.
For the aforementioned reasons, we have developed an opensource, finitevolume framework that is tailored for largescale studies of windfarminduced gravity waves and cluster wake–atmosphere interaction, with the objective of gaining sufficient understanding of the physics of atmospheric flow within and around wind plants, a grand challenge of modern wind energy according to Shaw et al. (2022). The new framework is called TOSCA (Toolbox fOr Stratified Convective Atmospheres) and exploits stateoftheart parallel libraries, such as Open MPI (Gabriel et al., 2004), PETSc (Balay et al., 2022), HYPRE (Falgout and Yang, 2002) and HDF5 (The HDF Group, 2006), for the parallel solution of partial differential equations and handling of intense I/O operations. TOSCA is specifically designed to enable LESs of large finite wind farms. Wind turbines can be modeled using the actuator line (Sørensen and Shen, 2002) and the actuator disk (Jimenez et al., 2007, 2008) models. As inlet–outlet boundary conditions produce a consistent and undesirable reflection of atmospheric gravity waves, we introduce a hybrid offline precursor–concurrentprecursor methodology which, coupled with periodic boundary conditions, limits artificial wave reflections while simultaneously reducing the computational cost associated with initializing the turbulent precursor. The concurrentprecursor method (Inoue et al., 2014) is to our knowledge not available in other finitevolume solvers, though it is extensively used in pseudospectral methods (Wu and PortéAgel, 2017; Allaerts and Meyers, 2017, 2018). For this reason, gravity wave studies to date have been only performed using the latter discretization technique, which does not allow for grid refinement in the pseudospectral directions. This forces a uniform grid resolution, leading to high cell counts. Conversely, the finitevolume method allows for grid stretching, enabling us to resolve larger domains with the same number of degrees of freedom while providing greater geometrical flexibility. Finally, TOSCA also features a sharpinterface immersed boundary method (IBM) based on Haji Mohammadi et al. (2019) that allows us to simulate moving objects and complex terrain features, but its validation will be covered in a followup paper.
The present paper is organized as follows. First, in Sect. 2, we describe the developed LES framework. Next, Sect. 3 presents comparisons with existing numerical and experimental studies to validate TOSCA's actuator models, the evolution of thermally stratified ABLs and wake interactions inside a periodic wind farm in neutral conditions. In Sect. 4, we compare results obtained from CNBL simulations using the newly developed velocity and temperature controlling techniques against the commonly used geostrophic forcing combined with a wind angle controller. In Sect. 5, we present the simulated flow field around a reference 100turbine finite wind farm immersed in a turbulent CNBL, highlighting TOSCA's ability to accurately predict gravity wave blockage effects. Finally, conclusions are outlined in Sect. 6.
TOSCA is a finitevolume code, formulated in generalized curvilinear coordinates, allowing it to also take as input nonCartesian structured meshes. The present section is organized as follows. The governing equations in Cartesian coordinates are reported in Sect. 2.1, while actuator models used to represent wind turbines in the domain are described in Sect. 2.2.
To provide a better flow of the paper, the numerical method, the governing equations in curvilinear coordinates actually solved in TOSCA and a brief overview of generalized curvilinear coordinates are reported in Appendix A, while the LES turbulence model in the curvilinear frame is detailed in Appendix B.
An overview of TOSCA's parallel efficiency is given in Appendix D, where we analyze the time per iteration with an increasing number of nodes and mesh elements on the Niagara (Niagara, 2018; Loken et al., 2010; Ponce et al., 2019) highperformance computer at the SciNet HPC Consortium. In addition, TOSCA has been used to run finite wind farm simulations on the whole Niagara cluster (2024 nodes, 40 cores per node) and on all Cascade nodes of the UBC ARC (Sockeye, 2023; 170 nodes, 40 cores per node) cluster, demonstrating its capability to handle massively parallel computations.
In order to run ABL simulations, we developed a novel methodology, described in Sect. 2.3, that enforces a desired hubheight wind speed while simultaneously avoiding inertial oscillations produced by the Coriolis force above the boundary layer. In addition, we show that disagreement exists between different computational fluid dynamics (CFD) codes in predicting the final mean potentialtemperature profile inside the boundary layer. In this regard, we propose the use of a mean temperature controller which maintains a prescribed average potentialtemperature profile, harmonizing the comparison of simulation results in future studies. Finally, Sect. 2.4 details our hybrid offline–concurrentprecursor methodology, which saves computational resources when performing the turbulence initialization in the precursor phase.
2.1 Governing equations
Governing equations correspond to mass and momentum conservation for an incompressible flow with Coriolis forces and Boussinesq approximation for the buoyancy term. The latter is calculated using the modified density ρ_{k}, evaluated by solving a transport equation for the potential temperature. These equations, expressed in Cartesian coordinates using tensor notation, read
where u_{i} is the Cartesian velocity; $p/{\mathit{\rho}}_{\mathrm{0}}$ is the kinematic pressure; θ is the potential temperature, defined as $\mathit{\theta}=T({p}_{\mathrm{0}}/p{)}^{R/{c}_{\mathrm{p}}}$ (T is the absolute temperature, R is the specific gas constant for dry air, c_{p} is the specific heat at constant pressure and p_{0} is the reference pressure); g_{i} is the gravitational acceleration vector; and Ω_{j} is the rotation rate vector at an arbitrary location on the planetary surface (defined as $\mathit{\omega}\mathrm{cos}\mathit{\varphi}\widehat{y}+\mathit{\omega}\mathrm{sin}\mathit{\varphi}\widehat{z}$, where ϕ is the latitude, in a local reference frame having $\widehat{z}$ aligned and opposite to the gravitational acceleration vector, $\widehat{x}$ tangent to Earth's parallels, and $\widehat{y}$ such that the frame is righthanded). Source terms f_{i}, ${s}_{i}^{v}$ and ${s}_{i}^{h}$ are body forces introduced by turbines and by vertical and horizontal damping regions, respectively. Moreover, the modified density ρ_{k} is defined as
where θ_{0} is a reference potential temperature, chosen as the ground temperature. Parameters ν_{eff} and κ_{eff} are the effective viscosity and thermal diffusivity, respectively. The former is the sum of the kinematic viscosity ν and the subgridscale viscosity ν_{t}, while the latter is sum between the thermal diffusivity $\mathit{\kappa}=\mathit{\nu}/Pr$ and the turbulent thermal diffusivity κ_{t}. Both ν_{t} and κ_{t} are defined in Appendix B, while the Prandtl number Pr is set to 0.7 in all simulations. The third term on the righthand side of Eq. (2) is a uniform horizontal pressure gradient that balances turbulent stresses and the Coriolis force, allowing the boundary layer to reach a statistically steady state. This term is commonly referred to as the velocity controller, and it is explained in Sect. 2.3.1.
2.2 Actuator models
To represent wind turbines, different models have been implemented. In TOSCA, they are referred to as the actuator line model (ALM), actuator disk model (ADM) and uniform actuator disk (UAD) model. Following Sørensen and Shen (2002), Sørensen et al. (2015), PortéAgel et al. (2010), and Jimenez et al. (2007), the first two models require detailed blade information (i.e., airfoils, twist and chord), while the UAD only requires the turbine thrust coefficient and general rotor information such as diameter and hub height (Jimenez et al., 2007, 2008). The idea behind actuator models is to represent the wind turbine as a distribution of points, each associated with a Lagrangian force. For the UAD model, the sum of forces from all points must be equal to the total wind turbine thrust, while the ALM and ADM in TOSCA additionally include rotor torque, as they also model blade rotation. Once the Lagrangian force at each point has been calculated, it is distributed to the surrounding mesh cells through a projection function. In TOSCA, a classical isotropic Gaussian projection is used, namely
where (x_{0}, y_{0}, z_{0}) is the position of the actuator point and ϵ is a tunable parameter, corresponding to the standard deviation of the Gaussian projection function. Note that while the projection function should integrate to unity to preserve each point force, this is never exactly possible, and the projection distance is cut when 99 % of the Gaussian volume has been taken into account. Moreover, the Gaussian widthtogrid size ratio should be larger than 2 in order to avoid large projection errors and numerical instabilities (MartínezTossas et al., 2015).
The definition of the turbine point mesh and the evaluation of the point force are different depending on the specific model. In the ALM, each rotor blade is represented by a line of points, which are physically rotated at each iteration, making it an unsteady model. In the ADM and UAD model, the number of points in the azimuthal (tangential) direction is not equal to the number of blades, and it is usually set to a high value. For both the ALM and ADM, the point force is calculated exploiting the blade element momentum theory (BEM; see Glauert, 1935). First, the radially varying velocity is estimated at each point, using information from the CFD mesh and the wind turbine angular velocity. This is known as velocity sampling, and different methods have been proposed (Churchfield et al., 2017). TOSCA samples the velocity at the actuator point, using nearestneighbor interpolation from the closest mesh cell. Next, the velocity magnitude and angle of attack are given as input to appropriate airfoil tables, which return lift and drag at the point location. Various airfoil tables are used along the blade radius because the airfoil type usually changes along the blade span, as does the operating Reynolds number. Lift and drag at each actuator point are then distributed to the surrounding CFD cells by convolution with the projection function. Conversely, for the UAD model, the blade loading is uniform, and the force is calculated by dimensionalizing the turbine thrust coefficient with the freestream velocity and the portion of rotor area belonging to each actuator point. In waked conditions, the concept of freestream velocity is not well defined. Hence, a common practice is to first average the wind velocity on the rotor disk and then use the momentum theory to infer the corresponding freestream velocity (Meyers and Meneveau, 2010). In our framework, since the ADM and ALM also account for blade rotation, they can be coupled with a rotor inertia and control system dynamics solver (pitch and angular velocity controllers), while nacelle yaw can be applied to any of the three models.
2.3 Controllers
This section reviews the current state of the art for velocity controllers in precursor simulations and presents a novel technique, which we refer to as geostrophic damping, which allows control of the hubheight velocity while avoiding inertial oscillations generated by the Coriolis force. Moreover, a simple temperature controller is also presented that maintains a constant average potentialtemperature profile throughout the precursor simulation.
2.3.1 Velocity controller
In the precursor simulation, the flow is usually driven by a uniform horizontal pressure gradient, which is related to the geostrophic wind components by the geostrophic balance at equilibrium:
where f_{c}=2Ω_{z} is known as the Coriolis parameter. Using the above equations to prescribe the driving pressure gradient does not give any control over velocity magnitude and direction at the wind turbine hub height. In fact, the latter will be a result of the turbulent stresses inside the boundary layer, which are not known a priori. However, being able to control these parameters is convenient in wind farm simulations as it allows the operation point of the turbines to be easily prescribed. To this end, Sescu and Meneveau (2014) and Allaerts and Meyers (2015) developed and tuned an algorithm that slowly rotates the flow in the domain, allowing us to control the wind direction at a specified height. Later, Stieren et al. (2021) used the same approach to impose dynamic wind direction changes. Besides the driving pressure gradient, evaluated using Eq. (6), the additional cross product ${\mathit{\u03f5}}_{ijk}{\mathit{\omega}}_{i}{u}_{j}{\widehat{x}}_{k}$ is added to the momentum equation's righthand side, where the angular frequency ω_{i} is calculated based on the angle difference at the reference height (see Allaerts and Meyers, 2015, for details on this procedure). Such a method, which we will refer to as the geostrophic controller, does not entirely solve the issue, as velocity magnitude at the hub height is still unknown a priori. Nevertheless, a different approach exists, available for example in SOWFA, which allows us to prescribe both magnitude and direction at a specified height h_{ref}. In particular, given a desired velocity u_{ref,i}, which should be maintained at h_{ref}, an error vector can be defined as the difference between the reference wind and the velocity sampled at the reference height, averaged over the homogeneous directions. At this point, a proportional–integral controller can be used to evaluate the driving pressure gradient (i.e., the third term on the righthandside of Eq. (2) such that the desired speed and angle are maintained at h_{ref}. This approach will be referred to as the pressure controller. In TOSCA both controller methods are implemented, and the driving pressure gradient in the second type of controller is evaluated as
where subscript i refers to the ith component, e_{P,i} is the proportional error, ${e}_{I,i}^{n}$ is the integral error evaluated at time step n, r is a relaxation factor, α is the proportional fraction of the controlling action, T is the time filter for the integral error, Δt is the time step size, and 〈⋅〉_{xy} denotes a spatial average along the homogeneous directions x and y. In the present study, we set r=0.7, α=0.8 and T=2 h.
On one hand, the pressure controller is more convenient for wind turbine simulations, as hub wind and direction can be directly specified. However, unlike the geostrophic controller, it does not provide knowledge of the geostrophic wind a priori, making it impossible to initialize the flow such that Eq. (6) is satisfied. An inconsistency in the initial condition produces inertial oscillations above the boundary layer, as the initial wind speed aloft differs from its equilibrium geostrophic value. This can be easily verified by noting that the unsteady form of Eq. (6) (see for example Stull, 2016), namely
represents an undamped linear oscillator with angular frequency f_{c}. In particular, if v≠V_{G} or u≠U_{G}, at any point during the simulation, inertial oscillations will be produced. In some wind energy applications, for example, when studying the formation of atmospheric gravity waves above the boundary layer, the physics of the problem strongly depends on the magnitude of the geostrophic wind. In such cases, results would be negatively impacted by these inertial oscillations, whose amplitude depends on the initial condition.
Nevertheless, being able to exactly define the wind speed and direction at a specified height is a desirable property of the pressure controller. For this reason, we developed a new methodology that allows us to remove these inertial oscillations, also enabling the use of the pressure controller in those cases where a steadystate geostrophic wind is preferred. First, we note that the system of Eq. (10) can be damped by introducing an additional term as follows:
where the coefficient α determines if the system is overdamped (α > 1), underdamped (α < 1) or critically damped α=1. With some manipulation, Eq. (11) can be rewritten as
These equations slightly differ from a conventional spring–mass–damper system in the additional coupling terms $\mathrm{2}\mathit{\alpha}{f}_{\mathrm{c}}^{\mathrm{2}}(v{V}_{G})$ and $\mathrm{2}\mathit{\alpha}{f}_{\mathrm{c}}^{\mathrm{2}}(u{U}_{G})$.
We observed that the presence of these terms enhances the damping action, halving the exponent of the decay rate, characterized by an efolding time of $\mathrm{1}/\left(\mathrm{2}\mathit{\alpha}{f}_{\mathrm{c}}\right)$. In order for the oscillation amplitude to reach less than 3 % of the initial value, a damping time ${T}_{\mathrm{3}\phantom{\rule{0.125em}{0ex}}\mathit{\%}}=\mathrm{ln}(\mathrm{100}/\mathrm{3})/\left(\mathrm{2}\mathit{\alpha}{f}_{\mathrm{c}}\right)$ is necessary. Note that, in order for the damping term to be evaluated, knowledge about the geostrophic wind components is still required. We deduce U_{G} and V_{G} from the driving pressure gradients imposed by the pressure controller (Eqs. 7, 8, 9) by means of the definition of the geostrophic wind speed (i.e., the geostrophic balance given by Eq. 6). In addition, we filter the obtained geostrophic components using a filter constant of $\mathrm{0.2}\mathit{\pi}/{f}_{\mathrm{c}}$ (corresponding to $\mathrm{1}/\mathrm{10}$ of the inertial oscillation period). Finally, we highlight that the damping action should start after the boundary layer has become fully developed, as the pressure gradient prescribed by the controller depends on the turbulent stresses if h_{ref} is inside the ABL. In our simulations, we start the damping action after almost one inertial period (${T}_{\mathrm{D}}\approx \mathrm{2}\mathit{\pi}/{f}_{\mathrm{c}}$), and we maintain the damping active for at least a time equal to T_{3 %}. In order for this geostrophic damping method not to affect the velocity inside the ABL, we smoothly bring the damping to zero below a certain height by multiplying the damping terms with the following function:
where H_{d} is the height where the damping has halved its strength. If the simulation models a capping inversion layer, we set H_{d}=H, where H is the capping inversion center, and Δ_{d}=Δ, where Δ is set as the capping inversion width.
2.3.2 Temperature controller
When running precursor CNBL simulations, the predicted ABL height, as well as the final value of potential temperature at the ground, depends on the mixing history experienced inside the boundary layer. This is in turn affected by the specific LES setup and the type of discretization used. Moreover, the impact of such code details is made even more noticeable by the fact that these simulations usually run for a very long time (of the order of $\mathrm{2}\mathit{\pi}/{f}_{\mathrm{c}}$). For example, SPWind, which employs a pseudospectral discretization in the horizontal directions and an energyconservative fourthorder advection scheme in the vertical, predicts less mixing than other pseudospectral codes, such as NCARLES (Pedersen et al., 2014) or WiRELES (Abkar, M. and PortéAgel, 2013), which use for example a secondorder central scheme in the vertical direction (this can be appreciated in Fig. 3). Besides, in finitevolume codes, like TOSCA or SOWFA, an upwindbiased advection scheme is usually preferred, as it stabilizes the numerical method but does not allow us to conserve mechanical energy. These considerations pose comparison issues among different codes, as their differences will have an impact for example on the final ABL height, inversion thickness and potentialtemperature jump and, in general, on the heating history of the boundary layer, ultimately affecting the successor solution, in which wind turbines are present. In particular, if a CNBL precursor is run with a certain initial potentialtemperature profile, this will evolve differently based on both the adopted simulation framework and the length of the precursor run, determining a discrepancy in the initial condition of the wind farm simulation.
In the present work, we propose to apply a potentialtemperature controller in the precursor simulation so that the successor can be exactly run with the intended temperature profile. This could be beneficial for example when making comparisons between different codes, as it ensures that the successor background potentialtemperature profile matches the precursor initial condition. In particular, we apply the following heightdependent source term on the righthand side of Eq. (3):
where $\overline{\mathit{\theta}}\left(h\right)$ is the desired vertical potentialtemperature profile, taken as the initial value of 〈θ(h)〉_{xy}, which is average of the potential temperature along the homogeneous directions at a given height. The parameter r is a relaxation coefficient that we set to 0.7. Note that a very similar method was used by Allaerts et al. (2020, 2023) to drive LESs using mesoscale models or observations.
2.4 Hybrid offline–concurrent precursor
Wind turbine wake recovery and thus power production are greatly influenced by background atmospheric turbulence. As a consequence, prescribing a physical turbulent inflow is necessary if real wind turbine operation is to be simulated. A commonly used approach is the socalled precursor–successor method (Churchfield et al., 2012b, a), where a first simulation of the sole ABL, without wind turbines, is run until turbulence reaches steadystate statistics. After this first phase, the latter is further progressed, and velocity and potential temperature are saved on a plane parallel to the inlet boundary, at each iteration, forming the inflow database. At this point, the simulation with the wind turbines (successor) is started and, at each time step, the inflow boundary condition is interpolated from the saved slices at the two closest times in the database so that precursor and successor time steps can take different values. The above methodology implies no periodicity of the domain in the streamwise direction and has proved to work extremely well for isolated turbine simulations, cases where the ABL height is not perturbed by objects located below it or in the absence of thermal stratification. On the contrary, when thermal stratification is present, atmospheric gravity waves can be triggered above the ABL, and a careful design of the simulation should prevent such waves from being reflected by the physical boundaries. In particular, the LES setup should be equipped with damping regions at the top, inlet and outlet, or if streamwise periodic boundary conditions are used, only one damping region in the streamwise direction is then required (Calaf et al., 2010; Inoue et al., 2014; Allaerts and Meyers, 2017). The latter, also known as the fringe region, must ensure that turbine wakes are prevented from being reintroduced into the domain by periodic boundaries and that an unperturbed turbulent inflow is reached again at the fringe exit. To achieve this, the desired flow that is used to compute the damping term should ideally contain timeresolved turbulent structures at every cell located in the fringe. The precursor is then advanced in sync with the successor, in a domain larger than or equal to the fringe region, so that velocity and temperature fields are available, at each time step and spatial location, to compute damping sources.
In order to spin up the wind farm simulation, we developed a threestep procedure (see Fig. 1), where both offline and concurrentprecursor methodologies are adopted. We believe that the latter method is necessary when dealing with windfarminduced gravity waves, as it allows us to avoid wave reflections while prescribing a timeresolved turbulent inflow at the same time. Since the concurrentprecursor domain should coincide with the successor – whose size is dictated by the gravity waves and wind farm – in both the spanwise and the vertical directions, it is usually oversized from the turbulence generation point of view. This leads to a considerable quantity of computational resources being consumed when starting up the unperturbed ABL. In fact, domains of a much smaller size are used in the literature when the sole ABL is of interest or when the inflow data are generated using the offline precursor technique.
In TOSCA, we exploit the flexibility of the finitevolume formulation by combining the two techniques. In particular, we initialize the ABL on what we refer to as the offline precursor domain, which can be arbitrarily defined in both the streamwise and the vertical directions. The only requirement is that the spanwise size of the successor is an exact multiple of the offline precursor domain. When turbulence has reached a statistically steady state, we save flow slices of velocity and potential temperature from this domain into a socalled inflow database. After this first phase, the concurrentprecursor and successor simulations are started, and streamwise inflow–outflow boundary conditions are used in the former for one flowthrough time. Inflow slices from the inflow database are periodized in the spanwise direction, while extrapolation is performed in the vertical direction. We note that it is extremely important that the flow above the inversion layer does not contain any periodic variations in time, as this would be noticed in the successor, at streamwise intervals equal to the concurrentprecursor domain length. For this reason, we average the offline precursor data at the 10 highest cells and slowly merge the instantaneous data to this average across such an interval. This removes even the smallest periodic content in the flow above the boundary layer, which is now characterized by a truly constant geostrophic wind. For instance, we weight the average and instantaneous velocities using two hyperbolic weighting functions (Eq. 13 is used for both, but a minus sign after unity is applied for the instantaneous velocity). H_{d} is set to the height of the fifth cell center from the offline precursor top boundary, while Δ_{d} is equal to the width of the 10 highest cells.
After the concurrent precursor has run for one flowthrough time, the turbulent inflow has reached the outlet, and streamwise boundary conditions are switched to periodic. At this point, the simulation is selfsustained, and we run the precursor and successor simultaneously for one successor flowthrough time so that gravity waves and wind turbine wakes are formed. For the simulation presented in Sect. 5, we ran the offline precursor for 10^{5} s, the concurrentprecursor spinup phase where we used inflow–outflow boundary conditions for 600 s, and the overall successor spinup for 5000 s (note that this phase can start in parallel with the previous one). Data were gathered from 105 000 to 120 000 s.
The drawback of such a method is that a spanwise periodicity is introduced into the concurrentprecursor and successor domains. If the Coriolis force is active, this can be broken everywhere except at hub height, where the flow is aligned with the x axis, by setting the offline precursor and the concurrentprecursor streamwise domain length to a different value. At the hub height, the larger turbulence structures might be locked in position if they span the whole domain length. Although they will eventually disappear, they result in slow convergence of flow averages at the hub height. This issue has been already observed in the past for example by Munters et al. (2016), who proposed to use shifted periodic boundary conditions in the concurrentprecursor domain.
Nevertheless, the proposed hybrid method, sketched in Fig. 1, is very convenient as it allows us to reduce the overall computational cost of the ABL spinup phase, where windfarminduced gravity waves are not yet present. In fact, this initial phase is run on a domain whose size is dictated by the current flow physics, rather than on quantities that will only become relevant at later simulation stages.
In this section, we validate the developed solver using three different benchmark cases. In Sect. 3.1, we simulate an NREL 5 MW reference wind turbine, operating in a uniform inflow equal to 8 m s^{−1}, and compare our results to those of MartínezTossas et al. (2015). In Sect. 3.2 we validate the ability of TOSCA to simulate conventionally neutral boundary layer (CNBL) evolution, comparing our results against data from different LES codes reported by Allaerts (2016). Finally, in Sect. 3.3, an infinite wind farm in a turbulent boundary layer without thermal stratification is compared to experimental and numerical data collected by Chamorro and PortéAgel (2011) and Stevens et al. (2018), respectively.
3.1 Isolated rotor in uniform inflow
In this validation case, we perform two simulations of the NREL 5 MW reference wind turbine (Jonkman et al., 2009) using the ADM and the ALM techniques, with a uniform inflow velocity of 8 m s^{−1}. Periodic boundary conditions are applied at the upper, lower and spanwise boundaries. At the outlet, a zero normal gradient on velocity outflow is specified. The domain is 10 rotor diameters in all directions, with the turbine rotor placed in the geometric center of the domain. The mesh is graded in all directions from a resolution of 16.8 m next to all boundaries to 2.1 m near the wind turbine. In particular, this fine region, where the mesh is uniform, extends 1 diameter upstream of the turbine and 5 diameters downstream. In the other two directions, it extends beyond the edge of the rotor for 1 diameter.
For this case, the standard Smagorinsky model was used, where we set the C_{s} coefficient of Eq. (B6) to 0.028224 (corresponding to the value of ${c}_{s}=\sqrt{{C}_{\mathrm{s}}}=\mathrm{0.168}$ used in MartínezTossas et al., 2015). The ALM has 63 points in the radial direction, while the ADM has 63 and 72 points in the radial and azimuthal directions, respectively. The rotational speed of the wind turbine is set to a constant value of 9.1552 rpm in all cases. The projection ϵ in Eq. (5) is set to 4.2 m. Both simulations are advanced for 300 s, after which data are averaged for the next 300 s. Figure 2 shows the normalized wind speed deficit at 1 and 4 downstream rotor diameters for ALM and ADM simulations performed both with TOSCA and by MartínezTossas et al. (2015). An excellent match can be observed at 1 diameter for both models, while at 4 diameters TOSCA predicts sightly higher deficits, especially for the ALM. This difference is due to an earlier breakdown of the bladetip vortices in the simulations of MartínezTossas et al. (2015), which were performed with OpenFOAM. As OpenFOAM is an unstructured code, we believe that nonhexahedral elements arising from the three successive refinement regions produce some small oscillations in the velocity, which is seen by the simulation as added turbulence intensity, determining an earlier breakdown of the bladetip vortices. This effect is not present in TOSCA, as the mesh is fully structured and smoothly graded from 16.8 to 2.1 m. In van der Laan et al. (2014), the same case is run without a turbulence model using EllipSys3D and SnS, and a higher maximum deficit than both TOSCA and MartínezTossas et al. (2015) has been observed at 2.5 diameters. Such a discrepancy between different codes in predicting turbine wake recovery is not observed when a precursor is used to prescribe the inflow, as wake mixing is guided by ABL turbulence instead of numerical oscillations. In Table 1 we report the aerodynamic power produced by the wind turbine as predicted by TOSCA with the two actuator models and that obtained by MartínezTossas et al. (2015).
The ADM matches well with results from MartínezTossas et al. (2015), while TOSCA's ALM predicts a slightly higher power. The reason for such a difference is presently unknown to the authors.
3.2 CNBL evolution
In this section, we validate TOSCA's ability to perform CNBL simulations by running two cases of a neutral boundary layer, developing against a stable background stratification with lapse rate of 1 and 10 K km^{−1}. In both validation cases, the geostrophic wind is G=10 m s^{−1}, the Coriolis parameter is set to ${f}_{\mathrm{c}}={\mathrm{10}}^{\mathrm{4}}$ s^{−1}, the surface roughness is z_{0}=0.01 m and the reference temperature is 290 K. The numerical domain size is 3 km × 3 km × 2 km, with 256^{3} grid points. This corresponds to a grid resolution of approximately 11.7 m × 11.7 m × 7.8 m. Results are compared with results from SPWind (Allaerts, 2016), WiRELES (Abkar, M. and PortéAgel, 2013) and NCARLES (Pedersen et al., 2014), which all use a pseudospectral horizontal discretization.
Velocity and temperature fields are initialized with a constant and linear profile equal to the geostrophic velocity and background stratification, respectively. Furthermore, sinusoidal perturbations are added to the velocity profile, below 100 m, with an amplitude of 0.1G and 12 periods in the x and y directions to trigger turbulent fluctuations. Simulations are advanced in time for 24 h, and results are averaged over the last hour. Periodic boundary conditions are applied in the horizontal directions, while a slip boundary condition is applied at the upper boundary. At the ground, the wall shear stress is prescribed through classic Monin–Obukhov similarity laws (Monin and Obukhov, 1954; Paulson, 1970; Etling, 1996), following the approach of Yang et al. (2017) to address the loglayer mismatch. The flow is driven using pressure gradients obtained from Eq. (6). We do not apply any velocity controller in these simulations so that the wind direction inside the boundary layer is free to change while the geostrophic wind remains aligned with the x axis. Figure 3 compares vertical profiles of velocity magnitude, horizontal wind direction, potential temperature and kinematic heat flux obtained from TOSCA, with profiles reported by Allaerts (2016), Abkar, M. and PortéAgel (2013) and Pedersen et al. (2014). Very good agreement is found in the horizontal wind direction and magnitude. Regarding temperature profiles, TOSCA is more aligned with NCARLES and WiRELES results, while SPWind predicts a slightly lower inversion layer and potential temperature at the ground. This highlights how turbulent mixing is predicted differently by the four codes. Heat flux profiles agree well below the inversion layer, with the exception that TOSCA predicts a more diffused kinematic heat flux profile above the inversion layer for the 10 K km^{−1} case. We do not have a clear explanation for such behavior.
In Table 2 quantitative parameters of the resulting ABLs are reported for the different cases and LES codes. The reference temperature θ_{0}, the capping inversion strength Δθ and the inversion width Δh are evaluated by a leastsquares fit of the resulting temperature profiles with the model proposed by Rampanelli and Zardi (2004). The ABL height is taken as the center of the capping inversion layer. Both vertical profiles of Fig. 3 and quantitative ABL parameters reported in Table 2 demonstrate that TOSCA is well aligned with results from Pedersen et al. (2014) and Abkar, M. and PortéAgel (2013) and is thus capable of conducting CNBL simulations.
In addition, we note that simulations performed using SPWind consistently predict less mixing than other codes, which is confirmed by the lower absolute value of the minimum kinematic heat flux q_{min}. This could be due to the fourthorder energyconservative scheme, which is adopted in SPWind simulations, while other codes employ second or thirdorder nonconservative advection schemes.
3.3 Infinite wind farm in neutral conditions
In this section, we run the same infinite wind farm simulation as has been conducted in Stevens et al. (2018), corresponding to the wind tunnel experiments performed by Chamorro and PortéAgel (2011). The scaled wind farm consists of 30 wind turbines, arranged in an aligned configuration with 3 columns and 10 rows. Spanwise and streamwise spacings are set to S_{y}=4 D and S_{x}=5 D, respectively, where D=0.15 m is the turbine diameter. The wind farm is made periodic in the spanwise direction by placing turbine columns 1 and 3 at a distance of ${S}_{y}/\mathrm{2}$ from the lateral boundaries, where periodic boundary conditions are applied. In Stevens et al. (2018), simulations are run with both the ALM and the nonrotating uniform ADM model, which we refer to as the uniform actuator disk model (UADM), and each wind farm row has a different C_{t} value. For this validation case, we did not attempt to use rotating actuator models (ADM or ALM), as C_{t} coefficients applied in Stevens et al. (2018) at some rows are higher than the value of C_{t,max} from their reported BEM calculations. Therefore, since it would not have been possible to match their exact angular velocity for some of the rows, we decided to opt for the UADM, where turbinespecific thrust calculation at the pth disk element is solely based on the thrust coefficient C_{t} and the freestream velocity U_{∞} as
In the above expression, dA_{p} is the disk area associated with the pth actuator disk point; ${\widehat{\mathit{e}}}_{t}$ is a vector normal to the rotor disk, pointing in the upstream direction; and U_{∞} is evaluated from the average disk velocity U_{disk} exploiting the momentum theory. For instance, ${U}_{\mathrm{\infty}}={U}_{\mathrm{disk}}/(\mathrm{1}a)$, where a is the induction factor, related to the thrust coefficient as ${C}_{t}=\mathrm{4}a(\mathrm{1}a)$. Note that Eq. (15) can be rewritten in an equivalent form by using U_{disk} and the diskbased thrust coefficient ${C}_{t}^{\prime}={C}_{t}/(\mathrm{1}a{)}^{\mathrm{2}}$ in place of U_{∞} and C_{t}. We use the latter formulation in the present case, as Stevens et al. (2018) reported the value of ${C}_{t}^{\prime}$ at each wind farm row.
The wind tunnel model used in the experiment by Chamorro and PortéAgel (2011) is the GWS/EP6030 turbine; it has a hub height of 0.125 m, an overhang of 0.03 m, a hub radius of 0.0075 m and a tower diameter of 0.01 m. The tower and nacelle have been modeled following the approach of Stevens et al. (2018), except for the projection function, which is given by Eq. (5). The value of ϵ has been set to 0.02452 for the tower and nacelle and to 0.03515625 for the rotor in order to closely match their approach. The tower is represented by 50 actuator points and is characterized by a drag coefficient of 0.68, while the nacelle consists of a single point where the force is calculated by dimensionalizing a drag coefficient of 4. Moreover, the rotor has been discretized using 20 radial points and 50 azimuthal points. To prescribe a turbulent inflow, Stevens et al. (2018) used the concurrentprecursor technique, as their code is pseudospectral. Since TOSCA is a finitevolume code and inflow–outflow conditions can be applied, we opted for the computationally cheaper offline precursor technique described in Sect. 2.4. The precursor domain is 1.8 m × 1.8 m × 0.675 m, with 129 × 129 × 145 cells in each direction in order to match their cell size. The flow is driven by the pressure controller described in Sect. 2.3.1, with a desired flow velocity u_{ref} of 3 m s^{−1} at the hub height. Potentialtemperature stratification is turned off so that the boundary layer height coincides with the domain size in the vertical direction z. We used periodic boundary conditions in the horizontal directions, while at the upper boundary a slip condition is applied. At the ground, we used the same similarity laws as Stevens et al. (2018), with an equivalent roughness height of 0.03 mm.
The precursor is run for 100 s (corresponding to ≈ 160 flowthrough times), after which we saved the inflow field at each time step for 300 s. In the successor, we apply the precalculated inflow and source terms from the precursor, linearly interpolating from the two closest available times. At the outlet, we use a zeronormalgradient condition applied to the velocity. The remaining boundaries are treated in the same manner as the precursor. The successor domain is 8.25 m × 1.8 m × 0.675 m, with 588 × 129 × 145 cells in each direction. The first row of the wind farm is located 5 diameters from the inlet boundary, matching the setup of Stevens et al. (2018). The successor is advanced in time for 300 s, and we start gathering data after one flowthrough time (≈ 3 s).
In Fig. 4 we show the velocity profiles, averaged across the wind farm columns, for rows 1, 4, 7 and 10, together with experimental data from Chamorro and PortéAgel (2011) and numerical results from Stevens et al. (2018). As can be noticed, TOSCA matches very well with both numerical and experimental data. In the upper portion of the velocity profile, for increasing wind turbine row, both TOSCA and results from Stevens et al. (2018) predict higher velocities than the experiment. This effect is given by wind farm area blockage in the numerical domain, which causes the flow to accelerate close to the upper boundary in order to conserve mass. In the experimental data, this is not observed, as the height of the wind tunnel test section was 1.7 m.
In Fig. 5, profiles of ${u}^{\prime}/{u}_{\mathrm{ref}}$ are reported for the same location as Fig. 4. Given the velocity time history at a point, we first evaluate ${u}^{\prime}{u}^{\prime}$ by averaging the square of the fluctuation history, obtained as the difference between the velocity signal and its average. Then u^{′} is obtained as the square root of ${u}^{\prime}{u}^{\prime}$. Results show that TOSCA is well aligned with results from Stevens et al. (2018), both predicting higher fluctuations than experiments in the topmost downwind part of the wind farm for the reason mentioned above. These results demonstrate that TOSCA accurately predicts turbine–wake interactions inside a wind farm, in terms of both the mean and the fluctuations, making it suitable for the simulation of wind turbines immersed in a turbulent boundary layer.
In this section, we present CNBL results obtained using the different velocity and temperature controllers described in Sect. 2.3. In particular, we compare case S2 from Allaerts and Meyers (2017) against results obtained from TOSCA using both pressure and temperature controllers at the same time (case PT), and pressure and geostrophic controllers with no temperature forcing (case P and G, respectively). A summary of the different cases with the relative controllers is given in Table 3. The simulations employ periodic boundary conditions in the horizontal directions and a slip boundary condition at the upper boundary. Classic Monin–Obukhov similarity theory is enforced at the ground.
Following Allaerts and Meyers (2017), the domain size is 9.6 km × 4.8 km × 1.5 km in the streamwise, spanwise and vertical directions, respectively, discretized using 320 × 320 × 300 cells in each direction. Case G is forced with a geostrophic wind speed of 12 m s^{−1}, matching the setup used by Allaerts and Meyers (2017) in case S2. Conversely, in the P and PT cases the pressure controller aims to maintain a wind speed of 10.871 m s^{−1} at h_{ref}=100 m, thus matching the hubheight wind speed obtained from case G. The Coriolis parameter f_{c} is set to 10^{−4}. In all cases, potential temperature has been initialized using the Rampanelli and Zardi (2004) model, the inputs of which are reported in Table 4.
The model parameter c which determines potentialtemperature smearing across the capping inversion is set to 0.33. For the velocity, we use a uniform log law to prescribe the initial condition; namely
In case G, where geostrophic damping is not applied, care must be paid in prescribing an initial geostrophic wind consistent with geostrophic forcing in order not to trigger inertial oscillations above the capping inversion.
For this reason, u^{*} is set to $\mathit{\kappa}G/\mathrm{ln}(H/{z}_{\mathrm{0}})$ in case G, while it is calculated as $\mathit{\kappa}{u}_{\text{ref}}/\mathrm{ln}({h}_{\text{ref}}/{z}_{\mathrm{0}})$ for cases P and PT. Note that, when the pressure controller is used, inertial oscillations cannot be avoided since geostrophic wind is not known a priori.
All three simulations are carried out for 10^{5} s (≈27.8 h), while Allaerts and Meyers (2017) run case S2 for 20 h, gathering statistics from 54 000 to 72 000 s (over the last 5 h of simulated time). Geostrophic damping in cases P and PT starts at ${T}_{\mathrm{D}}=\mathrm{6}\times {\mathrm{10}}^{\mathrm{4}}$ s (this value is close to the oscillation period of the geostrophic wind), and it will be later shown that the wind angle controller of case G stabilizes the wind angle at around 6.5×10^{4} s. Hence, we average flow statistics from cases G, P and PT from 92 800 to 100 000 s, while in Allaerts and Meyers (2017) results are averaged from 54 000 to 72 000 s. In Fig. 6, we report vertical profiles of velocity magnitude, direction, shear stress and potential temperature obtained from the four different cases.
As shown before in Sect. 3.2, TOSCA predicts more mixing than SPWind, used by Allaerts and Meyers (2017). This results in a higher inversion height for a given set of ABL parameters and can be observed by comparing cases G and S2, which feature the same wind angle controller but which differ in the profile of potential temperature obtained. This leads to an increased surface temperature predicted by TOSCA and a different wind veer profile between the two codes. Although we note that such differences are accentuated by the fact that statistics from SPWind are collected at an earlier time, i.e., the CNBL has grown by a lower extent, case S2 seems to be more aligned to case PT, where the average potentialtemperature profile is kept constant by the controller. The difference in mixing between the two codes also affects the average hubheight velocity, which differs by 0.33 m s^{−1} between case G and case S2. For cases P and PT such a parameter is an input, and it has been set according to results from case G. In Table 5, output quantities extracted from the four different simulations are reported, averaging flow statistics in the abovementioned time intervals. The capping inversion center H, ground temperature θ_{0}, inversion strength Δθ and inversion width Δh are calculated by fitting the Rampanelli and Zardi (2004) model in a leastsquares sense.
Figure 6, together with quantitative data reported in Table 5, demonstrates how the pressure controller with geostrophic damping (case P) almost exactly matches results obtained using the geostrophic controller (case G), predicting a geostrophic wind that only differs by 0.25 % with respect to G=12 m s^{−1}. Figure 6 also highlights how sensitive the ABL is to its heating history, since case PT – where the average θ profile is kept constant – predicts a lower geostrophic wind than cases G and P. In fact, it can be noticed from Fig. 7b how the inversion height is kept constant in case PT, while it grows in time during simulations P and G. For this reason, in the latter cases the ABL will experience a slightly higher amount of dissipation, which results in a small increase in the geostrophic wind if compared to case PT. Therefore, in simulations P, G and S2, the boundary layer is developing against a potentialtemperature profile that is slowly evolving, in turn affecting the mean velocity profile. This mechanism, which is of course physical, does not reproduce what happens in real life, where the boundary layer stability evolves following the timescale of the diurnal cycle instead. As a consequence, since such temperature drift is physical but arises from an idealization, we believe that fixing the average potentialtemperature profile would represent a more consistent method for conducting these types of idealized simulations. In addition, we also suggest driving the ABL with a pressure controller, which allows specifying the hubheight velocity, as the issue related to geostrophic inertial oscillations can be addressed using the proposed geostrophic damping method. This would allow precursor simulations to reach a truly statistically steady state, leading to better agreement between different codes applied to inflow conditions used for successor simulations.
Figure 7a shows the evolution of the geostrophic wind (components and magnitude), calculated as the spatial average in the homogeneous directions and at those cells where $h>H+\mathrm{\Delta}$, produced by cases P and G. It can be seen how the developed damping technique is able to stop inertial oscillations after a time T_{D}+T_{3 %}, reaching a geostrophic wind that only differs by 0.25 % from the simulation where the geostrophic controller has been applied (T_{3 %}=17 500; see Sect. 2.3.1 for definition). Moreover, in Fig. 7b we report wind angle and velocity magnitude horizontally averaged at the reference height, together with the height of the inversion center over time, evaluated by fitting the Rampanelli and Zardi (2004) model at each time step. It is evident that the pressure controller precisely maintains the wind at the desired speed and in the desired direction. Interestingly, it can be also noticed that the geostrophic controller produces small oscillations in the hubheight wind speed. These are inertial oscillations as well, but they are naturally damped by turbulence as they happen inside the boundary layer. Finally, looking at the evolution of the inversion layer height in Fig. 7b and at the final potentialtemperature profile in Fig. 6, it is clear that controlling the mean potential temperature prevents the boundary layer from growing indefinitely, preserving the initial capping inversion height and the initial value of potential temperature at the ground.
In this section, we present results from the simulation of a finitesize wind farm consisting of 100 NREL 5 MW wind turbines, aligned in 20 rows and 5 columns, with streamwise and spanwise spacing of 5 and 4.76 rotor diameters, respectively. We include thermal stratification to assess the effects of gravity wave blockage for a lapse rate of 1 K km^{−1}, a capping inversion centered at 500 m with a strength of 7.312 K. Given the large scale of the gravity waves, the numerical domain is set to 40 × 21 × 28 km in the streamwise, spanwise and vertical directions, respectively. All directions are graded to reach a mesh resolution of 30 × 12.5 × 10 m around the wind turbines. The hybrid offline–concurrentprecursor technique described in Sect. 2.4 has been used to spin up turbulence in the precursor, providing a timeresolved CNBL inflow for the successor simulation. This technique is combined with a Rayleigh damping layer and the advection damping technique (Lanzilao and Meyers, 2022a) to ensure low reflectivity of gravity waves from the top boundary and the fringe region exit. Further details on the successor–precursor meshes and simulations; CNBL parameters; and tuning of the fringe, Rayleigh and advection damping region coefficients are given in Appendix C.
Figure 8 shows hubheight instantaneous velocity and pressure contours around the wind farm. The gravity wave footprint inside the ABL can be clearly noticed in the pressure field, together with the smallscale pressure increase in front of each rotor. This effect is superimposed on the much larger pressure variation due to atmospheric gravity waves, which take place from the farm entrance to the exit.
Regarding the instantaneous velocity field, streamwise streaks generated by elongated turbulence structures can be appreciated. The large size of these structures is related to the high value of the prescribed equivalent roughness height z_{0} and to the fact that periodic boundary conditions artificially increase their length when they span the entire domain in the streamwise direction. If averages are gathered for a sufficient amount of time, these streaks are not expected to alter the simulation results from the wind farm power production standpoint.
At any given location, we define the perturbation value of a quantity as the difference between its successor time average and the precursor time average, evaluated at the same height. Figure 9 shows horizontal contours of pressure and temperature perturbations at the hub and inversion heights, respectively. An interesting aspect is that, due to the presence of the Coriolis force, the direction of propagation of interfacial waves in the inversion layer is not aligned with the wind farm streamwise symmetry axis.
For instance, the two trains of waves generated by the positive and negative inversion layer displacements, at the wind farm entrance and exit, respectively, have a spanwise offset, resulting in a much more complex interaction. Moreover, spurious wave interactions with their periodic images can also be noticed, but the spanwise size of the domain ensures that they happen far from and downstream of the wind farm. Nevertheless, we are developing a lateral fringe region, which is aimed at removing this effect, where the instantaneous desired flow is reconstructed from the concurrent precursor, allowing for a smaller spanwise domain size. At this time, we do not believe wave interactions from periodic images alter the gravity wave pattern in the region near and upstream of the wind farm, which is of primary interest in order to assess wind farm blockage.
Figure 10a shows the average hubheight velocity field, from which the effect of gravity waves on velocity can be assessed. In particular, positive perturbations are observed where negative pressure gradients are experienced and vice versa. The perturbation, averaged along the spanwise direction and within the region enclosed by the wind farm spanwise limits, is shown quantitatively in Fig. 10b. In particular, we record velocity reductions as high as 5.80 %, 2.15 % and 1.02 % at 2.5, 20 and 40 diameters upstream of the first turbine row. Another aspect that is noticeable from the average fields of Figs. 9 and 10 is the presence of highfrequency oscillations in the spanwise direction. These are caused by the turbulence streaks mentioned earlier in the present section, and, while we plan on eliminating this shortcoming by applying a spanwise offset in the streamwise periodicity to artificially break the locking in position of such structures (Munters et al., 2016), we note that their manifestation could be readily eliminated by averaging for a time period longer than 15 000 s.
Overall, these results indicate that the developed framework and methodology allow for conducting finitesize wind farm simulations, capturing gravity wave effects unaltered by spurious wave reflections from the fringe region or interaction from periodic images.
In the present paper, we introduced TOSCA, a new opensource LES framework for the simulation of large wind farms interacting with thermally stratified boundary layers. We validated TOSCA's wind turbine models and its ability to simulate the evolution of conventionally neutral boundary layers and to accurately predict the flow around infinite wind farms in neutral conditions. We presented a new controlling methodology for ABL precursors that allows us to prescribe a desired wind speed at a reference height – located inside the boundary layer – while at the same time avoiding velocity oscillations produced by the Coriolis force in the geostrophic region above the inversion height. This approach, if combined with a potentialtemperature controller, allows different codes to obtain CNBL inflow profiles that are characterized by the same potentialtemperature profile and hubheight wind speed. Conversely, using geostrophic forcing makes the hubheight velocity dependent on the amount of numerical dissipation specific to the adopted code, while the final temperature profile depends on both numerical dissipation and precursor simulated time. We also described a new methodology for simulating finitesize wind farms under atmospheric gravity wave effects. In particular, we introduced the hybrid offline precursor–concurrentprecursor method, where the offline technique is used on a small domain, in order to spin up ABL turbulence, while the concurrent method is adopted for the turbine simulation. In fact, we found that the concurrent precursor, combined with a fringe region, is a crucial element to simultaneously avoid spurious gravity wave reflections and provide a timeresolved turbulent inflow. The offline precursor data are used to start up the flow field in the concurrent precursor by means of spanwise periodization. The concurrentprecursor domain is usually bigger than required, as its size is determined by the successor domain that runs concurrently. Hence, being able to reach steadystate turbulent statistics on a smaller domain is indeed convenient, as it makes finite wind farm simulations less computationally intensive.
Finally, we demonstrated that TOSCA is able to simulate wind farm gravity wave interactions and largescale blockage effects. Specifically, for the CNBL simulated herein, we measured a velocity reduction of 5.80 % at 2.5 diameters upstream of the first row.
In the future, we will implement shifted periodic boundary conditions to obtain field statistics which are less dependent on the spanwise location, and we will address the heat flux mismatch above the inversion layer.
The adoption of generalized curvilinear coordinates allows for the computational mesh to follow terrain coordinates if required or to be stretched and deformed with the only condition that the indexing remains structured. We denote a set of generalized curvilinear coordinates as l_{i}, with $i=\mathrm{1},\mathrm{2},\mathrm{3}$, by which points in a threedimensional Euclidean space E^{3} may be defined. Cartesian coordinates are a special case of such a generalization and will be denoted as x_{i}, with $i=\mathrm{1},\mathrm{2},\mathrm{3}$. When using explicit notation, the three curvilinear directions will be identified by Greek symbols as ξ, η and ζ. With these definitions, and given the position vector r of a point P in Cartesian space, the covariant base vectors can be expressed as ${\mathit{g}}_{i}=\partial \mathit{r}/\partial {l}_{i}$ (with Cartesian components $({\mathit{g}}_{i}{)}_{j}=\partial {x}_{j}/\partial {l}_{i}$), while contravariant base vectors are given by g^{i}=∇l_{i} (with Cartesian components $({\mathit{g}}^{i}{)}_{j}=\partial {l}_{i}/\partial {x}_{j}$). As a result, the following relation holds between covariant and contravariant base vectors:
where J is the Jacobian of the transformation defining l_{i} in terms of x_{j}, i.e., the determinant of the matrix of partial derivatives $\partial {l}_{i}/\partial {x}_{j}$. It is required that J≠0, which is equivalent to asking that covariant base vectors are not coplanar. Note that they are usually neither unit vectors nor orthogonal to each other. Given a set of curvilinear coordinates l_{i}, with covariant base vectors g_{i} and contravariant base vectors g^{i}, it is possible to define the covariant and contravariant metric tensors through the scalar products:
where the repeated index implies summation. Metric tensors satisfy $J=\sqrt{det\left({g}_{ij}\right)}$ and ${J}^{\mathrm{1}}=\sqrt{det\left({g}^{ij}\right)}$.
The use of generalized curvilinear coordinates allows differential operators on any structured mesh to be expressed using a Cartesianlike discretization along the curvilinear directions, which are chosen to be the local structured grid lines. Moreover, the quantities
are equal to face area vectors if evaluated at cell faces. Contravariant base vectors components $\partial {l}_{i}/\partial {x}_{j}$ in Eq. (A3) are evaluated using Eq. (A1), i.e., from the covariant base vectors, which are easily obtained exploiting finite differences.
Contravariant components of any vector field u, a function of position r, can be expressed in terms of its Cartesian components as ${u}^{i}=\mathit{u}\cdot {\mathit{g}}^{i}$. If one instead uses ${\mathit{g}}^{i}/J$ (namely the face area vector along the ith direction S^{i}) and the relation is again evaluated at cell faces, contravariant fluxes are obtained as ${V}^{i}=\mathit{u}\cdot {\mathit{S}}^{i}$. In TOSCA, only the independent variables (positions) are transformed in curvilinear coordinates using the chain rule and integration by parts, while dependent variables are retained in Cartesian coordinates. This partial transformation avoids computing the Christoffel symbols of the second kind, which are cumbersome to evaluate numerically. Moreover, they would increase the requirements of smoothness of the computational mesh, as they involve secondorder derivatives of the transformation metrics. The momentum equation is finally dotted with the face area vectors, so it can be partially written in terms of contravariant fluxes as
Equation (A5) is used to solve for contravariant fluxes, which are staggered at cell faces, while pressure is located at cell centers. In contrast to a staggered formulation using a full transformation, where Cartesian velocity does not appear in the equations, in a partial transformation all Cartesian velocity components are required at each face center in order to discretize Eq. (A5) with the same accuracy. One alternative would be to solve all components of the momentum equation at each face center in order for the Cartesian velocity to be attainable without interpolation. Although this approach has been adopted in the literature (Maliska and Raithby, 1984), it triples the computational cost. In TOSCA, we follow the approach of Ge and Sotiropoulos (2007), where the momentum equation is first discretized at cell centers and then interpolated and solved at face centers in a staggered fashion. Cartesian velocity is subsequently reconstructed at cell centers by interpolating contravariant fluxes at the same location. With respect to a standard staggered formulation (e.g., in Cartesian coordinates), this procedure encompasses additional steps for interpolating the discretized momentum equation at face centers and for transforming the interpolated fluxes into the Cartesian velocity at cell centers (flux interpolation is required in either case). It should be noted that the overhead in computational cost is minimal because it only involves 1D interpolations along grid lines, for which a secondorder central scheme is used. Another slightly different approach (Rosenfeld et al., 1992) is to discretize the momentum equation in a staggered manner. This avoids interpolating the whole momentum righthand side at cell faces, but it requires interpolating contravariant fluxes instead. In either case, methods based on partially transformed equations involve an additional interpolation step (as contravariant fluxes and Cartesian velocity are defined at different locations). This imposes slightly tighter constraints on the time step value in order to keep the method stable. For this reason, we opted for an implicit treatment of advection and viscous terms in Eq. (A5). Specifically, we use the matrixfree Newton–Krylov solver implemented in PETSc (Balay et al., 2022), where the iterative Krylov subspace generalized minimum residual (GMRES) method (Saad and Schultz, 1986) is used to solve the linear system associated with each inner iteration. (See Knoll and Keyes, 2004, for a comprehensive review and application of matrixfree Newton–Krylov methods.) In addition, such a hybrid staggered–nonstaggered formulation facilitates the application of boundary conditions, which are prescribed for the Cartesian velocity using ghost cells.
Pressure–velocity coupling is provided using a secondorder fractional step method similar to that of van Kan (1986), where velocity is first guessed by solving for the contravariant fluxes, which are then projected into a divergencefree space by means of a pressure correction ϕ obtained by solving a Poisson equation. Potential temperature is subsequently solved using the new velocity field, with an implicit treatment of the righthand side. Time discretization uses a secondorder implicit scheme for both momentum and temperature equations. All derivatives are discretized using the secondorder central scheme, while the advection term in Eq. (A5) is discretized using a blend between central and QUICK (Leonard, 1979) schemes. The blending is such that QUICK is used in regions of almost uniform or slowly varying velocity, avoiding the oscillations produced by the central scheme in such regions.
To model the subgrid stresses, TOSCA uses the dynamic Smagorinsky model (Lilly, 1992; Germano et al., 1991), with Lagrangian averaging of the model coefficient C_{s} (Meneveau et al., 1996). The model has been recast in generalized curvilinear coordinates, similar to what was presented in Armenio and Piomelli (2000). The effect of unresolved scales in the momentum equation, after the filtering operation, appears in Cartesian coordinates through the equation
where the overbar indicates the filtering operation applied e.g. to a variable q as
where the integration is performed over the entire flow domain 𝒟 (Moeng, 1984). In curvilinear coordinates, Eq. (B1) reads
where
and assuming a linear eddy viscosity model,
where ${S}_{j}^{k}=\mathrm{1}/J\partial {l}_{k}/\partial {x}_{j}$ comprises the face area vectors, ${S}_{ij}=\frac{\mathrm{1}}{\mathrm{2}}(\partial {u}_{i}/\partial {x}_{j}+\partial {u}_{j}/\partial {x}_{i})$ is the symmetric part of the velocity gradient tensor and Δ is the cubic root of the local cell volume. Using the idea of Germano et al. (1991), a second filter, denoted as $\stackrel{\mathrm{\u0303}}{\cdot}$, can be applied which has $\stackrel{\mathrm{\u0303}}{\mathrm{\Delta}}=\mathrm{3}\mathrm{\Delta}$ in TOSCA, leading to the tensor
which accounts for the effect of the unresolved plus the smallest resolved scales. The Germano tensor, i.e., the contribution to the resolved stresses from the largest unresolved motions, is defined in generalized curvilinear coordinates by subtracting the tildefiltered Eq. (B3) from Eq. (B7):
Using Eqs. (B5) and (B6) to express $\stackrel{\mathrm{\u0303}}{{\mathit{\sigma}}_{j}^{k}}$ and ${T}_{j}^{k}$ gives
where in Eq. (B10) the approximation $\stackrel{\mathrm{\u0303}}{\stackrel{\mathrm{\u203e}}{{S}_{j}^{k}{S}_{ij}}}=\stackrel{\mathrm{\u0303}}{\stackrel{\mathrm{\u203e}}{{S}_{j}^{k}}}{\stackrel{\mathrm{\u0303}}{\stackrel{\mathrm{\u203e}}{S}}}_{ij}$ has been used. In fact, as the LES filter has a size of three mesh cells in each direction, $\stackrel{\mathrm{\u203e}}{{S}_{j}^{k}}$ (the face area vectors at the central cell) and $\stackrel{\mathrm{\u0303}}{\stackrel{\mathrm{\u203e}}{{S}_{j}^{k}}}$ (the filtered face area vectors within the box) are almost identical provided that mesh grading is smooth enough. Conversely, the equality holds exactly if a homogeneous filter and a uniform mesh are considered. Inserting Eqs. (B9) and (B10) into Eq. (B8) leads to
where
It is now possible to find C_{s} in a leastsquares sense as
Note that the above relation is not invariant with respect to rotation of the reference frame because it implicitly contains the face area vectors; hence tensors are no longer symmetric. Variables must then be transformed into physical space to find C_{s} as
where g_{ij} is the covariant metric tensor. Since the C_{s} coefficient oscillates in space, some sort of average is required. TOSCA follows the approach presented in Meneveau et al. (1996), where the numerator and denominator of Eq. (B13) are averaged along streamlines as
where $W\left(t\right)=\mathrm{1}/{T}_{\mathrm{s}}\mathrm{exp}(t/T)$ is a weighting function and T_{s} is a timescale defined as
The integrals of Eqs. (B14) and (B15) can be evaluated as
where $\mathit{\u03f5}=(\mathrm{\Delta}t/{T}_{\mathrm{s}})/(\mathrm{1}+\mathrm{\Delta}t/{T}_{\mathrm{s}})$. We use trilinear interpolation formulas to evaluate the integrals I_{GM} and I_{MM} at the x−uΔt position, and all quantities are evaluated at cell centers, including contravariant fluxes, which are linearly interpolated from the faces.
Regarding the potentialtemperature equation, subgrid fluxes are evaluated following the approach of Moeng (1984), i.e., through the definition of a thermal eddy diffusivity ${\mathit{\kappa}}_{t}={\mathit{\nu}}_{t}/P{r}_{t}$, where Pr_{t} is the turbulent Prandtl number, which depends on stability as
Note that, if the potentialtemperature gradient is locally stable, Pr_{t}→1, while for neutral or unstable cases, $P{r}_{t}=\mathrm{1}/\mathrm{3}$. This reflects the decrease in the mixing length scale under stable conditions (Schumann, 1991).
In this section, we describe in detail the setup of the finite wind farm case presented in Sect. 5. To avoid wave reflections from inflow–outflow boundaries, we adopt periodic boundary conditions and the concurrentprecursor technique. This also provides a suitable turbulent inflow, eliminating the wind farm wake readvected at the inlet by the periodic boundaries. To avoid wave reflections from the upper boundary, we use a Rayleigh damping layer, while lateral boundaries are periodic. Spanwise periodicity implies that gravity waves induced by the wind farm will interact with their periodic images, requiring the domain to be sufficiently large for these interactions to happen far from and downstream of the wind turbines. Moreover, we use the advection damping technique developed by Lanzilao and Meyers (2022a) to ensure that interactions between fringegenerated and physical waves are not advected downstream but remain trapped inside the advection damping region.
The size of the successor domain is 40 km × 21 km × 28 km in the streamwise, spanwise and vertical directions, respectively, discretized with 1554 × 1194 × 345 cells. All directions are graded to reach a mesh resolution of 30 m × 12.5 m × 10 m around the wind farm.
The wind farm has a rectangular planform, with 20 rows and 5 columns. The first row is located at x=0, and extends from 300 to 2700 m. This determines a lateral spacing of 600 m (4.76 D), while streamwise spacing is set to 630 m (5 D). Wind turbines are equipped with angular velocity and pitch controllers described in Jonkman et al. (2009). A very simple yaw controller is also added, which rotates each wind turbine independently with a uniform speed of 0.5^{∘} s^{−1} when flow misalignment exceeds 1^{∘}. The flow angle is calculated by filtering the wind velocity at a sampling point located 1 D upstream of the rotor center, using a time constant of 600 s. Turbines are modeled using the ADM, while tower and nacelle are not accounted for. The projection width ϵ is set to 18.75 m.
The concurrentprecursor mesh coincides with the portion of the successor domain which is located inside the fringe region. In particular, it is 5 km × 21 km × 28 km. The mesh resolution in the streamwise direction is 15 m, while in the spanwise and vertical directions, it is the same as the successor.
In order to save computational resources, we do not run the whole precursor simulation on the concurrentprecursor mesh, whose size is determined by the wind farm and gravity wave parameters. Instead, we perform the spinup phase on a 6 km × 3 km × 1 km domain, characterized by a resolution of 15 m × 15 m in the streamwise and spanwise directions. The vertical direction is discretized in the same manner as the successor in order to increase the resolution inside the capping inversion layer. This spinup phase is carried out for 10^{5} s, after which an inflow database is collected. The generated inflow database is then used to start up the solution in the concurrentprecursor and successor domains. This technique, which we refer to as the hybrid offline–concurrent precursor, is explained in Sect. 2.4. In the successor, after a spinup of 5000 s, corresponding to slightly more than one flowthrough time, data are gathered for 15 000 s.
The offline precursor simulation uses the pressure and temperature controllers described in Sect. 2.3, while in both the concurrentprecursor and the successor simulations, velocity is controlled using a constant source term, obtained by averaging the offline precursor source from 100 000 to 120 000 s. The temperature controller is retained in the concurrentprecursor simulation, but it is switched off in the successor so that the inversion height is free to be perturbed by the wind farm.
CNBL parameters used for the offline precursor are summarized in Table C2. They are calculated based on the sensitivity analysis performed in Allaerts and Meyers (2019). In particular, our objective is to choose a set of nondimensional parameters such that the capping inversion layer is strongly perturbed by the wind farm. This results in a capping inversion Froude number of Fr=0.94 and an internal wave parameter of P_{N}=3.02. These nondimensional groups are related to the physics and magnitude of interfacial waves inside the inversion layer and internal gravity waves above the ABL, respectively. The temperature profile is initialized according to Rampanelli and Zardi (2004).
In Fig. C1 we show vertical profiles of the wind speed magnitude, inflow angle, nondimensional shear stress and potential temperature, averaged over the last 15 000 s from the concurrentprecursor domain. It can be noticed how the pressure controller accurately maintains the desired wind speed and direction at h_{ref} and how the temperature controller removes the ground temperature shift observed in the previous sections by keeping the average profile constant in time. As a consequence, inversion height and strength are maintained equal to their initial values (Table C2), while the resulting friction velocity corresponds to 0.432 m s^{−1}.
Regarding the successor case, we followed the approach of Lanzilao and Meyers (2022a) to choose the damping layer and fringe coefficients. In particular, after a reflectivity study that employed a computationally cheap canopy model (not shown here), we found that a Rayleigh damping coefficient of ν_{RDL}=0.05 and a fringe damping coefficient of ν_{FR}=0.03 yielded minimal gravity wave reflectivity. The dominant vertical wavelength of the gravity waves is estimated as ${\mathit{\lambda}}_{z}=\mathrm{2}\mathit{\pi}G/N\approx \mathrm{11.8}$ km (Allaerts and Meyers, 2017), where N is the Brunt–Väisälä frequency ($N=\mathrm{5.72}\times {\mathrm{10}}^{\mathrm{3}}$ s^{−1} based on parameters listed in Table C2) and G=10.815 m s^{−1}. We ensure that at least one λ_{z} is contained in the Rayleigh damping layer by setting its width to 12 km. Regarding the advection damping technique developed by Lanzilao and Meyers (2022a), we observed that their guidelines in how to choose the length of the advection damping region did not apply to our case, which is characterized by a very strong inversion layer and a shallow boundary layer. In fact, we believe that a key parameter that needs to be tuned in order to avoid spurious gravity wave interactions is the length of the region where advection damping is applied after the fringe. This holds in particular for subcritical (Fr < 1) cases, where waves inside the capping inversion can propagate against the flow. Here, perturbations would be propagated upstream of the wind farm to the fringe exit, being suddenly forced to obey the precursor inflow inside the fringe region. Such a sharp change in the boundary layer displacement at the fringe exit induces spurious gravity waves which remain trapped at their streamwise location if horizontal advection of vertical velocity is turned off. Nevertheless, these waves would interact with physical waves from inside the domain, resulting in more spurious interactions. As a consequence, it is crucial to ensure that all spurious interactions generated by this mechanism are fully contained within the advection damping region and are not advected downstream.
We used the same damping functions as Lanzilao and Meyers (2022a), and in Table C3 their parameters are reported for our finite wind farm simulation.
In this section, we show TOSCA's strong and weak parallel performance by running CNBL simulations with an increasing number of Niagara nodes for three different mesh sizes. The simulation setup corresponds to the offline precursor described in Appendix C. The different meshes are evaluated by systematically doubling the number of elements in each direction, starting from 300 × 300 × 100 cells. As a consequence, they consist of 9 million, 72 million and 576 million elements in total. Table D1 reports the number of nodes for each run, which only consisted of 2 h of wallclock time.
Tests have been performed on Compute Canada's Niagara (2018) cluster, which consists of 2024 nodes, each with 40 Intel “Skylake” cores at 2.4 GHz or 40 Intel “Cascade Lake” cores at 2.5 GHz. Node interconnection consists of an EDR InfiniBand network, organized in a “Dragonfly+” topology with five dragonfly wings. Figure D1 shows the time per iteration as the node count increases for each of the CNBL meshes. TOSCA's strong scaling performance remains close to linear until roughly 25 000 cells per core are reached, which we identify as a reasonable tradeoff between efficiency and speed. TOSCA was also successfully run at scale on the entire Niagara cluster to simulate a finite wind farm on a mesh exceeding 1 billion elements, proving TOSCA's suitability for massively parallel computations.
We also highlight that the time per iteration does not reflect the actual speed at which the simulation advances in time, as the time step size depends on the numerical method. In implicit methods like the Newton–Krylov solvers employed by TOSCA, the computational cost of each time step depends on the time step size, whereas these quantities are unrelated in explicit methods. Nevertheless, implicit methods are able to advance in time with a Courant–Friedrichs–Lewy (CFL) number greater than 1 (we used 0.9 for these analyses), while explicit methods are usually limited to a value close to 0.5.
Regarding actuator models, their solution and I/O operations are also parallelized in TOSCA. Specifically, we define a sphere of cells, for each wind turbine, that contains all cells that the rotor can possibly intersect when yawing. Processors owning mesh cells belonging to the sphere are then grouped into turbinespecific subcommunicators, which are used to solve wind turbines simultaneously. Hence, provided that a sufficiently high core count is used, the wind turbine update time in TOSCA is independent of the number of wind turbines in the simulation, and each communicator can write turbine data to file simultaneously.
For the finite wind farm simulation presented in Sect. 5 (100 wind turbines), the turbine update time was less than 0.1 s. The individual turbine update time depends on the cell size and on the actuator pointprocessor ownership search (which processor in the communicator controls which actuator point). The latter is only triggered by a change in yaw for the ADM and uniform ADM models, while it has to be performed at every iteration for the ALM, as actuator points are physically rotating.
TOSCA is available at https://doi.org/10.17605/OSF.IO/Q4VAF (Stipa et al., 2023).
The dataset can be made available from the authors upon request.
Conceptualization: SS, DA, JB; methodology: SS; software: SS, AA; validation: SS; formal analysis: SS; investigation: SS; computational resources: JB; data curation: SS; writing – original draft preparation: SS; writing – review and editing: JB, DA; visualization: SS; supervision: JB, DA; project administration: JB; funding acquisition: JB. All authors have read and agreed to the published version of the paper.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
Computational resources provided by the Digital Research Alliance of Canada and Advanced Research Computing at the University of British Columbia are gratefully acknowledged.
This research has been supported by the Natural Sciences and Engineering Research Council of Canada (grant no. 556326).
This paper was edited by Sukanta Basu and reviewed by two anonymous referees.
Abkar, M. and PortéAgel, F.: The effect of freeatmosphere stratification on boundarylayer flow and power output from very large wind farms, Energies, 6, 2338–2361, https://doi.org/10.3390/en6052338, 2013. a, b, c, d
Ahsbahs, T., Nygaard, N. G., Newcombe, A., and Badger, M.: Wind Farm Wakes from SAR and Doppler Radar, Remote Sensing, 12, 462, https://doi.org/10.3390/rs12030462, 2020. a
Ainslie, J.: Calculating the flowfield in the wake of wind turbines, J. Wind Eng. Ind. Aerod., 27, 213–224, https://doi.org/10.1016/01676105(88)900372, 1988. a
Allaerts, D.: PhD Thesis: Largeeddy Simulation of Wind Farms in Conventionally Neutral and Stable Atmospheric Boundary Layers, Tech. rep., KU Leuven, https://lirias.kuleuven.be/1731391&lang=en (last access: 1 November 2022), 2016. a, b, c
Allaerts, D. and Meyers, J.: Large eddy simulation of a large windturbine array in a conventionally neutral atmospheric boundary layer, Phys. Fluids, 27, 065108, https://doi.org/10.1063/1.4922339, 2015. a, b
Allaerts, D. and Meyers, J.: Boundarylayer development and gravity waves in conventionally neutral wind farms, J. Fluid Mech., 814, 95–130, https://doi.org/10.1017/jfm.2017.11, 2017. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o
Allaerts, D. and Meyers, J.: Gravity Waves and WindFarm Efficiency in Neutral and Stable Conditions, Bound.Lay. Meteorol., 166, 269–299, https://doi.org/10.1007/s1054601703075, 2018. a, b
Allaerts, D. and Meyers, J.: Sensitivity and feedback of windfarminduced gravity waves, J. Fluid Mech., 862, 990–1028, https://doi.org/10.1017/jfm.2018.969, 2019. a, b
Allaerts, D., Quon, E., Draxl, C., and Churchfield, M.: Development of a Time–Height Profile Assimilation Technique for LargeEddy Simulation, Bound.Lay. Meteorol., 176, 329–348, https://doi.org/10.1007/s10546020005385, 2020. a
Allaerts, D., Quon, E., and Churchfield, M.: Using observational meanflow data to drive largeeddy simulations of a diurnal cycle at the SWiFT site, Wind Energy, 26, 469–492, https://doi.org/10.1002/we.2811, 2023. a
Armenio, V. and Piomelli, U.: A Lagrangian Mixed SubgridScale Model in Generalized Coordinates, Flow Turbul. Combust., 65, 51–81, https://doi.org/10.1023/A:1009998919233, 2000. a
Balay, S., Abhyankar, S., Adams, M. F., Benson, S., Brown, J., Brune, P., Buschelman, K., Constantinescu, E., Dalcin, L., Dener, A., Eijkhout, V., Faibussowitsch, J., Gropp, W. D., Hapla, V., Isaac, T., Jolivet, P., Karpeev, D., Kaushik, D., Knepley, M. G., Kong, F., Kruger, S., May, D. A., McInnes, L. C., Mills, R. T., Mitchell, L., Munson, T., Roman, J. E., Rupp, K., Sanan, P., Sarich, J., Smith, B. F., Zampini, S., Zhang, H., Zhang, H., and Zhang, J.: PETSc/TAO Users Manual, Tech. Rep. ANL21/39 – Revision 3.18, Argonne National Laboratory, https://doi.org/10.2172/1893326, 2022. a, b
Bastankhah, M. and PortéAgel, F.: A new analytical model for windturbine wakes, Renew. Energ., 70, 116–123, https://doi.org/10.1016/j.renene.2014.01.002, 2014. a
Bleeg, J., Purcell, M., Ruisi, R., and Traiger, E.: Wind Farm Blockage and the Consequences of Neglecting Its Impact on Energy Production, Energies, 11, 1609, https://doi.org/10.3390/en11061609, 2018. a
Branlard, E. and Gaunaa, M.: Cylindrical vortex wake model: right cylinder, Wind Energy, 524, 1973–1987, 2014. a
Branlard, E. and Meyer Forsting, A. R.: Assessing the blockage effect of wind turbines and wind farms using an analytical vortex model, Wind Energy, 23, 2068–2086, https://doi.org/10.1002/we.2546, 2020. a
Branlard, E., Quon, E., Forsting, A. R. M., King, J., and Moriarty, P.: Wind farm blockage effects: comparison of different engineering models, J. Phys. Conf. Ser., 1618, 062036, https://doi.org/10.1088/17426596/1618/6/062036, 2020. a
Breton, S., Sumner, J., Sørensen, J., Hansen, K., Sarmast, S., and Ivanell, S.: A survey of modelling methods for highfidelity wind farm simulations using large eddy simulation, Philos. T. R. Soc. A, 375, 20160097, https://doi.org/10.1098/rsta.2016.0097, 2017. a
Calaf, M., Meneveau, C., and Meyers, J.: Large eddy simulations of fully developed windturbine array boundary layers, Phys. Fluids, 22, 015110, https://doi.org/10.1063/1.3291077, 2010. a
Centurelli, G., Vollmer, L., Schmidt, J., Dörenkämper, M., Schröder, M., Lukassen, L. J., and Peinke, J.: Evaluating Global Blockage engineering parametrizations with LES, J. Phys. Conf. Ser., 1934, 012021, https://doi.org/10.1088/17426596/1934/1/012021, 2021. a
Chamorro, L. and PortéAgel, F.: Turbulent Flow Inside and Above a Wind Farm: A WindTunnel Study, Energies, 4, 1916–1936, https://doi.org/10.3390/en4111916, 2011. a, b, c, d
Churchfield, M., Lee, S., Michalakes, J., and Moriarty, P.: Numerical Study of the Effects of Atmospheric and Wake Turbulence on Wind Turbine Dynamics, J. Turbul., 13, N14, https://doi.org/10.1080/14685248.2012.668191, 2012a. a, b
Churchfield, M., Lee, S., Moriarty, P., Martínez Tossas, L., Leonardi, S., Vijayakumar, G., and Brasseur, J.: A LargeEddy Simulation of WindPlant Aerodynamics, https://doi.org/10.2514/6.2012537, ISBN 9781600869365, 2012b. a, b
Churchfield, M., Schreck, S., Martínez Tossas, L., Meneveau, C., and Spalart, P.: An Advanced Actuator Line Method for Wind Energy Applications and Beyond, in: AIAA SciTech 35th Wind Energy Symposium, 9–13 January 2017, Grapevine, Texas, https://doi.org/10.2514/6.20171998, 2017. a
Etling, D.: Modelling the vertical ABL structure, in: collection: Modelling of Atmospheric Flow Fields, edited by: Lalas, D. P. and Ratto, C. F., World Scientific, 45–86, https://doi.org/10.1142/9789814447164_0003, 1996. a
Falgout, R. D. and Yang, U. M.: hypre: A Library of High Performance Preconditioners, in: Computational Science – ICCS 2002, edited by Sloot, P. M. A., Hoekstra, A. G., Tan, C. J. K., and Dongarra, J. J., Springer Berlin Heidelberg, Berlin, Heidelberg, 632–641, ISBN 9783540477891, 2002. a
Fleming, P. A., Gebraad, P. M., Lee, S., van Wingerden, J.W., Johnson, K., Churchfield, M., Michalakes, J., Spalart, P., and Moriarty, P.: Evaluating techniques for redirecting turbine wakes using SOWFA, Renew. Energ., 70, 211–218, https://doi.org/10.1016/j.renene.2014.02.015, 2014. a
Gabriel, E., Fagg, G. E., Bosilca, G., Angskun, T., Dongarra, J. J., Squyres, J. M., Sahay, V., Kambadur, P., Barrett, B., Lumsdaine, A., Castain, R. H., Daniel, D. J., Graham, R. L., and Woodall, T. S.: Open MPI: Goals, Concept, and Design of a Next Generation MPI Implementation, in: Proceedings, 11th European PVM/MPI Users' Group Meeting, Budapest, Hungary, 97–104, https://doi.org/10.1007/9783540302186_19, 2004. a
Ge, L. and Sotiropoulos, F.: A numerical method for solving the 3D unsteady incompressible Navier–Stokes equations in curvilinear domains with complex immersed boundaries, J. Comput. Phys., 225, 1782–1809, https://doi.org/10.1016/j.jcp.2007.02.017, 2007. a
Germano, M., Piomelli, U., Moin, P., and Cabot, W. H.: A dynamic subgrid‐scale eddy viscosity model, Phys. Fluids AFluid, 3, 1760–1765, https://doi.org/10.1063/1.857955, 1991. a, b
Glauert, H.: Airplane Propellers, Springer Berlin Heidelberg, Berlin, Heidelberg, 169–360, https://doi.org/10.1007/9783642914874_3, ISBN 9783642914874, 1935. a
Gribben, B. J. and Hawkes, G. S.: A potential flow model for wind turbine induction and wind farm blockage, Systems and Engineering Technology, ]https://www.fnc.co.uk/windfarmblockage/ (last access: 1 April 2023), 2019. a
Haji Mohammadi, M., Sotiropoulos, F., and Brinkerhoff, J.: Moving least squares reconstruction for sharp interface immersed boundary methods, International Journal for Numerical Methods in Fluids, 90, 57–80, https://doi.org/10.1002/fld.4711, 2019. a
Haupt, S. E., Kosović, B., Berg, L. K., Kaul, C. M., Churchfield, M., Mirocha, J., Allaerts, D., Brummet, T., Davis, S., DeCastro, A., Dettling, S., Draxl, C., Gagne, D. J., Hawbecker, P., Jha, P., Juliano, T., Lassman, W., Quon, E., Rai, R. K., Robinson, M., Shaw, W., and Thedin, R.: Lessons learned in coupling atmospheric models across scales for onshore and offshore wind energy, Wind Energ. Sci., 8, 1251–1275, https://doi.org/10.5194/wes812512023, 2023. a
Inoue, M., Matheou, G., and Teixeira, J.: LES of a Spatially Developing Atmospheric Boundary Layer: Application of a Fringe Method for the Stratocumulus to Shallow Cumulus Cloud Transition, Mon. Weather Rev., 142, 3418–3424, https://doi.org/10.1175/MWRD1300400.1, 2014. a, b
Jensen, N.: A note on wind generator interaction, Tech. Rep. Risø M2411, https://orbit.dtu.dk/en/publications/anoteonwindgeneratorinteraction (last access: 1 April 2023), 1983. a
Jimenez, A., Crespo, A., Migoya, E., and Garcia, J.: Advances in largeeddy simulation of a wind turbine wake, J. Phys. Conf. Ser., 75, 012041, https://doi.org/10.1088/17426596/75/1/012041, 2007. a, b, c
Jimenez, A., Crespo, A., Migoya, E., and Garcia, J.: Largeeddy simulation of spectral coherence in a wind turbine wake, Environ. Res. Lett., 3, 015004, https://doi.org/10.1088/17489326/3/1/015004, 2008. a, b
Johlas, H. M., MartínezTossas, L. A., Churchfield, M. J., Lackner, M. A., and Schmidt, D. P.: Floating platform effects on power generation in spar and semisubmersible wind turbines, Wind Energy, 24, 901–916, https://doi.org/10.1002/we.2608, 2021. a
Jonkman, J., Butterfield, S., Musial, W., and Scott, G.: Definition of a 5MW Reference Wind Turbine for Offshore System Development, National Renewable Energy Laboratory (NREL), https://doi.org/10.2172/947422, 2009. a, b
Knoll, D. and Keyes, D.: Jacobianfree Newton–Krylov methods: a survey of approaches and applications, J. Computat. Phys., 193, 357–397, https://doi.org/10.1016/j.jcp.2003.08.010, 2004. a
Lanzilao, L. and Meyers, J.: An Improved FringeRegion Technique for the Representation of Gravity Waves in Large Eddy Simulation with Application to Wind Farms, Bound.Lay. Meteorol., 186, 567–593, https://doi.org/10.1007/s1054602200772z, 2022a. a, b, c, d, e, f
Lanzilao, L. and Meyers, J.: Effects of selfinduced gravity waves on finite windfarm operations using a largeeddy simulation framework, J. Phys.: Conf. Ser., 2265, 022043, https://doi.org/10.1088/17426596/2265/2/022043, 2022b. a
Larsen, G.: A simple wake calculation procedure. Tech. Rep. RisøM2760 Risø., Tech. rep., https://orbit.dtu.dk/en/publications/asimplewakecalculationprocedure (last access: 1 April 2023) 1988. a
Leonard, B.: A stable and accurate convective modelling procedure based on quadratic upstream interpolation, Comput. Meth. Appl. M., 19, 59–98, https://doi.org/10.1016/00457825(79)900343, 1979. a
Lilly, D. K.: A proposed modification of the Germano subgrid‐scale closure method, Phys. Fluids AFluid, 4, 633–635, https://doi.org/10.1063/1.858280, 1992. a
Loken, C., Gruner, D., Groer, L., Peltier, R., Bunn, N., Craig, M., Henriques, T., Dempsey, J., Yu, C.H., Chen, J., Dursi, L. J., Chong, J., Northrup, S., Pinto, J., Knecht, N., and Zon, R. V.: SciNet: Lessons Learned from Building a Powerefficient Top20 System and Data Centre, J. Phys. Conf. Ser., 256, 012026, https://doi.org/10.1088/17426596/256/1/012026, 2010. a
Lundquist, J., Duvivier, K., Kaffine, D., and Tomaszewski, J.: Costs and consequences of wind turbine wake effects arising from uncoordinated wind energy development, Nature Energy, 4, 25–34, https://doi.org/10.1038/s4156001802812, 2019. a
Maliska, C. R. and Raithby, G. D.: A method for computing three dimensional flows using nonorthogonal boundaryfitted coordinates, International Journal for Numerical Methods in Fluids, 4, 519–537, https://doi.org/10.1002/fld.1650040606, 1984. a
Maronga, B., Gryschka, M., Heinze, R., Hoffmann, F., KananiSühring, F., Keck, M., Ketelsen, K., Letzel, M. O., Sühring, M., and Raasch, S.: The Parallelized LargeEddy Simulation Model (PALM) version 4.0 for atmospheric and oceanic flows: model formulation, recent developments, and future perspectives, Geosci. Model Dev., 8, 2515–2551, https://doi.org/10.5194/gmd825152015, 2015. a
MartínezTossas, L. A., Churchfield, M. J., and Leonardi, S.: Large eddy simulations of the flow past wind turbines: actuator line and disk modeling, Wind Energy, 18, 1047–1060, https://doi.org/10.1002/we.1747, 2015. a, b, c, d, e, f, g, h, i
Meneveau, C., Lund, T. S., and Cabot, W. H.: A Lagrangian dynamic subgridscale model of turbulence, J. Fluid Mech., 319, 353–385, https://doi.org/10.1017/S0022112096007379, 1996. a, b
Meyers, J. and Meneveau, C.: Large Eddy Simulations of Large WindTurbine Arrays in the Atmospheric Boundary Layer, in: 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 4–7 January 2010, Orlando, FL, https://doi.org/10.2514/6.2010827, 2010. a
Min, M., Brazell, M., Tomboulides, A., Churchfield, M., Fischer, P., and Sprague, M.: Towards Exascale for Wind Energy Simulations, arXiv [preprint], https://doi.org/10.48550/ARXIV.2210.00904, 2022. a
Moeng, C.H.: A LargeEddySimulation Model for the Study of Planetary BoundaryLayer Turbulence, J. Atmos. Sci., 41, 2052–2062, https://doi.org/10.1175/15200469(1984)041<2052:ALESMF>2.0.CO;2, 1984. a, b
Monin, A. and Obukhov, A.: Basic laws of turbulent mixing in the surface layer of the atmosphere, Tr. Akad. Nauk SSSR Geophiz. Inst., 151, 163–187, 1954. a
Munters, W., Meneveau, C., and Meyers, J.: Shifted periodic boundary conditions for simulations of wallbounded turbulent flows, Phys. Fluids, 28, 025112, https://doi.org/10.1063/1.4941912, 2016. a, b
Niagara: Niagara, https://docs.alliancecan.ca/wiki/Niagara (last access: 1 April 2023), 2018. a, b
Niayifar, A. and PortéAgel, F.: Analytical Modeling of Wind Farms: A New Approach for Power Prediction, Energies, 9, 741, https://doi.org/10.3390/en9090741, 2016. a
Nygaard, N. G., Steen, S., Poulsen, L., and Pedersen, J. G.: Modelling cluster wakes and wind farm blockage, J. Phys. Conf. Ser., 1618, 062072, https://doi.org/10.1088/17426596/1618/6/062072, 2020. a
Nygaard, N. G., Poulsen, L., Svensson, E., and Grønnegaard Pedersen, J.: Largescale benchmarking of wake models for offshore wind farms, J. Phys.: Conf. Ser., 2265, 022008, https://doi.org/10.1088/17426596/2265/2/022008, 2022. a
OpenCFD: OpenFOAM – The Open Source CFD Toolbox – User's Guide – Version 6, OpenCFD Ltd., United Kingdom, https://www.openfoam.com/ (last access: 1 April 2023) 2018. a
Ørsted: Ørsted presents update on its longterm financial targets, https://orsted.com/en/companyannouncementlist/2019/10/1937002 (last access: 5 March 2023), 2019. a
Paulson, C. A.: The Mathematical Representation of Wind Speed and Temperature Profiles in the Unstable Atmospheric Surface Layer, J. Appl. Meteorol. Climatol., 9, 857–861, https://doi.org/10.1175/15200450(1970)009<0857:TMROWS>2.0.CO;2, 1970. a
Pedersen, J. G., Gryning, S.E., and Kelly, M.: On the Structure and Adjustment of InversionCapped Neutral Atmospheric BoundaryLayer Flows: LargeEddy Simulation Study, Bound.Lay. Meteorol., 153, 43–62, https://doi.org/10.1007/s105460149937z, 2014. a, b, c, d
Pedersen, J. G., Svensson, E., Poulsen, L., and Nygaard, N. G.: Turbulence Optimized Park model with Gaussian wake profile, J. Phys.: Conf. Ser., 2265, 022063, https://doi.org/10.1088/17426596/2265/2/022063, 2022. a
Ponce, M., van Zon, R., Northrup, S., Gruner, D., Chen, J., Ertinaz, F., Fedoseev, A., Groer, L., Mao, F., Mundim, B. C., Nolta, M., Pinto, J., Saldarriaga, M., Slavnic, V., Spence, E., Yu, C.H., and Peltier, W. R.: Deploying a Top100 Supercomputer for Large Parallel Workloads: The Niagara Supercomputer, in: Proceedings of the Practice and Experience in Advanced Research Computing on Rise of the Machines (Learning), PEARC '19, Association for Computing Machinery, New York, NY, USA, https://doi.org/10.1145/3332186.3332195, ISBN 9781450372275, 2019. a
PortéAgel, F., Lu, H., and Wu, Y.T.: A largeeddy simulation framework for wind energy applications, in: Fifth International Symposium on Computational Wind Engineering (CWE2010), 23–27 May 2010, Chapel Hill, North Carolina, USA, 2010. a
Rampanelli, G. and Zardi, D.: A Method to Determine the Capping Inversion of the Convective Boundary Layer, J. Appl. Meteorol., 43, 925–933, https://doi.org/10.1175/15200450(2004)043<0925:AMTDTC>2.0.CO;2, 2004. a, b, c, d, e, f, g
Rosenfeld, M., Kwak, D., and Vinokur, M.: Development of a fractionalstep method for the unsteady incompressible NavierStokes equations in generalized coordinate systems, NASA Technical Memorandum no. 103912, NASA, 1992. a
Saad, Y. and Schultz, M. H.: GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems, SIAM J. Sci. Stat. Comp., 7, 856–869, https://doi.org/10.1137/0907058, 1986. a
Schneemann, J., Rott, A., Dörenkämper, M., Steinfeld, G., and Kühn, M.: Cluster wakes impact on a fardistant offshore wind farm's power, Wind Energy Science, 5, 29–49, https://doi.org/10.5194/wes5292020, 2020. a
Schumann, U.: Subgrid lengthscales for largeeddy simulation of stratified turbulence, Theor. Comput. Fluid Dyn., 2, 279–290, 1991. a
Segalini, A.: An analytical model of windfarm blockage, J. Renew. Sustain. Ener., 13, 033307, https://doi.org/10.1063/5.0046680, 2021. a
Sescu, A. and Meneveau, C.: A control algorithm for statistically stationary Large Eddy Simulations of thermally stratified boundary layers, Q. J. Roy. Meteor. Soc., 140, 2017–2022, https://doi.org/10.1002/qj.2266, 2014. a
Shaw, W. J., Berg, L. K., Debnath, M., Deskos, G., Draxl, C., Ghate, V. P., Hasager, C. B., Kotamarthi, R., Mirocha, J. D., Muradyan, P., Pringle, W. J., Turner, D. D., and Wilczak, J. M.: Scientific challenges to characterizing the wind resource in the marine atmospheric boundary layer, Wind Energ. Sci., 7, 2307–2334, https://doi.org/10.5194/wes723072022, 2022. a
Smith, R. B.: Gravity wave effects on wind farm efficiency, Wind Energy, 13, 449–458, https://doi.org/10.1002/we.366, 2010. a
Sockeye: UBC ARC: Sockeye – Detailed Technical Specifications, https://arc.ubc.ca/sockeyetechspecs (last access: 1 April 2023), 2023. a
Sørensen, J. N. and Shen, W. Z.: Numerical Modeling of Wind Turbine Wakes, J. Fluids Eng., 124, 393–399, https://doi.org/10.1115/1.1471361, 2002. a, b
Sørensen, J. N., Mikkelsen, R. F., Henningson, D. S., Ivanell, S., Sarmast, S., and Andersen, S. J.: Simulation of wind turbine wakes using the actuator line technique, Philos. T. R. Soc. A, 373, 20140071, https://doi.org/10.1098/rsta.2014.0071, 2015. a
Stevens, R. J., MartínezTossas, L. A., and Meneveau, C.: Comparison of wind farm large eddy simulations using actuator disk and actuator line models with wind tunnel experiments, Renew. Energ., 116, 470–478, https://doi.org/10.1016/j.renene.2017.08.072, 2018. a, b, c, d, e, f, g, h, i, j, k, l
Stieren, A., Gadde, S. N., and Stevens, R. J.: Modeling dynamic wind direction changes in large eddy simulations of wind farms, Renew. Energ., 170, 1342–1352, https://doi.org/10.1016/j.renene.2021.02.018, 2021. a
Stipa, S., Ajay, A., and Brinkerhoff, J.: Toolbox fOr Stratified Convective Atmospheres (TOSCA), OSF [code], https://doi.org/10.17605/OSF.IO/Q4VAF, 2023. a
Stull, R.: Practical Meteorology: An Algebrabased Survey of Atmospheric Science, BC Open Textbook Collection, AVP International, University of British Columbia, ISBN 9780888651761, https://books.google.nl/books?id=xP2sDAEACAAJ (last access: 1 April 2023), 2016. a
The HDF Group: Hierarchical data format version 5, http://www.hdfgroup.org/HDF5 (last access: 1 April 2023), 2006. a
Troldborg, N. and Meyer Forsting, A.: A simple model of the wind turbine induction zone derived from numerical simulations, Wind Energy, 20, 2011–2020, https://doi.org/10.1002/we.2137, 2017. a
van der Laan, M. P., Storey, R., Sørensen, N., Norris, S., and Cater, J.: A CFD code comparison of wind turbine wakes, J. Phys. Conf. Ser., 524, 012140, https://doi.org/10.1088/17426596/524/1/012140, 2014. a
van Kan, J.: A SecondOrder Accurate PressureCorrection Scheme for Viscous Incompressible Flow, SIAM J. Sci. Stat. Comp., 7, 870–891, https://doi.org/10.1137/0907059, 1986. a
Wu, K. L. and PortéAgel, F.: Flow Adjustment Inside and Around Large FiniteSize Wind Farms, Energies, 10, 2164, https://doi.org/10.3390/en10122164, 2017. a, b
Yang, X. I. A., Park, G. I., and Moin, P.: Loglayer mismatch and modeling of the fluctuating wall stress in wallmodeled largeeddy simulations, Phys. Rev. Fluids, 2, 104601, https://doi.org/10.1103/PhysRevFluids.2.104601, 2017. a
 Abstract
 Introduction
 Methodology
 Validation
 CNBL simulations with different controllers
 Finite wind farm with thermal effects
 Conclusions
 Appendix A: Numerical procedure
 Appendix B: LES modeling
 Appendix C: Finite wind farm setup
 Appendix D: TOSCA parallel scaling
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.
 Article
(7346 KB)  Fulltext XML
 Abstract
 Introduction
 Methodology
 Validation
 CNBL simulations with different controllers
 Finite wind farm with thermal effects
 Conclusions
 Appendix A: Numerical procedure
 Appendix B: LES modeling
 Appendix C: Finite wind farm setup
 Appendix D: TOSCA parallel scaling
 Code availability
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References