Laplace's equation in spherical coordinates with Neumann b.c

Question

I am trying to find the temperature field in a semi-infinite solid on whose surface there is an isotherm spherical cap sunken by a length p. For example:

R = 10;
p = 3;
SphericalPlot3D[
 1/2 (Sqrt[2] Sqrt[Cos[2 θ] (R - p)^2 + 2 p R + R^2 - p^2] + 
    2 Cos[θ] (R - p)), {θ, Pi/2, 3/2 Pi}, {ϕ, 0, 
  2 Pi}]

The rest of the surface is adiabatic.

In a spherical coordinate system (physical convention) centered at the "center" of the cap, the PDE is: $$\frac{\partial }{\partial r}\left(r^2\frac{\partial T}{\partial r}\right)+\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial T}{\partial \theta }\right)=0$$ With boundary conditions in dimensionless form given by:

$T=0$, $r\to \infty$, this sets $T$ equal to the initial value far from the cap,

$T=1$, $r=\frac{1}{2} \left(\sqrt{2} \sqrt{-p^2+(1-p)^2 \cos (2 \theta )+2 p+1}+2 (1-p) \cos (\theta )\right)$, this imposes $T$ on the cap

$\frac{\partial T}{\partial \theta}\bigg| _{\theta=\pi/2}=0$, adiabatic condition

$\frac{\partial T}{\partial \theta}\bigg| _{\theta=\pi}=0$, symmetry condition.

I tried this:

p = 0.2;
    
boundaries = {-r + 1/2 (Sqrt[2] Sqrt[Cos[2 θ] (1 - p)^2 + 2 p + 1 - p^2] + 
      2 Cos[θ] (1 - p)), r - 100, -θ + Pi/2, θ - Pi}
Ω = ImplicitRegion[And @@ (# <= 0 & /@ boundaries), {r, θ}];
RegionPlot[Ω, PlotRange -> {{0, 3}, {1, 5}}]
NDSolveValue[{r^2 D[T[r, θ], {r, 2}] + 
    2 r D[T[r, θ], r] + D[T[r, θ], {θ, 2}] + 
    Cot[θ] D[T[r, θ], θ] == {NeumannValue[0., 
     boundaries[[3]] == 0], 
    NeumannValue[0., boundaries[[4]] == 0]}, {DirichletCondition[
    T[r, θ] == 1., boundaries[[1]] == 0.],
   DirichletCondition[T[r, θ] == 0., 
    boundaries[[2]] == 
     0.]}}, T, {r, θ} ∈ Ω]

But it does not seem to work. I have two Dirichlet conditions and two Neumann conditions, but I don't know if I inserted them in NDSolve in the right way.

Answer 1

Two issues here.

First of all, you've chosen 100 to approximate Infinity, which is way too large in this case. Something like 5 is OK:

p = 0.2; inf= 5;
boundaries = {-r + 1/2 (Sqrt[2] Sqrt[Cos[2 θ] (1 - p)^2 + 2 p + 1 - p^2] + 
               2 Cos[θ] (1 - p)), r - inf, -θ + Pi/2, θ - Pi};
<< NDSolve`FEM`
Ω = ToElementMesh@ImplicitRegion[And @@ (# <= 0 & /@ boundaries), {r, θ}];

ToElementMesh is added to help NDSolve analyzing the domain properly, otherwise the femcbtd warning will pop up, at least in v12.2. (Alternatively, DiscretizeRegion can be used in place of ToElementMesh, which is a bit slower. )

The next issue is, you haven't set NeumannValue correctly. If you read the Details section of NeumannValue and the FEM document carefully, you might notice NeumannValue is actually defined based on the formal form of a PDE. How can we check the formal form of a PDE? The new-in-12.2 NDSolve`FEM`GetInactivePDE does the work. (If you're not yet in v12.2, try the function in this post. )

seq = Sequence[{r^2 D[T[r, θ], {r, 2}] + 2 r D[T[r, θ], r] + 
     D[T[r, θ], {θ, 2}] + Cot[θ] D[T[r, θ], θ] == 
    0, {DirichletCondition[T[r, θ] == 1., boundaries[[1]] == 0.], 
    DirichletCondition[T[r, θ] == 0., boundaries[[2]] == 0.]}}, 
  T, {r, θ} ∈ Ω]

{state} = NDSolve`ProcessEquations@seq

GetInactivePDE@state

exprInBlueBox = -{{-r^2, 0}, {0, -1}} . Inactive[Grad][T[r, θ], {r, θ}]; 

The normal vector $\overset{\rightharpoonup }{n}=(0,1)$ at $\theta=\pi$, so the left hand side (LHS) of Neumann b.c. at $\theta=\pi$ is:

normalVector = {0, 1};
- exprInBlueBox . normalVector // Activate
(* - Derivative[0, 1][T][r, θ] *)

The normal vector $\overset{\rightharpoonup }{n}=(0,-1)$ at $\theta=\pi/2$, so the LHS of Neumann b.c. at $\theta=\pi/2$ is:

normalVector = {0, -1};
- exprInBlueBox . normalVector // Activate
(* Derivative[0, 1][T][r, θ] *)

Thus $\left.\frac{\partial T}{\partial \theta}\right| _{\theta=\pi/2}=0$ and $\left.\frac{\partial T}{\partial \theta}\right| _{\theta=\pi}=0$ are equivalent to zero NeumannValue in this case. Once again, according to the Details section of NeumannValue:

…not specifying a boundary condition at all is equivalent to specifying a Neumann 0 condition.

In other words, we don't need to explicitly set NeumannValue for your problem. So the problem can be solved with:

sol = NDSolveValue@seq

DensityPlot[sol[Sqrt[x^2 + z^2], ArcTan[z, Abs@x]], {x, -5, 5}, {z, -5, 0}, 
 AspectRatio -> Automatic, PlotPoints -> 100, PlotRange -> All]

You may further adjust value of inf and MaxCellMeasure option of ToElementMesh to see how the solution varies.

Answer 2

In a previous answer 240190, I showed how one could use anisotropic meshing to add a DirichletCondition at "infinity" for a 1D problem. In this answer, I shall extend the technique to a 2D problem.

Geometry description

In many FEM software packages, problems with spherical symmetry can be posed as an axisymmetric problem. Since it is easier for me to think in these terms, I will recast the problem.

As I understood the system, a spherical cap is embedded in a semi-infinite domain, as I have sketched below. The y-axis is the symmetry axis.

Helper functions

Mesh helper functions

I use some of the following helper functions to construct an anisotropic mesh based on connecting and extending edge segments. A structured Quad mesh can then be easily constructed using RegionProduct

(*Import required FEM package*)
Needs["NDSolve`FEM`"];
(*Define Some Helper Functions For Structured Meshes*)
pointsToMesh[data_] := 
  MeshRegion[Transpose[{data}], 
   Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
unitMeshGrowth[n_, r_] := 
 Table[(r^(j/(-1 + n)) - 1.)/(r - 1.), {j, 0, n - 1}]
meshGrowth[x0_, xf_, n_, r_] := (xf - x0) unitMeshGrowth[n, r] + x0
firstElmHeight[x0_, xf_, n_, r_] := 
 Abs@First@Differences@meshGrowth[x0, xf, n, r]
lastElmHeight[x0_, xf_, n_, r_] := 
 Abs@Last@Differences@meshGrowth[x0, xf, n, r]
findGrowthRate[x0_, xf_, n_, fElm_] :=(*Quiet@*)
 Abs@FindRoot[
    firstElmHeight[x0, xf, n, r] - fElm, {r, 0.00000001, 
     100000000/fElm}, Method -> "Brent"][[1, 2]]
meshGrowthByElm[x0_, xf_, n_, fElm_] := 
 N@Sort@Chop@meshGrowth[x0, xf, n, findGrowthRate[x0, xf, n, fElm]]
meshGrowthByElm0[len_, n_, fElm_] := meshGrowthByElm[0, len, n, fElm]
flipSegment[l_] := (#1 - #2) & @@ {First[#], #} &@Reverse[l];
leftSegmentGrowth[len_, n_, fElm_] := meshGrowthByElm0[len, n, fElm]
rightSegmentGrowth[len_, n_, fElm_] := 
 Module[{seg}, seg = leftSegmentGrowth[len, n, fElm];
  flipSegment[seg]]
reflectRight[pts_] := 
 With[{rt = ReflectionTransform[{1}, {Last@pts}]}, 
  Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
reflectLeft[pts_] := 
 With[{rt = ReflectionTransform[{-1}, {First@pts}]}, 
  Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
extendMesh[mesh_, newmesh_] := Union[mesh, Max@mesh + newmesh]

Model specific helper functions

Typically, the structured Quad mesh is used on rectangular domains. I use the following helper functions to map a square UV space mesh onto the curved domain.

Clear[β, γ, rcap, rinf, rl, capMesh]
β[R_, h_] := ArcCos[(R - h)/R]
γ[R_, h_, ρ_] := ArcCos[(R - h)/ρ]
rcap[R_, h_][u_] := Module[{angle = β[R, h], r = R},
  R {Sin[angle u], -Cos[angle u]}]
rinf[R_, h_, ρ_][u_] := 
 Module[{angle = γ[R, h, ρ], r = ρ},
  r {Sin[angle u], -Cos[angle u]}]
rl[R_, h_, ρ_][u_, v_] := 
 Module[{rc = rcap[R, h][u], ri = rinf[R, h, ρ][u]},
  (ri - rc) v + rc]
capMesh[R_, h_, ρ_][rh_, rv_] := 
 Module[{sqr, crd, inc, msh, mean, mrkrs, bmrkrs, pEle, pe, pm, pcrd,
   sdf, n, leIds, bcEle, z = {0, 0}, ex = {1, 0}, ey = {0, 1}, f, g},
  sqr = RegionProduct[rh, rv];
  crd = MeshCoordinates[sqr];
  inc = ( Delete[0] /@ MeshCells[sqr, 2]);
  mean = Mean /@ GetElementCoordinates[crd, #] & /@ {inc} // First;
  mrkrs = If[#2 > rf/ρ, 2, 1] & @@@ mean;
  msh = ToElementMesh["Coordinates" -> crd, 
    "MeshElements" -> {QuadElement[inc, mrkrs]}];
  pm = epm[msh] /. {0 -> 4};
  pe = epi[msh];
  pcrd = crd[[Flatten@pe]];
  sdf = Flatten@
     Position[SignedRegionDistance[#, pcrd], _?(Abs[#] < 0.0000001 &),
       1] &;
  g = (pm[[#1]] = First@#2) &;
  MapIndexed[g, 
   sdf /@ Table[
     TransformedRegion[Line[{z, ex}], 
      RotationTransform[i 90 °, 1/2 (ex + ey)]], {i, 0, 3}]];
  pEle = {PointElement[pe, pm]};
  bmrkrs = ebm[msh];
  n = ebn[msh];
  leIds = Range@Length@n;
  f = Function[{d}, Flatten@Position[n, _?(0.9999 < d . # &), 1]];
  g = (bmrkrs[[#1]] = First@#2) &;
  MapIndexed[g, 
   f /@ Table[RotationTransform[i 90 °][-ey], {i, 0, 3}]];
  bcEle = {LineElement[ebi[msh], bmrkrs]};
  crd = rl[R, h, ρ][#1, #2] & @@@ crd;
  inc = inc /. {{i_, j_, k_, l_} :> {l, k, j, i}};
  ToElementMesh["Coordinates" -> crd, 
   "MeshElements" -> {QuadElement[inc, mrkrs]}, 
   "BoundaryElements" -> bcEle, "PointElements" -> pEle]
  ]

Mesh construction

The following workflow constructs a mesh with an angular resolution of 1°. Radially, there are two segments. There is a fine mesh in the region of interest (defined as 5X the radius) and an infinite segment that extends 10,000X the region of interest.

(*Define geometric and meshing parameters*)
R = 1; h = 1/5; rf = 5 R; Rinf = 
 10000 rf; nelmr = 80; nelminf = 40; nelmang = 90;
Print["Angular discretization segment"]
segu = Subdivide[0, 1, nelmang];
ru = pointsToMesh@segu
Print["Mesh segment in the radial region of interest"]
segr = leftSegmentGrowth[rf, nelmr, rf/100];
pointsToMesh@segr
Print["Mesh segment infinite radial domain"]
seginf = meshGrowthByElm0[Rinf - rf, nelminf, Last@segr - segr[[-2]]];
reginf = pointsToMesh@seginf
Print["Combined radial mesh segment"]
rv = pointsToMesh@(#/Last[#] &@extendMesh[segr, seginf])
mesh = capMesh[R, h, Rinf][ru, rv];
Print["Full domain"]
Show[mesh["Wireframe"], Axes -> True]
Print["Zoomed region"]
Show[mesh["Wireframe"], PlotRange -> {{0, 2}, {-R + h, -2}}, 
 Axes -> True]

PDE set up and solution

In the Heat Transfer Verification Manual there are some helper functions to create a well-formed operator for axisymmetric heat transfer problems. The code is reproduced here:

Clear[HeatTransferModelAxisymmetric, TimeHeatTransferModelAxisymmetric]
HeatTransferModelAxisymmetric[T_, {r_, z_}, k_, ρ_, Cp_, 
  Velocity_, Source_] := 
 Module[{V, Q}, 
  V = If[Velocity === "NoFlow", 
    0, ρ*Cp*Velocity . Inactive[Grad][T, {r, z}]];
  Q = If[Source === "NoSource", 0, Source];
  1/r*D[-k*r*D[T, r], r] + D[-k*D[T, z], z] + V - Q]
TimeHeatTransferModelAxisymmetric[T_, TimeVar_, {r_, z_}, k_, ρ_,
   Cp_, Velocity_, Source_] := ρ*Cp*D[T, {TimeVar, 1}] + 
  HeatTransferModelAxisymmetric[T, {r, z}, k, ρ, Cp, Velocity, 
   Source]

After all the heavy lifting has been done to create the mesh, the construction and solution of the PDE are straightforward.

parms = {k -> 1, ρ -> 1, Cp -> 1, hc -> 10, Ta -> 0};
Γhot = 
  DirichletCondition[θ[r, z] == 1, ElementMarker == 1];
Γcold = 
  DirichletCondition[θ[r, z] == 0, ElementMarker == 3];
Γconv = 0;
parmop = HeatTransferModelAxisymmetric[θ[r, z], {r, z}, 
   k, ρ, Cp, "NoFlow", "NoSource"];
op = parmop /. parms;
pde = {op == Γconv, Γhot, \
Γcold};
Tfun = NDSolveValue[pde, θ, {r, z} ∈ mesh];

Now, we can construct some plots:

Plot[{Tfun[0, -z], Tfun[z, -R + h]}, {z, 0, 5}, PlotPoints -> 100, 
 PlotRange -> {0, 1.0}, PlotLegends -> "Expressions", 
 PlotLabel -> "Temperature along symmetry edges"]
uRange = MinMax[Tfun["ValuesOnGrid"]];
legendBar = 
  BarLegend[{"TemperatureMap", uRange}, 50, 
   LegendLabel -> Style["[°C]", Opacity[0.6`]]];
options = {PlotRange -> {{-2, 2}, {-2.8, -R + h}, uRange}, 
   ColorFunction -> ColorData[{"TemperatureMap", uRange}], 
   ContourStyle -> Opacity[0.5`], ColorFunctionScaling -> False, 
   Contours -> 10, AspectRatio -> Automatic,
PlotPoints -> 100, FrameLabel -> {"r", "z"}, 
   PlotLabel -> Style["Temperature Field: θ(r,z)", 18], 
   ImageSize -> 650};
Legended[ContourPlot[Tfun[Abs[r], z], {r, -2, 2}, {z, -2.8, 0}, 
  Evaluate[options]], legendBar]

As you can see in the first plot, at a distance of 5, the temperature has decayed about 90%.

A benefit of anisotropic meshing is that one can pose some stringent questions to the model with minimal computational cost. For example, suppose you had a requirement that you needed to know the distance where temperature decayed 99.99% of the spherical cap value. One can easily find that this occurs at a distance of about 4000 as shown below:

FindRoot[Tfun[0, -z] - 0.0001, {z, 100}]
(* {z -> 4024.02} *)

Convectively cooled top surface

It is straightforward to create a convectively cooled top surface (ElementMarker==2) using a Robin-type condition. From the previously defined parms, I defined a convective heat transfer coefficient of 10 and an ambient fluid temperature of 0°. To set up, we simply need to modify the NeumannValue.

Γconv = 
  NeumannValue[hc (Ta - θ[r, z]), ElementMarker == 2] /. parms;
pde = {op == Γconv, Γhot, \
Γcold};
Tfun = NDSolveValue[pde, θ, {r, z} ∈ mesh];

We can plot the solution as before:

Plot[{Tfun[0, -z], Tfun[z, -R + h]}, {z, 0, 5}, PlotPoints -> 100, 
 PlotRange -> {0, 1.0}, PlotLegends -> "Expressions", 
 PlotLabel -> "Temperature along symmetry edges"]
uRange = MinMax[Tfun["ValuesOnGrid"]];
legendBar = 
  BarLegend[{"TemperatureMap", uRange}, 50, 
   LegendLabel -> Style["[°C]", Opacity[0.6`]]];
options = {PlotRange -> {{-2, 2}, {-2.8, -R + h}, uRange}, 
   ColorFunction -> ColorData[{"TemperatureMap", uRange}], 
   ContourStyle -> Opacity[0.5`], ColorFunctionScaling -> False, 
   Contours -> 10, AspectRatio -> Automatic,
PlotPoints -> 100, FrameLabel -> {"r", "z"}, 
   PlotLabel -> Style["Temperature Field: θ(r,z)", 18], 
   ImageSize -> 650};
Legended[ContourPlot[Tfun[Abs[r], z], {r, -2, 2}, {z, -2.8, 0}, 
  Evaluate[options]], legendBar]

Comparison with another code

When possible, it is often conducive to compare the Mathematica results with another simulation code. To simulate boundary conditions at infinity, the FEM software COMSOL introduces an Infinite Element Domain (IED) concept below.

A large scaling factor (e.g., 1000) is applied to the equations in the IED.

As shown below, there is an excellent agreement between the Mathematica and COMSOL simulations. That should give one more confidence in the validity of this approach to solve the infinite domain problem.

Answer 3

One can also consider the 3D statement of the problem. Solution of a such linear problem is not so time consuming nowadays. For mesh generation let's take advantage of OpenCascadeLink procedures which are very useful for tessellation of domains with complex geometry. Let's $r$ is a radius of a cap and $R$ is "infinity" radius. The computational region is defined by a difference between spherical wedge of radius $R$ and a ball of radius $r$. Whereas the spherical wedge can be defined as intersection of rectangular cuboid and a ball of radius $R$.

Definition of computation domain and FE mesh generation

Needs["OpenCascadeLink`"]
Needs["NDSolve`FEM`"]

R = 5; (*infinity radius*)
r = 10; (*radius of a cap*)
p = 0.2;

shape1 = OpenCascadeShape[Ball[{0, 0, 0}, R]];
shape2 = OpenCascadeShape[
   Hexahedron[{{0, -R, -R}, {R, -R, -R}, {Sqrt[2] R, 
      Sqrt[2] R, -R}, {0, R, -R}, {0, -R, R}, {R, -R, R}, {Sqrt[2] R, 
      Sqrt[2] R, R}, {0, R, R}}]];
shape3 = OpenCascadeShape[Ball[{p - r, 0, 0}, r]];

intersection = OpenCascadeShapeIntersection[shape1, shape2];
difference = OpenCascadeShapeDifference[intersection, shape3];
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[
  difference];(*boundary mesh geteration*)

mesh = ToElementMesh[bmesh];(*FE mesh generation*)

groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;

Show[mesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]], 
 Axes -> True, AxesLabel -> {"x", "y", "z"}, 
 AxesStyle -> RGBColor[0, 0, 0], BaseStyle -> 14]

Numerical solution

In the code below we solve the problem by means of FEM low level routines. Boundary elements with ElementMarkers=4;5belong to the surface of cap whereas on "infinity" surface ElementMarkers=1;3. This is taken into account when implementing Dirichlet BC. The rest surface of computational domain is adiabatic.

nr = ToNumericalRegion[mesh];
vd = NDSolve`VariableData[{"DependentVariables", 
     "Space"} -> {{u}, {x, y, z}}];
sd = NDSolve`SolutionData["Space" -> nr];

pded = InitializePDECoefficients[vd, sd, 
   "DiffusionCoefficients" -> {{-IdentityMatrix[3]}}];
bcd = InitializeBoundaryConditions[vd, 
  sd, {{DirichletCondition[u[x, y, z] == 1, 
     ElementMarker == 4 || ElementMarker == 5], 
    DirichletCondition[u[x, y, z] == 0, 
     ElementMarker == 1 || ElementMarker == 3]}}]
md = InitializePDEMethodData[vd, sd];

dpde = DiscretizePDE[pded, md, sd];
dbc = DiscretizeBoundaryConditions[bcd, md, sd];
{load, stiffness} = Take[dpde["SystemMatrices"], 2];
DeployBoundaryConditions[{load, stiffness}, dbc];

res = LinearSolve[stiffness, load];

ufun = ElementMeshInterpolation[{mesh}, res];

Postprocessing

Temperature distributions along axes $y$ and $z$ should be the same

Show[
          
 Plot[ufun[x, 0, 0], {x, p, R}, PlotRange -> All, 
  PlotLegends -> {"along x axis"}, 
  PlotStyle -> {RGBColor[1, 0, 0], Thickness[0.005]}]
          ,
           
 Plot[{ufun[0, x, 0], ufun[0, 0, x]}, {x, Sqrt[r^2 - (p - r)^2], R}, 
  PlotRange -> All, PlotLegends -> {"along y axis", "along z axis"}, 
  PlotStyle -> {{RGBColor[0, 1, 0], 
     Thickness[0.005]}, {RGBColor[0, 0, 1], Thickness[0.005]}}],
         Frame -> True , FrameLabel -> {"Distance", "Temperature"}, 
 FrameStyle -> RGBColor[0, 0, 0], BaseStyle -> 18, ImageSize -> 600, 
 GridLines -> {{p, Sqrt[r^2 - (p - r)^2]}, None}, 
 GridLinesStyle -> {Dashed, RGBColor[0, 0, 0]}
      ]