% ________________________________________________________________________

%
%__________________________________________________________________________
%
%

% THE LIST OF SPECIES GRI 3.0
% O, O2, M, H, OH, H2, HO2, H2O2, CH, CO, CH2, HCO, CH2(S), CH3, CH2O, CH4, CO2, CH2OH, CH3O, CH3OH, C2H, C2H2, HCCO, C2H3, CH2CO, C2H4, C2H5, C2H6, H2O, N2, AR, C, HCCOH, N, NO, N2O, NO2, NH, HNO, NH2, NNH, CN, NCO, HCN, HOCN, HNCO, H2CN, HCNN, HCNO, NH3, CH2CHO, CH3CHO, C3H8, C3H7, (M)
cleanup;
clear all;


%%Enginedata
epsilon = 12; 
Bohrung = 85; %mm
r_KW = 52.5; %mm
f = 1500/60; %U/s
omega_Motor = 2*pi*f;
V_h = Bohrung^2/4*pi*2*r_KW/(1000^3); %[m^3]
V_c = V_h/(epsilon-1); %[m^3]
V_0 = V_h+V_c; %[m^3]
l_Pl = 240; %[mm]
lambda_s = r_KW/l_Pl;
A_Kolben = (Bohrung/1000)^2/4*pi; %[m^2]

m_Luft=450;         % kg/h Airmass
mH2OinL=0.2;        % kg/h Water in air




%% Gemischzusammensetzung
% Volumenanteile Luft (lt. ISO 2533, niedrigere Anteile vernachlässigt (siehe https://atemschutzlexikon.de/lexikon/a/atemluft/))
nL_N2 = 0.78084/(0.78084+0.209476+0.00934);
nL_O2 = 0.209476/(0.78084+0.209476+0.00934);
nL_Ar = 0.00934/(0.78084+0.209476+0.00934);
% Massenanteile Luft
M_Luft = (nL_N2*28+nL_O2*32+nL_Ar*40); %kg/kmol Molmasse Ansaugluft
ML_N2 = nL_N2*28/M_Luft;
ML_O2 = nL_O2*32/M_Luft;
ML_Ar = nL_Ar*40/M_Luft;

n_Luft=m_Luft/M_Luft;
nLf_H2O= mH2OinL/18;

n_Luft_feucht=n_Luft+nLf_H2O;

nLf_N2=(nL_N2*n_Luft)/(n_Luft+n_Luft_feucht);
nLf_O2=(nL_O2*n_Luft)/(n_Luft+n_Luft_feucht);
nLf_Ar=(nL_Ar*n_Luft)/(n_Luft+n_Luft_feucht);


Vdot_in = f*3600/2*V_h; %[m^3/h] %Volumenstrom Ansaugluft ("Ansaugen" nur jede zweite Umdrehung)
L_min =  8/ML_O2; %kg Luft/kg H2 mindestens notwendige Luftmasse für vollständige Verbrennung (32g/mol / (2*2g/mol) = 8)


n=0;
%%      Abfrage mit Variation Anfangstemperatur, Massenstrom, Lambda
T_abbruch=0;
i=0;
for k = 1:1 %(Lambda -)   
    for j = 1:1 %(Massenstrom kg/h) 200 - 1500
       
        for i = 1:1:3 %Variation des Wassermassenstroms durch Veränderung des Einspritzbeginns 
T_abbruch=0;
n=0;
 % Variation wateramount
 % water di
injector_open = (279-i*9)/180*pi;
injector_close = 270/180*pi;
injector_delta = mod(injector_close-injector_open, 4*pi);

injector_mass_rate= 0.0025;
injector_mass = (injector_delta*180/pi)/(f*360)*injector_mass_rate;


T_injector=293.15;

% Gemischzusammensetzung

m_H2_kgh = 5.2; %massflow hydrogen kg/h

Lambda = m_Luft/(m_H2_kgh*L_min);


mdot_in = m_Luft; % 760 %kg/h Luftmassenstrom Ansaugluft
ndot_in = mdot_in /(nLf_N2*28+nLf_O2*32+nLf_Ar*40+nLf_H2O*18); %kmol/h Stoffmenge Ansaugluft

%Stoffmengenanteil H2
% mdot_in_H2 = mdot_in/Lambda/L_min; %kg/h Massenstrom H2
mdot_in_H2=m_H2_kgh; % massflow hydrogen
ndot_H2 = mdot_in_H2/2; %kmol/h Stoffmenge H2
ndot_N2 = ndot_in * nLf_N2; %kmol/h Stoffmenge N2
ndot_O2 = ndot_in * nLf_O2; %kmol/h Stoffmenge O2
ndot_Ar = ndot_in * nLf_Ar; %kmol/h Stoffmenge Ar
ndot_H2O = ndot_in * nLf_H2O;

rhodot_in = (mdot_in+mdot_in_H2)/Vdot_in; %[kg/m^3] Dichte Ansaugluft

%Brennraumanfangsdruck
P0_Br=(2.3697)*10^5; %[Pa]
% P0_Br=1.3*10^5; %[Pa]

%Stoffmengenanteile Ansaugluft + Wasserstoff
n_H2 = ndot_H2 / (ndot_H2+ndot_N2+ndot_O2+ndot_Ar+ndot_H2O);
n_N2 = ndot_N2 / (ndot_H2+ndot_N2+ndot_O2+ndot_Ar+ndot_H2O);
n_O2 = ndot_O2 / (ndot_H2+ndot_N2+ndot_O2+ndot_Ar+ndot_H2O);
n_Ar = ndot_Ar / (ndot_H2+ndot_N2+ndot_O2+ndot_Ar+ndot_H2O);
n_H2O = ndot_H2O / (ndot_H2+ndot_N2+ndot_O2+ndot_Ar+ndot_H2O);

% Molanteile

Molanteile =  ['N2:' num2str(n_N2), ', O2:' num2str(n_O2), ', AR:' num2str(n_Ar),', H2:' num2str(n_H2),', H2O:' num2str(n_H2O)];
% Molanteile =  ['N2:' num2str(n_N2), ', O2:' num2str(n_O2), ', H2:' num2str(n_H2),', H2O:' num2str(n_H2O)];

% Startbedingungen Brennraum
T0_Br = 618; %[K] Anfangstemperatur Ansaugluft 

P_Br(1)=P0_Br;
T_Br(1)=T0_Br;

%% Reaktor
%Reaktoranfangszusammensetzung


    gas = GRI30;                                                                %Festlegen des Reaktionsmechanismus, Erstellen der reagierenden Komponente "gas"  
    set(gas,'T',T0_Br,'P',P0_Br,'X', Molanteile);                           %Zuweisen der Anfangsbedingungen/-zusammensetzung des Reaktionsgases im Reaktor
    r = Reactor(gas);                                                       %Erstellen eines Reaktors, gefüllt mit Komponente "gas"
 

 water =Water;

    set(water,'T',T_injector,'P',500000);

    rw = Reservoir(water);
    
    
   f_inj=MassFlowController(rw,r); % set connection between reservoir and reactor
   
   m_w=@(t)(mod(crank_angle(t,f)-injector_open,4*pi)< injector_delta)*injector_mass_rate; % define function for water inj timing

   C_P0 = cp_mole (gas);                                                   %Wärmekapazität bei konstantem Druck, Anfangswert
   C_V0 = cv_mole (gas);                                                   %Wärmekapazität bei konstantem Volumen, Anfangswert
   kappa = C_P0 / C_V0;                                                   %Isentropenexponent



a = IdealGasMix('air.cti');                                                 %Umgebungsgas, Luft
set(a,'P',101325, 'T', 573)                                                 %Anfangsbedingungen Umgebungsgas (Schätzwert für Kolbentemperatur, 300°C)
env = Reservoir(a);                                                         %Environment filled with gas
w = Wall;                                                                   %partition wall between reactor an environment
setArea(w, A_Kolben);                                                       %area for heat exchange
install(w,r,env);                                                           %Komponenten Wand, Reaktor, Umgebung installieren
setInitialVolume(r, V_0)                                                    %Anfangsvolumen Reaktor vor Kompression
setHeatTransferCoeff(w, 290) %lt.Kuchling: Taschenbuch der Physik           %W/m²/K Wärmeübergangskoeffizient Wand, Schätzwert für erzwungene Luftströmung

%% Schleife Kompression Reaktor
network = ReactorNet({r});                                                  %Erstellen eines Netzwerkes (notwendig für Aktivieren der Reaktorreaktion)
t = 0.0;   
dt = 1/(360*f);
t0 = cputime;
t_stop=(1/f)-dt;

Massenanteil_H2_vorher = massFraction(r, 'H2');%g 
Gesamtgewicht_Reaktanden_vorher_var= mass(r)*1000;
Masse_H2_vorher_var = Gesamtgewicht_Reaktanden_vorher_var*Massenanteil_H2_vorher;

 while  t<t_stop
     n=n+1;  
                                          %!Cantera-Fehlermeldungen bei zu großen Schritten
t = t + dt;  
% t=time(network);
advance(network,t);        
          
        
                          %Auslesen Massenanteil H2 vor Beginn Kompression       


                                                  %360° -> vollständige Motordrehung
    
    
                                %Start im UT 
    
%     f=MassFlowController(rw,r, mwtime);
    
% 
%     massFlowRate(f,0.013);

    C_P(n) = cp_mole (gas);                                                 %Auslesen der Gaskonstante im Reaktor 
    C_V(n) = cv_mole (gas);                                                 %Auslesen der Gaskonstante im Reaktor
    kappa = C_P(n) / C_V(n);                                                %Aktualisieren des Isentropenexponenten für nächsten Schleifendurchlauf

    C_P_Verlauf(n,k,j,i)= cp_mole (gas);  %g
    C_V_Verlauf(n,k,j,i)= cv_mole (gas); %g
    kappa_Verlauf (n,k,j,i)= C_P(n) / C_V(n);

    Kolbenweg = l_Pl+r_KW-r_KW*cos(crank_angle(t,f))-l_Pl*(1-lambda_s^2*sin(crank_angle(t,f))^2)^0.5; %Quelle: Vorlesung Motormechanik, Herleitung Massenkräfte
    Kolbenweg_t=l_Pl+r_KW-r_KW*cos(crank_angle(t+dt,f))-l_Pl*(1-lambda_s^2*sin(crank_angle(t+dt,f))^2)^0.5;


   f_inj_flow = Func('polynomial',0,m_w(t)); % define function behavior
   setFunction(f_inj,f_inj_flow); % set function for massflow controller

       V_Brennraum (n) = V_c+(Bohrung^2/4*pi*Kolbenweg)/1000^3;
       V_Brennraum (n+1) = V_c+(Bohrung^2/4*pi*Kolbenweg_t)/1000^3;
       
       Kolbengeschwindigkeit (n) = r_KW*omega_Motor*[sin(crank_angle(t,f))+lambda_s/2*sin(2*crank_angle(t,f))/(1-lambda_s^2*[sin(2*crank_angle(t,f))]^2)^0.5]/1000; %mm/s %Quelle: Vorlesung Motormechanik, Herleitung Massenkräfte
%     end
    P_Br(n+1) = P_Br(n) /(V_Brennraum(n+1)/V_Brennraum(n))^kappa;  %p_0 / (V_1/V_0)^k           %Berechnung Brennraumdruck
    T_Br(n+1) = T_Br(n) /(V_Brennraum(n+1)/V_Brennraum(n))^(kappa-1);                           %Berechnung Brennraumtemperatur
    f_pist = Func('polynomial',0,Kolbengeschwindigkeit (n));                     %Übergabe Kolbengeschwindigkeit an Cantera
    setVelocity(w,f_pist)                                                        %Übergabe Kolbengeschwindigkeit an Cantera
    KW(n) = crank_angle(t,f)*180/pi;                                                   %Auslesen für Plot
        

            

   
%%Auslesen Zustand/Zusammensetzung Reaktorinhalt                            %Lambda k 
                                                                            %Massenstrom j
                                                                            %Temperatur i

%             Gesamtgewicht_Reaktanden_vorher(k,j,i,n)= mass(r)*1000; %g      %Auslesen Masse H2 vor Beginn Kompression Reaktor (unabhängig von Reaktorzeitfortschritt) 
%             Gesamtgewicht_Reaktanden_vorher(n,k,j,i)= mass(r)*1000; %g
            Gesamtgewicht_Reaktanden_vorher(n,k,j,i)= Gesamtgewicht_Reaktanden_vorher_var;
%             Masse_H2_vorher(n,k,j,i) = Gesamtgewicht_Reaktanden_vorher(k,j,i,n)*Massenanteil_H2_vorher;%g %Auslesen Masse H2 vor Beginn Kompression Reaktor (unabhängig von Reaktorzeitfortschritt)
%             Masse_H2_vorher(n,k,j,i) = Gesamtgewicht_Reaktanden_vorher(n,k,j,i)*Massenanteil_H2_vorher;%g
            Masse_H2_vorher(n,k,j,i) = Masse_H2_vorher_var;
            Temperatur_Verlauf_KW(n,k,j,i) = temperature(r);                %Auslesen Reaktortemperatur (für jeden Reaktorzeitfortschritt, Schleifenfortschritt)
            Temperatur(i) = T0_Br;                                          %Reaktoranfangstemperatur, geändert durch Schleife, Parameter i
            Temperatur_Reaktor_end(k,j,i) = temperature(r);                 %Auslesen Reaktorendtemperatur nach Kompression im UT
            Anfangstemperatur_Brennraum(k,j,i)=T0_Br;                       %Zusatzinformation, nicht geplottet
            Temperaturerhoehung_Reaktor(k,j,i)= Temperatur_Reaktor_end(k,j,i) - Anfangstemperatur_Brennraum(k,j,i);%Zusatzinformation, nicht geplottet
            Massenstrom (j) = mdot_in;                                      %Auslesen Massenstrom geändert durch Schleife, Parameter j
            Luftverhaeltnis (k) = Lambda;                                   %Auslesen Luftverhältnis geändert durch Schleife, Parameter k
            Gesamtgewicht_Reaktanden(n,k,j,i)= mass(r)*1000; %g;            %Auslesen verbleibende Masse H2 im Reaktor (abhängig von Reaktorzeitfortschritt)
            Massenanteil_H2(n) = massFraction(r, 'H2');                     %Auslesen verbleibende Masse H2 im Reaktor (abhängig von Reaktorzeitfortschritt)
            Masse_H2(n,k,j,i) = Gesamtgewicht_Reaktanden(n,k,j,i)*(Massenanteil_H2(n));%g
            
            Massenanteil_H2O(n) = massFraction(r, 'h2o');
            Masse_H2O(n,k,j,i) = Gesamtgewicht_Reaktanden(n,k,j,i)*(Massenanteil_H2O(n));
            
            Anteil_H2_verbraucht (n,k,j,i) =  (Masse_H2_vorher(n,k,j,i) - Masse_H2(n,k,j,i))/Masse_H2_vorher(n,k,j,i); %massebezogen
            Molanteil_H2(n)= moleFraction(gas,'H2');                        %Auslesen verbleibender Molanteil H2 im Reaktor (abhängig von Reaktorzeitfortschritt)
            Molanteil_H2O(n)= moleFraction(gas,'H2O');                      %Auslesen verbleibender Molanteil H2O im Reaktor (abhängig von Reaktorzeitfortschritt)
            Molanteil_H2_min(k,j,i)= min(Molanteil_H2(n));                  %Auslesen minimaler Molanteil H2 im Reaktor (abhängig von Lambda/Massenstrom/Anfangstemperatur)
            Molanteil_H2O_max(k,j,i)= min(Molanteil_H2O(n));
            %Auslesen maximaler Molanteil H2O im Reaktor (abhängig von Lambda/Massenstrom/Anfangstemperatur)
            prod = netProdRates(gas);                                       %Variable zum Auslesen der Produktionsraten im Reaktor
            H2 = speciesIndex(gas, 'H2');                                   %Definition H2 zum Auslesen aus Cantera
            H2O = speciesIndex(gas, 'H2O');                                 %Definition H2O zum Auslesen aus Cantera
            Produktionsrate_H2(n,k,j,i) = prod(H2); %(kmol/m^3/s)           %Zusatzinformation, nicht geplottet
            Produktionsrate_H2O(n,k,j,i) = prod(H2O); %(kmol/m^3/s)         %Zusatzinformation, nicht geplottet
            min_Produktionsrate_H2(k,j,i)= min(Produktionsrate_H2(:,k,j,i));%Zusatzinformation, nicht geplottet
            max_Produktionsrate_H2O(k,j,i)= max(Produktionsrate_H2O(:,k,j,i));%Zusatzinformation, nicht geplottet
%             Matrix_P_Br(n,k,j,i) = P_Br(n);   %Alle Brennraumdrücke
            
            Matrix_P_Br(n,k,j,i) = pressure(r);

            mw_dot(n)=massFlowRate(f_inj,t);
            
            mw_inj_dur(i)=nnz(mw_dot(:));
            
            mw_rate(i)=(injector_mass*f/2)*3600;
            
            masse_H2O_Injektor(i)=injector_mass*1000;
            
            
    
for T_lauf = 1:max(size(Temperatur))
    alpha=1;
    h=1;
    EU10(T_lauf)=0;
while h==1
    if Anteil_H2_verbraucht(alpha,1,1,T_lauf) > 0.1
    EU10(T_lauf)=KW(alpha);
    h=0;
    else
        if alpha<max(size(Anteil_H2_verbraucht))
            alpha=alpha+1;
        else
            h=0;
            EU10(T_lauf)=NaN;
        end
    end
end

  alpha=1;
    h=1;
    
while h==1
    if Anteil_H2_verbraucht(alpha,1,1,T_lauf) > 0.01
    EU1(T_lauf)=KW(alpha)
    h=0;
    else
        if alpha<max(size(Anteil_H2_verbraucht))
            alpha=alpha+1;
        else
            h=0;
            EU1(T_lauf)=NaN;
        end
    end
end
end   
    
% if abs(min(EU10))>180
    T_abbruch=1;
% end
            
            

          

            end
        end
 
    

    end  
end