# Coupled Ion Exchange - Electrodialysis model using partial differential equations

15 views (last 30 days)
Zain Shabeeb on 9 Jul 2019
Commented: Zain Shabeeb on 17 Sep 2019
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.
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

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?
So I can't tell what makes sense here as initial condition for q.
Zain Shabeeb on 11 Jul 2019
I have spent some time analyzing my solution. The variation in the z direction seems fine as I am getting a reasonable result for DcDz. However, the x direction is giving strange results. For example DJDx, DphiDx. I have printed DJDx and DphiDx so that it shows their values when you run the code.
Could this be because of the way I have discretized DJDx and DphiDx? I used the upwind method that you provided in one of your earlier comments.
Data file:
%% 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^-1; %Film diffusion coefficient
v = 0.1; %Velocity of fluid
epsilon = 0.9; %Bed porosity
K = 0.001; %Equilibrium constant between solid phase and interface
u = 2; %Mobility
Itot = 5; %Total current
cHstar = 1; %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 = 0.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.00000001*ones(n,m); %Initial solid phase concentration matrix. Dependent on x and z directions.
phi0 = zeros(n,m); %Initial electric potential matrix. Dependent on x and z directions.
J0 = zeros(n,m); %Initial Solid phase flux matrix. Dependent on x and z directions.
y0 = [c0;reshape(q0,n*m,1)]; %Initial y vector appending together c0 and q0
[T,Y] = ode15s(@(t,y)EDIfun2(t,y,z,x,n,m,W,L,a,delta,D,v,epsilon,K,u,Itot,F,cHstar,qH,cFeed),t,y0);
%plot(z,Yc(:,:));
%Yc = Y(:,[1:n]);
%Yq = Y(:,[n+1:n+n*m]);
%Yqr = reshape(Yq(2,:),n,m);
Function file:
function DyDt=EDIfun2(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(n,m); %Solid phase concentration matrix. Dependent on x and z directions.
phi = zeros(n,m); %Potential matrix. Dependent on x and z directions.
J = zeros(n,m); %Solid phase flux matrix. Dependent on x and z directions.
cstar = zeros(n,m); %Cstar variable representing
DcDt = zeros(n,1); %Time differential vector of bulk concentration
DqDt = zeros(n,m); %Time differential matrix of solid phase concentration. Dependent on x and z directions.
DyDt = zeros(n+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(n,m); %Initial matrix for solid phase concentration change in x and z directions
DJDx = zeros(n,m); %Initializing matrix for solid phase flux change in x and z directions
DphiDx = zeros(n,m); %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 = reshape(y(n+1:n+n*m),n,m); %Initial vector for solid phase concetration allocated to y vector
DcDz(1) = 0; %Initial bulk concentration differential
DqDx(1,1) = 0; %Initial solid phase concentration differential
J(1,1) = 0; %Initial solid phase flux
DphiDx(1,1) = 0; %Initial potential gradient
for i = 2:n-1
for j = 2:m-1
DcDz(i) = (c(i)-c(i-1))/(z(i)-z(i-1)); %Discretization of bulk concentration with z direction
DqDx(i,j) = (q(i,j)-q(i,j-1))/(x(j)-x(j-1)); %Discretization of solid phase concentration with x and z direction
DphiDx(i,j) = Itot/(F*(1-epsilon)*W*L*u*q(i,j)); %Equation for DphiDx
J(i,j) = u*q(i,j)*DphiDx(i,j); %Equation for solid phase flux
DJDx(i,j) = (J(i,j)-J(i,j-1))/(x(j)-x(j-1)); %Discretization of solid phase flux with x and z direction
cstar(i,j) = (q(i,j)*cHstar)/(K*qH); %Equation for cstar
end
end
DcDz(n) = (c(n)-c(n-1))/(z(n)-z(n-1)); %Final bulk concentration differential
DqDx(n,m) = (q(n,m)-q(n,m-1))/(x(m)-x(m-1)); %Final solid phase concentration differential
DJDx(n,m) = (J(n,m)-J(n,m-1))/(x(m)-x(m-1)); %Final solid phase flux
DphiDx(n,m) = Itot/(F*(1-epsilon)*W*L*u*q(n,m)); %Final potential gradient
%% Differential equation with time
DcDt(1) = 0; %Initial bulk concentration differential with time
DqDt(1,1) = 0; %Initial solid phase concentration differential with time
for i = 2:n
for j = 2:m
DqDt(i,j) = a*D*(q(i,j)-(cstar(i,j))) - DJDx(i,j);
DcDt(i) = -v*DcDz(i) + (((1-epsilon)/epsilon)*(DqDt(i,j) + DJDx(i,j)));
end
end
DyDt = [DcDt;reshape(DqDt,n*m,1)]
DJDx
DphiDx
end
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
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: hbulusoy on 17 Sep 2019
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
Hi,
I did not solve this problem on MATLAB. However I believe it can be solved on GAMS software because it solves all equations simultaneously, whereas MATLAB looks at each line one by one. Good luck.
Zain