MATLAB: Solving a system of PDEs with non constant coefficients and dirichlet and neumann boundary conditions


I am trying to solve the following system of equations which deal with the diffusion of chemical species within an electrical double layer of thickness δ:
The diffusion coefficients, with are constant values as well as the reaction rate constants .
The concentration of the species are all dependent on the spacial x coordinate:
with .
For each equation we have one dirichlet condition at (which is a positive value) equal to the initial concentration of species i.
At each equation has a Neumann boundary condition:
I have rewritten the equations in the form:
which is the form that they need to be so I can use the PDE toolbox. I have managed to code the f matrix but I am unclear as to how to write the a matrix, which is a diagonal matrix (see attached images – sorry, couldn't insert them directly here…)
The code runs without issues, but the concentration values do not seem to change along x. I have set the initial conditions for the problem as well, equal to the initial concentrations of the species. I am not sure whether I have set the $a$ matrix correctly, but I believe that the f matrix is correct, since I followed the documentation. I tried to do the same with the a matrix, but I was a bit unclear as to how to proceed with it, so I just adapted the example that is given.
Is it possible to solve this kind of system using bvp4c (or bvp5c)? from what I understood of the documentation, these methods apply only to equations with boundary values that do not use derivatives [g(y(a),y(b))=0], right? is there a way to use Neumann boundary conditions with these problems?
Any help will be greatly appreciated!
Thank you all!
I am running the following code:
%Physical constants
R = 8.3145; %J mol^-1 K^-1
T = 298; %K
P = 100000; %Pa (1 bar)
F = 96485.3329; %s A mol^-1 - Faraday's constant
dl = 50*10^-6; %delta [m] - thickness of double layer
A = 10^-4; %m^2 - area of electrode
V = 25*10^-6; %m^3 - volume of cell
z = 2; %electrons flowing
j = [1 5 10 50].*10; %current densities in A m^2
%Reaction constants
kco2 = 0.84;
Kw = 1.00*10^-14; %mol^2 l^-2
global k1f, k1f = 2.23*10^3; %l mol^-1 s^-1

global k1r, k1r = 1.7*10^-5; %s^-1
K1 = k1f/k1r;
global k2f, k2f = 6.00*10^9; %l mol^-1 s^-1
global k2r, k2r = 7.06*10^4; %l mol^-1 s^-1
K2 = k2f/k2r;
ka = 1.8*10^-4; %l mol^-1
kplus = 1000; %mol m^-3 (1 M) - supporting electrolyte
%Diffusion constants
global Dco2, Dco2 = 1.67*10^-9; %m^2/s

global Dhco3, Dhco3 = 1.04*10^-9; %m^2/s
global Dco32, Dco32 = 8.09*10^-10; %m^2/s
global Doh, Doh = 5.27*10^-9; %m^2/s

%Determining the initial concentrations of the species
%CO2 concentration in the gas bubbles
co2g = P/(R*T); %mol m^-3
co2ini = co2g*kco2; %mol m^-3 -- Initial concentration of CO2
%coefficients of the [OH] equation after rearranging
ohcoef = [2*K1*K2*co2ini K1*co2ini+1 -kplus -Kw];
r = roots(ohcoef);
%pick correct root
for i = 1:length(r)
if (imag(r(i))==0)&&(real(r(i))>0)
ohini = r(i); %mol m^-3 -- Initial concentration of OH
co32ini = K1*K2*ohini^2*co2ini; %mol m^-3 -- Initial concentration of CO32-
hco3ini = K1*ohini*co2ini; %mol m^-3 -- Initial concentration of HCO3-
inivals = ["[CO2]" "[OH]" "[H]" "[CO32]" "[HCO3]";co2ini ohini hini co32ini hco3ini];%for visualization purpose

%Pseudo steady state concentration profiles in the boundary layer
global neqs, neqs = 4;
model = createpde(neqs); %4 equations in the system
height = 5*10^-6; %m = 5 microns
x1 = 0; y1 = -height/2;
x2 = dl; y2 = -height/2;
x3 = dl; y3 = height/2;
x4 = 0; y4 = height/2;
Rect = [3;4;x1;x2;x3;x4;y1;y2;y3;y4];
g = decsg(Rect);
geom = geometryFromEdges(model,g); %Creting the geometry
%Neumann boundary conditions
co2neug = j(k)/(z*F);
ohneug = j(k)/(z*F);
co32neug = 0;
hco3neug = 0;
iniconds = [co2ini; hco3ini; co32ini; ohini];
eqindex = [1 2 3 4];
dirichletconds = [co2ini, hco3ini, co32ini, ohini];
neumannconds_g = [co2neug, hco3neug, co32neug, ohneug];
%apply Dirichlet boundary condition to edge 2 (right side)
'u',dirichletconds,'EquationIndex', eqindex);
%apply Neumann boundary condition to edge 4 (left side)
ccoef = [Dco2;Dhco3;Dco32;Doh];
specifyCoefficients(model, 'm',0,'d',0,'c',ccoef,'a',@acoeffunction,'f', @fcoeffunction);
p = generateMesh(model,'Hmax',0.5*10^-6); %0.5micron mesh size
results = solvepde(model);
u = results.NodalSolution;
%since f is not constant, we must define it
function f = fcoeffunction(location,state)
global k1f
global k1r
global k2f
global k2r
global neqs; % Number of equations

nr = length(location.x); % Number of columns

f = zeros(neqs,nr); % Allocate f
% Now the particular functional form of f
f(1,:) = state.u(2,:).*k1r;
f(2,:) = state.u(3,:).*k2r + state.u(1,:).*state.u(4,:).*k1f;
f(3,:) = state.u(2,:).*state.u(4,:).*k2f;
f(4,:) = state.u(2,:).*k1r + state.u(3,:).*k2r;
function a = acoeffunction(location,state)
global k1f
global k1r
global k2f
global k2r
global neqs; % Number of equations
nr = length(location.x); % Number of columns
a = zeros(neqs,nr); % Allocate a
% Now the particular functional form of a
a(1,:) = state.u(4,:).*k1f;
a(2,:) = state.u(4,:).*k2f+k1r;
a(3,:) = k2r;
a(4,:) = state.u(1,:).*k1f + state.u(2,:).*k2f;

Best Answer

If your equations change in x direction only they are ODE (ordinary differential equations)
  • Is it possible to solve this kind of system using bvp4c (or bvp5c)?
So yes, it's possible to solve them using bvp4c
function du = f(t,u)
CCO2 = u(1);
CHCO3 = u(2);
CCO3 = u(3);
COH = u(4);
du(1:4,1) = u(5:8);
du(5,1) = eq1
du(5,1) = eq2
% ...
boundary conditions can be written as:
function res = fb(ua,ub)
% concentration of species at initial moments
res(1) = ub(1) - c1;
res(2) = ub(1) - c2;

% derivetives
res(5) = ua(2) - ...
res = res(:);