8 of 17 added 115 characters in body; added 1 character in body

Gravitational potential created by a thin disc using NDSolve

Context

I would like to compute the torque that a (thin) disc applies onto a ring. I.e. I 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.

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 and that the vertical gradient should be zero outside

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]

Mathematica graphics

Problems:

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

(ii) This is imposing fixed potential on the disk not fixed density

(iii) Strangely enough The code fails if I use

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

to impose a ring condition.

Question

How to implement the Potential created by a disc using NDSolve?

In a broader sense I think I am asking how can FEM methods deal with PDEs with boundaries at infinity?

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