Dynamic flow analysis using an OpenFOAM based CFD tool : Validation of Turbulence Intensity in a testing site

The presenting paper investigates on the validation of the turbulence intensity (TI) modeled by a CFD tool. Six meteorological masts, equipped with cup anemometers, have been used for the purpose. Three different turbulence closure schemes, which are the SST k-omega and the k-epsilon in two different configurations, have been tested. The flow analysis shows a qualitative agreement between measurements and models, which are capable to simulate the turning of the wind towards South when it comes from SSE. Furthermore, the simulations predict a zone of high turbulence in the northern part of the site that is confirmed by the local measurements. The scores for TI have been quantified by considering the observed directional frequencies in the validation analysis. For the testing site, the SST k-omega scheme achieves the best performance when using the TI definition which is representative of the longitudinal fluctuations of the velocity vector, against the other one, which considers the fluctuation of the horizontal vector. Lastly, the model errors have been used to correct the simulated values using two approaches; the analysis shows that, for the presented case, these correction methods do not always improve the accuracy of the simulations.


Introduction
Because large shares of wind power projects are nowadays located in areas with complex terrain, Computational Fluid Dynamics (CFD) models are often used to model the dynamic atmospheric flow at specific points where wind measurements are not available.These models are tested to capture effects on wind conditions like flow separation, recirculation zones, where the planned wind turbines must be kept away from.Turbulence is an important parameter to evaluate the suitability of the designed Wind Turbine Generator (WTG) through the fatigue loads calculations: the wind speed at any one time oscillates around the average value of a longer period and these fluctuations are transferred by the force of the wind into fluctuating loads on the WTG.The turbulence can be quantified by the Turbulence Kinetic Energy (TKE), which is calculated from a closure scheme in the CFD model.In the wind power industry it is usually measured by cup anemometers and expressed in terms of Turbulence Intensity (TI).Our work is focused on the validation of two used approximations which relate TKE to TI.The wind flow analysis is modeled by a CFD tool which uses a modified OpenFOAM [1] to solve the 3D Reynolds Averaged Navier-Stokes combined with a chosen turbulence closure scheme.The testing site ITM Web of Conferences consists of a plateau bounded by cliffs and a range of hills located in a coastal area of Chile, where six meteorological masts (met masts), equipped with cup anemometers, are used for the validation.The wind data present a good availability, quantified in two full years of concurrent records that allows avoiding any seasonal bias in the analysis.
Different settings are tested for the purpose: two different translations of TKE into TI and two turbulence closure schemes, the SST k-omega and the k-epsilon.For the k-epsilon two configurations are validated by using different model constants.All the results are quantified and compared in terms of Root Mean Square Error and BIA scores, which are computed in a "per sector" analysis that takes into account the local observed directional frequencies.
Eventually, two different approaches are used to correct the results from the score analysis.The presented techniques use the model errors in order to improve the accuracy of the simulated values of TI.

Turbulence definitions
For a turbulent flow, the turbulence can be quantified by TKE per unit mass, given by: where u v w are the fluctuations of the three wind components at the time t, according to the Reynolds decomposition [2].In the wind power industry, turbulence is defined as turbulence intensity (TI) and derived from wind speed and standard deviation, which are usually measured by cup anemometers, according to the IEC (International Electrotechnical Commission) [3].The target of our work is the validation of two approaches which relate TKE, obtained from the turbulence closure scheme, to the measured TI.The first method is expressed by the following equation: The other testing approximation is to represent TI for a given direction (x, for instance) like: which is often found in literature [4].In Eq. ( 2) and Eq.(3) Ū and ūh are respectively the 3D and the horizontal wind speed averaged over an interval of 10 minutes, the standard deviation and u and v are obtained from: where u and v refer to the horizontal components of the wind, i is the index of the single measurement inside the sample interval, which has the frequency of 1 Hz and N is the total number of the measurements taken over 10 minutes.The approach described by Eq. ( 2) is obtained by considering the fluctuations of the horizontal velocity vector, whereas Eq. (3) considers instead the longitudinal fluctuations of the velocity vector.In both cases the flow inclination is neglected, which corresponds to measurement sites located in gentle terrain.In addition, cup anemometers only capture turbulent motions up to a certain frequency, while potentially CFD can model all the scales of turbulence.This can in general origin a source of bias when comparing measurements taken over 10 minutes with CFD simulations.
There are open questions about the ability of the cup anemometers to measure the turbulence and about which of the two presented approximations is more appropriate [5].In a way Eq.( 2) and Eq.(3) act as a lower and upper bound when comparing the turbulence with the measurements [6].In our case the output of the turbulence closure scheme is obtained from steady-state simulations and compared (using Eq. ( 2) and Eq. ( 3)) to the yearly mean value derived from two years of measurements.In general, the use of full years of TI data allows avoiding seasonal bias [7] in further calculations that are important to establish the suitability of the designed WTG.
TI is function of wind speed but it tends towards an asymptotic value for high wind speeds.In our case, for each met mast, the asymptotic value of TI has been calculated from an exponential decay function that is used to fit the measured TI and wind speed data [6] in a "per sector" analysis.This is a reasonable approximation for model validation when running neutral atmospheric conditions.Under these conditions, in fact, also the simulated TI does not change much with wind speed.
Lastly, it must be highlighted that the measurements represent a wider range of velocities and temperature stabilities than the simulations, which are calculated, as said, using neutral condition and also for one only velocity.Furthermore, in the steady case run flow features, which have boundary conditions that change with time, are not taken into the account.

OF-Wind tool
The wind flow analysis is based on the results obtained from OF-Wind tool, which uses the general Open Foam [1] model to solve the Reynolds Averaged Navier-Stokes (RANS) equations, combined with a turbulence closure scheme when modeling the flow over the terrain.The main features of the tool and the used setting are described in the following subsections.

Meshing
The meshing is derived using the SnappyHexMesh utility [8], which allows building a 3-D structure containing hexahedra and spit-hexahedra cells that adjust to an input geometry given in triangulated surface format (stl).Figure 1 illustrates the grid generated for the testing case and displayed on two vertical planes, which are respectively aligned in NNW-SSE and W-E directions.The same mesh k-e 0.09 1.44 1.92 1.00 1.11 k-e_ mod 0.03 1.21 1.92 1.0 1.30 configuration is used for all the simulations.The finest used grid contains 8.9 M of elements; above the surface layers the side length of the cell is 12.5 m.Close to the ground, 10 layers are used with an expansion ratio of 1.25 that allows setting the first node height to 0.3 m (see also Sect.3.4.2).

Solver
In all the simulations Open Foam is run using the steady-state SimpleFoam [9] solver.Segregated Solver is used for pressure-velocity scheme.The other chosen numerical schemes [10] are illustrated in Table 1.
The following assumptions are considered in the model formulation: • Coriolis terms neglected in the momentum equations • Air density constant at 1.225 (kg/m 3 ) • Incompressible atmospheric flow • Neutral stable atmospheric condition.
Convergence was checked by monitoring points which have been placed at the six met mast locations, at 59 m above the ground level (AGL).In details, in all the run it was assured that the local flow variables (TKE and wind speed) reached constant values at these points; this matched with the other condition of having the root mean square errors of all the residuals less than 10 −5 .

Turbulence schemes
Two different turbulence closure schemes are used in our case: SST k-omega, which takes into account the transport of the turbulent shear stress [11], and the k-epsilon model [12].In the k-epsilon model, TKE and the turbulence eddy dissipation (TED), (illustrated at the inlet by Eq. ( 6) in the next section), satisfy their respective conservation equations, which are expressed in terms of model constants.In Table 1 the values of these constants are shown in the two configurations (here called k-e and k-e_mod) which will be used in our simulations.For k-e the constants are obtained from the model of Apsley and Castro [13], while for k-e_mod the used configuration is derived from Bechmann and Sørensen (2009) mentioned in [14].

Inlet
The inlet wind profile is described by Eq. ( 5): First Symposium on OpenFOAM in Wind Energy where z is the altitude, z 0 the roughness length, k the Von Karman constant (0.41), u * the friction velocity derived from a prescribed reference velocity (set at 8m/s) at a reference height (set at 80 m AGL) and from an inlet roughness (set at 0.03 m).TKE and TED are expressed in the following way: where c is a model constant [12].Zero gradient equation is employed for the pressure field.

Bottom
No-slip condition is set for the wind field and zero gradient condition is used for the pressure.The wall shear stress is computed by extended wall functions based on empirical formula for sand-grain roughness.The used approach [15,16] consists in a modification of the standard smooth law of the wall which, for the first cell node, can be written as: where: are respectively the dimensionless z coordinate calculated at the height of the first cell node, z p and the dimensionless roughness height; ν is the kinematic viscosity and k s is the roughness height.In Eq. ( 7) B is the smooth log constant and the function B measures the departure of the wall velocity from smooth conditions.In OpenFOAM and in a equilibrium boundary layer, assuming a fully rough regime, Eq. ( 7) can be approximated by: where E is a smooth constant and C s is a roughness constant, which is set to ensure first order matching between the law of the wall and the inlet profile conditions.When the first cell node is much bigger than the roughness length, that is: a relationship between k s and z o is found and written as: This last equation will be always used in our simulations.It has to be highlighted that an inconsistency of the rough wall function is present in the formulation even when Eq. ( 10) is satisfied [15,16].In addition, the used approach for the rough wall function does not limit k s according to the first cell size and the accuracy of the results depends on the location of the first cell [17].In general it's difficult to establish how this inconsistency can affect the results.In [18], an agreement between rough law of the wall and measurements seems to be found when z p is approximately greater than 0.4k s .This constraint will be considered in all the work.In details, the maximum roughness length in the simulated domain is equal to 3 cm, so that z p will be set to 0.3 m in all the used grids.

Other boundary conditions
Zero gradient condition is used for all the variables at the outlet, except for the pressure, where a fixed value of zero is set.Constant shear stress [12] is used for the wind vector at the top of the domain, where slip condition is used for the pressure.The top of the computational domain is set to 6500 m, which corresponds approximately to 10 times the differences between the maximum and the minimum of altitude in the whole domain.Lastly, cylinder type domains are used in all the cases and symmetry wall is used at the patches parallel to the flow.

Validation and results in the testing case
The region of interest is a coastal area consisting of a plateau bounded by cliffs to the West and a range of hills to the East.The available wind data, used for the validation, were recorded at six met masts at 59 m AGL. Figure 2 shows the location of the six met masts (M1-M6) and the altitude of the terrain.Two main directional sectors of the wind rose, SSE and South, have been calculated for the testing case, running three different angles, using a step of 10 degrees.This technique allows increasing the accuracy of the results by reducing the influence of localized effects arising in a single angle, which are not necessary representative of the average flow in the wind sector. Figure 3 illustrates the mesh convergence analysis for TI computed from Eq. ( 2) using SST k-omega.In this analysis, the side length of the cell above the surface layers has been progressively double (from 12 to 100 meters), whereas the first node location does not change and it has a distance of 0.3 m from the wall for all the four grids.In addition, compared to the finest mesh, the number of surface layers has been increased for the coarser grids in order to achieve the same vertical resolution at the ground.The results, which are obtained for an inlet wind direction of 150 degrees and shown at the met mast locations (at 59 m AGL), offer an idea of the discretization error when comparing the values obtained from the finest grid to the coarser ones.Figure 4 illustrates the simulated TI at 59 m AGL obtained from SST k-omega for an inlet wind direction of 150 degrees; the horizontal wind vector is also displayed.The graph highlights that, when the wind comes from SSE at M5, the flow is shifted towards south in the centre (M3) and in the north part of the site (M4), where higher turbulence is predicted; the same behaviour is also observed from the output of the other two closure schemes (not shown).This is proved from the measured data time series, shown in Figure 5, where a qualitative agreement with the simulated wind direction and TI is achieved.

04002-p.6
First Symposium on OpenFOAM in Wind Energy  As results, the derived wind rose (Fig. 6) for M5 exhibits higher frequency at SSE compared to the other two (M3 and M4), confirming the trend.Furthermore, the simulations of the flow at 190 degrees (not shown) exhibit a turning of the wind towards SSW at M4 in agreement with the higher frequency measured for this directional sector (Fig. 6c).
In order to quantify the results, two scores have been used [19]: the weighted BIAS, here called WBIAS, and the weighted Root Mean Square Error (WRMS).Both scores have been computed from First Symposium on OpenFOAM in Wind Energy the following equations, which take into account the observed directional frequencies: where TI is again the turbulence intensity, j is the index of the local simulated direction, l is the total number of directions (l = 6 in our case from 140 to 190 degrees), i is the met mast index, n is the total number of met masts (n = 6), S and M respectively refer to simulated and measured values and f is the number of the measured data.WBIAS offers a mean assessment about the overestimation (if positive sign) or underestimation (if negative sign) of the modeled values, while WRMSE quantifies the error committed by the model.The results for both scores are illustrated in Table 2.For the testing case ke mod model is the less accurate one.In details, this scheme overestimates with both definitions of TI and exhibits the largest error when using Eq. ( 2).The best score is achieved from SST k-omega with Eq. (3) (WRMSE equal to 2.3%).Interestingly, both SST-k-omega and k-e underestimate with Eq. ( 3) and overestimate with Eq. ( 2), proving that the two tested approximations bound the measured TI.
Figure 7 illustrates wbias obtained at the six met mats, proving that k-e_mod always overestimates with Eq. (2).

BIAS and Linear correction methods
In this section the attempt is to use the model errors in order to improve the accuracy of the simulation results.In this analysis we use the score obtained at met mast M5 to compute the final results.Two approaches are presented for the purpose.The first one, BIAS method, is based on Eq. ( 12): 04002-p.9 where again i is the met mast index and the other symbols were defined in section 4. The new simulated TI at the target met mast o, using score analysis from the i-th met mast, is calculated in the following way: where r refers to wbias method and TI s o indicates the original simulated value computed as: The second approach, here named Linear, is based on a simply linear function between simulated and measured couples of data.In our case the used equation is expressed by: Eq. ( 17) is then used to compute the new corrected TI by the following expression: where R refers to Linear method, and the other symbols were defined previously.Finally, we introduce the last score, the Absolute Error, AE [18], which is computed in order to compare three results at met mast o: the two modified values, calculated from the Eq. ( 15) and Eq. ( 18), with the original (not modified) ones.Indicating with S the original uncorrected value, the AE is then computed by Eq. ( 19): where: is the measured TI calculated as weighted average by directional frequency.
The results are shown in Table 3 for SST k-omega model when using TI definition from Eq. ( 3) and validation analysis from met mast M5.Except for o = M3 and, of course, o = M5 both methods do not improve the model score in the other cases.BIA method seems to work better compared to the linear one at met masts M1 and M4, whereas few differences between the two approaches are seen at the other three met masts.Finally, it must to be highlighted that also a method of best fit, which is carried out by considering five couple of simulated and measured data, does not improve the final accuracy at the target met mast.This is the case (not shown) when the parameters obtained from the linear regression analysis are used to correct the TI simulated by SST k-omega scheme at M4 and derived from Eq. (3).

04002-p.10
First Symposium on OpenFOAM in Wind Energy

Conclusions
A CFD tool has been used to simulate the dynamic of the flow over a hilly terrain.The simulations, which have been carried out for two directional sectors, SSE and South, using a step of 10 degrees, have shown that the model is capable to qualitatively reproduce a characteristic feature of the wind flow at the site.This has been proved from the experimental data, which show higher turbulence for the northern met masts and also the shift of the wind direction to South when the wind comes from SSE at the met mast M5.A quantitative analysis for TI has been carried out using two different definitions of this parameter and three turbulence models schemes.The SST k-omega scheme achieves the best performance when using the turbulence intensity definition which is representative of the longitudinal fluctuations of the velocity vector, against the second one which considers instead the fluctuation of the horizontal vector.The k-e closure scheme, using modified model constants, exhibits the largest error, whereas the standard model constants for k-e seem to offer more accurate results and comparable to the ones obtained from SST k-omega.For these last two schemes, in fact, BIAS analysis has shown that the two used approximation of TI can be considered as a lower and upper bound for the measurements.This behaviour is not seen for k-e_mod that overestimates with both definitions of TI.Lastly, two correction methods have been tested in order to reduce the model errors from the validation analysis.For the presented tests the first method, which is based on the BIAS score, exhibits better performance compared to the second one, where a linear function is applied.Anyway, both methods do not always improve the model score when using the validation analysis for one met mast.

Figure 1 .
Figure 1.Cells generated from snappymesh tool for the testing case and displayed on two vertical planes.Altitude in meters is also shown.

Figure 3 .
Figure 3. Convergence grid analysis for an inlet angle of 150 degrees using SST k-omega closure scheme: results of TI at the met mast locations (at 59 m AGL) obtained using different side lengths of the cell.

Figure 4 .
Figure 4. Turbulence Intensity at 59 m AGL (obtained from Eq. (2) and horizontal wind vector obtained from SST k-omega for an inlet wind direction of 150 degrees.

Figure 7 .
Figure 7. Wbias for TI obtained at the six met masts.

Table 1 .
Numerical schemes used in the solver.

Table 2 .
Model constants used in the k-epsilon scheme.

Table 3 .
Results of WBIAS and WRMSE for TI.

Table 4 .
Absolute Error (in %) for TI using met mast i = M5 computed from: uncorrected (AE s ), BIAS method (AE r ) and Linear method (AE R ).