FEMM Axisymmetric Interpolation Functions, Flux Density, and Element Matrices

David Meeker
08Oct2021

Introduction

The purpose of this note is to explain the somewhat unusual formulation used in FEMM's axisymmetric formulation. This is basically the method for first-order triangles considered in Section IV of Henrotte et al. The basic difference between this and other first-order methods is that this method interpolates radially in terms of \(r^2\), rather than in terms of just \(r\). The result is a formulation that is more physically consistent than many other competing methods, yielding an accurate solution close to \(r=0\) where other methods can break down.

Interpolation

Instead of the regular magnetic vector potential A, the nodal field values are in terms of nodal flux \(\phi\), where \(A\) and \(\phi\) are related by: \[ \phi = 2 \pi r A\] For ease of notation, define: \[s=r^2\] Now, consider a "unit triangle" that lies between the points \((p,q)= (0,0),(1,0),(0,1)\) as pictured in Figure 1.
 (image: https://www.femm.info/wiki/Images/files.xml?action=download&file=elm.png)
Figure 1: Unit triangle.

We can map the \((p,q)\) coordinates back to the physical \((s,z)\) coordinates usina standard linear interpolation scheme:
\[ \begin{eqnarray}
s & = & (s_1 - s_0) p + (s_2 - s_0) q + s_0 \\
z & = & (z_1 - z_0) p + (z_2 - z_0) q + z_0
\end{eqnarray}\] Nodal \(\phi\) is interpolated in a similar way as: \[ \phi = (\phi_1 - \phi_0) p + (\phi_2 - \phi_0) q + \phi_0\]
This relationship between element and physical coordinates can be inverted to yield:
\[ \left\{ \begin{array}{c} p \\ q \end{array} \right\}
= \frac{1}{d}
\left[ \begin{array}{cc} (z_2 - z_0) & -(s_2-s_0) \\ -(z_1-z_0) & (s_1-s_0) \end{array}\right]
\left\{ \begin{array}{c} s-s_0 \\ z-z_0 \end{array} \right\}
\] where \[ d = (s_1-s_0)(z_2 - z_0) - (z_1-z_0)(s_2-s_0)\]
We could write this more compactly as:
\[ \left\{ \begin{array}{c} p \\ q \end{array} \right\}
= \frac{1}{d}
\left[ \begin{array}{cc} b_1 & c_1 \\ b_2 & c_2 \end{array}\right]
\left\{ \begin{array}{c} s-s_0 \\ z-z_0 \end{array} \right\}
\] where \[ d = c_2 b_1 - b_2 c_1\] where the \(b's\) and \(c's\) are the corresponding entries in the previous matrices.

It is interesting to note that interpolating \(s\) rather than \(r\) produces an element with "warped" sides in the \(rz\) plane.

Flux Density

In terms of \(\phi\), the flux density can be written as: \[B_r = - \frac{1}{2 \pi r} \frac{\partial \phi}{\partial z} \;\;\;\; B_z = \frac{1}{2 \pi r} \frac{\partial \phi}{\partial r} \] To be easily integrated with the fashion in which we are interpolating the element, the derivatives can be re-written as in terms of \(s\) rather than \(r\) as: \[B_r = - \frac{1}{2 \pi \sqrt{s}} \frac{\partial \phi}{\partial z} \;\;\;\; B_z = \frac{1}{\pi } \frac{\partial \phi}{\partial s} \]
Writing out the derivatives of the interpolation function yields:
\[ \left\{ \begin{array}{c} \frac{\partial \phi}{\partial s} \\ \frac{\partial \phi}{\partial z} \end{array} \right\} =
\left\{ \begin{array}{c}
\frac{\partial \phi}{\partial q} \frac{\partial q}{\partial s} +\frac{\partial \phi}{\partial p} \frac{\partial r}{\partial s} \\
\frac{\partial \phi}{\partial q} \frac{\partial q}{\partial z} +\frac{\partial \phi}{\partial p} \frac{\partial r}{\partial z} \\
\end{array} \right\}
= \frac{1}{d}
\left[ \begin{array}{cc} b_1 & b_2 \\ c_1 & c_2 \end{array}\right]
\left\{ \begin{array}{c} \phi_1-\phi_0 \\ \phi_2-\phi_0 \end{array} \right\}
\] Substituting into the formulas for the derivatives yields:
\[ B_r = - \frac{1}{2 \pi d \sqrt{s}} \left( c_1 (\phi_1-\phi_0) + c_2 (\phi_1 - \phi_0) \right) \]
\[ B_rz= - \frac{1}{ \pi d } \left( b_1 (\phi_1-\phi_0) + b_2 (\phi_1 - \phi_0) \right) \]

Interpretation of the form of \(B\)

Now, we can observe that \(B_z\) is constant over the element, which is a reasonable behavior for a first-order element. The \(B_r\) flux becomes smaller as \(r\) increases, but the \(1/r\) scaling factor is consistent with a constant amount of radially directed flux being spread over a larger area as the radius increases. This form is similar to the standard 3-node triangle formulation for 2-D planar problems, in which flux is constant over each element.

Element Matrices

Using this interpolation, element matrices can then be derived. A Mathematica notebook that does the derivation, along with a PDF printout of the notebook, are linked below in the Files section. Defining the additional quantities:
\[ \begin{eqnarray} b_o & = & z_1 - z_2 \\
c_0 & = & s_2 - s_1 \end{eqnarray} \]
the element matrices are:
\[ M_z = \frac{2 \pi}{d \mu_z} \left[
\begin{array}{ccc}
b_0^2 r_0^2 & b_0 b_1 r_0 r_1 & b_0 b_2 r_0 r_2 \\
b_0 b_1 r_0 r_1 & b_1^2 r_1^2 & b_1 b_2 r_1 r_2 \\
b_0 b_2 r_0 r_2 & b_1 b_2 r_1 r_2 & b_2^2 r_2^2 \\
\end{array}
\right] \]
\[M_r = \frac{\pi}{2 d \mu_r s_{avg}}
\left[ \begin{array}{ccc}
c_0^2 r_0^2 & c_0 c_1 r_0 r_1 & c_0 c_2 r_0 r_2 \\
c_0 c_1 r_0 r_1 & c_1^2 r_1^2 & c_1 c_2 r_1 r_2 \\
c_0 c_2 r_0 r_2 & c_1 c_2 r_1 r_2 & c_2^2 r_2^2 \\
\end{array} \right]
\] where \[ s_{avg} = \frac{1}{2} \left( \frac{(s_0-s_1)(s_1-s_2)(s_2-s_0)}
{s_0 s_1 \log \frac{s_1}{s_0}+s_1 s_2 \log \frac{s_2}{s_1}+s_0 s_2 \log \frac{s_0}{s_2}} \right)
\] for a formulation where the nodal unknowns are the vector potential at the corners of the element.

For elements that are small compared to the radius, \(s_{avg} \) is closely approximated as \((s0+s_1+s_2)/3\). However, care must be taken to evaluate \(s_{avg}\) in the cases where any node is located at \(r=0\) or any of the nodes have the same exactly the same radius. The special cases are:
\[ \begin{eqnarray}
\left. s_{avg} \right|_{s_0 \rightarrow 0} & = & \frac{s_1-s_2}{2 \log \frac{s_1}{s_2}} \\
\left. s_{avg} \right|_{s_1 \rightarrow 0} & = & \frac{s_2-s_0}{2 \log \frac{s_2}{s_0}} \\
\left. s_{avg} \right|_{s_2 \rightarrow 0} & = & \frac{s_0-s_1}{2 \log \frac{s_0}{s_2}} \\
\left. s_{avg} \right|_{s_0 \rightarrow s_1 } & = & \frac{(s_1-s_2)^2}{2 \left( s_1-s_2 + s_2 \log \frac{s_2}{s_1} \right)} \\
\left. s_{avg} \right|_{s_1 \rightarrow s_2 } & = & \frac{(s_2-s_0)^2}{2 \left( s_2-s_0 + s_0 \log \frac{s_0}{s_2} \right)} \\
\left. s_{avg} \right|_{s_2 \rightarrow s_0 } & = & \frac{(s_0-s_1)^2}{2 \left( s_0-s_1 + s_1 \log \frac{s_1}{s_0} \right)} \\
\left. s_{avg} \right|_{s_0 \rightarrow 0, s_1 \rightarrow s_2} & = & \frac{s_2}{2} \\
\left. s_{avg} \right|_{s_1 \rightarrow 0, s_2 \rightarrow s_0} & = & \frac{s_0}{2} \\
\left. s_{avg} \right|_{s_2 \rightarrow 0, s_0 \rightarrow s_1} & = & \frac{s_1}{2}
\end{eqnarray} \]
For more than one node point on \(r=0\), the entire \(M_r\) matrix is zero.

References

F. Henrotte et al., "A new method for axisymmetric linear and nonlinear problems," IEEE Transactions on Magnetics, 9(2):1352-1355, March 1993. doi:10.1109/CEFC.1992.720582

Attachments
File Last modified Size
DeriveAxiElement.nb 2021-10-08 09:06 177Kb
DeriveAxiElement.pdf 2021-10-08 09:06 555Kb

Valid XHTML :: Valid CSS: :: Powered by WikkaWiki