The Finite Volume Method an Alternative Method for LF Electromagnetic Problems

To give an idea as well as more details on the FVM method, in this section we will study the FVM in both cases 2D with triangular element and 3D with prism volume. For both cases, starting from Maxwell’s equations [19], in terms of the magnetic vector potential A, the equation to solve is:

$\nabla \times \frac{1}{\mu} \nabla \times \mathbf{A}=\mathbf{J}_{\mathrm{s}}+\mathbf{J}_{i n d}$  (1)

$\mathbf{A}=0$ at infinity  (2)

μ is the magnetic permeability, J_s is the source current density and J_ind is the eddy current density. Eq. (2) describes the impressed Dirichlet condition at infinity.

2.1 2D case

In the 2D case the magnetic flux density is defined only in the Oxy plan and consequently the magnetic vector potential and the source current density have only one component as:

$\mathbf{A}=A_{z} \mathbf{k}$ (3)

$\mathbf{J}_{s}=J_{s z} \mathbf{k}$ (4)

where, k is the unit vector in the z direction. Thus, in this case Eq. (1) can be rewritten as:

$-\nabla \cdot\left(\frac{1}{\mu} \nabla A_{z}\right)=J_{s z} \mathbf{k}+\mathbf{J}_{i n d}$  (5)

The eddy current density, which is induced either by movement or temporal variation of the current source, is given by:

$\mathbf{J}_{\text {ind}}=\sigma\left(-\frac{\partial A_{z}}{\partial t} \mathbf{k}+\mathbf{v} \times \mathbf{B}\right)$  (6)

In Eq. (6) v is the velocity, B is the magnetic flux density and $\sigma$ is the electrical conductivity. In the 2D case; the velocity term can be rewritten as follows:

$\mathbf{v} \times \mathbf{B}=-\left(v_{x} \cdot \frac{\partial A_{z}}{\partial x}+v_{y} \cdot \frac{\partial A_{z}}{\partial y}\right) \mathbf{k}$  (7)

2.jpg

Figure 2. Triangular mesh scheme for 2D FVM method

v_x and v_y are the velocity components in the Oxy plan. By substituting Eq. (7) in Eq. (6), and by using the harmonic representation, Eq. (5) becomes:

$-\nabla \cdot\left(\frac{1}{\mu} \nabla A_{z}\right)+\sigma\left(j \omega A_{z}+v_{x} \cdot \frac{\partial A_{z}}{\partial x}+v_{y} \cdot \frac{\partial A_{z}}{\partial y}\right) \mathbf{k}=J_{s z} \mathbf{k}$  (8)

At this stage, Eq. (8) has to be solved by the 2D FVM method with triangular mesh. As described in Figure 2, in the case of triangular mesh scheme, each principal control volume Ω_p is represented by its principal node P and surrounded by three neighbourhood control volumes represented by their nodes K, L and M.

Subsequently, the FVM method involves a surface integration of Eq. (8) over the principal control volume Ω_p of node P as:

$\iint_{\Omega_{p}}-\nabla \cdot\left(\frac{1}{\mu} \nabla A_{z}\right) \mathrm{d} \Omega+\iint_{\Omega_{p}} \sigma \omega A_{z} \mathrm{jd} \Omega$

$+\iint_{\Omega_{p}} \sigma\left(v_{x} \cdot \frac{\partial A_{z}}{\partial x}+v_{y} \cdot \frac{\partial A_{z}}{\partial y}\right) \mathrm{d} \Omega=\iint_{\Omega_{p}} J_{s z} \mathrm{d} \Omega$  (9)

The first term in Eq. (9) can be rewritten as follows:

$\begin{align}  & \iint\limits_{{{\Omega }_{p}}}{-\nabla .\left( \frac{1}{\mu }\nabla {{A}_{z}} \right)}d\Omega =-\frac{1}{\mu }\oint\limits_{\partial {{\Omega }_{p}}}{\nabla {{A}_{z}}.\mathbf{n}.\mathbf{dl}} \\ & \begin{matrix}   \begin{matrix}   \begin{matrix}   {} & {}  \\\end{matrix} & {}  \\\end{matrix} & {} & {} & {}  \\\end{matrix}=-\frac{1}{\mu }\sum\limits_{i=1}^{3}{\nabla {{A}_{z}}.{{\mathbf{n}}_{i}}.{{d}_{ei}}} \\\end{align}$   (10)

n is the normal vector and d_ei is the i^th edge length with respect to i=1, 2, 3. The term $\nabla A_{z} \cdot \boldsymbol{n}_{i}$ can be approximated for each edge by a finite difference, for example for the edge d_e1:

$\left(\nabla A_{z} \cdot \mathbf{n}_{1}\right)_{d e 1}=\frac{A_{z}^{K}-A_{z}^{P}}{d_{K P}}$  (11)

By the way, Eq. (10) leads to an algebraic expression as:

$-\frac{1}{\mu} \sum_{i=1}^{3} \nabla A_{z} \mathbf{n}_{i} \cdot d_{e i}=-\frac{1}{\mu}\left[\frac{A_{z}^{K}-A_{z}^{P}}{d_{K P}} \cdot d_{e 1}\right.$$\left.+\frac{A_{z}^{L}-A_{z}^{P}}{d_{L P}} \cdot d_{e 2}+\frac{A_{z}^{M}-A_{z}^{P}}{d_{M P}} \cdot d_{e 3}\right]$  (12)

Eq. (12) can be rewritten as follows:

$-\frac{1}{\mu} \sum_{i=1}^{3} \nabla A_{z} \mathbf{n}_{i} \cdot d_{e i}=\frac{1}{\mu}\left(\frac{d_{e 1}}{d_{K P}}+\frac{d_{e 2}}{d_{L P}}+\frac{d_{e 3}}{d_{M P}}\right) \cdot A_{z}^{P}$$-\frac{1}{\mu}\left(\frac{d_{e 1}}{d_{K P}} \cdot A_{z}^{K}+\frac{d_{e 2}}{d_{L P}} \cdot A_{z}^{L}+\frac{d_{e 3}}{d_{M P}} \cdot A_{z}^{M}\right)$ (13)

Let’s assume a uniform distribution of the current density in the principal control volume Ω_p. Therefore, the source term is integrated as follows:

$\iint_{\Omega p} J_{sz} d \Omega=J_{szP} \cdot S_{\Omega p}$  (14)

The integration of all terms given by Eq. (9) leads in an algebraic equation that express the potential in the principal control volume in terms of potentials of its three neighbourhoods as:

$A_{z}^{P}=\frac{1}{C_{p}}\left(C_{\mathrm{k}} A_{z}^{K}+C_{1} A_{z}^{L}+C_{\mathrm{m}} A_{z}^{M}+J_{\mathrm{szP}} S_{\Omega p}\right)$  (15)

The coefficients C_p, C_k, C_l and C_m describe geometrical and physical properties associated with the principal control volume and its neighbourhoods. J_szP is the current density in the principal control volume. Rewriting Eq. (15) for all control volumes, thus leads to the following matrix equation:

$\left[\mathbf{M}_{2}\right]\left[A_{z}\right]=\left[J_{s z}\right]$  (16)

Compared to the FEM method, we point out that the FVM method leads to a sparser matrix M₂; this particular property are very significant regards to computational time and required storage memory.

2.2 3D case

The 3D standard FVM method featured by hexahedral control volume shape has been already demonstrated in previous works showing its effectiveness for electromagnetic devices. Unfortunately, the hexahedral element appears too cost in the case of complex geometry, such as for example for circular coils. To overcome the limitation of the standard FVM method and enable it to handle more complex geometry, we propose a 3D unstructured FVM method. In such alternative the simple hexahedral control volume is replaced by a prism control volume (Figure 3).

3.jpg

Figure 3. Prism control volume for 3D FVM method

Here, each control volume element D_p consists of a central node P and delimited by five faces. The neighbourhood control volumes of D_p are represented by their central nodes K, L and M in the Oxy plan, B and T according to the z direction. As in the 2D case, Eq. (1) must be integrated over the prism control volume:

$\int_{z}\left[\iint_{\Omega_{p}}\left(\nabla \times \frac{1}{\mu} \nabla \times \mathbf{A}\right) d \Omega_{p}\right] d z=\int_{z}\left[\iint_{\Omega_{p}}\left(\mathbf{J}_{\mathrm{s}}+\mathbf{J}_{i n d}\right) d \Omega_{p}\right] d z$  (17)

with:

$d \Omega_{p} d z=d D_{p}$  (18)

According to the Green-Ostrogradsky theorem, Eq. (17) becomes:

$-\iiint_{D_{p}}\left(\nabla \cdot\left(\frac{1}{\mu} \nabla \mathbf{A}\right)\right) \mathrm{d} D_{p}=\iiint_{D_{p}}\left(\mathbf{J}_{\mathrm{s}}+\mathbf{J}_{i n d}\right) \mathrm{d} D_{p}$  (19)

For the first term we have:

$-\iiint\limits_{Dp}{\nabla .(\frac{1}{\mu }\nabla A)d}{{D}_{p}}=-\mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}\limits_S \frac{1}{\mu } \nabla A.d\mathbf{s}$  (20)

Let consider for example the integration of the x component:

$-\mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}\limits_s \frac{1}{\mu } \nabla {{A}_{x}}.ds.\mathbf{n}=-\sum\limits_{i=1}^{5}{\frac{1}{{{\mu }_{i}}}\nabla {{A}_{x}}.d{{\mathbf{s}}_{i}}}$

$=-\left[\frac{S_{1}}{\mu_{1}} \frac{A_{x 1}-A_{x p}}{d l_{1} \sin \left(d e_{1}, d l_{1}\right)}+\frac{S_{2}}{\mu_{2}} \frac{A_{x 2}-A_{x p}}{d l_{2} \sin \left(d e_{2}, d l_{2}\right)}\right.$

$+\frac{S_{3}}{\mu_{3}} \frac{A_{x 3}-A_{x p}}{d l_{3} \sin \left(d e_{3}, d l_{3}\right)}+\frac{S_{4}}{\mu_{4}} \frac{A_{x 4}-A_{x p}}{d l_{4} \sin \left(d e_{4}, d l_{4}\right)}$

$\left.+\frac{S_{5}}{\mu_{5}} \frac{A_{x 5}-A_{x p}}{d l_{5} \sin \left(d e_{5}, d l_{5}\right)}\right]$  (21)

Let us define the first coefficient as:

$C_{1}=\frac{S_{1}}{\mu_{1}} \frac{1}{d l_{1} \sin \left(d e_{1}, d l_{1}\right)}$  (22)

Subsequently, all the coefficients in the previous algebraic equation can be rewritten as the the following general form:

$C_{i}=\frac{S_{i}}{\mu_{i}} \frac{1}{d l_{i} \sin \left(d e_{i}, d l_{i}\right)}$  (23)

The developpement of the x component of the first term leads to the following algebraic equation:

$-\left(C_{1}\left(A_{x 1}-A_{x p}\right)+C_{2}\left(A_{x 2}-A_{x p}\right)+C_{3}\left(A_{x 3}-A_{x p}\right)\right.$

$\left.+C_{4}\left(A_{x 4}-A_{x p}\right)+C_{5}\left(A_{x 5}-A_{x p}\right)\right)=\left[-\left(C_{1} A_{x 1}+C_{2} A_{x 2}\right.\right.$

$\left.+C_{3} A_{x 3}+C_{4} A_{x 4}+C_{5} A_{x 5}\right)$

$\left.+A_{x p}\left(C_{1}+C_{2}+C_{3}+C_{4}+C_{5}\right)\right]$       (24)

By assuming that the current density is constant within the control volume, the x component of the potential in the central node P is given by:

$A_{x p}=\frac{1}{C_{p}}\left(\sum_{i=1}^{5} C_{i} A_{x i}-C_{6} \frac{V_{x 2}-V_{x 1}}{d x}+C_{x s}\right)$  (25)

C_p, C₁, C₂, C₃, C₄, C₅and C₆, describe geometrical and physical properties associated with the principal control volume and its five neighbourhoods. V_xi is the electrical potential. Similarly, we can find both y and z components of the potential in the central node as:

$A_{y p}=\frac{1}{C_{p}}\left(\sum_{i=1}^{5} C_{i} A_{y i}-C_{6} \frac{V_{y 2}-V_{y 1}}{d y}+C_{y s}\right)$   (26)

$A_{z p}=\frac{1}{C_{p}}\left(\sum_{i=1}^{5} C_{i} A_{z i}-C_{6} \frac{V_{z 2}-V_{z 1}}{d z}+C_{z s}\right)$   (27)

These last three equations can be rewritten in a matrix form as:

$\left[\mathbf{M}_{3}\right][\mathbf{A}]=\left[\mathbf{J}_{s}\right]$  (28)

As in the 2D case, the 3D FVM method leads to sparse and symmetric matrix M₃; this feature is strongly significant regards to computational time and required storage memory especially for the 3D case. By the way, both 2D and 3D computer codes are developed in Matlab. To check the developed computer codes, two TEAM workshop problems and an experimental electromagnetic micro-actuator are analyzed in the next sections.

This paper presents the FVM as an alternative method for low frequency electromagnetic problems and in particular those require eddy current and electromagnetic force or torque computation. Both 2D and 3D computer codes are developed and tested. To examine the efficiency of both programs and showing the competitiveness of the FVM, both TEAM workshop problems No.28 and No.30 and an experimental electromagnetic micro-actuator are analyzed. One can see that numerical results obtained using the FVM are compared to the analytical or experimental results provided by these problems. Furthermore, the 3D FVM associated to the prism mesh can be a very promising for modelling problems without being constrained by their geometries.