Episode #125 of the Stack Overflow podcast is here. We talk Tilde Club and mechanical keyboards. Listen now
16 of 17 added 96 characters in body

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.

Mathematica graphics

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)

Mathematica graphics

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]

Mathematica graphics

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

Mathematica graphics

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[[1]]^2 + fD[[2]]^2] // FullSimplify, fD[[3]]};

Mathematica graphics

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

Mathematica graphics

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]

Mathematica graphics