Help regarding 1D compressible flow (Sod Shock tube)

selig5576 · July 25, 2016, 17:22

Hi all,

I am coding up a 1D compressible flow solver via Lax-Friedrichs method (for now), however have run into a when calculating certain quantities, such as velocity. In other words, I believe my implementation is correct, however something is not being properly updated per time step. I have attached my code for proper documentation.

Thoughts: My way of updating at each time step is to have U_current = U_new.

EDIT: Just for clarification, I am not sure what I am doing wrong as my plots for pressure, velocity, and density are not right

selig5576 · July 26, 2016, 11:40

Hi All,

I am still having issues with my 1D compressible solver. I have narrowed it down to the time step, however I am not sure what I'm doing wrong in the process of updating it. If someone point/help me in the right direction I would heavily appreciate it.

Code:

clear; 
close all; 
clc;

gamma = 1.4;
b = 0.5;
a = -0.5;

nx = 101;
dx = (b-a)/(nx-1); 
x = zeros(nx,1);

tMax = 0.15;
dt = 0.004;
nt = round(tMax/dt)+1;
t = zeros(nt,1);

r = zeros(nx,nt);
p = zeros(nx,nt);
u = zeros(nx,nt);

F1 = zeros(nx,nt);
F2 = zeros(nx,nt);
F3 = zeros(nx,nt);

U1 = zeros(nx,nt);
U2 = zeros(nx,nt);
U3 = zeros(nx,nt);

for j=1:nt
    for i=1:nx
        x(i) = a + (i-1)*dx;
        t(j) = (j-1)*dt;
    end
end


for i=1:nx
    if (x(i) <= 0)
        r(i,1) = 1;
        u(i,1) = 0;
        p(i,1) = 1;
    elseif (x(i) > 0)
        r(i,1) = 0.125;
        u(i,1) = 0;
        p(i,1) = 0.1;
    end
end   

U1(1) = 1;
U1(nx) = 0.125;

U2(1) = 0;
U2(nx) = 0;

U3(1) = 1;
U3(nx) = 0.1;

E = p./(gamma-1) + 0.5*r.*u.^2;

for j=1:nt
    for i=1:nx
        U1(i,j) = r(i,j);
        U2(i,j) = r(i,j).*u(i,j);
        U3(i,j) = E(i,j);
        
        F1(i,j) = r(i,j).*u(i,j);
        F2(i,j) = r(i,j).*u(i,j).^2 + p(i,j);
        F3(i,j) = u(i,j).*(E(i,j) + p(i,j));
    end

    for i=2:nx-1    
        U1(i,j+1) = 0.5*(U1(i+1,j) + U1(i-1,j)) - dt/(2*dx)*(F1(i+1,j) - F1(i-1,j));
        U2(i,j+1) = 0.5*(U2(i+1,j) + U2(i-1,j)) - dt/(2*dx)*(F2(i+1,j) - F2(i-1,j));
        U3(i,j+1) = 0.5*(U3(i+1,j) + U3(i-1,j)) - dt/(2*dx)*(F3(i+1,j) - F3(i-1,j));
    end
    
    r = U1;
    u = U2./U1;
    E = U3; 
    p = (gamma-1)*(E-0.5*r.*u.^2);
end
u
s1 = subplot(2,2,1); 
plot(x,r,'-b'); 
xlabel('x'); 
ylabel('r');
title('Density');

s2 = subplot(2,2,2); 
plot(x,u,'-b'); 
xlabel('x'); 
ylabel('u'); 
title('Velocity');

s3 = subplot(2,2,3);
plot(x,p,'-b'); 
xlabel('x'); 
ylabel('p');
title('Pressure');

I have attached an image as well.

agd · July 26, 2016, 12:29

I am not completely familiar with the language you are using, but it looks to me like you are overwriting values during your update that should not be overwritten until you complete the entire sweep. You construct u(i,j+1) and then use that value for u(i,j) in the next pass through the j-loop.

selig5576 · July 26, 2016, 12:39

Hi agd,

I am using Matlab for my language. In terms of

Code:

for j=1:nt
    for i=1:nx
        U1(i,j) = r(i,j);
        U2(i,j) = r(i,j).*u(i,j);
        U3(i,j) = E(i,j);
        
        F1(i,j) = r(i,j).*u(i,j);
        F2(i,j) = r(i,j).*u(i,j).^2 + p(i,j);
        F3(i,j) = u(i,j).*(E(i,j) + p(i,j));
    end

    for i=2:nx-1    
        U1(i,j+1) = 0.5*(U1(i+1,j) + U1(i-1,j)) - dt/(2*dx)*(F1(i+1,j) - F1(i-1,j));
        U2(i,j+1) = 0.5*(U2(i+1,j) + U2(i-1,j)) - dt/(2*dx)*(F2(i+1,j) - F2(i-1,j));
        U3(i,j+1) = 0.5*(U3(i+1,j) + U3(i-1,j)) - dt/(2*dx)*(F3(i+1,j) - F3(i-1,j));
    end
    
    r = U1;
    u = U2./U1;
    E = U3; 
    p = (gamma-1)*(E-0.5*r.*u.^2);
end

Are you saying that specigying the density, velocity, energy, and pressure, is causing issues? Should I make a separate loop for them in terms of space and time? A little more details to help: j represents time and i is space. I noticed that when I have r = U1, u = U2./U1, etc. I get "NaNs" sprinkled in the arrays of r, u, E, and p.

P.S. I'm new to CFD

agd · July 26, 2016, 14:00

OK - I see what your code is doing now. I'll look at it a little more and see if anything looks wrong. One thing you can try for debugging is to turn off the boundary conditions and see if your code will maintain a uniform set of conditions. Sometimes that can quickly point to a problem.

selig5576 · July 26, 2016, 14:49

Hi agd (and all),

So I took the boundary conditions out (commented it out) and nothing happens. I've looked at the arrays and I see a lot of "Nans." That being said, I removed:

Code:

    r = U1;
    u = U2./U1;
    E = U3; 
    p = (gamma-1)*(E-0.5*r.*u.^2);

and the Nans get removed, however all the values of the array are 0. As such, still no luck.

Below is the hyperbolic system (compressible euler) and some physical definitions I use for solving for total energy and pressure:

$\frac{\partial \rho}{\partial t} = \frac{\partial \rho u}{\partial x}$

$\frac{\partial \rho u}{\partial t} = \frac{\partial \rho u^{2} + p}{\partial x}$

$\frac{\partial E}{\partial t} = \frac{\partial u(E+p)}{\partial x}$

$E = \frac{p}{\gamma -1} + \frac{1}{2} \rho u.^{2}$

$p = (\gamma - 1)(E - \frac{1}{2} \rho u^{2})$

The scheme I'm using is the Lax-Friedrichs, i.e.

$u_{i}^{n+1} = \frac{1}{2} \left(u_{i+1}^{n} + u_{i-1}^{n}\right) - \frac{\Delta t}{2 \Delta x} \left(F_{i+1}^{n} - F_{i-1}^{n}\right)$

agd · July 27, 2016, 10:53

The conservation equations as you have written them are not correct. Is that simply an error in how you wrote them in your post, or did you code them up that way?

selig5576 · July 27, 2016, 11:27

Hi,

I coded them up two different ways to see if what you said was an issue, however the issue does not change.

Code:

clear; 
close all; 
clc;

gamma = 1.4;
b = 0.5;
a = -0.5;

nx = 101;
dx = (b-a)/(nx-1); 

tMax = 0.15;
dt = 0.004;
nt = round(tMax/dt)+1;

x = zeros(nx,nt);
t = zeros(nt,nt);

r = zeros(nx,nt);
p = zeros(nx,nt);
u = zeros(nx,nt);

F1 = zeros(nx,nt);
F2 = zeros(nx,nt);
F3 = zeros(nx,nt);

U1 = zeros(nx,nt);
U2 = zeros(nx,nt);
U3 = zeros(nx,nt);

for j=1:nt
    for i=1:nx
        x(i,j) = a + (i-1)*dx;
        t(i,j) = (j-1)*dt;
    end
end


for i=1:nx
    if (x(i) <= 0)
        r(i,1) = 1.0;
        u(i,1) = 0.0;
        p(i,1) = 1.0;
    elseif (x(i) > 0)
        r(i,1) = 0.125;
        u(i,1) = 0.0;
        p(i,1) = 0.1;
    end
end   

E = p./((gamma-1).*r) + 0.5*r.*u.^2;

U1(1) = 1;
U1(nx) = 0.125;

U2(1) = 0;
U2(nx) = 0;

U3(1) = 1;
U3(nx) = 0.1;

for j=1:nt
    for i=1:nx
        U1(i,j) = r(i,j);
        U2(i,j) = r(i,j)*u(i,j);
        U3(i,j) = r(i,j).*E(i,j);
        
        F1(i,j) = r(i,j)*u(i,j);
        F2(i,j) = r(i,j)*u(i,j).^2 + p(i,j);
        F3(i,j) = u(i,j)*(r(i,j).*E(i,j) + p(i,j));
    end

    for i=2:nx-1    
        U1(i,j+1) = 0.5*(U1(i+1,j) + U1(i-1,j)) - (dt/(2*dx))*(F1(i+1,j) - F1(i-1,j));
        U2(i,j+1) = 0.5*(U2(i+1,j) + U2(i-1,j)) - (dt/(2*dx))*(F2(i+1,j) - F2(i-1,j));
        U3(i,j+1) = 0.5*(U3(i+1,j) + U3(i-1,j)) - (dt/(2*dx))*(F3(i+1,j) - F3(i-1,j));
    end
    
    r = U1;
    u = U2./r;
    E = U3./r;
    p = (gamma-1).*r.*(E-0.5*r.*u.^2);
end

subplot(2,2,1); 
plot(x,r(:,26),'-b'); 
xlabel('x'); 
ylabel('r');
title('Density');

subplot(2,2,2); 
plot(x,u(:,26),'-b'); 
xlabel('x'); 
ylabel('u'); 
title('Velocity');

subplot(2,2,3);
plot(x,p(:,26),'-r'); 
xlabel('x'); 
ylabel('p');
title('Pressure');

The code aboves solves
$\frac{\partial \rho}{\partial t} + \frac{\partial \rho u}{\partial x} = 0$
$\frac{\partial \rho u}{\partial t} + \frac{\partial (\rho u.^{2} + p)}{\partial x} = 0$
$\frac{\partial \rho E}{\partial t} + \frac{\partial u (\rho E + p)}{\partial x} = 0$

Link: http://bender.astro.sunysb.edu/hydro...ible/Euler.pdf

The link I provided above is to illustrate the equations I'm solving not the method.

The part in which really stumps me (and maybe you can have wisdom on that) is why I keep getting NaNs in my arrays. I spent a couple hours debugging with no luck. I believe if I did not get "NaNs" in my arrays I would have a complete solution as it looks very close.

FMDenaro · July 27, 2016, 11:36

the RHS of your system must have the minus sign...

FMDenaro · July 27, 2016, 11:45

for i=1:nx U1(i,j) = r(i,j); U2(i,j) = r(i,j)*u(i,j); U3(i,j) = r(i,j).*E(i,j); F1(i,j) = r(i,j)*u(i,j); F2(i,j) = r(i,j)*u(i,j).^2 + p(i,j); F3(i,j) = u(i,j)*(r(i,j).*E(i,j) + p(i,j)); end

why do you mix the operators "*" and ".*" in the for cycle? they are different operations in matlab

selig5576 · July 27, 2016, 11:48

Hi,

Thanks for pointing that out. That was a mistake in terms of not putting the - sign, however in the code I already took that into account. In terms of wrongly using .* and * I have fixed that and is below.

Code:

clear; 
close all; 
clc;

gamma = 1.4;
b = 0.5;
a = -0.5;

nx = 101;
dx = (b-a)/(nx-1); 

tMax = 0.15;
dt = 0.004;
nt = round(tMax/dt)+1;

x = zeros(nx,nt);
t = zeros(nt,nt);

r = zeros(nx,nt);
p = zeros(nx,nt);
u = zeros(nx,nt);

F1 = zeros(nx,nt);
F2 = zeros(nx,nt);
F3 = zeros(nx,nt);

U1 = zeros(nx,nt);
U2 = zeros(nx,nt);
U3 = zeros(nx,nt);

for j=1:nt
    for i=1:nx
        x(i,j) = a + (i-1)*dx;
        t(i,j) = (j-1)*dt;
    end
end


for i=1:nx
    if (x(i) <= 0)
        r(i,1) = 1.0;
        u(i,1) = 0.0;
        p(i,1) = 1.0;
    elseif (x(i) > 0)
        r(i,1) = 0.125;
        u(i,1) = 0.0;
        p(i,1) = 0.1;
    end
end   

E = p./((gamma-1)*r) + 0.5*u.^2;

U1(1) = 1;
U1(nx) = 0.125;

U2(1) = 0;
U2(nx) = 0;

U3(1) = 1;
U3(nx) = 0.1;

for j=1:nt
    for i=1:nx
        U1(i,j) = r(i,j);
        U2(i,j) = r(i,j).*u(i,j);
        U3(i,j) = r(i,j).*E(i,j);
        
        F1(i,j) = r(i,j).*u(i,j);
        F2(i,j) = r(i,j).*u(i,j).^2 + p(i,j);
        F3(i,j) = u(i,j).*(r(i,j).*E(i,j) + p(i,j));
    end

    for i=2:nx-1    
        U1(i,j+1) = 0.5*(U1(i+1,j) + U1(i-1,j)) - (dt/(2*dx))*(F1(i+1,j) - F1(i-1,j));
        U2(i,j+1) = 0.5*(U2(i+1,j) + U2(i-1,j)) - (dt/(2*dx))*(F2(i+1,j) - F2(i-1,j));
        U3(i,j+1) = 0.5*(U3(i+1,j) + U3(i-1,j)) - (dt/(2*dx))*(F3(i+1,j) - F3(i-1,j));
    end
    
    r = U1;
    u = U2./r;
    E = U3./r;
    p = (gamma-1)*r.*(E-0.5*u.^2);
end

subplot(2,2,1); 
plot(x,r(:,26),'-b'); 
xlabel('x'); 
ylabel('r');
title('Density');

subplot(2,2,2); 
plot(x,u(:,26),'-b'); 
xlabel('x'); 
ylabel('u'); 
title('Velocity');

subplot(2,2,3);
plot(x,p(:,26),'-r'); 
xlabel('x'); 
ylabel('p');
title('Pressure');

With that said, I still get the same results.

agd · July 28, 2016, 11:45

Stupid question - are you sure that you are obeying the stability limit for the Lax-Friedrichs scheme?

selig5576 · July 28, 2016, 11:58

Hi agd,

Yes I am. I've made sure there is no instability aside the downfalls of the method.

FMDenaro · July 28, 2016, 12:04

Quote:

Originally Posted by agd

Stupid question - are you sure that you are obeying the stability limit for the Lax-Friedrichs scheme?

that's never a stupid question

he wrote an arbitrary dt and the cfl does not appear explicitly computed:

b = 0.5; a = -0.5; nx = 101; dx = (b-a)/(nx-1); x = zeros(nx,1); tMax = 0.15; dt = 0.004;

FMDenaro · July 28, 2016, 12:06

remember to estimate that dt/dx = 0.4 and the velocity should be considered as convective and sound velocity

selig5576 · August 2, 2016, 10:22

Hi agd and FMDenaro,

I have managed to get my 1d compressible flow solver to work. There were some things I did different:

1. I used the Local Lax Friedrichs which provided a more robust method. There were no small oscillations
2. The indexing I was using was running me into trouble, i.e. the time index
3. I adjusted my CFL condition

Thank you for the help

silent_hunter · January 20, 2020, 17:12

Hello!!Could you post the right code after your changes?

javeria · August 25, 2022, 05:11

hello i need lax friedrich matlab codes for 1D shock tube problem. Is there any one to help me?

javeria · August 25, 2022, 05:12

I need your codes

javeria · August 25, 2022, 05:16

hello selig can you post your codes after correction i need your codes

July 26, 2016, 12:29		#3
agd Senior Member Join Date: Jul 2009 Posts: 358 Rep Power: 19	I am not completely familiar with the language you are using, but it looks to me like you are overwriting values during your update that should not be overwritten until you complete the entire sweep. You construct u(i,j+1) and then use that value for u(i,j) in the next pass through the j-loop. selig5576 likes this.

July 26, 2016, 12:39	Thanks	#4
selig5576 Senior Member Selig Join Date: Jul 2016 Posts: 213 Rep Power: 11	Hi agd, I am using Matlab for my language. In terms of Code: for j=1:nt for i=1:nx U1(i,j) = r(i,j); U2(i,j) = r(i,j).u(i,j); U3(i,j) = E(i,j); F1(i,j) = r(i,j).u(i,j); F2(i,j) = r(i,j).u(i,j).^2 + p(i,j); F3(i,j) = u(i,j).(E(i,j) + p(i,j)); end for i=2:nx-1 U1(i,j+1) = 0.5(U1(i+1,j) + U1(i-1,j)) - dt/(2dx)(F1(i+1,j) - F1(i-1,j)); U2(i,j+1) = 0.5(U2(i+1,j) + U2(i-1,j)) - dt/(2dx)(F2(i+1,j) - F2(i-1,j)); U3(i,j+1) = 0.5(U3(i+1,j) + U3(i-1,j)) - dt/(2dx)(F3(i+1,j) - F3(i-1,j)); end r = U1; u = U2./U1; E = U3; p = (gamma-1)(E-0.5r.u.^2); end Are you saying that specigying the density, velocity, energy, and pressure, is causing issues? Should I make a separate loop for them in terms of space and time? A little more details to help: j represents time and i is space. I noticed that when I have r = U1, u = U2./U1, etc. I get "NaNs" sprinkled in the arrays of r, u, E, and p. P.S. I'm new to CFD

July 26, 2016, 14:49	Further debugging	#6
selig5576 Senior Member Selig Join Date: Jul 2016 Posts: 213 Rep Power: 11	Hi agd (and all), So I took the boundary conditions out (commented it out) and nothing happens. I've looked at the arrays and I see a lot of "Nans." That being said, I removed: Code: r = U1; u = U2./U1; E = U3; p = (gamma-1)(E-0.5r.*u.^2); and the Nans get removed, however all the values of the array are 0. As such, still no luck. Below is the hyperbolic system (compressible euler) and some physical definitions I use for solving for total energy and pressure: $\frac{\partial \rho}{\partial t} = \frac{\partial \rho u}{\partial x}$ $\frac{\partial \rho u}{\partial t} = \frac{\partial \rho u^{2} + p}{\partial x}$ $\frac{\partial E}{\partial t} = \frac{\partial u(E+p)}{\partial x}$ $E = \frac{p}{\gamma -1} + \frac{1}{2} \rho u.^{2}$ $p = (\gamma - 1)(E - \frac{1}{2} \rho u^{2})$ The scheme I'm using is the Lax-Friedrichs, i.e. $u_{i}^{n+1} = \frac{1}{2} \left(u_{i+1}^{n} + u_{i-1}^{n}\right) - \frac{\Delta t}{2 \Delta x} \left(F_{i+1}^{n} - F_{i-1}^{n}\right)$

July 28, 2016, 11:58	CFL condition	#13
selig5576 Senior Member Selig Join Date: Jul 2016 Posts: 213 Rep Power: 11	Hi agd, Yes I am. I've made sure there is no instability aside the downfalls of the method.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Issues on the simulation of high-speed compressible flow within turbomachinery	dowlee	OpenFOAM Running, Solving & CFD	11	August 6, 2021 07:40
Compressible Flow Pressure Outlet Back Flow	okstatecheme	OpenFOAM Running, Solving & CFD	1	March 21, 2018 09:11
rhoCentralFoam not reflecting shock in Shock Tube?	Astaria	OpenFOAM Running, Solving & CFD	5	March 4, 2012 04:07
compressible flow	maria teresa	FLUENT	1	September 7, 2007 17:58
Inviscid Drag at subsonic, subcritical Mach #	Axel Rohde	Main CFD Forum	1	November 19, 2001 13:19

July 26, 2016, 14:00		#5
agd Senior Member Join Date: Jul 2009 Posts: 358 Rep Power: 19	OK - I see what your code is doing now. I'll look at it a little more and see if anything looks wrong. One thing you can try for debugging is to turn off the boundary conditions and see if your code will maintain a uniform set of conditions. Sometimes that can quickly point to a problem.

July 27, 2016, 10:53		#7
agd Senior Member Join Date: Jul 2009 Posts: 358 Rep Power: 19	The conservation equations as you have written them are not correct. Is that simply an error in how you wrote them in your post, or did you code them up that way?

July 27, 2016, 11:36		#9
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,897 Rep Power: 73	the RHS of your system must have the minus sign...

July 27, 2016, 11:45		#10
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,897 Rep Power: 73	for i=1:nx U1(i,j) = r(i,j); U2(i,j) = r(i,j)u(i,j); U3(i,j) = r(i,j).E(i,j); F1(i,j) = r(i,j)u(i,j); F2(i,j) = r(i,j)u(i,j).^2 + p(i,j); F3(i,j) = u(i,j)(r(i,j).E(i,j) + p(i,j)); end why do you mix the operators "" and "." in the for cycle? they are different operations in matlab

July 28, 2016, 11:45		#12
agd Senior Member Join Date: Jul 2009 Posts: 358 Rep Power: 19	Stupid question - are you sure that you are obeying the stability limit for the Lax-Friedrichs scheme?

July 28, 2016, 12:06		#15
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,897 Rep Power: 73	remember to estimate that dt/dx = 0.4 and the velocity should be considered as convective and sound velocity

January 20, 2020, 17:12		#17
silent_hunter New Member ΑΤΤΙΚΗΣ Join Date: Jan 2020 Posts: 1 Rep Power: 0	Hello!!Could you post the right code after your changes?

August 25, 2022, 05:11		#18
javeria New Member javeria Join Date: Aug 2022 Posts: 21 Rep Power: 4	hello i need lax friedrich matlab codes for 1D shock tube problem. Is there any one to help me?

August 25, 2022, 05:12		#19
javeria New Member javeria Join Date: Aug 2022 Posts: 21 Rep Power: 4	I need your codes

August 25, 2022, 05:16		#20
javeria New Member javeria Join Date: Aug 2022 Posts: 21 Rep Power: 4	hello selig can you post your codes after correction i need your codes