# 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: 17 Rep Power: 13 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: 17 Rep Power: 13 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 Search this Thread Show Printable Version Email this Page Search this Thread: Advanced Search Display Modes Linear Mode Switch to Hybrid Mode Switch to Threaded 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 Off 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 19:54 libin Main CFD Forum 4 July 1, 2004 05:32 ger Main CFD Forum 1 September 24, 2003 08:41 AERO Main CFD Forum 2 January 19, 2000 15:36

All times are GMT -4. The time now is 11:07.