1 of 17

Potential created by a disc using NDSolve

Context

I would like to compute the torque that a razor thin disc applies onto a razor thin ring.

Attempt

I have defined a domain

dom = ImplicitRegion[0 <= x <= 1 && 0 <= 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

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