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!
r1 = k1[r, \[Phi]] nH[r] nOH[r]
etc? $\endgroup$