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$
    – andre314
    Sep 22, 2017 at 18:27
  • $\begingroup$ That's what I was thinking. My message was for "visitors" $\endgroup$
    – andre314
    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$
    – andre314
    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$
    – user21
    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$
    – user21
    Oct 17, 2017 at 0:57

4 Answers 4

Reset to default
5
$\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$
    – user21
    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$
    – Young
    Oct 1, 2017 at 17:32
  • $\begingroup$ I think there is something wrong with how I formed the differential equation. $\endgroup$
    – Young
    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$
    – xzczd
    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$
    – Ronald Monson
    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$
    – xzczd
    Oct 1, 2017 at 10:55
  • 1
    $\begingroup$ @xzczd see my edit. Nice ! $\endgroup$
    – andre314
    Oct 1, 2017 at 11:07
  • 1
    $\begingroup$ It helps if you define the boundRadius = 2; before the bound. $\endgroup$
    – user21
    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$
    – user21
    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$
    – andre314
    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$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.