Skip to content
Quanscient logo Quanscient Allsolve Docs

𝑯-πœ‘ formulation

Maxwell’s magnetoquasistatic model

The formulation is derived based on 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}\\ \nabla \times \boldsymbol{H} &= \boldsymbol{J}\\ \nabla \cdot \boldsymbol{B} &= \boldsymbol{0} \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 constituve relations as

B=μHE=ρJ,\begin{align} \boldsymbol{B}&=\mu \boldsymbol{H}\\ \boldsymbol{E}&=\rho \boldsymbol{J}, \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

The strong formulation can be obtained by expressing (1) and (3) with H\boldsymbol{H} using (2), (4) and (5), to obtain

βˆ‡Γ—(ΟΒ βˆ‡Γ—H)+βˆ‚t(ΞΌH)=0βˆ‡β‹…(ΞΌH)=0\begin{align} \nabla \times (\rho~\nabla \times \boldsymbol{H})+\partial_t(\mu\boldsymbol{H})&=\boldsymbol{0}\\ \nabla \cdot (\mu\boldsymbol{H})=\boldsymbol{0} \end{align}
Decomposition of H\boldsymbol{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.

β– \blacksquare 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.

β– \blacksquare 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}.

β– \blacksquare 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.

β– \blacksquare 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.

Variational formulation

To obtain the variational formulation of the problem, (6) and (7) are multiplied 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 integrated 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\\ \int_\Omega \left( \nabla \cdot \mu \boldsymbol{H}\right )\cdot\varphi'~\rm{d}V&=0 \end{align}

The final form of the variational 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}')\\ \left( \nabla \cdot \boldsymbol{B}\right) \cdot\varphi'&=\boldsymbol{B}\cdot (\nabla \varphi')-\nabla \cdot (\boldsymbol{B}\varphi') \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\\ \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 \end{align}

where the boundary integral terms are the Neumann boundary terms.

β– \blacksquare The more detailed formulation is obtained by substituting the decompositions of H\boldsymbol{H} and Hβ€²\boldsymbol{H}' into (12) and (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 &+\\ \int_{\Omega_c} \mu\partial_t(\boldsymbol{H}_c+\boldsymbol{H}_s - \nabla \varphi)\cdot (\boldsymbol{H}_c'+\boldsymbol{H}_s')~\rm{d}V&+\\ \int_{\Omega_{nc}} \mu\partial_t(\boldsymbol{H}_s - \nabla \varphi)\cdot \boldsymbol{H}_s'~\rm{d}V&+\\ \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&+\\ \int_{\partial \Omega_c} (\boldsymbol{n} \times \boldsymbol{E})\cdot(\boldsymbol{H}_c'+\boldsymbol{H}_s')~\rm{d}A&+\\-\int_{\partial\Omega}(\boldsymbol{n}\cdot \boldsymbol{B})\cdot \varphi' ~\rm{d}A=0 \end{align}

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