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

Question

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?

Answer 1

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.