- #1

- 624

- 11

## Main Question or Discussion Point

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 math

L=1.0

gamma = 1.4

N=1000

h=L/(N-1)

CFL=0.1

v=[[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

from matplotlib import pyplot as plt

import numpy as np

import copy

import random

import math

L=1.0

gamma = 1.4

N=1000

h=L/(N-1)

CFL=0.1

v=[[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*

vv

*[1]=u*rho #and third for e, and a vector for the volume of each grid*

vv

*[2]=e*

# calculate the maximum speed in the program

def cMax():

uMax=0

for i in range(N):

pressure = (gamma-1)*v# calculate the maximum speed in the program

def cMax():

uMax=0

for i in range(N):

pressure = (gamma-1)*v

*[0]*v**[2]*

speed=vspeed=v

*[1]/v**[0]*

c=np.sqrt(gamma*abs(pressure)/abs(vc=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 vif 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 vrho=v[i-1][0]

if v

*[0]> v[i-1][0]:*

rho=vrho=v

*[0]*

if vif v

*[0]==v[i-1][0]:*

rho=0

Frho=0

F

*[0]=rho*v**[1]/v**[0]*

#energy flux

if v#energy flux

if v

*[2]<v[i-1][2]:*

e=v[i-1][2]

if ve=v[i-1][2]

if v

*[2]>v[i-1][2]:*

e=ve=v

*[2]*

if vif v

*[2]==v[i-1][2]:*

e=0

Fe=0

F

*[2]=F**[0]*e*

#momentum flux

if v#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 vu=v[i-1][1]/v[i-1][0]

if v

*[1]/v**[0]>v[i-1][1]/v[i-1][0]:*

u=vu=v

*[1]/v**[0]*

if vif v

*[1]/v**[0]==v[i-1][1]/v[i-1][0]:*

u=0

Fu=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=vinitialise()

def update():

tau=CFL*h/cMax()

fluxCalculator()

for i in range(1,N-1):

rho=v

*[0]*

m=vm=v

*[1]*

e=ve=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]-Fpressure=(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]-Fe=e-tau/h*(F[i+1][2]-F

*[2])*

m=m-tau/h*(F[i+1][1]-Fm=m-tau/h*(F[i+1][1]-F

*[1])*

vv

*[0]=rho*

vv

*[1]=m*

vv

*[2]=e*

density=[]

position=np.arange(0,1+h/2,h)

for i in range(1000):

update()

for i in range(N):

density=density+[vdensity=[]

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()plt.plot(position,density)

plt.show()