%Konzentrations-, Geschwindigkeits- und Temperaturverteilung
%
%m-file zur Bestimmung von Konzentrations-, Geschwindigkeits- und
%Temperaturverteilung in Freistrahlen
%
%Autor: Carsten von der Heyde
%
%Variablen
%
clc
clear all
close all                       %Alle (offenen) Fenster schließen
syms y z                        %Deklaration Variablen
%Brennerleistung in MW
P=10;
%Gasdruck in mbarü  
p0=600;
%Gastemperatur in K
T0=293;
%Duesenbohrung in mm
d0=10.4;
%Freistrahlwinkel in °
Gamma=19;
%Primaergasanteil
Pri=0.41;
%Sekundaergasanteil
Sek=0.59;
%Lambda
lam=1.08;
%Luftfeuchte
phi=0.3;
%Berechnungspunkt='x-Koordinate in mm eingeben'
x=300;
%
%
%------------------------------
%Feste Parameter
%------------------------------
%Gaskonstanten [kJ/(kg K)]
R_L=287;
R_co2=188.9;
R_o2=259.8;
R_n2=296.8;
R_h20=462;
R_brst=449;
Rm=8.314472;                    %allg. Gaskonstante [J/(mol*K)]
TN=273.15;                      %Normtemperatur [K]
pN=1013.25;                     %Normdruck [mbar]
priGD=8;                        %Anzahl Primaergasdüsen
sekGD=8;                        %Anzahl Sekundärgasdüsen
k0=1;                           %Anfangskonzentration Brst
kf=1;                           %Korrekturfaktor für Düsenaustrittsgeometrie
ci=0.075;                       %Übertragungsfaktor Impuls
cq=0.087;                       %Übertragungsfaktor Wärme
cc=cq;                          %Übertragungsfaktor Konzentration
TL=293;                         %Verbrennugslufttemperatur[K]
ps=0.023392;                    %Wasserdampfpartialdruck [mbar]
pL=1013.25;                     %Luftdruck [mbar]
TRG=1473;                       %Rausgastemperatur [K]
pFR=pL;                         %Feuerraumdruck
cp0=1.7;                        %Wärmekapazität, isobar, Brst
cpL=1;                          %Wärmekapazität, isobar, Luft
%Molare isobare Wärmekapazität 0-1900°C [kJ/(kmol*K)]
cp_co2=54.1;                        
cp_h2o=43.5;
cp_n2=33.1;
cp_o2=35;
%Molmasse [kg/kmol]
M_co2=44.01;
M_h2o=18.01528;
M_n2=28.013;
M_o2=31.9988;
%Erdgas L - Zusammensetzung
ch4=0.82;
c2h6=0.033;
c3h8=0.006;
c4h10=0.003;
n2=0.126;
co2=0.012;
%-----------------------------------------------
%Mindestsauerstoffbedarf o_min
%-----------------------------------------------
o_min=(1+4/4)*ch4+(2+6/4)*c2h6+(3+8/4)*c3h8+(4+10/4)*c4h10
%-----------------------------------------------
%Mindestluftbedarf l_min
%-----------------------------------------------
l_min=o_min/0.21
%-----------------------------------------------
%Tatsaechlicher Luftbedarf l_tat und l_tatf
%-----------------------------------------------
l_tat=lam*l_min;
wL=phi*ps/(pL-phi*ps);
l_tatf=(wL+1)*l_min;
%-----------------------------------------------
%Norm- und Betriebsvolumenstrom
%-----------------------------------------------
VN_brst=P/31.8                          %Normvolumenstrom Brst [m³/s]
VB_brst=VN_brst*pN*T0/((p0+pN)*TN)      %Betreibsvolumenstrom Brst[m³/s]
VN_L=l_tatf*VN_brst                     %Normvolumenstrom Luft [m³/s]
VB_L=pN*VN_L*TL/(pL*TN)                 %Betreibsvolumenstrom Brst[m³/s]
%-----------------------------------------------
%Abgaszusammensetzung
%-----------------------------------------------
RG_co2=10*(ch4+c2h6+c3h8+c4h10)+co2
RG_h2o=28/2*(ch4+c2h6+c3h8+c4h10)+lam*wL*l_min
RG_n2=0.79*lam*l_min
RG_o2=0.21*(lam-1)+l_min
Rauchgasvolumen_spez=RG_co2+RG_h2o+RG_n2+RG_o2
VR_spez=Rauchgasvolumen_spez;
Rauchgasvolumenstrom=VN_brst*VR_spez
VR_N=Rauchgasvolumenstrom;
vR_N_tr=RG_co2+RG_n2+RG_o2;
VR_N_tr=vR_N_tr*VN_brst;
%-----------------------------------------------
%Raumanteile der Rauchgasbestandteile im VR
%-----------------------------------------------
r_co2=RG_co2*VN_brst/VR_N;
r_n2=RG_n2*VN_brst/VR_N;
r_o2=RG_o2*VN_brst/VR_N;
r_h2o=RG_h2o*VN_brst/VR_N;
%-----------------------------------------------
%Berechnung der Abgasdichte
%-----------------------------------------------
M_RG=r_co2*M_co2+r_n2*M_n2+r_o2*M_o2+r_h2o*M_h2o    %Molmasse des Rauchgases
roh_RGN=M_RG/22.4;                                  %Abgasdichte im Normzustand
R_RG=Rm/(M_RG/1000)                                 %Gaskonstante Rauchgas
roh_RG=pFR*10^2/(R_RG*TRG)                          %Abgasdichte bei Rauchgastemperatur
%-----------------------------------------------
%Berechnung von cp-Abgas 
%-----------------------------------------------
cp_rg1=r_co2*cp_co2+r_n2*cp_n2+r_o2*cp_o2+r_h2o*cp_h2o;       %[kJ/[kmol*K)]
cp_RG=cp_rg1/(22.4*roh_RGN);                                  %[kJ/(kg*K)]                  NOrmdichte oder Betriebsdichte??????
%------------------------------
%Berechnung div. Parameter
%------------------------------
w0=VB_brst*Sek/sekGD/(pi*(d0/1000)^2/4);
roh_brst=(p0+pN)*10^2/(R_brst*T0);                  %Brst-Dichte [kg/m³]
roh_L=(pL)*10^2/(R_L*TL)                            %Luft-Dichte [kg/m³]
GammaB=(Gamma*pi/180);                              %Umrechnung des Freistrahlwinkel [rad]
a=d0/(2*tan(GammaB/2));                             %Abstand Strahlursprungspunkt - KS [mm]
d1=2*(x+a)*tan(GammaB/2);                           %Strahldurchmesser an Stelle x [mm]
ymax=round(d1*0.5)+1;                                 %Obere Grenze y-Achse
ymin=round(-d1*0.5)-1;                                %Untere Grenze y-Achse
%------------------------------
%Berechnung der theoretischen Verbrennungstemperatur 
%------------------------------
m_brst=VN_brst*roh_brst;
m_L=VN_L*roh_L;
Hu=VN_brst*31800;
H_brst=m_brst*cp0*(T0-TN);
H_l=m_L*cpL*(TL-TN);
t_theo=(Hu+H_brst+H_l)/((m_brst+m_L)*cp_RG)+TN
%------------------------------
%Berechnung von Konzentrations-, Geschwindigkeits- und Temperaturverteilung
%------------------------------
w1=kf/(2*ci*(x/d0))*(roh_brst/roh_RG)^0.5*exp(-((1/ci^2)*(y./x)^2));
k1=kf*(roh_brst*w0/(roh_RG*w1*w0))/(4*cc^2*(x/d0)^2)*exp(-((y./x)^2)*1/cc^2);
T1=kf/(4*cq^2*(x/d0)^2)*(roh_brst/roh_RG)*(cp0/cp_RG)*(w0/(w1*w0))*exp(-((y./x)^2)*1/cq^2);
%------------------------------
%Berechnung der Massenzunahme
%------------------------------
m1=2*tan(GammaB/2)*(z/d0)*(roh_RG/roh_brst)^0.5;
%------------------------------
%Plotten der Ergebnisse
%------------------------------
 %
 figure(1)
% Punkte in xy-Ebene
pwx = ymin:1:ymax;
pwy = kf/(2*ci*(x/d0))*(roh_brst/roh_RG)^0.5*exp(-((1/ci^2)*(pwx./x).^2));
plot(pwx,pwy);
title('Geschwindigkeitsprofil');
xlabel('Strahldurchmesser [mm]');
ylabel('Geschwindigkeit [w1/w0]');
grid on;

% Rotieren um x-Achse
n=25;
t=(0:n)'*2*pi/n;
figure;
h1=surf(ones(n+1,1)*pwy,cos(t)*pwx,sin(t)*pwx); 
title('Geschwindigkeitsprofil');
ylabel('Strahldurchmesser [mm]');
xlabel('Geschwindigkeit [w1/w0]');
zlabel('Strahldurchmesser [mm]');
%
figure(2)
% Punkte in xy-Ebene
pkx = ymin:1:ymax;
pky = kf*roh_brst*w0/(roh_RG*w1.*w0)/(4*cc^2*(x/d0)^2)*exp(-((1/cc^2)*(pkx./x).^2));
plot(pkx,pky);
title('Konzentrationsverlauf');
xlabel('Strahldurchmesser [mm]');
ylabel('Konzentration [k1/k0]');
grid on;
% 
% Rotieren um x-Achse
n=25;
t=(0:n)'*2*pi/n;
figure;
h1=surf(ones(n+1,1)*pky,cos(t)*pkx,sin(t)*pkx); 
title('Konzentrationsprofil');
ylabel('Strahldurchmesser [mm]');
xlabel('Konzentration [k1/k0]');
zlabel('Strahldurchmesser [mm]');
%