%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%                                              %%%%%
%%%%%                   CASE I                     %%%%%
%%%%%      Two Heat Exchangers 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 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 wh1 wh2 Th1in Th2in UA1 UA2]

caseI    = [60 100  30 50  120 140  50 80];
caseII   = [60 100  30 50  120 140  120 200];
caseIII  = [60 100  15 15  120 140  120 200];
caseIV   = [60 100  50 50  120 140  50 80];
caseV    = [60 100  50 50  120 140  120 200];
caseVI   = [60 100  80 60  120 140  50 80];
caseVII  = [60 100  80 60  120 140  120 200];
caseVIII = [60 100  80 60  120 140  95 110];
HEN_F    = [147.0776 154477.15 91650.28 90371.74 ...
            248.07 248.07 110709 112203];
HEN_H    = [241.1040+0.8019 422096.96 159606.70 ... 
            147970.52 358.5740 138842 209356];

% Select case
casesel = caseI;

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

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

% Design parameters
UA1 	= casesel(7);       % [kW/K]
UA2 	= casesel(8);       % [kW/K]

% Number of iterations
iter=10000;

n = zeros(iter,1);
T1=n; T2=n; Th1=n; Th2=n; Tmix=n; e1=n; eh1=n; e2=n; eh2=n; 
C1=n; C2=n; NTU1=n; NTU2=n; U=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] = tempCalc(T0,w0,UA1,UA2,Th1in,wh1,Th2in,wh2,u);
     
    T1(i)=T(1); T2(i)=T(2); Th1(i)=T(3); Th2(i)=T(4); 
    Tmix(i)=T(5); e1(i)=HE(1); eh1(i)=HE(2); e2(i)=HE(3); 
    eh2(i)=HE(4); C1(i)=HE(5); C2(i)=HE(6); NTU1(i)=HE(7); 
    NTU2(i)=HE(8);
    
end

%% Results

% Self-optimizing variable
c = (T1-T0).^2./(Th1in-T0) - (T2-T0).^2./(Th2in-T0);

% Chen approximation
CHEN =  2.*T1.^2.*T2.^2.*Th2.*Th1in.*Th2in+8.*T1.*T2.^3.*Th2.*Th1.*Th1in-4.*T1.*Th2.*Th2in.*Th1in.^2.*T2.^2-4.*T1.*Th2.^2.*Th1.*Th1in.*T2.^2+4.*T1.*T2.^3.*Th2.*Th1in.^2+8.*T2.^2.*Th1.*Th2.*Th2in.*Th1in.^2-8.*T1.^2.*Th2.*Th1.*Th1in.*Th2in.^2-2.*T1.^2.*Th2.^2.*Th1.*Th2in.*Th1in+4.*T1.^2.*T2.*Th2.*Th1.^2.*Th2in+4.*T1.^2.*T2.*Th1.*Th1in.*Th2in.^2-2.*T1.^2.*T2.^2.*Th1.*Th2in.*Th1in+2.*T1.^3.*Th1.*Th2in.*T2.^2-4.*T1.^3.*T2.*Th1.*Th2in.^2+2.*T1.^3.*Th1.*Th2.^2.*Th2in+8.*T1.^3.*Th1.*Th2.*Th2in.^2-4.*T1.^2.*Th2.*Th1.^2.*Th2in.^2+2.*T1.^2.*T2.*Th1.^2.*Th2in.^2-T1.^2.*T2.^2.*Th1.^2.*Th2in+T1.^2.*T2.^2.*Th2.^2.*Th1in-2.*T1.^2.*Th2.*Th1in.*T2.^3-T1.^2.*Th2.^2.*Th1.^2.*Th2in-2.*T1.^2.*Th1.*Th2in.^3.*Th1in-2.*T1.*Th2.^2.*Th1in.^2.*T2.^2-2.*T2.^3.*Th1.^2.*Th2.*Th1in-8.*T2.^3.*Th1.*Th2.*Th1in.^2+T2.^2.*Th1.^2.*Th2.^2.*Th1in+4.*T2.^2.*Th1.*Th2.^2.*Th1in.^2+2.*T2.^2.*Th2.*Th1in.^3.*Th2in-T1.^2.*Th1.^2.*Th2in.^3+2.*T2.^2.*Th1.^2.*Th2.*Th1in.*Th2in-2.*T2.^3.*Th2.*Th1in.^3+T2.^2.*Th2.^2.*Th1in.^3+2.*T1.^3.*Th1.*Th2in.^3-8.*T1.^3.*T2.*Th1.*Th2.*Th2in+8.*T1.^2.*T2.*Th2.*Th1.*Th2in.*Th1in-8.*T1.*Th2.*Th1.*Th2in.*Th1in.*T2.^2;

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

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

GRAD = UA2.*UA1.*T1.*wh2.*wh1+UA1.*w2.*T1.*UA2.*wh1+UA2.*T0.*w1.*wh2.*UA1-UA1.*T0.*w2.*UA2.*wh1-UA2.*UA1.*T2.*wh2.*wh1-UA2.*T2.*wh2.*w1.*UA1+2.*UA2.*T0.*w1.*wh2.*wh1-2.*UA2.*T2.*wh2.*w1.*wh1-2.*UA1.*T0.*w2.*wh2.*wh1+2.*UA1.*w2.*T1.*wh2.*wh1;

SUM = [U T1 Th1 T2 Th2 Tmix c];
HE = [e1 eh1 e2 eh2 C1 C2 NTU1 NTU2];

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

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

% Information about the heat transfer (optimal split)
he=[e1(nr) eh1(nr) e2(nr) eh2(nr) C1(nr) C2(nr) ...
    NTU1(nr) NTU2(nr)]';

% Finding the self-optimizing split
[Ceq2,nr2]=min(abs(c));

split2=U(nr2);
T1m2=T1(nr2);
Th1m2=Th1(nr2);
T2m2=T2(nr2);
Th2m2=Th2(nr2);
Ceq2;
Tmixm2=Tmix(nr2);

% Information about the heat transfer
he2=[e1(nr2) eh1(nr2) e2(nr2) eh2(nr2) C1(nr2) C2(nr2) ...
    NTU1(nr2) NTU2(nr2)]';

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

% Looking for the 'T1-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);

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

Tr1=dT11./dT12;
Tr2=dT21./dT22;

% Calculating errors
[eAM1 eAM2] = errCalc(dT11, dT12, dT21, dT22);

%% Some plots
close all;

% Tend and controlled variable 
h = figure;

y1start=190; y1slutt=212; sprangy1=2;
y2start=-100; y2slutt=120; sprangy2=20;

[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)
axis(AX(1),[0 1 y1start y1slutt]);
axis(AX(2),[0 1 y2start y2slutt]);
set(AX(1),'YLim',[y1start y1slutt])
set(AX(1),'YTick',y1start:sprangy1:y1slutt)
set(AX(2),'YLim',[y2start y2slutt])
set(AX(2),'YTick',y2start:sprangy2:y2slutt)
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,Tr1,U,Tr2,'linewidth',2);
axis([0 1 0 0.1]);
xlabel('Split, u','fontsize',12);
ylabel('\theta_1/\theta_2','fontsize',12);
legend('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,eAM1,U,eAM2,'linewidth',2); 
xlabel('Split, u','fontsize',12);
ylabel('Error [%]','fontsize',12);
legend('HX1','HX2','fontsize',12);
axis([0 1 0 500]);
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 dT1opt dT2opt] = Tprof(T0,Th1m,Th2m,w1opt,w2opt,...
    wh1,wh2,UA1,UA2);

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

[z dT1soc dT2soc] = Tprof(T0,Th1m2,Th2m2,w1soc,w2soc,...
    wh1,wh2,UA1,UA2);

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

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