# How to implement finite difference method for one dimensional Navier-Stokes PDEs

I am trying to use backwards finite difference method to numerically solve a pair of partial differential equations:

1. $\frac{\partial \left(pv\right)}{\partial x}+\frac{\partial p}{\partial t}=0$
2. $\frac{\partial \left(pv^2\right)}{\partial x}+\frac{\partial \left(pv\right)}{\partial t}+c\frac{\partial p}{\partial x}=0$

with initial and boundary conditions:

1. $v\left(x,0 \right)=0$
2. $p\left(x,0 \right)=p_a$
3. $p\left(L,t \right)=p_a$
4. $\left\{\begin{matrix} if \, t<t_s \, then & p\left(0,t \right)=p_a \, and \, v\left(0,t \right)=0 \\ else & p\left(0,t \right)=p_r \end{matrix}\right.$

These are the Navier-Stokes equations for a one dimensional compressible inviscid isothermal fluid in a tube with one open end and step pressure signal at $t=t_s$ in the other end. gravity, radiation and conduction are negligible.

Using a backwards finite difference method we have:

1. ${p_t}_{i,j}=\frac{\left(p_{i,j}-p_{i,j-1}\right)}{\delta t}$
2. ${p_x}_{i,j}=\frac{\left(p_{i,j}-p_{i-1,j}\right)}{\delta x}$
3. ${v_t}_{i,j}=\frac{\left(v_{i,j}-v_{i,j-1}\right)}{\delta t}$
4. ${v_x}_{i,j}=\frac{\left(v_{i,j}-v_{i-1,j}\right)}{\delta x}$

putting these in the PDEs:

1. $v_{i,j}\frac{\left(p_{i,j}-p_{i-1,j}\right)}{\delta x}+ p_{i,j}\frac{\left(v_{i,j}-v_{i-1,j}\right)}{\delta x}+ \frac{\left(p_{i,j}-p_{i,j-1}\right)}{\delta t}=0$

2. $v_{i,j}^2\frac{\left(p_{i,j}-p_{i-1,j}\right)}{\delta x}+ v_{i,j}\frac{\left(p_{i,j}-p_{i,j-1}\right)}{\delta t}+ 2v_{i,j}p_{i,j}\frac{\left(v_{i,j}-v_{i-1,j}\right)}{\delta x}+ p_{i,j}\frac{\left(v_{i,j}-v_{i,j-1}\right)}{\delta t}+ c\frac{\left(p_{i,j}-p_{i-1,j}\right)}{\delta x}=0$

I have applied these in a very primitive python code:

import numpy
import math
from scipy.optimize import fsolve

l=1.0 #Length of the tube in meters
t=10.0 # length of the simulation in seconds
M=10 # spatial discretisation resolution
N=4 # temporal discretisation resolution
pr=200 #reservoir pressure pascal
pa=100 #ambient pressure pascal
Ta=300 #ambient temperature in kelvin
R=286.9 #Gas Constant for air

#step signal
ts=1 #in seconds,  after ts p start will be pr

dx=l/(M-1) #meters
dt=t/(N-1) #meters
p=numpy.ones((M,N))*pa
v=numpy.zeros((M,N))
p[0,(math.ceil(ts/dt)):N]=numpy.ones(N-(math.ceil(ts/dt)))*pr #assign the start boundry condition as step signal

def eqns(ivars, *data):
pij, vij=ivars
pij_1,pi_1j,vij_1,vi_1j=data
return (vij*(pij-pi_1j)/dx+pij*(vij-vi_1j)/dx+(pij-pij_1)/dt, 2*pij*vij*(vij-vi_1j)/dx-(vij**2)*(pij-pi_1j)/dx+vij*(pij-pij_1)/dt+pij*(vij-vij_1)/dt+R*Ta*(pij-pi_1j)/dx)

for i in range(1,M):
for j in range(1,N):
data=(p[i][j-1],p[i-1][j],v[i][j-1],v[i-1][j])
p[i][j], v[i][j]=fsolve(eqns,(p[i][j-1],v[i][j-1]), args=data)

print(v)
print(p)


I have three issues:

1. the code above runs without any errors, but the result is nothing but the initial condition!
2. I'm obviously not able to incorporate the boundary condition 3 at the end
3. I'm not able to calculate $v_{0,j}$s

and my questions are:

1. are my equations correct?
2. are my boundary conditions realistic and enough?
3. am I implementing backwards finite difference method correctly?
4. how to change my python code to solve the three issues above?
• Your equations are correct if you are able to add heat along the tube so that it keeps its temperature constant (I suppose c=Rg*T in your equations). However, I think this is difficult to achieve in practice. In that case, you should use the full Euler equations and consider that some expansion and shock waves may appear, and thus pressure and temperature will not be continuous along all the tube. Have a look at this paper that talks about the Riemann's problem (which is not exactly the same as yours, but similar): astro.uu.se/~hoefner/astro/teach/ch10.pdf – Manu Aug 7 '17 at 19:09
• In that paper talks about how to numerically solve this kind of problems. In all cases they used centered schemes or combine forward and backward finite differences, so the boundary conditions you impose at the two sides of the tube are effectively applied. – Manu Aug 7 '17 at 19:41
• Btw, how do you impose pr(t) at x=0? Is it by the motion of a piston? How is that motion? If its displacements are great, then the fluid domain changes and it should be necessary to redefine the mesh at every instant or to use a lagrangian formulation. – Manu Aug 7 '17 at 20:07
• @Manu thanks for your comments: 1. indead $c=\mathring{R} T_a$ 2. Do you mean the isothermal assumption makes the PDE insoluble? It actually makes sense because I can't solve the PDE using Mathematica either see this 3. I will go through the paper. thanks. – Foad Aug 7 '17 at 22:41
• 4. since I have posted this I have realised that finite difference has issues with boundary condition. I still need to learn how to tackle that. maybe central difference would help? 5. there is a valve at the beginning of the tube. it is not really a step signal, but we can assume that for now. – Foad Aug 7 '17 at 22:41

First assume that the flow is adiabatic and not isothermal, since the change of pressure at the left is a very fast process and thus the flow does not have time to exchange heat enough with the surroundings to be at its same temperature.

Then, the problem is governed by the Euler equations, which can be written in conservative form as \begin{equation*} \dfrac{\partial}{\partial t}\begin{bmatrix} \rho \\ \rho u \\ \frac{p}{\gamma-1} + \frac{1}{2}\rho u^2 \end{bmatrix} + \dfrac{\partial}{\partial x}\begin{bmatrix} \rho u \\ p + \rho u^2 \\ \frac{\gamma}{\gamma-1} pu + \frac{1}{2}\rho u^3 \end{bmatrix} = 0. \end{equation*}

Suppose we know the values of $p$, $\rho$ and $u$ at each point $i=1,\ldots,N$ of the grid at an instant $t^n$. With any of the (explicit) schemes shown in http://www.astro.uu.se/~hoefner/astro/teach/ch10.pdf, it is possible to compute $p$, $\rho$ and $u$ in the center points $i=2,\ldots,N-1$ at the following instant $t^{n+1}$.

The question now is how to impose the boundary conditions. It can be shown that, at each point, there are three magnitudes ---the so-called Riemann invariants--- that keep their value constant along some characteristics. These invariants are \begin{equation*} \zeta_1 = \dfrac{p}{\rho^\gamma}, \quad \zeta_2 = u + \dfrac{2c}{\gamma-1}, \quad \zeta_3 = u - \dfrac{2c}{\gamma-1}, \end{equation*} which propagate at speeds $u$, $u+c$ and $u-c$, respectively. Thus, in the case of a subsonic inlet (as ours), in which $0<u<c$, there are two characteristics coming from the outside of the tube ---which correspond to $\zeta_1$ and $\zeta_2$--- and one coming from the inside ---$\zeta_3$---. Since the characteristics can be seen as waves that propagate information, it is deduced that we must impose two boundary conditions at $x=0$ (pressure input and conservation of the stagnation enthalpy), whereas the other condition that determines the three unknowns ($p_1^{n+1},\rho_1^{n+1},u_1^{n+1}$) is obtained from the invariant that arrives from the inside of the tube ($\zeta_3$). According to http://inis.jinr.ru/sl/Simulation/Hirsch,_Numerical_Computation_of_Internal&External_Flows,1994/Hirsch,_Numerical_Computation_of_Internal&External_Flows,v2,1994/ch19.pdf, this condition is implemented by

• Imposing a compatibility equation. In the book, the compatibility equations have been linearized, but in the end is the same as imposing (in this case) that \begin{equation*} \dfrac{\partial \zeta_3}{\partial t} + (u-c) \dfrac{\partial \zeta_3}{\partial x} = 0. \end{equation*} Thus, the sought values of $p_1^{n+1},\rho_1^{n+1},u_1^{n+1}$ are the solutions of the following system of equations \begin{align*} p_1^{n+1}-p_r(t^{n+1}) &= 0, \\ c_p T_1^{n+1} + \dfrac{u^2}{2} - c_p T_r &= 0, \\ \dfrac{(\zeta_3)_1^{n+1} - (\zeta_3)_1^{n}}{\Delta t} + (u_1^n-c_1^n) \dfrac{(\zeta_3)_2^{n} - (\zeta_3)_1^{n}}{\Delta x} &= 0. \end{align*} Note that we use forward differencing in the last formula since the information comes precisely from the interior of the tube.
• Using an extrapolation technique. In this case, $p_1^{n+1},\rho_1^{n+1},u_1^{n+1}$ are the solutions of: \begin{align*} p_1^{n+1}-p_r(t^{n+1}) &= 0, \\ c_p T_1^{n+1} + \dfrac{u^2}{2} - c_p T_r &= 0, \\ (\zeta_3)_1^{n+1}-2(\zeta_3)_2^{n+1}+(\zeta_3)_3^{n+1} &= 0. \end{align*}

For the outlet, we have to consider two cases. If the flow is subsonic at the exit, we will have two characteristics coming from the inside of the tube (those corresponding to $\zeta_1$ and $\zeta_2$) and one coming from the outside. Thus, in this case we must impose only one boundary condition ($p_N^{n+1}=p_a$). The other two equations needed for determining $\rho_N^{n+1}$ and $u_N^{n+1}$ are obtained from ''conservation of the invariants'', that is, imposing \begin{align*} \dfrac{(\zeta_1)_N^{n+1} - (\zeta_1)_N^{n}}{\Delta t} + u_N^n \dfrac{(\zeta_1)_N^{n} - (\zeta_1)_{N-1}^{n}}{\Delta x} &= 0, \\ \dfrac{(\zeta_2)_N^{n+1} - (\zeta_2)_N^{n}}{\Delta t} + (u_N^n+c_N^n) \dfrac{(\zeta_2)_N^{n} - (\zeta_2)_{N-1}^{n}}{\Delta x} &= 0. \end{align*} (note that backward differences are used now) or \begin{align*} (\zeta_1)_{N-2}^{n+1}-2(\zeta_1)_{N-1}^{n+1}+(\zeta_1)_N^{n+1} &= 0, \\ (\zeta_2)_{N-2}^{n+1}-2(\zeta_2)_{N-1}^{n+1}+(\zeta_2)_N^{n+1} &= 0. \end{align*} However, if the flow is supersonic at the exit, we have three characteristics leaving the tube, so no information comes from the outside (choked tube) and therefore we cannot impose any boundary conditions. In that case, we should proceed in a similar way but imposing the values of the three Riemann invariants.

Some methods may present some spurious oscillations at the beginning. I'd suggest (i) not to impose a step $p_r(t)$ but a arctangent one (or similar) and (ii) if $v_1^{n+1}<0$ (or $v_N^{n+1}<0$) at some of the beginning instants, set it to zero or take $v_1^{n+1}=v_2^{n+1}$, otherwise the boundary conditions should change.

Take into account that these methods require $\Delta t/\Delta x$ to verify the CFL condition to be stable. You can also try implicit methods that are more stable but have to solve a big system of nonlinear equations at every instant. In those cases, it is important (from the performance point of view) to provide the solver an expression of the jacobian ---which is sparse, and this is also very important--- or, at least, its sparsity pattern. I do not know anything about python, but it seems for me that you were trying an implicit scheme without providing the solver any information about the jacobian.

• this is great. thanks a lot. it will take me a while to implement this though. Have you tried implementing it yourself? – Foad Aug 18 '17 at 21:24
• Yes, I have, but in Matlab. I can send you the codes if you are interested. Btw, I think $p_a$ should be $10^5\,Pa$, and not $100\,Pa$. – Manu Aug 19 '17 at 10:32
• how about sharing the code in a github snippet? – Foad Aug 19 '17 at 18:31
• – Manu Aug 23 '17 at 19:51
• Thanks @Manu , At the moment I'm doing lots of experiments to develop some empirical models for my problem. later I'll come back and use this. – Foad Aug 23 '17 at 20:52