Packed Bed Absorber - Axial Dispersion Model

2 次查看(过去 30 天)
% Initialisation
close all
clear all
clc
%System Set-up %
%Define Variables
D = 0.006; % Axial Dispersion coefficient
v = 1.061; % Superficial velocity
epsilon = 0.3; % Voidage fraction
k = 0.1; % Mass Transfer Coefficient
cFeed = 2; % Feed concentration
H=2;
L = 10; % Column length
t0 = 0; % Initial Time
tf = 50; % Final time
dt = 0.006; % Time step
dx=0.01;
z = [0:0.01:L]; % Mesh generation
t = [t0:dt:tf];% Time vector
n = numel(z); % Size of mesh grid
%Initial Conditions / Vector Creation
c0 = zeros(n,1);
c0(1) = cFeed;
q0 = zeros(n,1); % t = 0, q = 0 for all z,
y0 = [c0 ; q0]; % Appends conditions together
%ODE15S Solver
[T, Y] = ode15s(@(t,y) MyFun(t,y,z,n),t,y0);
function DyDt=MyFun(~, y, z, n)
% Defining Constants
D = 0.006; % Axial Dispersion coefficient
v = 1.061; % Superficial velocity
epsilon = 0.3; % Voidage fraction
k = 0.1; % Mass Transfer Coefficient
H=2;
dx=0.01;
% Variables being allocated zero vectors
c = zeros(n,1);
q = zeros(n,1);
DcDt = zeros(n,1);
DqDt = zeros(n,1);
DyDt = zeros(2*n,1);
DcDz = zeros(n,1);
D2cDz2 = zeros(n,1);
c = y(1:n);
q = y(n+1:2*n);
for i=2:n-1
DcDz(i) = (c(i+1)-c(i-1))/(2*dx);
D2cDz2(i)= (c(i-1)-(2*c(i))+c(i+1))/(dx^2) ;
end
% Calculate DcDz and D2cDz2 at z=L for boundary condition dc/dz = 0
% Standard setting for q
DqDt(1) = k*((H*c(1))-q(1));
DcDt(1) = 0.0;
% Set time derivatives in remaining points
for i=2:n
DqDt(i) = k*((H*c(i))-q(i));
%Equation: dc/dt = D * d2c/dz2 - v*dc/dz - ((1-e)/(e))*dq/dt
DcDt(i)=D*D2cDz2(i)-v*DcDz(i)-2.33*DqDt(i);
end
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end
I dont know whats wrong with my code. Any help will be appreciated Screenshot 2019-06-25 at 4.24.23 PM.png
BC C=2 @ t>0 @z=0
IC @ t=0 &Z= 0 to L ----> c=0 q=0
  1 个评论
Torsten
Torsten 2019-7-15
You only calculated DcDz and D2cDz2 for i=2,...,n-1. In your time-derivative calculation, you need DcDz and D2cDz2 for i=n which does not exist.

请先登录,再进行评论。

回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Sparse Matrices 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by