Dear Experts,

Please let me, in advance, apologies for my lengthy introduction to the problem. I also apologize if the code formatting is horrible. If anyone has any suggestions on how I can improve it, please let me know.

I have previously, with success, used Mathematica to model the reactant concentration profile inside catalyst pellets of different geometry, and to check for internal mass transport limitations. Previously the reactions have involved solely one reactant. I am now trying to apply the same set of equations to the oxidation reaction of orthoxylene to phthalic acid anhydride.

Initially (just to get going) I tried to express the concentration of Oxygen inside the catalyst pellet, as a function of the O-Xylene concentration, using stoichiometry. This might not be a good approach, but that is a problem for later. In effect, I ended up with a single second order differential equation. The rate of reaction I have defined as:

RateOfReaction[T_, Coxylen_] :=
 Module[{temp , kr, ka, cxyl, cO2},
  temp = T;
  cxyl = Coxylen;
  kr = 8.84*10^3*Exp[-(8170/temp)]*0.001;
  ka = 199*10^3*Exp[-(11880/temp)]*0.001;
  cO2 = CXylenS*(CO2S/CXylenS - 3/1*(CXylenS - cxyl)/CXylenS);

  (ka*kr*cO2*cxyl)/(ka*cO2 + 3*kr*cxyl)
  ]

The reaction is carried out at 600 Kelvin. The single second order differential equation becomes:

First@ NDSolve[{D[
     D[\[Xi][\[Zeta]], \[Zeta]], \[Zeta]] == -(1/\[Zeta])*
      D[\[Xi][\[Zeta]], \[Zeta]] + ( 
     Lp^2*(RateOfReaction[600, \[Xi][\[Zeta]]*CXylenS])*\[Rho]p)/(
     CXylenS*Deff), \[Xi]'[10^-9] == 0 , \[Xi][1] == 
    1}, {\[Xi]}, {\[Zeta], 10^-9, 1}]

At each instance where it states "10^-9", it has been use to replace 0, as to avoid division by zero in the fraction -(1/Zeta). Mathematica evaluates it without problems.

Now I remove the stoichiometric relation between the O-xylene and oxygen concentrations. My new rate expression becomes:

RateOfReaction2[T_, Coxylen_, CO2_] :=
 Module[{temp , kr, ka, cxyl, cO2},
  temp = T;
  cxyl = Coxylen;
  cO2 = CO2;
  kr = 8.84*10^3*Exp[-(8170/temp)]*0.001;
  ka = 199*10^3*Exp[-(11880/temp)]*0.001;

  (ka*kr*cO2*cxyl)/(ka*cO2 + 3*kr*cxyl)
  ]

In effect, I will now need to use two second order coupled differential equations. These become:

First@ NDSolve[{D[
     D[\[Xi]1[\[Zeta]], \[Zeta]], \[Zeta]] == -(1/\[Zeta])*
      D[\[Xi]1[\[Zeta]], \[Zeta]] + ( 
     Lp^2*(RateOfReaction2[
        600, \[Xi]1[\[Zeta]]*CXylenS, \[Xi]2[\[Zeta]]*
         CO2S])*\[Rho]p)/(CXylenS*Deff), 
   D[D[\[Xi]2[\[Zeta]], \[Zeta]], \[Zeta]] == -(1/\[Zeta])*
      D[\[Xi]2[\[Zeta]], \[Zeta]] + ( 
     Lp^2*(RateOfReaction2[
        600, \[Xi]1[\[Zeta]]*CXylenS, \[Xi]2[\[Zeta]]*
         CO2S])*\[Rho]p)/(CO2S*Deff), \[Xi]1'[10^-6] == 
    0, \[Xi]1[1] == 1, \[Xi]2'[10^-6] == 10^-6, \[Xi]2[1] == 
    1}, {\[Xi]1, \[Xi]2}, {\[Zeta], 10^-6, 1}]

I've used the samme approach above as previous, substituting some of the zeros by "10^-6", to circumvent division by zero. However, I still get the following errors:

Power::infy: Infinite expression 1/0. encountered. >>
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>
Power::infy: Infinite expression 1/0. encountered. >>
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>
Power::infy: Infinite expression 1/0.^2 encountered. >>
General::stop: Further output of Power::infy will be suppressed during this calculation. >>
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>
General::stop: Further output of Infinity::indet will be suppressed during this calculation. >>
NDSolve::ndnum: Encountered non-numerical value for a derivative at \[Zeta] == 1.`*^-6. >>

It may be that I have been looking at it so much that I can't find my own mistakes, but I hope someone else could hint me in the right direction.

I should quickly mention that zeta and xi are dimensionless variables, containing the dimensionless distance inside the pellet, and reactant concentration respectively.

Some defined variables I use are:

Lp = (7.5*10^-5)/2; (* Pellet Halfwidth*)
\[Rho]p = 2.3*10^6;  (* Pellet density*)
Deff = 1.4*10^-7; (*Effective diffusivity of O-xylene*)
CXylenS = 0.01*1; (*Surface concentration of O-xylene*)
CO2S = 0.99*1; (*Surface concentration of Oxygen*)
Ts = 600; (*Surface temperature*)

Best Regards