**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 cost of grid 3 is greater.

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

the coupling of the spatial centered scheme with a second order temporal scheme makes the

With further increase of Reynolds number for 1,000, there is an increase in the numerical instabilities in the three methods, being more pronounced in the Euler and Runge-Kutta methods. These instabilities were already expected once a turbulence model is not being used. Being the spatial discretization scheme, centered and without numerical diffusion, it is natural that the calculation becomes unstable, leading to divergence as seen through the

(a) (b)

Fig. 3. Time evolution of the drag (*Cd*) (a) and lift (*Cl*) (b) coefficients at Re=300, for the three

(a) (b)

Fig. 4. Time evolution of the drag (*Cd* ) (a) and lift (*Cl*) (b) coefficients, at Re=1,000, for the

It is important to appreciate that for a high Reynolds number the turbulence model is needed to ensure that the kinetic energy of turbulence is carried by the wave number or

method more stable (Ferziger, 2002).

time discretization methods.

three time discretization methods.

time evolution of the dynamic coefficients in Fig. (4).

Fig. 2. Time evolution of the drag coefficient for the three grid refinement, at Re=100. Euler (a), Adams-Bashforth (b) and Runge-Kutta (c).

For Re=100, it is noted that the results are identical both qualitatively and quantitatively for the three time discretization methods (Fig. (2)). Again it is illustrated that the transient regime, with instabilities, appears later for the coarse grid (grid 1). For this grid, the start time of the instabilities formation is 150, while for the grids 2 and 3, this time is 75. Another interesting fact is that the drag coefficient oscillations is more pronounced for grids 2 and 3. This is due to the fact that the vortices are formed closer to the cylinder.

Increasing the Reynolds number from 100 to 300, it is found that the flow becomes more unstable, appearing instabilities and the drag coefficient decreases. Such instabilities can be related to the centered scheme of spatial discretization, where for this Reynolds number, the nonlinear effects become important. Figure 3 shows the time evolution of the drag and lift coefficients for the Euler, Adams-Bashforth and Runge-Kutta methods. There is a small difference in the results obtained with Euler's method when compared with the others two. When observe the lift coefficient in Fig. (3b) we see that the oscillations amplitude for Euler's method is larger than the amplitude of the signal for the others methods. This shows that

(a) (b)

(c) Fig. 2. Time evolution of the drag coefficient for the three grid refinement, at Re=100. Euler

For Re=100, it is noted that the results are identical both qualitatively and quantitatively for the three time discretization methods (Fig. (2)). Again it is illustrated that the transient regime, with instabilities, appears later for the coarse grid (grid 1). For this grid, the start time of the instabilities formation is 150, while for the grids 2 and 3, this time is 75. Another interesting fact is that the drag coefficient oscillations is more pronounced for grids 2 and 3.

Increasing the Reynolds number from 100 to 300, it is found that the flow becomes more unstable, appearing instabilities and the drag coefficient decreases. Such instabilities can be related to the centered scheme of spatial discretization, where for this Reynolds number, the nonlinear effects become important. Figure 3 shows the time evolution of the drag and lift coefficients for the Euler, Adams-Bashforth and Runge-Kutta methods. There is a small difference in the results obtained with Euler's method when compared with the others two. When observe the lift coefficient in Fig. (3b) we see that the oscillations amplitude for Euler's method is larger than the amplitude of the signal for the others methods. This shows that

This is due to the fact that the vortices are formed closer to the cylinder.

(a), Adams-Bashforth (b) and Runge-Kutta (c).

the coupling of the spatial centered scheme with a second order temporal scheme makes the method more stable (Ferziger, 2002).

With further increase of Reynolds number for 1,000, there is an increase in the numerical instabilities in the three methods, being more pronounced in the Euler and Runge-Kutta methods. These instabilities were already expected once a turbulence model is not being used. Being the spatial discretization scheme, centered and without numerical diffusion, it is natural that the calculation becomes unstable, leading to divergence as seen through the time evolution of the dynamic coefficients in Fig. (4).

Fig. 3. Time evolution of the drag (*Cd*) (a) and lift (*Cl*) (b) coefficients at Re=300, for the three time discretization methods.

Fig. 4. Time evolution of the drag (*Cd* ) (a) and lift (*Cl*) (b) coefficients, at Re=1,000, for the three time discretization methods.

It is important to appreciate that for a high Reynolds number the turbulence model is needed to ensure that the kinetic energy of turbulence is carried by the wave number or

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

simulations at high Reynolds number since the turbulence modeling and the damping

(a) (b) Fig. 6. Instantaneous vorticity fields for Re=10,000, Adams-Bashforth; without damping

Firstly, are presented, simulations'results of flow over a pair of cylinders of the same

One of the main applications of this type of study is to obtain a better understanding of the flow around a set of risers, which is subject to shear flows by ocean currents. The flow interference over bluff bodies is responsible for changes in characteristics of the fluid load that acts on immersed bodies. Consequently, the study of cylinders pair even in twodimensional simulations has received considerable attention both from the standpoint of academic and industrial fields. In addition, flow over circular cylinders involve different fundamentals dynamic phenomena, such as boundary layer separation, shear layer

The configurations with a cylinders pair in tandem and side by side are the most extensively discussed in the literature (Sumner et al., 1999; Deng et al., 2006; Silva et al., 2009), although the form more general is the staggered configuration ( Sumner et al., 2008; Sumner et al., 2005; Silva et al., 2009). According to the literature, there are various interference regimes, which were based on flow visualization in experiments. The wake behavior of two cylinders can be classified into groups according to the pitch ratio between the cylinders centers by

Here, the two cylinders have equal diameters *d* and the distance center to center of the cylinders, is called *L*. The cylinder *A* is located upstream and cylinder *B* is located downstream of the inlet. In all simulated cases, the two cylinders are disposed such that the minimum distance from the surface of each cylinder to the end of the uniform grid region is 1.25*d* in the x direction and 2*d* in the *y* direction as shown in Fig. (7). The non-uniform grid

diameter, following by the results of rotating and rotationally-oscillating cylinder.

function is also applied to ensure the stability of the methodology.

**5. Applications of the immersed boundary methodology** 

**5.1 Flow around two circular cylinders in tandem configuration** 

development and vortex dynamic (Akbari & Price, 2005).

diameter (*P/D*) (Sumner et al., 2005).

function (a), with damping function (b).

cutoff frequency. The apparent convergence given by the Adams-Bashforth method can be misleading as will be seen in the item 4.3. It is noteworthy that the spatial centered schemes have no numerical viscosity, as in the case of upwind schemes, which are stable without turbulence model, even at high Reynolds numbers. The following are presented the simulations results with sub-grid Smagorinsky modeling, needed to ensure the stability of the methodology as previously commented.

### **4.3 Simulations with the sub-grid Smagorinsky modeling**

The motion equations are sufficient to model flows in any regime and for any value of Reynolds number. However, as the Reynolds number increases the energy spectrum associated with the flow becomes wider, making it necessary the use of grid extremely fine, which implies high computational costs. Thus, with the use of coarse grids, the grid filtering process will eliminate all high frequencies providing only the low frequencies, hence the restriction on its use, without additional turbulence modeling. It is observed in Fig. (5) that even for the most stable method, Adams-Bashforth, for high Reynolds number (Re=10,000) the calculation diverges without turbulence modeling.

Fig. 5. Time evolution of the drag coefficient (*Cd*) (a) and of the lift coefficient (*Cl*) (b); both without and with turbulence model, Re=10,000.

#### **4.4 Simulations with damping function**

Figure 6 shows the flow visualization, through the instantaneous vorticity fields, for the grid of 250x500 points, using the damping function in the outlet of the domain. It is noted for the case without damping function, Fig. (6a) that the wake vortex presents an unusual behavior for two-dimensional flows, which can lead to divergence of the calculations. With the damping function, Fig. (6b), the calculation becomes more stable even at greater times of simulations. The damping function aims to remove the vortices in the outlet of the calculation domain, thus enabling the application of the boundary condition of the developed flow. This function eliminates the input of mass at the domain outlet that occur due to the vortices rotation. As verified in the presented results the second order spatial centered scheme with the second order time discretization scheme may be perfectly used for

cutoff frequency. The apparent convergence given by the Adams-Bashforth method can be misleading as will be seen in the item 4.3. It is noteworthy that the spatial centered schemes have no numerical viscosity, as in the case of upwind schemes, which are stable without turbulence model, even at high Reynolds numbers. The following are presented the simulations results with sub-grid Smagorinsky modeling, needed to ensure the stability of

The motion equations are sufficient to model flows in any regime and for any value of Reynolds number. However, as the Reynolds number increases the energy spectrum associated with the flow becomes wider, making it necessary the use of grid extremely fine, which implies high computational costs. Thus, with the use of coarse grids, the grid filtering process will eliminate all high frequencies providing only the low frequencies, hence the restriction on its use, without additional turbulence modeling. It is observed in Fig. (5) that even for the most stable method, Adams-Bashforth, for high Reynolds number (Re=10,000)

(a) (b)

Fig. 5. Time evolution of the drag coefficient (*Cd*) (a) and of the lift coefficient (*Cl*) (b); both

Figure 6 shows the flow visualization, through the instantaneous vorticity fields, for the grid of 250x500 points, using the damping function in the outlet of the domain. It is noted for the case without damping function, Fig. (6a) that the wake vortex presents an unusual behavior for two-dimensional flows, which can lead to divergence of the calculations. With the damping function, Fig. (6b), the calculation becomes more stable even at greater times of simulations. The damping function aims to remove the vortices in the outlet of the calculation domain, thus enabling the application of the boundary condition of the developed flow. This function eliminates the input of mass at the domain outlet that occur due to the vortices rotation. As verified in the presented results the second order spatial centered scheme with the second order time discretization scheme may be perfectly used for

the methodology as previously commented.

**4.3 Simulations with the sub-grid Smagorinsky modeling** 

the calculation diverges without turbulence modeling.

without and with turbulence model, Re=10,000.

**4.4 Simulations with damping function** 

simulations at high Reynolds number since the turbulence modeling and the damping function is also applied to ensure the stability of the methodology.

Fig. 6. Instantaneous vorticity fields for Re=10,000, Adams-Bashforth; without damping function (a), with damping function (b).
