Simple Sheets
Download this Jupyter File
Constrained Eddy Current Problem in 2D with the Magnetic Vector Potential¶
\begin{align*} \int_\Omega \nu \nabla u \cdot \nabla u' + j\omega\sigma u u' \;d\Omega = 0 \end{align*} with $u = B_0 x$ or $u = B_0 y$ on the boundary $\Gamma_D = \Gamma_{left}\cup\Gamma_{right}$ or $\Gamma_D = \Gamma_{top}\cup\Gamma_{bottom}$
and the constraint for each sheet: \begin{align*} \int_{\Omega_{c, i} } J \;d\Omega = 0 \end{align*}
results in an extension of the base with lagrange multipliers
\begin{align*} \int_\Omega \nu \nabla u \cdot \nabla u' + j\omega\sigma (u+\lambda_i)(u'+\lambda_i') \;d\Omega = 0 \end{align*}
In [1]:
from netgen.occ import *
import numpy as np
from ngsolve import *
from ngsolve.webgui import Draw
maxh = 1/10
freq = 50
mu0 = 4e-7*np.pi
mu_Fe = mu0 *1000
mu_Air = mu0
sigma_Fe =2e6
sigma_Air = 1e-3
omega = freq*2*np.pi
delta = np.sqrt(2/(sigma_Fe*omega*mu_Fe))
In [2]:
Nsheets = 6
ff = 0.5
d = delta/2
order0 = 2
B0 = 1
use_symmetry = True
use_insulation=False
In [3]:
dFe = d*ff
d0 = d*(1-ff)
H_core = Nsheets*dFe + (Nsheets-1)*d0
W_core = H_core
W = 10*W_core
H = 10*H_core
wp = WorkPlane()
outer = wp.RectangleC(W, H).Face()
outer.name = "air"
outer.edges.Max(X).name = "right"
outer.edges.Min(X).name = "left"
outer.edges.Max(Y).name = "top"
outer.edges.Min(Y).name = "bottom"
if use_insulation:
insulation = wp.RectangleC(W_core, H_core).Face()
insulation.name = "insulation"
insulation.edges.Max(X).name = "insulation_right"
insulation.edges.Min(X).name = "insulation_left"
insulation.edges.Max(Y).name = "insulation_top"
insulation.edges.Min(Y).name = "insulation_bottom"
x_pos = - W_core/2
rec_sheets =[]
for i in range(Nsheets):
wp.MoveTo(x_pos, -H_core/2)
rec_sheets.append(wp.Rectangle(dFe, H_core).Face())
rec_sheets[-1].name = f"iron{i}"
rec_sheets[-1].edges.name = "interface"
rec_sheets[-1].edges.maxh = delta/10
x_pos += d
if use_symmetry:
wp.MoveTo(0, 0)
cutting = wp.Rectangle(W, H).Face()
cutting.edges.Min(X).name = "symmetry_left"
cutting.edges.Min(Y).name = "symmetry_bottom"
all_sheets = Glue(rec_sheets)
if use_insulation:
geo = Glue([cutting * (outer - insulation), cutting * (insulation - all_sheets), cutting * all_sheets])
else:
geo = Glue([cutting * (outer - all_sheets), cutting * all_sheets])
else:
geo = Glue([outer - Glue(rec_sheets), Glue(rec_sheets)])
meshRef = Mesh(OCCGeometry(geo, dim=2).GenerateMesh(maxh=delta))
settings = {"Objects":{"Wireframe":False}}
Draw(meshRef.MaterialCF({"iron.*":1, "air":2}, default=0), meshRef)
#Draw(x, meshRef)
print("domains", set(meshRef.GetMaterials()))
print("bnds", set(meshRef.GetBoundaries()))
print("penetration depth", delta)
domains {'iron5', 'iron3', 'air', 'iron4'}
bnds {'top', 'interface', 'symmetry_bottom', 'right', 'symmetry_left'}
penetration depth 0.0015915494309189533
In [4]:
mu = meshRef.MaterialCF({"iron.*":mu_Fe, "air":mu_Air, "insulation":mu_Air}, default=0)
sigma = meshRef.MaterialCF({"iron.*":sigma_Fe}, default=0)
nu = 1/mu
rho = 1/sigma
from myPackage import drawBndAll
#drawBndAll(meshRef, drawFunc=Draw, block=False)
In [5]:
excitation_orientation = "y"
# ------------------------------------------------------------------------------
# --- Excitation
# ------------------------------------------------------------------------------
def solveWithA(mesh, order0, omega):
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# +++ reference solution
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
print("order0", order0)
print("numSheets", Nsheets)
dir_A_sym = ""
if excitation_orientation == "x":
dir_A = "top|bottom"
if use_symmetry:
dir_A_sym = "|symmetry_bottom"
else:
dir_A = "left|right"
if use_symmetry:
dir_A_sym = "|symmetry_left"
VA = H1(mesh,order=order0, complex=True, dirichlet=dir_A + dir_A_sym)
VNum = []
for i in range(0 if not use_symmetry else int(Nsheets/2), Nsheets):
VNum.append(NumberSpace(mesh, complex=True, definedon=mesh.Materials(f"iron{i}")))
V = FESpace([VA] + VNum)
ndof = V.ndof
print(f"VA :{VA.ndof}")
print(f"ndof :{ndof}")
# Gridfunctions
sol_ref = GridFunction(V, "sol")
A_ref = sol_ref.components[0]
trials = V.TrialFunction()
tests = V.TestFunction()
uA = trials[0]
vA = tests[0]
# ------------------------------------------------------------------------------
# Matrix
# ------------------------------------------------------------------------------
with TaskManager():
# Bilinear form with
ah_ref = BilinearForm(V, symmetric=True)
# A:
ah_ref += nu*grad(uA) * grad(vA) * dx
# lagrange multipliers
for i in range(0 , Nsheets if not use_symmetry else int(Nsheets/2)):
domain_name = f"iron{i + (0 if not use_symmetry else int(Nsheets/2))}"
ah_ref += 1j*omega * sigma_Fe * (uA + trials[1+i]) * (vA + tests[1+i]) * dx(domain_name)
prec = Preconditioner(ah_ref, type="direct")
ah_ref.Assemble()
print(ah_ref.mat.AsVector().Norm())
f_ref = LinearForm(V)
f_ref.Assemble()
# ------------------------------------------------------------------------------
# ------ Solve It
# ------------------------------------------------------------------------------
with TaskManager():
if excitation_orientation == "x":
A_ref.Set(B0*y, BND, definedon=mesh.Boundaries(dir_A))
else:
A_ref.Set(B0*x, BND, definedon=mesh.Boundaries(dir_A))
# solve it
solvers.BVP(bf=ah_ref, lf=f_ref, gf=sol_ref, pre=prec, maxsteps=5)
gradA = grad(A_ref)
B = CF((gradA[1], -gradA[0]))
H = 1/mu * B
E = -1j*omega*(sum(sol_ref.components))
J = sigma * E
Q = Integrate(1/2 * InnerProduct(H, B), mesh).real
P = Integrate(1/2 * InnerProduct(E, J), mesh).real
if use_symmetry:
Q, P = Q*4, P*4
return A_ref, B, E, H, J, P, Q, sol_ref
In [6]:
A, B_A, E_A, H_A, J_A, P_A, Q_A, sol_A = solveWithA(meshRef, order0, omega)
print(f"eddy current losses\t{P_A:10.3f}")
print(f"reactive energy\t\t{Q_A:10.3f}")
order0 2 numSheets 6 VA :2039 ndof :2042 59228946.21955126 CG iteration 1, residual = 76.81714315775568 CG iteration 2, residual = 3.1962311498725343e-13 eddy current losses 0.251 reactive energy 746.039
In [7]:
settings.update({"camera":{"transformations":[{"type":"move", "dir":(0.3, 0.3, 1.000), "dist":1}]}})
scene = Draw(J_A.imag, meshRef, settings=settings)
Draw(B_A.real, meshRef, vectors=True, settings=settings, min = 0, max=4, autoscale=False )
Out[7]:
BaseWebGuiScene
T Formulation¶
In [8]:
def solveWithTPhi(mesh, order0, omega, nograds=False):
# ------------------------------------------------------------------------------
# --- materials
# ------------------------------------------------------------------------------
mu = mesh.MaterialCF({"air": mu_Air, "iron.*": mu_Fe}, default=0.001)
nu = 1/mu
sigma = mesh.MaterialCF({"air": sigma_Air, "iron.*": sigma_Fe}, default=0.001)
rho = 1/sigma
# ------------------------------------------------------------------------------
# --- fem
# ------------------------------------------------------------------------------
dir_T = "interface"
dir_Phi_sym = ""
if excitation_orientation == "x":
dir_Phi = "left|right"
if use_symmetry:
dir_Phi_sym = "|symmetry_left"
dir_T += "|symmetry.*"
else:
dir_Phi = "top|bottom"
if use_symmetry:
dir_Phi_sym = "|symmetry_bottom"
dir_T += "|symmetry.*"
fesT = HCurl(mesh, order=order0, dirichlet=dir_T, complex=True, definedon=mesh.Materials("iron.*"), nograds=nograds)
fesPhi = H1(mesh, order=order0+1, dirichlet=dir_Phi+dir_Phi_sym, complex=True)
fes = FESpace([fesPhi, fesT])
trials = fes.TrialFunction()
tests = fes.TestFunction()
uPhi, vPhi = trials[0], tests[0]
uT, vT = trials[1], tests[1]
sol = GridFunction(fes, "sol")
Phi = sol.components[0]
T = sol.components[1]
# ------------------------------------------------------------------------------
# --- fields
# ------------------------------------------------------------------------------
H = T - grad(Phi)
B = mu * H
J = curl(T)
E = rho * J
# ------------------------------------------------------------------------------
# --- formulation
# ------------------------------------------------------------------------------
ah = BilinearForm(fes, symmetric=True)
ah += rho * curl(uT) * curl(vT) * dx("iron.*")
ah += 1j * omega * mu * (uT - grad(uPhi)) * (vT - grad(vPhi)) * dx
ah += 1e-1 * uPhi * vPhi * dx("iron.*")
prec = Preconditioner(ah, type="direct")
f = LinearForm(fes)
ah.Assemble()
f.Assemble()
if excitation_orientation == "x":
Phi.Set( x * B0 * 1/mu0, BND, definedon=mesh.Boundaries(dir_Phi))
else:
Phi.Set( y * B0 * 1/mu0, BND, definedon=mesh.Boundaries(dir_Phi))
solvers.BVP(bf = ah, lf= f, pre=prec, gf=sol, maxsteps=10)
# sol.vec.data = ah.mat.Inverse(freedofs=fes.FreeDofs()) * f.vec
Q = Integrate(1/2 * InnerProduct(H, B), mesh).real
P = Integrate(1/2 * InnerProduct(E, J), mesh).real
if use_symmetry:
Q, P = 4*Q, 4*P
return T, Phi, B, E, H, J, P, Q, sol
T, Phi, B_T, E_T, H_T, J_T, P_T, Q_T, sol_T = solveWithTPhi(meshRef, order0, omega)
CG iteration 1, residual = 1308.3639121271974 CG iteration 2, residual = 1.3778666050097097e-07
In [9]:
Draw(J_T.imag, meshRef, settings=settings)
Draw(B_T.real, meshRef, vectors=True, settings=settings, min = 0, max=4, autoscale=False )
Out[9]:
BaseWebGuiScene
In [10]:
print(f"A-formulation eddy current losses\t{P_A:10.3f}")
print(f"reactive energy\t\t{Q_A:10.3f}")
print(f"ndof\t\t{sol_A.space.ndof}")
print(f"ndof A\t\t{A.space.ndof}")
print("")
print(f"TPhi-formulation eddy current losses\t{P_T:10.3f}")
print(f"reactive energy\t\t{Q_T:10.3f}")
print(f"ndof\t\t{sol_T.space.ndof}")
print(f"T ndof\t\t{T.space.ndof}")
print(f"Phi ndof\t\t{Phi.space.ndof}")
A-formulation eddy current losses 0.251 reactive energy 746.039 ndof 2042 ndof A 2039 TPhi-formulation eddy current losses 0.263 reactive energy 778.748 ndof 7194 T ndof 2669 Phi ndof 4525
In [ ]:
In [ ]: