- #1
Silviu
- 624
- 11
Hi! I have this code which should simulate a sod shock tube like this: https://en.wikipedia.org/wiki/Sod_shock_tube But when i plot density versus distance I just get something like a noisy signal. Any idea where I am doing something wrong?
from matplotlib import pyplot as plt
import numpy as np
import copy
import random
import mathL=1.0
gamma = 1.4
N=1000
h=L/(N-1)
CFL=0.1v=[[0 for x in range(3)] for x in range(N)]
F=[[0 for x in range(3)] for x in range(N)]
# give the values to each cell
def initialise():
for i in range(N):
if i<N/2:
rho=1
p=1
u=0
else:
rho=0.125
p=0.1
u=0
e=p/((gamma-1)*rho)
v[0]=rho #create a matrix with first column for density, second for density times speed
v[1]=u*rho #and third for e, and a vector for the volume of each grid
v[2]=e
# calculate the maximum speed in the program
def cMax():
uMax=0
for i in range(N):
pressure = (gamma-1)*v[0]*v[2]
speed=v[1]/v[0]
c=np.sqrt(gamma*abs(pressure)/abs(v[0]))
if uMax< (c+abs(speed)):
uMax=c+abs(speed)
return uMax
def fluxCalculator():
for i in range(1,N-1):
#density flux
if v[0]< v[i-1][0]:
rho=v[i-1][0]
if v[0]> v[i-1][0]:
rho=v[0]
if v[0]==v[i-1][0]:
rho=0
F[0]=rho*v[1]/v[0]
#energy flux
if v[2]<v[i-1][2]:
e=v[i-1][2]
if v[2]>v[i-1][2]:
e=v[2]
if v[2]==v[i-1][2]:
e=0
F[2]=F[0]*e
#momentum flux
if v[1]/v[0]<v[i-1][1]/v[i-1][0]:
u=v[i-1][1]/v[i-1][0]
if v[1]/v[0]>v[i-1][1]/v[i-1][0]:
u=v[1]/v[0]
if v[1]/v[0]==v[i-1][1]/v[i-1][0]:
u=0
F[1]=0.5*u*(F[i+1][0]+F[0])initialise()
def update():
tau=CFL*h/cMax()
fluxCalculator()
for i in range(1,N-1):
rho=v[0]
m=v[1]
e=v[2]
pressure=(gamma-1)*rho*e
e=e-tau/h*pressure*(v[i+1][1]/v[i+1][0]-m/rho)
m=m-tau/h*(pressure-(gamma-1)*v[i-1][0]*v[i-1][2])
rho=rho-tau/h*(F[i+1][0]-F[0])
e=e-tau/h*(F[i+1][2]-F[2])
m=m-tau/h*(F[i+1][1]-F[1])
v[0]=rho
v[1]=m
v[2]=e
density=[]
position=np.arange(0,1+h/2,h)
for i in range(1000):
update()
for i in range(N):
density=density+[v[0]]
plt.plot(position,density)
plt.show()
from matplotlib import pyplot as plt
import numpy as np
import copy
import random
import mathL=1.0
gamma = 1.4
N=1000
h=L/(N-1)
CFL=0.1v=[[0 for x in range(3)] for x in range(N)]
F=[[0 for x in range(3)] for x in range(N)]
# give the values to each cell
def initialise():
for i in range(N):
if i<N/2:
rho=1
p=1
u=0
else:
rho=0.125
p=0.1
u=0
e=p/((gamma-1)*rho)
v[0]=rho #create a matrix with first column for density, second for density times speed
v[1]=u*rho #and third for e, and a vector for the volume of each grid
v[2]=e
# calculate the maximum speed in the program
def cMax():
uMax=0
for i in range(N):
pressure = (gamma-1)*v[0]*v[2]
speed=v[1]/v[0]
c=np.sqrt(gamma*abs(pressure)/abs(v[0]))
if uMax< (c+abs(speed)):
uMax=c+abs(speed)
return uMax
def fluxCalculator():
for i in range(1,N-1):
#density flux
if v[0]< v[i-1][0]:
rho=v[i-1][0]
if v[0]> v[i-1][0]:
rho=v[0]
if v[0]==v[i-1][0]:
rho=0
F[0]=rho*v[1]/v[0]
#energy flux
if v[2]<v[i-1][2]:
e=v[i-1][2]
if v[2]>v[i-1][2]:
e=v[2]
if v[2]==v[i-1][2]:
e=0
F[2]=F[0]*e
#momentum flux
if v[1]/v[0]<v[i-1][1]/v[i-1][0]:
u=v[i-1][1]/v[i-1][0]
if v[1]/v[0]>v[i-1][1]/v[i-1][0]:
u=v[1]/v[0]
if v[1]/v[0]==v[i-1][1]/v[i-1][0]:
u=0
F[1]=0.5*u*(F[i+1][0]+F[0])initialise()
def update():
tau=CFL*h/cMax()
fluxCalculator()
for i in range(1,N-1):
rho=v[0]
m=v[1]
e=v[2]
pressure=(gamma-1)*rho*e
e=e-tau/h*pressure*(v[i+1][1]/v[i+1][0]-m/rho)
m=m-tau/h*(pressure-(gamma-1)*v[i-1][0]*v[i-1][2])
rho=rho-tau/h*(F[i+1][0]-F[0])
e=e-tau/h*(F[i+1][2]-F[2])
m=m-tau/h*(F[i+1][1]-F[1])
v[0]=rho
v[1]=m
v[2]=e
density=[]
position=np.arange(0,1+h/2,h)
for i in range(1000):
update()
for i in range(N):
density=density+[v[0]]
plt.plot(position,density)
plt.show()