The Stack Overflow podcast is back! Listen to an interview with our new CEO.

12 of 17 deleted 54 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 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)

## 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]


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 the ring-like boundary condition instead:

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


## Question

How to implement the 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?

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


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]