Solving PDE with unreasonable results

Question

I have a simplified version of my chemical reaction problem. Basically it solves the concentration (ca) and grain (g2) diameter. The ca function is depended on particle diameter R, while g2 function is depend on time t. The ca and g2 function are depended on each others. In this case, both ca and g2 should be a function of R and t. For my understanding, it should be clarified as PDE problem and can be solved using MethodofLines.

The code is given as:

del = 10^(-8);

Monitor[{casol, gsol} =

  NDSolveValue[   {  
    D[R^2*D[ca[R, t], R], R]/R^2 -
     ( 39.7*(g2[R, t])^2*ca[R, t])/(1 + g2[R, t]*(1-g2[R, t])) == 0,

    D[g2[R, t], t] == -ca[R, t]/(1 + g2[R, t]*(1 - g2[R, t])),

    ca[R, 0.0] == 1.0, 
    ca[1.0, t] == 1.0, 
   (D[ca[R, t], R] /. R -> del) == 0,

    g2[R, 0.0] == 1.0  },

  {ca, g2}, {R, del, 1}, {t, 0, 1} ,

   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       MinPoints -> 200}} ,


   EvaluationMonitor :> (monitor = 
      Row[{"t=", CForm[t], " gsol=", CForm[g2[1.0, t]]}])], 

 monitor] 

It can solve the problem with a few warnings but the result is not right. The g2 gives you a constant value of 1 with (R,t). This is not right. A simple analytic analysis for g2 equation at R=1 can prove g2 is a function of t, and should not be constant.

Also it seems the boundary/ initial conditions are not forced in the results. For example, if I try to Plot[casol[1.0, t], {t, 0, 1}]. It expect give you value of 1,as this is part of the boundary condition, but the result shows some unreasonable value.

Not sure what cause the problems? Really appreciate some one can provide some hints.

Answer 1

The system of equations can be solved using the MethodOfLines, but several options need to be added.The physical model is not finalized, since in the process of calculation for t>= 0.44 we have g2<0. You can calculate up to t = 0.49 after which an instability arises.

del = 10^(-8); tmax = 0.49;

Monitor[{casol, gsol} = 
  NDSolveValue[{D[R^2*D[ca[R, t], R], R]/
       R^2 - (397/10*(g2[R, t])^2*ca[R, t])/(1 + 
         g2[R, t]*(1 - g2[R, t])) == 0, 
    D[g2[R, t], t] == -ca[R, t]/(1 + g2[R, t]*(1 - g2[R, t])), 
    ca[R, 0] == 1, ca[1, t] == 1, Derivative[1, 0][ca][del, t] == 0, 
    g2[R, 0] == 1}, {ca, g2}, {R, del, 1}, {t, 0, tmax}, 
   Method -> {"IndexReduction" -> Automatic, 
     "EquationSimplification" -> "Residual", 
     "PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"TensorProductGrid", 
         "MinPoints" -> 137, "MaxPoints" -> 137, 
         "DifferenceOrder" -> "Pseudospectral"}}}, 
   EvaluationMonitor :> (monitor = 
      Row[{"t=", CForm[t], " gsol=", CForm[g2[1.0, t]]}])], monitor]

{Plot3D[casol[x, t], {x, del, 1}, {t, 0, .49}, PlotRange -> All, 
  AxesLabel -> Automatic, Mesh -> None, ColorFunction -> "Rainbow"], 
 Plot3D[gsol[x, t], {x, del, 1}, {t, 0, .49}, PlotRange -> All, 
  AxesLabel -> Automatic, Mesh -> None, ColorFunction -> "Rainbow"]}