Development and application of a mesh generator intended for unsteady vortex-lattice method simulations of wind turbines and wind farms

In the last decades, the unsteady vortex-lattice method (UVLM) has gained a lot of acceptance to study large onshore–offshore wind turbines (WTs). Furthermore, and due to the development of more powerful computers, parallelization strategies, and algorithms like the fast multipole method, it is possible to use vortex-based methods to analyze and simulate wind farms (WFs). However, UVLM-based solvers require structured meshes, which are generally very tedious to build using classical mesh generators, such as those utilized in the context of finite element methods (FEMs). Wind farm meshing is further complicated by the large number of design parameters associated with the wind turbine (coning angle, tilt angle, blade shape, etc.), farm layout, modeling of the terrain topography (for onshore WFs), and modeling of the sea level surface (for offshore WFs), which makes the use of FEM-oriented meshing tools almost inapplicable.

In the literature there is a total absence of meshing tools when it comes to building aerodynamic grids of WTs and WFs to be used along with UVLM-based solvers. Therefore, in this work, we present a detailed description of the geometric modeling and computational implementation of an interactive UVLM-oriented mesh generator, named UVLMeshGen, developed entirely in MATLAB^® and easily adaptable to GNU OCTAVE, for wind turbines and onshore–offshore wind farms. The meshing tool developed here consists of (i) a geometric processor in charge of designing and discretizing an entire wind farm and (ii) an independent module in charge of computing the kinematics for the entire WF. The output data provided by the UVLMeshGen consist of nodal coordinates and connectivity arrays, making it especially attractive and useful to be used by other flow potential solvers using vortices, sources and sinks, or dipoles/doublets, among others. The work is completed by providing a series of aerodynamic results related to WTs and WFs to show the capabilities of the mesh generator, without going into detailed discussions of wind turbine aerodynamics, which are not the focus of this paper. The meshing tool developed here is freely available under a Creative Commons Attribution 4.0 International License (Roccia, 2023).

1 Introduction

For some years now, wind energy has become one of the fundamental pillars on the world stage of renewable energy. This fact has been materialized by an increasing number of different wind turbine (WT) designs: going from small wind turbines like the Vestas V27 of 200 KW (Torabi, 2022) and moderately sized designs for onshore applications to the current large offshore WTs such as the Vestas V236-15.0MW prototype (Vestas, 2022) or the CSSC Haizhuang H260-18MW concept design with a rotor diameter of 260 m (CSSC Haizhuang, 2023).

One of the most important challenges of wind turbine technology is the accurate characterization of the WT loads under inflow conditions that trigger complex aerodynamic effects. Although the description of the flow surrounding a wind turbine has been a subject of interest for many years, the study of three-dimensional and expensive-to-model unsteady aerodynamics of WTs and wind farms is still an active field of research (Muñoz-Simón et al., 2022). Throughout the years, a wide variety of aerodynamic models for wind turbines have been proposed, implemented, verified, validated, and successfully applied. They range from basic approaches such as those based on blade element momentum (BEM) theory, widely spread through the industry for the initial design loops, to advanced high-fidelity models using computational fluid dynamics (CFD) techniques.

While classical and enhanced versions of BEM-based solvers have been found to provide good agreement with measurements and CFD simulations, they require a series of engineering corrections to model challenging unsteady aerodynamic phenomena of increasingly large WTs (Perez-Becker et al., 2020). The inaccuracies observed in BEM simulations, when compared to more sophisticated approaches, are a natural consequence of the underling theory behind the method (Hansen, 2015). Instead, high-fidelity CFD computations can capture more flow physics, thus providing a better prediction accuracy than BEM codes (Nigam et al., 2017; Liu et al., 2017). However, solving the full Navier–Stokes equations for three-dimensional unsteady flows with boundaries undergoing large and complex motion is by far a challenging and time-consuming task requiring typically 2–5 weeks on average (Terry, 2018).

As an intermediate option between the BEM and CFD approaches, we introduce the potential flow solvers, among which the so-called vortex-lattice methods (VLMs) represent a good alternative to assess the aerodynamic performance of different aeronautical–mechanical engineering applications. Its extension to the study of transient aerodynamic loads for slender lifting surfaces undergoing complex motions, the well-known nonlinear unsteady vortex-lattice method (UVLM), has proven to be a more than viable option, presenting an excellent trade-off between precision and computational cost (Verstraete et al., 2023). Furthermore, UVLM-based solvers have been continuously gaining ground in the context of those problems, in which free-wake methods become a necessity due to the geometric complexity and the presence of large displacement–rotations, such as morphing wings (Verstraete et al., 2015, 2019), flapping wings (Stanford and Beran, 2010; Roccia et al., 2013; Nguyen et al., 2016), rotorcraft (Wie et al., 2009; Colmenares et al., 2015; Lee et al., 2022), wind turbines (van Garrel, 2003; Gebhardt et al., 2010; Gebhardt and Roccia, 2014), and non-conventional wind energy devices (Abdelkefi et al., 2014; Beltramo et al., 2020; Roccia et al., 2020), among others.

Among the most promising UVLM-based solvers capable of performing aerodynamic analysis of wind turbines, we can mention OpenVOGEL (Hazebrouck, 2023), WinDS (WinDS, 2023), GSF-Aero (Verstraete et al., 2023), the general-purpose framework developed by Pérez Segura et al. (2020), and VLMSim (Verstraete et al., 2023). The latter one can handle aerodynamic studies of arbitrary onshore and offshore wind farms. Although UVLM frameworks are mid-fidelity simulation tools with an enormous potential to be used in the wind energy sector, they require structured meshes, which are very tedious to build by using classical meshing codes as those utilized in the context of the finite element method (FEM), such as Gmsh, CUBIT^®, MeshLab, and GID^®, among others.

In addition, wind turbines are characterized by a large number of design parameters, such as the coning angle, the tilt angle, multiple airfoils defining the blade, the twist angle, and pre-bent and pre-sweep shapes, etc. In this sense, a proper wind turbine meshing process should incorporate an easy way to handle such information. When considering wind farms, the meshing is further complicated by the need of including parameters associated with the farm layout, terrain topography (for onshore wind farms), and the description of the sea level surface for offshore wind farms. Another key point, and by no means less important, is the generation of the kinematics for the entire wind farm. This aspect includes everything from basic rotor kinematics and laws of motion for yaw and pitch (if any) to sea level surface kinematics (to simulate waves) and substructure motions for floating wind turbines (Sant and Cuschieri, 2016; Lee and Lee, 2019). On this basis, a versatile UVLM-oriented meshing tool for onshore and offshore wind farms must necessarily involve design plus discretization plus kinematic modules to provide all the data required by any standard UVLM engine.

To the best of our knowledge, there is to date no freely available UVLM-oriented mesh generator intended for arbitrary wind farms that allows for (i) designing wind park layouts, (ii) considering different wind turbines (with their own design parameters), (iii) including the terrain topography and/or the sea surface description, and (iv) computing the wind farm kinematics. The only attempts of which we are aware of are those individual efforts to mesh specific geometries and limited meshing tools already incorporated into UVLM codes such as OpenVOGEL or modules developed in-house at companies that are inaccessible to the wind energy community.

In this work, we present a detailed description of the geometric modeling and computational implementation of an interactive UVLM-oriented mesh generator for onshore and offshore wind farms, hereafter referred to as UVLMeshGen. The meshing tool, fully developed in the MATLAB^® language and easily adaptable to GNU OCTAVE, allows for the generation of structured and conformal aerodynamic grids of wind farms, including the terrain and/or sea level surface modeling. The structured mesh data provided by UVLMeshGen consist of nodal coordinates and connectivity arrays, similar in some way to classical FEM-oriented mesh generators. In this regard, such output data are not limited to only UVLM-based simulation frameworks, but they can be used by any potential flow solver relying on other singularities elements, such as vortices, sources and sinks, or dipoles/doublets, among others. Besides the geometric processor, UVLMeshGen has an additional independent module in charge of computing the kinematics for the entire wind farm based on user-defined input data. Furthermore, this meshing engine allows the addition of user-defined scripts or add-ons to post-process the aerodynamic grids and/or to extend their data structure with additional information to be used as input data in any kind of UVLM solver. All these features make the meshing tool presented here a valuable resource in larger projects and endeavors such as AVATAR (Advanced Aerodynamic Tools for Large Rotors, 2013–2017), which pursued, as the main goal, the assessment of different aerodynamic models for large (10 MW+) wind turbines (Schepers, 2015) or the CRC 1463 Offshore Megastructures (Integrated design and operation methodology for offshore megastructures, 2021–2024) targeting an integrative design and operation for offshore wind turbines (20 MW+) (Hannover, 2021). UVLMeshGen is freely available under a Creative Commons Attribution 4.0 International License (Roccia, 2023).

The remainder of this work is organized as follows: in Sect. 2, we introduce the geometric entities and the meshing process to obtain the aerodynamic grid of each component of a wind turbine. In Sect. 3, we describe the main aspects associated with the computational implementation of UVLMeshGen. In Sect. 4, we present a series of aerodynamic simulations to show the capabilities of the mesh generator, without entering into quantitative discussion about the aerodynamics of wind turbines. Finally in Sect. 5, we provide conclusions and future work to be addressed in a follow-up paper.

2 Geometric modeling

A wind turbine is characterized by a large number of parameters. A proper aerodynamic analysis of such mechanical systems by UVLM-based codes necessarily requires an accurate description of the wind turbine surfaces. On this basis, this section first presents the geometric object to be used to represent the surfaces of a wind turbine. This subsection is followed by a full description of the geometric modeling of each component of a wind turbine (e.g., hub, nacelle, blades, tower).

2.1 Geometric preliminaries

In the case of using boundary integral equation methods (BIEMs) to predict the aerodynamic forces and wake structures of very complex engineering systems, an accurate description of their boundaries (solid surfaces) is mandatory. In particular, the aerodynamic analysis of wind turbine farms by using BIEMs requires providing (i) precise data associated with the discretization of their boundaries (ground, rotors, towers, etc.); (ii) additional data according to the method adopted (collocation points, shedding zones, type of singularities, etc.); and (iii) kinematics of the system (positions and velocities over time).

In general, the main components of a wind turbine are the hub, the nacelle, the blades, the tower, and the ground where the turbine is located. For offshore WTs, we need to add the substructure, which can be fixed or floating depending on whether the WT is placed in shallow, moderately deep, or deep waters. On this basis, two different basic geometric entities can be identified from which all the discretized surfaces of a WT farm will be built. These are a hole plate (called GO₁) and a skinned surface (called GO₂) (see Fig. 1), which are obtained by following different geometric construction procedures. In view of this, and according to the geometric modeling spirit of the current work, it is necessary to make a distinction between them.

Figure 1 Basic geometric entities.

Here, both objects, GO₁ and GO₂, are built from the very beginning as discretized surfaces by using quadrilateral elements QEs (also called cells, panels, or boundary elements). Such simple geometrical elements were chosen because most of VLM- or PM-based codes rely on QE discretizations to represent the lifting and non-lifting surfaces. Next, we present a detailed description of how these objects are built through a discretized setting by using QE.

2.1.1 Object GO₁

This kind of object consists of a rectangular plate together with a circular or ellipsoidal hole. In order to mesh it with QE we make use of the FG-squircular mapping¹, which allows one to smoothly transform a circular domain $\mathcal{D}=\mathit{\left\{}\right(u,v)\in {\mathbb{R}}^{\mathrm{2}}\mid u^{\mathrm{2}}+v^{\mathrm{2}}\le r^{\mathrm{2}}\mathit{\}}$ into a square region $[-r,r]\times [-r,r]$ parameterized as $\mathcal{S}=\mathit{\left\{}\right(x,y)\in {\mathbb{R}}^{\mathrm{2}}\mid x^{\mathrm{2}}+y^{\mathrm{2}}-x^{\mathrm{2}}{y}^{\mathrm{2}}\le r^{\mathrm{2}}\mathit{\}}$ (Fernandez-Guasti, 1992). Knowing that 𝒟 and 𝒮 can be represented as a set of concentric circles and shrunken FG-squircles, we can establish a correspondence between the r disk and the 2r square region by mapping every circular contour in the interior of the disk to a squircular contour in the interior of the square (see Fig. 2).

Figure 2 Geometric object GO₁ – squircle shapes.

According to Fong (2021), an intermediate shape between a circle and a square can be represented by the following implicit equation:

$$\begin{array}{}\text{(1)}& x^{\mathrm{2}}+y^{\mathrm{2}}-{\displaystyle \frac{s^{\mathrm{2}}}{r^{\mathrm{2}}}}{x}^{\mathrm{2}}{y}^{\mathrm{2}}=r^{\mathrm{2}},\end{array}$$

where (x,y) is a set of Cartesian coordinates in ℝ², s is the so-called squareness parameter which allows the shape to be interpolated between a circle and a square, and r is the radius of the original circle. While s=0 generates a circle of radius r, s=1 produces a square with a side length of 2r. In turn, Eq. (1) can be slightly modified to allow an ellipse to be smoothly transformed into a rectangle. Such a representation is given as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{x^{\mathrm{2}}}{r_x^{\mathrm{2}}}}+{\displaystyle \frac{y^{\mathrm{2}}}{r_y^{\mathrm{2}}}}-{\displaystyle \frac{s^{\mathrm{2}}}{r_x^{\mathrm{2}}{r}_y^{\mathrm{2}}}}{x}^{\mathrm{2}}{y}^{\mathrm{2}}=\mathrm{1}.\end{array}$$

When s=0, we obtain the equation of an ellipse with semi-axes r_x, r_y; when s=1 we obtain the equation for a rectangle with sides 2r_x, 2r_y. In order to facilitate the computation of intermediate shapes, Eq. (2) is recast in parametric form through mapping $\mathrm{\Phi }\left(\mathit{\theta }\right):{\mathbb{R}}_{\mathit{\theta }}⟶(x,y)\in {\mathbb{R}}^{\mathrm{2}}$ as

$$\begin{array}{}\text{(3)}& \mathrm{\Phi }\left(\mathit{\theta }\right)=\left\{\begin{array}{c}x\left(\mathit{\theta }\right)\\ y\left(\mathit{\theta }\right)\end{array}\right\}=\left\{\begin{array}{c}\frac{r_x\text{sgn}\left(\cos \mathit{\theta }\right)}{s\sqrt{\mathrm{2}}\left|\sin \mathit{\theta }\right|}\sqrt{\mathrm{1}-\sqrt{\mathrm{1}-s^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathrm{2}\mathit{\theta }}}\\ \frac{r_y\text{sgn}\left(\sin \mathit{\theta }\right)}{s\sqrt{\mathrm{2}}\left|\cos \mathit{\theta }\right|}\sqrt{\mathrm{1}-\sqrt{\mathrm{1}-s^{\mathrm{2}}{\sin }^{\mathrm{2}}\mathrm{2}\mathit{\theta }}}\end{array}\right\},\end{array}$$

with $\mathit{\theta }\in {\mathbb{R}}_{\mathit{\theta }}=\mathit{\{}\mathit{\theta }\mid \mathrm{0}\le \mathit{\theta }\le \mathrm{2}\mathit{\pi }\mathit{\}}\setminus F_{\mathit{\theta }}$, $F_{\mathit{\theta }}=\mathit{\{}\mathrm{0},\frac{\mathit{\pi }}{\mathrm{2}},\mathit{\pi },\frac{\mathrm{3}}{\mathrm{2}}\mathit{\pi },\mathrm{2}\mathit{\pi }\mathit{\}}$ and provided that s≠0. On the one hand, when s=0, the shape corresponds to a circle or ellipse and therefore Eq. (3) reduces to the well-known parametric equation of a circle–ellipse. On the other hand, for θ∈F_θ, Eq. (3) becomes indeterminate and an alternative expression must be used. To avoid problems of a numerical nature, we define an interval $[\mathit{\theta }-\mathrm{TOL},\mathit{\theta }+\mathrm{TOL}]$ around each element of F_θ, where TOL is the allowed tolerance. A deeper look at the FG-Squircular mapping allows us to recognize that, for given r_x and r_y, the values of θ∈F_θ generate points on both the rectangle and its inscribed ellipse (see Fig. 2). Therefore, we can again use the parametric equations of the ellipse to map F_θ.

The object GO₁ is geometrically decomposed into a finite set of quadrilateral cells ${\mathcal{A}}_{{\mathrm{GO}}_{\mathrm{1}}}=\mathit{\left\{}{B}_k\mathit{\right\}}$ as

$${\mathcal{A}}_{{\mathrm{GO}}_{\mathrm{1}}}=\bigcup _{k\in E_{\mathrm{1}}}{B}_k,$$

where $E_{\mathrm{1}}=\mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },(N_r-\mathrm{1}\left)\right(N_c-\mathrm{1}\left)\mathit{\right\}}$, N_r is the number of intermediate shapes (including both inner and outer contours), and N_c is the number of elements along their tangential direction. The total number of cells (or panels) is then given by the cardinal of E₁, i.e., card(E₁). In Fig. 3 (left) we present a schematic of how the division along the radial and tangential directions are performed. In Fig. 3 (right) we show the final mesh of a typical GO₁ for N_r=4 and N_c=17, thus giving card(E₁)=48 panels.

Figure 3 GO₁ – (a) discretization procedure, (b) typical mesh.

2.1.2 Object GO₂

Surface generation in the context of computer-aided design (CAD) is typically done by using lofting or skinning processes. Although skinning can be understood as a sort of lofting (Ball, 1993), some differences have been introduced over the years. However, both processes are intended for passing a surface through a set of so-called cross-sectional curves. In this work, we adopted a specific skinning procedure, hereafter referred to as ruled skinning.

Figure 4 (a) Skinned body (C-S_k stands for cross-section k), (b) ruled surface, (c) representation of a skinned body as a set of ruled surfaces, and (d) typical mesh of a GO₂.

Ruled skinning provides the ability to skin a series of three or more profiles by placing ruled surfaces in between each section of profiles (see Fig. 4a). In turn, a ruled surface is defined by the property that through every point in the surface, there is at least one straight line which also lies in the surface. We can define a ruled surface more formally as a two-dimensional differentiable manifold constructed as the union of one parametric family of lines.

Definition 2.1: ruled surface

The following three definitions of a ruled surface are equivalent (Biran, 2019):

1. a surface such that through each point of it passes a straight line that is fully contained in the surface

2. a surface generated by the motion of a straight line

3. the set of a family of straight lines depending on a parameter that spans a set of real numbers.

Mathematically, a ruled surface can be described by

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\mathit{R}(u,v)& ={\mathit{C}}_{\mathrm{1}}\left(u\right)+v\mathit{r}\left(u\right),v\in \mathbb{R}\\ & =(\mathrm{1}-v){\mathit{C}}_{\mathrm{1}}\left(u\right)+v{\mathit{C}}_{\mathrm{2}}\left(u\right),\end{array}\end{array}$$

where ${\mathit{C}}_k:\mathbb{R}⟶{\mathbb{R}}^{\mathrm{3}}$ is a parameterization for the curve 𝒞_k⊂ℝ³. Any curve R(u₀,v) with fixed parameter u₀ is a generator line, the curve 𝒞_k is a directrix of the representation, and the vectors r(u)≠0 describe the directions of the generators. Alternatively, we can generate a ruled surface by starting with two non-intersecting curves C₁(u) and C₂(u) as directrices and get the line directions as $\mathit{r}\left(u\right)={\mathit{C}}_{\mathrm{2}}\left(u\right)-{\mathit{C}}_{\mathrm{1}}\left(u\right)$ (see Fig. 4b).

In the context of wind turbines, some components (hub, nacelle, blades), although easily generated by a skinning process for modeling purposes, they cannot be represented by ruled surfaces. As an example, let us consider the surface ℛ shown in Fig. 4c. It is not a ruled surface if considered as a whole (first definition in Definition 2.1). However, ℛ can be obtained as the union of three ruled surfaces. Furthermore, every GO₂ object is geometrically decomposed into a finite set of quadrilateral cells ${\mathcal{A}}_{{\mathrm{GO}}_{\mathrm{2}}}=\mathit{\left\{}{B}_k\mathit{\right\}}$ as follows:

$$\begin{array}{}\text{(5)}& {\mathcal{A}}_{{\mathrm{GO}}_{\mathrm{2}}}=\bigcup _{i\in R_{\mathrm{S}}}{\mathcal{A}}_{{\mathrm{GO}}_{\mathrm{2}}}^i,\text{and}{\mathcal{A}}_{{\mathrm{GO}}_{\mathrm{2}}}^i=\bigcup _{k\in E_{\mathrm{2},\mathrm{1}}}{B}_k,\end{array}$$

where $R_{\mathrm{S}}=\mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{rs}}\mathit{\}}$, N_rs is the number of ruled surfaces in what ${\mathcal{A}}_{{\mathrm{GO}}_{\mathrm{2}}}$ is decomposed, and E_2,i is a finite subset of ℕ. Then, the total number of QE used to discretize ${\mathcal{A}}_{{\mathrm{GO}}_{\mathrm{2}}}$ is calculated as ${\sum }_{i=\mathrm{1}}^{N_{\mathrm{rs}}}\text{card}\left(E_{\mathrm{2},i}\right)$. In Fig. 4d we present the mesh of a typical hub nose, represented by N_rs=9 and $\text{card}\left(E_{\mathrm{2},i}\right)=\mathrm{24}$, thus giving a total of 216 panels.

Finally, it should be stressed that any pair of cells belonging to GO₁ or GO₂ must meet the following requirements:

• If B_k∩B_j for k≠j is exactly one point, then it is a common vertex (node) of B_k and B_j.

• If B_k∩B_j for k≠j is a line, then it is a common facet of B_k and B_j (edge in two dimensions).

2.2 Wind turbine farm modeling

Here, we present a detailed description on how the surface of each component in a wind turbine farm is modeled in terms of the geometric objects already introduced in Sect. 2.1. Figure 5 provides a summary on which type of geometric object (GO₁ or GO₂) is involved in generating the discretized surface associated with each component of a wind turbine. From now on, we refer to mesh (or bound-vortex lattice) to make reference to a discretized surface.

Figure 5 Geometric modeling of a wind turbine: IEA Wind 15 MW offshore reference wind turbine.

Each wind turbine can be composed, at most, of six different components:

• Blades are responsible for capturing part of the energy available in the wind.

• The hub supports the blades and houses the pitch system.

• The nacelle houses the components that convert mechanical energy into electrical energy.

• The tower gives height to the rotor and supports the mass of the nacelle, hub, and blades.

• The ground represents either the terrain for onshore wind turbines or the sea surface for offshore wind turbines.

• The monopile holds the tower and the rest of wind turbine components above the sea floor.

The user can generate different wind turbine configurations by turning any of these components on or off. In the following subsections we will discuss in detail how to model each component and what parameters should be provided to generate them. All these parameters are specified through some options available in the main mesher script, as well as through a configuration ASCII file. This topic will be covered to some extent in Sect. 3 on computational implementation.

Table 1 Standard parameters of a wind turbine.

As usual, some of those parameters are related to the global configuration of the wind turbine; these are listed in Table 1. To complete the turbine setup, it is also required to indicate which components should be considered to build the wind turbine mesh (see Table 2).

Table 2 Wind turbine components.

Clearly, the number of blades cannot be increased indiscriminately. If this happens, there may be geometric interference between the hub surface and blade roots. Furthermore, as the number of blades is becoming larger and larger, it may happen that two blades are very close to each other, and therefore some geometrical interference may arise between them.

2.2.1 The tower and monopile

The tower and monopile components are essentially of the same kind. Both of them are ruled surfaces and therefore generated by GO₂-like entities. As input data are necessary to generate the aerodynamic mesh associated with the tower and monopile, we need to provide the diameter at the ends of the component, its length, and number of aerodynamic nodes along its longitudinal and tangential directions. Tables A1 and A2 in Appendix A provide a summary of the variables associated with both components.

2.2.2 The nacelle

The aerodynamic mesh of the nacelle 𝒩 is generated as the union of four sub-meshes: (i) the tower–nacelle connector ${\mathcal{N}}_{{\mathrm{GO}}_{\mathrm{2}}}^{\mathrm{1}}$, (ii) a curved patch ${\mathcal{N}}_{{\mathrm{GO}}_{\mathrm{1}}}^{\mathrm{2}}$, (iii) a cylindrical surface ${\mathcal{N}}_{{\mathrm{GO}}_{\mathrm{2}}}^{\mathrm{3}}$, and iv) the tail of the nacelle ${\mathcal{N}}_{{\mathrm{GO}}_{\mathrm{2}}}^{\mathrm{4}}$. All nacelle components are generated by GO₂-like entities with the exception of the curved patch, which is of type GO₁. Figure 6 shows an exploded schematic of a typical nacelle of a wind turbine.

Figure 6 Exploded view of a typical wind turbine nacelle.

Among the several parameters involved in the design of wind turbines, the cone angle, tilt angle, and pitch angle are directly related with the aerodynamic behavior of the rotor. In particular, the tilt angle is used to provide sufficient clearance between the rotor blades and the tower. Here, for meshing purposes, such tilt angle γ is defined as the angle between the longitudinal axis of the nacelle and the horizontal plane. This definition implies that the nacelle axis is, in general, not orthogonal to the longitudinal axis of the tower and, therefore, not orthogonal to the axis of the tower–nacelle connecting piece ${\mathcal{N}}_{{\mathrm{GO}}_{\mathrm{2}}}^{\mathrm{1}}$ either. In light of the above, the most complicated step to generate the aerodynamic mesh of the nacelle lies in the connection between the tower–nacelle connector and the nacelle itself, which is done through the curved patch ${\mathcal{N}}_{{\mathrm{GO}}_{\mathrm{1}}}^{\mathrm{2}}$. Such an object is obtained by following the next two steps in order:

1. Generate a typical flat GO₁-like object with appropriate dimensions.

2. Deform it into a curved surface as shown in Fig. 6. In other words, it means rolling the patch over a “virtual” cylinder representing the nacelle (see Fig. 7).

For the second step, we consider a continuous deformation 𝒯_n that maps each point r in a connected subset 𝒟⊂ℝ² to a point R=𝒯_n(r) on a surface in the three-dimensional space ℝ³. Mathematically, ${\mathcal{T}}_{\text{n}}:\mathcal{D}\subset {\mathbb{R}}^{\mathrm{2}}⟶{\mathbb{R}}^{\mathrm{3}}$ is given by

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\mathcal{T}}_{\text{n}}(x,y)& =x{\widehat{\mathit{E}}}_{\mathrm{1}}+R_{\mathrm{N}\mathrm{1}}\sin \left({\displaystyle \frac{\mathrm{1}}{R_{\mathrm{N}\mathrm{1}}}}y\right){\widehat{\mathit{E}}}_{\mathrm{2}}\\ & +R_{\mathrm{N}\mathrm{1}}\cos \left({\displaystyle \frac{\mathrm{1}}{R_{\mathrm{N}\mathrm{1}}}}y\right){\widehat{\mathit{E}}}_{\mathrm{3}},\end{array}\end{array}$$

where R_N1 is the radius of the cylinder representing the nacelle. It is straightforward to prove that the deformation map 𝒯_n is an isometric map. Therefore, it preserves the length of every possible arc of material points on the parametric domain 𝒟. This is true if and only if the Jacobian D𝒯_n on 𝒟 preserves the lengths of vectors in ℝ² in the sense that $\left|\left|D{\mathcal{T}}_{\text{n}}\right(\mathit{r})\right||=\left||\mathit{r}\right||$ for each r∈ℝ². Equivalently, D𝒯_n must obey ${\left(D{\mathcal{T}}_{\text{n}}\right)}^{T}D{\mathcal{T}}_{\text{n}}={\mathit{I}}_{\mathrm{2}}$, with I₂ being the identity linear transformation on ℝ² (Chen et al., 2018). Once the flat patch has been deformed into a curved patch, an affine transformation (translation and/or rotation) ${\mathcal{A}}_{\text{n}}:{\mathbb{R}}^{\mathrm{3}}⟶{\mathbb{R}}^{\mathrm{3}}$ is applied to assemble it with the rest of the nacelle. Thereby, the ${\mathcal{N}}_{{\mathrm{GO}}_{\mathrm{1}}}^{\mathrm{2}}$ object is obtained by means of the following map composition:

$$\begin{array}{}\text{(7)}& {\mathcal{N}}_{{\mathrm{GO}}_{\mathrm{1}}}^{\mathrm{2}}={\mathcal{A}}_{\text{n}}\circ {\mathcal{T}}_{\text{n}}\left(\mathcal{D}\right).\end{array}$$

Finally, we need to address the non-smooth connection between the curved patch and the tower–nacelle connector. Such a non-smoothness arises as a consequence of two reasons: (i) the deformation of the original circular hollow cavity (flat patch) into an elliptical cavity (curved patch) and (ii) the misalignment between the tower longitudinal axis and the hollow axis of the curved patch. For this purpose, we compute in advance the radii of an elliptical hollow in the original flat patch that will ensure a smooth connection between the tower and the deformed curved patch ${\mathcal{N}}_{{\mathrm{GO}}_{\mathrm{1}}}^{\mathrm{2}}$ (see Fig. 7). Such radii are given by

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{r}_x& ={\displaystyle \frac{R_{\mathrm{N}\mathrm{3}}}{\cos \mathit{\gamma }}},\\ r_y& =R_{\mathrm{N}\mathrm{1}}{\sin }^{-\mathrm{1}}\left({\displaystyle \frac{R_{\mathrm{N}\mathrm{3}}}{R_{\mathrm{N}\mathrm{1}}}}\right),\end{array}\end{array}$$

where R_N3 is the tower–nacelle connector radius (see Table A3). When γ=0, the radius r_x in Eq. (7) is directly the radius of the upper part of the tower (no tilt angle). However, the radius r_y is different from the upper tower radius because of the deformation of the flat patch. Table A3 in Appendix A provides a summary of the variables associated with the nacelle.

Figure 7 Exploded view of a typical wind turbine nacelle.

2.2.3 The hub

The aerodynamic mesh of the hub ℋ is generated as the union of a series of sub-meshes depending on the number of blades of the wind turbine: (i) blade–hub connectors ${\mathcal{H}}_{{\mathrm{GO}}_{\mathrm{2}}}^k$ for $k=\mathrm{1},\mathrm{\dots },N_{\mathrm{B}}$, (ii) curved patches ${\mathcal{H}}_{{\mathrm{GO}}_{\mathrm{1}}}^{N_{\mathrm{B}}+k}$ for $k=\mathrm{1},\mathrm{\dots },N_{\mathrm{B}}$, and (iii) the nose of the hub ${\mathcal{H}}_{{\mathrm{GO}}_{\mathrm{2}}}^{\mathrm{2}{N}_{\mathrm{B}}+\mathrm{1}}$. All hub components are generated by GO₂-like entities with the exception of the curved patches, which are of type GO₁. Figure 8 shows an exploded schematic of a typical hub for a three-blade wind turbine.

Figure 8 Exploded view of a typical three-blade wind turbine hub.

As before, the most complicated step to generate the aerodynamic mesh of the hub lies in the connection between the blade–hub connector and the hub itself, which is done through the curved patches. Those objects are generated by following a similar procedure as the nacelle:

1. Generate a typical flat GO₁-like object with appropriate dimensions.

2. Deform it into a curved surface (see Fig. 7).

3. Make N_B copies of this object and place them properly to assemble the rotor.

The deformation of a flat GO₁-like object into a curved one resembling a part of the hub's surface is performed by means of a continuous isometric mapping ${\mathcal{T}}_{\text{h}}:\mathcal{D}\subset {\mathbb{R}}^{\mathrm{2}}⟶{\mathbb{R}}^{\mathrm{3}}$ as follows:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\mathcal{T}}_{\text{h}}(x,y)& =x{\widehat{\mathit{E}}}_{\mathrm{1}}+R_{\mathrm{H}\mathrm{2}}\sin \left({\displaystyle \frac{\mathrm{1}}{R_{\mathrm{H}\mathrm{2}}}}y\right){\widehat{\mathit{E}}}_{\mathrm{2}}\\ & +R_{\mathrm{H}\mathrm{2}}\cos \left({\displaystyle \frac{\mathrm{1}}{R_{\mathrm{H}\mathrm{2}}}}y\right){\widehat{\mathit{E}}}_{\mathrm{3}},\end{array}\end{array}$$

where R_H2 is the radius of the cylinder representing the hub (without considering the blade–hub connectors). Once the flat patch has been deformed into a curved patch, we generated as many copies of the curved patch as the number of blades in the wind turbine. Then, an affine transformation ${\mathcal{A}}_{\text{h}}^k:{\mathbb{R}}^{\mathrm{3}}⟶{\mathbb{R}}^{\mathrm{3}}$ for $k=\mathrm{1},\mathrm{\dots },N_{\mathrm{B}}$ is applied to assemble each curved patch with the rest of the hub. Thereby, the ${\mathcal{H}}_{{\mathrm{GO}}_{\mathrm{1}}}^k$ object is obtained by means of the following map composition:

$$\begin{array}{}\text{(10)}& {\mathcal{H}}_{{\mathrm{GO}}_{\mathrm{1}}}^k={\mathcal{A}}_{\text{h}}^k\circ {\mathcal{T}}_{\text{h}}\left(\mathcal{D}\right),\text{for}k=\mathrm{1},\mathrm{\dots },N_{\mathrm{B}}.\end{array}$$

It should be noted that the dimension of the flat patches ${\mathcal{H}}_{{\mathrm{GO}}_{\mathrm{1}}}^k$ (before deformation) along the y coordinate (see Fig. 8) depends on the number of blades and the radius of the hub, i.e., $W_{\mathrm{Hy}}=\mathrm{2}{R}_{\mathrm{H}\mathrm{1}}\mathit{\pi }/N_{\mathrm{B}}$. Finally, the radii of the elliptical hollow in the original flat patch that will ensure a smooth connection between the blade–connector and the deformed curved patch ${\mathcal{H}}_{{\mathrm{GO}}_{\mathrm{1}}}^k$ are given by

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{r}_x& ={\displaystyle \frac{R_{\mathrm{H}\mathrm{3}}}{\cos \mathit{\beta }}},\\ r_y& =R_{\mathrm{H}\mathrm{2}}{\sin }^{-\mathrm{1}}\left({\displaystyle \frac{R_{\mathrm{H}\mathrm{3}}}{R_{\mathrm{H}\mathrm{2}}}}\right),\end{array}\end{array}$$

where β is the cone angle of the rotor, and R_H3 is the blade–hub connector radius. When β=0, the radius r_x in 11 is directly the radius of the blade–hub connector (no cone angle). Table A4 in Appendix A provides a summary of the variables associated with the hub.

2.2.4 The blade

The aerodynamic mesh of a blade ℬ is generated as the union of two sub-meshes: (i) the blade root ${\mathcal{B}}_{{\mathrm{GO}}_{\mathrm{2}}}^{\mathrm{1}}$ and (ii) the lifting surface of the blade ${\mathcal{B}}_{{\mathrm{GO}}_{\mathrm{2}}}^{\mathrm{2}}$. All blade components are generated by GO₂-like entities.

The mesh generation for wind turbine blades is a non-trivial process because it involves discretizing a three-dimensional surface whose shape changes along its longitudinal axis. Table 3 lists the geometric data that must be provided, as a function of the longitudinal coordinate, to build a three-dimensional wind turbine blade. To simplify what follows and avoid falling into excessive formalism, let us define B₀ as the set of design parameters of the blade. Figure 9 presents an example of how such parameters are usually defined along the longitudinal axis of the blade. A good representation of its surface necessarily requires specifying these parameters at “several” points along the blade.

Table 3 Blade geometry definition: user-input parameters.

Figure 9 Parameter definition along the blade.

On this basis, the blade is divided into N_B2−1 non-uniform intervals such that $[\mathrm{0},L_B]={\cup }_{i=\mathrm{1}}^{N_{\mathrm{B}\mathrm{2}}-\mathrm{1}}[z_i,z_{i+\mathrm{1}}]$ such that $z_{i+\mathrm{1}}>z_i$, where L_B is length of the blade, z₁=0, $z_{{\mathrm{N}}_{\mathrm{B}\mathrm{2}}}=L_B$, and N_B2 is the number of nodes along the blade (see Fig. 9). The geometric parameters at each coordinate z_i are interpolated, according to the provided data B₀, by using cubic splines. In the case of pre-bend and pre-sweep, we have two options, they are either provided by the manufacturer or they can be included by using Zuteck's formula (Larwood et al., 2014). The following equation is used to define the sweep and pre-bend:

$$\begin{array}{}\text{(12)}& x=x_{\mathrm{tip}}{\left({\displaystyle \frac{z-z_{\mathrm{0}}}{L_B-z_{\mathrm{0}}}}\right)}^{{\mathit{\xi }}_{\mathrm{1}}}\text{and}y=y_{\mathrm{tip}}{\left({\displaystyle \frac{z-z_{\mathrm{0}}}{L_B-z_{\mathrm{0}}}}\right)}^{{\mathit{\xi }}_{\mathrm{2}}},\end{array}$$

where y and x are the local distance from the elastic axis (or pitch axis) to the sweep/pre-bend curve, y_tip and x_tip are the distance from the pitch axis to the sweep/pre-bend curve at the blade tip, z is the local distance along the blade measured from the blade root, z₀ is the position of the beginning of the blade sweep/pre-bend, and ξ is the sweep/pre-bend exponent (see Fig. 10a). As mentioned above, the pre-bend and pre-sweep can be either automatically generated by using Zuteck's formula or manually specified by providing $x_k^s$ and $y_k^s$ at every station k.

Figure 10 (a) Parameters that define the wind turbine blade sweep and pre-bend and (b) sweep and pre-bend blade deformations.

The arc length for all curved blade shapes is equal to the original length of the straight blade, but with a slightly smaller rotor radius. The arc length of the blade should be kept the same to avoid blade extension, which will bias results towards longer blades that produce more power. To this end, we consider the position of an arbitrary point on the elastic axis of the blade at the reference configuration to be given by ${\mathit{r}}_{\mathrm{0}}=z{\widehat{\mathit{E}}}_{\mathrm{3}}$, while the position of the same point at the deformed configuration can be expressed as follows:

$$\begin{array}{}\text{(13)}& \mathit{r}=x\left(z\right){\widehat{\mathit{E}}}_{\mathrm{1}}+y\left(z\right){\widehat{\mathit{E}}}_{\mathrm{2}}+z{\widehat{\mathit{E}}}_{\mathrm{3}}+\mathit{u},\end{array}$$

where u is the displacement vector. For a given coordinate z along the blade, the shortening in the axial direction, u₃, due to the arc-length conservation can then be written as

$$\begin{array}{}\text{(14)}& \begin{array}{rl}s\left(z\right)& ={\int }_{\mathrm{0}}^z\left|\left|{\mathit{r}}^{\prime }\right(\mathit{\eta })\right||\mathrm{d}\mathit{\eta },\\ & ={\int }_{\mathrm{0}}^z\sqrt{x^{\prime \mathrm{2}}\left(\mathit{\eta }\right)+y^{\prime \mathrm{2}}\left(\mathit{\eta }\right)+(\mathrm{1}+u_{\mathrm{3}}^{\prime }(\mathit{\eta }){)}^{\mathrm{2}}}\mathrm{d}\mathit{\eta },\\ u_{\mathrm{3}}\left(z\right)& ={\int }_{\mathrm{0}}^z\left\{\sqrt{\mathrm{1}-\left[x^{\prime \mathrm{2}}\left(\mathit{\eta }\right)+y^{\prime \mathrm{2}}\left(\mathit{\eta }\right)\right]}-\mathrm{1}\right\}\mathrm{d}\mathit{\eta },\end{array}\end{array}$$

where η is a dummy integration variable, s is the arc-length coordinate along the curved blade, and $(\cdot {)}^{\prime }$ stands for derivative with respect to η. To obtain the shortening at a given blade section z (last expression in Eq. 14), we enforced the arc-length conservation by imposing that the length of the deformed blade must be equal to the original straight blade length, i.e., s(z)=z. The blade root is represented as z=0 in Eq. (14). The variables x and y represent the local distance from the pitch axis to the sweep and pre-bend curve, respectively. Then, the new nodal z coordinates associated with the initial partition [0,L_B] are obtained as $z_i^d=z_i-u_{\mathrm{3}}\left(z_i\right)$.

Once the new z coordinates are obtained, we perform a sanity analysis on all curved blades to check that the arc length is close enough to the length of the straight blade. Such an analysis is carried out by computing the difference between the original blade length and the curved blade length:

$$\begin{array}{}\text{(15)}& L_B^d={\int }_{\mathrm{0}}^{z_{\mathrm{c}}}\sqrt{\mathrm{1}+x^{\prime \mathrm{2}}\left(\mathit{\eta }\right)+y^{\prime \mathrm{2}}\left(\mathit{\eta }\right)}\mathrm{d}\mathit{\eta },\end{array}$$

where z_c is the tip radius of the curved blade. A script is used to automatically calculate z_c, perform the arc length sanity check, and prepare the blade geometry. In Fig. 10b, we present some examples for sweep and pre-bend blade configurations. Table A5 in Appendix A provides a summary of the variables associated with the blade.

2.2.5 The ground

The aerodynamic mesh of the ground 𝒢 is generated using only a GO₁-like object. As input data we need to provide (i) the ground–tower connection radius, (ii) the extent of the ground, namely, the side length of the square representing the ground area, and (iii) the number of aerodynamic nodes along its radial and tangential directions. Additionally, it is possible to generate an uneven ground by providing a user-defined function to compute the ground elevation or a scattered data set representing the ground elevation. The last option requires a fitting procedure (regression) to obtain the elevation on the ground aerodynamic mesh. This feature will be discussed in more detail in the computational implementation section. Table A6 in Appendix A provides a summary of the variables associated with the ground.

2.2.6 Wind turbine assembling

Once all the components for a given wind turbine configuration are generated, the next step is to assemble them to obtain the complete wind turbine mesh. The code internally calculates all the data required to assemble each turbine within the wind farm. The only information to be provided by the user is the initial rotor angle, the initial yaw angle, and initial pitch angles. Table 4 lists the data that must be supplied by the user to assemble the turbine. These data are located at the end of each turbine data sheet.

Table 4 Assembling parameters.

It should be stressed that the number of pitch angles to be provided must match the number of blades in the rotor, i.e., ${\mathit{\theta }}_{\mathrm{p}}\in {\mathbb{R}}^{N_{\mathrm{B}}}$, and their values can be different from each other. As an example, we present below how the initial angles must be specified to configure a three-blade rotor.

0 % Yaw0WT - Initial yaw angle [degree]
0 % Rot0WT - Initial rotor angle [degree]
4,4,4 % Pitch0WT - Initial pitch angle [degree]

3 Computational implementation

In this section we describe the main aspects associated with the computational implementation of the mesh generator code called UVLMeshGen developed at the University of Bergen (Norway) in collaboration with Universidad Nacional de Río Cuarto (Argentina). The full code is available in a GitHub repository (Roccia, 2023) for public access under a Creative Commons Attribution 4.0 International License.

This meshing tool is fully implemented in MATLAB^®, and it is intended for generating meshes for onshore and offshore wind farms consisting of horizontal-axis wind turbines. Among the main features of UVLMeshGen, we can mention the following:

• wind farm meshing,

• terrain meshing,

• specification of terrain topography,

• kinematic processor, and

• exporting files (Tecplot, UVLM, …).

UVLMeshGen provides meshing data into a series of structure variables by using typical FEM-oriented data arrays, such as nodal coordinate arrays and connectivity arrays. Next, we describe to some extent the organization of the code, description of the main script, important variables, and the features introduced above.

3.1 Code structure

The mesh generator is designed by following a procedural programming paradigm. The software contains two main blocks: (i) the geometric processor, which is responsible for generating the mesh of each wind turbine component and its assembling, and (ii) the kinematic processor, which is in charge for computing the kinematics of lifting and non-lifting surfaces. The exporting module, which is responsible for writing output files (Tecplot, input data for UVLM-based codes, etc.), can be added/modified by the user according to their needs. However, the code incorporates an option to export meshes in Tecplot format by default. In Fig. 11 we present a flowchart of UVLMeshGen.

Figure 11 UVLMeshGen flowchart.

The input files needed by UVLMeshGen are ASCII files and user-defined MATLAB functions. To keep all the input files clean and separate from the main code, all of them are organized in several folders.

The code starts by loading as many input files as different wind turbines the wind farm contains. These files must be located inside the folder “WTDataSheetFiles/”, and they gather the geometric definition of each wind turbine. Then, the code generates the meshes associated with each wind turbine component. Among them, the blade requires an additional input file containing its complete setup (see Table 3), which must be located inside the folder “BladeDataSheetFiles/”. The blade configuration file also makes reference to the airfoil distribution along its spanwise direction, thus requiring further information about airfoil coordinates. These data are stored inside the folder “Airfoils/”.

Once all the components for each wind turbine have been generated, the code continues with the wind farm terrain, depending on whether it is activated or not in the main script. Then, all wind turbine parts are assembled to shape the final wind farm. The aforementioned processes, namely the generation of wind turbine components, terrain and their final assembly into a wind park, make up the so-called geometric processor.

UVLMeshGen also has a kinematics module that is in charge of generating pitch, yaw, and rotor motions. This module is general enough to allow any kinematics to be prescribed by the user through custom user-defined functions. Such a feature is very useful to investigate the aerodynamic performance of wind turbines. To write the outputs, users can either select some of the file formats included by default or provide their own scripts to export data.

3.2 Main script

The code is executed through a main script, which contains a set of general options, such as wind farm layout, terrain elevation, wind turbine kinematics, and output files. In Table 5 we present a brief description of the variables to be specified in the main script.

Table 5 Main script variables – general options.

WTNames is a cell array variable of dimension N_WT×4. The first column contains the file names associated with the configuration of each wind turbine in the wind farm. The last three columns contain the $(x,y,z)$ coordinates of each wind turbine. As an example, let us consider a wind farm made up of four different wind turbines (N_WT=4). Under these assumptions, the variable WTNames could take the following values:

WTNames = {'DataSheet_DTU_10MW.DAT', 0.0, 0.0, 0.0
           'DataSheet_IEA_15MW.DAT', -100.0, 100.0, 0.0
           'DataSheet_Sandia01.DAT', -120.0, -120.0, 0.0
           'DataSheet_Sandia01.DAT',  100.0, -100.0, 0.0}.

EQWT_FLAG is a string variable which can take only two values: ON or OFF. If the value ON is used, the wind farm will consist of only one type of wind turbine. The setup of said turbine corresponds to the one specified in the first row of the variable WTNames. However the number of rows of WTNames must still be equal to the number of turbines considered in NumWT. This is because the wind turbines will be placed within the wind farm by extracting the coordinates from the variable WTNames.

GroundDivision is a cell array variable of dimension 2×1. The cell 𝙶𝚛𝚘𝚞𝚗𝚍𝙳𝚒𝚟𝚒𝚜𝚒𝚘𝚗{1,1} contains the coordinates of the patches into which the terrain will be divided along the x direction. In a similar way, the cell 𝙶𝚛𝚘𝚞𝚗𝚍𝙳𝚒𝚟𝚒𝚜𝚒𝚘𝚗{2,1} contains the coordinates of the patches into which the terrain will be divided along the y direction. The number of patches along the x and y directions are determined as $\text{dim}\left(\mathtt{\text{ GroundDivision}}\mathit{\{}\cdot ,\mathrm{1}\mathit{\}}\right)-\mathrm{1}$.

To ensure a smooth terrain surface with a regular discretization of elements, it is recommended to consider one patch per turbine in the wind farm. The dimensions of the patch must be large enough to accommodate the ground–tower coupling. Once the surface of the terrain is divided, the code will automatically determine where the turbines are located, it will generate the corresponding ground–tower coupling, and in those patches without turbines, the code will generate a rectangular grid. As an example, let us consider the wind farm introduced above (see variable WTNames). For this wind farm layout, a suitable terrain division could be as follows:

GroundDivision={[180, 60, -60, -180]
                [180, 60, -60, -180]},

from which it is clear that three patches were considered along both directions, i.e., nine patches in total. In Fig. 12 we present a schematic intended to explain how the terrain division is performed and the reference frame used.

Figure 12 Example of a wind farm terrain.

Ground_FLAG is a cell array variable of dimension 1×6. All the cells contained in this array are of type string, and they specify whether or not to generate the wind farm terrain. In addition, this variable allows configuring the terrain elevation and what options will be used to generate it. The reader is referred to Sect. 3.4 for a more detailed description of this feature.

Kinematics_FLAG is a cell array variable of dimension 1×5. All the cells contained in this array are of type string, and they specify whether or not kinematics associated with different wind turbine components are generated (e.g., yaw motion, pitch motion, and/or rotor motion). Furthermore, this variable allows configuring the time step to be used for the kinematics when considering heterogeneous wind farms. The reader is referred to Sect. 3.5 for a more detailed description of this feature.

OPT_FLAG is a cell array variable of dimension 1×2. All the cells contained in this array are of type string, and they specify whether or not the generated meshing data will be exported. This option consists of two fields as follows:

OPT_FLAG={'ON','MyOutput'},

where the first cell supports two values {𝙾𝙽,𝙾𝙵𝙵}, which allow us to indicate whether the data coming from the meshing procedure should be exported or not, and the second cell specifies the name of the user-defined script to be used for exporting purposes. The code has built-in three export options by default: (i) TecplotOutput, which allows us to export the assembled wind farm in Tecplot format, (ii) TecplotKinematics, which allows us to export the assembled wind farm and its kinematics in Tecplot format, and (iii) VLMSim, which allows us to write the input files for the aerodynamic solver VLMSim (Vortex-Lattice Method – Simulation) developed by the Group of Applied Mathematics, Universidad Nacional de Río Cuarto (Verstraete et al., 2023). The second option requires the Kinematics_FLAG to be ON, otherwise an error will occur. All user-defined export scripts must be located inside the folder “Output Scripts/”.

3.3 Output data structure

When UVLMeshGen is executed, it will generate four main structure variables where all meshing data are stored. These variables are WIND_TURBINE, CONNECT, GROUND_FARM, and KINEMATICS. The first variable contains only data associated with mesh coordinates and geometric dimensions related to the entire wind farm. The second structure variable contains all connectivity arrays associated with the entire wind farm discretization. The third variable contains the wind farm kinematics, namely, positions and velocities of the entire wind farm mesh for the stipulated simulation time grid. The last variable contains all data associated with the wind farm terrain.

• The WIND_TURBINE variable is indexed by the number of wind turbines within the wind park, i.e., WIND_TURBINE(i) for $i=\mathrm{1},\mathrm{\dots },N_{\mathrm{WT}}$. In Table B1 in Appendix B, we list the main fields associated with WIND_TURBINE structure. If any component of the wind turbine is disabled during the geometric modeling, the field associated with it is assigned the “empty value”.

• As before, the CONNECT variable is indexed by the number of wind turbine within the wind park. In Table B2 in Appendix B, we list the main fields associated with CONNECT structure.

• The GROUND_FARM variable contains only information regarding the wind farm terrain. There is no kinematics associated with the terrain since it is motionless, so its position is the same over time, and therefore the velocities at the control points are zero for all t. In Table B3 in Appendix B, we list the main fields associated with GROUND_FARM structure. The field GROUND_FARM.PATCH(i) for $i=\mathrm{1},\mathrm{\dots },N_{\mathrm{Gp}}$, where N_Gp is the number of patches into which the terrain was divided, contains nodal coordinates and connectivities, among other data.

• The KINEMATICS variable is indexed by the number of simulation time steps used by the kinematics processor, i.e., KINEMATICS(i) for $i=\mathrm{1},\mathrm{\dots },N_{\mathrm{\Delta }t}$, where N_Δt is the number of time steps. As fields, this variable contains (i)  KINEMATICS(i).GroundFarm.PATCH(k) for $k=\mathrm{1},\mathrm{\dots },N_{\mathrm{Gp}}$ and (ii) KINEMATICS(i).WT(j) for $j=\mathrm{1},\mathrm{\dots },N_{\mathrm{WT}}$. Although the terrain is motionless, KINEMATICS(i).GroundFarm stores the nodal coordinates, control point (CP) coordinates, and CP velocities (which are zero) over time. This decision on the storage of terrain data is based on a potential future implementation of “terrain motion” to simulate the sea surface and the effects of waves on offshore wind turbines. The fields associated with the second structure, KINEMATICS(i).WT(j), are listed in Table B4 in Appendix B.

3.4 Wind farm terrain processor

This module is in charge of generating the wind farm terrain. All necessary parameters are introduced via the Ground_FLAG cell array variable located in the main script. As previously mentioned, this variable consists of six cells as follows:

Ground_FLAG={'ON','ON','userfunction','MyGround',
'GroundData.DAT','poly23'},

where the first and second cells admit two values {𝙾𝙽,𝙾𝙵𝙵} indicating whether the wind farm terrain generation is activated or not and whether the terrain elevation feature is enabled or not, respectively. The third cell supports two different keywords {'𝚞𝚜𝚎𝚛𝚏𝚞𝚗𝚌𝚝𝚒𝚘𝚗','𝚎𝚡𝚝𝚎𝚛𝚗𝚊𝚕𝚍𝚊𝚝𝚊'}, which specify whether the terrain elevation will be generated via some user-defined function or whether it will be generated by fitting a scattered terrain data set. The fourth and fifth cells specify the name of the user-defined function and the name of the file containing the terrain data set, which should be placed inside the “Terrain Data/” folder. The mesher will use either the function or the data set according to the option specified in the third cell. The last cell specifies the fit type to use if externaldata is selected. Internally, the mesher uses the MATLAB^® intrinsic function “fit” to fit a surface to the data provided by the external file. Any fitting option accepted by fit can be specified as a valid option in Ground_FLAG.

The user-defined function for generating the terrain elevation must receive two inputs: (i) an array of dimension N_WF×3 (where WF denotes wind farm) containing the nodal coordinates of the flat wind farm terrain and (ii) the number of nodes of the wind farm terrain, N_WF. The terrain coordinates within the input array are organized as follows: x coordinates (first column), y coordinates (second column), and z coordinates (third column). The user can impose/calculate any elevation profile on the ground (coordinates z) as long as it does not present abrupt changes. As output, the function only needs to provide an array of coordinates of dimension N_WF×3, where the third column contains the z coordinates modified according to the elevation model proposed by the user. As an example, we provide below a very simple user-defined function that allows us to impose a terrain elevation profile.

function C = MyTerrainLevel (XYZ, NWF)}}
Ymin = min(XYZ(:,2));
Ymax = max(XYZ(:,2));
DD   = Ymax - Ymin;
C = XYZ;

A = [Ymin^2 Ymin 1; Ymax^2 Ymax 1; 2*Ymin 1 0];
F = [0; DD/6; 0];
Coef = A\textbackslash F;

for i = 1:NWF
 C(i,3) = Coef(1)*C(i,2)^2 + Coef(2)*C(i,2)^2 + Coef(3);
end

Figure 13 shows the elevation of the ground surface generated by using a parabolic and a sinusoidal profile. The ground surface corresponding to the parabolic profile was generated using the user-defined function introduced above. As can be seen, the elevation profile does not present any type of abrupt changes or jumps.

Figure 13 Example of wind farm terrain elevations.

When the option 'externaldata' is chosen, an external ASCII file containing $(x,y,z)$ coordinates of points on the ground must be provided. The amount/quality of the points considered should be large/good enough to ensure that the fitted surface renders a realistic elevation profile. Once the fitting process is done, the resulting surface equation will be used to find the elevation for a given nodal terrain coordinate (x,y). The data must be provided in three columns: x coordinates (first column), y coordinates (second column), and z coordinates (third column).

3.5 Kinematic processor

This module is in charge of generating the kinematics for the entire wind farm. As before, all necessary parameters are introduced via the Kinematic_FLAG cell array variable located in the main script. As previously mentioned, this variable consists of five cells as follows:

Kinematic_FLAG={'ON','RotorOFF','YawOFF','PitchOFF','min'},

where the first cell admits two values {'𝙾𝙽','𝙾𝙵𝙵'} indicating whether the wind farm kinematics is enabled or not. The second through fourth cells allow us to specify whether rotor, yaw, and pitch kinematics are enabled or not. Each of the following cells admits two values: {'𝚁𝚘𝚝𝚘𝚛𝙾𝙽','𝚁𝚘𝚝𝚘𝚛𝙾𝙵𝙵'}, {'𝚈𝚊𝚠𝙾𝙽','𝚈𝚊𝚠𝙾𝙵𝙵'}, and {'𝙿𝚒𝚝𝚌𝚑𝙾𝙽','𝙿𝚒𝚝𝚌𝚑𝙾𝙵𝙵'}.

When considering a wind farm, it may happen that it is composed of different types of turbines (heterogeneous wind farm), which can in turn result in very different aerodynamic grids (e.g., large differences in the size of aerodynamic panels). In UVLM-based codes it is customary to define characteristic magnitudes for computing force coefficients. Typically they are the characteristic density ρ_C, the characteristic length L_C, the characteristic velocity V_C, and the characteristic time T_C, which is obtained as $T_{\mathrm{C}}=L_{\mathrm{C}}/V_{\mathrm{C}}$. UVLMeshGen offers two options for setting the characteristic length. This can be either provided by the user or computed by using a default internal procedure. If the automatic option is selected, L_C is calculated based only on the discretization of the blade lifting surface as follows:

$$\begin{array}{}\text{(16)}& L_{\mathrm{C}}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^{N_{\mathrm{LS}}}{A}_i}{N_{\mathrm{LS}}}}},\end{array}$$

where N_LS is the number of panels on the blade lifting surface, and A_i is the surface area of the ith panel belonging to the lifting surface. Regarding the characteristic velocity, it is usually set to be the magnitude of the free-stream velocity, i.e., V_C=V_∞.

Under the above definitions, it is clear that the time step used by standard UVLM codes based on time-stepping schemes directly depends on the spatial discretization of the blade, i.e., how fine or coarse the aerodynamic grid is. This fact is particularly important for generating the kinematics of heterogeneous wind farms. In a scenario like this, each wind turbine will have its own time step, which greatly complicates the aerodynamic simulation of the wind farm. A simple way to overcome this problem is to define an unique time step for the entire wind farm, Δt_wf. Clearly, there are several ways to define such Δt_wf: (i) a minimum time step, (ii) a maximum time step, or (iii) an average among the time steps associated with different wind turbines. These three options are available in UVLMeshGen by specifying 'min', 'max', or 'average' in the last cell of Kinematic_FLAG. According to the option selected, Δt_wf is computed as follows:

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& \mathit{\left\{}\mathtt{\text{'min'}}\mathit{\right\}}⟶\mathrm{\Delta }{t}_{\mathrm{wf}}=min\mathit{\{}\mathrm{\Delta }{t}_{\mathrm{1}},\mathrm{\Delta }{t}_{\mathrm{2}},\mathrm{\dots },\mathrm{\Delta }{t}_{{\mathrm{N}}_{\mathrm{WT}}}\mathit{\}}\\ & \mathit{\left\{}\mathtt{\text{'max'}}\mathit{\right\}}⟶\mathrm{\Delta }{t}_{\mathrm{wf}}=max\mathit{\{}\mathrm{\Delta }{t}_{\mathrm{1}},\mathrm{\Delta }{t}_{\mathrm{2}},\mathrm{\dots },\mathrm{\Delta }{t}_{{\mathrm{N}}_{\mathrm{WT}}}\mathit{\}}\\ & \mathit{\left\{}\mathtt{\text{'average'}}\mathit{\right\}}⟶\mathrm{\Delta }{t}_{\mathrm{wf}}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{WT}}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{WT}}}\mathrm{\Delta }{t}_i,\end{array}\end{array}$$

where Δt_i is the time step associated with each wind turbine in the farm. In Table 6, we provide a summary of further kinematic variables that must be specified in each wind turbine sheet file (variable WTNames described in Sect. 3.2).

Table 6 Kinematic variables.

User-defined functions must be placed inside the “Kinematic Files/” folder. The mesher will use such M-Files functions to generate the kinematic laws for rotor, yaw, and pitch depending on whether such options are enabled in Kinematic_FLAG or not. Each of these functions must receive three inputs: (i) number of time steps, (ii) time step, and (iii) initial angle. For functions related to the rotor and yaw kinematics, the third argument is a scalar variable that specifies their initial angles. In the case of pitching, such a variable is a vector of dimension N_B containing the initial pitch of each blade in the rotor. As output, the function must provide two arrays of dimension: 1×N_Δt for rotor and yaw kinematics and N_B×N_Δt for pitch kinematics. The first output variable contains the angle time series, and the second one contains the time derivative of the angle time series. As an example, we provide below a very simple user-defined function which allows us to impose a harmonic yaw motion on a rotor wind turbine.

function [A, DA] = MyYaw (NTS, DT, A0)
A0 = A0 * pi / 180;
YH = 60 * pi / 180;
YW = 0.35;
T = linspace(0,NTS*DT,NTS+1);

A = A0 + YH * sin (YW * T);
DA = YH * YW * cos (YW * T);

It should be noted that all input angles are in degrees.

4 Numerical results

In this section, we present a series of results to show the capabilities of the mesh generator developed here to build arbitrary wind turbines and wind farm UVLM meshes. To this end, we use the UVLMeshGen along with the VLMSim to reproduce several standard results in the field of wind energy. Furthermore, we present some numerical simulations related to well-established wind turbine concepts such as the Sandia 100 m 13.2 MW wind turbine SNL100-00 and the DTU 10 MW reference wind turbine. Finally, we show the versatility and potential of the mesher through the generation of two different entire wind farms. Additionally, without falling into quantitative aspects about the aerodynamic loads generated on the blades, we present some qualitative snapshots of the wakes emanating from the wind park by using the method described in Appendix C.

All study cases presented in this work, the VLMSim solver was run on a desktop computer with an Intel(R) Core(TM) i7-8700 CPU at 3.20 GHz with 16 GB of RAM memory.

4.1 Solidity study

The rotor solidity is a dimensionless number commonly used for designing rotors, such as rotorcraft, propellers, and wind turbines. This is function of the aspect ratio and number of blades in the rotor thus providing a measure of how close a lifting rotor system is to an ideal actuator disk in momentum theory. Rotor solidity σ is defined as the fraction of the annular area in the control volume which is covered by the blades, i.e.,

$$\begin{array}{}\text{(18)}& \mathit{\sigma }\left(r\right)={\displaystyle \frac{c\left(r\right)N_{\mathrm{B}}}{\mathrm{2}\mathit{\pi }r}},\end{array}$$

where c(r) is the local chord and r is the radial position of the annular control area.

To analyze the influence of rotor solidity on the power output of a wind turbine, we consider 10 different configurations, where the number of blades is varied from 1 to 10. We select a blade similar to the SNL100-00 wind turbine. The rest of the parameters are as follows: air density ρ=1.29 kg m⁻³, free-stream velocity V_∞=13.0 m s⁻¹, wind turbine radius R=110 m, angular velocity Ω_r=7.44 RPM, pitch angle ${\mathit{\theta }}_{\mathrm{p}}=\mathrm{0}{}^{\circ }$, and rotor swept area A_r=πR². The blade was discretized into 10 elements along the chord and 40 elements along the span. The characteristic length and velocity are L_C=2 m and V_C=V_∞ respectively, which in turn gives a time step $\mathrm{\Delta }t=T_{\mathrm{C}}=L_{\mathrm{C}}/V_{\mathrm{C}}=\mathrm{0.1538}$ s. The rest of the simulation parameters are number of time steps N_Δt=500, bound-vortex-lattice cut-off δ_c=0.01, and wake cut-off δ_c=0.01 (see Appendix C). Using such parameters, the wake length obtained for all rotors considered in this subsection is approximately 4.5 times the wind turbine rotor diameter.

Figure 14 Power coefficient as a function of the rotor solidity.

As a reference power, we select the wind power crossing the rotor area, i.e.,

$$\begin{array}{}\text{(19)}& P_{\text{ref}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{V}_{\mathrm{\infty }}^{\mathrm{3}}{A}_r,\end{array}$$

which allows us to introduce the power coefficient as $C_{\mathrm{p}}=P/P_{\text{ref}}$. Figure 14 depicts the C_p as a function of the number of blades. VLMSim predicts a solution that is in complete agreement with theoretical results. The study shows that power coefficient increases rapidly at first with the number of blades. Although the power continues to grow, the curve presents a very noticeable asymptotic behavior for 10 blades. Such a trend exhibits a difference of approximately 30 % with respect to the well-known Betz limit for 10 blades. Let us remember that it establishes a theoretical limit for the power that can be extracted from the wind ($C_{\mathrm{p}}=\mathrm{16}/\mathrm{27}\approx \mathrm{0.5926}$), so that such a difference of around 30 % can be associated with a number of factors, such as finite number of blades, blade geometry, non-optimal rotor configuration, etc. In practice real wind rotors have maximum C_p values in the range of 25 %–45 % (Bedon et al., 2012).

4.2 Yaw study

Another key factor associated with the power extraction from the wind is the yaw angle of the wind turbine with respect to the free-stream velocity. Rotor yaw reduces the effective projected area exposed to wind flow, thus reducing the energy conversion efficiency of the turbine. Yaw occurs when the wind direction is not perpendicular to the rotor plane. As a consequence, the blade will experience a varying relative velocity and angle of attack, leading to even more unsteady aerodynamics phenomena.

It is clear that the effective velocity of the wind to be considered to estimated the output power of a yawed rotor is its projection on the rotor axis (see Fig. 15) (Gebhardt, 2012), which can be expressed as

$$\begin{array}{}\text{(20)}& V_{\mathrm{\infty },\mathrm{e}}=V_{\mathrm{\infty }}\cos \mathit{\gamma }\cos {\mathit{\theta }}_y,\end{array}$$

where the tilt angle γ was previously introduced, and θ_y is the yaw angle measured in a plane normal to the vertical (see Fig. 15). In addition, the output power of a wind turbine can be expressed as a function of the aerodynamic torque and the rotor angular velocity as follows:

$$\begin{array}{}\text{(21)}& P=\mathit{M}\cdot \mathbf{\Omega }=q_E{L}_{\mathrm{C}}^{\mathrm{3}}{C}_{\mathrm{M}}\mathrm{\Omega },\end{array}$$

where $q_E=\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{V}_{\mathrm{\infty },\mathrm{e}}^{\mathrm{2}}$ is the effective dynamic pressure, L_C is some characteristic length (e.g., L_C in Table 6), C_M is aerodynamic moment coefficient, and Ω the rotor angular velocity. Introducing Eq. (20) into Eq. (21) and dividing by the power at ${\mathit{\theta }}_y=\mathrm{0}{}^{\circ }$, we get the following expression for a normalized output power:

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{P}{P_{\mathrm{0}}}}={\displaystyle \frac{C_{\mathrm{M}}\left(\mathit{\xi }\right)}{C_{\mathrm{M}\left(\mathrm{0}\right)}}}{\cos }^{\mathrm{2}}{\mathit{\theta }}_y,\end{array}$$

where the ratio $C_{\mathrm{M}}\left(\mathit{\xi }\right)/C_{\mathrm{M}}\left(\mathrm{0}\right)$ is a nonlinear function of θ_y. However, it approaches unity for small values of yaw angle. Therefore, as a first approximation, it can be considered that the ratio $P/P_{\mathrm{0}}$ behaves like the function cos ²θ_y.

Figure 15 Free-stream direction.

To analyze the influence of yawed rotors on the power output of a wind turbine, we consider different values of θ_y ranging from −60 to 60^∘. For this study, we selected the SNL100-00 wind turbine operating under the same working conditions and spatial/temporal discretizations as in Sect. 4.1. The resulting output power is then normalized with respect to the power at ${\mathit{\theta }}_y=\mathrm{0}{}^{\circ }$.

Figure 16 Normalized power as a function of the yaw angle.

In Fig. 16 we show how the yaw angle affects the normalized output power. For angles lower than 15^∘, the power practically behaves as the function cos ²θ_y; i.e., nonlinearities associated with C_M(θ_y) are small. Beyond 15^∘, simulations diverge significantly from the function cos ²θ_y, thus meaning that nonlinearities become important.

4.3 Pitch study

One of the main control parameters in wind turbines is the pitch angle θ_p, which allows us to regulate the output power according to the environmental conditions. For flow velocities less than the rated wind speed, the blade pitch is modified in order to maximize the output power. Similarly, for flow velocities beyond the rated wind speed, the blade pitch is also modified in order to avoid excessive angular speeds or runaway events.

The pitch angle is defined as an inward rotation of the blade leading edge towards the center of rotation of the turbine; in other words, θ_p is the angle between the tip chord and the rotor plane (see Fig. 17). Moreover, the local pitch is usually expressed as a combination of the pitch angle and the twist angle, i.e., $\mathit{\phi }\left(z\right)={\mathit{\theta }}_{\mathrm{p}}+{\mathit{\theta }}^s\left(z\right)$. It should be stressed that V_total is the actual wind velocity hitting a blade section located at a distance z from the blade root.

Figure 17 Definition of the pitch angle.

To analyze the influence of the pitch angle on the output power produced for a wind turbine, we carried out a series of simulations considering the same rotor configuration as before and a pitch angle ranging from −6 to 12^∘. Figure 18 shows how the output power change as a function of the pitch angle. Note that the power is normalized with respect to the maximum power, $\stackrel{\mathrm{‾}}{P}=P\left({\mathit{\phi }}_p\right)/P_{\text{ref}}$.

As it can be observed, the power curve is shifted to the right reaching its maximum at ${\mathit{\theta }}_{\mathrm{p}}=\mathrm{4}{}^{\circ }$. For this case, a positive pitch angle is needed to get the maximum possible power. As expected, the pitch angle has a significant influence on output power, giving a glimpse of the great controllability associated with this parameter.

Figure 18 Normalized power as a function of the pitch angle.

4.4 Pre-bend and cone angle study

Besides the twist and distribution of airfoils along the blade, cone angle and blade pre-bend are other two important parameters defining the geometry of a wind turbine rotor. The cone angle β is defined as the angle between the rotor plane and the blade longitudinal axis (see Fig. 19). Its main function is to tilt the rotor so that the blades are no longer at right angles to the nacelle, thus preventing them from hitting the tower as the rotor turns in strong winds (Hau and von Renouard, 2006).

Pre-bending was also primarily conceived to achieve tower clearance, i.e., to make sure that there is sufficient distance between the blade tips and the tower during the wind turbine operation (see Fig. 19). Both methods, coned rotors and pre-bending, require adjustments to the nacelle design.

Figure 19 Definition of the cone angle and blade pre-bend.

Here we study how the cone angle and blade pre-bend affect power output independently of each other. To this end, we consider the following two different scenarios: (i) no blade pre-bend and a cone angle ranging from −12 to 12^∘ and (ii) a blade pre-bend characterized by a blade tip deflection ranging from −20 % to 20 % and a cone angle $\mathit{\beta }=\mathrm{0}{}^{\circ }$. For all cases, the pitch angle is set to zero, ${\mathit{\theta }}_{\mathrm{p}}=\mathrm{0}{}^{\circ }$. Here, pre-bending is included by using Zuteck's procedure, already described in Sect. 2.2.4, and available in UVLMeshGen. The parameters of the Zuteck formula are listed in Table 7.

Table 7 Zuteck's formula parameters.

The reader can find an explanation of such parameters in Appendix A, Table A5. Both tip deflection x₁ and cone angle range were calculated with the goal of having a tip blade deflection of 20 % of the rotor radius (i.e., blade length + hub radius = 106 m). For this study we selected again the SNL100-00 wind turbine operating under the same working conditions and spatial–temporal discretizations as in the previous subsections. For this study, the power is normalized with respect to the power obtained for the reference configuration, which is characterized by no pre-bending and cone angle $\mathit{\beta }=\mathrm{0}{}^{\circ }$.

Figure 20 Normalized power as a function of pre-bending and cone angle.

Figure 20 shows that blade pre-bending or cone angles affect the power output in a similar way. As it can be observed, the power curves are both shifted to the left reaching their maximum at $\mathit{\beta }=-\mathrm{6}{}^{\circ }$ and $x_{\mathrm{1}}=-\mathrm{6.9756}$ m. These findings suggest that those rotors having downwardly inclined blades provide higher performance than those with positive cone angles or positive blade pre-bending. However, such sorts of rotors are of no practical importance since in general cone angles and pre-bending are measures in order to achieve tower clearance. Although the results obtained in this study are original, they are only valid for the wind turbine considered here since there is not enough information available to draw general conclusions.

4.5 Aerodynamics of full wind turbines

In this subsection, we present some aerodynamic simulations of two well-known wind turbine concepts. One of them is the DTU 10 MW reference wind turbine developed by DTU Vindenergi (Institut for Vindenergi) (Bak et al., 2013). The second turbine considered here is the SNL100-00 13.2 MW wind turbine concept developed by Sandia National Laboratories (Griffith and Ashwill, 2011).

According to the technical report Bak et al. (2013), the DTU 10 MW wind turbine was designed for offshore siting for an IEC class 1A wind climate, which is characterized by a rated wind speed of 11.4 m s⁻¹, minimum rotor speed of Ω_min=6.0 RPM, and maximum rotor speed of Ω_max=9.6 RPM. In addition, the DTU 10 MW blades have pre-bends in order to guarantee tower clearance. In Table 8 we list the main geometric and kinematic parameters used to compute the power and wake evolution of the DTU 10 MW wind turbine. A complete geometric description of the wind turbine is available in Bak et al. (2013). Moreover, the reader can obtain the configuration files for generating the aerodynamic grid and kinematics for this study case in our UVLMeshGen GitHub repository (Roccia, 2023).

Table 8 DTU 10 MW wind turbine – geometric and kinematic parameters.

In Fig. 21 we present an aerodynamic simulation of the DTU 10 MW wind turbine obtained by using the VLMSim flow solver. The working parameters correspond to those specified in Table 8. Peaks appearing in the power curve arise as consequence of the blades passing through the tower influence region. According to the simulation parameters, the period of each revolution is 6.316 s, thus giving three peaks per revolution. Such a result is in total agreement with the fact of having the passage of three blades through the tower influence region per revolution.

Figure 21 Aerodynamic simulation of the DTU 10 MW wind turbine.

Such a phenomenon is recognized as an aerodynamically unsteady region, where the flow angle and velocities are significantly affected. As the upwind turbine blades pass through this region of velocity perturbations, the flow seen by the blade is directly modified, thus resulting in periodical drops of the lift forces and, therefore, a power curve with peaks. This finding is justified from a physical point of view, and its magnitude could be influenced by effects of numerical origin, for which more studies are required to understand this phenomenon.

Another aspect worth mentioning is the interaction of the wakes with the tower. As can be seen in Fig. 21, it may seem that the flow solver considers some sort of wake rupture as used by Gebhardt et al. (2010). However, VLMSim does not have such a capacity, but it can handle wake–body interference to some extent by using a vortex core-growth strategy based on the Lamb–Oseen model (Bhagwat and Leishman, 2002; Roccia et al., 2018). This add-on allows for diffusing vortex segment intensities to handle situations in which two or more vortex segments are getting extremely close. Despite the fact that both methods are radically different, they predict very similar behaviors in terms of the peaks observed in the aerodynamic loads (consequence of the tower shadow). Furthermore, it should be emphasized that the in-house solver VLMSim has been extensively verified and validated in Verstraete et al. (2023).

The second wind turbine adopted here is the Sandia 13.2 MW SNL100-00 with a 100 m blade length, which is based on the NREL 5 MW model (Jonkman et al., 2009). The cut-in, cut-out, and rated wind speeds are 3.0, 25, and 11.3 m s⁻¹, respectively. The maximum rotation rate of this variable speed machine is Ω_max=7.44 RPM. In Table 9 we list the main geometric and kinematic parameters used to compute the power and wake evolution of this wind turbine. Moreover, the reader can obtain the configuration files for generating the aerodynamic grid and kinematics for this study case in our UVLMeshGen GitHub repository (Roccia, 2023).

Table 9 Sandia 13.2 MW SNL100-00 – geometric and kinematic parameters.

In Fig. 22 we present the time series of the output power for the SNL100-00 wind turbine together with the spatial evolution of the wake. The power after 40 s of simulation is around 12.4 MW (steady state), this value being 6 % less than the theoretical power predicted in the technical report. This difference may be attributed to different sources. First of all, we are considering a constant pitch angle throughout the simulation, when in fact this machine operates with a variable pitch angle. Second, the report utilizes a version of the low-fidelity and well-known blade element momentum theory to predict the aerodynamic forces, which can lead to some differences when compared to more accurate methods such as mid-fidelity vortex-based approaches. However, the results included here are not aimed at an exhaustive aerodynamic study of the aerodynamic performance of turbines but to show the capacity of the UVLMeshGen to generate suitable UVLM meshes in a fast and versatile way.

Figure 22 Aerodynamic simulation of the Sandia 13.2 MW SNL100-00 MW wind turbine.

Finally, the peaks observed in the power curve have the same origin as those discussed above for the DTU 10 MW wind turbine; however their magnitudes are slightly larger. As previously stated, although the appearance of these peaks has a well-founded physical explanation, their magnitude may not be entirely correct and may be affected by numerical issues.

4.6 Wind energy farms

This subsection has a main goal to show the versatility and capacity of the meshing tool developed through the generation of two hypothetical wind farms. These examples give a glimpse the scalability power of the UVLMeshGen in generating aerodynamic grids of heterogeneous parks consisting of an arbitrary number of wind turbines, including the terrain modeling (for both onshore and offshore farms). The meshes obtained from UVLMeshGen are used as input for the VLMSim solver to carry out qualitative aerodynamic simulations of such wind farms.

4.6.1 Onshore wind farm

Here we focus on modeling an onshore wind farm consisting of four wind turbines, of which two turbines are of the DTU 10 MW type and the other two are of the SNL100-00 type. Figure 23 shows the layout of the wind park, and Table 10 lists the main geometric and kinematic parameters used to carry out an aerodynamic simulation of the entire wind farm.

Figure 23 Onshore wind farm layout.

Table 10 Onshore wind farm – geometric and kinematic parameters

As can be seen from the schematic, the grid used is relatively coarse in order to reduce the computational cost associated with the aerodynamic simulation. In UVLM-based solvers there are mainly two time-consuming processes: (i) the computation of the circulations, which requires solving a linear system of dimension N_pb×N_pb, and (ii) the convection of the wakes. Although for systems discretized into a large number of panels the solution of the linear system may take longer initially, as time evolves, the wake convection undoubtedly becomes the bottleneck of the aerodynamic simulation.

Figure 24 shows the aerodynamic simulation for the entire onshore wind farm after 100 time steps. Clearly, the wakes emanating from the DTU wind turbine develop more than those shed from the Sandia machine because of the different angular velocities. However, the time step used for the whole farm is the same and its value is determined according to what is explained in Sect. 3.5. As an example and to bring up the computational cost of simulating wind turbine farms, the simulation time for 100 time steps of the farm described above took 28 h on the desktop computer described at the beginning of this section.

Figure 24 Aerodynamic simulation of an onshore wind farm composed of four wind turbines.

4.6.2 Offshore wind farm

Finally, we present the modeling of an offshore wind farm consisting of nine IEA 15 MW wind turbines. This reference wind turbine is a Class IB direct-drive machine, with a rotor of 240 m and a fixed-bottom monopile support structure which was jointly developed between the National Renewable Energy Laboratory and the Technical University of Denmark (Gaertner et al., 2020). The terrain was modeled with a sinusoidal function to simulate an ocean wave profile, similar to the example shown in Fig. 13, to incorporate sea level as a boundary for the flow solver. Although this wave does not have movement (or kinematics), it is not the objective of this work to carry out an exhaustive study of offshore farms but rather to show the versatility of the meshing tool. Figure 25 shows the layout of the offshore wind park, and Table 11 lists the main geometric and kinematic parameters used to carry out an aerodynamic simulation of the entire wind farm.

Figure 25 Aerodynamic simulation of an offshore wind farm composed of nine wind turbines.

Table 11 Offshore wind farm – geometric and kinematic parameters

Figure 25 shows the aerodynamic simulation for the entire offshore wind farm after 200 time steps. As before, the grid used for this farm is coarse in order to reduce the computational cost associated with the aerodynamic simulation. Even with a coarse mesh, the total number of panels is around 35 000, which translates into an aerodynamic matrix of the order of 1×10⁹ elements. The solution of the linear system plus wake convection has increased the computational cost for this example from 28 h (onshore farm) to 96 h.

Although the aerodynamic simulations of wind turbine farms presented above are purely qualitative in nature, they do highlight the time-consuming problem associated with these types of studies. Although these cases have been run on a relatively outdated desktop computer, it is necessary to implement methods to speed up UVLM-based solvers. On this basis, it is essential to explore the use of the fast multipole technique in order to reduce from O(n²) to O(n) and O(nlog n) the number of operations performed evaluating the Biot–Savart law (Willis et al., 2007; Bogateanu et al., 2010; Deng et al., 2021), or the implementation of parallelization and vectorization algorithms in graphics processing units (GPU) using CUDA for example (Chabalko et al., 2013; Türkal et al., 2014).

5 Conclusions

In this article, we presented a detailed description of the geometric modeling and computational implementation of an interactive and versatile UVLM-oriented mesh generator for wind turbines and onshore–offshore wind farms. The meshing tool was developed entirely in MATLAB^® and easily adaptable to GNU OCTAVE. We also provided a full explanation of the input data needed by the tool, including tables where the reader can find a description of all the variables and their names in MATLAB^®. The output data provided by UVLMeshGen, nodal coordinates and connectivity arrays, were successfully used by the UVLM-based solver VLMSim. In addition, the meshing tool was robust when generating different configurations of WTs and WFs according to airfoil data, geometric parameters, terrain topography, wind farm layouts, and meshing requirements. UVLMeshGen has been tested with a large number of examples and has proven to be efficient in building aerodynamic grids of onshore–offshore wind turbines and wind farms. Furthermore, UVLMeshGen was intensively used to generate all the meshes for the aerodynamic simulations included in Sect. 4. Although such results were mainly conceived to show the capabilities of the developed meshing tool, the versatility of the mesh generator allowed us to investigate different rotor configurations, whose aerodynamic characteristics are not commonly found in the literature. Among these findings, coning angle and pre-bending were found to affect the output power in a similar way.

Although a freely available meshing tool such as UVLMeshGen may significantly contribute to the community focused on wind farm aerodynamic simulations, it still has important limitations that should be addressed in the future as a follow-up to this contribution. Among the most important improvements that can be made we identify (i) the meshing of different types of substructures (for offshore wind energy), (ii) the meshing of the blade considering its thickness, (iii) the kinematics of the substructure, and (iv) the kinematics of the sea surface (to simulate waves). Finally, we encourage the community to actively participate in this open project related to providing and improving meshing tools intended for potential flow solvers for the wind energy sector.

Appendix A: Wind turbine components

Table A1 Tower variables.

Table A2 Monopile variables.

Table A3 Nacelle variables.

Table A4 Hub variables.

Table A5 Blade variables.

Table A6 Ground variables.

Appendix B: Main variable structures

Table B1 WIND_TURBINE – main fields and substructure variables.

Table B2 CONNECT – main fields and substructure variables.

Table B3 GROUND_FARM – main fields and substructure variables.

Table B4 KINEMATICS(i).WT – main fields and substructure variables.

Appendix C: Unsteady vortex-lattice method

In this appendix, we present a brief review of the unsteady vortex-lattice method in order to highlight the relationship between the geometric modeling introduced in previous sections and the data needed by UVLM-based solvers when it comes to aerodynamic simulations of complex aeronautical/mechanical engineering applications – here wind energy farms in particular.

According to Preidikman (1998) and Katz and Plotkin (2001), in UVLM-based computational implementations, the continuous bound-vortex sheets are discretized into a lattice of short, straight vortex segments of circulation Γ(t). Such segments divide ∂ℬ into a finite number of elements B_k (also called panels or boundary elements). The wakes shed from the separation zones (trailing edges (TEs), wing tips or blade tips, and leading edges (LEs)) are also represented by vortex lines. In Fig. C1 we present a schematic representation of the vortex lattices for the hub–nacelle assembly of a wind turbine.

Figure C1 Vortex lattice associated with the hub–nacelle assembly.

Following Verstraete et al. (2023), the complete boundary of an aeronautical–mechanical system is geometrically decomposed into a finite set of boundary elements ${\mathcal{A}}_i=\mathit{\left\{}{B}_k^i\mathit{\right\}}$, such that

$$\begin{array}{}\text{(C1)}& \mathcal{A}=\bigcup _{i\in S_B}{\mathcal{A}}_i\text{and}{\mathcal{A}}_i=\bigcup _{k\in S_i}{B}_k^i,\end{array}$$

where $S_B=\mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{B}}\mathit{\}}$, N_B is the number of bodies, $S_i=\mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },N_{{\mathrm{pb}}_i}\mathit{\}}$, and $N_{{\mathrm{pb}}_i}$ is the number of panels associated with each aerodynamic subgrid 𝒜_i. Then, the total number of panels used to discretize the whole surface 𝒜 is calculated as $N_{\mathrm{pb}}={\sum }_{i=\mathrm{1}}^{N_{\mathrm{B}}}\text{card}\left(S_i\right)$.

It is well known that UVLM solvers strongly depend on the quality with which lifting and non-lifting surfaces are represented. Wind turbines, and even more so wind farms, are characterized by very complex geometries in general (rotor, blades, terrain topography, etc.), and therefore robust and precise meshing processes are needed to ensure a correct estimation of aerodynamic loads. In this sense, it has been found that the geometric entities GO₁ and GO₂, described in Sect. 2, allow us to generate all the grids and subgrids associated with sets 𝒜_i.

As mentioned above, the edges of each $B_k^i$ are represented by straight, finite vortex segments of circulation Γ(t), whose contribution to the velocity field is computed through a discrete version of the Biot–Savart law:

$$\begin{array}{}\text{(C2)}& \mathit{V}(\mathit{r},t)={\displaystyle \frac{\mathrm{\Gamma }\left(t\right)}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{\left({\mathit{r}}_{\mathrm{1}}\times {\mathit{r}}_{\mathrm{2}}\right)\left(\left||{\mathit{r}}_{\mathrm{1}}\right||+\left||{\mathit{r}}_{\mathrm{2}}\right||\right)}{\left||{\mathit{r}}_{\mathrm{1}}\right|\left|\left||{\mathit{r}}_{\mathrm{2}}\right|\right|\left(\left||{\mathit{r}}_{\mathrm{1}}\right|\left|\left||{\mathit{r}}_{\mathrm{2}}\right|\right|+{\mathit{r}}_{\mathrm{1}}\cdot {\mathit{r}}_{\mathrm{2}}\right)+{\left({\mathit{\delta }}_{\mathrm{c}}\left||\mathit{u}\right||\right)}^{\mathrm{2}}}},\end{array}$$

where r₁ and r₂ are the position vectors of the point where the velocity is being evaluated relative to the ends of the straight vortex segment, $\mathit{u}={\mathit{r}}_{\mathrm{1}}-{\mathit{r}}_{\mathrm{2}}$, and δ_c is a cut-off parameter, which is introduced to remove its singular kernel. Although introducing the term $\left({\mathit{\delta }}_{\mathrm{c}}\left||\cdot \right||\right)$ into Eq. (C2) is interpreted as essentially an ad hoc technique (Chorin, 1994; van Garrel, 2003), it has been proven to work satisfactorily well in practice.

C1 Aerodynamic influence coefficients

In UVLM-based codes, the non-penetration condition is generally imposed either at the geometric centers of each $B_k^i$ (the so-called control or collocation points, CPs) or at $\mathrm{3}/\mathrm{4}$ of the panel chord (Katz and Plotkin, 2001), resulting in a linear system of algebraic equations (usually with time-varying coefficients). The unknowns are the circulations around the individual bound vortex segments; however, the linear system can be rewritten in terms of vortex ring circulations G_j(t), thus reducing the size of the problem (Verstraete et al., 2023). Under these assumptions, the linear system takes the following form:

$$\begin{array}{}\text{(C3)}& \begin{array}{rl}\sum _{j=\mathrm{1}}^{N_{\mathrm{pb}}}& a_{ij}\left(t\right)G_j\left(t\right)+\left[{\mathit{V}}_{\mathrm{\infty }}+{\mathit{V}}_W({\mathit{r}}_i,t)\right.\\ & \left.-{\mathit{V}}_{\mathrm{S}}({\mathit{r}}_i,t)\right]\cdot {\widehat{\mathit{n}}}_i\left(t\right)=\mathrm{0},i=\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{pb}},\\ & \mathbf{A}\left(t\right)\mathit{G}\left(t\right)=\mathit{R}\mathit{H}\mathit{S}\left(t\right),\end{array}\end{array}$$

where a_ij(t) are the aerodynamic influence coefficients, ${\widehat{\mathit{n}}}_i$ is the unit vector normal at the ith control point, V_W is the velocity induced by the free-vortex lattice, V_S is the velocity of the solid, V_∞ is the free-stream velocity, A(t) is the aerodynamic influence matrix, G(t), and RHS(t) is the right-hand side, which collects the contributions of the wake, free-stream, and body velocities along the normal direction at each CP.

As it can be observed in Eq. (C3), solving for the unknown circulations requires knowing the body velocities V_S at each $B_k^i$ control point. Such velocities depend on the kinematics imposed on the rotor (including yaw and pitch motions, if any), substructure kinematics, and sea level surface motions for offshore wind turbines. Consequently, the position and velocity data provided by the kinematic module of the UVLMeshGen (described in Sect. 3.5) play a fundamental role in developing high-quality aerodynamic simulations of arbitrary onshore–offshore wind farms.

C2 Wake convection

In order to formalize the convection process, let us consider 𝒱_i for $i=\mathrm{1},\mathrm{\dots },N_W$ be a set of panels ${\mathcal{V}}_i=\left\{L_k^i\right\}$ representing the wake shed from shedding zones (SZs) of 𝒜_j∈𝒜, such that

$$\begin{array}{}\text{(C4)}& \mathcal{V}=\bigcup _{i\in S_W}{\mathcal{V}}_i\text{and}{\mathcal{V}}_i=\bigcup _{k\in W_i\left(t\right)}{L}_k^i,\end{array}$$

where $S_W=\mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },N_W\mathit{\}}$, N_W⩽N_B is the number of lifting surfaces, $W_i\left(t\right)=\mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\dots },N_{{\mathrm{pw}}_{\mathrm{i}}}(t\left)\mathit{\right\}}$, and $N_{{\mathrm{pw}}_{\mathrm{i}}}\left(t\right)$ is the number of vortex rings in 𝒱_i. Then, the total number of free-vortex rings at time t is determined as $N_{\mathrm{pw}}\left(t\right)={\sum }_{i=\mathrm{1}}^{N_W}\text{card}\left(W_i\right(t\left)\right)$.

Once the circulations G_j(t) are calculated, the wakes are convected to their new positions and new vortex segments are propagated into the free-vortex lattices. Because all the quantities involved in the convection are functions of time, the question of which instantaneous quantities to use in the approximation is raised. There are several options; for example, one can use the quantities that were calculated at the previous time step, the present time step, or their averaged values for the two time steps. In all cases except the first, iterations are needed, which increase the computational cost. Kandil et al. (1976) showed that explicit one-step methods are stable, and there are little differences in the computed results when compared with higher-order procedures. In this respect, here we use an explicit first-order method to propagate the wake:

$$\begin{array}{}\text{(C5)}& \begin{array}{rl}{\mathit{r}}_{\text{node}}& (t+\mathrm{\Delta }t)\approx {\mathit{r}}_{\text{node}}\left(t\right)+{\mathit{V}}_{\text{node}}\left(t\right)\mathrm{\Delta }t,\\ & \text{node}=\mathrm{1},\mathrm{\dots },N_{\mathrm{nw}}\left(t\right),\end{array}\end{array}$$

where the subscript “node” is introduced to refer to the corners of a vortex segment, N_nw(t) is the number of aerodynamic nodes in 𝒱, and Δt is the time step. The vector V_node(t) collects the contributions from all surface vortex rings $B_k^i$, all free-vortex rings $L_j^i$, and the free-stream velocity.

C3 Aerodynamic loads

Among the several procedures proposed in the literature for computing aerodynamic loads by using the UVLM, we can mention the Joukowski approach (Simpson et al., 2013; Lambert and Dimitriadis, 2017), the Katz approach (Katz and Plotkin, 2001; Lambert and Dimitriadis, 2017), and an alternative formulation based on the Katz method developed at Virginia Tech (VT) (Preidikman, 1998). Here we present a quick review of the VT approach for predicting aerodynamic forces.

The VT approach is similar to that proposed by Katz. It computes the pressure jump across the bound-vortex lattice by using the Bernoulli equation for unsteady flows:

$$\begin{array}{}\text{(C6)}& \begin{array}{rl}{\displaystyle \frac{D_{\mathrm{p}}}{\mathit{\rho }}}& =\left[{\left.({\partial }_t\mathit{\phi }+{\partial }_t\mathit{\psi })\right|}_U-{\left.({\partial }_t\mathit{\phi }+{\partial }_t\mathit{\psi })\right|}_L\right]\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{V}}_U\cdot {\mathit{V}}_U-{\mathit{V}}_L\cdot {\mathit{V}}_L\right),\end{array}\end{array}$$

where D_p is the pressure jump, ρ is the fluid density, ${\partial }_t(\cdot )$ stands for partial time derivative, ${\mathit{V}}_U={\left.(\mathrm{\nabla }\mathit{\phi }+\mathrm{\nabla }\times \mathbf{\Psi })\right|}_U$, ${\mathit{V}}_L={\left.(\mathrm{\nabla }\mathit{\phi }+\mathrm{\nabla }\times \mathbf{\Psi })\right|}_L$, ∇(⋅) is the Nabla operator in ℝ³, φ is a scalar potential, Ψ is a vector potential, and ψ is another scalar potential such as $\mathrm{\nabla }\times \mathbf{\Psi }=\mathrm{\nabla }\mathit{\psi }$. The component of the velocity field coming from the scalar potential is irrotational, while any vorticity contribution to it is captured by the vector potential component.

After some algebraic manipulations, the pressure jump for the k element in 𝒜 can be expressed as follows:

$$\begin{array}{}\text{(C7)}& D_{\mathrm{p},k}=\mathit{\rho }\left[{\mathit{V}}_k^{\mathrm{m}}-{\mathit{V}}_k\right]\cdot \mathrm{\Delta }{\mathit{V}}_k+\mathit{\rho }{\displaystyle \frac{G_k\left(t\right)-G_k(t-\mathrm{\Delta }t)}{\mathrm{\Delta }t}},\end{array}$$

where ${\mathit{V}}_k^{\mathrm{m}}={\mathit{V}}_{\mathrm{B},\mathrm{k}}+{\mathit{V}}_{W,k}+{\mathit{V}}_{\mathrm{\infty },k}$ is the “mean” velocity, which does not recognize the presence of the local vorticity, and ΔV_k represents the jump in the tangential velocity across B_k. Finally, the vector force on B_k, F_k, is obtained as the product of Eq. (C7) times the element area times the normal unit vector at CP_k:

$$\begin{array}{}\text{(C8)}& {\mathit{F}}_k=D_{\mathrm{p},k}^d{A}_k{\widehat{\mathit{n}}}_k.\end{array}$$

For more details about the theory as well as implementation aspects related to the UVLM, the reader is referred to Preidikman (1998); Roccia et al. (2013); Verstraete et al. (2023)