% This program uses implicit centered-difference technique to solve the 
% unsteady-state countercurrent heat exchanger problem in section 8.3.
%

global Pt Ps Vt Vs;

% simulation parameters
%
Pt = 0.2; Ps = 0.6;	% no values specified in the book
Vt = 2.5; Vs = 8;	% these numbers are assumed
L = 10;			% length of the heat exchanger
N = 20;                 % number of grid points
delta_x = L/N;		% divide the distance grid into N elements.
delta_t = delta_x/Vt;
time = 20;		% simulation time in minutes
iteration = time/delta_t;

% solve for the initial steady state temperature profile. In order to
% get temp. values on the N distance grid points, integration is 
% divided into N small steps.
% Note that the initial steady state problem is a split boundary value
% problem and superposition technique is used.
% 
% generate one particular solution:
%
x0 = 0.0; xf = L;
Tp0 = [45 50]';
for i = 2:1:N+1,
    xf = x0 + delta_x;
    [xp, Tp] = ode45('Model83', x0, xf, Tp0, 1.e-8, 1);
    Tpt(i) = Tp(length(Tp),1);
    Tps(i) = Tp(length(Tp),2);
    x(i) = xp(length(xp));
    Tp0 = [Tpt(i) Tps(i)]';
    x0 = xf;
end
x(1) = 0;
Tpt(1) = 45;
Tps(1) = 50;
%plot(x, Tpt, '-', x, Tps, '-.');
%title('Particular Solution');
%pause;

% generate two homogeneous solutions. In  this special case, the
% differential equations used for the particular solution and 
% homogeneous solutions are the same, which is defined in Model83.m.
%
x0 = 0.0; xf = L;
Th10 = [45 0]';
for i = 2:1:N+1,
    xf = x0 + delta_x;
    [xh1, Th1] = ode45('Model83', x0, xf, Th10, 1.e-8, 1);
    Th1t(i) = Th1(length(Th1),1);
    Th1s(i) = Th1(length(Th1),2);
    x(i) = xh1(length(xh1));
    Th10 = [Th1t(i) Th1s(i)]';
    x0 = xf;
end
x(1) = 0;
Th1t(1) = 45;
Th1s(1) = 0;
%plot(x, Th1t, '-', x, Th1s, '-.');
%title('Homogeneous Solution#1');
%pause;

x0 = 0.0; xf = L;
Th20 = [0 50]';
for i = 2:1:N+1,
    xf = x0 + delta_x;
    [xh2, Th2] = ode45('Model83', x0, xf, Th20, 1.e-8, 1);
    Th2t(i) = Th2(length(Th2),1);
    Th2s(i) = Th2(length(Th2),2);
    x(i) = xh2(length(xh2));
    Th20 = [Th2t(i) Th2s(i)]';
    x0 = xf;
end
x(1) = 0;
Th2t(1) = 0;
Th2s(1) = 50;
%plot(x, Th2t, '-', x, Th2s, '-.');
%title('Homogeneous Solution#2');
%pause;

% solve linear equation set to get the constants that satisfy the real
% boundary conditions.
%
A = [Th1t(1) Th2t(1);
     Th1s(N+1)  Th2s(N+1)];

B = [45-Tpt(1) 75-Tps(N+1)]';

const = A\B;

% compose the final solution.
%
Tt(1,:) = Tpt + const(1)*Th1t + const(2)*Th2t;
Ts(1,:) = Tps + const(1)*Th1s + const(2)*Th2s;

plot(x, Tt(1,:));
title('Tube side temperature');
xlabel('length, ft');
ylabel('temperature, deg F');
%print;
pause;
plot(x, Ts(1,:));
title('Shell side temperature');
xlabel('length, ft');
ylabel('temperature, deg F');
%print;
pause;

% define the constants used in the matrices:
%
at = 1/(2*delta_t) - Vt/(2*delta_x);
as = 0;
bt = -Pt/4;
bs = 1/(2*delta_t) + Vs/(2*delta_x) + Ps/4;
ct = 1/(2*delta_t) + Vt/(2*delta_x) + Pt/2;
cs = -Ps/2;
dt = -Pt/4;
ds = 1/(2*delta_t) - Vs/(2*delta_x) + Ps/4;
f1 = 1/(2*delta_t) + Vt/(2*delta_x) - Pt/2;
f2 = 1/(2*delta_t) - Vt/(2*delta_x);
f3 = 1/(2*delta_t) - Vs/(2*delta_x) - Ps/4;
f4 = 1/(2*delta_t) + Vs/(2*delta_x) - Ps/4;
%
% define the sparse a matrix
%
r1 = [1 1 1 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 6 6 6 6];
c1 = [1 2 3 1 2 3 2 3 4 5 2 3 4 5 4 5 6 7 4 5 6 7];
v1 = [bt ct dt bs cs ds at bt ct dt as bs cs ds at bt ct dt as bs cs ds];
r2 = [ 7 7 7 7 8 8 8 8 9 9 9 9 10 10 10 10 11 11 11 11 12 12 12 12];
c2 = [ 6 7 8 9 6 7 8 9 8 9 10 11 8 9 10 11 10 11 12 13 10 11 12 13];
v2 = [at bt ct dt as bs cs ds at bt ct dt as bs cs ds at bt ct dt];
v22 = [as bs cs ds];
r3 = [13 13 13 13 14 14 14 14 15 15 15 15 16 16 16 16 ];
c3 = [12 13 14 15 12 13 14 15 14 15 16 17 14 15 16 17];
v3 = [at bt ct dt as bs cs ds at bt ct dt as bs cs ds ];
r4 =[17 17 17 17 18 18 18 18 19 19 19 19 20 20 20 20];
c4 =[16 17 18 19 16 17 18 19 18 19 20 21 18 19 20 21];
v4 =[at bt ct dt as bs cs ds at bt ct dt as bs cs ds];
r5 =[21 21 21 21 22 22 22 22 23 23 23 23 24 24 24 24];
c5 =[20 21 22 23 20 21 22 23 22 23 24 25 22 23 24 25];
v5 =[at bt ct dt as bs cs ds at bt ct dt as bs cs ds];
r6 =[25 25 25 25 26 26 26 26 27 27 27 27 28 28 28 28];
c6 =[24 25 26 27 24 25 26 27 26 27 28 29 26 27 28 29];
v6 = v5;
r7 =[29 29 29 29 30 30 30 30 31 31 31 31 32 32 32 32];
c7 =[28 29 30 31 28 29 30 31 30 31 32 33 30 31 32 33];
v7 = v6;
r8 =[33 33 33 33 34 34 34 34 35 35 35 35 36 36 36 36];
c8 =[32 33 34 35 32 33 34 35 34 35 36 37 34 35 36 37];
v8 = v7;
r9 =[37 37 37 37 38 38 38 38 39 39 39 40 40 40];
c9 =[36 37 38 39 36 37 38 39 38 39 40 38 39 40];
v9 =[at bt ct dt as bs cs ds at bt ct as bs cs];
r =[r1 r2 r3 r4 r5 r6 r7 r8 r9];
c =[c1 c2 c3 c4 c5 c6 c7 c8 c9];
v =[v1 v2 v22 v3 v4 v5 v6 v7 v8 v9];

aa = sparse(r,c,v,40,40);

% define B.C.2:
%
Ttnew = 35;

% define B.C.3:
%
Tsold = 75;
%
% generate the right hand side vector
%
for j=1:iteration,
    B(1) = f1*Tt(j,1) + f2*Tt(j,2) - bt*(Ts(j,1)+Ts(j,2)) - at*Ttnew;
    B(2) = f3*Ts(j,1) + f4*Ts(j,2) - cs*Tt(j,1);

%    B(3) = f1*Tt(j,2) + f2*Tt(j,3) - bt*(Ts(j,2)+Ts(j,3));
%    B(4) = f3*Ts(j,2) + f4*Ts(j,3) - cs*Tt(j,2);
    for k=3:2:37,
	B(k) = f1*Tt(j,(k+1)/2) + f2*Tt(j,(k+1)/2+1) ...
	       - bt*(Ts(j,(k+1)/2)+Ts(j,(k+1)/2+1));
	B(k+1) = f3*Ts(j,(k+1)/2) + f4*Ts(j,(k+1)/2+1) - cs*Tt(j,(k+1)/2);
    end

    B(39) = f1*Tt(j,20) + f2*Tt(j,21) - bt*(Ts(j,21)+Tsold) - dt*Tsold;
    B(40) = f3*Ts(j,20) + f4*Ts(j,21) - ds*Tsold - cs*Tt(j,20);
%
% generate the solution vector
%
      TT = aa\B;
%
% update the right hand side
%
    Tt(j+1,1) = 35;
    for m=1:N,
        Tt(j+1,m+1) = TT(2*m);
	Ts(j+1,m) = TT(2*m-1);
    end
    Ts(j+1,m+1) = 75;
end

for n=1:5:iteration+1
    plot(x, Tt(n,:));
    hold on;
end
title('Dynamic Response of Tube Side to Step Disturbance');
xlabel('Length (ft)');
ylabel('Tube Side Water Temperature (degree F)');
print -deps fig812a.eps
pause;
hold off;

for n=1:5:iteration+1
    plot(x, Ts(n,:));
    hold on;
end
title('Dynamic Response of Shell Side to Step Disturbance');
xlabel('Length (ft)');
ylabel('Shell Side Water Temperature (degree F)');
print -deps fig812b.eps
hold off;

clear Tt, Tp;
