**The Numerical Simulation of Fluid Flow**

Editors

**Robert Castilla Anton Vernet**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editors* Robert Castilla Fluid Mechanics Universitat Politecnica de ` Catalunya Terrassa Spain

Anton Vernet Mechanical Engineering Universitat Rovira i Virgili Tarragona Spain

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Energies* (ISSN 1996-1073) (available at: www.mdpi.com/journal/energies/special issues/ Numerical Simulation Fluid Flow).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-2931-8 (Hbk) ISBN 978-3-0365-2930-1 (PDF)**

© 2022 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **About the Editors**

#### **Robert Castilla**

Robert Castilla graduated in Physics from the Universitat Autonoma de Barcelona (Spain) in ` 1991 and obtained a Ph.D. in Applied Physics from the Universitat Politecnica de Catalunya (Spain) ` in 2001. He received a position as Associate Professor in the Department of Fluid Mechanics in the Universitat Politecnica de Catalunya, in Terrassa, in 2011. He is also a member of the ` Steering Committee of CATMech, Center of Mechanical Technologies of Catalonia. His main interest is in Computational Fluid Dynamics with diverse applications: incompressible complex flows (fluid power components as pumps, valves, cylinders, etc.), compressible flows, turbulent flows, aeroacoustics, etc. The present projects are related to supersonic flow in vacuum ejectors and aeroacoustics in cars.

#### **Anton Vernet**

Anton Vernet graduated in Chemistry from the Universitat de Barcelona (Spain) in 1991 and obtained a Ph.D. in Chemical Engineering from the Universitat Rovira i Virgili (Spain) in 1997. After one year working as a researcher at the University of Western Ontario (Canada), he obtained a position as associate professor in the area of Fluid Mechanics at the Universitat Rovira i Virgili. He is currently working as a director of the Master in Computational Fluid Mechanics taught jointly between the Universitat Rovira i Virgili and the Universidad Internacional de la Rioja. He is part of the ECoMMFiT research group. His basic lines of research are focused on data analysis, pattern recognition, fuzzy logic and the analysis of images, mainly applied to the study of turbulent flows, heat transfer and biomedical flows. He also works experimentally in the development and application of PIV (Velocimetry by particle imaging) techniques to the study of turbulent flows.

## **Preface to "The Numerical Simulation of Fluid Flow"**

Almost any energy production process involves fluid flow. From the most obvious like air through the wind turbine blades or fuel flow in an internal combustion engine, to a secondary, though still essential, role, like lubrication in a mechanical power transmission hub. In the last decades, the simulation of fluid flows has been so relevant that CFD (Computational Fluid Dynamics) has become a discipline that is included in any textbook of Fluid Mechanics. The growth of computing capacity, summarized in Moore's law, and the development of numerical methodologies provide increasingly efficient and accurate simulations.

This book, which collects the articles published in a Special Issue of *Energies*, aims to be more focused on practical application of available methodologies and models rather than presentation of new numerical methods.

The selected papers for this publication cover a wide variety of topics related to fluid flow numerical simulation. They present problems related to aerodynamics, aeroacoustics and fluid-particle interaction. The ultimate goal has been to show the great diversity of methodologies and procedures that are being used in the numerical solution of flow problems. Some of the papers present innovative methodologies developed for a particular application. Chiu et al. present a solver that allows analyzing the flow-particle interaction. They compare the results obtained with those already published in cases such as Flow Past Circular Cylinders in Tandem or Lid-Driven Semi-Circular Cavity Flow, using them as a benchmark. Gorakifard et al. use the Lattice-Boltzman method to analyze an aeroacoustics problem in wind turbines. They describe the propagation of planar acoustic waves, including the temporal decay of a standing plane wave and the spatial decay of a planar acoustic pulse. The analysis of these specific benchmark problems has yielded qualitative and quantitative data on acoustic dispersion and dissipation, and their deviation from analytical results demonstrates the accuracy of the method. Yuan et al. also address the issue of wind turbines. In this case, they propose a new two-dimensional numerical method, which can facilitate and accelerate the design of the array. Bimbato et al. analyzes a turbulent boundary layer flows past a bluff body and the effect of surface roughness using the Lagrangian vortex method in association with Large-Eddy Simulation (LES). The problem of the effect of roughness is also discussed in the article by de Oliveira al. Xu and Wang use a finite difference (FD) method of monopole source to simulate the response of full-wave acoustic-logging in cave formations. The results of their study show that borehole full-waves have response characteristics to the caves, which are mainly reflected in the energy and the first arrival of waves. Finally, the open-source code OpenFoam is used by Guerrero and Castilla to evaluate the aerodynamic performance of a Formula 1 car under the influence of a turbulent wake.

> **Robert Castilla, Anton Vernet** *Editors*

## *Article* **A Fast Two-Dimensional Numerical Method for the Wake Simulation of a Vertical Axis Wind Turbine**

**Zheng Yuan <sup>1</sup> , Jin Jiang 2,\*, Jun Zang 3,\*, Qihu Sheng <sup>1</sup> , Ke Sun <sup>1</sup> , Xuewei Zhang <sup>1</sup> and Renwei Ji <sup>1</sup>**


**\*** Correspondence: jiangjin@jit.edu.cn (J.J.); j.zang@bath.ac.uk (J.Z.)

**Abstract:** In the array design of the vertical axis wind turbines (VAWT), the wake effect of the upstream VAWT on the downstream VAWT needs to be considered. In order to simulate the velocity distribution of a VAWT wake rapidly, a new two-dimensional numerical method is proposed, which can make the array design easier and faster. In this new approach, the finite vortex method and vortex particle method are combined to simulate the generation and evolution of the vortex, respectively, the fast multipole method (FMM) is used to accelerate the calculation. Based on a characteristic of the VAWT wake, that is, the velocity distribution can be fitted into a power-law function, a new correction model is introduced to correct the three-dimensional effect of the VAWT wake. Finally, the simulation results can be approximated to the published experimental results in the first-order. As a new numerical method to simulate the complex VAWT wake, this paper proves the feasibility of the method and makes a preliminary validation. This method is not used to simulate the complex three-dimensional turbulent evolution but to simulate the velocity distribution quickly and relatively accurately, which meets the requirement for rapid simulation in the preliminary array design.

**Citation:** Yuan, Z.; Jiang, J.; Zang, J.; Sheng, Q.; Sun, K.; Zhang, X.; Ji, R. A Fast Two-Dimensional Numerical Method for the Wake Simulation of a Vertical Axis Wind Turbine. *Energies* **2021**, *14*, 49. https://dx.doi.org/ 10.3390/en14010049

Received: 14 November 2020 Accepted: 20 December 2020 Published: 24 December 2020

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 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 (https://creativecommons.org/ licenses/by/4.0/).

**Keywords:** vertical axis wind turbine (VAWT); two-dimensional wake simulation; finite vortex method; vortex particle method; three-dimensional effect correction model of the wake

#### **1. Introduction**

Since the energy crisis and global warming are getting more serious, the development of renewable energy has received increasing attention. As common renewable energy, wind energy has been well developed due to its wide distribution and nonpollution. A review published by the Intergovernmental Panel on Climate Change (IPCC) shows that the contribution of wind power to the global electricity supply in 2050 will reach 13–14% [1]. There are two common types of wind turbines: horizontal axis wind turbine (HAWT) and vertical axis wind turbine (VAWT). Because of those advantages, such as no need for a yaw control mechanism and low manufacturing costs [2], VAWT has become a research hotspot in recent years. Although the power coefficient of a single VAWT is lower than that of a HAWT, however, some previous research showed that the performance of VAWTs in array configuration would improve significantly. The power coefficient (cp) will increase in a counter-rotating VAWT pairs configuration [3], and the power density of the VAWT farms is potentially more than that of the HAWT farms [4]. To make the VAWT array design more efficient, it is necessary to calculate the velocity field quickly. Therefore, a fast simulation method of the VAWT wake is proposed in this paper so that the wake effect can be measured quickly according to the velocity deficit.

Compared with the HAWT, the study on the wake of a VAWT is still limited. Moreover, in the limited experiments on the wake-field, most of them only focus on the near wake within 3 diameter distances [5–7]. In general, the VAWT array is designed based on the

velocity field of the whole wind farm; it is not enough that only near-field data are available, and the complete wake data are necessary. For this reason, there are three requirements that need to be met to be the benchmark test case. These requirements are: First, the velocity profiles from near-field to far-field need to be included in the experimental data (It is better to be as far as 10D downstream). Second, there is no blockage effect; third, the blade is a common type. According to the current research status of the wake experiments, only three experiments [8–10] provide the velocity profiles more than 10D downstream. However, two of them do not meet the other two requirements. In this view, the experiment conducted by WIRE wind tunnel [10] is the only relatively representative experiment so far, which satisfies these three requirements at the same time, and it is chosen as the benchmark test case in this paper. The details of this experiment are clearly reported in Table 2 (Section 3.3).

It should be noted that the scale of the VAWT model in this experiment is not big; the diameter of the rotor is 16.6 cm. However, in order to reduce the blockage effect of the wind tunnel or to obtain relatively farther velocity profiles under the limited measurement distance, many experimental models used in the VAWT wake studies are normally small [11,12]. However, according to the previous studies, the velocity distribution of the wake is quite similar, working on different scale ratios [13,14] and inflow velocities [15,16]. Although the scaling effect of the small model may cause some deviation, in the current, limited research, the benchmark test case selected in this paper is still the most representative. Moreover, only when more experiments of the VAWT wake are published in the future, more validation can be done at that time. However, it wouldn't affect the feasibility of this method in theory and the practicality in engineering.

In this paper, the numerical methods of the VAWT wake simulation are compared and analyzed in Section 2. According to the analysis of each method and the requirement for rapid simulation in array design, the new numerical method is proposed and introduced in detail in Section 3, which includes the finite vortex method, the vortex particle method, the corresponding convergence criteria, and the three-dimensional effect correction of the velocity field. Finally, the numerical results are compared and validated with the published PIV experiment, and the characteristics of the VAWT wake are analyzed in Section 4.

#### **2. Comparison and Analysis of Numerical Methods for Wake Simulation**

Because the current wake simulation methods of the vertical axis wind turbine (VAWT) are the same as that of the vertical axis tidal turbine (VATT), therefore, the general numerical methods of the vertical axis turbine (VAT) are analyzed in the following. Apart from the experimental method mentioned above, there are increasing numerical methods to simulate the VAT wake in recent years. Different from the traditional classification of the research methods [17], the simulation process in this paper will be divided into three steps, and each step can be completed by different sub-methods. These three steps and the corresponding sub methods are shown in the Table 1 below:

**Table 1.** Three steps of the wake simulation and the corresponding sub methods (Note: there is a label in bold before the method, which is the abbreviation of the method and will be used in the following analysis. The number of the label is the step of this process; the letter of the label is the initial of this method).


The process of the wake simulation can be divided into three steps. In the following examples, the whole simulation process can be represented by the combination of the label mentioned above: Li [27] used the experimental data of static airfoil in the wake research, the evolution of the wake was simulated in the form of the discrete vortex with an empirical decay function (The label of this method is 1T-2V-3VE). The actuator line model of Mendoza [28] used the experimental data modified by the dynamic stall model, and the force coefficient was projected onto the grid to simulate the evolution of the flow field (The label of this method is 1T-2F-3G). There is a hybrid method [32] which was used to calculate the forces on the Euler grid, then the vorticity source term was transformed into the vortex particle in a near-wall transition region, and the far-field was simulated in the Lagrangian framework (The label of this method is 1G-2V-3VA).

Based on the compromise between accuracy and computational efficiency, the vortex panel method was applied to calculate the blade force and vortex shedding because of its affordable computational costs. As a type of the vortex panel method, the finite vortex method [19] is more applicable in a wider range of working conditions due to its improved Kutta condition. However, the simulation of a VAWT wake is a challenging task; the wake evolution of the VAWT is much more complex than that of HAWT, such as the breakup and merging behavior of the vortex swarm. Moreover, there is also a blade-vortex interaction (BVI) problem in the VAWT wake; that is, the blade will encounter the shedding vortex at the downstream half revolution. Nevertheless, compared with the grid-based method, the vortex particle method will be more competent for this job because of its low numerical dissipation and easy-to-parallel. In addition, according to the characteristics of the VAWT wake-field, high vorticity is mainly distributed in the upper and lower rows of the wake instead of the whole wake area. In contrast, the whole domain needs to be calculated in the grid-based method, while only a small region where the vorticity is relatively high needs to be calculated in the vortex particle method. In other words, the calculation domain of the vortex particle method is smaller than that of the grid-based method. Thus, it can be concluded that the evolution of the VAWT wake can be calculated quickly and accurately by the combination of the finite vortex method and the vortex particle method. The corresponding label of this method proposed in this paper is 1 V-2 V-3VA.

To save computing costs, the two-dimensional model was often used in the preliminary array design [33–35]. Even in the atmospheric boundary layer, the wake can also be simplified as the two-dimensional Gaussian wake model [36]. In addition, the aspect ratio (height/diameter) of the VAWT in real engineering practice is usually big, so that more wind energy can be captured by increasing the height without increasing the land area. Therefore, the influence of the three-dimensional effect, which is caused by the blade end, will become smaller in the case of a large aspect ratio. The above reasons make the two-dimensional numerical method more suitable in the VAWT array design.

When the two-dimensional numerical model is used in the preliminary design of the VAWT array, both the blade force and the velocity field should be modified by the threedimensional effect correction model. There are some three-dimensional effect correction models for the force and performance [37,38], but there are few three-dimensional effect correction models for the VAWT wake. By comparing the three-dimensional and twodimensional wake simulation [39], it can be found that the length of the three-dimensional VAWT wake is shorter than that of the two-dimensional VAWT wake because velocity deficit will recover in the third dimension. In order to consider the three-dimensional effect in the velocity recovery, a new correction model of the VAWT wake is proposed, which is the projection coefficient from the two-dimensional to three-dimensional along the streamwise direction.

It is true that the turbulent evolution of the VAWT wake is a three-dimensional problem. However, the purpose of this paper is to propose a numerical method that can calculate the velocity distribution quickly and make the VAWT array design more efficient; the complex details of the three-dimensional turbulence evolution can be ignored. Therefore, it is not necessary to pay too much attention to the turbulent problems such as the dynamic stall and the turbulent inflow conditions because these turbulence details are not important for the preliminary array design and will be covered by the three-dimensional effect correction model.

In general, according to the tradeoff between accuracy and efficiency, a two-dimensional numerical method with necessary three-dimensional effect correction is proposed; it can calculate the wake-field quickly and relatively accurately, the first-order approximation of the velocity distribution is sufficient, which meets the requirement in the preliminary array design. More details of this method will be introduced in Section 3.

#### **3. Methodology**

The whole process of this new method is like this: First, the finite vortex method is used to calculate the blade load and vortex shedding. Then, the shedding vortex is discretized into the vortex particles to simulate the wake evolution. Finally, it is necessary to introduce a three-dimensional effect correction model into the two-dimensional velocity field so that the modified simulation results can be used in the preliminary array design. Details are given in the following subsections.

#### *3.1. Finite Vortex Method*

The forces and circulation of the blade can be calculated accurately by the vortex panel method, which was confirmed in the published papers [18,20]. The finite vortex method is one type of the vortex panel method [19,40]. A schematic of the vertical axis turbine in the reference coordinate is shown in Figure 1. A close-up view of the blade is shown in Figure 2. The motion of the blade can be described by the translational velocity *U* and angular velocity Ω relative to the local coordinate system *o*. The inflow velocity is *VA*. There are *N* sources and sinks distributed on the surface of the blade and a linear vortex lineγon the chord line. The schematic diagram is shown in Figure 3. *Ω* γ *Ω* γ

**Figure 1.** Reference coordinate system.

**Figure 2.** Local coordinate system.

*Ω*

γ

**Figure 3.** The schematic diagram of the finite vortex method. Subscript *i* is the serial number of source and sink singularities, *p* is the number of blades, *f* represents the attached vortex on the blade, *w* represents the shedding vortex in the wake, superscript *k* is the time step.

The total velocity potential *φ* satisfies the Laplace equation in the global coordinate system, as shown in Equation (1) below:

$$
\nabla^2 \phi(\mathbf{X}, \mathbf{Y}, t) = \mathbf{0} \tag{1}
$$

The surface of the blade satisfies no through-flow boundary condition, as shown in Equation (2), where *r* is the vector from the local coordinate system *o* to the blade surface, *nb* is the unit normal vector that points to the interior of the blade, the symbol *S<sup>b</sup>* denotes the surface of the blade.

$$\frac{\partial \phi}{\partial n}|\_{S\_b} = \left(\stackrel{\rightarrow}{\mathcal{U}} - \stackrel{\rightarrow}{V}\_A + \stackrel{\rightarrow}{\mathcal{\Omega}} \times \stackrel{\rightarrow}{r}\right) \cdot \stackrel{\rightarrow}{n\_b} \tag{2}$$

The symbol *p* is the field point; in the far-field, the induced velocity potential *ϕ* satisfies lim *<sup>p</sup>*→<sup>∞</sup> ∇*ϕ* = 0. At the initial time *t* = 0, the induced velocity potential satisfies the relationship shown in Equation (3):

$$
\nabla \varphi(p, 0) = 0 \tag{3}
$$

Some experimental studies [41,42] on the flow field showed that there is a finite pressure difference at the trailing edge under the unsteady conditions, which makes the streamline roll-up. In order to solve the stall condition at a high angle of attack, the improved Kutta condition set up a finite pressure difference at the trailing edge. Therefore it is called the finite vortex method [19,43]. The subscripts *u* and *d* represent the up and downsides of the trailing edge. The improved Kutta condition of the finite vortex model is shown in Equation (4):

$$p\_u - p\_d = \Delta p \tag{4}$$

In order to determine the value of this pressure difference, set a parameter *l* to satisfy the relation: <sup>∆</sup>*<sup>p</sup>* = *<sup>p</sup><sup>u</sup>* − *<sup>p</sup><sup>d</sup>* = −*lργ* (*k*) *<sup>w</sup>* , where *γ* (*k*) *<sup>w</sup>* is the circulation of the vortex shedding at the step *k*. By adjusting the value of *l*, the lift coefficient of numerical results would be kept consistent with the lift coefficient obtained by the experiment at different angles of attack [19]. The value of *l* is close to 0 at the small angle of attack, which indicates that the computational error of conventional Kutta condition (pressure difference at trailing edge is equal to 0) is not large at the small angle of attack. With the increase of angle of attack, the absolute value of *l* is also increasing, which indicates that the pressure difference of the trailing edge is also increasing gradually. Therefore, due to the existence of the correction factor *l* and the finite value of pressure difference at the trailing edge, the aerodynamic loads calculated by the finite vortex method can still be consistent with the experimental data at a high angle of attack. It also means that the finite vortex method can be more accurate than the conventional vortex panel method in the case of a high angle of attack.

According to the unsteady Bernoulli equation (as shown in Equation (5)), the pressure in the flow field is:

$$\frac{p - p\_{\infty}}{\rho} = -\frac{\partial \rho}{\partial t} + \left(\stackrel{\rightarrow}{\mathcal{U}} - \stackrel{\rightarrow}{V}\_{A} + \stackrel{\rightarrow}{\Omega} \times \stackrel{\rightarrow}{r}\right) \cdot \nabla \varphi - \frac{1}{2} (\nabla \varphi)^{2} \tag{5}$$

As the up and downsides of the trailing edge are very close to each other, the vector of the trailing edge can be written as Equation (6):

$$
\stackrel{\rightarrow}{r\_u^{\rightarrow}} \approx \stackrel{\rightarrow}{r\_d^{\rightarrow}} \approx \frac{1}{2} \left( \stackrel{\rightarrow}{r\_u^{\rightarrow}} + \stackrel{\rightarrow}{r\_d^{\rightarrow}} \right) = \stackrel{\rightarrow}{r\_{\text{TE}}} \tag{6}
$$

Therefore, after the linearization, the Kutta condition can be rewritten as Equation (7):

$$\begin{split} \frac{\partial(\varrho\_{u} - \varrho\_{d})}{\partial t} &= \frac{\partial \mathbf{r}\_{f}}{\partial t} = \left( \overset{\rightarrow}{\mathcal{U}} - \overset{\rightarrow}{V}\_{A} + \overset{\rightarrow}{\mathcal{D}} \times \overset{\rightarrow}{r\_{TE}} \right) \cdot \left( \nabla \varrho\_{u} - \nabla \varrho\_{d} \right) + \frac{1}{2} (\nabla \varrho\_{d} - \nabla \varrho\_{u}) \cdot \left( \nabla \varrho\_{d} + \nabla \varrho\_{u} \right) - \frac{\Delta p}{\rho} \\ &= \left( \overset{\rightarrow}{\mathcal{U}} - \overset{\rightarrow}{V}\_{A} + \overset{\rightarrow}{\mathcal{D}} \times \overset{\rightarrow}{r\_{TE}} - \overset{\rightarrow}{V\_{TE}} \right) \cdot \left( \nabla \varrho\_{u} - \nabla \varrho\_{d} \right) - \frac{\Delta p}{\rho} \end{split} \tag{7}$$

where mean disturbance velocity vector *VTE* of the trailing edge is shown in Equation (8):

$$
\stackrel{\rightarrow}{V}\_{TE}^{\rightarrow} = \frac{1}{2} (\nabla \varphi\_{\mu} + \nabla \varphi\_{d}) \tag{8}
$$

According to Kelvin's theorem, the total circulation does not change with time under the potential flow theory; in other words, the sum of circulation of the bound vortex (Γ<sup>f</sup> ) and free vortex (Γw) does not change with time. Hence, there is Equation (9). At different azimuthal angles, the lift and circulation of the airfoil will change. In each time step, the circulation variation of the bound vortex and that of the free vortex are equal in value, but the sign of them are opposite. It can be expressed as Equation (10). Here, the subscript *f* represents the circulation of the foil, and the subscript *w* represents the vortex in the wake, the superscript *k* represents the time step. Combining Equation (10) with Kelvin's theorem (Equation (9)), Equation (11) is obtained. Then, the circulation of the vortex shedding γ can be calculated by inserting Equation (11) into Equation (7) above. Moreover, the γ obtained here is exactly the circulation that will be discretized into the vortex particles. The detailed process is shown in the references [18–20].

$$\frac{d\Gamma\_{\text{total}}}{dt} = \frac{d\left(\Gamma\_f + \Gamma\_w\right)}{dt} = 0 \tag{9}$$

$$
\Gamma\_f^{(k)} - \Gamma\_f^{(k-1)} = -\gamma\_w^{(k)} \tag{10}
$$

$$\frac{d\Gamma\_f}{dt} = -\frac{d\Gamma\_w}{dt} = \frac{\Gamma\_f^{(k)} - \Gamma\_f^{(k-1)}}{\Delta t} = -\frac{\gamma\_w^{(k)}}{\Delta t} \tag{11}$$

#### *3.2. Vortex Particle Method*

The circulation of the shedding vortex is calculated by the finite vortex method in the potential flow. Then, it will be discretized into a vortex particle swarm to calculate the wake evolution in viscous flow. The vortex particle program is divided into six subsections as follows, and the definition of symbols in each equation is only applicable to the relevant subsection.

#### 3.2.1. Discretization of the Circulation

In this paper, the Vatistas model [44] is used as the vortex core model, which is a common smoothing function to remove the singularity, and it is also similar to the Gaussian smoothing function. The Vatistas model also has been applied in the research of a HAWT wake [13]. The regularization function *K* is shown in Equation (12): (where *ρ* = *r*/*r<sup>c</sup>* and *r<sup>c</sup>* is the radius of the vortex core):

$$K(\rho) = \frac{1}{\pi(\rho^4 + 1)^{\frac{3}{2}}} \tag{12}$$

Which satisfies the normalization condition, as shown in Equation (13):

$$2\pi \int \mathcal{K}(\rho)\rho d\rho = 1\tag{13}$$

The function of vorticity distribution is shown in Equation (14):

$$
\omega(r) = \frac{\Gamma}{r\_c^2} K(\rho) = \frac{\Gamma r\_c^4}{\pi (r^4 + r\_c^4)^{\frac{3}{2}}} \tag{14}
$$

#### 3.2.2. Vorticity Transport Equation

The vortex particle method is a viscous method, the governing equation as follows is the vorticity transport equation, which is obtained by taking the curl of the incompressible NS equation, as shown in Equation (15):

$$\frac{d\omega}{dt} = (\omega \cdot \nabla)V + \nu \nabla^2 \omega \tag{15}$$

It is worth mentioning that in the two-dimensional case, there is no stretching term, so the control equation is simplified as Equation (16):

$$\frac{d\omega}{dt} = \nu \nabla^2 \omega \tag{16}$$

It can be seen that the governing equation of the vortex particle method in the twodimensional case is much simpler than the NS equation of the grid-based method.

#### 3.2.3. The Calculation of the Velocity of Vortex Particles

The convection of the vortex particles is calculated by the Biot–Savart Equation. The induced velocity of the Vatistas model in Equation (17):

$$
\stackrel{\rightarrow}{V} = \frac{\stackrel{\rightarrow}{\Gamma}}{2\pi} \times \frac{\stackrel{\rightarrow}{r}}{\sqrt{r^4 + r\_c^4}} \tag{17}
$$

It can be seen that Equation (17) is very close to the induced velocity equation of a two-dimensional singular vortex, as shown in Equation (18):

$$V = \frac{\Gamma}{2\pi r} \tag{18}$$

In this paper, the fast multipole method (FMM) [45] is introduced to accelerate the calculation of the induced velocity of the vortex particles. Here, the induced velocity can be written in the form of induced velocity potential *ϕ*, *z*<sup>0</sup> is the target position, *z* is the position of the vortex particle, and *zc* is the centroid position of the block. The derivation of the Green function is shown in Equation (19):

$$G(z\_0, z) = -\frac{1}{2\pi} \log(z\_0 - z) = -\frac{1}{2\pi} \left[ \log(z\_0 - z\_c) + \log\left(1 - \frac{z - z\_c}{z\_0 - z\_c}\right) \right] \tag{19}$$

Applying the Taylor series expansion log(1 − *ξ*) = − ∞ ∑ *k*=1 *ξ k k* to the second logarithmic term on the right side of Equation (19), Equation (19) can be transformed into the following form (Equation (20)):

$$G(z\_0, z) = \frac{1}{2\pi} \sum\_{k=0}^{\infty} O\_k(z\_0 - z\_c) I\_k(z - z\_c) \tag{20}$$

There are two auxiliary functions (*Ik*(*z* − *zc*) and *Ok*(*z*<sup>0</sup> − *zc*)) defined here: ≥0

$$\mathbf{I}\_{k}(z - z\_{c}) = \frac{\left(z - z\_{c}\right)^{k}}{k!}, \text{ for } k \ge 0 \tag{21}$$

$$\text{O}\_k(\mathbf{z}\_0 - \mathbf{z}\_\varepsilon) = \frac{(k-1)!}{\left(\mathbf{z}\_0 - \mathbf{z}\_\varepsilon\right)^{k'}}, \text{ for } k \ge 1; \text{ and } \text{O}\_0(\mathbf{z}\_0 - \mathbf{z}\_\varepsilon) = -\log(\mathbf{z}\_0 - \mathbf{z}\_\varepsilon) \tag{22}$$

The final expression of velocity potential *ϕ* (Equation (23)) and accuracy criterion of order *p* (Equation (24)) can be written as follows: (Γ is the circulation of each vortex particle, *p* is the order of Taylor expansion, *N* is the number of the vortex particles, the error *ε* is set to 10 −5 in this paper, *R* is the diagonal length of the block. Moreover, the black dot here represents the scalar multiplication.) can be written as follows: (Γ *ε*  − ( ) ( )!

$$\phi = \frac{1}{2\pi} \sum\_{k=1}^{p} \left\{ \frac{(k-1)!}{(z\_0 - z\_c)^k} \cdot \left[ \sum\_{0}^{N} \frac{(z - z\_c)^k \cdot \Gamma}{k!} \right] \right\} - \frac{1}{2\pi} \ln(z\_0 - z\_c) \cdot \left[ \sum\_{0}^{N} \Gamma \right] \tag{23}$$

$$p = \ln\left[\varepsilon \cdot (2\pi) \cdot (1 - \frac{R}{|z\_0 - z\_c|})\right] / \ln(\frac{R}{|z\_0 - z\_c|}) - 1\tag{24}$$

There is an error criterion to judge whether the FMM can replace the direct calculation. The order of the FMM program is *NlogN*; by contrast, the order of the direct calculation method is *N*<sup>2</sup> [46,47]. In the FMM program, it is necessary to allocate the whole vortex particles into a tree structure and store them in the leaf nodes of the tree structure. The schematic diagram of the fast multipole method is shown in Figure 4.

**Figure 4.** The schematic diagram of how to allocate the vortex particles into the tree structure.

#### 3.2.4. The Calculation of the Vorticity of Vortex Particles

ω When the position of the vortex particles is updated, it is necessary to update the vorticity of the vortex particles. There are two common methods to simulate the viscous dissipation. One is called the random walk algorithm; the other is the particle strength exchange method (PSE). In this paper, the latter method (PSE) is used. In order to calculate the vorticity dissipation term, a transformation Equation (25) is introduced here to transform the second derivative form of the NS equation into the integral or summation form [48] (The symbol *f* in the following equation represents the function of vorticity ω).

$$
\nabla^2 f(\mathbf{x}) = \frac{2}{\sigma^2} \int [f(y) - f(\mathbf{x})] \eta\_\sigma(\mathbf{x} - y) dy \tag{25}
$$

 

 

where:

$$\eta\_{\sigma}(r) = \frac{1}{\sigma^2} \eta\left(\frac{r}{\sigma}\right) \text{ and } \eta(\rho) = -\frac{1}{\rho} \frac{dK(\rho)}{d\rho} \tag{26}$$

The governing equation is transformed into the summation form (Equation (27)):

$$\frac{d\omega\_P}{dt} = \frac{12\upsilon}{\pi} \sum\_{q} \frac{\sigma^4 r\_{pq}^2 \left(\Gamma\_q - \Gamma\_p\right)}{\left(r\_{pq}^4 + \sigma^4\right)^{5/2}}\tag{27}$$

Subscript *p* represents the target point, *q* is the vortex particles, and the symbol *σ* is the radius of the vortex particle.

#### 3.2.5. The Calculation of the Turbulent Viscosity Coefficient

It is worth mentioning that the viscosity coefficient *ν* in Equation (27) includes physical viscosity and turbulent viscosity, and the turbulent viscosity is often larger than the physical viscosity. As mentioned in Section 2, in the two-dimensional simulation, no matter what numerical method is used, the three-dimensional effect correction model is needed. However, the turbulence model is still necessary, so here is a brief introduction. The large eddy simulation (LES) in the Lagrangian frame [13] is introduced to calculate the turbulent viscosity in the flow field. The method proposed by Mansfield et al. [49] in Lagrangian large eddy simulation can establish the relationship between filter scale and vortex particle scale. The smoothing function *K* is related to the particle filter *G*; the relationship is shown in Equation (28):

$$\mathbf{G}\_{\Delta}(r) = \mathbf{K}\_{\Delta/c}(r) \tag{28}$$

where *c* is a constant, the filtering scale ∆ and the radios of the vortex particle *σ* satisfy the relationship: ∆ = *cσ*. The particle filter and circular filter meet the energy conservation equation in the wavenumber domain. By numerical calculation, the solution *c* = ∆/*σ* = 2.83 can be obtained. More numerical processes can be found in the work of Liu [13].

When the filtering scale ∆ is obtained, the turbulent viscosity can be calculated from the following Equation (29) according to the LES Smagorinsky model [50].

$$
\langle \nu\_{\rm sgs} \rangle = (c\Delta)^2 \langle \overline{\omega} \cdot \overline{\omega} \rangle^{1/2} \tag{29}
$$

#### 3.2.6. The Redistribution of the Vortex Particles

In order to ensure the convergence of the vortex particle method, a redistribution algorithm is in need [51]. Especially when the distance between particles is too large, the redistribution of particles can ensure a certain overlap rate between particles. There are some interpolation stencils, which distribute the vortex particle on the surrounding nodes with different weights. In this paper, the third-order interpolation function (Equation (30)) is used

$$\Lambda\_2(\lambda) = \begin{cases} 1 - \lambda^2 & , 0 \le \lambda < 0.5\\ \left(\lambda^2 - 3\lambda + 2\right) / 2 & , 0.5 \le \lambda < 1.5\\ 0 & , \text{others} \end{cases} \tag{30}$$

$$
\Lambda^{2D}\left(\frac{\chi}{\bar{h}}, \frac{y}{\bar{h}}\right) = \Lambda\left(\frac{\chi}{\bar{h}}\right)\Lambda\left(\frac{y}{\bar{h}}\right) \tag{31}
$$

where *λ* = *x*/*h*, *x* is the distance from the center of mesh to the old particles, *h* is the mesh size. Each old vortex particle will be replaced by nine new particles on the surrounding mesh nodes, as shown in Figure 5.


**Figure 5.** The schematic diagram of the redistribution of one vortex particle.

#### *3.3. Discrete Parameters and Relevant Convergence Criteria*

In this paper, the PIV experiment [10] is used as the benchmark test case for validation. The experimental parameters are shown in the Table 2 below.

**Table 2.** The basic information of the experiment.


The relevant experience of the discrete parameters and the corresponding convergence criteria will be introduced one by one:

The finite vortex method is not a time-consuming algorithm, and in order to ensure convergence accuracy, in this paper, there are 200 discrete elements on the blade surface. The time step of the finite vortex method is the time for 1 ◦ rotation, and it is a common discrete scale, which is proved by another paper [52]. The numerical parameters in the finite vortex method are conservative, and the accuracy of the program also has been verified by the previously published paper [18].

The diameter of the vortex core in the vortex particle method should be slightly larger than the chord length of VAT. In this case, the diameter of the vortex core is 1.2 times the chord length. This experience is also similar to the selection of the smooth length in the actuator line model [53].

It is worth mentioning that the vortex particle method is a mesh-free method. Therefore, there is no need to verify the independence of mesh size, and it is diffusion free in theory [54]. As long as the number of the vortex particles is enough, the wake evolution could be accurately simulated. The redistribution program in the vortex particle method redistributes each vortex particle to nine new particles through the interpolation function, which satisfies the conservation of the circulation from the zero-order moment to the second-order moment. In this way, the number of the vortex particles will be enough at the place where the vorticity is high, and the vortex particle methods will be automatically adaptive and will concentrate on the regions of physical interest [55]. Since the number of the vortex particles is adjusted adaptively by the redistribution program, the initial number of the vortex particles has little impact on the simulation of the wake-field. Therefore, it is not necessary to verify the independence of the number of vortex particles. In general, the numerical viscosity of the grid-based method will produce a big diffusive error in high Reynolds number. In contrast, the redistribution program does not bring discernible numerical dissipation. Redistribute the vortex particles in 2D simulations, which shows substantially lower errors than infrequent redistribution [46]. The redistribution program could provide a potential solution for the accurate simulation of the vorticity [56]. It can be seen that the redistribution is of great significance to the convergence of the vortex particle method.

The mesh size of the redistribution program is *h* = 0.0025 m. In order to ensure that there is a certain overlap rate between adjacent particles, the convergence criterion required that the vortex particle radius should be larger than the mesh size [46]. Hence, the radius of the vortex particle is *r*<sup>0</sup> = 0.0033 m.

According to the equation πD/c ≈17.4, the circumference is 17.4 times of the chord length. Therefore, it is set to discharge the vortex 18 times in one rotation period. That is 360◦/18 = 20◦ ; the time interval of the vortex shedding is the time for 20◦ rotation. The corresponding time step here is 0.0028 s. According to the stability criterion [46], the time step should be less than the reciprocal of the maximum vorticity. In this view, the maximum allowable vorticity under the convergence criterion should be 1/0.0028 s <sup>≈</sup>350 s−<sup>1</sup> . Although the value of vorticity changes with time due to the redistribution, the magnitude of vorticity is within ten in this case, which obviously satisfies the stability criterion.

It can be seen from the benchmark experiment that the blocking effect can be ignored. Therefore, the wake evolution in this paper is simulated under unbounded conditions, and the blocking effect has not been involved. In the vortex particle method, there is no need to set special boundary conditions for unbounded flow. This is a feature different from the grid method.

#### *3.4. The Three-Dimensional Effect Correction Model of Wake*

This new correction model is a projection function along the streamwise direction to correct the three-dimensional effect of the velocity field. As we all know, the threedimensional effect of the wake is that the third dimension will accelerate the recovery of the velocity deficit in the two-dimensional wake. By multiplying the two-dimensional velocity deficit field with a projection function, the influence of the third dimension on the velocity recovery can be equivalent.

The correction process is based on some empirical fitting functions of the VAWT wake. According to the analysis of PIV experiment results [57], the velocity deficit of the wake can

be fitted into the power-law function of the downstream position *X*: ∆*U U*<sup>0</sup> = *c*<sup>1</sup> *X Xt c*<sup>2</sup> + *c*3, where *X<sup>t</sup>* is the downstream transition location which is related to the dynamic solidity *σD*, and it can be obtained by this equation: *Xtransition* = *D*(4.78 − 4.93*σD*), where *D* is the diameter of the VAWT. These three coefficients, *c*<sup>1</sup> *c*<sup>2</sup> *c*3, are the empirical fitting coefficients which are influenced by solidity and tip speed ratio. Moreover, this form of the power-law function is also applicable to the HAWT wake and the two-dimensional or three-dimensional bluff body wake [57].

Therefore, the process of the three-dimensional effect correction proposed in this paper is based on the three-dimensional velocity distribution functions summarized by the previous study [57]. Then, we divided these three-dimensional velocity distribution functions by a two-dimensional velocity distribution function, and the proportion obtained by the division is the projection coefficient from the two-dimensional to the three-dimensional, which is exactly what we want to correct the three-dimensional effect of the wake-field. Finally, the two-dimensional velocity field is multiplied by the projection coefficient along the streamwise direction, and the two-dimensional simulation results considering the three-dimensional effect correction can be used in the preliminary array design.

In order to get the appropriate three-dimensional velocity distribution function which is close to the benchmark case in this paper, two models are selected from the previous study [57] because the parameters of these two models are relatively close to that of the benchmark case, and they were marked as VAWT\_1 and VAWT\_2, respectively. The velocity distribution of the three-dimensional wake models (VAWT\_1 and VAWT\_2) and the general form of the velocity distribution of the two-dimensional wake model are described in Table 3 below.


**Table 3.** The parameters of these three vertical axis wind turbine (VAWT) wake models.

It can be seen from these functions that the power-law exponent of the two-dimensional wake is −0.5, which is significantly larger than that of the other two three-dimensional wake models (−0.684 and −0.622). That is to say, due to the three-dimensional effect, the decay rate of the two-dimensional wake model is slower than that of the three-dimensional model.

In addition, it should be noted that the two models (VAWT\_1 and VAWT\_2) selected here are relatively close to the model in this paper; this can be seen in Figure 25 in the reference [57] that there is little difference in the power function between different VAT conditions. Therefore, there is no need to be too concerned about whether the selection of the experimental model is completely consistent with the validation model used in this paper.

In Figure 6, the curves formed by the square symbol () represent the power functions of the three-dimensional wake models (VAWT\_1 and VAWT\_2), the curve formed by the closed circle symbol (•) represents the general form of the power function of the twodimensional wake model. The functions of these three curves are described above. Then, these two power functions of the three-dimensional wake models (symbol ) are divided by the power function of the two-dimensional wake model (symbol •), and the ratio functions obtained (symbol ×) are the correction coefficient along the streamwise, which could be used to correct the velocity deficit from the two-dimensional to the three-dimensional.

deficit Δ ● ■ **Figure 6.** The correction coefficient along the streamwise direction. The X-axis is the downstream position normalized by the rotor diameter (*X/D*). The Y-axis represents two normalized parameters; one is the normalized velocity deficit ∆*U*/*U*<sup>0</sup> (symbol • and ) of the wake model, the other is the correction coefficient (symbol × and ) obtained by the division.

♦

♦

♦

♦

There are two other considerations: The first is to eliminate the unreasonable situation that the ratio is much greater than one within one diameter. Second, considering that the three-dimensional effect will make the power coefficient of VAT smaller than that of the two-dimensional model [33], and there is a positive correlation between the power coefficient and the velocity deficit, that is to say, the velocity deficit and power coefficient of the three-dimensional model are smaller than those of the two-dimensional model.

Hence, the correction coefficient should be slightly less than 1 in the near field. Based on these two considerations, combine the two ratio functions (symbol ×), what we get is the empirical coefficient (symbol) related to the streamwise position, and the combined function (symbol) can be approximately regarded as the final three-dimensional effect projection function which is used in this paper.

It can be found that the curve is a function of the distance along the streamwise direction; by multiplying the two-dimensional velocity field with the projection coefficient (symbol) along the streamwise direction, the wake-field corrected by the threedimensional effect can be obtained. Moreover, this is the whole process of the threedimensional effect correction for two-dimensional wake simulation in this paper.

In this paper, the power-law function is proposed to correct the three-dimensional effect, which will be a universal idea. Not only the wake of HAWT and VAWT but also the wake of the two-dimensional or three-dimensional bluff body conform to the form of the power-law function [57,58]. Moreover, the two-dimensional Gaussian power function was also used as the basic wake model for array layout design [35,59]. Even if the influence of the boundary layer on the wake is considered, the velocity distribution function could also be fitted into the form of the power-law function [60].

In fact, the correction model is widely used in real engineering practice to simplify the simulation, and it can be equivalent to the effect of some complicated turbulence details, such as dynamic stall and complex turbulent inflow conditions. Especially for the array design, the first-order approximation of the velocity distribution is enough to meet the demand, and the three-dimensional effect correction of the two-dimensional wake simulation is a good choice. The correction model proposed in this paper is derived from the power-law function of velocity distribution, which could make the numerical results in good agreement with the three-dimensional results. The validation of this numerical method is shown in Section 4.

#### **4. Validation and Discussion**

In order to verify the reliability of the new numerical method proposed in this paper, a PIV experiment is selected to compare with the simulation results. The parameters of this PIV experiment have been described in detail in Table 2 in Section 3.3.

The comparison results are shown in Figure 7. Relatively speaking, the average velocity profiles calculated by the proposed method match well with the PIV experimental data [10]. At least for the preliminary array design, it is an encouraging result to achieve at least first-order accuracy in the two-dimensional wake simulation. It also proves that the two-dimensional numerical method with a reasonable three-dimensional correction is an acceptable and practicable approach that can make the wake simulation faster and the array design easier. In Figure 8, the rotation direction of the VAWT is counterclockwise. According to the analysis of velocity distribution of a VAWT wake, the peak value of velocity deficit is found not on the centerline. This asymmetry can be seen in the velocity field (Figure 8).

It needs to be clarified here, in the preliminary array design, only the time-averaged velocity field needs to be considered, and the unsteady flow details can be ignored. In addition, since the validation data of the reference is time-averaged, and the benchmark test case did not provide the transient experimental data at different times. In order to keep consistent with the benchmark test case, the time-averaged velocity field was used in this paper.

In the vorticity field (Figure 9), due to the velocity deficit caused by the turbine and the conservation of momentum, the wake vortex will expand in the near field, and the expansion phenomenon of the stream tube in the numerical simulation is similar to the flow field observed by another PIV experiment [5]. It can also be seen that there is a skewness in the distribution of the wake vortex. According to the Magnus effect, a lateral force perpendicular to the inflow direction will be produced on a spinning object. In addition, according to Newton's third law, the forces acting on the turbine and the forces acting on the

flow field belong to action and reaction. In this case, the VAWT rotates counterclockwise and is subjected to a downward lateral force. Therefore, the wake is subjected to the upward force and will move upward. The similar phenomena of the wake skewness have also been found in some other VAT wake experiments [8,11]. It can be seen that there is a skewness in the VAT wake caused by the Magnus effect. Moreover, the faster the rotation, the greater the wake skewness. The explanation that wake skewness is related to the Magnus effect can also be found in other studies [5,9,61]. However, from the micro point of view, the generation of the wake skewness is caused by the wind-blade interaction, and it can also be explained from the macro point of view that the wake skewness is caused by the Magnus effect. Both explanations are correct. By comparison, there is usually no skewness in the HAWT wake, this consideration will make the array layout of VAT different from that of HAT. Moreover, more suggestions can be provided for the array arrangement of VATs, such as, the arrangement distance for the downstream VAT should be different considering the rotation direction of the upstream VAT.

Δ **Figure 7.** The comparison of average velocity profiles at six different cross-sections. The X-axis represents the sum of the two normalized parameters; one is the streamwise position, which is normalized by the turbine diameter (*X/D*), the other is the velocity deficit, which is normalized by the inflow velocity (∆*U*/*U*<sup>0</sup> ). The Y-axis is the lateral position, which is normalized by the turbine diameter (*Y/D*). The data of black curves are from a published paper [10]; the color curves are calculated by the numerical method proposed in this paper. Δ

Δ Δ **Figure 8.** The time-average velocity field of the VAWT wake. The X-axis is the streamwise position, which is normalized by the turbine diameter (*X/D*). The Y-axis is the lateral position, which is normalized by the turbine diameter (*Y/D*). The color bar is the velocity deficit, which is normalized by the inflow velocity (∆*U*/*U*<sup>0</sup> ). The blue dashed lines indicate that the velocity profiles at that location are used for comparison with the PIV experimental results.

'

**Figure 9.** The vorticity field at the near wake. The X-axis is the streamwise position which is normalized by the turbine diameter (*X/D*). The Y-axis is the lateral position which is normalized by the turbine diameter (*Y/D*). The color bar is the nondimensional vorticity which is normalized by the maximum vorticity in the field (*ω/ωmax*).

*ω/ω* There is another interesting phenomenon in the far wake in Figure 10. It can be seen from the following eight moments in one period (the period here is not the rotation period, but the period of the vortex street oscillation), due to the complex evolution, the vortex pairs in the far-field are similar to the Karman vortex street in the flow around a cylinder. This characteristic can also be found in another reference [62]. The similar wake evolution will occur in the far-field of a HAT wake as well [39]. It is still unknown whether this complex evolution is similar to the meandering of the HAWT wake. However, at least, this interesting phenomenon in these eight images shown here can lead to a new research direction in the future.

It is also worth mentioning that the vorticity of the vortex particle swarm is different at each time step, and the redistribution program makes the number and vorticity of the vortex particles change all the time. Therefore, the peak value of the vorticity field (*ωmax*) and the range of the color bar (*ω*) also change at each time step. In order to unify the color bar for observation, the vorticity field (Figures 9 and 10) is normalized by the maximum vorticity at that moment. The color bars of the wake-field are represented by the dimensionless vorticity (*ω/ωmax*). In general, although the redistribution program makes the range of vorticity field change all the time, for the vortex particle method, the whole field still satisfies the conservation condition, so it may not be necessary to pay too much attention to the absolute value of the vorticity at each moment.

However, the purpose of this paper is to propose a fast simulation method for the VAWT wake, so the complex wake evolution phenomenon will not be discussed in depth. The simulation results are calculated by the code written by the authors, and there is no commercial software to use this numerical method so far. In terms of the calculation efficiency of the algorithm, the vortex particle method takes less than 2 h on 4-core CPUs, and the number of the vortex particles is about a hundred thousand, when the focus is only on the flow field within 10D, as we all know, the calculation cost here is much less than that of the grid-based method.

**Figure 10.** *Cont.*

*ω/ω* **Figure 10.** The eight moments of vorticity field in one period at a wide view. The interval of each moment is approximately equal to one rotation period. The X-axis is the streamwise position, which is normalized by the turbine diameter (*X/D*). The Y-axis is the lateral position, which is normalized by the turbine diameter (*Y/D*). The color bar is the nondimensional vorticity, which is normalized by the maximum vorticity in the field (*ω/ωmax*).

#### **5. Conclusions**

*ω ω ω ω* In this paper, by combining the finite vortex method and the vortex particle method, a fast and relatively accurate two-dimensional numerical method is proposed to simulate the wake-field of the vertical axis wind turbine. Based on the power-law function of velocity distribution, a new correction model is proposed to correct the three-dimensional effect of the VAWT wake. Then, the feasibility of the numerical method is preliminarily verified by the PIV experiment. Finally, according to the flow field analysis, the skewness of the VAWT wake is discussed. In the far-field, the VAWT wake evolves into the Karman vortex street, which is similar to the flow around a cylinder.

Although the published experiments on the VAT wake, which can be used for the verification, are relatively limited, however, the numerical method proposed in this paper is feasible in theory and has been preliminarily verified. It has been proven that this new method can simulate the velocity distribution quickly and relatively accurately. The first-order approximation can be achieved by this two-dimensional numerical method, which meets the requirement for efficiency in the preliminary array design. In the future, with more experimental data published, the numerical method proposed in this paper will be systematically studied and improved. Nevertheless, we believe that it can be a promising approach to make the VAWT array design easier and more reasonable.

**Author Contributions:** Z.Y.: Investigation, methodology, writing the original draft preparation; J.J.: methodology, supervision; J.Z.: methodology, supervision; Q.S.: methodology; K.S.: methodology, validation; X.Z.: methodology; R.J.: Writing—review and editing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (No. 51879064, No. U1706227, No. 51979062).

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

#### **References**


## *Article* **Aerodynamic Study of the Wake Effects on a Formula 1 Car**

#### **Alex Guerrero and Robert Castilla \***

LABSON—Department of Fluid Mechanics, Universitat Politècnica de Catalunya, ES-08222 Terrassa, Catalunya, Spain; alex.guerrero.lorente@estudiantat.upc.edu

**\*** Correspondence: robert.castilla@upc.edu

Received: 7 September 2020; Accepted: 23 September 2020; Published: 5 October 2020

**Abstract:** The high complexity of current Formula One aerodynamics has raised the question of whether an urgent modification in the existing aerodynamic package is required. The present study is based on the evaluation and quantification of the aerodynamic performance on a 2017 spec. adapted Formula 1 car (the latest major aerodynamic update) by means of Computational Fluid Dynamics (CFD) analysis in order to argue whether the 2022 changes in the regulations are justified in terms of aerodynamic necessities. Both free stream and flow disturbance (wake effects) conditions are evaluated in order to study and quantify the effects that the wake may cause on the latter case. The problem is solved by performing different CFD simulations using the OpenFoam solver. The significance and originality of the research may dictate the guidelines towards an overall improvement of the category and it may set a precedent on how to model racing car aerodynamics. The studied behaviour suggests that modern F1 cars are designed and well optimised to run under free stream flows, but they experience drastic aerodynamic losses (ranging from −23% to 62% in downforce coefficients) when running under wake flows. Although the overall aerodynamic loads are reduced, there is a fuel efficiency improvement as the power that is required to overcome the drag is smaller. The modern performance of Ground Effect by means of vortices management represent a very unique and complex way of modelling modern aerodynamics, but at the same time notably compromises the performance of the cars when an overtaking maneuver is intended.

**Keywords:** Formula 1; Computational Fluid Dynamics (CFD); external aerodynamics; OpenFoam; snappyHexMesh; incompressible flow; Federation Internationale de l'Automobile (FIA); downforce; drag; vortex; wake

#### **1. Introduction**

For many years, it has been openly stated that F1 has lost much of its spectacular nature due to the difficulty of the cars in being able to follow each other closely for a long period of time. The sophisticated aerodynamics of these single-seater cars has compromised the chasing of the leading car (mostly due to the turbulent wake generation and clearly disturbed flow [1]). This way, it was found to be convenient to numerically analyse and quantify the actual loss of aerodynamic loads on a F1 car due to the immediate presence of a rival in front of it and then, establish a fair comparison with the situation of a free stream condition.

In recent years, a significant number of investigations has been performed, both conducted by the Federation Internationale de l'Automobile (FIA) or Formula One teams and by other sources of investigation. The works of Ravelli and Savini [2,3] have been taken as a reference, as they are considered to be a feasible approach to such a complex CFD problem as the very same CAD (Computer Aided Design) model is evaluated. The obtained results show interesting vorticity behaviours and the

reference data for the comparison of their numerical results under free stream flows are taken into account in this manuscript.

Besides, the works of Newbown, et al. [4] and Perry, et al. [5] have set a positive precedent on how to numerically evaluate a F1 car under wake flows, despite using a quite outdated CAD geometry. The methodology that they applied, consisting of varying the different distances between cars, is taken as a good inspiration on how to operate with the present model. However, it is believed that studying closer distances between the cars under the wake flows could potentially give a more complete definition of the study. That is why the closest distance that is evaluated between the two cars on the present manuscript is around a quarter of a car length, which is significantly smaller that the ones that are reflected on other publications.

Other authors, such as Ogawa et al. [6] or Senger et al. [7], have also contributed to the definition of the whole methodology applied this CFD study. Special attention is also placed on the research carried out by Larsson [8] that was backed by the BMW Sauber F1 team, where the description given of the wing interaction with the tire wake is crucial for the understanding of the on-set flow to the underbody.

The main contribution and goal of this paper is to evaluate, study, and numerically quantify the aerodynamic performance of a 2017 spec. adapted F1 car (see Figure 1) under free stream conditions and wake flows with the purpose of being able to argue whether the 2022 changes in the regulations are somehow justified in terms of aerodynamic necessities. In addition to, an original approach is given when including some energetic calculations, so as to see how the fact of running under wake flows can harm or benefit the overall power requirements and, therefore, affect the fuel efficiency of the cars.

This work can be considered to be contemporary study, as very few sources have extracted meaningful conclusions regarding the latest F1 aerodynamic package. Accordingly, the originality of this work recalls on the judgement that the authors may offer on the present F1 regulations, so as to contribute to a significant upcoming improvement of the category.

**Figure 1.** Ferrari SF70H, a 2017 spec. car. Extracted from Reddit [9].

The CFD methodology is primarily selected as other resources, such as the use of a wind tunnel or any other experimental solutions, are currently out of reach to deal with such a study (see Newbon et al. works in [10]). The amount of energy resources that a wind tunnel experimental execution would require (without mentioning the environmental impact of its construction) is also a fundamental point to opt for a more sensible approach. Moreover, as the CFD discipline involves a rather strict and accurate process to be able to deal with external aerodynamic problems, the methodology is

accepted in order to discern and evaluate the hypothesis established. It has to be commented, though, that CFD models are usually calibrated (see [11,12]), but for the present research, it was not possible. Accordingly, a future investigation regarding the usage of experimental data acquisition would validate this approach.

In the following Section 2, the Materials and Methods used along the research are described. In Section 3, the obtained results are presented with a proper discussion and, finally, Section 4 sums up the most relevant and significant conclusions.

#### **2. Materials and Methods**

Computational Fluid Dynamics is the branch of Fluid Mechanics that uses numerical methods and algorithms so as to solve and analyse the behaviour of the fluids. In order to benefit from CFD techniques, it is quite important to know as much as possible regarding the real problem that is intended to be simulated (this are the physical properties of the fluid, the boundary conditions as well as other variables). Moreover, CFD involves the design of a CAD model to be meshed and eventually solve a compendium of mathematical equations on it, as they are the fundamental tool to solve in numerical problems. However, some drawbacks that are inherent to CFD involve the model calibration (to ensure reliable results), complex meshing and usually powerful CPU requirements, among others.

Two different scenarios were initially contemplated within the evaluation of the present work. A free stream approach, where the behaviour of the F1 car is analysed under optimal conditions and a more challenging one, based on a conditioned situation under wake flows. As for the latter, different distances were evaluated in order to appreciate the changes and behaviour of the car in second place in direct comparison with the leading car.

#### *2.1. Geometry*

The CAD geometry employed was the PERRIN F1 car [13], as it was a realistic approach to a real F1 car because it was designed under the 2017 FIA regulations. This includes several small realistic details, such as winglets, vanes, vorticity generators and slots, being the smallest elements 1.5 mm thick. The length *L* of the car, which is later considered for the dimensions of the overall domain is around 5.350 m, while its wheelbase, required for the turbulent length scale measures 3.745 m. The geometry was exported to Parasolid format (.*x* \_ *t*.), so as to be processed by the SnappyHexMesh command. Figure 2 displays a generic view of the aforementioned F1 car.

**Figure 2.** Generic view of the F1 car CAD model used.

#### *2.2. Solver*

The simulations were executed using the OpenFoam [14] toolbox that solves the equations for unsteady incompressible flow of mass and momentum presented in Equation (1), (also known as Navier–Stokes equations [15]); as well as the continuity equation shown in (3).

$$
\rho \frac{\partial \mathbf{u}}{\partial t} + \rho (\mathbf{u} \cdot \nabla) \mathbf{u} = \rho \mathbf{g} - \nabla p + \mu \triangle \mathbf{u} \tag{1}
$$

where

$$
\triangle \equiv \nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} \tag{2}
$$

$$
\nabla \mathbf{u} = \mathbf{0} \tag{3}
$$

The discretisation of the domain was obtained applying the FVM (Finite Volume Method), as this method is locally conservative due to being based on a balance approach [16]. The gradient, divergence, and laplacian terms of the Navier–Stokes equations were discretised by means of the Gaussian schemes: at cell interfaces, the interpolation schemes were linear (upwind).

On the other hand, the simpleFoam algorithm (Semi-Implicit Method for Pressure-Linked Equations) was the one chosen, as it is appropriate for incompressible, steady, turbulent flows. The reason why a steady solver was applied on a flow field with unsteady characteristics resides in the nature of turbulence, which is, by definition, non-stationary. Accordingly, the presented results (wake vortices among others) are averaged on a medium flow, which is what simpleFoam solves. The goal was not to study transient phenomena (which could be the subject for future works), but these averages to evaluate the effect on average stability of the car.

The GAMG (Geometric Algebraic Multi Grid) solver was utilised for the pressure equation, while smoothSolver was the one that was selected for velocity and turbulence variables.

As for the execution, the initial 800 iterations were performed under a first order discretisation, while the remaining ones were performed using second order schemes.

For the turbulence model, the k-*ω* SST (k-*ǫ* is used in the outer region of and outside of the boundary layer and k-*ω* is used in the inner boundary layer [17]) was selected, since it offered a faster convergence and better results as compared to k-*ǫ* and Spart Allmaras.

#### *2.3. Domain and Mesh*

The fluid domain dimensions were inspired by other publications, such as Broniszewski's [18], who puts special attention to the back region of the domain in order to be able to capture the generation of the wake properly. Figure 3 shows the boundaries of the problem.

**Figure 3.** Fluid domain dimensions.

A three dimensional hexahedral mesh by means of the snappyHexMesh and blockMesh commands was generated. The height of the first cell at the solid surfaces is set at 0.01 mm with a layer expansion ratio of 1.3. The resulting average value of y+ is around 40, which involves the use of wall functions as opposed to a near wall treatment, which generally adopts y+ around 1 when solving for low Reynolds (Re) number [19].

Many refinement enclosures were specified along the geometry in order to be able to capture the effects of the wake and and other phenomena. Special attention was placed in the massive wake region as well as on the downforce generation zones, such as the wings and floor. The results of the mesh procedures may be checked in Figures 4–6.

**Figure 4.** Overall mesh and refinement enclosures of the free stream case.

**Figure 5.** Overall mesh and refinement enclosures under wake effects.

**Figure 6.** Mesh (**A**) front, (**B**) rear wing, (**C**) trimetric, and (**D**) cockpit.

A Grid Convergence Index (GCI) has been calculated in order to guarantee that the analysis of the results that were obtained was independent of the grid size. The GCI methodology is performed according to Celik et al. [20] and the study is performed by varying the level of definition of the three different meshes studied. The refinement levels used range from 5–6 in the coarsest, 6–7 in the intermediate, and 7–8 in the finest. The *φ* value analysed corresponds to the lift coefficient multiplied by the surface (*S* · *CL*). Table 1 presents the results obtained of the GCI.


**Table 1.** Grid Convergence Index result of the meshes studied.

The GCI values are in the asymptotic range of convergence, both *GC If ine* and *GC Icoarse*, as shown in Table 1. It is as well possible to note that the difference in *φ* values between the finest and intermediate meshes is considerably small thus, obtaining a rather small percentage of error. As for the computational resources, the finest mesh (which is formed by 19.2 million cells in the free stream case and around 40 million under wake flows) took an approximated time of 8 and 14 hours respectively to be generated. The grid size is found to be reasonable as it is almost identical to the one with which the results are later compared to (see Section 3). Hence, as computational resources do not suppose a huge compromise nor restriction, it was preferred to choose the finest mesh for the following study. The sacrifice in computational time and consumption was notably overcome by the superior amount of detail and definition that the finest mesh could deliver.

#### *2.4. Boundary Conditions*

The general boundary conditions established are as follows:


According to [21], it is possible to assign the car wheelbase, which represents the size of the largest eddy, as turbulent length scale *l*. All of these data are shown in Table 2.


#### **Table 2.** Setup of the simulation.

#### *2.5. Simulation Performance*

The simulations were executed in a cluster that was equipped with 32 cores and 64 Gb of RAM memory and the post-processing operations were carried out by means of an i7 six core laptop with 16 Gb of RAM.

Finally, it was considered to reach convergence when the pressure and velocity residuals were lower than 10 −3 , just as seen in Figure 7. The run time of the cases was around 10 hours under free stream conditions and 17 h under wake flows.

**Figure 7.** Residuals of the simulation

The obtained results are compared with the reference data proposed by Perrin [23], which are obtained through a CFD simulation using TotalSim, a professional enterprise highly experienced in CFD problems. This is taken as a reference as the boundary conditions and the problem definition coincide with the ones that are proposed in this study, although the solvers are different. Besides, reference data have a public background of more than 120 runs evaluated, so this may be considered to be a reliable source of validation.

The numerical prediction deals with the aerodynamic parameters: Downforce (S*CL*), Drag (S*CD*), overall aerodynamic efficiency *E* (*CL*/*CD*), and Front Balance *FB*, which is defined as the fraction between the downforce that is generated by the front axle and the total one.

#### **3. Results and Discussion**

#### *3.1. Free Stream Condition*

Table 3 shows that the results of the performed RANS simulation with k-*ω*SST model encounter a notable agreement with the reference data. The error found in both downforce and drag coefficient differs by 4.45% and 6.50%, respectively. This indicates that the overall simulation is acceptable, as the difference between the aerodynamic coefficients is found to be small enough to accept them and validate the reference datum.


**Table 3.** Comparison between the simulation and reference data.

On the other hand, it is found appropriate to indicate the relative contribution in terms of downforce and drag of the main components of the car. This way, it is possible to see the influence of each element and appreciate its aerodynamic efficiency, which is very useful in terms of redesign purposes and issues detection. The different components evaluated can be checked in Figure 8, while Table 4 shows the distribution of these aerodynamic forces.

**Figure 8.** Evaluated components of a F1.

**Table 4.** Relative downforce and drag contributions of the different elements.


The underbody, which is composed by the flat floor, the plank, and the diffuser, is responsible for the 60% of the total downforce generation. Following this trend, the rear and the front wing represent, respectively, around 35% and 23% of the overall downforce of the car. The bodywork, shaped as a wing profile, counterbalances these gains by producing lift as well as other elements, such as the front suspension, and both rear and front tires.

On the other hand, Table 4 shows that tires represent around 30% of the total drag, specially on rear's, as the wheels are not covered. This area is well influenced and dominated by big pressure losses in the wake, which leads to these undesired highly turbulent zones of the total drag of the car. Additionally, the front and the rear wing are also present, by being responsible for 13% and 20% of the total drag, respectively. Other elements, such as the underbody (15% of the drag) and the bodywork (10.70%), take an important contribution regarding the aerodynamic resistance.

In terms of aerodynamic efficiency, it is important to note that the underbody is, with no hesitation, the most efficient part of the car. This can be explained due to the use of Ground effect (despite being limited by the flat floor and the size of the diffuser). As opposed to that, the rear wing is known for its low aspect ratio, which helps to generate big downforce quantities, but it suffers from the production of induced drag. The efficiency of the front wing is somewhere between the underbody and the rear wing: the beneficial points of the rake and ground effect that are commented in [24] are counterbalanced by the high angle of attack of the flaps and the vortices generated at the tips.

Figure 9 shows the pressure distribution by means of the dimensionless coefficient (*CP*) around the car, both in terms of a plot and in a three dimensional visualisation. The upper view reflects high-pressure zones that are located in the nose (specially in the front wing) as well the rear wing due to its high curvature. Some stagnation areas are also found around the cockpit, where the pressure distribution is low and smooth. These results are in line with what it was previously commented, with the wings being one of the main generators of downforce due to its high angle of attack and curvature.

On the other hand, the underbody of the car shows low-pressure zones under the wings (as it was clearly expected by its nature and shape). Besides, the low-pressure zones along the floor and diffuser suggest that the car is working properly under the Ground effect. It is possible to see a smooth transition from a medium pressure zone to a low pressure region (meaning that the airflow is being accelerated), and finally an increase of pressure that returns the airflow in a lower velocity to the wake. However, the region in close proximity to the tires is affected by an increase of pressure, which implies the loss of the Ground effect benefits.

In general, the pressure distribution is rather smooth around the car, with no abrupt pressure gradients or unexpected transitions, just like that observed in [4].

**Figure 9.** Pressure coefficient plot and representation on the upper and under surfaces of the car.

Figure 10 pretends to illustrate to the reader the shape, position and prolongation of some three-dimensional vortical structures, as well as provide several details regarding the strength and rotational axis of such vortices.

31

**Figure 10.** Several axial vorticity contours: (**A**) YZ plane at 0.65 m in the X direction, (**B**) YZ plane at 3.65 m in the X direction, (**C**) YZ plane at 2.60 m in the X direction, and (**D**) XY plane at 0 m in the Z direction.

Figure 10A shows the multifaceted front wing: not only generating downforce, but also axial vorticity that tends to avoid the front tires and energises the flow downstream. Several vortical structures can be seen on the tip of the endplate (A1) or the winglet endplate (A2), among others. The rotation of the vortices—clockwise or counter clockwise—depends on the pressure field around them [25]. A quite interesting phenomena, as it is clearly the biggest vortex generated in the front wing, is the so-called Y250 vortex. This vortex is developed between the middle section of the wing and the multi flap surface, and it is aimed at recirculating the flow towards the underbody of the car (*inwash*).

Figure 10B, shows the rear part of the car, where several vortices are originated as a result of the wake of the spinning wheels and other devices. Special attention is placed in the Venturi vortices that are generated on the side of the diffuser due to the pressure gradient between the underbody and the outside. Additionally, the strakes of the diffuser generate small vortices that are coupled with an opposite rotating vortex due to the interaction with the ground boundary layer, similar to the ones observed in [26].

Figure 10C reflects the generation of vortices in the upper middle region of the car. The bargeboards and vanes play a special and important role here. The goal of the latter is to seal and canalise the flow over the bodywork, making sure that the flow keeps attached along the car (Coanda effect). Besides, the sealing vortices that are generated by the bargeboards tend to act as skirts, therefore preventing the underbody airflow to escape and maximise the Ground effect [27] and the diffuser efficiency.

Figure 10D displays a top view of the whole car to understand the behaviour of the overall vorticity.

Finally, Figure 11 shows a three dimensional representation of the streamlines of the vorticity. It is possible to see, in general traits, how the flow behaves around different areas of the car (Y250 vortex, tip of the wings, middle section) and how chaotic and turbulent the resulting wake looks like.

**Figure 11.** Streamlines of the axial vorticity.

#### *3.2. Under Wake flows*

Similarly, the evaluated results deal with the aerodynamic parameters: downforce (S*CL*), drag (S*CD*), overall aerodynamic efficiency *E* (*CL*/*CD*), and front balance *FB*.

However, four different situations are studied, which differ in the distance established between the 2 cars: 0.25 L, 0.5 L, 1 L, and 2 L, where *L* represents the length of one car (approximately 5.3 m). Figures 12–15 display the evaluated parameters of the second car (follower car) and Table 5 shows the percentage of change of those with respect to the first car (leading car).

**Figure 12.** Comparison between values of S*C<sup>L</sup>* under wake flows and under free stream conditions.

**Figure 13.** Comparison between values of S*C<sup>D</sup>* under wake flows and under free stream conditions.

**Figure 14.** Comparison between values of efficiency under wake flows and under free stream conditions.

**Figure 15.** Comparison between values of front balance under wake flows and under free stream conditions.

**Table 5.** Percentage of change of the aerodynamic coefficients of the second car as regards the first car.


The obtained results show that the reduction in the aerodynamic coefficients is clearly visible from an initial distance of two car lengths (approximately 10.6 m) to the closest case studied of 0.25 L (less than 1.5 m). The reduction of downforce ranges from a −23.5% to a very significant −62% in the worst case scenario. In a similar progression, the drag is reduced from a −14.2% to a −40% and, for this reason, so does the overall efficiency of the second car.

Besides the distinguished loss of downforce, the second car experiences a dramatic increase of front balance (FB) from +26% to 40%. This sudden increase on the front aerodynamic loads may presumably lead to experiencing oversteer (oversteer is caused when a car steers more than intended, thus losing the rear end) and safety issues while braking and on high-speed corners [28].

In general traits, it can be seen that, as the second car approaches and gets closer to the leading one, the aerodynamic loads are reduced, which worsens the performance of the car, but so does the drag. This is energetically a key point, as the power that is required to overcome the drag is smaller when the distance between the cars gets closer, which enables less fuel consumption for the car behind. Table 6 shows a hypothetical situation on the main straight of "Circuit de Barcelona", Catalunya, where the power and the energy that are required to overcome the different situations are evaluated (it has been assumed a distance of 1 km and a car speed of 50 m/s along the whole straight).

These same conclusions can be easily extrapolated to current road cars, although the conditions and data may differ notably, but not the overall conclusions.


**Table 6.** Energetic parameters required to overcome the air resistance.

Figure 16 is aimed at showing the obtained data in a visual representation, so that is possible to appreciate the rate of change in the various studied parameters.

**Figure 16.** Percentage change of the performance of the second car in respect of the leading car.

As mentioned, the increase of the front balance levels on the second car enables a short, but rather interesting discussion. It is known that the weight distribution of a F1 2017 specification car is around 45.5% on the front axle [29], so the car is not supposed to exceed this 45.5% of front balance, as it may lead to stability concerns. The Center of Pressure, which, by definition, is such where the total sum of pressure fields act on, should always remain behind the Center of Gravity. This can be explained, as the yawing moment of the aerodynamic forces counterbalances the steering of the driver and, therefore, stabilises the car.

On the other hand, if the Center of Pressure is ahead of the Center of Gravity, then the yawing moment increases the sideslip angle and produces instability. As the obtained results show, the front balance of the second car adopts always values that are greater than 45.5% (see Table 5). This implies not only a reduction of the stability and the performance of the car (slower laptimes and more degradation of the front tires), but also a more challenging approach when driving the car; major ease of spinning and safety issues.

Moreover, the study analyses the performance of the most relevant aerodynamic devices on the follower car: front wing, rear wing, and diffuser.

Table 7 reflects that the loss of downforce of the front wing starts to appear severely at a distance of 1 L up until a very critical −38%, when it reaches 0.25 L. On the other hand, the drag levels are found to be grater than the leading car at a distance of two car lengths, but a similar behaviour to the downforce is found as long as the distance is decreased. This could be explained due to the turbulent conditions that are found far away from the leading car. Moreover, as the second car approaches the leading one, it enters into a very strong wake region that initially hits the front wing and eventually influences all behaviour of the car.

**Table 7.** Percentage of change of the aerodynamic coefficients in the front wing of the second car as regards the first car.


This is precisely seen in Figure 17, where it is possible to see the behaviour commented earlier, as the low pressure zones of the front wing suffer an abrupt change and loss as the car enters into the closer wake region (cases 0.25 L and 0.5 L, where the tendency of the pressure distribution appears slightly inconsistent, although the numerical results do not). It is not at all elementary to try to find a reason why it occurs, but the most sensible explanation lies in the chaotic and unpredictable nature of turbulence. Besides, the larger pressure area located in 0.5 L is found in the very first tip of the front wing, which corresponds to the first element that is encountered by the airflow. In general, these changes are noticeably far less progressive than the rear wing's, which also evidences the premature front balance shift, but allows a more robust behaviour at greater distances.

As for the streamlines of the velocity, Figure 18 (2 L) shows that the magnitude of the velocity is somehow similar to the one of the free stream flow, although it presents low velocity areas on the central part of the wing. However, endplates and wing tip elements are still useful for redirecting the flow towards the rear. As the distance is halved, the velocity of the flow is reduced and the front wing loses its capacity to govern the airflow, but it is not until cases 0.5 L and 0.25 L that the front wing is fed by really low kinetic energy flow that leaves it notably inoperative. The effectiveness of the generation of vortices and redirection of the flow is insignificant, just as checked in Table 7 due to the strong wake region.

**Figure 18.** Streamlines of the velocity. Front wing of the second car at a distance of 0.25 L, 0.5 L, 1 L, and 2 L as compared to the free stream case.

Shifting to the rear wing, Table 8 portrays that the behaviour of the rear wing is notably different: the loss of downforce is very notable since the very first distance of 2 L and it keeps increasing as the slipstream distance gets closer. Similarly, the drag levels are also strongly reduced from the beginning and matching a comparable ratio that keeps the efficiency almost constant.

It can be stated then, that the rear wing is more affected under the wake effects than the front wing, as the latter seems to suffer less and only under close proximity. This somehow explains the increase of front balance (and forward change of the center of pressure) and its posterior decrease, as sketched in Figure 16.


**Table 8.** Percentage of change of the aerodynamic coefficients in the rear wing of the second car as regards the first car.

On the pressure side, from an early stage, one can see that the pressure loss in the rear wing is severely and gradually appreciated. The approach to the leading car distinctly damages the functioning of both its stabilising and suction purposes. Figure 19 sketches these comments.

**Figure 19.** Pressure coefficient distribution on the rear wing of the second car at a distance of 0.25 L, 0.5 L, 1 L, and 2 L as compared to the free stream case.

Similarly, under a normal regime, Figure 20 shows that the rear wing works perfectly with aligned flow that is redirected by other aerodynamic devices, but it is possible to see that, as the distance is reduced, the streamlines tend to deflect slightly inboard. The clear jump in terms of velocity magnitude is produced from case 1 L to 0.5 L, where the flow experiences a moderate deceleration as it enters into the strongest wake region. It could be speculated that the rear wing is more sensitive to the direction of the flow than the front wing due to its low aspect ratio, the pronounced curvature, and the massive endplates at both sides. Figure 20 (0.25 L) shows a mixture of velocity flows with very low speed that unexpectedly create an inwash phenomena at the rear end of the wing. This evidences that the whole aerodynamic package of the follower car is prominently disrupted under the wake conditions.

Finally, Table 9 reports that the diffuser is the device that suffers the most under the wake conditions, as the downforce loss that is in close proximity is around 70% and the drag is found to be reduced by around 57% at the same distance. However, the interesting conclusion is that the reduction of downforce starts from the very first beginning and keeps increasing as the distance is reduced.

This evidences that overtaking maneuvers are heavily influenced since the start (it is important not to forget that the diffuser is the greatest source of downforce generation, as seen in Table 4).

This generates an interesting and deep discussion, as it is encountered that the greatest way of generating aerodynamic loads under a free stream flow, is, at the same time, the worst one under wake flows. The whole conception and functioning of the underbody is found to be absolutely pointless; therefore, other aerodynamic paths and solutions should be evaluated if these losses are wished to be recovered.

On the pressure side, the results that are presented in Figure 21 again show that the whole low pressure zone is affected since the very first beginning, with moderate areas near the diffuser strakes with higher pressure values. However, the degradation of the performance is perfectly noticeable and again proves that the diffuser suffers excessively under wake flows until its contribution becomes almost negligible.

As for the streamlines of the velocity, Figure 22 describes how the underbody, and particularly the diffuser, is affected by the wake flow. In free stream conditions, the diffuser is fed by a high energised flow that is redirected by the front wing and guided around the flat floor. Nonetheless, the streamlines of the flow are not completely straight, as the vortices management allow for the control of the airflow around the underbody. Figure 22 (2 L) displays a quite non-disturbed behaviour of the airflow, as the management of it is still acceptable. When the distance is reduced, the performance of the underbody starts worsening due to the flow arriving more disturbed into the diffuser, hence not being able to

operate properly. Smaller distances, such as cases 0.5 L and 0.25 L, reflect an underbody region that is fed by a very low kinetic and rotational flow.

The main reason for the massive downforce losses that are reported in Table 9 is that the underfloor is notably sensitive under wake flows, as it works closer to the ground than the front wing. This means that the low energy (and highly-rotational) airflow may not be compressed around this small area, therefore experiencing a massive performance loss.

**Figure 20.** Streamlines of the velocity. Rear wing of the second car at a distance of 0.25 L, 0.5 L, 1 L, and 2 L as compared to the free stream case.



**Figure 21.** Pressure coefficient distribution on the diffuser of the second car at a distance of 0.25 L, 0.5 L, 1 L, and 2 L as compared to the free stream case.

**Figure 22.** Streamlines of the velocity. Diffuser of the second car at a distance of 0.25 L, 0.5 L, 1 L, and 2 L as compared to the free stream case.

Finally, Figure 23 shows the plots of Tables 7–9 in order to appreciate the aforementioned changes in the aerodynamic coefficients of the different elements.

**Figure 23.** Percentage of change on the performance of the second car.

On the wake side, these effects can easily be appreciated in Figure 24, where the second car is affected by a flow that is lower in terms of kinetic energy, as the wake that is generated by the leading car is released far away disturbing its follower. It is also clear that, as the second car gets closer, it inherently enters into a unique wake structure characterised by very low speed flow, that ranges from 0 to 10 m/s, therefore, resulting severely affected. It is seen that as the second car reduces the distance, its wake originates a separation region that enlarges and becomes evident as the distance is closed. At a large distance (2 L), the second car's wake adopts a needle shape, which somehow imitates the free stream natural wake, but this shape soon disappears at closer distances.

As for the leading car, it can be noted that its wake is not notably modified (in terms of shape and contours) by the presence of the follower car. The aerodynamic coefficients that are evaluated on the leading car only experience a low variation —less than 3%—when the distance between the two cars is set at 0.25 L.

**Figure 24.** Velocity contours on a top plane at a distance of (**A**) 0.25 L, (**B**) 0.5 L, (**C**) 1 L, and (**D**) 2 L.

Shifting now to the streamlines of the velocity, Figure 25 displays two different planes for each scenario. The first thing that one may notice is how the airflow is perfectly attached to the first car; from the front wing, going through the sidepods and bodywork and finally exiting the rear wing.

However, it is possible to also see that the exiting airflow on the rear-end of the leading car is somehow divided into two characteristic flows: the first one, located on the superior area, which adopts high-speed values due to being accelerated around the bodywork (low pressure zone and smooth behaviour), and the second one, which exits the diffuser upwards and is mainly a turbulent flow continuously undergoing changes in both magnitude and direction, as reported in [30]. The mixture of the previously mentioned flows is what originates such a chaotic wake region, as it is formed by the combination of multiple flows with various natures and velocities [31].

**Figure 25.** Streamlines of the velocity at a distance of (**A**) 0.25 L, (**B**) 0.5 L, (**C**) 1 L, and (**D**) 2 L.

#### **4. Conclusions**

The initial purpose of the present study was to establish whether the current F1 regulations required an urgent aerodynamic change in order to objectively improve the category. In other words, to test and validate, by means of CFD analyses, that the upcoming regulation changes are aerodynamically justified. Accordingly, the results gathered in the present manuscript may be useful to F1 teams, FIA administrations, and to anyone interested in racing aerodynamics.

The obtained results suggest that modern F1 cars are designed and well optimised to run under free stream flows (as the reference data coincides with the numerical results), but they suffer excessively when running under wake flows. Overall, the aerodynamic loads tend to be reduced when running under close proximity, ranging from 23% to a very significant 62% in the closest case.

In regards to the individual focus that is placed on the different aerodynamic devices, it has been found that the front wing experiences a sudden jump on downforce losses only when it enters into the closer wake region. On the other hand, the rear wing massively suffers from long distances, but its losses are way more linear and moderate. As for the diffuser, it is found that it represents the most affected aerodynamic device, as its performance is reduced from a considerable 25% to a huge 70%. This evidences that the conception of the diffuser and vortex management under the floor becomes critically compromised under wake flow conditions.

The modern performance of Ground Effect by means of vortices management represents a very unique and complex way of modelling modern aerodynamics, but, at the same time notably compromise the performance of the cars when an overtaking maneuver is intended. For this reason, it is possible to guarantee that the FIA changes in the current regulations is considered to be adequate, necessary, and justified in terms of aerodynamic necessities.

However, it is essential to comment that the present study presents some research limitations, as CFD methodology always requires an experimental validation. Accordingly, directions in further research point towards experimental validation (involving the use of wind tunnels with movable ground and spinning tires) in order to ensure that the data gathered reflect the reality properly.

**Author Contributions:** A.G. configured and performed the numerical simulations as well as the analysis of the results and the paper elaboration. R.C. supervised the study and reviewed the final version of the paper. All authors have read and agreed to the published version of the manuscript.

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

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

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **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* **Study of Surface Roughness E**ff**ect on a Blu**ff **Body—The Formation of Asymmetric Separation Bubbles**

#### **Alex Mendonça Bimbato 1,\* , Luiz Antonio Alcântara Pereira <sup>2</sup> and Miguel Hiroo Hirata <sup>3</sup>**


Received: 31 October 2020; Accepted: 18 November 2020; Published: 20 November 2020 -

**Abstract:** Turbulent flows around bluff bodies are present in a large number of aeronautical, civil, mechanical, naval and oceanic engineering problems and still need comprehension. This paper provides a detailed investigation of turbulent boundary layer flows past a bluff body. The flows are disturbed by superficial roughness effect, one of the most influencing parameters present in engineering applications. A roughness model, recently developed by the authors, is here employed in order to capture the main features of these complex flows. Starting from subcritical Reynolds number simulations (Re = 1.0 × 10<sup>5</sup> ), typical phenomena found on critical and supercritical flow regimes are successfully captured, like non-zero lift force and its direction change, drag crisis followed by a gradual increase on this force, and separation and stagnation points displacement. The main contribution of this paper is to present a wide discussion related with the temporal history of aerodynamic loads of a single rough circular cylinder capturing the occurrence of asymmetric separation bubbles generation. The formation of asymmetric separation bubbles is an intrinsic phenomenon of the physical problem, which is successfully reported by our work. Unfortunately, there is a lack of numerical results available in the literature discussing the problem, which has also motivated the present paper. Previous study of our research group has only discussed the drag crisis, without to investigate its gradual increase and the change on lift force direction. Our results again confirm that the Lagrangian vortex method in association with Large-Eddy Simulation (LES) theory enables the development of two-dimensional roughness models.

**Keywords:** bluff body aerodynamics; roughness model; boundary layer separation; vortex shedding; Lagrangian vortex method

#### **1. Introduction**

The technologic and science enhancement in many fields of knowledge is because of deep analysis and consequent comprehension of flow fundamentals around bodies of arbitrary shape. Many of the advances observed in ocean, civil and wind engineering have been possible because of studies conducted over time in the field of aerodynamics. The importance of knowing and mastering this subject can be easily recognized when it is observed that, most of bodies present in situations of practical interest for engineering, are exposed to air or water flows. These flows are typical examples of fluid-structure interaction problems, where the transition to turbulence is undeniably a complex phenomenon of nonlinear hydrodynamics. Such phenomenon arouses great scientific interest, has a

great impact on engineering applications (aeronautical, civil, mechanical, naval and oceanic), and has still its origin linked to the development of instabilities associated with the interaction of two shear layers of opposite signals; the latter are caused by separation points originated from the bluff body surface [1], like a circular cylinder, which is commonly studied because of its rich collection of flow phenomena. On the other hand, the transition to turbulence originates in the viscous wake region and, as the Reynolds number, Re, increases, it occurs more and more close to the body surface, passing through the two shear layers formed from the separation points (commonly defined by an angle θsep), until reaching the boundary layer region of the body.

In the literature, there are traditional works discussing transition to turbulence e.g., [2–10], where the surface roughness influence has also been included. The current scope of the research theme about "surface roughness effect" indicates few experimental and numerical studies e.g., [11–18]. The importance of the more relevant works will be discussed throughout the paper. The occurrence of the transition to turbulence into the boundary layer region is particularly important for practical applications, since the flow has a greater amount of momentum supporting the adverse pressure gradient; furthermore, the separation points move downstream of the body surface, causing the so-called drag crisis, which consists of a reduction in the drag force, that is opposed to body motion. The physics involved in this whole process is so complex, and under certain Reynolds number, the separation of the boundary layer is followed by a reattachment and a new separation again; the closed thin separated region was termed in the classical literature as a separation bubble.

According to Tani [2], the formation of separation bubbles causes low drag force values. Although the circular cylinder is a symmetrical body, the formation of separation bubbles does not occur simultaneously on the upper and lower surfaces of the body, which causes the appearance of an asymmetric pressure distribution and, as consequence, a non-zero lift force, contrary to what is expected to be obtained in any flow around a symmetrical bluff body.

The physical phenomenon above mentioned manifests at high Reynolds numbers, so that the presence of small disturbances in the flow will be amplified changing the bluff body aerodynamics. One of the types of disturbance most encountered in practical engineering problems refers to the surface roughness effect of a body, which affects heat exchangers efficiency, ship propellers performance, aerodynamic of sport materials, flow around offshore structures and wind turbines performance. The use of surface roughness effect on fluid flow requires to change the body geometry, being, therefore, classified as a passive method of vortex shedding control.

In Achenbach [3], the surface roughness effect is experimentally studied on the flow past a circular cylinder. He used a high-pressure wind tunnel, where high Reynolds number flows up to Re = 3 × 10<sup>6</sup> could be obtained. It was presented the classical C<sup>D</sup> × Re curve (*C*<sup>D</sup> is the drag coefficient) for the flow past a circular cylinder dividing it into four flow regimes: subcritical (Re < 2.0 × 10<sup>5</sup> − 5.0 × 10<sup>5</sup> ), critical (Re 2.0 × 10<sup>5</sup> − 5.0 × 10<sup>5</sup> ), supercritical (2.0 × 10<sup>5</sup> − 5.0 × 10<sup>5</sup> < Re ≤ 3.5 × 10<sup>6</sup> ) and transcritical (Re > 3.5 × 10<sup>6</sup> ); the latter is nowadays called pos-critical, as will be illustrated latter. According to him, the surface roughness effect does not affect the subcritical flow regime. On the other hand, as the Reynolds number increases, the drag force suddenly drops until reach a minimum value; this range of Reynolds number is within the critical flow regime. Further, when the Reynolds number exceeds the minimum value of the drag force, the drag force grows up again (supercritical flow regime) reaching a nearly constant value (transcritical flow regime). It was also observed that the drag force, within the transcritical flow regime, becomes bigger with increasing roughness. In regarding of the percentage of the skin friction (or viscous) component of the drag force, it was about 1% in the subcritical flow regime and 3% for the other three flow regimes.

Bearman [5] examined the flow past a smooth circular cylinder for the Reynolds number ranging from Re = 1.0 × 10<sup>5</sup> to 7.5 × 10<sup>5</sup> aiming to investigate the base pressure coefficient behavior depending on the Reynolds number. An interesting feature identified by him is that drag coefficient has almost the same value of the base absolute pressure value. Besides, according to the experimental data, the separation bubble formed on one side of the body provoked a discontinuity recorded at Re = 3.4 × 10<sup>5</sup> of the base pressure coefficient; during this phase, the lift coefficient, *C*L, calculated from the pressure distribution, was about 1.3. The separation bubble on the other side formed at Reynolds number flow of Re = 3.8 × 10<sup>5</sup> , and the lift coefficient was CL=0.0, being measured as soon as this happens and remains in this way, as the Reynolds number increases.

Guven et al. [7] presented in one same diagram the results obtained in their experimental tests with other results available in the literature [3,8,9] for circular cylinders with different roughness heights. It was showed scatter results with up to 60%, even for the same roughness and Reynolds numbers. This happens because, although the Reynolds number governs this kind of flow, it is so difficulty to guarantee that two experimental tests have exactly the same roughness surface; besides that, there are some influencing parameters which can cause a premature initiation of a flow regime or even its modification.

According to Zdravkovich [10], the above referred influencing parameters are: (a) free stream turbulence, described by its intensity and scale; (b) surface roughness effect, described by its mean relative surface roughness, ε ∗ /D<sup>∗</sup> , and its texture; (c) wall confinement (blockage effect); (d) aspect ratio; (e) presence of other boundary near the body; (f) end condition of the body; (g) body oscillation effect. The free stream turbulence and surface roughness influence are the most common disturbances in all engineering applications. The wall confinement and the aspect ratio are important influencing parameters that interfere in the comparison between experimental and numerical results, especially when the surface roughness effect is an influencing parameter too [17,18]. These features help us to understand the scatter in the experimental data presented by Guven et al. [7].

In a recent paper, Bimbato et al. [14] developed a Lagrangian vortex method using two-dimensional roughness model, which was based on the physics involved in turbulent boundary layer flows. The roughness surface effect was simulated without make changes on geometry of the body, what means that the investigated circular cylinder surface remains smooth. The key idea is that a point set strategically located close to the body surface injects momentum in its boundary layer aiming to represent surface roughness effect. Their strategy successfully shows that the separation of boundary layer is delayed. They studied the incompressible flow around both smooth and rough circular cylinders at Re = 1.0 × 10<sup>5</sup> , and it was only reported the drag crisis. Based on aerodynamic loads computations, it was concluded that their methodology, starting from a typical subcritical Reynolds number flow, was able to capture supercritical flow patterns. Overall, the drag coefficient computed just fall and it did not present the gradual increase, which is a typical feature of supercritical flows. Also, it was not presented results about the lift force, which would be interesting since a change on the lift force direction is expected [2].

The present paper, therefore, contributes to the literature using the same roughness model developed by Bimbato et al. [14] in order to capture physical details concerning the complex flow around a single rough circular cylinder. Put in other words, the main contribution here states on physical discussions involving the drag and lift forces connected with the occurrence of asymmetric separation bubbles generation. The generation of asymmetric separation bubbles is an intrinsic phenomenon to the physics of the problem [2]. In Section 4, it will be reported that our numerical results have captured that phenomenon, which has rarely been discussed by other numerical investigations. The results also confirm that two-dimensional Lagrangian vortex method with both LES and roughness models is able to capture not only the drag crisis, but also the non-zero lift force, which is typical of critical Reynolds number flows, and the gradual increase on drag coefficient, which is a typical feature of supercritical flows. Furthermore, our results show a delay in the boundary layer separation and also a stagnation point displacement.

Finally, it can be said that numerical simulations using Lagrangian vortex method have been employed to study problems in several areas of engineering. It is possible to mention a large number of studies, from the classics e.g., [19–22] to the most current ones e.g., [23,24]. Despite the high Reynolds numbers involved in practical engineering problems, the use of two-dimensional numerical simulations is necessary not only to develop and validate computational codes, but also to help understand the complex physical phenomena involved in such flows, and also to help conservative designs for engineering applications. A wide range of two-dimensional studies can be found in literature e.g., [14,17,18,25–31] which includes relevant contributions of our research group.

#### **2. Mathematical Formulation**

Figure 1 illustrates the two-dimensional fluid domain, Ω, to be studied, and defined by surface S=S1∪S<sup>2</sup> (S<sup>1</sup> refers to the circular cylinder boundary and *S<sup>2</sup>* is the boundary far from the body). In the same figure, U<sup>∗</sup> is the incompressible inlet flow, D<sup>∗</sup> defines the outer cylinder diameter, θ denotes upper angular position, x <sup>∗</sup> and y <sup>∗</sup> are the global coordinate system, and the symbol \* means dimensional quantities.

**Figure 1.** Definition of the fluid domain for the investigated problem, Ω.

The flow to be studied is assumed turbulent and, thus, a LES modeling is employed to disconnect the large eddies from the smaller ones. As practical result, the continuity and the momentum equations can be filtered resulting in Equations (1) and (2), such as, respectively:

$$\frac{\partial \overline{\mathbf{u}\_i^\*}}{\partial \mathbf{x}\_i^\*} = \mathbf{0} \tag{1}$$

$$\frac{\partial \overrightarrow{\mathbf{u}\_{i}^{\*}}}{\partial \mathbf{t}^{\*}} + \frac{\partial}{\partial \mathbf{x}\_{\cdot \cdot}^{\*}} (\overrightarrow{\mathbf{u}\_{i}^{\*}} \overrightarrow{\mathbf{u}\_{\cdot}^{\*}}) = -\frac{1}{\rho} \frac{\partial \overrightarrow{\mathbf{p}^{\*}}}{\partial \mathbf{x}\_{\cdot}^{\*}} + 2 \frac{\partial}{\partial \mathbf{x}\_{\cdot \cdot}^{\*}} [(\boldsymbol{\upsilon} + \boldsymbol{\upsilon}\_{\rm t}) \overrightarrow{\mathbf{S}\_{\overline{\mathbf{u}}}}] \tag{2}$$

In Equations (1) and (2), the *i-*th filtered velocity component of the flow field is identified by u ∗ i (being u ∗ i = u ∗ <sup>i</sup> + u ′∗ <sup>i</sup> with u ′∗ i representing the velocity fluctuation), ρ defines the density of fluid, p ∗ represents the pressure filtered field (p ∗ = p <sup>∗</sup> + p ′∗ ; p ′∗ is related to the pressure fluctuation), υ defines the fluid kinematic viscosity, υ<sup>t</sup> defines the local eddy viscosity coefficient, and S ∗ ij characterizes the deformation tensor of the filtered field [32].

The present methodology generates vortex blobs around the body surface during each time stepping aiming, thus, to solve the large eddy phenomena through the Equations (1) and (2). On the other hand, the small eddy phenomena are simulated using the concept of eddy viscosity coefficient (a typical LES approach). According to Lesieur and Métais [33], the local eddy viscosity coefficient at a point **x** ∗ in the flow field, and at an instant t ∗ , can be calculated via the local kinetic energy spectrum:

$$\text{tr}\_{\mathbf{t}} \begin{pmatrix} \mathbf{x}^\* \ \boldsymbol{\Delta}^{+\*} \end{pmatrix} = 0.105 \mathbf{C}\_{\mathbf{k}}^{-3/2} \boldsymbol{\Delta}^{+\*} \sqrt{\mathbf{F}\_2^\*(\mathbf{x}^\*, \boldsymbol{\Delta}^{+\*}, \mathbf{t}^\*)} \tag{3}$$

being C<sup>k</sup> = 1.4 the constant of Kolmogorov and F ∗ 2 (**x** ∗ , ∆ +∗ , t ∗ ) the local second-order velocity structure function of the filtered field; the latter is defined as follows [33]:

$$\overrightarrow{\mathbf{F}}\_{2}(\mathbf{x}^{\*}, \Delta^{+\*}, \mathbf{t}^{\*}) = \left\| \overleftarrow{\overline{\mathbf{u}}^{\*}(\mathbf{x}^{\*}, \mathbf{t}^{\*}) - \overline{\mathbf{u}}^{\*}(\mathbf{x}^{\*} + \mathbf{r}^{\*}, \mathbf{t}^{\*})} \right\|\_{|\mathbf{r}^{\*}| = \Delta^{+\*}}^{2} \tag{4}$$

In Equation (4), and for the three dimensional space, **u** ∗ (**x** ∗ , t <sup>∗</sup>) − **u** ∗ (**x** <sup>∗</sup> + **r** ∗ , t <sup>∗</sup>) represents the average speed differences between the center of a sphere of radius |**r** ∗ | = ∆ +∗ , located at **x** ∗ , and points placed on the sphere surface (∆ +∗ separates the scales that are computed by the numerical method of the ones that should be modeled). In this approach the small scales are assumed to be homogeneous and isotropic, and the point of the flow field where one wants to quantify the turbulent activity is assumed as the center of the sphere [33].

The impermeability boundary condition on the cylinder surface (at S<sup>1</sup> in Figure 1) is established such as:

$$
\overline{\mathbf{u}\_{\mathbf{n}}^{\*}} - \mathbf{v}\_{\mathbf{n}}^{\*} = \mathbf{0} \text{ on } \mathbf{S}\_{\mathbf{l}} \tag{5}
$$

being u ∗ n the fluid normal velocity and v ∗ n the normal velocity of the body surface.

The no-slip boundary condition on the cylinder surface (at S<sup>1</sup> in Figure 1) is imposed through the following equation:

$$
\overrightarrow{\mathbf{u}}\_{\tau}^{\*} - \mathbf{v}\_{\tau}^{\*} = \mathbf{0} \text{ on } \mathbf{S}\_{1} \tag{6}
$$

being u ∗ τ the fluid tangential velocity and v ∗ τ the tangential velocity of the body surface.

The far away boundary condition is given by:

$$\left| \overline{\mathbf{u}} \right| \to \mathbf{U}^\* \text{ on } \mathbf{S\_2} \text{ (Figure 1)}\tag{7}$$

It is very practical to make dimensionless all the quantities previously presented, aiming to obtain generality gain from the numerical simulations. Thus, in the above equations, D<sup>∗</sup> and U<sup>∗</sup> are assumed as length and velocity scales, respectively. As consequence, the non-dimensional time is defined by t = t <sup>∗</sup>U<sup>∗</sup> /D<sup>∗</sup> .

#### **3. The Lagrangian Vortex Method**

In order to solve numerically the mathematical problem presented on Section 2, the first step is to discretize the circular cylinder surface shown in Figure 1. This is done through the so-called panel method [34], which consists on discretizing any solid boundary in panels, over which singularities are distributed in order to guarantee one of the boundary conditions presented in Section 2 at the "pivotal point" (center point) of each panel (Figure 2). In the present approach, the cylinder surface is discretized and represented by NP flat panels with source distribution of constant density. The sources generation is necessary to ensure the impermeability boundary condition Equation (5) in association with the mass conservation of the problem.

**Figure 2.** Circular cylinder surface discretized by NP flat panels (in this example, NP = 6).

The Lagrangian vortex method is nothing more than a boundary layer model, which consists on discretizing a flow property (in this case, the vorticity filtered field) using vortex blobs (or vortex particle elements), such that [35]:

$$\overline{\boldsymbol{\varpi}}(\mathbf{x}, \mathbf{t}) = \sum\_{\mathbf{k}=1}^{\text{NV}} \frac{\Gamma\_{\mathbf{k}}}{\pi \,\sigma\_{0\_{\mathbf{k}}}^2} \exp\left(-\frac{|\mathbf{x}|^2}{\sigma\_{0\_{\mathbf{k}}}^2}\right) \tag{8}$$

In Equation (8), ω (ω = ∇ × **u**) represents the vorticity filtered field; NV is the total number of Lamb discrete vortices necessary to simulate the viscous wake; and Γ<sup>k</sup> defines the strength of the k-th vortex blob placed at a pivotal point of a flat panel. The discrete vortices generation is necessary to ensure the no-slip boundary condition Equation (6) in association with the global circulation conservation of the problem. In addition, σ<sup>0</sup> defines the Lamb vortex core, which is measured such as [36,37] (see also Figure 3):

$$
\sigma\_0 = 1.41421 \sqrt{\frac{\Delta t}{\text{Re}}} \tag{9}
$$

**Figure 3.** The boundary layer model for a smooth surface.

In Equation (9), the time stepping, ∆t, is computed from an estimate of the velocity and advective length of the flow.

As previously described, the potential problem is solved using the panel method through constant-density source distributions [34]. On the other hand, the vorticity effect is included into the problem by taking the curl of Equation (2), which, in association with Equation (1), results in the well-known vorticity transport equation [38]. In two dimensions, the mentioned equation is scalar and the pressure term is eliminated, such as:

$$\frac{\partial \overline{\boldsymbol{w}}}{\partial \mathbf{t}} + (\overline{\mathbf{u}} \cdot \boldsymbol{\nabla}) \overline{\boldsymbol{w}} = \frac{1}{\text{Re}\_{\mathbb{C}}} \mathbf{V}^{2} \overline{\boldsymbol{w}} \tag{10}$$

being Re<sup>c</sup> a local Reynolds number including the local eddy viscosity coefficient, and defined as [25]:

$$\text{Re}\_{\mathbf{c}\_{\parallel}}(\mathbf{t}) = \frac{\mathbf{U}^{\star}\mathbf{D}^{\star}}{\mathbf{v} + \mathbf{v}\_{\mathbf{f}\_{\parallel}}(\mathbf{t})} \tag{11}$$

Equation (10) simulates the vorticity filtered field dynamic, being solved to each vortex blob during each time stepping. Figure 3 sketches the typical generation process of vortex blobs used for a smooth surface. Each vortex blob is positioned tangent to the pivotal point of the panel to ensure the no slip boundary condition of the problem Equation (6), as already commented.

The present roughness model modifies the generation process of vortex particle elements without change the body surface. It is known that the roughness effect on flows is to stimulate the turbulence development. Thus, the mentioned generation process is altered, in order to add an instantaneous inertial effect to each nascent vortex blob. This additional inertial effect works as an injection of momentum into the boundary layer and, as consequence, the flow can support the unfavorable pressure gradient present in the boundary layer for longer, in agreement with the physics involved in this kind of complex flows. Since roughness and turbulence are intimately related, the momentum injection imposed by the present roughness model is obtained from turbulence modeling developed by Lesieur and Métais [33], which was originally implemented in a two-dimensional Lagrangian manner by Alcântara Pereira et al. [25], and with improvements made by Bimbato et al. [37].

Thus, the generation process of vortex blobs is modified considering the turbulent activity around "panel's shedding point" (Figure 2), what is not taken into account for smooth surfaces (Figure 3). The turbulent activity in the vicinity of the shedding point of the i-th panel is computed adopting a semicircle, with radius k**b**k = 2ε− σ0, and centered on i-th shedding point (ε = ε ∗ /D<sup>∗</sup> is the mean relative surface roughness). According to Bimbato et al. [14], a set of NR points is placed on the semicircle and the average speed differences, required to compute the second-order velocity structure function of the filtered field Equation (4), are calculated between the center of the semicircle (the shedding point) and the points on it (Figure 4). The second-order velocity function is obtained through the following adapted expression [14]:

$$\overline{\mathbf{F}}\_{\text{2}\_{\text{l}}}(\mathbf{t}) = \frac{1}{\text{NR}} \sum\_{\mathbf{w}=1}^{\text{NR}} \left\| \overline{\mathbf{u}}\_{\text{l}\_{\text{l}}}(\mathbf{x}\_{\text{l}}, \mathbf{t}) - \overline{\mathbf{u}}\_{\text{l}\_{\text{W}}}(\mathbf{x}\_{\text{l}} + \mathbf{b}, \mathbf{t}) \right\|\_{\text{W}}^{2} (1 + \varepsilon) \tag{12}$$

**Figure 4.** Turbulent activity computation in the vicinity of each shedding point.

In Equation (12), **u**t<sup>i</sup> (**x**i , t) is the total velocity on shedding point of the i−th panel (which is located at **x**<sup>i</sup> ) at time t, **u**t<sup>w</sup> (**x**<sup>i</sup> + **b**, t) is the total velocity at each one of the NR points, defined around the i−th shedding point (that is, on the semicircle of radius k**b**k = 2ε − σ0) at time t, and (1 + ε) is idealized to simulate the kinetic energy gain because of the roughness interference. Figure 4 gives us an idea of the turbulent activity computation around the shedding points, which is the velocity fluctuation in the boundary layer.

It is noted from Equation (12) that, at each time stepping, there is a different velocity fluctuation around the shedding point of the i−th panel. Therefore, there is also different local eddy viscosity coefficients associated with these velocity fluctuations, Equations (3) and (4). Thus, the eddy viscosity coefficient at i−th shedding point, at time t, is computed such as:

$$
\sigma\_{\mathbf{l}\_{\parallel}}(\mathbf{t}) = 0.105 \, \text{C}\_{\mathbf{k}}^{-3/2} \, \sigma\_{\mathbf{0}\_{\text{k}}} \sqrt{\overline{\mathbf{F}\_{\mathbf{2}\_{\parallel}}(\mathbf{t})}} \tag{13}
$$

In Equation (13), σ0<sup>k</sup> is the vortex core of the k−th vortex blob placed at i−th shedding point Equation (9).

Since the eddy viscosity is associated with each time and with each panel shedding point, a local Reynolds number must be computed at each time stepping, as defined in Equation (11). Therefore, Rec<sup>i</sup> (t) represents the Reynolds number modified locally (at i−th shedding point) at time t, and this modification only is computed if the surface roughness effect cannot be neglected, that is, for υt<sup>i</sup> (t) , 0.

Equation (9) shows how the Reynolds number of the flow and the Lamb vortex core are related. Thus, if the surface roughness effect is present in the boundary layer flow (if υt<sup>i</sup> (t) , 0), the Reynolds number is modified locally Equation (11). The numerical effect is to change the core radius, σ0, of each nascent vortex blob by satisfying Equation (14) instead of Equation (9). Therefore, the generation process of discrete vortex elements, illustrated in Figure 3, is now modified to include an additional inertial effect to each nascent vortex blob (it can also be interpreted by an increase of strength ∆Γ in Figure 5).

**Figure 5.** The boundary layer model to simulate a rough surface.

The new Lamb vortex core assumes the following form [14]:

$$
\sigma\_{0c\_k}(\mathbf{t}) = 1.41421 \sqrt{\frac{\Delta \mathbf{t}}{\mathrm{Re}} \bigg( 1 + \frac{\mathbf{v}\_{\mathbf{t}\_i}(\mathbf{t})}{\mathbf{v}} \bigg)}\tag{14}
$$

In Equation (14), σ0c<sup>k</sup> (t) is the vortex core of k−th vortex blob, placed at i−th shedding point at time t, and modified by surface roughness effect (compare Figure 3 with Figure 5). The roughness model is based only on physics of turbulent flows.

The procedure above described to represent the body surface (sources generation) and the vorticity field (creation of nascent vortex blobs) is mathematically defined by two different kinds of linear algebraic equation systems; their solution is iterative and they are solved using the method of least squares, and can be easily linked with the roughness model, when activated.

The set of discrete vortex elements, positioned at shedding points, is necessary to simulate the vorticity filtered field dynamics according to vorticity Equation (10). With this purpose, Chorin [19] idealized an algorithm that splits the advective-diffusive operator of Equation (10) such as, respectively:

$$\frac{\partial \overline{\boldsymbol{\omega}}}{\partial \mathbf{t}} = \frac{\partial \overline{\boldsymbol{\omega}}}{\partial \mathbf{t}} + (\overline{\mathbf{u}} \cdot \nabla) \overline{\boldsymbol{\omega}} = \mathbf{0} \tag{15}$$

$$\frac{\partial \overline{\boldsymbol{w}}}{\partial \mathbf{t}} = \frac{1}{\text{Re}\_{\mathbb{C}}} \mathbf{V}^{2} \overline{\boldsymbol{w}} \tag{16}$$

Equation (15) implies that the vortex cloud is advected as a set of material particles of the flow in a typical Lagrangian manner.

The solution of the advection problem Equation (15) is given by integrating each vortex blob path equation, and using an explicit Euler scheme, in the following form:

$$\frac{d\mathbf{x}\_{\mathbf{k}}}{dt} = \overline{\mathbf{u}}\_{\mathbf{k}}(\mathbf{x}\_{\mathbf{k}'}\mathbf{t}), \; \mathbf{k} = 1, \; \text{NV} \tag{17}$$

In Equation (17), **u**t<sup>k</sup> (**x**k, t) represents the velocity filtered field computed at position occupied by k−th vortex blob. Thus, the velocity field of the problem is composed by three contributions, i.e., (i) the incident flow (or the uniform flow) velocity, ui(**x**, t); the cylinder surface (using the source panels distribution), ub(**x**, t); the viscous wake (using the vortex cloud), uv(**x**, t).

For the incident flow, shown in Figure 1, its components take the form:

$$
\overline{\mathbf{u}\mathbf{i}\_1} = \mathbf{1} \tag{18}
$$

$$
\overline{\text{uċ}\_2} = 0 \tag{19}
$$

The velocity induced by NP flat panels with constant-density source distribution at position occupied by k−th vortex blob is computed in the following form [34]:

$$\overline{\mathbf{u}\mathbf{b}\_{\mathbf{k}}^{\mathrm{n}}}(\mathbf{x}\_{\mathbf{k}}\text{ t}) = \sum\_{\mathbf{i}=1}^{\mathrm{NP}} \sigma\_{\mathbf{i}} \mathbf{c}\_{\mathbf{k}\mathbf{i}}^{\mathrm{n}} [\mathbf{x}\_{\mathbf{k}}(\mathbf{t}) - \mathbf{x}\_{\mathbf{i}}], \; n = 1, \; \mathbf{2}; \; k = 1, \; \mathrm{NV} \tag{20}$$

In Equation (20), σ<sup>i</sup> is the source strength per panel length and c n ki[**x**k(t) <sup>−</sup> **<sup>x</sup>**<sup>i</sup> ] represents the n−th component of the velocity induced at the k−th vortex blob by the i-th source flat panel.

The third contribution for the velocity filtered field computation is because of the vortex cloud (i.e., the vortex-vortex interaction, obtained through the Biot-Savart law [35]), which takes the following form:

$$\mathbf{u}\overline{\mathbf{u}\mathbf{v}\_{\mathbf{k}}^{n}}\left(\mathbf{x}\_{\mathbf{k}},\mathbf{t}\right)=\sum\_{\mathbf{j}=1}^{\text{NV}}\Gamma\_{\mathbf{j}}\mathbf{c}\_{\mathbf{k}\mathbf{j}}^{n}[\mathbf{x}\_{\mathbf{k}}(\mathbf{t})-\mathbf{x}\_{\mathbf{j}}(\mathbf{t})]\_{\text{f}}\text{ }n=1\text{, 2: }k=1\text{, NV}\tag{21}$$

In Equation (20), Γ<sup>j</sup> is the strength of the j−th vortex blob and c n kjh **x**k(t) − **x**j(t) i represents the n−th component of the velocity induced at the k−th vortex blob by other vortex blob placed at **x**<sup>j</sup> .

Therefore, the velocity filtered field,**u**<sup>t</sup> (**x**, t), is obtained by the summation of Equations (18)–(20) and its computation is used:


Once advective equation Equation (15) has been solved, to properly simulate the vorticity field dynamics, it is necessary to solve, within the same time stepping, the diffusion equation (Equation (16)) to include molecular viscous effect. In the present numerical method, Equation (16) is solved through the random walk method [21], which was originally developed by Chorin [19]. The integral solution of the random walk method is [19]:

$$\overline{\boldsymbol{w}}(\mathbf{x},\mathbf{t}) = \int\_{-\infty}^{\infty} [\mathbf{G}(\mathbf{x},\mathbf{y},\mathbf{t}) - \mathbf{G}(\mathbf{x}-\mathbf{y},\mathbf{t})] \mathbf{f}(\mathbf{y}) d\mathbf{y} \tag{22}$$

being:

$$\mathbf{G}(\mathbf{x}, \mathbf{y}, \mathbf{t}) = \frac{1}{\sqrt{\frac{4 \pi \mathbf{t}}{\mathrm{Re}\_{\mathbb{C}}}}} \exp\left(\frac{-(\mathbf{x} - \mathbf{y})^2}{\frac{4 \mathbf{t}}{\mathrm{Re}\_{\mathbb{C}}}}\right) \tag{23}$$

Therefore, the diffusive displacement of the k−th vortex blob is computed, such as:

$$\varphi\_{\mathbf{k}}(\mathbf{t}) = \sqrt{\frac{4}{\text{Re}\_{\mathbf{c}}}} \ln(\frac{1}{\text{P}}) [\cos(2\,\pi\,\mathbf{Q})\_{\mathbf{x}} + \sin(2\,\pi\,\mathbf{Q})\_{\mathbf{y}}] \tag{24}$$

In Equation (24), P and Q are random numbers, which may assume values in the range: 0 < P < 1 and 0 < Q < 1.

The presence of the local Reynolds number (Rec) in Equation (23) indicates that turbulence manifestations must be computed into the diffusion process [25]. The essence of turbulence modeling used in this work consists on calculating the turbulent activity in points of the flow through the average speed differences existing around such points (Equation (4)). Instead of compute these differences between a center of a sphere and points on its surface, as established by Lesieur and Métais [33] for three-dimensional flows, Alcântara Pereira et al. [25] proposed to compute the average speed differences between a center of an annulus and points located between its internal and external radius. Both the center of the annulus and the points between its inner and outer radius are defined for each vortex blob presents in the computational domain simulating the vorticity field dynamics.

The adaptations proposed by Alcântara Pereira et al. [25] let to compute the second-order velocity structure function of the filtered field in the following form:

$$\overline{\mathbf{F}}\_{2\_{\mathbf{k}}} = \frac{1}{N} \sum\_{\mathbf{j}=1}^{N} \left\| \overline{\mathbf{u}}\_{\mathbf{k}}(\mathbf{x}\_{\mathbf{k}}) - \overline{\mathbf{u}}\_{\mathbf{i}}(\mathbf{x}\_{\mathbf{k}} + \mathbf{r}\_{\mathbf{j}}) \right\|\_{\mathbf{j}}^{2} \left( \frac{\sigma\_{0\mathbf{k}}}{\mathbf{r}\_{\mathbf{j}}} \right)^{2/3} \tag{25}$$

In Equation (25), **u**<sup>t</sup> is the total velocity induced on points of interest ( **u**t <sup>=</sup> ui <sup>+</sup> ub <sup>+</sup> uv), N is the number of discrete vortices inside the annulus and **r**<sup>j</sup> defines the distance between the center of annulus (vortex blob under analysis, i.e., the k−th particle) and the points located between the internal and the external radius of annulus (the j−th particle placed inside it).

A statistical analysis performed by Bimbato et al. [37] concluded that proper internal and external radius of annulus is, respectively: rint = 0.1σ0<sup>k</sup> and rext = 4.0σ0<sup>k</sup> (for smooth surfaces) or rint = 0.1σ0c<sup>k</sup> and rext = 4.0σ0c<sup>k</sup> (for rough surfaces).

The use of Lagrangian vortex method is justified by the ease of work with high Reynolds number simulations (there is no numerical instabilities, which are associated with Eulerian techniques). Furthermore, there is no mesh in the Lagrangian technique and it is kind simple to implement the

turbulence LES modeling using velocity differences instead of derivatives Equation (24), which must be closely related to the roughness model Equation (12).

It remains to calculate the aerodynamic coefficients (drag and lift forces). With this purpose, it is used the stagnation pressure definition, such as:

$$
\overrightarrow{\mathbf{Y}}^{\*} = \frac{\overrightarrow{\mathbf{P}}^{\*}}{\rho} + \frac{\overrightarrow{\mathbf{u}}^{\*2}}{2}; \overrightarrow{\mathbf{u}}^{\*} = \|\overrightarrow{\mathbf{u}}^{\*}\|\tag{26}
$$

In Equation (26), p ∗ defines the static pressure and u ∗ represents the velocity at any point into the fluid domain Ω (Figure 1). Then, the methodology to compute aerodynamic loads starts from a Poisson equation for the pressure, which is solved through the following integral formulation presented by Shintani and Akamatsu [39]:

$$\mathbf{H}\overline{\mathbf{Y}}\_{\mathrm{i}} - \int\_{\mathrm{S}\_{\mathrm{l}}} \overline{\mathbf{Y}} \nabla \overline{\mathbf{z}}\_{\mathrm{i}} \cdot \mathbf{e}\_{\mathrm{n}} \mathrm{d}\mathbf{S} = \iint\_{\Omega} \nabla \overline{\mathbf{z}}\_{\mathrm{i}} \cdot (\overline{\mathbf{u}} \times \overline{\mathbf{w}}) \mathrm{d}\Omega - \frac{1}{\mathrm{Re}} \int\_{\mathrm{S}\_{\mathrm{l}}} (\nabla \overline{\mathbf{z}}\_{\mathrm{i}} \times \overline{\mathbf{w}}) \cdot \mathbf{e}\_{\mathrm{n}} \mathrm{d}\mathbf{S} \tag{27}$$

In Equation (27), the static pressure can be computed at i−th calculation point, being H = 1.0 valid for the fluid domain and H = 1/2 valid for solid boundary, Ξ represents a fundamental solution of Laplace equation, and **e**<sup>n</sup> defines the unit vector normal to each flat panel used to discretize the cylinder surface.

Using the static pressure obtained from Equation (27), the drag and lift coefficients are evaluated, such as, respectively:

$$\mathbf{C\_D} = \sum\_{\mathbf{i}=1}^{\text{NP}} \mathbf{2} \left( \overline{\mathbf{p\_i}} - \mathbf{p\_{os}} \right) \Delta \mathbf{S\_i} \sin \beta\_{\mathbf{i}} = \sum\_{\mathbf{i}=1}^{\text{NP}} \mathbf{C\_{P\_i}} \Delta \mathbf{S\_i} \sin \beta\_{\mathbf{i}} \tag{28}$$

$$\mathbf{C}\_{\rm L} = -\sum\_{\rm i=1}^{\rm NP} \mathbf{2} \left( \overline{\mathbf{p}}\_{\rm i} - \mathbf{p}\_{\rm co} \right) \Delta \mathbf{S}\_{\rm i} \cos \beta\_{\rm i} = -\sum\_{\rm i=1}^{\rm NP} \mathbf{C}\_{\rm P\_{\rm i}} \Delta \mathbf{S}\_{\rm i} \cos \beta\_{\rm i} \tag{29}$$

In Equations (28) and (29), p<sup>∞</sup> is the reference pressure far from the solid boundary, ∆S<sup>i</sup> is the length of the i−th panel, and β<sup>i</sup> is the angle of the i−th panel.

Finally, the algorithm implemented essentially consists of eight steps in the following sequence: (i) simultaneous generation of source panels and nascent vortex blobs (the roughness model is used, when activated); (ii) velocity field computation at each vortex blob; (iii) static pressure computation at each pivotal point followed by integrated aerodynamic loads computation; (iv) advection of the vortex cloud; (v) vorticity diffusion (the LES modeling is used); (vi) reflection of vortex blobs, which eventually migrate into the cylinder; (vii) velocity field computation at each pivotal point to restore the boundary conditions of the problem for new generation of source panels and nascent vortex blobs; (viii) advance by time ∆t. In the last years, the authors of this paper have made an effort to develop the in-house code by using FORTRAN programming language.

#### **4. Results and Discussion**

This Section presents the numerical results obtained for the problem of the two-dimensional, incompressible and unsteady flow past a circular cylinder. The body surface was represented by NP = 300 source flat panels. That chosen number for the panels ensures the convergence for the potential solution of the problem and, as consequence, a refined discretization for the vorticity field, which is reflected on the aerodynamic loads computation. In the case of flow around smooth or rough cylinders for the Reynolds number value studied here (Re = 1.0 × 10<sup>5</sup> ), the drag coefficient is computed only by form (pressure) drag contribution. According to Achenbach [3] measurements, the friction drag contribution does not exceed 2–3% to the total drag. Thus, the friction drag computation has been omitted in this work.

In the present approach, it has been used NR = 21 points on each semicircle defined around each panel shedding point in order to compute the turbulent activity around it (Figure 4); this number is enough to obtain a reasonable value for the average speed differences Equation (12) in consonance with past study reported by Bimbato et al. [14].

In order to obtain good results, the time increment used is ∆t = 0.05 and an explicit Euler time-marching scheme is used Equation (17). All numerical simulations run until dimensionless time of t = 75, which is enough to reach a statistical equilibrium; the mean coefficient values (pressure distribution, drag and lift) are computed between 37.5 ≤ t ≤ 75.0. As previously presented by Oliveira et al. [18], the adopted criterion to compute aerodynamic loads has sufficiently been refined attaining the saturation state of the numerical simulations; this behavior is typical of a Lagrangian technique.

Table 1 presents a summary of some typical experimental results published in the literature and the numerical ones obtained using the present methodology. There is a remarkable scatter among the experimental results, even those referring to the smooth circular cylinder. The scatter among experimental results is because of influencing parameters, especially by the low aspect ratio and high blockage. The cylinder tested by Achenbach [3] presented an aspect ratio of L<sup>∗</sup> /D<sup>∗</sup> = 3.3 (L∗ is the cylinder length) and the wind tunnel blockage was D<sup>∗</sup> /B <sup>∗</sup> = 0.17 (B<sup>∗</sup> is the distance between tunnel walls). The cylinder tested by Achenbach and Heinecke [4] presented an aspect ratio of L ∗ /D<sup>∗</sup> = 6.75 and the wind tunnel blockage was D<sup>∗</sup> /B <sup>∗</sup> = 0.17, while Zhou et al. [12] used a configuration of cylinder with aspect ratio of L<sup>∗</sup> /D<sup>∗</sup> = 10. Thus, the comparison among experimental results and the present numerical results is not so fair, as previously discussed in Section 1. Unfortunately, there is a lack of numerical data of aerodynamic forces for the flow past a rough cylinder at upper subcritical Reynolds number of Re = 1.0 × 10<sup>5</sup> .


**Table 1.** Summary for the drag coefficient at Re = 1.0 × 10<sup>5</sup> .

Despite this, it is undeniable that the present numerical results have the same tendency as the experimental results of Achenbach [3] (Table 1). The roughness model used here is able to predict supercritical flow features from a subcritical Reynolds number flow simulation. This conclusion comes from the comparison between Figure 6a,b.

The present numerical results show that, for small superficial roughness (0.00 < ε ≤ 0.10 %), the drag coefficient suffers a very small reduction (about 2.9%), which occurs in the subcritical regime (Figure 6a). With the progressive increase in surface roughness (0.20% ≤ ε ≤ 0.45%), the roughness model injects momentum on the boundary layer in a manner that flows with higher Reynolds numbers are simulated; it is remarkable the drag coefficient sudden drops until reach a minimum value (Figure 6b), that is a typical feature of critical flows. The mean drag reduction is about 16.2%. For even greater roughness (ε > 0.45 %), the drag coefficient grows up gradually (about 4.5%), which is the main feature of the supercritical flow regime (Figure 6a).

**εε**

Δ

ε **Figure 6.** Flow regimes described by drag coefficient: (**a**) Classical representation of C<sup>D</sup> = C<sup>D</sup> (Re) curve for a cylinder (reproduced from [40] with permission; copyright © 2006 by World Scientific Publishing Co. Pte. Ltd.); (**b**) Form of C<sup>D</sup> = C<sup>D</sup> (ε) curve obtained by roughness model for a fixed Reynolds number (Re = 1.0 × 10<sup>5</sup> ).

ε ε ε Although the percentage drop in the drag coefficient measured by Achenbach [3] is slightly more than double that calculated in the present work, most of the features of transitional flows are qualitatively captured by present methodology. Table 2 shows that boundary layer separation point agree well with the drag coefficient behavior pointed out in Table 1 and Figure 6b, what means the separation point moves downstream with drag reduction (critical flow regime) and upstream with drag gradual increase (supercritical flow regime).


**Table 2.** Upper side angular position of boundary layer separation, θsep (Re = 1.0 × 10<sup>5</sup> ).

*Energies* **2020**, *13*, 6094

Furthermore, Figure 7 indicates that the lift coefficient,CL, oscillates regularly because of the periodic vortex shedding mechanism [41], and the lift force variation of the smooth cylinder is greater than that of the rough cylinder, which agrees with experimental tests conducted by Zhou et al. [12]; the dash-dot lines in Figure 7 indicates the mean amplitude of lift force. In the present work, no analysis regarding the vortex shedding frequency is presented, however, a detailed discussion can be easily recuperated [14,17,18].

ε ε **Figure 7.** Temporal history for drag and lift coefficients (Re = 1.0 × 10 5 ): (**a**) Smooth; (**b**) ε = 0.10%; (**c**) ε = 0.45%.

Figure 8 presents the pressure distributions around the cylinder surface at time instants represented by points A, B, C, D and E; these points are also shown in Figure 7. It can be noticed that the surface roughness effect causes irregular disturbances on instantaneous pressure distribution curves, which is reflected on time history of drag and lift forces.

**Figure 8.** Pressure distributions for particular time instants of the flow past a circular cylinder (Re = 1.0 × 10 5 ): (**a**) Smooth; (**b**) ε = 0.10%; (**c**) ε = 0.45%.

The mentioned disturbances occur because of inertial effect manifesting on the flow, which is provoked by roughness model. Thus, the boundary layer presents more momentum supporting the unfavorable pressure gradient, as expected, and, as consequence, the flow separation is delayed. This was shown in the joint analysis of Figure 6b and Table 2; now it can also be verified with the analysis of Figure 9.

**Figure 9.** Separation bubble identification (Re = 1.0×10 5 ): (**a**) Mean lift coefficient; (**b**) Stagnation point.

In Figure 9a, the non-zero lift force for the smooth cylinder (C<sup>L</sup> = 0.021) and for the less roughened cylinders (C<sup>L</sup> = 0.035 for ε = 0.10% and C<sup>L</sup> = 0.005 for ε = 0.20%) can be attributed to numerical rounding errors (panels method and vortex cloud contributions). It can be noted that the stagnation point of these three tests is located at θstag = 0.6 ◦ ; this happens because the discretization of the cylinder using the panel method starts with θ = 0 ◦ (Figure 1) and the panel pivotal point is placed at θ = 0.6 ◦ . It is important to emphasize that all the computations are done on pivotal points (Figure 2). So, this is the reason why θstag , 0 for the three cases mentioned.

On the other hand, it is notorious not only the non-zero lift force for the rougher cylinders but its direction change (C<sup>L</sup> = +0.100 and θstag = −1.8 ◦ for ε = 0.45%; C<sup>L</sup> = −0.132 and θstag = +1.8 ◦ for ε = 0.70%), see Figure 9a,b. The boundary layer transition is the only explanation of this occurrence on a symmetrical bluff body, like the circular cylinder. The involved physics indicates that first an asymmetric flow is caused because the boundary layer becomes turbulent at one side of the cylinder; furthermore, a separation bubble and a non-zero lift force can be identified. On the other side of the cylinder is generated another separation bubble, when the boundary layer becomes turbulent on that side. According to Kamiya et al. [6], a fully symmetric flow is not reach as soon as the second separation bubble is formed. They measured C<sup>L</sup> = 0.2 and θstag = +2 ◦ when the first separation bubble was formed, followed by C<sup>L</sup> = 0.0 and θstag = 0 ◦ at higher Reynolds number (Re ≥ 7.0 × 10 5 ). The critical flow regime is so unstable, and the lift force changing direction was also captured by Kamiya et al. [6] by using a smooth cylinder.

One more feature of the complex transitional flow is captured by the present methodology. According to Bearman [5], the drag coefficient is approximately equal to the measured base pressure coefficient for three representative Reynolds numbers, i.e., Re = 2.0 × 10 5 (no separation bubble), Re = 3.7 × 10 5 (one separation bubble) and Re = 4.0 × 10 5 (two separation bubbles). Figure 10 shows the comparison between drag and base pressure coefficients, calculated by the present Lagrangian approach.

**Figure 10.** Comparison between computed results of base pressure and mean drag coefficients (Re = 1.0 × 10 5 ).

It is remarkable that the present methodology agrees with the observations given by Bearman [5] (using the L ∗ /D<sup>∗</sup> = 12 and D<sup>∗</sup> /B ∗ = 0.06 configuration) for the critical flow regime. As expected, there is a disagreement between the computed values of base pressure and drag coefficients for the supercritical flow regime, because of increasing in drag coefficient.

In addition, Figure 11 presents the time-averaged pressure distribution for the flow around a cylinder with different surface roughness influences, where the experimental result was reported by Blevins [42] for smooth cylinder. It can be observed only small differences in the base pressure between the less rough cylinder and the smooth cylinder. As can be identified, the numerical result for the rough cylinder (ε = 0.45%) clearly shows a higher increase in the base pressure, which physically agrees with the higher drag reduction (Figure 6b).

**Figure 11.** Time-averaged pressure distribution around the circular cylinder surface (Re = 1.0 × 10 5 ).

Finally, Figure 12 illustrates the viscous wake developed downstream the cylinder in the end of numerical simulations (at t = 75). The blue points in the same figure represent instantaneous vortex blob distributions; they indicate that the vortical structures for ε = 0.45% are narrower than the

other ones, which is in accordance to what is expected from critical flow regimes. That conclusion is supported by upper and lower separation points behavior for ε = 0.45%, which are 93◦ (Table 2) and 85.8◦ , respectively. <sup>ε</sup>

ε ε ε ε **Figure 12.** Final position of the viscous wake at t = 75 (Re = 1.0 × 10<sup>5</sup> ): (**a**) Smooth; (**b**) ε = 0.10%; (**c**) ε = 0.20%; (**d**) ε = 0.45%; (**e**) ε = 0.70%.

#### **5. Conclusions**

This paper describes a two-dimensional Lagrangian vortex method blended with both LES and roughness models to capture changes in a bluff body aerodynamics. It can be observed that the roughness model employed in the present work is classified as a passive method of vortex shedding control. The passive control method implies to modify the body surface, which is reflected on its aerodynamic loads behaviour. In fact, the roughness model presented here acts as a passive method modifying the boundary layer flow by an injection of momentum into it (which occurs on rough surfaces), but without changing the discretization of the body surface.

In the present approach, the vorticity field is discretized and represented by a cluster of vortex blobs, which are generated during each time stepping from the body surface. The potential theory is adapted to compute the circular cylinder surface, which is discretized using flat panels with constant-density source distribution over them. The aerodynamic loads are calculated through an integral formulation derived from a Poisson equation for the pressure, where the instantaneous vorticity field contributes to the computation. The chosen number of flat panels ensures the convergence of the solution for the potential velocity field and, as consequence, establishes the number of nascent vortex blobs. As the non-dimensional time runs, the number of vortex blobs increases and thus satisfies the vorticity transport equation forming the von-Kármán vortex street. In an upper subcritical Reynolds flow, the drag coefficient amplitudes are significantly lesser then the lift ones (Figure 7a). This behaviour indicates that the numerical simulations have captured the vortex shedding mechanism properly; it is important to observe that the saturation state of the present numerical simulations was attained and it is supported by previous discussions given by Oliveira et al. [18], since the numerical approach is the same.

In regarding of the rough circular cylinder, the numerical results show that the present methodology is able to capture, even using two-dimensional simulations, many attributes of this complex flow such as: (i) the drag crisis followed by a gradual increase on the drag force; (ii) the non-zero lift force, since according to Zdravkovich [10] the low and high pressure on two sides of the cylinder interchange in different runs, because the side at which the separation is turbulent switches from one side to another occasionally; and (iii) the occurrence of asymmetric separation bubbles generation. The asymmetries are part of the physics involved in critical Reynolds number flows around a circular cylinder; that behavior was previously described by Zdravkovich [10], whom dedicated all his research to investigate different aspects of flows around cylinders.

In order to obtain the results showed here, the present work ran simulations with roughness heights of ε = 0.00% (smooth case), ε = 0.10%, ε = 0.20%, ε = 0.45% and ε = 0.70%. The roughness heights of ε = 0.10% and ε = 0.20% were performed to show that their effect does not change the flow patterns. Further, the tests jumped for ε = 0.45% (+0.25%) and then for ε = 0.70% (+0.25%). The critical (ε = 0.45%) and upper transition (ε = 0.70%) Reynolds number flows are so unstable, that it could be difficult to physically explain the results on a higher resolution area around ε = 0.45%.

Although the lift force behavior is an important feature of the critical/supercritical flow regimes, studies that focus on this force are scarce, which valorizes the present numerical results. The results for both predictions of boundary layer separation and stagnation point displacement also help us to report the sensitivity of the roughness model. Furthermore, each experimental study available in the literature uses a different technique to represent the surface roughness effect (wires, spheres, small holes, among others), and each technique is subject to different influencing parameters (free stream turbulence, aspect ratio, blockage and surface texture), which produces a great scatter in the measurements. As consequence, in most of situations it is not appropriate to compare two different experimental results or even an experimental and a numerical one.

Finally, the authors also intend to adapt the present vortex code to study other practical engineering situations, such as problems involving wake interference [43], ground effect mechanisms [17,18,27,28], vortex-induced vibrations [13,26] and thermal effect [44,45]. The present paper also aims to lay the theoretical foundation of the present numerical method for further algorithmic extensions to simulations of such fluid flows in three dimensions.

**Author Contributions:** Conceptualization, A.M.B., L.A.A.P. and M.H.H.; data curation, A.M.B.; formal analysis, A.M.B., L.A.A.P. and M.H.H.; funding acquisition, L.A.A.P.; investigation, A.M.B., L.A.A.P. and M.H.H.; methodology, A.M.B., L.A.A.P. and M.H.H.; project administration, A.M.B. and L.A.A.P.; resources, A.M.B. and L.A.A.P.; software using FORTRAN, A.M.B. and L.A.A.P.; validation, A.M.B.; visualization, A.M.B.; supervision, L.A.A.P.; writing—Original draft, A.M.B.; and writing—review and editing, A.M.B. and L.A.A.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by FAPEMIG, grant number APQ-02175-14.

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

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 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*

## **Control and Suppression of Vortex Shedding from a Slightly Rough Circular Cylinder by a Discrete Vortex Method**

**Marcos André de Oliveira <sup>1</sup> , Paulo Guimarães de Moraes <sup>1</sup> , Crystianne Lilian de Andrade <sup>1</sup> , Alex Mendonça Bimbato <sup>2</sup> and Luiz Antonio Alcântara Pereira 1,\***


Received: 31 July 2020; Accepted: 26 August 2020; Published: 31 August 2020

**Abstract:** A discrete vortex method is implemented with a hybrid control technique of vortex shedding to solve the problem of the two-dimensional flow past a slightly rough circular cylinder in the vicinity of a moving wall. In the present approach, the passive control technique is inspired on the fundamental principle of surface roughness, promoting modifications on the cylinder geometry to affect the vortex shedding formation. A relative roughness size of ε\*/*d*\* = 0.001 (ε\* is the average roughness and *d*\* is the outer cylinder diameter) is chosen for the test cases. On the other hand, the active control technique uses a wall plane, which runs at the same speed as the free stream velocity to contribute with external energy affecting the fluid flow. The gap-to-diameter varies in the range from *h*\*/*d*\* = 0.05 to 0.80 (*h*\* is the gap between the moving wall and the cylinder bottom). A detailed account of the time history of pressure distributions, simultaneously investigated with the time evolution of forces, Strouhal number behavior, and boundary layer separation are reported at upper-subcritical Reynolds number flows of *Re* = 1.0 × 10<sup>5</sup> . The saturation state of the numerical simulations is demonstrated through the analysis of the Strouhal number behavior obtained from temporal history of the aerodynamic loads. The present work provides an improvement in the prediction of Strouhal number than other studies no using roughness model. The aerodynamic characteristics of the cylinder, as well as the control of intermittence and complete interruption of von Kármán-type vortex shedding have been better clarified.

**Keywords:** bluff body; roughness model; Venturi effect; suppression hybrid control; Lagrangian description

#### **1. Introduction**

In the literature, "bluff body" is defined as being a structure that when immerse in a fluid flow will present significant proportion of its surface generating separated flow. This idea is also associated with characteristics of the flow around the body, especially the two shear layers of opposite signals formed from the separation points of the body [1]. Since 1900, numerous analytical, numerical, and experimental investigations have been conducted to study bluff body aerodynamics on either two-dimensional or axisymmetric shapes; the two most common bodies being the circular cylinder and the sphere. The more relevant results have contributed to develop important researches in aerospace/aeronautical, civil, marine, mechanical, and computer engineering.

The present paper aims to contribute with more discussions concerning the surface roughness effect on bluff body aerodynamics. An important practical engineering application of the flow around cylindrical structures is that of the fluid–elastic interaction between the flow and the structure exciting the body into flow-induced vibrations (FIV) [2–6]. Undoubtedly, the surface roughness effect must be used to control this non-linear hydrodynamic phenomenon. This is one of important motivations for future extension of our present research, which will also include the study of slender body aerodynamics [7]. Another one is planned to investigate heat transfer behavior in mixed convection of a fluid by a temperature particles method [8].

In the context of the present work, the literature has reported studies of surface roughness effect on flow past a bluff body with low frequency [2,5,6,9–19]. It is important to mention that numerical analysis of rough bluff body aerodynamics near a ground plane is very limited, which valorizes the recent methodology proposed by Alcântara Pereira et al. [10] and classified such as "hybrid control technique of vortex shedding" [20].

In fact, control and suppression of vortex shedding from a bluff body has been considered as one of the most important research areas in the field of aerodynamics and hydrodynamics applications, such as vibration of pipelines, interactions of currents and wave with offshore structures, suspension bridges and chimneys near tall buildings.

There are different situations where the vortex shedding may cease, and, of particular interest here, is that of a circular cylinder in the vicinity of a ground plane. To the best of our knowledge, the ground effect is governed by three mechanisms [21–28]: (i) The wake interference because of the intertwine of the body wake and the boundary layer formed on the ground, the latter is less influent despite several intensive studies reported so far. (ii) The three-dimensional effect, which presents a momentum transfer in the axial direction of the body leading to a lower drag force value as compared to the two-dimensional results. (iii) The blockage effect (or Venturi effect), which contributes to appearing big and small peaks during the temporal evolution of the drag curve. These peaks have been identified within the large-gap regime, i.e., *h*\*/*d*\* > 0.40, by Bimbato et al. [22] using a discrete vortex method implemented with Large Eddy Simulation (LES) theory. The large-gap regime is characterized by the presence of strong vortical structures, which are generated at the rear part of a bluff body [25]. In Section 3.2, these peculiar peaks' behavior on the drag curve will be again studied, now using the present numerical method, to contribute on discussion of the roughness model sensitivity. The applicability of the Venturi effect in aerodynamic models has been discussed in the literature [29,30].

Within the scenario above described, Roshko et al. [26] reported the aerodynamic forces behavior for the smooth circular cylinder in the vicinity of a fixed wall using a wind tunnel at high Reynolds number flow of *Re* = 2.0 × 10<sup>4</sup> . As the cylinder came close to the ground plane, it was reported that the drag force rapidly decreased; on the contrary, the lift force increased.

Zdravkovich [27] also reported aerodynamic forces behavior for the smooth circular cylinder near a fixed ground at high Reynolds numbers flows in the range of 4.8 × 10<sup>4</sup> < *Re* < 3.0 × 10<sup>5</sup> . As the gap-to-diameter ratio *h*\*/*d*\*, reduced to less than the thickness of the boundary layer δ\*/*d*\**(*δ\* is the boundary layer thickness) on the ground plane, it was identified a rapid decrease in the drag force. The drag variation was dominated by *h*\*/δ\* and the ratio *h*\*/*d*\* was less dominant. It was observed that the state of the boundary layer could interfere in the lift force, although the boundary layer thickness was less influent.

In other relevant research, Zdravkovich [28] investigated the drag force behavior of the smooth circular cylinder near a ground plane. The wall plane ran at the same speed as the free stream velocity in an upper-subcritical Reynolds number flow of *Re* = 2.5 × 10<sup>5</sup> into the critical flow regime. The experimental results showed contrast to all previous studies using fixed ground, e.g., [26,27], because practically no boundary layer was generated from the moving wall. Interestingly, the expected decrease of the drag force, as the ratio *h*\*/*d*\* decreased, did not appear. That behavior was attributed to the non-existence of the boundary layer formed on the ground or the high Reynolds number, or any other influencing factors, such as surface texture and structural vibration.

Nishino [25] reported experimental results of the flow around the smooth circular cylinder with an aspect ratio of 8.33 in a wind tunnel. Two upper-subcritical Reynolds numbers of *Re* = 0.4 × 10<sup>4</sup> and 1.0 × 10<sup>5</sup> were investigated during his experiments. The cylinder with and without end plates configurations were near and parallel to a wall plane running at the same speed as the free stream velocity. The moving wall effect eliminated the less influent effects of boundary layer formed on the ground and, therefore, the experimental results contributed to clarify the fundamental mechanisms of ground effect in more details. The experimental study also produced new perceptions into the physics of ground effect, and still nowadays serves as database to support numerical investigations. Of importance for the present work, is that, for the cylinder with end plates, the oil flow patterns were observed to be approximately two-dimensional. In contrast to the cylinder near a fixed ground configuration [26], the drag force rapidly decreased, as the ratio *h*\*/*d*\*decreased to less than 0.50 and became constant for *h*\*/*d*\*of less than 0.35.

Bimbato et al. [21,22] implemented an algorithm of the discrete vortex method with LES modeling to study the two-dimensional flow around the circular cylinder near a moving wall. Their numerical strategy was to represent the ground plane motion using a plane wall with no vorticity generation on it. The numerical results for an upper-subcritical Reynolds number of *Re* = 1.0 × 10<sup>5</sup> presented a good agreement with the experimental results reported by Nishino [25] using the cylinder with end plates configuration. The authors concluded that the Venturi effect almost completely suppressed the vortex shedding from the cylinder placed closer to the ground plane. Furthermore, the drag force decreased as consequence of the suppression.

In a recent paper, Alcântara Pereira et al. [10] proposed a hybrid control technique of vortex shedding, combining passive and active controls, to study the flow past the rough circular cylinder in the vicinity of a moving wall at upper-subcritical Reynolds number of *Re* = 1.0 × 10<sup>5</sup> . They successfully associated the methodologies developed by Bimbato et al. [11,21,22] focusing on the effect of higher relative roughness sizes, namely, ε\*/*d*\* = 0.0045 and 0.007, on flow dynamics of the cylinder at small-gap regime, which was identified in past investigation by Bimbato et al. [22] at *h*\*/*d*\* < 0.20. The small-gap regime is characterized by vortex shedding suppression. Their results shown an anticipation of the vortex shedding suppression when using ε\*/*d*\* = 0.007 at *h*\*/*d*\* = 0.10. Of importance, the Strouhal number completely vanished (*St* = 0.0) and it was observed the formation of two nearly parallel shear layers of opposite signals at the rear part of the cylinder. It is interesting to comment that Bimbato et al. [22] did not capture the complete interruption of vortex shedding, even at *h*\*/*d*\* = 0.05. This flow behavior will be investigated later in Section 3.2 to also contribute on discussion of the roughness model effect on it.

It is important to observe that Alcântara Pereira et al. [10] only computed the form component of the drag force. According to the discussions of Achenbach [9], the form (or pressure) component dominates the drag force on the smooth cylinder contributing more than 98% of the total drag force. On the other hand, the skin friction (or viscous) component of the drag force is responsible for the remaining 1–2%. In the literature, it is expected neither component of the drag force of the rough cylinder can be neglected. However, Achenbach [9] investigated the viscous drag force from a sand-roughened surface (ε\*/*d*\* = 0.0011 and 0.0045) and concluded that it contributed about 2–3% of the total drag force. That result showed a slight increase of the viscous drag force as compared as that of the smooth cylinder. Alcântara Pereira et al. [10] have reported differences above 10% between the smooth cylinder and other rough when integrating only static pressure, and, therefore, as discussed by them, the surface roughness effect has been captured by their numerical approach.

The main contribution of the present paper is shown that the methodology proposed by Alcântara Pereira et al. [10] also captures important changes in the flow dynamics of the slightly rough circular cylinder in the vicinity of a moving wall. Thus, the effect of small relative roughness size, namely, ε\*/*d*\* = 0.001, on flow dynamics around the cylinder is investigated and compared to smooth cylinder configurations. The chosen roughness size can be found, for example, in support columns of large offshore floating structures after a few years of operation. Put in other words, the main goal here is to

report understanding of the influence of surface roughness for control of vortex shedding frequency reduction and wake destructive behavior behind the rough cylinder placed near a moving wall at upper-subcritical Reynolds number of *Re* = 1.0 × 10<sup>5</sup> . Our highlight results are for temporal history of aerodynamic loads, Strouhal number behavior and separation point prediction. The focus is to track the Strouhal number behavior until it vanishes. Overall, the results are found agree with the physics expected for this kind of vortical flow. The Section 4 will summarize the key findings of this study. The focus is to track the Strouhal number behavior until it vanishes. Overall, the results are found agree with the physics expected for this kind of vortical flow. The Section 4 will summarize the key findings of this study. The rough circular cylinder aerodynamics has been reasonably reported in the literature [2,9,11,12,14,16,18], however, very little attention has been paid to the problem of the flow around the rough cylinder near a moving wall [10]. In recent past works [10,11], it has been reported that the effect of two-dimensional roughness model is much more sensitive than single turbulence modeling.

*Energies* **2020**, *13*, x FOR PEER REVIEW 4 of 24

temporal history of aerodynamic loads, Strouhal number behavior and separation point prediction.

. Our highlight results are for

The rough circular cylinder aerodynamics has been reasonably reported in the literature [2,9,11,12,14,16,18], however, very little attention has been paid to the problem of the flow around the rough cylinder near a moving wall [10]. In recent past works [10,11], it has been reported that the effect of two-dimensional roughness model is much more sensitive than single turbulence modeling. In general, numerical simulations of high Reynolds number flows around two-dimensional bluff bodies over predict aerodynamic forces behavior. However, the results are very important for applications of conservative designs in practical engineering problems, where higher integrated loads are computed, specially the drag force, in association with accurate vortex shedding

In general, numerical simulations of high Reynolds number flows around two-dimensional bluff bodies over predict aerodynamic forces behavior. However, the results are very important for applications of conservative designs in practical engineering problems, where higher integrated loads are computed, specially the drag force, in association with accurate vortex shedding frequencies. In the last years, our research group has made an effort to develop the in-house code for future extension to three-dimensional flows, and this research is integral part of the project. frequencies. In the last years, our research group has made an effort to develop the in-house code for future extension to three-dimensional flows, and this research is integral part of the project. **2. Theory and Numerical Method**  *2.1. Physical Modeling* 

with a free stream velocity, *U*\*, at infinity. The fluid is Newtonian with constant kinematic viscosity,

#### **2. Theory and Numerical Method** Figure 1 illustrates the smooth circular cylinder immersed in a semi-infinity fluid domain, ,

#### *2.1. Physical Modeling .* The flow is assumed to be unsteady, incompressible, and two-dimensional. The fluid domain can

Figure 1 illustrates the smooth circular cylinder immersed in a semi-infinity fluid domain, Ω, with a free stream velocity, *U*\*, at infinity. The fluid is Newtonian with constant kinematic viscosity, ν. The flow is assumed to be unsteady, incompressible, and two-dimensional. The fluid domain can be identified by a surface *S* = *S*1∪*S*2∪*S*3, being *S*<sup>1</sup> the cylinder surface, *S*<sup>2</sup> the moving wall surface and *S*<sup>3</sup> the far away boundary. The surface *S* is required to establish the boundary conditions of the physical problem. The location of the separation points of the flow for the top (open) and the bottom (gap) sides of the cylinder are defined by θ + *sep* and θ − *sep*, respectively. The blockage effect is captured by reducing the gap-to-diameter ratio, *h*\*/*d*\*. be identified by a surface *S* = *S*1*S*2*S*3, being *S*1 the cylinder surface, *S*2 the moving wall surface and *S*3 the far away boundary. The surface *S* is required to establish the boundary conditions of the physical problem. The location of the separation points of the flow for the top (open) and the bottom (gap) sides of the cylinder are defined by + *sep* and − *sep*, respectively. The blockage effect is captured by reducing the gap-to-diameter ratio, *h*\*/*d*\*. To non-dimensionalize all the quantities of the problem, it is chosen the cylinder diameter, *d*\*, as scale length. The dimensionless time is defined by *t*\**U*\*/*d*\*. In the general formulation, the symbol \* is used to identify dimensional quantities.

**Figure 1.** Schematic representation of the semi-infinity fluid domain. **Figure 1.** Schematic representation of the semi-infinity fluid domain.

The Reynolds number is defined in the following form: *U d Re* (1) To non-dimensionalize all the quantities of the problem, it is chosen the cylinder diameter, *d*\*, as scale length. The dimensionless time is defined by *t*\**U*\*/*d*\*. In the general formulation, the symbol \* is used to identify dimensional quantities.

 The Reynolds number is defined in the following form:

$$\text{Re} = \frac{\text{LT}^\* d^\*}{\nu} \tag{1}$$

temporal series of the lift force curve.

as:

According to Bearman [1], the presence of two shear layers of opposite signals is primarily responsible for vortex shedding, and the body merely modifies the mechanism by allowing feedback between the viscous wake and the shedding of circulation at the separation points (Figure 1). As consequence, the interaction between the two shear layers as a function of the free stream velocity, *U*\*, and the body diameter, *d*\*, is the key factor to define the rate at which vortical structures are cyclically detached at the rear part of the body. Thus, the dimensionless Strouhal number is defined as: \* \* \* *d St f U* (2) being the frequency of vortex shedding, *f*\*, related with the scales *d*\* and *U*\*. The Strouhal number for the smooth cylinder with no wall confinement has been experimentally measured in the literature and found close to *St* ≈ 0.19 for an upper-subcritical Reynolds number of *Re* = 1.0 × 10<sup>5</sup> [31]. In Equation (2), the quantity *f*\* is originally obtained from the

*Energies* **2020**, *13*, x FOR PEER REVIEW 5 of 24

between the viscous wake and the shedding of circulation at the separation points (Figure 1). As consequence, the interaction between the two shear layers as a function of the free stream velocity, *U*\*, and the body diameter, *d*\*, is the key factor to define the rate at which vortical structures are

$$St = f^\* \frac{d^\*}{\mathcal{U}^\*} \tag{2}$$

being the frequency of vortex shedding, *f*\*, related with the scales *d*\* and *U*\*. *2.2. Introduction of Source Singularity Elements and Nascent Vortex Blobs* 

The Strouhal number for the smooth cylinder with no wall confinement has been experimentally measured in the literature and found close to *St* ≈ 0.19 for an upper-subcritical Reynolds number of *Re* = 1.0 × 10<sup>5</sup> [31]. In Equation (2), the quantity *f*\* is originally obtained from the temporal series of the lift force curve. In the present numerical method, the smoothed cylinder and moving wall surfaces (Figure 1) are treated by an integral formulation of the potential component of the flow [32]. These surfaces are discretized and represented by flat panels with distribution of source singularity elements with constant density. Each flat panel has a center point, named pivotal point, where the impermeability

boundary condition must be satisfied. The impermeability condition is imposed on each pivotal

#### *2.2. Introduction of Source Singularity Elements and Nascent Vortex Blobs* point and it establishes equality between the normal velocity component of a fluid particle and the

In the present numerical method, the smoothed cylinder and moving wall surfaces (Figure 1) are treated by an integral formulation of the potential component of the flow [32]. These surfaces are discretized and represented by flat panels with distribution of source singularity elements with constant density. Each flat panel has a center point, named pivotal point, where the impermeability boundary condition must be satisfied. The impermeability condition is imposed on each pivotal point and it establishes equality between the normal velocity component of a fluid particle and the normal velocity component of each pivotal point. Of numerical importance is the fact that the impermeability condition and the mass conservation of the problem are guaranteed by the source elements generation in each temporal step. normal velocity component of each pivotal point. Of numerical importance is the fact that the impermeability condition and the mass conservation of the problem are guaranteed by the source elements generation in each temporal step. In addition to this, the discrete vortex method engages in to discretize spatially the vorticity field using an instantaneous vortex blobs collection, as illustrated in Figure 1. The vortex blobs are represented by Lamb vortex elements, each one presenting a distribution of vorticity, <sup>σ</sup><sup>0</sup> ς , (commonly called the cut-off function), a circulation strength, Γ, a core size, <sup>0</sup> *σ* , and a spatial position, *x* [33,34]. The no-slip boundary condition is imposed on each pivotal point and it establishes equality between the tangential velocity component of a fluid particle and the tangential

In addition to this, the discrete vortex method engages in to discretize spatially the vorticity field using an instantaneous vortex blobs collection, as illustrated in Figure 1. The vortex blobs are represented by Lamb vortex elements, each one presenting a distribution of vorticity,ςσ<sup>0</sup> , (commonly called the cut-off function), a circulation strength, Γ, a core size, σ0, and a spatial position, *x* [33,34]. The no-slip boundary condition is imposed on each pivotal point and it establishes equality between the tangential velocity component of a fluid particle and the tangential velocity component of each pivotal point. Of numerical importance is the fact that the no-slip condition and the global circulation conservation of the problem are guaranteed by the vortex blobs generation in each temporal step. The moving wall effect dispenses this boundary condition on the ground plane, being the vortex blobs only generated on the cylinder surface [21,22,25]. Figure 2 illustrates as a vortex blob is introduced into the fluid domain. In Lagrangian manner, the vorticity generated from a flat panel stay concentrated inside blob vortex core, σ0. velocity component of each pivotal point. Of numerical importance is the fact that the no-slip condition and the global circulation conservation of the problem are guaranteed by the vortex blobs generation in each temporal step. The moving wall effect dispenses this boundary condition on the ground plane, being the vortex blobs only generated on the cylinder surface [21,22,25]. Figure 2 illustrates as a vortex blob is introduced into the fluid domain. In Lagrangian manner, the vorticity generated from a flat panel stay concentrated inside blob vortex core, <sup>0</sup> *σ* . To satisfy the boundary conditions above mentioned are necessary two different kinds of system of linear algebraic equations, which are solved iteratively using the method of least squares. The simultaneous generation of sources elements and vortex blobs can also be coupled with the roughness model, when activated; details of the roughness model used in simulations of interest will be presented in Section 2.4.

**Figure 2.** Shedding of a vortex blob during one time stepping (*coi* is the pivotal point and *epsi'* defines the "smooth shedding point" location of the vortex blob). **Figure 2.** Shedding of a vortex blob during one time stepping (*co<sup>i</sup>* is the pivotal point and *eps<sup>i</sup> '* defines the "smooth shedding point" location of the vortex blob).

To satisfy the boundary conditions above mentioned are necessary two different kinds of system of linear algebraic equations, which are solved iteratively using the method of least squares. The simultaneous generation of sources elements and vortex blobs can also be coupled with the

roughness model, when activated; details of the roughness model used in simulations of interest will be presented in Section 2.4.

#### *2.3. Discrete Vortex Method with LES Modeling*

The motion of each vortex blob is governed by the vorticity transport equation, which is obtained by taking the curl of the Navier–Stokes equations [33]. Chorin [35] proposed an algorithm that splits the vorticity transport equation to separately solve the advection and diffusion problems. Alcântara Pereira et al. [36] originally presented the solution of the two problems including LES modeling into the two-dimensional discrete vortex method, such as, respectively:

$$\frac{D\overline{\omega}}{Dt} = \frac{\partial \overline{\omega}}{\partial t} + (\overline{\mathfrak{u}} \cdot \nabla)\overline{\omega} = 0 \tag{3}$$

$$\frac{\partial \overline{\boldsymbol{\omega}}}{\partial t} = \left(\frac{1}{\text{Re}} + \boldsymbol{\nu}\_t^\*\right) \nabla^2 \overline{\boldsymbol{\omega}} \tag{4}$$

where

$$\nu\_t^\* = \frac{\nu\_t}{\overline{U^\* d^\*}}\tag{5}$$

represents the local eddy viscosity coefficient and ω = ∇ × *u* defines the vorticity scalar field.

The vorticity generated from each flat panel in a time stepping is regarded "free" to undergo advection and diffusion process satisfying Equations (3) and (4), respectively. The vortex blobs transport by advection (Equation (3) is computed through the following expression:

$$\mathbf{x}\_{i}^{t+\Delta t} = \mathbf{x}\_{i}^{t} + \boldsymbol{\mu}(\mathbf{x}\_{i})^{t+a\Delta t}\Delta t \ = \mathbf{x}\_{i}^{t} + \left[ (\mathbf{K} \* \overline{\boldsymbol{\omega}})(\mathbf{x}) + \overline{\boldsymbol{\omega}}\overline{\boldsymbol{b}\_{i}^{n}}(\mathbf{x}\_{i}) + 1 \right]^{t+a\Delta t} \Delta t \tag{6}$$

where *u*(*xi*) *t*+α ∆*t* represents the velocity vector of the filtered field, and α is the temporal integrating parameter, such that 0 ≤ α ≤ 1, where α = 0 and α = 1 define an explicit and an implicit Euler scheme, respectively. In this work, is adopted an explicit Euler scheme.

In Equation (6), the velocity vector is computed at the point occupied by the *i*th vortex blob according to the Biot–Savart law (vortex–vortex interaction), panel method (vortex–panel interaction) and free stream velocity (vortex–mainstream interaction) contributions, respectively. It is remarkable that the Lagrangian manner [37] dispenses the need to explicitly treat advective derivatives into the Equation (3).

Equation (4) is solved using to the random walk method [35], where is imposed a displacement for each vortex blob in the following form [36]:

$$
\zeta\_i(t) = \sqrt{\frac{4\Delta t}{\text{Re}\_{c\_i}}} \ln \left(\frac{1}{p}\right) \left[\cos(2\pi Q)\_x + \sin(2\pi Q)\_y\right] \tag{7}
$$

being *p* and *Q* random numbers generated between 0.0 and 1.0, and *Re<sup>c</sup>* the local Reynolds number modified by the local eddy viscosity coefficient, such as [36]:

$$Re\_{\mathbf{c}\_l}(t) = \frac{\mathcal{U}^\* d^\*}{\nu + \nu\_{l\_l}(t)}\tag{8}$$

In Equation (7), *p* and *Q* define two random displacements with mean equal to zero and a variance given by twice the product of the kinematic viscosity and the time. The vorticity diffusion process through the random displacements of vortex blobs simulates the viscosity effect.

In the present formulation, the local turbulence effect is simulated during the vorticity diffusion process. Therefore, the local eddy viscosity coefficient computation is necessary to include the effect of the small scales through the concept of differences of velocity between vortex blobs [36]. The support to the turbulence modeling success is that each vortex blob needs to move with the local fluid velocity

*Energies* **2020**, *13*, 4481 function of the filtered field [38] for each vortex blob that constructs the viscous wake. Alcântara

diffusion process.

induced at each vortex blob (see *<sup>t</sup> α Δt*

in a Lagrangian manner to simulate the vorticity advection process, and the velocity induced at each vortex blob (see *u*(*xi*) *t*+α ∆*t* in Equation (6)) is calculated before the vorticity diffusion process. two-dimensions, such as: 2 3 

Pereira et al. [36] proposed an adaption to compute the second-order velocity structure function in

The average velocity differences are required to calculate the second-order velocity structure

*Energies* **2020**, *13*, x FOR PEER REVIEW 7 of 24

fluid velocity in a Lagrangian manner to simulate the vorticity advection process, and the velocity

*i u*

The average velocity differences are required to calculate the second-order velocity structure function of the filtered field [38] for each vortex blob that constructs the viscous wake. Alcântara Pereira et al. [36] proposed an adaption to compute the second-order velocity structure function in two-dimensions, such as: 2 0 1 2 1 *j j j j j k k N j F N t t k k u x u x r r* (9)

$$\overline{F}\_{2\_j} = \frac{1}{N} \sum\_{k=1}^{N} \left\| \overline{u}\_{l\_j}(\mathbf{x}\_j) - \overline{u}\_{l\_j}(\mathbf{x}\_j + \mathbf{r}\_k) \right\|\_{k}^{2} \left( \frac{\sigma\_0}{r\_k} \right)^{2/3} \tag{9}$$

*x* in Equation (6)) is calculated before the vorticity

where *j* defines the position of *j*th vortex blob, *u<sup>t</sup>* defines the velocity filtered field computed at each vortex blob, *N* characterizes a special group of vortex blobs inside a circular crown idealized around the *j*th vortex blob under analysis, and *r<sup>k</sup>* measures the distance between the *j*th and *k*th vortex blobs, the latter necessarily belonging to that special group (for more details, please see Bimbato et al. [11]). vortex blobs, the latter necessarily belonging to that special group (for more details, please see Bimbato et al. [11]). After the solution of Equation (9), the local eddy viscosity coefficient computation is obtained through the following formula [38]:

After the solution of Equation (9), the local eddy viscosity coefficient computation is obtained through the following formula [38]: 3 2 0 2 ( ) 0.105 ( ) *i k i t k t C F t* (10)

$$\nu\_{l\_l}(t) = 0.105 \text{C}\_{\text{k}}^{-3/2} \sigma\_{\text{0}} \sqrt{\overline{F}\_{2\_l}(t)} \tag{10}$$

being *C<sup>k</sup>* = 1.4 the Kolmogorov constant. necessary to stabilize the numerical solution of the problem; furthermore, it also provides a basis for

It is important comment that use of two-dimensional LES-based turbulence modeling is necessary to stabilize the numerical solution of the problem; furthermore, it also provides a basis for a future three-dimensional turbulence modeling. Bimbato et al. [21,22] validated the LES modeling used in this work. a future three-dimensional turbulence modeling. Bimbato et al. [21,22] validated the LES modeling used in this work. As already mentioned in Section 2.2, the smoothed cylinder surface is represented by source flat panels [32], being that each one also produces one vortex blob of strength Γ*i* at every time stepping.

As already mentioned in Section 2.2, the smoothed cylinder surface is represented by source flat panels [32], being that each one also produces one vortex blob of strength Γ*<sup>i</sup>* at every time stepping. Over a period of time, the viscous boundary layer develops to take on the form of a cluster of several thousand vortex blobs, as can be seen from the computer output shown in Figure 3. The separation that occurs on the cylinder surface originates an unsteady flow with the presence of Von Kármán large-scale vortices downstream the body. The proximity of a moving wall (blockage effect) will certainly interfere in any way with that vortex formation regime. Over a period of time, the viscous boundary layer develops to take on the form of a cluster of several thousand vortex blobs, as can be seen from the computer output shown in Figure 3. The separation that occurs on the cylinder surface originates an unsteady flow with the presence of Von Kármán large-scale vortices downstream the body. The proximity of a moving wall (blockage effect) will certainly interfere in any way with that vortex formation regime.

**Figure 3.** Starting flow around the smooth circular cylinder with no wall confinement. **Figure 3.** Starting flow around the smooth circular cylinder with no wall confinement.

It is remarkable that with the Lagrangian tracking of vortex blobs, one need not take into account the far away boundary conditions, i.e., at *S*3 in Figure 1. In addition, the computations are It is remarkable that with the Lagrangian tracking of vortex blobs, one need not take into account the far away boundary conditions, i.e., at *S*<sup>3</sup> in Figure 1. In addition, the computations are only concentrated on regions containing significant vorticity, which are the regions of high activities of the flow. The Lagrangian formulation of the discrete vortex method indeed dispenses a grid for the spatial discretization of the interest domain. Thus, numerical instabilities associated to high Reynolds number flows do not require special care in contrast to the Eulerian schemes. In order to take care of both advection and diffusion of the vorticity, one makes use of an advection-diffusion splitting

algorithm [35]; according to it the advection of each vortex blob is carried out independently of the diffusion (Equations (6) and (7)). In other cases, when vortex blobs migrate to the interior of a solid surface, they are reflected from their paths. independently of the diffusion (Equations (6) and (7)). In other cases, when vortex blobs migrate to the interior of a solid surface, they are reflected from their paths. *2.4. The Roughness Model* 

*Energies* **2020**, *13*, x FOR PEER REVIEW 8 of 24

only concentrated on regions containing significant vorticity, which are the regions of high activities of the flow. The Lagrangian formulation of the discrete vortex method indeed dispenses a grid for the spatial discretization of the interest domain. Thus, numerical instabilities associated to high Reynolds number flows do not require special care in contrast to the Eulerian schemes. In order to

splitting algorithm [35]; according to it the advection of each vortex blob is carried out

#### *2.4. The Roughness Model* The surface roughness effect is associated with the vortex blobs creation process, such that the

The surface roughness effect is associated with the vortex blobs creation process, such that the circulation strength, Γ*smooth*, of each nascent vortex blob is increased by the amount of ∆Γ. The value of ∆Γ is defined by the turbulent activity around a "smooth shedding point" (Figure 2) and accounting by changing the vorticity of the nascent vortex blob. Adjacent to each flat panel, used to represent the cylinder surface, is simulated an inertial effect promoted by body roughness. circulation strength, Γ*smooth*, of each nascent vortex blob is increased by the amount of ΔΓ. The value of ΔΓ is defined by the turbulent activity around a "smooth shedding point" (Figure 2) and accounting by changing the vorticity of the nascent vortex blob. Adjacent to each flat panel, used to represent the cylinder surface, is simulated an inertial effect promoted by body roughness. The key idea of the roughness model capturing changes in the vortex shedding frequency

The key idea of the roughness model capturing changes in the vortex shedding frequency (Strouhal number) is sketched in Figure 4a, where *co<sup>i</sup>* represents a pivotal point and *eps<sup>i</sup>* (*t*) defines the location of a vortex blob, after its vorticity to be changed by surface roughness effect. The numerical effect of the roughness model is to change the core size, σ0, of each nascent vortex blob, using the following equation [21]: (Strouhal number) is sketched in Figure 4a, where *coi* represents a pivotal point and *epsi*(*t*) defines the location of a vortex blob, after its vorticity to be changed by surface roughness effect. The numerical effect of the roughness model is to change the core size, <sup>0</sup> *σ* , of each nascent vortex blob, using the following equation [21]:

$$
\sigma\_{0c\_k}(t) = 1.41421 \sqrt{\frac{\Delta t}{\text{Re}} \left( 1 + \frac{\nu\_{l\_i}(t)}{\nu} \right)}\tag{11}
$$

being ∆*t* the temporal step estimated from velocity scale of the flow and *Re* the Reynolds number of the flow, as defined in Equation (1). being Δ*t* the temporal step estimated from velocity scale of the flow and *Re* the Reynolds number of the flow, as defined in Equation (1).

**Figure 4.** (**a**,**b**) Shedding of a vortex blob during one time stepping (*coi* is the pivotal point and *epsi*(*t*) defines the "rough shedding point" location of the vortex blob). **Figure 4.** (**a**,**b**) Shedding of a vortex blob during one time stepping (*co<sup>i</sup>* is the pivotal point and *eps<sup>i</sup>* (*t*) defines the "rough shedding point" location of the vortex blob).

In Equation (11), the local eddy viscosity coefficient computation newly depends on the solution of Equation (10). Therefore, the average velocity differences, required to calculate the second-order velocity structure function of the filtered field and to simulate surface roughness effect, are evaluated between the center of a semicircle with radius of *b* = 2*ε* − <sup>0</sup> *σ* and rough points, both placed near to each pivotal point, as follow [11]: In Equation (11), the local eddy viscosity coefficient computation newly depends on the solution of Equation (10). Therefore, the average velocity differences, required to calculate the second-order velocity structure function of the filtered field and to simulate surface roughness effect, are evaluated between the center of a semicircle with radius of *b* = 2ε − σ<sup>0</sup> and rough points, both placed near to each pivotal point, as follow [11]:

$$\overline{F}\_{2\_i}(t) = \frac{1}{NR} \sum\_{w=1}^{NR} \left\| \overline{\boldsymbol{u}}\_{l\_i}(\mathbf{x}\_i, t) - \overline{\boldsymbol{u}}\_{l\_W}(\mathbf{x}\_i + \mathbf{b}, t) \right\|\_{w}^2 (1 + \varepsilon) \tag{12}$$

where *u<sup>t</sup>* is the instantaneous velocity filtered field, *NR* defines twenty-one rough points adjacent to each flat panel (Figure 4b), and (1 + ε) characterizes the kinetic energy gain associated with the average roughness effect, ε = ε <sup>∗</sup>/*d* ∗ (for more details, please see Bimbato et al. [11]).

#### *2.5. Aerodynamic Loads*

As previously commented in Section 2.3, the only governing equation in the discrete vortex method is the vorticity transport equation of the filtered field. It should be remarked that the pressure is missing in the formulation, because the pressure term is eliminated when is applied the curl in the

Navier–Stokes equations. The pressure filtered field is then recuperated by taking the divergence operator to the Navier–Stokes equations. The procedure starts with the stagnation pressure definition, such as:

$$
\overrightarrow{Y}^\* = \frac{\overrightarrow{p}^\*}{\rho} + \frac{\overrightarrow{u}^2}{2}, \overrightarrow{u}^\* = \|\overrightarrow{u}^\*\|\tag{13}
$$

where *p* ∗ defines the static pressure, ρ is the fluid density and *u* ∗ represents the velocity.

Thus, a Poisson equation is derived, and the static pressure is obtained using an integral formulation, such as [39,40]:

$$H\overline{Y}\_p - \int\_{S\_1 + S\_2} \overline{Y} \nabla \overline{\Sigma}\_p \cdot \mathbf{e}\_n d\mathcal{S} \ = \iint\_{\Omega} \nabla \overline{\Sigma}\_p \cdot (\overline{u} \times \overline{\omega}) d\Omega - \frac{1}{Re} \int\_{S\_1} \left( \nabla \overline{\Sigma}\_p \times \overline{\omega} \right) \cdot \mathbf{e}\_n d\mathcal{S} \tag{14}$$

where *p* is the any point to compute pressure, *H* = 1 for points of the fluid domain, *H* = 0.5 for pivotal points, Ξ represents the fundamental solution of Laplace equation, and *e<sup>n</sup>* defines the unit vector normal to each solid boundary.

Finally, the drag and lift coefficients are calculated as, respectively [21]:

$$\mathcal{C}\_{D} = \sum\_{p=1}^{NP} \mathbf{2} (\overline{p}\_p - p\_{\infty}) \Delta S\_p \text{sen}\beta\_p = \sum\_{p=1}^{NP} \mathcal{C}\_{Pp} \Delta S\_p \text{sen}\beta\_p \tag{15}$$

$$\mathcal{C}\_{L} = -\sum\_{p=1}^{NP} \mathbf{2} \{ \overline{p}\_p - p\_{\infty} \} \Delta S\_p \cos \beta\_p = -\sum\_{p=1}^{NP} \mathcal{C}\_{P\_p} \Delta S\_p \cos \beta\_p \tag{16}$$

where *NP* represents the total number of flat panels, *p*<sup>∞</sup> represents the reference pressure at *S*<sup>3</sup> (Figure 1), ∆*S<sup>p</sup>* defines a flat panel length, and β*<sup>p</sup>* defines a flat panel angle.

#### *2.6. Computational Sequence for Solution of the Navier–Stokes Equations*

The numerical method described above is implemented to run sequentially according to the followings steps: (i) simultaneous generation of source panels and nascent vortex blobs (including the roughness model); (ii) calculation of the velocity vector at the point occupied by each vortex blob; (iii) calculation of surface pressure distribution and hence drag and lift on the cylinder; (iv) advection of the vortex blobs; (v) diffusion of the vortex blobs (including the LES modeling); (vi) reflection of the vortex blobs that migrate into the cylinder or ground plane; (vii) calculation of the velocity induced by vortex cloud at pivotal points; and (viii) advance by time ∆*t*.

#### **3. Results and Discussion**

#### *3.1. Simulation Setup*

The chosen upper-subcritical Reynolds number for all computations was *Re* = 1.0 × 10<sup>5</sup> , which allows comparisons with experimental results presented by Nishino [25], when possible. Test cases were previously performed for the smooth circular cylinder aiming to find suitable values for the following parameters [21,22]: dimensionless time, *t* = 75; dimensionless time step, ∆*t* = 0.05; number of nascent vortex blobs during each time step, *M* = 300; and initial displacement of each nascent vortex blob adjacent to the flat panel, *eps<sup>i</sup>* ' = σ0*<sup>i</sup>* = 0.001, as illustrated in Figure 2.

The dimensionless time step was estimated according to ∆*t* = *2*π*K*/*N*, being 0 < *K* < 1 and *N* = 300 [21,22] and it depends on the accuracy of the explicit Euler scheme; in this procedure *N* represents the number of flat panels adopted to discretize the cylinder surface [32]. All obtained values were computed in the interval 37.50 ≤ *t* ≤ 75.00 to compute time-averaged results. The chosen dimensionless time step was found suitable to compute aerodynamic loads with accuracy reducing the final time of the simulations.

In Figure 1, the wall plane length was fixed as *L*\*/*d*\* = 10 and discretized using 950 flat panels. A horizontal distance of *Ll*\*/*d*\* = 3 allows to identify the cylinder front stagnation point, where the reference starts from the origin of inertial frame of reference placed at (*x*\*/*d*\*; *y*\*/*d*\*) = (0.0; 0.0). That horizontal distance was previously investigated, being the length of *Ll*\*/*d*\* = 3 found suitable to capture the wall confinement (blockage) effect [21,22]. The blockage effect was captured when the gap-to-diameter ratio, *h*\*/*d*\*, was reduced from 0.80 to 0.05 (Tables 1 and 2). The relative roughness size of ε\*/*d*\* = 0.001 was chosen for the test cases, when the roughness model was activated. All the numerical results with no roughness model were identified by ε\*/*d*\* = 0.000 in Section 3.2.



**Table 2.** Summary of data for the Strouhal number and position of separation for the circular cylinder in the vicinity of a moving wall (*Re* = 1.0 × 10<sup>5</sup> ).


#### *3.2. Circular Cylinder in the Vicinity of a Moving Wall*

This section presents simultaneous measurements of integrated aerodynamic loads and surface pressure distributions for the circular cylinder, which are essentials to support all bellow discussions. The main objective is to report that the roughness model associated with blockage effect captures the full interruption of vortex shedding from the slightly rough cylinder surface placed closer to a moving wall. Thus, the results of drag force reduction, positive lift force, Strouhal number behavior, and location of the separations points of the flow around the cylinder will successfully support the analyses in a very good physical sense.

Table 1 presents experimental and numerical results of the aerodynamic force coefficients for the cylinder at upper-subcritical Reynolds number flows of *Re* = 1.0 × 10<sup>5</sup> . The experimental results are for the smooth cylinder case at different gap-to-diameter ratios, *h*\*/*d*\*, being the uncertainties in the drag, *CD*, and lift, *CL*, coefficients of ±0.016 and ±0.011, respectively, with 95% confidence [25]. The goal is to compare them to our numerical results, also presented in Table 1.

The experimental results presented in Table 1 and identified as "using end plates" were obtained for the ratios of *ye*\*/*d*\* = 0.00 and 0.40. In the experiments of Nishino [25], the length *ye*\*was defined as the distance from the edge of the cylinder to the bottom border of the end plate. In the experimental study, the bottom border of the end plate was placed below the cylinder, being the edge of the cylinder its bottom side (θ + *sep* = 270◦ in Figure 1). *Energies* **2020**, *13*, x FOR PEER REVIEW 12 of 24

With reference to Table 1, the experiments of Nishino [25] revealed that the drag force increases with the use of end plates for *h*\*/*d*\* ≥ 0.45, i.e., the flow becomes more two-dimensional. The use of end plates with *ye*\*/*d*\* = 0.40 revealed the flow closest to a two-dimensional pattern. The use of a pair of end plates, especially at high Reynolds numbers flows, was justified by Nishino [25] since the effect of the end condition of the cylinder cannot usually be large enough in practical investigations. As comparison, the numerical results of drag force for the smooth cylinder (ε\*/*d*\* = 0.000) present a very good agreement with the end plates configuration at *ye*\*/*d*\* = 0.40, being the difference of the drag coefficient between them around 10%. When the small-gap regime is identified for *h*\*/*d*\* < 0.20, the drag force significantly reduces, and this behavior occurs because of the surface roughness effect. Lei et al. [23] observed that the critical drag behavior cannot be accurately determined, since experimental and numerical investigations are carried out using discrete gap-to-diameter ratios, and the vortex shedding suppression manifests as ratio *h*\*/*d*\* is gradually reduced. In general, the numerical results show that the lift force for the rough cylinder is slightly lower as compared as the smooth cylinder; that behavior agrees with observations [2]. shear layer of opposite signal (Figure 7b). As consequence, a low pressure region is created at the cylinder upper surface (Figure 6a). The instant represented by point C in Figure 5a indicates a minimum value of the lift coefficient, where a counter-clockwise vortical structure detaches from the smooth cylinder lower surface and moves toward the near wake. That vortical structure can be visualized in Figure 7c and it creates a low pressure region at the cylinder lower surface (Figure 6a). The same counter-clockwise vortical structure grows and starts to attract the shear layer of opposite signal (Figure 7d), the latter is feeding the clockwise vortical structure causing its detachment. It is important to observe that a new upper vortical structure will be born on the upper surface and will start to grow attracting the lower shear layer; the latter also will feed the lower vortical structure causing its detachment. The complete incorporation to the near wake of the lower and upper vortical structures is revealed at instants identified by points B and D in Figure 7b,d, respectively. In Figure 6a, the instant defined by point E represents the same event previously described for the point A.

Table 2 summarizes results of the Strouhal number, St, and separation point prediction, θ + *sep* and θ − *sep* (Figure 1), for the same study cases shown in Table 1. There are no experimental data available of Strouhal number for the flow around smooth and rough cylinders near a moving wall. The experimental results for position of separation of the flow past the cylinder with no end plates [25] were included for comparison purposes. Nishino [25] did not report experimental results of the separation points prediction for the configuration of the cylinder using end plates at *ye*\*/*d*\* = 0.40. The mechanism above reported is cyclic repeating alternatively on the upper and lower sides of the smooth cylinder surface. Thus, an unsteady flow with the presence of von Kármán-type vortices takes place downstream of the cylinder (Figure 8a). That phenomenon agrees to the classical vortex shedding mechanism of the smooth cylinder with no ground effect [41]. In Figure 8a, the viscous wake downstream takes the form of ''mushroom-type" vortical structures, being that the blockage effect will destroy them far from the cylinder. For the large-gap regime, antisymmetrical

As illustration, the temporal history of the drag and lift coefficients of the smooth cylinder and other rough, both at *h*\*/*d*\* = 0.50 and 0.05, is shown in Figure 5. perturbations are captured from the near wake region and they are felt near the cylinder surface. These perturbations intrinsically relate to the Von Kármán large-scale vortices formation mode.

**Figure 5.** *Cont.*

**Figure 5.** Temporal history of the drag and lift coefficients for the circular cylinder in the vicinity of a moving wall (*Re* = 1.0 × 10<sup>5</sup> ). **Figure 5.** Temporal history of the drag and lift coefficients for the circular cylinder in the vicinity of a moving wall (*Re* = 1.0 × 10<sup>5</sup> ).

Figure 6a,b present instantaneous pressure coefficient, *CP*, distributions related to instants represented by points A–E, marked in Figure 5a,b, for the cylinder at *h*\*/*d*\* = 0.50. In Figure 6a,b, the separation angle on upper (open) side of the cylinder is identified as θ≡θ + *sep*.

In Figure 5a, the instant represented by point A characterizes a maximum value of the lift coefficient, where a clockwise vortical structure detaches from the smooth cylinder upper surface and moves toward the near wake (Figure 7a). That vortical structure grows and starts to attract the shear layer of opposite signal (Figure 7b). As consequence, a low pressure region is created at the cylinder upper surface (Figure 6a).

The instant represented by point C in Figure 5a indicates a minimum value of the lift coefficient, where a counter-clockwise vortical structure detaches from the smooth cylinder lower surface and moves toward the near wake. That vortical structure can be visualized in Figure 7c and it creates a low pressure region at the cylinder lower surface (Figure 6a). The same counter-clockwise vortical structure grows and starts to attract the shear layer of opposite signal (Figure 7d), the latter is feeding the clockwise vortical structure causing its detachment. It is important to observe that a new upper vortical structure will be born on the upper surface and will start to grow attracting the lower shear layer; the latter also will feed the lower vortical structure causing its detachment. The complete incorporation to the near wake of the lower and upper vortical structures is revealed at instants identified by points B and D in Figure 7b,d, respectively. In Figure 6a, the instant defined by point E represents the same event previously described for the point A.

The mechanism above reported is cyclic repeating alternatively on the upper and lower sides of the smooth cylinder surface. Thus, an unsteady flow with the presence of von Kármán-type vortices takes place downstream of the cylinder (Figure 8a). That phenomenon agrees to the classical vortex shedding mechanism of the smooth cylinder with no ground effect [41]. In Figure 8a, the viscous wake downstream takes the form of "mushroom-type" vortical structures, being that the blockage effect will destroy them far from the cylinder. For the large-gap regime, antisymmetrical perturbations are captured from the near wake region and they are felt near the cylinder surface. These perturbations intrinsically relate to the Von Kármán large-scale vortices formation mode.

**Figure 6.** Instantaneous pressure distributions for the circular cylinder in the vicinity of a moving wall at *h*\*/*d*\* = 0.50 (*Re* = 1.0 × 10<sup>5</sup> ). **Figure 6.** Instantaneous pressure distributions for the circular cylinder in the vicinity of a moving wall at *h*\*/*d*\* = 0.50 (*Re* = 1.0 × 10<sup>5</sup> ).

Gerrard [41] stated that a vortical structure shedding from bluff body surface continues to grow (see, e.g., the upper vortical structure in Figure 7b), being fed by circulation originated of the connected shear layer and that, when it is strong enough, draws the opposing shear layer across the near wake. In this mechanism, the approximation of the shear layers of opposite signals is able to cut off further supply of circulation to the growing vortical structure, the latter, then, sheds and moves downstream the cylinder.

In this work, the mechanism of vortical structures formation at the rear part of the cylinder with no wall confinement [41] has also been identified for the smooth cylinder at *h*\*/*d*\* = 0.50 (Figure 7a–d). Therefore, it can be concluded that the Venturi effect really redraws the lower vortical structure shed downstream. Figure 7c gives us an idea of the Venturi effect acting on a nascent vortical structure and deforming it.

The Venturi effect also contributes creating two different highest peaks for the drag coefficient curve, which is synchronized with the lift coefficient curve (Figure 5a). The explication for this interesting aerodynamic force behavior is that while the upper vortical structure finds total freedom to grow at the rear part of the cylinder until to be incorporated by the viscous wake (Figure 7b), leading to a bigger value in the drag coefficient (see approximately the drag coefficient value by projecting point C in Figure 5a), the developing of the lower vortical structure is affected by the Venturi effect (Figure 7c). This second event reflects the smaller peak in the drag curve (see approximately the drag coefficient value by projecting point D in Figure 5a). *Energies* **2020**, *13*, x FOR PEER REVIEW 15 of 24

**Figure 7.** Near wake patterns for the circular cylinder in the vicinity of a moving wall (ε\*/*d*\* = 0.000; *Re* = 1.0 × 10<sup>5</sup> ).

**Figure 7.** Near wake patterns for the circular cylinder in the vicinity of a moving wall (*ε*\*/*d*\* = 0.000;

(see, e.g., the upper vortical structure in Figure 7b), being fed by circulation originated of the connected shear layer and that, when it is strong enough, draws the opposing shear layer across the near wake. In this mechanism, the approximation of the shear layers of opposite signals is able to cut off further supply of circulation to the growing vortical structure, the latter, then, sheds and moves

*Re* = 1.0 × 10<sup>5</sup>

downstream the cylinder.

).

(**d**) *h*\*/*d*\* = 0.05 for *ε*\*/*d*\* = 0.001

**Figure 8.** Final positions of clusters of vortex blobs for the circular cylinder in the vicinity of a moving wall at *t* = 75 (*Re* = 1.0 × 10<sup>5</sup> ). **Figure 8.** Final positions of clusters of vortex blobs for the circular cylinder in the vicinity of a moving wall at *t* = 75 (*Re* = 1.0 × 10<sup>5</sup> ).

In this work, the mechanism of vortical structures formation at the rear part of the cylinder with no wall confinement [41] has also been identified for the smooth cylinder at *h*\*/*d*\* = 0.50 (Figure 7a–d). Therefore, it can be concluded that the Venturi effect really redraws the lower vortical structure shed downstream. Figure 7c gives us an idea of the Venturi effect acting on a nascent vortical structure and deforming it. The Venturi effect also contributes creating two different highest peaks for the drag coefficient curve, which is synchronized with the lift coefficient curve (Figure 5a). The explication for this interesting aerodynamic force behavior is that while the upper vortical structure finds total freedom to grow at the rear part of the cylinder until to be incorporated by the viscous wake (Figure 7b), leading to a bigger value in the drag coefficient (see approximately the drag coefficient value by projecting point C in Figure 5a), the developing of the lower vortical structure is affected by the Venturi effect (Figure 7c). This second event reflects the smaller peak in the drag curve (see approximately the drag coefficient value by projecting point D in Figure 5a). The study case for the smooth cylinder at *h*\*/*d*\* = 0.50 forecasts the occurrence of a (mean) stagnation point near the frontal part of the cylinder, where, at that position, the pressure coefficient is around 1.0 at *θ* ≈ 0.6 (Figure 6a). It is important to observe that the front stagnation point of the The study case for the smooth cylinder at *h*\*/*d*\* = 0.50 forecasts the occurrence of a (mean) stagnation point near the frontal part of the cylinder, where, at that position, the pressure coefficientis around 1.0 at <sup>θ</sup> <sup>≈</sup> 0.6◦ (Figure 6a). It is important to observe that the front stagnation point ofthe cylinder with no wall confinement is located at <sup>θ</sup>*stag* <sup>=</sup> 0.6◦ because the discretization of thecylinder surface using the panel method [32] initiates at <sup>θ</sup> <sup>=</sup> <sup>0</sup> ◦ , with the first pivotal point located at θ = 0.6◦ . The upper (open) separation angular position is identified at about θ + *sep* ≈ 79◦ (Table 2). The experimental result for the smooth cylinder with no wall confinement at *Re* = 1.0 × 10<sup>5</sup> is predictedaround <sup>θ</sup> <sup>≈</sup> <sup>82</sup>◦ [31]. The separation angle at bottom (gap) side is about θ − *sep* ≈ 91◦ (Table 2). This resultis consistent with the expected physics for the problem because the wall confinement (blockage) effectreally changes downstream the separation angle at cylinder bottom side. The experimental resultof Nishino [25] for the cylinder with no end plates at *<sup>h</sup>*\*/*d*\* <sup>=</sup> 0.40 is predicted to occur about <sup>θ</sup> − *sep* ≈ 92◦ (Table 2). It is important to observe that a lower rear pressure (Figure 6a), identified for thecylinder in ground effect, reflects a higher value for the drag coefficient (*C<sup>D</sup>* <sup>≈</sup> 1.474 in Table 1) as compared to the experimental value reported by Blevins [31], that is *C<sup>D</sup>* ≈ 1.2, with 10% uncertainty. As additional information, a numerical result available in the literature [10] for smooth cylinder withno wall confinement is about *<sup>C</sup><sup>D</sup>* <sup>≈</sup> 1.198.

$$0.6^{\circ}$$

because the discretization of the cylinder

The numerical result of the mean lift coefficient for the smooth cylinder at *h*\*/*d*\* = 0.50 is predicted around *C<sup>L</sup>* ≈ 0.104 (Table 1). The appearance of lift force pointing away from the moving wall is because of the viscosity effect, which definitively contributes to move the cylinder front stagnation point downstream and, as consequence, it is created an additional positive circulation with increasing of the lift force. This change on the front stagnation point is captured through the temporal history of points A–E in Figure 6a, in which a slight change can be identified (i.e., *C<sup>p</sup>* is not more equal to 1.0 at θ = 0.6◦ ). For the smooth cylinder with no wall confinement, the lift force oscillates around mean value of zero [10].

The numerical result of the Strouhal number for the cylinder at *h*\*/*d*\* = 0.50 seems insensitive to blockage effect, being predicted around *St* ≈ 0.204 (Table 2). This conclusion is supported by the experimental value available in the literature for the smooth cylinder with no wall confinement at upper-subcritical flow regime of *Re* = 1.0 × 10<sup>5</sup> , that is around *St* ≈ 0.19 [31], and also with 10% uncertainty. In accordance, the vortices formation mechanism [41] is not delayed, which justifies the Strouhal number value does not change. The latter conclusion is substantiated by defining the period *t<sup>E</sup>* − *t<sup>A</sup>* corresponding to the detachment of a pair of vortical structures from cylinder surface and connected to each other by a vortex sheet, which rotate in opposite directions until they be completely incorporated into the viscous wake. The period computed is of *t<sup>E</sup>* − *t<sup>A</sup>* = 61.15 − 56.15 = 5.0 for the smooth cylinder at *h*\*/*d*\* = 0.50 (Figure 6a) and of *t<sup>E</sup>* − *t<sup>A</sup>* = 4.7 for the smooth cylinder with no wall confinement [10], being the difference between them around 6.0%.

The saturation state of a typical numerical simulation at dimensionless time of *t* = 75 can be demonstrated through of the difference around 2%between the Strouhal number, obtained from the inverse of the period *t<sup>E</sup>* − *t<sup>A</sup>* ≈ 0.20 (Figure 6a), and other one of *St* ≈ 0.204, which was computed from a Fast Fourier Transformation of the lift curve between 37.5 ≤ *t* ≤ 75 (Figure 5a).

On the other hand, the temporal history of the drag and lift coefficients for the smooth cylinder at *h*\*/*d*\* = 0.05 can be seen in Figure 5c. The drag reduction is about 40.5% as compared as the smooth cylinder at *h*\*/*d*\* = 0.50 (Table 1). Figure 7e–h sketch the near wake pattern for the smooth cylinder at *h*\*/*d*\* = 0.05, being the instants defined by respective points A–D in Figure 9a. The Strouhal number reduces to *St* ≈ 0.107 (Table 2), which characterizes intermittency on the von Kármán large-scale vortex formation mode (Figure 8b).

It is of great importance for the present work that the moving wall control demonstrates the efficiency of the roughness model to reduce the drag force of the rough cylinder, when the passive control technique of vortex shedding is activated, for the chosen relative roughness size of ε\*/*d*\* = 0.001 (Figure 5b,d can be compared). The hybrid control technique is therefore able to reduce the drag force of the smooth cylinder at *h*\*/*d*\* = 0.50 around 62.1% as compared as the rough cylinder at *h*\*/*d*\* = 0.05 (from *C<sup>D</sup>* ≈ 1.474 to 0.559 in Table 1). It is interesting to comment that Alcântara Pereira et al. [10] also identified a strong reduction on the drag force when using the relative roughness size of ε\*/*d*\* = 0.007 at *h*\*/*d*\* < 0.20 (small-gap regime). In their numerical experiment, the drag force of the rougher cylinder at *h*\*/*d*\* = 0.05 reduced around 60.2% as compared as the smooth cylinder near a moving wall at *h*\*/*d*\* = 0.50. In both numerical studies, the higher drag reduction is because of the passive control technique promoting the full interruption on the von Kármán large-scale vortices formation mode.

*Energies* **2020**, *13*, x FOR PEER REVIEW 19 of 24

**Figure 9.** Instantaneous pressure distributions for the circular cylinder in the vicinity of a moving wall at *h*\*/*d*\* = 0.05; *Re* = 1.0 × 10<sup>5</sup> . **Figure 9.** Instantaneous pressure distributions for the circular cylinder in the vicinity of a moving wall at *h*\*/*d*\* = 0.05; *Re* = 1.0 × 10<sup>5</sup> .

Figure 10a–h can be accompanied in sequence to better understand both interferences of surface roughness and wall confinement, which combined, completely destroy the orderly von Kármán vortex street. Instead, the Venturi effect interferes redrawing the negative shear layer parallel to ground plane behind the rough cylinder at *h*\*/*d*\* = 0.05. Figures 8d and 10e–h sketch the near wake The rough cylinder at *h*\*/*d*\* = 0.05 configuration is strongly sensitive to interfere to the Venturi effect, which controls the smooth cylinder aerodynamics when submitted to ground effect. The surface roughness effect also participates to control the flow dynamics of the cylinder. The lift coefficient reduces of *C<sup>L</sup>* ≈ 0.531 (smooth case) for *C<sup>L</sup>* ≈ 0.313 (rough case), a reduction of 41.1% (Table 1).

pattern of the rough cylinder at *h*\*/*d*\* = 0.05, being the instants defined by points A–D in Figure 9b. Figure 5b shows that the surface roughness effect interferes on the orderly behavior of the big and small peaks previously identified in the drag curve for the smooth cylinder at *h*\*/*d*\* = 0.50 (Figure 5a). Now, in Figure 5b, there are no more single big and small peaks as the dimensionless time runs from *t* = 40 on, approximately. In Figure 5b, it is until difficult to identify the big and small peaks in the drag curve because of noise increasing originated from roughness effect. The period *t<sup>E</sup>* − *t<sup>A</sup>* is of 5.0 in Figure 6b, and of *t<sup>E</sup>* − *t<sup>A</sup>* = 4.7 for the cylinder with no wall confinement [10]. Once again, the difference between them is around 6.0%, and the suppression of vortex shedding cannot be promoted when using ε\*/*d*\* = 0.001 at *h*\*/*d*\* = 0.50 (Figure 8c).

Figure 9b presents instants randomly chosen and represented by points A–E for instantaneous pressure distributions of the rough cylinder at *h*\*/*d*\* = 0.05. The simulation also predict the appearing of a (mean) stagnation point near the cylinder frontal part, where, at that position, the pressure coefficient, *Cp*, is not more 1.0 at θ = 0.6◦ (Figure 9b). In Figure 5d, is more difficult to identify points A–E as compared as the Figure 5a.

The Strouhal number completely vanishes, *St* = 0.0, for the cylinder with roughness model at *h*\*/*d*\* = 0.05 (Table 2). Bimbato et al. [22] also identified a decrease of the Strouhal number for the smooth cylinder at *h*\*/*d*\* = 0.05, however, in their numerical experiment, the Strouhal number did not vanish (*St* ≈ 0.080). In the present simulation, the correspondent computed value was of *St* ≈ 0.107 (Table 2). An important conclusion is that for the smooth cylinder placed too close of the moving wall, i.e., at *h*\*/*d*\* = 0.05, the full interruption of von Kármán-type vortex shedding is not captured only using the active control technique by moving wall (Figure 7e–h).

The Strouhal number behavior was reported by Buresti and Lanciotti [12] for smooth and rough cylinders near a fixed ground at Reynolds numbers in the range from *Re* = 8.5 × 10<sup>4</sup> to 3.0 × 10<sup>5</sup> . The boundary layer thickness on the ground was about δ\*/*d*\* = 0.1 at cylinder location. For the flow around the smooth cylinder, within the subcritical and critical regimes (*Re* < 1.9 × 10<sup>5</sup> ), the critical gap-to-diameter ratio, *h*\*/*d*\*, was identified at 0.40, and the Strouhal number was estimated around of 0.20 for all ratios *h*\*/*d*\* greater than 0.40. The same results were obtained within the subcritical and critical regimes (in this situation, *Re* < 1.4 × 10<sup>5</sup> because of the surface roughness effect) also for the rough cylinder. Although there is a lack of experimental data of Strouhal number for the flow around the smooth and rough cylinders near a moving wall, our numerical results agree basically with an experimental result of *St* ≈ 0.20 for the cylinder near a fixed ground at large-gap regime. For the small-gap regime at *h*\*/*d*\* < 0.20, the Strouhal number decreases, and the boundary layer separation is delayed because of the combined effects of surface roughness and wall confinement (Table 2).

Figure 10a–h can be accompanied in sequence to better understand both interferences of surface roughness and wall confinement, which combined, completely destroy the orderly von Kármán vortex street. Instead, the Venturi effect interferes redrawing the negative shear layer parallel to ground plane behind the rough cylinder at *h*\*/*d*\* = 0.05. Figures 8d and 10e–h sketch the near wake pattern of the rough cylinder at *h*\*/*d*\* = 0.05, being the instants defined by points A–D in Figure 9b.

The difference of drag reduction between the study cases at *h*\*/*d*\* = 0.05 is explained because the base pressure increases for the rough cylinder (Figure 11b). It can be identified a lesser increase in the base pressure of the smooth cylinder (Figure 11a), which explains the difference about 36.3% in the drag force between them (Table 1).

Further numerical investigation will be carried out elsewhere to fully understand the relationship of the instantaneous surface pressure behavior with the mutual interaction between the two layers of opposite signals, when the complete interruption of vortex shedding is anticipated, for the rough cylinder at *h*\*/*d*\* = 0.05. Figures 8d and 10e–h give us some hints. Some small vortical structures observed in Figure 8d has been formed because of the two shear layers of opposite signals injecting vorticity at the rear part of the cylinder. As captured in our animations, the advection of vorticity of signals positive and negative creates those small vortical structures.

*Energies* **2020**, *13*, x FOR PEER REVIEW 20 of 24

**Figure 10.** Near wake patterns for the slightly rough circular cylinder in the vicinity of a moving wall **Figure 10.** Near wake patterns for the slightly rough circular cylinder in the vicinity of a moving wall (ε\*/*d*\* = 0.001; *Re* = 1.0 × 10<sup>5</sup> ).

The difference of drag reduction between the study cases at *h*\*/*d*\* = 0.05 is explained because the base pressure increases for the rough cylinder (Figure 11b). It can be identified a lesser increase in the base pressure of the smooth cylinder (Figure 11a), which explains the difference about 36.3% in

(*ε*\*/*d*\* = 0.001; *Re* = 1.0 × 10<sup>5</sup>

the drag force between them (Table 1).

).

*Energies* **2020**, *13*, x FOR PEER REVIEW 21 of 24

**Figure 11.** Time-averaged pressure distributions for the circular cylinder in the vicinity of a moving wall (*h*\*/*d*\* = 0.05; *Re* = 1.0 × 10<sup>5</sup> ). **Figure 11.** Time-averaged pressure distributions for the circular cylinder in the vicinity of a moving wall (*h*\*/*d*\* = 0.05; *Re* = 1.0 × 10<sup>5</sup> ).

#### Further numerical investigation will be carried out elsewhere to fully understand the **4. Conclusions**

relationship of the instantaneous surface pressure behavior with the mutual interaction between the two layers of opposite signals, when the complete interruption of vortex shedding is anticipated, for the rough cylinder at *h*\*/*d*\* = 0.05. Figures 8d and 10e–h give us some hints. Some small vortical structures observed in Figure 8d has been formed because of the two shear layers of opposite signals injecting vorticity at the rear part of the cylinder. As captured in our animations, the advection of vorticity of signals positive and negative creates those small vortical structures. **4. Conclusions**  The present numerical study was addressed for the control and suppression of vortex shedding from a slightly rough bluff body in the vicinity of moving wall. A discrete vortex method implemented with a hybrid control technique of vortex shedding was employed. The association of active and passive control techniques of vortex shedding to study the two-dimensional flow past the The present numerical study was addressed for the control and suppression of vortex shedding from a slightly rough bluff body in the vicinity of moving wall. A discrete vortex method implemented with a hybrid control technique of vortex shedding was employed. The association of active and passive control techniques of vortex shedding to study the two-dimensional flow past the rough circular cylinder placed too close to a moving wall (i.e., for the gap spacing of *h*\*/*d*\* = 0.05 in Figure 1) successfully captured the complete interruption of von Kármán-type vortex shedding (Figure 8d). The philosophy of our research line is to attain supercritical Reynolds number flow patterns starting from subcritical flows [11]. Achenbach [9] presented a classical *C<sup>D</sup>* × *Re* diagram (being *C<sup>D</sup>* the mean drag coefficient) for the circular cylinder, and divided it into four flow regimes, i.e., subcritical (*Re* < 2.0 × 10<sup>5</sup> − 5.0 × 10<sup>5</sup> ), critical (*Re* 2.0 × 10<sup>5</sup> − 5.0 × 10<sup>5</sup> ), supercritical (2.0 × 10<sup>5</sup> − 5.0 × 10<sup>5</sup> < *Re* ≤ 3.5 × 10<sup>6</sup> ) and transcritical (*Re* > 3.5 × 10<sup>6</sup> ), being the latter nowadays called post-critical. The numerical results reported in the present work have been substantiated by experimental data of the flow past

rough circular cylinder placed too close to a moving wall (i.e., for the gap spacing of *h*\*/*d*\* = 0.05 in Figure 1) successfully captured the complete interruption of von Kármán-type vortex shedding the smooth cylinder in the vicinity of a moving wall at upper-subcritical Reynolds number flow of *Re* = 1.0 × 10<sup>5</sup> [25]. The key findings of the present study are summarized below.


**Author Contributions:** Conceptualization, M.A.d.O., C.L.d.A., and L.A.A.P.; data curation, M.A.d.O., P.G.d.M., and C.L.d.A.; formal analysis, C.L.d.A., A.M.B., and L.A.A.P.; funding acquisition, L.A.A.P.; investigation, M.A.d.O., P.G.d.M., C.L.d.A., and A.M.B.; methodology, C.L.d.A., A.M.B., and L.A.P.P.; project administration, A.M.B. and L.A.A.P.; resources, C.L.d.A. and L.A.A.P.; software using FORTRAN, C.L.d.A. and A.M.B.; validation, M.A.d.O., P.G.d.M., and C.L.d.A.; visualization, M.A.d.O., P.G.d.M., and C.L.d.A.; supervision, L.A.A.P.; writing—original draft, L.A.A.P.; and writing—review and editing, A.M.B. and L.A.A.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported by the FAPEMIG through the research project APQ-02175-14.

**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* **Analysis of Aeroacoustic Properties of the Local Radial Point Interpolation Cumulant Lattice Boltzmann Method**

**Mohsen Gorakifard <sup>1</sup> , Clara Salueña 1,\* , Ildefonso Cuesta <sup>1</sup> and Ehsan Kian Far <sup>2</sup>**


**Abstract:** The lattice Boltzmann method (LBM) has recently been used to simulate wave propagation, one of the challenging aspects of wind turbine modeling and simulation. However, standard LB methods suffer from the instability that occurs at low viscosities and from its characteristic lattice uniformity, which results in issues of accuracy and computational efficiency following mesh refinement. The local radial point interpolation cumulant lattice Boltzmann method (LRPIC-LBM) is proposed in this paper to overcome these shortcomings. The LB equation is divided into collision and streaming steps. The collision step is modeled by the cumulant method, one of the stable LB methods at low viscosities. In addition, the streaming step, which is naturally a pure advection equation, is discretized in time and space using the Lax–Wendroff scheme and the local radial point interpolation method (RPIM), a mesh free method. We describe the propagation of planar acoustic waves, including the temporal decay of a standing plane wave and the spatial decay of a planar acoustic pulse. The analysis of these specific benchmark problems has yielded qualitative and quantitative data on acoustic dispersion and dissipation, and their deviation from analytical results demonstrates the accuracy of the method. We found that the LRPIC-LBM replicates the analytical results for different viscosities, and the errors of the fundamental acoustic properties are negligible, even for quite low resolutions. Thus, this method may constitute a useful platform for effectively predicting complex engineering problems such as wind turbine simulations, without parameter dependencies such as the number of points per wavelength *Nppw* and resolution *σ* or the detrimental effect caused by the use of coarse grids found in other accurate and stable LB models.

**Keywords:** local radial point interpolation cumulant LBM; aeroacoustics; dispersion; dissipation; wind turbine

#### **1. Introduction**

Evidence of early sailing boats on the Nile and of Persian pumps and mills from the first century B.C. shows humans have been interested in Wind Energy since ancient times [1]. In general, a wind turbine is defined as a device which converts the wind's kinetic energy into electrical energy [2,3]. It plays a key role on producing intermittent renewable energy and implementing a strategy to lower costs and reducing the reliance on fossil fuels [4,5]. Wind turbines have unique aerodynamic and aeroacoustic behavior that makes their prediction most challenging [6,7], particularly their simulation needs an enormous number of grid points or cells, and long enough time samples [8]. Researchers and centers such as the National Renewable Energy Laboratory (NREL) and the National Wind Technology Center (NWTC) have initiated multi-year programs on aeroacoustic wind turbine modeling [9] to develop efficient and appropriate computational aeroacoustic (CAA) implementations. Among particular issues specified to wind turbine problems, the propagation of sound is always a significant computational challenge [10,11]. With this

**Citation:** Gorakifard, M.; Salueña, C.; Cuesta, I.; Far, E.K. Analysis of Aeroacoustic Properties of the Local Radial Point Interpolation Cumulant Lattice Boltzmann Method. *Energies* **2021**, *14*, 1443. https://doi.org/ 10.3390/en14051443

Academic Editor: Robert Castilla

Received: 11 January 2021 Accepted: 2 March 2021 Published: 6 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 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 (https:// creativecommons.org/licenses/by/ 4.0/).

aim, different numerical approaches were developed in the field of computational aeroacoustics. Tam [12] and Wells et al. [13] proposed popular numerical schemes such as compact and non-compact optimized schemes like the high-order compact difference scheme [14] and the dispersion-relation-preserving (DRP) scheme [15]. Cheong et al. proposed grid-optimized dispersion-relation-preserving (GODRP) schemes for aeroacoustic simulations [16]. Due to the huge cost of CAA simulations, hybrid methods using two sets of equations, one for the flow and another one for the acoustic disturbance field were developed [17].

Direct aerostatic simulations are computer-intensive due to the small ratio between sound pressure and pressure variation as a whole, the spreading of acoustic fields over a large area, and the time-consuming nature of traditional methods [18]. As an illustration, the direct numerical simulation of waves using Navier–Stokes equations requires schemes of fifth-order accuracy in space and fourth order accuracy in time [19,20]. Therefore, the lattice Boltzmann method (LBM), an explicit time marching scheme [21], has been widely used as an alternative to simulate sound wave propagation due to its kinetic nature, relative simplicity of implementation and parallelization. The LBM as a mesoscopic method uses probability density functions (probability of finding particles within a certain range of velocities at a certain range of locations at a given time) to model the momentum distribution in discrete space, thereby economizing computer resources [22]. Buick et al. [23] and Dellar et al. [24] studied sound wave propagation using LBM and achieved acceptable results. Bres et al. [25] and Gorakifard et al. [26] presented the dissipation and dispersion of acoustic waves using the BGK-LBM and the cumulant LBM, respectively. Furthermore, a regularized method for the BGK-LBM [27] and the recursive and regularized LBM (LBMrrBGK) [28,29] have been developed to model wave propagation.

One of the drawbacks of the LBM is lattice uniformity, originated from symmetric lattice velocity models such as the square and cubic lattice meshes in 2D/3D simulations [30]. Lattice uniformity causes the streaming step to occur at uniform neighboring grid points. Thus, the LBM is only applied to uniform meshes, whereas issues of accuracy and computational efficiency mainly affect simulations of problems that require non-uniform meshes. For example, numerical simulations with curved and irregular boundaries, common in wind turbine simulations, encounter difficulties when fitting the grids to the boundaries or adapting complex computational domains. Grid refinement schemes or adaptive LBM can help to simulate curved and moving boundaries more accurately [31]. Wood [32] used refinement in LBM simulations to analyze wind energy and utilized adaptive LBM for moving boundaries [33]. However, these schemes are associated with higher computational costs and even additional perturbations in acoustical problems [34].

Generally, the methods commonly used along with the LBM on non-uniform meshes include three distinct categories [35]. The first is the interpolation-supplemented LBM (IS-LBM) [36,37]. This method adds an interpolation step to the collision and streaming steps of the LBM. Two major drawbacks of the IS-LBM are its inefficiency due to the need to interpolate at each time-step, and the appearance of negative particle distributions [38]. The second method is the combination of the LBM with the finite difference method (FDLBE) [39], finite volume methods (FVLBE) [40–42], or finite element methods (FELBE) [43–45] used to stabilize the computation. The third method is the Taylor-series expansion and least-squares-based lattice Boltzmann method (TLLBM) [46,47] instead of direct interpolation. These methods, the implementation of which is somewhat simple, use continuous distribution functions in physical space.

Although the above numerical schemes show robustness for complex problems, they are affected by the inherent shortcomings of using meshes in numerical methods, such as the enormous cost of generating meshes, the low level of stress accuracy in fluid–structure interaction simulations (FSI) [48], obstacles in the adaptive analysis, and limitations in simulations of physical phenomena with singular, or nearly singular behavior. Mesh-free or meshless methods have been devised to eliminate mesh-related problems [49]. The MFree method based on Liu's definition [50] is "a method that generates a system of algebraic

equations for the entire problem domain without consideration of a predefined mesh". This means that the method needs a set of scattered nodes inside the problem domain and on the boundaries of the domain called field nodes. In addition, the relationship between the nodes for the interpolation or approximation of the unknown field variables is not required [49]. Some of these well-known methods are the local point interpolation method (LPIM) [51], the local radial point interpolation method (LRPIM) [52], and the meshless local Petrov–Galerkin method (MLPG) [53,54].

The coupling of LBM and MFree methods has achieved acceptable results in some applications [55–58]; however, this approach is at an early stage of development and must be improved to address instability in low viscosities and high Re numbers. A key point in the stability and the accuracy of these methods is the collision operator, which is not remarkable in early LB methods such as the BGK model and the multi-relaxation times (MRT) model [59], both of which also violate the principle of Galilean invariance. Therefore, using a suitable collision operator to simulate complex engineering problems accurately has been recognized as a necessity. The advent of the more stable cumulant LBM [60,61] could dramatically improve MFree-LB methods and contribute to their advancement as powerful numerical tools for complex simulations.

The aim of this paper is to study the capability of the local radial point interpolation cumulant LBM (LRPIC-LBM) to simulate the propagation of planar acoustic waves, including the temporal decay of a standing plane wave and the spatial decay of a planar acoustic pulse of Gaussian shape and calculate the deviation from theoretical results, and to determine whether this method might be useful for wind turbine problems. The LB equation is deconstructed into collision and streaming parts. The collision step is performed by means of the cumulant method. The streaming step, which represents a pure advection equation, is discretized first in time using the Lax–Wendroff scheme, and then in space using the local radial point interpolation method (RPIM). Thus, this paper is arranged in the following sections; Section 2 is devoted to a brief summary of our new LB method. Section 2.1 presents the collision part referred to as cumulant collision step. Section 2.2 is devoted to the second most important part of the LBM, the streaming step. The Mfree shape function construction (RPIM shape function), time discretization (Lax–Wendroff scheme), and space discretization (LRPIM) are discussed in Sections 2.2.1–2.2.3. Section 3 reports the results of the LRPIC-LBM. Particularly, the planar standing wave and planar pulse wave propagation results are discussed in Sections 3.1 and 3.2.

#### **2. The Lattice Boltzmann Method**

The lattice Boltzmann method (LBM) is obtained from the kinetic theory of gases. In the LBM, a key point for modeling the momentum distribution in discrete space is the use of probability density functions [22]. The lattice Boltzmann equation without an external force is

$$f\_i(\mathbf{x} + e\_{\text{xi}}c\Delta t, \mathbf{y} + e\_{\text{yi}}c\Delta t, \mathbf{t} + \Delta t) - f\_i(\mathbf{x}, \mathbf{y}, \mathbf{t}) = \Omega\_{\text{i}} \tag{1}$$

where *f<sup>i</sup>* , Ω*<sup>i</sup>* , and *c* are the particle distribution function, the collision operator, and the lattice speed, respectively.

In general, the LBM consists of collision and streaming steps. In the local radial point interpolation cumulant LBM (LRPIC-LBM), the cumulant method is used for the collision part and the local radial point interpolation method (RPIM) is used for the streaming parts.

#### *2.1. Collision Step*

The cascade method has been proposed [62,63] as a way to overcome the instability problems and modeling artifacts of previous LB methods. However, it is hindered by the effects of lower order moments over higher order moments. The cumulant method, presented by Seeger to solve the Boltzmann Equation [64,65], effectively resolves these issues. Cumulants can be efficiently generated from central moments. The cumulant method for solving the LBM is briefly described in this section.

The probability density function (PDF) [22] is

$$f(\boldsymbol{\xi}, \boldsymbol{\eta}) = \sum\_{\boldsymbol{i}\boldsymbol{j}} f(\boldsymbol{\xi}\_{i\boldsymbol{\nu}} \boldsymbol{\eta}\_{\boldsymbol{j}}) \delta(\boldsymbol{\xi} - \boldsymbol{\xi}\_{\boldsymbol{i}}) \delta(\boldsymbol{\eta} - \boldsymbol{\eta}\_{\boldsymbol{i}}) \tag{2}$$

where *f* is a probability mass function (PMF) and *ξ*, *η* are discrete random variables with ranges *R<sup>ξ</sup>* = {*ξ*1, *ξ*2, . . .}, *R<sup>η</sup>* = {*η*1, *η*2, . . .} based on microscopic velocities in *x*, *y* directions.

Using the moment-generating function, *M*(Ξ, H) = ∑ *ij f*(*ξ<sup>i</sup>* , *η<sup>j</sup>* )*e* <sup>Ξ</sup>*ξie* H*η<sup>j</sup>* , the moments can be determined without any discontinuity issues as [66]

$$\mu'\_{\mathfrak{S}^m \mathfrak{h}^n} = \left. \frac{\partial^m \mathfrak{S}^n}{\partial \Xi^m \partial \mathbf{H}^n} M(\Xi, \mathbf{H}) \right|\_{\Xi = \mathbf{H} = 0} \tag{3}$$

where H, Ξ are the normalized wave numbers.

It is a good idea to shift the moment-generating function into the moving reference frame of the fluid to reduce the Galilean invariant issues in the collision process. Therefore, using the central moment-generating function, *M*ı(Ξ, H) = *e* <sup>−</sup>Ξ*u*/*c*−H*v*/*cM*(Ξ, H), the central moments can be defined as [67]

$$
\mu\_{\xi^{m}\eta^{n}} = \left. \frac{\partial^{m}\partial^{n}}{\partial\Xi^{m}\partial\mathbf{H}^{n}} \widehat{\mathcal{M}}(\Xi,\mathbf{H}) \right|\_{\Xi=\mathbf{H}=\mathbf{0}}.\tag{4}
$$

The cumulants are calculated using the logarithm form of the moment-generating function [68] as follows:

$$\varepsilon\_{\mathfrak{Z}^m \mathfrak{Y}^n} = \left. \frac{\mathfrak{Y}^m \mathfrak{Y}^n}{\mathfrak{Y}^m \mathfrak{Y}^n} \ln(\mathsf{M}(\mathsf{Z}, \mathsf{H})) \right|\_{\mathfrak{Z} = \mathsf{H} = 0}. \tag{5}$$

Each cumulant relaxes with an individual relaxation rate [69].

$$
\omega^\*\_{\xi^m \eta^n} = \omega\_{\xi^m \eta^n} + \omega\_{\xi^m \eta^n} (\mathcal{c}^{\eta}\_{\xi^m \eta^n} - \mathcal{c}\_{\xi^m \eta^n}),
\tag{6}
$$

where *c eq ξmη <sup>n</sup>* are the cumulants of the equilibrium state.

The sound speed and the kinematic viscosity are defined as *c<sup>s</sup>* = ∆*x*/( √ 3∆*t*) and *υ* = ∆*tc*<sup>2</sup> *s* (1/*ω* − 0.5), respectively.

#### *2.2. Streaming Step*

To model the streaming step, the second most important part of the LBM, a pure advection equation is normally solved from a Lagrangian approach within uniform structured meshes with CFL numbers equal to one. However, considering the Eulerian perspective for the calculation can effectively resolve this step when meshes are non-uniform and unstructured. The pure advection equation is

$$\frac{\partial f\_{\mathbf{i}}}{\partial t} + c\_{\mathbf{i},\alpha} \frac{\partial f\_{\mathbf{i}}}{\partial x\_{\alpha}} = 0,\tag{7}$$

One alternative is a semi-discrete formulation with time and space derivatives discretized separately. Thus, Equation (7) is discretized using the explicit Lax–Wendroff scheme in time, followed by the local radial point interpolation method (LRPIM) in space. In addition, it is important to approximate the unknown field functions using trial (shape) functions as an approximate solution for the partial differential equation.

2.2.1. MFree Shape Function Construction—Radial Point Interpolation Shape Functions

Radial point interpolation method (RPIM) shape functions were developed to circumvent the singularity problem arising in the point interpolation method (PIM). The RPIM interpolation augmented with polynomials is

$$f^{h}(\mathbf{x},t) = \sum\_{i=1}^{n} R\_i(\mathbf{x})a\_i(t) + \sum\_{j=1}^{m} p\_j(\mathbf{x})b\_j(t) = \mathbf{R}^T(\mathbf{x})a(t) + p^T(\mathbf{x})b(t) \tag{8}$$

where *R<sup>i</sup>* (*x*) is a radial basis function (RBF), and *p<sup>j</sup>* (*x*) is monomial in the coordinate space *x <sup>T</sup>* = [*x*, *y*]: *R <sup>T</sup>* <sup>=</sup> [*R*1(*x*) *<sup>R</sup>*2(*x*) · · · *<sup>R</sup>n*(*x*)] and *<sup>p</sup> <sup>T</sup>* <sup>=</sup> [*p*1(*x*) *<sup>p</sup>*2(*x*) · · · *<sup>p</sup>m*(*x*)]. Parameters *<sup>n</sup>* and *m* are the number of RBFs and polynomial basis functions. Variables *a<sup>i</sup>* and *b<sup>j</sup>* are time dependent unknown coefficients. It should be noted that the independent variable in RBF *Ri* (*x*) is the distance between the point of interest *x* and a node at *x<sup>i</sup>* .

Different radial basis functions (RBF), and their characteristics have been studied extensively in the meshless RPIM. In this paper, the multi-quadrics (MQ) function is used as

$$R\_i = \left(r\_i^2 + (\alpha\_c d\_c)^2\right)^q \tag{9}$$

where *α<sup>c</sup>* = 1.0, *q* = 1.03, and *d<sup>c</sup>* = 3.0.

To determine the *n* + *m* unknown coefficients in Equation (8), some specific constraint equations and the Kronecker delta function property are dictated. These constraints are

$$\sum\_{i=1}^{n} p\_j(\mathbf{x}\_i) a\_i(t) = \mathbf{P}\_m^T \mathbf{a}(t) = 0, \qquad j = 1, 2, \dots, m \tag{10}$$

where

$$\mathbf{P}\_m^T = \begin{bmatrix} 1 & 1 & \cdots & 1 \\ \mathbf{x}\_1 & \mathbf{x}\_2 & \cdots & \mathbf{x}\_n \\ \mathbf{y}\_1 & \mathbf{y}\_2 & \cdots & \mathbf{y}\_n \\ \vdots & \vdots & \ddots & \vdots \\ p\_m(\mathbf{x}\_1) & p\_m(\mathbf{x}\_2) & \cdots & p\_m(\mathbf{x}\_n) \end{bmatrix} \tag{11}$$

Thus, the approximation function can be obtained as

$$f^{\hbar}(\mathbf{x}, t) = \sum\_{i=1}^{n} \phi\_i(\mathbf{x}) f\_i(t) = \Phi(\mathbf{x}) \mathbf{F}(t) \tag{12}$$

where *F* is a vector containing the nodal values of the distribution function and **Φ** is a vector containing the first n components of the **Φ**˜ vector

$$\boldsymbol{\Phi} = \begin{bmatrix} \mathbf{R}^T \boldsymbol{\mathcal{P}}^T \end{bmatrix} \mathbf{G}^{-1} \tag{13}$$

where

$$\mathbf{G} = \begin{bmatrix} \mathbf{R}\_0 \mathbf{P}\_m \\ \mathbf{P}\_m^T \mathbf{0} \end{bmatrix}, \quad \mathbf{R}\_0 = \begin{bmatrix} R\_1(r\_1) & R\_2(r\_1) & \cdots & R\_n(r\_1) \\ R\_1(r\_2) & R\_2(r\_2) & \cdots & R\_n(r\_2) \\ \vdots & \vdots & \ddots & \vdots \\ R\_1(r\_n) & R\_2(r\_n) & \cdots & R\_n(r\_n) \end{bmatrix}, \quad r\_k = \sqrt{(\mathbf{x}\_k - \mathbf{x}\_l)^2 + \left(y\_k - y\_l\right)^2}. \tag{14}$$

#### 2.2.2. Semi-Discrete Formulation—Time Discretization

The Taylor series expansion of the particle distributions is

$$f\_i^{n+1} = f\_i^n + \Delta t \frac{\partial f\_i}{\partial t} \Big|^n + \frac{\Delta t^2}{2} \frac{\partial^2 f\_i}{\partial t^2} \Big|^n + O(\Delta t^3), \tag{15}$$

where here *n* refers to the time step. Substituting the time derivatives in terms of the *t* derivatives up to second order results in the time discretization of Equation (7) based on the Lax–Wendroff scheme,

$$f\_{i}^{n+1} = f\_{i}^{n} - \Delta t c\_{i,\alpha} \frac{\partial f\_{i}^{n}}{\partial \mathbf{x}\_{a}} + \frac{\Delta t^{2}}{2} c\_{i,\alpha} c\_{i,\beta} \frac{\partial^{2} f\_{i}^{n}}{\partial \mathbf{x}\_{a} \partial \mathbf{x}\_{\beta}}.\tag{16}$$

#### 2.2.3. Semi-Discrete Formulation—Space Discretization

The local radial point interpolation method (LRPIM) was developed to avoid the side effects of using global background cells in the global weak-form. In this method, the numerical integration is performed within the local domain consisting of a set of distributed nodes. The LRPIM is based on the RPIM shape functions with the delta function property. The main advantage of the LRPIM is the excellent interpolation stability of RBFs.

MFree local weak-form methods use the weak form of the problem obtained from the method of weighted residuals (MWR). The weighted residual statement of Equation (16) on the local domain Ω*<sup>I</sup>* of point *I* bounded by Γ*<sup>I</sup>* is posed as

$$\int\_{\Omega\_{l}} w\_{l} f\_{l}^{n+1} \, \mathrm{d}\Omega = \int\_{\Omega\_{l}} w\_{l} f\_{l}^{n} \, \mathrm{d}\Omega - \Delta t \int\_{\Omega\_{l}} w\_{l} c\_{i,\mathfrak{a}} \frac{\partial f\_{i}^{n}}{\partial \mathbf{x}\_{\mathfrak{a}}} \, \mathrm{d}\Omega + \frac{\Delta t^{2}}{2} \int\_{\Omega\_{l}} w\_{l} c\_{i,\mathfrak{a}} c\_{i,\mathfrak{b}} \frac{\partial^{2} f\_{i}^{n}}{\partial \mathbf{x}\_{\mathfrak{a}} \partial \mathbf{x}\_{\mathfrak{b}}} \, \mathrm{d}\Omega \tag{17}$$

where *w<sup>I</sup>* is the local weight function of node *I* considered as

$$w\_I(r) = \begin{cases} \frac{2}{3} - 4r^2 + 4r^3 & r \le 0.5\\ \frac{4}{3} - 4r + 4r^2 - \frac{4}{3}r^3 & 0.5 < r \le 1\\ 0 & r > 1 \end{cases} \tag{18}$$

where *r* = |*x* − *x<sup>i</sup>* |/*d* max and *d max* is the radius of the compact support. Equation (17) is applied to all nodes in the problem domain.

To work with a continuous approximate solution, it is necessary to decrease the differentiation requirements of the unknown in the weighted residual statement by employing integration by parts in Equation (17),

$$\int\_{\Omega\_{l}} w\_{l} f\_{l}^{n+1} \, \mathrm{d}\Omega = \int\_{\Omega\_{l}} w\_{l} f\_{l}^{n} \, \mathrm{d}\Omega - \int\_{\Omega\_{l}} \left( \Delta t w\_{l} c\_{i,a} \frac{\partial f\_{i}^{n}}{\partial x\_{a}} + \frac{\Delta t^{2}}{2} c\_{i,a} c\_{i,\theta} \frac{\partial w\_{l}}{\partial x\_{\beta}} \frac{\partial f\_{i}^{n}}{\partial x\_{a}} \right) \mathrm{d}\Omega + \frac{\Delta t^{2}}{2} \int\_{\Gamma\_{l}} w\_{l} c\_{i,a} c\_{i,\theta} \frac{\partial f\_{i}^{n}}{\partial x\_{a}} n\_{\beta} \, \mathrm{d}\Gamma \tag{19}$$

where Γ*<sup>I</sup>* is the boundary of the local domain Ω*<sup>I</sup>* and *n<sup>β</sup>* is the unit outward normal vector. Substitution of the approximate solution given in Equations (12) into the weak form

given by Equation (19) leads to

$$\sum\_{I=1}^{N\_I} \mathbf{M}\_{II} f\_{i,I}^{n+1} = \sum\_{I=1}^{N\_I} \left[ \mathbf{M}\_{II} + \mathbf{K}\_{i,II} \right] f\_{i,I}^{n} \tag{20}$$

where *MI J* and *Ki*,*I J* are the nodal mass and stiffness matrix, respectively, defined as

$$\mathbf{M}\_{\rm II} = \int\_{\Omega\_{\rm I}} w\_{\rm I} \boldsymbol{\Phi}\_{\rm I} \, \mathrm{d}\Omega \tag{21}$$

$$\mathbf{K}\_{i,I} = -\int\_{\Omega\_{I}} \left(\Delta t w\_{I} + \frac{\Delta t^{2}}{2} c\_{i,\theta} \frac{\partial w\_{I}}{\partial \mathbf{x}\_{\beta}}\right) c\_{i,a} \frac{\partial \Phi\_{I}}{\partial \mathbf{x}\_{a}} d\Omega + \frac{\Delta t^{2}}{2} \int\_{\Gamma\_{I}} w\_{I} c\_{i,a} \frac{\partial \Phi\_{I}}{\partial \mathbf{x}\_{a}} c\_{i,\theta} n\_{\beta} d\Gamma \tag{22}$$

Thus, the global equation system for all nodes in the entire domain is obtained as

$$\mathbf{M}f\_i^{n+1} = [\mathbf{M} + \mathbf{K}\_i]f\_i^n \tag{23}$$

where *M*, *K*, and *f<sup>i</sup>* are the global mass matrix, stiffness matrix, and particle distribution vector, respectively. This system has N equations with N unknowns which is solved separately for each direction.

To numerically evaluate the area and the curve integrals in Equations (21) and (22) the Gauss quadrature scheme is used as follows

$$\mathbf{M}\_{I\bar{I}} = \sum\_{k=1}^{n\_{\bar{\mathcal{K}}}} \tilde{w}\_k w\_I(\mathbf{x}\_k) \boldsymbol{\Phi}\_I(\mathbf{x}\_k) \left| f^{\Omega\_I} \right| \tag{24}$$

$$\mathbf{K}\_{i,II} = -\sum\_{k=1}^{n\_{\mathcal{S}}} \tilde{w}\_{k} \left( \Delta t w\_{I}(\mathbf{x}\_{k}) + \frac{\Delta t^{2}}{2} c\_{i,\emptyset} \frac{\partial w\_{I}}{\partial \mathbf{x}\_{\beta}} \bigg|\_{\mathbf{x}\_{k}} \right) c\_{i,\mathfrak{a}} \left. \frac{\partial \Phi\_{I}}{\partial \mathbf{x}\_{\mathbf{z}}} \right|\_{\mathbf{x}\_{k}} \left| I^{\Omega\_{I}} \right| + \frac{\Delta t^{2}}{2} \sum\_{k=1}^{n\_{\mathcal{S}}^{2}} \tilde{w}\_{k} w\_{I}(\mathbf{x}\_{k}) c\_{i,\mathfrak{a}} \left. \frac{\partial \Phi\_{I}}{\partial \mathbf{x}\_{\mathbf{z}}} \right|\_{\mathbf{x}\_{k}} c\_{i,\emptyset} n\_{\beta} \left| I^{\Gamma\_{I}} \right| \tag{25}$$

where *n<sup>g</sup>* and *n b <sup>g</sup>* are the total number of Gauss points in the quadrature domain and boundaries, *<sup>w</sup>*e*<sup>k</sup>* is the Gauss weight factor for Gauss point *x<sup>k</sup>* , *J* <sup>Ω</sup>*<sup>I</sup>* and *J* <sup>Γ</sup>*<sup>I</sup>* are the Jacobian matrix for the domain and boundary integrations, respectively.

#### **3. Results and Discussion**

One of the complicated phenomena that has recently received major interest from researchers using the LBM is wind turbine aeroacoustics [32,33,70,71], which can be directly simulated without additional computational cost. In this section, our aim is to demonstrate the standard analysis procedure using the local radial point interpolation cumulant lattice Boltzmann method (LRPIC-LBM) to predict acoustic properties for benchmark cases. Thus, numerical studies are conducted for the propagation of planar acoustic waves, concentrating on numerical dissipation and dispersion. Considering the total pressure equation as *p* = *p*<sup>0</sup> + *p* ′ , where *p*<sup>0</sup> is the atmospheric pressure and *p* ′ is the acoustic pressure, a lossy wave equation is [72]

$$
\left(1 + \tau\_s \frac{\partial}{\partial t}\right) \nabla^2 p' = \frac{1}{c\_s^2} \frac{\partial^2 p'}{\partial t^2} \tag{26}
$$

where

$$
\tau\_{\rm s} = \frac{1}{\rho\_0 c\_{\rm s}^2} \left(\frac{4}{3}\eta + \eta\_B\right).
$$

where *η* is the coefficient of shear viscosity, and *η<sup>B</sup>* is the coefficient of bulk viscosity. A quantitative assessment of the method's two setups, including the temporal decay of a standing plane wave in a periodic domain and the spatial decay of a propagating planar acoustic pulse of Gaussian shape for regular and irregular nodal distributions (shown in Figure 1a,b) are discussed below. It should be noted that the default nodal arrangement is a regular distribution; the base units are in the LB system.

(**b**)

**0 5 10 15 20 25 30 35 40 45 50**

**Figure 1.** Nodal arrangement for the propagation of planar acoustic waves in LRPIC-LBM simulations: (**a**) Regular nodes. (**b**) Irregular nodes.

#### *3.1. Planar Standing Wave*

**0**

**10**

As an initial case study, we performed a temporal analysis on a standing plane acoustic wave in a periodic domain. The dissipation and dispersion relations based on the temporal analysis of Equation (26) are [25]

$$c\_T = c\_s \sqrt{1 - \left(\frac{k\nu}{c\_s}\right)^2} \tag{27a}$$

$$
\omega\_T = k^2 \nu \tag{27b}
$$

where *k* is the wave number. The assumptions for this set-up are presented in Table 1. They were chosen as in reference [26] for ease of comparison.


**Table 1.** Parameters for the planar standing wave

The acoustic pressure at time *t* is [26]

$$p'(\mathbf{x}, y, t) = A \exp[-\alpha\_T t] \sin[k(\mathbf{x} - c\_T t)]\tag{28}$$

It should be noted that the results of the temporal analysis can be considered as a function of the number of points per wavelength *Nppw* = *λ*/∆*x* or the non-dimensional wave-number *k*∆*x* = 2*π*/*Nppw*.

In accordance with the concepts discussed in [26] to study the accuracy of the propagation of waves using the local radial point interpolation cumulant lattice Boltzmann method (LRPIC-LBM), various viscosities and different resolutions (i.e., the number of points per wavelength) were studied. Figure 2 shows the acoustic pressure time history for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ( ∆*x* 2 ∆*t* ) with *Nppw* = 12 points per wavelength. The analytical result is represented by a solid black line, whereas the cumulant LBM and the LRPIC-LBM are shown with a red dashed line and a blue five-pointed star, respectively. The deviations of the numerical phase speed and the temporal dissipation rate from the theoretical values are less than one percent for the BGK LBM [25] and the cumulant LBM [26] for resolutions lower than 12 points per wavelength. However, the LRPIC-LBM exhibits even better behavior in predicting the acoustic pressure of the analytical values.

**Figure 2.** Acoustic pressure ( *kg* ∆*x*∆*t* 2 ) versus time step (∆*t*) for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ( ∆*x* 2 ∆*t* ) with *Nppw* = 12 points per wavelength.

The acoustic pressure time history for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>4</sup> ( ∆*x* 2 ∆*t* ) with *Nppw* = 12 points per wavelength is presented in Figure 3. It shows the comparison between the analytical solution, the cumulant LBM and the LRPIC-LBM numerical solution at low viscosities. One of the most problematic issues in the BGK LBM is the low viscosity limit, which makes the solution unstable. However, the cumulant LBM, with phase speed and temporal dissipation rate errors of less than 1 percent, as in reference [26], did not present difficulties at low viscosity. Although both approaches were successful, the new method exhibited better performance, closely following the theoretical result for the same viscosity value.

**Figure 3.** Acoustic pressure ( *kg* ∆*x*∆*t* 2 ) versus time step (∆*t*) for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>4</sup> ( ∆*x* 2 ∆*t* ) with *Nppw* = 12 points per wavelength.

An important characteristic highlighted by some researchers [25,26] is that these numerical deviations are only a function of *Nppw* and are independent of other parameters such as frequency and viscosity. They found that the errors of the BGK [25] and the cumulant LBM [26] are about 7 percent for *Nppw* = 4. However, the acoustic pressure time history for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ( ∆*x* 2 ∆*t* ) with *Nppw* = 4 points per wavelength illustrated in Figure 4 for the analytical solution, the cumulant LBM and the LRPIC-LBM reveals that the deviation is less than 2 percent for the LRPIC-LBM, with ∆*t* = 0.25. Thus, the LRPIC-LBM is much more successful in predicting theoretical results even with a low number of points per wavelength.

**Figure 4.** Acoustic pressure ( *kg* ∆*x*∆*t* 2 ) versus time step (∆*t*) for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ( ∆*x* 2 ∆*t* ) with *Nppw* = 4 points per wavelength with ∆*t* = 0.25.

The choice of the time step size has an influence on the accuracy and stability of the solution. A smaller time step leads to more accurate results, especially in a hyperbolic partial differential Equation (PDE). To minimize the phase speed and dissipation rate errors, the time step was reduced. Figure 5 shows the acoustic pressure time history for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ( ∆*x* 2 ∆*t* ) with *Nppw* = 4 points per wavelength and ∆*t* = 0.1. The LRPI-CLBM replicates the analytical results with negligible errors. Therefore, this method makes it possible to predict wave motion more accurately with no dependency on the number of points per wavelength *Nppw*.

**Figure 5.** Acoustic pressure ( *kg* ∆*x*∆*t* 2 ) versus time step (∆*t*) for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ( ∆*x* 2 ∆*t* ) with *Nppw* = 4 points per wavelength with ∆*t* = 0.1.

Although the results of the propagation of acoustic waves for regular nodal distributions were good, it is important to adequately estimate accuracy when considering irregular nodes (Figure 1b). The acoustic pressure time history for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ( ∆*x* 2 ∆*t* ) with *Nppw* = 12 points per wavelength is presented in Figure 6. The results show that the local radial point interpolation cumulant lattice Boltzmann method (LRPIC-LBM) with irregular nodal distributions again very closely reproduced the analytical acoustic pressure.

**Figure 6.** Acoustic pressure ( *kg* ∆*x*∆*t* 2 ) versus time step (∆*t*) for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ( ∆*x* 2 ∆*t* ) with *Nppw* = 12 points per wavelength with irregular nodal distributions.

#### *3.2. Planar Pulse Wave*

Next, we studied a planar pulse wave by replacing the plane wave with a Gaussian shape planar pulse, initially located at the center of the domain. The dissipation and dispersion relationships derived from the spatial analysis of Equation (26) are [25]

$$c\_S = \sqrt{2}c\_s \sqrt{\frac{1 + \left(\omega \tau\_s\right)^2}{\sqrt{1 + \left(\omega \tau\_s\right)^2} + 1}}\tag{29a}$$

$$\alpha\_S = \frac{\omega}{\sqrt{2}c\_s} \sqrt{\frac{\sqrt{1 + (\omega \tau\_s)^2} - 1}{1 + (\omega \tau\_s)^2}} \tag{29b}$$

The variables and parameters for this case are presented in Table 2. A planar pulse emits from the origin throughout the domain, where periodic boundary conditions are imposed.



To assess the accuracy of the LRPIC-LBM in simulating the pulse wave emission, we proceeded as in reference [26]. Different viscosities and resolutions (which are related to the choice of *σ*) were studied. Figure 7 depicts the acoustic pressure time history for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ( ∆*x* 2 ∆*t* ) with *σ* = 0.1. The cumulant LBM results (dashed red line) are compared to the LRPIC-LBM (solid black line). As stated in [25,26] the intensity loss at any location is proportional to the distance of propagation, regardless of the precise location. Thus, the data were extracted from the center of the domain, and at 5, 11, and 17 nodes apart from the center. For resolutions up to *σ* = 0.1, deviations between the cumulant LBM and the LRPIC-LBM results were less than 1 percent.

**Figure 7.** Acoustic pressure ( *kg* ∆*x*∆*t* 2 ) versus time step (∆*t*) for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ( ∆*x* 2 ∆*t* ) with *σ* = 0.1: the cumulant LBM (dashed red line), LRPIC-LBM (solid black line).

The acoustic pressure time history for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>4</sup> ( ∆*x* 2 ∆*t* ) with *σ* = 0.1 is shown in Figure 8. It presents the comparison between the cumulant LBM and the LRPIC-LBM at low viscosity. While the standard BGK LBM creates instabilities and noisy results, less than 1 percent deviation was found between the cumulant LBM and this method, with stable behavior at low viscosities. In addition, as in the case of the temporal results shown in Figure 3, reducing the viscosity in the considered range did not substantially impact the results.

Like the parameter *Nppw* introduced in the temporal analysis of the first set up, *σ* is another effective parameter used in spatial analysis which affects the accuracy of the LBM results.The acoustic pressure time history is depicted in Figure 9 for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ( ∆*x* 2 ∆*t* ), with *σ* = 0.06. It shows that the deviations between the cumulant LBM and the LRPIC-LBM results increase after the reduction of *σ*. The wiggling observed in the results of the cumulant LBM is due to the reduction of the number of nodes inside the pulse, which brings the accuracy of the results into question. However, the LRPIC-LBM graphs are smooth and unaffected by the lesser number of nodes. To better assess the accuracy of the LRPIC-LBM compared to the cumulant LBM, the data shown in Figure 9 were compared with the analytical results. The Fourier transform of the pressure time history of the passing wave yields the Fourier coefficient of pressure "*p* ′ (*x*, *ω*) = "*p* ′ (*x*, 0) exp[−*αSx*] exp<sup>î</sup> *iω x cS* ó which is the solution to Equation (26) [25]. Figure 10 shows the ratio "*p* ′ (*x*, *ω*)/"*p* ′ (*x*, 0) as a function of angular frequency for the analytical solution. This figure shows that the LRPIC-LBM is more successful than the cumulant LBM at predicting theoretical results with *σ* = 0.06. Thus, this method can model wave propagation more accurately even at smaller resolutions.

**Figure 8.** Acoustic pressure ( *kg* ∆*x*∆*t* 2 ) versus time step (∆*t*) for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>4</sup> ( ∆*x* 2 ∆*t* ) with *σ* = 0.1: the cumulant LBM (dashed red line), LRPIC-LBM (solid black line).

**Figure 9.** Acoustic pressure( *kg* ∆*x*∆*t* 2 ) versus time step (∆*t*) for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ( ∆*x* 2 ∆*t* ) with *σ* = 0.06: the cumulant LBM (dashed red line), LRPIC-LBM (solid black line).

Figure 11 illustrates the acoustic pressure time history for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ( ∆*x* 2 ∆*t* ) with *σ* = 0.1 considering irregular nodes (Figure 1b). The cumulant LBM and the LRPIC-LBM results are represented by a dashed blue line and a solid black line, respectively. The data are extracted as in Figure 7. The results show that the LRPIC-LBM with irregular nodal distributions reproduced the behavior observed previously.

**Figure 10.** The ratio "*p* ′ (*x*, *ω*)/"*p* ′ (*x*, 0) as a function of angular frequency ( <sup>1</sup> ∆*t* ) for the analytical solution, the cumulant LB and the LRPIC-LB methods.

**Figure 11.** Acoustic pressure ( *kg* ∆*x*∆*t* 2 ) versus time step (∆*t*) for *<sup>ν</sup>* <sup>=</sup> 1.0 <sup>×</sup> <sup>10</sup>−<sup>2</sup> ( ∆*x* 2 ∆*t* ) with *σ* = 0.1 considering irregular nodal distributions: the cumulant LBM (dashed blue line), LRPIC-LBM (solid black line).

#### **4. Conclusions**

This paper presents a numerical study of the propagation of acoustic waves, one of the challenging issues occurring in wind turbine simulations, that includes the temporal decay of a standing plane wave and the spatial decay of a planar acoustic pulse, using the local radial point interpolation cumulant lattice Boltzmann method (LRPIC-LBM). The LB equation is divided into collision, modeled by the cumulant method, and streaming, discretized in time and space, using the Lax–Wendroff scheme and the local radial point interpolation method (RPIM). The LRPIC-LBM results were compared with the cumulant LBM and the analytical solutions. Both methods showed a similar acoustic pressure time history, and the deviations for the phase speed and the dissipation rate were minor for high number of points per wavelength *Nppw* and resolution *σ*. In addition, they showed that reduced viscosity does not affect the stability of either LB method due to the intrinsic characteristics of the cumulant method. Unlike LB methods such as the BGK LBM and the cumulant LBM, the time history for the acoustic pressure and the phase speed and dissipation rate predicted by LRPIC-LBM showed considerably smaller errors for low *Nppw* and *σ* due to the construction of the meshless method itself. Moreover, the LRPIC-LBM with irregular nodal distributions reproduces the same propagation of acoustic waves obtained with regular nodal distributions.

In summary, the freedom to scatter nodes based on problem conditions and the occurrence of sharp gradients plus the accuracy obtained from the RPIM along with the good stability and simplicity achieved by the cumulant LBM may provide an adequate platform with which to model wind turbine problems.

**Author Contributions:** Conceptualization: M.G., C.S., I.C. and E.K.F.; methodology, M.G., C.S., I.C. and E.K.F.; coding, M.G., C.S., I.C. and E.K.F.; validation, M.G., C.S., I.C. and E.K.F.; formal analysis, M.G., C.S., I.C. and E.K.F.; investigation, M.G., C.S., I.C. and E.K.F.; resources, M.G., C.S., I.C. and E.K.F.; data curation, M.G., C.S., I.C. and E.K.F.; writing—original draft preparation, M.G.; writing review and editing, C.S., I.C. and E.K.F.; visualization, M.G., C.S., I.C. and E.K.F.; supervision, C.S. and I.C.; project administration, C.S. and I.C.; funding acquisition, C.S. and I.C. All authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work has been funded by the Spanish Ministry of Economy, Industry and Competitiveness Research National Agency (under project DPI2016-75791-C2-1-P), by FEDER funds and by Generalitat de Catalunya - AGAUR (under project 2017 SGR 01234).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


## *Article* **Numeric Simulation of Acoustic-Logging of Cave Formations**

#### **Fanghui Xu and Zhuwen Wang \***

College of Geoexploration Science and Technology, Jilin University, Changchun 130026, China; xufh19@mails.jlu.edu.cn

**\*** Correspondence: wangzw@jlu.edu.cn; Tel.: +86-135-7894-7784

Received: 1 July 2020; Accepted: 22 July 2020; Published: 31 July 2020

**Abstract:** The finite difference (FD) method of monopole source is used to simulate the response of full-wave acoustic-logging in cave formations. The effect of the cave in the formation of borehole full-waves was studied. The results show that the radius of cave is not only linearly related to the first arrival of the compressional wave (P-wave), but also to the energy of the shear wave (S-wave). The converted S (S–S wave) and P-waves (S–P wave) are formed when the S-wave encounters the cave. If the source distance is small, the S–S and S–P waves are not separated, and the attenuation of the S-wave is not large, due to superposition of the converted waves. The S–P wave has been separated from the S-wave when the source distance is large, so the attenuation of the S-wave increases. The amplitude of the P and S–waves changes most when the distance of the cave to the borehole wall reaches a certain value; this value is related to the excitation frequency. The amplitude of the Stoneley wave (ST wave) varies directly with the radius of cave. If the radius of the cave is large, the energy of ST wave is weak. The scattered wave is determined by the radius and position of the cave. The investigation depth of a monopole source is limited. When the distance of the cave to the borehole wall exceeds the maximum investigation depth, the borehole acoustic wave is little affected by the cave. In actual logging, the development of the cave can be evaluated by using the first arrival of the P-wave and the energy of the S and ST waves.

**Keywords:** cave formation; P-waves; S-waves; Stoneley wave; scattered wave

#### **1. Introduction**

Caves are widely distributed in carbonate and igneous reservoirs, causing heterogeneity and making it difficult to evaluate the reservoirs by conventional logging. Based on the observation of cores, the reservoir spaces of igneous and carbonate formations can be divided into primary pores and secondary pores, according to their genesis [1]. Secondary pores contain caves. These can be divided into matrix-dissolution pores, intragranular-dissolution pores and apricot-kernel dissolution pores in igneous formations [2]. According to the genesis, carbonate formations have porous caves, intergranular caves and fractured caves. On the basis of the size of the cave, they may fall into middle and small-type caves and large caves. Generally, caves with a diameter greater than 10 mm can be regarded as large caves [3]. Because dissolution pores and caves often transform primary pores, making them better connected and improving the reservoir property of the rock, they are important reservoir spaces [4]. In fact, the caves in the formation can be defined as voids, washouts or large pores of secondary porosity in a more general. The existence of caves greatly changes the speed and energy of acoustic-logging [5]. Therefore, it is very important to grasp the effect of caves on borehole acoustic wave propagation, which plays an important role in the detection and evaluation of caves. Currently, the effect of fractured and gas-bearing reservoirs on borehole acoustics has been investigated [6–8]. However, there are few studies on the acoustic-logging response of caves. Heterogeneous bodies such

as fractures and caves, scatter seismic-wave propagation [9]. Elastic-wave scattering in heterogeneous media can be divided into "velocity" and "wave impedance" types [10]. Due to the fact that the scattered wave of a borehole dipole source is mainly a SH shear wave [11–13], the dipole source has been used to study the scattering effect in heterogeneous formation [14]. The relationship between S-wave wavelength and the size of cave can be obtained by dipole acoustic-logging [15]. According to the different filling materials, caves can be divided into four types: unfilled, sand-mud-filled, breccia filled, and calcite-filled [16]. The filling degree of a cave can slow P and S-waves. Studies have shown that different filling degrees of caves has a good corresponding relationship with acoustic simulation curves, and the results have certain application effect [17]. In cave formation, the slowness of the P-wave increases and the density of well-logging decreases. The resistivity is generally low in conventional logging [18]. Acoustic-logging and electrical-imaging logging can evaluate caves comprehensively; ST wave energy of acoustic-logging reduces due to caves [3].

The above studies show that it is feasible to use acoustic-logging to evaluate near-borehole caves. Based on actual logging data, the results have qualitatively shown the effects of cave on conventional logging, such as P-wave slowness. The response of component waves to caves has been studied. There is no detailed regularity of the size and position of cave on borehole acoustic waves. Monopole is a common source in acoustic-logging [19,20], and the component waves with different frequencies formed by monopole contain the information of specific formations. It is difficult to analyze scattered-wave fields of monopole acoustic-logging because of their complexity. Numeric simulation in cave reservoirs by monopole sources has not been reported. In this study, the finite difference (FD) method of the elastic wave equation is used to simulate a monopole full-wave acoustic-logging response of a fluid-filled cave in a homogeneous formation [21,22]. On this basis, the results of a numeric simulation are used to analyze actual logging data [23].

In the process of simulating a wave field inside and outside of a borehole by FD, matrix calculations are mainly carried out. It is necessary to divide small grids to meet the calculation requirements of cave, which greatly increases the time costs of the numeric simulation. Therefore, computer-unified device architecture (CUDA) parallel-computing was used to improve the calculation efficiency [24–27]. The essence of CUDA is to change the computing mode from central processing unit (CPU) serial to graphics-processing unit (GPU) parallel.

#### **2. Theory and Methods**

#### *2.1. FD Method for Elastic Wave Equation*

In order to study the issue more simply, we use the elastic wave equation in the 2D-cylindrical coordinate system. A 2D-model simulation has advantages over a 3D-model simulation in terms of calculation times. A 2D calculation of a cave can also show the main characteristics of the problem studied [15]. The staggered-grid finite difference method is used to simulate the propagation of elastic waves in our simulation [6]. The wave equation of velocity stress in isotropic media in cylindrical coordinates is as follows:

$$\left[\pi^{n+1}\_{r(i+\frac{1}{2},j+\frac{1}{2})}-\pi^{n}\_{r(i+\frac{1}{2},j+\frac{1}{2})}\right]\wedge\Lambda t = \left[\lambda\_{(i+\frac{1}{2},j+\frac{1}{2})}+2\mu\_{i+\frac{1}{2},j+\frac{1}{2}}\right]\left[L^{1}\_{\mathbb{R}}\widetilde{\nu}^{n+\frac{1}{2}}\_{r(i,j+\frac{1}{2})}\right] + \lambda\_{(i+\frac{1}{2},j+\frac{1}{2})}\left[L^{1}\_{\mathbb{Z}}\widetilde{\nu}^{n+\frac{1}{2}}\_{z(i+\frac{1}{2},j)} + \sigma\_{\mathbb{R}}\widetilde{\nu}^{n+\frac{1}{2}}\_{r(i,j+\frac{1}{2})}\right] \tag{1}$$

$$
\begin{bmatrix}
\tau^{n+1}\_{\varpi(i+\frac{1}{2},j+\frac{1}{2})} - \tau^{n}\_{zz(i+\frac{1}{2},j+\frac{1}{2})}
\end{bmatrix} / \Delta t = \lambda\_{(i+\frac{1}{2},j+\frac{1}{2})} \begin{bmatrix}
\iota^{1}\_{\varpi}\nu^{n+\frac{1}{2}}\_{r(i,j+\frac{1}{2})} + \frac{\sigma\_{r}\sigma^{n+\frac{1}{2}}}{r\_{i+\frac{1}{2}}}
\end{bmatrix} + \begin{bmatrix}
\lambda\_{(i+\frac{1}{2},j+\frac{1}{2})} + 2\mu\_{i+\frac{1}{2},j+\frac{1}{2}}
\end{bmatrix} \begin{bmatrix}
\iota^{1}\_{Z}\nu^{n+\frac{1}{2}}\_{z(i+\frac{1}{2},j)}
\end{bmatrix} \tag{2}
$$

$$\left[\boldsymbol{\pi}\_{0\varnothing(i+\frac{1}{2}j+\frac{1}{2})}^{n+1} - \boldsymbol{\pi}\_{0\varnothing(i+\frac{1}{2}j+\frac{1}{2})}^{n}\right] \wedge \boldsymbol{\Lambda} = \lambda\_{(i+\frac{1}{2}j+\frac{1}{2})} \left[\boldsymbol{L}\_{\boldsymbol{R}}^{1}\boldsymbol{\nu}\_{r(i+\frac{1}{2})}^{n+\frac{1}{2}} + \boldsymbol{L}\_{\boldsymbol{Z}}^{1}\boldsymbol{\nu}\_{r(i+\frac{1}{2}j)}^{n+\frac{1}{2}}\right] + \left[\lambda\_{(i+\frac{1}{2}j+\frac{1}{2})} + 2\mu\_{i+\frac{1}{2},j+\frac{1}{2}}\right] \left[\boldsymbol{\sigma}\_{r}\boldsymbol{\nu}\_{r(i+\frac{1}{2})}^{n+\frac{1}{2}} / r\_{i+\frac{1}{2}}\right] \tag{3}$$

$$\left[\boldsymbol{\tau}\_{rz(i,j)}^{n+1} - \boldsymbol{\tau}\_{rz(i,j)}^{n}\right] / \Delta t = \mu\_{(i-\frac{1}{2}, j-\frac{1}{2})}^{H} \left[L\_{\mathbb{R}}^{1} \boldsymbol{\upsilon}\_{z(i-\frac{1}{2}, j)}^{n+\frac{1}{2}} + L\_{Z}^{1} \boldsymbol{\upsilon}\_{r(i, j-\frac{1}{2})}^{n+\frac{1}{2}}\right] \tag{4}$$

$$\begin{split} \left[ \sigma\_{i} \rho\_{(i-\frac{1}{2},j+\frac{1}{2})} \middle| \upsilon^{n+\frac{1}{2}}\_{r(i,j+\frac{1}{2})} - \upsilon^{n-\frac{1}{2}}\_{r(i,j+\frac{1}{2})} \right] / \Delta t &=& L^{1}\_{\mathbb{R}} \tau^{n}\_{rr(i-\frac{1}{2},j+\frac{1}{2})} + L^{1}\_{\mathbb{Z}} \tau^{n}\_{rz(i,j)} \\ &+ \left[ \sigma\_{i} \pi^{n}\_{rr(i-\frac{1}{2},j+\frac{1}{2})} - \sigma\_{i} \pi^{n}\_{\theta\theta(i-\frac{1}{2},j+\frac{1}{2})} \right] / r\_{i} \end{split} \tag{5}$$

$$
\sigma\_l \boldsymbol{\sigma}\_{(i+\frac{1}{2},j-\frac{1}{2})} \Big| \boldsymbol{v}\_{z(i+\frac{1}{2},j)}^{n+\frac{1}{2}} - \boldsymbol{v}\_{z(i+\frac{1}{2},j)}^{n-\frac{1}{2}} \Big| / \Delta t = \boldsymbol{L}\_{\mathbb{R}}^{1} \boldsymbol{\tau}\_{z\mathbf{z}(i,j)}^{n} + \boldsymbol{L}\_{\mathbb{Z}}^{1} \boldsymbol{\tau}\_{z\mathbf{z}(i+\frac{1}{2},j-\frac{1}{2})}^{n} + \sigma\_l \boldsymbol{\tau}\_{r\mathbf{z}(i,j)}^{n} / r\_{i+\frac{1}{2}} \tag{6}
$$

In the equations mentioned above, the subscripts *i* and *j* in brackets refer to spatial index of nodes in *x* and *z* directions, respectively, the superscript *n* represents the time discrete index, τ*rr*, τθθ and τ*zz* refer to normal stress in *x*, θ and *z* directions, τ*rz* is shear stress, ∆*t* is the time step, *L* 1 *R* and *L* 1 *Z* are second-order difference operators of space, and the parameters ρ and µ represent the density and shear moduli of the medium. On the medium interface, the arithmetic mean value of ρ and the harmonic mean value of µ are employed; on the solid–liquid interface, µ is set to 0 [6,7]. σ*<sup>i</sup> f* (*i*,*j*) and σ*<sup>j</sup> f* (*i*,*j*) are arithmetic average operators and defined as:

$$\sigma\_i f\_{(i,j)} = 0.5 \left[ f\_{(i+1,j)} + f\_{(i,j)} \right] \tag{7}$$

$$
\sigma\_{\dot{\jmath}} f\_{(i,\dot{\jmath})} = 0.5 \left[ f\_{(i,\dot{\jmath})} + f\_{(i,\dot{\jmath}+1)} \right] \tag{8}
$$

µ *H* (*i*,*j*) is the harmonic mean and defined as

$$\mu\_{(i,j)}^H = 4/\left[1/\mu\_{(i+\frac{1}{2},j+\frac{1}{2})} + 1/\mu\_{(i+\frac{1}{2},j-\frac{1}{2})} + 1/\mu\_{(i-\frac{1}{2},j+\frac{1}{2})} + 1/\mu\_{(i-\frac{1}{2},j-\frac{1}{2})}\right] \tag{9}$$

To reduce the numeric dispersion, the space grid needs to satisfy the following condition:

$$
\Delta \mathbf{r} \le v\_{\min} / 10 f\_{\max} \tag{10}
$$

where *vmin* is the minimum velocity of numeric model, and *fmax* is the maximum frequency of the source [6,7]. The stability condition of the second-order FD is

$$v\_{\text{max}} \Delta t \left[ \frac{1}{\Delta r^2} + \frac{1}{\Delta z^2} \right]^{\frac{1}{2}} < 1 \tag{11}$$

where *vmax* is the maximum velocity of the numeric model [6,7]. The monopole source function is

$$s(t) = \begin{cases} \frac{1}{2} \left[ 1 + \cos\frac{2\pi}{T\_c} \left( t - \frac{T\_c}{2} \right) \right] \cos 2\pi f\_0 \left( t - \frac{T\_0}{2} \right), 0 \le t \le T\_c\\ 0, t > T\_c \end{cases} \tag{12}$$

where *f*<sup>0</sup> denotes the center frequency of the source, *t* represents the time, and *T<sup>c</sup>* is the acoustic source pulse width, *T<sup>c</sup>* = 2/ *f*0.

#### *2.2. The Principle of CUDA Parallel-Computing*

The CUDA parallel-computing model of FD method of elastic wave equation is established; the algorithm implementation flow is shown in Figure 1. During the initialization of the calculation model, the CPU first allocates storage space in the GPU memory for the velocities, stresses and medium parameters. The kernel functions of stress, velocity fields and source are then designed. The GPU is controlled by the CPU to complete the step-by-step calculations of velocity and stress fields. It was found that for the homogeneous medium, the calculation time of the conventional program was 10,571.5 s; the calculation time of the parallel program was only 704.7 s through the simulation experiment. CUDA parallel-computing is used as a computing method in this study, so it is not introduced in detail here.

() = {

= 2/<sup>0</sup>

1 2

[1 + cos

2 ( − 2

)] 2<sup>0</sup> ( −

>

0,

<sup>0</sup>

0 2 ),

0 ≤ ≤

**Figure 1.** Graphics-processing unit (GPU) parallel-computing progress based on computer unified device architecture (CUDA) platform.

#### *2.3. Borehole Model and Non-Splitting Perfect-Matched Layer Absorbing Boundary Condition*

r\_z ∆ = ∆ = 2 mm ∆ = 0.2 μs Figure 2 shows a borehole model surrounded by an elastic formation. The cave was filled with water and the formation was homogeneous. Perfect-matched layers (PML) were added at the top, bottom and right side of the model on the r\_z plane. The model dimension was 4 m × 2 m, the grid step was ∆*r* = ∆*z* = 2 mm, and the time step was ∆*t* = 0.2 µs. A non-splitting perfect-matched layer (NPML) approach [28] was used to mitigate the ghost reflection at the artificial boundary of the model.

The elastic parameters of each medium are shown in Table 1. Alteration usually occurs in carbonate and igneous reservoirs. The rock density of this kind of reservoir is generally between 2.5–2.8 g/cm<sup>3</sup> . In this study, 2.65 g/cm<sup>3</sup> was selected as the density of the formation [20]. The P and S-wave velocity were selected according to the principle that the P-wave velocity is generally 1.6–1.9 times that of the S-wave velocity [20]. The selection of these physical parameters is similar in numeric simulations of fractures and cave formations [7,15,29]. If the physical parameters of calcite and dolomite are selected in the numeric simulation, it will be of more practical significance.

**Table 1.** Physical parameters of different media.


In order to verify the correctness of the results of the FD method, the propagation of the ST wave in homogeneous elastic formation was simulated by using the FD and the RAI method [29]. The formation parameters are shown in Table 1. The center frequency of the source was 2 kHz. Figure 3 shows the comparison between the waveform calculated by the FD and RAI method. The waveforms of the two methods matched well and it showed that the FD method was effective. In addition, there was a little difference in waveform amplitude between the two methods, which was caused by discrete Fourier transform (DFT) in RAI method and finite difference dispersion in FD method. The maximum error of the two methods was 1.2%.

g/cm<sup>3</sup>

–

**Figure 2.** Schematic diagram of a cave formation in a borehole. The radius of the cave is *R*, the distance of the cave to the borehole wall is *L*, the radius of the borehole is *a*, *a* = 10 cm, and the source and receiver are placed in the borehole. The source is located under the cave; the vertical distance from the cave is 2 m. Perfect-matched layer (PML) is an absorbing boundary region. *r* is the radial coordinate and *z* is the axial coordinate of the borehole.

–

–

**Figure 3.** Stoneley (ST) waveforms propagating in a borehole surrounded by a homogeneous formation. The solid blue curve is a waveform obtained by the FD method; the dashed black curve is a waveform obtained by real-axis integration (RAI) method. The solid black curve indicates error.

–

#### **3. Numeric Simulation Results**

#### *3.1. Borehole Acoustic Wave and Wave Field Snapshot of Cave Formation*

Figure 4 shows the full-wave waveform of cave formation. In this study, a monopole source is selected. The P and S-wave cannot be excited when the center frequency of the source is 2 kHz. There is only a ST wave in the borehole. The ST wave is unchanged because the cave in the formation cannot directly affect it in Figure 4a. The P, S and pseudo-Rayleigh waves appear in the wave train (Figure 4b–d) with the increase of the source frequency. Scattering waves are formed when the P and S-waves encounter the water-filling cave in the formation. Borehole acoustic waves are also changed because of the superposition with scattered waves [30]. In order to suppress pseudo-Rayleigh wave, the center frequency of source is usually 8–12 kHz.

– **Figure 4.** Borehole acoustic wave formed by different excitation frequency. Center frequency of source in (**a**–**d**) are 2 kHz, 5 kHz, 8 kHz and 10 kHz, respectively. Radius of the cave is 2 cm and the distance of the cave to the borehole wall is 2 cm.

– – – – — – – The center frequency of monopole source was 10 kHz in this study. Figure 5. shows the wave field of acoustic wave at different times in cave formation. The acoustic wave propagated along the borehole or formation to the PML boundary of the model. There was no artificial boundary reflection, and this proves that the NPML was effective. The P-wave encountered the cave at 0.9 ms. Because the scattered wave formed by P-wave was weak, only the original wave field excited by source can be seen in Figure 5a. The S–P and S–S waves were formed when an S-wave encountered a cave due to energy conversion. The S–P wave can be observed in Figure 5b beside the original wave field formed by the source at 1.3 ms. The scattered-wave field formed by the S–wave can be clearly observed in Figure 5c—including S–P and S–S waves. It should be noted that the amplitude of the S-wave decreased significantly after it passed through the cave in Figure 5c (the color of the S-wave became lighter).

**Figure 5.** Wave-field snapshot of acoustic waves in a cave formation. (**a**–**c**) Wave fields at 0.9, 1.3 and **Figure 5.** Wave-field snapshot of acoustic waves in a cave formation. (**a**–**c**) Wave fields at 0.9, 1.3 and 1.6 ms, respectively. Bar scale represents the amplitude of the wave**.**

*3.2. The E*ff*ect of the Cave on the Borehole* P-Wave

Figure 6a shows the effect of different radii of cave on P-waves.

**Figure 6.** Effect of the cave on P-wave. (**a)** shows the effect of the cave radius on P-wave. (**b)** shows the effect of the cave position on P-wave. *R* represents the radius of the cave. *L* represents the distance of the cave to the borehole wall. Center frequency of source is 10 kHz. P-wave recorded by the receiver is 0.4 m above the cave.


–

**Figure 7.** Effect of cave and excitation frequency on the P-wave. (**a**,**b**) Excitation frequency is 10 kHz. (**c**) Excitation frequency is 7 kHz. *R* represents the radius of the cave. *L* represents the distance of the cave to the borehole wall. P-wave recorded by the receiver is 0.4 m above the cave.


 – Because the cave was fully filled with water and the wave propagation speed in the fluid was less than that in the formation, the first arrival of the P-wave was delayed. The larger the radius of the cave, the longer was the P-wave propagation time in the fluid, and the later the first arrival time of the P-wave. In this study, a single cave was simulated. There may be multiple caves of different sizes in an actual formation, which will cause more delays on the first arrival of the P-wave [31]. P-wave recorded by the receiver was the superposition of scattered wave and "gliding" P-wave refracted on the borehole wall. The scattered wave formed by the cave changed with different distances to the borehole wall (*L*), so the amplitude of the P-wave may be different when the scattered wave is superimposed with it. In this study, the P-wave amplitude was reduced slightly, and it may increase if the source distance was changed because the scattered wave also varied with different source distance.

#### *3.3. The E*ff*ect of the Cave on Borehole S-Wave*

Figure 8a shows the effect of the radius of the cave on the S-wave.


**Figure 8.** Effect of the radius and position of cave on S-wave. (**a)** shows the effect of the cave radius on S-wave. (**b)** shows the effect of the cave position on S-wave. *R* represents the radius of the cave and *L* represents the distance of the cave to the borehole wall. Center frequency of source is 10 kHz. S-wave recorded by the receiver 0.4 m above the cave.

2 When the S-wave encounters the water-filled cave, the water inhibits the propagation of the S-wave and a considerable part of the energy of the S-wave is converted into the compressed energy (the energy of S–P wave) of the fluid [31]. The compression energy may be converted into S-wave and other mode waves, and these converted waves are recorded by the receiver finally. Therefore, the S-wave energy recorded by the receiver above the cave will be partially attenuated due to scattering and conversion.

= 1 −

#### 1 *3.4. The Relationship between the Attenuation Coe*ffi*cient of S-Wave and the Cave*

<sup>1</sup> <sup>2</sup> It is found that the existence of a cave can make the S-wave energy decay through simulation and analysis. In this study, the attenuation coefficient A of the S-wave under the action of the cave is defined as

$$A = 1 - \frac{A\_2}{A\_1} \tag{13}$$

– – – – – – – where *A*<sup>1</sup> and *A*<sup>2</sup> are the amplitude of the S-wave for a homogeneous formation and a cave formation in the same source distance, respectively. The dotted line is a fitting linear equation in Figure 9a, and there is a positive linear correlation between the attenuation coefficient and the radius of the cave. The slope of linear equation of attenuation coefficient is variable when the distance of cave is different. The large slope indicates that the attenuation coefficient is sensitive to the radius of the cave. When *L* = 4 cm, the attenuation coefficient is the most sensitive to the cave radius in our model.

– When the S-wave propagates in the formation and run into the water-filled cave, it will form S–P and S–S waves. The S–P and S–S waves recorded by the receiver near the cave are superposed together, so the attenuation coefficient is smaller as shown in Figure 9b. The S–P and S–S waves are separated gradually with the increase of the source distance. When the S–P leaves the S-wave, it takes its energy away from the S-wave, so the attenuation coefficient increases [31]. The effect of the cave on the S-wave is reduced with the increase of the source distance, so the attenuation coefficient decreases gradually. In our model, the S–S wave was separated completely when the source distance is 2.6 m as can be seen from Figure 9b.

**Figure 9.** Effect of cave and source distance on attenuation coefficient of the S-wave. (**a)** shows the effect of the cave radius on attenuation coefficient. (**b**) shows the effect of the source distance on attenuation coefficient. *R* represents the radius of the cave and *L* represents the distance of the cave to the borehole wall. Center frequency of source is 10 kHz.

#### *3.5. The E*ff*ect the of Cave on Stoneley Wave*

– – The P and S-wave which can propagate in the formation are formed when the center frequency of source is high. At this time, the scattered wave propagates to the borehole and affects the ST wave. Figure 10a is the wave train in homogeneous elastic formation. Figure 10b–f represent the wave train in the formation with cave of different radii. It can be seen that when the radius of the cave is small (Figure 10b,c), the attenuation of the ST wave is weak compared with Figure 10a. When the radius of the cave is large (Figure 10d–f), the amplitude of the ST wave has a great attenuation. Therefore, the attenuation of the ST wave not only corresponds to fractures, but also to caves. – –

– – **Figure 10.** Effect of the radius of the cave on ST wave. (**a)** shows the wave train in homogeneous elastic formation. (**b**–**f)** represent the wave train in the formation with cave of different radii. *R* represents the radius of the cave. The distance of the cave to the borehole wall is 2 cm. Center frequency of source is 10 kHz.

Figure 11a shows acoustic waveforms. In order to study the relationship between the ST wave and the cave radius, it is necessary to calculate the energy of the ST wave. The full-wave frequency spectrum is obtained by Fourier transform in Figure 11b; the low-pass filter is used to obtain ST wave. The energy of the ST wave is obtained by summing the square of the corresponding value of each frequency. The frequency of the ST wave is less than 4 kHz according to Figure 11b. The cutoff frequency of the low-pass filter is 4 kHz. The black line in Figure 11a is the low-pass filtered waveform. The S-wave and pseudo-Rayleigh wave have been suppressed to a certain extent. However, due to the complexity of acoustic wave and the limitation of Fourier transform [32], S and pseudo-Rayleigh waves are difficult to be completely suppressed. Figure 11c represents the relationship between ST wave energy and cave radius at different source distances. The energy of the ST wave gradually decreases with the increase of cave radius at different source distances.

**Figure 11.** Calculation of ST wave energy and its relationship with cave radius. (**a**) shows the original and the filtered waves. (**b**) shows the frequency spectrum of original wave. (**c**) shows the effect of cave radius on the energy of ST wave. Center frequency of source is 10 kHz. *R* represents the radius of the cave. *z* is the axial coordinate.

#### *3.6. The E*ff*ect of the Cave on the Scattered Wave*

– – Figure 12a refers to the acoustic wave and the scattered wave. The wave propagation in the borehole of homogeneous formation is regarded as a direct wave, and the wave trains of cave formation in Figure 12a subtract the direct waves to separate the scattered waves. A Fourier transform was used to get the frequency spectrum of the scattered wave in Figure 12b. The peaks with larger amplitude are mainly S–S wave, and the dominant frequency is about 9–12 kHz. The peaks with smaller amplitude are mainly converted P-waves. We added the square of the value of each frequency to get the energy of the scattered wave. Taking the distance and radius of the cave as the abscissa and scattered wave energy as the ordinate, the relationship between the scattered wave energy and the cave was obtained in Figure 13.

The scattered wave recorded at 0.4 m below the cave (the source distance was 1.6 m) was selected. The cave with a diameter less than 0.25 times the wavelength of the S-wave was approximated as a point scatterer; otherwise the cave itself could form multiple-scattering [30]. In our model, a cave with a radius less than 2 cm could be regarded as a point-scatterer and could form multiple-scattering with a radius greater than 3 cm.

Figure 13a shows the relationship between the energy of scattered wave and the position of cave. When the cave can be regarded as a point-scatterer (*R* ≤ 2 cm), the energy of the scattered wave has a negative correlation with the distance of the cave. The energy of scattered wave is very small when the distance is more than 24 cm and reaches the maximum when *L* = 2 cm if the cave itself can form multiple-scattering (*R* > 3 cm). The scattered wave energy gradually decreases when *L* > 2 cm (*R* > 3 cm).

Figure 13b shows the relationship between the energy of scattered wave and the radius of the cave. The energy of scattered wave raises with the increase of cave radius if *R* ≤ 2 cm. If *R* > 3 cm, the energy of the scattered wave changes sharply when the cave is close to the borehole wall (*L* = 2 cm) and is relatively stable when L > 2 cm with the increase of cave radius.

– –

**Figure 12.** Full-wave acoustic-logging in cave formation. (**a**) shows the scattered wave in the cave formation. (**b**) shows the frequency spectrum of the scattered wave. Radius of cave is 2 cm. *L* represents the distance of the cave. Center frequency of source is 10 kHz. ≤

**Figure 13.** The relationship between the energy of the scattered wave and the cave. (**a**) shows the effect of cave position on the energy of the scattered wave. (**b**) shows the effect of cave radius on the energy of the scattered wave. *R* represents the radius of the cave*. L* represents the distance of the cave. Center frequency of source is 10 kHz.

#### **4. Discussion**

The numeric simulation of borehole acoustic waves in cave formations is a theoretical calculation; its ultimate purpose is to provide interpretation and reference for actual data. The results in our study show that borehole full-waves have response characteristics to the caves, which are mainly reflected in the energy and the first arrival of waves.

The porosity-frequency spectrum is calculated from the conductivity curve of electrical-imaging logging [33] and follows the following rules.


…

(3) When the primary pores and secondary pores are well developed, but heterogeneous, the porosity-distribution spectrum shows multiple peaks. The peaks move to the direction of large porosity.

Cumulative porosity curves PS3, PS5, PS7 . . . PS30 represent pore components with porosity values below 0.03, 0.05, 0.07 . . . 0.3 [34]. The porosity of electrical imaging can be obtained by filling the profile with different colors; the contribution of the different sizes of pores to the total porosity can be intuitively analyzed. According to the method in Figure 11b mentioned above, the energy of the ST wave is calculated. The acoustic full-wave was decomposed into several intrinsic mode function (IMF) components by empirical mode decomposition (EMD). P and S-wave were distributed in the highest frequency IMF1 component [35]. The energy of IMF1 component can be approximately regarded as the energy of the S-wave because the amplitude of the P-wave is very small. The slowness of the P-wave is the reciprocal of velocity. The first arrival delay of the P-wave indicates that the velocity is low, then the slowness is large.

Figure 14 shows the logging results of H-well in the igneous formation of Liaohe. Area **a** is basalt. Cores 1–2 show no alteration. The color of electrical-imaging logging is bright, indicating that the conductivity is low. The porosity-distribution spectrum shows a single peak moving forward, and the porosity of electrical imaging shows that small pores are mainly developed in the formation. The energy of the S and ST waves is high and stable. The density is large, and the porosity and slowness of the P-wave are small. All the above logging data indicate that this area is a tight formation. Cores 3–4 show altered basalt in area **b**. Dark spots appear in the electrical-imaging logging static image, which indicates that as conductivity increases, and caves develop. The porosity-distribution spectrum widens and moves in the direction of large porosity. It was found that there are mainly medium and large pores from the electrical imaging porosity. The energy of the S and ST waves was attenuated, and the slowness of the P-wave increased. This is in good agreement with the theoretical results of the cave formation simulated. – –

– – **Figure 14.** Comprehensive interpretation of logging data in Liaohe igneous formation. (**a**) Depth of the well; (**b**) GR is natural gamma and CAL is caliper; (**c**). FMI IMAGE represents electrical-imaging logging static image; (**d**) FDIST is the porosity-distribution spectrum; (**e**) PSMIN–PSMAX are the porosity of electrical imaging; (**f**–**h**) TFWV, PSWAVE and STWAVE are the waveforms of full-wave, P, S and ST waves, respectively; (**i**) SENERGY represents the energy of the S-wave. STENERGY is the energy of ST wave; (**j**) DEN is the density of the formation. POR is the porosity by conventional logging. DTC is the slowness of the P-wave. Figures 1 and 2 on the right are the core slice and core of the tight formation. Figures 3 and 4 are the core slice and core of the altered formation.

Figure 15 shows the effect of caves with different development degree on borehole acoustic waves. The attenuation of the S and ST waves becomes stronger, and the slowness of the P-wave becomes larger with an increasing number and size of caves.

– **– Figure 15.** Comprehensive interpretation of logging data of cave formations with different development degrees. (**a**) Depth of the well; (**b**) GR is natural gamma and CAL is caliper; (**c**). FMI IMAGE represents electrical-imaging logging static image; (**d**) FDIST is the porosity-distribution spectrum; (**e**) PSMIN–PSMAX are the porosity of electrical imaging; (**f–h**) TFWV, PSWAVE and STWAVE are the waveforms of full-wave, P, S and ST waves, respectively; (**i**) SENERGY represents the energy of the S-wave. STENERGY is the energy of ST wave; (**j**) DEN is the density of the formation. POR is the porosity by conventional logging. DTC is the slowness of the P-wave.

The theoretical results of this study can be applied to practical data, but there are also some shortcomings.


#### **5. Conclusions**

The following important conclusions can be drawn:


**Author Contributions:** F.X. is the main author of this manuscript and this work was conducted under the advisement of Z.W.; Z.W. is the supervision, project administration. Z.W. reviewed the manuscript and made contributions to its structure. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, Grant Number: 41874135.

**Conflicts of Interest:** The authors declare no conflicts 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* **An Immersed Boundary Method Based Improved Divergence-Free-Condition Compensated Coupled Framework for Solving the Flow–Particle Interactions**

**Pao-Hsiung Chiu <sup>1</sup> , Huei Chu Weng 2,\* , Raymond Byrne <sup>3</sup> , Yu Zhang Che <sup>4</sup> and Yan-Ting Lin 5,\***

	- <sup>3</sup> Centre for Renewables and Energy, School of Engineering, Dundalk Institute of Technology, Dundalk,
	- Tel.: +886-3-2654311 (H.C.W.); +886-3-4711400 (ext. 3356) (Y.-T.L.)

**Abstract:** A flow–particle interaction solver was developed in this study. For the basic flow solver, an improved divergence-free-condition compensated coupled (IDFC<sup>2</sup> ) framework was employed to predict the velocity and pressure field. In order to model the effect of solid particles, the differentially interpolated direct forcing immersed boundary (DIIB) method was incorporated with the IDFC<sup>2</sup> framework, while the equation of motion was solved to predict the displacement, rotation and velocity of the particle. The hydrodynamic force and torque which appeared in the equations of motion were directly evaluated by fluid velocity and pressure, so as to eliminate the instability problem of the density ratio close to 1. In order to effectively evaluate the drag/lift forces acting on the particle, an interpolated kernel function was introduced. The present results will be compared with the benchmark solutions to validate the present flow–particle interaction solver.

**Keywords:** immersed boundary method; quasi multi-moment method; incompressible Navier– Stokes equation; dispersion-relation-preserving; flow–structure interaction

#### **1. Introduction**

The offshore wind farm is a wind turbine installed in an offshore area where there is a strong wind field for a long period of time. It can avoid the noise, shading and visual obstruction caused by the wind turbine's operation in the land area. Offshore wind farms can also provide more stable wind energy and low turbulence effects. However, it requires consideration of the effects of wind turbulence on the wind turbine and its supporting structure. Recently, attention has also been paid to the offshore floating wind turbine. How to efficiently simulate the interaction between the floating structure and the water/air interface becomes an important research topic.

The coupling analysis of fluids and solids in practical engineering applications usually requires the consideration of complex geometries. These complex geometries are usually in a stationary or transient motion at high Reynolds number flows, such as offshore stationary wind turbines and marine systems. These fluid–structure problems are typically solved by the traditional body-fitted-grid. In order to generate high-quality mesh, a lot of efforts in terms of manpower are always required to refine the mesh at the high curvature region. The computational cost is then increased, especially when dealing with the moving boundary problems. In addition, the quality of the grid is also very crucial for simulation analysis.

**Citation:** Chiu, P.-H.; Weng, H.C.; Byrne, R.; Che, Y.Z.; Lin, Y.-T. An Immersed Boundary Method Based Improved Divergence-Free-Condition Compensated Coupled Framework for Solving the Flow–Particle Interactions. *Energies* **2021**, *14*, 1675. https://doi.org/10.3390/en14061675

Academic Editors: Robert Castilla, Anton Vernet and Antonio F. Miguel

Received: 23 December 2020 Accepted: 13 March 2021 Published: 17 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 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 (https:// creativecommons.org/licenses/by/ 4.0/).

In recent decades, the immersed boundary method has been shown to have great potential for modeling the flow–particle interaction problem. It is achieved by virtue of incorporating introduced momentum forcing terms in the equations of motion to predict the displacement, rotation and velocity of the solid object [1]. The main advantage of employing an immersed boundary method for this kind of problem is that the regular fixed Cartesian grid can still be used for simulating the complex time-varying geometries.

However, how to effectively evaluate the drag/lift of the coupling interface and the particle velocity is still a research topic to be refined. Uhlmann [2] proposed incorporating the regularized delta function approach to model the particulate flows, and the improved methods were then proposed by Kempe and Fröhlich [3] and Breugem [4]. Zhang and Gay [5] also proposed a finite element method based on an immersed boundary for fluid– structure interactions. Alternatively, Luo et al. [6] simulated the fluid–particle interactions by direct-forcing based modified immersed boundary method. To further improve the accuracy of velocity in moving solid, Wang et al. [7] introduced the multi-direct forcing method. Yang and Stern [8] also demonstrated a sharp interface-based immersed boundary method by virtue of the direct-forcing approach for solving particulate flows, while Horng et al. [9] exhibited a prediction-correction based direct-forcing immersed boundary projection method for solving fluid–particle interactions.

With the success of finite volume and finite element frameworks, scientists are also motivated to employ an immersed boundary based fluid–particle solver on the Lattice– Boltzmann (LB) framework. Wu and Shu [10] demonstrated a boundary condition-enforced immersed boundary LB scheme to solve particulate flows, while Hao and Zhu [11] revealed an implicit immersed boundary method on the Lattice–Boltzmann framework for fluid– structure interactions. Wang et al. [12] later evaluated three LB schemes to understand their performance when modeling the particulate flows. Zhou and Fan [13] also presented a LB immersed boundary scheme based on the work of Breugem [4]. Zhang et al. [14] further exhibited the particulate immersed boundary method (PIBM) to speed up the calculation of the particulate flow simulations. Coclite et al. [15] proposed kinematic and dynamic forcing strategies for predicting the transport of inertial capsules, and later extend to model the deformable inertial capsules [16].

Due to the fact that the convergence of the Lattice-Boltzmann framework is based on equilibrium equations, each fluid must have its own balance equation for the convergence of the Lattice–Boltzmann method. Furthermore, the definitions of boundary conditions in the Lattice–Boltzmann method is still developing for fitting the physical properties of interfaces. Moreover, some previous studies still only applied the prescribed motion of solid particle for the fluid–solid interaction simulation.

To alleviate the issue of evaluating drag and lift forces, Chiu and Poh [17] have recently successfully incorporated the improved divergence-free-condition compensated coupled (IDFC<sup>2</sup> ) framework with the direct forcing immersed boundary (DIIB) method [18] for solving the flows with prescribed-motion time-varying geometries. The spurious force oscillation (SFO) can be efficiently alleviated and the calculation is relatively simple. This motivated us to incorporate the IDFC<sup>2</sup> framework with equation of particle motion for simulating the fluid–particle interactions. The ultimate goal for this study is to develop the solver to model the interaction between the floating wind turbine and water/air interface.

In this study, the framework will be extended for simulating the fluid–particle interactions. To prevent the instability problem when the density ratio is close to 1, the drag and lift forces will be directly evaluated from fluid velocity and pressure. An interpolated kernel function to accurately and effectively evaluate the drag/lift forces will be proposed. Details of calculating the equation of motion as well as near-wall treatment when the solid particle is approaching the wall will also be addressed.

This paper is organized as follows. Section 2 presents the governing equations for fluid flows and solid particles. Section 3 present the DIIB method and the methodology for fluid–particle interactions. Section 4 presents the IDFC<sup>2</sup> based flow–particle interaction framework. Verification studies are presented in Section 5. This is followed by presenting the simulated results to show the applicability for the proposed framework. Finally, some concluding remarks are drawn in Section 6.

#### **2. Governing Equations**

#### *2.1. Incompressible Navier–Stokes Equations*

The non-dimensional primitive-variable Navier–Stokes equations for incompressible flows can be expressed as the following momentum equations and continuity equation:

$$\frac{\partial \stackrel{\smile}{u}}{\partial t} + \nabla \cdot \left(\stackrel{\smile}{u} \stackrel{\smile}{u}\right) = \frac{1}{Re} \nabla^2 \stackrel{\smile}{u} - \nabla p + \stackrel{\smile}{f} \tag{1}$$

$$\nabla \cdot \stackrel{\rightharpoonup}{u} = \mathbf{0} \tag{2}$$

where ⇀ *u* is the velocity vector, *p* is the pressure field, and *Re* is the Reynolds number. ⇀ *f* is external force vectors due to particle motions, and is modeled by the present immersed boundary based flow–particle interaction framework.

#### *2.2. Equations of Motion for Solid Particle*

directly evaluated in this study:

The equations of motion for the rigid body can be written as

$$
\rho\_S V\_s \frac{d\overline{U}}{dt} = (\rho\_S - \rho\_f) V\_s \overline{\hat{\mathbf{g}}} + \overline{\hat{F}\_f} \tag{3}
$$

$$\frac{d\vec{X}}{dt} = \vec{\mathcal{U}}\tag{4}$$

$$I\_S \frac{d\vec{\omega}}{dt} = \vec{\psi}\_f \tag{5}$$

$$\frac{d\vec{\Omega}}{dt} = \vec{\omega} \tag{6}$$

where the linear velocity vector ( ⇀ *V*) and particle position vector ( ⇀ *X*) are due to linear motion, while *ω* and Ω are the angular velocity and orientation due to the angular motion. *ρ<sup>S</sup>* and *ρ<sup>f</sup>* are the density for the solid particle and fluid, respectively. *V<sup>s</sup>* is the volume of the solid particle, ⇀ *g* is the gravity, and *I<sup>S</sup>* is the inertial moment of the solid particle. ⇀ *F<sup>f</sup>* and ⇀ *ψ<sup>f</sup>* are the hydrodynamic force and torque acting on the solid particle, respectively.

*d*

As stated in [2,3], when we evaluated the ⇀ *F<sup>f</sup>* and ⇀ *ψ<sup>f</sup>* by the summation of forcing terms ⇀ *f* together with the relation of momentum balance over the corresponding fluid domain [2], it will lead to an instability issue when the density ratio close to 1, or (*ρ<sup>S</sup>* − *ρ<sup>f</sup>* ) ≈ 0. This is due to the fact that the term (*ρ<sup>S</sup>* − *ρ<sup>f</sup>* ) will appear in the denominator when solving the equations of motion. In order to avoid this issue, hydrodynamic force and torque are

$$
\stackrel{\rightarrow}{F}\_f = \rho\_f \int \underline{\sigma} \cdot \underline{n} \, dA \tag{7}
$$

$$
\overrightarrow{\psi}\_f = \rho\_f \int r \times (\underline{\underline{\sigma}} \cdot \underline{n}) \, dA \tag{8}
$$

It is noted that by using the present approach, (*ρ<sup>S</sup>* − *ρ<sup>f</sup>* ) will be never appear in the denominator, so that there will be no instability problems even when the density ratio id close to 1. In Section 3.1, the present study will introduce methods to accurately and effectively evaluate the two above equations. Furthermore, the implicit midpoint temporal scheme was employed to solve Equations (3)–(6) in the present study.

#### **3. Differentially Interpolated Direct Forcing Immersed Boundary (DIIB) Method**

In this section, the DIIB method is employed to model the solid body. Here, we will briefly describe the fundamental idea of the DIIB method.

Assuming there is a solid object *S* immersed in the Eulerian grid; in order to model the effect of solid object, the direct forcing term ⇀ *f* is introduced in Equation (1) to satisfy the velocity condition ( ⇀ *V*) of the solid object. With the semi-implicit discretization, Equation (1) can be split into the following two equations: ⃑ ⃑ ∗− ⃑ 1

∗

$$\frac{\stackrel{\rightharpoonup}{u}^\* - \stackrel{\rightharpoonup}{u}^n}{\Delta t} = -\nabla \cdot \left(\stackrel{\rightharpoonup}{u}\stackrel{\rightharpoonup}{u}\right)^\* = \frac{1}{Re}\nabla^2 \stackrel{\rightharpoonup}{u}^\* - \nabla p \tag{9}$$

$$\frac{\stackrel{\rightharpoonup}{V} - \stackrel{\rightharpoonup}{u}^\*}{\Delta t} = \stackrel{\rightharpoonup}{f} \tag{10}$$

It is noted that in practice, the boundary points of the solid object do not always fall into the Eulerian grid points. The direct forcing method will then evaluate the modified velocity on some target Eulerian grid points, followed by obtaining the corresponding momentum forcing terms based on Equation (10) to model the effect of solid object with velocity condition ⇀ *V*. ⃑ 

For the DIIB method, the velocity (*uS*) at Eulerian grid point which is inside and close to the solid boundary, will be modified to satisfy the velocity condition *uIB* at the solid boundary by virtue of interpolated velocity at fluid domain *u<sup>f</sup>* and Taylor series expansion: = 2 −

$$
u\_S = 2\boldsymbol{u}\_{IB} - \boldsymbol{u}\_f \tag{11}$$

Taking Figure 1 as an example, in general *u<sup>F</sup>* does not always lie on the grid point. The evaluation of *u<sup>f</sup>* is done by "advect" *u<sup>f</sup>* from point Q to point A [18]: ∂ + (⃑ ∙ ∇) = 0 ⃑

$$\frac{\partial \mu\_f}{\partial \tau} + \left(\stackrel{\rightharpoonup}{n} \nabla \right) \mu\_f = 0 \tag{12}$$

 **Figure 1.** Schematic of the direct forcing immersed boundary (DIIB) method: (**a**) evaluation of *uS*; and (**b**) evaluation of *u<sup>f</sup>* .

In the above, ⇀ *n* is the unit normal vector of the solid object. In practice, Equation (12) is solved by first order upwind with single artificial time step ∆*τ* = *APQ*/ ⇀ *n* = 2*AP*/ ⇀ *n* . The reader can refer to [18] for details.

To alleviate the spurious force oscillations (SFOs) and to obtain smooth pressure, the velocity modification method [17] was employed in this study. The idea behind this method is to modify the velocity inside the velocity to the solid object moving velocity

when solving the pressure correction equation. It has been shown in [17] that this method successfully suppresses the SFO for the problems with prescribed motions and will then be implemented for dealing with the problem for which the motions are obtained by solving equations of motion.

#### *3.1. Evaluation of Drag and Lift Forces*

When employing the direct forcing immersed boundary method, it is suggested to evaluate forces *F* by the following equation for calculating drag and lift forces:

$$
\underline{F} = \int \underline{\sigma} \cdot \underline{n} \, dA \tag{13}
$$

where *σ* is the summation of shear and pressure stress, *n* is unit normal vector and *A* is surface area. In this study, a four-point piecewise discrete delta function based interpolated kernel is employed [17]:

$$F(X,Y) = \sum F(\mathbf{x}, \mathbf{y}) \phi'\left(\frac{\mathbf{x} - \mathbf{X}}{h}\right) \phi'\left(\frac{\mathbf{y} - \mathbf{Y}}{h}\right) \tag{14}$$

$$
\phi'\left(\frac{\mathbf{x} - \mathbf{X}}{h}\right) = \phi\left(\frac{\mathbf{x} - \mathbf{X}}{h}\right) / \sum \phi\left(\frac{\mathbf{x} - \mathbf{X}}{h}\right) \tag{15}
$$

$$
\phi'\left(\frac{y-\chi}{h}\right) = \phi\left(\frac{y-\chi}{h}\right) / \sum \phi\left(\frac{y-\chi}{h}\right) \tag{16}
$$

$$\phi(r) = \begin{cases} \frac{1}{8} \left( 3 - 2|r| + \sqrt{1 + 4|r| + 4r^2} \right), |r| \le 1 \\\frac{1}{8} \left( 5 - 2|r| - \sqrt{-7 + 12|r| - 4r^2} \right), 1 \le |r| \le 2 \\\ 0 & |r| \ge 2 \end{cases} \tag{17}$$

when the solid particle is approaching the wall, the interpolated kernel will be simplified to the following one-point piecewise discrete delta function [19]:

$$\phi(r) \begin{cases} 1 - |r|\_{\prime} & |r| \le 1 \\ 0 & |r| \ge 1 \end{cases} \tag{18}$$

#### **4. Improved Divergence-Free-Condition Compensated Coupled (IDFC<sup>2</sup> ) Framework**

The idea behind the IDFC<sup>2</sup> framework is to discretize momentum equations by virtue of cell-center and cell-face velocity [17]. With the present approach, we can not only obtain the accurate velocity vector (*u*, *v*), but also the fully coupled pressure field *p*. The following will first introduce the original IDFC framework, then briefly describe how to mitigate to IDFC<sup>2</sup> framework.

#### *4.1. Derivation of IDFC Framework*

With the 2D grid structure as shown in Figure 2, the derivation of the IDFC framework starts from the semi-discrete momentum equation. Without loss of generality, we take velocity *u* at the cell-center P (*up*) as an example:

$$\frac{\boldsymbol{u}\_{P}^{\*} - \boldsymbol{u}\_{P}^{n}}{\Delta t} + (\frac{\partial \boldsymbol{u} \boldsymbol{u}}{\partial \mathbf{x}})\_{P} + (\frac{\partial \boldsymbol{v} \boldsymbol{u}}{\partial \boldsymbol{y}})\_{P} = \frac{1}{\text{Re}} (\frac{\partial^{2} \boldsymbol{u}}{\partial \mathbf{x}^{2}} + \frac{\partial^{2} \boldsymbol{u}}{\partial \boldsymbol{y}^{2}})\_{P} - (\frac{\partial p^{n}}{\partial \mathbf{x}})\_{P} \tag{19}$$

**Figure 2.** Schematic of the present framework.

 (

 The idea for the IDFC method is that discretizing the momentum equation for the cell-face *u<sup>e</sup>* should have the same form as we discretize the momentum equation at the cell-center:

$$\frac{\partial u^\*\_{\varepsilon} - u^n\_{\varepsilon}}{\Delta t} + (\frac{\partial uu}{\partial \mathbf{x}})\_{\varepsilon} + (\frac{\partial vu}{\partial y})\_{\varepsilon} = \frac{1}{\text{Re}} (\frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\partial^2 u}{\partial y^2})\_{\varepsilon} - (\frac{\partial p^n}{\partial \mathbf{x}})\_{\varepsilon} \tag{20}$$

 <sup>∗</sup>− ∆ + ( ) + ( ) = 1 ( 2 <sup>2</sup> + 2 2 ) − ( ) However, obtaining the velocity derivatives on the cell-face is not as easy as on the cell-center. The linear-averaged approximation is then employed to evaluate the cell-face velocity derivatives on the cell-face *e* [20]:

$$(\frac{\partial \mu u}{\partial \mathbf{x}})\_{\varepsilon} \approx (\frac{\partial \hat{u}u}{\partial \mathbf{x}})\_{\varepsilon} = \frac{1}{2} (\left(\frac{\partial \mu u}{\partial \mathbf{x}}\right)\_{E} + \left(\frac{\partial \mu u}{\partial \mathbf{x}}\right)\_{P}) \tag{21}$$

$$(\frac{\partial vu}{\partial y})\_\epsilon \approx (\frac{\partial v\mathbf{\hat{u}}}{\partial y})\_\epsilon = \frac{1}{2}(\left(\frac{\partial vu}{\partial y}\right)\_E + \left(\frac{\partial vu}{\partial y}\right)\_P) \tag{22}$$

$$\frac{1}{Re}(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2})\_\varepsilon \approx \frac{1}{Re}(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2})\_\varepsilon = \frac{1}{2}(\frac{1}{Re}\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right)\_E + \frac{1}{Re}\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right)\_P) \tag{23}$$
 
$$\text{and} \quad \text{i.e.,} \quad \text{i.e.,} \quad \text{(24)} \quad (\text{23}): \quad \text{i.e.,} \quad \text{ii.,} \quad \text{iii.} \quad \text{i.e.,} \quad \text{i.e.,} \quad \text{i.e.,} \quad \text{iv.} \quad \text{i.e.,} \quad \text{iv.}$$

( ) ≈ ( ) = 2 (( ) + ( ) ) 1 2 2 1 2 2 ̂ 1 1 2 2 1 By using Equations (21)–(23) into Equation (20), *u* ∗ *e* can then be obtained. *u* ∗ *<sup>w</sup>*, *v* ∗ *n* , *v* ∗ *s* can also be obtained in a similar way. In order to fulfill the divergence-free condition, the Poisson equation for pressure correction is then solved with cell-face velocity:

$$\left(\frac{\partial^2 p'}{\partial x^2} + \frac{\partial^2 p'}{\partial y^2}\right)\_P = \frac{1}{\Delta t} (\frac{u\_e^\* - u\_w^\*}{\Delta x} - \frac{v\_n^\* - v\_s^\*}{\Delta y}) \tag{24}$$
 
$$\begin{array}{ccccccccc} \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \end{array}$$

It is followed by updating the cell-face velocity by virtue of the predicted pressure correction:

+1 =

+1 =

+1 =

+1 =

+ ′

$$
\mu\_{\varepsilon}^{n+1} = \mu\_{\varepsilon}^{\*} - \Delta t \left(\frac{\partial p'}{\partial \mathcal{X}}\right)\_{\varepsilon} \tag{25}
$$

 2 <sup>2</sup> + 2 2 ) )

<sup>∗</sup>

<sup>∗</sup> ∗

$$\mu\_w^{n+1} = \mu\_w^\* - \Delta t (\frac{\partial p'}{\partial \mathfrak{x}})\_w \tag{26}$$

$$v\_n^{n+1} = v\_n^\* - \Delta t (\frac{\partial p'}{\partial y})\_n \tag{27}$$

$$v\_s^{n+1} = v\_s^\* - \Delta t (\frac{\partial p'}{\partial y})\_s \tag{28}$$

<sup>∗</sup> − ∆(

<sup>∗</sup> − ∆(

<sup>∗</sup> − ∆(

′ )

′ )

′ ) while the pressure is updated as *p <sup>n</sup>*+<sup>1</sup> = *p <sup>n</sup>* + *p* ′ . Finally, the cell-center velocity is updated by again employing the linear-averaged approximation for the cell-center derivatives:

$$\frac{\partial u\_P^\* - u\_P^n}{\Delta t} + (\frac{\partial u \hat{u}}{\partial \mathbf{x}})\_P + (\frac{\partial v \hat{u}}{\partial y})\_P = \frac{1}{Re} (\frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\partial^2 u}{\partial y^2})\_P - (\frac{\partial p^n}{\partial \mathbf{x}})\_P \tag{29}$$

where:

$$\left(\frac{\partial \hat{u}u}{\partial \mathbf{x}}\right)\_P = \frac{1}{2}(\left(\frac{\partial u \Omega u}{\partial \mathbf{x}}\right)\_\varepsilon + \left(\frac{\partial u \Omega u}{\partial \mathbf{x}}\right)\_w) \tag{30}$$

$$(\frac{\partial v\partial u}{\partial y})\_P = \frac{1}{2}(\left(\frac{\partial v\partial u}{\partial y}\right)\_\varepsilon + \left(\frac{\partial v\partial u}{\partial y}\right)\_w) \tag{31}$$

$$\frac{1}{Re}(\frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\partial^2 u}{\partial y^2})\_P = \frac{1}{2}(\frac{1}{Re}\left(\frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\partial^2 u}{\partial y^2}\right)\_\varepsilon + \frac{1}{Re}\left(\frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\partial^2 u}{\partial y^2}\right)\_w) \tag{32}$$

It is noted that by virtue of Equations (30)–(32), the introduced IDFC forcing terms can efficiently stabilize the calculations [20].

#### *4.2. Derivation of IDFC<sup>2</sup> Framework*

It is known from [17] that the IDFC method may lead to non-physical velocity oscillations due to the fact that the cell-face and cell-center velocity is not strongly coupled when evaluating the velocity derivatives. To resolve this issue, we modified Equation (23) to include the contribution of cell-face velocity:

$$\underbrace{1}\_{\text{Re}} (\frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\partial^2 u}{\partial y^2})\_\varepsilon}\_{\varepsilon} = \frac{1}{\text{Re}} (\mathbb{C} \left( \frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\partial^2 u}{\partial y^2} \right)\_\varepsilon + (1 - \mathbb{C}) \left( \frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\partial^2 u}{\partial y^2} \right)\_\varepsilon \tag{33}$$

In the above, *C* is the newly introduced parameter which is utilized as coupling cell and face velocity. In this study, *C* was chosen as 1/2 to ensure the strong coupling.

The approximation equation for cell-face velocity *u* ∗ *e* is then rewritten as

$$\frac{\mathbf{u}\_{\varepsilon}^{\star} - \mathbf{u}\_{\varepsilon}^{\mathrm{n}}}{\Delta t} + (\frac{\partial \mathbf{u} \mathbf{\hat{u}}}{\partial \mathbf{x}})\_{\varepsilon} + (\frac{\partial \mathbf{v} \mathbf{\hat{u}}}{\partial \mathbf{y}})\_{\varepsilon} = \frac{1}{\mathrm{Re}} (\underbrace{\frac{\partial^{2} \mathbf{u}}{\partial \mathbf{x}^{2}} + \overbrace{\partial^{2} \mathbf{u}}^{2}}\_{\varepsilon})\_{\varepsilon} - (\frac{\partial p^{n}}{\partial \mathbf{x}})\_{\varepsilon} \tag{34}$$

The corresponding cell-center velocity is expressed as

$$\frac{u\_P^\* - u\_P^n}{\Delta t} + \left(\frac{\partial u\Omega}{\partial \mathbf{x}}\right)\_P + \left(\frac{\partial v\Omega}{\partial y}\right)\_P = \frac{1}{\text{Re}} \widehat{\left(\frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\partial^2 u}{\partial y^2}\right)\_P} - \left(\frac{\partial p^n}{\partial \mathbf{x}}\right)\_P \tag{35}$$

$$\underbrace{1}\_{\text{Re}} (\underbrace{\partial^2 u}\_{\text{Re}} + \overbrace{\partial y^2}^{\text{Re}})\_P = \frac{1}{2} (\underbrace{1}\_{\text{Re}} \left( \overbrace{\partial x^2}^{\text{\text{\textdegree}}} + \overbrace{\partial y^2}^{\text{\textdegree}} \right)\_\text{\textdegree} + \frac{1}{\text{Re}} \left( \overbrace{\partial x^2}^{\text{\textdegree}} + \overbrace{\partial y^2}^{\text{\textdegree}} \right)\_w \tag{36}$$

The overall algorithm for the present DIIB–IDFC<sup>2</sup> framework is summarized as follows. The linear and angular velocity, position and rotation of the solid objects are first solved by the equations of motion, where the hydrodynamic force and torque are obtained from the fluid velocity and pressure. To compute the summation of the hydrodynamic force and torque, the solid surface is first divided by a collection of Lagrangian points, and then interpolates the fluid velocity and pressure on each Lagrangian point. Finally, the linear and angular velocity of solid objects are then employed as the boundary condition for the solid–fluid interface. For the sake of completeness, we also plotted the algorithm in Figure 3. In this study, a finite volume-based, second order accurate dispersion-relation preserving upwinding scheme [20] was utilized to discretize the convection terms, while the central difference was employed for other derivative terms. The second order semi-implicit

Gear method was employed as the temporal scheme. For more details of the numerical scheme, people can refer to [20].

**Figure 3.** Overall algorithm for the present framework.

#### **5. Results**

#### *5.1. Taylor-Couette Flow*

– ) ) The 2-D Taylor–Couette flow problem has been chosen as the first validation problem. For the present problem, the radius for outer circular (*RO*) and inner circular (*RI*) are 0.4 and 0.2. The non-slip boundary condition is employed for the outer circular, and a constant rotating angular velocity *ωTC* = 3 is set as a boundary condition for the inner circular. The exact solution for this kind of flow setting can be derived as

$$u(x,y) = -K\left(\frac{R\_O^2}{r^2} - 1\right)y\tag{37}$$

$$v(x,y) = K\left(\frac{R\_O^2}{r^2} - 1\right)\mathbf{x} \tag{38}$$

$$p(x,y) = \mathcal{K}^2(\frac{r^2}{2} - \frac{R\_O^2}{2r^2} - R\_2^2 \log\left(r^2\right))\tag{39}$$

= √

$$K = \frac{\omega\_{TC} R\_I^2}{R\_O^2 - R\_I^2} \tag{40}$$

$$r = \sqrt{x^2 - y^2} \tag{41}$$

<sup>2</sup> − 2

Δ = Δ = Δℎ = 1/10 1/20

 = 2 <sup>2</sup> − 2 For the present study, four different mesh sizes, namely ∆*x* = ∆*y* = ∆*h* = 1/10, 1/20, 1/40 and 1/80 with *Re* = 500 are conducted for the validation study in a unit square domain. The inner and outer cylinder are modeled with DIIB method, as shown in Figure 4.

–

1/40 1/80 = 500

–

= 3

(, ) = − (

(, ) = (

2 ( 2

2 −

=

= √

(, ) =

 2

<sup>2</sup> − 1)

<sup>2</sup> − 1)

<sup>2</sup> − <sup>2</sup> 2

2

log(

2 ))

 2

 2

2

 2

> <sup>2</sup> − 2

<sup>2</sup> −

) )

Δ = Δ = Δℎ = 1/10 1/20

– **Figure 4.** Schematic of the Taylor–Couette flow problem.

From Table 1 and Figure 5, this shows that the rate of convergence for *u*, *v* and *p* are 2.033, 2.033 and 2.044, which match very well with the proposed accuracy order. The validation of the present solver was then confirmed. 

**Table 1.** The error L-2 norms for the Taylor–Couette flow problem. –


– **Figure 5.** The rate of convergence plot for the Taylor–Couette flow problem. (**a**) *u*; (**b**) *v*; and (**c**) *p*.

#### *5.2. Lid-Driven Semi-Circular Cavity Flow*

 = 1 = 0.5 = 0.25 = 1000 To further understand the performance of the present framework, the lid-driven semicircular cavity flow problem was then investigated. In a 1 × 0.5 rectangular domain, there is a semi-circular cavity with a radius of 0.5 modeled by immersed boundary method, and the lid is driven with velocity *ulid* = 1, as shown in Figure 6. Comparisons of cutting-line velocity at *x* = 0.5 and *y* = 0.25 with the benchmarking solutions of Glowinski et al. [21] and Ding et al. [22] were made with *Re* = 1000 and 100 × 50 mesh. This shows excellent agreement with Figure 7. To further show the accuracy of the present IDFC 2 framework, we also plotted the numerical results which was obtained by the fifth order upwinding schemebased DFC framework [23]. As can be seen in Figure 7, the present results match better than

the reference results. For the sake of completeness, the predicted contours of pressure and vorticity are also plotted in Figure 8 to show the applicability of the present framework.

= 1

= 0.5 = 0.25

–

= 1000

–

Δℎ <sup>−</sup> <sup>−</sup> <sup>−</sup> Δℎ <sup>−</sup> <sup>−</sup> <sup>−</sup> Δℎ <sup>−</sup> <sup>−</sup> <sup>−</sup> Δℎ <sup>−</sup> <sup>−</sup> <sup>−</sup>

**Figure 6.** Schematic of the lid-driven semi-circular cavity flow problem.

**Figure 7.** Comparisons of the velocity with the referenced solutions at the cutting lines for the lid-driven semi-circular cavity flow problem.

**Figure 8.** The contour plots for the lid-driven semi-circular cavity flow problem: (**a**) pressure; and (**b**) vorticity.

#### *5.3. Flow Past Circular Cylinders in Tandem*

200 ∆ = ∆ = /20

–

<sup>∞</sup> = 1 In this subsection, the flow past circular cylinders in tandem problem will be simulated to show that the present framework can produce accurate solutions for time transient flows. For this problem the computational domain is set as a 45 × 20 rectangular domain at the inlet with velocity *u*∞ = 1 from the left boundary, while the other three boundaries are set as convective boundary conditions. Inside the computational domain, there are two

<sup>∞</sup> = 1

= 1

= 1

=

=

200 ∆ = ∆ = /20

–

circular cylinders with diameter *d* = 1 in tandem. The centroid of the left one is located at (9.5, 0), while the right one is located at (9.5 + D, 0), where D is the distance between two circular centroids. In this study, two scenarios of D = 4 and D = 5 with *Re* = 200 and mesh size ∆*x* = ∆*y* = *d*/20 were investigated. The details of the present settings are plotted in Figure 9. The simulated vorticity and pressure contours at t = 400 are plotted in Figures 10 and 11 for D = 4 and D = 5, respectively. As shown in Figures 10 and 11, the present framework can correctly produce the vortex shedding behavior, while the pressure can remain smooth. The predicted results were then used to evaluate the drag coefficient (*CD*), lift coefficient (*CL*) and Strouhal number (*St*), and compare with the referenced solutions ([24–27]). From Table 2, good agreements can be observed between predicted results and benchmarking solutions. The drag and lift coefficient are also plotted in Figures 12 and 13 with respect to time. It is confirmed that the present framework for transient flows is applicable and accurate. 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2

1 2 1 2

 = 200 1 2

 = 200 1 2

 = 200 1 2

**Figure 9.** Schematic of the flow past circular cylinders in tandem problem.

**Figure 10.** The instantaneous contour plots for the flow past circular cylinders in tandem problem, D = 4: (**a**) vorticity; and (**b**) pressure.

**Figure 11.** The instantaneous contour plots for the flow past circular cylinders in tandem problem, D = 5: (**a**) vorticity; (**b**) pressure.


**Table 2.** Comparisons of *CD*, *C<sup>L</sup>* and *St* for the flow past circular cylinders in tandem problem with *Re* = 200. Note that *CD*1 is the left circular cylinder, while *CD*2 is the right circular cylinder.

1 2 **Figure 12.** The time history plots for the flow past circular cylinders in tandem problem, D = 4. Note that *CD*1 is the left circular cylinder, while *CD*2 is the right circular cylinder: (**a**) *CD*; and (**b**) *CL*. 1 2

1 2 **Figure 13.** The time history plots for the flow past circular cylinders in tandem problem, D = 5. Note that *CD*1 is the left circular cylinder, while *CD*2 is the right circular cylinder: (**a**) *CD*; and (**b**) *CL*.

( = −8~24, = −8~8)

( = −8~24, = −8~8)

), 0 ≤ ≤ 16

), 0 ≤ ≤ 16

), 0 ≤ ≤ 16

), 0 ≤ ≤ 16

− 16, 0 ≤ ≤ 32

− 16, 0 ≤ ≤ 32

32 − , 0 ≤ ≤ 32

32 − , 0 ≤ ≤ 32

, = 0,0

, = 0,0

() = {

() = {

4 sin ( 4

16 − 4 sin ( 4

4 sin ( 4

16 − 4 sin ( 4

() = {

() = {

−2

−2

(,

(,

∆ = ∆ = /40 ∆ = 10

∆ = ∆ = /40 ∆ = 10

= 16,1.5)

= 16,1.5)

#### *5.4. Two Circular Cylinders Moving towards Each Other in Quiescent Flow*

In the computational domain ( *x* = −8 ∼ 24, *y* = −8 ∼ 8), there are two circular cylinders with a unit diameter (d = 1), the lower initially located at (*x<sup>l</sup>* , *y<sup>l</sup>* = 0, 0), while the upper one located at (*xu*, *y<sup>u</sup>* = 16, 1.5), and then moving toward each other along the *x* direction with the prescribed motions until *t* = 32:

$$\mathbf{x}\_{l}(t) = \begin{cases} \frac{4}{\pi} \sin\left(\frac{\pi t}{4}\right), & 0 \le t \le 16\\ t - 16, & 0 \le t \le 32 \end{cases} \tag{42}$$

$$\chi\_{\iota}(t) = \begin{cases} 16 - \frac{4}{\pi} \sin\left(\frac{\pi t}{4}\right), & 0 \le t \le 16\\ 32 - t, & 0 \le t \le 32 \end{cases} \tag{43}$$

Based on the above prescribed motions, two cylinders will move periodically in the *x* axis until *t* = 16, then start to move towards each other. In this study, a no-slip boundary condition is applied to all boundaries, and the mesh size and time step size are chosen as ∆*x* = ∆*y* = *d*/40 and ∆*t* = 10 −2 . From Figures 14 and 15, it shows that the present framework can produce smooth pressure and vorticity. To validate the accuracy of the predicted results, we plotted Figure 16 to compare the time-history drag (CD) and lift (CL) coefficients of the upper cylinder with the referenced solutions [9,28] which shows excellent agreements. For the sake of completeness, we also plotted Figure 17 to show the effect of velocity modification method. This clearly shows that with the velocity modification method, SFO is well suppressed.

 **Figure 14.** The predicted pressure contour plots for the problem of two cylinders moving towards each other: (**a**) *t* = 4; (**b**) *t* = 16; (**c**) *t* = 24; and (**d**) *t* = 32.

 **Figure 15.** The predicted vorticity contour plots for the problem of two cylinders moving towards each other problem. (**a**) *t* = 4; (**b**) *t* = 16; (**c**) *t* = 24; and (**d**) *t* = 32. 

 **Figure 16.** Comparisons of the time history plots of the upper cylinder for the problem of two cylinders moving towards each other: (**a**) *CD*; and (**b**) *CL*.

 **Figure 17.** Comparisons of the effect of the velocity modification method for the problem of two cylinders moving towards each other: (**a**) *CD*; and (**b**) *CL*.

#### *5.5. Free-Falling Circular Cylinder in Quiescent Flow*

() () − The free-falling circular cylinder problem will be simulated in this subsection. In a 2 × 6 rectangular box, a circular cylinder located at (1,4) with a radius r = 0.125 is free falling from rest, as shown in Figure 18a. The densities of solid (*ρs*) and fluid *ρf* are 1.25 and 1, while the non-dimensional fluid viscosity and gravity for the investigated problem are 0.1 and (0, −980). The cylinder positions and vorticity contours at different times are plotted in Figure 18. In order to validate the present framework, comparisons of the time

∆ =

/)

−4

∆ = 1/128 ∆ = 10

−

∆ = 1/128 ∆ = 10

−4

history of cylinder positions and velocity are also plotted at Figure 19. It can be clearly shown that the preset results obtained by the proposed framework and methodology with ∆*x* = ∆*y* = 1/128 and ∆*t* = 10 <sup>−</sup><sup>4</sup> agree well with the reference benchmark results [29,30]. For the sake of completeness, CD, C<sup>L</sup> and *ω* are plotted in Figure 20a, while the effect of velocity modification method is plotted in Figure 20b. From Figure 20, it shows that the present framework can also well suppress for the free-falling cylinder problem. The pressure contours at t = 0.2 are also plotted to further show the benefit of employing velocity modification method. From Figure 21a, the smooth pressure is obtained, while Figure 21b clearly shows oscillated pressure near the downward cylinder interface. /)

() ()

∆ =

**Figure 18.** The predicted results for the free-falling circular cylinder problem: (**a**) the cylinder positions at different times; (**b**) vorticity contours at t = 0.4; (**c**) vorticity contours at t = 0.6; and (**d**) vorticity contours at t = 0.8.

**Figure 19.** Comparisons of the results for the free-falling circular cylinder problem with benchmark results [29,30]: (**a**) velocity; and (**b**) centric y location.

 **Figure 20.** The time history plots for the free-falling cylinder problem: (**a**) *CD*, *C<sup>L</sup>* and *ω*; and (**b**) comparison of the velocity modification method for *CL*.

**Figure 21.** The predicted pressure contours at t = 0.2 for the free-falling cylinder problem: (**a**) with velocity modification method; and (**b**) without velocity modification method.

To show the advantage of the present framework, different density ratios (*ρs*/*ρ<sup>f</sup>* ), 1.05 and 1.01, were also investigated in this study. Figure 22 clearly shows that the present framework can correctly and consistently predict solid velocity and position. The applicability and efficiency of the present framework for solving the free-falling cylinder problem is thus confirmed.

**Figure 22.** Comparisons of the results for the free-falling circular cylinder problem with benchmark results [25,26]: (**a**) velocity; and (**b**) centric y location.

#### *– – 5.6. Drafting–Kissing–Tumbling (DKT) Problem of Two Free-Falling Circular Cylinders in Quiescent Flow*

– – (, ) = (0~6, −1~1) − 5 × 10 −3 3Δ Δ(= Δ/) Δ Δ = Δ = 1/256 Δ = 5 × 10 −5 ≤ 0.17 > Finally, the drafting–kissing–tumbling (DKT) problem was investigated in this subsection. In a computational domain of (*x*, *y*) = ( 0 ∼ 6, −1 ∼ 1), there are two cylinders with a radius r = 0.125 located at (1, 0.001) and (1.5, −0.001). The density ratio and viscosity ratio was chosen as 1.5 and 0.01, while the gravity is (980, 0). To model the DKT scenario in this study, the repulsion force model [29] was employed in this study, with the choice of a stiffness parameter as 5 × 10 <sup>−</sup><sup>3</sup> and effective force range as 3∆*x*. Due to the large repulsion force, the sub-time stepping method is employed to prevent the divergence issue when solving the equations of motions. The idea is to introduce a sub-time step ∆*ts*(= ∆*t*/*s*), which is *s* times smaller than the original time step ∆*t*, and solves the sub-time step *s* times. In practice, the fluid solver costs the most of the computational time, so the additional cost due to the sub-time stepping can be ignored. In this study, the sub-time step *s* is set as 20. The predicted pressure and vorticity contours with <sup>∆</sup>*<sup>x</sup>* = <sup>∆</sup>*<sup>y</sup>* = 1/256 and <sup>∆</sup>*<sup>t</sup>* = <sup>5</sup> × <sup>10</sup> <sup>−</sup><sup>5</sup> at t = 0.15, 0.2, 0.25 and 0.3 are plotted in Figures 23 and 24. The comparisons of solid velocity and position with the referenced solution [2] are then made. From Figures 25 and 26, good agreements can be seen at *t* ≤ 0.17, (drafting and first kissing stage). For *t* > 0.17 (the later kissing and tumbling stage), it still shows reasonable agreements. As indicated by [2], the discrepancy should be due to strong instability and the choice of collision parameter and

0.17

*– –*

5 × 10 −5

0.17

– –

Δ/) Δ

5 × 10

the initial horizontal offset of two cylinders. It is noted that there is no drafting for the present framework if the horizontal offset is set to zero.

≤ 0.17 >

(, ) = (0~6, −1~1) −

Δ(=

−3 3Δ

Δ = Δ = 1/256 Δ =

**Figure 23.** The predicted pressure contour plots for the drafting–kissing–tumbling (DKT) problem: (**a**) *t* = 0.15; (**b**) *t* = 0.2; (**c**) *t* = 0.25; and (**d**) *t* = 0.3. – – 

 **Figure 24.** The predicted vorticity contour plots for the DKT problem: (**a**) *t* = 0.15; (**b**) *t* = 0.2; (**c**) *t* = 0.25; and (**d**) *t* = 0.3. 

– –

– –

 – – – – **Figure 25.** Comparisons of the x velocity (*uc*) and x position (*xc*) of two cylinders with referenced results for the drafting–kissing–tumbling (DKT) problems of two free-falling particles: (**a**) *uc*; and (**b**) *xc*.

– –

 – – **Figure 26.** Comparisons of y velocity (*vc*) and y position (*yc*) of two cylinders with referenced results for the drafting–kissing–tumbling (DKT) problems of two free-falling particles. (**a**) *vc*; (**b**) *yc* .

– –

#### **6. Conclusions**

In this study, we proposed a framework to model the flow–particle interaction problem. The solid object is modeled by the DIIB method, while the fluid velocity and pressure are solved by the IDFC <sup>2</sup> method. To prevent the instability problem of a density ratio close to 1, the hydrodynamic force and torque which appeared in the equations of motion are directly evaluated with modified interpolation kernel function. The methodology of evaluating the drag/lift forces and the treatment of a circular cylinder approaching the wall have also been addressed. To validate the applicability and accuracy, problems of one and two particles were investigated. To make the framework more stable, the idea of a sub-time stepping method was also employed for DKT problems of two particles. From the investigating validation/benchmarking problems, the simulated results reveal the applicability, accuracy and efficiency of the present framework.

**Author Contributions:** P.-H.C. developed the numerical model, performed the numerical simulations, analyzed the results, plotted the figures, as well as wrote the manuscript. H.C.W. analyzed the results and wrote/revised the manuscript. Y.Z.C. review and editing. R.B. review and editing. Y.-T.L. performed the numerical simulations, analyzed the results, as well as wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** Yan-Ting Lin would like to thank the Ministry of Science and Technology, Republic of China for the funding support under project No. MOST 108-3116-F-042A-006.

**Acknowledgments:** Raymond Byrne acknowledge the Research Office at Dundalk Institute of Technology, Ireland.

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

#### **References**


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

*Energies* Editorial Office E-mail: energies@mdpi.com www.mdpi.com/journal/energies

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-0365-2930-1