6
$\begingroup$

I am trying to solve a complex transient reaction-diffusion problem numerically in 2D. Let's take a simple system to start, say

$$ A+B\rightarrow C \\ C \rightarrow A+B $$

with reaction rate $k_1$ and $k_2$, respectively.

Looking at a radially symmetric system $(r,z)$ with no flux boundary conditions, $[A]$ being generated within the domain with arbitrary spatial distribution, and fixed initial concentrations, I put together the following in Mathematica 10.1:

(* Geometric Params *)
r1 = 0; r2 = 50; z1 = 0; z2 = 0.2; rb = 1;

(* Diffusion Coefficients *)
dA = 1; dB = 1; dC = 1;

(* Rate Constants *)
k1 = 1; k2 = 2;

(* Generation Rate *)
generation = 10;

(* Helper Function for Operator *)
operator[f_, d_] := D[f, t] - d Laplacian[f, {r, z}];



sol = NDSolveValue[
  {
   operator[cA[t, r, z], dA] == -k1 cA[t, r, z] cB[t, r, z] + 
     k2 cC[t, r, z] + generation (1 - 1/(1 + Exp[-20 r + 20]))
   ,operator[cB[t, r, z], dB] == -k1 cA[t, r, z] cB[t, r, z] + 
     k2 cC[t, r, z]
   ,operator[cC[t, r, z], dC] == 
    k1 cA[t, r, z] cB[t, r, z] - k2 cC[t, r, z] 


   ,DirichletCondition[cA[t, r, z] == 0, True]
   ,DirichletCondition[cB[t, r, z] == 0, True]
   ,DirichletCondition[cC[t, r, z] == 0, True]


   ,cA[0, r, z] == 0
   ,cB[0, r, z] == 100
   ,cC[0, r, z] == 0

   }, {cA, cB, cC}, {t, 0, 100}, {r, r1, r2}, {z, z1, z2}
  , Method -> {"MethodOfLines", "TemporalVariable" -> t}
  ]

This gives me an error stating that "Nonlinear coefficients are not supported in this version of NDSolve".

Is there a work around to solve a system like this?

I tried playing around with some of the low level FEM solver functions, but get similar errors along the way.

I really appreciate the help!

$\endgroup$
5
$\begingroup$

Replace

,DirichletCondition[cA[t, r, z] == 0, True]
,DirichletCondition[cB[t, r, z] == 0, True]
,DirichletCondition[cC[t, r, z] == 0, True]

by

, DirichletCondition[cA[t, r, z] == 0, z == z1 || z == z2]
, DirichletCondition[cB[t, r, z] == 0, z == z1 || z == z2]
, DirichletCondition[cC[t, r, z] == 0, z == z1 || z == z2]
, DirichletCondition[cA[t, r, z] == 0, r == r1 || r == r2]
, DirichletCondition[cB[t, r, z] == 0, r == r1 || r == r2]
, DirichletCondition[cC[t, r, z] == 0, r == r1 || r == r2]

in order to specify boundary conditions in r and z, but not in t, which already are specified by cA[0, r, z] == 0, etc.

Even these boundary conditions generate a warning message,

NDSolveValue::ibcinc: Warning: boundary and initial conditions are inconsistent. >>

because cB[0, r, z] == 100 intersects other boundaries where cB ==0, but I do not believe this to be a serious issue. A sample plot for cB is

Plot3D[sol[[2]][t, r, .1], {t, 0, 100}, {r, r1, r2}, PlotRange -> All]

enter image description here

sol[[1]] and sol[[3]] are zero.

$\endgroup$

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.