Lax-wendroff scheme for Shock tube problem.

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

 November 26, 2011, 10:54 Lax-wendroff scheme for Shock tube problem. #1 New Member   Manu Kamin Join Date: Nov 2011 Posts: 17 Rep Power: 14 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.

 November 26, 2011, 14:10 #2 New Member   Manu Kamin Join Date: Nov 2011 Posts: 17 Rep Power: 14 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.