Magnetostatics problem.

[Manuel Pineda]

Hi all, 

I'm trying to solve a magneto statics problem based on the magnetostatics subdomain problem that is on the fenics tutorial (page 99), but a different geometry.

 I should be having a magnetic field with maximum 0.8-0.9 T, instead I'm having values of 6.4 T. 

Parameters:

air permeability = 4*pi*1e-7  

iron permeability = 1000

copper permeability = 1.

J_1 = J_2= 1 A.

In the near future, I'll need to use the magnetic field to determine forces and displacement in a strong coupled problem as this image explains:

N_1, N_2 are copper windings.

Also I'm having this warning:

WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.

 And N_2 winding isn't responding but I'm solving one issue at time since I'm pretty new programming with fenics and mathematics.

Thanks for reading and any advise will be appreciated 

from __future__ import print_function

from fenics import *

from dolfin import *

from mshr import *

from math import sin, cos, pi

 

lx = 8 

ly = 8 

n = 3

r = 0.1 

sy = 5 

sx = 5 

esp = 3*r

#define geometries

domain = Rectangle(Point(0, 0), Point(lx,ly))

inf_arm = Rectangle(Point(1, 2), Point(7,1))

rig_arm = Rectangle(Point(1, 6), Point(2,2)) 

sup_arm = Rectangle(Point(1, 7), Point(7,6))

lef_arm = Rectangle(Point(6, 6), Point(7,2))

armor = sup_arm + inf_arm + rig_arm + lef_arm

a_Sn1 = [Circle(Point(0.8,sy+i/n),r) for i in range(n)]

a_En1 = [Circle(Point(2.2,sy+i/n),r) for i in range(n)]

a_Sn2 = [Circle(Point(sx+i/n,2.2),r) for i in range(n)]

a_En2 = [Circle(Point(sx+i/n,0.8),r) for i in range(n)]

domain.set_subdomain(1, armor)

for (i, wire) in enumerate(a_Sn1):

    domain.set_subdomain(2 + i, wire)

for (i, wire) in enumerate(a_En1):

    domain.set_subdomain(2 + n + i, wire)

for (j, wire) in enumerate(a_Sn2):

    domain.set_subdomain(10 + j, wire)

for (j, wire) in enumerate(a_En2):

    domain.set_subdomain(10 + n + j, wire)    

    

    

mesh = generate_mesh(domain, 64)

V = FunctionSpace(mesh, 'P', 1)

bc = DirichletBC(V, Constant(0), 'on_boundary')

markers = MeshFunction('size_t', mesh, 2, mesh.domains())

dx = Measure('dx', domain=mesh, subdomain_data=markers)

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] = 4*pi*1e-7  # perm vacio

        elif self.markers[cell.index] == 1:

            values[0] = 1000   #perm acero 

        else:

            values[0] = 1 #perm cobre  

mu = Permeability(markers, degree=1)

J_s = Constant(1)   

J_e = Constant(-1)

A_z = TrialFunction(V)

v = TestFunction(V)

a = (1 / mu)*dot(grad(A_z), grad(v))*dx

N1_s = sum(J_s*v*dx(i) for i in range(2, 2 + n))

N1_e = sum(J_e*v*dx(i) for i in range(2 + n, 2 + 2*n))

N2_s = sum(J_s*v*dx(j) for j in range(10, 2 + n))

N2_e = sum(J_e*v*dx(j) for j in range(10 + n, 2 + 2*n))

L = N1_s + N1_e + N2_s + N2_e

A_z = Function(V)

solve(a == L, A_z, bc)

W = VectorFunctionSpace(mesh, 'P', 1)

B = project(as_vector((A_z.dx(1), -A_z.dx(0))), W)

#plot(A_z)

plot(B)

vtkfile_A_z = File('coupled/potential.pvd')

vtkfile_B = File('coupled/field.pvd')

vtkfile_A_z << A_z

vtkfile_B << B

#interactive()

import matplotlib.pyplot as plt

plt.show()

Greetings,

Manuel Pineda.