% =========================================================================
% Script for solving an fsolve model of a heat exchanger at the Brobekk
% plant in Oslo. The results are hot and cold flow temperatures, and a
% profile of transferred heat.
%
% H. Smedsrud, 2007
%==========================================================================

clc
clear all
close all

% Parameters
Th_in = 180;
Tc_in =62;
FB = 250;
FV = 250;
N = 10;
Q_max = 16;
deltaTlm = 40;
UA = Q_max/deltaTlm;
UA_cell = UA/N;

% Correcting the heat transfer number
UA_cell = UA_cell*0.682731;     % Tc_in = 62; FV = 250
%UA_cell = UA_cell*0.560221;     % Tc_in = 62; FV = 450

par = [Th_in;Tc_in;UA_cell;FB;FV];

% Initial values
x0 = [Th_in-((1/N)*(Th_in-126));
      Th_in-((1.9/N)*(Th_in-126));
      Th_in-((2.8/N)*(Th_in-126));
      Th_in-((3.7/N)*(Th_in-126));
      Th_in-((4.6/N)*(Th_in-126));
      Th_in-((5.5/N)*(Th_in-126));
      Th_in-((6.4/N)*(Th_in-126));
      Th_in-((7.3/N)*(Th_in-126));
      Th_in-((8.2/N)*(Th_in-126));
      126;
      Q_max*0.9/N;
      Q_max*0.8/N;
      Q_max*0.7/N;
      Q_max*0.6/N;
      Q_max*0.5/N;
      Q_max*0.4/N;
      Q_max*0.3/N;
      Q_max*0.2/N;
      Q_max*0.1/N;
      Q_max*0.05/N;
      Tc_in+((1/N)*(Th_in-126));
      Tc_in+((2/N)*(Th_in-126));
      Tc_in+((3/N)*(Th_in-126));
      Tc_in+((4/N)*(Th_in-126));
      Tc_in+((5/N)*(Th_in-126));
      Tc_in+((6/N)*(Th_in-126));
      Tc_in+((7/N)*(Th_in-126));
      Tc_in+((8/N)*(Th_in-126));
      Tc_in+((9/N)*(Th_in-126));
      130];

% Solving
options = optimset('TolFun',1e-10);
[x fval exitflag] = fsolve(@Brobekk_HE_fsolve_fun,x0,[],par);
Th = x(1:10)
Tc = x(21:30)
Q = x(11:20)
Q_tot = sum(x(11:20))
exitflag

% Plotting
n = 1:10;
Tc = [x(30);x(29);x(28);x(27);x(26);x(25);x(24);x(23);x(22);x(21)];

% Make subplot 1
figure(1)
subplot(221)
plot(n,Th,'r',n,Tc,'b')
xlabel('Cell number')
ylabel('Temperature [°C]')
axis([1 10 50 200])
title('A')

subplot(222)
plot(n,Q,'k')
xlabel('Cell number')
ylabel('Transferred heat [MW]')
axis([1 10 1.3 1.9])
title('B')

% Make subplot 2
load 'Th_450';
load 'Tc_450';
figure(2)
subplot(221)
plot(n,Th,'r',n,Tc,'b')
xlabel('Cell number')
ylabel('Temperature [°C]')
axis([1 10 50 200])
title('A')

load 'Q_450';
subplot(222)
plot(n,Q,'k')
xlabel('Cell number')
ylabel('Transferred heat [MW]')
axis([1 10 1.3 1.9])
title('B')