Flow in a large wind field with multiple actuators in the presence of constant vorticity

Study of wind farms is an area of active research. Researchers have proposed simpliﬁed wind farm models which deﬁne the wake structure in a wind farm and how they aﬀect the performance of the wind turbines. Interestingly these models do not take into account an important aspect of ﬂuid ﬂow i.e. the ﬂuid-structure interaction (FSI) between the turbines and the wind, which has an important role to playl. This motivated researchers to implement numerical analysis tools to model the geometry of the wind turbines in computational ﬂuid dynamics (CFD) based models of a wind farm in order to better understand the wake structure and study the performance of the wind farms. But, modelling the complex geometry of the blades and the turbines makes these models computationally expensive. In this paper, we propose an FSI methdology which can simplify the blade resolving CFD models and eliminate the requirement for modelling these complex geometries during preliminary engineering phase. As an example, we present simulations of up to three back-to-back wind turbines and compare the results with those obtained from wind engineering software, FLORIS. The proposed methodology demonstrates how the approach can be used to develop a relatively less computationally expensive wind farm model. The approach formulated in this paper follows an intermediate way between the analytical wind farm models and CFD models by introducing modiﬁcations to one of the most basic wind farm models (Jensen’s model) and using it to develop a simpliﬁed CFD model using the Decomposed Immersed Interface Method strategy.


I. INTRODUCTION
The classical actuator disc theory 1 considers the horizontal axis wind turbine (HAWT) as a circular disc capable of extracting kinetic energy from the wind. This theory implements the law of conservation of momentum and Bernoulli's principle to compute the velocity in the wake of the wind turbine. However, this theory is valid for a single wind turbine subjected to uniform steady flow. Thus, over the years researchers have proposed different wind farm models which are based on this basic model.
One of the earliest work is Jensen's top hat model for wake expansion and velocity deficit based on the law of conservation of momentum 2 . With this model, Jensen demonstrated wake structure and energy outputs for 10 back-to-back wind turbines (or actuators) and for turbines placed in a circle. Subsequently, this model was implemented to calculate the output from wind farms 3 . Softwares like FLORIS 4 , WindPRO 5 , WAsP 6 , WindSIM 7 and a few others use Jensen's model extensively. As such, the current work is focussed primarily on Jensen's model.
Other models like Infinite Wind Farm Boundary Layer (IWFBL) model 8 or another IWFBL model based on mixing-length theory 9 exist, but some of their parameters are difficult to determine. More recently, a model for an array of wind turbines was proposed by splitting the wake into three regimes, the first is turbines' expoa) Electronic mail: basub@tcd.ie sure to multiple wake flow linking the flow deficit to wake expansion, the next corresponds to wake expansion in the vertical direction only after wakes have merged and finally, flow in balance with boundary layer for an infinitely large wind farm 10 . Researchers have even proposed a Gaussian wake model, the results of which were validated with Large Eddy Simulation (LES) and wind tunnel experiments of a miniature wind turbine as well as LES simulation of wake of a Vestas V80-2MW turbine 11 . Another work models wind farm based on super Gaussian shape function 12 . The Gauss curl hybrid model 13 incorporates the effects of wake steering in large wind farms. Apart from these, a few other models have also been proposed over time like a curled wake model for wind turbine under yaw condition 14 and two new models based on the Park model by including turbulence explicitly and another which considers blockage contribution from individual wind turbine to model the entire wind farm and hence, compute its energy output 15 .
With the advancement in computational technology, researchers started looking at wind farms numerically from late 1990s. Few such works on numerical analysis of wind turbine wakes using methods like Blade Element Momentum and actuator line have been carried out in the recent past [16][17][18] . However, an in-between approach between full CFD models and simple analytical models to suit the preliminary engineering needs is not something that has been explored in detail.
FSI problems demand extensive computational resources. This is due to the fact that geometric model of such problems using methods like finite element or finite volume require creating a complex mesh all around the body immersed in the fluid. The more irregu-lar the body, the more complex and finer the mesh gets, thereby demanding more computational cost. A full three-dimensional (3D) computational fluid dynamics (CFD) model of the entire wind farm which might consist of hundreds of turbines is computationally expensive both in terms of processing time as well as the memory and the number of processors required. During the preliminary engineering analysis and design stage when detailed site data is not available, investing in such computational resource is not worthwhile. Nonetheless one or more arrangement of the wind turbines in the wind farm is still required to be investigated in order to demonstrate the pros and cons of each arrangement from the perspective of power production. This is where the analytical two-dimensional (2D) models mentioned earlier become useful. The primary aim of these models is that during this preliminary analysis stage, they give an insight into the approximate power output of a wind farm without the need to develop the full 3D CFD model. However, an interesting fact about these models is that none of them include the effect of the presence of blades. Therefore, the current work aims at proposing an approach so that these models can be further refined to include the effect, the blades might have on the flow field, so that the preliminary estimate of power production from a wind farm using these models are more realistic. We intend to achieve the same by not going into complexities of the aero-elasticity of the blades, which of course provide useful information 19,20 , but rather trying to include its effect geometrically by considering a zone where the behaviour of wind is different compared to the rest of the domain.
It has been shown how complex geometries can be avoided and yet an acceptable solution can be arrived at by using reasonable computational resources if special FSI methodologies are implemented. 21 Immersed Boundary (IB) Method 22,23 or Immersed Interface Method (IIM) 24 are two such notable FSI methods. Over the years, researchers have put considerable effort to improve these methods in order to increase their numerical efficacy and reduce their dependence on extensive computational resources which has led to variants like Decomposed Immersed Interface Method (DIIM) 25 and applied to elliptic partial differential equations (PDEs) or Explicit Jump Immersed Interface Method (EJIIM) 26 for boundary value problems on irregular domain.
The main contribution of this paper is in proposing an FSI based simplified CFD approach suitable for modelling a large wind farm which can aid in providing a more realistic power output of a wind farm during preliminary engineering analysis stage. A wind farm model is coupled with DIIM to simulate the wake structure for back-to-back wind turbines in the presence of constant vorticity in 2D. The inclusion of constant vorticity is another advancement of the model over irrotational flows.

II. SETTING UP THE MODEL
Like most fluid flow problems, the atmospheric boundary layer (ABL) also has a mean and turbulent component. The work is based on the fact that the mean velocity in atmospheric boundary layer is predominantly responsible for the power production from wind turbines and the contribution from turbulence when compared to the mean velocity is not that much significant. However, turbulence does affect the fatigue life performance of the blades, which however is not the focus here. The Reynolds number in the present case is very large i.e. Re = vD/ν ∼ 10 2 × 10/10 −5 = 10 8 and as such viscosity has no role to play unless the scales of turbulence to be resolved are small enough to investigate the local effects on the blades for detailed analysis. As pointed out in literature review, one of the popular wind farm model to estimate power production is Jensen's model. However, it does not consider the effects of either vorticity or viscosity and is based on the mean velocity of ABL.
In this study, we therefore focus on the mean wind velocity field only and propose a model to include the effects of vorticity and compare the velocity fields of this proposed model with that of Jensen's model obtained using FLORIS. The velocity components parallel (axis X 1 ) and perpendicular (axis X 3 ) to the earth surface are denoted by V i ; i = 1, 3 respectively. It can be shown using the equation of continuity ∂V i /∂X i = 0 that for a 2D steady incompressible flow with constant vorticity, ω = ∂V 3 /∂X 1 − ∂V 1 /∂X 3 , the following holds true. 27 The two mean velocity components are now decoupled resulting in two Laplace equations. The simulation of the two components provide a complete description of the flow field in 2D.
B. Defining jumps using modifications to Jensen's model Actuators of equal radii, r 0 are dealt with in this paper. Here, r denotes the overall radius of the expanded wake at any typical section and r i,j denotes radius of the wake at actuator j generated by actuator i. For any typical actuator A n ; n = {0, 1, 2, ...}, f i and w i represent sections just in front and behind the actuator and V 1,i , V 1,i,d and V w,i represent velocity upstream, at and downstream of actuator respectively.V i is the weighted velocity of air entrained from the previous overlapping wakes before approaching Consider the case presented in Figure (1) where back to back actuators are shown with overlapping wake effects. Jensen used the law of conservation of momentum and proposed a linear wake expansion model for the velocity deficit assuming r = r 0 + αx, where, α(= 0.1) is an entrainment constant. Thus, when applied at w 0 and f 1 , one obtains where, U ∞ is the mean free stream velocity outside the expanded wake, and r 0,1 = r 0 + αx 0 . From the basic principles of the classical actuator disc theory, it is known that the velocity variation, −aU ∞ induced by the actuator disc on the streamwise component of the flow when superimposed onto the free stream velocity changes a is the axial induction factor. With these values of r 0,1 and V w,0 , V 1,1 comes out to be The optimum value of a = 1/3 enables an actuator to achieve its maximum power output (Betz limit). So, Jensen based the derivations for the value of a = 1/3 i.e. V w,0 = U ∞ /3. Experimental measurements by 31 when used to calibrate α produced a value of 0.070. However, Jensen considered α = 0.1 as the assumption V w,0 = U ∞ /3 is fairly uncertain. A similar approach by considering sections w 1 and f 2 would give an equation Jensen claimed thatV 2 ≈ U ∞ since the spacing between the actuators is usually large andV 2 is a(r 0 /αx 0 ) 2 U ∞ less than U ∞ ; a factor which can be neglected. Accordingly, for A n , the velocity at the face of the actuator can be defined by For the simulations presented, a is assumed to be the same across all the actuators. However, if the spacing reduces Jensen's model cannot be applied straightaway as the assumptionV 2 ≈ U ∞ is no longer valid and needs to be modified. For this case, it is reasonable to assumé V 2 ≈ V 1,1 just behind A 1 . From eq.(4), it can be observed that the free stream velocity field on the downwind actuator is 1 − 2ak times the free stream velocity field of upwind actuator. In a similar way, it can be considered that with respect to A n+1 onwards, the free stream velocity with respect to A n is changed from This implies From the perspective of DIIM, the wind turbines can be considered to be the regions of discontinuities in the domain of the wind bounded by the interface, Γ where the change in velocity occurs. Usually the root width of blades of such large turbines can range from 3m to 4m. Thus, a zone of width of the order of the root width of the blade is considered within which the free stream velocity changes as defined by eq.(8) and eq.(9). The boundary defined by a rectangle of height equal to the diameter of the actuator and width equal to the root width of the blade with its centre located at the hub height from the MSL defines Γ. It is assumed that there is no change in the derivative of the V 1 in a direction normal to Γ. The jumps for A n can hence be defined as 1. Case 1: Actuators are spaced far enough such that 2. Case 2: A n−1 and A n are spaced such thatV n ≈ V 1,n−1 In both the cases [V 1 ] t,n = [V 1 ] b,n = 0. As will be seen from the simulation results that this approach of using Jensen's model to evaluate the jumps show a good resemblance with the theoretical values. Iterative approaches require appropriate stopping criteria to both detect when the solution has been resolved to a required level of accuracy and also to prevent an onset of divergent behavior due to issues of stability or floating-point arithmetic. For many iterative methods, there are worst-case bounds on the number of iterations based on eigenvalue distribution in the case of a normal matrix (which includes the case of, e.g., symmetric matrices) and additionally non-orthogonality of the eigenvectors in the case of non-normal matrices.
The simulations use BiCGStab method combines two different approaches (BiCG and restarted GMRES) and is still amenable to some of this analysis. However, the DIIM approach introduces adjustments to the residual which can be interpreted as the action of a dynamically induced implicit "flexible" pre-conditioner. Thus, we do not have an implicit or explicit representation of the system matrix with which to apriori analyze expected convergence behavior or stability. We adopt instead a more heuristic approach by studying the residual to ascertain if sufficient convergence has been achieved and to detect any divergence due to instability.
Accordingly, the following criteria has been set for convergence.
where, ||r|| 2 and ||Φ|| 2 are the L2 norms of the residual and the solution vectors respectively. Some supplemental iterations (approx. 100 − 200) are performed, and it is observed that the residual has stabilized when these criteria are satisfied.
One additional issue one must be cognizant of is that these convergence criteria are related to the discretized problem of algebraic equations which approximate the physical model in question. Our convergence tolerances therefore also should not be more stringent than the er- ror induced by the process of modeling reality or the discretization process. Finite difference method is used for discretization. For the purpose of validation, we consider a commonly studied elliptic PDE problem, ∇ 2 φ = 0 as shown in Figure ( The analytical solution of this problem is given by eq.(10).
(10) For the purpose of numerical simulation, a mesh size of 40 × 40 is considered. Figure (3) shows a plot of the numerical solution of φ obtained using the iterative solver. Values at a few selected points are also shown in table I. It can be observed that the numerical and analytical results bear a good resemblance, thereby validating the performance of the solver.

III. RESULTS AND DISCUSSIONS
We now look into the results of the simulations carried out for back-to-back actuators. The challenge is that (2). The value of a = 1/3 is considered throughout. Though both V 1 and V 3 represent the full flow field, in this paper, we have simulated the horizontal velocity component, V 1 only, as the power output from an energy extracting actuator disc is dependent mainly on this component. Though the approach is valid for multiple actuators, in this paper, we present the study for upto three back-to-back actuator discs each of diameter, = 150 m with a hub height of 149.3 m. The spacing between adjacent actuators is considered to be 5 . A zone of 3 is considered in front of the first actuator to allow for the free flow of the wind towards the actuator. An additional zone of 5 is considered downstream of the last actuator following which it is assumed that the initial velocity profile is regained. The height of the domain is fixed by considering a zone of 5 above the topmost point of the disc. Accordingly, the domain size of 2700 m × 975 m is considered. Node to node spacing along the X 1 and X 3 axes are considered to be X 1 = 1 m and X 3 = 2 m respectively, which is a fine enough resolution for such a large domain.
It is to be noted that the velocity profile changes sign as X 3 → 0 m with a discontinuity at X 3 = 0 m. One way to circumvent this issue is to assume the values of V 1 near to the surface as 0. However, doing so would introduce a discontinuity in the boundary condition which in turn might incur instabilities in the numerical simulation. Instead, the model can be developed ignoring this small zone, the effect of which is practically negligible in such a large domain. Figure (

A. Single actuator
First the simulation is run for a single actuator located at X 1 = 448.3m. From eq.(8), the jumps conditions are The plot of the velocity field incorporating the effects of FSI is shown in Figure (5). Figure (6) shows that the shear profile of the wind is maintained for a reasonable distance. As expected, this profile starts getting affected approximately from X 1 = 325m and the effect becomes more prominent closer to the actuator. At X 1 = 449m, that is just after the jump at the interface occurs, drop in the velocity is observed as expected. The value of the velocity averaged over the height, after the drop is around 10m/s which is nearly equal to the theoretical value of   at X 1 = 453m is observed to be approximately 5.5m/s which is nearly equal to the average theoretical value of 4.64m/s. These results validate the approach of DIIM to model FSI in the context of the classical actuator disc theory. However, an interesting fact is observed from Figure (5) i.e. the presence of stationary points in the velocity field and the formation of vortex. Such stationary points are generated by the presence of vorticity. In fact, these points are going to affect the amount of energy extracted from the downstream actuators, a parameter not considered by existing wind farm models.

B. Two actuators
For the case with more than one actuator, the classical actuator disc theory cannot be applied straightaway as U ∞ = V 1 (X 3 ) is no longer valid. This is where Jensen's model can be put to use. At this stage, two parameters need to be validated first:-1. Average value of V 1 obtained numerically in the wake of A 0 compared with Jensen's model 2.V 1 → U ∞ as the distance from actuator A 0 increases For this purpose, numerical integration using Simpson's rule was carried out to determine the average value of the velocity within the expanded wake at X 1 = {650 m, 950 m, 1198 m, 1300 m, 1500 m, 1900 m}. The coordinates of the lower and upper limits of the expanded wake i.e. the limits within which the numerical integration is carried out to determine the average velocity is shown in table II. Jensen based his work on a freely expanding wake. So, for larger distances from actuator, X 3,l tends to become negative which is physically inconsistent. Hence, a correction to the average velocity obtained from Jensen's model is required. This is done by conserving the mass within the expanded wake.
where, ρ is the density of air, V J is the velocity for freely expanding wake, V JC is the corrected velocity. As V JC → U ∞ , this correction is no longer required as it can be assumed that the wake effect has died out and the original velocity has been regained. From the results presented in table II, it can be observed that (i) average value of V 1 in the wake of A 0 computed numerically shows a good resemblance with the values obtained from Jensen's model. (ii) as the distance from the actuator increases, the value of the average velocity increases such that at X 1 = 1500 m, the value (14.02 m/s) nearly reaches the average value of U ∞ = 13.40 m/s, thus, validating the fact that Jensen's model can be considered to assign the jumps for A 1 . We now look into two cases such that A 1 is placed at a distance where a) average value ofV 1 ≈ U ∞ ; b)average value ofV 1 ≈ V 1,1 (7) shows the velocity field for two actuators considering a = 1/3 for both A 0 and A 1 . A 0 is placed at X 1 = 1500.3 m considering the fact that V 1,1 ≈ U ∞ , an assumption consistent with Jensen's model for downstream actuators. Results show a nice convergence.
2. Case-2: A1 at X3 = 1198.3m Next we look into a case when A 1 lies within the wake of A 0 and try to examine the axial induction factor. It is observed that for a = 1/3 that the convergence of the numerical solution is only upto the order of 10 −4 and the residual does not stabilize. This means that it might not be possible to have a = 1/3 inside the wake and the power production will be even less. Accordingly, a needs to be modified to an acceptable value. Eq. (9) gives this value of modified axial induction factor as a(1 − 2ak) 2−1 = 0.27. Contour plot for this modified value of a is shown in Figure ( 8). As expected, the residual stabilizes with this value of a.

C. Three actuators
Now that we have established the applicability of both classical actuator disc theory and Jensen's model in the context of FSI, we now try to model a case with three back-to-back actuators in order to check the performance of DIIM for a greater number of actuator discs. As per eq.(9), the modified axial induction factors for this case are 0.27 for A 1 and 0.23 for A 2 . A stable solution is obtained as expected. Zones of low velocity are observed above A 1 and A 2 actuators signifying deficit of wind energy locally. As mentioned earlier, one of the well-recognised softwares to understand the preliminary wake structure is FLORIS. In this section, we compare the results of implementing DIIM with respect to a wind farm model modelled using FLORIS. Figure (10) shows the wake structure for each of the four cases dealt with in the previous sections. It can be clearly seen that none of the plots show the effect of vorticity. The reason being, though Jensen's model provides an outline of the expected wake structure, it does not consider the impact vorticity has on the velocity field in any way i.e. whether the flow is rotational or irrotational is not reflected in the velocity field. Further, the model itself does not take into account the effect the presence of blades will have on the velocity field. Our proposed approach, on the other hand is able to capture these details by using eq.(1) to incorporate vorticity coupled with the simplified FSI approach (DIIM) to incorporate the effect of the presence of the geometry of the blades for preliminary analysis. This can clearly be observed by comparing Figure ( We now compare the change in velocity profile for the case with one actuator. The change in velocity profile using DIIM both upstream and downstream of the actuator, A 0 is shown in Figure ( has on the upstream velocity at X 1 = 400 m is clearly not reflected when FSI effects are ignored in FLORIS model. Further, the velocity profile just upstream and downstream (i.e. at X 1 = 449 m and at X 1 = 453 m respectively) of the actuator are seen to be identical in FLORIS. However, because of the effects of vorticity and FSI, the magnitude of velocity and hence, the thrust at X 1 = 449 m is observed to be on the higher side compared to the case when these effects are ignored.

IV. CONCLUSION
In this work, we have proposed an FSI strategy using DIIM in the presence of vorticity and demonstrated how this can be implemented to simplify complex CFD wind farm models. Modifications to Jensen's model have been proposed to compute the appropriate jump conditions for the purpose of FSI analysis. We have theoretically computed values using actuator disc theory and modifications to Jensen's model and found that the results obtained numerically using FSI strategies are in close resemblance with the theoretical value. Presence of stationary points in the fluid are detected in the contour plots signifying the presence of vorticity and the impact of FSI. This fact also illustrates the importance of FSI in the wind farm models and how this interaction influences the generation of vortices. Though the simulations are carried out for a fixed value of axial induction factor, the approach can still be used for cases where the axial induction factor varies from turbine to turbine. We expect that other wind farm models like the Gaussian models can be similarly implemented alongwith DIIM.
The results obtained have been compared with FLORIS and it is found that the velocity deficit in the wake using FLORIS and without considering any vortic-ity are significantly different. The energy production of a wind farm which is a function of square of velocity, therefore is considerably affected, if FSI effects are ignored. Ideally in any atmospheric boundary layer, vortices are inherent. Hence, the proposed model can give us a better approximation to the power output of a wind farm during preliminary engineering analysis without spending excessive computational resources.
Though the application of DIIM has been been demonstrated for 2D elliptic PDEs, it can be extended to 3D elliptic PDEs as well. 25 In the context of wind farms, the approach presented in this paper when extended to 3D will then be able to give a more realistic overview of the flow field. Further, the rectangular internal boundary which is representative of the actuator disc is approximate and was considered to demonstrate that DIIM can be applied for the case of wind farms. The approach is valid for any shape of the internal boundary which is clearly reflected from the validation problem presented in this work. Thus, in future, a detailed study can be performed with a more accurate geometrical shape of the blades.