Comparison between 2.5-D and 3-D realistic models for wind field adjustment

Comparison between 2.5-D and 3-D realistic models for wind field adjustment L. Ferraguta,∗, R. Montenegrob, G. Monterob, E. Rodr´ ıguezb, M.I. Asensioa, J.M. Escobarb aUniversity of Salamanca, University Institute on Fundamental Physics and Mathematics, Spain. bUniversity of Las Palmas de Gran Canaria, University Institute of Intelligent Systems and Numerical Applications in Engineering, Spain. Abstract In previous works, many authors have widely used mass consistent models for wind field simulation by the finite element method. On one hand, we have developed a 3-D mass consistent model by using tetrahedral meshes which are simultaneously adapted to complex orography and to terrain roughness length. In addition, we have included a local refinement strategy around several measurement or control points, significant contours, as for example shorelines, or numerical solution singularities. On the other hand, we have developed a 2.5D model for simulating the wind velocity in a 3-D domain in terms of the terrain elevation, the surface temperature and the meteorological wind, which is consider as an averaged wind on vertical boundaries. Using the meteorological wind as datum, the 2.5-D model provides a 3-D local wind modified by topography and thermal gradients on the surface by solving only a 2-D optimal control problem where the boundary condition is the control. In this case, the finite element discretization consists on a triangular mesh adapted to the terrain topography and roughness length. In both models, the wind field adjusts to several wind speed measurements at several points in the 3-D domain and eventually to an average wind flux on the boundary. In this paper we introduce several advances in the 2.5-D and 3-D wind models and we compare their results on a region located in the Province of Lugo (Spain) with realistic data that have been provided by the company Desarrollos E´ olicos S.A. (DESA). In order to obtain the best adjustment of models results to the measurements, the main parameters governing the models are estimated by using genetic algorithms with a parallel implementation. Key words: Wind field modelling, mass consistent models, parameter estimation, genetic algorithms, adaptive mesh refinement, finite element method. PACS: Preprint submitted to Elsevier 19 April 2010 Manuscript Click here to view linked References 1 Introduction The society has been becoming aware of environmental problems (climate change) and nowadays it appreciates the use of renewable energies. Along last years the use of wind power for producing electric energy has augmented considerably. So companies of this sector are requesting more and more sophisticated tools that allow them to face the competitive and demanding market. Wind models are interesting tools to the study of several problems related to the atmosphere, such as, the effect of wind on structures, pollutant transport [26], fire spreading [19], wind farm location, etc. Diagnostic models are not used to make forecasts through integrating conservative relations [30]. Therefore they are also called kinetic models [12]. These models generate wind fields that satisfy some physical conditions. If mass conservation law is the only imposed restriction, we are defining a mass consistent model [13,35]. The relative simplicity of diagnostic models makes them attractive from the practical point of view, since they do not require many input data and may be easily used. Pennel [29] checked that, in some cases, improved mass consistent models such as NOABL and COMPLEX obtained better results than other dynamic models which are more complex and expensive. However, we have to take into account that a simple mass consistent models neither consider thermal effects nor those due to pressure gradients. As a consequence, problems like sea breezes can not be simulated with these models unless such effects are incorporated into the initial wind data using observations in selected locations [11,28]. So, 3-D mass consistent models have been designed for simulating the effects of the orography on the steady average wind (i.e., average wind in intervals from 10 minutes to 1 hour) by imposing the incompressibility condition for the air. There exists a wide range of diagnostic models which have been used by scientists in problems of meteorology and air pollution. The primitive 2-D mass consistent models did not consider the terrain orography and the vertical profile of the wind. They built an horizontal interpolated wind field taking into account only the distance from the mesh nodes to the measurement stations and then they solved the two-dimensional elliptic problem arising from the discretization in a plane; see [37] for a 2-D adaptive finite element model with a mixed formulation. Nowadays, in problems defined over complex terrain, it is possible to use high quality adaptive 3-D meshes of the studied domains. However, most of the existing models use to work with uniform meshes. This strategy is impracticable for problems with complex terrain since the size of the elements must be very small in order to capture the digital information of the terrain elevation. Moreover, in this case there would be regions where such small element size would ∗Corresponding author. Email address: [email protected](L. Ferragut). URL: http://web.usal.es/ferragut(L. Ferragut). 2 not be necessary. This would finally lead us to larger linear systems of equations and higher computational cost for solving them. Our first 3-D models (see [23,24]) had some of these limitations since the meshes used for defining the terrain surface were uniform. In [26,27], we presented a new 3-D finite element model that uses adaptive unstructured meshes of tetrahedrons with elements of small size where it is necessary but maintaining greater elements where such level of discretization is not required. The resulting 3-D mesh were constructed by using a refinement/derefinement process to adapt a 2-D mesh to the terrain surface, a vertical spacing function to locate nodes in the air and a Delaunay triangulation algorithm [6]. This adaptive mesh generation technique was introduced in [20,21,25]. The generated mesh has a higher density of nodes near the terrain surface, where we need more precision, or close to significant contours that has an important role in the simulation, like shorelines or roughness length contours [7,8]. In a postprocess, the 3-D mesh is smoothed and, if necessary, untangled by using the algorithm in [7] in order to improve the mesh quality. In addition, a local refinement procedure was proposed for improving the numerical solution [10,16]. Finally, although mass consistent models are widely used, they are often criticized because their results strongly depend on some governing parameters. These parameters are generally approximated using empirical criteria. Our 3-D wind model includes a tool for the parameter estimation based on genetic algorithms [17,27,33]. We have also introduced in the 3-D model a characterization of the atmospheric stability that is carried out by means of the experimental measures of the intensities of turbulence. In addition, since several measures are often available at the same vertical line, we have constructed a least square adjustment of such measures for developing a vertical profile of wind velocities from an optimum friction velocity. The main goal of this paper is to compare and show the ability, as well as establish its limits, of two models providing a wind field with a minimum of experimental wind data in a real and practical case, more precisely in the case of wind farms where the user needs to have a short term prediction of the wind from the last meteorological measures and predictions. The space scale of the meteorological prediction is usually too large to cover the engineering requirements in a wind farm. These models provide a solution that can satisfy these requirements and represent an improvement from the existing mass consistent codes. We compare results of the 3-D mass consistent model with a 2.5-D vertical diffusion wind field model. If the significant phenomena that we want to simulate occurs in a zone where the horizontal dimensions are much larger than the vertical one, then an asymptotic approximation of the primitive Navier-Stokes equations can be obtained as in the model developed in [1]. The most relevant aspect of this asymptotic approach is that provides a three–dimensional velocity wind field that verifies the incompressibility condition in the air layer, governed by a two–dimensional 3 equation, so that it can be coupled with temperature surface distribution in order to take into account thermal effects as sea breezes. In addition, the terrain elevation is also considered by the model. The validity of this 2.5-D model has the following limits: The nonlinear terms are neglected and we assume that the air temperature linearly decreases with the height. On the contrary the model takes into account buoyancy forces, slope effects and mass conservation. The 2.5-D wind model presented in this article is an adaptation of the wind model proposed in [1], such that the data can be given on several points located in the domain. These data should be obtained from experimental measurements or meteorological predictions. In order to check the accuracy of the results and the efficiency of the two models with realistic data, the company Desarrollos E´ olicos S.A. (DESA) has provided us with technical support about digital terrain elevation maps related to orography and roughness length, as well as measurements of wind and turbulence intensity in several anemometers located in Lugo (Spain). The outline of this paper is as follows. In section 2, we describe the general notation common to both models. In section 3, we present the 2.5-D vertical diffusion wind field model, including the asymptotic equations and the adjustment of point data by solving an optimal control problem. In section 4, we summarize the 3-D mass consistent model and we introduce a technique for inserting new information about wind measures at different heights and turbulence intensities. The 2-D adaptive discretization of the terrain surface and the 3-D mesh generation procedure is presented in section 5, including an adaptive strategy to capture the orography and roughness information simultaneously and additional local refinements in different regions of the terrain. Several ideas about the parameter estimation of both models, with a parallel implementation of genetic algorithms, are summarized in section 6. A comparison between the results of the 2.5-D and 3-D models, in a realistic case for an episode along a day, is presented in section 7. Finally, we summarize the conclusions of this work and the topics that need further research. 2 Notation Let us consider the three–dimensional domain Ω={(x, z): x∈ω, h(x)< z < δ} representing the air layer of the studied region. In the 2.5-D model, we assume that the height δis small in front of the width and that the height of the terrain surface at point x,h(x), is smaller than δ. We note that previous assumption is not necessary in the 3-D model. We decompose the boundary of Ωinto ∂Ω = S∪A∪L, where S={(x, z) : x∈ ω, z =h(x)}is the terrain surface, A={(x, z) : x∈ω, z =δ}is the air upper 4 horizontal boundary and L={(x, z) : x∈∂ω, h(x)< z < δ}is the air lateral vertical boundary. For the 2.5-D model, let ω⊂ ℜ2be a two–dimensional normalized bounded domain, representing the projection of the three–dimensional terrain surface S. We denote by (x, z)any point of the three–dimensional domain Ω, and by (x)any point of the bi–dimensional domain ω. Then we denote by an index xzthe three– dimensional operators, and by an index xthe bi-dimensional operators. We use small letters for the two–dimensional problem, and capital letters for the three–dimensional problem. U= (U1, U2, U3)denotes the air velocity field. For the 2.5-D model, we distinguish the vertical velocity from the horizontal one denoting W=U3,V= (U1, U2).P is the potential, Tis the temperature, τthe time and vmis the meteorological wind. Depending on the context the symbol Urepresents instantaneous velocity or mean velocity. Finally Nis the inner unit normal vector field to ∂Ω, and n= (n1, n2)is the inner unit normal vector field to ∂ω. 3 Vertical Diffusion Wind Field Model in 2.5-D In this section we present the 2.5-D wind model which includestemperature effects. An asymptotic analysis gives a three dimensional convective model governed by a two dimensional equation. This model adjusts a three dimensional velocity wind field in a thin layer under the influence of the orography and temperature distribution, so that it can be coupled with two dimensional fire spread simulation models [1]. 3.1 Asymptotic equations The air velocity U= (U1, U2, U3)and the potential Psatisfy the Navier–Stokes equations. On one hand, the momentum equation reads ∂τU+U· ∇xzU−1 Re∆xzU+∇xzP=λ′Te3in Ω (1) where Re is the Reynolds number, λ′is an expansion factor and e3= (0,0,1). 5 On the other hand, the air compressibility is also neglected, so that ∇xz·U= 0 in Ω (2) Boundary conditions are U·N= 0,∂U ∂Ntan =ζUon S(3) U3= 0, ∂zU1=∂zU2= 0,on A(4) UL= (vm,0) on L(5) where ζis the friction coefficient and the subscript tan denotes the tangential component, vmis the meteorological wind, that we assumed to be known, horizontal, non depending on zand with a null total flux through the lateral boundary L, that is, ∂zvm= 0,Z∂ω(δ−h)vm·nds = 0 (6) We complete these equations with the initial condition U|τ=0 =U0(7) where U0is the initial velocity, that we assume to be known. Equations (1) to (7) are well posed. We distinguish the vertical velocity from the horizontal one denoting W=U3, V= (U1, U2), and we define the horizontal flux at a point x∈ωand time τby V(τ, x) = Zδ h(x) V(τ, x, z)dz (8) The incompressibility and the fact that the air does neither cross Snor A, that is, U·N= 0, involve that the horizontal flux is also incompressible, then ∇x·V= 0 in ω(9) Using the fact that thickness δof the considered air layer is small compared with its width and assuming that the wind is not too strong, more precisely δ′2Re ≪1 where δ′=δ width of the layer , then preserving only the dominant terms and re-scaling P, equations (1) and (2) give in Ω −∂2 zzV+∇xP= 0 (10) ∂zP=λT (11) ∇x·V+∂zW= 0 (12) 6 where λ=λ′Re. Reynolds number can be absorbed and does not appear in (10) by choosing properly the scaling of the pressure P. Conditions (3), (4) and (5) particularly give (V, W )·N= 0, ∂zV=ζV,on S(13) W= 0, ∂zV= 0 on A(14) V·n= (δ−h)vm·non ∂ω (15) See [4] for a complete derivation of this vertical diffusion model in a similar case. Equations (10) to (15) are well posed: For given Tand vm, there exists a unique solution (V, W, P )(up to an additive constant for P). For more details about this convection asymptotic model see [1]. Equations (10) and (11), together with conditions (13) and (14), provide V(x, z) = m(x, z)∇xp(x) + n(x, z)∇xˆ t(x)(16) where m(x, z) = 1 2z2−δz −1 2h2(x) + (δ+ξ)h(x)−ξδ n(x, z) = −1 24z4+1 6δz3−1 3δ3z+1 24h4(x)−1 6h3(x)(δ+ξ) +1 2ξδh2(x) + 1 3δ3h(x)−1 3ξδ3 being ξ=1 ζthe inverse of the friction coefficient ζand ˆ ta re-scaled temperature related to the terrain surface temperature t=t(x)by ˆ t(x) = λt(x) δ−h(x). We are assuming that the air temperature linearly decreases with the height, T(τ, x, z) = t(τ, x)δ−z δ−h(x). The potential p(x)satisfies the following boundary problem −∇x(a∇xp) = ∇x(b∇xˆ t) in ω(17) a∂p ∂n=−b∂ˆ t ∂n+ (δ−h)vm·non ∂ω (18) where a=a(x) = 1 3(δ−h(x))2(3ξ+δ−h(x)) and b=b(x) = 1 30(δ−h(x))22δ2(2δ+5ξ)−2δ(δ−5ξ)h(x)−(3δ+5ξ)h2(x)+h3(x) Finally, the vertical velocity Wcan be obtained from equation (12). 7 3.2 Adjustment of point data by solving an optimal control problem Usually the wind flux is not known on the boundary, then problem (17)-(18) can not be solved directly because in general we do not know vmon the boundary. Alternately we know the direction and the intensity of the wind at the points where the meteorological stations are placed. In this subsection we reformulate the former problem so that the required data values are the wind measures at some given points. In order to simplify notation, as in the following we are related only to the bidimensional problem, we avoid the index xof the differential operators. Let v= (δ−h)vm·n, the wind flux on the boundary ∂ω,v∈ V the space of functions verifying R∂ω v= 0. We are going to formulate the former problem as an optimal control problem. Given Nexperimental measurements of the horizontal wind velocity Vi, i = 1, ..., N, at Ngiven points Pi= (xi, zi), i = 1, ...N, we search for the value of v∈ V such that the value V(xi, zi)given by the expression in (16) are as close as possible to the experimental value Vi. This is an optimal control problem where •v∈ V is the control. •Problem (17),(18) are the state equations. •We will choose as cost function J(v) = 1 2 N X i=1 kV(xi, zi)−Vik2+α 2Z∂ω v2 that is J(v) = 1 2 N X i=1     m(xi, zi)∇p(xi) + n(xi, zi)∇ˆ t(xi)−Vi    2 +α 2Z∂ω v2(19) In practice, we can use instead of (19) the following regularized functional J(v) = 1 2 N X i=1 Zωρǫ,i(x)    m(x, zi)∇p(x) + n(x, zi)∇ˆ t(x)−Vi    2 +α 2Z∂ω v2(20) where ρǫ,i is a suitable smoothing function given for example by ρǫ,i(x) = 1 ǫ2ρ(x−xi ǫ) 8 ρ(x) = Me−1 1−||x||2for ||x|| <1 0 for ||x|| ≥ 1 for a small ǫand Msuch that Rρǫ,i(x)dx= 1 The optimal control problem to be solved is posed as follow: Find u∈ V such that J(u) = inf v∈V J(v)(21) The solution uof the optimal control problem (21) is characterized by J′(u) = 0. Using the general optimal control theory [15], and introducing the adjoin state, then the problem (21) is characterized by the following three equations relating p , q and u: •Zω a∇p(u)· ∇ϕ+1 αZ∂ω qϕ =−Zω b∇ˆ t· ∇ϕ∀ϕ(22) • Zωa∇q(u)· ∇ψ − N X i=1 Zωρǫ,i(m∇p(u) + n∇ˆ t−Vi)m· ∇ψ= 0 ∀ψ(23) • u=−1 αqon ∂ω (24) There exist a unique solution of the problem (9). Moreover, problem (20),(21), has a unique solution [15]. 3.3 A modified optimal control problem If a good estimation of the wind flux v≈u∗is known on the boundary then we can modify the cost function in (20) as follows J(v) = 1 2X iZω ρǫ,i(x)    m(x, zi)∇p(x) + n(x, zi)∇q(x)−Vi    2 +α 2Z∂ω(v−u∗)2(25) Note that we cannot impose at the same time the value of the wind flux on the boundary and the value of the solution at several givenpoints, as once the wind flux 9 We note that, in the case of the 2.5-D model, the problem is solved by only considering an adaptive 2-D triangulation of the rectangular region which is studied. 6 Parameters Estimation with Genetic Algorithms Genetic algorithms (GAs) are optimisation tools based on the natural evolution mechanism [2,17,36]. They produce successive trials that have an increasing probability to obtain a global optimum. This work is based on the model developed by Levine [14]. It is a standard genetic algorithm code (pgapack library), with string real coding. In the numerical experiments with the 2.5-D model, we look for optimal values of the quadratic adjustment parameters a0, a1and a2of the friction coefficient ζin terms of the roughness of the terrain, i.e., ζ=a0+a1z0+a2z2 0. We search for the optimum of the linear parameter a0in [1,10], the first order parameter a1in [0,5] and the second order parameter a2in [−0.05,0.05]. In the numerical experiments with the 3-D model, we look for optimal values of α,ε,γand γ′. Specifically, the so called stability parameter α=α1 α2determines the rate between horizontal and vertical wind adjustment. For α >> 1flow adjustment in the vertical direction predominates, while for α << 1flow adjustment occurs primarily in the horizontal plane (see equation (30)). Thus, the selection of αallows the air to go over a terrain barrier or around it. We search the optimum in [10−2,102]. The second parameter to be estimated is the weighting coefficient ε (0 ≤ε≤1) involved in the horizontal interpolation of wind measurements (see equation (36)). For ε→1, the importance of the horizontal distance from each point to the measurement stations is greater, while ε→0signifies more importance of the height difference between each point and the measurement stations. The parameter γis related to the height of the planetary boundary layer (see equation (41)). There exist different versions of where to search for this parameter. The interval [0.15,0.4] considered in our simulations includes all the proposed search spaces. Finally, the parameter γ′appears in the computation of the mixing height for stable atmosphere (see equation (42)). Several authors have proposed that the value of γ′should be searched in the surroundings of 0.4. More details and references about the discussion of these parameters can be found in [27]. We propose to minimise the following fitness function which is defined as the average relative error of the wind velocities given by the model with respect to the measures at the reference stations, F=1 Nr Nr X i=1 ||Vi−V(xi, zi)|| ||Vi|| (53) 16 where V(xi, zi)is the horizontal wind velocity obtained by the model at the location of station i,Viis the horizontal measured wind and Nris the number of reference stations. Note that in the case of the 2.5-D model, F=F(a0, a1, a2)and in the case of the 3-D model, F=F(α, ε, γ, γ′). 7 Wind Simulations in a Realistic Episode In order to compare the results of the two wind field models, we have considered a simulation using realistic wind data that have been supplied by DESA in several measurement points for an episode along the March 21, 2003, see [22]. The first step is to discretize the studied domain, the second is to estimate the main parameters of the model and, then, apply the wind model using the estimated values. Next, the wind velocity is checked in the control points. We present several applications to show the improvements carried out in our wind model. All experiment were run on a XEON precision 530, except the parameter estimation problem which was solved using a cluster of PCs. 7.1 Surface mesh adaption to orography and roughness The studied three-dimensional domain Ωis located in a region of Lugo, Spain, at 43Nof latitude and it is defined by four points of UTM coordinates A(609980, 4799020),B(626000,4799020),C(626000,4813040) and D(609980,4813040), respectively. The upper boundary Aof Ωhas been taken at a height δ= 4000 m in the 3-D model and δ= 1080 min the 2.5-D model. A digital elevation map was provided by DESA on a quadrilateral grid of element size 20 ×20 m. The X axis corresponds to East direction and the Yone to North. Thus, we are working with a region of 16020 ×14020 m. The minimum and maximum terrain heights are 420 mand 1020 m, respectively. Figure 1 represents a color map of the heights of the terrain. The measurement stations and the control points have been approximately plotted, such that from North to South we can find E243, E208, E212, E242, E206 and E283. Tables 3 and 4 contain their coordinates, respectively. The height of all measurements and control points to the ground are indicated in Table 5 in parentheses. All these points are located close to the top of the hills. Roughness is an essential factor on the atmospheric stratification, and therefore, on the characteristics of the resulting wind profile. Figure 2 shows the roughness length of the terrain which were supplied by DESA. We remark that some stations and control points are closed to contours of the roughness. In this case, the roughness length values are 0.03 m,0.05 m,0.08 m,0.3mand 0.8m. 17 Station UTM-E UTM-N Height E206 615396 4805218 924.8 E208 616917 4807256 945.0 E212 617423 4806382 895.0 Table 3 Coordinates in mof the measurement stations. Control point UTM-E UTM-N Height E242 618290 4806136 873.2 E243 616629 4808235 947.0 E283 617473 4804111 849.0 Table 4 Coordinates in mof the control points. Fig. 1. Elevation map (m) of the studied region in Lugo. From North to South, we can see the stations or control points E243, E208, E212, E242, E206 and E283. Starting from a regular mesh of the rectangular region with element size of 1× 1km approximately, five global refinements are carried out using 4-T Rivara’s algorithm [32]. With this number of refinement steps, we obtain a mesh with an element size about 31 m. In order to improvethe discretization near the stations and control points, five additional local refinements are applied inside six circles with centre at the stations and control points, respectively, and diameter 200 m. This 18 Fig. 2. Roughness length map (m) of the studied region in Lugo with the station and control points. produces a local element size about 1m. Once we have interpolated the height and the roughness length in the nodes of these refined two-dimensional mesh, we use the derefinement algorithm [9,31] described in section 5.1 with εh= 10mand εr= 0.01m, keeping in any case the nodes located inside the six circles. In Figure 3 we can see the resulting triangulation of the terrain surface. The corresponding threedimensional mesh, see Figure 4, contains 102662 nodes and 515812 tetrahedrons. 7.2 Parameter estimation along a day We have taken into account Table 2 for determining the stability class from the available turbulence intensity values. For the studied day, we have obtained neutral conditions. So, for the 3-D model, we must estimate the stability parameter α, the weighting parameter εrelated to the horizontal interpolation of wind velocities and the parameter γinvolved in the computation of the planetary boundary layer; see, e.g., [27]. The estimation has been carried out each hour (24 computations). We have applied genetic algorithms to solve these parameter estimation problems, where the fitness functions (see equation (53)) are defined in terms of the relative velocity errors obtained by the model at the measurement stations. It is evident that, in order to avoid spurious solutions, more than 20 repetitions for parameter setting of each hour should be done. This fact would obviously imply an important increas19 Fig. 3. Triangulation of the terrain simultaneously adapted to orography and roughness corresponding to the studied region in Lugo. Fig. 4. Adaptive 3-D mesh corresponding to the studied region in Lugo. 20 ing on the computational cost of the parameter estimation process, even more if we take into account that each evaluation of the fitness function supposes the resolution of a finite element problem (in this case, about one hundred thousand unknowns). Nevertheless, from our previous experience in this kind of wind simulation problems, we have observed that just one computation is enough for reaching a good solution. In Figure 5 we can see the evolution of the values of the three parameters of the 3-D model along the episode. The values of εare practically constant and approximately equal to 1. This means that only the horizontal distance has effect on the horizontal interpolation. This result is agreed with the orographic characteristics of the studied domain. Likewise thevalues obtained for γare closed to0.15, that is, the lower limit for this parameter which is related to low planetary boundary layers. However, the stability parameter αvaries in the interval 8-20. This range of values makes the wind predominantly flow more over the obstacles than around them. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 00:00 06:00 12:00 18:00 23:50 α/ 20 ε γ Fig. 5. Results of the estimation of α,εand γ(3-D model parameters) along the studied episode (March 21, 2003). In Figure 6 we see the evolution of the values of the three parameters of the 2.5-D model along the episode. These parameters givethe relationship between the roughness and the friction coefficient. The variation of these parameters with the meteorological conditions can be explain by certain hidden nonlinearity of the model, this means that the friction coefficient depends on the solution. As we can see in Figure 6, the variation of the parameters is not too high, so in practice we could assume constant values for a0,a1and a2and we would obtain similar results. 7.3 Comparison of model results with empirical data Once the main parameters are estimated, we start the wind modelling along the selected episode using the obtained values. For March 21, 2003, only measures from two control points, E242 and E283, were available. Figures 7 and 8 show the 21 Fig. 6. Results of the estimation of a0,a1and a2(2.5-D model parameters) along the studied episode (March 21, 2003). wind speeds obtained with the models and the reference values measured at the control points E242 and E283, respectively. More details of the errors of computed winds with respect to the measured wind may be seen in Table 5. We remark that the average errors at the measurement stations are small as expected. The average error for the 3-D model is 27.24% at control point E242 and 4.94% at E283. In addition, the average for the 2.5-D model is 1.05% and 43.28%, respectively. Then, the 2.5-D model obtain better results close to measurements points. However, the 3-D model is more accurate far from measurements points. Specifically this effect can be observed at the open boundaries of the domain, where the asymptotic model cannot fit very well the boundary conditions. Ideally the three graphics in figures 7 and 8 should be coincident. Of course the measured wind, as many experimental measurements, are obtained by sampling a Gaussian distribution. On average the measured wind verify mass conservation. With both models we obtain a wind field that verify mass conservation. 8 Conclusions We have presented two models for the wind field adjustment comparing the results by means of an example with real data. Both models are mass consistent models. The 3-D model needs a major use of empirical laws (initial interpolation, logarithmic profile, Pasquill stability class). Instead the 2.5-D model is a physical model since it is an asymptotic approximation of Navier-Stokes equations. The 3-D model is applicable in very general orographies since it admits all kinds of irregularities and discontinuities. The 2.5-D model is an asymptotic model, it is valid only for not very abrupt orographies. 22 0 5 10 15 20 00:00 02:00 04:00 06:00 08:00 10:00 12:00 14:00 16:00 18:00 20:00 22:00 00:00 m/s time Measured wind 3-D Computed wind 2.5-D Computed wind Fig. 7. Comparison of the wind velocities measured at the control station E242 (March 21, 2003). 0 5 10 15 20 00:00 02:00 04:00 06:00 08:00 10:00 12:00 14:00 16:00 18:00 20:00 22:00 00:00 m/s time Measured wind 3-D Computed wind 2.5-D Computed wind Fig. 8. Comparison of the wind velocities measured at the control station E283 (March 21, 2003). Both models give very good results in the measurements points as expected. In the control points we observe that the asymptotic model can provide very good results, even better that the 3D model in some cases, when the hypotheses of the model are satisfied. We have used a technique for constructing tetrahedral meshes which are simultaneously adapted to the terrain orography and the roughness length. The use of our refinement/derefinement process in the 2-D mesh corresponding to the terrain surface allows us to obtain meshes that are accurately adapted to different functions as well as are locally refined around several points. These characteristics of the generated meshes are very important in the wind simulation since, on the one hand, the quality of the representation of both orography and roughness is critical for obtaining accurate results with the two models, and on the other hand, the local refinement at the stations and control points is essential for inserting the wind data of the stations or recovering such data at any required point. 23 Stations and Average Average % Maximum Minimum Model control points measured computed average absolute absolute wind wind error error error E206 (49 m) 15.37 15.50 0.81 % 0.46 0.01 3-D 15.37 14.82 3.62 % 1.15 0.23 2.5-D E208 (15 m) 8.57 8.98 4.74 % 1.25 0.00 3-D 8.57 9.13 6.52 % 2.45 0.01 2.5-D E208 (30 m) 9.25 9.92 7.21 % 1.36 0.05 3-D 9.25 9.42 1.82 % 1.41 0.00 2.5-D E212 (15 m) 8.46 8.44 0.20 % 0.63 0.00 3-D 8.46 8.18 3.33 % 1.89 0.01 2.5-D E212 (30 m) 9.02 9.85 9.25 % 1.60 0.31 3-D 9.02 8.46 6.17 % 2.16 0.01 2.5-D E242 (40 m) 8.40 10.69 27.24% 5.09 0.09 3-D 8.40 8.49 1.05 % 2.54 0.01 2.5-D E283 (49 m) 13.62 12.95 4.94 % 3.04 0.02 3-D 13.62 7.73 43.28 % 9.88 1.34 2.5-D Table 5 Error of the computed wind at stations and control points. Some improvements for the 3-D model have been carried out in the construction of the initial wind based on the horizontal interpolation of wind measures and vertical extrapolation in stratified atmosphere. The optimization of the friction velocity for several measures in the same tower allows to minimize the differences between the constructed vertical profile of wind and the measures. However, though such differences are small, further research is needed in order to construct new wind profiles that exactly satisfy all the available measures of wind velocities. In addition, the inclusion of observations of turbulence intensities has made the model to be able of automatically updating the suitable wind profile as function of the corresponding stability class. The periodic updating of the main parameters of the models has proved to be fundamental for reducing the errors of the computed wind. However, further considerations should be taken into account in future works for a better performance of the models. For example, a finer map of roughness, a more sophisticated interpolation of wind velocities, a better approximation of the friction coefficient and a greater number of measurement stations well distributed over the studied region will help to reduce the errors of the models. In order to obtain an accurate wind 24 field in zones with very steep slopes, the mesh should be adapted to the contour lines, since a change in the direction of edges in the mesh may strongly affect the computed wind. In short, when the terrain is very rugged the 3-D model is recommended, howeverif the asymptotic assumptions are verified the 2.5-D model provides a good solution and has the advantage of incorporating thermal effects if required. Acknowledgement The work has been partially supported by Secretar´ ıa de Estado de Universidades e Investigaci´ on of the Ministerio de Educaci´ on y Ciencia and of the Ministerio de Ciencia e Innovaci´ on of the Spanish Government and FEDER, grant contracts: CGL2007-65680-C03-01,CGL2007-65680-C03-03,CGL2008-06003-C03-01and CGL2008-06003-C03-03, UNLP08-3E-010, the Junta de Castilla y Le´ on, grant number SA124A08 and the Instituto Tecnol´ ogico de Canarias. The authors are also grateful to Ignacio L´ainez and Antonio Ruiz for their technical support provided under the scope of the collaboration agreement signed by the University of Las Palmas de Gran Canaria and Desarrollos E´olicos (DESA) in January 2005. The wind data correspond to a DESA’s wind farm located in Lugo. References [1] Asensio MI, Ferragut L, Simon J (2005) A convection model for fire spread simulation. Appl Math Letters 18:673–677. [2] B¨ack T, Fogel DB, Michalewicz Z (1997) Handbook of evolutionary computation. Oxford Univ. Press, New York-Oxford. [3] Businger JA, Arya SPS (1974) Heights of the Mixed Layer in the Stable, Stratified Planetary Boundary Layer. Adv Geophys 18A:73–92. [4] Bresch D., Lemoine J., Simon J. (1999) A vertical diffusion model for lakes. Siam J. Math Anal 30:603-622. [5] Davis L (1991) Handbook of genetic algorithms. Van Nostrand Reinhold. [6] Escobar JM, Montenegro R (1996) Several aspects of three-dimensional Delaunay triangulation. Adv Eng Soft 27(1/2):27–39. [7] Escobar JM, Rodr´ıguez E, Montenegro R, Montero G, Gonz´alez-Yuste JM (2003) Simultaneous untangling and smoothing of tetrahedral meshes. Comp Meth Appl Mech Eng 192:2775–2787. 25