Is this finite difference approach correct?

Question

I am solving incompressible 2D Navier-Stokes equations with zero y-component velocity. The geometry is a simple 2D pipe of a length $L$ and diameter $W$ and there is only two boundary conditions:

1. Non-slip on the pipe walls: $u = 0$
2. Flow velocity at inlet $u = U_{max}$

and so, the two equations I have are:

$$ \frac{\partial u}{\partial x} = 0$$ $$ \frac{\partial P}{\partial x} = \mu \frac{\partial^2 u}{\partial y^2}$$

using finite difference discretization I came up with the following:

1.continuity equation: $$ u_{i,j} = u_{i+1,j} $$ 2. Pressure using forward difference: $$P_{i,j} = P_{i+1,j} - \frac{\mu \triangle x}{(\triangle y) ^2}(u_{i,j+1} + u_{i,j-1} - 2u_{i,j})$$ 3. Pressure using backward difference: $$P_{i,j} = P_{i-1,j} + \frac{\mu \triangle x}{(\triangle y) ^2}(u_{i,j+1} + u_{i,j-1} - 2u_{i,j})$$ 4. Manipulating equation (2) to get an equation for $u$: $$u_{i,j} = \frac{u_{i,j+1}}{2} + \frac{u_{i,j-1}}{2} - (\frac{P_{i+1,j} - P_{i,j}}{\frac{2\mu \triangle x}{(\triangle y) ^2}}) $$

I shall then use Gauss–Seidel method to solve the previous equations, but since I made two finite difference equation for u and P using the same partial differential equation I can't tell whether this is correct or not as I am getting very weird results.

So, is my approach correct?

Answer 1

There are two aspects to your question I think.

1) Do your equations match the physical problem you're trying to model?

2) Do your finite difference equations converge to the continuous ones as dx and dt approach zero?

First, I'd like to address question 1, since this affects question 2.

Question 1: Assuming you're looking for a general (developing) flow solution, this image is a helpful illustration:

As you can see, the entrance region (ER) extends to a critical point where the flow reaches a fully developed (FD) state. In the FD region, there's only 1 non-zero component of velocity, however, in ER there are multiple components.

The governing equation (for the entire domain) are momentum and continuity:

$\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \bullet \nabla \mathbf{u} = - \nabla p + \frac{1}{Re} \nabla^2 \mathbf{u}$

$\nabla \bullet \mathbf{u} = 0$

According to here, this simplifies to (assuming $u_{\theta} = 0, \frac{\partial}{\partial \theta} = 0$, which also implies that we're assuming a laminar solution)

r: $\frac{\partial u_r}{\partial t} + u_r \frac{\partial u_r}{\partial r} + u_z \frac{\partial u_r}{\partial z} = -\frac{\partial p}{\partial r} + \frac{1}{Re} \left[ \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_r}{\partial r} \right) - \frac{u_r}{r^2} + \frac{\partial^2 u_r}{\partial z^2} \right]$

z: $ \frac{\partial u_z}{\partial t} + u_r \frac{\partial u_z}{\partial r} + u_z \frac{\partial u_z}{\partial z} = -\frac{\partial p}{\partial z} + \frac{1}{Re} \left[ \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_z}{\partial r} \right) + \frac{\partial^2 u_z}{\partial z^2} \right] $

Continuity: $\frac{1}{r} \frac{\partial (r u_r)}{\partial r} +\frac{\partial u_z}{\partial z} = 0$

With this, we may apply general inlet BCs as you suggested:

$u_r = u_{r,inlet}(r,0)$

$u_z = u_{z,inlet}(r,0)$

Likely outlet BCs might be

$u_r = 0 $ (more fewer unknowns than $\frac{\partial u_r}{\partial z} = 0$)

$\frac{\partial u_z}{\partial z} = 0$

Question 2: Since there are many different approaches to solve equations in question 1, I will just briefly comment on some approaches rather than writing finite difference equations.

You may discretize in space and time and then apply a solution method (you suggested Gauss-Seidel). You may use staggered variables, where velocity and pressure are located on the cell face and center respectively, or you may use a collocated scheme, but then you will need to compute fluxes for the advection term in order to avoid pressure checkerboard phenomena. Typically, 2nd-order central difference schemes are applied to these equations for spatial derivatives, unless higher order discretization is needed. Since it seems you've written your equations for steady state, I'll assume that you're not interested in transient behavior and suggest using 1st-order (explicit Euler) time marching.