I am trying to use backwards finite difference method to numerically solve a pair of partial differential equations:
- $\frac{\partial \left(pv\right)}{\partial x}+\frac{\partial p}{\partial t}=0$
- $\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:
- $v\left(x,0 \right)=0$
- $p\left(x,0 \right)=p_a$
- $p\left(L,t \right)=p_a$
- $\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:
- ${p_t}_{i,j}=\frac{\left(p_{i,j}-p_{i,j-1}\right)}{\delta t}$
- ${p_x}_{i,j}=\frac{\left(p_{i,j}-p_{i-1,j}\right)}{\delta x}$
- ${v_t}_{i,j}=\frac{\left(v_{i,j}-v_{i,j-1}\right)}{\delta t}$
- ${v_x}_{i,j}=\frac{\left(v_{i,j}-v_{i-1,j}\right)}{\delta x}$
putting these in the PDEs:
$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$
$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:
- the code above runs without any errors, but the result is nothing but the initial condition!
- I'm obviously not able to incorporate the boundary condition 3 at the end
- I'm not able to calculate $v_{0,j}$s
and my questions are:
- are my equations correct?
- are my boundary conditions realistic and enough?
- am I implementing backwards finite difference method correctly?
- how to change my python code to solve the three issues above?