Solving a chemical system using NDSolve

Question

I have to solve a set of coupled differential equations and am not sure how to proceed. It is a model of a chemical process where I'm interested to find the species concentration as $r$ varies.
For a chemical reaction given by $A + B \rightarrow AB + hv$,

$\displaystyle \frac{d(AB)}{dr}=k*n(A)*n(B)$, is the rate of production of the species $AB$

So far I have built the following code

Te[ϕ_]:= 2200 + 440*Cos[ϕ] 
Ra[ϕ_]:= 2 + 0.25*Cos[ϕ + Pi]

t[r_,ϕ_]:= Te[ϕ]*(1 - Sqrt[1 - (Ra[ϕ]/r)^2])^(1/4)

Note that, for my case, $k_{i}$ also changes with ϕ (to make it simple to solve, I can set ϕ manually from 0 to 2π)

k1[r_,ϕ_]:= (6.99*10^-14)*((t[r,ϕ]/300)^2.8)*Exp[-1950/t[r,ϕ]] 
k2[r_,ϕ_]:= (1.59*10^-11)*((t[r,ϕ]/300)^1.2)*Exp[-9610/t[r,ϕ]]
k17[r_,ϕ_]:= (3.14*10^-13)*((t[r,ϕ]/300)^2.7)*Exp[-3150/t[r,ϕ]] 
k18[r_,ϕ_]:= (2.05*10^-12)*((t[r,ϕ]/300)^1.52)*Exp[-1736/t[r,ϕ]]
k62[r_,ϕ_]:= (1.77*10^-11)*Exp[178/t[r,ϕ]]
k63[r_,ϕ_]:= (1.85*10^-11)*((t[r,ϕ]/300)^0.95)*Exp[-8571/t[r,ϕ]] 
k94[r_,ϕ_] := (1.65*10^-12)*((t[r,ϕ]/300)^1.14)*Exp[-50/t[r,ϕ]] 

Now I should build the rate equations, based on the differential equantion I pointed in the begining. And there is where I'm not sure how to proceed. The Wolfram Documentation has a good example about how to work with chemical systems, but there the $k_i$ are constants which make it much more simple to solve. There, the rate equations would be write:

r1 = k1 nH[r] nOH[r];
r2 = k2 nH[r] nH2O[r];
r17 = k17 nH2[r] nO[r];
r18 = k18 nH2[r] nOH[r];
r62 = k62 nO[r] nOH[r];
r63 = k63 nO[r] nH2O[r];
r94 = k94 nOH[r] nOH[r];

As the $k's$ in my case are dependent on $r$ and $ϕ$, should I write it like

r1[r_,ϕ_] = k1[r,ϕ]*nH[r]*nOH[r];

for each rate equation?

Just as an example, I wrote the code to solve the system where $k_{i}$ are constants:

r1 = k1 nH[r] nOH[r];
r2 = k2 nH[r] nH2O[r];
r17 = k17 nH2[r] nO[r];
r18 = k18 nH2[r] nOH[r];
r62 = k62 nO[r] nOH[r];
r63 = k63 nO[r] nH2O[r];
r94 = k94 nOH[r] nOH[r];
eqns = {nH'[r] == (r17 + r18 + r62) - (r1 + r2), 
   nOH'[r] == (r2 + r17 + r63) - (r1 + r18 + r62 + r94), 
   nO'[r] == (r1 + r94) - (r62 + r63), 
   nH2'[r] == (r1 + r2) - (r17 + r18), nH2O'[r] == (r18 + r94) - (r2)};


eqEqn = {nH[r] + nOH[r] + nO[r] + nH2[r] + nH2O[r] + nCONS[r] == 
    0.50001459};
ic = {nH[2] == 10^-7, nOH[2] == 4.9*10^-7, nO[2] == 10^-5, 
   nH2[2] == 5*10^-1, nH2O[2] == 4*10^-6, nCONS[2] == 0};

params = {k1 -> 7.1143*10^-13, k2 -> 2.41635*10^-12, 
   k17 -> 1.744057*10^-11, k18 -> 2.00061*10^-11, 
   k62 -> 1.91656*10^-11, k63 -> 2.70681*10^-12, 
   k94 -> 1.594279*10^-11 };

sol = NDSolve[{eqs, igEqs, ic} /. params, {nH, nOH, nO, nH2, 
    nH2O}, {r, 2, 8}];

and it worked pretty fine. How could I rewrite this code to work with $k$ changing with $r$? Moreover, is it possible to relate the independent variable $r$ that I used both in the equations and NDSolve? For example, in the NDSolve I've used $r=2$, as the first step. The program will use it to evaluate the previous equations

t[2,ϕ_]:= Te[ϕ]*(1 - Sqrt[1 - (Ra[ϕ]/2)^2])^(1/4);
k1[2,ϕ_]:= (6.99*10^-14)*((t[2,ϕ]/300)^2.8)*Exp[-1950/t[2,ϕ]] 
k2[2,ϕ_]:= (1.59*10^-11)*((t[2,ϕ]/300)^1.2)*Exp[-9610/t[2,ϕ]]
k17[2,ϕ_]:= (3.14*10^-13)*((t[2,ϕ]/300)^2.7)*Exp[-3150/t[2,ϕ]] 
k18[2,ϕ_]:= (2.05*10^-12)*((t[2,ϕ]/300)^1.52)*Exp[-1736/t[2,ϕ]]
k62[2,ϕ_]:= (1.77*10^-11)*Exp[178/t[2,ϕ]]
k63[2,ϕ_]:= (1.85*10^-11)*((t[2,ϕ]/300)^0.95)*Exp[-8571/t[2,ϕ]] 
k94[2,ϕ_] := (1.65*10^-12)*((t[2,ϕ]/300)^1.14)*Exp[-50/t[2,ϕ]]  

and use the results on NDSolve? If don't, is there some other way to appoach this problem?
I hope I made myself clear. Thank you in advance!

Answer 1

I think you just have to rewrite your "r" equations.

I took your code:

Te[ϕ_] := 2200 + 440*Cos[ϕ]
Ra[ϕ_] := 2 + 0.25*Cos[ϕ + π]

t[r_, ϕ_] := Te[ϕ]*(1 - Sqrt[1 - (Ra[ϕ]/r)^2])^(1/4)

k1[r_, ϕ_] := (6.99*10^-14)*((t[r, ϕ]/300)^2.8)*
  Exp[-1950/t[r, ϕ]]
k2[r_, ϕ_] := (1.59*10^-11)*((t[r, ϕ]/300)^1.2)*
  Exp[-9610/t[r, ϕ]]
k17[r_, ϕ_] := (3.14*10^-13)*((t[r, ϕ]/300)^2.7)*
  Exp[-3150/t[r, ϕ]]
k18[r_, ϕ_] := (2.05*10^-12)*((t[r, ϕ]/300)^1.52)*
  Exp[-1736/t[r, ϕ]]
k62[r_, ϕ_] := (1.77*10^-11)*Exp[178/t[r, ϕ]]
k63[r_, ϕ_] := (1.85*10^-11)*((t[r, ϕ]/300)^0.95)*
  Exp[-8571/t[r, ϕ]]
k94[r_, ϕ_] := (1.65*10^-12)*((t[r, ϕ]/300)^1.14)*
  Exp[-50/t[r, ϕ]]

Here I only changed your code the way your ks depend on r and ϕ, like you defined them above:

r1 = k1[r, ϕ] nH[r] nOH[r];
r2 = k2[r, ϕ] nH[r] nH2O[r];
r17 = k17[r, ϕ] nH2[r] nO[r];
r18 = k18[r, ϕ] nH2[r] nOH[r];
r62 = k62[r, ϕ] nO[r] nOH[r];
r63 = k63[r, ϕ] nO[r] nH2O[r];
r94 = k94[r, ϕ] nOH[r] nOH[r];

With your equations and boundary conditions:

eqns = {nH'[r] == (r17 + r18 + r62) - (r1 + r2), 
   nOH'[r] == (r2 + r17 + r63) - (r1 + r18 + r62 + r94), 
   nO'[r] == (r1 + r94) - (r62 + r63), 
   nH2'[r] == (r1 + r2) - (r17 + r18), nH2O'[r] == (r18 + r94) - (r2)};

eqEqn = {nH[r] + nOH[r] + nO[r] + nH2[r] + nH2O[r] + nCONS[r] == 
   0.50001459}; ic = {nH[2] == 10^-7, nOH[2] == 4.9*10^-7, 
  nO[2] == 10^-5, nH2[2] == 5*10^-1, nH2O[2] == 4*10^-6, 
  nCONS[2] == 0};

This gives at least a numerical solution. Note, that I defined ϕ to be 2π (you mentioned you can change it manually):

sol = NDSolve[{eqns, eqEqn, ic} /. ϕ -> 2π, {nH, 
   nOH, nO, nH2, nH2O}, {r, 2, 8}]

In addition: If you not want to change ϕ manually, just generate a table over your different ϕ (note: I didnt get a solution for all of them):

solTable = Table[
   NDSolve[{eqns, eqEqn, ic} //. ϕ -> ϕi, {nH, nOH, nO, nH2,
      nH2O}, {r, 2, 8}],
   {ϕi, 0, 2π, π/10}];

(* delete empty entrys *)
shortSolTable = DeleteCases[solTable, {}];

To visualise your solutions you can use:

Table[nH[r] /. shortSolTable[[i, 1]], {i, 1, Length[shortSolTable]}];
Plot[%, {r, 2, 8}, PlotLegends -> Automatic,PlotLabel -> "nH[r,ϕ]"]