%KK5.m

function f = KK5(t,X,flag,u,d,T_BKs,h_DRs,h_KTs,T_DEs,u_0,N_OH,N_FD,N_K,N_EK,N_FV,T_OH_RG_0,T_OH_D_0,T_FD_RG_0,T_EK_RG_0,T_FV_RG_0)
global W V Y YY tt U
    
    %Pådrag og tilstander
    F_CH4 = u(1);           %mol/s
    %F_LUFT = u(2);           %mol/s
    %w_8 = u(3);             %kg/s
    %w_15 = u(4);            %kg/s
    w_12 = u(5);            %kg/s
    w_21 = u(6);            %kg/s
    X_BK = X(1:5);         %n_CH4 er her utelatt
    T_OH_RG = X(6:6+N_OH-1);     %Temperaturer i de 3 cellen på varm side i OH
    T_OH_D = X(6+N_OH:6+2*N_OH-1);     %Temperaturer i de 3 cellen på kald side i OH
    T_FD_RG = X(6+2*N_OH:6+2*N_OH+N_FD-1); %K, temp. i de 3 cellene på varm side i FD
    T_FD_H2O = X(6+2*N_OH+N_FD);        %Temperatur i damp fra drum, vannet i FD er konstant lik uttemp.
    h_DR = X(6+2*N_OH+N_FD+1);   %m, væskenivå i drum
    T_K_cw = X(6+2*N_OH+N_FD+2:6+2*N_OH+N_FD+2+N_K-1); %K, kjølevannstemperatur i cellene i K
    T_K_D = X(6+2*N_OH+N_FD+2+N_K); %Kondenseringstemp. i Kondensator
    h_KT = X(6+2*N_OH+N_FD+2+N_K+1); %m, væskenivå i KT
    T_DE = X(6+2*N_OH+N_FD+2+N_K+2); %K, temp. i deaerator
    T_EK_RG = X(6+2*N_OH+N_FD+2+N_K+3:6+2*N_OH+N_FD+2+N_K+3+N_EK-1);   %K, temp. i de 3 cellene, varmside, EK
    T_EK_H2O = X(6+2*N_OH+N_FD+2+N_K+3+N_EK:6+2*N_OH+N_FD+2+N_K+3+2*N_EK-1); %K, temp. i de 3 cellene, kaldside, EK
    T_FV_RG = X(6+2*N_OH+N_FD+2+N_K+3+2*N_EK:6+2*N_OH+N_FD+2+N_K+3+2*N_EK+N_FV-1); %K, temp. i de 3 cellene, varmside, FV
    T_FV_H2O = X(6+2*N_OH+N_FD+2+N_K+3+2*N_EK+N_FV:6+2*N_OH+N_FD+2+N_K+3+2*N_EK+2*N_FV-1); %K, temp. i de 3 cellene, kaldside, FV

%--------------------------------------------------------------------------------------------------------------    
    %Konstanter
    R = 8.3145;                     %J/(K*mol) fra SICD [8]
    Vm = 23.690E-3;                 %m3/mol, molart volum ved 1 atm og 288.7 K
    eta_KOMPR = 0.9;                %Polytropisk kompressorvirkningsgrad (fra Bolland (høstprosjekt))
    eta_GT = 0.9;                   %Adiabatisk turbinvirkningsgrad (fra Bolland (høstprosjekt))
    eta_DT = 0.9;                   %Adiabatisk turbinvirkningsgrad (fra Bolland (høstprosjekt))
    eta_p = 0.75;                   %Adiabatisk pumpevirkningsgrad (fra Bolland (høstprosjekt))
    delta_rH_ref = -803E3;          %J/mol, forbrenningsvarme fra delta_fH i SICD, for T_0=298 K
    delta_rcp = 11;                 %J/(K*mol), fra cp-verdier i SICD
    delta_vapH_ref = 2.4444E6;      %J/kg, fra SICD [8], fordampningsentalpi ved T_ref=298K 
    T_ref = 298;                    %K, referansetemperatur for entalpievalueringer
    
    Mm = [32.0; 28.0; 44.0; 18.0; 16.0];    %g/mol  Mm =[O2; N2; CO2; H20; CH4] verdier fra SICD [8]
    cp = [29; 29; 37; 34; 36];              %J/(K*mol) cp =[O2; N2; CO2; H20(g); CH4] verdier fra [8]
    cp_D = 1888;                            %J/kg,K, cp for vanndamp, [8]
    cp_s = 473;                             %J/kg,K, cp for stål, verdi hentet fra Geankoplis [13]
    cp_l = 4180;                            %J/kg,K, cp for H2O(l), [6]
    rho_s = 7801;                           %kg/m3, tetthet av stål, [13]
    rho_l = 1000;                           %kg/m3, tetthet av vann, [8]
    x_LUFT = [0.21; 0.79; 0; 0; 0];         % mol% (Luft-sammensetning)
    x_CH4 = [0; 0; 0; 0; 1];                % mol% (Metanføde-sammensetning)
    
    cp_LUFT = cp'*x_LUFT;
    gamma_LUFT = cp_LUFT/(cp_LUFT-R);   %adiabateksponent for luft
    gamma_CH4 = cp(5)/(cp(5)-R);        %adiabateksponent for metan
    gamma_D = cp(4)/(cp(4)-R);          %adiabateksponent for damp
    
    V_BK = 0.4;                         %m3 BK-volum
    
    d_tube = 0.05;      %m, rørdiameter, verdi hentet fra [2]
    dd_tube = 0.005;    %m, tykkelse av rørvegg, verdi hentet fra [2]
    A_DE = 36;          %m2, grunnflate av deaerator, fra [2]
    A_DR = 10;          %m2, grunnflate av drum, fra [2]
    A_EK = 15500;       %m2, varmeoverføringsareal i EK, verdi hentet fra [2]
    A_FD = 21000;       %m2, varmeoverføringsareal i fordamper, fra [2]
    A_FV = 2300;        %m2, varmeoverføringsareal i forvarmer, fra [2]
    A_K = 5000;         %m2, varmeoverføringsareal i kondensator, fra [2]
    A_KT = 1;           %m2, grunnflate av kondensatortanken, [2]
    A_OH = 12000;       %m2, varmeoverføringsareal i OH, verdi hentet fra [2]
    
    U_EK_0 = 100;   %W/K,m2 varmegjennomgangstall for ekonomiser ved kjente betingelser, fra [2]
    U_FD_0 = 100*3;   %W/K,m2 varmegjennomgangstall for fordamper ved kjente betingelser, fra [2]
    U_FV_0 = 100;   %W/K,m2 varmegjennomgangstall for forvarmer ved kjente betingelser, fra [2]
    U_K = 4000*3;     %W/K,m2 verdi hentet fra [2]
    U_OH_0 = 80*3;    %W/K,m2, varmegjennomgangstall for OH ved kjente betingelser, verdi hentet fra [2]
    w_4_0 = 500;	%Røykgassmengde ved kjent varmegjennomgangstall for OH, ref [2]
    w_20_0 = 80;    %kg/s, [2]
    w_21_0 = 1.5;   %kg/s, [2]
    m_H2O_tot = 65000;  %kg, totalmasse for vann i de 3 tankene (KT, DR, DE)
    w_cw = 10000;   %kg/s, [2]
    T_cw = 283;     %K, kjølevannstemp. inn på kondensator [2]
    T_FD_H2O_0 = 584;   %K, Referanseverdi fra [2] for å beregne P_DR_0
    T_20_0 = 793;   %K, Referanseverdi fra [2] for å beregne rho_DT_0 og w_DT_0
    m=0.57;			%Eksponent for beregning av U_OH fra U_OH_0, w_4 og w_4_0, ref [2] og [7]
    
    P_4 = 1.013E5;      %Pa, trykk i strøm 4 ut av GT
        
    %Priser
    P_el = 0.15;                        %NOK/kWh, (verdi fra Bolland)
    P_NG = 0.70;                        %NOK/Sm3, (verdi fra Bolland)
    
%-----------------------------------------------------------------------------------------------  
%Regulering 

%Beregning av reguleringsavvik
    e_T_BK = T_BKs - X_BK(5);     
    e_h_DR = h_DRs - h_DR;
    e_h_KT = h_KTs - h_KT;
    e_T_DE = T_DEs - T_DE;
    
    e_T_BK_int = X(35);    %Integrert reguleringsavvik 
    e_h_DR_int = X(36);    %Integrert reguleringsavvik 
    e_h_KT_int = X(37);    %Integrert reguleringsavvik 
    e_T_DE_int = X(38);    %Integrert reguleringsavvik
    
if 1==1
%Temperatur i BK
Kc_T_BK = 0.06;
tau_T_BK = 0.004;
F_LUFT_0 = u_0(2);

F_LUFT = F_LUFT_0 - Kc_T_BK*(e_T_BK + e_T_BK_int/tau_T_BK);  %mol/s, dersom T_BK > T_BKs ==> reduseres F_CH4
    if F_LUFT >= 2*F_LUFT_0
        F_LUFT = 2*F_LUFT_0;
    elseif F_LUFT <= 0
        F_LUFT = 0;
    end
end

if 1==1
%Væskenivå i drum
Kc_h_DR = 200;
tau_h_DR = 200; %s
w_15_0 = u_0(4);    %kg/s

w_15 = w_15_0 + Kc_h_DR*(e_h_DR + e_h_DR_int/tau_h_DR); %kg/s, dersom h_DR > h_DRs ==> reduseres innstrømmen w_15 
    if w_15 >= 2*w_15_0
        w_15 = 2*w_15_0;
    elseif w_15 <= 0
    w_15 = 0;
    end
end

if 1==1
%Væskenivå i kondensatortank
Kc_h_KT = 100;
tau_h_KT = 40;  %s
w_8_0 = u_0(3); %kg/s

w_8 = w_8_0 - Kc_h_KT*(e_h_KT + e_h_KT_int/tau_h_KT);   %kg/s, dersom h_KT > h_KTs ==> øker utstrømmen w_8
    if w_8 >= 2*w_8_0
        w_8 = 2*w_8_0;
    elseif w_8 <= 0
    w_8 = 0;
    end
end

if 1==1
%Temperatur i deaerator
Kc_T_DE = 45.5;
tau_T_DE = 4000;  %s
F_CH4_0 = u_0(1);   %mol/s

F_CH4 = F_CH4_0 + Kc_T_DE*(e_T_DE + e_T_DE_int/tau_T_DE); %mol/s, dersom T_DE > T_DEs ==> reduseres F_CH4
    if F_CH4 >= 2*F_CH4_0
        F_CH4 = 2*F_CH4_0;
    elseif F_CH4 <= 0
        F_CH4 = 0;
    end
end

%-----------------------------------------------------------------------------------------------
    %Brennkammer, BK (her beregnes trykket ut av kompressorene)
    %Antar all hold-up i GT skjer i brennkammeret.
    %Antar forbrenningen skjer momentant og fullstendig, og 100% omsetning av CH4
    
    n_BK = [X_BK(1:4); 0];              %mol i BK av komponentene, n_CH4 settes her =0
    T_BK = X_BK(5);                 %K 
    
    n_BKtot = sum(n_BK);            %mol i BK totalt
    x_BK = n_BK/n_BKtot;            % mol% (sammensetning i BK)
    
    p_BK = n_BKtot*R*T_BK/V_BK;     %Pa 
    
    ny = [-2; 0; 1; 2; -1];          %reaksjonskoeffisienter, [O2; N2; CO2; H2O; CH4] 
%------------------------------------------------------------------------------------------------    
    %Her beregnes flow gjennom kompressorene og temp. ut av kompressorene
        %Metankompressor
            %Flow og Temp_ut
            Ws_CH4 = F_CH4*cp(5)*d(3)*((p_BK/d(1))^((gamma_CH4-1)/gamma_CH4)-1)/eta_KOMPR; %J/s
            %F_CH4 = (u(1)*eta_KOMPR)/(cp(5)*d(3)*((p_BK/d(1))^((gamma_CH4 - 1)/gamma_CH4) - 1));   %mol/s
            F_2 = F_CH4;                                            %mol/s (ingen hold-up i kompressor)
        
        if F_CH4 > 0    
            T_2 = d(3) + Ws_CH4/(F_CH4*cp(5));     %K, temp. i strøm 2, d(3)=T_CH4
        else T_2 = d(3);
        end
        
        if 1==1    
            if Ws_CH4 <= 0      %Hvis p_BK<P_CH4, antas det at CH4 ekspanderes gjennom metankompr uten at det tas ut arbeid,
                Ws_CH4 = 0;     %og temperaturen endres som gitt fra ideell gasslov, Tut=Tinn*Put/Pinn
                T_2 = d(3)*p_BK/d(1);
            end
        end
        
        %Luftkompressor
            %Flow og Temp_ut
            Ws_LUFT = F_LUFT*cp_LUFT*d(4)*((p_BK/d(2))^((gamma_LUFT-1)/gamma_LUFT)-1)/eta_KOMPR; %J/s
            %F_LUFT = (u(2)*eta_KOMPR)/(cp_LUFT*d(4)*((p_BK/d(2))^((gamma_LUFT - 1)/gamma_LUFT) - 1)); %mol/s
            F_1 = F_LUFT;                                           %mol/s (ingen hold-up i kopressor)
              
            T_1 = d(4) + Ws_LUFT/(F_LUFT*cp_LUFT); %K, temp. i strøm 1, d(4)=T_LUFT
%---------------------------------------------------------------------------------------------------            
    %Turbin (her beregnes massestrøm ut fra BK)
        p_3_0 = 3.0E6;            %Pa, verdi fra høstprosjektet
        rho_3_0 = 5.75;         %kg/m3, verdi beregnet fra høstprosjekt samt oppholdstid (Bolland)
        w_3_0 = (1.73E6)/3600;  %kg/s, verdi hentet fra høstprosjektet
        p_3 = p_BK;             %Pa
        rho_3 = (n_BK'*Mm)/(V_BK*1000);    %kg/m3    
        w_3 = w_3_0*sqrt((p_3*rho_3)/(p_3_0*rho_3_0));   %kg/s, (fra Hyllseth (perm))
        w_4 = w_3;  %kg/s, Antar ingen hold-up i turbinen
%----------------------------------------------------------------------------------------------------        
    %Brennkammer
        delta_rH_T_BK = delta_rH_ref + delta_rcp*(T_BK - T_ref);      %J/mol (CH4)
    
        F_3 = (w_3*1000)/(Mm'*x_BK);                                %mol/s
        %Materialbalanse   
        DnDt = F_LUFT*x_LUFT + F_CH4*x_CH4 - F_3*x_BK + F_CH4*ny;   %mol/s
        
        %Energibalanse
        DT_BKDt = (F_1*cp_LUFT*(T_1 - T_BK) + F_2*cp(5)*(T_2 - T_BK) + (-delta_rH_T_BK)*F_CH4)/(n_BK'*cp);
%-----------------------------------------------------------------------------------------------------------        
%Turbin (her beregnes arbeidet som tas ut i turbinen)
        gamma_3 = (cp'*x_BK)/((cp'*x_BK) - R);      %adiabateksponent
        Wsrev_GT = F_3*(cp'*x_BK)*T_BK*((P_4/p_3)^((gamma_3 - 1)/gamma_3) - 1); %J/s, reversibelt arbeid
        Ws_GT = Wsrev_GT*eta_GT;                    %J/s, arbeid som tas ut i GT
        
        cp_4 = cp'*x_BK;                            %J/(K*mol)
        T_4 = T_BK + Ws_GT/(F_3*cp_4);              %K

%---------------------------------------------------------------------------------------------------------------------     
%Drum og fordamper

%Fysiske data for fordamperen
V_FD_D = d_tube*A_FD/4;         %m3, volum av fordamper på dampsiden
V_FD_RG = d_tube*A_FD/2;		%m3, dobbeltvolum i forhold til dampsiden [2]
m_FD_s = A_FD*dd_tube*rho_s;    %kg, masse av stål i rørvegg (men hva med drumveggen?? kanskje dd_tube er tilpasset?)

U_FD = U_FD_0*(w_4/w_4_0)^m;     %W/K,m2 varmegjennomgangstall for fordamper

delta_vapH_T_FD_H2O = delta_vapH_ref + (cp_D - cp_l)*(T_FD_H2O - T_ref); %J/kg

%Dampsiden
T_16 = T_EK_H2O(1,1);   %K, temp. i vann ut av EK
m_FD_H2O = (A_DR*h_DR + V_FD_D*0.5)*rho_l;  %kg, hold-up av H2O i fordamper+drum
                                            %antar at halve fordamperen er fylt med vann (neglisjerer damp)

%Røykgassiden
n_FD_RG = zeros(N_FD,1);
i = [1:N_FD];
    n_FD_RG(i,1) = P_4*(V_FD_RG/N_FD)./(R*T_FD_RG_0(i,1)); %mol, hold-up røykgass i fordamper
    Q_FD(i,1) = U_FD*(A_FD/N_FD)*(T_FD_RG(i,1)-T_FD_H2O); %J/s, tilført varme, slipper delta_T_lm pga antar celler med perfekt mix (T=T_ut)

Q_FD_tot = sum(Q_FD(:,1)); %J/s, (kan være denne må endres??)

%---------------------------------------------------------------------------------------------------------------------
% Beregner strøm ut av drum, dvs. kritisk strømning gjennom dampturbinen
T_20 = T_OH_D(N_OH);    %K, damptemp. ut fra OH
P_DR = trykk(T_FD_H2O); %Pa, metningstrykk ved damptemp. ut av fordamper/drum
rho_DT = Mm(4)*P_DR/(R*T_20*1000);  %kg/m3, tetthet av dampen inn på dampturbin
P_DR_0 = trykk(T_FD_H2O_0); %Referansetrykk for å benyttet "choked nozzel"-likningen
rho_DT_0 = Mm(4)*P_DR_0/(R*T_20_0*1000);  %kg/m3,referanseverdi for tetthet av dampen inn på dampturbin

w_DT = (w_20_0-w_21_0)*sqrt(P_DR*rho_DT/(P_DR_0*rho_DT_0));         %kg/s, "choked nozzel"-likning
w_19 = w_DT + w_21; %kg/s
%F_19 = w_19/Mm(4)*1000; %mol/s

%--------------------------------------------------------------------------------------------------------------------
%Pumpe, P2, mellom deaerator og ekonomiser
P_DE = trykk(T_DE); %Pa, trykk i deaerator beregnet på bakgrunn av temp. i deaerator
Ws_P2 = w_15*(P_DR-P_DE)/(rho_l*eta_p);

if w_15 > 0
    T_15 = T_DE + Ws_P2/(w_15*cp_l);        %K, temp. i strøm 15 ut av P2
else T_15 = T_DE;
end
%--------------------------------------------------------------------------------------------------------------------
%Drum og fordamper
%Diff.likninger
DT_FD_RGDt = zeros(N_FD,1);
DT_FD_RGDt(1,1) = (F_3*cp_4*(T_OH_RG(N_OH) - T_FD_RG(1,1)) - Q_FD(1,1))/(n_FD_RG(1,1)*cp_4);    %K/s
i = [2:N_FD];
    DT_FD_RGDt(i,1) = (F_3*cp_4*(T_FD_RG(i-1,1) - T_FD_RG(i,1)) - Q_FD(i,1))./(n_FD_RG(i,1)*cp_4);   %K/s

DT_FD_H2ODt = (w_15*cp_l*(T_16 - T_FD_H2O) - w_19*delta_vapH_T_FD_H2O + Q_FD_tot)/(m_FD_H2O*cp_l + m_FD_s*cp_s);  %K/s

Dh_DRDt = (w_15 - w_19)/(rho_l*A_DR); %m/s

%-------------------------------------------------------------------------------------------------------------------
%Overheter
    %Fysiske data for OH
    V_OH_D = A_OH*d_tube/4;         %m3, volum inni rørene
    V_OH_RG = A_OH*d_tube/2;        %m3, dobbelt volum utenfor rørene, ref.[2]
    m_s = A_OH*dd_tube*rho_s;       %kg, masse av stål i rørvegg

U_OH = U_OH_0*(w_4/w_4_0)^m; %W/K,kg Varmegjennomgangstall for OH ved w_4    

%Røykgass-siden, varm side
n_OH_RG = zeros(N_OH,1);
Q_OH = zeros(N_OH,1);
i = [1:N_OH];
    n_OH_RG(i,1) = P_4*(V_OH_RG/N_OH)./(R*T_OH_RG_0(i,1)); %mol, hold-up røykgass i hver av cellene
    Q_OH(i,1) = U_OH*A_OH/N_OH*(T_OH_RG(i,1)-T_OH_D(i,1)); %J/s, tilført varme, slipper delta_T_lm pga antar celler med perfekt mix (T=T_ut)

%Dampsiden, kald side
n_OH_D = zeros(N_OH,1);
i = [1:N_OH];
    n_OH_D(i,1) = P_DR*(V_OH_D/N_OH)./(R*T_OH_D_0(i,1)); %mol, hold-up damp i hver av cellene

%Diff.likninger
DT_OH_RGDt = zeros(N_OH,1);
DT_OH_DDt =  zeros(N_OH,1);

DT_OH_RGDt(1,1) = (F_3*cp_4*(T_4 - T_OH_RG(1,1)) - Q_OH(1,1))/(n_OH_RG(1,1)*cp_4); %Tar ikke med "ut"-leddet pga at T=T_ut (perfekt mix i hver celle)
i = [2:N_OH];
    DT_OH_RGDt(i,1) = (F_3*cp_4*(T_OH_RG(i-1,1) - T_OH_RG(i,1)) - Q_OH(i,1))./(n_OH_RG(i,1)*cp_4);

DT_OH_DDt(N_OH,1) = (w_19*cp_l*(T_FD_H2O - T_OH_D(N_OH,1)) + Q_OH(N_OH,1))/(n_OH_D(N_OH,1)*cp(4) + m_s/N_OH*cp_s); %Tar heller ikke her med "ut"-leddet
i = [1:N_OH-1];
    DT_OH_DDt(i,1) = (w_19*cp_l*(T_OH_D(i+1) - T_OH_D(i,1)) + Q_OH(i,1))./(n_OH_D(i,1)*cp(4) + m_s/N_OH*cp_s); %Tar heller ikke her med "ut"-leddet

%----------------------------------------------------------------------------------------------------------------
%Dampturbin, DT
P_23 = trykk(T_K_D);    %Beregner metningstrykket til dampen ved temperaturen på varm side i kondensatoren
Ws_DT = w_DT*cp_D*T_OH_D(1,1)*((P_23/P_DR)^((gamma_D-1)/gamma_D)-1)*eta_DT; %J/s, reelt arb. fra DT (neg.)

delta_vapH_T_K_D = delta_vapH_ref + (cp_D - cp_l)*(T_K_D - T_ref); %J/kg

T_23_utenkondens = T_OH_D(1,1) + Ws_DT/(w_DT*cp_D); %K, teoretisk temp. i strøm 23 uten kondensering
x_kondens = cp_D*(T_K_D - T_23_utenkondens)/delta_vapH_T_K_D; %Massebrøk av strøm 23 som er kondensert

%----------------------------------------------------------------------------------------------------------------
%Pumpe, P1, mellom kondensatortank og forvarmer
Ws_P1 = w_8*(P_DE-P_23)/(rho_l*eta_p);
%w_8 = Ws_P1*eta_p*rho_l/(P_DE-P_23); %kg/s, massestrøm ut av P1
if w_8 > 0
    T_9 = T_K_D + Ws_P1/(w_8*cp_l); %K, temp. i strøm 9 ut av P1
else T_9 = T_K_D;
end
%------------------------------------------------------------------------------------------------------------
%Kondensator
%Fysiske data for K
V_K_D = d_tube*A_K/4; %m3, volum av dampsiden i kondensatoren
V_K_cw = d_tube*A_K/2;  %m3, volum av kjølevannsiden dobbelt så stort som dampsiden
m_K_s = A_K*dd_tube*rho_s;  %kg, tilpasset stålmasse i rørveggen i kondensatoren, [2]

%Hold-up
m_cw = V_K_cw*rho_l;    %kg, masse av kjølevannet (cooling water)
m_K_H2O = (A_KT*h_KT + 0.5*V_K_D)*rho_l; %kg, antar at halvparten av volumet på dampsiden er fylt med vann 
                                            %og neglisjerer massen av dampen
i = [1:N_K];
    Q_K(i,1) = U_K*A_K/N_K*(T_K_D - T_K_cw(i,1));   %J/s, varme overført fra dampsiden til kjølevannsiden i celle i

Q_K_tot = sum(Q_K(:,1));  %J/s, totalt overført varme i kondensator

%Diff.likninger
DT_K_cwDt  = zeros(N_K,1);
i = [1:N_K-1];
    DT_K_cwDt(i,1) = (w_cw*cp_l*(T_K_cw(i+1,1) - T_K_cw(i,1)) + Q_K(i,1))./(m_cw/N_K*cp_l);  %K/s

DT_K_cwDt(N_K,1) = (w_cw*cp_l*(T_cw - T_K_cw(N_K,1)) + Q_K(N_K,1))/(m_cw/N_K*cp_l); %K/s

DT_K_DDt = (w_DT*(1-x_kondens)*delta_vapH_T_K_D - Q_K_tot)/(m_K_H2O*cp_l + m_K_s*cp_s);   %K/s

Dh_KTDt = (w_DT - w_8)/(rho_l*A_KT);    %m/s, nivåendring i kondensatortanken

%----------------------------------------------------------------------------------------------------------------
%Ekonomiser
%Fysiske data
V_EK_H2O = A_EK*d_tube/4;   %m3, volum av vannsiden i ekonomiseren
V_EK_RG = 2*V_EK_H2O;       %m3, volum av røykgassiden i EK, dobbelt volum i forhold til vannsiden
m_EK_s = A_EK*dd_tube*rho_s;    %kg, tilpasset stålmasse i rørveggen i EK

U_EK = U_EK_0*(w_4/w_4_0)^m;    %W/K,kg Varmegjennomgangstall for EK ved w_4

%Røykgassiden, varm side
n_EK_RG = zeros(N_EK,1);
Q_EK = zeros(N_EK,1);
i = [1:N_EK];
    n_EK_RG(i,1) = P_4*(V_EK_RG/N_EK)./(R*T_EK_RG_0(i,1)); %mol, hold-up røykgass i hver av cellene i EK
    Q_EK(i,1) = U_EK*A_EK/N_EK*(T_EK_RG(i,1)-T_EK_H2O(i,1));    %J/s, tilført varme, slipper delta_T_lm pga antar celler med perfekt mix (T=T_ut)

%Vannsiden
m_EK_H2O = V_EK_H2O*rho_l;  %kg, hold-up vann i hele EK

%Diff.likninger
DT_EK_RGDt = zeros(N_EK,1);
DT_EK_H2ODt = zeros(N_EK,1);

DT_EK_RGDt(1,1) = (F_3*cp_4*(T_FD_RG(N_FD,1)-T_EK_RG(1,1)) - Q_EK(1,1))/(n_EK_RG(1,1)*cp_4);
i = [2:N_EK];
    DT_EK_RGDt(i,1) = (F_3*cp_4*(T_EK_RG(i-1,1)-T_EK_RG(i,1)) - Q_EK(i,1))./(n_EK_RG(i,1)*cp_4);

i = [1:N_EK-1];
    DT_EK_H2ODt(i,1) = (w_15*cp_l*(T_EK_H2O(i+1,1) - T_EK_H2O(i,1)) + Q_EK(i,1))/((m_EK_H2O*cp_l+m_EK_s*cp_s)/N_EK);

DT_EK_H2ODt(N_EK,1) = (w_15*cp_l*(T_15 - T_EK_H2O(N_EK,1)) + Q_EK(N_EK,1))/((m_EK_H2O*cp_l+m_EK_s*cp_s)/N_EK);

%------------------------------------------------------------------------------------------------------------
%Forvarmer, FV
%Fysiske data
V_FV_H2O = A_FV*d_tube/4;   %m3, volum av vannsiden i forvarmeren
V_FV_RG = 2*V_FV_H2O;       %m3, volum av røykgassiden i FV, dobbelt volum i forhold til vannsiden
m_FV_s = A_FV*dd_tube*rho_s;    %kg, tilpasset stålmasse i rørveggen i FV

U_FV = U_FV_0*(w_4/w_4_0)^m;    %W/K,kg Varmegjennomgangstall for FV ved w_4

%Røykgassiden, varm side
n_FV_RG = zeros(N_FV,1);
Q_FV = zeros(N_FV,1);
i = [1:N_FV];
    n_FV_RG(i,1) = P_4*(V_FV_H2O/N_FV)./(R*T_FV_RG_0(i,1));    %mol, hold-up røykgass i hver av cellene i FV
    Q_FV(i,1) = U_FV*A_FV/N_FV*(T_FV_RG(i,1) - T_FV_H2O(i,1));%J/s, tilført varme, slipper delta_T_lm pga antar celler med perfekt mix (T=T_ut)

%Vannsiden, kald side
w_10 = w_8 + w_12;  %kg/s, massestrøm 10
if w_10 > 0
    T_10 = (T_9*w_8 + T_FV_H2O(1,1)*w_12)/w_10; %K, beregnes fra (w*cp*T)_inn = (w*cp*T)_ut, cp_l er konstant og forkortet bort,
else T_10 = T_9;
end
m_FV_H2O = V_FV_H2O*rho_l;  %kg, hold-up vann i hele FV

%Diff.likninger
DT_FV_RGDt = zeros(N_FV,1);
DT_FV_H2ODt = zeros(N_FV,1);

DT_FV_RGDt(1,1) = (F_3*cp_4*(T_EK_RG(N_EK,1)-T_FV_RG(1,1)) - Q_FV(1,1))/(n_FV_RG(1,1)*cp_4);
i = [2:N_FV];
    DT_FV_RGDt(i,1) = (F_3*cp_4*(T_FV_RG(i-1,1)-T_FV_RG(i,1)) - Q_FV(i,1))./(n_FV_RG(i,1)*cp_4);

i = [1:N_FV-1];
    DT_FV_H2ODt(i,1) = (w_10*cp_l*(T_FV_H2O(i+1,1) - T_FV_H2O(i,1)) + Q_FV(i,1))/((m_FV_H2O*cp_l+m_FV_s*cp_s)/N_EK);

DT_FV_H2ODt(N_FV,1) = (w_10*cp_l*(T_10 - T_FV_H2O(N_FV,1)) + Q_FV(N_FV,1))/((m_FV_H2O*cp_l+m_FV_s*cp_s)/N_EK);

%-------------------------------------------------------------------------------------------------------------
%Deaerator, DE

%Hold-up
m_DE_H2O = m_H2O_tot - rho_l*(h_DR*A_DR + h_KT*A_KT);    %kg, hold-up i DE blir her beregnet fra tot.massebalanse for de 3 tankene
h_DE = m_DE_H2O/(rho_l*A_DE);               %m, væskenivå er da avhengig av dynamikken i de andre tankene

%Energibetrakninger
delta_vapH_T_DE = delta_vapH_ref + (cp_D - cp_l)*(T_DE - T_ref); %J/kg, fordampningsentalpi ved T_DE
delta_H_D = cp_D*(T_20 - T_DE) + delta_vapH_T_DE;   %J/kg, entalpibidrag fra dampen i strøm 21(avkjøling av damp og kondensering)

%Diff.likninger
DT_DEDt = (w_21*delta_H_D + w_8*cp_l*(T_FV_H2O(1,1)-T_DE))/(m_DE_H2O*cp_l);

%------------------------------------------------------------------------------------------------------------
%Netto el-produksjon
 sigma_el = Ws_GT + Ws_DT + Ws_CH4 + Ws_LUFT + Ws_P1 + Ws_P2;                  %J/s (neg.)
       
        
%Kostnadsfunksjon, legges inn for å sjekke optimaliseringen
 J = -3600*((P_el/3600)*(-sigma_el/1000) - P_NG*Vm*F_CH4);  %NOK/s
  
  
%Returnerer       
f1 = DnDt(1:4);     %tilstand 1-4
f2 = DT_BKDt;       %tilstand 5
f3 = DT_OH_RGDt;    %tilstand 6-8
f4 = DT_OH_DDt;     %tilstand 9-11
f5 = DT_FD_RGDt;    %tilstand 12-14
f6 = DT_FD_H2ODt;   %tilstand 15
f7 = Dh_DRDt;       %tilstand 16
f8 = DT_K_cwDt;     %tilstand 17-19
f9 = DT_K_DDt;      %tilstand 20
f10 = Dh_KTDt;      %tilstand 21
f11 = DT_DEDt;      %tilstand 22
f12 = DT_EK_RGDt;   %tilstand 23-25
f13 = DT_EK_H2ODt;  %tilstand 26-28
f14 = DT_FV_RGDt;   %tilstand 29-31
f15 = DT_FV_H2ODt;  %tilstand 32-34

De_T_BKDt = T_BKs - T_BK;   %Reguleringsavvik integreres med tid ved å legges inn i f, tilstand 35
De_h_DRDt = h_DRs - h_DR;   %Reguleringsavvik integreres med tid ved å legges inn i f, tilstand 36
De_h_KTDt = h_KTs - h_KT;   %Reguleringsavvik integreres med tid ved å legges inn i f, tilstand 37
De_T_DEDt = T_DEs - T_DE;   %Reguleringsavvik integreres med tid ved å legges inn i f, tilstand 38

f = [f1; f2; f3; f4; f5; f6; f7; f8; f9; f10; f11; f12; f13; f14; f15; De_T_BKDt; De_h_DRDt; De_h_KTDt; De_T_DEDt];

%Globale variable  
%V = [J; sigma_el; p_BK; rho_3; uX; T_4; w_4; F_3; cp_4; T_OH_RG(3)];     %Denne tas med bare for å sjekke disse verdiene
V = [J; sigma_el; Ws_P1; Ws_P2; Ws_CH4; Ws_LUFT; p_BK; rho_3; F_CH4; F_LUFT; w_8; w_15; w_12; w_21; X; T_4; h_DE; x_kondens; P_23; P_DE; P_DR]; %samme som Z
W = [T_BK; h_KT; h_DR];      %målinger/utganger
Y = [F_CH4, F_LUFT, w_8, w_15, J];
YY = [YY; Y];
tt = [tt; t];
U = [F_CH4; F_LUFT; w_8; w_15; w_12; w_21; X(1:34)];    %Pådrag og tilstander sendes tilbake som global variabel for evt å lagres
if 1==0
    disp(sprintf('Name T F Q'));
    for i=1:3
       disp(sprintf('T_OH_RG: %6.4f %6.4f',T_OH_RG(i),w_4));
    end
    for i=1:3
       disp(sprintf('T_FD_RG: %6.4f %6.4f',T_FD_RG(i),w_4));
    end
    for i=1:3
       disp(sprintf('T_EK_RG: %6.4f %6.4f',T_EK_RG(i),w_4));
    end
    for i=1:3
       disp(sprintf('T_FV_RG: %6.4f %6.4f',T_FV_RG(i),w_4));
    end
    disp(' ');
    for i=1:1
       disp(sprintf('T_K_D: %6.4f %6.4f',T_K_D(i),w_DT));
    end
    for i=3:-1:1
       disp(sprintf('T_FV_H2O: %6.4f %6.4f %6.4f',T_FV_H2O(i),w_8,w_12));
    end
    for i=1:1
       disp(sprintf('T_DE: %6.4f %6.4f',T_DE(i),w_15));
    end
    for i=3:-1:1
       disp(sprintf('T_EK_H20: %6.4f %6.4f',T_EK_H2O(i),w_15));
    end
    for i=1:1
       disp(sprintf('T_FD_H20: %6.4f %6.4f',T_FD_H2O(i),w_15));
    end
    for i=3:-1:1
       disp(sprintf('T_OH_D: %6.4f %6.4f',T_OH_D(i),w_15));
    end
    for i=1:1
       disp(sprintf('T_21: %6.4f %6.4f',T_OH_D(1),w_21));
    end
    for i=1:3
       disp(sprintf('Q_OH: %6.4f %6.4f %6.4f',Q_OH(i)/1E6));
    end
    for i=1:3
       disp(sprintf('Q_FD: %6.4f %6.4f %6.4f',Q_FD(i)/1E6));
    end
    for i=1:3
       disp(sprintf('Q_EK: %6.4f %6.4f %6.4f',Q_EK(i)/1E6));
    end
    for i=1:3
       disp(sprintf('Q_FV: %6.4f %6.4f %6.4f',Q_FV(i)/1E6));
    end
    for i=1:3
       disp(sprintf('Q_K: %6.4f %6.4f %6.4f',Q_K(i)/1E6));
    end
    disp(sprintf('Ws_P1 Ws_P2 Ws_DT: %6.4f %6.4f %6.4f',Ws_P1,Ws_P2,Ws_DT));
    disp('');
    disp(sprintf('P_23: %6.4f',P_23/1e5));
    disp(sprintf('P_DE: %6.4f',P_DE/1e5));
    disp(sprintf('P_DR: %6.4f',P_DR/1e5));
    pause
end
return