%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%                                              %%%%%
%%%%%                   CASE II                    %%%%%
%%%%%        Two Heat Exchangers in Series         %%%%%
%%%%%              and One in Parallel             %%%%%
%%%%%                                              %%%%%
%%%%%          Daniel Greiner Edvardsen            %%%%%
%%%%%                 IKP, NTNU                    %%%%%
%%%%%               June 14, 2011                  %%%%%
%%%%%                                              %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Info: This m-file calculates the hot and cold outlet 
% temperatures for two heat exchangers in series and 
% one in parallel for a given set of conditions. 
% Further, the optimum is found (optimal split) that 
% results in the maximum end temperature. The SOC split 
% is also found and compared with the optimal split. 

clc;
clear all;
close all;

%% Defining parameters

% Cases evaluated. Vector parameters: 
% [T0 w0 wh11 wh1 wh2 Th11in Th1in Th2in UA11 UA1 UA2]
caseI    = [60 100  30 50 20   120 140 140  50 80 65];
caseII   = [60 100  30 50 40   120 140 140  100 180 60];
caseIII  = [60 100  25 50 25   70 160 125   30 70 33];
caseIV   = [60 100  150 120 70 140 120 141  90 150 65];
caseV    = [60 100  30 70 40   110 75 125   50 80 65];
caseVI   = [60 100  30 50 40   120 140 80   50 80 65];

%Select case
casesel = caseI;

%%%%%%%%%%%%%%% Operation parameters %%%%%%%%%%%%%%%%%%%%
T0 	= casesel(1);		% Feed stream temperature [K] 	
w0	= casesel(2);		% m*Cp [kW/K]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%% Utility parameters %%%%%%%%%%%%%%%%%%%%%%
wh11    = casesel(3); % Hot stream 1 Heat Capacity  [kW/K]
wh1 	= casesel(4); % Hot stream 11 Heat Capacity [kW/K]
wh2 	= casesel(5); % Hot stream 2 Heat Capacity  [kW/K]
Th11in	= casesel(6); % Hot stream 1 Temperature    [kW/K]
Th1in	= casesel(7); % Hot stream 11 Temperature   [kW/K]
Th2in 	= casesel(8); % Hot stream 2 Temperature    [kW/K]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%% Design parameters %%%%%%%%%%%%%%%%%%%%
UA11 	= casesel(9);   % [kW/K]
UA1 	= casesel(10);  % [kW/K]
UA2 	= casesel(11);  % [kW/K]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Number of iterations
iter=1000;

n = zeros(iter,1);
T11=n; T2=n; T1=n; Th11=n; Th2=n; Th1=n; Tmix=n; 
e11=n; eh11=n; e2=n; eh2=n; e1=n; eh1=n; C11=n; 
C2=n; C1=n; NTU11=n; NTU2=n; NTU1=n; U=n; CON=n; 
CHEN=n; CHEN2=n; CON2=n;

%% T-out, HE and CONSTRAINT 
for i=1:iter
    
    u = i/iter;
    U(i)=u;
    
    % Calculating outlet temperatures and info about HEs 
    % Only u is changing
    [T HE] = tempCalc2(T0,w0,UA11,UA2,UA1,Th11in,Th2in,...
        Th1in,wh11,wh2,wh1,u);
    
    T11(i)=T(1); T2(i)=T(2); T1(i)=T(3); Th11(i)=T(4); 
    Th2(i)=T(5); Th1(i)=T(6); Tmix(i)=T(7); e11(i)=HE(1);
    eh11(i)=HE(2); e2(i)=HE(3); eh2(i)=HE(4); C11(i)=HE(5);
    C2(i)=HE(6); NTU11(i)=HE(7); NTU2(i)=HE(8);
    
end

%% Results

% Self-optimizing variable
c = (T2-T0).^2./(Th2in-T0) - (((Th1in-T1)./(Th11in-T0)-1).*(T11-T0).^2./(Th1in-T11)+(T1-T0).^2./(Th1in-T11));

%No exergy loss: T1 = T2
EQ = T1-T2;

% Gradient, GRAD
w1=w0.*U;
w2=w0.*(1.-U);

GRAD = -4.*UA2.*T2.*wh1.*w1.^2.*wh11.*wh2-2.*UA2.*T2.*w1.^2.*UA1.*wh11.*wh2+2.*w1.*T1.*wh1.*UA1.*wh11.*w2.*UA2+4.*w1.*T1.*wh1.*UA1.*wh11.*w2.*wh2+2.*T1.*UA1.*wh1.*w1.*wh11.*w2.*UA2+4.*T1.*UA1.*wh1.*w1.*wh11.*w2.*wh2+2.*UA2.*w1.*T1.*wh1.*UA1.*wh11.*wh2+2.*UA2.*T1.*UA1.*wh1.*w1.*wh11.*wh2-2.*UA2.*T2.*w1.*wh1.*UA1.*wh11.*wh2-2.*UA2.*UA11.*T2.*wh1.*w1.^2.*wh2-UA2.*T2.*w1.^2.*UA1.*UA11.*wh2+2.*UA11.*w1.*T1.*wh1.*wh11.*w2.*UA2+4.*UA11.*w1.*T1.*wh1.*wh11.*w2.*wh2+2.*UA2.*UA11.*w1.*T1.*wh1.*wh11.*wh2-2.*UA2.*UA11.*T2.*w1.*wh1.*wh11.*wh2+UA1.*UA11.*T1.*wh1.*wh11.*w2.*UA2+UA11.*w1.*T1.*wh1.*UA1.*w2.*UA2+2.*UA11.*T1.*UA1.*wh1.*w1.*w2.*wh2+UA11.*T1.*UA1.*wh1.*w1.*w2.*UA2+2.*UA11.*w1.*T1.*wh1.*UA1.*w2.*wh2+UA2.*UA11.*w1.*T1.*wh1.*UA1.*wh2+UA2.*UA11.*T1.*UA1.*wh1.*w1.*wh2-UA2.*UA11.*T2.*w1.*wh1.*UA1.*wh2+UA11.*w1.*T1.*UA1.*wh11.*w2.*UA2+UA2.*UA11.*w1.*T1.*UA1.*wh11.*wh2+2.*UA11.*w1.*T1.*UA1.*wh11.*w2.*wh2-UA2.*UA11.*T2.*w1.*UA1.*wh11.*wh2+2.*UA1.*UA11.*T1.*wh1.*wh11.*w2.*wh2+UA1.*UA2.*UA11.*T1.*wh1.*wh11.*wh2-UA1.*UA2.*UA11.*T2.*wh1.*wh11.*wh2;

SUM = [U T11 Th1 T2 Th2 Tmix c];
HE = [e11 eh11 e2 eh2 C11 C2 NTU11 NTU2];

% Finding optimal split
[Tmixm,nr]=max(Tmix);

split=U(nr);
T11m=T11(nr); T2m=T2(nr); T1m=T1(nr); Th11m=Th1(nr); 
Th2m=Th2(nr); Th1m=Th1(nr); Ceq=c(nr); Tmixm;

Tm=[T11m T2m T1m Th11m Th2m Th1m Tmixm]';

% Information about the heat transfer (optimal)
he=[e11(nr) eh11(nr) e2(nr) eh2(nr) C11(nr) C2(nr) 
    NTU11(nr) NTU2(nr)]';

% 'theta1'/'theta2 (in optimum)'
Tr=[(Th11m-T0)/(Th11in-T11m) 
    (Th2m-T0)/(Th2in-T2m)
    (Th1m-T11m)/(Th1in-T1m)];

% Calculating the temperature differences 
% on each side of the HEs
dT111=Th1-T0; dT112=Th11in-T11; dT21=Th2-T0; 
dT22=Th2in-T2; dT11=Th1-T11; dT12=Th1in-T1;

Tr11=dT111./dT112;
Tr2=dT21./dT22;
Tr1=dT11./dT12;

% Calculating errors
[eAM11 eAM2 eAM1] = errCalc2(dT111, dT112, dT21, dT22, ...
    dT11, dT12);

int=round(split*iter);
eTOT=sqrt(eAM11(int)^2+eAM2(int)^2+eAM1(int)^2);

% Finding the self-optimizing split
[Ceq2,nr2]=min(abs(c));
split2=U(nr2);
T11m2=T11(nr2); Th11m2=Th1(nr2); T2m2=T2(nr2); 
Th2m2=Th2(nr2); T1m2=T1(nr2); Th1m2=Th1(nr2); 
Ceq2; Tmixm2=Tmix(nr2);

Tm2=[T11m2 T2m2 T1m2 Th11m2 Th2m2 Th1m2 Tmixm2]';

% Information about the heat transfer (SCOC)
he2=[e11(nr2) eh11(nr2) e2(nr2) eh2(nr2) C11(nr2) C2(nr2) 
    NTU11(nr2) NTU2(nr2)]'

% Looking for the 'Gradient split'
[Ceq3,nr3]=min(abs(GRAD));
split3=U(nr3);

% Looking for the 'T11-T2 split'
[Ceq4,nr4]=min(abs(EQ));
split4=U(nr4);

% Looking for the 'CHEN split'
[Ceq5,nr5]=min(abs(CHEN));
split5=U(nr5);

% Deviation from optimal split
best = abs([split2,split3,split4,split5]'-split);

%% Plots
close all;

% Tend and controlled variable 
h = figure;
[AX,H1,H2] = plotyy(U, Tmix, U,c); 
set(get(AX(1),'Ylabel'),'String','T_{end} [ \circC]',...
    'fontsize',12) 
set(get(AX(2),'Ylabel'),'String',...
    'Controlled variable, c [ \circC]','fontsize',12)
set(H1,'linewidth',2)
set(H2,'linewidth',2)
xlabel('Split, u','fontsize',12);
hold on; 
gridxy(split2,'Color','r','LineStyle','--','LineWidth',2); 
hold on;
grid on;
print(h,'-depsc','temp.eps');

%Error plots
h = figure;
plot(U,Tr11,U,Tr1,U,Tr2,'linewidth',2);
axis([0 1 0 10]);
xlabel('Split, u','fontsize',12);
ylabel('\theta_1/\theta_2','fontsize',12);
legend('HX11','HX1','HX2','fontsize',12);
hline([1/1.4 1.4 1],{'magenta','magenta','black:'},...
    {'LB','UB',''})
hold on; 
gridxy(split,'Color','r','LineStyle','--','LineWidth',2);
print(h,'-depsc','approx.eps');

h = figure;
plot(U,eAM11,U,eAM1,U,eAM2,'linewidth',2); 
xlabel('Split, u','fontsize',12);
ylabel('Error [%]','fontsize',12);
legend('HX11','HX1','HX2','fontsize',12);
axis([0 1 0 30]);
hold on; 
gridxy(split,'Color','r','LineStyle','--','LineWidth',2);
print(h,'-depsc','error.eps');

% Temperature profiles

%Optimal profiles
w1opt=split*w0;
w2opt=(1-split)*w0;

[z dT11opt dT2opt dT1opt] = Tprof2(T0,T11m,Th11m,Th2m,...
    Th1m,w1opt,w2opt,wh11,wh2,wh1,UA11,UA2,UA1,iter);

%SOC profiles
w1soc=split2*w0;
w2soc=(1-split2)*w0;

[z dT11soc dT2soc dT1soc] = Tprof2(T0,T11m2,Th11m2,...
    Th2m2,Th1m2,w1soc,w2soc,wh11,wh2,wh1,UA11,...
    UA2,UA1,iter);

% Optimal
h = figure;
plot(z,dT11opt,z,dT1opt,z,dT2opt,'linewidth',2);
xlabel('Length','fontsize',12);
ylabel('\DeltaT [ \circC]','fontsize',12);
legend('HX11','HX1','HX2','fontsize',12);
grid on;
print(h,'-depsc','optprof.eps');

% SOC
h = figure;
plot(z,dT11soc,z,dT1soc,z,dT2soc,'linewidth',2);
xlabel('Length','fontsize',12);
ylabel('\DeltaT [ \circC]','fontsize',12);
legend('HX11','HX1','HX2','fontsize',12);
grid on;
print(h,'-depsc','socprof.eps');
