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.


NeumannValue
, related: mathematica.stackexchange.com/q/224812/1871 2. How can $\theta=3\pi/2$ in spherical coordinates? What convention do you follow? $\endgroup$