How to input a novel boundary condition for a coupled PDE system
112 次查看(过去 30 天)
显示 更早的评论
Hi all,
I am trying to simulate a coupled PDE system with a non typical boundary conditions using the pde1dm solver which is based on the pdepe solver in MATLAB.
My coupled PDE system is as such.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752574/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752579/image.png)
The system is solved for two time periods
For first time period ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752739/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752739/image.png)
For s at the L.H.S we have
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752584/image.png)
For s at the R.H.S we have
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752589/image.png)
For species m at the L.H.S we have the to set m value to a constant
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752594/image.png)
For species m at the R.H.S we have the no flux condition
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752599/image.png)
For second time period ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752744/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752744/image.png)
For s at the L.H.S and R.H.S we have the same such as
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752584/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752589/image.png)
Also for the m R.H.S we have the same boundary condition.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752599/image.png)
But, for the species
we have the L.H.S boundary condition for timr period
. It is this boundary condtiuon that I am seeking help with. How can I go about implementing this on the L.H.S for
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752614/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752749/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752749/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1752944/image.png)
function [c25, c2, vode] = pde1dm_PS_OCP_v1(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_r)
% Piecewise PDE Implementation
% Pulsing Simulation
% Constants Used in the Nernst Equation
n = 2; % No. of Electrons Involved
F = 96485.3329; % C/mol
R = 8.314; % J/(mol.K) or CV/(mol.K)
Temp = 273.15+20; % K
N = 21;
m = 0;
Area = 0.0011; % Geometric Area (cm^2)
m_layer = 4.13E-05; % m_layer (mol/cm^3)
D_m = 8.14E-08; % cm^2/s (De of polymer)
k = n*F*Area*l*m_layer;
% DIMENSIONLESS PARAMETERS
chi = linspace(0, 1, N);
tau = linspace(0, tau_M1, N); % Dimensionless Time
% FIRST POTENTIAL STEP
E0 = 0.22;
E_set = .350;
epsilon1 = ((n*F)/(R*Temp))*(E_set-E0); % Dimensionless Potential
c_M_ox = 1/(1+exp(-epsilon1)); % Mox BC
ic_arg_1 = {@(x)ones(size(N)) ; @(x)zeros(size(N))};
IC = @(x)PDE_ND_PS_EK_IC(x, ic_arg_1);
BC = @(xl, yl, xr, yr, t)PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox);
sol1 = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r), ...
IC, BC, chi, tau);
c1 = sol1(:, :, 1); % Substrate Conc.
c2 = sol1(:, :, 2); % Mox Conc.
% OPEN CIRCUIT POTENTIAL
tau2 = linspace(tau_M1, tau_M2, N); % Dimensionless Time
ic_arg_1 = {@(x)interp1(chi, sol1(N, :, 1), x) ; @(x)interp1(chi, sol1(N, :, 2), x)};
IC2 = @(x)PDE_ND_PS_EK_IC(x, ic_arg_1);
BC2 = @(xl, yl, xr, yr, t, v)PDE_PSw_EK_BC_2(xl, yl, xr, yr, v);
xOde = [1, .9].';
ode_IC = @() ode_ICFun(tau_M1);
opts.vectorized='on';
% sol2 = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), ...
% IC2, BC2, chi, tau2);
[sol2, vode] = pde1dm(m, @(x, t, c, DcDx)PDE_ND_PS_EK(x, t, c, DcDx, kappa, eta, gamma, mu, D_r), ...
IC2, BC2, chi, tau2, @odeFunc, ode_IC, xOde(2));
% Concentration Profiles c(i, j, k)(Solutions)
c14 = [c1; sol2(:, :, 1)]; % Substrate Conc.
c25 = [c2; sol2(:, :, 2)]; % Mox Conc.
c36 = 1-c25; % Mred Conc.
mox = abs(sol2(:, 1, 2));
mred = 1-mox;
E_PS = E_set.*ones(1,N);
for counter2 = 1:N
E_OCP(counter2) = E0 + (((R*Temp)/(n*F).*(log(mox(counter2,:)./mred(counter2,:)))));
end
E_array = [E_PS, E_OCP];
figure(2);
plot(chi, c25(N,:));
xlabel('\chi'); ylabel('m_{ox}');
function [cc, ff, ss] = PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r)
% S; Mox;
cc = [ D_r; 1];
ff = [ 1; 1].*DcDx;
S_kin = ((gamma/eta)*(kappa^2)*c(1)*c(2))./...
(gamma.*c(2).*(1+(mu.*c(1)))+c(1));
M_kin = ((kappa^2).*c(1).*c(2))./...
(gamma.*c(2).*(1+(mu.*c(1)))+c(1));
ss = [-S_kin; -M_kin];
end
function c0 = PDE_ND_PS_EK_IC(x, ic_arg_1)
% Substrate; Mox;
c0 = [ic_arg_1{1}(x).'; ic_arg_1{2}(x).'];
end
% % % % % % % % % % % % % % % % For 0<tau_m<tau_1
function [pl, ql, pr, qr] = PDE_PS_EK_BC(xl, yl, xr, yr, c_M_ox)
% ---------------------
% | | |
% | | |
% | | |
% | | |
% | | |
% | | |
% | | |
% -------|-------------
% pl pr
% Substrate; Mediator;
pl = [0 ; yl(2)-c_M_ox];
ql = [1 ; 0];
pr = [yr(1)-1 ; 0];
qr = [0 ; 1];
end
% % % % % % % % % % % % % % % % For tau_1<tau_m<tau_2
function [pl, ql, pr, qr] = PDE_PSw_EK_BC_2(xl, yl, xr, yr, t, v, vdot)
pl = [0 ; 0];
ql = [1 ; 1];
pr = [yr(1)-1 ; 0];
qr = [0 ; 1];
end
function f = odeFunc(x, t, c, DcDx, v, vdot)
f = (DcDx*(1E-03./(v(1)*(1-v(1)))))-vdot(1);
end
function v0=ode_ICFun(tau_0)
v0 = ones(tau_0);
end
end
0 个评论
采纳的回答
Torsten
2024-8-14,23:44
xOde = [0 1].'; % Define coupling points for ODE definition
function [pl, ql, pr, qr] = PDE_PSw_EK_BC_2(xl, yl, xr, yr, t, v, vdot)
...
pl(2) = yl(2) - v;
ql(2) = 0;
...
end
function f = odeFunc(x, t, c, DcDx, v, vdot)
% Assuming first index in DcDx and c is equation, second index is coupling point in xOde vector
% (might be vice versa, look up in the User's Guide of pde1dm)
f = DcDx(2,1) - k/(c(2,2)*(1-c(2,2))*vdot
end
function v0=ode_ICFun(tau_0)
v0 = m_ox %@t=0,chi=0
end
18 个评论
Torsten
2024-8-15,22:12
编辑:Torsten
2024-8-15,22:12
The reason for your problems is obvious: c2(N,1) = 1.
That means that you divide by 0 in the expression f=DcDx(2) - 8/(c(2)*(1-c(2)))*vdot(1) .
In my opinion, there is something wrong in the mathematical formulation of your problem or in the choice of the starting values of c2 from phase 1 in phase 2.
更多回答(1 个)
Bill Greene
2024-8-20,13:30
编辑:Bill Greene
2024-8-20,13:53
Here is my code to solve this problem with pde1dm based on my understanding of the description and included code. I'm not at all sure I understand the problem definition and since I know nothing about the physics of the problem, I have no idea whether this solution is actually correct.
function matlabAnswers_8_15_2024
kappa=1; eta=1; gamma=1; mu=1; tau_M1=1; tau_M2=2; l=1; D_r=1;
pde1dm_PS_OCP_v1(kappa, eta, gamma, mu, tau_M1, tau_M2, 1, 1);
end
function [c25, c14, mox, vode] = pde1dm_PS_OCP_v1(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_r)
% Piecewise PDE Implementation
% Pulsing Simulation
% Constants Used in the Nernst Equation
n = 2; % No. of Electrons Involved
F = 96485.3329; % C/mol
R = 8.314; % J/(mol.K) or CV/(mol.K)
Temp = 273.15+20; % K
Nx = 21;
Nt=15;
m = 0;
Area = 0.0011; % Geometric Area (cm^2)
m_layer = 4.13E-05; % m_layer (mol/cm^3)
D_m = 8.14E-08; % cm^2/s (De of polymer)
k = n*F*Area*l*m_layer;
% DIMENSIONLESS PARAMETERS
chi = linspace(0, 1, Nx);
tau1 = linspace(0, tau_M1, Nt); % Dimensionless Time
tau2 = [linspace(tau_M1+1e5*eps, tau_M2, Nt)];
% FIRST POTENTIAL STEP
E0 = 0.22;
E_set = .350;
epsilon1 = ((n*F)/(R*Temp))*(E_set-E0); % Dimensionless Potential
c_M_ox = 1/(1+exp(-epsilon1)); % Mox BC
IC = @(x) PDE_ND_PS_EK_IC(x);
BC = @(xl, yl, xr, yr, t,v,vDot) PDE_PS_EK_BC(xl, yl, xr, yr, t,v,vDot, c_M_ox, tau_M1, k);
pdef= @(x, t, c, DcDx) PDE_ND_PS_EK(chi, tau1, c, DcDx, kappa, eta, gamma, mu, D_r);
xOde=0;
odeIcF= @() ode_ICFun(IC, chi);
[sol1,sol1Ode] = pde1dm(m, pdef, IC, BC, chi, tau1, @odeFunc, odeIcF, xOde);
% figure; plot(tau1, sol1Ode, tau1, sol1(:,1,2));
% title 'mox at left end as a function of time'; grid;
IC2=@(x) icFun2(x, chi, sol1);
odeIcF2= @() ode_ICFun(IC2, chi);
[sol2,sol2Ode] = pde1dm(m, pdef, IC2, BC, chi, tau2, @odeFunc, odeIcF2, xOde);
sol=[sol1; sol2];
solOde=[sol1Ode; sol2Ode];
tau=[tau1 tau2];
s = sol(:, :, 1); % Substrate Conc.
mox = sol(:, :, 2); % Mox Conc.
figure; plot(tau, mox(:,1));
title 'mox at left end as a function of time'; grid;
figure; plot(tau, mox(:,end)); grid;
title 'mox at right end as a function of time';
figure; plot(chi, s(end,:));
title 's at final time'; grid;
figure; plot(chi, mox(end,:));
title 'mox at final time'; grid;
end
function [cc, ff, ss] = PDE_ND_PS_EK(chi, tau, c, DcDx, kappa, eta, gamma, mu, D_r)
% S; Mox;
cc = [ D_r; 1];
ff = [ 1; 1].*DcDx;
S_kin = ((gamma/eta)*(kappa^2)*c(1)*c(2))./...
(gamma.*c(2).*(1+(mu.*c(1)))+c(1));
M_kin = ((kappa^2).*c(1).*c(2))./...
(gamma.*c(2).*(1+(mu.*c(1)))+c(1));
ss = [-S_kin; -M_kin];
end
function c0 = PDE_ND_PS_EK_IC(x)
% Substrate; Mox;
c0 = [1 0]';
end
function c0 = icFun2(x, chi, sol1)
s1Final=squeeze(sol1(end,:,:));
c0=interp1(chi, s1Final, x)';
end
function [pl, ql, pr, qr] = PDE_PS_EK_BC(xl, yl, xr, yr, t,v,vDot, c_M_ox, tau_M1, k)
% ---------------------
% | | |
% | | |
% | | |
% | | |
% | | |
% | | |
% | | |
% -------|-------------
% pl pr
% Substrate; Mediator;
if t <= tau_M1
pl = [0 ; yl(2)-c_M_ox];
ql = [1 ; 0];
else
moxR=yr(2);
dMoxDTau=vDot;
p=-k/(moxR*(1-moxR))*dMoxDTau;
pl = [0 ; p];
ql = [1 ; 1];
end
pr = [yr(1)-1 ; 0];
qr = [0 ; 1];
end
function F=odeFunc(t,v,vdot,x,u)
F=v-u(2);
end
function v0=ode_ICFun(icf,chi)
x0=chi(1);
c0=icf(x0);
v0 = c0(2);
end
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1756354/image.png)
1 个评论
Torsten
2024-8-20,13:56
OP changed his problem description from
moxR=yr(2);
to
moxL=yl(2);
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!