1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
| %boucle temporelle
for i=2:N
temp0=Tp(i-1,j);
temp1=0;
temp2=Tp(i-1,j);
rayon0=r(i-1,j);
test=0;
rayon=r(i-1,j);
%boucle de résolution de l'équation différentielle
while (abs(temp2-temp1)>prec)&(test==0)
%%%%calcul des caractéristique non constantes%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if rayon<0
rayon=0;
temp2=Tg;
test=1;
else
%flux de masse [kg/m².s]
mu=((M_v/(2*pi*R_GP*temp2))^0.5)*Pref*exp(delta_H*(temp2-Tref)/(R_GP*temp2*Tref));
rayon=rayon0-mu/rho_s;
%capacacité calorifique des suies [J/kg.K]
Cp_p=(16.86+0.00477*temp2-85400/(temp2*temp2))/0.012;
%emmisivité
eps=48*pi*rayon*n_r*n_i/(lambda*(4*n_r*n_r*n_i*n_i+(2+n_r*n_r-n_i*n_i)*(2+n_r*n_r-n_i*n_i)));
%capacité calorifique gaz
somme2=0;
dTemp=(temp2-Tg)/Ni;
for ka=1:Ni
tampon=28.58+0.00377*temp2-(5e4)/(temp2*temp2);
tampon=tampon-R_GP;
somme2=somme2+tampon*dTemp;
end
integrale=somme2/M_v;
%coefficient de transfert radiatif
h_r=sigma*(Te*Te+temp2*temp2)*(Te+temp2);
%%%%calcul du temps d'integration suivant%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
temp1=temp2;
tampon=dt/(3*rayon*rho_s*Cp_p);
temp2=(temp0+tampon*(eps*(q(i,j)/4+h_r*Te)-mu*(delta_H/M_s+alpha*(-R_GP*Tg/(2*M_v)+integrale))))/(1+tampon*(eps*h_r+alpha*mu*R_GP/(2*M_v)));
end
end
Tp(i,j)=temp2;
r(i,j)=rayon; |
Partager