Numerical solution of the incompressible Navier–Stokes equations in primitive variables and velocity–vorticity formulation

Applied Mathematics and Computation 222 (2013) 575–588

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Numerical solution of the incompressible Navier–Stokes equations in primitive variables and velocity–vorticity formulation V.C. Loukopoulos a,⇑, G.T. Messaris a, G.C. Bourantas b a b

Department of Physics, University of Patras, Patras, 26500 Rion, Greece Department of Medical Physics, School of Medicine, University of Patras, Patras, GR 26500 Rion, Greece

a r t i c l e

i n f o

Keywords: Navier–Stokes equations Primitive variables Velocity–Vorticity Numerical solution LFDM Staggered grid

a b s t r a c t A finite-difference method is presented for the numerical solution of the Navier–Stokes equations of motion of a viscous incompressible fluid in two dimensions in primitive-variables and velocity–vorticity formulation. For the case of primitive-variables, using a staggered grid and introducing an auxiliary function of the coordinate system and considering the form of the initial equation on lines passing through the nodal point (x0, y0) and parallel to the coordinate axes, we can separate it into two parts that are finally reduced to ordinary linear differential equations, one for each dimension. Discretization of these equations leads to a system of linear equations in n-unknowns which is solved by an iterative technique and the method converges rapidly giving satisfactory results. For the pressure variable we consider a pressure Poisson equation with suitable Neumann type boundary conditions. In case of velocity–vorticity formulation similar procedure is used, for a collocated grid. Numerical results up to Reynolds number 5000, confirming the accuracy of the proposed method, are presented for configurations of interest, such as Poiseuille flow, liddriven cavity flow in square domain and flow in a backward-facing step in the presence of a pressure gradient. Ó 2013 Elsevier Inc. All rights reserved.

1. Introduction The Navier–Stokes equations for incompressible viscous flow are coupled partial differential equations of the velocity and pressure. Two are the most popular approaches for the solution of these equations: the velocity–vorticity approach [1–3] and the primitive-variable approach [4,5]. The velocity–vorticity approach requires the transformation of the Navier–Stokes equations into equations of derived variables (velocity and vorticity components). Pressure and velocity components can be easily determined once these derived variables are known. An advantage of solving the transformed equations is that there are two and three equations for two- and three-dimensional problems, respectively, which are fewer than the number of equations for the corresponding problems in the primitive-variables approach. However, imposition of boundary conditions in the velocity–vorticity approach may require substantial effort, especially for three-dimensional problems, which is a major disadvantage of this strategy [6,7]. For two-dimensional problems, there is only one vorticity component and the two velocity components can be substituted by the stream function. Thus, the velocity–vorticity approach is generally known as the stream function-vorticity approach [8,9]. The two equations of stream function and vorticity can be reduced to a fourth-order partial differential equation ⇑ Corresponding author. E-mail address: [email protected] (V.C. Loukopoulos). 0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.08.002

576

V.C. Loukopoulos et al. / Applied Mathematics and Computation 222 (2013) 575–588

of stream function. This implies that it should be possible to reduce the four Navier–Stokes equations for three dimensional problems to two equations in two variables. Usually, this approach is followed when the physical problem presents symmetry around an axis. Unfortunately, the velocity–vorticity approach does not allow a straightforward reduction of equations. In the primitive-variables approach, the Navier–Stokes equations are to be solved for primitive variables namely pressure and velocity components. The number of equations is three and four for two- and three-dimensional problems, respectively. The imposition of boundary conditions in this approach is more straightforward. Nevertheless, there are different kinds of limits when primitive variables formulation is used. One of these limits consists of the compatibility condition known as the inf– sup condition or the LBB condition. The consequence of not satisfying this condition is the appearance of severe spatial oscillations in the pressure field, usually called spurious pressure modes. Consequently, two sets of elements or grids are necessary for methods such as the Galerkin-based FEM methods or FDM in order to satisfy the LBB condition [9]. The solution of Navier–Stokes equations remains an attractive research theme for many researchers [1–5,7–16]. Chantasiriwan in [7] has proposed an alternative approach for solving the incompressible Navier–Stokes equations, where the governing equations are reduced to higher-order partial differential equations of the remaining velocity components. This approach reduces the number of the Navier–Stokes equations for two- and three-dimensional problems by two without requiring the transformation of the Navier–Stokes equations into equations of derived variables like the vorticity–velocity approach. The number of equations is reduced by eliminating pressure and one of the velocity components from the governing equations. As a result, there is only one equation for two-dimensional problems and two equations for three-dimensional problems. Imposition of boundary conditions in the proposed approach is simple because the variables to be solved for are still primitive variables. Johnston and Liu [5] have studied incompressible Navier–Stokes equations based on a primitive variable formulation. The pressure was treated explicitly in time, completely decoupling the computation of the momentum and kinematic equations. The authors propose that their method consists an efficient Navier–Stokes solver and the key to the schemes was a Neumann boundary condition for the pressure Poisson equation which enforces the incompressibility condition for the velocity field. In the present work we present a numerical scheme for the solution of the incompressible Navier–Stokes equations based on both formulations, that is, a primitive-variables formulation in which the incompressibility constraint has been replaced by a pressure Poisson equation and a velocity–vorticity formulation. In both cases, the discretization of the Navier–Stokes equations result in a system of linear algebraic equations, where the matrix of the coefficients of the unknown functions is diagonally dominant, ensuring thus the convergence of the numerical algorithm. The discretization of the equations of the steady flow is obtained using a finite difference scheme of second order accuracy and suitable mathematical transformations. This method is called LFDM (Linearized Finite Differences Method). The LFDM has been used for the solution of the unsteady Navier–Stokes equations in stream function-vorticity formulation for every system of coordinates [16], as well as for the solution of the steady Navier–Stokes equations, which can include the terms of external forces BFD flows (Biomagnetic Fluid Dynamics) [17] and natural convection problems considering the Boussinesq formulation [18]. In addition, LFDM method has been used to simulate turbulent flows [19] and the three dimensional Convection–Diffusion problem for a Cartesian system of coordinates [20]. Parallel to the construction and development of numerical methods such as finite element method (FEM), finite volume methods (FVM), boundary element method (BEM) and Meshless (or meshfree) methods, the finite differences method (FDM) remain a useful tool for the internal construction of some of them. For example, finite differences can be part of meshless methods [21,22]. In recent years meshless methods have emerged as a class of effective numerical methods which are capable of avoiding the difficulties encountered in conventional computational mesh based methods, such as, meshing complex geometries, mesh distortion due to large deformation and remeshing in moving boundary problems [23]. The paper is organized as follows. In Section 2, we formulate the Navier–Stokes equations, a Poisson equation for the pressure and set the Neumann type boundary conditions for the pressure. Moreover the velocity–vorticity formulation is presented. In Section 3, we reformulate the equations in the LFDM framework. The LFDM has been used successfully for the solution of the Navier–Stokes equations either in stream function-vorticity formulation [16,17] or in velocity–vorticity formulation for a cylindrical coordinate system [24]. In the present work, we extend and apply the above mentioned numerical technique for the solution of the Navier–Stokes equations in primitive-variables formulation for first time and for this reason the procedure’s steps are presented analytically. In Section 4, we study three benchmark problems, that is, the Poiseuille flow, the lid-driven cavity flow in square domain and the flow in a backward-facing step in the presence of a pressure gradient. In Section 5, we show results for the case of velocity–vorticity formulation for high values of Reynolds number. Thus, the accuracy of the proposed scheme is demonstrated yet at high values of the Reynolds number up to Re = 5000. Finally, in Section 6, we gather our conclusions. 2. Equations of the problem 2.1. Primitive variables formulation The two-dimensional incompressible Navier–Stokes equations (with the absence of body forces), in nondimensional form, are the following nonlinear partial differential equations

V.C. Loukopoulos et al. / Applied Mathematics and Computation 222 (2013) 575–588

u

! @u @u @p 1 @ 2 u @ 2 u ; þt ¼ þ þ @x @y @x Re @x2 @y2

! @t @t @p 1 @ 2 t @ 2 t ; u þt ¼ þ þ @y Re @x2 @y2 @x @y @u @ t þ ¼ 0; @x @y

577

ð1Þ

ð2Þ

ð3Þ

where Re is the Reynolds number. The equation of the pressure can be derived from (1) and (2) after differentiation with respect to x and y, respectively. Eqs. (1) and (2) can be written as

" !# @p @u @u 1 @ 2 u @ 2 u ¼ u þt  þ ¼ G1 ðx; yÞ; @x @x @y Re @x2 @y2 " !# @p @t @t 1 @2t @2t ¼ u þt  þ ¼ G2 ðx; yÞ: @y @x @y Re @x2 @y2

ð4Þ

ð5Þ

Finally, the equation of the pressure takes the form

@ 2 p @ 2 p @G1 @G2 þ ¼ þ ¼ Sðx; yÞ: @x2 @y2 @x @y

ð6Þ

Eq. (6) is a Poisson equation and the boundary conditions for the pressure are of the Neumann type and are defined from Eqs. (4) and (5). We notice that the numerical solution of Eq. (6) is acceptable if Green’s theorem is verified, that is

Z Z

r2 pdxdy ¼

Z Z

S dx dy ¼

Z

@p dC @n

ð7Þ

where C is the boundary of the domain X, op/on = rp  n, and n is the unit outward normal vector on the boundary C. 2.2. Velocity–vorticity formulation ^ the associated  ¼ xk  ) and x For two-dimensional problems, with (u, t) being the components of the velocity vector (u  ¼ru  , the vorticity transport equation, in non dimensional form, can be written as vorticity vector, given by x

  rx ¼ u

1 2 r x; Re

ð8Þ

while the velocity components satisfy the equations

r2 u ¼ 

@x @x and r2 t ¼ : @y @x

ð9Þ

Thus the solution of the vorticity transport equation, in combination with the velocity Poisson equations with reference to boundary conditions, provides the velocity and the vorticity distribution all over the domain. The boundary conditions are  ¼ ðr  u ¼u  C and x  ÞC , where C is the boundary of the spatial domain. Moreover, the boundary values of vorticity are calu culated by using x ¼ @u  @@xt along with suitable finite differences expressions at each side and corners of the spatial domain, @y such as forward, central and backward differences. Eq. (8) is solved using LFDM described in details in Section 3.1. Eqs. (9) are Poisson type equations and are solved in a fast and efficient manner by using a simple SOR (Successive Over Relaxation) method, after its discretization with central finite differences scheme. 3. Numerical solution 3.1. The LFDM method The LFDM (Linearized Finite Differences Method) is described in the present section. We consider a uniform grid of mesh points with size Dx and Dy in x and y direction, in which the grid lines are parallel to the coordinate axes at each point (x, y). Let (x0, y0) be a nodal point and (x0 + Dx, y0), (x0, y0 + Dy), (x0  Dx, y0), (x0, y0  Dy), be neighboring points on the grid lines passing through this nodal point. All quantities at this set of grid points are denoted by the subscripts 0, 1, 2, 3, and 4, respectively, Fig. 1. Eq. (6) is a Poisson equation and can be solved fast and efficiently with the use of a simple SOR (Successive Over Relaxation) method, after its discretization with central finite differences scheme. The main effort is focused on solving the other

578

V.C. Loukopoulos et al. / Applied Mathematics and Computation 222 (2013) 575–588

2

(i, j +1)

(i -1, j)

(i +1, j) 0

3

1 (i, j)

4

(i, j -1)

Fig. 1. Five points stencil.

equations of the system, namely (1) and (2). When the difference equations are associated with matrices which are diagonally dominant, then they are amenable to solution by iterative methods (such as SOR). Thus, our aim is to make the matrices of the coefficients of the unknowns diagonally dominant. The difficulty in using the standard central-difference numerical schemes is due to the presence of the non linear terms

u

@u @u þt @x @y

and u

@t @t þt @x @y

which cause diagonal dominance to be lost in the matrix associated with the approximating finite-difference equations. To avoid this difficulty Eq. (1) (and Eq. (2)), after nondimensionalization, can be split into two equations, namely,

@2u @u @p  u Re ¼ Re þ Qðx; yÞ @x2 @x @x

ð10Þ

@2u @u  t Re ¼ Qðx; yÞ; @y2 @y

ð11Þ

where Re ¼ uL m is the Reynolds number of the flow, u a characteristic velocity, L a characteristic length, m the kinematic viscosity of the flow, Q is an arbitrary unknown function of x, y and we postulate that each equation is valid in the neighbor of (x0, y0). In particular (10) is valid along the line 103 and (11) is valid along the line 204, Fig. 1. Eq. (11) is transformed locally for x0  Dx 6 x 6 x0 þ Dx, keeping y equal to y0 and setting

uðx; y0 Þ ¼ Cðx; y0 Þe½cðx;y0 Þ ;

where cðx; y0 Þ ¼ 

1 2

Z

Re uðx; y0 Þdx

ð12Þ

x0 x

and C(x, y0) is an unknown function of x and y0. In the same way (11) is transformed locally for y0  Dy 6 y 6 y0 þ Dy, keeping x equal to x0 and setting

tðx0 ; yÞ ¼ Hðx0 ; yÞe½qðx0 ;yÞ ; where qðx0 ; yÞ ¼ 

1 2

Z

y

Re tðx0 ; yÞdy

ð13Þ

y0

and H(x0, y) is an unknown function of x0 and y. Thus, substitution of (12) to (10) and (13) into (11) gives the following system of differential equations

!   @2C 1 @u Re2 u2 @p ecðx;y0 Þ ; þ Re  þ Q ðx; y Þ C ¼ Re 0 @x2 2 @x @x 4 ! @2H 1 @ t Re2 t2 H ¼ Q ðx0 ; yÞeqðx0 ;yÞ ; þ Re  @y2 2 @y 4

ð14Þ

ð15Þ

in which the non linear terms have been eliminated and partial derivatives are actually ordinary derivatives. Next we replace in Eqs. (14) and (15) all partial derivatives by central–difference approximations at the point (x0, y0), eliminate Q(x0, y0) between the two new equations and substitute C and H by their appropriate expressions from (12) and (13) at the corresponding points around (x0, y0). In this way we obtain from Eqs. (14) and (15) the following equation:

!   Dx2 @u Re2 u20 Dx2 Dy2 @ t Re2 t20 Dx2 @p 2 u0 ¼ Dx2 Re u1 e þ u3 e þ k u2 e þ k u4 e þ 2 þ Re    2k þ Re   @x 0 @x 2 4 2 @y 0 4 c1

c3

2

q2

2

q4

ð16Þ

V.C. Loukopoulos et al. / Applied Mathematics and Computation 222 (2013) 575–588

579

where k ¼ DDyx. However, the matrix of the coefficients of the unknowns associated with (16) is not necessarily diagonally dominant, which is a prerequisite for the convergence of the iterative procedure. Diagonal dominance is obtained by expanding the exponential terms in Taylor series at the point (x0, y0) and keeping the sufficient number of terms so that the order of the truncation error is conserved. In particular we expand us(x, y) = Reu(x, y) in a Taylor series at the point (x0, y0) in the direction of x increasing, so that Eq. (12) can be integrated to give c(x, y0) in powers of (x  x0). If we set in this latter equation successively x = x0 + Dx and x ¼ x0  Dx we obtain c1 and c3, respectively, in powers of Dx. The values of c1 and c3 are used in order to take expressions of ec1 and ec3 in the form of Taylor series, which will be substituted in the first two terms of the left– hand side of Eq. (16). Thus, for the quantities us ðx; yÞ, c1, c3, ec1 and ec3 we take the following expressions u00 ðx ;y Þ

us ðx; y0 Þ ¼ us ðx0 ; y0 Þ þ u0s ðx0 ; y0 Þðx  x0 Þ þ s 2!0 0 ðx  x0 Þ2 þ    ;    s 1 @ 2 us Dx2  12 Dx3 þ    ; c1 ¼  12 us0 Dx  14 @u @x 0 @x2 0   1 @ 2 us s Dx2 þ 12 ð @x2 Þ0 Dx3 þ    ; c3 ¼ 12 us0 Dx  14 @u @x 0     @u ec1 ¼ 1  12 Re u0 Dx  Re Dx2 þ 18 Re2 u20 Dx2 þ O Dx3 ; 4 @x 0 @u   Dx2 þ 18 Re2 u20 Dx2 þ O Dx3 : ec3 ¼ 1 þ 12 Re u0 Dx  Re 4 @x 0 We follow the same procedure to deduce likewise expressions for the other two terms of the left–hand side member of Eq. (16). In this way and using the equation of continuity, Eq. (16) or (1) takes the form

a1 u1 þ a2 u2 þ a3 u3 þ a4 u4  a0 u0 ¼ Re

@p 2 Dx ; @x

ð17Þ

where the coefficients are:

h    i a1 ¼ 1  12 Reu0 Dx þ 18 Re2 u20 Dx2 þ O Dx3 ; a2 ¼ 1  12 Ret0 Dy þ 18 Re2 t20 Dy2 þ O Dy3 k2 ; h    i a3 ¼ 1 þ 12 Reu0 Dx þ 18 Re2 u20 Dx2 þ O Dx3 ; a4 ¼ 1 þ 12 Ret0 Dy þ 18 Re2 t20 Dy2 þ O Dy3 k2 ;     a0 ¼ 2 þ 14 Re2 u20 Dx2 þ 14 Re2 t20 Dx2 þ 2k2 þ O Dx3 þ O Dy3 and the suffix 0 denotes that the corresponding factors have been evaluated at (x0, y0). Following the same procedure for the Eq. (2) we obtain

a1 t1 þ a2 t2 þ a3 t3 þ a4 t4  a0 t0 ¼ Re

@p 2 Dx ; @y

ð18Þ

Eqs. (17) and (18) are consistent to the initial Eqs. (1) and (2) and are valid at every nodal point of the grid. The matrix of the coefficients of the unknowns in (17), (18) is diagonally dominant, such that it makes the algorithm to be convergent and stable, apart from being accurate and efficient. We notice that the same procedure is followed for the case of vorticity equation, Eq. (8). 3.2. Staggered grid In regard to the grid, when central differences are used for the discretization of the incompressible Navier–Stokes equations, with the use of a collocated grid, the resulting difference equations are of the form that, when presented with the nonsensical velocity and pressure distributions, will tend to perpetuate these distributions. A fix of previous problem is attained if we maintain the central differencing but stagger the grid [25]. A staggered grid is illustrated in Fig. 2. The pressure is calculated at the solid grid points and the velocities and stream function are calculated at the open grid points. As far as, the velocities, the stream function and pressure are concerned, we calculate their values at the rest points of the grid as mean values of successively grid points. Also, using fictitious points for the velocities and pressure, the use and calculation of central differences is possible at every point of the grid. Actually, we use the nonconservation form of the Navier–Stokes equations and not the conservation form. This explains why it is necessary to estimate the values of the unknowns’ functions at any open or solid grid point. When conservation form of the Navier–Stokes is used, first order partial derivatives are appeared in equations and using central differences approximation the value of the derivative at point (i, j) is estimated from the values at points (i + 1, j) and (i  1, j) or the values at points (i, j + 1) and (i, j  1), Fig. 1. For the case of nonconservation form and specifically for the present method it is required to know the values of the unknown functions at the five points of the stencil, during the iterative procedure. For our computations different types of staggered grid were used, as they are presented in Fig. 2. Also, the components of the velocity were calculated either at the same grid points or not. Among the various types I to IV, Type II was the most effective in our computations. In this type, the components of the velocity were calculated at the corners of the corresponding square cell and the pressure at the center of the cell.

580

V.C. Loukopoulos et al. / Applied Mathematics and Computation 222 (2013) 575–588

Type I

Typ e II

Type III

Typ e IV

Fig. 2. Staggered grid. Solid grid points for the pressure and open grid points for the velocities (and the stream function).

Adopting the staggered grid of Type II, it is clear that, for the pressure, points of the stencils that correspond to different internal points of the grid, next to the boundaries, are located out of the computational domain. Because of this we use fictitious points outside of the computational domain for the pressure and the velocity components. These values of the velocity are necessary for the estimation of G1(x, y) and G2(x, y), the functions on the right side of Eqs. (4) and (5), which are used for the solution of (6) and the calculation of the Neumann boundary conditions. On the lower (and upper) boundary it is u(i, j) = 0, that is @uði;jÞ ¼ 0 and from the continuity equation it is deduced that t(i, j  1) = t(i, j + 1). Similarly, for the left @x ði;jÞ (and right) boundary is t(i, j) = 0, that is @ t@y ¼ 0 and from the continuity equation it result that is u(i  1, j) = u(i + 1, j). 3.3. The LFDMPLE algorithm The procedure that we develop for the calculation of the flow field has been given the name LFDMPLE which stands for Linearized Finite Differences Method for Pressure-Linked Equations. For the case of primitive variables formulation, the sequences of the important operations in the order of their execution are: 1. Make initial guesses at the interior points of the computational domain for the velocity (u, t) and the pressure field p, as well as for the boundary conditions for the pressure. 2. Solve the momentum equations, such as Eqs. (1) and (2) and calculate a new estimation, considering p as a known quantity. 3. Calculate a new estimation for p by solving Eq. (6) once, considering u and t known. 4. Return to step  2, and repeat the whole procedure until the criterion convergence is satisfied. The criterion convergence s   5 used is max1  FFsþ1ðx;yÞ  6 10 , where s is the number of the iterations performed and F stands for the unknown functions ðx;yÞ u, t, and p. Additionally, for the velocity–vorticity formulation similar steps are followed.

V.C. Loukopoulos et al. / Applied Mathematics and Computation 222 (2013) 575–588

581

4. Numerical examples for primitive variables formulation 4.1. Model problem 1: Poiseuille flow between two parallel plates We consider the flow of an incompressible and viscous fluid between two parallel plates in the presence of a pressure gradient (dp/dx – 0). Flow is assumed to be fully developed as it passes the inlet at x = 0. The boundary conditions of the problem are

Inflow ðx ¼ 0; 0 6 y 6 hÞ : Outflow ðx ¼ L; 0 6 y 6 hÞ :

t ¼ 0:   t ¼ 0Þor @x@ ¼ 0 : u ¼ 0; t ¼ 0: u ¼ 0; t ¼ 0:

u ¼ uðyÞ;

ðu ¼ uðyÞ;

Upper plate ðy ¼ h; 0 6 x 6 LÞ : Lower plate ðy ¼ 0; 0 6 x 6 LÞ :

The pressure at the outlet is set to zero. The boundary conditions for the pressure equation on the other boundaries are of the Neumann type:



  @p 1 ¼  AðVÞ  r2 V  n @n Re C

Also, the unknown quantities (u, t) at the exit of the channel can be assumed independent of x, (o/ox = 0). Thus, the outflow conditions for all quantities are determined from the interior grid points of the computational domain see Fig. 3. As far as the stream function W, the boundary conditions are easily implemented using u = oW/oy, t = oW/ox, since the velocity components are known (fully developed flow at the entrance and non-slip conditions on the plates). Thus, the value at the entrance is W(0, y) = 2y2  (4/3)y3. In Fig. 4 results are shown for Reynolds number 1000, the number of nodes being 140 in x-direction and 20 in y-direction. Additional testing of the code was conducted for denser grids such as 150  25, 170  30 etc., ensuring grid independency. The contours of stream function, u-component of velocity and pressure verify the theory. This means that, u-component of velocity retain the parabolic profile along the entire flow field and the stream functions lines are skew lines parallel to x-axis. Our results were compared qualitatively with corresponding results, obtained using a mixed finite element solver [26]. 4.2. Model problem 2: Lid-driven cavity flow in square domain The second model problem is the lid-driven cavity flow in a square domain. It is one of the most widely used benchmark problems to test steady state incompressible fluid dynamics codes. The boundary conditions of the problem on a [0,1]  [0,1] are given by

Lower boundary u ¼ 0; Upper boundary u ¼ 1;

t ¼ 0 on y ¼ 0; 0 6 x 6 1 t ¼ 0 on y ¼ 1; 0 6 x 6 1

t ¼ 0 on x ¼ 0; 0 6 y < 1 Right boundary u ¼ 0; t ¼ 0 on x ¼ 1; 0 6 y < 1 Left boundary u ¼ 0;

Left, right and down walls are fixed. The top wall is assumed to be moving with a given velocity (u = 1, for the present study). Results are compared with analogous results obtained by Lewis et al. [27], Chinchapatman et al. [28], Ghia et al. [29] and the agreement is obvious. The former developed a time marching algorithm and solutions were produced using the characteristic-based-split (CBS) scheme and the characteristic Galerkin procedure, the second used an RBF (Radial Basic Function)

Fig. 3. Fictitious points.

582

V.C. Loukopoulos et al. / Applied Mathematics and Computation 222 (2013) 575–588

1 0.5 0

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

16

18

20

0

2

4

6

8

10

12

14

16

18

20

1 0.5 0

1 0.5 0

Fig. 4. Contours of stream function, u-component of velocity and pressure for Re = 1000, L = 20, h = 1.

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fig. 5. Contours of the stream function and pressure for Stokes flow, Re  0.

meshless method, whereas the latter applied a coupled strongly implicit and multigrid method.In Figs. 5–8 results are shown for Reynolds number Re  0 (Stokes flow), 100, 400 and 1000, the number of nodes being 100 in x-direction and 100 in ydirection for Re  0 and 100. The number of grid’s points was increased to 125 in x-direction and 100 in y-direction for Re = 400 and 1000. Additional testing of the code has been conducted for denser grids, ensuring the convergence of the method. In Fig. 9 the current numerical solutions are compared with results of previous work [29] and in this way verify the accuracy of the numerical scheme. 4.3. Model problem 3: Flow in a backward-facing step The third problem presented is the flow over a backward-facing step. Fluid is flowing in a channel of width h downstream of the origin and h/2 upstream of the origin, separated by a backward-facing step. The flow is assumed to be fully developed as it passes the inlet at x = 0. The problem domain is the channel starting at the inlet and extending downstream a distance L long enough for the flow to again become fully developed. The problem has been studied extensively in [30,31]. The selection of an outflow boundary condition that allows an accurate numerical simulation of a flow problem to be achieved was the main question in [30]. The Navier–Stokes equations

583

V.C. Loukopoulos et al. / Applied Mathematics and Computation 222 (2013) 575–588 1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.5

0.6

0.7

0.8

0.9

1

Fig. 6. Contours of the stream function and pressure for Re = 100.

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0

0.1

0.2

0.3

0.4

Fig. 7. Contours of the stream function and pressure for Re = 400.

provide a widely applicable mathematical description for the analysis of fluid motion. Individual flow problems are differentiated by specific choices of fluid properties, the geometry of interest and a suitable set of boundary and initial conditions for the flow variables. Generally, the boundary conditions are the most difficult to define (precisely) since it is difficult to isolate a fluid system from the effects of its environment. This is particularly true for uncontained flows, which require some rational specification of how the fluid enters and leaves the domain of interest. In the mentioned research paper the downstream distance L is taken to be 30h in order for the flow to be fully developed. This completes the specification of the problem. Nevertheless, the previous set of questions is beyond the scope of the present study.

584

V.C. Loukopoulos et al. / Applied Mathematics and Computation 222 (2013) 575–588

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fig. 8. Contours of the stream function and pressure for Re = 1000.

0.6

1.0

LFDM Ghia et al.

0.4

v-velocity (y=0.5)

0.8

Y

0.6

0.4

0.2

0.2 0.0 -0.2 -0.4

LFDM Ghia et al.

0.0 -0.4

-0.2

0.0

0.2

0.4

0.6

0.8

1.0

-0.6

1.2

0.0

u-velocity (x=0.5)

(a)

0.2

0.4

0.6

0.8

1.0

X

(b)

Fig. 9. Square lid-driven cavity problem. (a) The u – velocity on the vertical section, for Re = 1000. (b) The t – velocity on the horizontal section y = 0.5, for Re = 1000. Results of Ghia et al. [29] are compared with the current numerical solutions.

For our computations the boundary conditions of the problem are

t ¼ 0 on y ¼ 0:5; 0 6 x 6 L: t ¼ 0 on y ¼ 0:5; 0 6 x 6 L: u ¼ 12yð1  2yÞt ¼ 0 on x ¼ 0; 0 6 y 6 h=2: u ¼ 0; t ¼ 0 on x ¼ 0; h=2 6 y 6 0: @u   ¼ 0 or u ¼ 3 14  y2 ; t ¼ 0 on x ¼ L; h=2 6 y 6 h=2: @x

Lower plate u ¼ 0; Upper plate Inflow Step Outflow

u ¼ 0;

The results were obtained either with parabolic exit profile of the u-component of the velocity or using the boundary condition @u=@x ¼ 0. We compare our results with the results of the numerical experiments [32–34]. Figs. 10–12 present contours of the stream function, u component of the velocity and pressure for Re equal to 50, 100 and 200. For Reynolds number 50, the number of nodes being 150 in x-direction and 20 in y-direction. Additionally running of the code was executed for more dense grids as 200  25, 200  30 etc., ensuring in this way grid independency. Moreover, for Reynolds number 100 and 200, grids as 220  30, 250  30 was required to ensure grid independency. It can be seen from the figures that a

585

V.C. Loukopoulos et al. / Applied Mathematics and Computation 222 (2013) 575–588 0.5 0 -0.5 0

5

10

15

20

25

30

5

10

15

20

25

30

5

10

15

20

25

30

0.5 0 -0.5 0 0.5 0 -0.5 0

Fig. 10. Contours of stream function, u-component of the velocity and pressure for Re = 50, L = 30, h = 1 and step = 0.5.

0.5 0 -0.5 0

5

10

15

20

25

30

5

10

15

20

25

30

5

10

15

20

25

30

0.5 0 -0.5 0 0.5 0 -0.5 0

Fig. 11. Contours of stream function, u-component of the velocity and pressure for Re = 100, L = 30, h = 1 and step = 0.5.

0.5 0 -0.5 0

5

10

15

20

25

30

5

10

15

20

25

30

5

10

15

20

25

30

0.5 0 -0.5

0

0.5 0 -0.5 0

Fig. 12. Contours of stream function, u-component of the velocity and pressure for Re = 200, L = 30, h = 1 and step = 0.5.

recirculation zone is formed downstream of the step face. For all values of Reynolds number our solutions agree very well with respect to the vortex shape and the length of the circulation region for stream function contours, as well as to the range of values of the stream function and the pressure. In particular, the recirculation zone details obtained by the present method are compared with those obtained by using global Radial Basis Functions (RBF) [32], the h-p finite elements [33] and the collocation meshless method [34]. The numerical results obtained with LFDM, regarding the length of the recirculation and the value of wmin are in agreement with previous obtained numerical resulats, as shown in Table 1.

586

V.C. Loukopoulos et al. / Applied Mathematics and Computation 222 (2013) 575–588

Table 1 Backward-facing step for Re = 200, length of recirculation region, primary vortex strength and location.

Length of recirculation

wmin wmin location

h-p finite elements [32]

RBF [33]

MPC [34]

LFDM

2.67 0.0331 (1.0021,0.2030)

2.72 0.0315 (1.333,0.2167)

2.64 0.0331 (1.0001,0.2000)

2.66 0.0330 (1.0011,0.2070)

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.5

0.6

0.7

0.8

0.9

1

Fig. 13. Contours of the stream function and the vorticity for Re = 3200.

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0

0.1

0.2

0.3

0.4

Fig. 14. Contours of the stream function and the vorticity for Re = 5000.

587

V.C. Loukopoulos et al. / Applied Mathematics and Computation 222 (2013) 575–588 0,6

1,0

LFDM Ghia et al.

0,4

0,8

v-velocity (y=0.5)

0,2

Y

0,6

0,4

0,2

0,0

-0,2

-0,4

LFDM Ghia et al.

0,0 -0,4

-0,2

0,0

0,2

0,4

0,6

0,8

-0,6 0,0

1,0

0,2

0,4

0,6

0,8

1,0

X

u-velocity (x=0.5)

(b)

(a) 0,6

1,0

LFDM Ghia et al.

0,4

0,8

v-velocity (y=0.5)

0,2

Y

0,6

0,4

0,2

0,0

-0,2

-0,4

LFDM Ghia et al.

0,0 -0,4

-0,2

0,0

0,2

0,4

0,6

0,8

1,0

-0,6 0,0

0,2

0,4

0,6

u-velocity (x=0.5)

X

(c)

(d)

0,8

1,0

Fig. 15. Square lid-driven cavity problem. (a) The u – velocity on the vertical section x = 0.5, for Re = 3200. (b) The t – velocity on the horizontal section y = 0.5, for Re = 3200. (c) The u – velocity on the vertical section x = 0.5, for Re = 5000. (d) The t – velocity on the horizontal section y = 0.5, for Re = 5000. Results of Ghia et al. [29] are compared with the current numerical solutions.

5. Numerical examples for velocity–vorticity formulation and high values of Reynolds number The lid-driven cavity flow problem in a square domain is solved in velocity–vorticity formulation, where, we are able to derive results for high values of the Reynolds number. Thus, in Figs. 13–15, results are presented for Re = 3200 and 5000. The current numerical solutions are compared with results of previous work [29] and in this way verify the accuracy of the numerical scheme. The number of nodes being 170 in x-direction and 150 in y-direction for Re = 3200. The number of grid’s points was increased to 225 in x-direction and 200 in y-direction for Re = 5000. 6. Conclusions In the present article we extend and apply the LFDM numerical method for the approximate solution of the Navier–Stokes equations of motion of a viscous incompressible fluid in two dimensions in the primitive-variables and velocity–vorticity formulation. Three representative examples for configurations of interest, demonstrate the accuracy of the proposed method. Numerical results up to Reynolds number 5000, were achieved. Comparison of primitive variables and velocity–vorticity formulation reveal that results for high values of Reynolds number (Re > 1000) can be derived for the second. The numerical results presented herein confirm the important characteristics of the method: the numerical results of the three model problems show that the method can simulate fluid flow problems accurately and efficiently in any formulation of the Navier–Stokes equations. The method is sufficiently flexible to handle a variety of boundary conditions, including Dirichlet, Neumann etc. The application of the method is simple.

588

V.C. Loukopoulos et al. / Applied Mathematics and Computation 222 (2013) 575–588

Although results are presented for two dimensional flows, we estimate that the previous mentioned technique can also be used for the solution of the three dimensional flows, as well as the solution of unsteady flow, in the presence of external forces, in any system of coordinates. LFDM method could possibly be improved, combined with other techniques. This could be done for example, if the equations of unsteady flow were solved, in which the variable of time can operate as a relaxation factor, or, if the LFDM method were associated with velocity and pressure correction methods, such as SIMPLE or SIMPLER [6] and on grid generation techniques for the mapping of the physical problem and the solution of more complex flows. Finally, another factor, which should be investigated, on the accuracy of a numerical solution, is the type of the order of the numerical boundary conditions used in computation. References [1] D.C. Lo, K. Murugesan, D.L. Young, Numerical solution of three-dimensional velocity–vorticity Navier–Stokes equations by finite difference method, Int. J. Numer. Meth. Fluids 47 (2005) 1469–1487. [2] D.L. Young, S.K. Yang, T.I. Eldho, Solution of the Navier–Stokes equations in velocity–vorticity form using a Eulerian–Lagrangian boundary element method, Int. J. Numer. Meth. Fluids 34 (2000) 627–650. [3] M.O.L. Hansen, J.N. Sørensen, W.Z. Shen, Vorticity–velocity formulation of the 3D Navier–Stokes equations in cylindrical co-ordinates, Int. J. Numer. Meth. Fluids 41 (2003) 29–45. [4] H. Johnston, J.-G. Liu, Finite difference schemes for incompressible flow based on local pressure boundary conditions, J. Comput. Phys. 180 (2002) 120– 154. [5] H. Johnston, J.-G. Liu, Accurate, stable and efficient Navier–Stokes solvers based on explicit treatment of the pressure term, J. Comput. Phys. 199 (2004) 221–259. [6] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, 1980. [7] S. Chantasiriwan, An alternative approach for numerical solutions of the Navier–Stokes equations, Int. J. Numer. Meth. Eng. 69 (2007) 1331–1344. [8] S.J. Liao, Higher-order stream function-vorticity formulation of 2D steady-state Navier–Stokes equations, Int. J. Numer. Meth. Fluids 15 (1992) 595– 612. [9] Y. Kim, D.W. Kim, S. Jun, J.H. Lee, Meshfree point collocation method for the stream-vorticity formulation of 2D incompressible Navier–Stokes equations, Comput. Meth. Appl. Mech. Eng. 196 (2007) 3095–3109. [10] D.C. Lo, D.L. Young, K. Murugesan, An accurate numerical solution algorithm for 3D velocity–vorticity Navier–Stokes equations by the DQ method, Commun. Numer. Meth. Eng. 22 (2006) 235–250. [11] C. de Falco, R. Sacco, G. Scrofani, Stabilized 3D finite elements for the numerical solution of the Navier–Stokes equations in semiconductors, Comput. Meth. Appl. Mech. Eng. 196 (2007) 1729–1744. [12] E.H. van Brummelen, H.C. Raven, B. Koren, Efficient numerical solution of steady free-surface Navier–Stokes flow, J. Comput. Phys. 174 (2001) 120–137. [13] Qi S. Zhang, Global solutions of Navier–Stokes equations with large L2 norms in a new function space, Adv. Differ. Equ. 9 (2004) 587–624. [14] S.A. Salem, On the numerical solution of the incompressible Navier–Stokes equations in primitive variables using grid generation techniques, Math. Comput. Appl. 11 (2006) 127–136. [15] J.S.B. Gajjar, N.A. Azzam, Numerical solution of the Navier–Stokes equations for the flow in a cylinder cascade, J. Fluid Mech. 520 (2004) 51–82. [16] V.C. Loukopoulos, A numerical technique for the solution of the Navier–Stokes equations of unsteady flow, Comput. Meth. Appl. Mech. Eng. 195 (2006) 534–550. [17] V.C. Loukopoulos, E.E. Tzirtzilakis, Biomagnetic channel flow in spatially varying magnetic field, Int. J. Eng. Sci. 42 (2004) 571–590. [18] V.C. Loukopoulos, Taylor vortices in annular heated spherical flow at medium and large aspect ratios, Phys. Fluids 17 (2005) 0181081–0181084. [19] E.E. Tzirtzilakis, M. Xenos, V.C. Loukopoulos, N.G. Kafoussias, Turbulent biomagnetic fluid flow in a rectangular channel under the action of a localized magnetic field, Int. J. Eng. Sci. 44 (2006) 1205–1224. [20] N. Niakas, V.C. Loukopoulos, C. Douskos, An accurate numerical solver on second order PDEs with variable coefficients in three dimensions, Appl. Math. Comput. 204 (2008) 50–57. [21] G.R. Liu, J. Zhang, H. Li, K.Y. Lam, B.B.T. Kee, Radial point interpolation based finite difference method for mechanics problems, Int. J. Numer. Meth. Eng. 68 (2006) 728–754. [22] M. Cheng, G.R. Liu, A novel finite method for flow simulation, Int. J. Numer. Meth. Eng. 39 (2002) 1161–1178. [23] X. Jin, G. Li, N.R. Aluru, Positivity conditions in meshless collocation methods, Comput. Meth. Appl. Mech. Eng. 193 (2004) 1171–1202. [24] N. Niakas, V.C. Loukopoulos, C. Douskos, LFDM method on the Navier–Stokes equations in three-dimensional flow in cylindrical coordinates on the von Kármán problem, in: AIP Conf. Proc., March 19, 1108, 2009, pp. 242–251. [25] J.D. Anderson, Computational Fluid Dynamics, The Basics with Applications, McGraw-Hill, Inc., 1995. [26] H.C. Elman, D.J. Silvester, A.J. Wathen, Finite Elements and Fast Iterative Solvers with Applications in Incompressible Fluid Dynamics, Oxford university press, 2008. [27] R.W. Lewis, P. Nithiarasu, K.N. Seetharamu, Fundamentals of the finite elements method for heat and fluid flow, John Wiley & Sons, Ltd., 2004. [28] P.P. Chinchapatnam, K. Djidjeli, P.B. Nair, Radial basis function meshless method for the steady incompressible Navier–Stokes equations, Int. J. Comput. Math. 84 (2007) 1509–1521. [29] U. Ghia, K.N. Ghia, C.T. Shin, High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method, J. Comput. Phys. 48 (1982) 387–411. [30] D.K. Gartling, A test problem for outflow boundary step conditions-flow over a backward-facing, Int. J. Numer. Meth. Fluids 11 (1990) 953–967. [31] P.M. Gresho, D.K. Gartling, J.R. Torczynski, K.A. Cliffe, K.H. Winters, T.J. Garratt, A. Spence, J.W. Goodrich, Is the steady viscous incompressible twodimensional flow over a backward-facing step ar Re = 800 stable?, Int J. Numer. Meth. Fluids 17 (1993) 501–541. [32] E. Barragy, Parallel finite element methods and iterative solution techniques for viscous incompressible flows. PhD thesis, University of Texas at Austin, 1993. [33] P.P. Chinchapatnam, K. Djidjeli, K. Nair, Radial basis function meshless method for the steady incompressible Navier–Stokes equations, Int. J. Comput. Math. 84 (2007) 1509–1526. [34] G.C. Bourantas, E.D. Skouras, V.C. Loukopoulos, G.C. Nikiforidis, Meshfree point collocation schemes for 2D steady state incompressible Navier–Stokes equations in velocity–vorticity formulation for high values of Reynolds numbers, CMES: Comput. Model. Eng. Sci. 1451 (2010) 1–33.