Magnetomechanical static coupled problem compared with commercial software

1

Hi everyone,

I’m analyzing the magneto mechanical static 2D problem of a ferromagnetic bar influenced by an electromagnet, the image below describes the geometry:

image_2022-01-24_192049
image_2022-01-24_192049797×723 21.2 KB

All materials are assumed isotropic and the ferromagnetic material is considered as linear. Later, I’ll try to implement the nonlinearities of the material and analyze the modal response (also seeking for advise in these matters).

When comparing with a commercial software (COMSOL) I’m getting a lot of difference in the displacement results.

Thanks for reading and any help is appreciated.

here’s the code:

from __future__ import print_function
from fenics import *
from dolfin import *
import numpy as np
from mshr import *
from math import sin, cos, pi
import matplotlib.pyplot as plt
import time


w = 0.2   
h = 3.5  
eh = 0.1  

Ab = w * h  #coil area

domain = Rectangle(Point(0, 0), Point(8, 8))

#geometry

inf_arm = Rectangle(Point(1, 2), Point(6, 1))
lef_arm = Rectangle(Point(1, 6), Point(2, 2))
sup_arm = Rectangle(Point(1, 7), Point(6, 6))
arm = Rectangle(Point(6.1, 7), Point(7.1, 1))
core = sup_arm + inf_arm + lef_arm 

coil1 = Rectangle(Point(0.7, 5.7), Point(0.9, 2.2))
coil2 = Rectangle(Point(2.1, 5.7), Point(2.3, 2.2))

domain.set_subdomain(1, core)
domain.set_subdomain(2, coil1)
domain.set_subdomain(3, coil2)
domain.set_subdomain(4, arm)


mesh = generate_mesh(domain, 64)


mu_0 = 4 * pi * 1e-7 
mu_ra = 1500  
mu_rc = 1


class Permeability(UserExpression): 
    def __init__(self, markers, **kwargs):
        super(Permeability, self).__init__(**kwargs)  
        self.markers = markers

    def eval_cell(self, values, x, cell):
        if self.markers[cell.index] == 0:
            values[0] =  mu_0 #vac perm
        elif self.markers[cell.index] == 1:
            values[0] = mu_ra * mu_0   # iron perm
        elif self.markers[cell.index] == 4:
            values[0] = mu_ra * mu_0
        else:
            values[0] = mu_rc * mu_0   


mu = Permeability(markers, degree=1)
perm = project(mu, FunctionSpace(mesh, "DG", 0))

# current (A)
I_s = Constant(1)
I_e = Constant(-1)

# magnetostatic

VE = FunctionSpace(mesh, 'P', 2)
A_z = TrialFunction(VE)
v_E = TestFunction(VE)
aE = (1 / mu) * dot(grad(A_z), grad(v_E)) * dx


N1_s = I_s * v_E * dx(2)
N1_e = I_e * v_E * dx(3)

LE = (1 / Ab) * (N1_s + N1_e)  


bcE = DirichletBC(VE, Constant(0), 'on_boundary')


Az = Function(VE)
solve(aE == LE, Az, bcE)

W = VectorFunctionSpace(mesh, 'DG', 0)
B = project(as_vector((Az.dx(1), -Az.dx(0))), W)

Bx = Az.dx(1)
By = -Az.dx(0)

Hx = Bx / mu
Hy = By / mu

Fx = mu * Hx * Hy * eh
Fy = (0.5 * mu) * (Hy ** 2 - Hx ** 2) * eh
dF = project(as_vector((Fx, Fy)), W)


f = dF

E = Constant(212E6)
nu = Constant(0.3)
lmbda = (E * nu) / ((1 + nu) * (1 - 2 * nu))
miu = E / 2 / (1 + nu) 
Vi = VectorFunctionSpace(mesh, 'Lagrange', 2)

def eps(v):
    return sym(grad(v))


def sigma(v):
    return lmbda * tr(eps(v)) * Identity(2) + 2.0 * miu * eps(v)


class Solid(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (6.1, 7.1)) and between(x[1], (1.0, 7.0))


subdomain_marker = MeshFunction("size_t", mesh, mesh.geometry().dim(), 0)
Solid().mark(subdomain_marker, 1)


class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1e-14
        return between(x[0], (6.1, 7.1)) and x[1] == 1  
        

boundary_marker = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1, 0)
arm_bound = Boundary()
Boundary().mark(boundary_marker, 2)


du = TrialFunction(Vi)
u_ = TestFunction(Vi)
a = inner(sigma(du), eps(u_)) * dx(4)
fm = inner(f, u_) * dx(4)

start = time.time()

bc = DirichletBC(Vi, Constant((0, 0)), boundary_marker, 2)

u = Function(Vi, name="Displacement")
solve(a == fm, u, bc)
2

I’m comparing with comsol in order to validate the results of the code, so I’m expecting to have values approximated to the comsol results.

3

The code that you have supplied does not run, using the latest release of dolfin.
I have done several modifications to the code to make it executable. Please carefully look through the code to make sure you understand the differences

from dolfin import *
import numpy as np
from mshr import *
from math import sin, cos, pi
import matplotlib.pyplot as plt
import time


w = 0.2
h = 3.5
eh = 0.1

Ab = w * h  # coil area

domain = Rectangle(Point(0, 0), Point(8, 8))

# geometry

inf_arm = Rectangle(Point(1, 2), Point(6, 1))
lef_arm = Rectangle(Point(1, 6), Point(2, 2))
sup_arm = Rectangle(Point(1, 7), Point(6, 6))
arm = Rectangle(Point(6.1, 7), Point(7.1, 1))
core = sup_arm + inf_arm + lef_arm

coil1 = Rectangle(Point(0.7, 5.7), Point(0.9, 2.2))
coil2 = Rectangle(Point(2.1, 5.7), Point(2.3, 2.2))

domain.set_subdomain(1, core)
domain.set_subdomain(2, coil1)
domain.set_subdomain(3, coil2)
domain.set_subdomain(4, arm)


mesh = generate_mesh(domain, 64)


mu_0 = 4 * pi * 1e-7
mu_ra = 1500
mu_rc = 1


class Permeability(UserExpression):
    def __init__(self, markers, **kwargs):
        super(Permeability, self).__init__(**kwargs)
        self.markers = markers

    def eval_cell(self, values, x, cell):
        if self.markers[cell.index] == 0:
            values[0] = mu_0  # vac perm
        elif self.markers[cell.index] == 1:
            values[0] = mu_ra * mu_0   # iron perm
        elif self.markers[cell.index] == 4:
            values[0] = mu_ra * mu_0
        else:
            values[0] = mu_rc * mu_0


markers = MeshFunction("size_t", mesh, mesh.topology().dim(), mesh.domains())

mu = Permeability(markers, degree=1)
perm = project(mu, FunctionSpace(mesh, "DG", 0))

# current (A)
I_s = Constant(1)
I_e = Constant(-1)

# magnetostatic
dxC = Measure("dx", domain=mesh, subdomain_data=markers)
VE = FunctionSpace(mesh, 'P', 2)
A_z = TrialFunction(VE)
v_E = TestFunction(VE)
aE = (1 / mu) * dot(grad(A_z), grad(v_E)) * dxC


N1_s = I_s * v_E * dxC(2)
N1_e = I_e * v_E * dxC(3)


LE = (1 / Ab) * (N1_s + N1_e)


bcE = DirichletBC(VE, Constant(0), 'on_boundary')


Az = Function(VE)
solve(aE == LE, Az, bcE)

W = VectorFunctionSpace(mesh, 'DG', 0)
B = project(as_vector((Az.dx(1), -Az.dx(0))), W)

Bx = Az.dx(1)
By = -Az.dx(0)

Hx = Bx / mu
Hy = By / mu

Fx = mu * Hx * Hy * eh
Fy = (0.5 * mu) * (Hy ** 2 - Hx ** 2) * eh
dF = project(as_vector((Fx, Fy)), W)


f = dF

E = Constant(212E6)
nu = Constant(0.3)
lmbda = (E * nu) / ((1 + nu) * (1 - 2 * nu))
miu = E / 2 / (1 + nu)
Vi = VectorFunctionSpace(mesh, 'Lagrange', 2)


def eps(v):
    return sym(grad(v))


def sigma(v):
    return lmbda * tr(eps(v)) * Identity(2) + 2.0 * miu * eps(v)


class Solid(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (6.1, 7.1)) and between(x[1], (1.0, 7.0))


subdomain_marker = MeshFunction("size_t", mesh, mesh.geometry().dim(), 0)
Solid().mark(subdomain_marker, 1)


class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1e-14
        return between(x[0], (6.1, 7.1)) and x[1] == 1


boundary_marker = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1, 0)
arm_bound = Boundary()
Boundary().mark(boundary_marker, 2)


du = TrialFunction(Vi)
u_ = TestFunction(Vi)
a = inner(sigma(du), eps(u_)) * dxC(4)
fm = inner(f, u_) * dxC(4)
start = time.time()

bc = DirichletBC(Vi, Constant((0, 0)), boundary_marker, 2)
A = assemble(a, keep_diagonal=True)
b = assemble(fm)
A.ident_zeros()
bc.apply(A, b)
u = Function(Vi, name="Displacement")
solve(A, u.vector(), b)
XDMFFile("u.xdmf").write_checkpoint(u, "u", 0)
4

Hi,
Thank you!, yes I realize I omitted the markers definition. And you added the lines where the matrix is defined.
cheers