An Introduction of Central Difference Scheme Stability for High Reynolds Number

The purpose of this book is to introduce researchers and graduate students to a broad range of applications of computational simulations, with a particular emphasis on those involving computational fluid dynamics (CFD) simulations. The book is divided into three parts: Part I covers some basic research topics and development in numerical algorithms for CFD simulations, including Reynolds stress transport modeling, central difference schemes for convection-diffusion equations, and flow simulations involving simple geometries such as a flat plate or a vertical channel. Part II covers a variety of important applications in which CFD simulations play a crucial role, including combustion process and automobile engine design, fluid heat exchange, airborne contaminant dispersion over buildings and atmospheric flow around a re-entry capsule, gas-solid two phase flow in long pipes, free surface flow around a ship hull, and hydrodynamic analysis of electrochemical cells. Part III covers applications of non-CFD based computational simulations, including atmospheric optical communications, climate system simulations, porous media flow, combustion, solidification, and sound field simulations for optimal acoustic effects.

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.

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.

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.

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: is the viscous force and p F [N] is the pressure force.

Turbulence model
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.

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: where t  is the turbulent viscosity.

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: where  is the characteristic length, ij S is the strain rate, s C 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.

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.

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, 1 n u  1 n v  , obtained from the respective Navier-Stokes equations, satisfy the mass conservation at time 1 n  .

www.intechopen.com
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.

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.

Euler method
It is a first-order method for solving transient problems. With this method the time can be approximated by: 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 n i P and n i F 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.

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:

Second order Runge-Kutta method
It is a single stage method, ie, to determine , one needs only the information available at the previous time nn i j u,u . 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

Problem description
Stability analysis of the second order spatial centered scheme with the time discretization schemes is performed by two-dimensional simulations of incompressible flows around a stationary circular cylinder. The rectangular domain is chosen to be d 15 x d 30 with the cylinder located at 5 . 16 cylinder diameters from the inlet as illustrated in Fig. (1). The time step used in all simulations is 0.0001 s.

Analyses of the grid refinement
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)  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 refinement for each of the time discretization methods.

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.
www.intechopen.com 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 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). 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 www.intechopen.com 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.

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.  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 simulations at high Reynolds number since the turbulence modeling and the damping function is also applied to ensure the stability of the methodology.

Applications of the immersed boundary methodology
Firstly, are presented, simulations'results of flow over a pair of cylinders of the same diameter, following by the results of rotating and rotationally-oscillating cylinder.

Flow around two circular cylinders in tandem configuration
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 development and vortex dynamic (Akbari & Price, 2005). 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 diameter (P/D) (Sumner et al., 2005).
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.25d in the x direction and 2d in the y direction as shown in Fig. (7). The non-uniform grid www.intechopen.com region is composed by 600x300 points, the Reynolds number equal to 39,500 and the pitch ratio equal to L/D=2. Fig. 7. Illustrative scheme of the distance from the cylinder surface to the uniform region boundaries. Figure 8 shows the flow visualization through the instantaneous vorticity field after the flow has reached steady state. It is noted that the shear layers originated from the surface of the upstream cylinder surrounding the downstream cylinder, forming a single wake behind the cylinder B. It is also noted, that the vortex wake oscillates around the symmetry line of the domain. The interaction between the two layers occurs only behind the downstream cylinder, which is within the wake of the upstream cylinder. For this case, the '2S' vortex shedding mode compose the wake. It is important to appreciate that for this pitch ratio and geometrical configuration, the two cylinders behave as a single body. According to Naudascher & Rockwell (1994) no detectable vortex shedding behind the upstream cylinder occur, for L/D<3.8. Also according to these authors, as the spacing between the cylinders increases, vortex shedding occur in the upstream cylinder with a frequency that gradually approaches to the frequency for a stationary cylinder. Deng et al. (2006), in they work at low Reynolds number (Re=220), concluded that for two-dimensional simulations, each cylinder will produce a vortex wake only for L/D  4.0, with no vortex shedding between the cylinders for L/D  3.5. They also affirmed that even in threedimensional flows, for this configuration and L/D  3.5, the flow is equal to the twodimensional.   Fig. (9b) shows the time evolution of the lift coefficients. It is verified that the drag coefficient on the cylinder B is considerably smaller than the cylinder A, with mean close to zero. This can be understood by the fact that the cylinder B is inside of the upstream cylinder wake. The fluctuations of the lift coefficient of the two cylinders have zero mean, as shown in Fig. (9b). The amplitude obtained for the cylinder B is approximately seven times greater than the amplitude of cylinder A. The absence of vortices behind the upstream cylinder minimizes the lift fluctuations. Note also that the both fluctuations are in phase, Fig. (9b). This is consistent, once the vortices that are formed and transported induce forces on both cylinders simultaneously.

Flow around a rotating cylinder
The flow dynamics around a rotating cylinder is different from that observed for a stationary cylinder. The rotation of a cylinder in a uniform viscous flow modifies the vortices configuration and probably has an effect on flow-induced oscillations. As the cylinder rotates the flow is accelerated in one side and decelerated in the other side. This can be attributed to viscous effects injected by the cylinder on the flow. Therefore, the pressure on the accelerated side becomes smaller than the pressure at the decelerated side resulting on a lift force, transverse to the flow. In recent years more attention has been given to control the wake formed behind the cylinder, especially in order to suppress the vortices with the use of active or passive controls. The rotating motion of an immersed body can suppress partially or totally the vortex shedding process, so that the wake separation on one side of the body, be displaced from the axis of vertical symmetry.

Comparison of results
Aiming to compare the present results with the literature, simulations were carried out at low Reynolds numbers, which are 60, 100 and 200. For this simulations, the grid is composed by 200x125 points, refined over the cylinder (twenty grids per diameter) to ensure good accuracy in the results. The rotating moviment is imposed clockwise around its axis and is achieved by the imposition of the velocity components in each Lagrangean point. Figure 10 shows the amplitude of the drag and lift coefficients in function of the specific rotation  (the ratio of the tangential velocity and free-stream velocity) compared with the numerical results of Kang et al. (1999), for Re=60 and Re=100. Note that the drag coefficient amplitude, Fig. (10a) increases until given  and then decreases, reaching a near-zero amplitude. Note also that the amplitudes increase with the Reynolds number and the rotation in which the amplitude decreases is different for each Reynolds number. For Re=60, the amplitude of the drag is reduced for >1.0 and for Re=100 and Re=200, the reduction occur for >1.5. On the other hand, the amplitude values of the lift coefficient, Fig. (10b), shows small variations for   1.0, for all Reynolds numbers and then decreases, tending to zero. As observed, there was good agreement between the present results with those of Kang et al. (1999).

Flow over a rotationally-oscillating circular cylinder
For the stationary cylinder at low Reynolds numbers, it is known that the vortex wake is aligned and symmetrical about the central axis of the flow. The behavior is not verified when the cylinder is subjected to rotationally-oscillating moviment around its own axis. The mutual interaction between cylinder moviment and the adjacent fluid modifies the pattern wake of the flow through the acceleration and deceleration of the flow around the cylinder. Thus, there is a transition between different vortex shedding modes as the relationship between oscillation frequency and the vortex shedding frequency for the stationary cylinder varies for the same amplitude A . Commonly, some authors present two different flow regimes, being the no lock-in regime and the lock-in regime (Cheng et al. 2001a(Cheng et al. , 2001b. According to Löhner & Tuszynski (1998), the flow around a rotationally-oscillating cylinder is a forced oscillator form, or a nonlinear system, that in some cases, can become chaotic.
Here, the rotationally-oscillating cylinder is started impulsively from rest and the tangential velocity on the cylinder is given by the expression: where  is the angular velocity, A is the oscillation amplitude, R is the cylinder radius, f c is the oscillation or imposed frequency and t is the physical time. The simulations were performed for Reynolds number equal to 1,000, the non-uniform grid is composed by 400x125 points and the turbulence model and damping function in the outlet of the domain were applied.

Different vortex shedding modes
In Fig. (11) the flow visualizations are presented, through the instantaneous vorticity fields for the dimensionless time equal to 200, at different amplitude values and frequency ratios. Figure (11a) corresponds to the stationary cylinder. Figures (11b) and (11c) Fig. (11a), corresponding to the stationary cylinder, as already mentioned, there is the classical Von Kárman Street, represented by the classical '2S' vortex shedding mode. This mode indicates the generation of a positive vortex in one side of the cylinder and a negative vortex on the other side, at each oscillation cycle, forming a single vortex wake with displaced vortices around the symmetry line of the flow. In Fig. (11b), /1 . 0 5 co ff  , the vortex wake is similar to pattern wake ('2S' mode), however, the vortices are presented more rounded and with smaller longitudinal and transversal spacing between them when compared with Fig. (11a). Increasing the frequency ratio to /2 . 5 co ff  and keeping the amplitude 1 A  , Fig. (11c), there is a new vortex shedding mode called 'P+S'. This mode corresponds to a pair of vortices and single vortex composing the wake. Pairs of vortices having opposite signs are located at the inferior side of the central line of the flow, while the single vortices are released at the superior side of the cylinder. For /0 . 5 co ff  and A=2 it is also observed a new vortex shedding mode called '2P ', which corresponds to pairs of vortices of opposite signs along the wake. Keeping the same oscillation amplitude and increasing the frequency ratio to /2 . 5 co ff  , Fig. (11e), it is noted the same vortex shedding mode of the previous case, Fig (11d). Interesting to note, in this case, that the pairs are disposed symmetrically about the centerline of the flow forming a cone-shaped wake. Increasing the amplitude to A=3, and taking /0 . 5 co ff  again, a new vortex shedding mode is obtained, called '2C ', as quoted in Williamson & Jauvtis (2004). It is noteworthy that the '2C' mode is not taken by other authors for the case of circular cylinder in rotationallyoscillating moviment. According to Williamson & Jauvtis (2004) this mode was obtained for pivoted cylinders. For /2 . 5 co ff  , Fig. (11g), there is a new standard of vortex emission, in which the double vortex wake near the cylinder, composed by vortices of the same sign in each row, after a given distance away from the cylinder are coupled to form a single wake. The double wake length decreases with increasing the frequency ratio. In Fig. (11h), corresponding to /6 . 0 co ff  the instability caused by the cylinder oscillation is limited to a region near the cylinder, while far from the immersed body, the vortices reorient themselves to form the stable Von Kármán Street. Occurs, therefore, a vortex-vortex interaction of the   The concomitant use of second order temporal schemes with the spatial centered scheme is crucial for the stability of the methodology. The Adams-Bashforth temporal scheme presented more stable than the second order Runge-Kutta scheme. As the Reynolds number is increased the methodology showed to be unstable for all second-order temporal discretization schemes. This result is expected once the centered scheme has no numerical diffusion. Thus, it is concluded that for high Reynolds number, the use of turbulence modeling for the energy transfer process between the largest and smallest scales of turbulence is needed. It is important to appreciate that without the modeling and numerical diffusion the kinetic energy of the physical instabilities accumulates on the cutoff frequency and the simulation diverges. The cutoff frequency is determined by the mesh discretization.

St Pow er Spect rum
It is known that the use of developed flow boundary condition at the outlet of the domain is common in the literature. However when there are physical instabilities, which must leave the domain, there may be problems of numerical stability, especially when using centered spatial schemes. This is due to the fact that the physical instabilities carry spurious information from the outside of the domain to inside. The result is also the divergence of the simulations. To solve this problem the use of a damping function is essential to ensure stability for higher values of Reynolds number. Aiming to illustrate the applicability of the Immersed Boundary method used togheter the second order spatial centered scheme and second order temporal discretization scheme, simulations were carried out with a circular cylinder pairs, rotating cylinder and rotationally-oscillating cylinder. For the rotating cylinder case, the results showed good agreement with literature data. It was found that the rotation has greater influence on the amplitude of the drag coefficient than on the amplitude of the lift coefficient. It's worth noting that with increasing rotation the amplitude of the dynamic coefficients tends to null, as expected, once the vortex shedding process decreases. For simulations with rotationally oscillating cylinder is analyzed the influence of amplitude and frequency ratio in the vortex shedding modes, as well as in the vortex shedding frequency. It is observed different vortex shedding modes when fixed the oscillation amplitude and varies the frequency ratios. It is important to appreciate the 2C mode obtained in this study once this mode is not found in the literature for rotatinally-oscillating cylinder and it is worth mentioning that, according to Williamson & Jauvtis (2004) the 2C mode is obtained for pivoted cylinder. It is also obtained for the amplitude and frequency ratios considered the lock-in regime, whose range increases as the oscillation amplitude increases.

Acknowledgments
The authors are deeply gratefully to the following organizations: Minas Gerais State Agency FAPEMIG for the continued support to their research work, especially through the postdoctorate scholarship granted to A.R. da Silva; Brazilian Research Council -CNPq for the financial support to their research activities; CAPES Foundation from the Brazilian Ministry of Education and the School of Mechanical Engineering College of the Federal University of Uberlândia, Brazil.