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

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 .

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
    end
end  
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 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;
end

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;
end

Answer 1

[darova]

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(:);

[Carlos Lobo]

Thanks a lot!