|
[Sponsors] |
August 29, 2008, 03:01 |
HLL Riemann Shock Tube Matlab Problem
|
#1 |
Guest
Posts: n/a
|
Hi There.
I am trying to write a simple code in MATLAB for an air-air shock tube using Godunov's method and the HLL Flux. I am using the ideal gas equation of state, and also assuming constant specific heats so we have P=RT*rho and e=(5/2)RT. My present code is having problems and I think it may be me misunderstanding the Godunov/HLL Flux. Please, could someone look at my code and help me out. See below. Thanks for the help clear all clc format long %----------------Constants and Initial Conditions-------------------------- Dx=0.001; %Length of the Cells Dt=1*10^(-7); %Timestep N=1000; %Number of cells = 2N T=300; %Number of Timesteps R=287.05; %Specific Gas Constant for air P_A=zeros(2*N,1); %Air Pressure of shock tube T_A=zeros(2*N,1); %Air Temperature of shock tube rho_A=zeros(2*N,1); %Air density of shock tube u_A=zeros(2*N,1); %Air velocity profile of shock tube c_A=zeros(2*N,1); %Air sound velocity profile of shock tube Ta=20; %Atmospheric Temperature in Celcius Pa=101325; %Atmospheric Pressure in Pascals rho_a=(Pa)/(R*(Ta+273.15));%density of air (Ideal Gas Law) u_a=0; %initial shock tube air velocity c_a=331.3*sqrt(1+Ta/273.15); %Air sound speed e_a=5/2*R*(Ta+273.15); %Air Internal Energy PD=50000; % Pressure difference in the tube %-----------Setting the initial condition in the tube---------------- for i=1:N P_A(i)=Pa; T_A(i)=Ta; rho_A(i)=rho_a; u_A(i)=u_a; c_A(i)=c_a; P_A(N+i)=Pa-PD; T_A(N+i)=Ta; rho_A(N+i)=P_A(N+i)/(R*(T_A(N+1)+273.15)); u_A(N+i)=u_a; c_A(N+i)=c_a; end e_A=5/2.*R.*(T_A+273.15); %-------------------------------------------------------------------------- %----------------Creating initial set of conserved variables--------------- for j=1:2*N Q_A(1,j)=rho_A(j); Q_A(2,j)=rho_A(j)*u_A(j); Q_A(3,j)=rho_A(j)*(e_A(j)+0.5*u_A(j)^2); end Q_OLD=Q_A; AIRPRESSURE=P_A; %--------------------HLL RIEMANN SOLVER---------------------------------- for timestep=1:T time=T-timestep Q_NEW(:,1)=Q_OLD(:,1); Q_NEW(:,2*N)=Q_OLD(:,2*N); for i=2:2*N-1 %--------------Calculating F(i+1/2)-------------------------- SL=min([u_A(i)-c_A(i),u_A(i+1)-c_A(i+1)]); SR=max([u_A(i)+c_A(i),u_A(i+1)+c_A(i+1)]); RHO_L=Q_A(1,i); U_L=Q_A(2,i)/RHO_L; E_L=Q_A(3,i)/RHO_L; P_L=AIRPRESSURE(i); RHO_R=Q_A(1,i+1); U_R=Q_A(2,i+1)/RHO_R; E_R=Q_A(3,i+1)/RHO_R; P_R=AIRPRESSURE(i+1); if SL>=0 F_A_PLUS=[RHO_L*U_L;RHO_L*U_L^2+P_L;U_L*RHO_L*E_L+U_L*P_L]; marker=1; else if SR<=0 F_A_PLUS=[RHO_R*U_R;RHO_R*U_R^2+P_R;U_R*RHO_R*E_R+U_R*P_R]; marker=2; else F_A_PLUS=(SR.*[RHO_L.*U_L;RHO_L.*U_L.^2+P_L;RHO_L.*E_L.*U_L+U_L.* P_L]-SL.*[RHO_R.*U_R;RHO_R.*U_R.^2+P_R;RHO_R.*E_R.*U_R+U_R.* P_R]+SL.*SR.*([RHO_R;RHO_R.*U_R;RHO_R.*E_R]-[RHO_L;RHO_L.*U_L;RHO_L.*E_L]))./(SR-SL); end end F(:,i)=F_A_PLUS; end %end cell iteration F(:,1)=F(:,2); F(:,2*N)=F(:,2*N-1); %------------------------------------------------------------------ %----------------Updating Vector of Conserved Variables--------------- for i=2:2*N-1 Q_NEW(:,i)=Q_OLD(:,i)+(Dt/Dx).*(F(:,i-1)-F(:,i)); end %---------------------------------------------------------------------- %------------------Updating Primitive variables----------------------- rho_A=Q_NEW(1,; u_A=Q_NEW(2,./rho_A; e_A=Q_NEW(3,./rho_A-0.5.*u_A.^2; T_A=(2/5).*e_A./R-273.15; c_a=331.3.*sqrt(1+Ta./273.15); AIRPRESSURE=(2/5).*rho_A.*e_A; %--------------------------------------------------------------------- clear Q_OLD Q_OLD=Q_NEW; clear Q_NEW end %end timestep |
|
July 7, 2009, 13:44 |
|
#2 |
Member
Nishant Kumar
Join Date: Jun 2009
Posts: 32
Rep Power: 17 |
Hi,
I am looking into kind of similar problem. Just wanted to know, were you able to crack out the problem? Thanks |
|
May 20, 2016, 02:10 |
|
#3 |
Member
LUQILIN
Join Date: May 2016
Location: Harbin, China
Posts: 35
Rep Power: 10 |
You should learn how to debug.
I suggest to give an linear initial condition to test your code. Then flux values at left and right of each interface should be same. Try to find out the errors. I go through the code quickly. Maybe something is wrong with SL ans SR: use abs(u)-c and abs(u)+c instead. |
|
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
shock tube!!! | showwa | FLUENT | 12 | May 20, 2012 19:14 |
1D Riemann problem for real gas | Nishu | Main CFD Forum | 0 | July 9, 2009 10:38 |
shock tube validation | AB | Siemens | 0 | November 21, 2004 18:43 |
About gas-water two-phase Riemann problem | Ma, Dong-Jun | Main CFD Forum | 4 | April 25, 2001 06:48 |
Shock tube source code | sencal | Main CFD Forum | 1 | December 10, 2000 03:27 |