Mesh-adaptive simulations of horizontal-axis turbine arrays using the actuator line method

Numerical models of the flow and wakes due to turbines operating within a real-scale offshore wind farm can lead to a prohibitively large computational cost, particularly when considering blade-resolved simulations. With the introduction of turbine parametrizations such as the actuator disk (AD) or the actuator line (AL) models, this problem has been partially addressed, yet the computational cost associated with these simulations remains high. In this work, we present an implementation and validation of an AL model within the mesh-adaptive three-dimensional fluid dynamics solver, Fluidity, under a unsteady Reynolds-averaged Navier-Stokes based turbulence modelling approach. A key feature of this implementation is the use of mesh optimization techniques, which allow for the automatic refinement or coarsening of the mesh locally according to the resolution needed by the fluid flow solver. The model is first validated against experimental data from wind tunnel tests. Finally, we demonstrate the benefits of mesh-adaptivity by considering flow past the Lillgrund offshore wind farm.


Introduction
The actuator line model (ALM), introduced almost 15 years ago by Sørensen and Shen [1], has been an increasingly popular approach for the parametrisation of horizontal-axis wind turbines within computational fluid dynamics (CFD) solvers. With the ALM a rotor is represented by a number of rotating lifting lines corresponding to the turbine's blades, for which tabulated airfoil polar data is used to compute the forces. The coupling between the turbine and the underlying fluid flow solver is achieved by projected the negative blade force back to the fluid with the aid of compactly supported momentum source terms. ALM has been proved to provide accurate solutions for both the near and the far wake fields, particularly when it is coupled with transient flow solvers. The accuracy of the model has been tested in a number of studies starting from [1] and followed by [2,3] amongst many others. More recently, a number of enhancements have been proposed in order to reproduce a momentum source which better approaches the actual blade behaviour [4,5], as well as the idea of coupling ALM with a dynamic stall model aiming to capture the response of a blade element to unsteady loading [6].
Turbine parametrisations models (TPMs) such as the actuator line model (ALM) and the actuator disc model (ADM) exhibit a large number of advantages compared to blade-resolved simulations, both in terms of their respective computational efficiency but also as far as their implementation within a computational fluid dynamics (CFD) solver is concerned. First, by using TPMs the number of degrees of freedom needed by the fluid solver is significantly reduced since the boundary layer of the individual blades is no longer required to be resolved. Second, the introduction of the momentum source to represent the motion of the blades circumvents the need to use either a rotating or an overlapping mesh strategy to capture the motion of the turbine rotor. These two factors have rendered the use of TPMs a computationally affordable alternative approach for the modelling of large-scale wind farms. Thus, the actuator disk model (ADM) was used by [7,8,9] to model the wake field and predict the power output of a large-scale wind farm while [10,11,12] undertook ALM simulations to solve for the wake field as well as to obtain statistics for blade loads. In all of these studies either a uniform mesh or a block-mesh strategy was employed within the presented simulations. Inherently, such an approach requires a priori knowledge of the wake length and width, or otherwise the assignment of a relatively large volume to be a refined "wake region". However, as wind farms are getting larger and as ever more realistic simulations are required, including changes in the wind (and therefore wake) direction, this leads to the requirement that the "refined wake regions" must be expanded to cover a larger proportion of the domain. This approach is not optimal for the discretisation of the domain as the overall computational problem size can be significantly increased even through only a small increase in the size of this pre-refined region. Therefore, some sort of flexible mesh adaptivity procedure (e.g. the dynamic mesh optimization approach used here) employed during the course of a simulation is an attractive approach to consider as we proceed towards longer and more realistic simulations of offshore wind farms. This particular need has already been expressed in previous studies. For instance, Churchfield et al. [12] mentions that ". . . adaptive mesh refinement would be useful in providing higher resolution only where necessary, but may incur a run-time penalty in performing the refinement and the processor load balance . . . ". In a similar note, Nilsson et al. [9] also used different resolution grids and pointed out that the simulations using the most refined grids were abandoned due to limitations on the available computational resource.
In an effort to resolve the issues related to the optimum use of the computational resources, we present here the implementation and validation of an actuator line model which employs dynamic mesh optimization techniques. The optimization of the mesh is achieved through a strategy which allows control over both the numerical error and mesh size at run time [13,14,15]. A similar approach was first undertaken by [16] and more recently by [17] who both developed turbine models (actuator volume and actuator disk based respectively) and combined them with adaptive mesh methods to resolve the flow field. The former also undertook a real-scale wind farm simulation by considering the case of the Lillgrund offshore wind farm [18]. The present AL model is implemented within the open source code Fluidity [14,19] which is a general purpose finite element solver with the additional ability to conduct optimization based dynamic adaptivity of the underlying unstructured mesh. To model the fluid flow we have adopted here an unsteady Reynoldsaveraged Navier-Stokes (uRANS) based approach, combined with the k-ω SST turbulence model [20]. Before proceeding to the large-scale simulations, the accuracy of the new ALM implementation in Fluidity is investigated through comparisons with a series of wind tunnel tests [21,22] for the power and wake of one and two turbines operating in line. The model shows very good agreement with the rotor's thrust and power coefficients predictions as well as the wake field. These comparisons give us confidence that the model can predict the wake characteristics with high accuracy when real-world scale wind farms are considered. To demonstrate the efficiency of the mesh optimization approach we compare the results from an adaptive-mesh simulation with those from a static pre-refined mesh simulation for the Lillgrund offshore wind farm and compare to data from the literature [9].
The paper starts with section 2 which introduces the mesh-adaptive fluid solver and section 3 which discuss the turbine parametrisation and its coupling to the fluid solver, respectively. In Section 4, the model's two mesh approaches (fixed vs. adaptive) are validated against laboratory scale experimental data [21,22]. Next in section 5 wake prediction from the Lillgrund offshore wind farm are presented using an upscaled version of the "Blind Tests" turbine. The mesh optimization techniques are presented from the point of view of the same solver (Fluidity) and its contribution to increasing accuracy and reducing computational cost are finally discussed thereafter in section 6.

Unsteady RANS formulation of the governing equations
The fluid flow is modelled using the unsteady Reynolds-averaged Navier-Stokes (uRANS) equations, in which the velocity is decomposed into mean and fluctuating (turbulent) components. These equations take the form where u is the time-averaged component of velocity, u is the fluctuating velocity component, p is the mean pressure, µ is the dynamic viscosity of the fluid, and F T is a momentum source term computed at each time step from the AL model. The term −∇ · (ρu ⊗ u ) is a residual term from the application of the time-averaging operator on the outer product of the fluctuating velocity components, called the Reynolds stress tensor τ R . The presence of the Reynolds stress tensor in (2) introduces a number of additional unknowns and leads to a lack of closure in the system of equations. To address this closure problem we employ the Boussinesq approximation in which the Reynolds stresses are related to the time-averaged turbulence kinetic energy and strain-rate tensor as where k = 1 2 u · u is the time-averaged turbulence kinematic energy, µ T is the eddy viscosity and I is the unit tensor. Here, we make use of the standard k-ω SST turbulence model to compute the time-averaged turbulence kinetic energy and eddy viscosity proposed by Menter [20]. The closure model is a blend of two turbulence models, combining the classic k-ω model of Wilcox [23] to capture the flow separation in near-wall regions and the classic kmodel [24] for the regions of free shear flow. In the k-ω SST model the turbulent energy dissipation is transformed to the turbulent frequency ω ∝ /k and the two transport equations are expressed as and ρ ∂ω ∂t where is the eddy viscosity and F 1 is the blending function taking the value 1 for the Wilcox k-ω model and 0 for the standard k-model. The tildered quantitỹ P k denotes a limiting turbulence kinetic energy production and is given bỹ which is applied to prevent the build-up of turbulent energy in stagnation regions as proposed more recently by Menter [20]. The closure coefficients σ k , σ ω , α, β and β * which are selected for the analysis are found through an interpolation from the blending function value F 1 of the form where, η 1 , η 2 denote two dummy coefficients for the two limiting cases of F 1 = 0 and F 2 = 1 respectively. The limiting coefficients that are selected for this analysis are σ k1 = 0.85, σ k2 = 1.0, σ ω1 = 0.5, σ ω2 = 0.856, α 1 = 5/9, α 2 = 0.44, β 1 = 0.075, β 2 = 0.0828 and β * = 0.09. The details for the blending function F 1 can be found in the appendix of Abolghasemi et al. [17] and the references therein.

Numerical implementation
The system of the governing equations together with the additional scalar transport equations required for the turbulence modelling have been implemented within the open-source code Fluidity [19]. Fluidity is a general purpose three-dimensional, unstructured mesh, finite element/control volume based PDE solver [25,13,15] with the ability to make use of dynamic mesh adaptivity. For our analysis, the ensemble-averaged continuity and momentum equations are discretized using mixed finite elements for which we use a linear discontinuous basis functions for the velocity and a continuous quadratic basis functions for the pressure field over a tetrahedral element (P1 DG -P2).This scheme is known to be LBB stable [26], and it behaves well for advection-dominated problems [27]. In addition, the formulation uses a slope limiter [28] to ensure a robust solution for the velocity and pressure fields. The k-ω SST model employs a control volume discretisation [15] which is shown to prevent an oscillatory behaviour of the turbulent kinetic energy k, while for time marching the second-order accurate Crank-Nicolson scheme is used and is combined with an additional explicit sub-cycling approach used for the advection operators [19].

Mesh optimization
The underlying unstructured mesh of the fluid domain is subject to optimization-based adaptivity algorithms which are used in order to improve the quality of the mesh as well as provide higher or lower resolution at locations which are identified by the solver. For instance, introducing the actuator lines' momentum source to the fluid domain will create a requirement for a particular edge length over an assigned volume. This is achieved here with the use of a scalar field which identifies the region that the rotating actuator lines will occupy during the simulation (usually a disk of diameter which slightly larger than the diameter of the rotor). At the same time two additional fluid properties, the velocity vector field and the turbulence kinetic energy (a scalar field) are also used to guide the mesh-optimization process. While a detailed description of the optimization process is beyond the scope of this paper, we have included a short description of the process herein. The present mesh-optimization strategy is guided by the Hessian of selected prognostic and diagnostic fields [29]. The method considers a series of scalar fields φ (including potentially the individual components of the velocity field) from which an interpolation theory based error estimate on an element ∆ is derived. Next, we may define an element-valued Hessian H and |H| the positive-definite metric formed by taking the absolute value of the eigenvalues of H. Finally, we define the metric tensor field where u is a desired 'interpolation error' defined by the user. A metric is formed for each of the selected adaptivity fields (velocity, TKE and actuator line scalar field) which are then superimposed to form a single metric. In addition, in order to restrain the number of elements to a reasonable amount we impose a maximum number of mesh nodes, a maximum and a minimum edge lengths as well as an element gradation parameter. A series of mesh optimization operations are then performed in order to obtain a mesh of equilateral tetrahedra with unit edge lengths when measures with respect to the metric. Finally, mesh to mesh interpolation is used to transfer solution data from the old to the new mesh. The process is generally repeated every specified number of time steps. More information about the particular implementation used here can be found in [29,13,14,15] and the references therein.

Actuator Line Model
The actuator line model implementation presented as part of this work is based on the work of Sorensen and Shen [1]. We start the implementation of our ALM formulation by substituting the rotor's blades with moving (rotating) actuator lines which are discretised into blade elements (foil sections with finite span-wise length). The analysis of each blade element is conducted independently, assuming the local incident fluid velocity U = (U x , U y , U z ) and solid body velocity U b = (U bx , U by , U bz ), as shown in figure 1. The incident fluid velocity U is directly extracted from the fluid solver Fluidity by evaluating the finite element solution at the midpoint of every blade element. The solid body velocity is computed from the motion of the actuator lines, and can either be known a priori or computed at run-time. For instance, a horizontal axis turbine rotating with a constant prescribed rotational velocity field having non-flexible blades will experience a solid body velocity U b = Ω × r on a blade element where Ω is the rotor's rotational velocity, and r is the vector between the blade cross-section mid-point and the rotor's root. The normal and tangential components of the local blade velocity are computed as where n and t are the local element's normal and tangential (chord-wise) unit vectors respectively. We may now calculate the relative velocity magnitude from and the effective angle of attack can be calculated from α E = tan −1 U rn /U rt .
Having computed the relative velocity of each blade element, the hydrodynamic forces and moments on each blade section will be evaluated in terms of the normal and tangential forces, and the pitching moment coefficients C n , C t , C M s 25 defined as as well as the blade element area A and chord length c. The details on calculating the force coefficients are presented below 3.1.

Static foil Data
The force and moment coefficients of equations (12) are calculated using tabulated static lift, drag and pitching moment coefficients (C L , C D , C M 25 ) as functions of the local chord Reynolds number Re c = U r c/ν (where ν is the fluid's kinematic viscosity) and the effective angle of attack α E . These coefficients are generally extracted from experimental data or computational models and in many cases can be a priori corrected to account for threedimensional effects. In the present ALM, look-up procedures are used to obtain static foil coefficients from predefined tables, while corrections for rotational and three-dimensional effects are made by applying a correction factor to the lift coefficient (see section 3.2). The final lift, drag and pitching moment coefficients are then used to compute the respective normal and tangential force coefficients as follows,

Blade end effects
As mentioned before, the static foil data were obtained from two-dimensional numerical or experimental procedures and therefore will have to be corrected in order to account for blade end effects. According to Prandlt [30] in a real rotor the circulation near the blade tip tends exponentially to zero. To reproduce such effects in our ALM computations, we consider the tip loss correction model of Shen et al. [31]. In their approach the correction coefficient is computed as where the parameter g is given by B is the number of the rotor blades, λ the rotor's tip speed ratio, R the rotor radius, r the distance of a particular element from the tip, φ is the angle formed between the axis of rotation and the direction of the element's local relative velocity, and c 1 , c 2 are two empirically obtained coefficients given the values c 1 = 0.125 and c 2 = 21, respectively. Alternatively, the blade end effects measured can be accounted by extracting the lift and drag coefficients from blade-resolved simulations as it was recently suggested by [32].

Tower modelling
Finally, the tower behind the turbine will also have an effect on both the near and the far wake fields. Schumann et al. [33] showed that the tower wake is entrained into the rotor wake thanks to the swirling components of the velocity. As a result the velocity profiles in the near wake region exhibit an asymmetry. The method adopted herein to model the wake of the tower is the same as in [34] and [35] -the tower is also represented by an actuator line with elements corresponding to cylindrical blade elements. The drag coefficients for the cylinder are assumed to be constant during the simulations and the value C D = 1 is assigned. On the other hand, the lift coefficient C L varies with time, in an effort to reproduce the von Karman street behind the cylinder, via where f = 0.2 × U ∞ /D tower is the Strouhal number. An alternative method would be to consider the exact geometry of the tower with a boundary fitted mesh, and retain the actuator lines for the blades only as suggested by [4].

Actuator line to mesh interpolation
The final computed forces on the actuator line nodes are subsequently projected onto the fluid mesh and represented through momentum source terms. To ensure a smooth transition from a concentrated point AL force to the mesh, a smoothing interpolation functions needs to be employed. Traditionally in AL models, the momentum source terms F T are computed by multiplying the local force f AL with the smoothing function and adding a minus sign to account for the action-reaction principle. For a particular choice of smoothing interpolation function this may be written as Here, |r| is the distance of the mesh point from the AL node and is a smoothing parameter selected after taking into account the mesh size, drag force and chord size: where, V elem is the volume of the element in which the actuator line node lies, c the element's chord size and C D the drag coefficient of the respective blade element. This selection process was demonstrated in [5] to yield a smoothing parameter which minimizes the error between the point force and the flow field past an airfoil for linearised 2D conditions. In a deviation from the traditional interpolation method, we shall project the forces using the radial basis function (RBF) methodology. The RBF approach has long been used for solving multi-variable scattered data problems [36]. To demonstrate the method, let us assume that a single actuator line consists of a set of points X = {x 1 , . . . , x N } in the three-dimensional space and let F T := f AL on X, then where y is the coordinate vector of any fluid mesh node, Φ j = η(|y − x|) is computed from the Gaussian smoothing kernel and s j is an interpolation coefficient. In order to find the s j a linear system of equations needs to be solved by specifying that the interpolation conditions hold on the actuator line nodes. In this work, the RBF method proved to provide a better interpolation scheme for the unstructured tetrahedral mesh and therefore was used for all simulations.

Model Validation
To validate the newly implemented actuator line model we use data from previous wind tunnel experiments. These are experimental data results from a series of "Blind Test" workshops organized by the Norwegian Research Centre for Offshore Wind Technology (NoWiTech) and the Norwegian Centre for Offshore Wind Energy (NoCOWE) and which took place in the wind tunnel facilities of the Norwegian University of Science and Technology (NTNU) [21,22]. The reported data has been used as a benchmark solution for the validation of many turbine parametrisation models [34,37] and thus it was selected for the validation of the present model as well. We will refer to the two papers that reported the data as "blind tests" or when mentioned individually as "blind test 1" (BT1) and "blind test 2" (BT2). The wind tunnel facility used for the two blind tests is 11.14 m long, 2.71 m wide, and 1.8 m high. In (BT1) a single turbine with diameter D = 0.894 m with a hub height equal to H hub = 0.817 m is placed at a distance of 2D from the wind tunnel inlet. The turbine has three blades which consist of 14% NREL S826 airfoils and is designed so that its optimum operation tip speed ratio is equal to λ = 6. For this case, power and thrust coefficients are reported by [21] for a large range of tip speed ratios (5 to 11.5) while wake statistics including the stream-wise velocity deficit and the turbulence kinetic energy are reported for λ = 6. On the other hand, BT2 involves wake predictions from two turbines operating in-line. Two similarly designed turbines with slightly different diameters D 1 = 0.944 m and D 2 = 0.894 m and the subscripts 1 and 2 correspond to the "Front" and the "Rear" turbines respectively) are placed in the wind tunnel. The front turbine is located at 2D 2 = 1.788 m from the tunnel's entrance while the rear is placed at a distance 3D 2 = 2.686 m from the first as shown in figure 2. The experimental study reports three different operational scenarios for the two turbines, a tip speed ratio λ 1 = ΩR/U ∞ = 6 for the front, and λ 2 = 4, 7 and 2.5 for the rear. In both BT1 and BT2 the free-stream velocity in the wind tunnel is taken to be 10 m s −1 and a low ambient turbulence intensity of approximately I = 0.3% is assumed for both cases. In the computations we use the turbulence intensity I, the free stream velocity and the length scale, and calculate the initial and boundary conditions for the turbulence kinetic energy k = 1.3 × 10 −3 m 2 /s 2 , and a turbulent frequency ω = 0.5 s −1 .

Mesh convergence study
For simplicity, we have chosen the design power coefficient as a representative quantity for our mesh size and time step convergence study. The mesh convergence study involves six set-up cases in which we choose between the two mesh types (uniform fixed, adaptive) and vary the minimum element edge length h = 0.1 m, 0.075 m and 0.05 m. Here, all adaptive simulations were conducted using an adaptation period T adapt = 2.5 s. Table 1 shows the six set-up cases along with the average number of elements used by the adaptive meshes during the simulations. In addition, we have conducted a time step convergence study in which three time step sizes ∆t = 0.02 s, 0.01 s and 0.005 s, which correspond to 25, 50 and 100 time steps per rotor revolution respectively, are used for the higher resolution cases (h = 0.05 m) for both the fixed and adaptive mesh cases. To accurately predict the static hydrodynamic coefficients for the 14% NREL S826 section we use the experimental data obtained by [34] which include data spanning Reynolds number from  with the information about the sections used for the supporting structure (tower). Finally, the turbine performance and thrust coefficients are computed based on the uniform upstream velocity U ∞ as well as the the nominal radius R = D/2 via where P is the turbine's extracted power computed from the rotor's torque magnitude P = QΩR. To obtain time-averaged values for the performance and thrust coefficients which have reached converged values we run long enough simulations to allow the turbines to undertake more than 100 revolutions (approximately 8 seconds in actual time). Results from the convergence study are shown in figure 4. On the left of figure 4 the relative error in the predicted power coefficient is plotted against the number of elements used in the simulations and for which the minimum time step ∆t = 0.005 s is used. We observe that for a given minimum element edge length, the two approaches yield very similar power coefficient relative errors. However, the adaptive mesh method will always require a lower number of elements due to the optimum and flexible use of the underlying mesh within the computational domain. Therefore, the minimum edge length used is the primary factor determining the accuracy of the model, while adaptivity is used primarily to reduce the overall number of degrees of freedom, and therefore the associated computational cost. Later we will be able to demonstrate the advantage of mesh-adaptivity, by undertaking large-scale wind farm simulations. In that case the substantially fewer number of elements required during the spin-up period of the wake development leads to an important saving in CPU time. On the other hand, the temporal convergence study (right-hand side of figure 4) shows that as we reduce the time step (and thus increase the number of time steps used during one rotor revolution), the accuracy of the model converges to the element-based maximum obtainable accuracy. Again, the plotted figures correspond to simulations using the minimum edge length h = 0.05 m. Therefore, the element edge length and the time step used in all simulations presented hereafter, will be based upon this preliminary convergence study. We should also notice, that once we have selected an element edge length for our analysis, the magnitude requirement on the time step is entirely dictated by the actuator line method and not the stability of the fluid solver, which allows for a far more relaxed condition assuming a Courant number of near unity. Finally, we should mention that the smoothing parameter varies with the element's edge length and for all these simulations is taken equal to 2.5 times the edge length, which also satisfies equation (18).

Blind Test 1 and 2
Having established that the element edge length determines the accuracy of our model we will present the wake predictions from BT1 for all six set-up cases of table 1using the smallest time step ∆t = 0.005 s, while the power and thrust coefficients and the wake predictions of BT2 using only the fine mesh-adaptive results, in order to maintain clarity. We start with the power and thrust coefficients from the blind tests. For the wake field, we present three horizontal profiles downstream of the single turbine at x/D = 1, 3 and 5 for BT1 and the three horizontal profiles downstream of the rear turbine at x/D = 1, 2.5 and 4 for only the first scenario (λ 2 = 4) from BT2. Figures 6 and 7 show the mean stream-wise velocity and the turbulence kinetic energy (TKE) for BT1, while figures 8 and 9 show the mean stream-wise velocity and the stream-wise turbulent stresses u u for BT2. All quantities are time-averaged after a spin-up period which is taken to be approximately 1 sec and have been normalised by the upstream velocity U ∞ and presented  as a function of the normalised horizontal distance y/R. Next we consider the second Blind Test (BT2). The power and thrust coefficient predictions from the ALM for the second turbine for the three scenarios are shown in table 2. Similarly for resolving the wake field of the turbines operating inline only the refined adaptive case is used. For the sake of brevity we present the wake predictions only from scenario 1. We should notice here that in order to obtain the Reynolds stress u 2 we make use of the isotropy turbulence relation k = 3u 2 /2. Such an assumption is not appropriate when the turbulence stresses are highly anisotropic. This may explain the large discrepancies in figure 9 for x/D = 1. However, we may in general accept that it provides a very good estimate for the other two downstream profiles. The results from both BT1 and BT2 provide us with confidence that the proposed adaptive methodology is both faster and more accurate than the fixed mesh one, assuming that the same solver is used with both approaches. However,  applying mesh adaptivity in laboratory scale (wind tunnel) experiments cannot fully demonstrate its potential. This is due to the high-blockage present in the wind tunnel experiments and which necessitates the use of the finest mesh resolution almost everywhere in the domain. To better demonstrate the effect and potential benefits of mesh adaptivity we present in the next sections simulations from the Lillgrund offshore wind farm and discuss some options that can further reduce the computational time.

Parametrisation of the wind farm
The Lillgrund offshore wind farm is located between Copenhagen in Denmark and Malmö in Sweden and consists of 48 Siemens SWT93-2.3 MW (Siemens Wind Power, Hamburg, Germany) turbines. The farm harnesses the south-westerly winds coming from a wind fetch as large as 50 km and has a layout made up from 8 rows and 8 columns named from 1 to 8 and A to H respectively as shown in figure 10. For our simulations, the coordinate axis system (x-y) has been rotated so that it aligns with the statistically dominant wind direction SW-NE (43 • / 223 • ). To parametrise the turbines making up the wind farm we need to make a number of assumptions as their blade geometry and airfoil characteristics are not publicly available. The only information available to us is the turbines' rotor radius R = 46.5 m, the hub height H hub = 65 m as well as the thrust and power coefficients as a function of the wind speed [39]. For the rest of the turbine parameters we shall follow the approach of Nilsson et al. [9] who considered a downscaled version of the conceptual National Renewable Energy Laboratory (NREL) 5 MW turbine as presented by [40] and confirmed its suitability by comparing the thrust and power (through torque) output for different wind speeds. For our simulations however, we consider an up-scaled version of the "Blind Test" turbine which we will need to scale up by approximately a factor of 100. To also confirm the validity of our choice, we present in figure 11 the thrust and power coefficients for velocities varying from 5 m/s to 11 m/s similar with [9]. So far in section 4.2, we have expressed both the thrust and power coefficient as a function of the tip speed ratio λ. Given that no active control of the rotor or the blade is taken into account, we will assume that a stand-alone turbine will operate optimally by adjusting its angular velocity in accordance with the optimum tip speed ratio λ = 6 as shown in figure 5. The results of this assumption are shown in figure 11, which will be assumed a reasonable approximation for the purpose of our simulations.

Mesh strategies
Turning now to the mesh strategy selected for the simulations, we begin by introducing the two approaches of fixed (pre-refined) versus dynamically adaptive meshing. As mentioned before, the problem at hand is inherently multi-scale in which a number of length scales need to be resolved. These scales span from individual turbine wake's turbulence length scale, usually determined by the modelling approach (LES, uRANS), to the far larger atmospheric meso-scales usually in the order of O(200 m). With this in mind, To compare the computational cost of the two approaches we consider the case of an upstream wind speed U 0 = 8 m/s aligned with the x-coordinate axis (223 • ). The level of turbulence intensity is chosen to be as high as TI = 5.7% while all turbines operate with the optimum tip speed ratio as described earlier in this section. The simulations were performed for a case with a pre-refined region (PRR) in which the finest resolution of h =20 m was pre-assigned in an area of 4 km × 6.5 km × 0.16 km (see fig. 12), and three adaptive simulations in which different frequencies of mesh adaptations (T adapt = 2.5 s, 5 s and 10 s) were applied. It should be noted here that the selection of the pre-refined region (PRR) was performed through a "trial and error" approach as the width and the height of the individual turbine wakes were not known a priori. Therefore our goal was to create a domain which has the refined region within the wake region but does not extend excessively beyond this purpose. An alternative approach of course would be to assume no a priori knowledge and simply use the minimum edge length over the entire domain, but this would of course result in a huge problem size. For mesh adaptive simulations the frequency with which mesh adaptations are performed in time is selected from a larger number of user-defined parameters available to us such as the elements' maximum aspect ratio, the minimum/maximum element ratio etc. This option is considered in order to understand its impact on reducing the overall CPU time in the present uRANS set-up. The selection of mesh adaptivity frequency as a key parame-ter is motivated by the computational cost (runtime penalty) associated with it, particularly during the interpolation process. Indeed, mesh adaptations can be seen as a potential bottleneck step in the computation procedure, as they require information to be passed from all processors and the application of relatively costly interpolation methods to be used (e.g. Galerkin projection) in order to transfer the information from the previous mesh to the new one. The application of a Galerkin projection method is necessitated here by the use of a discontinuous function space, as the recent study of Farrell and Maddison [41] showed that a consistent interpolation is not suitable for discontinuous fields. In addition, other parameters such as the element edge length size h, have been excluded as in section 4 we showed that there is a direct correlation between the size of the underlying mesh and the model's accuracy. Thus, changing the frequency with which these adaptations occur is the only remaining parameter which can be varied in order to reduce the overall CPU time while maintaining a similar accuracy for the turbine performance and wake predictions. The latter will also be re-assessed later in this section. Moving on to the simulation set-up, we apply free slip conditions on all boundaries as well as a uniform velocity U 0 = 8 m/s profile and a uniform turbulence intensity TI = 5.7 % (by setting the turbulence kinetic energy to k = 0.3 m 2 /s 2 ) to the initial field and at the flow inlet. To obtain a quasi-steady solution for the far wake field, we time-march the solution with a constant time step of ∆t = 0.1 s, which corresponds to more than 100 time steps within a rotor revolution for all turbines while also remains well bellow a unity Courant number and therefore satisfies the stability criteria for the fluid solver. The solution is time-marched for a total number of 4000 time steps which corresponds to more than 3 times the time required by the uniform flow to travel the full length of the computational domain considered here. The solutions obtained at the final time step are shown in figures 13 and 14, for both the fixed mesh and the adaptive mesh using the smallest adaptation period T adapt = 2.5 s, respectively. The plots include the underlying mesh, the stream-wise velocity u x and the turbulence intensity T I shown in a horizontal cross-section slice at hub height H hub = 65 m. In addition to the wake field, the power estimates normalised by the median of the production of the front row turbines P md similar to Nilsson et al. [9], are presented compared against the measured data extracted from [39] in figure  15. Again the power production of the model P is considered at the final time step. Looking at the model's accuracy, we begin by noticing that the two solu-  tions obtained by the adaptive and the fixed mesh simulations exhibit very similar behaviour both qualitatively and quantitatively as is shown in figures 13 and 14. An important feature in both simulations is the dissipation of the inlet TKE which was due to the assignment of a turbulence frequency required for the dissipation of the turbine wakes. The dissipation process also created a moving front for the (TKE) which disappears after about 2000 time steps. The latter also affected the adaptive mesh since the TKE field was used within the definition of the metric controlling the mesh optimization procedure. This TKE inlet front had an important impact on the mesh elements and will be discussed in more detail in the next subsection. Finally, the power production trends are also found to agree well with the measured  Figure 15: Results from the fixed-mesh (black line) and the mesh-adaptive (red line) simulations for turbines in Row A-H compared against the measured data (symbols) of [39]. The power production is normalised with the average of the first row turbines as reported by [39].
ones as is shown in figure 15. The only major discrepancy observed in the plots (Rows A-H) is the over-prediction of the relative power of the second turbine. This result is believed to be related to the operation strategy adopted in the simulations for the individual turbines, which is based on adjusting the angular velocity of the turbine to the optimum tip speed ratio by using the rotor-averaged velocity as the reference value, and not taking into account control strategies that may be used in the wind farm operation such as active torque or blade pitch control. Similar trends were observed by [9] who adopted a similar approach for the operation, although the discrepancies between experimental values and their LES results appear to be much smaller. Nevertheless, the power production as predicted by the three adap-tive simulations at the final time level are essentially identical and therefore only one of the three is plotted in figure 15. The mean relative error over all turbines between the CFD simulations' power prediction and the measured data amounts to 7.38% for all three adaptive-mesh simulations while 8.02% for the fixed-mesh ones. This result, confirms our initial hypothesis that varying the frequency of mesh adaptations will not affect the accuracy of the model predictions in the uRANS approach employed here. On the other hand, a small difference between the fixed-mesh and the adaptive-mesh simulations power predictions is observed.

Computational efficiency
Returning to the main objective of our investigation which is to examine the computational efficiency of the two proposed mesh strategies. It was hypothesized that the computational cost of the adaptive simulations should be lower than that for the fixed pre-refined mesh simulations, as well as the computational cost becoming smaller as the adaptation frequency is increased while the same or similar accuracy is achieved when the same edge lengths are used. To test our hypothesis, we examined the development of the wakes within the Lillgrund offshore wind farm, and conducted simulations using both a pre-refined mesh and three adaptive meshes with different adaptation periods (T adapt = 2.5 s, 5 s and 10 s). The simulations were run in parallel using MPI on a total of 80 processing cores (4 nodes each of 20 cores) on the cx1 cluster at Imperial College London. Each case was given a total wall clock time of 96 hours in which the pre-refined simulations required approximately 86 hours to complete. We should mention here that the absolute value of the required CPU time would differ for different algorithms or code implementations. For this reason, we present the adaptive simulation CPU time normalised by the pre-refined one. Nevertheless, we begin by presenting the evolution of the fluid mesh elements with time as shown in figure 16. First, we observe that the number of elements follow a similar trend for all three adaptive cases. The simulations begin with a spin-up of the mesh in the vicinity of the turbines and a moving front starting from the inlet. The latter, is a result of the dissipation of the turbulent kinetic energy in the inlet and will start disappearing after about 50 s. After this spin-up period and between 50 s to 100 s the number of elements are slightly reduced due to the "dissipation" of the moving front after which the element counts starts to increase again with a constant rate after 150 s. Finally, the number of elements becomes constant after about 380 s at which time a steady state solution has been reached for the far field. The discussed mesh evolution is also depicted in the horizontal and vertical mesh slices shown in figures 17 and 18. In addition, from figure 16 we also observe that the curves of the three adaptive simulations remain well below the pre-refined simulation line until end of the simulations. Inherently, we may argue that since the adaptive-mesh simulations always use a smaller number of elements than the fixed-mesh ones as well as Fluidity generally scaling approximately linearly with the number of elements, then the adaptive-mesh simulations should require a much lower CPU time. However, as discussed previously, an additional cost "penalty" is present due to the mesh adaptation procedure itself. In addition, there is also potential for unequal balancing of the elements across the MPI processors (although a dynamic load balancing step is incorporated within the parallel mesh optimization procedure), as well as the fixed number of cores utilised throughout the simulation not necessarily being optimal from a parallel scalability perspective. Thus, the resulting CPU time is significantly impacted by the frequency of these mesh adaptations. Starting with the pre-refined simulations, a total of 6880 CPU hours where required for the simulation to complete, while 5641.6, 4974.24 and 4747.19 CPU hours were required for the adaptive simulations using mesh adaptation periods of (T adapt = 2.5 s, 5 s and 10 s), respectively. Firstly, we notice that the CPU time is reduced by 18% through the use of mesh adaptivity even when making the most frequent mesh adaptation operations. Second, by reducing the frequency of mesh adaptations additional reductions in the CPU time are achieved: 27.7 and 31 % respectively. We should notice however, that the selection of 80 processors while it can be an optimum in terms of scalability in the case of pre-refined region, it is not always for the adaptive ones. This is due to the fact that during the spin-up period, a relatively small number of elements will be distributed across a large number of processors and the computational cost will be dominated by the MPI communications between the partitions. Thus an additional step which could be taken to optimise the CPU usage would be to start the simulations with a smaller number of processors (e.g. 20) and as the overall problem size grows start and stop the simulation (using checkpoints) in order to gradually increase the number of processors to 80. It should be further emphasised that these cost savings are in the context of a static pre-refined mesh where some a priori knowledge of the flow dynamics was used to generate that mesh. In the case that this knowledge was lacking, or with a time-varying wind field, this 'optimal mesh generation' would be more complex to achieve and the adaptive mesh performance increase would be expected to be more substantial.

Discussion
We have presented the implementation and validation of a mesh-adaptive actuator line (AL) model which can be used to optimise the number of elements/cells that are required to resolve the wake field within a real-scale wind farm. Our actuator line model implementation is based on the original model of [1] with additional models being used to account for blade end effects and the impact of the tower shadow on the near wake field. The fluid flow is resolved using a uRANS formulation of the governing equations and the k-ω SST turbulence closure model. The developed Fluidity-ALM model was validated against laboratory scale data before attempting full-scale wind farm simulations. For the wind farm simulations we selected the Lillgrund offshore wind farm for which numerous previous studies and performance measurements are available. The effect of mesh-adaptivity on reducing the computational cost and accuracy is also discussed.
More specifically, for the developed model a detailed mesh convergence study was undertaken for both the fixed-mesh and the adaptive mesh simulations showing that for a given edge length h the two mesh strategies (fixed and adaptive) will yield similar results for the integrated rotor characteristics such as the power or thrust coefficient. Moreover, as the element's edge length is reduced, the model predictions for the power coefficient C P converge to the experimental reference values. Looking at the wake solutions on the other hand, mesh-adaptivity provides better estimates for the wake field, particularly in regions where large shear (velocity gradients) exist. This is inherently due to the ability of the mesh optimization algorithm to identify these regions and allow elements with large aspect ratios to align with the underlying velocity gradient. A fixed mesh does not have such an ability and therefore a better solution at these regions can be obtained only by further refining the mesh locally. The validation of the model using the two Blind Tests (BT1 & BT2) shows that both the near and far wake field can be accurately predicted given that enough resolution is used. More satisfactory however, are the predictions from the second validation case, BT2. Deep array wake modelling has been a great challenge for many wake models and it is a crucial step towards accurately predicting offshore wind farm power output. The results from the comparison with BT2 give us confidence that, although only two turbines in-line are used, we are capable of obtaining high-fidelity wake solutions and closely predict the power curtailment with an error being estimated as low as 9.92% for scenario 1 of BT2.
Looking at the Lillgrund offshore wind farm simulation results, the power predictions from each row (A-H) agree well with measured data from [39]. The largest discrepancy between the model predictions and the measured data is presented for the second turbine in each row. This is systematically appearing in all simulations, and is believed to be due to the operational scenario chosen for the second turbine. Nevertheless, the two modelling approaches using the fixed and the adaptive-mesh strategies, yield almost identical results for the power coefficients re-affirming our conclusions from the mesh convergence study which states that the turbine predictions are not dependent on the mesh strategy but rather the edge-length used at the location of the rotor. As far as the wake field is concerned, although there is lack of data for the same of comparison, the accuracy of the power coefficients predictions suggests that the model was able capture the magnitude of the wake-deficits along each row, at least in a time-averaged sense. Nonetheless, as far as the mesh-adaptive simulations are concerned, adaptivity was found to be a favourable choice leading to a reduction in the overall computational cost while maintaining the same accuracy. Finally, by changing the period with which mesh adaptation was invoked, T adapt , we were able to further reduce the overall computational cost without compromising the model's accuracy.
Finally, we should also mention the limitations of our modelling approach. In the current study, we investigated the potential of mesh adaptivity to reduce the computational costs of large-scale wind farm simulations by comparing two meshing strategies, and ignored a number of parameters/phenomena that will in general affect the operation of the wind farm. To begin with, a uniform inlet velocity with slip conditions on both the bottom and the top of the domain were considered, instead of a no-slip condition at the bottom, allowing for the development of a boundary layer. In addition, buoyancy effects and atmospheric boundary layer (ABL) stability were also ignored for the sake of simplicity. However, a number of studies [10,42,12,11] have suggested that wake recovery is highly dependent on the stability of the ABL and therefore our model will also need to be tested under such conditions. However, by conducting our investigation from a point of view of a uRANS solver, we are unable to capture many of the key effects occurring within the ABL such as the sporadic flow bursts the intermittent turbulence loading on the turbines or the effect of small upwind eddies driven by the temperature difference. For such simulations higher-accuracy LES simulations should be used which can capture smaller time scales and resolve the turbulence kinetic energy within the arrays. Inherently, a large-eddy simulation solver which uses mesh-adaptivity will impose additional constraints on the frequency of mesh adaptations and the elements' aspect ratios, while a far smaller edge lengths and time step will be required. All these additional factors put the efficiency of mesh adaptivity in question and additional studies will be required to assess its applicability. Creech et al. [18] studied the Lillgrund offshore wind farm using the Fluidity solver with mesh-adaptivity and an LES formulation and found good agreement with the SCADA (Supervisory control and data acquisition) system data. However, their study did not examine the efficiency of mesh-adaptivity. Finally, the actuator line model chosen for the parametrisation of the turbines within the fluid domain, while it exhibits a number of advantages particularly in modelling the near wake field, was also found to be particularly computationally expensive for the nature of our simulations. A recent study by Deskos et al. [43] showed that the actuator line method can be two orders of magnitude more computationally demanding, largely due to the need to use a smaller time step for the rotor, while blockage effects may also be under-predicted.
Future planned work will include an examination of large scale offshore and onshore wind farms using large eddy simulation solvers (LES) where the atmospheric boundary layer (ABL) will also be resolved. Mesh adaptivity should also find value in resolving the boundary layer over an uneven terrain and studying the interaction of terrain-generated turbulence with wind turbines wakes.