Mesh-adaptive simulations of horizontal-axis turbine arrays using the actuator line method
RReceived 21 October 2017; Revised 01 June 2018; Accepted 14 June 2018DOI: xxx/xxxx
RESEARCH ARTICLE
Mesh-adaptive simulations of horizontal-axis turbine arraysusing the actuator line method
Georgios Deskos* | Matthew D. Piggott
Department of Earth Science and Engineering,South Kensington Campus, Imperial CollegeLondon, SW7 2AZ, UK
Correspondence *G. Deskos, Department of Earth Science andEngineering, South Kensington Campus,Imperial College London, SW7 2AZ, UK.E-mail : [email protected] InformationThis research was supported by theEngineering and Physical Sciences ResearchCouncil (EPSRC), Grant Numbers:EP/R007470/1 and EP /L000407/1.
Abstract
Numerical models of the flow and wakes due to turbines operating within a real-scale offshorewind farm can lead to a prohibitively large computational cost, particularly when consideringblade-resolved simulations. With the introduction of turbine parametrisations such as the actua-tor disk (ADM) or the actuator line (ALM) models this problem has been partially addressed, yetthe computational cost associated with these simulations remains high. In this work we presentan implementation and validation of an ALM within the mesh-adaptive 3D fluid dynamics solver,Fluidity, under a uRANS-based turbulence modelling approach. A key feature of this implemen-tation is the use of mesh optimization techniques which allow for the automatic refinement orcoarsening of the mesh locally according to the resolution needed by the fluid flow solver. Themodel is first validated against experimental data from wind tunnel tests. Finally, we demonstratethe benefits of mesh-adaptivity by considering flow past the Lillgrund offshore wind farm.
KEYWORDS: actuator line model; mesh optimization; uRANS turbulence models; Lillgrund offshore wind farm
Turbine parametrisation models (TPMs) such as the actuator line model (ALM) and the actuator disk model (ADM) exhibit a large number of advan-tages compared to blade–resolved simulations, both in terms of their respective computational efficiency but also as far as their implementationwithin a CFD solver is concerned. First, by using TPMs the number of degrees of freedom needed by the fluid solver is significantly reduced sincethe boundary layer of the individual blades is no longer required to be resolved. Second, the introduction of the momentum source to representthe 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.Hence, the ADM was used by (1, 2, 3, 4, 5) to model the wake field and predict the power output of operating offshore wind farms (e.g. Lillgrund,Horns Rev) while (6, 7, 8) undertook ALM simulations to solve for the wake field as well as to obtain statistics for blade loads. In all of these studieseither a uniform mesh or a block-mesh strategy was employed within the presented simulations. Such an approach requires a priori knowledge ofthe wake length and width, or otherwise a large volume of the computational domain to be a assigned as the refined “wake region”. Moreover, asmore realistic simulations are required for utility–scale wind farms, including changes in the wind (and therefore wake) direction, the “refined wakeregions” will need to be expanded to a greater extend in order to provide the required resolution of such simulations. Undoubtedly, this approachis not optimal for the discretization of the domain, as even a moderate expansion (in the order of a few decades of metres) of the wake region cansignificantly increase the number of degrees of freedom (DoF). Inherently, some sort of flexible mesh adaptivity procedure (e.g. the dynamic meshoptimization approach used here) which is employed during the course of the simulation would be an attractive approach to consider for the abovedescribed problem. This particular need has already been expressed in previous studies. For instance, Churchfield et al. (8) mentions that “...adaptivemesh refinement would be useful in providing higher resolution only where necessary, but may incur a run-time penalty in performing the refinement and the a r X i v : . [ phy s i c s . f l u - dyn ] O c t DESKOS AND PIGGOTT processor load balance ... ”. In a similar note, Nilsson et al. (4) also used different resolution grids and pointed out that the simulations using the mostrefined grids were abandoned due to limitations on the available computational resource.Combining turbine parametrisation models with mesh–optimisation techniques was first considered by (9) and more recent implementations canalso be found in (10) and (11), each employing a different turbine parametrisation model (actuator volume and actuator disk and actuator line based,respectively). Mesh–adaptivity has been used in conjunction with blade–resolved simulations by (12, 13, 14) who employed an unstructured over-lapping mesh strategy to obtain the near–body solution and either unstructured or Cartesian mesh–adaptive solvers for the off–body flow. Suchdual mesh approaches have been effectively used to calculate wind turbine and rotorcraft wakes and it has been postulated that the computationalefficiency of the approach is thanks to the off–body dynamic mesh–adaptivity solver. Indeed, a significant amount of computational resources canbe saved by using an optimal local refinement or coarsening of the mesh while at the same time maintaining the desired levels of solution accuracy.In the context of large–scale wind farm simulations, mesh–adaptivity was used by (15) and (16) to simulate the wakes developed within the Lillgrundoffshore wind farm. Both studies demonstrated the ability of a mesh–adaptive solver (either unstructured or Cartesian–based) to be coupled withturbine parametrisation models or overlapping mesh strategies, respectively, and to be used as a high–fidelity, multiscale wind farm modelling tool.In the above mentioned mesh–optimization/adaptivity algorithms and studies the obtained solutions have not been compared with traditionalstatic mesh solutions using the same solver in order to provide a rigorous estimate of either the potential accuracy gains and/or reductions incomputational effort. Such questions are important in shedding light on the efficiency and accuracy of mesh–adaptive solvers, and more specificallyon their applicability to wind energy research. To partially address these questions, we present herein the implementation and validation of an ALMwhich employs dynamic mesh optimization techniques. The optimization of the mesh is achieved through a strategy which allows control over boththe numerical error and mesh size at run time (17, 18, 19). Both the ALM and mesh–adaptivity approach are developed within the open source codeFluidity (18, 20) which is a general purpose unstructured mesh based finite element solver. In addition, the fluid flow is modelled in this work usingan unsteady Reynolds-averaged Navier-Stokes (uRANS) based approach, combined with the k – ω SST turbulence model (21). Before proceedingto the large-scale simulations, the accuracy of the new ALM implementation in Fluidity (combined with the uRANS configuration) is investigatedthrough comparisons with a series of wind tunnel tests (22, 23) for the power and wake of one and two turbines operating in line. The model showsvery good agreement with the rotor’s thrust and power coefficients predictions as well as the wake field. These comparisons give us confidence thatthe model can predict the wake characteristics with high accuracy when real-world scale wind farms are considered. To demonstrate the efficiencyof the mesh optimization approach we compare the results from an adaptive-mesh simulation with those from a static pre-refined mesh simulationfor the Lillgrund offshore wind farm and data from the literature (24).The paper starts with section 2 which introduces the mesh-adaptive fluid solver and section 3 which discusses the implementation of the turbineparametrisation. In Section 4, the two mesh approaches (fixed vs. adaptive) are validated against the experimental data of (22, 23) while in section 5simulation of the Lillgrund offshore wind farm are undertaken. The mesh optimization techniques are presented from the point of view of the samesolver (Fluidity) and its contribution to increasing accuracy and reducing computational cost are finally discussed thereafter in section 6.
The wind turbine wakes are modelled using the unsteady Reynolds-averaged Navier-Stokes (uRANS) equations, in which the velocity is decomposedinto mean and fluctuating (turbulent) components. Within uRANS the governing equations take the form ∇ · u = 0 , (1) ρ ∂ u ∂t + ρ u · ∇ u = −∇ p + µ ∇ u − ∇ · ( ρ u (cid:48) ⊗ u (cid:48) ) + F T , (2)where u is the time-averaged component of velocity, u (cid:48) is the fluctuating velocity component, p is the mean pressure, µ is the dynamic viscosityof the fluid, and F T is a momentum source term computed at each time step from the ALM. The term −∇ · ( ρ u (cid:48) ⊗ u (cid:48) ) is a residual term from theapplication of the time-averaging operator on the outer product of the fluctuating velocity components, called the Reynolds stress tensor τ R . Thepresence of the Reynolds stress tensor in (2) introduces a number of additional unknowns and therefore the Boussinesq approximation is adoptedto provide closure to the system of equations. That is, the Reynolds stresses are related to the time-averaged turbulence kinetic energy k = u (cid:48) · u (cid:48) and strain-rate tenso r via −∇ · ( ρ u (cid:48) ⊗ u (cid:48) ) = τ R = − kρ I + µ T (cid:18) ∇ u + ( ∇ u ) T (cid:19) . (3) ESKOS AND PIGGOTT where µ T is the eddy viscosity and I the unit tensor. To compute k we make use of the standard k − ω SST model proposed by Menter (21) whichrequires the solution of two additional transport equations: ρ ∂k∂t + ρ u · ∇ k = ∇ · (cid:18)(cid:0) µ + µ T σ k (cid:1) ∇ k (cid:19) + ˜ P k − ρβ ∗ kω, (4)and ρ ∂ω∂t + ρ u · ∇ ω = s./b ∇ · (cid:18)(cid:0) µ + µ T σ ω (cid:1) ∇ ω (cid:19) + (cid:18) ραµ T (cid:19) ˜ P k − ρβω + 2 (cid:18) − F (cid:19) ρσ ω ω ∇ k ∇ ω, (5)where µ T = ρ kω , (6)is the eddy viscosity and F a blending function defined in (21). The tildered quantity ˜ P k denotes a limiting turbulence kinetic energy productiongiven by ˜ P k = min( P k , ρβ ∗ ω ) , (7)which is applied to prevent the build-up of turbulent energy in stagnation regions (25). Finally, the closure coefficients σ k , σ ω , α , β and β ∗ areselected by linear interpolation using the blending function value F . Further information for the model implementation can be found in theappendix of (10) and the references therein. The system of governing equations, including the additional scalar transport equations required for the turbulence modelling, have been discretizedwithin the open-source code Fluidity (20). Fluidity is a general purpose three-dimensional, unstructured mesh, finite element/control volume basedPDE solver (26, 17, 19) with the ability to make use of optimization based anisotropic mesh adaptivity. For our analysis, the continuity and momen-tum equations are discretized using mixed finite elements for which piecewise-linear discontinuous basis functions are used to represent velocity,while continuous piecewise-quadratic basis functions are used for pressure over tetrahedral elements (the so-called P1 DG –P2 element pair). Thisscheme is known to be a Ladyzhenskaya-Babuska-Brezzi (LBB) stable combination (27) and to perform well for advection-dominated problems (28).The formulation also uses a slope limiter (29) to ensure a robust solution for the velocity and pressure fields. The k – ω SST model makes use of a con-trol volume based discretization (19) with flux limiters to help prevent oscillatory behaviour of the turbulent kinetic energy k . For time marching thesecond-order accurate Crank-Nicolson scheme is used and is combined with an additional explicit sub-cycling approach for momentum advection(20).The underlying unstructured tetrahedral mesh is also subject to optimization-based adaptivity algorithms which are used in order to improvethe quality of the mesh as well as provide higher or lower resolution at locations which are identified by the solver. For instance, the introductionof the ALM momentum source will create a requirement for a particular edge length over the assigned rotor volume. This is achieved here throughthe specification of a scalar field which identifies the region that the rotating actuator lines will occupy during the simulation. At the same time twoadditional fluid properties, the velocity vector field and the turbulence kinetic energy (a scalar field) are also used to guide the mesh-optimizationprocess. This is achieved through the derivation of a metric tensor field. If we consider a single scalar field φ which we want to adapt our mesh tooptimally resolve, we form a metric tensor, M φ , by first computing the Hessian, H φ , of that scalar field and defining M φ = 1 (cid:15) φ | ¯ H φ | , (8)where (cid:15) φ is a user-defined weight for field φ ( (cid:15) φ can in some sense be considered a requested error – a smaller value leading to a larger M φ andconsequently a finer mesh). The Hessian encodes information about the curvature of the scalar field and includes both spatial and direction infor-mation; we desire finer mesh resolution at location and in directions where curvature is high, and coarser resolution where the solution is close tolinear. | ¯ H φ | indicates that we are interested in the magnitudes of curvatures when deciding on optimal mesh resolution, and do not care about thesign. Once metric tensor fields have been calculated for all the scalar fields (for the velocity vector we consider each component separately) we wishto adapt to, a final metric is obtained through superimposition of the individual metrics (30). At this stage additional constraints on the total numberof elements in the calculation, and/or maximum and minimum edge length, and the maximum rate at which edge lengths can vary in space can all beincorporated. The metric can then be used to measure the length of vectors – primarily the edge lengths of elements. A perfectly optimized meshin physical space is defined as one which is made up of unit length edges in metric space. The inhomogeneous and directionally dependent natureof the Hessians and hence the metrics thus leading to a mesh which is variable in both space and direction (i.e. is potentially anisotropic). An opti-mization functional is defined which measures how well the current mesh achieves this goal, and a series of topological operations are performedon the current mesh to improve this agreement. These operations include edge collapse and splitting, and face to face and face to edge swapping(30, 19). Finally, conservative mesh to mesh interpolation is used to transfer solution data from the old to the new mesh (31, 32). The entire meshoptimization process is conducted every user-defined number of time period, termed the adaptation period T adapt . Further information on the meshoptimization process, including its parallelisation, may be found in (30, 17, 18, 33, 19) and the references therein. DESKOS AND PIGGOTT
Our actuator line model (ALM) follows the approach of (34, 35) in which the turbine blades are represented by rotating virtual lines – the actuatorlines (AL). Point forces are computed along each AL, (at each blade element’s midpoint) using the relative velocity U extracted from the fluid solverby evaluating the globally-defined finite element solution at these points, the solid body velocity U b of each point, and the lift and drag coefficientsobtained from look-up tables using airfoil data for the blade element’s respective profile. Extra care has been taken for the elements near the bladetip for which the tip loss correction model of (36) is used. The tower behind the turbine is also incorporated by adopting a model similar to (37) – thetower is represented by an actuator line with element having a cylindrical cross-section of constant drag coefficient C D = 1 and a time–dependentlift coefficient C L which is “tuned” in order to reproduce the Von Kármán street behind the cylinder, C L = A sin(2 πft ) , (9)where A is the amplitude of the “dynamic lift force” and f = 0 . × U ∞ /D tower is the Strouhal number, based on the uniform velocity and thediameter of the tower D tower . For all simulations presented in this work the hub and nacelle of the turbines are not modelled. The final ALM forcesare projected onto the fluid mesh and represented in the governing equations via the momentum source term F T . To ensure a smooth transitionfrom a concentrated point AL force f AL to the source term F T , a smoothing interpolation function is used, F T = − (cid:15) π / exp (cid:18) − | r | (cid:15) (cid:19) f AL . (10)where | r | is the distance of the mesh point from the AL node and (cid:15) is a smoothing parameter selected after taking into account the mesh size, dragforce C D , chord size c and V elem is the volume of the element in which the AL node lies, (cid:15) = max (cid:20) c , (cid:112) V elem , cC D (cid:21) , (11)consistent with the recommendations of (38). The newly implemented ALM is validated using data from a series of “Blind Test” workshops organized by NoWiTech and NoCOWE and which wasobtained from the wind tunnel facilities of NTNU (22, 23). We will refer to the two papers that reported the data as “blind tests” or when mentionedindividually as “blind test 1” (BT1) and “blind test 2” (BT2). The wind tunnel facility used for the two blind tests is .
15 m long, .
72 m wide, and . high. In (BT1) a single turbine with rotor diameter D = .
894 m , hub height equal to H hub = .
817 m is placed at a distance of D from the windtunnel inlet. The turbine has three blades consisting of % NREL S826 airfoils, a tower (support structure) with diameter D tower = .
11 m and isdesigned to operate optimally for a tip speed ratio of λ = . For BT1, power and thrust coefficients are reported by (22) for a large range of tip speedratios ( to . ) while wake statistics including the stream-wise velocity deficit and the turbulence kinetic energy for λ = 6 . On the other hand,BT2 involves wake predictions from two similar turbines (with slightly different diameters D = .
944 m and D = .
894 m where the subscripts1 and 2 correspond to the “Front” and the “Rear” turbines respectively) operating in- line. The front turbine is located at D = .
788 m from thetunnel’s entrance while the rear is placed at a distance D = .
686 m from the first as shown in figure 1 . BT2 reports three operational scenarios,in which the front turbine operates with λ = Ω R/U ∞ = , and the rear with λ = , and . . In both BT1 and BT2 the mean free–stream velocitywas measured to be
10 m s − whereas the vertical profile of the ambient turbulence intensity found to be nearly uniform with a value of I ≈ . % .For our computations, we also assume a uniform velocity and TKE profiles with initial and inlet conditions U =
10 m s − and k = . × − m / s respectively while the turbulent frequency is also considered to be uniform and equal to ω = . − . Finally, the tower model is enabled for all thewind tunnel simulations with the lift coefficient given a value of A = . , consistent with the study of (37). The design power coefficient from the single turbine experiments (BT1) is used as a representative quantity for our mesh size and time step con-vergence studies. For the mesh convergence study we consider six set-up cases by assuming two mesh types (uniform fixed, adaptive) and threeminimum element edge lengths h = . , .
075 m and .
05 m as shown in table 1 , while the same time step ∆ t = .
005 s is used for all simulations.We should note here that the adaptive simulations are also restricted by a maximum edge length of . . Additionally, all adaptive simulations wereconducted using an adaptation period T adapt = . . On the other hand, the time step convergence study considers three time step sizes ∆ t = .
02 s , .
01 s and .
005 s , corresponding to approximately , and time steps per rotor revolution, and using the finest mesh ( h = .
05 m ) for boththe fixed and adaptive mesh. The lift and drag coefficients for the
NREL S826 section are taken from (39) and include data spanning a Reynoldsnumber from × to × and angles of attacks from − to as shown in figure 2 . The turbine performance and thrust coefficients are ESKOS AND PIGGOTT D X=1D X=3D X=5D L = 11 . m ≃ . D xz zy W = 2 . m ≃ D X=03D2D X=1D X=2.5D X=4DX=0 xzBT1BT2 L = 11 . m ≃ . DH = 1 . m ≃ D FIGURE 1
Schematic representation from the two “Blind Tests” set ups.
TABLE 1
Table of simulation set-up cases.Case Min. Mesh size h [m] Average. Num. Elem. Mest Type Tag01 0.1 234,098 Fixed FC02 0.075 561,611 Fixed FM03 0.05 1,699,635 Fixed FF04 0.1 169,939 Adaptive AC05 0.075 351,688 Adaptive AM06 0.05 698,236 Adaptive AFcomputed based on the uniform upstream velocity U ∞ and the nominal radius R = D/ via C T = T ρπR U ∞ , (12) C P = P ρπR U ∞ . (13)To obtain mean values for the power and thrust coefficients we run long enough simulations and allow the turbines to undergo more than 100revolutions. In figure 3 the relative error in the predicted power coefficient is plotted against the average number of elements used in the simula-tions. We observe that for a given minimum element edge length, the two approaches (fixed and adaptive mesh) yield very similar relative errors.However, the adaptive mesh method required a smaller number of elements due to the optimum and flexible use of the underlying mesh within thecomputational domain. Therefore, within the ALM/uRANS configuration considered here, the minimum edge length in the mesh can be considered DESKOS AND PIGGOTT − − − . . . C L − − . . . . . C D Re = 4 · Re = 6 · Re = 8 · Re = 1 · Re = 1 . · Re = 1 . · Re = 2 · Re = 3 · Re = 4 · FIGURE 2
Lift and drag coefficients as a function of the angle of attack (AoA) (39) used for the model validation. . . . . . P model − P exp P exp fixed meshadaptive mesh
25 50 75 100 1250 %5 %10 %15 %20 % Time steps per revs fixed meshadaptive mesh
FIGURE 3
Convergence study with respect to spatial and temporal resolution.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 further demonstrate the potential advantages of mesh–adaptivity by undertakinglarge-scale wind farm simulations. In that case the substantially fewer elements required during the spin-up period of the wake development leadsto an important saving in CPU time. On the other hand, the temporal convergence study (right-hand side of figure 3 ) shows that as we reduce thetime step (and thus increase the number of time steps used during one rotor revolution), the accuracy of the model converges to the element-basedmaximum obtainable accuracy. Again, the plotted figures correspond to simulations using the minimum edge length h = .
05 m . Therefore, the ele-ment edge length and the time step used in all simulations presented hereafter, will be based upon this preliminary convergence study. We shouldalso note that once we have selected an element edge length for our analysis, the magnitude requirement on the time step is entirely dictated bythe ALM and not the stability of the fluid solver, which allows for a far more relaxed condition assuming a Courant number of near unity. Lastly, thesmoothing parameter (cid:15) varies with the element’s edge length and for all simulations hereafter is taken equal to . times the edge length, which alsosatisfies equation (11). Having established that the element edge length determines the accuracy of our model we will present the wake predictions from BT1 for all sixset-up cases described in table 1 using the smallest time step ∆ t = .
005 s , while results for the power and thrust coefficients as well as thewake predictions of BT2 are presented only the fine mesh-adaptive, in order to maintain clarity. Starting with the wake profiles, we present three
ESKOS AND PIGGOTT − − . − − . . . . y/R − u / U ∞ x/D = 1 − − . − − . . . . y/R − u / U ∞ x/D = 3 − − . − − . . . . y/R − u / U ∞ x/D = 5 FIGURE 4
Blind Test 1: Horizontal mean stream-wise velocity profiles at x/D = , and . The plotted lines correspond to: fixed coarse (FC) ( ),fixed medium (FM) ( ), fixed fine (FF) ( ), adaptive coarse (AC) ( ), adaptive medium (AM) ( ), adaptive fine (AF) ( ) and the symbols( ) to the experimental values reported by (22). − − . − − . . . − − − y/R T K E / U ∞ x/D = 1 − − . − − . . . − − − y/R T K E / U ∞ x/D = 3 − − . − − . . . − − − y/R T K E / U ∞ x/D = 5 FIGURE 5
Blind Test 1: Horizontal TKE profiles at x/D = 1 , and . Line colouring the same as in figure 4 .horizontal profiles downstream of the single turbine at x/D = , and for BT1 and the three horizontal profiles downstream of the rear turbineat x/D = , . and for only the first scenario ( λ = ) from BT2. Figures 4 and 5 show the mean stream-wise velocity and the TKE for BT1,while figures 6 and 7 show the mean stream-wise velocity and the stream-wise turbulent stresses u (cid:48) u (cid:48) for BT2. All quantities are time-averagedafter a spin-up period which is taken to be approximately , and have been normalised by the upstream velocity U ∞ and presented as a functionof the normalised horizontal distance y/R . Similarly, to resolve the wake field in BT2 only the refined adaptive case was used, and for brevity wepresent the wake predictions only from scenario 1. We should also note that in order to obtain the Reynolds stress u (cid:48) we make use of the isotropyturbulence relation k = 3 u (cid:48) / . Such an assumption is not appropriate when the turbulence stresses are highly anisotropic, which may explain thelarge discrepancies in figure 7 for x/D = 1 . However, better estimates are obtained for the other two downstream profiles. Flow anisotropy wasalso found to affect the predictions of the wake for BT1. This is an inherent inability of all turbulence models, particularly for the estimation of TKE(40). Looking at integrated rotor quantities such as the power and thrust coefficients, an overall good agreement is observed between the ALMpredictions an the wind tunnel measurements. Results are shown for both BT1 and BT2 in figure 8 . The only large discrepancy that can be observedin figure 8 is for scenario 2 of BT2 (rear turbine operating with λ = ). In that case, the power coefficient exhibits a discrepancy of 457.14 % incomparison with the measurements which can be attributed primarily to the ALM limitation rather than the modelling of the wake. Indeed, (23)reported this scenario as the most challenging one for turbine parametrisation models. This is due to a non-uniform span-wise pressure experiencedby the blades of the rear turbine, necessitating the use of blade-resolved simulations to accurately predict the lift and drag coefficients of theindividual blade elements. Still, the results from both BT1 and BT2 give us confidence that the proposed adaptive methodology is both faster andmore accurate than the fixed mesh one. However, applying mesh adaptivity in laboratory scale (wind tunnel) experiments cannot demonstrate itsfull potential. This is due to the high-blockage created in the wind tunnel which necessitates the use of the finest mesh resolution almost everywherein the domain. To better demonstrate the potential benefits of mesh adaptivity we present in the next sections simulations for the Lillgrund offshorewind farm and discuss some options that can further reduce the computational time. DESKOS AND PIGGOTT − − . − − . . . . y/R u / U ∞ x/D = 1 − − . − − . . . . y/R u / U ∞ x/D = 2 . − − . − − . . . . y/R u / U ∞ x/D = 4 FIGURE 6
Blind Test 2: Horizontal mean stream-wise velocity profiles at x/D = , . and . Only the adaptive fine ( ) solution and theexperimental data ( ) are shown. − − . − − . . . . · − y/R u u / U ∞ x/D = 1 − − . − − . . . . · − y/R u u / U ∞ x/D = 2 . − − . − − . . . . · − y/R u u / U ∞ x/D = 4 FIGURE 7
Blind Test 2: Stream-wise turbulent stress profiles at x/D = , . and . Only the adaptive fine ( ) solution and the experimentaldata ( ) are shown. The Lillgrund offshore wind farm is located near the southern coast of Sweden and has recently attracted the interest of the research communityas a benchmark for numeical model validation (8, 4, 15, 41). Model validation is achieved through comparison with supervisory control and dataacquisition (SCADA) system measurements and it primarily tests the ability of the model to predict power losses along a row of turbines for differentwind directions. Part of the SCADA data are accessible via an online technical report (24) while the complete information used for the validationcan be extracted from the validation studies just mentioned, and in particular from the studies of (4, 15). The wind farm consists of 48 SiemensSWT93-2.3 MW (Siemens Wind Power, Hamburg, Germany) distributed in 8 rows (A to H) as shown in figure 9 . For our simulations, we considerthe Southwestern statistically dominant wind direction ( ◦ / ◦ ) and we align the coordinate axis system ( x - y ) with the rows A, B etc resulting ina layout where the aligned row turbines have a 4.3 rotor diameter spacing. To parametrise the turbines making up the wind farm we need to makea number of assumptions as their exact blade geometry and airfoil characteristics are not publicly available. The available information does includethe turbines’ rotor radius R = . , the hub height H hub =
65 m as well as the thrust and power coefficients as a function of wind speed (24). Forthe rest of the turbine parameters we adopt the approach of (4) who considered a downscaled version of the conceptual NREL 5 MW turbine aspresented by (42) and confirmed its suitability by comparing the thrust and power (through torque) output for different wind speeds. Here, however,we consider an up-scaled version of the “Blind Test” turbine, scaled up by approximately a factor of . To confirm the validity of our choice, wepresent in figure 10 the thrust and power coefficients for velocities varying from / s to
11 m / s , similar with (4). So far in section 4.2, we haveexpressed both the thrust and power coefficient as a function of the tip speed ratio λ . As an active control strategy, we assume that the turbine willsimply adjust its angular velocity in accordance with the optimum tip speed ratio, which was found to be around λ = as shown in figure 8 . Thisassumption is justified in part by the good agreement between the upscaled “Blind Test” turbine and the Siemens SWT-93 manufacturer’s curve ESKOS AND PIGGOTT . . . λC P . . . . . λC T Exp. FrontExp. RearFluidity (AF) FrontFluidity (AF) Rear
FIGURE 8
Power and thrust coefficient curves as computed by the present model using the adaptive fine mesh and as measured by (22, 23) plottedagainst the tip speed ratio λ . y [ k m ] ABCDEFGHXY
FIGURE 9
Layout of the Lillgrund offshore wind farm. Rows A, B, etc. are aligned with the x –axis of the coordinate system.shown in figure 10 . Here, it is worth emphasizing that a recent investigation by (43) showed that when wake predictions are made the need formesh resolution is primarily driven by the velocity deficit, which inherently relates to the thrust coefficient. Therefore any turbine parametrisationthat accurately captures the thrust force would result in a similar computational mesh. For the simulations presented here two mesh strategies are investigated, namely the fixed (but pre-refined, i.e. variable resolution) and fully dynam-ically adaptive meshing. Considering that the problem at hand is inherently multi-scale and a number of length scales, spanning from the individualturbine wake’s turbulence length scale to the far larger atmospheric meso-scales usually in the order of O (
200 m ) , need to be resolved, the twoselected meshing strategies are the following: DESKOS AND PIGGOTT . . . . U (m/s) P o w e r c o e ffi c i e n t , C P Siemens SWT-93Upscaled “Blind Test” turbine . . . . . U (m/s) T h r u s t c o e ffi c i e n t , C T Siemens SWT-93Upscaled “Blind Test” turbine
FIGURE 10
Comparison of between the real power and thrust coefficient ( C P and C T ) and the up-scaled “Blind Test” turbine as a function of theupstream wind speed.
10 km k m k m k m Adaptive region
FIGURE 11
Schematic representation of the two mesh strategies considered for the Lillgrund offshore wind farm test case. The pre-refined regionis centred around the middle of the y -axis while it starts from the domain inlet and ends . before its outlet.1. a large domain of dimensions
10 km ×
10 km × is considered with a meso-scale resolution O (
200 m ) everywhere expect for an innerregion with prescribed resolution using an edge length equal to h =
20 m (in the subsequent discussions this is referred to as the pre-refinedregion – PRR) and,2. a large domain of the above mentioned dimensions and meso-scale initial resolution everywhere without the a priori chosen inner regionrefinement, but with mesh-adaptivity enabled and operating over the whole domain and a minimum allowed edge length of h =
20 m specified.The simulations were performed for a case with a pre-refined region (PRR) assigned in an area of × . × .
16 km (see fig. 11 ), and threeadaptive simulations in which different frequencies of mesh adaptations ( T adapt = . , and
10 s ) were applied. The selection of the pre-refinedregion (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 covering the wake region but does not extend excessively beyond this point.An alternative approach (which would arguably represent a fairer comparison to the adaptive approach) would be to assume no a priori knowledgeand simply use the minimum edge length over the entire domain, but this would of course result in a huge problem size. For the mesh adaptive
ESKOS AND PIGGOTT simulations the adaptation period is selected from a larger number of user-defined parameters available to us such as the elements’ maximumaspect ratio, the minimum/maximum element ratio etc. This option is considered in order to understand its impact on reducing the overall CPUtime in the present uRANS set-up. The selection of mesh adaptivity frequency as a key parameter is motivated by the computational cost (runtimepenalty) associated with it, particularly during the interpolation process. Indeed, mesh adaptations can be seen as a potential bottleneck step in thecomputations, as load-rebalancing represents a data migration overhead and the application of relatively costly interpolation methods to be used(e.g. Galerkin projection) to transfer the information from the previous mesh to the new one. The application of a Galerkin projection method isnecessitated here by the use of a discontinuous function space, as the recent study of (32) showed that a consistent interpolation is not suitablefor discontinuous fields. For the present simulations, we consider the wind speed to be constant and therefore the required frequency of the meshadaptation is not expected to vary with time, except for when the wakes have converged. However, if longer simulations are to be undertakenin which the wind speed is expected to vary with time, the frequency of the adaptations may be adjusted during the course of the simulation tobetter capture the wake dynamics. It should be noted of course that a varying wind direction would also necessitate a larger PRR in the fixed meshcase, leading to substantial additional computational costs. Finally, other parameters such as the element edge length size h , have been excludedas 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 thefrequency with which these adaptations occur is a key remaining parameter which can be varied in order to reduce the overall CPU time whilemaintaining 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, for the purpose of these simulations we assume uniform initial and inlet mean velocity and TKE profiles( U = − and k = .
31 m / s which corresponds to a turbulence intensity I = . % ). On all other boundaries, we apply free-slip velocity andzero gradients for all other quantities ( ∂/∂n ( u, v, w, k, ω ) = ). The selected inlet and boundary conditions considered herein are not representa-tive of the levels of shear that the turbine wakes will experience within a realistic atmospheric boundary layer. However, as a first approach, we willdemonstrate the benefits of mesh adaptivity by imposing a uniform incident velocity and TKE profile, and therefore based on the selected adapta-tion metrics (velocity and TKE) not need to refine outside the wake region. The accuracy of our choice is also discussed later, when the numericalresults are compared with the observed SCADA data and some discrepancies are observed. It should be remembered that the focus of this studyis to compare the PRR solutions with the mesh-adaptive ones, and the observed data is added only to show the overall performance of the model.Future studies will consider more complex scenarios in which the mesh will also be further controlled in a manner that allows for anisotropic gra-dation in the vertical direction only, outside the wake region. Nevertheless, to obtain a quasi-steady solution for the far wake field we time-marchthe solution with a constant time step of ∆ t = . and the solution for a total number of time steps. The final wake solutions are shown infigures 12 and 13 , for both the fixed and adaptive mesh (using the smallest adaptation period T adapt = . ) simulations. In addition, the powerestimates normalised by the median of the production of the front row turbines P md (similar to quantities considered by (4)) are presented and com-pared against the measured SCADA data in figure 14 . The two approaches’ (adaptive and fixed mesh) results exhibit very similar behaviour bothqualitatively and quantitatively as observed in figures 12 and 13 . An important feature in both solutions is the dissipation of the inlet TKE dueto the assignment of a turbulence frequency required for the dissipation of the turbine wakes. The increased dissipation creates a moving front forthe inlet TKE which however disappears after about time steps. This impacted on the adaptive mesh since the TKE field was used within thedefinition of the metric controlling mesh optimization. The impact of the TKE inlet front is discussed in more detail in the next subsection. Next, thepower production trends are also found to agree well with measurements as shown in figure 14 . The only major discrepancy observed in the plots(Rows A–H) is the over-prediction of the relative power of the second turbine. This discrepancy may be attributed to our selection of uniform inflowconditions, as the amount of computed shear and wake asymmetry will be underestimated. Similar trends were also observed by (8, 4, 15) whoadopted a log–law profile, although the discrepancies between the experimental values and their LES results appear to be much smaller. Neverthe-less, the power production as predicted by the three adaptive simulations at the final time level are essentially identical and therefore only one ofthe three is plotted in figure 14 . The mean relative error over all turbines between the CFD simulations’ power prediction and the measured dataamounts to . % for all three adaptive-mesh simulations while it takes the value . % for the fixed-mesh ones. This result confirms our initialhypothesis that varying the frequency of mesh adaptations will not affect the accuracy of the model predictions in the case of the uRANS approachemployed here. On the other hand, a small difference between the fixed-mesh and the adaptive-mesh simulations’ power predictions is observed;this can be attributed to the ability of the adaptive simulations to better resolve the wake field by applying an optimum element aspect ratio, as wasalso shown for the wind tunnel tests. DESKOS AND PIGGOTT
FIGURE 12
Fixed mesh simulations: From left to right a horizontal cross-section at hub height of the underlying mesh, the stream-wise velocity u x and the turbulent intensity T I are shown.
FIGURE 13
Adaptive mesh simulations ( T adapt = . ): From left to right a horizontal cross-section at hub height of the underlying mesh, thestream-wise velocity u x and the turbulent intensity T I are shown.
Returning to a key objective of our investigation, which is to examine the computational efficiency of the two proposed mesh strategies, it washypothesized that the computational cost of the adaptive simulations should be lower than that for the fixed pre-refined mesh simulations, as wellas the computational cost becoming smaller as the adaptation period is increased. In addition, the same or similar accuracy should be achievedwhen the same minimum edge length is used. To test our hypothesis, we conducted simulations using one fixed and three adaptive meshes. Thesimulations 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.The pre-refined simulations required approximately 86 hours to complete which is considered as a reference value. We should mention here thatthe absolute value of the required CPU time would differ for different algorithms or code implementations. For this reason, we present the CPUtimes for the adaptive mesh simulations normalised by the CPU time for the pre-refined mesh case. Here the term “CPU time” refers to the wall–clock time multiplied by the number of processing cores. This CPU time is dominated by the actual run time of the jobs and auxiliary computational
ESKOS AND PIGGOTT . . z/R R e l a t i v e p o w e r P / P m d . . z/R R e l a t i v e p o w e r P / P m d . . z/R R e l a t i v e p o w e r P / P m d . . z/R R e l a t i v e p o w e r P / P m d . . z/R R e l a t i v e p o w e r P / P m d . . z/R R e l a t i v e p o w e r P / P m d . . z/R R e l a t i v e p o w e r P / P m d . . z/R R e l a t i v e p o w e r P / P m d FIGURE 14
Normalised power computed by the fixed-mesh (black line) and the mesh-adaptive (red line) simulations for turbines in Row A–H andcompared against the measured data (dots) of (24)procedures such as the decomposition and re–partition of the mesh were found to result in trivial computational times in comparison. With this inmind, we start by presenting the temporal evolution of the number of fluid mesh elements in figure 15 ; the number of elements are seen to follow asimilar 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 startingfrom the inlet. The latter, is a result of the dissipation of the turbulent kinetic energy in the inlet and disappears after about
50 s . After this spin-upperiod, 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 elementcounts start to increase again with a constant rate after
150 s . Finally, the number of elements becomes constant after approximately
380 s at whichtime a steady state solution has been reached for the far field. A similar trend for the count of elements/DoF was also obtained by (16). In theirsimulations of the full Lillgrund wind farm they observed a strong linear spin–up curve due to wake transients, and subsequently a flattened peakwhere the wakes have started interacting with each other, and small variations in the mesh due to wake–wake and turbine–wake interactions werefound. Such fluctuations are not observed however in the dynamic mesh evolution of the present simulations, which can be attributed to the low–pass temporal filtering incurred by the uRANS equations. On the right hand side of the same figure, we have plotted the evolution of the relativepower (
P/P md ) for all the back-row turbines against the elapsed time, for one mesh-adaptive case ( T adapt = . ) and the PRR simulation. It can beobserved that the relative power of the individual turbines converges to near its final value long before
380 s . The evolution of the adaptive mesh ispresented via horizontal and vertical slices through the domain in figures 16 and 17 . Note that slices through completely unstructured tetrahedralmeshes do result in complex polygons and slivers, but these images do convey the spin up of the mesh as it resolves the developing wakes. Zoomedpictures of the vertical profiles are also shown in 17 to demonstrate the element gradation (particularly the higher aspect ratio of the elements) DESKOS AND PIGGOTT . . . · Time Elapsed (seconds) N u m b e r o f e l e m e n t s PRR simulationsAdaptive simulations T adapt = 2 . T adapt = 5 sAdaptive simulations T adapt = 10 s . . . . . R e l a t i v e p o w e r P / P m d PRR simulationsAdaptive simulations T adapt = 2 . FIGURE 15
Left: Number of tetrahedral elements used for the underlying mesh as a function of the elapsed simulation time, Right: Convergence ofthe relative power of the rear row of turbines – individual lines for each turbine are plotted here.
FIGURE 16
From left to right: A horizontal slice at hub height ( H =
65 m ) through the adaptive mesh at ,
100 s and
400 s , and having in total
100 760 , and elements in the computational domain, respectively.near the edge of the wake field. Element–wise, the plots (figure 15 ) of the three element counts for the adaptive mesh simulations remain wellbelow that for the pre-refined simulation line until the end of the simulations and thus we may argue that since the adaptive-mesh simulationsalways use a smaller number of elements than the fixed-mesh ones, and that the Fluidity solver scales approximately linearly with the number ofelements, then the adaptive-mesh simulations should require a much lower CPU time. For the adaptive mesh simulations, however, additional CPU“penalties” are imposed due to the various stages of the mesh adaptation procedure itself. Thus, the resulting CPU time can be significantly impactedby the frequency of these mesh adaptations. In addition, there is also potential for unequal balancing of the elements across the MPI processes(although a dynamic load balancing step is incorporated within the parallel mesh optimization procedure). The fact that a fixed number of processingcores is utilised throughout the simulation, and that this number is unlikely to be optimal from a parallel scaling perspective given the significantchanges in element count during the course of the adaptive mesh simulations, can also impact on computational efficiency gains. The actual CPUrequired by our numerical simulations confirm our initial hypothesis. For the pre-refined simulations a total of CPU hours were required, while . , . and . CPU hours were required for the adaptive mesh simulations using mesh adaptation periods of ( T adapt = . , and
10 s ), respectively. From these results, we may make two observations. First, by switching to the adaptive approach the CPU time is reduced by % even when making relatively frequent mesh adaptation operations. Second, by reducing the frequency of mesh adaptations additional reductionsin the CPU time are achieved: to . and % respectively. However, we should mention that the selection of processors while it can be optimal ESKOS AND PIGGOTT FIGURE 17
Left: A vertical slice along the x -axis passing through the centre of the domain ( y = ) for the adaptive mesh at ,
100 s and
400 s (from top to bottom) having in total
100 760 , and elements in the computational domain, respectively, and zoomed–in snapshotsfrom the final adaptive mesh ( t =
400 s ) at different locations. Right: A 3D view of the computational mesh and stream wise velocity u x at t =
100 s using a crinkle slice passing through hub height. Iso–contours of the magnitude of the turbine source term F T are also shown to identify the locationof each turbine. The iso–contour level is chosen small enough such that all turbines are visible. T=2.5 sec T=5 sec T=10 sec00 . . . . . C P U T i m e a d a p t C P U T i m e P RR Case A (1 run)Case B (2 runs)Case C (4 runs)
FIGURE 18
CPU time required by the adaptive simulations in the three scenarios (A,B and C), normalised by the CPU time required for the PRRsimulation.in terms of scalability in the case of the pre-refined region, it is not always so for the adaptive ones. This is due to the fact that during the spin-upperiod, a relatively small number of elements will be distributed across a large number of processors and the computational cost will be dominatedby the MPI communications between the partitions. In order to further optimize the CPU usage of the adaptive simulations, we have re-run thesame simulations but this time starting the parallel computation on a smaller number of processing cores (Nproc= ), and as the overall problemsize grows we start and stop the simulation (using checkpoints) and increase the number of cores to the final as the element counts increase.We should emphasize here that prior knowledge of the final mesh from the original adaptive simulations using processing cores was essential inbetter designing the decomposition of the domain when checkpoints were used. Here we present three cases, the first one (Case A) is the originalcase in which, the simulations were run on 80 processors for the whole time period, the second one (Case B) starts with processors, then afterthe first spin–up time ( t =
100 s ) the simulation is stopped and restarted with , while the last one (Case C) uses , , and at the periods to , to , to and to , respectively. The cumulative CPU time for the three cases is shown in figure 18 . It is observed thatcases B and C result in a smaller overall CPU time than case A. In addition, while in case C the decomposition/repartition is performed four times,this does not lead to smaller CPU times. Instead, case B seems to behave better, although by a small amount. Based on this, we may infer that if amoderately large number of elements ( (cid:39)
25 000 ) is used per processor after the first spin-up time ( t =
100 s ), changing the number of processorsmore often does not lead to significant changes in the observed CPU times. This might be due to the fact that Fluidity exhibits excellent scalabilityproperties (both strong and weak) when a small number of processors is used, and enough elements are distributed to each processor. DESKOS AND PIGGOTT
In this work, we presented the implementation and validation of a uRANS-based, mesh-adaptive ALM which is able to optimize the number ofelements/cells that are required to resolve the wake field within a full-size wind farm. Our ALM implementation is based on the original model of(34) 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 flowis resolved using a uRANS formulation of the governing equations and the k – ω SST turbulence closure model. The developed Fluidity-ALM modelwas validated against wind tunnel measurements before attempting full-scale wind farm simulations. For the wind farm simulations we selected theLillgrund offshore wind farm for which numerous previous studies and observed data are available. The effect of mesh-adaptivity on reducing thecomputational cost was also examined with particularly emphasis being given on processors’ load balancing, especially during the transient periodwhere the number of elements follows a linear trend.More specifically, for the developed model a detailed mesh convergence study was undertaken for both the fixed and the adaptive meshapproaches showing that for a given edge length h the two mesh strategies yield similar results for the integrated rotor characteristics such as thepower or thrust coefficient. Moreover, as the element’s edge length is reduced, the model predictions for the power coefficient C P converge to theexperimental reference values. Looking at the wake solutions on the other hand, mesh-adaptivity provides better estimates for the wake field, par-ticularly in regions where large shear (velocity gradients) exist. This is due to the ability of the mesh optimization algorithm to identify these regionsand allow elements with large aspect ratios to align with the underlying gradients in solution fields. A fixed mesh does not have such an ability andtherefore a better solution at these regions can be obtained only by further refining the mesh locally. The validation of the model using the two BlindTests (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 modelsand it is a crucial step towards accurately predicting offshore wind farm power output. The results from the comparison with BT2 give us confi-dence that, although only two turbines in-line are used, we are capable of obtaining high-fidelity wake solutions for the back row turbines, whenboth operate at peak conditions (optimal tip speed ratio).Looking at the Lillgrund offshore wind farm simulation results, the power predictions from each row (A–H) agree well with measured data from(24). The largest discrepancy between the model predictions and the measured data is observed on the second turbine in each row. This systemati-cally appears in all simulations, and is believed to be due to the torque chosen for the second turbine. Nevertheless, the two modelling approachesusing the fixed and the adaptive-mesh strategies yield almost identical results for the power coefficients. This re-affirms our conclusions from themesh convergence study which states that the turbine predictions are not dependent on the mesh strategy but rather the edge-length used at thelocation of the rotor. As far as the wake field is concerned, although there is lack of data for comparison, the accuracy of the power coefficients pre-dictions suggests that the model is able capture the magnitude of the wake-deficits along each row, at least in a time-averaged sense. In addition,mesh adaptivity was found to be a favourable choice leading to a reduction in the overall computational cost while maintaining the same accu-racy. Moreover, by changing the adaptation period T adapt , we were able to further reduce the overall computational cost without compromising themodel’s accuracy. An additional reduction of the computational cost was also observed when the adaptive simulations were initialised on a smallernumber of CPUs, with this number increased in response to the spin–up of the computational mesh. The adaptation period was considered to beconstant during each individual simulation. Again, this is not necessarily an optimal approach, since the frequency of the adaptations should alwaysbe based on the state of the dynamics. However, such an approach requires relatively complex error measure designs which are beyond the scopeof the present work. For an effective variable-frequency mesh-adaptive approach a “goal-based” a posteriori error measure would be good to con-sider. This seeks to generate the optimal mesh at every instance purely for maximising the accuracy in a user-defined “goal” (e.g. power or thrust),and thus resolution is not wasted at locations and times where it does not contribute to this goal.Regarding the limitations of the present approach, we should begin by discussing the validity of our choice to use a uniform inlet velocity with slipconditions on both the bottom and the top of the domain instead of a log-law profile and a rough wall model at the bottom. Such a simplification ofthe flow conditions were found to affect the ability of the model to predict the power output in the large scale simulations, particularly for the secondturbine of each row. Future formulations of the mesh–adaptive solver should consider a boundary layer, although extra care should be taken to avoidextensive refinement near the bottom of the domain via for example using a vertically variant mesh gradation technique. In addition, by adoptinga uRANS framework, many of the higher frequency interactions between the flow and the turbines were ignored and thus limiting the informationthat can be extracted from the simulations. To capture such effects, turbulence–resolving simulations need to be undertaken. Creech et al. (15)studied the Lillgrund offshore wind farm using the Fluidity solver (although using different discretisation options) with an LES formulation andmesh adaptivity and found good agreement with SCADA data. However, their study did not examine the efficiency of mesh-adaptivity, e.g. throughcomparisons with fixed-mesh simulations. Inherently, an LES solver which employs mesh-adaptivity techniques will impose additional constraintson the frequency of mesh adaptations and the elements’ aspect ratios, while a far smaller edge lengths and time step will be required. All theseadditional factors put the efficiency of coupling mesh adaptivity with an LES solver and conducting wind farm simulations in question and thereforeadditional studies will be needed to assess its applicability. On the other hand, coupling a uRANS solver with mesh adaptivity for wind energy ESKOS AND PIGGOTT problems may be seen as more appropriate, particularly for regional scale simulations that focus for example on the interaction of adjacent windfarms (44), or when the wind farm simulations are used within an adjoint–based optimisation algorithm (45).Based on the above, the mesh–adaptive uRANS framework was shown to be an appropriate and efficient tool for large scale wind farm sim-ulations. The observed computational efficiency of the present simulations suggests that the same approach may be applied to more complexconfigurations such as those of onshore wind turbines placed over an uneven terrain and that mesh–adaptivity could be used to solve for theensemble–averaged flow quantities (mean velocity, TKE) while requiring substantially fewer CPU hours. ACKNOWLEDGEMENTS
The authors would like to thank Dr James Percival for advice given on this work, and acknowledge funding from Imperial College London’s EnergyFutures Lab and the EPSRC (grant numbers EP/R007470/1 and EP/L000407/1). In addition, several insightful comments and meaningful sugges-tions provided by the anonymous reviewers helped us improve the quality of the presentation of the results. Support for parallel computations wasprovided by Imperial College’s Research Computing Service.
References [1] Ammara Idriss, Leclerc Christophe, Masson Christian. A Viscous Three-Dimensional Differential/Actuator-Disk Method for the AerodynamicAnalysis of Wind Farms.
Journal of Solar Energy Engineering.
Journal of Physics: Conference Series.
Physics ofFluids.
Wind Energy.
Renewable Energy.
Physics of Fluids.
Journal of Turbulence.
Wind Energy.
Journal of Fluids and Structures.
Energies.
Surveys in Geophysics. DESKOS AND PIGGOTT [17] Piggott Matthew D., Pain CC, Gorman GJ, Power PW, Goddard AJH. h, r, and hr adaptivity with applications in numerical ocean modelling.
OCEAN MODELLING.
International Journal for Numerical Methods in Fluids.
Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences.
AIAA Journal.
RenewableEnergy.
Renewable Energy.
Assessment of the Lillgrund wind farm: power performance, wake effects. report: LG Pilot Rep., Vattenfall Vindkraft AB, Stock-holm.; 2009. Available at https://corporate.vattenfall.se/globalassets/sverige/om-vattenfall/om-oss/var-verksamhet/vindkraft/lillgrund/assessment.pdf.[25] Menter F. R., Kuntz M., Langtry R.. Ten years of industrial experience with the sst turbulence model. In: :625-632; 2003.[26] Pain C.C., Piggott M. D., Goddard A.J.H., et al. Three-dimensional unstructured mesh ocean modelling.
Ocean Modelling.
Finite Element Methods for Flow Problems . John Wiley & Sons, Ltd; 2005.[28] Cotter Colin J., Ham David A., Pain Christopher C., Reich Sebastian. {LBB} stability of a mixed Galerkin finite element pair for fluid flowsimulations.
Journal of Computational Physics.
Journal of Computational and AppliedMathematics.
Computer Methods in Applied Mechanics and Engineering.
Computer Methods in Applied Mechanics and Engineering..
Computer Methods in AppliedMechanics and Engineering.
Ocean Modeling in an Eddying Regime , AGU 2008 (pp. 383–408).[34] Sørensen Jens, Shen Wen Zhong. Numerical Modeling of Wind Turbine Wakes.
Journal of Fluids Engineering.
Wind Energy.
RenewableEnergy.
Wind Energy.
Turbulent Flows . Cambridge University Press; 2000.[41] Eriksson O, Lindvall J, Breton S-P, Ivanell S. Wake downstream of the Lillgrund wind farm - A Comparison between LES using the actuator discmethod and a Wind farm Parametrization in WRF. In: ; 2015.[42] Jonkman J., Butterfield S., Musial W., Scott G..
Definition of a 5-MW Reference Wind Turbine for Offshore System Development. : NationalRenewable Energy laboratory; 2009.[43] Deskos G, Abolghasemi M. A, Piggott M. D. Wake predictions from two turbine models using mesh-optimisation techniques. In: Lewis A., ed.
Proceedings of the Twelfth European Wave and Tidal Energy Conference , EWTEC; 2017; University College Cork, Ireland. ISSN: 2309-1983.
ESKOS AND PIGGOTT [44] Hansen K S, Réthoré P-E, Palma J, et al. Simulation of wake effects between two wind farms. Journal of Physics: Conference Series.
Wind Energy Science.
How cite this article:
Deskos G., and M. D. Piggott (2018), Mesh-adaptive simulations of horizontal-axis turbine arrays using the actuator linemethod,
Wind Energy ,2018;00:1–6