Computed B field too high for a solenoid

[themagpipers]

I would like to compute the B field of a solenoid (and then compute other derived quantities such as the inductance, etc.). So far, I am getting a B field's value that's too high compared to the theoretical one, by a factor of roughly 7. Note that I am aware of the default units for lengths (mm), B fields (milli Tesla) and H fields ( kA /m), which all differ from SI units.

Since there is no solenoid object in magpylib, I decided to create one such object using a parametrization of current.Line()

Here is my code:

import magpylib as mag
import numpy as np

# Vacuum permeability.
µ0 = 1.256637e-6

# Length (height) of the coil, in mm.
L = 150
n_winding = 400

# Resolution of the coil. The right int specifies how many points are going to be
# used for every loop of the solenoid.
t_res = n_winding * 10
t = np.linspace(0, L, t_res)
ω = 2 * np.pi * n_winding / L
radius = 20

# Parametric definition of the vertices making up the coil.
vertices = np.asarray([radius * np.cos(ω * t), radius * np.sin(ω * t), t]).T
current = 1.0

# Define the solenoid as a magpylib object.
solenoid = mag.current.Line(current, vertices)

# Move it to the center.
solenoid.move((0,0,-L/2))

# Define the dimensions of the box in millimeters
dx, dy, dz = 100, 100, 300

# Define the number of cells in each dimension
nx, ny, nz = 11, 11, 15

# Create the grid
grid = pv.ImageData(
    dimensions=(nx+1, ny+1, nz+1),  # Dimensions of the box. That's the number of points in each direction.
    spacing=(dx/nx, dy/ny, dz/nz),  # Spacing between cells.
    origin=(-dx/2, -dy/2, -dz/2),  # The bottom left corner of the data set
)

grid["B"] = solenoid.getB(grid.points)
grid["H"] = solenoid.getH(grid.points)
print(f'Maximum computed B field: {np.max(grid["B"]) / 1e3} T.')
print(f"Expected B field's magnitude inside solenoid: {µ0 * n_winding * current} T.")

# Optional code below to visualize the B field, to make sure the B max computed is not an isolated bad point.
pl = pv.Plotter()
pl.add_mesh(grid.outline(), color="k")
mag.show(solenoid, canvas=pl, backend="pyvista")
# Add the B field vectors to the scene
pl.add_mesh(grid.glyph(scale="B", orient="B", factor=np.max(grid["B"])*1))
pl.show()
pl.close()

So, for a coil 15 cm long, with a radius of 2 cm, 400 windings, the expected B field in the solenoid is about 5e-4 T (as given by the formula shown in the code). However, the computed value is 3.3e-3 T. This is significantly higher than the theoretical value. I get the same problem if I refine the grid where the B field is computed, and also if I refine the solenoid by having a much higher point density.

I would be glad to know if I have done something wrong, or if anyone has an idea about the reason I am getting such a high B field computed compared to the theoretical solution.

[OrtnerMichael]

hi @themagpipers

You forgot the 1/L in the expected field - this will solve the problem.

Also: While your parametrization of the coil is nicely done, you also have 4000 individual Line segments. This is not a high number and for a circular coil you should probably choose t_res=n_winding * 100. I guess you chose 10 because it makes no big difference and because you realize that a large t_res makes the computation slow.

I suggests, that you use instead the Loop class. In this case you have a perfect circle, and only n_winding computations internally. Something like this:

import numpy as np
import magpylib as mag

current = 1.0
L = 150
n_winding = 400
radius = 20

coil = mag.Collection()
for z in np.linspace(-L/2, L/2, n_winding):
    coil.add(mag.current.Loop(current, 2*radius, (0,0,z)))

print(coil.getB((0,0,0))[2])

gives pretty much the same result as your computation - and is probably a bit faster, especially when computing for larger arrays of positions.

You can also use the magnet Cylinder class to describe the field of solenoids. Here the outer cylinder wall functions like a current sheet.

I will add a documentation example of all this - I think this is a really nice example reproducing known physics limits !

br,
michael

[themagpipers]

Thank you very much @OrtnerMichael ! I totally overlooked this factor, despite checking a couple of times. I now get a computed field lower than the one from the analytical formula, which makes sense, I guess (5 mT vs 4.32 mT for a given configuration). This is in much better agreement with the theory.

I am aware I have a very high number of Line segments, I was thinking to use a collection (but apparently this wouldn't change anything related to the speed of the calculations). I used many different resolutions (and similarly for the grid). I didn't know I could use a Loop to create a solenoid. But indeed, I have hit a performance wall quite early. Not only does my code takes several minutes to run when the grid contains as little as 10 k elements, but the RAM usage goes over 32 GB and I have to abort the process (as far as I remember, I wouldn't face this problem with finite elements with a mesh with that many elements).

I will try your suggestion to use a Collection of Loops. Even though it looks like technically this is different from a solenoid, this may give very similar results, but possibly not quite the same either.

Thanks for all.