Dear all,
I want to calculate the electromagnetic wave number of an infinite long conductor of 1.5 mm radius with finite conductivity in the frequency range of 0-100KHz. The mesh size around the conductor is small compared with the wavelength. Can anyone give some suggestions on how to solve such problem?
Dear all,
I calculated the propagation wave number of the electromagnetic wave problem at the frequency of 80000 Hz, but got an inaccurate solution 0.0017806972731223838-0.001756968505041909j and the accurate one should be 0.0016463654919279397-0.0017340090317261698j. However, the code can give the correct result at 2MHz frequency. I don’t know if it’s caused by the large mesh size range which span 5 order of magnitude. The code is here:
import pygmsh
import numpy as np
import gmsh, meshio
from mpi4py import MPI
from slepc4py import SLEPc
from ufl import curl, dx, FiniteElement, grad, inner, Measure, TestFunctions, TrialFunctions
from dolfinx import Constant, Function, FunctionSpace, geometry, cpp, DirichletBC
#from dolfinx.fem import locate_dofs_topological, DofMapRestriction, assemble_matrix
from dolfinx.fem import locate_dofs_topological, assemble_matrix
from scipy.constants import c, epsilon_0
from dolfinx.io import XDMFFile
from petsc4py import PETSc
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
if prune_z:
out_mesh.prune_z_0()
return out_mesh
with pygmsh.geo.Geometry() as geom:
core_r = 0.0015
max_r = 1900
inner_lcar = 0.00007
outer_lcar = 40.1
theta = np.pi / 10
p0 = [0, 0]
p1 = [core_r, 0]
p2 = [core_r * np.cos(theta), core_r * np.sin(theta)]
p3 = [max_r * np.cos(theta), max_r * np.sin(theta)]
p4 = [max_r, 0]
point_coords = [p0, p1, p2, p3, p4]
lcars = [inner_lcar, inner_lcar, inner_lcar, outer_lcar, outer_lcar]
points = []
for p, lcar in zip(point_coords, lcars):
points.append(geom.add_point([p[0], p[1]], lcar))
line1 = geom.add_line(points[0], points[1])
inner_arc = geom.add_circle_arc(points[1], points[0], points[2])
line2 = geom.add_line(points[2], points[0])
inner_loop_curves = [line1, inner_arc, line2]
line3 = geom.add_line(points[2], points[3])
outer_arc = geom.add_circle_arc(points[3], points[0], points[4])
line4 = geom.add_line(points[4], points[1])
outer_loop_curves = [inner_arc, line3, outer_arc, line4]
inner_curve_loop = geom.add_curve_loop(inner_loop_curves)
outer_curve_loop = geom.add_curve_loop(outer_loop_curves)
inner_surface = geom.add_plane_surface(inner_curve_loop)
outer_surface = geom.add_plane_surface(outer_curve_loop)
geom.add_physical(outer_arc, 'Dirichlet')
geom.add_physical(inner_surface, 'Conductor')
geom.add_physical(outer_surface, 'Air')
mesh = geom.generate_mesh(2)
gmsh.write("mesh.msh")
msh = meshio.read("mesh.msh")
# Create and save one file for the mesh, and one file for the facets
triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mt.xdmf", line_mesh)
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf_infile:
mesh = xdmf_infile.read_mesh(name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)
mvc_subdomain = xdmf_infile.read_meshtags(mesh, "Grid")
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf_infile:
mvc_boundaries = xdmf_infile.read_meshtags(mesh, "Grid")
V = FiniteElement('N1curl', mesh.ufl_cell(), 2)
Q = FiniteElement('P', mesh.ufl_cell(), 3)
W = FunctionSpace(mesh, V*Q)
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
omega = 8e4 #frequency
sigma_core = 8678
sigma_corona = 5e-6
k0_sqr = (omega/c)**2
print(k0_sqr)
S = FunctionSpace(mesh, ('DG', 0))
e_r = Function(S)
x = S.tabulate_dof_coordinates()
print(np.unique(mvc_subdomain.values))
for i in range(x.shape[0]):
if mvc_subdomain.values[i] == 3:
e_r.vector.setValueLocal(i, 1-1j*sigma_core/(omega*epsilon_0))
else:
e_r.vector.setValueLocal(i, 1)
theta_sqr = k0_sqr * 1.1
print(theta_sqr)
theta_sqr = -4e-7-7e-6j
#print(theta_sqr)
a = 1/theta_sqr*(inner(curl(u), curl(v)) - k0_sqr*e_r*inner(u, v))*dx
b = inner(grad(p), v)*dx + inner(u, v)*dx + inner(grad(p), grad(q))*dx + inner(u, grad(q))*dx - k0_sqr*e_r*inner(p, q)*dx
c = b + a
electric_wall = Function(W)
with electric_wall.vector.localForm() as bc_local:
bc_local.set(0)
print(mvc_boundaries.values)
bndry_facets = np.where(mvc_boundaries.values == 1)[0]
bdofs = locate_dofs_topological(W, 1, bndry_facets)
bc = DirichletBC(electric_wall, bdofs)
#B = assemble_matrix(b, bcs=[], restriction=(restriction_W_Gamma, restriction_W_Gamma))
B = assemble_matrix(b, [bc])
B.assemble()
#C = assemble_matrix(c, bcs=[], restriction=(restriction_W_Gamma, restriction_W_Gamma))
C = assemble_matrix(c, [bc])
C.assemble()
eps = SLEPc.EPS()
eps.create()
eps.setOperators(B, C)
eps.setDimensions(100)
eps.solve()
n = eps.getConverged()
print(n)
xr, xi = B.getVecs()
for j in range(n):
lam = eps.getEigenpair(j, xr, xi)
if (((1-1/lam)*theta_sqr)**0.5).real > 0.0002 and abs((((1-1/lam)*theta_sqr)**0.5).imag) < 1:
print(lam)
print(((1-1/lam)*theta_sqr)**0.5, 'wave number')
x_r = Function(W)
xr.copy(x_r.vector)
x_i = Function(W)
xi.copy(x_i.vector)
xt_r, xz_r = x_r.split()
xt_i, xz_i = x_i.split()
with XDMFFile(MPI.COMM_WORLD, "Et_r.xdmf", 'w',
encoding=XDMFFile.Encoding.HDF5) as file:
file.write_mesh(mesh)
file.write_function(xt_r)
with XDMFFile(MPI.COMM_WORLD, "Ez_r.xdmf", 'w',
encoding=XDMFFile.Encoding.HDF5) as file:
file.write_mesh(mesh)
file.write_function(xz_r)
with XDMFFile(MPI.COMM_WORLD, "Et_i.xdmf", 'w',
encoding=XDMFFile.Encoding.HDF5) as file:
file.write_mesh(mesh)
file.write_function(xt_i)
with XDMFFile(MPI.COMM_WORLD, "Ez_i.xdmf", 'w',
encoding=XDMFFile.Encoding.HDF5) as file:
file.write_mesh(mesh)
file.write_function(xz_i)
break
Your help is appreciated.