How to input a novel boundary condition for a coupled PDE system
Hi all,
I am trying to simulate a coupled PDE system with a non typical boundary conditions.
My coupled PDE system is as such.
The system is solved for two time periods
For first time period 0<tau_M1<tau
For s at the L.H.S we have
For s at the R.H.S we have
For species m at the L.H.S we have the to set m value to a constant
For species m at the R.H.S we have the no flux condition
For second time period tau<tau_M1<tau_2
For s at the L.H.S and R.H.S we have the same such as
For m we have
For the species we have the L.H.S boundary condition for timr period (tau<tau_M1<tau_2). It is this boundary condtiuon that I am seeking help with. How can I go about this.
function [c25, c2, vode] = pde1dm_PS_OCP_v1()
% 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;
kappa=1;
eta=1;
gamma=1;
mu=1;
tau_M1=1;
tau_M2=2;
l=1;
D_r=1;
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 = .450;
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);
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, .99);
% 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];
% This is the case finder function
% Video(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_m, N, chi, c14, c25, c36, E_array)
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
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
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 = ((1E-03./(v(1)*(1-v(1))))*(DcDx))-vdot(1);
end
function v0=ode_ICFun(tau_0)
v0 = ones(tau_0).*.9;
end
function myVideo = Video(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_m, N, chi, c14, c25, c36, E_array)
% Initialize Video
G = figure(1);
myVideo = VideoWriter(sprintf(‘κ%.2f γ%.2f η%.2f µ%.2f’, kappa, gamma, eta, mu));
myVideo.FrameRate = 10; % can adjust this, 5 – 10 works well for me
myVideo.Quality = 100;
open(myVideo);
color = ‘red’;
u = uicontrol(G, ‘Style’,’slider’,’Position’,[0 40 10 360],…
‘Min’,0,’Max’,N*2,’Value’,0);
tspan1 = linspace(0, round(((tau_M1*(l^2))/D_m)), N);
tspan2 = linspace(round(((tau_M1*(l^2))/D_m)), round(((tau_M2*l^2)/D_m),2), N);
tspan = [tspan1,tspan2];
for ii = 2:(N*2)
% Plot Species Conc. in the Layer
subplot(2,1,1);
yyaxis left
plot(chi,c25(ii,:),chi,c36(ii,:), ‘LineWidth’, 2);
ylabel(‘M’)
ylim([0,1]);
yyaxis right
plot(chi, c14(ii,:), ‘LineWidth’, 2);
ylabel(‘S’)
ylim([0,1]);
xlabel(‘chi’);
legend(‘M_{ox}’,’M_{red}’, ‘S’);
title(‘Normalized Concentration’);
subplot(2,1,2);
addpoints(animatedline(tspan,E_array,’marker’, ‘.’, ‘markersize’, 6, ‘color’, color,…
‘linestyle’, ‘–‘, ‘MaximumNumPoints’, 1),tspan(ii),E_array(ii));
ylim([0,0.5]);
xlim([0, tspan(end)]);
hold on;
title([‘Potential = ‘, num2str(E_array(ii))]);
% text(8, 8,{”,”})
drawnow
xlabel(‘t / s’); ylabel(‘E / V’);
u.Value = ii;
uicontrol(‘Style’,’Edit’,’Position’,[0,00,40,40], …
‘String’,num2str(tspan(ii),3));
pause(0.001)
M = getframe(G);
writeVideo(myVideo, M);
end
end
endHi all,
I am trying to simulate a coupled PDE system with a non typical boundary conditions.
My coupled PDE system is as such.
The system is solved for two time periods
For first time period 0<tau_M1<tau
For s at the L.H.S we have
For s at the R.H.S we have
For species m at the L.H.S we have the to set m value to a constant
For species m at the R.H.S we have the no flux condition
For second time period tau<tau_M1<tau_2
For s at the L.H.S and R.H.S we have the same such as
For m we have
For the species we have the L.H.S boundary condition for timr period (tau<tau_M1<tau_2). It is this boundary condtiuon that I am seeking help with. How can I go about this.
function [c25, c2, vode] = pde1dm_PS_OCP_v1()
% 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;
kappa=1;
eta=1;
gamma=1;
mu=1;
tau_M1=1;
tau_M2=2;
l=1;
D_r=1;
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 = .450;
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);
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, .99);
% 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];
% This is the case finder function
% Video(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_m, N, chi, c14, c25, c36, E_array)
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
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
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 = ((1E-03./(v(1)*(1-v(1))))*(DcDx))-vdot(1);
end
function v0=ode_ICFun(tau_0)
v0 = ones(tau_0).*.9;
end
function myVideo = Video(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_m, N, chi, c14, c25, c36, E_array)
% Initialize Video
G = figure(1);
myVideo = VideoWriter(sprintf(‘κ%.2f γ%.2f η%.2f µ%.2f’, kappa, gamma, eta, mu));
myVideo.FrameRate = 10; % can adjust this, 5 – 10 works well for me
myVideo.Quality = 100;
open(myVideo);
color = ‘red’;
u = uicontrol(G, ‘Style’,’slider’,’Position’,[0 40 10 360],…
‘Min’,0,’Max’,N*2,’Value’,0);
tspan1 = linspace(0, round(((tau_M1*(l^2))/D_m)), N);
tspan2 = linspace(round(((tau_M1*(l^2))/D_m)), round(((tau_M2*l^2)/D_m),2), N);
tspan = [tspan1,tspan2];
for ii = 2:(N*2)
% Plot Species Conc. in the Layer
subplot(2,1,1);
yyaxis left
plot(chi,c25(ii,:),chi,c36(ii,:), ‘LineWidth’, 2);
ylabel(‘M’)
ylim([0,1]);
yyaxis right
plot(chi, c14(ii,:), ‘LineWidth’, 2);
ylabel(‘S’)
ylim([0,1]);
xlabel(‘chi’);
legend(‘M_{ox}’,’M_{red}’, ‘S’);
title(‘Normalized Concentration’);
subplot(2,1,2);
addpoints(animatedline(tspan,E_array,’marker’, ‘.’, ‘markersize’, 6, ‘color’, color,…
‘linestyle’, ‘–‘, ‘MaximumNumPoints’, 1),tspan(ii),E_array(ii));
ylim([0,0.5]);
xlim([0, tspan(end)]);
hold on;
title([‘Potential = ‘, num2str(E_array(ii))]);
% text(8, 8,{”,”})
drawnow
xlabel(‘t / s’); ylabel(‘E / V’);
u.Value = ii;
uicontrol(‘Style’,’Edit’,’Position’,[0,00,40,40], …
‘String’,num2str(tspan(ii),3));
pause(0.001)
M = getframe(G);
writeVideo(myVideo, M);
end
end
end Hi all,
I am trying to simulate a coupled PDE system with a non typical boundary conditions.
My coupled PDE system is as such.
The system is solved for two time periods
For first time period 0<tau_M1<tau
For s at the L.H.S we have
For s at the R.H.S we have
For species m at the L.H.S we have the to set m value to a constant
For species m at the R.H.S we have the no flux condition
For second time period tau<tau_M1<tau_2
For s at the L.H.S and R.H.S we have the same such as
For m we have
For the species we have the L.H.S boundary condition for timr period (tau<tau_M1<tau_2). It is this boundary condtiuon that I am seeking help with. How can I go about this.
function [c25, c2, vode] = pde1dm_PS_OCP_v1()
% 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;
kappa=1;
eta=1;
gamma=1;
mu=1;
tau_M1=1;
tau_M2=2;
l=1;
D_r=1;
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 = .450;
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);
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, .99);
% 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];
% This is the case finder function
% Video(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_m, N, chi, c14, c25, c36, E_array)
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
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
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 = ((1E-03./(v(1)*(1-v(1))))*(DcDx))-vdot(1);
end
function v0=ode_ICFun(tau_0)
v0 = ones(tau_0).*.9;
end
function myVideo = Video(kappa, eta, gamma, mu, tau_M1, tau_M2, l, D_m, N, chi, c14, c25, c36, E_array)
% Initialize Video
G = figure(1);
myVideo = VideoWriter(sprintf(‘κ%.2f γ%.2f η%.2f µ%.2f’, kappa, gamma, eta, mu));
myVideo.FrameRate = 10; % can adjust this, 5 – 10 works well for me
myVideo.Quality = 100;
open(myVideo);
color = ‘red’;
u = uicontrol(G, ‘Style’,’slider’,’Position’,[0 40 10 360],…
‘Min’,0,’Max’,N*2,’Value’,0);
tspan1 = linspace(0, round(((tau_M1*(l^2))/D_m)), N);
tspan2 = linspace(round(((tau_M1*(l^2))/D_m)), round(((tau_M2*l^2)/D_m),2), N);
tspan = [tspan1,tspan2];
for ii = 2:(N*2)
% Plot Species Conc. in the Layer
subplot(2,1,1);
yyaxis left
plot(chi,c25(ii,:),chi,c36(ii,:), ‘LineWidth’, 2);
ylabel(‘M’)
ylim([0,1]);
yyaxis right
plot(chi, c14(ii,:), ‘LineWidth’, 2);
ylabel(‘S’)
ylim([0,1]);
xlabel(‘chi’);
legend(‘M_{ox}’,’M_{red}’, ‘S’);
title(‘Normalized Concentration’);
subplot(2,1,2);
addpoints(animatedline(tspan,E_array,’marker’, ‘.’, ‘markersize’, 6, ‘color’, color,…
‘linestyle’, ‘–‘, ‘MaximumNumPoints’, 1),tspan(ii),E_array(ii));
ylim([0,0.5]);
xlim([0, tspan(end)]);
hold on;
title([‘Potential = ‘, num2str(E_array(ii))]);
% text(8, 8,{”,”})
drawnow
xlabel(‘t / s’); ylabel(‘E / V’);
u.Value = ii;
uicontrol(‘Style’,’Edit’,’Position’,[0,00,40,40], …
‘String’,num2str(tspan(ii),3));
pause(0.001)
M = getframe(G);
writeVideo(myVideo, M);
end
end
end pde MATLAB Answers — New Questions
