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?(adsbygoogle = window.adsbygoogle || []).push({});

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

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[i][0]*e

#momentum flux

if v[i][1]/v[i][0]<v[i-1][1]/v[i-1][0]:

u=v[i-1][1]/v[i-1][0]

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

u=v[i][1]/v[i][0]

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

u=0

F[i][1]=0.5*u*(F[i+1][0]+F[i][0])

initialise()

def update():

tau=CFL*h/cMax()

fluxCalculator()

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

rho=v[i][0]

m=v[i][1]

e=v[i][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[i][0])

e=e-tau/h*(F[i+1][2]-F[i][2])

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

v[i][0]=rho

v[i][1]=m

v[i][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[i][0]]

plt.plot(position,density)

plt.show()[/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i][/i]

**Physics Forums | Science Articles, Homework Help, Discussion**

Dismiss Notice

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Python Hydrodynamics simulation

Tags:

Have something to add?

Draft saved
Draft deleted

**Physics Forums | Science Articles, Homework Help, Discussion**