5
$\begingroup$

I'm trying to numerically solve Poisson's equation for the following scenario: The potential inside a cylinder of radius R=1 and height H=2 with uniform charge density(which I'll set to 1). Poisson's equation, in cylindrical form, is apparently $$\frac{\partial}{\partial r}(r\frac{\partial V}{\partial r})+\frac{\partial^2 V}{\partial z^2} = -\frac{\rho}{\varepsilon}.$$

In this case the charge density is obviously a bunch of dirac deltas, which Mathematica cannot handle. How should I numerically solve the equation then?

I assume referencing this discussion on a similar problem in magnetostatics will be helpful. However, I was unable to understand the spherical boundary trick and I'm looking for some help there as well.

Thanks in advance!

Edit: I should note that I know an analytic solution is possible; in fact, I have constructed an approximate solution using the following code(edited; I made some typos in the first version):

(*Initialize*) R = 1; H = 2; (*Analytic solution:*) U[r_, z_] := 1/(2 Pi) (NIntegrate[R/Sqrt[ R^2 + r^2 - 2 R r Cos[\[Theta]] + (\[Zeta] - z)^2], {\[Theta], 0, 2 \[Pi]}, {\[Zeta], -H, H}] + NIntegrate[\[Rho]/ Sqrt[\[Rho]^2 + r^2 - 2 \[Rho] r Cos[\[Theta]] + (H - z)^2], {\[Rho], 0, R}, {\[Theta], 0, 2 \[Pi]}] + NIntegrate[\[Rho]/ Sqrt[\[Rho]^2 + r^2 - 2 \[Rho] r Cos[\[Theta]] + (H + z)^2], {\[Rho], 0, R}, {\[Theta], 0, 2 \[Pi]}]); (*Plotting the field over a discrete set of points:*) u = Table[U[r, z], {r, 0.05, R-0.05, 0.05}, {z, -H+0.1, H-0.1, 0.1}]; field = Table[{u[[nr]][[nz]] - u[[nr + 1]][[nz]], u[[nr]][[nz]] - u[[nr]][[nz + 1]]}, {nr, 1, 18}, {nz, 1, 38}]; ListVectorPlot[field]

Yields the following output:Output

Improve this question
$\endgroup$
8
  • 5
    $\begingroup$ Why should the uniform charge density become a bunch of dirac's ? $\endgroup$ Commented Aug 7 at 15:38
  • $\begingroup$ Take the side of the cylinder for example; the charge density should be $$\sigma \delta (r-R)~,~~|z|\le H.$$ That's what I mean by deltas. I am fully aware that Neumann conditions are most probably a better choice, but I have no idea how to deploy them in this scenario. Looking for some help here as well. @UlrichNeumann $\endgroup$ Commented Aug 8 at 4:01
  • $\begingroup$ I should add that I know an analytic approach is possible, but I tried plotting it on mathematica and it won't work, I'm assuming it's because I needed to invoke NIntegrate and that slowed down the calculation. I was able to analytically calculate the potential over a set of discrete points as a rough approximation, but I expect that there is a better solution. $\endgroup$ Commented Aug 8 at 4:11
  • $\begingroup$ Isn't the charge density a boundary condition at r==R? What is V, perhaps a Charge distribution? $\endgroup$ Commented Aug 8 at 6:43
  • 1
    $\begingroup$ You have seen the ElectrostaticsPDEComponent and the accompanying monograph? $\endgroup$ Commented Aug 8 at 7:23

1 Answer 1

Reset to default
3
$\begingroup$

modified (charged cylinder surface)

Assuming the charge distribution proposed by QP we can solve the problem numerically, but we have to use an approximation for DiracDelta

eps = .001 dirac = Function[x, Exp[-(x^2/(2 eps))]/Sqrt[2 Pi eps]]

To solve the pde we need additional DirichletCondition.

R = 1; H = 2; v = NDSolveValue[{D[V[r, z], {z, 2}] +1/r D[r D[V[r, z], r], r] == -dirac[R - r] Boole[-H < z < H] - (dirac[z - H] + dirac[z + H])Boole[0 < r < R], DirichletCondition[V[r, z] == 0, r == 5 R]}, V, {r, 0, 5 R}, {z, -3 H, 3 H}, Method -> {"FiniteElement"}] Plot3D[v[r, z], {r, 0, 5 R}, {z, -3 H, 3 H}, PlotRange -> {0, All}, AxesLabel -> {"r", "z", "V[r,z]"}, MeshFunctions -> (#3 &), MaxRecursion -> 4, PlotPoints -> 50]

enter image description here

VectorPlot[Evaluate@Grad [v[r, z], {r, z}], {r, 0, 5 R}, {z, -3 H, 3H}]

enter image description here

Hope it helps!

Improve this answer
$\endgroup$
10
  • $\begingroup$ I know this might be dumb to ask, but it seems to me that there's a problem with your solution: why would you assume that the potential on the axis is a constant, which is clearly not the case? $\endgroup$ Commented Aug 9 at 15:01
  • $\begingroup$ @JeffGiff With NeumannValue==0 at r==R && z==H,-H the solution doesn't depend on z . What are the boundary conditions you want to simulate? $\endgroup$ Commented Aug 9 at 16:07
  • $\begingroup$ I intend to simulate a cylinder of finite length as stated above, so I expect the solution to depend on z. $\endgroup$ Commented Aug 10 at 13:15
  • $\begingroup$ @JeffGiff Sorry for the delay. Have a look at my modified answer, now with finite cylinder and new DirichletCondition. $\endgroup$ Commented Aug 11 at 15:11
  • $\begingroup$ That was very helpful! I should note, however, that you forgot to take into account the top and bottom faces of the cylinder, but that doesn't really add much to the solution anyways. $\endgroup$ Commented Aug 13 at 4:16

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.