# Is this finite difference approach correct?

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?

• No, that doesn't make sense. You seem to be missing an equation. – David Ketcheson Apr 14 '16 at 8:52
• @DavidKetcheson I have another equation: $$\frac{\partial P}{\partial y} = 0$$ But I don't think adding it to the system of equation is going to make any difference, or is it? – Algo Apr 14 '16 at 9:13
• @Algo, You are right, you do not need the extra pressure equation since you are resolving the continuity equation which implicitly is an equation for pressure. It looks like you are using a collocated grid which is prone to an oscillatory pressure distribution and an unrealistic velocity field. This can be remedied using a staggered grid for the pressure and velocity variables. – namu Apr 15 '16 at 0:14
• As @namu said, your considering the momentum equation, with the continuity constraint (for incompressible flow). You also must specify boundary conditions, this is important and missing from your question. – Charles Apr 15 '16 at 17:57
• Also, you cannot say you're solving a 2D NS equation with one of the remaining components of velocity being zero, this is 1D NS. – Charles Apr 15 '16 at 18:00

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.