Postprocessing - calculation of eddy current loss

[martinbaun]

Hi,

I´ve got a question about postprocessing. I've run a magneto-harmonic 2d Simulation with linear triangle elements. Now I want to calculate the eddy current losses

$P = Re ( \frac{1}{2 \sigma} \int \int_\Delta J J^* dx dy ) $

In a book I found an equation for the integration with linear triangular elements

$P_e = \frac{\Delta}{12 \sigma} ( |J_i|^2 + |J_j|^2 + |J_k|^2 ) + Re( \frac{1}{2} ( J_i J_j^* + J_i J_k^* + J_k J_i^* +J_j J_k^* + J_j J_i^* + J_k J_j^* ) )$

with the current density $J = j \omega \sigma A$ and vector potential A given at the three nodes $i, j, k$

This works well but I wonder if there is a more simpler way to do this with skfem - even if you want to use other elements.

[kinnala]

Yes, it is possible to calculate such integrals using @Functional decorator. What is exactly your unknown? Is it $\vec{A}$ defined as ElementVector(ElementTriP1())?

[martinbaun]

That sounds good.
The vectorpotential $\vec A$ (in 2D only in z-axis $A_z$) ist my unkown and therefor the solution of the solver. Here is a small example:

import skfem
from skfem.io.json import from_dict
from skfem.helpers import dot, grad
import numpy as np
import matplotlib.pyplot as plt
plt.ion()

omega = 2 * np.pi * 0.1  # frequencies
sigma = 1e6              # electrical conductivity
mu0 = np.pi * 4e-7       # magnetic field constant

@skfem.BilinearForm(dtype=np.complex64)
def poisson_cmplx(u, v, w):
    return dot(grad(u), grad(v) / mu0) + 1j * omega * sigma * u

m = skfem.MeshTri.init_circle(4)
basis = skfem.Basis(m, skfem.ElementTriP1())

S = skfem.asm(poisson_cmplx, basis)
b = basis.point_source(np.array([0, 0]))
A = skfem.solve(*skfem.condense(S, b, D=basis.get_dofs()))

plt.figure(1)
plt.tricontourf(m.p[0,:], m.p[1,:], m.t.T, A.real  )
plt.axis('equal')
plt.colorbar()

[martinbaun]

I´ve tried the following code but i got only a single value. Seems wrong... Any hints?

J = 1j*omega*sigma*A

@skfem.LinearForm
def J2fun(u, _):
    return  1/(2*sigma)*(u*u.conj()).real

x = skfem.asm(J2fun, basis)
J2 = x @ J
ecloss = J2.real

[kinnala]

That is integral over the whole domain.

@skfem.Functional
def J2(w):
return (w.u*w.u.conj()).real

J2.elemental(J2, u=x)

[martinbaun]

Thanks, it's working :-)

For those who are interested, here is a small example:

import skfem
from skfem.io.json import from_dict
from skfem.helpers import dot, grad
import numpy as np
import matplotlib.pyplot as plt
plt.ion()

omega = 2 * np.pi * 0.01  # frequencies
sigma = 1e6              # electrical conductivity
mu0 = np.pi * 4e-7       # magnetic field constant

@skfem.BilinearForm(dtype=np.complex64)
def poisson_cmplx(u, v, w):
    return dot(grad(u), grad(v) / mu0) + 1j * omega * sigma * u

m = skfem.MeshTri.init_circle(4, smoothed=True)
basis = skfem.Basis(m, skfem.ElementTriP1())
S = skfem.asm(poisson_cmplx, basis)
b = basis.point_source(np.array([0, 0]))
A = skfem.solve(*skfem.condense(S, b, D=basis.get_dofs()))

# eddy current loss in every element
@skfem.Functional
def J2(w):
    return  1/(2*sigma)*(w.u*w.u.conj()).real

J = 1j*omega*sigma*A
ecloss = J2.elemental(basis, u=J)

# Element area (integrate over '1')
@skfem.Functional
def area_fun(w):
    return  1
area = area_fun.elemental(basis)

plt.figure(1)
plt.clf()
plt.subplot(2,1,1)
plt.title('eddy current loss')
plt.tripcolor(m.p[0,:], m.p[1,:], m.t.T, ecloss  )
plt.axis('equal')
plt.colorbar()
plt.subplot(2,1,2)
plt.title('eddy current loss density')
plt.tripcolor(m.p[0,:], m.p[1,:], m.t.T, ecloss/area  )
plt.axis('equal')
plt.colorbar()
plt.tight_layout()