# Lax-wendroff scheme for Shock tube problem.

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

 November 26, 2011, 11:54 Lax-wendroff scheme for Shock tube problem. #1 New Member   Manu Kamin Join Date: Nov 2011 Posts: 16 Rep Power: 7 Sponsored Links 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 12:40.

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

 Thread Tools Display Modes Linear 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 On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post manukamin Main CFD Forum 3 March 22, 2016 02:02 shahpar73 CFX 1 September 28, 2010 18:54 libin Main CFD Forum 4 July 1, 2004 04:32 ger Main CFD Forum 1 September 24, 2003 07:41 AERO Main CFD Forum 2 January 19, 2000 15:36