H-𝜑 formulation

What you’ll learn

• Domain decomposition — splitting conducting and non-conducting regions for the H-𝜑 formulation.
• Field decomposition — how $\boldsymbol{H}$, $\varphi$, and cohomology source fields $\boldsymbol{H}_s$ fit together.
• Boundary conditions — Neumann terms and constraints for $\boldsymbol{H}_c$ and $\varphi$.
• Weak formulation — the coupled variational equations used for HTS and electric motor simulations.

This page describes 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. The weak formulation is stated up front; the full derivation is at the end of this page.

Weak formulation

The weak form solved by Allsolve consists of two coupled variational equations:

$$\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{1}\\ \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{2} \end{align*}$$

where $\boldsymbol{H}$ is decomposed into conducting and non-conducting parts, $\varphi$ is the magnetic scalar potential, and the boundary integral terms are the Neumann boundary terms.

See Formulation derivation for the full strong-to-weak derivation, domain decomposition, and the expanded form used in Allsolve.

Formulation derivation

Maxwell’s magnetoquasistatic model

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

$$\begin{align*} \nabla \times \boldsymbol{E} +\partial_t \boldsymbol{B} &= \boldsymbol{0} \tag{3}\\ \nabla \times \boldsymbol{H} &= \boldsymbol{J} \tag{4}\\ \nabla \cdot \boldsymbol{B} &= \boldsymbol{0} \tag{5} \end{align*}$$

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

$$\begin{align*} \boldsymbol{B}&=\mu \boldsymbol{H} \tag{6}\\ \boldsymbol{E}&=\rho \boldsymbol{J}, \tag{7} \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

Express $\text{(3)}$ and $\text{(5)}$ in terms of $\boldsymbol{H}$ using $\text{(4)}$, $\text{(6)}$, and $\text{(7)}$:

$$\begin{align*} \nabla \times (\rho~\nabla \times \boldsymbol{H})+\partial_t(\mu\boldsymbol{H})&=\boldsymbol{0} \tag{8}\\ \nabla \cdot (\mu\boldsymbol{H})&=\boldsymbol{0} \tag{9} \end{align*}$$

Decomposition of H

The $\boldsymbol{H}-\varphi$ formulation is obtained by first decomposing the whole modeling domain $\Omega$ into conducting domain $\Omega_{c}$ and non-conducting domain $\Omega_{nc}$ such that $\Omega=\Omega_{c} \cup \Omega_{nc}$.

Then, $\boldsymbol{H}$ is decomposed as

$$\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 $\boldsymbol{H}_c$ is the magnetic field strength in the conducting domain. The cohomology source field $\boldsymbol{H}_s$ is used, for example, to impose the total current to the High Temperature Superconducting (HTS) tape. Moreover, $\boldsymbol{H}_c$ and $\boldsymbol{H}_s$ are interpolated with edge elements and $\varphi$ with nodal elements.

• Note that $\varphi$ and $\boldsymbol{H}_s$ have degrees of freedom only in $\Omega_{nc}$ but their support is the whole $\Omega$. This allows the coupling of $\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 $\Omega_{nc}$.

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

• $\boldsymbol{H}_c$ should be constrained to zero at parts of the $\partial \Omega_c$ where no current is desired to flow through.
• At parts of $\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 $\partial \Omega_{nc}$ where magnetic flux density is desired to flow perpendicularly through the surface.
• At parts of $\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.

Deriving the weak form

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

$$\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

$$\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{10}\\ \int_\Omega \left( \nabla \cdot \mu \boldsymbol{H}\right )\cdot\varphi'~\rm{d}V&=0 \tag{11} \end{align*}$$

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

$$\begin{align*} (\nabla \times \boldsymbol{E})\cdot \boldsymbol{H}'&=\nabla \cdot (\boldsymbol{E} \times \boldsymbol{H}')+\boldsymbol{E}\cdot (\nabla \times \boldsymbol{H}') \tag{12}\\ \left( \nabla \cdot \boldsymbol{B}\right) \cdot\varphi'&=\boldsymbol{B}\cdot (\nabla \varphi')-\nabla \cdot (\boldsymbol{B}\varphi') \tag{13} \end{align*}$$

to obtain the coupled equations stated at the top of this page. The more detailed formulation is obtained by substituting the decompositions of $\boldsymbol{H}$ and $\boldsymbol{H}'$ into $\text{(1)}$ and $\text{(2)}$ to obtain

$$\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

[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