7
$\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$

1 Answer 1

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 and acknowledge you have read our privacy policy.