Boundary condition error with a partial differential equation

Questions
1

Hallo everyone :slight_smile: ,

I’m new to Mathematica so I need some help solving this steady state mass transfer problem in a pipe.
Due to simplisity I set all factors to one and the R as well. Here is my partial differential equation

NDSolve[{2*(1 - (r/R)^2)*
    D[\[Rho][r, z], {z}] == (r*D[\[Rho][r, z], {r, 2}] + 
      D[\[Rho][r, z], {r}])/r, 
  D[\[Rho][0, z], {r}] == 0, \[Rho][R, z] == 0, \[Rho][r < R, 0] == 
   1}, \[Rho], {r, 0, R}, {z, 0, 10}, Reals]

I get the error that my furst boundary condition is aways true and I can’t fix this error.
Do you have an Idea how to fix it. I want to plot the solution as well.
Thanks
Slavi

2

The usual cause of that warning is that something (probably ? in your case) is already defined.
You might try

??
?z
?r

to make sure none has a value. Or to just remove any definitions,

Clear[?, z, r]
3

To solve the PDE using NDSolve, you will have to formulate the problem as follows:

R = 1 + 10^-15;
eqns = (2 (1 - (r/R)^2)) D[rho[r, z], z] == (r D[rho[r, z], r, r] + D[rho[r, z], r])/r ;
bcs = {rho[R, z] == 0, Derivative[1, 0][rho][10^-19, z] == 0, rho[r, 0] == 1};

sol = NDSolveValue[{eqns, bcs}, rho[r, z], {r, 10^-19, R}, {z, 0, 10}, Method -> {"PDEDiscretization" -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 10^4}}}]

A couple of points need to be made here. There are singularities at r=0 and (r/R) = 1. So, in order to avoid the left singulariy, the boundary conditions are offset by a tiny bit. The second point is associated with the initial condition rho[z,0] = 1. This is inconsistent with the BC rho[R,z]=0. This is why the additional options have been added in NDSolve to give additional weight to the BCs.

The solution can be plotted as:

Plot3D[sol, {r, 0, R}, {z, 0, 0.5}, PlotRange -> All]

image.jpg
image.jpg345×270 30.4 KB

4

Paritosh, this is very interesting. Is it possible to plot this solution?

5

Paritosh, thank you very much. You probably use Mathematica 9 that’s why it didn’t work with the plot. But I plot it this way.

R = 1 + 10^-15;
eqns = (2 (1 - (r/R)^2)) D[rho[r, z], 
     z] == (r D[rho[r, z], r, r] + D[rho[r, z], r])/r;
bcs = {rho[R, z] == 0, Derivative[1, 0][rho][10^-19, z] == 0, 
   rho[r, 0] == 1};

sol = NDSolve[{eqns, bcs}, rho[r, z], {r, 10^-19, R}, {z, 0, 10}]  (*NDSolve instead of NDSolveValue *)
Plot3D[rho[r, z] /. sol, {r, 0, R}, {z, 0, 0.5}, PlotRange -> All] (* I added rho[r, z] /.sol here *)

I wonder how can I solve the problem, if I have to plot this over the hole pipe. The radius is not going from 0 to R, but from -R to R. It sould be like this plot but with additional mirrored plot of this across the z axis. In this case I’ll have somewere problem whit the zero.

6

If you are using V8, then you would have to change the syntax a little bit as follows:

R = 1 + 10^-15;
eqns = (2 (1 - (r/R)^2)) D[rho[r, z],z] == (r D[rho[r, z], r, r] + D[rho[r, z], r])/r;
bcs = {rho[R, z] == 0, Derivative[1, 0][rho][10^-19, z] == 0,rho[r, 0] == 1};

sol = Quiet@NDSolve[{eqns, bcs}, rho, {r, 10^-19, R}, {z, 0, 10},Method -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> {True,"ScaleFactor" -> 10^4}}] ;

If you want to plot the solution from {-R, R}, then you can just mirror the solution.

Show[Plot3D[Evaluate[rho[r, z] /. sol], {r, 0, R}, {z, 0, 0.5},PlotRange -> All] ,
Plot3D[Evaluate[rho[-r, z] /. sol], {r, -R, 0}, {z, 0, 0.5}]]

8170image.jpg
8170image.jpg357×269 27.4 KB