 |
BELFEM 0.9.0
Berkeley Lab Finite Element Framework
|
|
BELFEM 0.9.0
Berkeley Lab Finite Element Framework
|
Date: 2026-04-09 Purpose: Derive the surface bilinear form for a thin resistive contact layer (rint) collapsed to a zero-thickness interface in the H-formulation, and clarify what this buys us relative to the existing PENTA6TS thin-shell treatment. Module: fem/maxwell Status: Theory note. No implementation yet. Related: Schnaubelt et al. 2023 §III
The rint layer sits between the silver coating and the YBCO film in a REBCO tape stack. Its role is to model the electrical contact resistance between the two conductors.
| Property | Value |
|---|---|
| Resistivity rho_rint | ~1e-3 Ohm*m |
| Thickness h_rint | ~0.1 um = 1e-7 m |
| Contact resistance per unit area R_ct = rho * h | ~1e-10 Ohm*m^2 |
| Through-thickness stiffness rho / h | ~1e4 Ohm |
The dominant physics is through-thickness current transfer: current flows from Ag through rint into YBCO (normal to the layer plane). In-plane currents within the rint are negligible because the layer is too thin and too resistive for significant tangential current flow.
Starting from the H-formulation weak form in the rint volume Omega_rint (thickness h, resistivity rho):
The first term (mass) is O(mu_0 * h) ~ 1e-13 and negligible for thin layers. The stiffness term dominates.
Following Schnaubelt 2023 §III (Equations 6-8), the thin-shell approximation (TSA) collapses the rint volume to a surface Gamma (the rint midplane). The H field is discretized through the virtual thickness using first-order Lagrange basis functions:
where w in [0, h], Xi_0(w) = (h - w)/h, Xi_1(w) = w/h, and H^-, H^+ are the tangential H traces on the bottom (Ag) and top (YBCO) sides of the interface.
The curl through the thickness gives the tangential current:
These are O(1/h) and dominate the curl. Integrating the stiffness through the thickness with the 1D Lagrange FE matrices:
The resulting surface integral is:
where [[H_t]] = H^- - H^+ is the tangential H jump across the interface.
This is the contact-impedance term. It scales as rho / h.
The tangential curl (in-plane derivatives) gives the normal current component. Integrating through the thickness:
The resulting surface integral scales as rho * h ~ 1e-10, which is 14 orders of magnitude smaller than the through-thickness term and is negligible for the contact resistance case.
The mass contribution scales as mu_0 * h ~ 1e-13 and is negligible.
For the rint application, retaining only the dominant through-thickness term:
This is the simplified surface bilinear form for the collapsed rint interface. The no-consistency-term assumption rests on treating this as a natural (Robin-type) boundary condition derived from the energy functional. Juntunen & Stenberg 2009 ("Nitsche's method for general boundary conditions") supports this for scalar elliptic problems, and the physical derivation (§3.1) suggests it extends to the H(curl) case. A formal proof for BELFEM's specific H(curl) conductor-conductor interface is still missing. Either a proof or a numerical convergence study confirming optimal order is an implementation gate. Until one of the two exists, treat the absence of consistency terms as a plausible inference, not a settled result.
Two quantities with reciprocal thickness scaling must not be confused:
Both the bulk discretization and the collapsed interface produce the same matrix scaling:
The scaling is the same, so the TSA does not change the conditioning ratio. For rint with rho = 1e-3 and h = 1e-7, the matrix entry scales as 1e4 in both forms, and the conditioning ratio against adjacent Ag (rho/h ~ 2e-3) remains ~10^7.
The surface form does not improve conditioning. Its advantages are:
| Benefit | Bulk rint layer | TSA surface interface |
|---|---|---|
| DOFs in rint | 6 edge_h per element (3 bottom + 3 top) | 0 (collapsed to interface) |
| Coupling interfaces | 2 (rint-Ag + rint-YBCO) | 1 (Ag-YBCO directly) |
| Nitsche penalty parameter | Required (alpha must bridge rint-Ag contrast) | Not needed (the impedance is the coupling) |
| Matrix extreme values | Entire rint block has K ~ 1e4 | Only interface DOFs have ~ 1e4 contribution |
| Conditioning ratio | ~10^7 (rint block vs Ag block) | ~10^7 (interface DOFs vs Ag interior DOFs) |
| Localization | Extreme values distributed across all rint DOFs | Extreme values localized to interface DOFs only |
The key practical advantage is the third table row: eliminating the separate penalty parameter. The current h_ghost() kernel in mt_maxwell_h.cpp already bounds alpha via a regularized harmonic mean with k_reg, so the interface does not produce runaway off-diagonal entries. The regularized harmonic mean is a stabilization choice that trades accuracy for conditioning. The TSA contact-impedance form eliminates the issue at the source: the physical impedance rho/h is the coupling coefficient, exactly right by construction, with no tunable parameter and no accuracy-conditioning tradeoff.
The existing PENTA6TS thin-shell model (the h_picard / h_newton_* kernels driving the thin-shell branch of calculator::MaxwellData) is the N=1 TSA. It assembles the full surface bilinear form: mass (§3.3), in-plane stiffness (§3.2), and through-thickness stiffness (§3.1). For the rint application, only the through-thickness term (§3.1) is numerically significant.
The difference between the existing thin-shell model and the collapsed interface:
h_ghost() in mt_maxwell_h.cpp weakly enforces tangential H continuity:
The penalty form is: alpha * [[H_t]] . [[v_t]] plus SIPG consistency/adjoint terms.
The contact-impedance interface enforces a different condition:
The coefficient is the physical contact impedance rho / h, not a tunable stabilization parameter. Consistency and adjoint terms remain an open question (see §3.4 above); current evidence suggests they are not needed, but this has not been proven for the H(curl) case.
The matrix forms look structurally similar because both are coefficient * [[H_t]] . [[v_t]], with these differences:
The contact-impedance kernel operates on an interface sideset between two conductor blocks (Ag and YBCO). "No DOFs in rint" means no rint-owned block DOFs because the rint block does not exist. The H+ and H- traces at the interface still carry edge_h DOFs from the Ag and YBCO conductor blocks. These DOFs must remain distinct, rather than merged, so the jump [[H_t]] = H_m - H_s is representable. This is the same distinct-trace requirement as in Schnaubelt 2023 §III and Alves 2024.
At each integration point on the interface, the kernel assembles:
This is assembled into the stiffness matrix K() (not M()), because it represents resistive dissipation in the contact layer. The mass-matrix contribution is negligible for thin rint and is omitted.
This is structurally identical to the penalty part of h_ghost(), except that it omits the consistency and adjoint terms (the Dm/Ds curl-operator terms in h_ghost()). A new kernel function h_contact_impedance() in mt_maxwell_h.cpp would be cleaner than reusing h_ghost() with modified alpha.
The buffer layer (phi formulation) and the rint contact impedance solve different problems:
| Aspect | Buffer (phi formulation) | Rint (contact impedance) |
|---|---|---|
| Physics | True insulator, J = 0 everywhere | Resistive contact, J_through ≠ 0 |
| DOFs | Nodal phi (scalar potential) | None (collapsed to interface) |
| Bulk operator | mu_0 * grad(phi) . grad(v) (Laplace) | None |
| Interface coupling | Static condensation (edge-to-node) | (rho/h) * [[H_t]] . [[v_t]] (impedance) |
| Current transfer | Zero | Finite, determined by rho/h |
| Cut topology | Rerouted through buffer node (see buffer_cut_topology.md) | Not affected for the local same-loop case (rint sits between two conductors in the same current loop). In NI-coil geometries where radial contact currents change the current partition between turns, Schnaubelt 2023 shows an additional free cut coefficient IC2 is needed; that case is outside our current scope. |
The two can coexist in the same tape stack:
The buffer splits the conductor into two halves (requiring cut rerouting). The rint provides the resistive bridge for current transfer between Ag and YBCO (if present). Without rint, the YBCO half is only connected to the terminals via the cut constraint.
The implementation is more than adding one kernel: the kernel itself is small (§9.2), but the Maxwell machinery does not yet provide the support code for a collapsed interface (§9.1), and that support code dominates the effort.
A new domain type InterfaceTsCond already exists in en_DomainType.hpp:62, but it is not wired through any dispatch tables. The following sites all need a new case or branch:
| Site | File:line | What's missing |
|---|---|---|
| String parser | en_DomainType.cpp:97 | No "contact impedance" or "interfacetscond" string entry |
| Topology production | cl_Topology.cpp:270 | Not produced by sideset_type(master, slave) |
| Domain parsing | cl_FEM_Domain.cpp:28 | Not handled in the input-section dispatch |
| IWG dispatch | cl_IWG_Maxwell.cpp:558 | Not dispatched to any kernel function |
| DOF allocation | cl_Maxwell_FieldList.cpp:411 | Not given DOFs in the sideset DOF table |
Additionally:
| Component | Effort | Risk |
|---|---|---|
| Plumbing (§9.1): domain type wiring, topology, DOFs | 1-2 days | Medium (many files, all small) |
| Plumbing (§9.1): interface-property path design + implementation | 1 day | High (new design decision) |
| Plumbing (§9.1): ThinShellFactory skip-extrusion for rint | Half a day | Medium |
| Kernel (§9.2): h_contact_impedance() | Half a day | Low (simple bilinear form) |
| Consistency-term derivation (§3.4 gate) | 1 day | Medium (may not be needed, but must be checked) |
| Testing: analytical test problem + tape-stack validation | 1 day | Low |
| Total | 4-6 days |
The working value ρ = 1e-3 Ω·m used in the conditioning estimates above sits at the extreme high end of the measured range. The actual values depend critically on which interface is being modeled.
This is the metallurgical bond between the silver cap and the YBCO film, the "rint" in our tape model. Measured contact resistance per unit area R_ct:
| Source | R_ct [Ω·m²] | Notes |
|---|---|---|
| Hayasaka & Ito 2019 | 7e-12 to 3.7e-11 | Contact-probing current transfer length method, temperature-dependent |
| Soldered overlap measurements | 2.5e-12 to 5e-12 | Indium-soldered overlap |
| As-supplied coated conductors | ~2.5e-12 | ~25 nΩ·cm², Wimbush & Strickland data |
| Cu-stabilized REBCO lap joints | ~3.6e-12 | ~36 nΩ·cm², dominated by tape interface itself |
Typical range: R_ct ≈ 1e-12 to 1e-10 Ω·m²
Converting to bulk resistivity for a modeled layer thickness h = 0.1 μm:
This is mechanical contact between tape surfaces under winding pressure, a completely different physical interface. Typical values R_ct ≈ 1e-9 to 1e-7 Ω·m² (10-100 μΩ·cm²), i.e. 100-10000× higher than the internal Ag/YBCO interface. This is the application where the contact-impedance TSA (Schnaubelt 2023) is essential.
The conditioning ratio ρ/h against adjacent Ag (ρ_Ag/h_Ag ≈ 2e-3):
| Scenario | ρ_rint [Ω·m] | ρ/h | Contrast vs Ag | Assessment |
|---|---|---|---|---|
| Typical internal contact | 1e-4 | 1e3 | ~5e5 | MUMPS handles this routinely |
| Moderate contact | 1e-3 | 1e4 | ~5e6 | Challenging, may need regularization |
| NI coil T2TCL | 1e-1 to 1 | 1e6 to 1e7 | ~5e8 to 5e9 | Requires TSA or resistivity cap |
The resistivity cap at ρ_eff = 1e-4 Ω·m is both Alves 2022b's literature limit for what the pure H-formulation can handle and the physical value for typical REBCO tapes. It coincides with the measured Ag/YBCO interface resistivity for a typical coated conductor (R_ct ≈ 1e-11 Ω·m², h ≈ 0.1 μm → ρ ≈ 1e-4 Ω·m).
This means:
A resistivity of 1e-3 Ω·m corresponds to R_ct ≈ 1e-10 Ω·m², the upper extreme of the measured range. This occurs in:
For these cases, the conditioning ratio rises to ~5e6, which is more challenging but may still be manageable with the existing regularized h_ghost() approach.
The collapsed-interface formulation (§9) is worth the 4-6 day implementation effort only for:
For the current work (single-tape transport current at 77 K): the resistivity cap at ρ_eff = 1e-4 Ω·m is sufficient and physically correct for typical REBCO tapes. The existing h_ghost() with regularized harmonic mean handles the conditioning. For the internal Ag/YBCO contact, the regularized h_ghost() path may already solve the ghost problem.
For future work (NI coils, quench, multi-tape stacks): the contact-impedance formulation becomes necessary. The theory in §2-7 of this note and the implementation scope in §9 provide the starting point.
Immediate action: verify that the regularized h_ghost() at ρ = 1e-4 Ω·m produces the correct current-transfer behavior in the single-tape test case. If convergence is clean and the physics is right, the collapsed-interface implementation (§9) can be deprioritized indefinitely.
1.18.0