**2. Mathematical modeling**

The mathematical model that describes the flow consists of a set of coupled differential equations representing the physical phenomenon for which we want the solution. The literature has shown that only a fraction of practical problems can be resolved due to the complexity of the equations. Thanks to high-performance computers and numerical methods, the solution of several problems becomes possible. The formulations of the Immersed Boundary Method and Virtual Physical Model are briefly presented.

#### **2.1 Immersed boundary method**

The Immersed Boundary Method (Peskin, 1977) along with the Virtual Physical Model (Lima e Silva et al., 2007) are used to solve two-dimensional, incompressible, isothermal and transient flows. It is based on the motion equations plus a force term which model the interface. Thus, it becomes necessary to use two formulations: one for the fluid (Eulerian fixed mesh) and another for the immersed interface (Lagrangean mesh). These meshes are geometrically independent and coupled through the force term.

#### **2.1.1 Mathematical formulation for the fluid**

The Navier-Stokes equations, Eq. (1), and the continuity equation, Eq. (2), for a Newtonian fluid, are presented below using the tensor form:

$$\frac{\partial u\_i}{\partial t} + \frac{\partial \left(u\_i u\_j\right)}{\partial \mathbf{x}\_j} = -\frac{1}{\rho} \left. \frac{\partial p}{\partial \mathbf{x}\_i} + \frac{\partial}{\partial \mathbf{x}\_j} \right[ \nu \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right) \bigg] + f\_{i,j} \tag{1}$$

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

where [kg/m3] and [m2/s] are respectively the specific mass and the kinematic viscosity, properties that characterize the fluid. The variables of interest are represented by: the velocity components *ui* [m/s], the pressure *p* [N/m2] and the components of the Eulerian force acting on the interface *<sup>i</sup> f* [N/m3]. The force term only exists in Eulerian points coincident or close to the Lagrangean mesh, being zero for the remaining points of the calculation domain. This term is calculated by distributing the components of the Lagrangean interfacial force vector, using a distribution function (Peskin & McQueen, 1994):

$$f\left(\mathbf{x}\right) = \sum\_{k} D\_{ij}\left(\mathbf{x} - \mathbf{x}\_{k}\right) F\left(\mathbf{x}\_{k}\right) \Delta S^{2}\left(\mathbf{x}\_{k}\right) \tag{3}$$

where *f x* [N] is the Eulerian force vector, *x* [m] and *<sup>k</sup> x* [m] are respectively the position vectors of Eulerian and Lagrangean points, *S* [m] is the arc length centered on each Lagrangean points, *F x <sup>k</sup>* [N] is the interfacial force calculated by the IBM and *Dij* [m-2] is a interpolation/distribution function, which properties are the same of the Gaussian function.

#### **2.1.2 Mathematical formulation for the immersed interface – Virtual Physical Model (VPM)**

The VPM allows the calculation of Lagrangean force based on physical interaction of the fluid and immersed solid surface in the flow. This model is based on applying the balance of momentum quantity over the fluid particles located at the Lagrangean points. The equation that determines this force is expressed as:

$$F(\mathbf{x}\_{k'},t) = \underbrace{\rho \frac{\partial \mathbf{V}(\mathbf{x}\_{k'},t)}{\partial t}}\_{\mathbf{F}\_{\mathbf{z}}} + \underbrace{\rho \nabla \left(\mathbf{V}(\mathbf{x}\_{k'},t)\mathbf{V}(\mathbf{x}\_{k'},t)\right)}\_{\mathbf{F}\_{\mathbf{z}}} - \underbrace{\mu \nabla \left[\nabla \mathbf{V}(\mathbf{x}\_{k'},t) + \nabla^{T} \mathbf{V}(\mathbf{x}\_{k'},t)\right]}\_{\mathbf{F}\_{\mathbf{z}}} + \underbrace{\nabla p(\mathbf{x}\_{k'},t)}\_{\mathbf{F}\_{\mathbf{z}}} \tag{4}$$

where *Fa* [N] is the acceleration force, *Fi* [N] is the inertial force, *Fv* [N] is the viscous force and *Fp* [N] is the pressure force.

#### **2.2 Turbulence model**

62 Computational Simulations and Applications

The Immersed Boundary Method (IBM), due to their capability to deal with problems of complex and mobile interfaces, becomes attractive, especially in cases involving large displacements. In the modeling process of physical problems, the equations that govern the physics of the problem appear naturally. These models can range from those involving only one differential equation to those involving a system of differential equations, which can be fully coupled. However, in most cases, exact solutions can not be obtained and numerical methods appear as a tool to solve these problems. The Immersed Boundary method is used here with the Virtual Physical Model in order to simulate two-dimensional incompressible flows over stationary, rotating and rotationally-oscillating circular cylinders. Different time discretization methods are used: first order Euler scheme and the second-order Adams-Bashforth and Runge-Kutta schemes. The sub-grid Smagorinsky model and a damping function are also used. Considering the existence of a mistaken view about the mentioned numerical methods, their stability analyses are made in the present work. The results are

compared with numerical and experimental results obtained from the literature.

Immersed Boundary Method and Virtual Physical Model are briefly presented.

geometrically independent and coupled through the force term.

**2.1.1 Mathematical formulation for the fluid** 

fluid, are presented below using the tensor form:

The mathematical model that describes the flow consists of a set of coupled differential equations representing the physical phenomenon for which we want the solution. The literature has shown that only a fraction of practical problems can be resolved due to the complexity of the equations. Thanks to high-performance computers and numerical methods, the solution of several problems becomes possible. The formulations of the

The Immersed Boundary Method (Peskin, 1977) along with the Virtual Physical Model (Lima e Silva et al., 2007) are used to solve two-dimensional, incompressible, isothermal and transient flows. It is based on the motion equations plus a force term which model the interface. Thus, it becomes necessary to use two formulations: one for the fluid (Eulerian fixed mesh) and another for the immersed interface (Lagrangean mesh). These meshes are

The Navier-Stokes equations, Eq. (1), and the continuity equation, Eq. (2), for a Newtonian

*j ij ji*

0 *<sup>i</sup> i u x*

viscosity, properties that characterize the fluid. The variables of interest are represented by: the velocity components *ui* [m/s], the pressure *p* [N/m2] and the components of the Eulerian force acting on the interface *<sup>i</sup> f* [N/m3]. The force term only exists in Eulerian

[m2/s] are respectively the specific mass and the kinematic

*i*

, (2)

, (1)

*i j* 1 *<sup>j</sup> i i*

*u u u u p u + = + +f t x <sup>ρ</sup> xx xx*

**2. Mathematical modeling** 

**2.1 Immersed boundary method** 

where

[kg/m3] and

Turbulence is one of the most challenging problems of modern physics and is among the most complex and beautiful phenomena in nature. Due to several practical implications for many sectors, the number of research related to understanding and controlling these flows has increased. The turbulence effects can be modeled and simulated since emprirical correlations and diagrams up to modern methodology of numerical simulation.

#### **2.2.1 Turbulence equations**

It is known that even for flows controled by moderate Reynolds numbers, it is not possible to solve directly all frequencies present in a turbulent flow. Reynolds (1894) proposed a decomposition process of the Navier-Stokes equations in a mean and floating part in order to solve the turbulent flow. The decomposition process of the scales yielded two groups of equations for the turbulence, the first being called Reynolds Averaged Navier-Stokes equations, and another called the filtered Navier-Stokes equations (Smagorinsky, 1963). After applying the filtering and the decomposition process and applying the definitions in Eqs. (1) and (2), we obtain the following equation:

$$\frac{\partial \overline{u}\_{i}}{\partial t} + \frac{\partial \left(\overline{u}\_{i}\overline{u}\_{j}\right)}{\partial \mathbf{x}\_{j}} = -\frac{1}{\rho} \frac{\partial p^{\*}}{\partial \mathbf{x}\_{i}} + \frac{\partial}{\partial \mathbf{x}\_{j}} \left[\nu\_{\text{eff}} \left(\frac{\partial \overline{u}\_{i}}{\partial \mathbf{x}\_{j}} + \frac{\partial \overline{u}\_{j}}{\partial \mathbf{x}\_{i}}\right)\right] + f\_{i} \tag{5}$$

An Introduction of Central Difference Scheme Stability for High Reynolds Number 65

The spatial discretization is performed using the second order centered finite difference scheme and the time discretization with the first order Euler method, Adams-Bashforth and

It is presented, then, a brief description of the time discretization methods used here,

It is a first-order method for solving transient problems. With this method the time

*u u <sup>=</sup> <sup>f</sup> u ,u P +F <sup>Δ</sup><sup>t</sup>*

where *f* includes advective and diffusive terms of the motion equation. The index *n* , is related to time and *Δt* is the time step. This method is easy to implement, but requires the use of small time step to ensure the stability of the solution. The terms *<sup>n</sup> Pi* and *<sup>n</sup> Fi* are the pressure gradient and force field, respectively in the *i* direction. The term *u* is an estimate

It is a multi-point method, where the velocity fields at the current time is obtained using information from two previous time instants. In other words, the advective and diffusive terms in *n* and 1 *n* are needed for the calculations in time 1 *n +* . Multipoint methods are easy to be implemented and require only an evaluation of the derivatives by time step, making them relatively inexpensive. The main disadvantage of these methods is that as they require information on previous points, they can not be started by themselves. For this purpose, we use the Euler method for initial calculus. For *i* component of the estimated

*i j i j i i*

(9)

*i i nn n n n n*

It is a single stage method, ie, to determine *n+*<sup>1</sup> *ui* , one needs only the information available at the previous time *n n u,u <sup>i</sup> <sup>j</sup>* . In this method or in the superior orders the function in one or more additional points should be calculated. The first step until the middle of interval can be regarded as a predictor step, based on the explicit Euler method, which is accompanied by a correction to the end of the range. In summary, this method needs of information calculated only on the last time. Moreover, it requires the calculation of the function

*<sup>i</sup> <sup>j</sup> f u ,u* twice and thus consumes more time. For the *i* component of the estimated

*u u <sup>=</sup> <sup>f</sup> u ,u <sup>f</sup> u ,u P +F <sup>Δ</sup><sup>t</sup>*

<sup>1</sup> 3 1 1 1 2 2

*i j i i*

(8)

 *n+ n* <sup>1</sup> *i i nn n n*

Runge-Kutta, both of second order.

**3.2 Time discretization methods** 

**3.2.1 Euler method** 

*t* 

derivative *ui*

*n n*

Predictor step:

already making an analogy with the motion equation.

can be approximated by:

of the velocity inherent of the coupling method used.

**3.2.2 Second order Adams-Bashforth method** 

velocity, this method can be represented by:

**3.2.3 Second order Runge-Kutta method** 

velocity, this method can be represented by:

*n+ n*

$$p^\* = \overline{p} + \frac{2}{3}\rho \text{ } k \text{ and } \nu\_{\text{cf}} = \nu + \nu\_t \tag{6}$$

where *<sup>t</sup>* is the turbulent viscosity.

#### **2.2.2 Sub-grid modeling and large eddy simulation**

The sub-grid models are suitable for the calculation of turbulent viscosity. The sub-grid Smagorinky model used here is based on the assumption that the production of sub-grid turbulent stress is equal to the dissipation. The turbulent viscosity is a function of strain rate and the length scale and is expressed as:

$$\boldsymbol{\nu}\_t = \left(\mathbb{C}\_S \boldsymbol{\ell}\right)^2 \sqrt{2\overline{\mathbf{S}}\_{ij}\overline{\mathbf{S}}\_{ij}}\tag{7}$$

where is the characteristic length, *Sij* is the strain rate, *Cs* is the Smagorinsky constant. Large eddy simulation allows us to obtain three-dimensional and transient results using the motion equations, as well as to simulate flows at high Reynolds numbers with the use of refined meshes. Like any methodology, the sub-grid model has some disadvantages, such as adjusting the constant in accordance with the problem, deficiency in modeling phenomena involving energy transfer from small scales to larger scales and disability in the calculation of viscosity near the walls, which may require the use of wall laws.

#### **3. Numerical methodology**

It is important to appreciate that numerical analysis of a two-dimensional flow is possible since that determine the values of the interest variables at discrete points. The result of a discretization process are finite difference equations, which are written for each point in the domain that we want to solve. After solving these equations, the approximated solution of the problem is found. As the number of grid points becomes large, the solution of discretized equations tends to the exact solution of the corresponding differential equation.

#### **3.1 Fractional step method**

The Fractional Step Method (Chorin, 1968) with displaced meshes for the coupling between the pressure and velocity fields is used here. This arrangement allows greater facility on discretization, without the need of mean calculus, because the velocity components are positioned on the face of control volume. Moreover, this arrangement provides more stability in the pressure-velocity coupling.

As the flow is incompressible, the pressure is no longer a function of specific mass, that is constant, ie, is not a function of the thermodynamic pressure of the fluid. The Fractional Step Method is a non-iterative method, where, from the force, velocity and pressure fields of the previous iteration, we estimate the velocity components fields. With these estimated fields, we calculate the pressure correction, through the solution of a linear system, by Modified Strongly Implicit Procedure (MSI) (Schneider & Zedan, 1981). The pressure acts as a Lagrange multiplier in minimization problems. The importance of the Poisson equation for pressure correction is that it makes the connection between the equations of motion and continuity. Provides values of *p* that allow that the values of velocities components, *<sup>n</sup>* <sup>1</sup> *u*

*<sup>n</sup>* <sup>1</sup> *v* , obtained from the respective Navier-Stokes equations, satisfy the mass conservation at time 1 *n* .

The spatial discretization is performed using the second order centered finite difference scheme and the time discretization with the first order Euler method, Adams-Bashforth and Runge-Kutta, both of second order.

#### **3.2 Time discretization methods**

It is presented, then, a brief description of the time discretization methods used here, already making an analogy with the motion equation.

#### **3.2.1 Euler method**

64 Computational Simulations and Applications

 and *ef = +* 

The sub-grid models are suitable for the calculation of turbulent viscosity. The sub-grid Smagorinky model used here is based on the assumption that the production of sub-grid turbulent stress is equal to the dissipation. The turbulent viscosity is a function of strain rate

where is the characteristic length, *Sij* is the strain rate, *Cs* is the Smagorinsky constant. Large eddy simulation allows us to obtain three-dimensional and transient results using the motion equations, as well as to simulate flows at high Reynolds numbers with the use of refined meshes. Like any methodology, the sub-grid model has some disadvantages, such as adjusting the constant in accordance with the problem, deficiency in modeling phenomena involving energy transfer from small scales to larger scales and disability in the calculation

It is important to appreciate that numerical analysis of a two-dimensional flow is possible since that determine the values of the interest variables at discrete points. The result of a discretization process are finite difference equations, which are written for each point in the domain that we want to solve. After solving these equations, the approximated solution of the problem is found. As the number of grid points becomes large, the solution of discretized equations tends to the exact solution of the corresponding differential equation.

The Fractional Step Method (Chorin, 1968) with displaced meshes for the coupling between the pressure and velocity fields is used here. This arrangement allows greater facility on discretization, without the need of mean calculus, because the velocity components are positioned on the face of control volume. Moreover, this arrangement provides more

As the flow is incompressible, the pressure is no longer a function of specific mass, that is constant, ie, is not a function of the thermodynamic pressure of the fluid. The Fractional Step Method is a non-iterative method, where, from the force, velocity and pressure fields of the previous iteration, we estimate the velocity components fields. With these estimated fields, we calculate the pressure correction, through the solution of a linear system, by Modified Strongly Implicit Procedure (MSI) (Schneider & Zedan, 1981). The pressure acts as a Lagrange multiplier in minimization problems. The importance of the Poisson equation for pressure correction is that it makes the connection between the equations of motion and continuity. Provides values of *p* that allow that the values of velocities components, *<sup>n</sup>* <sup>1</sup> *u <sup>n</sup>* <sup>1</sup> *v* , obtained from the respective Navier-Stokes equations, satisfy the mass conservation

*<sup>t</sup>* (6)

<sup>2</sup> <sup>2</sup> *tS ij ij <sup>ν</sup> = C SS* (7)

\* 2 <sup>3</sup> *<sup>p</sup> p k* 

of viscosity near the walls, which may require the use of wall laws.

*<sup>t</sup>* is the turbulent viscosity.

and the length scale and is expressed as:

**3. Numerical methodology** 

**3.1 Fractional step method** 

at time 1 *n* .

stability in the pressure-velocity coupling.

**2.2.2 Sub-grid modeling and large eddy simulation** 

where 

> It is a first-order method for solving transient problems. With this method the time derivative *ui t* can be approximated by:

$$\frac{\tilde{\boldsymbol{u}}\_{i}^{n+1} - \boldsymbol{u}\_{i}^{n}}{\Delta t} = \boldsymbol{f}\left(\boldsymbol{u}\_{i}^{n}, \boldsymbol{u}\_{j}^{n}\right) - P\_{i}^{n} + F\_{i}^{n} \tag{8}$$

where *f* includes advective and diffusive terms of the motion equation. The index *n* , is related to time and *Δt* is the time step. This method is easy to implement, but requires the use of small time step to ensure the stability of the solution. The terms *<sup>n</sup> Pi* and *<sup>n</sup> Fi* are the pressure gradient and force field, respectively in the *i* direction. The term *u* is an estimate of the velocity inherent of the coupling method used.

#### **3.2.2 Second order Adams-Bashforth method**

It is a multi-point method, where the velocity fields at the current time is obtained using information from two previous time instants. In other words, the advective and diffusive terms in *n* and 1 *n* are needed for the calculations in time 1 *n +* . Multipoint methods are easy to be implemented and require only an evaluation of the derivatives by time step, making them relatively inexpensive. The main disadvantage of these methods is that as they require information on previous points, they can not be started by themselves. For this purpose, we use the Euler method for initial calculus. For *i* component of the estimated velocity, this method can be represented by:

$$\frac{\tilde{u}\_{i}^{n+1} - u\_{i}^{n}}{\Delta t} = \frac{3}{2} f\left(u\_{i}^{n}, u\_{j}^{n}\right) - \frac{1}{2} f\left(u\_{i}^{n-1}, u\_{j}^{n-1}\right) - P\_{i}^{n} + F\_{i}^{n} \tag{9}$$

#### **3.2.3 Second order Runge-Kutta method**

It is a single stage method, ie, to determine *n+*<sup>1</sup> *ui* , one needs only the information available at the previous time *n n u,u <sup>i</sup> <sup>j</sup>* . In this method or in the superior orders the function in one or more additional points should be calculated. The first step until the middle of interval can be regarded as a predictor step, based on the explicit Euler method, which is accompanied by a correction to the end of the range. In summary, this method needs of information calculated only on the last time. Moreover, it requires the calculation of the function *n n <sup>i</sup> <sup>j</sup> f u ,u* twice and thus consumes more time. For the *i* component of the estimated velocity, this method can be represented by: Predictor step:

An Introduction of Central Difference Scheme Stability for High Reynolds Number 67

For these simulations three grids are used, which are shown in the Tab. (1), along with the three time discretization schemes. It is observed through the mean values of drag coefficients (Table 2), the similarity of results when different time discretization methods were considered and the same grid refinement. Considering the various refinements, it is noted that with the coarser grid the destabilization of the flow occurs more slowly. With the grid refinement, which filters the instabilities of high frequency, the transition of the flow is faster. It is also observed that with the grid refinement from the grid 2 to grid 3, the mean values of drag coefficients are approximately the same, which leads to the independence of the results for finer mesh than 125x250. The Sthouhal number, obtained by Fast Fourier Transform (FFT) of the lift coefficient signal is also shown in Tab. (2) for Reynolds number

Grid Points number Method

1 62x124 Euler

Re=100

Method Adams Euler R-K Adams Euler R-K Adams Euler R-K

*Cd* 1.41 1.41 1.41 1.38 1.39 1.38 1.38 1.38 1.38 *St* 0.12 0.12 0.12 0.15 0.15 0.15 0.15 0.16 0.16

Note that the mean values of the drag coefficient decreases with grid refinement for the three methods. No significant difference is observed when passing from the intermediate to the most refined grid, as mentioned previously. These results are also visualized through the time evolution of the drag coefficient, Fig. (2), which presents the different grid

**4.2 Stability of the time discretization schemes increasing the Reynolds number**  For this analysis, simulations are carried out with the different time discretization methods mentioned and Reynolds numbers of 100, 300 and 1,000. For these simulations the grid is composed by 125x250 points, once, as analysed, the grid refinement did not alter significantly the inherent characteristics of the flow as the drag coefficient. Moreover, the

Grid 1 2 3

Table 2. Mean values of drag coefficients and Strouhal number for the three time

discretization methods and different grids.

cost of grid 3 is greater.

refinement for each of the time discretization methods.

2 125x250 3 250x500

Table 1. Grids used for the three time discretization schemes, Re=100.

Adams-Bashforth Runge-Kutta (R-K)

**4.1 Analyses of the grid refinement** 

100.

$$\frac{\tilde{\mu}\_{i}^{n\star}\frac{1}{2} - \mu\_{i}^{n}}{\frac{\Delta t}{2}} = f\left(\mu\_{i}^{n}, \mu\_{j}^{n}\right) - P\_{i}^{n} + F\_{i}^{n} \tag{10}$$

Corrector step:

$$\frac{\tilde{\boldsymbol{u}}\_{i}^{n+1} - \boldsymbol{u}\_{i}^{n}}{\Delta t} = f\left(\tilde{\boldsymbol{u}}\_{i}^{n+} \frac{\mathbf{1}}{2}, \tilde{\boldsymbol{u}}\_{j}^{n+} \frac{\mathbf{1}}{2}\right) - P\_{i}^{n} + F\_{i}^{n} \tag{11}$$
