Imposing a Global Current Constraint in a Maxwell-Ampère Problem (FEniCSx)

variational formulation
1

Hello everyone,

I am solving a time-dependent Maxwell-Ampère equation in a 3D cylindrical geometry using Nédélec elements of the first kind in FEniCSx. My computational domain consists of two concentric cylinders:

  • Inner cylinder: Represents the conductor.
  • Outer cylinder: Represents the air domain.

The PDE I am solving is:

\nabla \times (\rho \nabla \times \mathbf{H}) = -\mu_0 \frac{\partial \mathbf{H}}{\partial t}​

where:

  • \mathbf{H} is the magnetic field,
  • \rho is the resistivity (different between the two materials),
  • \mu_0​ is the permeability.

Boundary Conditions (Already Implemented and Working)

  • Dirichlet BC on the lateral surface of the outer cylinder → Imposed separately.
  • Homogeneous Neumann BC on the top and bottom surfaces → Naturally arises from the weak form due to symmetry.

Outstanding Issue: Imposing a Global Current Constraint

I need to enforce a global constraint on the total flux of the current density \mathbf{J} = \nabla \times \mathbf{H} over the top and bottom surfaces of the conductor (facet tags ft = 5 and 6), ensuring that:

\int_{S} \mathbf{J} \cdot \mathbf{n} \, dS = \int_{S} \nabla\times\mathbf{H} \cdot \mathbf{n} \, dS= I.

This condition should prescribe a fixed total current I over these surfaces while still allowing for the natural distribution of \mathbf{H}.

From a theoretical perspective, this is often handled using Lagrange multipliers or similar approaches. I am looking for advice on the best way to implement this constraint in FEniCSx, considering that I have the facet tags identifying the relevant surfaces.

Current Implementation

Below is my current code, where I have already implemented the geometry, function spaces, weak formulation, and boundary conditions. However, the current constraint is still missing. Any suggestions on how to properly enforce this condition in FEniCSx would be greatly appreciated!

import gmsh
from dolfinx.io import gmshio
from mpi4py import MPI
import pyvista as pv
from dolfinx.fem import functionspace, Function, locate_dofs_topological, dirichletbc, Expression
import numpy as np
from ufl import as_vector, rot, dot, TestFunction, TrialFunction, Measure, sqrt
from dolfinx.fem.petsc import LinearProblem
from dolfinx.plot import vtk_mesh


######## MESH GENERATION #########
length = 10
R_cyl = 1
R_ext = 5

gmsh.initialize()

cyl = gmsh.model.occ.addCylinder(0,0,-length/2,0,0,length,R_cyl)
ext_cyl = gmsh.model.occ.addCylinder(0,0,-length/2,0,0,length,R_ext)
gmsh.model.occ.synchronize()

model_dim_tags = gmsh.model.occ.fragment([(3, ext_cyl)], [(3, cyl)])
gmsh.model.occ.synchronize()

# Add physical tag for the volume (internal cylinder->tag=0, sphexternal cylinderere->tag=1)
volume_entities = [model_dim_tags[0][0][1]]
gmsh.model.addPhysicalGroup(3, volume_entities, tag=0) # internal cylinder
volume_entities = [model_dim_tags[0][1][1]]
gmsh.model.addPhysicalGroup(3, volume_entities, tag=1) # external cylinder
gmsh.model.occ.synchronize()

# assign facets tags
ext_boundary = gmsh.model.getBoundary([model_dim_tags[0][1]])
int_boundary = gmsh.model.getBoundary([model_dim_tags[0][0]])
gmsh.model.occ.synchronize()

# create physical groups of surfaces
for tag, boundary in enumerate(ext_boundary + int_boundary):
    gmsh.model.addPhysicalGroup(2, [boundary[1]], tag=tag)
gmsh.model.occ.synchronize()

# Generate mesh
gmsh.model.mesh.generate(3)

_ = gmsh.write("issue.msh")

gmsh.finalize()

domain, ct, ft = gmshio.read_from_msh("issue.msh", MPI.COMM_WORLD, 0, gdim=3)


######## FUNCTION SPACES AND FUNCTIONS #########
N = functionspace(domain, ('N1curl',1)) # Nedelec of the first kind
DG0 = functionspace(domain, ('DG',0))   # Discrete Galerkin, degree 0
DG1 = functionspace(domain, ('DG',1))   # Discrete Galerkin, degree 0

H = Function(N)  # Magnetic field on Nadelec elements


######## MATERIAL PROPERTIES #########
rho = Function(DG0)
rho.x.array[:] = 100 # air (external cylinder)
rho.x.array[ct.find(0)] = 1e-6 # Conductor (internal Cylinder)

mu0 = 4*np.pi*1e-7


######## DIRICHELET BOUNDARY CONDITIONS #########
domain.topology.create_connectivity(1,2)
f2n = domain.topology.connectivity(1,2)
f2n_mat = np.vstack([f2n.array[range(i,f2n.array.shape[0],3)] for i in range(3)]).T
boundary_edges = list(set(f2n_mat[ft.find(1),:].flatten()))
external_boundary = locate_dofs_topological(N, entity_dim=2, entities=boundary_edges)

Hext = Function(N)
bc = dirichletbc(Hext, external_boundary)


######## INITIAL CONDITION AND TIME DISCRETIZATION #########
dt = 2e-4
t_end = 100*dt
n_samples = int(t_end/dt)
u_n = Function(N)
u_n.x.array[:] = 0


######## WEAK FORMULATION #########
dx = Measure("dx", domain=domain, subdomain_data=ct) # for volume integrals
ds = Measure("ds", domain=domain, subdomain_data=ft) # for surface integrals

u = TrialFunction(N)
v = TestFunction(N)

lhs =  dt*dot(rot(v), rho*rot(u))*dx() + dot(v, mu0*u)*dx
rhs =  dot(v, mu0*u_n)*dx
res = lhs - rhs


######## TIME DOMAIN SIMULATION #########
H_vec = np.zeros([H.x.array.shape[0], n_samples])
for t in range(n_samples):
    tt = t*dt
    print(f't: {tt}s')

    # Time varying bcs
    Hext_expr = Expression(H-H+as_vector([0,0,0.1/mu0*np.sin(2*np.pi*50*tt)]), N.element.interpolation_points())
    Hext.interpolate(Hext_expr)

    # solve the linear problem at instant t
    problem = LinearProblem(a=lhs, L=rhs, u=H, bcs=[bc])
    problem.solve()

    # save the solution at instant t
    H_vec[:,t] = H.x.array[:]
    # update previous iteration field
    u_n.x.array[:] = H.x.array[:]

Thanks in advance for your help!

2

I would probably use a real function space, as for instance done in

2 Likes