% Data passed from / to TRNSYS 
% ---------------------------- 
%
% trnTime (1x1)        : simulation time 
% trnInfo (15x1)       : TRNSYS info array
% trnInputs (nIx1)     : TRNSYS inputs 
% trnStartTime (1x1)   : TRNSYS Simulation Start time
% trnStopTime (1x1)    : TRNSYS Simulation Stop time
% trnTimeStep (1x1)    : TRNSYS Simulation time step
% mFileErrorCode (1x1) : Error code for this m-file. It is set to 1 by TRNSYS and the m-file should set it to 0 at the
%                        end to indicate that the call was successful. Any non-zero value will stop the simulation
% trnOutputs (nOx1)    : TRNSYS outputs  

% trnInputs
% ---------
%
% trnInputs(1) : ta_inlet
% trnInputs(2) : tf_inlet
% trnInputs(3) : tw_inlet
% trnInputs(4) : RHa_inlet
% trnInputs(5) : RHf_inlet
%
% trnOutputs
% ----------
%
% trnOutputs(1) : ta_outlet
% trnOutputs(2) : tf_outlet
% trnOutputs(3) : tw_outlet
% trnOutputs(4) : RHa_outlet
% trnOutputs(5) : RHf_outlet


mFileErrorCode = 100    % Beginning of the m-file 

%--------------------------------------------------------------------------
%----PARAMETERS------------------------------------------------------------ 
%--------------------------------------------------------------------------
%dimensies warmtewisselaar
L=2;                %lengte warmtewisselaar in m
W=1;                %breedte warmtewisselaar in m
H=1;                %hoogte warmtwisselaar in m
b_a=0.0045;         %breedte tussen de platen in m
b_f=0.0045;         %breedte tussen de platen in m

debiet=7100;        %totaal luchtdebiet in m³/h
ro=1.2;             %densiteit van lucht in kg/m³

%definieer constanten
zeta=6.1*10^(-6);   %dampcapaciteit in kg/kg/Pa
N=100;              %aantal knopen waarin temperaturen worden berekend
deltax=W*L/N;       %lengte interval in m*breedte in m
Ca=1008;            %specifieke warmte extractielucht J/kgK
Cf=1008;            %specifieke warmte pulsielucht J/kgK
hev=2500000;        %verdampingswarmte in J/kg

k_air=0.025;        %warmtegeleiding van lucht in W/mK
Pr=0.707;           %Prandl getal
nu=1.5*10^(-5);     %kinematische viscositeit in m²/s
d_p=0.0001;         %dikte van de platen in m    
k_p=0.040;          %warmtegeleidingscoefficient van de wisselaar in W/mK (polypropyleen)
d_w=0.00005;        %dikte waterfilm in m
k_w=0.606;          %warmtegeleidingscoefficient water in W/mK < uit paper

%berekening overgangscoefficienten en debieten per laag
aantal=H/(b_a+b_f+d_p+d_w);                 %aantal laagjes resp. in de twee stromen
Q=debiet/aantal;                            %debiet per laag in m³/h
Ga=(Q/2)*ro/3600;                           %half debiet per kanaal in kg/s
Gf=(Q/2)*ro/3600;    %orde 2.5*10-5         %half debiet per kanaal in kg/s

D_a=2*W*b_a/(W+b_a);                        %hydraulische diameter van 1 laagje: [4*sectie/natte omtrek] in m
v_a=Q/(3600*b_a*W);   % orde 2m/s           %luchtsnelheid per luchtlaag in m/s
Re_a=v_a*D_a/nu; % orde 600                 %Re berekend ahv hydraulische diam en debiet per laagje

D_f=2*W*b_f/(W+b_f);
v_f=Q/(3600*b_f*W);
Re_f=v_f*D_f/nu;

%Nusselt ook berekenen op basis van hydraulische diameter
if   Re_a<2300;
     Nu_a=((3.66^3)+(1.61^3)*(Re_a*Pr*D_a/L))^(1/3);      %laminair      
else zeta_a=(1.82*log(Re_a)-1.64)^(-2);
     Nu_a=((zeta_a/8)*(Re_a-1000)*Pr)/(1+12.7*((zeta_a/8)^(0.5))*((Pr^(2/3))-1));    %turbulent
end

if   Re_f <2300;
     Nu_f=((3.66^3)+(1.61^3)*(Re_f*Pr*D_f/L))^(1/3);    
else zeta_f=(1.82*log(Re_f)-1.64)^(-2);
     Nu_f=((zeta_f/8)*(Re_f-1000)*Pr)/(1+12.7*((zeta_f/8)^(0.5))*((Pr^(2/3))-1));  
end

%boek Heat Transfer ; A. BEJAN. p. 302 eq. (6.43)
%Nu (D) = h*D/k met D daarin de hydraulische diameter

ha=Nu_a*k_air/D_a;
hf=Nu_f*k_air/D_f;
%ha=Nu_a*k_air/(b_a/2);  %warmteovergangscoefficient aan zijde extractielucht
%hf=Nu_f*k_air/(b_f/2);  %warmteovergangscoefficient aan zijde pulsielucht
alfa_a=ha;              %W/m²K
alfa_f=1/((d_p/k_p)+(d_w/k_w)+(1/hf));  %+weerstand warmtewisselaar en waterfilm W/m²K

%berekening dampovergangscoefficient via Lewis - relatie
Pair=101325;            %luchtdruk Pa
beta=alfa_a/(Ca*Pair);  %in kg/s/Pa/m²      %orde 2*10-7

%iteraties
a=0.5;        %relaxatie-factor


mFileErrorCode = 110    % After setting parameters

%--------------------------------------------------------------------------
% --- Process Inputs-------------------------------------------------------
% -------------------------------------------------------------------------
nI = trnInfo(3);
nO = trnInfo(6);

ta_inlet = trnInputs(1);%inlaat temperatuur extractielucht in °C
tf_inlet = trnInputs(2);%inlaat temperatuur pulsielucht in °C
tw_inlet = trnInputs(3); %inlaat temperatuur waterfilm in °C
RHa_inlet = trnInputs(4);%inlaat relatieve vochtigheid extractielucht in %
RHf_inlet = trnInputs(5);%inlaat relatieve vochtigheid pulsielucht in %

mFileErrorCode = 120    % After processing inputs

%--------------------------------------------------------------------------
% --- First call of the simulation: initial time step (no iterations) -----
% -------------------------------------------------------------------------
% (note that Matlab is initialized before this at the info(7) = -1 call, but the m-file is not called)

if ( (trnInfo(7) == 0) & (trnTime-trnStartTime < 1e-6) )  
    % This is the first call (Counter will be incremented later for this very first call)
    iCall = 0;

    % This is the first time step
    iStep = 1;

    mFileErrorCode = 130    % After initialization
end

%--------------------------------------------------------------------------
% --- All iterative calls -------------------------------------------------
% -------------------------------------------------------------------------

% --- If this is a first call in the time step, increment counter ---

if ( trnInfo(7) == 0 )
    iStep = iStep+1;
end
mFileErrorCode = 150;   % After reading inputs 

%--------------------------------------------------------------------------
% --- Calculations --------------------------------------------------------
%--------------------------------------------------------------------------
%EERSTE RUN MET BEGINWAARDEN    

pa_inlet=(RHa_inlet/100)*exp(65.8094-(7066.27/(273.15+ta_inlet))-5.976*log(273.15+ta_inlet));
pf_inlet=(RHf_inlet/100)*exp(65.8094-(7066.27/(273.15+tf_inlet))-5.976*log(273.15+tf_inlet));

M=zeros(3*N,3*N);   %matrix met coefficienten 
T=zeros(3*N,1);     %matrix met onbekende temperaturen
D=zeros(3*N,1);     %matrix met constanten

%coefficienten matrix M
aa=zeros(N,1);
ba=zeros(N,1);
ca=zeros(N,1);
aw=zeros(N,1);
bw=zeros(N,1);
cw=zeros(N,1);
af=zeros(N,1);
bf=zeros(N,1);
cf=zeros(N,1);
%coefficienten matrix D
da=zeros(N,1);
dw=zeros(N,1);
df=zeros(N,1);
%onbekende temperaturen in matrix T
ta=zeros(N,2);  %3kolommen; eerste kolom oude waarde, tweede kolom nieuwe
tw=zeros(N,2);  %3kolommen; eerste kolom oude waarde, tweede kolom nieuwe
tf=zeros(N,2);  %3kolommen; eerste kolom oude waarde, tweede kolom nieuwe
%onbekende dampdrukken 
pa=zeros(N,1);
pw=zeros(N,1);
%startwaarden
da(1,1)=-Ga*Ca*ta_inlet;  %eerste coefficient matrix D  
df(N,1)=-Gf*Cf*tf_inlet;  %laatste coefficient matrix D
tw(1,1)=tw_inlet;
ta(1,1)=ta_inlet;
tf(N,1)=tf_inlet;
pw(1,1)=exp(65.8094-(7066.27/(273.15+tw(1,1)))-5.976*log(273.15+tw(1,1)));
pa(1,1)=(Ga*zeta*pa_inlet+beta*deltax*pw(1,1))/(Ga*zeta+beta*deltax); 


%VOOR ALLE CONTROLEVOLUMES !!!!

for i=1:N
%schatten temperatuur = inlaatcondities in elke knoop van de warmtewisselaar
ta(i,1)=ta(1,1);
tw(i,1)=tw(1,1);
tf(i,1)=tf(N,1); 

%berekening dampdruk wateropp uitgaande van tw_inlet
pw(i,1)=exp(65.8094-(7066.27/(273.15+tw(i,1)))-5.976*log(273.15+tw(i,1))); 
end

%berekening dampdruk extractielucht
for i=2:N
pa(i,1)=(Ga*zeta*pa(i-1,1)+beta*deltax*pw(i,1))/(Ga*zeta+beta*deltax);
end

for i=1:N
%%WARMTESYSTEEM
%berekening coefficienten matrix M (veronderstelling deltax constant)
aa(i,1)=Ga*Ca;
ba(i,1)=-(Ga*Ca+deltax*alfa_a);
ca(i,1)=deltax*alfa_a;
aw(i,1)=deltax*alfa_a;
bw(i,1)=-deltax*(alfa_a+alfa_f);
cw(i,1)=deltax*alfa_f;
af(i,1)=deltax*alfa_f;
bf(i,1)=-(Gf*Cf+deltax*alfa_f);
cf(i,1)=Gf*Cf;
dw(i,1)=-hev*beta*deltax*(pa(i,1)-pw(i,1));       

%plaatsing elementen in matrix M
        M(3*i-2,3*i-2)=ba(i,1);
        M(3*i-2,3*i-1)=ca(i,1);
        M(3*i-1,3*i-2)=aw(i,1);
        M(3*i-1,3*i-1)=bw(i,1);
        M(3*i-1,3*i)=cw(i,1);
        M(3*i,3*i-1)=af(i,1);
        M(3*i,3*i)=bf(i,1);
%plaatsing constanten in matrix D
        D(3*i-1,1)=dw(i,1);
end

D(1,1)=da(1,1);
D(3*N,1)=df(N,1);

%aa(1) en cf(N) moet hij niet plaatsen!
for i=2:N
    M(3*i-2,3*i-5)=aa(i,1);           
end         
for i=1:N-1
    M(3*i,3*i+3)=cf(i,1);
end

%M*T=D
A=inv(M);
T=A*D;

for i=1:N
    %uitlezen matrix T
    ta(i,2)=T(3*i-2,1);
    tw(i,2)=T(3*i-1,1);
    tf(i,2)=T(3*i,1);
end

taout=ta(N,2); %ter controle: waarden na eerste run
twout=tw(N,2);
tfout=tf(1,2);  

%--------------------------------------------------------------------------
%ITEREREN SYSTEEM TOT CONVERGENTIEVOORWAARDE VOLDAAN
%--------------------------------------------------------------------------
q=0;
%eerste kolom = oude waarde
while   ((abs(ta(N,2)-ta(N,1)))> 0.01 | (abs(tf(1,2)-tf(1,1)))> 0.01 | (abs(tw(N,2)-tw(N,1)))> 0.01 | (q<=50))
        %((abs(ta(N,2)-ta(N,1))> 0.1)| (q<=10) ) %oude - nieuwe waarde => convergentievoorwaarde !!
        %schuif "oude nieuwe" waarde naar "oude" waarde: shift kolom 2 ->1

        q=q+1 

        for k=1:N
        ta(k,1)=ta(k,1)+a*(ta(k,2)-ta(k,1));
        tw(k,1)=tw(k,1)+a*(tw(k,2)-tw(k,1));
        tf(k,1)=tf(k,1)+a*(tf(k,2)-tf(k,1));
        
        pw(k,1)=exp(65.8094-(7066.27/(273.15+tw(k,1)))-5.976*log(273.15+tw(k,1))); 
        end

        pa(1,1)=(Ga*zeta*pa_inlet+beta*deltax*pw(1,1))/(Ga*zeta+beta*deltax);   
        for k=2:N
        pa(k,1)=(Ga*zeta*pa(k-1,1)+beta*deltax*pw(k,1))/(Ga*zeta+beta*deltax);
        end

D=zeros(3*N,1);
M=zeros(3*N,3*N);
T=zeros(3*N,1);

        for k=1:N     
        %%WARMTESYSTEEM
        %berekening coefficienten
        aa(k,1)=Ga*Ca;
        ba(k,1)=-(Ga*Ca+deltax*alfa_a);
        ca(k,1)=deltax*alfa_a;
        aw(k,1)=deltax*alfa_a;
        bw(k,1)=-deltax*(alfa_a+alfa_f);
        cw(k,1)=deltax*alfa_f;
        af(k,1)=deltax*alfa_f;
        bf(k,1)=-(Gf*Cf+deltax*alfa_f);
        cf(k,1)=Gf*Cf;
        dw(k,1)=-hev*beta*deltax*(pa(k,1)-pw(k,1));       
        
            %plaatsing elementen in matrix M
            M(3*k-2,3*k-2)=ba(k,1);
            M(3*k-2,3*k-1)=ca(k,1);
            M(3*k-1,3*k-2)=aw(k,1);
            M(3*k-1,3*k-1)=bw(k,1);
            M(3*k-1,3*k)=cw(k,1);
            M(3*k,3*k-1)=af(k,1);
            M(3*k,3*k)=bf(k,1);
   
        %plaatsing constanten in matrix D
        D(3*k-1,1)=dw(k,1);
        end
        
        D(1,1)=da(1,1);     
        D(3*N,1)=df(N,1);

        %aa(1) en cf(N) moet hij niet plaatsen! 
        for k=2:N
        M(3*k-2,3*k-5)=aa(k,1);           
        end         
        for k=1:N-1
        M(3*k,3*k+3)=cf(k,1);
        end    
    
        %M*T=D
        B=inv(M);
        T=B*D;

        %uitlezen matrix T :: nieuwe waarden
        for k=1:N
        ta(k,2)=T(3*k-2,1);
        tw(k,2)=T(3*k-1,1);
        tf(k,2)=T(3*k,1);
        end 
        
    ta(N,2);
    tw(N,2);
    tf(1,2);
end

%outputs
ta_outlet=ta(N,2);
tf_outlet=tf(1,2);
tw_oulet=tw(N,2);

pa_outlet=pa(N,1);

%berekening RH outputs
pa_sat=exp(65.8094-(7066.27/(273.15+ta_outlet))-5.976*log(273.15+ta_outlet));
RHa_outlet=100*pa_outlet/pa_sat
pf_sat=exp(65.8094-(7066.27/(273.15+tf_outlet))-5.976*log(273.15+tf_outlet));
RHf_outlet=100*pf_inlet/pf_sat


%--------------------------------------------------------------------------
% --- Set outputs----------------------------------------------------------
%--------------------------------------------------------------------------
trnOutputs(1) = ta_outlet;
trnOutputs(2) = tf_outlet;
trnOutputs(3) = tw_outlet;
trnOutputs(4) = RHaf_outlet;
trnOutputs(5) = RHf_outlet;


mFileErrorCode = 0; % Tell TRNSYS that we reached the end of the m-file without errors
return