Coupled Ion Exchange - Electrodialysis model using partial differential equations
Show older comments
Dear All,
I am trying to model a coupled Ion Exchange - Electrodialysis model. This is similar to the adsorption model. The extra term in this model is that of transport of ions due to the electric field. The first snapshot I have attached shows the original model. The second snapshot shows the rearranged form of the model to match the adsorption model (without axial dispersion). In this second snapshot you can see the extra term containing 'phi', the potential gradient.
In the adsorption model, the 2 independent variables are time and one direction. The 2 dependent variables are bulk concentration and solid phase concentration. In this model, the 3 independent variables are time and 2 directions. 1 direction (x) is the direction the ions move under the influence of the electric field. The second direction (z) is the direction that the bulk fluid moves in. The 3 dependent variables are the bulk concentration, solid phase concentration and the electric potential.
The sample code for the adsorption model is given below (Taken from https://uk.mathworks.com/matlabcentral/answers/385756-adsorption-modelling-solving-pde-axial-dispersion-model. Coded by Ashley Gratton, Torsten). Please let me know how I can change this for my model which has 3 independent and 3 dependent variables.
Thank you in advance.
function main
% Initialisation
close all
clear all
clc
%System Set-up %
%Define Variables
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
b = 2.5*10^-5; % Langmuir parameter
qs = 2*10^4; % Saturation capacity
cFeed = 10; % Feed concentration
L = 1; % Column length
t0 = 0; % Initial Time
tf = 1000; % Final time
dt = 0.5; % Time step
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, this makes sense to me
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 = 3*10^-8; % Axial Dispersion coefficient
v = 1*10^-3; % Superficial velocity
epsilon = 0.4; % Voidage fraction
k = 3*10^-5; % Mass Transfer Coefficient
b = 2.5*10^-5; % Langmuir parameter
qs = 20000; % Saturation capacity
% 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);
zhalf = zeros(n-1,1);
DcDz = zeros(n,1);
D2cDz2 = zeros(n,1);
c = y(1:n);
q = y(n+1:2*n);
% Interior mesh points
zhalf(1:n-1)=(z(1:n-1)+z(2:n))/2;
for i=2:n-1
DcDz(i) = ((z(i)-z(i-1))/(z(i+1)-z(i))*(c(i+1)-c(i))+(z(i+1)-z(i))/(z(i)-z(i-1))*(c(i)-c(i-1)))/(z(i+1)-z(i-1));
D2cDz2(i) = (zhalf(i)*(c(i+1)-c(i))/(z(i+1)-z(i))-zhalf(i-1)*(c(i)-c(i-1))/(z(i)-z(i-1)))/(zhalf(i)-zhalf(i-1));
end
% Calculate DcDz and D2cDz2 at z=L for boundary condition dc/dz = 0
DcDz(n) = 0;
D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*(c(n)-c(n-1))/(z(n)-z(n-1));
% Set time derivatives at z=0
% DcDt = 0 since c=cFeed remains constant over time and is already set as initial condition
% Standard setting for q
DqDt(1) = k*( ((qs*b*c(1))/(1 + b*c(1))) - q(1) );
DcDt(1) = 0.0;
% Set time derivatives in remaining points
for i=2:n
%Equation: dq/dt = k(q*-q) where q* = qs * (b*c)/(1+b*c)
DqDt(i) = k*( ((qs*b*c(i))/(1 + b*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) - ((1-epsilon)/(epsilon))*DqDt(i);
end
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end
34 Comments
Zain Shabeeb
on 9 Jul 2019
Edited: Zain Shabeeb
on 9 Jul 2019
Torsten
on 10 Jul 2019
When looking at your model, the following questions arose:
c_i^b only depends on z ?
c_i^hat depends on z and x ?
d_phi_dx is only a single number (thus not dependent on x or z) ?
Boundary conditions for c_i^b at z = 0 ?
Boundary conditions for c_i^hat at x = 0 ?
How large is "n" in d_phi_dx ?
Zain Shabeeb
on 10 Jul 2019
Torsten
on 10 Jul 2019
ci^hat (solid phase concentration) changes only with x direction (please let me know if your understanding of the model suggests otherwise).
Since c_i^hat depends on c_i^b, c_i^hat must depend on z and x.
I was given this relationship of dphidx from literature. I understand that it is dependent on the x direction.
In your notation, I don't see a dependence on x for dphidx. It is a single number since all the space-dependent contributions are integrated out.
Zain Shabeeb
on 10 Jul 2019
Torsten
on 10 Jul 2019
How can I proceed with modelling this system?
First by clarifying the questions - maybe by contacting the authors.
Zain Shabeeb
on 10 Jul 2019
Your model is too much different from the one you want to take as reference. If I were you, I'd start programming your model from scratch.
If z_i (i=1,...,Z) is a discretization for the interval [0;1] and x_j (j=1,...,X) is a disretization for the interval [0;d], you will have to solve for c_k^b at z=z_i (i=1,...,Z) and c_k^hat at z=z_i and x=x_j(i=1,...,Z;j=1,...,X).
Discretize dc_k^b/dz at z=z_i and dc_k^hat/dx with a first-order upwind scheme.
The integral in the computation of dphidx as a single value can be evaluated using the trapezoidal rule, e.g.
This gives Z + X*Z ordinary differential equations that can be solved using ODE15S.
Good luck.
Zain Shabeeb
on 10 Jul 2019
Zain Shabeeb
on 10 Jul 2019
Zain Shabeeb
on 10 Jul 2019
Torsten
on 10 Jul 2019
Should read
c_k^b(i-1)
and
c_k^hat(i,j-1)
thus the solution variables at z_(i-1) resp. (z_i,x_(j-1)).
Zain Shabeeb
on 10 Jul 2019
Edited: Zain Shabeeb
on 10 Jul 2019
Zain Shabeeb
on 10 Jul 2019
Torsten
on 10 Jul 2019
With open questions I mean especially the dphidx-term. I doubt that it is constant for all positions.
Zain Shabeeb
on 10 Jul 2019
Torsten
on 11 Jul 2019
Count the number of variables in the call list to EDIfun2 and in the receive list. They do not coincide.
Zain Shabeeb
on 11 Jul 2019
Zain Shabeeb
on 11 Jul 2019
I don't understand why you can't identify such obvious programming errors on your own.
In EDIfun2, you set
q = y(n+1:n+n*m);
and then use q in
DqDx(i,j) = (q(i,j)-q(i,j-1))/(x(j)-x(j-1));
Don't you see the problem ?
Please don't immediately post every single error - start thinking and use MATLAB's debugger.
Zain Shabeeb
on 11 Jul 2019
Zain Shabeeb
on 11 Jul 2019
q=0 at t=0, and you divide by q when calculating DphiDx.
Again: Print Dydt in EDIfun2 until it contains reasonable values and use the debugger.
Note that according to your model equations, DphiDx does not depend on x and z . In your code from above, it does because you replace the integrand by its value in (z(i),x(j)).
Zain Shabeeb
on 11 Jul 2019
Why should DphiDx be a matrix of 0's ? If you start with q=1, it is a single number different from 0. And if it is a single number that does not depend on z and x, it will affect DJdx.
Note further that you use DqDt(i,j) in your loop before you define it.
Zain Shabeeb
on 11 Jul 2019
Torsten
on 11 Jul 2019
I defined DqDt and DcDt in the function file at the top. Here is the section of code which shows where i defined it:
This is only an initialization of the arrays.
What I mean is
for i = 2:n
for j = 2:m
DcDt(i) = -v*DcDz(i) + (((1-epsilon)/epsilon)*(DqDt(i,j) + DJDx(i,j)));
DqDt(i,j) = a*D*(q(i,j)-(cstar(i,j))) - DJDx(i,j);
end
end
You use DqDt(i,j) in the calculation of DcDt(i) before you define it as
DqDt(i,j) = a*D*(q(i,j)-(cstar(i,j))) - DJDx(i,j);
Zain Shabeeb
on 11 Jul 2019
Torsten
on 11 Jul 2019
If I have q = 0 then I cant divide it in the DphiDx equation. What would you suggest for this problem?
I don't know about the physical background of your problem.
So I can't tell what makes sense here as initial condition for q.
Zain Shabeeb
on 11 Jul 2019
Zain Shabeeb
on 13 Jul 2019
Hi,
I changed my code to have q vary only in the x direction because I did not see any variation of q in the z direction with the previous code.
But my problem remains. I have tried to debug it but I just don't understand why DcDz is reasonable in the z direction but DqDx is giving a strange result. My code is given below. I think the problem is in the discretization of DqDx. I have pasted my code below and the graph of DcDz versus z as well as DqDx versus x.
Code:
function main
clear all
%% Defining constants
W = 0.1; %Width of column
L = 1; %Length of column
a = 2; %Surface area to volume ratio
delta = 0.00001; %Film thickness
D = 3*10^-10; %Film diffusion coefficient
v = 0.1; %Velocity of fluid
epsilon = 0.4; %Bed porosity
K = 10; %Equilibrium constant between solid phase and interface
u = 0.01; %Mobility
Itot = 1; %Total current
F = 96485.3329; %Faraday constant
cHstar = 3; %Interfacial concentration of counter ion (hydrogen)
qH = 2; %Solid phase concentration of counter ion (hydrogen)
cFeed = 0.01; %Feed concentration
%% Initializing Vectors for Independent variables
t0 = 0; %Initial time
tf = 10; %Final time
dt = .1; %Time step
t = [0:dt:tf]; %Time vector
z = [0:0.1:L]; %Mesh generation for z direction
x = [0:0.01:W]; %Mesh generation for x direction
n = numel(z); %Number of elements in z
m = numel(x); %Number of elements in x
%% Initialization of Dependent variables
c0 = zeros(n,1) ; %Initial bulk concentration vector
c0(1) = cFeed; %Initial bulk concentration at z = 0
q0 = 0.00001*ones(m,1); %Initial solid phase concentration matrix. Dependent on x and z directions.
J0 = zeros(m,1); %Initial Solid phase flux matrix. Dependent on x and z directions.
y0 = [c0;q0]; %Initial y vector appending together c0 and q0
[t,y] = ode15s(@(t,y)EDIfun3(t,y,z,x,n,m,W,L,a,delta,D,v,epsilon,K,u,Itot,F,cHstar,qH,cFeed),t,y0);
%Yc = Y(:,[1:n]);
%plot(z,Y(:,[1:11]));
%Yq = Y(:,[n+1:n+n*m]);
%Yqr = reshape(Yq(2,:),n,m);
function DyDt=EDIfun3(t,y,z,x,n,m,W,L,a,delta,D,v,epsilon,K,u,Itot,F,cHstar,qH,cFeed);
%% Variables being allocated zero vectors
c = zeros(n,1); %Bulk concentration vector
q = zeros(m,1); %Solid phase concentration matrix. Dependent on x and z directions.
J = zeros(m,1); %Solid phase flux matrix. Dependent on x and z directions.
cstar = zeros(m,1); %Cstar variable representing
DcDt = zeros(n,1); %Time differential vector of bulk concentration
DqDt = zeros(m,1); %Time differential matrix of solid phase concentration. Dependent on x and z directions.
DyDt = zeros(n+m,1); %Time differential vector of concatenated c and q
%% Discretizing in x and z directions
DcDz = zeros(n,1); %Initial vector for bulk concentration change with z direction
DqDx = zeros(m,1); %Initial matrix for solid phase concentration change in x and z directions
DJDx = zeros(m,1); %Initializing matrix for solid phase flux change in x and z directions
DphiDx = zeros(m,1); %Initializing matrix for potential gradient change in x and z directions
c = y(1:n); %Initial vector for bulk concentration allocated to y vector
q = y(n+1:n+m); %Initial vector for solid phase concetration allocated to y vector
DcDz(1) = 0; %Initial bulk concentration differential
DqDx(1) = 0; %Initial solid phase concentration differential
J(1) = 0; %Initial solid phase flux
DphiDx(1) = 0; %Initial potential gradient
for i = 2:n-1
DcDz(i) = (c(i)-c(i-1))/(z(i)-z(i-1)); %Discretization of bulk concentration with z direction
end
for j = 2:m-1
DqDx(j) = (q(j)-q(j-1))/(x(j)-x(j-1)); %Discretization of solid phase concentration with x and z direction
DphiDx(j) = -Itot/(F*(1-epsilon)*W*L*u*q(j)); %Equation for DphiDx
J(j) = -u*q(j)*DphiDx(j); %Equation for solid phase flux
DJDx(j) = (J(j)-J(j-1))/(x(j)-x(j-1)); %Discretization of solid phase flux with x and z direction
cstar(j) = (q(j)*cHstar)/(K*qH); %Equation for cstar
end
DcDz(n) = (c(n)-c(n-1))/(z(n)-z(n-1)); %Final bulk concentration differential
DqDx(m) = (q(m)-q(m-1))/(x(m)-x(m-1)); %Final solid phase concentration differential
DJDx(m) = (J(m)-J(m-1))/(x(m)-x(m-1)); %Final solid phase flux
DphiDx(m) = -Itot/(F*(1-epsilon)*W*L*u*q(m)); %Final potential gradient
%DphiDx(m) = 0;
%% Differential equation with time
DcDt(1) = 0; %Initial bulk concentration differential with time
DqDt(1) = 0; %Initial solid phase concentration differential with time
for i = 2:n
for j = 2:m
DqDt(j) = a*D*(q(j)-(cstar(j))) - DJDx(j);
DcDt(i) = -v*DcDz(i) + (((1-epsilon)/epsilon)*(DqDt(j) + DJDx(j)));
end
end
DyDt = [DcDt;DqDt]
DJDx
DphiDx
-v*DcDz
((1-epsilon)/epsilon)*(DqDt + DJDx)
cstar
q
DqDx
DqDt
y(:,:)
end
figure;
subplot(3,3,1)
plot(x,y(:,[12:22]));
title('Solid concentration vs x direction')
xlabel('x');
ylabel('Solid phase concentration');
subplot(3,3,2)
plot(z,y(:,[1:11]));
title('Liquid concentration vs z direction')
xlabel('z');
ylabel('Liquid phase concentration');
subplot(3,3,3)
plot(x,DqDx);
title('DqDx vs x direction')
xlabel('x');
ylabel('DqDx');
subplot(3,3,4)
plot(x,DphiDx);
title('DphiDx vs x direction')
xlabel('x');
ylabel('DphiDx');
subplot(3,3,5)
plot(t,y(:,[12:22]));
title('q vs t')
xlabel('t');
ylabel('q');
subplot(3,3,6)
plot(x,J);
title('J vs x direction')
xlabel('x');
ylabel('J');
subplot(3,3,7)
plot(x,DJDx);
title('DJDx vs x direction')
xlabel('x');
ylabel('DJDx');
subplot(3,3,8)
plot(x,cstar);
title('cstar vs x direction')
xlabel('x');
ylabel('cstar');
subplot(3,3,9)
plot(z,DcDz);
title('DcDz vs z direction')
xlabel('z');
ylabel('DcDz');
end
DcDz vs z:
DqDx vs x:
Answers (1)
hbulusoy
on 17 Sep 2019
0 votes
Did you have any luck solving your problem? I'm very interested in this code. I have been trying to figure out to code this model in Python for a while but I am also not getting very good results.
1 Comment
Zain Shabeeb
on 17 Sep 2019
Categories
Find more on Structural Mechanics in Help Center and File Exchange
on 9 Jul 2019
on 17 Sep 2019
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
