% program principal
clear all
clc
format long
gama=1.4;
L=1;
nr_cel=100;
dx=L/nr_cel;
tf=0.0005;
CFL=1;
l_ref=dx/2;
% prealocarea memoriei
U1new=zeros(1,nr_cel);
U2new=zeros(1,nr_cel);
U3new=zeros(1,nr_cel);
ro=zeros(1,nr_cel);
u=zeros(1,nr_cel);
p=zeros(1,nr_cel);
Mach=zeros(1,nr_cel);
% definirea conditiilor initiale
[U1,U2,U3] = conditii_initiale(gama,nr_cel);
% afisarea contitiilor initiale
for i=1:nr_cel
[ ro(i),u(i),p(i)] = cons_to_prim( gama,U1(i),U2(i),U3(i) );
c = vit_sunet( gama,p(i),ro(i) );
Mach(i) = u(i)/c;
end
x = dx/2:dx:L-dx/2;
subplot(2,2,1), h1=plot(x,ro,'YDataSource','ro');
title ('densitate');
xlabel('x')
ylim([0 1.3])
subplot(2,2,2), h2=plot(x,u,'YDataSource','u'); title ('viteza'); xlabel('x')
ylim([-50 350])
subplot(2,2,3), h3=plot(x,p,'YDataSource','p'); title('presiune');
xlabel('x')
ylim([0 1.2e+5])
subplot(2,2,4), h4=plot(x,Mach,'YDataSource','Mach'); title('Mach');
xlabel('x')
ylim([-0.2 1.2])
t=0;
while t<tf
% pasul de timp
dt = pas_timp( l_ref,gama,nr_cel,U1,U2,U3 );
% definirea conditiilor la limita (perete solid la capetele tubului)
[ Us1,Us2,Us3,Ud1,Ud2,Ud3 ] = conditii_limita( nr_cel,gama,U1,U2,U3 );
% calculul fluxurilor
[ F1,F2,F3 ] = fluxuri( nr_cel,gama,U1,U2,U3,Us1,Us2,Us3,Ud1,Ud2,Ud3 );
% calculul variabilelor conservative la momentul actual
for i=1:nr_cel
U1new(i) = U1(i) - CFL*dt/dx*( F1(i+1) - F1(i) );
U2new(i) = U2(i) - CFL*dt/dx*( F2(i+1) - F2(i) );
U3new(i) = U3(i) - CFL*dt/dx*( F3(i+1) - F3(i) );
end
% calculul variabilelor primitive si a nr-ului Mach
for i=1:nr_cel
[ ro(i),u(i),p(i)] = cons_to_prim( gama,U1new(i),U2new(i),U3new(i) );
c = vit_sunet( gama,p(i),ro(i) );
Mach(i) = u(i)/c;
end
% afisarea rezultatelor
refreshdata(h1,
drawnow
refreshdata(h2,
drawnow
refreshdata(h3,
drawnow
refreshdata(h4,
drawnow
'caller')
'caller')
'caller')
'caller')
% update
U1 = U1new;
U2 = U2new;
U3 = U3new;
t = t + dt;
display(t)
end
function [Ui1,Ui2,Ui3] = conditii_initiale(gama,nr_cel)
format long
% functia defineste conditiile initiale ale problemei tubului de soc
ro4=1;
u4=0;
p4=10^5;
ro1=0.125;
u1=0;
p1=10^4;
% prealocare memorie
Ui1=zeros(1,nr_cel);
Ui2=zeros(1,nr_cel);
Ui3=zeros(1,nr_cel);
for i=1:nr_cel/2
Ui1(i) = ro4;
Ui2(i) = ro4*u4;
Ui3(i) = ro4 * ( 1/(gama-1) * p4/ro4 + u4^2/2);
end
for i=(nr_cel/2 + 1):nr_cel
Ui1(i) = ro1;
Ui2(i) = ro1*u1;
Ui3(i) = ro1 * ( 1/(gama-1) * p1/ro1 + u1^2/2);
end
end
function [ Us1,Us2,Us3,Ud1,Ud2,Ud3 ] = conditii_limita( nr_cel, gama ,U1, U2,
U3 )
format long
% functia defineste conditiile la limita ale problemei tubului de soc
% (perete solid la capetele tubului)
[ ro1,u1,p1 ] = cons_to_prim( gama, U1(1),U2(1),U3(1) );
[ ron,un,pn ] = cons_to_prim( gama, U1(nr_cel),U2(nr_cel),U3(nr_cel) );
ros = ro1;
us = -u1;
ps = p1;
rod = ron;
ud = -un;
pd = pn;
Us1 = ros;
Us2 = ros*us;
Us3 = ros * ( 1/(gama-1) * ps/ros + us^2/2);
Ud1 = rod;
Ud2 = rod*ud;
Ud3 = rod * ( 1/(gama-1) * pd/rod + ud^2/2);
end
function [ Us1,Us2,Us3,Ud1,Ud2,Ud3 ] = conditii_limita( nr_cel, gama ,U1, U2,
U3 )
format long
% functia defineste conditiile la limita ale problemei tubului de soc
% (perete solid la capetele tubului)
[ ro1,u1,p1 ] = cons_to_prim( gama, U1(1),U2(1),U3(1) );
[ ron,un,pn ] = cons_to_prim( gama, U1(nr_cel),U2(nr_cel),U3(nr_cel) );
ros = ro1;
us = -u1;
ps = p1;
rod = ron;
ud = -un;
pd = pn;
Us1 = ros;
Us2 = ros*us;
Us3 = ros * ( 1/(gama-1) * ps/ros + us^2/2);
Ud1 = rod;
Ud2 = rod*ud;
Ud3 = rod * ( 1/(gama-1) * pd/rod + ud^2/2);
end
function [ F1,F2,F3 ]= fluxuri( nr_cel,gama,U1,U2,U3,Us1,Us2,Us3,Ud1,Ud2,Ud3 )
format long
% functia calculeaza fluxurile pe interfetele ce marginesc celulele
% prealocarea memoriei
F1=zeros(1,nr_cel+1);
F2=zeros(1,nr_cel+1);
F3=zeros(1,nr_cel+1);
for i=1:nr_cel+1
% definirea starilor L si R
if i==1
UR(1) = U1(i);
UR(2) = U2(i);
UR(3) = U3(i);
UL(1) = Us1;
UL(2) = Us2;
UL(3) = Us3;
elseif i==nr_cel+1
UL(1) = U1(nr_cel);
UL(2) = U2(nr_cel);
UL(3) = U3(nr_cel);
UR(1) = Ud1;
UR(2) = Ud2;
UR(3) = Ud3;
else
UL(1) = U1(i-1);
UL(2) = U2(i-1);
UL(3) = U3(i-1);
UR(1) = U1(i);
UR(2) = U2(i);
UR(3) = U3(i);
end
% calculul fluxurilor la interfete
flux = schema_roe(gama,UL,UR);
F1(i) = flux(1);
F2(i) = flux(2);
F3(i) = flux(3);
end
end
function dt = pas_timp( l_ref,gama,nr_cel,U1,U2,U3 )
format long
% functia calculeaza pasul de timp
% prealocare memorie
lambda=zeros(1,nr_cel);
for i=1:nr_cel
[ ro,u,p ] = cons_to_prim( gama,U1(i),U2(i),U3(i) );
c = vit_sunet( gama,p,ro );
lambda(i) = u + c;
end
v = max(lambda);
dt = l_ref/v;
end
function p=presiune(gama,U1,U2,U3)
format long
% functia calculeaza presiunea in functie de valorile vectorului
% variabilelor conservative;
p =(gama-1)*(U3-0.5*U2^2/U1);
end
function [flux_num]=schema_roe(gama,UL,UR)
% variabile primitive
% starea l
[ rol,ul,pl ] = cons_to_prim( gama,UL(1),UL(2),UL(3) );
cl = vit_sunet( gama,pl,rol );
hl = cl*cl/(gama-1) + 0.5*ul*ul;
% starea r
[ ror,ur,pr ] = cons_to_prim( gama,UR(1),UR(2),UR(3) );
cr = vit_sunet( gama,pr,ror );
hr = cr*cr/(gama-1) + 0.5*ur*ur;
% valorile medii roe
rt = sqrt(ror/rol);
ro = rt*rol;
u = (ul + rt*ur)/(1 + rt);
h = (hl + rt*hr)/(1 + rt);
c = sqrt( (gama - 1)*( h - 0.5*u*u ) );
% amplitudinile undelor
dro = ror - rol;
du = ur - ul;
dp = pr - pl;
ldu = [ 0.5*(dp-ro*c*du)/(c^2) ,
dro - dp/(c*c) , 0.5*(dp+ro*c*du)/(c^2) ];
% valorile absolute ale vitezelor undelor
ws = [ abs(u-c) , abs(u) , abs(u+c) ];
% corectia de entropie (Harten); doar pentru campurile neliniare
eps = 0.001;
if ( ws(2) < eps )
ws(2) = 0.5 * ( ws(2)*ws(2)/eps + eps );
end
if ( ws(3) < eps )
ws(3) = 0.5 * ( ws(3)*ws(3)/eps + eps );
end
% vectori proprii la dreapta
r = [ 1 , 1 , 1 ;...
u - c , u , u + c;...
h - u * c , u^2/2 , h + u * c];
% termenul disipativ: |A|*(UR-UL)
diss = zeros(1,3);
for i=1:3
for j=1:3
diss(i) = diss(i) + ws(j)*ldu(j)*r(i,j);
end
end
% fluxul fizic
fl = [rol * ul , rol * ul * ul + pl , rol * ul * hl];
fr = [ror * ur , ror * ur * ur + pr , ror * ur * hr];
% fluxul numeric
flux_num = zeros(1,3);
for i=1:3
flux_num(i) = 0.5 * (fl(i) + fr(i) - diss(i));
end
end
function c = vit_sunet( gama,p,ro )
format long
% functia calculeaza viteza sunetului
c=sqrt(gama*p/ro);
end