I am physicist from the first education, so it is apparently my field. As it follows from my experience in 3D FEM testing with application to the magnetic field calculation there is a problem with equation $\nabla \times (\nu \nabla \times \vec {A})=\vec {j}$. Therefore we prefer another form of this equation, for example, this one $\nabla \nu \times (\nabla \times \vec {A})+\nu \nabla \times \nabla \times\vec {A} =\vec {j}$ (ungauged form).
Then if we have Coulomb gauge $\nabla.\vec {A}$, it automatically turns to $\nabla \nu \times (\nabla \times \vec {A})-\nu \nabla ^2\vec {A} =\vec {j}$ (Coulomb gauge). Now we can compare two form using mesh from Tim Laska answer (thanks to him), and function appro from xzczd answer (thanks to him as well). Let check Coulomb gauge first:
u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]}; appro =
With[{k = 1. 10^4}, Tanh[k #]/2 + 1/2 &];
\[Nu]1 = Simplify`PWToUnitStep@
PiecewiseExpand@If[x <= a && y <= a && z <= a, 1/(\[Mu]r), 1] /.
UnitStep ->
appro;(*permeability depending on iron cube in mesh*)\
\[CapitalGamma]d = {DirichletCondition[ux[x, y, z] == 0, y == 0],
DirichletCondition[ux[x, y, z] == -fluxDensity*b/2, y == b],
DirichletCondition[uy[x, y, z] == 0, x == 0],
DirichletCondition[uy[x, y, z] == fluxDensity*b/2, x == b],
DirichletCondition[uz[x, y, z] == 0,
y == b || y == 0 || x == 0 || x == b || z == 0 || z == b]};
\[CapitalGamma]n = {0, 0, 0};
op7 = Cross[Grad[\[Nu]1, {x, y, z}], Curl[u, {x, y, z}]] - \[Nu]1*
Laplacian[u, {x, y, z}];(*Coulomb gauged*){mvpAx, mvpAy, mvpAz} =
NDSolveValue[{op7 == {0, 0, 0}, \[CapitalGamma]d}, {ux, uy,
uz}, {x, y, z} \[Element] mesh];
Visualization
Now let check ungauged form
u = {ux[x, y, z], uy[x, y, z], uz[x, y, z]}; appro =
With[{k = 2. 10^4}, ArcTan[k #]/Pi + 1/2 &];
\[Nu]1 = Simplify`PWToUnitStep@
PiecewiseExpand@If[x <= a && y <= a && z <= a, 1/(\[Mu]r), 1] /.
UnitStep ->
appro;(*permeability depending on iron cube in mesh*)\
\[CapitalGamma]d = {DirichletCondition[ux[x, y, z] == 0, y == 0],
DirichletCondition[ux[x, y, z] == -fluxDensity*b/2, y == b],
DirichletCondition[uy[x, y, z] == 0, x == 0],
DirichletCondition[uy[x, y, z] == fluxDensity*b/2, x == b],
DirichletCondition[uz[x, y, z] == 0,
y == b || y == 0 || x == 0 || x == b || z == 0 || z == b]};
\[CapitalGamma]n = {0, 0, 0};
op7 = Cross[Grad[\[Nu]1, {x, y, z}], Curl[u, {x, y, z}]] - \[Nu]1*
Laplacian[u, {x, y, z}]; op8 =
Cross[Grad[\[Nu]1, {x, y, z}], Curl[u, {x, y, z}]] + \[Nu]1*
Curl[Curl[u, {x, y, z}], {x, y, z}];(*Coulomb gauged*){mvpAx,
mvpAy, mvpAz} =
NDSolveValue[{op8 == {0, 0, 0}, \[CapitalGamma]d}, {ux, uy,
uz}, {x, y, z} \[Element] mesh];
It looks reasonable but pay attention how we play with k and with Tanh[](Coulomb gauge) and ArcTan[] (ungauged form).
For reference we can compare 3 solution for problem of coil magnetic field first considered by N. Demerdash, T. Nehl and F. Fouad, "Finite element formulation and analysis of three dimensional magnetic field problems," in IEEE Transactions on Magnetics, vol. 16, no. 5, pp. 1092-1094, September 1980. doi: 10.1109/TMAG.1980.1060817. This solution explained without code on https://physics.stackexchange.com/questions/513834/current-density-in-a-3d-loop-discretising-a-model/515657#515657
We have to calculate the vector potential and magnetic field of a rectangular coil with a current of 20A. The number of turns = 861. Inner cross section is 10.42cm×10.42cm, outer cross section is 15.24cm×15.24cm, coil height is 8.89cm. Here we show code for Closed Form Solution Algorithm (CFSA), BEM (Integral) and Mathematica FEM. CFSA code:
h = 0.0889; L1 = 0.1042; L2 = 0.1524; n = 861 (*16AWG wire*); J0 = \
20(*Amper*); j0 = 20*n/(h*(L2 - L1)/2); mu0 = 4 Pi 10^-7; b0 = j0 mu0;
bx[a_, b_, x_, y_, z_] :=
z/(Sqrt[(-a + x)^2 + (-b + y)^2 +
z^2] (-b + y + Sqrt[(-a + x)^2 + (-b + y)^2 + z^2])) - z/(
Sqrt[(a + x)^2 + (-b + y)^2 +
z^2] (-b + y + Sqrt[(a + x)^2 + (-b + y)^2 + z^2])) - z/(
Sqrt[(-a + x)^2 + (b + y)^2 +
z^2] (b + y + Sqrt[(-a + x)^2 + (b + y)^2 + z^2])) + z/(
Sqrt[(a + x)^2 + (b + y)^2 +
z^2] (b + y + Sqrt[(a + x)^2 + (b + y)^2 + z^2]))
by[a_, b_, x_, y_, z_] :=
z/(Sqrt[(-a + x)^2 + (-b + y)^2 +
z^2] (-a + x + Sqrt[(-a + x)^2 + (-b + y)^2 + z^2])) - z/(
Sqrt[(a + x)^2 + (-b + y)^2 +
z^2] (a + x + Sqrt[(a + x)^2 + (-b + y)^2 + z^2])) - z/(
Sqrt[(-a + x)^2 + (b + y)^2 +
z^2] (-a + x + Sqrt[(-a + x)^2 + (b + y)^2 + z^2])) + z/(
Sqrt[(a + x)^2 + (b + y)^2 +
z^2] (a + x + Sqrt[(a + x)^2 + (b + y)^2 + z^2]))
bz[a_, b_, x_, y_,
z_] := -((-b + y)/(
Sqrt[(-a + x)^2 + (-b + y)^2 +
z^2] (-a + x + Sqrt[(-a + x)^2 + (-b + y)^2 + z^2]))) - (-a +
x)/(Sqrt[(-a + x)^2 + (-b + y)^2 +
z^2] (-b + y + Sqrt[(-a + x)^2 + (-b + y)^2 + z^2])) + (-b + y)/(
Sqrt[(a + x)^2 + (-b + y)^2 +
z^2] (a + x + Sqrt[(a + x)^2 + (-b + y)^2 + z^2])) + (a + x)/(
Sqrt[(a + x)^2 + (-b + y)^2 +
z^2] (-b + y + Sqrt[(a + x)^2 + (-b + y)^2 + z^2])) + (b + y)/(
Sqrt[(-a + x)^2 + (b + y)^2 +
z^2] (-a + x + Sqrt[(-a + x)^2 + (b + y)^2 + z^2])) + (-a + x)/(
Sqrt[(-a + x)^2 + (b + y)^2 +
z^2] (b + y + Sqrt[(-a + x)^2 + (b + y)^2 + z^2])) - (b + y)/(
Sqrt[(a + x)^2 + (b + y)^2 +
z^2] (a + x + Sqrt[(a + x)^2 + (b + y)^2 + z^2])) - (a + x)/(
Sqrt[(a + x)^2 + (b + y)^2 +
z^2] (b + y + Sqrt[(a + x)^2 + (b + y)^2 + z^2]))
da = (L2 - L1)/15/2;
dh = h/26/2; a = b = L1/2;
Bz[x_, y_, z_] :=
Sum[bz[a + da (i - 1), b + da (i - 1), x, y, z + dh j], {i, 1,
16}, {j, -26, 26, 1}] +
Sum[bz[a, b, x, y, z + dh j], {j, -6, 6,
1}];
Code for BEM (Integral)
reg = RegionDifference[
ImplicitRegion[-L2/2 <= x <= L2/2 && -L2/2 <= y <= L2/2 && -h/2 <=
z <= h/2, {x, y, z}],
ImplicitRegion[-L1/2 <= x <= L1/2 && -L1/2 <= y <= L1/2 && -h/2 <=
z <= h/2, {x, y, z}]];
j[x_, y_, z_] := Boole[{x, y, z} \[Element] reg]
jx[x_, y_, z_] := If[-y <= x <= y || y <= -x <= -y, Sign[y], 0]
jy[x_, y_, z_] := -jx[y, x, z]
Bx1[X_?NumericQ, Y_?NumericQ, Z_?NumericQ] :=
b0/(4 Pi) NIntegrate[
j[x, y, z] jy[x, y,
z] (Z - z)/(Sqrt[(x - X)^2 + (y - Y)^2 + (z - Z)^2])^3, {x, y,
z} \[Element] reg] // Quiet
By1[X_?NumericQ, Y_?NumericQ,
Z_?NumericQ] := -b0/(4 Pi) NIntegrate[
j[x, y, z] jx[x, y,
z] (Z - z)/(Sqrt[(x - X)^2 + (y - Y)^2 + (z - Z)^2])^3, {x, y,
z} \[Element] reg] // Quiet
Bz1[X_?NumericQ, Y_?NumericQ, Z_?NumericQ] :=
b0/(4 Pi) NIntegrate[
j[x, y, z] (jx[x, y, z] (Y - y) -
jy[x, y,
z] (X - x))/(Sqrt[(x - X)^2 + (y - Y)^2 + (z - Z)^2])^3, {x,
y, z} \[Element] reg] // Quiet
Code for FEM
eq1 = {Laplacian[A1[x, y, z], {x, y, z}] == -j[x, y, z] jx[x, y, z],
Laplacian[A2[x, y, z], {x, y, z}] == -j[x, y, z] jy[x, y, z]};
{Ax1, Ay1} =
NDSolveValue[{eq1,
DirichletCondition[{A1[x, y, z] == 0, A2[x, y, z] == 0},
True]}, {A1, A2}, {x, y, z} \[Element]
ImplicitRegion[-2 L2 <= x <= 2 L2 && -2 L2 <= y <=
2 L2 && -2 L2 <= z <= 2 L2, {x, y, z}]];
B = Evaluate[Curl[{Ax1[x, y, z], Ay1[x, y, z], 0}, {x, y, z}]];
Now we compute and visualize data
lst1 = Table[{z1, -b0 B[[3]] /. {x -> 0, y -> 0,
z -> z1}}, {z1, -.3, .3, .01}];
lst2 = Table[{z1, Bz[0, 0, z1] mu0 20/(4 Pi)}, {z1, -.3, .3, .01}];
lst3 = Table[{z1, -Bz1[0, 0, z1]}, {z1, -.3, .3, .01}];
{Region[reg],
Show[ListLinePlot[lst2, PlotStyle -> Orange, Frame -> True,
Axes -> False],
ListPlot[{lst1, lst2, lst3}, Frame -> True,
FrameLabel -> {"z", "\!\(\*SubscriptBox[\(B\), \(z\)]\)"},
PlotLegends -> {"FEM", "CFSA", "Integral"}]]}
