# Gravitational potential created by a thin disc using FEM and NDSolve

## Context

I would like to compute the torque that a (thin) disc applies onto a ring. I.e. I would like to try to understand what is the impact of this outer ring on the inner disc in the simulation below. For this I would like to compute the gravitational potential generated by a (razor) thin disc and a ring. So the abstraction is the following (seen from 9 different angles) Once I know how to compute the potential it should become straightforward to compute the torque one feature applies on the other. Here I want to use FEM here for flexibility, which I will use when I will account for a more realistic abstraction of the problem. (e.g. exponential surface density profile in disc).

## Attempt

I have defined a domain

dom = ImplicitRegion[0 <= x <= 1 && -1 <= y <= 1, {x, y}];


and the Laplace operator

op = -Laplacian[u[x, y], {x, theta, y}, "Cylindrical"];


I impose the edge condition that the potential should be 1 on the disc

edge = DirichletCondition[u[x, y] == 1, 0 <= x <= 1/2 && y == 0];


When I solve

uD = NDSolveValue[{op == 0, edge},  u, {x, y} \[Element] dom]


I get

StreamPlot[-{D[#, x], D[#, y]} &@uD[x, y] // Evaluate, {x, y} \[Element] dom,AspectRatio->2] Problems:

(i) The outer box imposes a (box like) symmetry which is not in the sought solution

(ii) Strangely enough The code fails if I use the ring-like boundary condition instead:

{DirichletCondition[u[x, y] == 1, 1/2 <= x <= 3/4 && y == 0]};


## Question

How to compute the gravitational potential created by a disc (and a ring) using FEM in NDSolve?

In a broader sense I think I am asking how can FEM methods deal with PDEs with boundaries at infinity? I am guessing that one strategy might be to move the boundary sufficiently far away and increasing sampling within the inner region?

Note that my attempt above is imposing fixed potential on the disk not fixed density. I am not sure this is important or not, but ideally (to compare to the analytical solution below) fixing density would be better.

PostScriptum

I have found this (nice!) blog which provides me with an analytic solution as follows

PhiDiskData =
WolframAlpha[
"electric potential of a charged disk", {{"Result", 1},
"Input"}] // ReleaseHold;
PhiDisk = PhiDiskData /.  QuantityVariable[a_, _] -> a /. {
Q -> Pi R^2, "ElectricConstant" -> 1};
Phi = PhiDisk /. { x -> r Cos[Theta], y -> r Sin[Theta]} //
Simplify[#, Assumptions -> {r > 0}] & Clear[fD]; fD =
FullSimplify[-D[PhiDisk, {{x, y, z}}]  /. x^2 + y^2 -> r^2,
Element[z, Reals] && r > 0]  /. {x -> r Cos[Theta], y -> r Sin[Theta]};
fD = -{Sqrt[fD[]^2 + fD[]^2] // FullSimplify, fD[]}; So that

phiN = (Phi /. { Theta -> 0, R -> 1/2}); pl1 =
ContourPlot[Evaluate[phiN], {r, 0, 2}, {z, -2, 2}, Exclusions -> {},
Contours -> 15,ColorFunction -> (ColorData["RedBlueTones"][1 - #] &),
Epilog -> {Thickness[0.02], Line[{{0, 0}, { 1/2, 0}}]},
FrameLabel -> {r, z}, PlotRange -> All, AspectRatio -> 2];
pl2 = StreamPlot[(fD /. R -> 1/2), {r, 0, 2}, {z, -2, 2},
AspectRatio -> 2, StreamStyle -> White];
pl3 = Show[pl1, pl2, PlotRange -> {{0, 1.5}, {-0.5, 1}},
AspectRatio -> 1]


yields So my question amounts to finding this solution numerically.

Note that the analytic solution works nicely for rings as well (if defined as the difference between two discs.)

phiN = (Phi /. { Theta -> 0,
R -> 1}) - (Phi /. { Theta -> 0, R -> 1/2}); pl1 =
ContourPlot[Evaluate[phiN], {r, 0, 2}, {z, -2, 2}, Exclusions -> {},
Contours -> 15,ColorFunction -> (ColorData["RedBlueTones"][1 - #] &),
Epilog -> {Thickness[0.02], Line[{{1/2, 0}, { 1, 0}}]},
FrameLabel -> {r, z}, PlotRange -> All, AspectRatio -> 2];
pl2 = StreamPlot[(fD /. R -> 1) - (fD /. R -> 1/2), {r, 0, 2}, {z, -2,
2},AspectRatio -> 2, StreamStyle -> White];
pl4 = Show[pl1, pl2, PlotRange -> {{0, 1.5}, {-0.5, 1}},
AspectRatio -> 1] • I don't know if it's me, but there are a lot of things I don't understand. Among these things : – andre314 Jan 20 '18 at 16:18
• – andre314 Jan 20 '18 at 16:26

I am not too familiar with gravitational potentials but here are some suggestions:

1. I would gauge them to be 0 at infinity and to satisfy Neumann conditions on solid matter.

2. Moreover when specifying boundary conditions with inequalities, you should avoid < and > because the boundary conditions will quite likely not be set correctly, if the boundary is given by the these inequalties with < and > replaced by =.

Afterwards, I gain this, which seems quite plausible to me:

Needs["NDSolveFEM"];
dz = 1/16; reg0 =
RegionDifference[Cylinder[{{0, 0, -dz}, {0, 0, dz}}, 1],
Cylinder[{{0, 0, -dz}, {0, 0, dz}}, 1/2]];
reg = RegionDifference[Ball[{0, 0, 0}, 2], reg0];
mesh = ToElementMesh[reg, MaxCellMeasure -> 0.01,
"MaxBoundaryCellMeasure" -> 0.025];
sol = NDSolveValue[{
Laplacian[u[x, y, z], {x, y, z}] ==
NeumannValue[1, 1/2^2 <= x^2 + y^2 <= 1 && -dz <= z <= dz],
DirichletCondition[u[x, y, z] == 0, x^2 + y^2 + z^2 >= 2^2]
}, u, {x, y, z} \[Element] mesh
];
Off[InterpolatingFunction::dmval]
GraphicsRow[
{
ContourPlot[sol[x, 0, z], {x, -2, 2}, {z, -2, 2}, AspectRatio -> 1,
Contours -> 20],
ContourPlot[sol[x, y, 0], {x, -2, 2}, {y, -2, 2}, AspectRatio -> 1,
Contours -> 25]
}
] The stream lines of the gravitational force:

plr = DiscretizeRegion[
ImplicitRegion[1/2^2 <= x^2 <= 1 && -dz < z < dz, {x, z}],
PlotTheme -> "Minimal"];
g = Show[
StreamPlot[-{D[sol[x, y, z], x], D[sol[x, y, z], z]} /. y -> 0 //
Evaluate, {x, -1, 1}, {z, -1, 1}, StreamPoints -> Fine
],
plr] Finally, you could try to apply transformations of following type to your PDE (and the thin disc region) in order to map the sphere at infinity to a sphere of finite radius R:

Φ = X \[Function] X /(R - Sqrt[X.X]);
Ψ = Y \[Function] R Y/(1 + Sqrt[Y.Y]);


Note that you have to transform the Laplacian and also the Neumann conditions; that might be not too convenient...

It might be even more realistic to model the ring as a mass density ρ (for example constant on the ring) instead of cutting the rings out of the computational domain. So you would have to put it on the right hand side of the Laplacian instead of the NeumannValue. This could look like this:

reg = Ball[{0, 0, 0}, 2];
mesh = ToElementMesh[reg, MaxCellMeasure -> 0.001];
ρ = Function[{x, y, z}, Boole[1/2^2 <= x^2 + y^2 <= 1 && -dz <= z <= dz]];
sol = NDSolveValue[{
Laplacian[u[x, y, z], {x, y, z}] == ρ[x, y, z],
DirichletCondition[u[x, y, z] == 0, x^2 + y^2 + z^2 >= 2^2]
}, u, {x, y, z} ∈ mesh
];
Off[InterpolatingFunction::dmval]
g = GraphicsRow[
{
ContourPlot[sol[x, 0, z], {x, -2, 2}, {z, -2, 2}, AspectRatio -> 1,
Contours -> 20],
ContourPlot[sol[x, y, 0], {x, -2, 2}, {y, -2, 2}, AspectRatio -> 1,
Contours -> 25]
},
ImageSize -> Large
] The stream lines of the gravitational force are not perpendicular to the boundary of the disk anymore, but that cannot be expected with a mass density:

plr = DiscretizeRegion[
ImplicitRegion[1/2^2 <= x^2 <= 1 && -dz < z < dz, {x, z}],
PlotTheme -> "Minimal"];
g = Show[
StreamPlot[-{D[sol[x, y, z], x], D[sol[x, y, z], z]} /. y -> 0 //
Evaluate, {x, -2, 2}, {z, -2, 2}, StreamPoints -> Fine
],
plr] This is in quite a good accordance with the previous result. Note that the used mesh is coarser than the one before, in particular in the region around the ring. One has to fine tune it, e.g. with a suitable choice of a MeshRefinementFunction.

• Dear Henrik, could you try f[x_, y_, z_] = -Laplacian[sol[x, y, z], {x, y, z}]; g = GraphicsRow[{ContourPlot[f[x, 0, z], {x, -2, 2}, {z, -2, 2}, AspectRatio -> 1, Contours -> 2, PlotRange -> All], ContourPlot[f[x, y, 0], {x, -2, 2}, {y, -2, 2}, AspectRatio -> 1, Contours -> 2, PlotRange -> All]}, ImageSize -> Large] and possibly comment on the inaccuracy of the result? i – chris Feb 27 '18 at 17:06
• Dear chris, that's a typical phenomenom for the finite element method with $C^0$ (and not $C^1$) elements: The discrete solutions $u_h$ converge in Sobolov space $H^{-1}(\varOmega)$ to the solution $u$ of the continuous problem when the mesh size $h$ approached $0$. That means that $u_h \to u$ in $L^2(\varOmega)$ and $Du_h \to Du$ in $L^2(\varOmega)$. But that does not mean that $\Delta u_h \to u$ in $L^2$ or even $L^\infty$. – Henrik Schumacher Feb 27 '18 at 18:10
• One has only convergence of $\Delta u_h$ to $\Delta u$ in $H^{-1}(\varOmega)$, that means in distributional sense: $\|\int_\varOmega \Delta_h(u_h - u) \, \varphi \, dx\| \colon= \|\int_\varOmega \langle D(u_h - u) \, D \varphi \rangle \, dx\| \leq C \, \varrho(h) \, (\|\varphi\|_{L^2}+\|D\varphi\|_{L^2})$ with $\varrho(h) \to 0$. In a nutshell: when integrated against a sufficiently smooth function $\varphi$, $\Delta_h u_h$ will converge to $\Delta_h u$. – Henrik Schumacher Feb 27 '18 at 18:10

The method outlined by Henrik Schumacher can also be using while explicitly relying on cylindrical symmetry (when applicable), so as to reach better accuracy.

## Laplace equation with a given boundary

    Needs["NDSolveFEM"];
dz = 1/32; reg0 = Rectangle[{{1/2, -dz}, {1, dz}}];
reg1 = RegionDifference[Disk[{0, 0}, 4],
Rectangle[{1/2, -dz}, {1, dz}]];
reg = RegionIntersection[reg1, ImplicitRegion[x > 0, {x, y}]];
mesh = ToElementMesh[reg, MaxCellMeasure -> 0.01,
"MaxBoundaryCellMeasure" -> 0.025];
Show[mesh["Wireframe"], Frame -> True, PlotRange -> All] sol = NDSolveValue[{-Laplacian[u[x, z], {x, y, z}, "Cylindrical"] ==
NeumannValue[1, 1/2^2 <= x^2 <= 1 && -dz <= z <= dz],
DirichletCondition[u[x, z] == 0, x^2 + z^2 >= 2^2]} // Evaluate,
u, {x, z} \[Element] mesh];
Off[InterpolatingFunction::dmval]


We can plot it as

plr = Region@ImplicitRegion[1/2^2 <= x^2 <= 1 && -dz < z < dz, {x, z}];
Clear[sol1]; sol1[x_, z_] = sol[Abs[x], z];
plc = ContourPlot[sol1[x, z], {x, -2, 2}, {z, -2, 2}, Contours -> 20,
AspectRatio -> 1, ColorFunction -> "Heat"];
pls = Show[
StreamPlot[{D[sol[Abs[x], z], x], D[sol[Abs[x], z], z]} /.
Abs'[x] -> If[x >= 0, 1, -1] // Evaluate, {x, -2, 2}, {z, -2,
2}, StreamStyle -> Gray], plr];
Show[plc, pls] ## Poisson equation with a given surface density

Similarly, adapting the 3D case to the cylindrical symmetry

reg = Disk[{0, 0}, 2];
reg = RegionIntersection[reg, ImplicitRegion[x > 0, {x, y}]];
mesh = ToElementMesh[reg, MaxCellMeasure -> 0.001];
\[Rho] = Function[{x, z},
Boole[1/2^2 <= x^2 <= 9/4 && -dz <= z <= dz]];
sol = NDSolveValue[{-Laplacian[u[x, z], {x, \[Theta], z},
"Cylindrical"] == \[Rho][x, z] // Evaluate,
DirichletCondition[u[x, z] == 0, x^2 + z^2 >= 2^2]},
u, {x, z} \[Element] mesh];
Off[InterpolatingFunction::dmval]
plc = ContourPlot[sol1[x, z], {x, -2, 2}, {z, -2, 2}, Contours -> 20,
AspectRatio -> 1, ColorFunction -> "Heat"];
pls = Show[
StreamPlot[{D[sol[Abs[x], z], x], D[sol[Abs[x], z], z]} /.
Abs'[x] -> If[x >= 0, 1, -1] // Evaluate, {x, -2, 2}, {z, -2,
2}, StreamStyle -> Gray], plr];
Show[plc, pls] We can check the accuracy of the solution as follows First compute the Laplacian of the solution

f[x_, z_] = -Laplacian[sol[x, z], {x, \[Theta], z},];


then

 {Plot[{f[x, 0.], \[Rho][x, 0.]}, {x, 0, 2}],
Plot[{f[3/4, x], \[Rho][3/4, x]}, {x, -1, 1}, PlotRange -> All]} and

{Plot[{f[x, 0.] - \[Rho][x, 0.]}, {x, 0, 2}],
Plot[{f[3/4, x] - \[Rho][3/4, x]}, {x, -1, 1}, PlotRange -> All]} • Nice plots! Btw.: The more I think about it: Using Neumann conditions (as well as Dirichlet conditions) on the boundary of solid matter is just wrong: Imagine a planet with uniform density but with mountains... and snow. Imagine that you could not have a sleigh ride down the hill for the gravitational force is always perpendicular to the underground. What a poor universe that would be! – Henrik Schumacher Feb 27 '18 at 9:10
• @HenrikSchumacher sure: it is not meant to be realistic but act as a basis for a decomposition of the interaction. – chris Feb 27 '18 at 9:27

I'm very sorry to spoil this party but solving PDE in the whole volume is a bad idea if all that you want is a force on one body. Use green function $V(\vec{x})=-\underset{\mathbb{R}^{3}}{\int}\frac{G}{|\vec{x}-\vec{r}|}\rho(\vec{r})dr$ to get potential and force exactly in points of the second body and calculate all you need.

• agree. That's what I mean here – andre314 Feb 27 '18 at 19:17
• The problem is that the frontier of the solid matter doesn't correspond to Neumann or Dirichlet conditions. The fem problem posed in this way can't give the right solution (independently of Mathematica). – andre314 Feb 27 '18 at 19:53
• I have asked a moderator to reopen the chatroom about this question, so that everybody can write messages. – andre314 Feb 27 '18 at 20:14
• One downvote for @Vselvood ! I don't understand. What's happening ? – andre314 Feb 27 '18 at 20:35
• @andre lol it's frozen again – Vsevolod A. Feb 28 '18 at 5:10