How to pass solution of pdefun1 as IC-2 for pdefun2
2 次查看(过去 30 天)
显示 更早的评论
my pde is
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/585181/image.png)
Given below is my effort so far...
function pdepe_philip_v3b
% 1-b Extend the above calculation to model a potential step to any chosen
% potential assuming that the ratio of the surface concentrations of Os2+
% and Os3+ obeys the Nernst equation both before and after the potential
% step. Now we want to further add a potential step which will make it a
% double potential step.
global D n F R
n = 1; % No. of Electrons Involved
F = 96485.3329; % sA/mol
R = 8.3145; % kgcm^2/s^2.mol.K
D = 5e-06; % cm^2/s
Area = 1; % cm^2
Os_bulk = 1e-05; % mol/cm^3
N = 10;
m = 0; % Cartesian Co-ordinates
% logarithmized xmesh for semi infinite boundary condition
x = logspace(log10(0.000006), log10(0.06), N); x(1) = 0;
% Time Span-1
t1 = linspace(0, 10, N); % s
E = 0.8;
sol1 = pdepe(m, @pdev3bpde, @(x)pdev3bic(x,Os_bulk), @(xl, ul, xr, ur, t)...
pdev3bbc(xl, ul, xr, ur, t, E, Os_bulk), x, t1);
c1 = sol1(:, :, 1);
I_num = ( n*F*Area) * c1(1,:) .* sqrt(D./(t1*pi));
% I Flux For Time Span-1
figure(1);
plot(t1, I_num, 'r--', 'LineWidth',1);
xlabel('t (s)'); ylabel('I (A/cm^2s)');
hold on
% Time Span-2
t2 = linspace(10, 20, N); % s
c_new = (c1(N,:))';
E = -1.6;
sol2 = pdepe(m, @pdev3bpdeNEW, @(x)pdev3bicNEW(x,c_new),@(xl, ul, xr, ur, t)...
pdev3bbcNEW(xl, ul, xr, ur, t, E, c_new), x, t2);
end
%% pdepe Function
%%
function [a, f, s] = pdev3bpde(x, t, u, DuDx)
global D
a = 1;
f = D*DuDx;
s = 0;
end
%% pdepe Function
%%
function [a, f, s] = pdev3bpdeNEW(x, t2, u, DuDx)
global D
a = 1;
f = D*DuDx;
s = 0;
end
%% Initial Condition-1
%%
function u0 = pdev3bic(x, Os_bulk)
u0 = Os_bulk;
end
%% Initial Condition-2
%%
function u0 = pdev3bicNEW(x, c_new)
u0 = c_new.*size(x);
end
%% Boundry Condition-1
%%
function [pl, ql, pr, qr] = pdev3bbc(xl, ul, xr, ur, t, E, Os_bulk)
global n F R
E0 = 0.2;
alpha = exp((E-E0).*((n*F)/(R*298.15)));
pl = ul - (Os_bulk./(1+alpha));
ql = 0;
pr = ur - Os_bulk;
qr = 0;
end
%% Boundry Condition-2
%%
function [pl, ql, pr, qr] = pdev3bbcNEW(xl, ul, xr, ur, t, E, Os_bulk)
global n F R
E0 = 0.2;
alpha = exp((E-E0).*((n*F)/(R*298.15)));
pl = ul - (Os_bulk./(1+alpha));
ql = 0;
pr = ur - Os_bulk;
qr = 0;
end
%%
%%
采纳的回答
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!