WhitneyAVSolver unphysical B field when currents "close" to steel

[tpgillam]

We're finding that we get unphysical magnetic fields when we have current densities "near" to regions of high magnetic permeability.

To demonstrate this I've setup a small toy example of a 3x1x1 cuboid of 'steel' (relative permeabilty of 5000), surrounded by a closed coil model. The gap between the steel and the coil region is set to be:
- 0.0 ("solenoid_0") (in this case there should be a conformal boundary between the coil and steel)
- 0.1 ("solenoid_1")
- 0.5 ("solenoid_5")

In the first two cases, we see vortices in the B field within the steel region. This doesn't make sense, as the current density should be zero there (and the material is linear, so we should have curl B = curl H = 0, i.e. no vortices). The third example behaves sensibly again.

Other notes on the setup:
- I'm using the CoilSolver to get the current density in the coil. AFAICT this is working great, and the resulting J field seems sensible when I inspect in paraview.
- The mesh is coarse at the edge of the airbox, and gets finer towards the bodies in the centre.


- The .zip files contain the full setup, including meshes & solver log & output VTU file, but for quite a coarse mesh to keep file sizes down.

(Seems like I can only upload three attachments, so the following will hopefully appear in replies)
- The .png files are screenshots from paraview with a much finer mesh, showing the B field.
- The case.sif fies are the same in all examples, so I've attached that separately in case it is useful.


Can anyone spot anything we might be doing wrong with the setup? Or have ideas of things we can try to improve the situation?

Thanks!
Tom

[tpgillam]

Screenshots from paraview mentioned above:

airgap = 0.0

solenoid_0.png (346.01 KiB) Viewed 289 times

airgap = 0.1

solenoid_1.png (425.89 KiB) Viewed 289 times

airgap = 0.5

solenoid_5.png (712.76 KiB) Viewed 289 times

[tpgillam]

And an example case.sif:

[tpgillam]

Just to add, I've tried a couple more things:
- use a uniform mesh size (i.e. just as fine at the edge of the airbox as on the surface of the material regions)
- use a hexahedral rather than tetrahedral mesh

In both cases I see the same vortex pattern as in the images above. I wasn't necessarily expecting that to have changed anything, but just trying out anything I can think of!

[raback]

Hi

Do you have a reference solution?

Last year around this time we looked at an even more simple case that intuitively didn't make sense. It was a case of cocentric cylinders with a jump in reluctivity. In the end we could find no problem in Elmer. We compared to an analytical model where it was easy to forget the BC between the cylinders. When enforcing continuity of H_t over the interface agreement with Elmer was spot on.

-Peter

[tpgillam]

Hi Peter,

Thanks very much for the reply. Unfortunately we don't have a theoretical or other reference right now.

I've been poking around a little more at the `jfix` calculation. I think this might be the source of the issue, since actually the $g_s$ that's being calculated right now does have non-negligible terms within the air & box.

Attached is a plot I made in paraview for the airgap=0.1 case, showing the gradient of the `jfix` field (which should be $\nabla Q=g_u - g_s$ if I understand the docs correctly). It's zoomed in on one corner of the coil:

jfix_gradient.png (522.85 KiB) Viewed 255 times

I think what I want to do is to add B.C.s that ensure that $\nabla Q = 0$ outside of the coil body (since it is closed then I would have thought it should be possible to construct currents that satisfy this. But I'm not sure about the details of where the source currents are defined, and shape functions used etc., so maybe this is tricky)

However, the only lever I've found is to set `jfix = 0` on certain boundaries -- if I do this on either the surface of the box, or on the surface of the coil, then the resultant WhitneyAV solve diverges. Similarly if I turn off the `Fix Input Currents` step, then it diverges.

Do you know if it's possible to constrain `jfix` in the way that we'd like here?

[raback]

Hi,

Ok, maybe the Jfix potential has a support that extends to the steel. On the other hand you cannot have a solution to the problem unless the divergence of current density is very closely zero.

The alternative is to use elemental coil current in CoilSolver and request the "Fix Input Current Density" in that solver (and have false at the other). For elemental currents this may be done elementwise without introducing a potential with such a wide support.

You can also say "Use Elemental CoilCurrent = True" for the AV solver and skip the body force (or add all components).

-Peter

[tpgillam]

Hi Peter,

Thanks very much for the pointer -- I've been experimenting with this option this morning and it helps a lot! With a little fiddling (*), I now get sensible B fields out:

solenoid_0_happier.png (522.52 KiB) Viewed 226 times

Tom

(*) initially after making this change I found that the WhitneyAV solve diverged. By adding some filleting to the winding mesh, to smooth out the corners, and slightly decreasing the mesh size, I got convergence again. I'm not _too_ surprised that forcing current to go around a sharp corner might cause some issues.

However I did find some slightly surprising failures on the way, for example the first non-linear iteration would converge nicely, then the second iteration -- which should be trivial since the material is linear -- diverges. I'm guessing there was something a bit nasty that meant that the scaling step, or similar, was causing issues.

[raback]

Hi

Great! You almost had me worried. Maybe you could share the final outcome at the elmer-elmag repo.

-Peter

[tpgillam]

Thanks for the suggestion! I made a PR here: https://github.com/ElmerCSC/elmer-elmag/pull/21

Tom