kinetics_crank_nicolson.asv

有限差分法（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
a =1; % nondimensionless x
b =5; % nondimensionless time t
n =50; % discrete dx
m =50;  % discrete dt
h=1/(n-1);
k=b/(m-1);
% input parameters
L=170; % cm
uv=120; % mL/min
R=5; % cm
Eb=0.34; % external porosity
t_bi=L*Eb/(uv/(pi*R*R/4)); % nondimensionless time -> time

c1=10;

Pe=30;
St=10;

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

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

U=CrankNicolson(@ic,c1,theta,fai,yita,keli,n,m);

t=0:dt:b;
out=[t'*t_bi,U(:,end)/c1];

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