|
[Sponsors] |
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. |
|
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
lax friedrich scheme for shock tube problem. | manukamin | Main CFD Forum | 3 | March 22, 2016 01:02 |
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 |