H-𝜑 formulation

What you’ll learn

• Domain decomposition — splitting conducting and non-conducting regions for the $\boldsymbol{H}$-$\varphi$ 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 $\boldsymbol{H}$-$\varphi$ 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 \frac{1}{\sigma}(\nabla \times \boldsymbol{H})\cdot (\nabla \times \boldsymbol{H}')\ d\Omega + \int_\Omega \frac{\partial (\mu \boldsymbol{H})}{\partial t}\cdot\boldsymbol{H}'\ d\Omega + \int_{\Gamma_c} (\boldsymbol{n} \times \boldsymbol{E})\cdot\boldsymbol{H}'\ d\Gamma_c &=0 \tag{1}\\ \int_\Omega (\mu \boldsymbol{H})\cdot (\nabla \varphi')\ d\Omega - \int_{\Gamma}(\boldsymbol{n}\cdot \boldsymbol{B})\cdot \varphi'\ d\Gamma &=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} + \frac{\partial \boldsymbol{B}}{\partial t} &= \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{J} &= \sigma \boldsymbol{E}, \tag{7} \end{align*}$$

where $\mu$ is the magnetic permeability and $\sigma$ is the electric conductivity. Note that, $\mu$ and $\sigma$ 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 \left(\frac{1}{\sigma}~\nabla \times \boldsymbol{H}\right) + \frac{\partial (\mu\boldsymbol{H})}{\partial t} &= \boldsymbol{0} \tag{8}\\ \nabla \cdot (\mu\boldsymbol{H}) &= \boldsymbol{0} \tag{9} \end{align*}$$

Decomposition of $\boldsymbol{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 $\Gamma_{nc}$.

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

• $\boldsymbol{H}_c$ should be constrained to zero at parts of $\Gamma_c$ where no current is desired to flow through.
• At parts of $\Gamma_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 $\Gamma_{nc}$ where magnetic flux density is desired to flow perpendicularly through the surface.
• At parts of $\Gamma_{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 (\nabla \times \boldsymbol{E})\cdot\boldsymbol{H}'\ d\Omega + \int_\Omega \frac{\partial (\mu \boldsymbol{H})}{\partial t}\cdot\boldsymbol{H}'\ d\Omega &= 0 \tag{10}\\ \int_\Omega (\nabla \cdot \mu \boldsymbol{H})\cdot\varphi'\ d\Omega &= 0 \tag{11} \end{align*}$$

For the first equation, using the divergence of a cross product we can rewrite the first term

$$\begin{align*} \int_\Omega \nabla \cdot (\boldsymbol{E} \times \boldsymbol{H}')\ d\Omega + \int_\Omega \boldsymbol{E} \cdot (\nabla \times \boldsymbol{H}')\ d\Omega + \int_\Omega \frac{\partial (\mu \boldsymbol{H})}{\partial t}\cdot\boldsymbol{H}'\ d\Omega &= 0. \tag{12} \end{align*}$$

Applying the divergence theorem on the divergence term we obtain

$$\begin{align*} \int_{\Gamma} (\boldsymbol{E} \times \boldsymbol{H}') \cdot \boldsymbol{n}\ d\Gamma + \int_\Omega \boldsymbol{E} \cdot (\nabla \times \boldsymbol{H}')\ d\Omega + \int_\Omega \frac{\partial (\mu \boldsymbol{H})}{\partial t}\cdot\boldsymbol{H}'\ d\Omega &= 0. \tag{13} \end{align*}$$

Making use of the scalar triple product rule on the boundary term and rearranging results in

$$\begin{align*} \int_\Omega \boldsymbol{E} \cdot (\nabla \times \boldsymbol{H}')\ d\Omega + \int_\Omega \frac{\partial (\mu \boldsymbol{H})}{\partial t}\cdot\boldsymbol{H}'\ d\Omega + \int_{\Gamma_c} (\boldsymbol{n} \times \boldsymbol{E}) \cdot \boldsymbol{H}'\ d\Gamma &= 0 \tag{14} \end{align*}$$

and by substituting $\boldsymbol{E} = \frac{1}{\sigma}(\nabla \times \boldsymbol{H})$ into the first term we derive the first coupled equation stated at the top of this page.

For the second equation, using the Leibniz rule for a nabla operator on $\text{(11)}$ we get

$$\begin{align*} \int_\Omega (\mu \boldsymbol{H}) \cdot (\nabla \varphi')\ d\Omega - \int_\Omega \nabla \cdot (\boldsymbol{B}\varphi')\ d\Omega &= 0. \tag{15} \end{align*}$$

Applying the divergence theorem on the divergence term we obtain the second coupled equation 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{(8)}$ and $\text{(9)}$. For the first equation we have

$$\begin{align*} &\int_{\Omega_c} \frac{1}{\sigma}(\nabla \times \boldsymbol{H}_c+\nabla \times \boldsymbol{H}_s)\cdot (\nabla \times \boldsymbol{H}_c'+\nabla \times \boldsymbol{H}_s')\ d\Omega \\ &+ \int_{\Omega_c} \mu\frac{\partial}{\partial t}(\boldsymbol{H}_c+\boldsymbol{H}_s - \nabla \varphi)\cdot (\boldsymbol{H}_c'+\boldsymbol{H}_s')\ d\Omega \\ &+ \int_{\Omega_{nc}} \mu\frac{\partial}{\partial t}(\boldsymbol{H}_s - \nabla \varphi)\cdot \boldsymbol{H}_s'\ d\Omega \\ &+ \int_{\Gamma_c} (\boldsymbol{n} \times \boldsymbol{E})\cdot(\boldsymbol{H}_c'+\boldsymbol{H}_s')\ d\Gamma =0 \tag{16} \end{align*}$$

For the second equation we obtain

$$\begin{align*} &\int_{\Omega_c} \mu (\boldsymbol{H}_c+\boldsymbol{H}_s-\nabla\varphi)\cdot (\nabla \varphi')\ d\Omega \\ &+ \int_{\Omega_{nc}} \mu (\boldsymbol{H}_s-\nabla\varphi)\cdot (\nabla \varphi')\ d\Omega \\ &- \int_{\Gamma}(\boldsymbol{n}\cdot \boldsymbol{B})\cdot \varphi'\ d\Gamma =0 \tag{17} \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