Electro‑thermal simulation of a bimetal – problem with imposed current on top surface (mesh from Gmsh, units in mm)

General
1

Dear FEniCS community,

I am developing an electro‑thermal solver for a bimetal strip in a circuit breaker. I generate the mesh with Gmsh (dimensions in mm ), then convert it to XDMF using meshio . The geometry is a 3D volume with two distinct contact areas on the top face (see attached sketch):

sketch
sketch1197×616 41.4 KB

  • Entrance (current injection): a small rectangle (26.5 × 4.5 mm²)
  • Outlet (voltage reference): a circle of radius 2.75 mm

I want to impose a total current (e.g. 10 A) on the rectangular face, while the circular face is set to zero electric potential (ground). The rest of the boundary is electrically insulated.

What works: The mesh loads correctly, and the thermal part seems to run, but the electrical part does not conserve current. The computed current through the entrance is always much lower than the imposed value, and the sum of entrance + outlet currents is not zero. Consequently the Joule power is far too low and the temperature rise is negligible (or sometimes unrealistically high if I mis‑tune properties).

What I have tried:

  1. Direct marking of physical groups in Gmsh – I assigned physical surfaces to the rectangle (tag 42) and the circle (tag 43). When I import the boundary mesh into FEniCS, the number of facets is very low (e.g. 86 for the rectangle, 262 for the circle) and the computed area does not match the geometric area (the rectangle should be ~9.9 mm², but I get ~8.1 mm²).
    Using these tags, the current is not conserved.
  2. Defining subdomains directly in FEniCS – I tried to bypass the Gmsh markers by writing SubDomain classes with geometric tolerances (e.g. checking x[2] for the top face and the rectangle/circle inequalities). The facet count improves but the current is still far from the imposed value.
  3. Unit consistency – My mesh coordinates are in mm . I have tried two approaches:
  • Keep everything in mm: σ in S/mm, k in W/(mm·K), h in W/(mm²·K), and compute areas in mm².
  • Scale the mesh to meters (mesh.coordinates()[:] *= 1e-3 ) and use SI units.
    Both give similar discrepancies, so I suspect the issue is not purely unit‑related.
from fenics import *
import sys

TAG_SURF_ENTREE = 42  
TAG_SURF_SORTIE = 43   

print("1. Chargement du maillage...")
mesh = Mesh()
with XDMFFile("bilame_volume.xdmf") as infile:
    infile.read(mesh)


boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
boundaries.set_all(0)


class SurfaceEntree(SubDomain):
    def inside(self, x, on_boundary):
        # Rectangle : x entre 26.5 et 28.7, y entre 1.5 et 6.0, z = 0.6
        tol = 1e-3
        if not (on_boundary and near(x[2], 0.6, tol)):
            return False
        return (26.5 - tol <= x[0] <= 28.7 + tol) and (1.5 - tol <= x[1] <= 6.0 + tol)

class SurfaceSortie(SubDomain):
    def inside(self, x, on_boundary):
        # Cercle : centre (2.75, 3.075), rayon 2.75, z = 0.6
        tol = 1e-3
        if not (on_boundary and near(x[2], 0.6, tol)):
            return False
        r2 = (x[0] - 2.75)**2 + (x[1] - 3.075)**2
        return r2 <= (2.75 + tol)**2

entree = SurfaceEntree()
entree.mark(boundaries, TAG_SURF_ENTREE)
sortie = SurfaceSortie()
sortie.mark(boundaries, TAG_SURF_SORTIE)
print("Nb facettes entrée :", sum(1 for _ in SubsetIterator(boundaries, TAG_SURF_ENTREE)))
print("Nb facettes sortie :", sum(1 for _ in SubsetIterator(boundaries, TAG_SURF_SORTIE)))
# Vérification rapide
print("Facettes tag 42 (entrée rect.) :", sum(1 for _ in SubsetIterator(boundaries, TAG_SURF_ENTREE)))
print("Facettes tag 43 (sortie cercle) :", sum(1 for _ in SubsetIterator(boundaries, TAG_SURF_SORTIE)))

ds = Measure('ds', domain=mesh, subdomain_data=boundaries)


V_space = FunctionSpace(mesh, 'P', 1)
T_space = FunctionSpace(mesh, 'P', 1)


sigma = Constant(5.86e4)     # S/mm (5.86e7 S/m) 
k = Constant(0.1121)        # W/mm.K (25 W/m.K)
T_ambiant = Constant(293.15)  # K


print("2. Résolution électrique (courant entrant sur le rectangle, masse sur le cercle)...")
I_value = 10.0  # A

u = TrialFunction(V_space)
v = TestFunction(V_space)

a_elec = sigma * dot(grad(u), grad(v)) * dx
L_elec = Constant(0.0) * v * dx


bc_elec_ground = DirichletBC(V_space, Constant(0.0), boundaries, TAG_SURF_SORTIE)
bcs_elec = [bc_elec_ground]


Area_entree = assemble(Constant(1.0) * ds(TAG_SURF_ENTREE))
if Area_entree <= DOLFIN_EPS:
    raise RuntimeError("Surface d'entrée vide ou mal marquée.")
print(f"Aire de l'entrée (rectangle) : {Area_entree:.2f} mm² (théorique 9.9 mm²)")


jn = Constant(I_value / Area_entree)



L_elec -= jn * v * ds(TAG_SURF_ENTREE)

Potentiel = Function(V_space, name="Potentiel_Electrique")
solve(a_elec == L_elec, Potentiel, bcs_elec)


n = FacetNormal(mesh)
J = - sigma * grad(Potentiel)
I_calc_entree = assemble(dot(J, n) * ds(TAG_SURF_ENTREE))
I_calc_sortie = assemble(dot(J, n) * ds(TAG_SURF_SORTIE))
Area_sortie = assemble(Constant(1.0) * ds(TAG_SURF_SORTIE))
V_entree_moy = assemble(Potentiel * ds(TAG_SURF_ENTREE)) / Area_entree
R_calc = V_entree_moy / I_value
Q_Joule = sigma * dot(grad(Potentiel), grad(Potentiel))
P_Joule = assemble(Q_Joule * dx)

Questions:

  1. Is there a robust way to impose a total current on a surface in FEniCS when the mesh is in mm ? Do I need to manually scale the facet normals or the integration measure?

  2. Does the sign of the Neumann term change when the mesh is not unit‑scaled? I have read that ds includes the correct area element, so the expression jn * v * ds should be dimensionally correct regardless of units.

  3. Are there any known pitfalls with meshio conversion of 2D physical groups from Gmsh? Perhaps the facet‑to‑cell mapping is lost, leading to wrong area computation.

Thank you very much for any insight!