2-D Iterative Scheme Pressure Correction Based Navier-Stokes Equation using Nodal Integral Method for Square Cylinder

2-D Iterative Scheme Pressure Correction Based Navier-Stokes Equation using Nodal Integral Method for Square Cylinder

1. 1 Mr. Vaibhav V. Bondge 2. Mr. Raghava Raju Khandabhattu

   HCL Technologies

   Abstract:- This paper deals with solving the Navier-Stokes equations arising during the flow across solid objects. The research work presented here deals with the numerical modeling of such flows by the application of the Nodal Integral Method (NIM) developed specifically for this purpose. The NIM-based schemes have been used to solve partial differential equations and are known to have very high accuracy compared to conventional numerical schemes. The scheme is implemented using a SIMPLE (Semi-Implicit Method for Pressure-Linked Equations)-like algorithm for pressure and velocity correction instead of solving the exact pressure Poisson equation. The objective of this paper is to model the unsteady flow of incompressible fluid across a rectangular prism using simple collocated grid arrangement for the Navier-Stokes equations. The flow is modelled numerically and flow visualizations have been generated for increasing values of Reynolds Number for a fluid at respective kinematic viscosity.

   Keywords: Navier-Stokes Equations, Nodal Integral Method, Square Cylinder, Poisson Equation, Kinematic Viscosity.

   1. INTRODUCTION

      Over the past several years many computational models have been developed in evaluating flow properties across solid objects. In current paper, Nodal Integral Method (NIM) has been employed to solve the Navier-Stokes equations arising during the flow across solid objects. As the nodal schemes are coarse-mesh methods developed to efficiently solve sets of linear and nonlinear partial differential equations (PDEs) for various applications, a significant saving in computational effort is expected without sacrificing accuracy.

      These schemes approximately satisfy the governing differential equations on finite-size bricklike elements called nodes or cells that are obtained by discretizing the space of variables. The idea of NIM is to generate the shape functions within the cells, i.e., cell analytical solution for the primitive variables of PDEs, by converting them to ordinary differential equations (ODEs).

      NIM schemes have been developed and applied to solve fluid flow and heat transfer problems. In earlier implementations of this method, schemes were developed by clubbing the nonlinear convection terms of Navier-Stokes equations into inhomogeneous terms of ODEs generated from N-S equations [2]. This leads to an inefficient generation of shape functions for velocity components that captures only the diffusion process and not convection. Moreover, use of the continuity equation to develop cell analytical solution leads to asymmetry in local solutions or shape functions of primitive variables (velocity components) in a cell [1]. The shortfalls of these schemes were addressed by the Modified Nodal Integral Method (MNIM) proposed by Rizwan-uddin [2].

      Moreover, since the MNIM is inherently implicit, a transient solution to get to the steady state requires significant effort. Therefore, in the present work a steady-state formulation of the N-S equation by a NIM approach is developed. It should be noted that the steady-state NIM formulation cannot be trivially obtained from the corresponding transient formulation.

      In the current formulation of NIM for N-S equations, instead of solving the exact Poisson equation for pressure, an approach similar to SIMPLE is used to correct the pressure at each iterative loop. Thus it takes full advantage of the SIMPLE philosophy, which is to solve for a velocity distribution using the momentum equations in such a way that it will eventually converge to the distribution which satisfies the continuity equation. This goal is achieved by having appropriate pressure correction iterations replace the exact pressure Poisson equation as used by Rizwan-uddin [2].

   2. NUMERICAL DETAILS

This section will give the details of modeling and the numerical procedure for solving the equations obtained.

1. Steps in NIM procedure

   1. Discretize the domain of interest into sub-domains called nodes or cells as shown in the Figure 1a

   2. Carry out the transverse integration procedure (TIP) for each PDE with the dependent variable (Ø). This transverse integration over the cell dimensions of all but one independent variable reduces each PDE to a set of ODEs. The dependent variables in these ODEs are referred to as transverse-averaged variables. These transverse-averaged variables are represented over the respective faces of a cell as shown in Figure 1c

      Figure 1: Discretization of domain for NIM. (a) Schematic of global (x, y) space, and its division into computational elements.

      (b) Local coordinate system with the node (i, j). (c) Transverse-averaged quantities on the surfaces in two adjacent nodes in x direction.

   3. Obtain the cell local solution of these ODEs.

   4. By using the continuity of these transverse-averaged variables (and their derivatives for second-order ODEs) on cell boundaries obtain a set of discrete equations.

   5. Apply the constraint conditions to get the full discretized equations.

2. Governing equations

Two-dimensional, steady-state, incompressible, isothermal Navier-Stokes equations in primitive variable form are

u + v = 0 (1)

x y

2 2

2 2

uu +v v ( u + u + 1 p + g (x,y) = 0 (2)

x y

x2

y2)

x x

2 2

2 2

uu +v v ( v + v

+ 1 p + g (x,y) = 0 (3)

x y

x2

y2)

y y

where gx (x, y) and gy (x, y) are body force terms such as gravity. Instead of solving the exact pressure Poisson equation, a correction formula similar to the SIMPLE algorithm is used to update the pressure along with correction of velocity by solving the momentum Equations (2)(3), respectively. The elliptical form of the pressure equation in terms of pressure correction p is given by

2 2

2 2

( p + p ) = (u + v ) (4)

x2

y2

x y

where, u* and v* are estimated velocity components in x and y directions, respectively.

2.2.1 Transverse integration procedure and the set of ODEs.

Applying the transverse integration operator, (1/2ai)

over the (i, j)th cell to Eqs. (2), (3),

and (4) yields, respectively,

1

1

1 2px(y)

y2

= (y) (5)

v ux(y) 2ux(y) = (y) (6)

c y y2 2

v vx(y) 2vx(y) = (y) (7)

c y y2 3

where, the cell-specific subscripts (i, j) on dependent variables are omitted. Terms not explicit in Eqs. (5)(7) are lumped into the right hand side as pseudo-source terms, such as

2

2

= 1 ai [ (u + v ) p]dx (8)

1 2ai

ai

x y

x2

The inhomogeneous pseudo-source terms in Eqs. (5)(7) are then expanded in Legendre polynomials and truncated at the zeroth order, which leads to consistent second-order numerical scheme. The above process reduces Eqs. (5)(7) to

1 2px(y)

y2

= s1xy (9)

v ux(y) 2ux(y) = s xy (10)

c y y2 2

v vx(y) 2vx(y) = s xy (11)

c y y2 3

Similarly, applying the transverse integral operator (1/2ai) and truncating the pseudo-source terms at the zeroth order

yields the aboe similar equations

2.2.2. Local solutions.

The set of ODEs generated using the transverse integration process are solved locally within each cell. The local solution of Eq.

(9) for x-averaged pressure update (y) is quadratic.

2

2

px (y) = (s xy y + C y + C ) (12)

1 2 1 2

Similarly, the local solution for transverse integrated velocities (y) and (y) are of constant +linear + exponential form. The local solutions of Eq. (10) and Eq. (11) are, respectively,

C vcy y C

ux(y) = 3 e( ) + s xy + 4

(13)

vc

C vcy

2 vc vc

y C

vx(y) = 5 e( ) + s xy + 6

(14)

vc 3 vc vc

where, C1C6 are constants of integration. These constants of integration are eliminated by transverse-averaged values of primitive variables at the boundary of the cell (i, j) as boundary conditions.

After eliminating the constants of integration, the resulting expressions for Eqs. (12)(14) are

px

i,j

(y) =

s1xyi,j 2

( y2 – b2i,j) +

1

2bi,j

( px

i,j

– px

i,j1

) y + pxi,j + pxi,j1

(

(

2

) (15)

{v (ux

ux

) 2b

sxy }

Revi,j ( 2 )

vc(i,j) y

sxy

( Revi,j x

x )

e

e

ux

(y) = { c(i,j) i,j i,j1 i,j 2 i,j e } e(

) + 2 i,j + e

u i,j1 u i,j

e

e

i,j

bi,j s2xy

Rev

i,j

i,j

(e + 1)

vc(i,j)

( Revi,j1)

vc(i,j)

( Revi,j1)

+

+

i,j

i,j

i,j

Rev

vc(i,j) (e 1)

(16)

{v (vx

vx

) 2b

sxy }

Revi,j ( 2 )

vc(i,j) y

sxy

( Revi,j x

x )

vx

i,j

(y) = { c(i,j) i,j i,j1 i,j 3 i,j e } e(

Revi,j

Revi,j

vc(i,j) (e 1)

) + 3 i,j + e

vc(i,j)

v i,j1 v i,j

i,j

i,j

Rev

(e 1)

bi,j s2xy

Rev

i,j

i,j

(e + 1)

+

+

i,j

i,j

i,j

Rev

vc(i,j) (e 1)

(17)

where the local or cell Reynolds number in the y direction is defined as

Rev

i,j

= 2bi,jvc(i,j)

(18)

Applying a similar process to the terms of y-averaged primitive variables in the x direction. The local Reynolds number in the x direction is defined as

Reu

i,j

= 2ai,juc(i,j)

(19)

To eliminate the psudo-source terms we apply constraint equations. Six discrete algebraic equations per cell, derived are in

terms of 12 unknowns,

(i,j), (i,j), (i,j),

(i,j),

(i,j),

(i,j), (i,j), (i,j), (i,j), (i,j), (i,j) and

1 2 3 4 5

6

6

(i,j), Thus six more equations are needed for closure. These equations are obtained by using six constraint equations. The

first three constraint equations are obtained by ensuring that the pressure Poisson and the momentum equations are satisfied over each cell in an integral sense. Thus, applying the cell-averaged operator

1 ai,j bi,j

4ai,jbi,j

ai,j

bi,j

on Eqs. (2), (3), and (4) yields the relations between pseudo-source terms. These are

(i,j) + (i,j) + f 1(i,j) = 0 (20)

1 1

(i,j) + (i,j) + f 2(i,j) = 0 (21)

2 2

(i,j) + (i,j) + f 3(i,j) = 0 (22)

3 3

respectively where,

f 1(i,j)

ui,jyui1,jy

= (

2ai,j

vi,jyvi1,jy

+ )

2bi,j

(23)

f = 1

pi,jypi1,jy

(24)

2(i,j)

( 2ai,j )

f = 1

pi,jxpi,j1x

(25)

3(i,j)

( 2b

)

i,j

These constraint equations ensure that the final numerical scheme is conservative. f1 will work as source term for the discretized pressure Poisson. f2 and f3 will work as source terms in the discretized x-momentum and y-momentum equations, respectively. The last three constraint equations are obtained by imposing the condition that the cell-averaged variables be unique, independent of the order of integration, i.e.

= 1

bi,j

dy = 1

ai,j

dx = (26)

,

2bi,j

bi,j

,

2ai,j

ai,j

,

,

= 1

bi,j

dy = 1

ai,j

dx = (27)

,

2bi,j

bi,j

,

2ai,j

ai,j

,

,

= 1

bi,j

dy = 1

ai,j

dx = (28)

,

2bi,j

bi,j

,

2ai,j

ai,j

,

,

The pseudo-source terms are eliminated by applying all six constraints equations to six discrete equations. The final set of two discrete algebraic equations for pressure is

F11 – F12 + F13 + F11 ( + ) + F15 ( + ) 1 +1 1 +1 1+1 = F16 f1(i,j) + F17 f1(i,j+1)	(29)
F21 F22 + F23 + F24 ( + ) + F25 ( + ) 1 +1 1 +1 1+1 = F26 f1(i,j) + F27 f1(i+1,j)	(30)

Similarly the two discrete x-momentum equation equations are

F32 + F33 + F34 + F35 F36 F37 ) 31 1 +1 1 1+1 +1 = F38 f2(i,j) F39 f2(i,j+1)	(31)
F41 F42 + F43 + F44 + F45 F46 F47 ) 1 +1 1 +11 +1 = F48 f2(i,j) F49 f2(i+1,j)	(32)

Due to similarity in the discretization process, two discrete y-momentum equations can be obtained by replacing variable u by v and source term f2 by f3 from the discrete x-momentum equations, i.e., Eqs. (31) and (32). Thus a set of six discrete algebraic equations per cell, two for pressure and four for momentum equations, for six transverse averaged primitive variables, is obtained. The coefficients, F11, F12, etc., are functions of ai, bj, u, Reui,j, and Revi,j. Note that the corresponding transient formulation has eight equations per cell (six for momentum and two for pressure).

1. Solver Details

   1. Algorithm for iterative procedure

The steps involved in the solution procedure are as follows.

1. Noting that the grid is collocated (as shown in Figure 1), guess the value of the transverse average pressures (), () at

   appropriate grid surfaces. Similarly, set the values of transverse average velocities [ (), () appropriate grid surfaces. Superscript n is the sequence number of the nonlinear iterative loop.

   () and

   ()] at

2. For given nth-level pressure (), () the discretized momentum equations are solved for transverse-average velocities using current pressure in their source term. The discretized pressure corrections for (px, py) are solved in the same loop after

   the momentum equations using the current available components of velocity in their source terms. The above equations are solved using the Gauss-Seidel method. Thus, both velocity and pressure are corrected one after the other in the internal loop.

3. The criterion used for convergence is based on mass residue term f1, which is also a source term for the pressure-correction equations. Convergence is achieved when f1 (which is also the nodal discretized form of the continuity equation) approaches a defined tolerance limit.

4. If convergenc is not achieved, the velocity and pressure are updated as follows:

   (u)n+1 = (1 v) (u)n + v (u)n+1	(30)
   (v)n+1 = (1 v) (v)n + v (v)n+1	(31)
   (p)n+1 = (p)n + p (p)n+1		(32)

   where, v and p are relaxation factors for velocity and pressure, respectively. From extensive numerical experimentation it was found that for low nonlinearity in the momentum equation, i.e., for low Reynolds (or Rayleigh) number, the values of v, p could be chosen between 0.8 and 1.0, but for high Reynolds (or Rayleigh) number, v is as low as 0.5 and p could be 0.5.

5. Repeat steps 24 until convergence is achieved. The error tolerance limit () used for the convergence criteria is based on the root-mean-squared value of mass residue f1 in an absolute sense.

   2.3.2 Boundary conditions

   In most cases free slip boundary condition is commonly used which involves vortex shedding because the cross-stream outlet boundary conditions is a very important issue. An outflow boundary (i) permits the flow to exit the domain with a smooth discharge of vortices, (ii) has a minimum effect on the flow inside the domain, and (iii) has a negligible effect on the near body flow. For finite difference and finite volume discretization, Neumann boundary condition (NBC) and convective boundary condition (CBC) are the two most important boundary conditions but here we have put up the Neumann boundary condition at the top and bottom boundary as well as outlet of the computational domain which include the gradient of the flow variable while Dirichlet boundary condition applied at the inlet where value of flow variables are known.

   Figure 2: Boundary conditions at inlet, outlet, top and bottom boundary of the computational domain including square cylinder placed in domain.

   In the above Figure 2 a solid square cylinder is placed inside the unconfined boundary (H = L = 1, B = 0.2) having constant inlet velocity of unit magnitude in x direction while magnitude of velocity component in y direction at inlet is zero. At top and bottom boundary U-velocity, V- velocity and pressure gradients are equal to zero. At the outlet of computational domain, velocity and pressure gradient in x direction are equal to zero.

   1. RESULTS AND DISCUSSION

      A square cylinder model was developed by considering a fluid of various kinematic viscosity and then it is compared with numerical models developed by Norberg [4] and Robichaux [5] as it was not possible to get the values from the paper, contour plots are been compared and are similar to the obtained result. Various results are obtained for different Reynolds number using MATLAB program. The figures below show the variation in vortex shedding characteristics with increasing values of the Reynolds number.

      Re=20 Re=100

      Re=200 Re=1000

      Figure 3: Contours obtained for various Reynolds numbers

      The above results shown in Figure 3 are displayed at steady state. Initially, the flow is highly laminar and the flow does not show any separation. Eventually as the Reynolds number reaches Re=100, eddy formation begins on the downstream side. As Reynolds number reaches a critical value, flow separation begins and tendency of vortex formation increases. Thus, for high kinematic viscosity, the tendency of the fluid to stick to the no-slip surfaces increases. Hence, a clear flow separation does not occur till higher values of Reynolds number are attained.

   2. CONCLUSION

From the above studies, we have concluded that the flow separation and vortex shedding phenomena are largely governed by the kinematic viscosities of the fluid at particular value of Reynolds number. Also, for a fluid, formation of eddies increases with increase in Reynolds number. Further investigations may be conducted into investigating such flows across series of square cylinder placed or parallel square cylinder placed, changing the boundary conditions, arbitrary shape to evaluate flow pattern characteristics across such obstacles.

Nomenclature

a, b = Half of the size of the discretized cell or node in the x and y directions, respectively

A = Coefficient in discretized coupling equation for transverse-averaged velocity components B = Term defined in expanded form of coefficients equations.

f = Source term for discretized pressure Poisson and momentum equations F = Ccoefficients for discertized pressure Poisson and momentum equations g = Volumetric body force term for momentum equations

n = Number of cells in domain of interest p = Pressure correction transerve averaged

Re, Ra = Reynolds and Rayleigh number, respectively

Reu, Rev = Cell Reynolds number in x and y directions, respectively S = Pseudo-source terms

u = Velocity component in x direction

u , v = Transverse-averaged u and v velocity components, respectively uc, vc = Cell-averaged u and v velocity components, respectively

u ,v = Estimated velocity components in x and y directions, respectively

v = Velocity component in y direction = underrelaxation factor

= thermal diffusivity of fluid

= Coefficient of thermal expansion = Kinematic viscosity

Ø = Arbitrary dependent variable = Density

Subscripts

c = Cell-averaged quantity

i, j = Index in x and y directions, respectively

Superscripts

x = x-averaged quantity y = y-averaged quantity

xy = cell-averaged quantity

Appendix

F11 = [ 1 3b ]

2b 2(a2+b2) i,j+1

F13 = [ 1 3b ]

2b 2(a2+b2) i,j

F12 = {[ 1 + 3b ] + [ 1 + 3b ] }

2b 2(a2+b2) i,j 2b 2(a2+b2) i,j+1

F14 = [ 3b ] , F15 = [ 3b ]

2(a2+b2) i,j 2(a2+b2) i,j+1

2

2

2

2

F = [ a b ] , F = [ b a ]

16 (a2+b2)

17

i,j

(a2+b2)

i,j

F31 = [ 1 ( Rev

) + 2b ( 1

1 ) 1 ( 1 1 )]

2b eRev1

eRev1

Rev

B Rev

eRev1

i,j+1

F33 = [ 1 ( Rev

) + 2b ( 1

1 ) 1 ( 1 1 )]

2b 1 eRev

1 eRev

Rev

B Rev

1 eRev

i,j

F32 = [ 1 ( Rev

) + 2b ( 1

1 ) 1 ( 1 1 )]

+ [ 1 ( Rev

) + 2b ( 1

1 ) 1 ( 1

2b eRev1

1 )]

eRev1

Rev

B Rev

eRev1

i,j+1

2b 1 eRev

1 eRev

Rev

B Rev

1 eRev

i,j

F34 = [ 2b ( 1

1 ) 1 ( 1

1 )]

1 eRev

Rev B

1 eReu

Reu

i,j

F35 = [ 2b ( 1 1 ) 1 ( 1 1 )]

1 eRev

Rev

B Reu

1 eReu

i,j

F36 = [ 2b ( 1

1 ) 1 ( 1

1 )]

1 eRev

Rev B

1 eReu

Reu

i,j+1

F37 = [ 2b ( 1 1 ) 1 ( 1 1 )]

1 eRev

Rev

B Reu

1 eReu

i,j+1

F38 = [ 2b ( 1 1 ) 1 E]

1 eRev

Rev B

i,j

F39 = [ 2b ( 1 1 ) 1 E]

Reu

Reu

2

2

Rev

Rev

1 eRev

Rev B

i,j+1

2

2

B = {2a [ 1

(e + 1) 2

] + 2b [ 1

(e + 1) 2 ] }

Reu

Reu

Reu

eReu 1

Reu2

Rev

eRev 1

Rev2

i,j

E =

E =

2a2 1

{ [

(e + 1) 2 ]}

Reu

eReu 1

Reu2

i,j

Since the formulation is symmetric in nature, the expressions for coefficient from F41 to F49 could be generated by replacing a by b, Reu by Rev, i by j, and vice-verse.

REFERENCES

1. Neeraj Kumar, Suneet Singh and Doshi J. B., Iterative scheme for Navier-Stokes equations using NIM, Numerical Heat Transfer B, vol. 62, 264288, 2012.

2. Toreja A. J. and Rizwan-uddin, Hybrid numerical methods for convection-diffusion problems in arbitrary geometries, Comput. Fluids, vol. 32, no. 6, pp. 835872, 2003.

3. Patankar S. V. Numerical Heat Transfer and Fluid Flow, Taylor and Francis Publishers 2005, p. 126-130.

4. Sohankar A., Norberg C. Davidson L., Numerical simulation of unsteady low Reynolds number flow around rectangular cylinders at incidence, J. Wind Engg. Ind. Aerodyn., vol. 69, p. 189-201, 1997.

5. Robichaux J. , Balachandar S. and Vanka S. P. Two-dimensional Floquet instability of the wake of square cylinder, Phys. Fluids, vol. 11 p. 560- 578, 1999.