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

Lax-wendroff scheme for Shock tube problem.

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

Reply
 
LinkBack Thread Tools Display Modes
Old   November 26, 2011, 10:54
Default Lax-wendroff scheme for Shock tube problem.
  #1
New Member
 
Manu Kamin
Join Date: Nov 2011
Posts: 16
Rep Power: 4
manukamin is on a distinguished road
Hey all

I'm stuck yet again. I'm once again not able to debug the matlab code for shock-tube problem using the lax-wendroff scheme. Can anyone please help me out ?

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;



x =-10:delx:10;


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




mom(1,1:c+1)=rho(1,1:c+1).*u(1,1:c+1);
e(1,1:c+1)=p(1,1:c+1)./((gamma-1)*rho(1,1:c+1));
et(1,1:c+1)=e(1,1:c+1)+0.5*u(1,1:c+1).^2;


Q(1,1:c+1)=rho(1,1:c+1);
Q(2,1:c+1)=rho(1,1:c+1).*u(1,1:c+1);
Q(3,1:c+1)=rho(1,1:c+1).*et(1,1:c+1);

F(1,1:c+1)=Q(1,1:c+1).*u(1,1:c+1);
F(2,1:c+1)=Q(2,1:c+1).*u(1,1:c+1)+ p(1,1:c+1);
F(3,1:c+1)=Q(3,1:c+1).*u(1,1:c+1)+ p(1,1:c+1).*u(1,1:c+1);


for n=2:1:3

for i=2:1:c

uavplus=0.5*(u(n-1,i+1)+ u(n-1,i));
uavminus=0.5*(u(n-1,i-1)+ u(n-1,i));

etavplus=0.5*(et(n-1,i+1)+ et(n-1,i));
etavminus=0.5*(et(n-1,i-1)+ et(n-1,i));

Aplus(1,1)=0;
Aplus(1,2)=1;
Aplus(1,3)=0;
Aplus(2,1)=.5*(gamma-3)*(uavplus)^2;
Aplus(2,2)=(3-gamma)*(uavplus);
Aplus(2,3)=gamma-1;
Aplus(3,1)=-gamma*uavplus*etavplus + (gamma-1)*(uavplus)^3;
Aplus(3,2)=gamma*etavplus - 3/2*(gamma-1)*(uavplus)^2;
Aplus(3,3)=gamma*uavplus;


Aminus(1,1)=0;
Aminus(1,2)=1;
Aminus(1,3)=0;
Aminus(2,1)=.5*(gamma-3)*(uavminus)^2;
Aminus(2,2)=(3-gamma)*(uavminus);
Aminus(2,3)=gamma-1;
Aminus(3,1)=-gamma*uavminus*etavminus + (gamma-1)*(uavminus)^3;
Aminus(3,2)=gamma*etavminus - 3/2*(gamma-1)*(uavminus)^2;
Aminus(3,3)=gamma*uavminus;


Qn(1:3,i)=Q(1:3,i)-0.5*lambda*(F(1:3,i+1)-F(1:3,i-1))+ .5*lambda^2*(Aplus*(F(1:3,i+1)-F(1:3,i)) - Aminus*(F(1:3,i)-F(1:3,i-1)));

rho(n,i)=Qn(1,i);
mom(n,i)=Qn(2,i);
et(n,i)=Qn(3,i)/Qn(1,i);
u(n,i)=mom(n,i)/rho(n,i);
p(n,i)=(gamma-1)*rho(n,i)*(et(n,i)-0.5*u(n,i)^2);
e(n,i)=et(n,i)-0.5*u(n,i)^2;

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);
et(n,1)=et(1,1);
et(n,c+1)=et(1,c+1);
u(n,1)=u(1,1);
u(n,c+1)=u(1,c+1);
p(n,1)=p(1,1);
p(n,c+1)=p(1,c+1);
e(n,1)=e(1,1);
e(n,c+1)=e(1,c+1);

Q(1:3,2:c)=Qn(1:3,2:c);
%Qn(1:3,:)=0;

F(1,1:c+1)=Q(1,1:c+1).*u(n,1:c+1);
F(2,1:c+1)=Q(2,1:c+1).*u(n,1:c+1)+ p(n,1:c+1);
F(3,1:c+1)=Q(3,1:c+1).*u(n,1:c+1)+ p(n,1:c+1).*u(n,1:c+1);

Qn(1:3,:)=0;

end





Manu

Last edited by manukamin; November 26, 2011 at 11:40.
manukamin is offline   Reply With Quote

Old   November 26, 2011, 14:10
Default
  #2
New Member
 
Manu Kamin
Join Date: Nov 2011
Posts: 16
Rep Power: 4
manukamin is on a distinguished road
I think there is no problem with my code. In the sense there are no bugs. My doubt is that my formulation of the equation itself might be wrong. Could anyone please tell me if my formulation is correct?

My vector equation is:

U(n,i)=U(n-1,i) -lambda/2*(F(n-1,i+1)-F(b-1,i-1)) +lambda^2/2*(A i+1/2*(F(i+1)-F(i) - A i-1/2*(F(i)-F(i-1)));

where U and F are primary and flux vectors. and A is the matrix dF/dU.
manukamin 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
lax friedrich scheme for shock tube problem. manukamin Main CFD Forum 2 November 3, 2013 02:21
problem with second order backward Euler scheme shahpar73 CFX 1 September 28, 2010 18:54
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 14:36


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