MATLAB: Adsorption Modelling – Solving PDE – Axial Dispersion Model

Question

Dear all,

I am trying to develop a simple model for an adsorption column. I have attached a file that contains the equations that I am trying to solve.

To solve the system, I am trying to solve the axial dispersion model (neglecting pressure drop, and hence, the dv/dz term) alongside the Linear Driving Force (LDF) model, which gives an expression for dq/dt in the axial dispersion PDE (see attached file). I am also using the Langmuir equation to give an expression for q* (equilibrium concentration) in the LDF equation.

Unknown variables are q and c, all other variables are known (i.e. D, v, epsilon etc.( / are the independent variables (i.e. time, t, and length, z).

I have tried using the Method of Lines to solve this system (again, see attached file), but I am aware that the system I have developed is flawed.

Any suggestions on how to approach this problem would be appreciated. I would preferably like to solve the system using either a PDE solver or ode45 applied to the method of lines.

My code so far is as follows:

Data File:

%%Constants
D = 3*10^-8;    % Axial Dispersion coefficient
v = 1*10^-3;    % Superficial velocity
epsilon = 0.4;  % Voidage fraction
k = 3*10^-5;    % Mass Transfer Coefficient
bH = 2.5*10^-5;  % Langmuir parameter for H2
L = 1;          % Column length
%%MoL Coefficients & Mesh Generation
N = 20;
dz = L/N;
z = [0:dz:L];
% MoL Coefficients from Discretisation
alph = ((D/(dz^2)) + ((v)/(2*dz)));
bet = ((-2*D)/(dz^2));
gamm = ((D/(dz^2)) - ((v)/(2*dz)));
%%Initial / Boundary Conditions
% Initial condition for ODE 1:
c0 = 0;         % C at t = 0 for all z
% Boundary Condition for ODE 1:
cFeed = 10;     % C at z = 0 for all t
qs = (1+bH*cFeed)/(bH*cFeed);         % Conc. of adsorbed phase at surface saturation
% Initial Condition for ODE 2:
q0 = qs*((bH*cFeed)/(1+bH*cFeed))-1;         % Initial concentration in adsorbed phase
% Time Boundaries
t0 = 0;
tf = 10000;
% Initial equilibrium conc.
qeq0 = qs*((bH*c0)/(1+bH*c0));
dq0 = k*qeq0;

Solving using ODE45:

close all
clear all
clc
% Call Data File
DataFile;
tspan = [0:10000];
[t,c] = ode45('cprime', tspan, c0);

Cprime file

function [ dcdt ] = cprime(~, c)
% Call data file to define constants
DataFile;
t=t0;
q = q0;
[dqdt] = qprime(t, q, k, qs, bH, c);
%%Create Tri-diagonal Matrix A
A = zeros(N,N);
for i = 2 : N-1         
    A(i, i-1)   = alph;
    A(i, i)     = bet;
    A(i, i+1)   = gamm;
end
A(1,1) = bet;
A(1,2) = gamm; 
A(N,N-1) = alph+gamm;
A(N-1, N) = gamm;
A(N,N) = bet;
b = ones(N,1);
b(1) = cFeed*alph-((1-epsilon)/(epsilon)*(k*(qeq-q)));
dcdt = (A*c + b);
end

qprime file:

function [dqdt] = qprime(~, q, k, qs, bH, c)
DataFile;
qeq = qs*((bH*c)/(1+bH*c));
dqdt = k*(qeq-q);
end

I am really struggling to figure out how to structure my code / approach these equations.

Any help would be muchly appreciated!

Thanks!

Answer 1

The code I provided here for a very similar problem might help:

https://de.mathworks.com/matlabcentral/answers/371313-error-in-solving-system-of-two-reaction-diffusion-equations

Best wishes

Torsten.