- Is it possible to solve this kind of system using bvp4c (or bvp5c)?
Solving a system of PDEs with non constant coefficients and dirichlet and neumann boundary conditions
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
0 commentaires
Réponse acceptée
darova
le 19 Mai 2020
If your equations change in x direction only they are ODE (ordinary differential equations)
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(:);
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Boundary Conditions dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!