# HLL Riemann Shock Tube Matlab Problem

 Register Blogs Members List Search Today's Posts Mark Forums Read August 29, 2008, 04:01 HLL Riemann Shock Tube Matlab Problem #1 Luke F 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, 14:44 #2 Member   Nishant Kumar Join Date: Jun 2009 Posts: 32 Rep Power: 14 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, 03:10 #3 Member   LUQILIN Join Date: May 2016 Location: Harbin, China Posts: 35 Rep Power: 7 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 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 showwa FLUENT 12 May 20, 2012 20:14 Nishu Main CFD Forum 0 July 9, 2009 11:38 AB Siemens 0 November 21, 2004 19:43 Ma, Dong-Jun Main CFD Forum 4 April 25, 2001 07:48 sencal Main CFD Forum 1 December 10, 2000 04:27

All times are GMT -4. The time now is 13:34.