View source for Magnetodynamics with cohomology conditions

{{metamodelGetDP|indheat}}

{{mytexdefs}}

== Additional information ==

The model represents an induction heating eddy current problem that utilizes the homology and cohomology solver of Gmsh <ref name=Pellika2013 />. In this example application, alternating current fed into the inductor coil generates eddy currents in a workpiece inside the coil, causing the workpiece to heat up. To run the example, open '''indheat.pro''' with Gmsh. 

The electromagnetic modeling aspects of the problem turn out to be subtle. The so-called $A-V$ formulation of the problem is straightforward to implement, but it results in a large linear system that might be difficult to solve. In the $T-\Omega$ formulation of the problem, the same accuracy is achieved with a smaller linear system, but its implementation involves so-called thick-cuts, or source fields, that aren't discussed much in the finite element curriculum and may be difficult to produce. Here, we call them cohomology basis functions, and generate them using the cohomology solver implemented in Gmsh. The role of the cohomology basis functions is first discussed in the $T-\Omega$ formulation and then for completeness, the cohomology functions are also expoited in the $A-V$ formulation to take an another perspective on their role in the electromagnetic modeling.

=== Problem definition ===

==== The domain ====

[[File:indheatdomain.png|200px|thumb|alt=Domain|Topology of the modeling domain.]]

Let $M \subset \mathbb{R}^3$ and let $\partial M = S_1 \cup S_2$ so that $\partial S_1 = \partial S_2 = S_1 \cap S_2$ denote the 3D modeling domain and its 2D boundary that is decomposed in two parts. Furthermore, the domain $M$ is decomposed in a conducting subdomain $M_c$ and a non-conducting subdomain $M_a$ so that $M = M_c \cup M_a$ and $M_c \cap M_a = \partial M_c \cap \partial M_a$ hold. We assume that $M$ is connected and has no holes nor voids, i.e. its [http://en.wikipedia.org/wiki/Betti_number Betti numbers] are $b_0(M)$ = 1 and $b_1(M) = b_2(M) = 0$.

The conducting subdomain $M_c$ consists of two parts: the workpiece being heated is a hollow cylinder around which the inductor coil wraps. The terminals of the inductor coil extend to the boundary $\partial M$ of the domain.

==== The boundary value problem ====

We solve time harmonic magnetoquasistatic approximation to the [[Maxwell's equations]] in $M$:

\begin{align} 
\Curl{\vec{h}} &= \vec{j}, & \Curl{\vec{e}} + i \omega \vec{b} &= 0 , \\
\Div{\vec{j}} &= 0 , & \Div{\vec{b}} &= 0, \\
\end{align}
with constitutive equations
\begin{align} 
\vec{b} &= \mu \vec{h}  , &
\vec{e} &= \rho \vec{j} .
\end{align}

[[File:indheatdomaincond.png|200px|thumb|alt=Boundary conditions|Applied boundary conditions.]]

In order to treat the domain $M$ as a circuit element, we apply magnetic isolation at its boundary $\partial M$. This is achieved by the boundary condition 
\begin{align} 
\vec{b} \cdot \vec{n} = 0 \textrm{ on } \partial M.
\end{align}

In the non-conducting domain $M_a$ we have no currents: 
\begin{align}
\Curl{\vec{h}} = 0 \textrm{ in } M_a,
\end{align}
and the current cannot pass from the conducting subdomain $M_c$ to the non-conducting subdomain $M_a$. Also, the current cannot pass through the boundary $\partial M$ to the non-conducting subdomain $M_a$. Thus, the conditions
\begin{align} 
\vec{j} \cdot \vec{n} = 0 \textrm{ on } \partial M_a \cap \partial M_c \textrm{ and on } M_a \cap \partial M = S_2
\end{align}
are required to hold at the interfaces. The terminals on $M_c \cap \partial M$ which connect the inductor coil to an external circuit are modeled as equipotential surfaces, i.e. as perfectly conducting surfaces. Thus, we apply a boundary condition 
\begin{align} 
\vec{e} \times \vec{n} = 0 \textrm{ on } M_c \cap \partial M = S_1
\end{align}
on them.

[[File:indheatdomaincoh.png|200px|thumb|alt=Cohomology conditions|Curves and a surface over which the cohomology conditions are specified.]]

Note that the boundary $\partial M$ is now decomposed on two parts according to the boundary condition applied on it: $\partial M = S_1 \cap S_2$. On $S_1 = M_c \cap \partial M$ we have the condition $\vec{e} \times \vec{n} = 0$ and on $S_2 = M_a \cap \partial M$ we have the condition $\vec{j} \cdot \vec{n} = 0$. We can either drive the problem by the net current $I_s$ through the terminals, or by the net voltage difference $V_s$ between the terminals:
\begin{align} 
\int_{\zeta} \vec{e} \cdot \mathrm{d} \vec{l} &= V_s, \\
\int_{\Sigma} \vec{j} \cdot \vec{n} \, \mathrm{d} a  = \int_{\Sigma} \Curl{\vec{h}} \cdot \vec{n} \, \mathrm{d} a = \int_{\partial \Sigma} \vec{h} \cdot \mathrm{d} \vec{l} &= I_s,  \\
\end{align}
where $\zeta$ is a curve on $S_2$ between the terminals on $S_1$ and $\Sigma \subset M$ is a surface that isolates the terminals. The boundary $\partial \Sigma$ also lies on the part $S_2$ of the boundary of the domain. In addition, the fields $\vec{e}$ and $\vec{j}$ are linked by the constitutive equation $\vec{e} = \rho \vec{j}$. All these considerations hint a duality between the electric field $\vec{e}$ and the current density $\vec{j}$ and between the voltage and the current conditions.

A [http://en.wikipedia.org/wiki/Equivalence_class class] $[\zeta]$ of curves on $S_2$ between the terminals on $S_1$ belong to so-called relative [http://en.wikipedia.org/wiki/Simplicial_homology homology] space $H_1(S_2, \partial S_2)$, while the class $[\partial \Sigma]$ of curves belong to a homology space $H_1(S_2)$. The integral conditions on fields $\vec{e}$ and $\vec{j}$ (or $\vec{h}$) fixing $V_s$ and $I_s$ are called [http://en.wikipedia.org/wiki/De_Rham_cohomology cohomology] conditions, since they consider integrals over an element of a homology space.

==== $T-\Omega$ potential formulation ====

In this formulation we solve for the magnetic field $\vec{h}$ in $M$ to obtain the current density $\vec{j} = \Curl{\vec{h}}$ in $M_c$. Let us 
denote by $N(M_a)$ the set of nodes in the mesh of $M_a$ and by $E(M_c \setminus (\partial M_a \cap \partial M_c))$ the set of edges in $M_c$ that are not on the interface $\partial M_a \cap \partial M_c$ of conducting and non-conducting subdomains. 

Denote by $\mathsf{n}_i$ nodal shape functions associated with the nodes of the mesh and by $\mathsf{e}_i$ edge elements associated with the edges of the mesh. By $\mathsf{E}_i$ we denote so-called 1-cohomology basis function, which is an integer combination of edge elements on $M$:
$$
\mathsf{E}_i = \sum_j z^j_i \mathsf{e}_j
$$
where the integer coefficient vector $\vec{z}_i$ is produced by the cohomology solver of Gmsh.

[[File:indheati1.png|200px|thumb|alt=Cohomology basis function|Edges that have non-zero coefficient in a cohomology basis function $\mathsf{E}_1 \in H^1(M_a)$ whose coefficient $I_1$ specifies the current in the inductor.]]
[[File:indheati2.png|200px|thumb|alt=Cohomology basis function|Edges that have non-zero coefficient in a cohomology basis function $\mathsf{E}_2 \in H^1(M_a)$ whose coefficient $I_2$ specifies the current in the workpiece.]]

Now, the unknown magnetic field $\vec{h}$ is of the form 
$$
\vec{h} = \sum_{i \in N(M_a)} \phi_i \,\Grad{\mathsf{n}_i} + \sum_{i \in E(M_c \setminus (\partial M_a \cap \partial M_c))} t_i \,\mathsf{e}_i + I_1 \mathsf{E}_1 + I_2 \mathsf{E}_2,
$$
where the cohomology basis functions $\mathsf{E}_1$ and $\mathsf{E}_2$ are produced by the computation of the cohomology space $H^1(M_a)$ in Gmsh. That is, the function space from which we look for the approximate solution of $\vec{h}$ is spanned by nodal shape functions $\mathsf{n}_i$ of $M_a$, edge elements $\mathsf{e}_i$ of $M_c \setminus (\partial M_a \cap \partial M_c)$, and by the cohomology basis functions $\mathsf{E}_i$ of the cohomology space $H^1(M_a)$. Those functions shall also serve as the test functions in the following weak formulation.

The weak formulation is
\begin{align}
-\int_M \mu \,\vec{h} \cdot \vec{h}' \, \mathrm{d} V &= 0 \quad \forall \, \vec{h}' \\
\int_{M_c} \rho \,\Curl{\vec{h}} \cdot \Curl{\vec{h}'} \, \mathrm{d} V + i \omega \int_M \mu \,\vec{h} \cdot \vec{h}' \, \mathrm{d} V &= -
\int_{S_2} \vec{e} \times \vec{h}' \cdot \vec{n} \, \mathrm{d} a \quad \forall \, \vec{h}',
\end{align}
where 
\begin{align}
\vec{h}' \in \{ \Grad{\mathsf{n}_i}, \mathsf{e}_i, \mathsf{E}_1, \mathsf{E}_2\} 
& \textrm{ and } &
\int_{S_2} \vec{e} \times \vec{h}' \cdot \vec{n} \, \mathrm{d} a = 
\begin{cases} V_i & \textrm{ when } \vec{h}' = \mathsf{E}_i \\ 
0 & \textrm{ otherwise. }
\end{cases}
\end{align}
By plugging in the expression for $\vec{h}$ one can construct a linear system from which the unknown coefficients $\phi_i$, $t_i$ and $V_i$ or $I_i$ can be solved. Note that either $V_i$ or $I_i$ from both pairs $(V_1, I_1)$ and $(V_2, I_2)$ must be fixed in order to get a unique solution.

==== $A-V$ potential formulation ====

In this formulation we solve for the vector potential $\vec{a}$ of the magnetic flux density $\vec{b}$ in $M$ and for the curl-free part $\boldsymbol{\epsilon}$ of the electric field $\vec{e}$ in $M_c$ to obtain the current density $\vec{j} = \rho^{-1} \vec{e} = \rho^{-1} (i \omega \vec{a} + \boldsymbol{\epsilon})$ in $M_c$.

The unknown magnetic vector potential $\vec{a}$ is of the form
$$
\vec{a} = \sum_{i \in E(M  \setminus \partial M)} a_i \mathsf{e}_i.
$$
That is, the function space from which we look for the approximate solution of $\vec{a}$ is spanned by edge elements $\mathsf{e}_i$ of $M \setminus \partial M$. Those functions shall also serve as the test functions in the following weak formulation.

[[File:indheatv1.png|200px|thumb|alt=Cohomology basis function|Edges that have non-zero coefficient in a cohomology basis function $\mathsf{E}_1 \in H^1(M_c, S_1)$ whose coefficient $V_1$ specifies the voltage across the inductor.]]
[[File:indheatv2.png|200px|thumb|alt=Cohomology basis function|Edges that have non-zero coefficient in a cohomology basis function $\mathsf{E}_2 \in H^1(M_c, S_1)$ whose coefficient $V_2$ specifies the voltage across the workpiece.]]

The unknown curl-free part $\boldsymbol{\epsilon}$ of the electric field is of the form
$$
\boldsymbol{\epsilon} = \sum_{i \in N(M_c \setminus S_1)} \varphi_i \Grad{\mathsf{n}_i} + V_1 \mathsf{E}_1 + V_2 \mathsf{E}_2,
$$
Where the cohomology basis functions $\mathsf{E}_1$ and $\mathsf{E}_2$ are produced by the computation of the relative cohomology space $H^1(M_c, S_1)$ in Gmsh. That is, the function space from which we look for the approximate solution of $\boldsymbol{\epsilon}$ is spanned by nodal shape functions $\mathsf{n}_i$ of $M_c \setminus S_1$, and by the cohomology basis functions $\mathsf{E}_i$ of the cohomology space $H^1(M_c, S_1)$. Those functions shall also serve as the test functions in the following weak formulation.

The weak formulation is
\begin{align}
\int_M \mu^{-1} \Curl{\vec{a}} \cdot \Curl{\vec{a}}' \, \mathrm{d} V + \int_{M_c} \rho^{-1} \boldsymbol{\epsilon} \cdot \vec{a}' \, \mathrm{d} V + i \omega \int_{M_c} \rho^{-1} \vec{a} \cdot \vec{a}' \, \mathrm{d} V &= 0 \quad \forall \, \vec{a}' \\
\int_{M_c} \rho^{-1} \boldsymbol{\epsilon} \cdot \boldsymbol{\epsilon}' \, \mathrm{d} V + i \omega \int_{M_c} \rho^{-1} \vec{a} \cdot \boldsymbol{\epsilon}' \, \mathrm{d} V &= -
\int_{S_1} \vec{j} \varphi' \, \mathrm{d} a  \quad \forall \, \boldsymbol{\epsilon}'
\end{align}
where 
\begin{align}
\vec{a}' \in \{\mathsf{e}_i\}, \quad \boldsymbol{\epsilon}' \in \{ \Grad{\mathsf{n}_i}, \mathsf{E}_1, \mathsf{E}_2\},
& \textrm{ and } &
\int_{S_1} \vec{j} \varphi' \, \mathrm{d} a = 
\begin{cases} I_i & \textrm{ when } \boldsymbol{\epsilon}' = \mathsf{E}_i \\ 
0 & \textrm{ otherwise. }
\end{cases}
\end{align}
By plugging in the expressions for $\vec{a}$ and $\boldsymbol{\epsilon}$ one can construct a linear system from which the unknown coefficients $a_i$, $\varphi_i$, and $V_i$ or $I_i$ can be solved. Note that either $V_i$ or $I_i$ from both pairs $(V_1, I_1)$ and $(V_2, I_2)$ must be fixed in order to get a unique solution.

==== Cohomology conditions ====

In both formulations, if the cohomology basis function $\mathsf{E}_1$ corresponds to the net current $I_s$ or the net voltage $V_s$ in the inductor coil, $I_1 = I_s$ and $V_1 = V_s$ hold. If cohomology basis function $\mathsf{E}_2$ corresponds to the net current or voltage in the workpiece, one should set $V_2 = 0$. That is, one drives the problem by setting either the source current $I_s$ or the source voltage $V_s$, and the other quantity together with the net current $I_2$ in the workpiece are solved as part of the finite element solution.

The implementation below produces such cohomology basis functions $\mathsf{E}_1$ and $\mathsf{E}_2$ as described above.

<!--

== Implementation == 

=== Indheat.geo: problem geometry and cohomology computation in Gmsh ===

{{GetDPTutorial|Magnetodynamics/GMSH_GETDP/indheat.geo|height=20em}}

=== Indheat.pro: weak formulation in GetDP ===

{{GetDPTutorial|Magnetodynamics/GMSH_GETDP/indheat.pro|height=20em}}

=== How to use ===

All the files (.geo and .pro) must be located in the same directory. 

; Meshing the domain and computing the cohomology:
Go to the directory and then type:
 gmsh indheat.geo -3

After the mesh is built, a file "indheat.msh" should have been created in the directory.

; Solving the problem with GetDP:
In a Terminal, type (in the right directory)
 getdp indheat.pro -solve MagDynTOComplex -pos MagDynTO
to solve with $T-\Omega$ formulation, or 
 getdp indheat.pro -solve MagDynAVComplex -pos MagDynAV
to solve with $A-V$ formulation.

; Showing the result

Open the file "jTO.pos" or "jAV.pos" with Gmsh by typing "gmsh jTO.pos" or "gmsh jAV.pos" in a terminal in the right directory.

=== Result ===
[http://onelab.info/files/getdp/Magnetodynamics/GMSH_GETDP/indheatanim.gif Animation in the time domain]
and real and imaginary part of the current density $\vec{j}$ in the conducting subdomain $M_c$, here computed using the $T-\Omega$ formulation:

[[File:indheatjr.jpg|500px|frameless|alt=Current real part|Real part of the current density in the conducting regions.]]
[[File:indheatji.jpg|500px|frameless|alt=Current imaginary part|Imaginary part of the current density in the conducting regions.]]

-->

== References ==

<references>

<ref name=Pellika2013> M. Pellikka, S. Suuriniemi, L. Kettunen and C. Geuzaine, [http://geuz.org/gmsh/doc/preprints/gmsh_homology_preprint.pdf Homology and cohomology computation in finite element modeling]. SIAM Journal on Scientific Computing 35(5), pp. 1195-1214, 2013.</ref>

</references>

{{metamodelfooter|indheat}}