% Kjørefil
%--------------------------------------------------------------------------
% Dette er del av hovedoppgaven til Jørgen Bauck Jensen våren 2003. 
% Denne filen definerer nødvendige data for å løse det dynamiske problemet
% gitt i filen model.m. Løseren som brukes er ode15s, og visualisering av 
% resultater gjøres med filen Plotter.m. I tillegg benyttes filen utfall.m
% for å fortelle løseren når endringer inntreffer.
%--------------------------------------------------------------------------
clear all
global u P w Tventil y
start = cputime;

% Generelle parametre 
%--------------------------------------------------------------------------
d.XXXXXX_Generelle_Parametre_XXXXXX = [];
h=0;
d.R = 8.314;
d.n = 20;
d.T0 = 298;
d.P0 = 1.013e5;
d.P0_sat_NG1A = 8.7e5;
d.T0_sat_NG1A = -24+273;
d.DHvapm_sat_NG1A = 15e3;
d.P0_sat_NG1B = 3.3e5;
d.T0_sat_NG1B = -55+273;
d.DHvapm_sat_NG1B = 15e3;

% Molvekter $kg/mol$
%--------------------------------------------------------------------------
d.XXXXXXXXXXXX_Molvekt_XXXXXXXXXXXX = [];
d.MW_NG1A = [34,17,34,29,21]/1000;
d.MW_NG1B = d.MW_NG1A;
d.MW_SW1 = [34,18]/1000;
d.MW_SW2 = d.MW_SW1;
d.MW_tank = [34,34,34]/1000;

% Cp verdier 
%--------------------------------------------------------------------------
% molar $J/(K mol)$
d.XXXXXXX_Molare_Cp_verdier_XXXXXXX = [];
d.Cpgm_NG1A = [65,66,65,68,53];
d.Cpgm_NG1B = d.Cpgm_NG1A;
d.Cpgm_SW1 = [65,0];
d.Cpgm_SW2 = d.Cpgm_SW1;
d.Cpgm_tank = [65,65,65];
d.Cplm_NG1A = [94,0,94,0,0];
d.Cplm_NG1B = d.Cplm_NG1A;
d.Cplm_SW1 = [0,76];
d.Cplm_SW2 = [94,76];
d.Cplm_tank = [0,0,94];
% spesifikk $J/(K kg)$
d.XXXXX_Spesifikke_Cp_verdier_XXXXX = [];
d.Cpg_NG1A = d.Cpgm_NG1A./d.MW_NG1A;
d.Cpg_NG1B = d.Cpgm_NG1B./d.MW_NG1B;
d.Cpg_SW1 = d.Cpgm_SW1./d.MW_SW1;
d.Cpg_SW2 = d.Cpgm_SW2./d.MW_SW2;
d.Cpg_tank = d.Cpgm_tank./d.MW_tank;
d.Cpl_NG1A = d.Cplm_NG1A./d.MW_NG1A;
d.Cpl_NG1B = d.Cplm_NG1B./d.MW_NG1B;
d.Cpl_SW1 = d.Cplm_SW1./d.MW_SW1;
d.Cpl_SW2 = d.Cplm_SW2./d.MW_SW2;
d.Cpl_tank = d.Cplm_tank./d.MW_tank;
d.Cps_w = [900,900,900,900];

% Tetthet $kg/m^3$
%--------------------------------------------------------------------------
d.XXXXXXX_Tetthet_for_vaeske_XXXXXX = [];
d.rhol_NG1A = [487,0,487,0,0];
d.rhol_NG1B = d.rhol_NG1A;
d.rhol_SW1 = [0,1000];
d.rhol_SW2 = [487,1000];
d.rhol_tank = [0,0,487];

% Fordampingsvarme
%--------------------------------------------------------------------------
% molar $J/mol$
d.XXXXX_Molar_fordampingsvarme_XXXX = [];
d.DHvapm_NG1A = [17000,9000,17000,16000,12000];
d.DHvapm_NG1B = d.DHvapm_NG1A;
d.DHvapm_SW1 = [17000,44000];
d.DHvapm_SW2 = d.DHvapm_SW1;
d.DHvapm_tank = [17000,17000,17000];
% spesifikk $J/kg$
d.XXX_Spesifikk_fordampingsvarme_XX = [];
d.DHvap_NG1A = d.DHvapm_NG1A./d.MW_NG1A;
d.DHvap_NG1B = d.DHvapm_NG1B./d.MW_NG1B;
d.DHvap_SW1 = d.DHvapm_SW1./d.MW_SW1;
d.DHvap_SW2 = d.DHvapm_SW2./d.MW_SW2;
d.DHvap_tank = d.DHvapm_tank./d.MW_tank;

%--------------------------------------------------------------------------
% Anleggsdata
%--------------------------------------------------------------------------

% Trykk Pa
%--------------------------------------------------------------------------
d.XXXXXXXX_Stasjonaere_trykk_XXXXXX = [];
d.P_NG1A = [8.7,62,20.7,20,38]*1e5;
d.P_NG1B = [3.3,62,20.7,20,38]*1e5;
d.P_SW1 = [21,1]*1e5;
d.P_SW2 = [20.8,1]*1e5;
d.P_tank = [2.6,7.3,20.7]*1e5;
P = [];

% Stasjonære temperaturer K
%--------------------------------------------------------------------------
d.XXXX_Stasjonaere_temperaturer_XXX = [];
load Startvektor
d.T_inn_NG1A = [T(1),12+273,T(5*d.n+12),12+273,12+273];
d.T_ut_NG1A = [-24+273,T(2*d.n+1),T(3*d.n+1),T(4*d.n+1),T(5*d.n+1)];
d.T_inn_NG1B = [T(5*d.n+13),d.T_ut_NG1A(2:5)];
d.T_ut_NG1B = [-55+273,T(7*d.n+13),T(8*d.n+13),T(9*d.n+13),T(10*d.n+13)];
d.T_inn_SW1 = [76,6]+273;
d.T_inn_SW2 = [T(5*d.n+2),6+273];
d.T_tank = [T(5*d.n+10),T(5*d.n+11),T(5*d.n+12)];
d.T_ut_SW1 = [T(5*d.n+2),T(5*d.n+3)];
d.T_ut_SW2 = [T(5*d.n+4),T(5*d.n+5)];
Tventil=[0,0];

% Volumer $m^3$
%--------------------------------------------------------------------------
d.XXXXXXXXXX_Totale_volum_XXXXXXXXX = [];
d.V_NG1A = [14,7,7,5,5];
d.V_NG1B = [14,10,6,4,4];
d.V_SW1 = [11.2,5.2];
d.V_SW2 = [113.7,68.8];
d.V_tank = [82.4,95.3,207.3];

% Varmeoverføringskoeffesienter $J/K$
%--------------------------------------------------------------------------
d.X_Varmeoverforingskoeffesienter_X = [];
d.hA_NG1A = [42000000,15000000,20000000,9000000,9000000];
d.hA_NG1B = [30000000,15000000,20000000,9000000,9000000];
d.UA_SW1 = 370000;
d.UA_SW2 = 7800000;

% Stasjonært holdup kg 
%--------------------------------------------------------------------------
d.XXXXXXX_Stasjonaere_holdup_XXXXXX = [];
d.ml_NG1A = [0.9,0,1,0,0].*d.V_NG1A.*d.rhol_NG1A;
d.mg_NG1A = [0.1,1,0,1,1].*d.P_NG1A.*d.V_NG1A.*d.MW_NG1A./(d.R.*d.T_inn_NG1A);
d.ml_NG1B = [0.9,0,1,0.5,0.5].*d.V_NG1B.*d.rhol_NG1B;
d.mg_NG1B = [0.1,1,0,0.5,0.5].*d.P_NG1B.*d.V_NG1B.*d.MW_NG1B./(d.R*d.T_inn_NG1B);
d.ml_SW1 = [0,1].*d.V_SW1.*d.rhol_SW1;
d.mg_SW1 = [1,0].*d.P_SW1.*d.V_SW1.*d.MW_SW1./(d.R*d.T_ut_SW1);
d.ml_SW2 = [0.9,1].*d.V_SW2.*d.rhol_SW2;
d.mg_SW2 = [0.1,0].*d.P_SW2.*d.V_SW2.*d.MW_SW2./(d.R*d.T_inn_SW2);
d.ml_tank = [0,0,0.44].*d.V_tank.*d.rhol_tank;
d.mg_tank = [1,1,0.56].*d.P_tank.*d.V_tank.*d.MW_tank./(d.R*d.T_tank);
d.m_w = [90000,70000,4200,63200];

% Stasjonære strømmer $kg/s$
%--------------------------------------------------------------------------
d.XXXXXX_Stasjonaere_strommer_XXXXX = [];
d.w_NG1A = [350000,627101,623000,434260,463642]/3600;
d.w_NG1B = [273000,627101,273000,434260,463642]/3600;
d.w_SW1 = [623000,1000000]/3600;
d.w_SW2 = [623000,15000000]/3600;
d.w_tank = [273000,623000,623000]/3600;
w = [];

% Regulatorparametre
%--------------------------------------------------------------------------
u.n=1; % Normalisert kompressorturtall
u.z = [0.5,0.5]; % Ventilåpninger
y.s = [d.P_tank(1),d.T_ut_NG1A(2),d.T_ut_NG1B(2)]; % Setpunkt

% Kompressor- og ventilparametre
%--------------------------------------------------------------------------
d.XXXXXX_Kompressorparametre_XXXXXX = [];
d.eff = [0.75,0.75]; % Effektivitetsfaktor
d.Q_komp = d.w_tank(1:2).*d.T_tank(1:2).*d.R./(d.P_tank(1:2).*d.MW_tank(1:2)); % Volumstrøm $m^3/s$
d.Cv = [d.w_NG1A(1)/(u.z(1)*sqrt(d.P_NG1A(3)-d.P_NG1A(1))),d.w_NG1B(1)/(u.z(2)*sqrt(d.P_NG1B(3)-d.P_NG1B(1)))];

% Strømningsparametre $m\cdot s \cdot \sqrt{K}$
%--------------------------------------------------------------------------
d.XXXXXXXX_Ventilparametre_XXXXXXXX = [];
d.k_flow_SW1 = d.w_SW1(1)/sqrt((d.P_SW1(1)-d.P_SW2(1))*(d.P_SW1(1)/d.T_ut_SW1(1)+d.P_SW2(1)/d.T_inn_SW2(1)));
d.k_flow_SW2 = d.w_SW2(1)/sqrt((d.P_SW2(1)-d.P_tank(3)));
d.k_flow_NG1A = d.w_NG1A(1)/sqrt((d.P_NG1A(1)-d.P_tank(2))*(d.P_NG1A(1)/d.T_inn_NG1A(1)-d.P_tank(2)/d.T_tank(2)));

d.k_flow_NG1B = d.w_NG1B(1)/sqrt((d.P_NG1B(1)-d.P_tank(1))*(d.P_NG1B(1)/d.T_inn_NG1B(1)-d.P_tank(1)/d.T_tank(1)));

%--------------------------------------------------------------------------
% Integrerer det dynamsike problemet i model.m med ode15s
%--------------------------------------------------------------------------

% Initialverdier
%--------------------------------------------------------------------------
load Startvektor
T0=T';
I = [0;0;0]; % For regulatorene
T0 = [T0;I];
d
tspan = [0 12000];
options = odeset('event',@utfall);
[t,T] = ode15s(@model,tspan,T0,options,d);
Tidsforbruk = cputime-start
y.s = [d.P_tank(1),d.T_ut_NG1A(2),d.T_ut_NG1B(2)]; % Setpunkt
P=[];w=[];u.z=u.z(1,:);u.n=u.n(1,:);Tventil=[];
last = length(t);
for i=1:last
    model(t(i),T(i,:)',d);
end

% Plotter resultatene
%--------------------------------------------------------------------------
Plotter(t/60,T,d);