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

lax friedrich scheme for shock tube problem.

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

Reply
 
LinkBack Thread Tools Display Modes
Old   November 24, 2011, 15:09
Default lax friedrich scheme for shock tube problem.
  #1
New Member
 
Manu Kamin
Join Date: Nov 2011
Posts: 16
Rep Power: 5
manukamin is on a distinguished road
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.

Can someone please help me debug it?

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
File Type: txt laxfriedrich_shocktube.txt (1.9 KB, 13 views)
manukamin is offline   Reply With Quote

Old   November 25, 2011, 01:15
Default
  #2
New Member
 
Manu Kamin
Join Date: Nov 2011
Posts: 16
Rep Power: 5
manukamin is on a distinguished road
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.
manukamin is offline   Reply With Quote

Old   November 3, 2013, 03:21
Default
  #3
New Member
 
Alpesh
Join Date: Oct 2013
Posts: 1
Rep Power: 0
alpesh_iitb is on a distinguished road
Hi,
can you tell me what is lambda here?
also can you post the updated code?
alpesh_iitb is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Problem simulating wave propagation in a tube cKnoop OpenFOAM Running, Solving & CFD 6 August 22, 2011 15:28
Problem with using icoFoam for flow in a straight tube edificelp OpenFOAM 4 February 11, 2011 04:27
problem about numerical scheme in LES. libin Main CFD Forum 4 July 1, 2004 04:32
Shock tube problem ger Main CFD Forum 1 September 24, 2003 07:41
Difference between Lax and MacCormack methods on solutions of 1D problem. AERO Main CFD Forum 2 January 19, 2000 15:36


All times are GMT -4. The time now is 09:21.