CFD Online Discussion Forums - Help regarding 1D compressible flow (Sod Shock tube)

Page 1 of 3

Show 40 post(s) from this thread on one page

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- Main CFD Forum (https://www.cfd-online.com/Forums/main/)

- - Help regarding 1D compressible flow (Sod Shock tube) (https://www.cfd-online.com/Forums/main/175183-help-regarding-1d-compressible-flow-sod-shock-tube.html)

selig5576

July 25, 2016 16:22

Help regarding 1D compressible flow (Sod Shock tube)

1 Attachment(s)

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 10:40

Still having issues

1 Attachment(s)

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 11: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 11:39

Thanks

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 13: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 13:49

Further debugging

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 09: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 10:27

Equations

1 Attachment(s)

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 10:36

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

FMDenaro

July 27, 2016 10: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 10: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 10:45

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

selig5576

July 28, 2016 10:58

CFL condition

Hi agd,

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

FMDenaro

July 28, 2016 11:04

Quote:

Originally Posted by agd (Post 611793)

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 11:06

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

selig5576

August 2, 2016 09:22

Thanks

1 Attachment(s)

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 16:12

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

javeria

August 25, 2022 04:11

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

javeria

August 25, 2022 04:12

I need your codes

javeria

August 25, 2022 04:16

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

All times are GMT -4. The time now is 13:03.

Page 1 of 3

Show 40 post(s) from this thread on one page