TEAM Problem 13

2. TEAM Problem 13#

2.1. Getting started#

The script dependencies are:

ngsolve
numpy
scipy
matplotlib
the .vol file of the geometry ./team13_mesh.vol

If the geometry file is erroneous, execute the python script \(\texttt{./geometry.py}\), which generates the necessary \(\texttt{.vol}\) file.

python3 geometry.py [-fullProblem True/False]

2.2. Problem Description#

This script solves the T.E.A.M problem 13 (3-D Non-Linear Magnetostatic Model). The problem description can be found under the link.

The problem is based on the Maxwell equations with the corresponding boundary conditions

\[\begin{split} \begin{align} \nabla \times \textbf{H} &= \textbf{J}, \qquad &&\textbf{H}\times\textbf{n} = \textbf{0}\\ \nabla \cdot \textbf{B} &= 0, &&\textbf{B}\cdot\textbf{n} = 0 \end{align} \end{split}\]

for magnetostatic problems. The second equation allows the introduction of a magnetic vector potential \(\textbf{A}\) as

\[ \begin{equation} \textbf{B} = \nabla \times \textbf{A}. \end{equation} \]

For the weak formulation

\[ \begin{equation} \int_\Omega \frac{1}{\mu(|\textbf{B}|)}\nabla\times \textbf{A}\cdot \nabla\times\textbf{v}\; d\Omega = \int_{\Omega_c} \textbf{J}\cdot \textbf{v}\;d\Omega \end{equation} \]

is used, which considers the non-linear material relationship

\[ \begin{equation} \textbf{H} = \frac{1}{\mu(|\textbf{B}|)}\textbf{B}. \end{equation} \]

Here, \(\textbf{A}\) and \(\textbf{v}\) are the trial and the test function, \(\Omega\) is the domain of interest and \(\Omega_c\) is the domain of the coil.

For a unique solution, an additional regularisation term with a small \(\varepsilon > 0\) has been added:

\[ \begin{equation} \int_\Omega \frac{1}{\mu(|\textbf{B}|)}\nabla\times \textbf{A}\cdot \nabla\times\textbf{v}\; d\Omega + \int_\Omega\varepsilon\;\textbf{A}\cdot\textbf{v}\;d\Omega= \int_{\Omega_c} \textbf{J}\cdot \textbf{v}\;d\Omega \end{equation} \]

The magnetic flux density of the solution is presented in Fig. 1.

title

Fig. 1: Illustration of the the magnetic flux density in the steel channels and magnetic field lines.

2.3. Coding#

2.3.1. Imports#

Importing the packages \(\textit{ngsolve, netgen, numpy, ... }\) enables all functionalities of this script. Additionally, the number of used threads is reduced and the permeability of vacuum \(\mu_0\) is set.

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.csg import *
from ngsolve.internal import *

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import time
from geometry import generateGeometry, getMeshingParameters
%matplotlib widget
SetNumThreads(10)

mu0 = 4*np.pi*1e-7

2.4. Setting Arguments#

First of all, the simulation parameters need to be set. The space order defines the used order of the trial function and the test function within the finite element space. The variable \(\texttt{I}_{ON}\) defines the value of Ampere-turns in the coil. Considering the description, the value for the Ampere-turns is set to \( 1,000 \) or \( 3,000 \).

space_order = 2   
I0N = 1000                  # Ampere-turns 1000 or 3000

2.5. Meshing#

Next, the available mesh is loaded. The command \(Curve(5)\) determines the roundness of the mesh and is obligatory for the provided geometry with its curved corners.

if True:
    geo = generateGeometry(maxh_iron = 100, fullProblem=False)
    mp = getMeshingParameters(maxh_global = 100)
    ngmesh = geo.GenerateMesh(mp=mp) # global maxh
    mesh = Mesh(ngmesh)
else:
    mesh = Mesh("./team13_mesh.vol")
    # mesh = Mesh("./team13_mesh_full.vol") # activate this to simulate the full problem
mesh.Curve(5)

# check domains
val = {"corner_right_back":1 , "corner_left_back":1, "corner_left_front":1, "corner_right_front":1, "brick_front":1, "brick_back":1, "brick_left":1, "brick_right":1, "iron":2, "air":1.5}
domains = mesh.MaterialCF(val, default=None)

To check the correct assignment of the domains \(\{air, iron, coil\}\), all corresponding subdomains are coloured accordingly.

title

Fig. 2: illustration of the geometry In Fig. 2 the coil is blue, the steel channels are red and the air is transperent. For a meaningful display a clipping plane was inserted.

2.5.1. Finite Element Space#

The finite element space \(H(curl, \Omega)\) is selected and homogeneous Dirichlet boundary conditions are prescribed on the far boundary. Additionally

the solution, i.e. magnetic vector potential \(\textbf{A}\), is defined as \(\texttt{sol}\),
trial and test function are introduced,
the \(\textit{CoefficientFunction}\) \(\texttt{B}\) for the magnetic flux density is determined and
\(\texttt{B}\) is drawn.

# create fe space
fes = HCurl(mesh, order=space_order, dirichlet="outer", nograds=True)

# magnetic vector potential as GridFunction
sol = GridFunction(fes, "A")
sol.vec[:] = 0

# Test- and Trialfunction
u = fes.TrialFunction()
v = fes.TestFunction()

# magnetic flux density
B = curl(sol)

2.5.2. Material Parameters#

The non-linear magnetisation curve of the material is defined in the description .

H_KL = [ -4.47197834e-13, 1.60000000e+01, 3.00000000e+01, 5.40000000e+01\
    , 9.30000000e+01, 1.43000000e+02, 1.91000000e+02, 2.10000000e+02\
, 2.22000000e+02, 2.33000000e+02, 2.47000000e+02, 2.58000000e+02\
, 2.72000000e+02, 2.89000000e+02, 3.13000000e+02, 3.42000000e+02\
, 3.77000000e+02, 4.33000000e+02, 5.09000000e+02, 6.48000000e+02\
, 9.33000000e+02, 1.22800000e+03, 1.93400000e+03, 2.91300000e+03\
, 4.99300000e+03, 7.18900000e+03, 9.42300000e+03, 9.42300000e+03\
, 1.28203768e+04, 1.65447489e+04, 2.07163957e+04, 2.55500961e+04\
, 3.15206135e+04, 4.03204637e+04, 7.73038295e+04, 1.29272791e+05\
, 1.81241752e+05, 2.33210713e+05, 2.85179674e+05, 3.37148635e+05\
, 3.89117596e+05, 4.41086557e+05, 4.93055518e+05, 5.45024479e+05\
, 5.96993440e+05, 6.48962401e+05, 7.00931362e+05, 7.52900323e+05\
, 8.04869284e+05, 8.56838245e+05, 9.08807206e+05, 9.60776167e+05\
, 1.01274513e+06, 1.06471409e+06, 1.11668305e+06, 1.16865201e+06\
, 1.22062097e+06, 1.27258993e+06, 1.32455889e+06, 1.37652785e+06\
, 1.42849682e+06, 1.48046578e+06, 1.53243474e+06, 1.58440370e+06\
, 1.63637266e+06, 1.68834162e+06, 1.74031058e+06, 1.79227954e+06\
, 1.84424850e+06, 1.89621746e+06, 1.94818643e+06, 2.00015539e+06\
, 2.05212435e+06, 2.10409331e+06, 2.15606227e+06, 2.20803123e+06\
, 2.26000019e+06]

B_KL = [  0.00000000e+00, 2.50000000e-03, 5.00000000e-03, 1.25000000e-02\
, 2.50000000e-02, 5.00000000e-02, 1.00000000e-01, 2.00000000e-01\
, 3.00000000e-01, 4.00000000e-01, 5.00000000e-01, 6.00000000e-01\
, 7.00000000e-01, 8.00000000e-01, 9.00000000e-01, 1.00000000e+00\
, 1.10000000e+00, 1.20000000e+00, 1.30000000e+00, 1.40000000e+00\
, 1.50000000e+00, 1.55000000e+00, 1.60000000e+00, 1.65000000e+00\
, 1.70000000e+00, 1.75000000e+00, 1.80000000e+00, 1.80000000e+00\
, 1.86530612e+00, 1.93061224e+00, 1.99591837e+00, 2.06122449e+00\
, 2.12653061e+00, 2.19183673e+00, 2.25714286e+00, 2.32244898e+00\
, 2.38775510e+00, 2.45306122e+00, 2.51836735e+00, 2.58367347e+00\
, 2.64897959e+00, 2.71428571e+00, 2.77959184e+00, 2.84489796e+00\
, 2.91020408e+00, 2.97551020e+00, 3.04081633e+00, 3.10612245e+00\
, 3.17142857e+00, 3.23673469e+00, 3.30204082e+00, 3.36734694e+00\
, 3.43265306e+00, 3.49795918e+00, 3.56326531e+00, 3.62857143e+00\
, 3.69387755e+00, 3.75918367e+00, 3.82448980e+00, 3.88979592e+00\
, 3.95510204e+00, 4.02040816e+00, 4.08571429e+00, 4.15102041e+00\
, 4.21632653e+00, 4.28163265e+00, 4.34693878e+00, 4.41224490e+00\
, 4.47755102e+00, 4.54285714e+00, 4.60816327e+00, 4.67346939e+00\
, 4.73877551e+00, 4.80408163e+00, 4.86938776e+00, 4.93469388e+00\
, 5.00000000e+00]
bh_curve = BSpline (2, [0]+list(B_KL), list(H_KL)) # [0] + is needed!

# create figure
plt.figure(1, figsize=[12, 9])
plt.clf()
plt.plot(H_KL, B_KL, '.-r')
plt.xlim(0, 2000)
plt.ylim(0, 2.12)
plt.grid()
plt.title("Magnetisation curve")
plt.xlabel("Magnetic Field Strength H  in  A/m")
plt.ylabel("Magnetic Flux Density B in T")
font = {'size'   : 16}
matplotlib.rc('font', **font)
plt.show(block=False)

The non-linear problem is solved by minimising the the functional

\[ \begin{equation} F(\textbf{A}) = \int_{\Omega}w(|\nabla\times \textbf{A}|)\;d\Omega - \int_{\Omega_c}\textbf{J}\cdot \textbf{A}\;d\Omega. \end{equation} \]

Therefore the magnetic energy density

\[ \begin{equation} w(B) = \int_0^B H(B')\;dB' \end{equation} \]

has to be set

energy_dens = bh_curve.Integrate()    # to be minimised

2.5.3. Stiffness Matrix and Dirichlet Boundaries#

The stiffness matrix is set on all regions and the regularisation term mentioned in the section description is added.

a = BilinearForm(fes, symmetric=True)

a += SymbolicBFI(1/mu0 * curl(u)*curl(v), definedon=~mesh.Materials("iron"))
a += SymbolicEnergy(energy_dens(sqrt(1e-12+curl(u)*curl(u))), definedon=mesh.Materials("iron")) # 1e-12 ... to avoid 0
a += SymbolicBFI(1e-1*u*v)  # regularisation

# preconditioner
c = Preconditioner(a, type="direct", inverse="sparsecholesky")

The boundary condition

\[ \begin{equation} \textbf{B}\cdot\textbf{n} = 0 \end{equation} \]

is prescribed with the magnetic vector potential by

\[\begin{split} \begin{align} \nabla\times\textbf{A} \cdot\textbf{n} &= 0\\ \textbf{A} &= \textbf{0}\qquad\text{on }\Gamma_{outer}. \end{align} \end{split}\]

2.5.4. Biot-Savart Field and Neumann Boundaries#

In this section, the impressed current density in the coil is modelled. This has to be done carefully, since the orientation of the current has to be defined in the curved corners consistently.

Without details:

A = 2500*1e-6       
J0 = I0N/A
# +++ bricks +++
J_brick_back = [1, 0]
J_brick_front = [-1, 0]
J_brick_left = [0, 1]
J_brick_right = [0, -1]

# +++ corners +++
# right back
x_right_back = x - 0.050
y_right_back = y - 0.050
r_right_back = (x_right_back**2 + y_right_back**2)**(1/2)
J_corner_right_back = [1/r_right_back * y_right_back, -1/r_right_back * x_right_back]

# left back
x_left_back = x + 0.050
y_left_back = y - 0.050
r_left_back = (x_left_back**2 + y_left_back**2)**(1/2)
J_corner_left_back =  [1/r_left_back * y_left_back, -1/r_left_back * x_left_back]


# left front
x_left_front = x + 0.050
y_left_front = y + 0.050
r_left_front = (x_left_front**2 + y_left_front**2)**(1/2)
J_corner_left_front = [1/r_left_front * y_left_front, -1/r_left_front * x_left_front]

# right front
x_right_front = x - 0.050
y_right_front = y + 0.050
r_right_front = (x_right_front**2 + y_right_front**2)**(1/2)
J_corner_right_front = [1/r_right_front * y_right_front, -1/r_right_front * x_right_front]

# J
val={   "corner_right_back":J_corner_right_back, "corner_left_back":J_corner_left_back, \
        "corner_left_front":J_corner_left_front, "corner_right_front":J_corner_right_front,\
        "brick_back":J_brick_back, "brick_left":J_brick_left, \
        "brick_front":J_brick_front, "brick_right":J_brick_right}    

J = J0 * CoefficientFunction([val[mat][0] if mat in val.keys() else 0 for mat in mesh.GetMaterials()])*CoefficientFunction((1,0,0)) + \
    J0 * CoefficientFunction([val[mat][1] if mat in val.keys() else 0 for mat in mesh.GetMaterials()])*CoefficientFunction((0,1,0))

Draw(J, mesh, "J")

BaseWebGuiScene

# figure
Nx = 30
Ny = 30
xi = np.linspace(-0.100, 0.100, Nx)
yi = np.linspace(-0.100, 0.100, Ny)

X, Y= np.meshgrid(xi, yi)
J_T = np.zeros([Nx, Ny, 2])

for i in range(Nx):
    for j in range(Ny):
        tmp = J(mesh(xi[i], yi[j], 0))
        J_T[i, j, 0] = tmp[0]
        J_T[i, j, 1] = tmp[1]
J_T[:, :, 0] = np.transpose(J_T[:, :, 0])
J_T[:, :, 1] = np.transpose(J_T[:, :, 1])
plt.figure(1, figsize=[10, 10])
plt.quiver(X, Y, J_T[:, :, 0], J_T[:, :, 1])
plt.xlabel("x coordinate / m")
plt.ylabel("y coordinate / m")
plt.title("Biot-Savart Field")
plt.xlim([-0.1, 0.1])
plt.ylim([-0.1, 0.1])
font = {'size'   : 16}
matplotlib.rc('font', **font)
plt.show(block=False)

The absolute value of the current density \(|\textbf{J}|\) is constant in the whole volume of the coil. The \(\textit{CoefficientFunction}\) \(\texttt{J}\) is used for the right-hand side of the weak formulation.

f = LinearForm(fes)
f += SymbolicLFI(J*v)

2.5.5. Solving the Problem#

The problem is iteratively solved by finding the minimum energy until an error falls bellow a specified limit.

t_simulation = time.time()
with TaskManager():
    f.Assemble()
    
    err = 1
    it = 1
    
    # create memories
    au = sol.vec.CreateVector()
    r = sol.vec.CreateVector()
    w = sol.vec.CreateVector()
    sol_new = sol.vec.CreateVector()

    while err > 1e-10:
        print ("nonlinear iteration", it)
        it = it+1

        # calculate current energy
        E0 = a.Energy(sol.vec) - InnerProduct(f.vec, sol.vec)
        print ("Energy old = ", E0)
        
        # solve linearized problem
        a.AssembleLinearization(sol.vec)
    
        a.Apply (sol.vec, au)
        r.data = f.vec - au

        # inverse of a
        inv = CGSolver (mat=a.mat, pre=c.mat)
        w.data = inv * r

        # calculate error
        err = InnerProduct (w, r)
        print ("err = ", err)


        sol_new.data = sol.vec + w
        # calculate new energy
        E = a.Energy(sol_new) - InnerProduct(f.vec, sol_new)
        print ("Enew = ", E)
        tau = 1
        while E > E0:
            tau = 0.5*tau
            sol_new.data = sol.vec + tau * w
            E = a.Energy(sol_new) - InnerProduct(f.vec, sol_new)
            print ("tau = ", tau, "Enew =", E)

        sol.vec.data = sol_new

        Redraw()
t_simulation = time.time() - t_simulation

print("simulation time %.3lf seconds" % t_simulation)

nonlinear iteration 1
Energy old =  2.1995519996926788e-13

err =  0.11291605303265723
Enew =  -0.06557866942113591
nonlinear iteration 2
Energy old =  -0.0655786694211361

err =  0.06287993017975554
Enew =  -0.09656272397297888
nonlinear iteration 3
Energy old =  -0.09656272397297892

err =  0.0003362050162805447
Enew =  -0.09675120222185078
nonlinear iteration 4
Energy old =  -0.09675120222185087

err =  4.407869894150019e-05
Enew =  -0.09677620250656811
nonlinear iteration 5
Energy old =  -0.09677620250656793

err =  5.506078322708533e-06
Enew =  -0.0967792369943546
nonlinear iteration 6
Energy old =  -0.0967792369943547

err =  4.0380443997320635e-07
Enew =  -0.09677945062189076
nonlinear iteration 7
Energy old =  -0.09677945062189065

err =  1.3447197699280657e-08
Enew =  -0.09677945738229421
nonlinear iteration 8
Energy old =  -0.0967794573822942

err =  1.348348115731841e-11
Enew =  -0.09677945738903328
simulation time 54.393 seconds

2.5.6. Final Step#

To show the final solution on the surface of the iron the finite element space \(H^1(\Omega)\) is defined on \(iron\) only. Further, the solution is projected in to that space.

# draw Bnorm only on iron
fesBnorm = H1(mesh, definedon=mesh.Materials("iron"))
B_norm = GridFunction(fesBnorm, "|B|")
B_norm.Set(B.Norm())
Draw(B_norm)

# some viewing options
ngsolve.internal.viewoptions.clipping.dist=-0.167

2.6. Results#

The goal of the T.E.A.M. problem 13 is to simulate the problem and to compare some characteristic values of the solution with measured values provided.