### Using "project" to compute the magnetic field from the potential field

134

views

0

Hello everybody,

I have a problem with the differentiation of the potential vector Ato compute the magnetic field B.

I have made the tutorial "magnetostatics" to compute the potential vector generated by straight wires, in this tutorial the instruction:

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

is used to compute B components as the partial derivatives of A.

I have then tried to solve my own simple magnetostatics problem, namely compute the distribution of the magnetic field generated by a single wire like the (Biot Savart law).

THe potential vector A is computed but when I try to get the field B, it seems to be computed but none plot appears and even in ParaView nothing appears.

Which error have I made with "dx", "project" ?

HEre below is my own script

from __future__ import division

from fenics import *

import numpy as np

import matplotlib.pyplot as plt

from mshr import *

# Problem data

I = 100 # [A] current intensity

r = 0.01 # [m] wire radius

A = np.pi*pow(r,2) # [m2] wire cross sectional area

mo = 4*np.pi*1e-7 # [N/A2] air magnetic permeability

mu_r = 1 # [-] relative magnetic permeability

mu = mo*mu_r

# FE solver data

R = 1 # [m] mesh radius - mesh boundary

N_elem = 60 # [-] number of elements created from radius mesh

domain = Circle(Point(0,0),R)

mesh = generate_mesh(domain,100)

J = Constant(I/A) # [A/m2] current density in the wire

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

Az = TrialFunction(V)

v = TestFunction(V)

def boundary(x,on_boundary):

return on_boundary

Az_D = Constant(0)

bc = DirichletBC(V,Az_D,boundary)

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

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

L = J*v*dx

Az = Function(V)

solve(a==L,Az,bc)

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

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

plot(mesh)

interactive()

plot(Az)

interactive()

plot(B)

interactive()

vtkfile = File('Spira_B_solution.pvd')

vtkfile << B

vtkfile_Az = File('Spira_Az_solution.pvd')

vtkfile_Az << Az

Thank you in advance for your help

Diego

I have a problem with the differentiation of the potential vector Ato compute the magnetic field B.

I have made the tutorial "magnetostatics" to compute the potential vector generated by straight wires, in this tutorial the instruction:

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

is used to compute B components as the partial derivatives of A.

I have then tried to solve my own simple magnetostatics problem, namely compute the distribution of the magnetic field generated by a single wire like the (Biot Savart law).

THe potential vector A is computed but when I try to get the field B, it seems to be computed but none plot appears and even in ParaView nothing appears.

Which error have I made with "dx", "project" ?

HEre below is my own script

from __future__ import division

from fenics import *

import numpy as np

import matplotlib.pyplot as plt

from mshr import *

# Problem data

I = 100 # [A] current intensity

r = 0.01 # [m] wire radius

A = np.pi*pow(r,2) # [m2] wire cross sectional area

mo = 4*np.pi*1e-7 # [N/A2] air magnetic permeability

mu_r = 1 # [-] relative magnetic permeability

mu = mo*mu_r

# FE solver data

R = 1 # [m] mesh radius - mesh boundary

N_elem = 60 # [-] number of elements created from radius mesh

domain = Circle(Point(0,0),R)

mesh = generate_mesh(domain,100)

J = Constant(I/A) # [A/m2] current density in the wire

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

Az = TrialFunction(V)

v = TestFunction(V)

def boundary(x,on_boundary):

return on_boundary

Az_D = Constant(0)

bc = DirichletBC(V,Az_D,boundary)

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

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

L = J*v*dx

Az = Function(V)

solve(a==L,Az,bc)

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

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

plot(mesh)

interactive()

plot(Az)

interactive()

plot(B)

interactive()

vtkfile = File('Spira_B_solution.pvd')

vtkfile << B

vtkfile_Az = File('Spira_Az_solution.pvd')

vtkfile_Az << Az

Thank you in advance for your help

Diego

Community: FEniCS Project

### 1 Answer

0

Please login to add an answer/comment or follow this question.