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.

Using physical convention $\theta \leq \pi$, you're right. The question was corrected according to your comment. Thank you. — umby, Apr 28 '21 at 11:17

xzczd · Accepted Answer · 2021-04-30T14:29:25.627

12

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.

edited Apr 30 '21 at 14:29

answered Apr 29 '21 at 06:32

xzczd

65,995
9
163
468

2

Great. Small suggestion, you could use NDSolveFEMToElementMesh in place of DiscretizeRegion as that would return a second order mesh with curved elements and should approximate the region and solution better. But this is probably a minor thing. – user21 Apr 29 '21 at 11:46
2

@user21 The mesh quality seems to be similar in this case, but ToElementMesh turns out to be faster. Edited. Thx for the suggestion. – xzczd Apr 29 '21 at 12:02
@xzczd Many thx. Knowing the formal form of a PDE is fundamental to set Neumann BC. If I understand, according to the documentation (can u check pls?), in my problem $a=f=0$, $\alpha=\gamma=0$, $\beta=\left{c,cot\theta\right}$ and $c=\left{-r^{2},-1\right}$. Then to set $\frac{\partial T}{\partial \theta}=0$ on the boundary3, I have to write: == NeumannValue[-1, boundaries[[3]] == 0.]. While, to impose the same b.c. on both boundary 3 and 4 implies: == NeumannValue[-1, boundaries[[3]] == 0.]+NeumannValue[-1, boundaries[[4]] == 0.]. Unfortunately, I use v11 so no NDSolveFEMGetInactivePDE. – umby Apr 30 '21 at 10:44
@xzczd, I beg your pardon, I meant $\frac{\partial T}{\partial \theta}=1$, formatting problems. In this case, is it correct what I wrote above? Many thanks. – umby Apr 30 '21 at 12:34
To use GetInactivePDE I need to create a .m function? In this case I can name it GetInactivePDE.m and put it into the packages folder? – umby Apr 30 '21 at 13:04
1

@umby You can create a .m of course, but it's not necessary, just copy and execute the code in a notebook is OK. In principle, $\left.\frac{\partial T}{\partial \theta}\right| _{\theta=\pi/2}=0$ should be interpreted to NeumannValue[1, boundaries[[3]] == 0.], because the normal vector is {0,-1} here. – xzczd Apr 30 '21 at 13:13
@user21 Strange, in v12.2 if I set NeumannValue[-1, boundaries[[4]] == 0.] in the code above, Plot[D[sol[r, th], th] /. th -> Pi // Evaluate, {r, 0.2, inf}, PlotRange -> All] shows the b.c. isn't the expected one. (Looks like something close to $\left.\frac{\partial T}{\partial \theta}\right| _{\theta=\pi}=0.4$). NeumannValue[1, boundaries[[3]] == 0.] works as expected: Plot[D[sol[r, th], th] /. th -> Pi/2 // Evaluate, {r, 0.2, 5}, PlotRange -> {0, 2}]. Have I done anything wrong? – xzczd Apr 30 '21 at 13:17
@user21 Ah, I see, for a well-behaved $T$, the only reasonable choice of Neumann Condition at $\theta=\pi$ is Neumann zero condition, because $\theta=\pi$ is the removable singularity in spherical coordinates, thus the solution should be smooth there. If a non-smooth solution is needed, one can set the b.c. at e.g. $\theta=0.99 \pi$ and make the mesh dense enough e.g. MaxCellMeasure -> 0.0001. – xzczd Apr 30 '21 at 13:42
@xzczd, great you figured this out before I even started ;-) - That's how we like it ;-) – user21 Apr 30 '21 at 13:48
@xzczd, Hopefully, I've got it. If I want $\frac{\partial T}{\partial \theta}=1$ at both $\frac{\pi}{2}$ and $\pi$, then: == NeumannValue[1, boundaries[[3]] == 0.]+NeumannValue[-1, boundaries[[4]] == 0.]. Ok? – umby Apr 30 '21 at 14:11
1

@umby Yes, but as mentioned in the comment above, NeumannValue[-1, boundaries[[4]] == 0.] is ill-posed. – xzczd Apr 30 '21 at 14:30
It is just an example to see if I understood. Many thanks to xzczd and to user 21. – umby Apr 30 '21 at 14:37
I am ashamed to ask, but how can I format a Mathematica code in comments; it is reported in the online help: indent four spaces or Ctrl + K, but this does not seem to work. – umby May 07 '21 at 13:46
@umby You just need the backtick ` , for more info read this: https://mathematica.stackexchange.com/editing-help#comment-formatting – xzczd May 07 '21 at 14:04
Thank you @xzczd, ALt + 96 on my keyboard! – umby May 07 '21 at 14:28
@umby Oh, you don't have this key on your keyboard? Usually it's below the ESC button. – xzczd May 07 '21 at 14:45
No, unfortunately I don't, but ALT + 96 is ok, otherwise I have to change the keyboard setting with Keyboard Layout Creator, to associate to AltGr + ' = `. – umby May 07 '21 at 14:49
@xzczd, I do not know if I can continue the discussion about my problem in the comments or start another question. I mean to present the complete version of my problem. I am afraid of making mistakes in using comments in contravention of community rules. Actually I don't think I could ask for the things I asked above a while ago. Or maybe I have to delete these kinds of comments not related to the original question? – umby May 07 '21 at 14:57
@umby If the complete version of the problem differs enough from the current one, it's better to start a new question. (Of course you should make the question not too localized, otherwise it may be considered as "question requires either advice from Wolfram support or the services of a professional consultant" and closed. ) While it might be considered as "too chatty", it's OK to ask about formatting issue in the comment, nobody is punished for this AFAIK. – xzczd May 07 '21 at 15:10

Tim Laska · Answer 2 · 2021-05-05T02:51:02.100

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.

Thanks for the comparative analysis @Tim Laska! Does the IED concept imply using a special elements or it is just increasing of element length far from heat source? — Oleksii Semenov, May 05 '21 at 13:06
@OleksiiSemenov You are welcome! In COMSOL, the IED is either mapped in 2D (quads) or swept in 3D (wedges or hexas) in the infinite direction. I do not think you can apply an IED to an unstructured mesh. I am not quite sure what COMSOL is doing under the hood. It is possible that they are scaling the elements invisible to the user or possibly they are scaling the thermal conductivity in a piecewise fashion to achieve the same effect. — Tim Laska, May 05 '21 at 13:51
@OleksiiSemenov From the [COMSOL documentation] (https://doc.comsol.com/5.5/doc/com.comsol.help.comsol/comsol_ref_definitions.12.115.html#1875595) applies coordinates stretching of the IED not dissimilar to what I have done. In my approach, the first element of the IED matches the width of the outer layer of the finite domain and stretches the elements to provide a gradual transition out to "infinity". — Tim Laska, May 05 '21 at 19:45
Thanks for the link @TimLaska. I also found some papers devoted to this aspect: Zienkiewicz, O. C., C. Emson, and P. Bettess. "A novel boundary infinite element." International Journal for Numerical Methods in Engineering 19.3 (1983): 393-404. and Marques, J. M. M. C., and D. R. J. Owen. "Infinite elements in quasi-static materially nonlinear problems." Computers & structures 18.4 (1984): 739-751. — Oleksii Semenov, May 06 '21 at 11:42

score 8 · Answer 3 · answered May 04 '21 at 09:19

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]}
      ]