Hi
I had a go at editing the code given above. I know there are many mistakes in this. The main one I see is the discretization in the x direction. My system contains an x direction and I don't know how to discretize in the x direction for my problem. Please let me know how I can get this code to run and if it is even possible to use this code for my problem. I have added 'EDIT' comments where I have changed the code.
Thanks in advance.
% Initialisation
close all
clear all
clc
%System Set-up %
%Define Variables
a = 0.4;
delta = 0.0001;
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
u = 0.1;
Itot = 4;
F = 96485.3329;
cHstar = 2*10^-5;
qH = 2*10^-4;
cFeed = 10; % Feed concentration
L = 1; % Column length
W = 1; % EDIT: Column Width
t0 = 0; % Initial Time
tf = 1000; % Final time
dt = 0.5; % Time step
z = [0:0.05:L]; % Mesh generation
x = [0:0.05:W]; % EDIT: Mesh generation in x-direction
t = [t0:dt:tf];% Time vector
n = numel(z); % Size of mesh grid in z direction
m = numel(x); % EDIT: Size of mess grid in x direction
%Initial Conditions / Vector Creation
c0 = zeros(n,1);
c0(1) = cFeed;
q0 = zeros(m,1); % t = 0, q = 0 for all z, this makes sense to me. EDIT: zeros(m,1) instead of zeros (n,1)
phi0 = zeros(m,1); % EDIT: Electric field vector initial conditions. zeros(m,1) instead of zeros (n,1)
J0 = zeros(m,1); % EDIT: Solid flux vector initial conditions. J = u.q.DphiDt. zeros(m,1) instead of zeros (n,1)
y0 = [c0 ; q0 ; phi0 ; J0]; % Appends conditions together
%ODE15S Solver
[T, Y] = ode15s(@(t,y) EDIfun(t,y,z,x,n),t,y0);
plot(T,Y);
function DyDt=EDIfun(~, y, z, x, n)
% Defining Constants
a = 0.4;
delta = 0.0001;
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
u = 0.1;
Itot = 4;
F = 96485.3329;
% Variables being allocated zero vectors
c = zeros(n,1);
q = zeros(m,1); % EDIT: zeros(m,1) instead of zeros (n,1)
phi = zeros(m,1); % EDIT: Electric field vector
J = zeros(m,1); % EDIT: Solid phase flux vector
DcDt = zeros(n,1);
DqDt = zeros(m,1); % EDIT: zeros(m,1) instead of zeros (n,1)
DphiDt = zeros(m,1); % EDIT: Electric field change with time vector
DJDt = zeros(m,1); % EDIT: Solid phase flux with time vector
DyDt = zeros(4*n,1); % EDIT: DyDt has 4 parts instead of 2
%zhalf = zeros(n-1,1); % EDIT: No need for this line because there is no
%D2CDz2.
DcDz = zeros(n,1);
DphiDx = zeros(m,1); % EDIT: added DphiDx vector
DJDx = zeros(m,1); % EDIT: added DJDx vector
c = y(1:n);
q = y(n+1:2*n);
phi = y(2*n+1:3*n); %EDIT: added phi to y vector
J = y(3*n+1:4*n); %EDIT: added J to y vector
% Interior mesh points
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));
end
%EDIT: added another for loop to discretize in the x direction. I dont know
%how to discretize in the x direction.
J(1) = u*q(1)*DphiDx(1);
DphiDx(1) = Itot/(F*(1-epsilon)*W*L*u*q(1));
DcDz(n) = 0;
DJDx(m) = 0;
DphiDx(m) = 0;
J(m) = 0;
cstar(1) = (q(1)*cHstar)/(K*qH);
for i = 2:m-1
J(i) = u*q(i)*DphiDx(1);
DJDx(i) = (J(i)-J(i-1))/m;
DphiDx(i) = Itot/(F*(1-epsilon)*W*L*u*q(i));
end
% 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
% EDIT: changed the DqDt equation
DqDt(1) = a*D*((c(1)-cstar(1))/delta) - DJDx(1);
DcDt(1) = 0.0;
% Set time derivatives in remaining points
for i=2:n
DqDt(i) = (a*K*(c(i)-cstar(i))/delta) - DJDx;
DcDt(i) = - v*DcDz(i) - (((1-epsilon)/(epsilon))*(DqDt(i)+DJDx(i)));
end
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt;DphiDt;DJDt]; %EDIT: concatenated for all 4 dependent variables
end