# Solving heat equation on a cylinder with insulated ends and convective boundary conditions

I am trying to solve heat equation on a cylinder whose ends are thermally insulated and its circular face is exposed to convection. Therefore I have Neumann boundary condition on all faces of the cylinder. Here is my code:

NDSolveValue[{Laplacian[u[x, y, z], {x, y, z}] == NeumannValue[0, z == 0] +
NeumannValue[0, z == 1] + NeumannValue[1 - u[x, y, z],
x^2 + y^2 == 1]}, u, {x, y, z} \[Element] Cylinder[{{0, 0, 0}, {0, 0, 1}}, 1]]


However I get the error:

NDSolveValue::bcnop: No places were found on the boundary where Coordinate was True, so BoundaryCondition[{Robin,{1,1},{CompiledFunction[{10,11.,5568},{_Real,_Real,_Real},{{3,0,0},{3,0,1},{3,0,2},{3,2,0}},{{{{1.}},{3,2,0}}},{0,0,3,0,1},{{1}},Function[{x,y,z},{{1.}},Listable],Evaluate],CompiledFunction[{10,11.,5568},{_Real,_Real,_Real},{{3,0,0},{3,0,1},{3,0,2},{3,2,0}},{{{{-1.}},{3,2,0}}},{0,0,3,0,1},{{1}},Function[{x,y,z},{{-1.}},Listable],Evaluate]},1},Coordinate,CompiledFunction[{10,11.,5568},{_Real,_Real,_Real},<<5>>,Evaluate],NeumannValue[1-u,x^2+y^2==1]] will effectively be ignored.


My question: what does No places were found on the boundary where Coordinate was True mean? As far I can see I have specified the faces of the cylinder correctly in the NeumannValue boundary conditions. Any help is much appreciated.

• Use NeumannValue[1 - u[x, y, z], x^2 + y^2 >= 0.9999]. Oct 10, 2020 at 12:40
• @flinty NeumannValue[1 - u[x, y, z], x^2 + y^2 >= 0.99] works for me, but any number greater than 0.99 doesn't work and I get the same error as before. What's going on here?
– Deep
Oct 10, 2020 at 12:47
• I suspect it could be due to the discretization of the cylinder. Oct 10, 2020 at 12:49
• @flinty Thanks a lot. How can I refine the discretization?
– Deep
Oct 10, 2020 at 12:53

To avoid having to tweak the boundary condition, load the finite elements package and make an actual mesh:

<< NDSolveFEM

mesh = ToElementMesh @ Cylinder[{{0, 0, 0}, {0, 0, 1}}, 1];

NDSolveValue[
{
Laplacian[u[x, y, z], {x, y, z}] ==
+ NeumannValue[0, z == 0]
+ NeumannValue[0, z == 1]
+ NeumannValue[1 - u[x, y, z], x^2 + y^2 == 1]
},
u,
Element[{x, y, z}, mesh]
]


By default the mesh will be second-order, and perhaps this is why it is able to handle the curved boundary properly. It seems that ToElementMesh is able to handle curved boundaries much better than the default discretisation method used by NDSolveValue.

• Thanks a lot. Your solution works. Do you know how I can make the mesh finer?
– Deep
Oct 11, 2020 at 4:38
• @Deep Take a look at the Element Mesh Generation tutorial and the options for ToElementMesh. Most useful to you will be the options MaxCellMeasure and "MaxBoundaryCellMeasure". Oct 11, 2020 at 12:23

The OP indicated that they are on version 11.0 that does not include the OpenCascadeLink. I do not have version 11, so I do not know if this works, but it does not depend on OpenCascade. Note that the curved surface came out as ElementMarker==3 in this case.

Needs["NDSolveFEM"]
c1 = Cylinder[{{0, 0, 0}, {0, 0, 1}}, 1];
bmesh = ToBoundaryMesh[c1];
groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp
bmesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]
mesh = ToElementMesh[bmesh];
mesh["Wireframe"]
ufun = NDSolveValue[{Laplacian[u[x, y, z], {x, y, z}] ==
NeumannValue[1 - u[x, y, z], ElementMarker == 3]},
u, {x, y, z} ∈ mesh];
SliceContourPlot3D[
ufun[x, y, z], "CenterPlanes", {x, y, z} ∈
Cylinder[{{0, 0, 0}, {0, 0, 1}}, 1], PlotLegends -> Automatic]


You can use OpenCascadeLink define the geometry and it will create ElementMarkers to the faces that you may refer to in your boundary condition specification. This will avoid guessing at what discretization is required when the object or scale changes.

Here is an example. Note that the $$\color{Red}{Red\ Surface}$$ corresponds to the curved face and is ElementMarker==1.

Needs["OpenCascadeLink"]
Needs["NDSolveFEM"]
cyl = OpenCascadeShape[c1 = Cylinder[{{0, 0, 0}, {0, 0, 1}}, 1]];
groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp
bmesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]
mesh = ToElementMesh[bmesh];
mesh["Wireframe"]
ufun = NDSolveValue[{Laplacian[u[x, y, z], {x, y, z}] ==
NeumannValue[1 - u[x, y, z], ElementMarker == 1]},
u, {x, y, z} ∈ mesh];
SliceContourPlot3D[
ufun[x, y, z], "CenterPlanes", {x, y, z} ∈
Cylinder[{{0, 0, 0}, {0, 0, 1}}, 1], PlotLegends -> Automatic]


• Unfortunately I am unable to load OpenCascadeLink. I am using Mathematica 11.0.
– Deep
Oct 11, 2020 at 4:17
• @Deep Look at the update and see if it works. I do not have v11, so I cannot test. Oct 11, 2020 at 15:21
• Thank you very much.
– Deep
Oct 12, 2020 at 3:21

I think it has to do with discretizing the region. Consider:

NDSolveValue[{Laplacian[u[x, y, z], {x, y, z}] ==
NeumannValue[0, z == 0] + NeumannValue[0, z == 1] +
NeumannValue[1 - u[x, y, z], 0.999 <= x^2 + y^2 <= 1.001]}
, u, {x, y, z} \[Element] Cylinder[{{0, 0, 0}, {0, 0, 1}}, 1]]


This produces your error. However, if we soften up the condition x^2+y^2==1 a bit, then it works:

NDSolveValue[{Laplacian[u[x, y, z], {x, y, z}] ==
NeumannValue[0, z == 0] + NeumannValue[0, z == 1] +
NeumannValue[1 - u[x, y, z], 0.99 <= x^2 + y^2 <= 1.01]}
, u, {x, y, z} \[Element] Cylinder[{{0, 0, 0}, {0, 0, 1}}, 1]]
(*InterpolatingFunction[{{-1., 1.}, {-1., 1.}, {0., 1.}}, <>]*)
`