 Open Access
 Authors : Vaibhav V. Bondge , Raghava Raju Khandabhattu
 Paper ID : IJERTV10IS080226
 Volume & Issue : Volume 10, Issue 08 (August 2021)
 Published (First Online): 08092021
 ISSN (Online) : 22780181
 Publisher Name : IJERT
 License: This work is licensed under a Creative Commons Attribution 4.0 International License
2D Iterative Scheme Pressure Correction Based NavierStokes Equation using Nodal Integral Method for Square Cylinder

1 Mr. Vaibhav V. Bondge 2. Mr. Raghava Raju Khandabhattu
HCL Technologies
Abstract: This paper deals with solving the NavierStokes 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 NIMbased 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 (SemiImplicit Method for PressureLinked 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 NavierStokes 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: NavierStokes Equations, Nodal Integral Method, Square Cylinder, Poisson Equation, Kinematic Viscosity.

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 NavierStokes equations arising during the flow across solid objects. As the nodal schemes are coarsemesh 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 finitesize 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 NavierStokes equations into inhomogeneous terms of ODEs generated from NS 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 Rizwanuddin [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 steadystate formulation of the NS equation by a NIM approach is developed. It should be noted that the steadystate NIM formulation cannot be trivially obtained from the corresponding transient formulation.
In the current formulation of NIM for NS 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 Rizwanuddin [2].

NUMERICAL DETAILS

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

Steps in NIM procedure

Discretize the domain of interest into subdomains called nodes or cells as shown in the Figure 1a

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 transverseaveraged variables. These transverseaveraged 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) Transverseaveraged quantities on the surfaces in two adjacent nodes in x direction.

Obtain the cell local solution of these ODEs.

By using the continuity of these transverseaveraged variables (and their derivatives for secondorder ODEs) on cell boundaries obtain a set of discrete equations.

Apply the constraint conditions to get the full discretized equations.


Governing equations
Twodimensional, steadystate, incompressible, isothermal NavierStokes 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 cellspecific subscripts (i, j) on dependent variables are omitted. Terms not explicit in Eqs. (5)(7) are lumped into the right hand side as pseudosource terms, such as
2
2
= 1 ai [ (u + v ) p]dx (8)
1 2ai
ai
x y
x2
The inhomogeneous pseudosource terms in Eqs. (5)(7) are then expanded in Legendre polynomials and truncated at the zeroth order, which leads to consistent secondorder 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 pseudosource 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 xaveraged 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 transverseaveraged 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 yaveraged 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 psudosource 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 cellaveraged operator
1 ai,j bi,j
4ai,jbi,j
ai,j
bi,j
on Eqs. (2), (3), and (4) yields the relations between pseudosource 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 xmomentum and ymomentum equations, respectively. The last three constraint equations are obtained by imposing the condition that the cellaveraged 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 pseudosource 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 xmomentum 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 ymomentum equations can be obtained by replacing variable u by v and source term f2 by f3 from the discrete xmomentum 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).

Solver Details

Algorithm for iterative procedure

The steps involved in the solution procedure are as follows.

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

For given nthlevel pressure (), () the discretized momentum equations are solved for transverseaverage 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 GaussSeidel method. Thus, both velocity and pressure are corrected one after the other in the internal loop.

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

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.

Repeat steps 24 until convergence is achieved. The error tolerance limit () used for the convergence criteria is based on the rootmeansquared 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 crossstream 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 Uvelocity, 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.

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 noslip surfaces increases. Hence, a clear flow separation does not occur till higher values of Reynolds number are attained.

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 transverseaveraged 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 = Pseudosource terms
u = Velocity component in x direction
u , v = Transverseaveraged u and v velocity components, respectively uc, vc = Cellaveraged 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 = Cellaveraged quantity
i, j = Index in x and y directions, respectively
Superscripts
x = xaveraged quantity y = yaveraged quantity
xy = cellaveraged 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 viceverse.
REFERENCES

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

Toreja A. J. and Rizwanuddin, Hybrid numerical methods for convectiondiffusion problems in arbitrary geometries, Comput. Fluids, vol. 32, no. 6, pp. 835872, 2003.

Patankar S. V. Numerical Heat Transfer and Fluid Flow, Taylor and Francis Publishers 2005, p. 126130.

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. 189201, 1997.

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