Lax-wendroff scheme for Shock tube problem.
 

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

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.