kinetics_crank_nicolson.m

有限差分法（crank_nicholson finite difference ) 求解色谱模型

function kinetics_crank_nicolson
% 使用Crank-Nicolson有限差分方法求解一维动态传热模型
%
% calling the function U = CrankNicolson(f,c1,c2,a,b,c,n,m) to solve the heat equation: 
%   du(x,t)/dt = α * d2u(x,t)/dx2 over R={(x,t):0<=x<=a,0<=t<=b}
%   with u(x,0)=f(x), for 0<=x<=a, and u(0,t)=c1,u(a,t)=c2,for 0<=t<=b.
% Input  - f=u(x,0) as a string 'f'
%           - c1=u(0,t) and c2=u(a,t)
%           - a and b right end points of [0,a] and [0,b]
%           - c the constant in the heat equation
%           - n and m number of grid points over [0,a] and [0,b]
% Output    - U solution matrix
%
%   Author: LIU jianqi
%   Copyright 2006,   
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2006/09/11 $

% Increment of r and z
clc
tic
a =1; % nondimensionless x
b =50; % nondimensionless time t
n =1000; % discrete dx
m =200;  % discrete dt
h=a/(n-1);
k=b/(m-1);
% input parameters
L=170; % cm
uv=120; % mL/min
R=5.2; % cm
Eb=0.345; % external porosity
t_bi=L*Eb/(uv/(pi*R*R/4)); % nondimensionless time -> time

c1=100;

Pe=50;
St=1;

porosity=0.345;
F=(1-porosity)/porosity;
N_a=-1/Pe;
N_b=F*St;
lamda=k*St;
alpha=N_a*k/(2*h*h);
beta=k/(4*h);
gama=N_b*k;
clear N_a N_b

theta=(alpha-beta)/(2*alpha-1);
fai=(alpha+beta)/(2*alpha-1);
yita=gama/(2*alpha-1);
keli=(2*alpha+1)/(2*alpha-1);

U=CrankNicolson(@ic,c1,theta,fai,yita,keli,lamda,a,b,n,m);
toc
t=0:k:b;
plot(t'*t_bi,U(:,end))
out=[t'*t_bi,U(:,end)];

% ------------------------------------------------------------------
function f = ic(x)
f = 0;