maelstrom.maxwell

\[\DeclareMathOperator{\div}{div} \DeclareMathOperator{\curl}{curl}\]

The equation system defined here are largely based on [CCG+97], with the addition of the material flux \(u\) of the liquid.

Given an electric field \(E\), a magnetic field \(B\), and a material flux \(u\), the current density \(J\) in the material is given by

\[J = \sigma (E + u \times B).\]

where \(\sigma\) is the electrical conductivity.

This leads to

\[\curl(\sigma^{-1} J - u\times B + \text{i} \omega A) = 0.\]

Assuming that \(B\) is given by the potential \(\phi\), \(B = \curl(\phi e_{\theta})\), we end up with

\[\begin{split}u \times B &= u \times \curl(\phi e_{\theta}) \\ &= u \times \left( - \frac{\text{d}\phi}{\text{d}z} e_r + \frac{1}{r} \frac{\text{d}(r\phi)}{\text{d}r} e_z \right) \\ &= - u_z \frac{\text{d}\phi}{\text{d}z} e_{\theta} + u_{\theta} \frac{\text{d}\phi}{\text{d}z} e_z - u_r \frac{1}{r} \frac{\text{d}(r\phi)}{\text{d}r} e_{\theta} + u_{\theta} \frac{1}{r} \frac{\text{d}(r\phi)}{\text{d}r} e_r.\end{split}\]

(Note that \(\curl\) is taken in cylindrial coordinates.)

Following Chaboudez, this eventually leads to the equation system

\[\begin{split}\begin{cases} - \div\left(\frac{1}{\mu r} \nabla(r\phi)\right) - \sigma u_\theta \div\left(\frac{1}{r}\nabla(r\phi)\right) + \sigma \left\langle (u_r, u_z)^T, \frac{1}{r}\nabla(r\phi) \right\rangle + \text{i} \sigma \omega \phi = \frac{\sigma v_k}{2\pi r} \quad\text{in } \Omega,\\ n\cdot\left(- \frac{1}{\mu r} \nabla(r\phi)\right) = 0 \quad\text{on }\Gamma \setminus \{r=0\}\\ \phi = 0 \quad\text{ for } r=0. \end{cases}\end{split}\]

The differential operators are interpreted like 2D for \(r\) and \(z\).

For the weak formulation, the volume elements \(2\pi r\,\text{d}x\) are used. This corresponds to the full 3D rotational formulation and also makes sure that at least the diffusive term is nice and symmetric. Additionally, it avoids dividing by \(r\) in the convections and the right hand side.

Here with no convection in direction \(\theta\):

\[\int_\Omega \div\left(\frac{1}{\mu r} \nabla(r u)\right) (2\pi r v) + \langle b, \nabla(r u)\rangle 2\pi v + \text{i} \sigma \omega u 2 \pi r v = \int_\Omega \sigma v_k v.\]

maelstrom.maxwell.build_system(V, dx, Mu, Sigma, omega, f_list, f_degree, convections, bcs)

Build FEM system for

\[\div\left(\frac{1}{\mu r} \nabla(r\phi)\right) + \left\langle u, \frac{1}{r} \nabla(r\phi)\right\rangle + \text{i} \sigma \omega \phi = f\]

by multiplying with \(2\pi r v\) and integrating over the domain and the preconditioner given by [KL12].

maelstrom.maxwell.compute_joule(Phi, voltages, omega, Sigma, Mu, subdomain_indices)

See, e.g., equation (2.17) in [CCG+97].

In a time-harmonic approximation with

..math:

\begin{align}
A &= \Re(a exp(\text{i} \omega t)),\\
B &= \Re(b exp(\text{i} \omega t)),
\end{align}

the time-average of \(A\cdot B\) over one period is

..math:

\overline{A\cdot B} = \frac{1}{2} \Re(a \cdot b^*)

see http://www.ece.rutgers.edu/~orfanidi/ewa/ch01.pdf. In particular,

..math:

\overline{A\cdot A} = \frac{1}{2} \|a\|^2.

Consequently, we can compute the average source term over one period as

..math:

s = \frac{1}{2} \|j\|^2 / \sigma = \frac{1}{2} \|E\|^2 \sigma.

(Not using \(j\) avoids explicitly dividing by \(\sigma\) which is 0 at nonconductors.)

maelstrom.maxwell.compute_lorentz(Phi, omega, sigma)

In a time-harmonic discretization with quantities

\[\begin{split}\begin{align} A &= \Re(a \exp(\text{i} \omega t)),\\ B &= \Re(b \exp(\text{i} \omega t)), \end{align}\end{split}\]

the time-average of \(A\times B\) over one period is

\[\overline{A\times B} = \frac{1}{2} \Re(a \times b^*),\]

see http://www.ece.rutgers.edu/~orfanidi/ewa/ch01.pdf. Since the Lorentz force generated by the current \(J\) in the magnetic field \(B\) is

\[F_L = J \times B,\]

its time average is

\[\overline{F_L} = \frac{1}{2} \Re(j \times b^*).\]

With

\[\begin{split}J &= \Re(\exp(\text{i} \omega t) j e_{\theta}),\\ B &= \Re\left( \exp(i \omega t) \left( -\frac{\text{d}\phi}{\text{d}z} e_r + \frac{1}{r} \frac{\text{d}(r\phi)}{\text{d}r} e_z \right) \right),\end{split}\]

we have

\[\begin{split}\overline{F_L} &= \frac{1}{2} \Re\left(j \frac{d\phi^*}{dz} e_z + \frac{j}{r} \frac{d(r\phi^*)}{dr} e_r\right)\\ &= \frac{1}{2} \Re\left(\frac{j}{r} \nabla(r\phi^*)\right)\\\end{split}\]

In the workpiece, we can assume

\[j = -\text{i} \sigma \omega \phi\]

which gives

\[\begin{split}\begin{align*} \overline{F_L} &= \frac{\sigma\omega}{2r} \Im\left( \phi \nabla(r \phi^*) \right)\\ &= \frac{\sigma\omega}{2r} \left( \Im(\phi) \nabla(r \Re(\phi)) -\Re(\phi) \nabla(r \Im(\phi)) \right) \end{align*}\end{split}\]

maelstrom.maxwell.compute_potential(coils, V, dx, mu, sigma, omega, convections, verbose=True, io_submesh=None)

Compute the magnetic potential \(\Phi\) with \(A = \exp(\text{i} \omega t) \Phi e_{\theta}\) for a number of coils.

maelstrom.maxwell.get_voltage_current_matrix(phi, physical_indices, dx, Sigma, omega, v_ref)

Compute the matrix that relates the voltages with the currents in the coil rings. (The relationship is indeed linear.)

This is according to [KP02].

The entry \(J_{k,l}\) in the resulting matrix is the contribution of the potential generated by coil \(l\) to the current in coil \(k\).

maelstrom.maxwell.prescribe_voltage(A, b, coil_rings, voltage, v_ref, J)

Get the voltage coefficients \(c_l\) with the total voltage prescribed.

maelstrom.maxwell.solve(V, dx, Mu, Sigma, omega, f_list, convections, f_degree=None, bcs=None, tol=1e-12, verbose=False)

Solve the complex-valued time-harmonic Maxwell system in 2D cylindrical coordinates

\[\div\left(\frac{1}{\mu r} \nabla(r\phi)\right) + \left\langle u, \frac{1}{r} \nabla(r\phi)\right\rangle + i \sigma \omega \phi = f\]

with finite elements.

Parameters:

• V – function space for potentials
• dx – measure
• Mu (dictionary) – mu per subdomain
• Sigma (dictionary) – sigma per subdomain
• omega (float) – current frequency
• f_list – list of right-hand sides for each of which a solution will be computed
• convections (dictionary) – convection, defined per subdomain
• bcs – Dirichlet boundary conditions
• tol (float) – relative solver tolerance
• verbose (boolean) – solver verbosity

Return type:

list of functions