# How to code Flux Boundary Conditions for a PDE using ODE15s and method of lines

16 views (last 30 days)
Nicholas Vecchiarello on 21 Sep 2022
Answered: Torsten on 21 Sep 2022
I have a PDE in spherical coordinates that I've discretized spatially using the method of lines for solving with ODE15s. I am having trouble figuring out how to represent the flux boundary conditions:
1.) At r=0, ∂c/∂r = 0
2.) at r=R, De ∂c/∂r = kf (C − c)
My code looks like this: Thanks for your help!
function [dy] = Kinetics(t,y)
kf=1.8*10^-3; %Mass transfer coefficient (cm/s)
De=8*10^-8;
R=20*10^-4;
ka=2*10^-4;
Keq=72;
Vm=40;
V=40;
n=100;
qm=45;
Ep=0.6;
cp=y(1:n);
Qp=y(n+1:2*n);
Cp=y(2*n+1);
dCp=zeros(1,1);
dcp=zeros(n,1);
dQp=zeros(n,1);
Vl_stage=(1/n);
%%At r=0 Boundary Conditions
dcp(1)=(cp(1)-cp(2))/Vl_stage; %This is where I tried to implement Boundary Condition at r=0
dQp(1)=((ka)*(((qm-Qp(1))*cp(1))-(Qp(1)/Keq)));
for i=2:n-1
Ri=(Vl_stage*(i)/R);
Ri1=(Vl_stage*(i+1)/R);
dcp(i)=(De/Ep)*(1/Ri1^2)*(((Ri1^2*(cp(i+1)-cp(i)))-(Ri^2*(cp(i)-cp(i-1))))/(Vl_stage^2))-(dQp(i)/Ep);
dQp(i)=((ka)*(((qm-Qp(i))*cp(i))-(Qp(i)/Keq)));
end
%%At r=R Boundar Conditions
for i=n
dcp(n)=(kf/De)*(Cp(1)-cp(n))/Vl_stage; %This is where I tried to implement Boundary Condition at r=R
dQp(n)=((ka/Ep)*(((qm-Qp(n))*cp(n))-(Qp(n)/Keq)));
dCp(1)=-(3*Vm/(R*V))*kf*(Cp(1)-cp(n));
end
dy=[dcp;dQp;dCp];
end
And I am calling with [t,y] = ode15s(@Kinetics,[0,1],[horzcat(zeros(1,200),10)]);
If it helps, here is the original system of equations:

Torsten on 21 Sep 2022
One of many possibilities:
dc/dr = 0
c2 - c1 = 0
->
dc1/dt = dc2/dt
De*dc/dr = kf*(C-c)
De*(cn-c(n-1))/Vl_stage = kf*(C-cn)
cn = (kf*C + De*c(n-1)/Vl_stage)/(kf+De/Vl_stage)
->
dcn/dt = (kf*(dC/dt) + De*(dc(n-1)/dt)/Vl_stage)/(kf+De/Vl_stage)

### Categories

Find more on Boundary Conditions in Help Center and File Exchange

R2022b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!