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


64
views
0
4 weeks ago by
Diego  
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
Community: FEniCS Project

1 Answer


0
4 weeks ago by
Diego  
I am sorry I was wrong.

The problem has been solved.

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

Similar posts:
Search »