Hello,
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 
.
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.
(which is a positive value) equal to the initial concentration of species i.
At 
each equation has a Neumann boundary condition:
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^-1global k2r, k2r = 7.06*10^4; %l mol^-1 s^-1K2 = 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/sglobal Dco32, Dco32 = 8.09*10^-10; %m^2/sglobal 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
endend 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
%{
pdegplot(model,'EdgeLabels','on')xlim([0,dl])%}
%Neumann boundary conditions
k=1;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)
applyBoundaryCondition(model,'dirichlet','Edge',2,... 'u',dirichletconds,'EquationIndex', eqindex);%apply Neumann boundary condition to edge 4 (left side)
applyBoundaryCondition(model,'neumann','Edge',4,... 'g',neumannconds_g);ccoef = [Dco2;Dhco3;Dco32;Doh];specifyCoefficients(model, 'm',0,'d',0,'c',ccoef,'a',@acoeffunction,'f', @fcoeffunction);setInitialConditions(model,iniconds);p = generateMesh(model,'Hmax',0.5*10^-6); %0.5micron mesh size
%pdeplot(model);
results = solvepde(model);u = results.NodalSolution;pdeplot(model,'XYData',u(:,1))%since f is not constant, we must define it
function f = fcoeffunction(location,state)global k1fglobal k1rglobal k2fglobal k2rglobal 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;endfunction a = acoeffunction(location,state)global k1fglobal k1rglobal k2fglobal k2rglobal neqs; % Number of equationsnr = length(location.x); % Number of columnsa = 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;end
Best Answer