Skip to content
Quanscient logo Quanscient Allsolve Docs

H-πœ‘ formulation

This page derives the H-πœ‘ mixed formulation used for magnetoquasistatic problems that combine conducting and non-conducting domains. The formulation decomposes the magnetic field into edge-element and scalar-potential parts, enabling efficient solution of superconductor and electric motor simulations.

Maxwell’s magnetoquasistatic model

Section titled β€œMaxwell’s magnetoquasistatic model”

The formulation is derived from the magnetoquasistatic approximation of Maxwell’s equations:

βˆ‡Γ—E+βˆ‚tB=0βˆ‡Γ—H=Jβˆ‡β‹…B=0\begin{align*} \nabla \times \boldsymbol{E} +\partial_t \boldsymbol{B} &= \boldsymbol{0} \tag{1}\\ \nabla \times \boldsymbol{H} &= \boldsymbol{J} \tag{2}\\ \nabla \cdot \boldsymbol{B} &= \boldsymbol{0} \tag{3} \end{align*}

The electric field strength E\boldsymbol{E}, the electric current density J\boldsymbol{J}, the magnetic flux density B\boldsymbol{B} and the magnetic field strength H\boldsymbol{H} are paired by the constitutive relations as

B=μHE=ρJ,\begin{align*} \boldsymbol{B}&=\mu \boldsymbol{H} \tag{4}\\ \boldsymbol{E}&=\rho \boldsymbol{J}, \tag{5} \end{align*}

where μ\mu is the magnetic permeability and ρ\rho is the electric resistivity. Note that, μ\mu and ρ\rho can depend on different quantities such as electromagnetic fields and temperature.

Strong formulation

Section titled β€œStrong formulation”

Express (1)\text{(1)} and (3)\text{(3)} in terms of H\boldsymbol{H} using (2)\text{(2)}, (4)\text{(4)}, and (5)\text{(5)}:

βˆ‡Γ—(ΟΒ βˆ‡Γ—H)+βˆ‚t(ΞΌH)=0βˆ‡β‹…(ΞΌH)=0\begin{align*} \nabla \times (\rho~\nabla \times \boldsymbol{H})+\partial_t(\mu\boldsymbol{H})&=\boldsymbol{0} \tag{6}\\ \nabla \cdot (\mu\boldsymbol{H})&=\boldsymbol{0} \tag{7} \end{align*}

Decomposition of H\boldsymbol{H}

Section titled β€œDecomposition of H\boldsymbol{H}H”

The Hβˆ’Ο†\boldsymbol{H}-\varphi formulation is obtained by first decomposing the whole modeling domain Ξ©\Omega into conducting domain Ξ©c\Omega_{c} and non-conducting domain Ξ©nc\Omega_{nc} such that Ξ©=Ξ©cβˆͺΞ©nc\Omega=\Omega_{c} \cup \Omega_{nc}.

Then, H\boldsymbol{H} is decomposed as

H={Ξ©c:Hcβˆ’βˆ‡Ο†+HsΞ©nc:βˆ‡Ο†+Hs\boldsymbol{H}=\begin{cases} \Omega_c: &\boldsymbol{H}_c-\nabla \varphi+\boldsymbol{H}_s\\ \Omega_{nc}: &\nabla \varphi + \boldsymbol{H}_s \end{cases}

where Ο†\varphi is a scalar field and Hc\boldsymbol{H}_c is the magnetic field strength in the conducting domain. The cohomology source field Hs\boldsymbol{H}_s is used, for example, to impose the total current to the High Temperature Superconducting (HTS) tape. Moreover, Hc\boldsymbol{H}_c and Hs\boldsymbol{H}_s are interpolated with edge elements and Ο†\varphi with nodal elements.

  • Note that Ο†\varphi and Hs\boldsymbol{H}_s have degrees of freedom only in Ξ©nc\Omega_{nc} but their support is the whole Ξ©\Omega. This allows the coupling of H\boldsymbol{H} through the boundary of the conducting domain.
  • The scalar field Ο†\varphi should be made unique by constraining the value of Ο†\varphi to zero at a point on the boundary of Ξ©nc\Omega_{nc}.

To impose the boundary conditions for the vector field Hc\boldsymbol{H}_c:

  • Hc\boldsymbol{H}_c should be constrained to zero at parts of the βˆ‚Ξ©c\partial \Omega_c where no current is desired to flow through.
  • At parts of βˆ‚Ξ©c\partial \Omega_c, where current density is desired flow pass perpendicularly through the surface, the Neumann boundary term is zero.

To impose the boundary conditions for the scalar field Ο†\varphi:

  • Ο†\varphi should be constrained to zero at parts of the βˆ‚Ξ©nc\partial \Omega_{nc} where magnetic flux density is desired to flow perpendicularly through the surface.
  • At parts of βˆ‚Ξ©nc\partial \Omega_{nc}, where zero magnetic flux density is desired to pass through, the Neumann boundary term should be equal to zero.
    • The Neumann boundary term is explained in the following section.

Weak formulation

Section titled β€œWeak formulation”

To obtain the weak (variational) formulation, multiply (6)\text{(6)} and (7)\text{(7)} by the test function fields:

Hβ€²={Ξ©c:Hcβ€²+βˆ‡Ο†β€²+Hsβ€²Ξ©nc:βˆ‡Ο†β€²+Hsβ€²\boldsymbol{H}'=\begin{cases} \Omega_c: &\boldsymbol{H}_c'+\nabla \varphi'+\boldsymbol{H}_s'\\ \Omega_{nc}: &\nabla \varphi' + \boldsymbol{H}_s' \end{cases}

and integrate over the whole modeling domain Ξ©\Omega to obtain

∫Ω(βˆ‡Γ—E)β‹…Hβ€²+βˆ‚t(ΞΌH)β‹…Hβ€²Β dV=0∫Ω(βˆ‡β‹…ΞΌH)β‹…Ο†β€²Β dV=0\begin{align*} \int_\Omega \left( \nabla \times \boldsymbol{E}\right )\cdot\boldsymbol{H}' +\partial_t(\mu\boldsymbol{H})\cdot\boldsymbol{H}'~\rm{d}V&=0 \tag{8}\\ \int_\Omega \left( \nabla \cdot \mu \boldsymbol{H}\right )\cdot\varphi'~\rm{d}V&=0 \tag{9} \end{align*}

The final form of the weak formulation used in Allsolve is obtained using the Stokes’ theorem and the identities

(βˆ‡Γ—E)β‹…Hβ€²=βˆ‡β‹…(EΓ—Hβ€²)+Eβ‹…(βˆ‡Γ—Hβ€²)(βˆ‡β‹…B)β‹…Ο†β€²=Bβ‹…(βˆ‡Ο†β€²)βˆ’βˆ‡β‹…(BΟ†β€²)\begin{align*} (\nabla \times \boldsymbol{E})\cdot \boldsymbol{H}'&=\nabla \cdot (\boldsymbol{E} \times \boldsymbol{H}')+\boldsymbol{E}\cdot (\nabla \times \boldsymbol{H}') \tag{10}\\ \left( \nabla \cdot \boldsymbol{B}\right) \cdot\varphi'&=\boldsymbol{B}\cdot (\nabla \varphi')-\nabla \cdot (\boldsymbol{B}\varphi') \tag{11} \end{align*}

to obtain

∫Ωρ(βˆ‡Γ—H)β‹…(βˆ‡Γ—Hβ€²)+βˆ‚t(ΞΌH)β‹…Hβ€²Β dV+βˆ«βˆ‚Ξ©c(nΓ—E)β‹…Hβ€²Β dA=0∫Ω(ΞΌH)β‹…(βˆ‡Ο†β€²)Β dVβˆ’βˆ«βˆ‚Ξ©(nβ‹…B)β‹…Ο†β€²Β dA=0\begin{align*} \int_\Omega \rho(\nabla \times \boldsymbol{H})\cdot (\nabla \times \boldsymbol{H}') +\partial_t(\mu\boldsymbol{H})\cdot\boldsymbol{H}'~\rm{d}V+\int_{\partial \Omega_c} (\boldsymbol{n} \times \boldsymbol{E})\cdot\boldsymbol{H}'~\rm{d}A &=0 \tag{12}\\ \int_\Omega (\mu \boldsymbol{H})\cdot (\nabla \varphi')~\rm{d}V-\int_{\partial\Omega}(\boldsymbol{n}\cdot \boldsymbol{B})\cdot \varphi' ~\rm{d}A&=0 \tag{13} \end{align*}

where the boundary integral terms are the Neumann boundary terms. The more detailed formulation is obtained by substituting the decompositions of H\boldsymbol{H} and Hβ€²\boldsymbol{H}' into (12)\text{(12)} and (13)\text{(13)} to obtain

∫Ωcρ(βˆ‡Γ—Hc+βˆ‡Γ—Hs)β‹…(βˆ‡Γ—Hcβ€²+βˆ‡Γ—Hsβ€²)Β dV+∫ΩcΞΌβˆ‚t(Hc+Hsβˆ’βˆ‡Ο†)β‹…(Hcβ€²+Hsβ€²)Β dV+∫ΩncΞΌβˆ‚t(Hsβˆ’βˆ‡Ο†)β‹…Hsβ€²Β dV+∫ΩcΞΌ(Hc+Hsβˆ’βˆ‡Ο†)β‹…(βˆ‡Ο†β€²)Β dV+∫ΩncΞΌ(Hsβˆ’βˆ‡Ο†)β‹…(βˆ‡Ο†β€²)Β dV+βˆ«βˆ‚Ξ©c(nΓ—E)β‹…(Hcβ€²+Hsβ€²)Β dA+βˆ’βˆ«βˆ‚Ξ©(nβ‹…B)β‹…Ο†β€²Β dA=0\begin{align*} \int_{\Omega_c} \rho(\nabla \times \boldsymbol{H}_c+\nabla \times \boldsymbol{H}_s)\cdot (\nabla \times \boldsymbol{H}_c'+\nabla \times \boldsymbol{H}_s')~\rm{d}V &+ \tag{14}\\ \int_{\Omega_c} \mu\partial_t(\boldsymbol{H}_c+\boldsymbol{H}_s - \nabla \varphi)\cdot (\boldsymbol{H}_c'+\boldsymbol{H}_s')~\rm{d}V&+ \tag{15}\\ \int_{\Omega_{nc}} \mu\partial_t(\boldsymbol{H}_s - \nabla \varphi)\cdot \boldsymbol{H}_s'~\rm{d}V&+ \tag{16}\\ \int_{\Omega_{c}} \mu (\boldsymbol{H}_c+\boldsymbol{H}_s-\nabla\varphi)\cdot (\nabla \varphi')~\rm{d}V+\int_{\Omega_{nc}} \mu (\boldsymbol{H}_s-\nabla\varphi)\cdot (\nabla \varphi')~\rm{d}V&+ \tag{17}\\ \int_{\partial \Omega_c} (\boldsymbol{n} \times \boldsymbol{E})\cdot(\boldsymbol{H}_c'+\boldsymbol{H}_s')~\rm{d}A&+ \tag{18}\\-\int_{\partial\Omega}(\boldsymbol{n}\cdot \boldsymbol{B})\cdot \varphi' ~\rm{d}A&=0 \tag{19} \end{align*}

References

Section titled β€œReferences”

[1] Lahtinen, V., Stenvall, A., Sirois, F. et al. A Finite Element Simulation Tool for Predicting Hysteresis Losses in Superconductors Using an H-Oriented Formulation with Cohomology Basis Functions. J Supercond Nov Magn 28, 2345–2354 (2015)

[2] J. Ruuskanen et al., β€œModeling Eddy Current Losses in HTS Tapes Using Multiharmonic Method,” in IEEE Transactions on Applied Superconductivity, vol. 33, no. 5, pp. 1-5, Aug. 2023, Art no. 5900605

[3] Pellikka, Matti, et al. β€œHomology and cohomology computation in finite element modeling.” SIAM Journal on Scientific Computing 35.5 (2013): B1195-B1214.

[4] Accelerating superconductivity simulations with the H-𝝓 formulation