function PDEs2DS_CrankNicolson
% 用有限差分求解固定床反應(yīng)器二維擬均相穩(wěn)態(tài)模型（二維穩(wěn)態(tài)PDE方程組）
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2002/05/02 $

global n F G rc dz M M1 b1 b2 C2 C3 C4
% Increment of r and z
n = 5+1;            % node number in r direction is 5 
m = 100;
dr = 0.01;          % m
dz = 0.0585;        % m

% Parameters
Ramda = 0.45;       % W/(m K)
G = 2500;           % kg/(m^2 h)
Cp = 2.18;          % kJ/(kg K)
dH = 140e3;         % J/mol
rho = 1440;         % kg/m^3
a = 0.05;           % radius of reactor, m
u0c0 = 0.069/(pi*a^2);  % u0*c0 obtained from u0*c0*A=0.069 kmol/h
Cp1 = 1.0e3;        % 煙道氣比熱，J/(kg K)
F = 130/3600;       % kg/s

% Equation coefficient（方程的系數(shù)）
a1 = Ramda/(G*Cp*1000)*3600;    % m
b1 = dH*rho/G/Cp;       % Note of unit accordance：dH*rho/G/Cp rc = [K/m]，rc=[kmol/(kg h)]
a2 = 0.000427;          % a2 = D2/mu
b2 = rho/u0c0;
           
C1 = 2*pi*a*Ramda/(F*Cp1);      % formula (12)
C2 = 4*dr/(dz*C1);              % formula (33)
M = a1*dz /dr^2;                % formula (17)
M1 = a2*dz /dr^2;               % formula (19)
C3 = M*C2/2*(1+0.5/(n-1))-(1+M);  
C4 = M*C2/2*(1+0.5/(n-1))-(1-M);

% Initializing variables, at the inlet of the reactor
% (在反應(yīng)器入口處即z=0或j=1時初值)
T(1:n-1) = 873;
T(n) = 893;
x(1:n) = 0;
X0 = [T;x];
Tresult(1,:) = X0(1,:);
xresult(1,:) = X0(2,:);

% Calculate F(i,j),G(i,j) and rc(i,j) at z=0 (initial value)
% 計算對應(yīng)于j=1時的F(i,1)和G(i,1)及反應(yīng)速度rc
FGrc(T,x);

for j=2:m
    X = fsolve(@TxEquations,X0);    % Solve nonlinear equations
    T = X(1,:);
    x = X(2,:);
    Tresult(j,:) = T
    xresult(j,:) = x
    xa(j,:) = sum(2*xresult(j,2)+4*xresult(j,3)+6*xresult(j,4)...
        +8*xresult(j,5)+5*xresult(j,6))/25 
    if xa(j)>0.45       % When the average conversion > 0.45, stop calculating and exit the loop
        break
    end
    
    % Calculate F(i,j),G(i,j) and rc(i,j) for the next iteration
    % 計算F(i,j)和G(i,j)及反應(yīng)速度rc(i,j)供下次迭代使用
    FGrc(T,x)
   
end
z = dz.*[0:j-1];
r = dr.*[0:n-1];
% When the  conversion rate of ethylbenzene is 45%, the reactor length required is:
z = z'
L = spline(xa,z,0.45)

% Plot the results
surf(r,z,Tresult)   % 反應(yīng)管軸徑向溫度分布
xlabel('r (m)')
ylabel('z (m)')
zlabel('T (K)')
figure
plot(z,xa)          % 平均轉(zhuǎn)化率沿管長的分布圖
xlabel('z (m)')
ylabel('x_a_v')
figure
surf(r,z,xresult)   % 軸徑向平均轉(zhuǎn)化率分布
xlabel('r (m)')
ylabel('z (m)')
zlabel('x')

% --------------------------------------------------------------------------
function f = ReactionRate(T,x)
% Calculate the reaction rate（計算反應(yīng)速度）
k = 0.027*exp(0.021*(T-773));
f = 15100*exp(-11000./T).*((1-x)./(11+x)-1.2*x.^2./k./(11+x).^2);

% --------------------------------------------------------------------------
function f = FGrc(T,x)
% Calculate F(i), G(i)
global F G rc n dz M M1 b1 b2 C2 C3 C4
rc = ReactionRate(T,x);         % Calculate the reaction rate

F(1) = ((1-2*M)*T(1)+2*M*T(2)-b1*dz/2*rc(1))/(1+2*M);       % formula (28)
G(1) = ((1-2*M1)*x(1)+2*M1*x(2)+b2*dz/2*rc(1))/(1+2*M1);    % formula (29)
i = (2:5);
    var1 = 1-0.5./(i-1);
    var2 = 1+0.5./(i-1);
    F(i) = ( M/2*( var1.*T(i-1)+var2.*T(i+1) )  + (1-M)*T(i) -b1*dz/2*rc(i) )/(M+1);    % formula (22)
    G(i) =( M1/2*( var1.*x(i-1)+var2.*x(i+1) )  + (1-M1)*x(i)+b2*dz/2*rc(i) )/(M1+1);   % formula (23)
F(n) = ( -M*T(n-1) +C4*T(n)+b1*dz/2*rc(n) )/ C3;            % formula (35) 
G(n) = ( M1*x(n-1)+(1-M1)*x(n)+b2*dz/2*rc(n) )/(1+M1);      % formula (31)

% --------------------------------------------------------------------------
function f = TxEquations(X)
% Define the linear equation system composed of (20)、(21) 、(26)、(27)、(30) and (34).
global n F G rc dz M M1 b1 b2 C2 C3
T = X(1,:)     
x = X(2,:)     
fT(1) = F(1) + (2*M*T(2)-b1/2*dz*rc(1) )/(1+2*M) - T(1);        % formula (26)
fx(1) = G(1) + (2*M1*x(2)+b2/2*dz*rc(1) )/(1+2*M1) - x(1);      % formula (27)
for i=2:n-1
    var1(i) = (1-0.5/(i-1));
    var2(i) = (1+0.5/(i-1));
    fT(i) = F(i)+ ( M/2*(var1(i)*T(i-1)+var2(i)*T(i+1))-b1/2*dz*rc(i) )/(M+1) - T(i);   % formula (20)
    fx(i) = G(i)+ ( M1/2*(var1(i)*x(i-1)+var2(i)*x(i+1))+b2/2*dz*rc(i) )/(M1+1) - x(i); % formula (21)
end
fT(n) = F(n) + ( -M*T(n-1)+b1/2*dz*rc(n) )/C3 -T(n);            % formula (34)
fx(n) = G(n) + ( M1*x(n-1)+b2/2*dz*rc(n) )/(1+M1) - x(n);       % formula (30)
f = [fT;fx];