# lax friedrich scheme for shock tube problem.

 Register Blogs Members List Search Today's Posts Mark Forums Read November 24, 2011, 15:09 lax friedrich scheme for shock tube problem.
#1
New Member

Manu Kamin
Join Date: Nov 2011
Posts: 17
Rep Power: 13 Hi

I have written a matlab code for shocktube problem using the lax friedrich scheme. However there seems to be a bug in it, due to which it is not giving me the expected results for plots of density, pressure and velocity.

Here's the code:

clear all;
close all;
clc;

R=287;
%Cv=717;
gamma=1.4;
rhol=1;
Ul=0;
Pl=100000;

rhor=0.125;
Ur=0;
Pr=10000;

c=50;
delx=20/c;
delt=.0004278;
lambda=.001069;

for i=1:1:c+1
x(i)=-10+(i-1)*delx;
end

for i=1:1:c+1

if (x(i)<0)

rho(1,i)=rhol;
u(1,i)=Ul;
p(1,i)=Pl;

else
rho(1,i)=rhor;
u(1,i)=Ur;
p(1,i)=Pr;
end
end

for i=1:1:c+1

mom(1,i)=rho(1,i)*u(1,i);
%t(1,i)=p(1,i)/(rho(1,i)*R);
e(1,i)=p(1,i)/((gamma-1)*rho(1,i));
energ(1,i)=(e(1,i)+.5*u(1,i)^2);
energrho(1,i)=energ(1,i)*rho(1,i);
end

for n=2:1:30

for i=2:1:c

rho(n,i)= .5*(rho(n-1,i-1)+ rho(n-1,i+1))-.5*lambda*(u(n-1,i+1)*rho(n-1,i+1)-u(n-1,i-1)*rho(n-1,i-1));
mom(n,i)= .5*(mom(n-1,i-1)+ mom(n-1,i+1))-.5*lambda*(u(n-1,i+1)*mom(n-1,i+1)+p(n-1,i+1)-p(n-1,i-1)-u(n-1,i-1)*mom(n-1,i-1));
energrho(n,i)= .5*(energrho(n-1,i-1)+energrho(n-1,i+1))-.5*lambda*((u(n-1,i+1)*energrho(n-1,i+1)+p(n-1,i+1)*u(n-1,i+1) - u(n-1,i-1)*energrho(n-1,i-1)+p(n-1,i-1)*u(n-1,i-1)));

energ(n,i)=energrho(n,i)/rho(n,i);
u(n,i)=mom(n,i)/rho(n,i);
e(n,i)=energ(n,i)-.5*u(n,i)^2;
%t(n,i)=e(n,i)/Cv;
p(n,i)=(gamma-1)*rho(n,i)*e(n,i);

end

rho(n,1)=rho(1,1);
rho(n,c+1)=rho(1,c+1);
mom(n,1)=mom(1,1);
mom(n,c+1)=mom(1,c+1);
e(n,1)=e(1,1);
e(n,c+1)=e(1,c+1);
%t(n,1)=t(n,2);
%t(n,1)=e(n,1)/Cv;
%t(n,c+1)=t(n,c);
%t(n,c+1)=e(n,c+1)/Cv;
p(n,1)=p(1,1);
%p(n,1)=rho(n,1)*R*t(n,1);
p(n,c+1)=p(1,c+1);
%p(n,c+1)=rho(n,c+1)*R*t(n,c+1);
u(n,1)=u(1,1);
%u(n,1)=mom(n,1)/rho(n,1);
u(n,c+1)=u(1,c+1);
%u(n,c+1)=mom(n,c+1)/rho(n,c+1);

end

%t=plot(x,p(30,: ),'or');
%axis([-10 10 0 100000])

%t=plot(x,u(30,: ),'or');
%axis([-10 10 -150 450])

t=plot(x,rho(25,: ),'or');
axis([-10 10 0 1.1])

Thanks
Manu
Attached Files laxfriedrich_shocktube.txt (1.9 KB, 34 views)   November 25, 2011, 01:15 #2 New Member   Manu Kamin Join Date: Nov 2011 Posts: 17 Rep Power: 13 Hi never mind.. I fixed it.. I was not updating the value of energy at boundary points and it was hence taking default value of zero which is not correct.   November 3, 2013, 03:21 #3 New Member   Alpesh Join Date: Oct 2013 Posts: 1 Rep Power: 0 Hi, can you tell me what is lambda here? also can you post the updated code?   March 22, 2016, 02:02 #4 New Member   Yi Han Join Date: Oct 2013 Location: Laramie WY Posts: 15 Rep Power: 11 Hi, can you upload the updated code? Thanks so much!  Thread Tools Search this Thread Show Printable Version Email this Page Search this Thread: Advanced Search Display Modes Linear Mode Switch to Hybrid Mode Switch to Threaded Mode 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 OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules Similar Threads Thread Thread Starter Forum Replies Last Post cKnoop OpenFOAM Running, Solving & CFD 6 August 22, 2011 16:28 edificelp OpenFOAM 4 February 11, 2011 04:27 libin Main CFD Forum 4 July 1, 2004 05:32 ger Main CFD Forum 1 September 24, 2003 08:41 AERO Main CFD Forum 2 January 19, 2000 15:36

All times are GMT -4. The time now is 06:16.