1
$\begingroup$

I'm working on solving an equation to estimate production in a spherical catalyst pellet. A working example is given below:

delta = 0.01;
k=1;
t0 = 260 + 273;
rt = 8.314;
presList = {p[x], p2[x], p3[x], p4[x]};
presBCList = 80*{9/100, 75/100, 0.17/100, 0.07/100};
{a1, a2, a3, a4, a5} = {4224.40, 0.00001446, 3453.38, 24.18, 96.629};

p2[x_] = presBCList[[2]] - 3 (presBCList[[1]] - p[x]);
p3[x_] = presBCList[[3]] + 1 (presBCList[[1]] - p[x]);
p4[x_] = presBCList[[4]] + 1 (presBCList[[1]] - p[x]);

The reaction speed is given by

r[x_] = k ((a1*p[x]*
   p2[x]*(1 - 1/a2*(p3[x]*p4[x])/(p2[x]^3*p[x])))/(1 + 
   a3*p4[x]/p2[x] + a4*Sqrt[p2[x]] + a5*p4[x])^3);

Next setting up NDSolve with no flux at pellet center BC and a given flux at the surface (x=1):

 diffEqList = {(Derivative[2][p][x] + 2/x Derivative[1][p][x]) - 
 r[x] == 0};
 bcList = {Derivative[1][p][delta] == 0, 
 Derivative[1][p][1] == -(p[1] - presBCList[[1]])};
 sol = NDSolve[Join[diffEqList, bcList], p[x], {x, delta, 1}];

The equation is immediatly solved when k=1 to k=11. However, when k>= 11.95, NDSolve fails.

Power::infy: Infinite expression 1/0. encountered.

Yet looking at the graphs, I can not see what quantity may be going to 0. As far as I can tell, the solutions should not be qualitatively different when k is increased. Changing the value of delta to increase it seems to help a little, but not much. Can this equation be solved with a better method ?

$\endgroup$
  • 1
    $\begingroup$ As well as Chris's answer, are you expanding the spherical coordinates around x=0 to get the boundary condition there, or just moving it delta away from 0? i.e. r'(delta) = r'(0) + delta r''(0) + ... $\endgroup$ – KraZug May 15 '17 at 12:03
  • $\begingroup$ Just moved the delta away from 0. Would that help ? $\endgroup$ – Whelp May 15 '17 at 12:34
  • $\begingroup$ It is probably not essential for this set of equations, but it is something to consider if you use higher order differential equations. $\endgroup$ – KraZug May 15 '17 at 12:39
  • $\begingroup$ I have seen in the literature that this problem can be solved using the orthogonal collocation method. Is this method implemented in MMA ? $\endgroup$ – Whelp May 15 '17 at 13:30
4
$\begingroup$

NDSolve is using the shooting method for this boundary value problem. To help get it on the right track, you can manually give initial guesses. As noted by @Szabolcs in this answer "this is documented in the Advanced Numerical Differential Equation Solving tutorial".

This seems to help with your problem:

k = 15;
r[x_] = k ((a1*p[x]*p2[x]*(1 - 1/a2*(p3[x]*p4[x])/(p2[x]^3*p[x])))/(1 + a3*p4[x]/p2[x] + a4*Sqrt[p2[x]] + a5*p4[x])^3);

diffEqList = {(Derivative[2][p][x] + 2/x Derivative[1][p][x]) - r[x] == 0};
bcList = {Derivative[1][p][delta] == 0, Derivative[1][p][1] == -(p[1] - presBCList[[1]])};
sol = NDSolve[Join[diffEqList, bcList], p, {x, delta, 1},
  Method -> {"Shooting", "StartingInitialConditions" -> {p[delta] == 6, p'[delta] == 0}}];

works fine for me.

$\endgroup$
  • $\begingroup$ This seems to help a lot indeed. A basic question but how do you get values or plot p'[x] ? I'm not able to find the right combination of commands to do it. Plot[Evaluate[p'[x] /. sol], {x, delta, 1}] does not seem to work. $\endgroup$ – Whelp May 15 '17 at 12:57
  • $\begingroup$ Better to make the unknown in NDSolve p, not p[x]. I'll edit my answer to reflect this. Then Plot[p'[x] /. sol[[1]], {x, delta, 1}] works. $\endgroup$ – Chris K May 15 '17 at 13:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.