CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

HLL Riemann Shock Tube Matlab Problem

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   August 29, 2008, 03:01
Default 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

  Reply With Quote

Old   July 7, 2009, 13:44
Default
  #2
Member
 
Nishant Kumar
Join Date: Jun 2009
Posts: 32
Rep Power: 17
Nishu is on a distinguished road
Hi,

I am looking into kind of similar problem.
Just wanted to know, were you able to crack out the problem?

Thanks
Nishu is offline   Reply With Quote

Old   May 20, 2016, 02:10
Default
  #3
Member
 
LUQILIN
Join Date: May 2016
Location: Harbin, China
Posts: 35
Rep Power: 10
LUQILIN is on a distinguished road
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.
LUQILIN is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

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 Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


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


All times are GMT -4. The time now is 10:58.