CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Help regarding 1D compressible flow (Sod Shock tube)

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   July 25, 2016, 17:22
Default Help regarding 1D compressible flow (Sod Shock tube)
  #1
Senior Member
 
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10
selig5576 is on a distinguished road
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
Attached Files
File Type: txt LF1DEuler.txt (1.7 KB, 10 views)

Last edited by selig5576; July 26, 2016 at 12:44.
selig5576 is offline   Reply With Quote

Old   July 26, 2016, 11:40
Default Still having issues
  #2
Senior Member
 
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10
selig5576 is on a distinguished road
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.
Attached Images
File Type: png PostProcess.png (61.5 KB, 27 views)
selig5576 is offline   Reply With Quote

Old   July 26, 2016, 12:29
Default
  #3
agd
Senior Member
 
Join Date: Jul 2009
Posts: 349
Rep Power: 18
agd is on a distinguished road
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.
agd is offline   Reply With Quote

Old   July 26, 2016, 12:39
Default Thanks
  #4
Senior Member
 
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10
selig5576 is on a distinguished road
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
selig5576 is offline   Reply With Quote

Old   July 26, 2016, 14:00
Default
  #5
agd
Senior Member
 
Join Date: Jul 2009
Posts: 349
Rep Power: 18
agd is on a distinguished road
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.
agd is offline   Reply With Quote

Old   July 26, 2016, 14:49
Default Further debugging
  #6
Senior Member
 
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10
selig5576 is on a distinguished road
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)
selig5576 is offline   Reply With Quote

Old   July 27, 2016, 10:53
Default
  #7
agd
Senior Member
 
Join Date: Jul 2009
Posts: 349
Rep Power: 18
agd is on a distinguished road
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?
agd is offline   Reply With Quote

Old   July 27, 2016, 11:27
Default Equations
  #8
Senior Member
 
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10
selig5576 is on a distinguished road
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.
Attached Images
File Type: png PostProcess.png (36.1 KB, 21 views)
selig5576 is offline   Reply With Quote

Old   July 27, 2016, 11:36
Default
  #9
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,760
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
the RHS of your system must have the minus sign...
FMDenaro is offline   Reply With Quote

Old   July 27, 2016, 11:45
Default
  #10
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,760
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
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
FMDenaro is offline   Reply With Quote

Old   July 27, 2016, 11:48
Default
  #11
Senior Member
 
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10
selig5576 is on a distinguished road
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.
selig5576 is offline   Reply With Quote

Old   July 28, 2016, 11:45
Default
  #12
agd
Senior Member
 
Join Date: Jul 2009
Posts: 349
Rep Power: 18
agd is on a distinguished road
Stupid question - are you sure that you are obeying the stability limit for the Lax-Friedrichs scheme?
agd is offline   Reply With Quote

Old   July 28, 2016, 11:58
Default CFL condition
  #13
Senior Member
 
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10
selig5576 is on a distinguished road
Hi agd,

Yes I am. I've made sure there is no instability aside the downfalls of the method.
selig5576 is offline   Reply With Quote

Old   July 28, 2016, 12:04
Default
  #14
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,760
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
Quote:
Originally Posted by agd View Post
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 is offline   Reply With Quote

Old   July 28, 2016, 12:06
Default
  #15
Senior Member
 
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,760
Rep Power: 71
FMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura aboutFMDenaro has a spectacular aura about
remember to estimate that dt/dx = 0.4 and the velocity should be considered as convective and sound velocity
FMDenaro is offline   Reply With Quote

Old   August 2, 2016, 10:22
Default Thanks
  #16
Senior Member
 
Selig
Join Date: Jul 2016
Posts: 213
Rep Power: 10
selig5576 is on a distinguished road
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
Attached Images
File Type: png Shock.png (39.7 KB, 32 views)
selig5576 is offline   Reply With Quote

Old   January 20, 2020, 17:12
Default
  #17
New Member
 
ΑΤΤΙΚΗΣ
Join Date: Jan 2020
Posts: 1
Rep Power: 0
silent_hunter is on a distinguished road
Hello!!Could you post the right code after your changes?
silent_hunter is offline   Reply With Quote

Old   August 25, 2022, 05:11
Default
  #18
New Member
 
javeria
Join Date: Aug 2022
Posts: 21
Rep Power: 3
javeria is on a distinguished road
hello i need lax friedrich matlab codes for 1D shock tube problem. Is there any one to help me?
javeria is offline   Reply With Quote

Old   August 25, 2022, 05:12
Default
  #19
New Member
 
javeria
Join Date: Aug 2022
Posts: 21
Rep Power: 3
javeria is on a distinguished road
I need your codes
javeria is offline   Reply With Quote

Old   August 25, 2022, 05:16
Default
  #20
New Member
 
javeria
Join Date: Aug 2022
Posts: 21
Rep Power: 3
javeria is on a distinguished road
hello selig can you post your codes after correction i need your codes
javeria is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


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


All times are GMT -4. The time now is 05:12.