% This program uses the centered-difference technique to solve the 
% unsteady-state steam heat exchanger problem (explicit) in section 8.2.
%

global Cp rho vz Tsteam D U;

% simulation parameters
%
Cp = 1.0;
D = 0.5;
rho = 62.4; 
vz = 3;
U = 0.2;
TT0 = 75;
Tsteam = 250;
L = 12;
N = 40;			% divide the distance grid into N elements.
delta_z = L/N;	
delta_t = delta_z/vz;


% the initial steady state temp. can also be solved analytically

for I = 2:(N+1),
    T(1,I) = Tsteam + (TT0-Tsteam)*exp(-4*U*(I-1)*delta_z/(rho*Cp*D*vz));
end

T(1,1) = 75;
z = [0:delta_z:12]';
plot(z, T(1,:), '+');
title('Dynamic Response of Steam Heat Exchanger to Step Disturbance');
xlabel('Reactor Length (ft)');
ylabel('Water Temperature (degree F)');
pause;
hold on;

% define the step change in inlet temp.
%
for K = 2:(N+1),
    T(K,1) = 60;
end

% constants for iteration:
%
a = 1/(2*delta_t);
b = vz/(2*delta_z);
c = 2*U/(rho*Cp*D);

% loop for time iteration
%
for i = 1:N,

    % loop for distance iteration
    %
    for n = 1:N,
	T(i+1, n+1) = (2*c*Tsteam + (b-a)*T(i+1,n) + (-b+a)*T(i,n+1) ...
		      + (b+a-c)*T(i,n))/(a+b+c);
    end
end

plot(z, T(3,:), ':');
plot(z, T(11,:), '-.');
plot(z, T(21,:), '-');
plot(z, T(31,:), '--');
plot(z, T(41,:), 'o');
text(0.5, 92, '+++   t=0')
text(0.5, 90, '...   t=0.2 sec')
text(0.5, 88, '-.-   t=1.0 sec')
text(0.5, 86, '___   t=2.0 sec')
text(8, 65, '- - -   t=3.0 sec')
text(8, 63, 'ooo   t=4.2 sec')

hold off;
print -deps fig85.eps;
