HLL Riemann Shock Tube Matlab Problem
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 |
Hi,
I am looking into kind of similar problem. Just wanted to know, were you able to crack out the problem? Thanks |
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. |
All times are GMT -4. The time now is 15:56. |