1
$\begingroup$

I am trying to solve a simple PDE in the cylindrical coordinates:

$$ \frac{\partial c}{\partial t} = D \bigg(\frac{1}{r}\frac{\partial c}{\partial r}+\frac{\partial c^2}{\partial r^2}\bigg)$$

The initial and boundary conditions are as follows:

$$ c[t=0,r] = {0,\hspace{2mm} r<=r_0 ;\hspace{2mm} C_b, \hspace{2mm} r> r_0}$$ $$ c[t,r=\infty] = C_b$$ $$ \frac{\partial c}{\partial r}[t,r=0] = 0$$

I wrote the following code in Mathematica:

Needs["NDSolve`FEM`"];

ClearAll[DiffCoeff, Cb, c1, r, t, pde, radialDerivative, \
secondRadialDerivative, boundaryConditions, initialCondition, sol, \
AvgConcentration, r0, tmax, inf]

DiffCoeff = 1; Cb = 1; r0 = 1; tmax = 10; inf = 100;

(*Define the PDE with NeumannValue for the flux condition*)
radialDerivative = 1/r*Derivative[1, 0][c1][t, r];
secondRadialDerivative = Derivative[2, 0][c1][t, r];
pde = Derivative[0, 1][c1][t, r] - 
    DiffCoeff*(radialDerivative + secondRadialDerivative) == 
   NeumannValue[0, r == 0];

(*Boundary Condition using DirichletCondition*)
boundaryConditions = 
 DirichletCondition[c1[t, r] == Cb, 
  r == inf];
(*Initial condition:c1(0,r)=0 for r<=r0 and Cb r>r0*)
initialCondition = c1[0, r] == If[r <= r0, 0, Cb];
(*Solve the PDE*)
sol = NDSolveValue[{pde, boundaryConditions, initialCondition}, 
  c1, {t, 0, tmax}, {r, 0, inf}]

This solves with a warning:

NDSolveValue::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.

Now, I was trying to solve the warning by smoothening the mesh with the following modification to the code:

(*Create a custom mesh*)
mesh = ToGradedMesh[{{Line[{{0}, {2*r0}}], <|
     "MaximalDistance" -> 0.1|>}, {Line[{{2*r0}, {inf}}]}}]
MeshRegion[mesh]

sol = NDSolveValue[{pde, boundaryConditions, initialCondition}, 
  c1, {t, 0, tmax}, {r} \[Element] mesh]

And this gives an error which I would like to solve:

NDSolveValue::femnotime: The differential equation cannot be solved as a time dependent equation as specified, most likely because the initial conditions given at (t==0.`) are not sufficient to define an initial value problem. As a consequence the differential equation will be solved as a time independent equation.

NDSolveValue::fememrc: The ranges {{0.,300.},ElementMesh[{{0.},{0.1},{0.2},{0.3},{0.4},{0.5},{0.6},{0.7},{0.8},{0.9},<<431>>},{LineElement[{{1,2,222},{2,3,223},{3,4,224},{4,5,225},{5,6,226},{6,7,227},{7,8,228},{8,9,229},{9,10,230},{10,11,231},<<210>>},{0,0,0,0,0,0,0,0,0,0,<<210>>}]},{PointElement[{{1},{221}}]},{PointElement[{{1},{221}}]}]} cannot be combined to a region. Please specify a combined region.

Now I am open to refining the original mesh in the NDSolve Function itself, maybe using the MaxCellMeasure function. But I want the mesh to be non-uniform. Will MaxCellMeasure allow me to do that?

$\endgroup$

1 Answer 1

1
$\begingroup$

You have the derivatives the wrong way round. Try this:

(*Define the PDE with NeumannValue for the flux condition*)
radialDerivative = 1/r*Derivative[0, 1][c1][t, r];
secondRadialDerivative = Derivative[0, 2][c1][t, r];
pde = Derivative[1, 0][c1][t, r] - 
    DiffCoeff*(radialDerivative + secondRadialDerivative) == 
   NeumannValue[0, r == 0];
$\endgroup$
3
  • $\begingroup$ Thank you so much, I will die in shame today :'( $\endgroup$ Oct 10, 2023 at 20:43
  • $\begingroup$ Can you once explain why my initial and boundary condition are not consistent at r = 100? The boundary condition at r = 100 is the Dirichlet boundary condition, (c1 = Cb). The initial condition at r = 100 is also (c1 = Cb). $\endgroup$ Oct 12, 2023 at 13:52
  • $\begingroup$ @Brownian_Motion, you are right I have miss read that and I have removed the statement about that. Note, however, that with the If statement for r=100 you have introduced a jump condition. At, say, at r=99 it's 0 and then it jumps to Cb. Most likely this will numerically not be an issue. For more a physically more correct model a smooth transition for 0 to Cb would be good. And that transition should be representable on the mesh. But again, this is probably minor. Happy FEM-ing. $\endgroup$
    – user21
    Oct 16, 2023 at 7:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.