Welcome to maelstrom’s documentation!

Relevant publications are [KP03], [DAG89], [Chu08], [DDKS12].

maelstrom.heat

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

maelstrom.stabilization

Stabilization techniques for PDEs with dominating convection. The classical article about SUPG is [BH82]; for an overview of methods, see [JK07], [JK08], and [BGS04].

maelstrom.stabilization.supg(mesh, convection, diffusion, element_degree)

For each cell, this function return the expression

..math:

\begin{align*}
\tau &= \frac{h}{2\|b\|}
\left(\frac{1}{\tanh Pe} - \frac{1}{Pe}\right)\\
& = \frac{h^2}{4\varepsilon} \frac{1}{Pe}
\left(\frac{1}{\tanh Pe} - \frac{1}{Pe}\right)
\end{align*}

with the element diameter in the direction of the convection vector \(b\) and the Péclet number \(Pe = \frac{\|b\| h}{2\varepsilon}\); see (3) in [JK08]. Note that \(\tau\) does not have a singularity for \(\|b\|=0\) since

..math:

\frac{1}{\tanh Pe} - \frac{1}{Pe} = \frac{1}{3}Pe + O(Pe^3)

for \(Pe\approx 0\). This Taylor expansion (with a few more terms) is made use of in the code.

maelstrom.stokes

Numerical solution schemes for the Stokes equation in cylindrical coordinates.

maelstrom.stokes.F(u, p, v, q, f, r, mu, my_dx)
maelstrom.stokes.stokes_solve(up_out, mu, u_bcs, p_bcs, f, my_dx=<Mock name='mock.dx' id='139668395699504'>)

maelstrom.navier_stokes

Numerical solution schemes for the Navier–Stokes equation in cylindrical coordinates,

\[\DeclareMathOperator{\div}{div}\]
\[\begin{split}\begin{align*} &\rho \left(\frac{du}{dt} + (u\cdot\nabla)u\right) = -\nabla p + \mu \left( \frac{1}{r} \div(r \nabla u) - e_r \frac{u_r}{r^2} \right) + f,\\ &\frac{1}{r} \div(r u) = 0, \end{align*}\end{split}\]

cf. https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations#Representations_in_3D. In the weak formulation, we consider integrals in pseudo 3D, resulting in a weighting with \(2\pi r\) of the equations. (The volume element is \(2\pi r\,\text{d}x\).)

The order of the variables is taken to be \((r, z, \theta)\). This makes sure that for planar domains, the \(x\)- and \(y\)-coordinates are interpreted as \(r\), \(z\).

Note that [LMW+12] contains a chapter with an extensive comparison of solution methods for the transient Navier–Stokes equations for many different problems. It is found that IPCS is the best (fastest, most accurate) solution method.

An overview of projection methods for incompressible flow can be found in [GMS06] and [Joh16].

class maelstrom.navier_stokes.IPCS(time_step_method='backward euler')

Bases: object

Incremental pressure correction scheme; for details see [GMS06].

order = {'pressure': 1, 'velocity': 1}
step(dt, u, p0, W, P, u_bcs, p_bcs, rho, mu, f, verbose=True, tol=1e-10, my_dx=<Mock name='mock.dx' id='139668395699504'>)
maelstrom.navier_stokes.compute_pressure(P, p0, mu, ui, u, my_dx, p_bcs=None, rotational_form=False, tol=1e-10, verbose=True)

Solve the pressure Poisson equation

\[\begin{split}\begin{align} -\frac{1}{r} \div(r \nabla (p_1-p_0)) = -\frac{1}{r} \div(r u),\\ \text{(with boundary conditions)}, \end{align}\end{split}\]

for \(\nabla p = u\).

The pressure correction is based on the update formula

\[\begin{split}\frac{\rho}{dt} (u_{n+1}-u^*) + \begin{pmatrix} \text{d}\phi/\text{d}r\\ \text{d}\phi/\text{d}z\\ \frac{1}{r} \text{d}\phi/\text{d}\theta \end{pmatrix} = 0\end{split}\]

with \(\phi = p_{n+1} - p^*\) and

\[ \frac{1}{r} \frac{\text{d}}{\text{d}r} (r u_r^{(n+1)}) + \frac{\text{d}}{\text{d}z} (u_z^{(n+1)}) + \frac{1}{r} \frac{\text{d}}{\text{d}\theta} (u_{\theta}^{(n+1)}) = 0\]

With the assumption that u does not change in the direction \(\theta\), one derives

\[\begin{split}- \frac{1}{r} \div(r \nabla \phi) = \frac{1}{r} \frac{\rho}{dt} \div(r (u_{n+1} - u^*))\\ - \frac{1}{r} \langle n, r \nabla \phi\rangle = \frac{1}{r} \frac{\rho}{dt} \langle n, r (u_{n+1} - u^*)\rangle\end{split}\]

In its weak form, this is

\[\int r \langle\nabla\phi, \nabla q\rangle \,2 \pi = - \frac{\rho}{dt} \int \div(r u^*) q \, 2 \pi - \frac{\rho}{dt} \int_{\Gamma} \langle n, r (u_{n+1}-u^*)\rangle q \, 2\pi.\]

(The terms \(1/r\) cancel with the volume elements \(2\pi r\).) If the Dirichlet boundary conditions are applied to both \(u^*\) and \(u_n\) (the latter in the velocity correction step), the boundary integral vanishes.

If no Dirichlet conditions are given (which is the default case), the system has no unique solution; one eigenvalue is 0. This however, does not hurt CG convergence if the system is consistent, cf. [vdV03]. And indeed it is consistent if and only if

\[\int_\Gamma r \langle n, u\rangle = 0.\]

This condition makes clear that for incompressible Navier-Stokes, one either needs to make sure that inflow and outflow always add up to 0, or one has to specify pressure boundary conditions.

Note that, when using a multigrid preconditioner as is done here, the coarse solver must be chosen such that it preserves the nullspace of the problem.

maelstrom.navier_stokes.compute_tentative_velocity(time_step_method, rho, mu, u, p0, dt, u_bcs, f, W, my_dx, tol)

Compute the tentative velocity via

\[\rho (u_0 + (u\cdot\nabla)u) = \mu \frac{1}{r} \div(r \nabla u) + \rho g.\]
maelstrom.navier_stokes.compute_velocity_correction(ui, p0, p1, u_bcs, rho, mu, dt, rotational_form, my_dx, tol, verbose)

Compute the velocity correction according to

\[U = u_0 - \frac{dt}{\rho} \nabla (p_1-p_0).\]

References

[BGS04]Pavel B. Bochev, Max D. Gunzburger, and John N. Shadid. Stability of the supg finite element method for transient advection–diffusion problems. Computer Methods in Applied Mechanics and Engineering, 193(23-26):2301–2323, June 2004. URL: http://dx.doi.org/10.1016/j.cma.2004.01.026, doi:10.1016/j.cma.2004.01.026.
[BH82]Alexander N. Brooks and Thomas J.R. Hughes. Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations. Computer Methods in Applied Mechanics and Engineering, 32(1-3):199–259, September 1982. URL: http://dx.doi.org/10.1016/0045-7825(82)90071-8, doi:10.1016/0045-7825(82)90071-8.
[CCG+97]C. Chaboudez, S. Clain, R. Glardon, D. Mari, J. Rappaz, and M. Swierkosz. Numerical modeling in induction heating for axisymmetric geometries. IEEE Transactions on Magnetics, 33(1):739–745, 1997. URL: http://dx.doi.org/10.1109/20.560107, doi:10.1109/20.560107.
[Chu08]Paul K. Chu. Semiconductor Materials & Physics. City University of Hong Kong, 2008. URL: http://personal.cityu.edu.hk/~appkchu/AP4120/2.PDF.
[DAG89]J.J. Derby, L.J. Atherton, and P.M. Gresho. An integrated process model for the growth of oxide crystals by the czochralski method. Journal of Crystal Growth, 97(3-4):792–826, October 1989. URL: http://dx.doi.org/10.1016/0022-0248(89)90583-6, doi:10.1016/0022-0248(89)90583-6.
[DDKS12]W. Dreyer, P. -É. Druet, O. Klein, and J. Sprekels. Mathematical modeling of czochralski type growth processes for semiconductor bulk single crystals. Milan Journal of Mathematics, 80(2):311–332, July 2012. URL: http://dx.doi.org/10.1007/s00032-012-0184-9, doi:10.1007/s00032-012-0184-9.
[GRS07]Christian Grossmann, Hans-Görg Roos, and Martin Stynes. Numerical Treatment of Partial Differential Equations. Springer Berlin Heidelberg, 2007. URL: http://dx.doi.org/10.1007/978-3-540-71584-9, doi:10.1007/978-3-540-71584-9.
[GMS06]J.L. Guermond, P. Minev, and Jie Shen. An overview of projection methods for incompressible flows. Computer Methods in Applied Mechanics and Engineering, 195(44-47):6011–6045, September 2006. URL: http://dx.doi.org/10.1016/j.cma.2005.10.010, doi:10.1016/j.cma.2005.10.010.
[Joh16]Volker John. Finite Element Methods for Incompressible Flow Problems. Springer International Publishing, 2016. URL: http://dx.doi.org/10.1007/978-3-319-45750-5, doi:10.1007/978-3-319-45750-5.
[JK07]Volker John and Petr Knobloch. On spurious oscillations at layers diminishing (sold) methods for convection–diffusion equations: part i – a review. Computer Methods in Applied Mechanics and Engineering, 196(17-20):2197–2215, March 2007. URL: http://dx.doi.org/10.1016/j.cma.2006.11.013, doi:10.1016/j.cma.2006.11.013.
[JK08]Volker John and Petr Knobloch. On spurious oscillations at layers diminishing (sold) methods for convection–diffusion equations: part ii – analysis for and finite elements. Computer Methods in Applied Mechanics and Engineering, 197(21-24):1997–2014, April 2008. URL: http://dx.doi.org/10.1016/j.cma.2007.12.019, doi:10.1016/j.cma.2007.12.019.
[KP02]O. Klein and P. Philip. Correct voltage distribution for axisymmetric sinusoidal modeling of induction heating with prescribed current, voltage, or power. IEEE Transactions on Magnetics, 38(3):1519–1523, May 2002. URL: http://dx.doi.org/10.1109/20.999125, doi:10.1109/20.999125.
[KP03]Olaf Klein and Peter Philip. Transient numerical investigation of induction heating during sublimation growth of silicon carbide single crystals. Journal of Crystal Growth, 247(1-2):219–235, January 2003. URL: http://dx.doi.org/10.1016/s0022-0248(02)01903-6, doi:10.1016/s0022-0248(02)01903-6.
[KL12]Michael Kolmbauer and Ulrich Langer. A robust preconditioned minres solver for distributed time-periodic eddy current optimal control problems. SIAM Journal on Scientific Computing, 34(6):B785–B809, January 2012. URL: http://dx.doi.org/10.1137/110842533, doi:10.1137/110842533.
[LMW+12]Anders Logg, Kent-Andre Mardal, Garth N. Wells, and others. Automated Solution of Differential Equations by the Finite Element Method. Springer, 2012. ISBN 978-3-642-23098-1. doi:10.1007/978-3-642-23099-8.
[vdV03]Henk A. van der Vorst. Iterative Krylov Methods for Large Linear Systems. Cambridge University Press, 2003. URL: http://dx.doi.org/10.1017/cbo9780511615115, doi:10.1017/cbo9780511615115.

Indices and tables