8
$\begingroup$

Bug introduced in 11.2.0 and fixed in 11.3.0

The system is a hollow cylinder (thin solenoid) with a current density $\text{J}$ and I'm looking to solve the magnetic potential ($\text{A}$) inside the radius of the cylinder using FEA. I'm struggling with the boundary conditions for the open ends (and in general). Suggestions are appreciated. 

u0 = 4 Pi 1*^-7; j = 100; (*current density parallel to z-axis*) solution = NDSolveValue[{ -Div[(1/u0) Grad[u[r, z], {r, th, z}, "Cylindrical"], {r, th, z}, "Cylindrical"] == NeumannValue[j, r == 1], DirichletCondition[u[r, z] == 0, z^2 == 1]}, u, {r, 0, 1}, {z, -1, 1}, Method -> {"PDEDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxBoundaryCellMeasure" -> 0.01}}}] ContourPlot[solution[r, z], {r, z} ∈ solution["ElementMesh"], PlotRange -> All, Contours -> 19, ImageSize -> Medium]

To compute the magnetic density field ($\text{B}$):

bField = Curl[solution[r, z], {r, th, z}, "Cylindrical"]

Update:

I've done additional research and the following is an approach by FEMM to solve such problems. The cross section is still in cylindrical coordinates but the open boundary is spherical.

bound = Table[{boundRadius Sin[theta], boundRadius Cos[theta]}, {theta, 0, Pi, Pi/100}]; boundIndex = Partition[Last[FindShortestTour[bound]], 2, 1]; coil = {{1, -1}, {1, 1}}; coilIndex = {{Length[bound] + 1, Length[bound] + 2}}; boundRadius = 2; reg = ImplicitRegion[r^2 + z^2 < boundRadius^2 && r >= 0, {r, z}]; bmesh = ToBoundaryMesh["Coordinates" -> Join[bound, coil], "BoundaryElements" -> {LineElement[Join[boundIndex, coilIndex]]}]; mesh = ToElementMesh[bmesh]; Show[ RegionPlot[reg, Epilog -> Point[bound], AspectRatio -> 2, ImageSize -> Medium], bmesh["Wireframe"] ]

enter image description here

I need to apply the following mixed boundary condition to the spherical surface:

$\frac{1}{\mu_0}\frac{\partial A}{\partial r}+\frac{A}{\mu_0 R} = 0$

I have the following failed attempt (NDSolveValue spits out bcnop warning and the internal b.c. is ignored):

u0 = 4 Pi 1*^-7; j = 100; solution = NDSolveValue[{-Div[(1/u0) Grad[u[r, z], {r, th, z}, "Cylindrical"], {r, th, z}, "Cylindrical"] - j == NeumannValue[0, r == 0] - NeumannValue[u[r, z]/(boundRadius u0), r^2 + z^2 == boundRadius^2], DirichletCondition[u[r, z] == u0 j, r == 1 && z^2 <= 1]}, u, {r, z} ∈ mesh, Method -> {"PDEDiscretization" -> {"FiniteElement"}}] ContourPlot[solution[r, z], {r, z} ∈ solution["ElementMesh"], PlotRange -> All, Contours -> 19, ImageSize -> Medium, AspectRatio -> 2]
Improve this question
$\endgroup$
12
  • 3
    $\begingroup$ Just a precison for beginners in magnetostatic : this problem can be solved anatically. See here $\endgroup$ Commented Sep 22, 2017 at 18:27
  • $\begingroup$ That's what I was thinking. My message was for "visitors" $\endgroup$ Commented Sep 22, 2017 at 18:30
  • $\begingroup$ The problem of open ends (=boundaries at infinity) arrises very often in magnetostatics and electrostatics. Example : In electrostatics in 2 dimensions, on may use a conformal mapping to map infinity to finite boundaries, and then use FEM. Otherwise, FEM solvers have special complex methods for dealing with some "infinite mesh elements" (?). I'm afraid that's not a simple affair. Would be happy to be wrong. $\endgroup$ Commented Sep 22, 2017 at 20:03
  • 2
    $\begingroup$ OK, @xzczd there is a regression in V11.2 that triggered the message and that will need to get fixed. Never the less I think the wire should have some extension. $\endgroup$ Commented Oct 1, 2017 at 12:07
  • 1
    $\begingroup$ @xzczd, I just merged a fix for this problem, the next release will work as expected. Thanks for digging this out. $\endgroup$ Commented Oct 17, 2017 at 0:57

4 Answers 4

Reset to default
6
$\begingroup$

This is fixed in version 11.3:

Needs["NDSolve`FEM`"] boundRadius = 2; bound = Table[{boundRadius Sin[theta], boundRadius Cos[theta]}, {theta, 0, Pi, Pi/100}]; boundIndex = Partition[Last[FindShortestTour[bound]], 2, 1]; coil = {{1, -1}, {1, 1}}; coilIndex = {{Length[bound] + 1, Length[bound] + 2}}; reg = ImplicitRegion[r^2 + z^2 < boundRadius^2 && r >= 0, {r, z}]; bmesh = ToBoundaryMesh["Coordinates" -> Join[bound, coil], "BoundaryElements" -> {LineElement[Join[boundIndex, coilIndex]]}]; mesh = ToElementMesh[bmesh]; Show[RegionPlot[reg, Epilog -> Point[bound], AspectRatio -> 2, ImageSize -> Medium], bmesh["Wireframe"]]

enter image description here

Show[Graphics[ GraphicsComplex[mesh["Coordinates"], Point[Flatten[ElementIncidents[mesh["PointElements"]]]]]]]

enter image description here

Note, how the interior points are now (again) part of the boundary mesh. And the solution:

u0 = 4 Pi 1*^-7; j = 100; solution = NDSolveValue[{-Div[(1/u0) Grad[u[r, z], {r, th, z}, "Cylindrical"], {r, th, z}, "Cylindrical"] - j == NeumannValue[0, r == 0] - NeumannValue[u[r, z]/(boundRadius u0), r^2 + z^2 == boundRadius^2], DirichletCondition[u[r, z] == u0 j, r == 1 && z^2 <= 1]}, u, {r, z} \[Element] mesh]; ContourPlot[solution[r, z], {r, z} \[Element] solution["ElementMesh"], PlotRange -> All, Contours -> 19, ImageSize -> Medium, AspectRatio -> 2, PlotPoints -> 50]

enter image description here

Improve this answer
$\endgroup$
8
+50
$\begingroup$

Well, I'm not quite familiar with electromagnetism and it's not immediately clear to me how to compare the numeric solution with the analytic solution, but I've circumvented the bcnop warning anyway, so let me post an answer.

The idea is creating a very narrow slit in the domain to approximate the wire i.e. the LineElement in your failed attempt. Notice I've also enlarge boundRadius because after some trial I noticed boundRadius = 2 is too small:

boundRadius = 4; u0 = 4 Pi 1*^-7; j = 100; eps = 10^-4; reg = ImplicitRegion[ r^2 + z^2 < boundRadius^2 && r >= 0 && ! (1 - eps < r < 1 + eps && z^2 < 1), {r, z}]; solution = NDSolveValue[{-Div[(1/u0) Grad[u[r, z], {r, th, z}, "Cylindrical"], {r, th, z}, "Cylindrical"] - j == -NeumannValue[u[r, z]/(boundRadius u0), r^2 + z^2 == boundRadius^2], DirichletCondition[u[r, z] == u0 j, 1 - eps < r < 1 + eps && z^2 < 1]}, u, {r, z} ∈ reg, Method -> {FiniteElement, MeshOptions -> MaxCellMeasure -> 0.001}] interestedRadius = 2; interestedreg = ImplicitRegion[ r^2 + z^2 < interestedRadius^2 && ! (1 - eps < r < 1 + eps && z^2 < 1), {r, z}]; ContourPlot[ Evaluate@PiecewiseExpand@If[r > 0, solution[r, z], solution[-r, z]], {r, z} ∈ interestedreg, PlotRange -> All, PlotPoints -> 100, Contours -> 19]

Mathematica graphics

Improve this answer
$\endgroup$
10
  • $\begingroup$ I think this is the way to go - also from a physics point of view. You may want to create the mesh manually like in the question for very 'thin' wires. $\endgroup$ Commented Oct 1, 2017 at 8:41
  • $\begingroup$ I eventually obtained a result close to your graph but when you take the Curl of solution to obtain the magnetic field it doesn't match the analytical solution. $\endgroup$ Commented Oct 1, 2017 at 17:32
  • $\begingroup$ I think there is something wrong with how I formed the differential equation. $\endgroup$ Commented Oct 1, 2017 at 17:39
  • $\begingroup$ @Young Er…sorry, but the $A$ i.e. u in your equation is a scalar, how can you calculate its Curl? $\endgroup$ Commented Oct 2, 2017 at 2:50
  • 1
    $\begingroup$ The use of PiecewiseExpand seems unnecessary; replace Evaluate@PiecewiseExpand@If[r > 0, solution[r, z], solution[-r, z]] by Evaluate@solution[Abs@r, z] $\endgroup$ Commented Oct 7, 2021 at 11:09
6
$\begingroup$

Extended comment

The result of the OP's code are not the same on Mathematica 11.0 and Mathematica 11.2

Here is the code I have used :

<<NDSolve`FEM` bound = Table[{boundRadius Sin[theta], boundRadius Cos[theta]}, {theta, 0, Pi, Pi/100}]; boundIndex = Partition[Last[FindShortestTour[bound]], 2, 1]; coil = {{1, -1}, {1, 1}}; coilIndex = {{Length[bound] + 1, Length[bound] + 2}}; boundRadius = 2; reg = ImplicitRegion[r^2 + z^2 < boundRadius^2 && r >= 0, {r, z}]; bmesh = ToBoundaryMesh["Coordinates" -> Join[bound, coil], "BoundaryElements" -> {LineElement[Join[boundIndex, coilIndex]]}]; mesh = ToElementMesh[bmesh]; (*Show[ RegionPlot[reg, Epilog -> Point[bound], AspectRatio -> 2, ImageSize -> Medium], bmesh["Wireframe"] ]*) u0 = 4 Pi 1*^-7; j = 100; solution = NDSolveValue[{-Div[(1/u0) Grad[u[r, z], {r, th, z}, "Cylindrical"], {r, th, z}, "Cylindrical"] - j == NeumannValue[0, r == 0] - NeumannValue[u[r, z]/(boundRadius u0), r^2 + z^2 == boundRadius^2], DirichletCondition[u[r, z] == u0 j, r == 1 && z^2 <= 1]}, u, {r, z} \[Element] mesh, Method -> {"PDEDiscretization" -> {"FiniteElement"}}] ContourPlot[solution[r, z], {r, z} \[Element] solution["ElementMesh"], PlotRange -> All, Contours -> 19, ImageSize -> Medium, AspectRatio -> 2]

Mathematica 11.2

enter image description here

Mathematica 11.0 - first evaluation (on a fresh kernel)

It doesn't work :

enter image description here

Mathematica 11.0 - second evaluation (= reevaluation of the same code)

enter image description here

Edit

same as before, but with PlotPoints->50added to ContourPlot:

enter image description here

Improve this answer
$\endgroup$
8
  • $\begingroup$ Will the last result be improved if you add e.g. PlotPoints->50 to ContourPlot? $\endgroup$ Commented Oct 1, 2017 at 10:55
  • 1
    $\begingroup$ @xzczd see my edit. Nice ! $\endgroup$ Commented Oct 1, 2017 at 11:07
  • 1
    $\begingroup$ It helps if you define the boundRadius = 2; before the bound. $\endgroup$ Commented Oct 1, 2017 at 11:54
  • 1
    $\begingroup$ OK, thanks Andre (@xzczd) this is a regression in 11.2 and I'll see to get it fixed for the next release. $\endgroup$ Commented Oct 1, 2017 at 12:05
  • 1
    $\begingroup$ A precision for visitors : the comment of @user21 above (boundRadius =2 before the bound) is destinated to fix the problem of the failing first evaluation on Mathematica version 11.0 $\endgroup$ Commented Oct 1, 2017 at 12:08
4
$\begingroup$

Extended comment:

Something like this

epsi = 10^-1; coil = {{1 - epsi, -1}, {1 - epsi, 1}, {1 + epsi, 1}, {1 + epsi, -1}}; coilIndex = {{1, 2}, {2, 3}, {3, 4}, {4, 1}} + Length[bound]; boundRadius = 2; reg = ImplicitRegion[r^2 + z^2 < boundRadius^2 && r >= 0, {r, z}]; Show[ bmesh["Wireframe"], Graphics[{Blue, Point[{1, 0}]}], Graphics[ GraphicsComplex[ bmesh["Coordinates"], {Red, Point[Flatten[ElementIncidents[bmesh["PointElements"]]]]}]] ]
Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.