Skip to content

Magnetic PDE

Maxwell's Equations

The set of Maxwell's equations is based on the experiments of Ampére, Gauss and Faraday. The novelty of Maxwell was to introduce the displacement current density in Ampére's law, which enabled him to describe the propagation of electromagnetic waves. The full set of this equations can be written in local form as \begin{equation} \label{eq:Ampere} \textbf{curl}( \mathbf{H}) = \mathbf{J} + \frac{\partial \mathbf{D}}{\partial t}, \end{equation}

\begin{equation} \label{eq:Faraday} \textbf{curl} (\mathbf{E}) = -\frac{\partial \mathbf{B}}{\partial t}, \end{equation}
\begin{equation} \label{eq:Gauss} \textbf{div} (\mathbf{D}) = q_e, \end{equation}
\begin{equation} \textbf{div} (\mathbf{B}) = 0. \label{eq:solenoidalB} \end{equation}

For the global form, we consider the arbitrary volume, depicted in the figure below and apply Gauss- and Stokes-theorem to (\ref{eq:Ampere})-(\ref{eq:solenoidalB}).

acoustic

Volume domain for the global form of Maxwell's equations

This leads to \begin{equation} \label{eq:Ampereglobal} \int_{\partial \Gamma} \mathbf{H}\cdot \mathbf \tau_{\Gamma} = \int_{\Gamma} \mathbf{J} \cdot \mathbf n_{\Gamma} + \int_{\Gamma} \frac{\partial \mathbf{D}}{\partial t} \cdot \mathbf n_{\Gamma} , \end{equation}

\begin{equation} \label{eq:Faradayglobal} \int_{\partial \Gamma} \mathbf{E} \cdot \mathbf \tau_{\Gamma} = - \int_{\Gamma} \frac{\partial \mathbf{B}}{\partial t} \cdot \mathbf n_{\Gamma} , \end{equation}
\begin{equation} \label{eq:Gaussglobal} \int_{\partial \Omega} \mathbf{D} \cdot \mathbf n_{\Omega} = \int_{\Omega} q_e, \end{equation}

\begin{equation} \int_{\partial \Omega} \mathbf{B} \cdot \mathbf n_{\Omega} = 0, \label{eq:solenoidalBglobal} \end{equation} which can be identified as the law of Ampére with Maxwell's contribution of the displacement current density (\ref{eq:Ampereglobal}), the law of Faraday (\ref{eq:Faradayglobal}), the law of Gauss (\ref{eq:Gaussglobal}) and the global solenoidal \mathbf B-field (\ref{eq:solenoidalBglobal}).

The equations above are additionally equipped with three constitutive laws, where isotropy of the electric conducticity \gamma, electric permittivity \epsilon and magnetic permeability \mu, is assumed \begin{equation} \mathbf{J} = \gamma (\mathbf{E} + \mathbf{v} \times \mathbf{B}), \label{eq:J} \end{equation} \begin{equation} \mathbf{D} = \epsilon \mathbf{E}, \end{equation} \begin{equation} \mathbf{B} = \mu \mathbf{H}. \label{eq:BH} \end{equation}

In the general case, the magnetic permeability \mu respectively reluctivity \nu=\frac{1}{\mu} is not constant for all values of \mathbf{B}, as it depends on the flux density \nu=\nu(|\mathbf{B}|).

The electric field intesity can then be split into a purely solenoidal part \mathbf E_s and an impressed one \mathbf E_{\textrm i}, which leads to the following splitting of the total current density (\ref{eq:J})

\begin{equation} \mathbf{J} = \underbrace{\gamma \mathbf{E}_{\textrm i}}_{\mathbf J_{\textrm i}} + \gamma (\mathbf{E}_{\textrm s} + \mathbf{v} \times \mathbf{B}), \end{equation}

in which \gamma \mathbf{E}_{\textrm i} can be considered as an impressed (prescribed) current density \mathbf J_{\textrm i}.

An important fact, which follows from the law of Ampére and (\ref{eq:J}), is that the total current is divergence-free \textbf{div} (\mathbf J) = 0.

Currently, only the low frequency regime of the electromagnetic spectrum is implemented in openCFS, which means that the time derivative of the displacement current density is neglected and therefore wave propagation is not part of the consideration.

PDE Types

Depending on if the weak form is discretized with nodal (classical Lagrange) elements or with Nédélec edge elements, we have different PDE's implemented:

  • Discretization with nodal elements:
<pdeList>
    <magnetic>
        ...
    </magnetic>
    ...
</pdeList>
  • Discretization with Nédélec edge elements:
<pdeList>
    <magneticEdge>
        ...
    </magneticEdge>
    ...
</pdeList>

However for most purposes, the edge element version is way more versatile but only 3D problems can be considered. If a 2D example has to be considered, an extraction of the 2D surface into the depth direction with one element is recommended, resulting in a quasi-2D problem.

Boundary conditions and loads

  • Essential (Dirichlet) boundary conditions on \Gamma_{\textrm D}, where \mathbf A \times \mathbf n and \mathbf A' \times \mathbf n is set to a zero. If these values are set to zero, we are in accordance with the considered function space V because it has zero tangential trace on the boundary. The fact that \mathbf A \times \mathbf n = 0 and $\mathbf A' \times \mathbf n = $ also results in
\begin{equation} \label{eq:binncurlazero} \mathbf B \cdot \mathbf n = \textbf{curl} \mathbf A \cdot \mathbf n = 0 \end{equation}

Boundary conditions

In openCFS this type of BC can be specified by

<bcsAndLoads>
    <fluxParallel name="Surface_1"/>
    <fluxParallel name="Surface_2"/>
</bcsAndLoads>
  • Natural (Neumann) boundary conditions on \Gamma_{\textrm N}, where \nu \textbf{curl} \mathbf A \times \mathbf n is set to a certain value, corresponding to a surface current density.

not implemented yet

  • Surface impedance boundary conditions (SIBC), where a linear relation between \mathbf E \times \mathbf n and \mathbf H \times \mathbf n is prescribed at the boundary, which can be used for surface excitation. This method presents an elegant way to drastically decrease the number of degrees of freedom in coil and inductor regions but it is only valid if the skin depth tends to zero, respectively if the frequency tends towards infinity.

partially implemented

Loads

In openCFS you can define these loads (list not completed):

  • Magnetic excitation
<fluxDensity name="V_magnet">
    <comp dof="z" value="1">
</fluxDensity> 

Coils

Coils are used in openCFS to model current-carrying conductors which produce a magnetic field. Here is a simple example how to use it as a xml-input, be aware that they must be specified after <bcsAndLoads> in <coilList>. It is recommended to prescibe the current density as current value (and use a wireCrossSection of 1). E.g. for a coil with a coilarea of "A", "n" windings and a wire current of "I" specify the current as shown:

<coilList>
    <coil id = " coil1 "> <! -- reference id for coil results -- >
        <source type = " current " value = "I*n/A" />
        <part id = " 1 "> <! -- may have several parts -- >
            <regionList>
                <region name = "coilwire" />
            </regionList>
            <direction>
                <analytic>
                    <comp dof = "z" value = " 1.0 " />
                </analytic>
            </direction>
            <wireCrossSection area = "1" />
        </part>
    </coil>
</coilList>

Potentials and Formulations

According to the solenoidal (divergence-free/no magnetic monopoles) \mathbf B-field (\ref{eq:solenoidalB}), there exists a vector potential \mathbf A, called the magnetic vector potential, defined by

\begin{equation} \label{eq:bequalcurlA} \mathbf B = \textbf{curl}( \mathbf A). \end{equation}

Plugging that into (\ref{eq:Faraday}), we obtain

\begin{equation} \label{eq:curlEpA} \textbf{curl} \left( \mathbf E + \frac{\partial \mathbf A}{\partial t} \right) = 0. \end{equation}

Since the kernel of the curl operator consists of all gradients of scalar functions, the gradient of a potential, called the electric scalar potential V can be added to the electric field intensity without changing (\ref{eq:curlEpA})

\begin{equation} \label{eq:Eprimegauge} \mathbf E = -\frac{\partial \mathbf A}{\partial t} -\textbf{grad} V. \end{equation}

The PDE-formulation containing both \mathbf A and V can be accessed via

<magneticEdge formulation="A-V">

Building on that, there are several ways how to fix the gauges, based on additional (constraint) equations. There is a large manifold of possible gauges, whereas only the most important ones for computational electromagnetics are mentioned here.

Coulomb gauge

This gauge states that

\begin{equation} \textbf{div} \mathbf A = 0. \end{equation}

The huge benefit of this gauge is, that if we are only interested in the magnetic flux denstity \mathbf B, there is no need to additionally gauge the formulation because of (\ref{eq:bequalcurlA}). If however the magnetic vector potential itself is observed, gauging techniques must be considered. For the magnetostatic case, we add a small ficticous conductivity (relaxation parameter), which does not effect the solution dramatically but we can ensure uniqueness of \mathbf A. If there is no definition of the relaxation parameter, the default value is 1\cdot 10^{-6}.

Temporal-/Weyl-gauge -- Classical A Formulation

The more important gauge, regarding computational electromagnetics is the temporal- or Weyl-gauge, e.g. see p.48 in 1, which sets the electric scalar potential to zero V=0. Doing this, results in the classical \mathbf A formulation, defined by the following PDE-formulation

<magneticEdge formulation="A">

Using the Weyl-gauge, results in the following description of the electric field intensity

\begin{equation} \label{eq:intermediate1} \mathbf E = -\frac{\partial \mathbf A}{\partial t}. \end{equation}

Combining (\ref{eq:intermediate1}), (\ref{eq:Ampere}), (\ref{eq:Faraday}) and the constitutive law (\ref{eq:J}) results in the well known curl-curl equation

\begin{equation} \label{eq:curl1} \textbf{curl} \left( \nu \textbf{curl}( \mathbf{A} ) \right) + \gamma \frac{\partial \mathbf{A}}{\partial t} - \gamma(\mathbf v \times \textbf{curl}( \mathbf A) ) = \mathbf{J}_{\textrm i}. \end{equation}

Weak Form

The considered function space for the classical A formulation is the curl-conforming space, defined by

\begin{equation} \label{eq:defVspace} V:=\lbrace \mathbf A \in H(\textbf{curl},\Omega): \textrm{tr}_\tau (\mathbf A)=0\, \textrm{on}\, \partial \Omega \rbrace. \end{equation}

By multiplying the strong form (\ref{eq:curl1}) with all test functions \mathbf A'\in V and integrating over the domain \Omega, we obtain

\begin{equation} \int_\Omega \left( \textbf{curl} \left( \nu \textbf{curl} \mathbf{A} \right) + \gamma \frac{\partial \mathbf{A}}{\partial t} \right) \cdot \mathbf A' = \int_\Omega \mathbf{J}_{\textrm i} \cdot \mathbf A', \quad \forall \mathbf A' \in V \end{equation}

and using the integration by parts for the first term, results in the final weak form

\begin{equation} \label{eq:weakcurl1} \int_\Omega \left( \nu \textbf{curl} \mathbf A \cdot \textbf{curl} \mathbf A' + \gamma \frac{\partial \mathbf A }{\partial t} \cdot \mathbf A' \right) = \int_{\partial \Omega} (\nu \textbf{curl} \mathbf A \times \mathbf n) \cdot \mathbf A' + \int_\Omega \mathbf{J}_{\textrm i} \cdot \mathbf A', \quad \forall \mathbf A' \in V, \end{equation}

where the first term on the right hand side defines the boundary conditions. By identifying the expression in the brackets as

\begin{equation} \nu \textbf{curl} \mathbf A \times \mathbf n = \mathbf H \times \mathbf n, \end{equation}

different boundary conditions can be stated.

Material

These material parameters are specified in the material-xml file, e.g. for copper

<material name="copper">
    <magnetic>
        <electricConductivity>
        <linear>
            <isotropic>5e7</isotropic>
        </linear>
        </electricConductivity>
        <magneticPermeability>
        <linear>
            <isotropic> 1.256637061E-06 </isotropic> 
        </linear>
        </magneticPermeability>
    </magnetic>
</material>

An example of a field dependent magnetic permeanbility is defined as follows

  <material name="iron">
    <magnetic>
      <electricConductivity>
        <linear>
          <isotropic>
            <real>1.0E5</real>
          </isotropic>
        </linear>
      </electricConductivity>
      <permeability>
        <linear>
          <isotropic>
            <real>  9.17290e-4 </real>
          </isotropic>
        </linear>
        <nonlinear>
          <isotropic>
            <dependency></dependency>
            <approxType>smoothSplines</approxType>
            <measAccuracy>   0.05      </measAccuracy>
            <maxApproxVal>   3.0       </maxApproxVal>
            <dataName>bh1.fnc</dataName>
          </isotropic>
        </nonlinear>
      </permeability>
    </magnetic>
  </material>

Analysis Types

Depending on the temporal setting, the formulation (\ref{eq:curl1}) can be used for

  • Static: \partial / \partial t (\cdot)= 0
<analysis>
  <static/>
</static>
  • Transient: \partial / \partial t (\cdot) \neq 0
<analysis>
  <transient>
        <numSteps>100</numSteps>
        <deltaT>1e-3</deltaT>
  </transient>
</analysis>
  • Harmonic: \partial / \partial t (\cdot) = j \omega (\cdot)
<analysis>
    <harmonic>
        <numFreq>10</numFreq>
        <startFreq>10</startFreq>
        <stopFreq>3000</stopFreq>
        <sampling>linear</sampling>
    </harmonic>
</analysis>

Postprocessing Results

Based on the solution A, different postprocessing results can be defined in openCFS, described in the following.

As in every other PDE in openCFS, the results are diveded into element-/ node-/ surface-element-/ surface-region and region-results, whereas we have to distinguish between the nodal- and edge-version.

Results for magEdgePDE (edge element version)

Element Results

  • Magnetic Vector Potential A (primary solution)
<elemResult type="magPotential">
  • Time Derivative of Magnetic Vector Potential
    • transient case: \begin{equation} \frac{\partial \mathbf A}{\partial t} \end{equation}
    • harmonic case: \begin{equation} j \omega \mathbf A \end{equation}
<elemResult type="magPotentialD1">
  • Magnetic Flux Density B \begin{equation} \mathbf B = \textbf{curl} (\mathbf A) \end{equation}
<elemResult type="magFluxDensity">
  • Magnetic Field Strength H \begin{equation} \mathbf H = \nu \mathbf B = \nu \textbf{curl} (\mathbf A) \end{equation}
<elemResult type="magFieldIntensity">
  • Electric Field Intensity
    • transient case: \begin{equation} \mathbf E = \frac{\partial \mathbf A}{\partial t} \end{equation}
    • harmonic case: \begin{equation} \mathbf E = j \omega \mathbf A \end{equation}
<elemResult type="elecFieldIntensity">
  • Magnetic Eddy Current Density
    • transient case: \begin{equation} \mathbf J_{\textrm{eddy}}= \gamma \mathbf E = \gamma \frac{\partial \mathbf A}{\partial t} \end{equation}
    • harmonic case: \begin{equation} \mathbf J_{\textrm{eddy}}= \gamma \mathbf E = \gamma j \omega \mathbf A \end{equation}
<elemResult type="magEddyCurrentDensity">
  • Magnetic Eddy Power Density (Joule Loss Power Density)
    • transient case: \begin{equation} P_{\textrm{Joule}} = \mathbf J_{\textrm{eddy}} \cdot \mathbf E = \mathbf J_{\textrm{eddy}} \cdot \frac{\partial \mathbf A}{\partial t} =\gamma \frac{\partial \mathbf A}{\partial t} \cdot \frac{\partial \mathbf A}{\partial t} \end{equation}
    • harmonic case: Averaged over one period of the base excitation frequency \begin{equation} \tilde{P}_{\textrm{Joule}} = \frac{1}{2} \gamma \omega^2 \mathbf A \cdot \mathbf A^* = \frac{1}{2} \gamma \omega^2 |\mathbf A| \end{equation}
<elemResult type="magJouleLossPowerDensity">
  • Total Current Density \begin{equation} \mathbf J_{\textrm{total}} = \mathbf J_{\textrm{impressed}} + \mathbf J_{\textrm{eddy}} \end{equation}
<elemResult type="magTotalCurrentDensity">
  • Lorentz Force Density \begin{equation} \mathbf F_{\textrm{Lorentz}} = \mathbf J_{\textrm{total}} \times \mathbf B = \mathbf J_{\textrm{total}} \times \mathbf B \end{equation}
<elemResult type="magForceLorentzDensity">

Surface Element Results

  • Magnetic Force Maxwell Density on Surface
<surfElemResult type="magForceMaxwellDensity">

Surface Region Results

  • Magnetic Eddy Current on Surface \begin{equation} J = \int_{A} \mathbf J_{\textrm{eddy}} \cdot \mathbf n\, dS \end{equation}
<surfRegionResult type="magEddyCurrent">
  • Magnetic Flux on Surface \begin{equation} \Phi = \int_{A} \mathbf B \cdot \mathbf n\, dS \end{equation}
<surfRegionResult type="magFlux">
  • Magnetic Force (Virtual Work Principle) on Surface
<surfRegionResult type="magForceVWP">
  • Magnetic Force Maxwell on Surface
<surfRegionResult type="magForceMaxwell">

Region Results

  • Magnetic Energy Implemented as \begin{equation} E_{\textrm{mag}} = \frac{1}{2} \sum_{e=1}^{n_e} \int_{\Omega^e} \mathbf A^{\textrm{T}} \cdot \mathbf K^e \cdot \mathbf A, \end{equation} where n_e is the number of elements in the considered region, \mathbf A the element solution vector and \mathbf K^e the element stiffness matrix. This corresponds to \begin{equation} E_{\textrm{mag}} = \frac{1}{2} \int_{\Omega} \mathbf B \cdot \mathbf B, \end{equation}
<regionResult type="magEnergy">
  • Magnetic Eddy Power \begin{equation} E_{\textrm{eddy}} = \int_{\Omega} \mathbf J_{\textrm{eddy}} \cdot \mathbf E \end{equation}
<regionResult type="magEddyPower">
  • Joule Loss Power
    • transient case: The same as Magnetic Eddy Power \begin{equation} E_{\textrm{Joule}} = \int_{\Omega} P_{\textrm{Joule}} = \int_{\Omega} \mathbf J_{\textrm{eddy}} \cdot \mathbf E \end{equation}
    • harmonic case: Averaged over one period of the base excitation frequency \begin{equation} E_{\textrm{Joule}} = \int_{\Omega} \tilde{P}_{\textrm{Joule}} \end{equation}
<regionResult type="magJouleLossPower">
  • Lorentz Force \begin{equation} F_{\textrm{Lorentz}} = \int_{\Omega} \mathbf F_{\textrm{Lorentz}} \end{equation}
<regionResult type="magForceLorentz">

Results for magneticPDE (nodal element version)

Tutorials

Permanent Magnet

References


  1. Miklos Kuczmann. Potential Formulations in Magnetics: Applying the Finite Element Method. Laboratory of Electromagnetic Fields, Szechenyi Istvan University, Gyor, Hungary, 2009. URL: http://maxwell.sze.hu/docs/C4.pdf