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