function PDEs2DS_Collocation
% 用對稱正交配置法求解固定床反應器二維擬均相穩態模型（二維穩態PDE方程組）
% (只在r方向無因次化)
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/03/07 $

clear all
clc

global A B a a1 a2 b1 b2 C1 C2 Rn T0 L N

% P3(x^2)對于圓柱對稱（a=2）的配置常數Aji和Bji
A = [-3.359794  5.2924315   -3.1010284  1.1683909
     -1.3980385 -1.5627540  4.3197367   -1.3589422
     0.69721650 -3.6766754  -1.1267583  4.1062172
     -1.2266754 5.4010626   -19.174383  15];
B = [-15.881426    19.636380   -5.2811862  1.5262327
     11.151861     -34.497415  29.235709   -5.890155
     -3.5405872    34.512110   -99.621159  68.649637
     -33.869987    136.21969   -252.37970  150];
  
Rn = [0.29763730  0.63989598  0.88750181  1];   % 全部配置點

% 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.0*1e3;  % CP', J/(kg K)
F = 130/3600;   % kg/s

% Equation coefficient（方程的系數）
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)
T0 = 873;
L = 1.053;
N = 3;              % 內配置點個數

y0 = [T0 T0 T0 893 0 0 0];  % y=[T1 T2 T3 T4 x1 x2 x3]
[z,y] = ode45(@Euqations,[0 L],y0);
z
T = y(:,1:N+1)
x = y(:,N+2:2*N+1);
for i = 1:length(z)
    xb(i) = - sum(A(N+1,1:N).* x(i,:)/A(N+1,N+1));
end
x = [x xb']

% 求沿管長的平均轉化率xa(i)
for i=1:length(z)
    xn = x(i,:);
    xa(i) = quadl(@func,0,1,[],[],Rn,xn)/(1^2/2)
end
L = spline(xa,z,0.45)

% Plot the results
surf(Rn*a,z,T)      % 反應管軸徑向溫度分布
xlabel('r (m)')
ylabel('z (m)')
zlabel('T (K)')
figure
plot(z,xa)          % 平均轉化率沿管長的分布圖
xlabel('z (m)')
ylabel('x_a_v')
figure
surf(Rn*a,z,x)      % 軸徑向平均轉化率分布
xlabel('r (m)')
ylabel('z (m)')
zlabel('x')

% ------------------------------------------------------------------
function dydx = Euqations(x,y)
global A B a a1 a2 b1 b2 C1 Rn T0 L N
T = y(1:N+1);
x = y(N+2:2*N+1);
rc = ReactionRate(T(1:N),x);
for i = 1:N
    dTdR(i) = a1/a^2* sum( (B(i,:)+A(i,:)./Rn(i)) .* T' ) - b1*rc(i);  
end

dTdR(N+1) = C1/a*sum(A(N+1,:).* T')

xb = - sum(A(N+1,1:N).* x'/A(N+1,N+1));
for i=1:N
    dxdR(i) = a2/a^2*( sum((B(i,1:N)+A(i,1:N)./Rn(i)).*x(1:N)') ...
        + (B(i,N+1)+A(i,N+1)./Rn(i)).*xb ) + b2*rc(i);
end
dydx = [dTdR dxdR]';

% ------------------------------------------------------------------
function f = ReactionRate(T,x)      % 計算反應速度
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 y = func(R,Rn,xn)
x = spline(Rn,xn,R);
y = R.*x;