**Numerical Investigation of Vortex Induced Vibration for Submerged Floating Tunnel under Di**ff**erent Reynolds Numbers**

#### **Ruijia Jin 1,2, Mingming Liu 3,\*, Baolei Geng 1, Xin Jin 3, Huaqing Zhang <sup>1</sup> and Yong Liu <sup>2</sup>**


Received: 9 October 2019; Accepted: 3 January 2020; Published: 7 January 2020

**Abstract:** A 2D numerical model was established to investigate vortex induced vibration (VIV) for submerged floating tunnel (SFT) by solving incompressible viscous Reynolds average Navier-Stokes equations in the frame of Abitrary Lagrangian Eulerian (ALE). The numerical model was closed by solving SST *k*-ω turbulence model. The present numerical model was firstly validated by comparing with published experimental data, and the comparison shows that good achievement is obtained. Then, the numerical model is used to investigate VIV for SFT under current. In the simulation, the SFT was allowed to oscillate in cross flow direction only under the constraint of spring and damping. The force coefficients and motion of SFT were obtained under different reduced velocity. Further research showed that Reynolds number has not only a great influence on the vibration amplitude and 'lock-in' region, but also on the force coefficients on of the SFT. A large Reynolds number results in a relatively small 'lock-in' region and force coefficient.

**Keywords:** Navier-Stokes equation; SST *k*-ω turbulence model; vortex-induced vibration (VIV); Arbitrary Lagrangian Eulerian (ALE) method; finite element method (FEM)

#### **1. Introduction**

Submerged floating tunnel (SFT) is a new type of traffic structure that crosses the strait, bay and lake. It is usually suspended more than 30 m below the water. The SFT has a large internal space, which is sufficient to meet the requirements of roads and even railways. For some fjords with harsh natural conditions, due to environmental conditions and technical constraints, traditional spanning methods (such as: cross-sea bridges, immersed tunnels) are not feasible, and SFT offers the possibility of crossing. Since the SFT is always in a deep-water depth, whose location is more the half wave length of normal wave, the normal wave has less influence on it. In addition, due to severe natural environment in the fjord, the flow velocity is usually fast, which has a greater influence on the SFT.

In view of the coupled analysis of flow and SFT, many scholars have done the relevant researches. Mai [1] considered the fluid-structure interaction effect, and studied the effects of surface velocity, tunnel section form and support form on the dynamic response of the SFT. It is found that the surface velocity would significantly affect the response displacement of the SFT, but it did not affect the stress distribution along the axial direction. Wang [2] analysed the variation law of the load on the SFT structure under the lateral lift force of the flow with the submerged depth, water depth, flow velocity and section size. Long [3] studied the dynamic response of SFT in different Buoyancy Weight Ratios (BWRs), and proposed the optimal range of BWRs for SFT under flow loading. The empirical lift formula based on Morison's formula was used in the above studies when discussing the flow load,

but the effect of SFT on the flow field was not considered in detail. In order to find the optimal section of SFT, Luo et al. [4,5] compared the flow field distribution and force acting on the fixed SFT with different cross-section forms by large eddy numerical simulation. It was found that the ear-shaped (Figure 1) SFT structure had smaller lift coefficient and drag coefficient, which was the most reasonable cross-section shape, followed by circle, ellipse, hexagon and rectangular.

**Figure 1.** A sketch of ear-shape SFT structure.

Since the SFT is suspended in the sea via the mooring system, it will move under the flow action. The periodically varying lift makes the SFT with elastically support vibrate perpendicular to the inlet flow direction, that is, 'Vortex-Induced Vibration' (VIV). When the vortex shedding frequency is close to the natural frequency of the structure, the phenomenon of resonance or lock-in occurs, which reflects the complex interaction between the fluid and the structure. Therefore, when frequency lock-in occurs in the SFT under the flow action, the fatigue damage of the structure will be significantly increased, which will have a negative impact on the safety of the project. Many experimental researches and numerical simulations on the VIV problem were carried out, such as Morse and Williamson [6], Govardhan and Williamson [7], Yan et al. [8], Luo et al. [9], Zheng et al. [10]. VIV experiment of rigid cylinder with elastically support under wind load was successfully studied by Feng [11]. Williamson and Khalak [12,13] and Govardhan et al. [14] performed VIV experiments on rigid columns with low mass ratio and elastic support in the wave tank, which became the verification test for many subsequent numerical simulations. Lu and Dalton [15] numerically studied the cylindrical VIV problem with cross-flow motion in the case of Reynolds number Re = 13,000, and the model used large eddy simulation to close the turbulence equation. Dong and Karniadakis [16] used direct numerical simulation (DNS) to study the force vibration of cylinders with a Reynolds number of 10,000. For the vibration analysis of SFT, Ge et al. [17] and Kang et al. [18] applied the von der Pol equation to simulate the VIV of the SFT in flow, and studied the influence of the spacing of the anchor chain on the vibration amplitude of the tunnel. Su and Sun [19] utilised the wake oscillation model to simulate the VIV of the SFT. It was found that the vortex-induced vibration resonance occurred and the axial stress of the structure increased significantly. Chen et al. [20] proposed a simplified theoretical model for vibration analysis of the coupled SFT tube-cable system under wave and current.

In general, the empirical formulas are employed in most research about the VIV of SFT, and the study based on computational fluid dynamic (CFD) are limited on a specific Reynolds number. For this reason, the investigation of VIV based on CFD in a serious of Reynolds number under a specific engineering background should be proceeded here. The vibration phenomenon of SFT under different Reynolds number is different because the size of the structure and flow velocity vary greatly in different situations. Based on the above background, this study aims to analyse the coupled motion of flow and SFT in the range of Reynolds number from 1000 to 100,000, and study the influence of different Reynolds numbers on the vibration of SFT as well as force coefficient under the flow action, so as to provide reference value for practical engineering.

#### **2. Numerical Model**

For underwater SFT, the flow around the structure is usually turbulent due to the large scale of the structure itself and the relatively fast flow velocity. At present, the simulation of turbulence can be approximated by direct numerical simulation (DNS) or by using a suitable turbulence model, and the turbulence model is used to simulate the problem in this paper. At the same time, since the length of the SFT structure is much larger than the section size, we can simulate this problem as a two-dimensional (2D) flow-structure interaction problem approximately. Although it is well known that the flow is indeed three-dimensional (3D) effect for large Reynolds number, however, 3D simulation cost a lot of computational resources. Therefore, a 2D numerical model was adopted in this study. Although 2D numerical model will overpredict the numerical results, it still can reveal the relationship between reduced velocity, vibration amplitude and force coefficient. In addition, many other researchers have also adopted 2D numerical model to solve similar problems (Lu and Dalton [15], Dong and Karniadakis [16]).

#### *2.1. Governing Equation and Turbulence Model*

The two-dimensional incompressible Reynolds-Averaged Navier-Stokes equations are adopted to describe the turbulence flow of incompressible viscous fluid. The governing equations in the Arbitrary Lagrangian-Eulerian (ALE) frame can be written as (Liu et al. [21])

$$\frac{\partial \overline{u\_i}}{\partial t} + (u\_j - u\_j^{\text{uv}}) \frac{\partial \overline{u\_i}}{\partial \mathbf{x}\_j} = -\frac{1}{\rho} \frac{\partial p}{\partial \mathbf{x}\_i} + \frac{\partial}{\partial \mathbf{x}\_j} [2\nu S\_{ij} - \overline{u\_i' u\_j'}] \tag{1}$$

$$\frac{\partial u\_i}{\partial \mathbf{x}\_i} = 0 \tag{2}$$

where *x*<sup>1</sup> = *x*, *x*<sup>2</sup> = *y* are the horizontal and vertical coordinates, respectively, *ui* is the fluid velocity in the *xi*-direction, *t* is the time, *u<sup>m</sup> <sup>j</sup>* is the velocity of moving grid in the *xj*-direction, *p* is the pressure, ρ is the fluid density, υ is the kinematic viscosity of the fluid, *Sij* is the mean strain rate tensor with *Sij* <sup>=</sup> ∂*ui*/∂*xj* + ∂*uj*/∂*xi* /2. The Reynolds stress term in Equation (1) reads can be expressed as

$$
\overline{u\_i' u\_j'} = \nu\_l \left(\partial u\_i / \partial x\_j + \partial u\_j / \partial x\_i\right) + \frac{2}{3} k \delta\_{ij} \tag{3}
$$

where υ*<sup>t</sup>* is the turbulent eddy viscosity, *k* is the turbulence kinetic energy and δ*ij* is the Kronecker operator.

In order to close the governing equations, the Shear Stress Transport (SST) *k*-ω turbulence model (Menter [22]; Menter et al. [23]) is adopted. The parameters in the equation have been widely accepted and successfully applied, which has shown good performance in simulating the boundary layer flows with significant adverse pressure gradient. The governing equation of the SST *k*-ω turbulence model can be written as follows:

$$\frac{\partial \mathbf{k}}{\partial t} + (\boldsymbol{\mu}\_{\dot{\boldsymbol{\beta}}} - \boldsymbol{\mu}\_{\dot{\boldsymbol{\beta}}}^{\text{m}}) \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{\dot{\boldsymbol{\beta}}}} = \frac{\partial}{\partial \mathbf{x}\_{\dot{\boldsymbol{\beta}}}} \Big[ (\boldsymbol{\nu} + \sigma\_{k} \boldsymbol{\nu}\_{t}) \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{\dot{\boldsymbol{\beta}}}} \Big] + P\_{k} - \boldsymbol{\beta}^{\star} \boldsymbol{\alpha} \mathbf{k} \tag{4}$$

$$\frac{\partial \boldsymbol{\omega}}{\partial t} + (\boldsymbol{u}\_{\boldsymbol{j}} - \boldsymbol{u}\_{\boldsymbol{j}}^{\rm m}) \frac{\partial \boldsymbol{\omega}}{\partial \mathbf{x}\_{\boldsymbol{j}}} = \frac{\partial}{\partial \mathbf{x}\_{\boldsymbol{j}}} \Big[ (\boldsymbol{\nu} + \boldsymbol{\sigma}\_{\boldsymbol{\omega}} \boldsymbol{\nu}\_{\boldsymbol{l}}) \frac{\partial \boldsymbol{\omega}}{\partial \mathbf{x}\_{\boldsymbol{j}}} \Big] + aS^{2} - \beta \boldsymbol{\omega}^{2} + 2(1 - F\_{1}) \boldsymbol{\sigma}\_{\boldsymbol{\omega}2} \frac{1}{\boldsymbol{\omega}} \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{\boldsymbol{j}}} \frac{\partial \boldsymbol{\omega}}{\partial \mathbf{x}\_{\boldsymbol{j}}} \tag{5}$$

where *Pk* is the production of turbulent kinetic energy and the related parameters in Equations (4) and (5) are calculated as follows

$$\upsilon\_{t} = \frac{\alpha\_{1}k}{\max(\alpha\_{1}\omega, \Omega F\_{2})}$$

$$P\_{k} = \min\left[\upsilon\_{t}\frac{\partial u\_{j}}{\partial x\_{i}}\left(\frac{\partial u\_{j}}{\partial x\_{i}} + \frac{\partial u\_{i}}{\partial x\_{j}}\right) \mathbf{1} \mathbf{0} \boldsymbol{\beta}^{\*} k \boldsymbol{\omega}\right]$$

$$F\_{1} = \tanh\left\{ \left(\min\left[\max\left(\frac{\sqrt{k}}{\boldsymbol{\beta}\,^{\*}\boldsymbol{\omega}\boldsymbol{\beta}^{\*}}, \frac{500\boldsymbol{\nu}}{\boldsymbol{y}\,^{\*2}\boldsymbol{\omega}}\right) \frac{4\rho\boldsymbol{\sigma}\_{\boldsymbol{\alpha}\boldsymbol{\beta}}\boldsymbol{k}}{D\_{k\boldsymbol{\omega}}\mathbf{y}^{\*2}}\right]\right)^{4} \right\}$$

*Water* **2020**, *12*, 171

where Ω is the absolute value of vorticity, *y*\* is the distance to the nearest solid wall, and the parameters *F*<sup>2</sup> and *Dk*<sup>ω</sup> are

$$D\_{\mathbf{k}\boldsymbol{\omega}} = \max\Big(2\sigma\_{\boldsymbol{\omega}}\varrho \frac{1}{\omega} \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{\boldsymbol{j}}} \frac{\partial \boldsymbol{\omega}}{\partial \mathbf{x}\_{\boldsymbol{j}}}, 10^{-10}\Big) \\ F\_{2} = \tanh\Big\{ \left[ \max\Big(\frac{2\sqrt{k}}{\boldsymbol{\beta}\*\boldsymbol{\omega}\boldsymbol{y}\*}, \frac{500\boldsymbol{\nu}}{\boldsymbol{y}\*^{2}\boldsymbol{\omega}} \Big) \right]^{2} \Big\}$$

By using the blending function *F*1, the following parameters can be calculated, i.e.,

$$
\sigma\_k = F\_1 \sigma\_{k1} + (1 - F\_1) \sigma\_{k2}; \\
\sigma\_{a\flat} = F\_1 \sigma\_{a\flat 1} + (1 - F\_1) \sigma\_{a\flat 2}.
$$

$$
\alpha = F\_1 \alpha\_1 + (1 - F\_1)\alpha\_2; \\
\beta = F\_1 \beta\_1 + (1 - F\_1)\beta\_2;
$$

The model constant in the SST *k*-ω model are listed in Table 1.

**Table 1.** Parameters of SST *k*-ω turbulent model.


When the flow field and the pressure field are obtained, the fluid force acting on the structure can be obtained by integrating the surface pressure and the viscous shear force over the body surface. The dimensionless drag coefficient *CD* and the lift coefficient *CL* are respectively

$$\mathcal{C}\_{D} = -\int\_{0}^{2\pi} p \cos\theta d\theta - \frac{1}{\text{Re}} \int\_{0}^{2\pi} \left(\frac{\partial v}{\partial \mathbf{x}} - \frac{\partial u}{\partial y}\right) \sin\theta d\theta \tag{6}$$

$$\mathcal{C}\_{L} = -\int\_{0}^{2\pi} p \sin \theta d\theta + \frac{1}{\text{Re}} \int\_{0}^{2\pi} \left(\frac{\partial v}{\partial \mathbf{x}} - \frac{\partial u}{\partial y}\right) \sin \theta d\theta \tag{7}$$

#### *2.2. Motion Response of SFT*

For the vibration problem of the SFT under the flow action, since it is necessary to ensure the anchor cable is always in the elastic range, the whole system can be simplified as a mass-damping-spring system. In this paper, only the vibration response of the SFT in cross flow direction is considered, and its motion equation can be expressed as follows.

$$
\dot{x}m\ddot{y} + c\dot{y} + ky = F\_y \tag{8}
$$

where *m*, *c*, and *k* are the mass, damping and stiffness of SFT, respectively. *Fy* is the fluid force of SFT in cross flow direction which is determined by the flow equation.

Utilising the relationship of structural dynamics, *c*/*m* = 4πξ *fn*, *k*/*m* = (2π*fn*) 2 , the ratio of the mass of the SFT to the mass of water discharged from SFT is defined *<sup>m</sup>*<sup>∗</sup> = 4*m*/πρ*D*2. Therefore, the expression can be obtained as follow

$$
\ddot{y} + 4\pi\xi f\_n \dot{y} + \left(2\pi f\_n\right)^2 y = \frac{4F}{\pi \rho D^2 m^\*} \tag{9}
$$

where ξ is the damping ratio of structure, *fn* is the natural frequency and *m*\* is mass ratio.

The following dimensionless relationship can be defined further

$$\ddot{Y} = \frac{\ddot{y}D}{\mathcal{U}^2}, \dot{Y} = \frac{\dot{y}}{\mathcal{U}}, Y = \frac{y}{\mathcal{D}}, F\_{\text{fl}} = \frac{f\_{\text{fl}}D}{\mathcal{U}} \tag{10}$$

At the same time according to the definition of lift coefficient, *Fy* = 0.5ρ*U*<sup>2</sup>*DCL*, Equation (9) can be transformed in dimensionless form

$$
\ddot{Y} + 4\pi\xi F\_n \dot{Y} + (2\pi F\_n)^2 Y = \frac{2C\_L}{\pi m^\*} \tag{11}
$$

By introducing the definition of the reduced velocity *Ur* = *U*/*fnD*, the motion equation of the SFT can also be expressed as a dimensionless equation in the form of reduced velocity.

$$
\ddot{Y} + \frac{4\pi\xi}{lI\_r}\dot{Y} + \left(\frac{2\pi}{lI\_r}\right)^2 Y = \frac{2C\_L}{\pi m^\circ} \tag{12}
$$

The lift coefficient at the right end of the above formula has been given before, and the dynamic response of the SFT can be calculated by the above formula.

#### *2.3. Calculation Model and Boundary Conditions*

The calculation model and boundary conditions are shown in Figure 2. Let the origin of the coordinate be at the initial centre of the cylinder, and the dimensionless cylinder diameter *D* = 1. The dimensionless velocity *u* = 1, *v* = 0 are set up at inlet. Symmetrical boundary conditions are applied at side wall ∂*u*/∂*y* = 0, *v* = 0. The outlet velocity boundary condition is ∂*ui*/∂*t* + *c*∂*ui*/∂*xi* = 0, where *c* is local average flow velocity. non-slip boundary conditions applied to the cylindrical surface *u* = *dx*/*dt*, *v* = *dy*/*dt*. In the calculation, in order to ensure the first layer of the grids is in the viscous boundary layer, the distance between the surface of the circular cylinder and the first layer of the grids is less than 0.05% *D*. In addition, the choice of boundary layer can be employed by the method of Palm et al. [24]. In the calculation, the pressure at outlet *p* = 0 and the pressure boundary condition ∂*p*/∂**n** = 0 is applied at other boundary conditions, **n** is unit normal vector pointing out the fluid domain. At initial time, the velocity and pressure in the fluid domain are set up to zero, i.e., the initial velocity field satisfies the continuous equation.

**Figure 2.** Sketch definition of computational domain and boundary conditions.

#### **3. Numerical Dispersion and Grid Update**

The convection-diffusion equation is solved by streamline upwind/petrov-galerkin finite element method in this paper (Brooks and Hughes [25]), and this method has been applied in the solution of impressive flow problem successfully (Mochida and Murakami [26], Kim et al. [27], Guilmineau and Queutey [28]). The distribution method is utilised in the time integration of the momentum equation. First, the pressure term is neglected, and the intermediate velocity of convection and diffusion term is considered; then the pressure equation is solved to calculate the pressure of the next time step; finally, the pressure gradient term is considered to correct the flow field. Streamline upwind method is used to predict the velocity.

The Newmark-β method is used to solve the motion equation of the structure. Given the displacement, velocity and acceleration at an initial moment, the time step Δ*t*, the parameters β and γ are selected, and then the equivalent stiffness is formed. The effective load at time *t* + Δ*t* is obtained, and the displacement at the time *t* + Δ*t* can be solved. For time advance, according to CFL (Courant-Friedrichs-Lewy) conditions, the following dynamic time steps are adopted:

$$
\Delta t = \mathcal{C}\_s \min \left( \sqrt{S\_c} / |u\_c| \right) \tag{13}
$$

where *Sc* is the mesh area, *ue* is the flow velocity at grid centre, min indicates the minimum in the computational domain, *Cs* is the safe coefficient, *Cs* = 0.2. Due to the reciprocating motion of the SFT under the flow action, the dynamic grid method based on ALE is used to simulate the fluid-structure coupling problem. In this paper, the mesh in the computational domain is assumed to be an elastic one, as shown in Figure 3, in order to achieve the purpose of adapting to the movement of the grid boundary nodes and internal nodes. The motion and deformation of the mesh can be obtained by solving the governing equation of linear elastodynamics (Johnson and Tezduyar [29]). The mesh updating method make the displacement of the mesh nodes more uniform and improve the stability of numerical calculation. In addition, the possibility of mesh distortion can be reduced by controlling the elastic modulus of the computing element. Specifically, the balance length of the mesh is equal to the length of the mesh itself at the initial moment. When the two ends of the joint move relative to each other, the mesh will be stretched or compressed, correspondingly. The mesh still satisfies Hooke's law, so the total force vector of any node *i* is

$$\mathbf{F}\_{i} = \sum\_{j=1}^{\nu\_{i}} \alpha\_{ij} (\mathbf{s}\_{j} - \mathbf{s}\_{i}) \tag{14}$$

where **F***<sup>i</sup>* is the total force vector on node *i*, α*ij* is the mesh stiffness between the nodes, υ*<sup>i</sup>* represents the number of nodes connected to node *i*, *j* = [*i*, υ*i*], δ*<sup>i</sup>* and δ*<sup>j</sup>* are displacement vector of node *i* and *j*. In order to avoid collisions of mesh nodes, the following expression is usually used to calculate the mesh stiffness

$$\alpha\_{i\bar{j}} = \frac{1}{\sqrt{\left(\mathbf{x}\_{\bar{j}} - \mathbf{x}\_{i}\right)^{2} + \left(y\_{\bar{j}} - y\_{i}\right)^{2}}}\tag{15}$$

where **x***i*, **x***<sup>j</sup>* are the position vectors of nodes *i* and *j*, that is, the value of α*ij* is considered to be the reciprocal of the side length. This mesh updating method has also been widely used in the VIV research (Tang et al. [30], Lu et al. [31]).

**Figure 3.** Illustration of the elastic mesh.

#### **4. Model Validation**

In order to obtain reliable numerical results, this paper firstly uses the free vibration cylinder problem under Re = 30,000 and *Ur* = 6.00 as an example to verify the mesh convergence of the numerical model. Four different meshes are considered in Table 2. From Table 2, it can be seen that the numerical results under the four meshes are very close, indicating that the numerical results have converged under the current grid density. Considering the computational efficiency, the latter numerical calculation takes mesh 3 (Figure 4) as the benchmark. In the table, *Y*max denotes the maximum vibration amplitude of the cylinder, Dismin denotes the minimum distance between the circular cylinder and the first layer gird, *C<sup>M</sup> <sup>D</sup>* is the mean drag coefficient, *<sup>C</sup>RMS <sup>L</sup>* is root mean square (RMS) of lift coefficient.


**Table 2.** Comparisons of numerical results with different mesh solutions.

**Figure 4.** Sketch of mesh 3 in the model validation.

In order to further verify the reliability of the numerical model in this paper, the coupling analysis of uniform flow and fixed cylinder is verified firstly. The diameter of cylinder is *D* = 1.0, the Reynolds number Re = 10,000. The calculation results containing the mean drag coefficient *C<sup>M</sup> <sup>D</sup>* , amplitude of lift coefficient *C<sup>A</sup> <sup>L</sup>* and strouhal number *St* are compared with experiment results conducted by Gopalkrishnan [32] and the numerical results of Dong and Karmiadakis [16], Zhao et al. [33], Song [34], which are shown in Table 3.


**Table 3.** Calculation results comparisons of fixed single cylinder.

From the comparisons in Table 3, the calculation results of *CM <sup>D</sup>* , *CA <sup>L</sup>* and *St* are almost consistent with those of other scholars, which proves the accuracy of the numerical model in the case of high Reynolds number. The problem of flow around a fixed cylinder is actually only a unilateral hydrodynamic calculation problem. It does not involve the motion response calculation of the circular cylinder itself, so it is not a true fluid-structure coupling problem. The vibration of cylinder with elastic support under the flow action involves fluid-structure coupling problem. There have been many experimental results on the VIV of rigid cylinders with elastic supports under the flow action, such as Khalak and Williamson [12], which have done systematic experimental studies. In order to facilitate comparison with the experimental results, the same calculation parameters as in the tests of Khalak and Williamson [12] were used. The Reynolds number Re = 12,000, the mass ratio *m*\* = 2.4, and the mass damping ratio *m*\*ξ = 0.013, and the maximum dimensionless vibration amplitudes versus reduced velocity are calculated in the paper.

From the comparison results in Figure 5, it can be seen that the maximum dimensionless vibration amplitude is close to 1 and the 'lock-in' region is from *Ur* = 4.0~10.0. In addition, the upper branch and lower branch can also be described in the numerical model clearly as shown in the experiment. Therefore, the results in the numerical model agree well with the experimental results of Khalak and Williamson [12]. It is illustrated that the model established in this paper can be used to investigate the fluid-structure coupling problem with high Reynolds number.

**Figure 5.** Calculation result comparisons of cylinder with spring and damping.

#### **5. Example Analysis**

Based on the above numerical model, this paper calculates the motion of the SFT under different constraint stiffness and different Reynolds numbers. The Reynolds number is calculated from 1000 to 100,000. Firstly, the time history curves of the cross-flow direction under the conditions of Reynolds number 50,000, mass ratio *m*\* = 2.5, damp ξ = 0.007 and reduced velocity *Ur* = 2.0, 5.0 and 12.0, respectively, are introduced, as shown in Figure 6. Then, the Fast Fourier Transform (FFT) is utilised in the lift force coefficient time history, and the result is displayed in Figure 7.

From Figure 5, it can be found that when reduced velocity are 2.0 and 12.0, the vibration amplitude in cross flow direction of the SFT is small, and when the reduced velocity is 5.0, the vibration amplitude is large. Subsequently, the flow field is analysed in the case of larger and smaller vibration amplitudes, as shown in Figures 8 and 9.

**Figure 6.** Time histories of submerged floating tunnel (SFT) motion under different reduced velocity: (**a**) *Ur* = 2; (**b**) *Ur* = 5; (**c**) *Ur* = 12 (Re = 50,000, *m*\* = 2.5, ξ = 0.007).

**Figure 7.** Fourier transform of SFT lift force coefficient when reduced velocity *Ur* = 2.0 (Re = 50,000, *m*\* = 2.5, ξ = 0.007).

From Figure 6, it can be seen that the vortex shedding frequency is about 0.23 Hz. According to the definition of *Ur* = *U*/*fnD*, the natural frequencies (*fn*) of SFT in the case here are 0.5 Hz, 0.2 Hz and 0.083 Hz, respectively. Therefore, it can be seen that when the reduced velocity of the SFT are 2.0 and 12.0, the vortex shedding frequency is far away from the natural frequency, and the VIV of the SFT is not obvious. The wake pattern is shown in Figure 9, and the wake shape is also regular in one vibration period. However, when the reduced velocity is 5.0 by adjusting the spring stiffness, the frequency of vortex shedding is close to the natural frequency. Under the action of flow lift force, a large VIV phenomenon occurs in the SFT. The vibration amplitude is larger, even reaching 0.7 times of the outer diameter of the SFT, whose wake pattern is shown in Figure 8. It can be seen that the wake has a long 'tail' after being separated, and the wake shape is irregular. Next, we compare the coupling effect of the flow and the SFT versus different reduced velocity under this Reynolds number, and the maximum dimensionless vibration amplitude of the SFT has been statistically obtained, as shown in Figure 10.

**Figure 8.** Vibration mode and wake vortex shedding mode of SFT in one cycle when reduced velocity *Ur* = 5.0 (Re = 50,000, *m*\* = 2.5, ξ = 0.007).

**Figure 9.** Vibration mode and wake vortex shedding mode of SFT in one cycle when reduced velocity *Ur* = 2.0 (Re = 50,000, *m*\* = 2.5, ξ = 0.007).

**Figure 10.** Vibration amplitude of SFT versus reduced velocity under Re = 50,000, *m*\* = 2.5 and ξ = 0.007.

It can be seen from the calculation results that when the reduced velocity is from 4.0 to 10.0, the structure is 'locked' under the flow action, while at other reduced velocity, the vibration amplitude of the structure is small. Then the VIV of the SFT under different Reynolds numbers are compared, and the calculation results are shown in Figure 11. From the results of VIV of SFT under different Reynolds numbers, it can be seen that Reynolds number has not only a great influence on the vibration amplitude, but also on the 'lock-in' region. In general, the lower the Reynolds number is, the larger the amplitude is out of the 'lock-in' region. The minimum vibration amplitude in the 'lock-in' region is about 0.4*D*, while the maximum vibration amplitude in the 'lock-in' region can reach to 0.8*D*. It can also be seen that larger Reynolds number leads to narrow 'lock-in' region in Figure 10. Therefore, when the size of the SFT is small or the inlet flow velocity is slow, VIV 'lock-in' phenomenon is more likely to occur for SFT because of lower Reynolds number.

**Figure 11.** Vibration amplitude of SFT versus natural frequency in different Reynolds numbers.

Figure 12 is the force coefficient on SFT versus reduced velocity at different Reynolds numbers. It can be seen that reduced velocity has a greater influence on the mean drag coefficient and RMS of lift coefficient. When the Reynolds number is low, the mean drag coefficient and RMS of lift coefficient of the SFT are relatively large from 1000 to 10,000. As the Reynolds number increases, the mean drag coefficient and RMS of lift coefficient become smaller. Therefore, when the size of the SFT is small or the inlet flow velocity action on the structure is slow, the force coefficient is large, and when the size is large or the flow velocity is fast, the mean drag coefficient and lift force coefficient of the structure are small.

**Figure 12.** Non-dimensional force coefficient on SFT versus natural frequency at different Reynolds numbers: (**a**) mean drag coefficient; (**b**) RMS of lift coefficient.

#### **6. Conclusions**

Based on the FEM solution of incompressible viscous Reynolds average Navier-Stokes equations, combining the frame of Abitrary Lagrangian Eulerian, through accurate computational fluid dynamic numerical simulation, the vortex-induced vibration problems of submerged floating tunnel with different Reynolds numbers are studied. Main conclusions are as follows:

Firstly, the analysis of uniform flow and fixed single cylinder proves the accuracy of the model in the case of high Reynolds number. Then, by simulating the vortex-induced vibration of a cylinder at high Reynolds number and comparing with other scholars' experimental results, it is proved that the model established in this paper can be used to study the fluid-structure coupling problem at high Reynolds number.

Through the research of vortex-induced vibration of SFT under the flow action, the force coefficient and motion of the SFT versus different reduced velocity at different Reynolds numbers are analysed. The results show that the Reynolds number has not only a great influence on the vibration amplitude and 'lock-in' region, but also on the force coefficient on the SFT. When the Reynolds number is low, the 'lock-in' region, the mean drag coefficient and RMS of lift coefficient of the SFT are relatively large. As the Reynolds number increases, the 'lock-in' region, the mean drag coefficient and RMS of lift coefficient become smaller. Therefore, when the size of the SFT is small or the flow velocity action on the structure is slow, the force coefficient and 'lock-in' region are relatively large, while when the size is large or the flow velocity is fast, the force coefficient and 'lock-in' region are relatively small.

**Author Contributions:** Formal analysis, X.J.; data curation, B.G.; writing—original draft preparation, R.J.; writing—review and editing, M.L. and Y.L.; funding acquisition, H.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors would like to acknowledge the support from State Key Laboratory of Coastal and Offshore Engineering Projects Program (Grant No. LP1912), National Natural Science Foundation of China (Grant No. 51809133) and China Postdoctoral Science Foundation (Grant No. 2019M652479).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Aeroelastic Performance Analysis of Wind Turbine in the Wake with a New Elastic Actuator Line Model**

#### **Ziying Yu 1, Zhenhong Hu 1,\*, Xing Zheng 1, Qingwei Ma 1,2 and Hongbin Hao <sup>1</sup>**


Received: 8 March 2020; Accepted: 21 April 2020; Published: 26 April 2020

**Abstract:** The scale of a wind turbine is getting larger with the development of wind energy recently. Therefore, the effect of the wind turbine blades deformation on its performances and lifespan has become obvious. In order to solve this research rapidly, a new elastic actuator line model (EALM) is proposed in this study, which is based on turbinesFoam in OpenFOAM (Open Source Field Operation and Manipulation, a free, open source computational fluid dynamics (CFD) software package released by the OpenFOAM Foundation, which was incorporated as a company limited by guarantee in England and Wales). The model combines the actuator line model (ALM) and a beam solver, which is used in the wind turbine blade design. The aeroelastic performances of the NREL (National Renewable Energy Laboratory) 5 MW wind turbine like power, thrust, and blade tip displacement are investigated. These results are compared with some research to prove the new model. Additionally, the influence caused by blade deflections on the aerodynamic performance is discussed. It is demonstrated that the tower shadow effect becomes more obvious and causes the power and thrust to get a bit lower and unsteady. Finally, this variety is analyzed in the wake of upstream wind turbine and it is found that the influence on the performance and wake flow field of downstream wind turbine becomes more serious.

**Keywords:** elastic actuator line model; OpenFOAM; NREL 5 MW wind turbine; aeroelastic performance

#### **1. Introduction**

With the improvement of wind power technology and the demand of high-power generation, the target of wind turbine design turns to large scale and offshore [1–7]. In 2009, the National Renewable Energy Laboratory (NREL) in America defined a 5 MW reference wind turbine for offshore system, in which the rotor diameter is 126 m [8]. Technical University of Denmark described a 10 MW reference wind turbine whose rotor diameter is 178.3 m [9] in 2013. As the blade of wind turbine gets longer, it becomes less stiff, more susceptible, and easily deformable, which will lead to increased fatigue damage and reduced production. When simulating the aeroelastic performances of a floating offshore wind turbine or a wind farm, there are many challenges to solve.

There are mainly three methods to study the aerodynamic performances of wind turbines. The first one is the Blade Element Momentum (BEM) theory, which combines the blade element theory and momentum theory. It has high efficiency and is widely used in the industrial application, but the information of flow field is not considered. The second one is the Computational Fluid Dynamics (CFD) method, which calculates the velocity and pressure fields by solving the Navier-Stokes equations. It can obtain quite accurate results and is gradually used in recent years, e.g., Tran et al. [10] analyzed the loading of 5 MW offshore wind turbine, Miao et al. [11] investigated the wake characteristics of

double wind turbines, Liu et al. [12] established a fully coupled model for simulating floating offshore wind turbine. Nevertheless, this research usually takes consuming time to run on a computer, because of the refinement mesh in the wake region and motion of the platform. The last one is Vortex lattice method, which can achieve the velocity distribution behind the rotor and calculate faster than CFD method. In this method, the vortex core model, using the empirical equations, is considered to avoid the unrealistic result and the dissipation of the wake structure. Compared with BEM theory, it pays more attention to the wake region rather than performance on the blade (both of them are considered in the CFD method). However, it may gain inaccurate results in far wake area and is seldom applied [13].

The structural dynamic characteristics are studied mainly by three methods, modal approach, Finite Element Method (FEM), and multi-body dynamics method. In the modal approach, the description of deformation can be made as a linear superimposition by some physical realistic modes. This method can reduce the number of DOF (Degree Of Freedom) and compute very fast [14]. FEM divides the structure into finite elements, in which a shape function is used to approximate the deformation [15]. However, its number of DOF is much bigger than that of former method and it will cost more computational time. In the multi-body dynamic method, every rigid part is connected by springs and hinges. This method is cheaper than FEM, but more expensive than the modal approach [14].

The aeroelastic performances of wind turbine contain the aerodynamic properties and structural responses. Multi combination methods of aerodynamic performances and structural dynamic characteristics method can solve this complex problem and have their own advantages and disadvantages. Fatigue, Aerodynamics, Structures, and Turbulence (FAST) is developed by the National Renewable Energy Laboratory (NREL), which combined the BEM and modal approach to simulate the aeroelastic performances. In a recent version, the multi-body dynamics method is also added into FAST. Ferede et al. [16] gave a framework for aeroelastic wind turbine blade analysis with BEM and multi-body dynamics method. Li et al. [17] coupled CFD and multi-body dynamics method by overset technology predict the aerodynamic performances. Ageze et al. [18] and Heinz et al. [19] used CFD and FEM to solve the fluid–structure interaction (FSI) in a wind turbine simulation.

Since BEM cannot predict wake behavior and CFD is a high computational cost method, the actuator line method (ALM), which combines advantages of BEM and CFD, was introduced by Shen and Sørensen in 2002 [20]. In ALM, the blade is replaced by a series of aerodynamic loading and this loading which acts as the body forces is added into the source item of Navier-Stokes equations. Thus, ALM is very effective and can be widely used in wind farm simulation [21,22]. Moreover, Meng et al. [23] proposed the elastic actuator line (EAL) model firstly, in which ALM and a finite difference structural model are used to analyze the wake-induced fatigue. Ma et al. [13] combined ALM and FEM to discuss the wake characteristics. The aim of the current research is developing a rapid technique to predict the aeroelastic performances. Therefore, a beam solver, which is used in the wind turbine blade design [24], is combined with ALM to accomplish this task.

In this paper, the concept of ALM is reviewed in brief, and the technique about how to add the structural solver method into turbinesFoam in OpenFOAM [25] is given. The NREL 5 MW baseline wind turbine is studied under the uniform inlet wind speed of 8 m/s and 11.4 m/s. The aerodynamic properties, including power and thrust, and the tip displacement are calculated and compared with related research and NREL's reference data, which can validate the new elastic actuator line model (EALM). Then, the difference of the aerodynamic performances with and without aeroelastic is discussed. Moreover, this variety is also studied in the wake of upstream wind turbine. It is to be observed that the foundation of the wind turbine is fixed in this research, and the floating one combined with structural dynamic characteristics of wind turbine blade is too complex and requires more validation which needs further study.

#### **2. Theoretical Model**

#### *2.1. Actuator Line Model*

In the actuator line model, the blades of wind turbine are replaced by lines with the aerodynamic force distributed on it, as shown in Figure 1. The lines are discretized into actuator sections, in which there are different airfoils, chords and structure twists. The loading in each section, which is usually named as the aerodynamic force, can be calculated by Equation (1).

$$f = (L, D) = \frac{1}{2} \rho l l^2 c (\mathbb{C}\_l \overrightarrow{\text{e}\_L} + \mathbb{C}\_d \overrightarrow{\text{e}\_D}) \tag{1}$$

where *f* is the aerodynamic force; *L* and *D* are lift and drag, respectively; ρ is the air density, and ρ = 1.225 kg/m3 in this paper; *U* is the relative velocity of the blade element in actuator sections; *c* is the chord length; *Cl* and *Cd* are the coefficients of lift and drag, which can be looked up from airfoil data tables derived from physical experiment; <sup>→</sup> *eL* and <sup>→</sup> *eD* represent the unit direction vectors of *L* and *D*.

Besides, the aerodynamic loads are concentrated on the aerodynamic center of the blade element, which are located in 0.25 chord length. To describe the influence of wind turbine blade on the flow field, the aerodynamic loads should be dispersed on the grid point. There are many distributional ways, and three-dimensional Gaussian distribution ηε is used in this paper, whose function is shown in Equation (2). This equation is a smooth function and can avoid numerical singularity.

$$\eta\_{\varepsilon}(d) = \frac{1}{\varepsilon^3 \pi^{3/2}} \exp[-\left(\frac{d}{\varepsilon}\right)^2] \tag{2}$$

where ε is a constant to adjust the strength of the distribution formula, and it is two times of local grid scale in this study according to Shives [26]; and *d* is the distance between the force point and the grid point. Therefore, the forces *f*ε on the nearby mesh can be calculated by Equation (3) and added into the momentum equation to solve the flow field as the body forces.

$$f\_{\ell} = f \times \eta\_{\ell}(d) \tag{3}$$

$$\frac{\partial V}{\partial t} + V \cdot \nabla V = -\frac{1}{\rho} \nabla p + \nu \nabla^2 V + f\_\varepsilon \tag{4}$$

where *V* and *p* are the wind speed and pressure in the flow; υ is the kinematic viscosity, and it is set as 1.5 <sup>×</sup> <sup>10</sup>−<sup>5</sup> <sup>m</sup>2/s ; *<sup>f</sup>*<sup>ε</sup> is the source item calculated through Equation (3).

Moreover, the Prandtl–Glauert tip loss model is used to the consideration that the velocity is zero at the tip. The function of tip loss model is shown as Equation (5).

$$F\_{tip} = \frac{2}{\pi} \arccos[\exp(-\frac{B(R-r)}{2r\sin\phi})] \tag{5}$$

where *Ftip* is the aerodynamic corrected force at the tip; *B* is the blades number; *R* is the radius of the blades; *r* is the distance between the root of blades and the action point of aerodynamic force; and φ is the angle between *U* and the rotor plane [27]. Tower and hub effect are also considered in this research by similar technology, and their airfoils are cylinder.

**Figure 1.** Blades represented by actuator line and discretized into actuator sections.

#### *2.2. Rotating Beam Solver*

The blades of wind turbine are considered as rotating variable cross-section projecting beams in this study. In this section, an approach of solving these beams in actual design is introduced [24]. Although this method is faulty and not suitable for structural solution, it still could catch some phenomena in vibration.

Figure 2 shows the diagram of the blade deflection, in which the blades have been substituted by actuator lines. In the picture, direction 0 is out of the rotating plane and called flapwise. Similarly, direction 1 is in the rotating plane and called edgewise. In this research, the subscript 0 means the component of the physical quantity in 0 direction and the subscript 1 represents component of the corresponding parameter in 1 direction. The diagrammatic sketch of the forces on the blade, which is considered as a cantilever beam, is given in Figure 3. According to the Newton's second law of a microelement, the sheer stresses *T*0, *T*<sup>1</sup> and the bending moments *M*0, *M*<sup>1</sup> of each blade element can be calculated from Equation (6) to Equation (9).

$$\frac{dT\_0}{dx} = -p\_0(\mathbf{x}) - m(\mathbf{x})\mathbf{g}\sin\theta + m(\mathbf{x})\ddot{\mathbf{y}}\_0(\mathbf{x})\tag{6}$$

$$\frac{dT\_1}{d\mathbf{x}} = -p\_1(\mathbf{x}) - m(\mathbf{x})\mathbf{g}\cos\theta + m(\mathbf{x})\ddot{\mathbf{y}}\_1(\mathbf{x})\tag{7}$$

$$\frac{dM\_1}{d\mathbf{x}} = T\_0 + N\_0 \tag{8}$$

$$\frac{d\mathcal{M}\_0}{d\mathbf{x}} = -T\_1 + N\_1 \tag{9}$$

where *p*(*x*) is the aerodynamic force distribution; *m*(*x*) is the blade mass distribution; *g* is the gravitational acceleration, and it is set as 9.8 m/s<sup>2</sup> in this paper; <sup>θ</sup> is the shaft tilt angle, and it is 5 degree here; .. *y* is the displacement *y* to the second derivative of time, and it can be calculated by prediction and correction of modified Euler method, given from Equation (10) to Equation (13); *N*0, *N*<sup>1</sup> are the component of centrifugal force in 0 direction and 1 direction.

$$
\tilde{y}\_{n+1} = 2y\_n - y\_{n-1} + \ddot{y}\_n dt dt \tag{10}
$$

$$
\overline{\ddot{y}}\_{n+1} = h(\overline{y}\_{n+1}, t\_{n+1}) \tag{11}
$$

$$y\_{n+1} = 2y\_n - y\_{n-1} + 0.5(\bar{y}\_n + \overline{\dot{y}}\_{n+1})dtdt\tag{12}$$

$$
\ddot{y}\_{n+1} = h(y\_{n+1}, t\_{n+1}) \tag{13}
$$

where the subscript *<sup>n</sup>* <sup>+</sup> 1 represents at time *tn*+1; the subscript *<sup>n</sup>* represents at time *tn*; *<sup>y</sup>* is predicted value; function *h* can be calculated from transposition of structure equations, as shown in Equation (15).

$$\mathcal{M}(\mathbf{x})\ddot{\mathbf{y}} + \mathcal{K}(\mathbf{x})\mathbf{y} = F(\mathbf{x}, t) \tag{14}$$

$$
\bar{y} = \frac{F(x,t)}{M(x)} - \frac{K(x)y}{M(x)} = g(y,t) \tag{15}
$$

where Equation (14) is the structure equations without damping; *M*(*x*) is mass of blade; *K*(*x*) is the stiffness of blade; *F*(*x*) is the external loading.

According to Figure 4, the bending moments *M*<sup>0</sup> and *M*<sup>1</sup> can be transformed into principal axes direction by Equation (16) and Equation (17).

$$M\_{11} = M\_1 \cos \beta - M\_0 \sin \beta \tag{16}$$

$$M\_{22} = M\_1 \sin \beta - M\_0 \cos \beta \tag{17}$$

where *M*<sup>11</sup> and *M*<sup>22</sup> are the bending moments on the first principal axis and the second principal axis; β is the structure twist angel.

**Figure 2.** Coordinate system of the blade deflection.

**Figure 3.** The diagrammatic sketch of the forces on the cantilever beam.

**Figure 4.** The blade element in the local coordinate system.

Based on the beam theory, the curvature of principal axes κ<sup>11</sup> and κ<sup>22</sup> can be obtained from Equation (18) and Equation (19).

$$\kappa\_{11} = \frac{M\_{11}}{EI\_1} \tag{18}$$

$$\kappa \kappa \mathbf{2} = \frac{M\_{22}}{EI\_2} \tag{19}$$

where *EI* is the stiffness of the blade element.

Equation (18) and Equation (19) are converted back to direction 0 and direction 1 through the formulas as following:

$$\kappa \wp = -\kappa\_{11}\sin\beta + \kappa\_{22}\cos\beta \tag{20}$$

$$\kappa\_1 = \kappa\_{11}\cos\beta + \kappa\_{22}\sin\beta\tag{21}$$

Therefore, the angular deformation θ and deflection *y* can be calculated by:

$$\frac{d\theta\_1}{d\mathbf{x}} = \kappa\_1 \tag{22}$$

$$\frac{d\theta\_0}{d\mathbf{x}} = \kappa\_0 \tag{23}$$

$$\frac{dy\_0}{dx} = -\theta\_1\tag{24}$$

$$\frac{dy\_1}{d\mathbf{x}} = \theta\_0 \tag{25}$$

The root of the blade is clamped, so the boundary condition at root is:

$$
\partial\_0 = 0, \partial\_1 = 0, \underline{y}\_0 = 0, \underline{y}\_1 = 0 \tag{26}
$$

The tip is free end, so the boundary condition at tip is:

$$T\_0 = 0, T\_1 = 0, M\_0 = 0, M\_1 = 0 \tag{27}$$

#### *2.3. New Elastic Actuator Line Model*

In this section, the ALM is improved into a new elastic actuator line model (EALM) based on turbinesFoam library, which is developed by Bachant et al. [25]. This library uses ALM to simulate wind and marine hydrokinetic turbines in OpenFOAM. Its interpolation, Gaussian projection, and vector rotation functions are all adapted from NREL's Simulator for Off/Onshore Wind Farm Applications (SOWFA).

In EALM, the body forces are calculated by traditional ALM and the blade deflection is computed by rotating beam solver, which is defined in Section 2.2. The computation process of EALM is given in Figure 5, in which the difference between EALM and ALM is marked by dashed rectangle. The part of structure solver is added into the actuatorLineSource class, and the actuator point is changed in the actuatorLineElement class. From Figure 5, it can be found that this combination of ALM and structure model is one-way coupling. Compared with the research by Meng et al. [23], the part of aerodynamics solver is similar, but the way of dealing with structural solver is different (see Section 2.2). The advantages of this model are that the EALM allows large time step when the position of actuator line changes every time and computes long terms of wind turbine working. This technology will be further improved and used in the simulation of the floating offshore wind turbine in the sea, which needs more simulation time to keep the floating foundation stability under waves.

**Figure 5.** The computational flow chart of new elastic actuator line model.

#### *2.4. Computational Wind Turbine Model*

In this study, the aerodynamic performances, which are power and thrust, and the structural responses, which are represented mainly by the blade tip displacement, will be compared with different research using varieties methods to validate EALM. According to the theory above, all these three physical quantities are obtained from the aerodynamic forces along the blade. Different cases use varieties method to achieve these forces. Only if the blade properties, including the aerodynamic and structural properties, and the computational condition such as wind speed are all the same, these three quantities could make equivalent comparisons. Thus, the computational condition will be given next. The model used in this research is NREL 5 MW wind turbine, and its gross properties are listed in Table 1. The details of the blade structural and aerodynamic properties can be referenced in Jonkman's technical report [8]. To be exact, the shaft tilt angle is considered in the structural part, which will make the external force of blade increased because of the gravity component.


**Table 1.** Properties of NREL 5 MW baseline wind turbine (NREL: National Renewable Energy Laboratory).

The simulation physical domains used in this paper are the same as Yu et al. [28]. Their three-dimensional sketches are shown in Figures 6 and 7. In these pictures, A is the outer mesh region and B is the refined mesh region, where the grid is refined by 4 levels. The mesh refined ratio is 0.5 between each level. The mesh independence test has been completed in that research and 1.5 m is chosen as the minimum size of grid. Besides, it is confirmed that the wind turbine gets a better working status at the wind speed of 8 m/s and its corresponding tip speed ratio is 7.55. As a consequence of this model, the cases studied in present research are mainly in two situations, 8 m/s and 11.4 m/s. To analyze the aeroelastic performance in the wake, double NREL 5 MW wind turbines set in a line are studied, and its simulation domain is shown in Figure 7. The main contents of this research are focused on the downstream wind turbine.

Compared to the work of Yu et al. [28], the parameter settings are all the same expect that the LES (Large Eddy Simulation) model is used instead of RANS (Reynolds Average Navier-Stokes). That is because the LES model will get more accurate results about the vorticities than RANS and the influence of the wake flows to elastic blade is studied in this research. In addition, the wind turbine blade deflection solver is added in ALM, named EALM.

**Figure 6.** The simulation domain of single wind turbine.

**Figure 7.** The simulation domain of double wind turbines.

#### **3. Verification and Analysis**

In this section, the mesh independence and uncertainty analysis of actuator point number are given first. Then, the aerodynamic performance and the structural responses are compared to validate EALM. The aerodynamic performance contains the power and the thrust. The structural responses are represented mainly by the blade tip displacement.

#### *3.1. Mesh Independence and Uncertainty Analysis*

Four levels of grids are tested to prove the mesh independence under the rated wind speed, and the non-dimension coefficients of power and thrust are compared in Table 2. They are defined as power coefficient *Cp* and thrust coefficient *Ct* as below:

$$\mathcal{C}\_{\mathcal{P}} = \frac{P}{0.5 \rho V^3 \pi R^2} \tag{28}$$

$$\mathcal{C}\_{l} = \frac{T}{0.5 \rho V^{3} \pi R^{2}} \tag{29}$$

where *P* is the mechanical power; ρ is the air density; *V* is the wind speed in the flow; *T* is the thrust on the blades. It shows that the results drop slow when the grid levels up especially from level 3 to level 4. Therefore, the level 3 mesh is used in the following research.


**Table 2.** Mesh independence test of power and thrust coefficient.

According to Shives' research [26], the number of actuator point has the rule that the maximum of the distance between the adjacent point should not more than the size of blade region. Different actuator point numbers are also tested in Table 3 and Figure 8. Besides, the mesh size and computational condition are kept as constant during this test. To improve the accuracy the uncertainty coefficient *U*<sup>ξ</sup> is calculated as follow steps.

(1) Calculate the difference of aerodynamic coefficients between neighbor level, which are defined as εζ<sup>21</sup> and εζ32;

(2) The convergence ratio *R*<sup>ξ</sup> can be computed by Equation (30).

$$R\_{\xi} = \frac{\varepsilon\_{\zeta\xi2}}{\varepsilon\_{\zeta21}}\tag{30}$$

(3) The power coefficient convergence ratio is 0.125 and the thrust coefficient convergence ratio is 0.167. They are both located in the interval which greater than 0 and less than 1. So, the Richardson extrapolation method is used to get the order of accuracy *p*<sup>ξ</sup> and the estimated value of error δReξ<sup>1</sup> in Equation (31) and Equation (32).

$$p\_{\xi} = \frac{\ln(\varepsilon\_{21}/\varepsilon\_{32})}{\ln(r\_{\xi})} \tag{31}$$

$$\delta^\*\_{\text{Re}\xi1} = \frac{\varepsilon\_{\xi32}}{r\_{\xi}^{p\_{\xi}} - 1} \tag{32}$$

where *r*<sup>ξ</sup> is refinement ratio of actuator point number.

(4) The uncertainty coefficient *U*<sup>ξ</sup> can be compute by Equation (33).

$$\mathcal{LL}\_{\xi} = \left( \left| \mathbb{C}\_{\xi} \right| + \left| 1 - \mathbb{C}\_{\xi} \right| \right) \left| \delta^\*\_{\text{Re}\xi 1} \right| \tag{33}$$

where *C*<sup>ξ</sup> is correction factor and defined by Equation (34).

$$\mathbf{C}\_{\xi} = \frac{r\_{\xi}^{p\_{\xi}} - 1}{r\_{\xi}^{2} - 1} \tag{34}$$

Following the calculation steps above, the uncertainty coefficient of C*<sup>p</sup>* and C*<sup>t</sup>* can be achieved. They are 0.129% D and 0.126% D, where D is the reference data from the NREL technical report [8]. Because both of them are less than 1% D, that proves the results are credible. Through the verification above, the level 2 actuator point number will be adopted, and it will obtain reliable results.

**Table 3.** Relationship between actuator point number and aerodynamic coefficients.

**Figure 8.** Relationship between actuator point number and aerodynamic coefficient, (**a**) power coefficient, and (**b**) thrust coefficient.

In addition, one of this method's advantages is that in the part of the aerodynamic solver it will cost less computational resources than CFD. The BEM theory is not considered here because the wake flow cannot be achieved in this way. The comparison of the computational cost between ALM and CFD method are shown in Table 4. It can be concluded that ALM uses less grids and computes faster than the CFD method under similar accuracy.

**Table 4.** Comparison of the computational cost between actuator line method (ALM) and computational fluid dynamics (CFD) method.


ALM: Actuator Line Method; CFD: Computational Fluid Dynamics; OpenFOAM: Open Source Field Operation and Manipulation, a free, open source CFD software package released by the OpenFOAM Foundation, which was incorporated as a company limited by guarantee in England and Wales; Star CCM+: Siemens Digital Industries Software, 5800 Granite Parkway, Suite 600, Plano, TX, USA.

#### *3.2. Comparison of the Power and the Thrust*

In this part, the results of aerodynamic performance are compared with 8 cases, in which 7 cases are the previous research and 1 case is the present one. Case 1 is the simulation of EALM and it is the present study. Case 2 is from the official technical report even though some data are given by NREL's FAST code [8]. In this report, the blade mode is derived by the mode's program and then passes a best-fit polynomial to get the equivalent polynomial representations of the mode shapes needed by FAST. Case 3 is simulated by HAWC2 (Horizontal Axis Wind turbine simulation Code 2nd generation) [29], which is an in-house nonlinear aeroelastic model developed by Technical University of Denmark. BEM is used as aerodynamic model and the multi-body dynamics method (MBD) is its structural model, in which each body is a linear Timoshenko beam element. The data of Case 4 are calculated by Li et al. [17] using CFD and MBD. This approach uses a dynamic overset CFD code for aerodynamics and MBD code for motion responses. The coupled way is done by exchanging the information between the fluid and structure solver in explicit form. The results of Case 5 are given by Jeong et al. [30] with BEM. In their research, the FSI (Fluid–Structure interaction) applies a strong coupling method which is using a first order implicit-explicit coupling scheme. Case 6 is set up in the research of Ponta et al. [31] by BEM and dynamic rotor deformation model, where the effects of rotor deformation are incorporates in the computation of aerodynamic loads. Yu et al. [32,33] combined CFD and CSD (Computational Structural Dynamics) in Case 7. The coupling methodology in this research is made by adopting the delta-airload loose-coupling technology. The last Case 8 is the results of the actuator line finite-element beam method (ALFBM) by Ma et al. [13]. The process of this research is three parts: the CFD solver given by OpenFOAM, the aerodynamic solver calculated by ALM, and the structural solver by finite-element beam method. The coupled steps are reading the velocity in the flow field from *t*–Δ*t* and adding it to source term, which is computed by aerodynamic solver and structural solver. To clearly compare the varieties, the detailed information about the solved method on the aerodynamic performance and the structural responses in different cases are given in Table 5. In this research, BEM theory cannot describe the detailed flow field, and the CFD method may cost too much time to calculate. Therefore, ALM is chosen to resolve aerodynamic performance. Compared to ALFBM, this research concentrates on the variation of the aeroelastic characteristics in wake flow caused by upstream wind turbine.



ALM: Actuator Line Method; BEM: Blade-Element Momentum; MBD: Multi-Body Dynamics; CFD: Computational Fluid Dynamics; CSD: Computational Structural Dynamics; FEM: Finite Element Method.

Table 6 lists the results of 8 cases with velocities of 8 m/s and 11.4 m/s, and their mean output power and thrust differences are described in a histogram to visualize disparity of these cases, as shown in Figures 9 and 10. According to the results, the NREL's reported results are chosen as the reference data. It can be seen that all the results of power are higher than the reference data, except Case 6. All the results of thrust are lower than the FAST result. The big difference of thrust results are found from Case 5 to Case 8. The main reason for this is that the tip loss correction is not considered, and it has a great influence on aerodynamic performance. Case 4 has unsatisfied power result and there is no data in Case 3 under the rated situation. Besides, Case 1 has the most approximate result to the reference data among all studies. Its power difference is less than 50 kW and thrust difference is not more than 50 kN.


**Table 6.** Comparison of the power and the thrust with different research.

Case 1: Elastic Actuator Line Model (present study); Case 2: NREL's Fatigue, Aerodynamics, Structures, and Turbulence (FAST) code; Case 3: Horizontal Axis Wind turbine simulation Code 2nd generation results; Case 4: research by Li et al. [17]; Case 5: research by Jeong et al. [30]; Case 6: research by Ponta et al. [31]; Case 7: research by Yu et al. [32,33]; Case 8: research by Ma et al. [13].

**Figure 9.** Comparison of different cases results in 8 m/s: (**a**) power; and (**b**) thrust Case 1: Elastic Actuator Line Model (present study); Case 2: NREL's FAST code; Case 3: Horizontal Axis Wind turbine simulation Code 2nd generation results; Case 4: research by Li et al. [17]; Case 5: research by Jeong et al. [30]; Case 6: research by Ponta et al. [31]; Case 7: research by Yu et al. [32,33]; Case 8: research by Ma et al. [13].

**Figure 10.** Comparison of different cases results in 11.4 m/s: (**a**) power; and (**b**) thrust Case 1: Elastic Actuator Line Model (present study); Case 2: NREL's FAST code; Case 3: Horizontal Axis Wind turbine simulation Code 2nd generation results; Case4: research by Li et al. [17]; Case 5: research by Jeong et al. [30]; Case 6: research by Ponta et al. [31]; Case 7: research by Yu et al. [32,33]; Case 8: research by Ma et al. [13].

To further examine the accuracy of EALM under different wind speed conditions, power and thrust are computed and compared with Case 2, which are shown in Figure 11. Besides, their relative errors are given in Table 7. The errors are smaller when the wind speed is close to the rated speed, which is 11.4 m/s, than those when the wind speed is below 8 m/s. This phenomenon is caused by the coefficients of lift and drag *Cl* and *Cd* in Equation (1). These two parameters are referenced from airfoil data tables and these tables are achieved by physical experiment which is given by NREL official report [8]. In that report, the airfoil data is only obtained under nearly rated wind speed. Therefore, the lift and drag of a blade element can get more accurate if the wind speed closed to 11.4 m/s. On the contrary, the result of the force in the blade element will not be satisfactory if the wind speed is far away from 11.4 m/s. In case of application that the airfoil data corresponds to the wind speed, the result will be perfect, which will be improved in future research. From results of Picture 11, the predict results show good agreements with NREL's reference data. In addition, the differences between NREL

and EALM are getting small with the inlet wind velocity increasing whether the output power or thrust. Especially in the rated situation, the relative errors of both power and thrust are less than 5%.

**Figure 11.** Comparisons of the aerodynamic performance results between elastic actuator line model (EALM) and NREL's FAST code at different wind speeds: (**a**) thrust; and (**b**) power.

**Table 7.** Comparisons of the power and the thrust relative error between elastic actuator line model (EALM) and NREL's FAST code at different wind speeds.


According to the validation above, it could be concluded that EALM can predict the power and thrust accurately. This conclusion can be predictable because the aerodynamic calculations of EALM are based on blade element theory. In this theory, the lift and drag of blade element are achieved from two-dimensional airfoil data, which is obtained from physical experiments.

#### *3.3. Comparison of the Tip Displacement*

In this section, the structural responses of blade tip are compared in 5 cases. Case 1 is the present research and obtained by EALM. Case 2 is from the official technical report and its data are given by NREL's FAST code [8]. Case 3 is the results of CFD and multi-body dynamics method by Li et al. [17]. Case 4 is set up by Yu et al. [32,33], in which CFD and CSD (Computational Structural Dynamics) are combined to calculate the coupled problem. The last Case 5 is using the actuator line finite-element beam method (ALFBM) by Ma et al. [13]. Detailed information in different cases can refer to Table 5.

Similarly, Table 8 lists the tip displacement results of 5 cases with 8 m/s and 11.4 m/s, and their differences are provided in the form of histograms, which are shown in Figures 12 and 13. Additionally, the NREL's reported results are chosen as the reference data. In these pictures, 0 direction means out of the rotor plane and 1 direction is in the rotor plane. The results of Case 3 and Case 5 are larger than those in Case 2, while the results of Case 4 are smaller. At the same time, the tip displacement of Case 1 is the nearest with the reference data in 0 direction and the farthest away from that reference data in 1 direction among 5 cases. It indicates that the tip displacement calculated by the beam solver, used in actual design, is similar to the modal approach, used to achieve structural dynamic characteristics, in 0 direction. Besides, the means used in this research could only predict approximate results in 1 direction. Even so, this beam solver can be still utilized to provide structural responses because of little deflection in edgewise.


**Table 8.** Comparison of tip deflection with different research.

Case 1: Elastic Actuator Line Model (present study); Case 2: NREL's FAST code; Case 3: research by Li et al. [17]; Case 4: research by Yu et al. [32,33]; Case 5: research by Ma et al. [13].

**Figure 12.** Comparison of different cases tip displacement results in 8 m/s: (**a**) 0 direction; and (**b**) 1 direction Case 1: Elastic Actuator Line Model (present study); Case 2: NREL's FAST code; Case 3: research by Li et al. [17]; Case 4: research by Yu et al. [32,33]; Case 5: research by Ma et al. [13].

**Figure 13.** Comparison of different tip displacement results in 11.4 m/s: (**a**) 0 direction; and (**b**) 1 direction Case 1: Elastic Actuator Line Model (present study); Case 2: NREL's FAST code; Case 3: research by Li et al. [17]; Case 4: research by Yu et al. [32,33]; Case 5: research by Ma et al. [13].

The tip displacement and its relative error in varied wind speeds are compared with Case 2 in Figure 14 and Table 9. The predict tendency of the tip displacement in 0 direction is almost the same as the reference data from Figure 13a, and its relative error is less than 5%. That indicates once again that the result of beam solver used in this study is approximate to that of the modal approach in 0 direction. Although the predicted value is not satisfied in high wind velocity in 1 direction, the relative differences are still not more than 15% except in the velocity of 10 m/s.

**Figure 14.** Comparison of the tip displacement results between elastic actuator line model (EALM) and NREL's FAST code at different wind speeds: (**a**) in 0 direction; and (**b**) in 1 direction.

**Table 9.** Comparison of the tip displacement relative errors between elastic actuator line model (EALM) and NREL's FAST code at different wind speeds.


To represent the exactitude of blade tip deformation in different positions when considering the rotating motion, the azimuthal variations of tip displacement are compared with FAST code in Figures 15 and 16. In addition, the Pearson simplified correlation coefficient *r* is introduced to describe the fit degree of two curves. It is a measure of the linear correlation between two variables in statistics. It has a value from −1 to +1, where 1 is total positive linear correlation, 0 is no linear correlation, and −1 is total negative linear correlation. Moreover, the larger the absolute value of the correlation coefficient is, the stronger proximity of two variables becomes. In this study, the coefficient *r* is rearranged in the formulas like Equation (35).

$$r = \frac{n\sum\_{i=1}^{n} x\_i y\_i - \sum\_{i=1}^{n} x\_i \sum\_{i=1}^{n} y\_i}{\sqrt{n\sum\_{i=1}^{n} x\_i^2 - \left(\sum\_{i=1}^{n} x\_i\right)^2} \sqrt{n\sum\_{i=1}^{n} y\_i^2 - \left(\sum\_{i=1}^{n} y\_i\right)^2}}\tag{35}$$

where *n* is sample size; *xi* and *yi* are the individual sample points indexed with *i*, and they represent the data from FAST and EALM in this research.

**Figure 15.** Comparison of azimuthal variations of blade tip deformations between elastic actuator line model (EALM) and NREL's FAST code at 8 m/s: (**a**) in 0 direction; and (**b**) in 1 direction.

**Figure 16.** Comparison of azimuthal variations of blade tip deformations between elastic actuator line model (EALM) and NREL's FAST code at 11.4 m/s: (**a**) in 0 direction; and (**b**) in 1 direction.

From the Figure 15, the results agree well with each other. Moreover, it can be seen that all the Pearson correlation coefficients are not less than 99% from Table 10, which means the curves of two methods have strong relativity. Consequently, it could be concluded that EALM can calculate reliable structural responses.


**Table 10.** Comparison of blade tip displacement Pearson simplified correlation coefficient between elastic actuator line model (EALM) and NREL's FAST code.

#### **4. Aeroelastic Performance Analysis**

#### *4.1. Influence on Aerodynamic Performance*

This section will discuss the influence of elastic blade on the aerodynamic performance. Figures 17 and 18 depict the power and thrust azimuthal history with traditional actuator line model (ALM) and elastic actuator line model (EALM) at 8 m/s and 11.4 m/s. The ALM results are achieved by hiding the part of rotating beam solver in EALM, and the parameters setting are all the same. The curves fluctuate every 60 degrees, which is caused by the tower shadow effect. When the structural deformation is considered, the mean power reduces about 11.38 kW in 8 m/s and 11.43 kW in 11.4 m/s. Besides, the averaged thrust decreases approximately 1.08 kN in 8 m/s and 0.35 kN in 11.4 m/s. It indicates that the influence of elasticity on mean power and thrust is small in the rated condition. As the blade passes through tower, the reduction of output power turns into 15.33 kW in 8 m/s and 31.45 kW in 11.4 m/s. The decrement of thrust changes into 1.32 kN in 8 m/s and 0.96 kN in 11.4 m/s. The distance between tower and blade gets closer due to the blade deformation. Hence, the tower shadow effect becomes more serious and the decrement of power and thrust is larger. Because the blade displacement in high wind speed is large, the tower effect is the most serious in the rated wind speed and its reduction is more than that in 8 m/s.

**Figure 17.** Comparisons of the aerodynamic performance results between traditional actuator line model (ALM) and elastic actuator line model (EALM) at 8 m/s: (**a**) power; and (**b**) thrust.

**Figure 18.** Comparisons of the aerodynamic performance results between traditional actuator line model (ALM) and elastic actuator line model (EALM) at 11.4 m/s: (**a**) power; and (**b**) thrust.

The stabilization of the normal and tangential force is present by standard deviation, as shown in Figure 19. The standard deviation s is defined as follows:

$$s = \sqrt{s^2} = \sqrt{\frac{\sum\_{i=1}^{n} (x\_i - \overline{x})^2}{n-1}} \tag{36}$$

where *n* is sample size; *xi* is the individual sample points indexed with *I*; and *x* is expectation, which is the mean value of sample. From the figure, the standard deviation of the loads on blade increases a little when the elasticity of the blades is considered.

**Figure 19.** Stabilization of forces on blade: (**a**) normal force; and (**b**) tangential force.

Figure 20 shows the wake structure comparison with EALM and ALM at 8 m/s. The wake vorticities are described by Q criterion, in which Q is calculated by the second invariant of the velocity gradient tensor and are stained by velocity field. The tip and hub vorticities are clear to be seen behind the rotor. The velocity in the internal surface of vorticities is smaller than the outside one. That is because the wind turbine extract momentum from the incoming airflow passing through it, and the wind energy is turned into mechanical energy, which causes the wind velocity after the rotor decreased.

**Figure 20.** The wake structures in terms of Q field at 8 m/s: (**a**) by EALM; and (**b**) by ALM.

The blade deflection is abbreviated and described in Figure 21. The deformation is too small to be discovered in 1 direction from left side view. Besides, it focuses on blade tip in 0 direction from front and vertical side view.

**Figure 21.** The abbreviated sketches of blade deflection: (**a**) front side view; (**b**) left side view; (**c**) vertical side view; and (**d**) southwest isometric side view.

#### *4.2. Wake Flows Analysis*

To analyze the influence of elastic blade in the wake flows, an array of two NREL 5 MW wind turbines is studied (see Figure 7) with EALM. The results are extracted after the upstream wind turbine rotating 40 revolutions, and the corresponding simulation time is around 260 s. According to research [28], the rotor speed of downstream wind turbine is set as 8.11 rpm to maximize its output power.

The comparison of mean output power among EALM, ALM, and Jha et al. [34] is presented in Table 11, in which the power generated by downstream wind turbine reduces a lot and only about 40% of upstream wind turbine because of the wake flow effects. In addition, the elastic blade makes the mean power of downstream wind turbine decrease about 14.6 kW, which is bigger than 11.38 kW of upstream wind turbine. Furthermore, the power and thrust azimuthal history of downstream wind turbine are depicted in Figure 22. When the blade passes through tower, the output power reduces 24.83 kW and the thrust decreases 7.01 kN. In contrast with the reduction of 8 m/s in Section 4.1, which are 15.33 kW and 1.32 kN, respectively, the tower shadow effect becomes more serious. This phenomenon is caused by the instability of the flow filed around the blade. The velocity and turbulence intensity in the wake of upstream wind turbine are changed into unstable. In the flow field of the downstream wind turbine blade, the variety of blade deflection with time will further disturb them. That will impact the relative velocity of the blade element. Thus, it can be concluded that the decrement of power and thrust caused by elastic blade are getting larger in the wake flows according to the comparisons above.



**Figure 22.** Comparison of downstream wind turbine aerodynamic performance results between Table 8. m/s: (**a**) power; and (**b**) thrust.

The stabilization of the normal and tangential force on downstream wind turbine blade are compared again in Figure 23. In contrast with Figure 19, the standard deviation of the loads on blade in the wake caused by upstream wind turbine getting larger when the elasticity of the blades is considered. It indicates that the influence of elastic blade becomes more serious when it is in the wake of upstream wind turbine. That will aggravate the fatigue loads and should be paid much more attention during corresponding studies.

**Figure 23.** Stabilization of forces on downstream wind turbine blades: (**a**) normal force; and (**b**) tangential force.

Figure 24 illustrates the process of double wind turbines array vorticity development. The wake expands as expected and the tip vorticity turns into continuum gradually in Figure 24a. Then the vorticity breakdown and the smaller-scale turbulence appears with distance increasing from upstream wind turbine (referring to Figure 24b). When the downstream wind turbine meets with this complicated wake, where the turbulence levels raise because of upstream wind turbine wakes and the momentum has not recovery completely, its tip vorticities break down earlier. That makes the wakes behind downstream wind turbine more complex, which can be seen in Figure 24c.

**Figure 24.** The wake structure development of double wind turbine aligned: (**a**) initial status without disturb; (**b**) developing status; and (**c**) interference status.

Figure 25 shows the development of wake velocity field correspondingly. The phenomenon of velocity deficit is easy to observe behind the rotor. Before the downstream wind turbine disturbed from the wakes generated by upstream wind turbine, the velocity field is stable relatively. As the interference status appears, the distribution of velocity field in nearby regions of downstream wind turbine is in a muddle condition. That will increase the fatigue loads on blade of downstream wind turbine.

*Water* **2020**, *12*, 1233

**Figure 25.** The wake velocity field development of double wind turbine aligned: (**a**) initial status without disturb; (**b**) developing status; and (**c**) interference status.

#### **5. Conclusions and Discussions**

A new elastic actuator line model, which combines the traditional actuator line model and a beam solver used in the blade design, is proposed to study the aeroelastic performance of wind turbine efficiently. The validation is set up through comparing the result of the power, the thrust, and the tip deflection with different research. It is found that this new model can obtain acceptable prediction, including aeroelastic performance and wake fields. Besides, the tip displacement in 0 direction is the nearest with the reference data from official technical report given by NREL's FAST code among these studies. Furthermore, the elasticity of the blades can aggravate tower shadow effect and has the possibility of wind turbine instability. The fatigue loads on the blades also become more serious. When the wake effect generated by upstream wind turbine is considered, the impact above will get more obvious.

However, this new elastic actuator line model is incomplete and needs to be improved. For example, the modal shape and natural frequency results are not yet considered. The purpose of this study is to propose an approach for investigating the aeroelastic performance rapidly. Further study is still required to refine the solver of structural dynamic characteristics and wave interactions [35].

**Author Contributions:** Z.Y. made the computations, conducted data analysis; Z.H. made the comparison of calculation results; X.Z. gave the suggested ideas and did the proof reading; Q.M. guided the whole project; H.H. gave the data analysis of double wind turbine wake flows. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** This research work was supported by the National Natural Science Foundation of China (Nos. 51879051; 51739001; 51639004). The fourth author also thanks the Chang Jiang Visiting Chair Professorship scheme of the Chinese Ministry of Education, hosted by HEU.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Water* Editorial Office E-mail: water@mdpi.com www.mdpi.com/journal/water

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-03936-803-7