MagHeatPDE

This work description illustrates all the possible combinations of electromagnetic-thermal coupled (induction-heating) problems.

The emphasis is put on the combinations of nonlinearities, e.g., nonlinear BH-curve for the electromagnetic problem in frequency domain, coupled with temperature-dependent thermal material parameters in a steady state problem.

Since there are differences in the implemented features between 2D and 3D setups, they are also included in this description.

Induction-Heating Introduction

Induction-heating simulations always involves two different PDEs:

• the electromagnetic problem, briefly and not exhaustively described here MagneticPDE,
• the heat conduction (thermal) problem, briefly and not exhaustively described here HeatConductionPDE, which can be coupled in a lot of different ways regarding their temporal domain (transient magnetics - steady state heat, harmonic magnetics - transient heat, ...).

The reason for this vast amount of coupling possibilities are the timescale differences between the magnetic- and the thermal field (qualitatively shown in the Figure down below). In the magnetic problem, the typical order of magnitude for the timescales are ms (milliseconds) - corresponds to an excitations in the kHz frequency regime), whereas the thermal problem has timeconstants in the s-region, so several orders of magnitude between them. If both PDEs are coupled in time-domain (transient electromagnetic - transient heat conduction), the simulation needs thousands of timesteps of the magnetics problem, without a significant change in the thermal field, which is far from efficient.

In the following, we will discuss the different coupling possibilities together with the corresponding nonlinearities, depicted in the Figure below.

Coupling Magnetic-and Heat PDE

Depending on the kind of simulation (list can be found below), the coupling between the MagneticPDE and the HeatConductionPDE can be established via the Joule losses (Eddy current losses) in time domain (\ref{eq:13}) or as usually in frequency domain (\ref{eq:17}).

To give an example where to use which coupling combination, the following simulations will be considered.

• "$\mu s$ heating"
• "classical" transient heating
• steady state heat solution

For $\mu s$ heating, where the heating process is in the time range of the electromagnetic problem, time varying Joule losses are necessary to couple both PDE's with each other. Depending on the simulation requirements, this kind of induction heating can be performed with the combinations found in 1a, 2a, 3a or 4a in the lists below.

In the "classical" transient heating where the heating object is excited with a time varying magnetic field in the kHz range and the heating process takes up to several seconds, it makes no sense to couple both PDE's in time domain as mentioned above. Hence, the coupling between the two PDE's will be established via the periodic averaged joule losses in frequency domain. This means that the magnetic PDE is solved in frequency and the heat PDE in time domain. In the list below, the different coupling combination for this kind of simulation can be found. At this point it should be mentioned that not all combinations of nonlinearity and temporal coupling are implemented, e.g., there is no option to couple nonlinear magnetics in frequency-domain (multiharmonic) with temperature dependent BH-curves.

Warning

Not all combinations are possible. It depends on the selected non-linearities and temporal couplings.

To obtain the steady state heat solution immediately, the term $\frac{\partial T}{\partial t}$ in the heat PDE will be neglected. The coupling between the PDE's is again done via the periodically averaged joule losses in frequency domain.

List of possible combinations for 3D setups

The list of possible combinations, which will also be discussed in detail in the next sections, is as follows:

1. Coupling: Pure forward coupling: no adaption of $B$ or $\sigma$ according to temperature ($\sigma = \textrm{const.}$, $B=B(H)$)
   1. Temporal combination: transient magnetics - transient heat
      1. linear magnetics - linear heat
      2. nonlinear magnetics ($B=B(H)$)) - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
      3. linear magnetics - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
      4. nonlinear magnetics ($B=B(H)$)) - linear heat
   2. Temporal combination: magnetics in frequency-domain - transient heat
      1. linear magnetics - linear heat
      2. linear magnetics - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
      3. nonlinear magnetics (multiharmonic) - linear heat
      4. nonlinear magnetics (multiharmonic) - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
   3. Temporal combination: magnetics in frequency-domain - steady state heat
      1. linear magnetics - linear heat
      2. linear magnetics - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
      3. nonlinear magnetics (multiharmonic) - linear heat
      4. nonlinear magnetics (multiharmonic) - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
2. Coupling: Forward coupling with adapted electric conductivity $\sigma = \sigma(T)$
   1. Temporal combination: transient magnetics - transient heat
      1. linear magnetics - linear heat
      2. nonlinear magnetics ($B=B(H)$)) - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
      3. linear magnetics - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
      4. nonlinear magnetics ($B=B(H)$)) - linear heat
3. Coupling: Forward coupling with temperature-dependent BH-curves $B = B(H,T)$ and $\sigma=\textrm{const.}$
   1. Temporal combination: transient magnetics - transient heat
      1. nonlinear magnetics ($B=B(H,T)$)) - linear heat
      2. nonlinear magnetics ($B=B(H,T)$)) - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
4. Coupling: Forward coupling with temperature-dependent BH-curves $B = B(H,T)$ and $\sigma=\sigma (T)$
   1. Temporal combination: transient magnetics - transient heat
      1. nonlinear magnetics ($B=B(H,T)$)) - linear heat
      2. nonlinear magnetics ($B=B(H,T)$)) - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)

List of possible combinations for 2D setups

The list of possible combinations, which will also be discussed in detail in the next sections, is as follows:

1. Coupling: Pure forward coupling: no adaption of $B$ or $\sigma$ according to temperature ($\sigma = \textrm{const.}$, $B=B(H)$)
   1. Temporal combination: transient magnetics - transient heat
      1. linear magnetics - linear heat
      2. nonlinear magnetics ($B=B(H)$)) - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
      3. linear magnetics - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
      4. nonlinear magnetics ($B=B(H)$)) - linear heat
   2. Temporal combination: magnetics in frequency-domain - transient heat
      1. linear magnetics - linear heat
      2. linear magnetics - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
   3. Temporal combination: magnetics in frequency-domain - steady state heat
      1. linear magnetics - linear heat
      2. linear magnetics - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
2. Coupling: Iterative coupling with adapted electric conductivity $\sigma = \sigma(T)$
   1. Temporal combination: transient magnetics - transient heat
      1. linear magnetics - linear heat
      2. nonlinear magnetics ($B=B(H)$)) - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
      3. linear magnetics - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
      4. nonlinear magnetics ($B=B(H)$)) - linear heat
3. Coupling: Iterative coupling with temperature-dependent BH-curves $B = B(H,T)$ and $\sigma=\textrm{const.}$
   1. Temporal combination: transient magnetics - transient heat
      1. nonlinear magnetics ($B=B(H,T)$)) - linear heat
      2. nonlinear magnetics ($B=B(H,T)$)) - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)
4. Coupling: Iterative coupling with temperature-dependent BH-curves $B = B(H,T)$ and $\sigma=\sigma (T)$
   1. Temporal combination: transient magnetics - transient heat
      1. nonlinear magnetics ($B=B(H,T)$)) - linear heat
      2. nonlinear magnetics ($B=B(H,T)$)) - nonlinear heat ($\rho(T), c_{\textrm p}(T), \lambda (T)$)

Eddy current and hysteresis losses in ferromagnetic materials

To decompose the eddy current and hysteresis losses in ferromagnetic materials, the pointing's vector can be used.

$$\begin{equation} \begin{split} \label{eq:7} P^{av}=\frac{1}{T} \int_{T} \frac{dW}{dt} dt = \frac{1}{T} \int_{T} \int_{\Omega} \nabla \cdot \textbf{S} d \Omega ' dt \\ = \frac{1}{T} \int_{T} \int_{\Omega} \nabla \cdot (\textbf{E} \times \textbf{H}) d \Omega ' dt \end{split} \end{equation}$$

By exploiting the two Maxwell equations (eddy current case)

$$\begin{equation} \label{eq:8} \nabla \times \textbf{H} = \textbf{J} = \textbf{J}_i + \textbf{J}_e + \gamma \underbrace{(\textbf{v}\times\textbf{B})}_{\text{assume}\,\textbf{v} = 0} \end{equation}$$

$$\begin{equation} \label{eq:9} \nabla \times \textbf{E}_s = -\frac{\partial\textbf{B}}{\partial t} \end{equation}$$

$\textbf{J}_e$ refers to the eddy current density

as well as the constitutive law and a vector identity,

$$\begin{equation} \label{eq:10} \textbf{J}_e = \gamma \textbf{E}_s; \,\,\, \nabla \cdot ( \textbf{E}_s \times \textbf{H})=\textbf{H}\cdot\nabla\times\textbf{E}_s -\textbf{E}_s \cdot \nabla \times \textbf{H} \end{equation}$$

one ends in the following relation:

$$\begin{equation} \begin{split} \label{eq:11} P^{av} = \frac{1}{T} \int_{\Omega} \int_{T} \textbf{H}\cdot \frac{\partial \textbf{B}}{\partial t} dt d \Omega ' + \frac{1}{T} \int_{\Omega} \int_{T} \textbf{E}_s \cdot (\textbf{J}_i + \textbf{J}_e)dt d\Omega ' = \\ \frac{1}{T} \int_{\Omega} \int_{B_1(t)}^{B_2(t)} \textbf{H}\cdot d \textbf{B} d \Omega ' + \frac{1}{T} \int_{\Omega} \int_{T} \textbf{E}_s\cdot \textbf{J} dt d\Omega ' = P_{hyst}^{av} + P_{e}^{av} \end{split} \end{equation}$$

Warning

Hysteresis losses are not taken into account in the magnetic heat coupling.

For the induction heating, hysteresis losses will be neglected, as it is assumed that only soft magnetic materials will be utilized. Soft magnetic materials are characterized by a narrow hysteresis curve. Compared to the eddy current losses, the hysteresis losses are order of magnitudes lower. For the harmonic balancing method, the hysteresis curve is not valid. As only soft magnetic materials are assumed to be used, the hysteresis loop can be idealized as a single curve (anhyteretic curve). The induced error by idealizing the curve is negligible for induction heating.

Joule losses in time domain

The heat sources in the HeatConductionPDE, are the eddy current losses in the control volume.

$$\begin{equation} \label{eq:12} P(t) = P_{diss}(t) = \int_{\Omega} \dot{q}_d d\Omega = P_e(t) = \int_{\Omega} \textbf{E}_s \cdot \textbf{J} d\Omega ' \end{equation}$$

$$\begin{equation} \label{eq:13} \dot{q}_d = \textbf{E}_s \cdot (\textbf{J}_i + \textbf{J}_e) = \frac{d\textbf{A}}{dt} \cdot \textbf{J}_i + \gamma \frac{d\textbf{A}}{dt}\cdot \frac{d \textbf{A}}{dt} \end{equation}$$

The first term in (\ref{eq:13}) represents the heat sources inside the conductor due to the impressed current density and the second term represent the heat sources due to the eddy currents in the heating object.

For the purposes of induction heating, the focus will be on the heat distribution in the workpiece and thus, the joule losses in the conductor will not be further discussed. The joule losses in the time domain are then finally modeled with the following equation.

$$\begin{equation} \label{eq:18} \dot{q}_d = \gamma \frac{d\textbf{A}}{dt}\cdot \frac{d \textbf{A}}{dt} \end{equation}$$

Joule losses in frequency domain

For classical induction heating, where the heating process takes seconds and the steady state solution of the magnetic field is reached in $\mu s$, joule losses in time domain are not an efficient method to couple the PDE's with each other. Thus, the periodic averaged joule losses $\bar{q}_d$ in frequency domain will be considered in the heat PDE.

To transform the eddy current losses into frequency domain, the multi-harmonic ansatz will be needed.

$$\begin{equation} \label{eq:14} \gamma \frac{d\textbf{A}}{dt}\cdot \frac{d \textbf{A}}{dt} \rightarrow \gamma \left( \sum_{k=-N}^{N} j k \omega \hat{\textbf{A}}_k\cdot e^{jk\omega t} \cdot \sum_{l=-N}^{N} j l \omega \hat{\textbf{A}}_l\cdot e^{jk\omega t}\right) \end{equation}$$

In a next step, period averaging is performed, to result in the total averaged power losses.

$$\begin{equation} \begin{split} \bar{q}_{d} = \frac{1}{T}\int_{T} \dot{q}_{d,1}(t)dt = \\ \frac{\gamma}{T}\int_{T} \frac{d\textbf{A}}{dt}\cdot \frac{d \textbf{A}}{dt} dt = & \,\omega \gamma \frac{j^2}{T}\int_{T} \left( \sum_{k=-N}^{N} \sum_{l=-N}^{N} k l \hat{\textbf{A}}_k\cdot \hat{\textbf{A}}_l\cdot e^{j(k+l)\omega t}\right) \end{split} \label{eq:15} \end{equation}$$

By exploiting the property (\ref{eq:2}) of an integrated harmonic e-function over a period, equation (\ref{eq:15}) simplifies to

$$\begin{equation} \label{eq:16} \bar{q}_{d} = \omega \gamma \frac{j^2}{T} \frac{T}{2} \left( \sum_{k=-N}^{N} -k^2 \hat{\textbf{A}}_k\cdot \hat{\textbf{A}}_{-k}\right) \end{equation}$$

Since the solution of $\textbf{A}(\textbf{x},t)$ is periodic in time, the Fourier coefficients have to be conjugate-symmetric $\hat{\textbf{A}}_{-k} = \bar{\hat{\textbf{A}}}_k$, where the bar over the vector potential denotes the conjugate complex. The final expression for the periodic averaged joule losses reads as:

$$\begin{equation} \label{eq:17} \bar{q}_{d} = \omega \gamma \frac{1}{2}\left( \sum_{k=-N}^{N} k^2 |\hat{\textbf{A}}_k|^2\right) \end{equation}$$

For the linear case, N = 0.

Time harmonic FEM

In the case where only the steady-state solution of the eddy current problem is needed, the (multi)-harmonic Ansatz can be used to solve for the magnetic vector potential $\textbf{A}$. A solution can be found by applying the following steps:

• transform the magnetic problem in the frequency domain
• solve for the harmonic(s) of the magnetic vector potential
• transform it back to get a periodic steady state solution of the magnetic vector potential

The coupling between the heat and magnetic PDE is not significantly impacted by the use of the steady state solution, since the transient process will be lost in the steady state solution due to the different timescales in the magnetic and heat PDE, as previously mentioned.

The magnetic PDE in $\textbf{A}$ formulation can be derived from the magneto-quasistatic subset of the maxwells equations. A more detailed description of the deviation can be found here, MagneticPDE.

$$\begin{equation} \label{eq:eddycurrent} \textbf{curl} \left( \nu (||\textbf{curl}( \mathbf{A}||_2)\textbf{curl}( \mathbf{A}) ) \right) + \gamma \frac{\partial \mathbf{A}}{\partial t} = \mathbf{J}_{\textrm i}. \end{equation}$$

The motional emf term will be neglected for all further considerations.

Linear Case

In this subsection, the reluctivity $\nu$ will be considered as constant. This leads to some simplifications of the previous stated eddy current equation (\ref{eq:eddycurrent}). With a harmonic Ansatz the equation reads as:

$$\begin{equation} \label{eq:lineareddycurrent} \textbf{curl} \left( \nu \textbf{curl}( \hat{\mathbf{A}}) ) \right) + \gamma j\omega \hat{\mathbf{A}} = \hat{\mathbf{J}}{\textrm i}. \end{equation}$$

The big advantage of the linear case is that it can be solved for different frequencies independently, which is not the case for nonlinear materials where couplings between the different frequencies occur.

The weak formulation of (\ref{eq:lineareddycurrent}) is given by

$$\begin{equation} \label{eq:weakform_linear} \int_{\Omega} \nu \textbf{curl}(\hat{\mathbf{A}}) \cdot \textbf{curl}(\hat{\mathbf{A}'}) + j\omega \gamma \hat{\mathbf{A}} \cdot \hat{\mathbf{A}'} d\Omega = \int_{\Omega} \hat{\mathbf{J}'} \cdot \hat{\mathbf{A}'} d\Omega, \hspace{1cm} \forall \hat{\mathbf{A}'} \in V \end{equation}$$

where $\hat{\mathbf{A}'}$ denotes the applied test function and V the considered function space (more details can be found here, MagneticPDE). To discretize the system, one can perform an edge finite element discretization, which approximate the vector potential in the following way:

$$\begin{equation} \label{eq:edge_discretization} \hat{\mathbf{A}} \approx \sum_{n=1}^{m_{eq}} \textbf{N}_n \Psi_n \end{equation}$$

$\Psi_n$ denotes the corresponding degree of freedom, namely the line integral of the magnetic vector potential along the n-th edge, and $\textbf{N}_n$ the edge shape function associated with the n-th edge. Inserting this into the weak form (\ref{eq:weakform_linear}) results in:

$$\begin{equation} \label{eq:discretized_weakform} \sum_{a = 1}^{m_{eq}} \sum_{b = 1}^{m_{eq}} \Psi_b \left( \Psi_a\int_{\Omega} \nu \textbf{curl}(\textbf{N}_a) \cdot \textbf{curl}(\textbf{N}_b) d\Omega + j \omega \Psi_a \int_{\Omega} \gamma \textbf{N}_a \cdot \textbf{N}_b d\Omega - \int_{\Omega} \textbf{N}_b \cdot \hat{\textbf{J}_i} d\Omega \right) = 0 \end{equation}$$

Since the solution of $\Psi_b$ is not from interest, equation (\ref{eq:discretized_weakform}) can be divided by $\Psi_b$. In compact matrix notation the final system of equation has the following form:

$$\begin{equation} \label{eq:discretized_weakform_matrix_notation} (\textbf{K}(\nu) + j\omega \textbf{M}\cdot) \underline{\Psi} = \hat{\underline{J}}_i \end{equation}$$

The vector $\underline{\Psi}$ contains all degrees of freedom and the matrices $\textbf{K}$ and $\textbf{M}$ are the assembled mass and stiffness matrix from each element.

Nonlinear Case

Since the reluctivity $\nu(B) = \nu(||\textbf{curl}( \mathbf{A}||_2)) $ is a solution dependent quantity in the general case, it is not sufficient to use only one harmonic (fourier transform for one frequency). This can only be done for linear materials where $\nu(B) = \nu = constant$. Hence, a fourier series expansion for the magnetic vector potential and the solution dependent reluctivity is required. This expansion can be written in the following manner:

$$\begin{equation} \label{eq:fourmagvec} \textbf{A}((\textbf{x},t) = \Re \left( \sum_{k=0}^{N} \hat{\textbf{A}}_k(\textbf{x})\cdot e^{jk\omega t}\right) = \sum_{k=-N}^{N} \hat{\textbf{A}}_k(\textbf{x})\cdot e^{jk\omega t} \end{equation}$$

$$\begin{equation} \label{eq:fourrel} \nu(||\textbf{curl}( \mathbf{A}||_2)) = \sum_{m=-M}^{M} \hat{\nu}_m(\textbf{x})\cdot e^{jm\omega t}, \hspace{1cm} M \leq N \end{equation}$$

The boundary N refers to the maximum number of harmonics which are used to extend the vector potential $\textbf{A}$ and M is the maximum number of harmonics which will be considered for the reluctivity $\nu(t)$.

To construct the nonlinear system in the frequency domain, the fourier expansions (\ref{eq:fourmagvec}) and (\ref{eq:fourrel}) have to be inserted into the PDE for the eddy current case (\ref{eq:eddycurrent}).

$$\begin{equation} \begin{split} \textbf{curl} \sum_{m=-M}^{M} \hat{\nu}_m(\textbf{x})\cdot e^{jm\omega t} \textbf{curl} \sum_{k=-N}^{N} \hat{\textbf{A}}_k(\textbf{x})\cdot e^{jk\omega t} = \\ \hat{\textbf{J}} e^{j \omega t} - j \omega \gamma \sum_{k=-N}^{N} k \hat{\textbf{A}}_k(\textbf{x})\cdot e^{jk\omega t} \end{split} \label{eq:fourinsertPDE} \end{equation}$$

Multiplying the equation (\ref{eq:fourinsertPDE}), by $e^{−iωnt}$ for n $\in$ [−N, ..., N ] and integrating over one period $T = 2\pi/\omega$, results in

$$\begin{equation} \label{eq:1} \begin{split} \int_{0}^{\frac{2 \pi}{\omega}} \textbf{curl} \sum_{m=-M}^{M}\sum_{k=-N}^{N} \hat{\nu}_m(\textbf{x})\cdot e^{j(m+k-n)\omega t} \textbf{curl} \sum_{k=-N}^{N} \hat{\textbf{A}}_k(\textbf{x})\cdot e^{j(k-n)\omega t} dt = \\ \int_{0}^{\frac{2 \pi}{\omega}} \hat{\textbf{J}} e^{j(1-n) \omega t} - j \omega \gamma \sum_{k=-N}^{N} k \hat{\textbf{A}}_k(\textbf{x})\cdot e^{j(k-n)\omega t} dt \end{split} \end{equation}$$

By exploiting the property

$$\begin{equation} \label{eq:2} \int_{T} e^{j(k-n)\omega t} dt = \delta_{kn}=\begin{cases} \frac{2 \pi}{\omega} = T & \text{ if } k = n\\ \;\;0 & \text{ if } k \neq n \end{cases} \end{equation}$$

$$\begin{equation} \label{eq:3} \int_{T} e^{j(k+m-n)\omega t} dt = \delta_{(k+m),n}=\begin{cases} \frac{2 \pi}{\omega} =T & \text{ if } k+m = n\\ \;\;0 & \text{ if } k+m \neq n \end{cases} \end{equation}$$

one ends up in the strong formulation of the eddy current case in frequency domain:

$$\begin{equation} \begin{split} \textbf{curl} \left( \begin{array}{cccccc} \hat{\nu}_0 & \hat{\nu}_{-1} & \cdots & \hat{\nu}_{-M}& & \\ \hat{\nu}_1 & \hat{\nu}_0 & \hat{\nu}_{-1} & & \ddots & \\ \vdots & & \ddots & & & \hat{\nu}_{-M} \\ \hat{\nu}_M & & & & & \vdots \\ & \ddots & & & \hat{\nu}_{0} & \hat{\nu}_{-1}\\ & & \hat{\nu}_M & \cdots & \hat{\nu}_{1} & \hat{\nu}_{0} \\ \end{array}\right) \left( \begin{array}{cccccc} \textbf{curl} \hat{\textbf{A}}_{-N} \\ \vdots\\ \textbf{curl} \hat{\textbf{A}}_0\\ \textbf{curl} \hat{\textbf{A}}_1\\ {\vdots}\\ \textbf{curl} \hat{\textbf{A}}_N \\ \end{array}\right) \\ + j\omega \gamma \left( \begin{array}{cccccc} -N & & & & & \\ & \ddots & & & & \\ & & 0 & & & \\ & & & 1 & & \\ & & & & \ddots & \\ & & & & & N \\ \end{array}\right) \left( \begin{array}{cccccc} \hat{\textbf{A}}_{-N} \\ \vdots\\ \hat{\textbf{A}}_0\\ \hat{\textbf{A}}_1\\ {\vdots}\\ \hat{\textbf{A}}_N \\ \end{array}\right) = \left( \begin{array}{cccccc} \hat{\textbf{J}}_{-N} \\ \vdots\\ \hat{\textbf{J}}_0\\ \hat{\textbf{J}}_1\\ {\vdots}\\ \hat{\textbf{J}}_N \\ \end{array}\right) \end{split} \label{eq:xxx} \end{equation}$$

In the next step, one has to derive the weak form of this system. This can be done in a very similar way described in MagneticPDE. For the discretization of the vector potential $\hat{\textbf{A}}$, the same Ansatz as in the linear case can be used.

The discretized weak form of system (\ref{eq:xxx}) is then the given by

$$\begin{equation} \begin{split} \label{eq:6} \left(\begin{bmatrix} \textbf{K}(\hat{\nu}_0) & \textbf{K}(\hat{\nu}_{-1}) & \cdots & \textbf{K}(\hat{\nu}_{-M})& & \\ \textbf{K}(\hat{\nu}_1) & \textbf{K}(\hat{\nu}_0) & \textbf{K}(\hat{\nu}_{-1}) & & \ddots & \\ \vdots & & \ddots & & & \textbf{K}(\hat{\nu}_{-M}) \\ \textbf{K}(\hat{\nu}_M) & & & & & \vdots \\ & \ddots & & & \textbf{K}(\hat{\nu}_{0}) & \textbf{K}(\hat{\nu}_{-1})\\ & & \textbf{K}(\hat{\nu}_M) & \cdots & \textbf{K}(\hat{\nu}_{1}) & \textbf{K}(\hat{\nu}_{0}) \\ \end{bmatrix}\right. \\ + j\omega \gamma \left.\begin{bmatrix} -N\textbf{M} & & & & & \\ & \ddots & & & & \\ & & 0 & & & \\ & & & \textbf{M} & & \\ & & & & \ddots & \\ & & & & & N\textbf{M} \\ \end{bmatrix}\right) \begin{bmatrix} \underline{\Psi}_{-N} \\ \vdots\\ \underline{\Psi}_{0}\\ \underline{\Psi}_{1}\\ {\vdots}\\ \underline{\Psi}_{N} \\ \end{bmatrix} = \begin{bmatrix} \hat{\underline{J}}_{-N} \\ \vdots\\ \hat{\underline{J}}_0\\ \hat{\underline{J}}_1\\ {\vdots}\\ \hat{\underline{J}}_N \\ \end{bmatrix} \end{split} \end{equation}$$

where the vector $\underline{\Psi}_{-N}$ contains all degree of freedom for the k-th harmonics. $\textbf{K}(\hat{\nu}_k)$ is the stiffness Matrix and $\textbf{M}$ the mass matrix, which do not change from iteration to iteration. These matrices are huge matrices, which contain all spatial degrees of freedom.

In contrast to the linear case (\ref{eq:discretized_weakform_matrix_notation}) where we can investigate different frequencies independent of each other (no couplings), the nonlinear case results in a more complex system, since the different harmonics are coupled with each other. Due to the nonlinearity of $\nu(\textbf{B})$, the system has to be solved repeatedly for each adaptation of $\hat{\nu}_k$ until the change in the solution disappears from one iteration to the other. With the help of the harmonic balancing method, this task can be solved. There are different ways to perform harmonic balancing, but for this example the approach consist of using an alternating time-frequency scheme depicted in figure below. The superscript $i$ in the following section denotes the iteration counter, where the subscript represent the harmonic number k.

Harmonic balancing: alternating time-frequency scheme

At first, the system will be solved only for the linear case $\nu(t) = \nu_0(t) = constant$. This results in the fourier coefficients $\underline{\Psi}_k^0$ of the magnetic vector potential. By performing an inverse fourier transform, the magnetic vector potential $\textbf{A}$ and furthermore the magnetic flux density $\textbf{B} = \textbf{curl}\textbf{A}$ are obtained. Now, the reluctivity $\nu(t)$ can be evaluated for the obtained magnetic field. The fourier transform of the time dependent reluctivity has now multiple harmonics $\hat{\nu}_k^1$ in contrast to the linear reluctivity which has only the linear part $\hat{\nu}_0^0$. With the multi-harmonic reluctivity the system can be solved again which results into new fourier coefficients $\underline{\Psi}_k^1$. The circle is then repeated $i$ times until no change in the vector potential $\textbf{A}$ is detected from one to the next iteration (convergence is achieved).

Tutorials

InductionHeating