Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Python Hydrodynamics simulation

  1. Apr 18, 2016 #1
    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[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]
     
  2. jcsd
  3. Apr 18, 2016 #2

    Mark44

    Staff: Mentor

    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.
    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: Apr 19, 2016
  4. Apr 19, 2016 #3

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    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
    Code (Python):
    [code=python] and ##[\ /code]##  appear without the browser taking action
    o0)???
     
  5. Apr 19, 2016 #4

    Mark44

    Staff: Mentor

    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Hydrodynamics simulation
  1. Quantum Simulation (Replies: 5)

  2. Physics simulation (Replies: 5)

  3. Simulating the world (Replies: 14)

Loading...