A Mathematical Model for Wind Velocity Field Reconstruction and Visualization Taking into Account the Topography Influence

Abstract

In this paper, we propose a global modelling for vector field approximation from a given finite set of vectors (corresponding to the wind velocity field or marine currents). In the modelling, we propose using the minimization on a Hilbert space of an energy functional that includes a fidelity criterion to the data and a smoothing term. We discretize the continuous problem using a finite elements method. We then propose taking into account the topographic effects on the wind velocity field, and visualization using a free library is also proposed, which constitutes an added value compared to other vector field approximation models.

1. Introduction

Vector field approximation has many applications, such as to predict wind turbine production; in oceanography, to study marine currents; and more generally in a computer sciences framework. As introduced in [1], in order to approximate a vector field, several approaches have been developed: a finite element method interpolation (see [2]), PDE-based methods, kriging methods, a Lagrange interpolation method and spline and Rational Basis Function (RBF) approximations (see [3,4,5,6,7,8]). These approaches have drawbacks, particularly when a small number of data are available and the approximation’s result is qualitatively insufficient. In this work, we are precisely in a case where the number of data from anemometers is considered low in comparison to the large study area. It is, therefore, necessary to propose a mathematical model allowing for this type of data to be processed via a robust energy functional minimization. The major methodological contribution of this work consists in the modelling using D^m splines, as well as the contribution of adding the topography effect into the numerical results; most of the methods proposed in the literature do not integrate this important aspect, but because there are little available data, it is important to include as much information in the model as possible to generate a realistic wind field, and taking the topography into account is quite simple and brings significant added value. Other methods (like [3,4,5]) only focus on the modelling and the mathematical aspects of the approximation without adding more information than the input dataset. To our knowledge, our approach is the first to offer such a global framework.

In this paper, we first give the considered vector field approximation model, using a D^m spline operator, rigorously introduced in [1,7,9,10,11]. The dataset consists of a finite set of vectors (xi, yi, wi)_i, where (xi, yi) locates the point in 2D and wi = (Ui, Vi) ∈ R² gives the direction and speed of the wind at location (xi, yi). A minimization problem is introduced, leading to a variational problem whose solution is the searched for wind vector field. We give the discretization using a classic finite element method. We then give details on how to integrate into the model the effect of the topography on the obtained wind velocity field; this last part greatly improves upon previous models (given in [1]), where topographic effects are not taken into account. We show the effects of this topographic influence on the synthetic dataset given by an explicit function. We then give numerical examples for a real dataset, including a specific tool for visualization. A global review of this complete approximation framework is given in Figure 1.

2. Mathematical Modelling

The mathematical modelling of our approximation problem is constructed using a D^m spline operator as follows: For all $v\in {\mathrm{H}}^{m+1}\left(\mathsf{\Omega }\right)$, we introduce the energy functional consisting of two terms. The first one is the data fidelity criterion, while the second one is the smoothing parameter:

$$J_{\epsilon }\left(v\right)={〈\rho \left(v-w\right)〉}^2+{\epsilon {\left|v\right|}^2}_{m+1,\mathsf{\Omega },R}$$

(1)

where Ω is an open subset corresponding to the studied zone. w = (w₁, w₂, …, w_N) $\in {\left({\mathrm{R}}^2\right)}^N$ is the vector field dataset, $\rho$ is a linear operator linked to the dataset, ${{\left|•\right|}^2}_{m+1,\mathsf{\Omega },IR}$ is the usual semi-norm of the usual Sobolev space H^m+1(Ω) as defined in [1], $〈•〉$ is the Euclidean norm in R² and ε is a smoothing parameter generally equal to 10⁻⁶ in many applications (see [1,4,7,8] for more details). We recall that H^m+1(Ω) is a space of functions belonging to L²(Ω) and their (m + 1) derivatives. We also introduce an ordered set of N points (x_i)_i of Ω, where we suppose as known the value of the wind velocity field. The linear operator $\rho$ is given by $\rho (v)=(v(x_1),v(x_2),\dots v(x_N))\in {\left({\mathrm{R}}^2\right)}^N.$ We use the D^m spline (see Gout et al. [1] and Arcangéli et al. [7] for a complete study of this approximation operator) approximation framework to solve this problem. We call ${\sigma }_{\epsilon }^d$ the smoothing spline on Ω relative to $\rho$, which is the unique solution of the minimization problem:

$$\left\{\begin{array}{c} \mathrm{find} {\sigma }_{\epsilon }\in {\mathrm{H}}^{m+1}(\mathsf{\Omega }), \mathrm{such} \mathrm{that} \mathrm{for} \mathrm{any} v \in { \mathrm{H}}^{m+1}(\mathsf{\Omega }):\\ J_{\epsilon }\left({\sigma }_{\epsilon }\right)\le J_{\epsilon }\left(v\right).\end{array}\right.$$

(2)

We can use the Lax–Milgram theorem to establish the uniqueness of this minimization problem, since the solution ${\sigma }_{\epsilon }^d$ of this minimization problem is the solution of the following variational problem

$$\left\{\begin{array}{c}\mathrm{find} {\sigma }_{\epsilon }\in {\mathrm{H}}^{m+1}(\mathsf{\Omega }), \mathrm{such} \mathrm{that} \mathrm{for} \mathrm{any} v \in { \mathrm{H}}^{m+1}(\mathsf{\Omega }):\\ 〈\rho {\sigma }_{\epsilon }^d,\rho v〉+\epsilon {\left({\sigma }_{\epsilon }^d,v\right)}_{m+1,\mathsf{\Omega },R}=〈w,\rho v〉.\end{array}\right.$$

(3)

where ${\left(•,•\right)}_{m+1,\mathsf{\Omega },R}$ denotes the semi-norm of ${\mathrm{H}}^{m+1}\left(\mathsf{\Omega }\right)$. To apply the Lax–Milgram theorem, we recall that all the hypotheses of this theorem are satisfied, since ${\mathrm{H}}^{m+1}\left(\mathsf{\Omega }\right)$ is a Hilbert space, and we also have as follows:

• $a(u,v)=〈\rho u,\rho v〉+\epsilon {\left(u,v\right)}_{m+1,\mathsf{\Omega },R}$ is a bilinear form, being the sum of scalar products.
• a(u,v) is continuous on ${\left({\mathrm{H}}^{m+1}\left(\mathsf{\Omega }\right)\right)}^2$ because $\left|a(u,v\right|)\le \max (1,\epsilon ){‖u‖}_{m+1,\mathsf{\Omega }}{‖v‖}_{m+1,\mathsf{\Omega }}$, using the Cauchy–Schwarz inequality and the norm equivalence between ${\left(〈\rho v,\rho v〉+{\left(v,v\right)}_{m+1,\mathsf{\Omega },R}\right)}^{1/2}$ and ${‖v‖}_{m+1,\mathsf{\Omega }}$.
• a(v,v) is elliptic on ${\mathrm{H}}^{m+1}\left(\mathsf{\Omega }\right)$ since $a(v,v)\ge \min (1,\epsilon ){‖v‖}^{2}m+1,\mathsf{\Omega }.$
• $〈w,\rho v〉$ is a continuous linear form.

We now propose a discretization of the variational problem using a finite element discretization (see [1,9,12] for more details on such a discretization). We recall that the main idea of the finite element method is to replace the (Hilbert) space ${\mathrm{H}}^{m+1}\left(\mathsf{\Omega }\right)$ used to define the variational Equation (3) by a finite dimensional subspace V_h. Of course, we have ${\mathrm{V}}_h\subset {\mathrm{H}}^{m+1}\left(\mathsf{\Omega }\right)$. The functions belonging to V_h. are piecewise polynomials, and the bases of the functions for the space V_h are constructed such that they have small support. For any real h > 0, let T_h be a triangulation of Ω by n-simplices or n-rectangles K with diameter h_K ≤ h. We classically approximate the space H^m+1(Ω) by the space V_h,, a finite dimensional space included in H^m+1(Ω) and admitting a polynomial basis of polynomial functions ${\left({\mathsf{\Phi }}_j\right)}_j$. We write the solution of ${\sigma }_{\epsilon ,h}^d$ on the basis of V_h as

$${\sigma }_{\epsilon ,h}^d={\displaystyle \sum _{j=1}^{\dim V_h}{\beta }_j{\mathsf{\Phi }}_j}\hspace{1em}:\hspace{0.33em}\hspace{0.33em}{\beta }_j\in R.$$

(4)

Note that the polynomials ${\left({\mathsf{\Phi }}_j\right)}_j$ are given, since they are computed following the chosen generic finite element. The generic finite element we choose here is the Bogner–Fox–Schmit (BFS) rectangle of class C¹, where a function of V^h. is completely determined by its four values (value, values of the two first derivatives and value of the twist derivative) at a nodal point (see Appendix A). The choice of the BFS finite element is due to their capability to easily tessellate rectangular domains and to guarantee a final approximation of class C¹.

From (4), we have to find ${\left({\beta }_j\right)}_j$ in order to find the solution ${\sigma }_{\epsilon ,h}^d$ of our approximation problem. We can now give the discretization of problem (3) using (4):

$$\begin{array}{l}Find\hspace{0.33em}{\left({\beta }_j\right)}_j\in \hspace{0.33em}{R}^{\dim V_h}\hspace{0.33em}such\hspace{0.33em}\hspace{1em}that\hspace{1em}\forall k=1,\dots ,\dim V_h,\\ {\displaystyle \sum _{i=1}^N}〈{\displaystyle \sum _{j=1}^{\dim V_h}{\beta }_j.}{\mathsf{\Phi }}_j,v_h(x_i)〉+\epsilon {\left({\displaystyle \sum _{j=1}^{\dim V_h}{\beta }_j.}{\mathsf{\Phi }}_j,v_h\right)}_{m+1,\Omega ,R^n}={\displaystyle \sum _{i=1}^N〈w_i,v_h(x_i)〉}.\end{array}$$

(5)

Equation (5) leads to the following linear system, taking as the test function $v_h$ all the basis functions ${\mathsf{\Phi }}_k$: k = 1,…, dim V_h

$${\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=1}^{\dim V_h}{\beta }_j.}{\mathsf{\Phi }}_j,{\mathsf{\Phi }}_k(x_i)}+\epsilon {\displaystyle \sum _{j=1}^{\dim V_h}{\beta }_j.}{\left({\mathsf{\Phi }}_j,{\mathsf{\Phi }}_k\right)}_{m+1,\Omega ,R^n}={\displaystyle \sum _{i=1}^N〈w_i,{\mathsf{\Phi }}_k(x_i)〉}.$$

(6)

We finally have to solve the following linear system to find the unknown real values ${\left({\beta }_j\right)}_j$:

$$\begin{gathered} \left({\left(A\right)}^{T}A+\epsilon R\right)\beta ={\left(A\right)}^{T}w, \\ \mathrm{with} A={\left[{\mathsf{\Phi }}_j(x_i)\right]}_{1\le i\le N,1\le j\le \dim V_h} \mathrm{and} R={\left[{({\mathsf{\Phi }}_j,{\mathsf{\Phi }}_j)}_{m+1,\Omega ,R}\right]}_{1\le i,j\le \dim V_h}. \end{gathered}$$

(7)

For the numerical simulation, we take m = 1, and we use the Bogner–Fox–Schmit finite element with the basis function as a polynomial of degree 3 (see [7] for more details). The modelling we have introduced in this section permits the approximation of a wind velocity field on all Ω from a finite set of data given, for instance, by several anemometers (as illustrated in the numerical section of this work—Figure 10 and Figure 11).

3. Taking into Account the Topography

In Section 2, we proposed a framework to approximate a wind velocity field from a finite set of measures. It is, of course, well known that the topography plays an important role in wind field velocity variations. Obstacles modify air flows due to pressure forces (see Figure 2). The wind slows down upstream of an obstacle, and accelerates downstream of it.

Since the approach introduced in Section 2 does not take into account physical considerations, we propose here a way to approximate the topography’s influence by post processing the approximated wind field obtained from the model given in Section 2. More precisely, let us consider a wind vector field on N points as

W’ = (w’₁,…, w’_N),

where each wind vector w’_i belongs to R x R, with given coordinates (x_i, y_i) for each w’_i, according to the topographic configuration around it as follows: W’_i = c_i(θ) w_i, where W’_i is the adjusted wind vector and c_i(θ) is the coefficient computed from the topographic configuration at point (x_i, y_i), depending on the wind direction θ. This approximation holds for local topographic effects. It cannot take into account large-scale effects, such as Venturi effects in valleys or straits. To compute the topographic coefficient c_i for a given wind direction θ at point (x_i, y_i), we used the formulas given in parts 1–4 of [13], depending on the slope Φ

$$\left\{\begin{array}{ll}1,& \Phi <0.05\\ 1+2s\Phi ,& 0.05<\Phi <0.3\\ 1+0.6s,& \Phi >0.3\end{array}\right.$$

(8)

where s is the characteristic coefficient of the obstacle, depending on its features (see

Table 1. Features used to compute the characteristic coefficient s of an obstacle (source: [13]).

Variables	Definition
s	Orographic location factor
Φ	Upwind slope H/Lu in the wind direction (see Figure 3)
Le	Effective length of the upwind side
Lu	Length of the upwind side
Ld	Length of the downwind side
H	Effective height of the obstacle
x	Horizontal distance between point (x,y) and the top of the obstacle
z	Height of the considered point (x,y)

and Figure 3).

The effective length L_e is computed as follows (type of slope Φ = H/L_u):

$$L_e=\left\{\begin{array}{c}{L}_u,\hspace{1em}0.05<\mathsf{\Phi }<0.3\hfill \\ {\displaystyle \frac{H}{0.3}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{1em}\mathsf{\Phi }>0.3.\end{array}\right.\hspace{0.33em}$$

(9)

We also have to compute the value of the orographic location factor s used in (8). As shown in [13], the value of s is related to the ratio H/L_e. More precisely, for an upwind section, for ranges $-1.5\le {\displaystyle \frac{x}{L_u}}\le 0\hspace{0.33em}\hspace{0.17em}\mathrm{and}\hspace{0.17em}\hspace{0.17em}0\le {\displaystyle \frac{z}{L_e}}\le 2\hspace{0.17em}$ we take $s=A\exp \left(B{\displaystyle \frac{x}{L_u}}\right),\hspace{0.17em}$ where

$$A=0.1552{\left({\displaystyle \frac{z}{L_e}}\right)}^4-0.8575\hspace{0.17em}{\left({\displaystyle \frac{z}{L_e}}\right)}^3+1.8133{\left({\displaystyle \frac{z}{L_e}}\right)}^2-1.9115\left({\displaystyle \frac{z}{L_e}}\right)+1.0124,$$

$\mathrm{and}\hspace{0.33em}B=0.3542{\left({\displaystyle \frac{z}{L_e}}\right)}^2-1.0577\left({\displaystyle \frac{z}{L_e}}\right)+2.6456.$ Note that when ${\displaystyle \frac{x}{L_u}}\le -1.5\hspace{0.17em}\hspace{0.17em}\mathrm{or}\hspace{0.17em}\hspace{0.17em}2\le {\displaystyle \frac{z}{L_e}}\hspace{0.17em}$ we take s = 0.

For a downwind section, as shown in [13], we take $s=A{\left(\log {\displaystyle \frac{x}{L_e}}\right)}^2+B\log {\displaystyle \frac{x}{L_e}}+C,\hspace{0.17em}$ with

$$A=-1.342{\left(\log {\displaystyle \frac{z}{L_e}}\right)}^3-0.822\hspace{0.17em}{\left(\log {\displaystyle \frac{z}{L_e}}\right)}^2+0.4609\log {\displaystyle \frac{z}{L_e}}-0.0791,$$

$$B=-1.0196{\left(\log {\displaystyle \frac{z}{L_e}}\right)}^3-0.891{\left(\log {\displaystyle \frac{z}{L_e}}\right)}^2+0.5343\log {\displaystyle \frac{z}{L_e}}-0.1156,$$

$$\mathrm{and}\hspace{0.17em}\hspace{0.33em}C=0.803{\left(\log {\displaystyle \frac{z}{L_e}}\right)}^3+0.4236{\left(\log {\displaystyle \frac{z}{L_e}}\right)}^2-0.5738\log {\displaystyle \frac{z}{L_e}}+0.1606.$$

To compute the topographic coefficients of a domain D = [0, 1]² × R², we consider a regular grid D_h of D of step h. Then, for each point (x_i, y_i) = (ih, jh) in D_h, we compute a coefficient for each wind direction θ. We split up the compass wind into eight directions θ_j from 0 to 360°, by steps of 45 degrees (see Figure 4).

Once the collection of topographic coefficients is computed using (8), we use it to adjust the approximated wind field, selecting coefficients according to the direction of each wind vector and using them on the obtained approximating wind velocity field.

4. Numerical Examples

In this section, we give several numerical examples, including the computation of the topographic coefficients (using a given function f) and the approximation of a vector field from a finite set of vectors giving the direction and speed of the wind.

4.1. Computation of the Topographic Coefficients

In order to illustrate the proposed methodology on synthetic data given by an explicit given function f, we simulate an obstacle (hills) in domain D using the basic 2D function f defined as follows:

$$\begin{array}{c}f(x,y)={\displaystyle \frac{3}{4}}\exp \left(-{\displaystyle \frac{1}{4}}{\left(-9x-2\right)}^2-{\displaystyle \frac{1}{4}}{\left(9y-2\right)}^2\right)\\ +{\displaystyle \frac{3}{4}}\exp \left(-{\displaystyle \frac{1}{49}}{\left(-9x+1\right)}^2-{\displaystyle \frac{1}{4}}{\left(9y+1\right)}^2\right)\\ +{\displaystyle \frac{1}{2}}\exp \left(-{\displaystyle \frac{1}{4}}{\left(-9x-7\right)}^2-{\displaystyle \frac{1}{4}}{\left(9y-3\right)}^2\right)\\ -{\displaystyle \frac{1}{5}}\exp \left(-{\displaystyle \frac{1}{4}}{\left(-9x-4\right)}^2-{\displaystyle \frac{1}{4}}{\left(9y-7\right)}^2\right)\end{array}$$

(10)

We define the discretized domain as D_h, with h = 1/n. The obstacle is obtained by computing f for every couple point (x_i, y_i), where i, j = 0, …, n. For n = 40, we obtain the following obstacle (see Figure 5), and we give on this image the computed topographic coefficients obtained using (10) (considering an arbitrary wind direction indicated by the red arrow).

The simulated topographical data given by the function f in (10) are then used to compute the topographic coefficients as described in the previous subsection. For each of the eight wind directions, we can plot the color map of the computed topographic coefficients on the obstacle. The associated colors go from dark, for zones where the topography slows down the wind flow (c < 1), to white, for zones where the topography accelerates the wind flow (c > 1). For instance, we have plotted color maps for situations where the wind comes from the northeast (see Figure 6), west (see Figure 7) and southeast (see Figure 8). For each figure, the arrow indicates the direction of the wind.

Note that the function given in (10) permits illustrating the topographic effect in cases where we consider hills and ridges (as in Figure 3). We have tested the topographic effect on more vertical hills and the results were satisfying. For the case of cliffs or escarpments (Figure 9), the computation is slightly different but the reasoning is analogous, and the values of the corresponding orographic parameters are given in Figure 10. An improvement could be to propose calculations to take into account the influence of buildings (vertical walls of buildings with a rectangular plan, the influence of the angle of roofs, etc.) or of vegetation (trees, etc.), especially if we want to reconstruct the wind on a micro-scale.

We now give several numerical experiments on real datasets. The experiments are performed on a 2.21 GHz Athlon with 1.00 GB of RAM.

We focus on the data for wind vector fields acquired in northwest France; the dataset takes into account eight weather stations (Meteo France, Figure 10, Figure 11 and Figure 12).

4.2. Numerical Simulations of the Global Algorithm

In this subsection, from a set of six velocity wind data, we give the approximation obtained by the method in Section 2, and we then compute the topographic coefficients of the studied zone using the method given in Section 3.

Here is some information about the numerical examples:

• Dataset: six anemometers located at six airports giving the direction and speed of the wind, see Figure 10 and Figure 11;
• Parameter $\epsilon$ = 0.000001;
• Generic finite element: Bogner–Fox–Schmit of class C¹ (See Appendix A);
• Studied domain: [3500, 6000] × [2.44, 2.62];
• Meshing: 4 × 4 rectangles and 3 × 3 rectangles. The results are given in Figure 12.

The choice of the finite element meshing is crucial, and it must be linked to the number of data we have in the input. For a grid of 3 × 3 rectangles, we have 9 rectangles, 16 nodes and, as we have four basis functions per node with the BFS finite element of class C¹, the dimension of the space V_h is equal to 64 (while it is 100 with a 4 × 4 meshing, leading to 16 rectangles and 25 nodes). As we do not have a large amount of data, we choose a low number of rectangles in our mesh. If we choose a finer grid, the approximation error increases.

We compute the quadratic error given by the following quotient

$$Quad\_Error={\left({\displaystyle \frac{{{\displaystyle \sum _{i=1}^N〈{\sigma }_{\epsilon ,h}^d(a_i)-w_i〉}}_2^2}{{\displaystyle \sum _{i=1}^N{〈w_i〉}_2^2}}}\right)}^{1/2},$$

where ${〈•〉}_2$ denotes the Euclidean scalar norm. In all our different tests, the quadratic errors is of 10⁻³ and 10⁻⁴ orders, which is considered as very good in the context of vector field approximation. We then compute this obtained approximation using the topography of the considered zone (see Figure 7 and Figure 13 for the result).

In order to show this method on more complicated datasets, we consider the wind conditions over 90 h; we have the value of the wind vector field at each Meteo France station every 3 h (total of 30 datasets). We apply the previous method for each time step. We then obtain the approximated wind velocity field over the 90 h. We have to propose a way to visualize such datasets.

5. Visualization

To obtain a simulation on time using a free library, we first propose using Matplotlib using Python. The following code was developed at INSA Rouen Normandie by the authors (and thanks to H. Merelle from the Applied Math. Department for his help). The main advantage of this code is that it gives a complete framework from the input (dataset) to the numerical simulation, including the approximation using the spline functions, finite element methods and the topography influence.

Algorithm for visualization using Matplolib [14].

Here is the list of files necessary for correct processing and the different steps of the proposed method linked to the flowcharts given in Figure 1:

-

Initial Input: dataset (xi, yi, (Ui, Vi))i.

-

Definition of Ω—meshing of Ω with rectangles (as we use rectangular finite elements): the number of subdivisions is linked to the number of data; in the examples here, it is 3 × 3 or 4 × 4 subdivisions in x and y.

-

D^m spline approximation: the output is the evaluation of the vector field on each point of a fine grid of Ω.

-

Computation of topography effect on the vector field: output.txt file.

-

Script_visualization.py

-

The “output.txt” file (in the same folder).

The purpose of the program is to visualize a vector flow from text files, with the possibility of adding a background (topography, etc.).

Data conditions:

-

For an animation:

• The “output.txt” file is of the form

  X1 Y1 U1_1 V1_1 U1_2 V1_2…

  X2 Y2 U2_1 V2_1 U2_2 V2_2…

  with X1 and Y1 being the first coordinates, followed by U1_1 V1_1 U1_2 V1_2…; the different sizes of the vectors are a function of time.

• To execute in a terminal under Ubuntu, we use the Python script_visualization.py, with the following instructions:

  ○

  The title: it represents the file name (when exporting) and the title of the figure.

  ○

  The size of the vector arrows: the bigger they are, the smaller the vectors appear.

  ○

  The number of images: if your output file is of the form “Animation”, in this case you will have the following question, “Enter the number of frames per second”; it determines the frame rate per second.

  ○

  For the background: “O” for accept or “N” otherwise.

  ○

  For the name of the image you must give the file extension: here is a non-exhaustive list of usable formats: [name].png, [name].jpg, [name].jpeg and [name].gif.

  ○

  To display the result: This command is only used to show you the result. The result is still saved even if you do not display it.

  ○

  Data output: For an animation, you can find the animation in the folder in the form [title].gif, and for a fixed image, you can find the rendering in the folder in the form [title].png.

In order to show a numerical simulation, we give a simulation for all of the Normandy region (wind velocity field using the Meteo France dataset) using this Python script and using Matplotlib.

Examples of the obtained visualizations on a sequence (time) of an approximated velocity vector field (test in the Normandy region, France) and marine current (Seine River at Rouen, France) are given in Figure 14 and Figure 15.

6. Discussion and Future Directions

About the approximation method using D^m splines, note that a theoretical study of the error in the approximation method following the used finite element mesh is a work in progress. To do that, we use previous results obtained from smoothing spline approximations from a finite set of points, as performed in [10,11]. Another development, linked to the approximation part, will be to propose an automatic meshing of the domain and different choices for the used generic finite element (based on triangles, etc.).

Another goal consists in improving both the modelling and visualization. About the modelling, the goal will be to include new kind of datasets; nowadays, it is possible to obtain wind datasets from Lidar located on wind turbines (see Figure 16). This dataset gives the wind velocity field with a specific geometry: along a spiral. It makes the computation much more difficult because it requires a specific finite elements meshing, which makes the process much less automatic.

We also plan to add a smooth visualization based on texture using the Matplotlib library. Another crucial point consists in the effect of using a smooth visualization of the flow using streamlines. Vector field data are produced by scientific experiments and numerical simulations, which are now widely used to study complex dynamic phenomena, using a robust method to visualize steady flow field with both line representations and textures. In [15], the authors specify that “a streamline is a line tangential to the vector field at any point. Covering an image with a set of streamlines is a very good way to visualize the flow features” (see Figure 17 and [15] for more details).

In order to show this method on more complicated datasets, we considered the wind conditions over 90 h. We have the value of the wind vector field at each Meteo France station every 3 h (total of 30 datasets equivalent to the one we show in Section 4).

We then computed the obtained vector field using the modelling proposed in this work with the help of B. Jobard [15], and we obtained the movie given in [16]. This result is smooth and promising. But improvements have to be made to propose a tool able to treat the whole process with the same software, and to maybe try other approximation methods (like the one in [17], and to mix this approximation/visualization tool with an image processing framework, or the one in [18,19,20], using radial basis functions; this is ongoing work). Moreover, it is also crucial to develop an algorithm with which to approximate more complex datasets, like Lidar ones (instead of anemometers) to compute wind velocity fields.

7. Conclusions

In this work, we successfully proposed a global tool, from vector field approximation to visualization. We proposed a method to obtain a visualization of a vector field from a sparse dataset, after computing its numerical approximation using a mathematical model using energy minimization and finite elements for the discretization. Note that we also integrated the topography effect into the modelling of a wind velocity field approximation method. To our knowledge, this is the first global approach for such numerical simulations from a dataset with few data.

As stressed in Section 6, several developments should occur in the future in order to improve this global approach. Many potential applications exist, from velocity wind approximation for wind turbine energy modelling, current simulation for modelling the morphodynamics of coastal zones and control theory for vehicle navigation (cars, submarines, etc.).