5 of 17 added 161 characters in body

Gravitational potential created by a razor thin disc using NDSolve

Context

I would like to compute the torque that a (razor 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 need to compute the gravitational potential generated by a (razor thin) disc and or 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 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?

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