Hydrodynamics simulation

  • #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[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()
 

Answers and Replies

  • #2
33,306
4,998
Please repost your code, using a [code=python] tag at the top and a [/code] tag at the bottom.

You might notice that what is displayed above switches to italics in the middle somewhere. This is because your code contains v[i][0]. The first pair of brackets is interpreted by browsers as a sign that the following text should be in italics. Each [i] in the code is "gobbled up" and doesn't appear.
But when i plot density versus distance I just get something like a noisy signal. Any idea where I am doing something wrong?
Given the amount of code, this is probably not enough for us to go on. If you can identify a section of code that you think is misbehaving, that would be very helpful.
 
Last edited:
  • Like
Reactions: BvU
  • #3
BvU
Science Advisor
Homework Helper
2019 Award
13,050
3,024
Please repost your code, using a [code=python] tag at the top and a [/code] tag at the bottom.
Hehe, I wanted to point that out to Silviu in a private conversation, so his (her) thread would remain unanswered, but never mind. Out of curiosity: How can you let the tags
Python:
[code=python] and ##[\ /code]##  appear without the browser taking action
o0)???
 
  • #4
33,306
4,998
Hehe, I wanted to point that out to Silviu in a private conversation, so his (her) thread would remain unanswered, but never mind. Out of curiosity: How can you let the tags
Python:
[code=python] and ##[\ /code]##  appear without the browser taking action
o0)???
I change the color of the first character in each tag. Black works fine -- this causes the browser to not interpret the tag as it usually would.
 

Related Threads on Hydrodynamics simulation

  • Last Post
Replies
4
Views
9K
Replies
3
Views
535
Replies
2
Views
8K
Replies
4
Views
3K
Replies
2
Views
2K
Replies
6
Views
4K
  • Last Post
Replies
4
Views
877
  • Last Post
Replies
7
Views
817
  • Last Post
Replies
4
Views
241
  • Last Post
Replies
10
Views
951
Top