Hi folks! After a couple of serious and not-so-serious lifestyle-esque posts, it’s time to get back to work. Mine, to be exact. If you’ve read the previous instalments in this series, you know that I’m currently working on a faster alternative to time-stepping FEM. If you haven’t, I suggest you start from the beginning. Today, it’s the time to discuss some implementation issues.


In the previous post, I described the basics of the harmonic balance method. The time-dependent vector potential \avec(t) is written as a truncated Fourier series

\avec(t) \approx \sum\limits_{k=1}^{2P+1} \avec_k \timefunction_k(t),

where \timefunction_k (t) is a convenience notation for either \cos or \sin. This sum approximation is then substituted into the governing equation

\stiffnessmatrix(t) \avec(t) = \loadvector(t),

and the unknown coefficients \avec_k are then solved with the Galerkin’s method.

This gives us a big block matrix equation, with the matrix blocks \stiffnessmatrix_{rc} equal to

\stiffnessmatrix_{rc} = \langle \timefunction_r(t) \stiffnessmatrix(t) \timefunction_c(t) \rangle,

where the angle brackets \langle \rangle correspond to calculating the time-average. From this matrix system, the coefficients \avec_k can then be solved.


However, this leaves us with the not-so-trivial problem of actually computing these blocks. As I described earlier, the stiffness matrix \stiffnessmatrix(t) can be divided into constant and time-dependent parts, i.e.

\stiffnessmatrix(t) = \stiffnessmatrix_0 + \Sag(t).

Obviously, the static component \stiffnessmatrix_0 is no problem to handle. However, this cannot be said for the time-dependent component \Sag.

Indeed, the most typical approach to model the rotation of the machine is to change the finite element mesh in the air-gap, where the stator and rotor are moving past each other. In this approach (and most others, for that matter), the entries of \Sag(t) will experience sharp, discontinuous changes as the elements of the mesh change. Thus, calculating the time-integral

\langle \timefunction_r(t) \stiffnessmatrix(t) \timefunction_c(t) \rangle

accurately could be quite complicated indeed.

However, and luckily, this is a problem I have already solved. Remember earlier when I spoke about calculating the Fourier series of the matrix \Sag(t)? In other words, I already know how to write \Sag(t) as

\Sag(t) \approx \sum\limits_{k=1}^{N} \Sag_k \timefunction_k(t),

and how to compute the coefficients \Sag_k.

Number of required terms

Of course, having to compute a huge number  N of the terms would not be very practical, would it. I mean, that would be required to present the highly discontinuous entries of \Sag accurately, right?

Well, yes and no. Accurate representation of \Sag would indeed need quite a many terms. However, we are not interested in an accurate series approximation of \Sag – we are interested in calculating the matrix blocks

\langle \timefunction_r(t) \stiffnessmatrix(t) \timefunction_c(t) \rangle!

Thus, let’s now substitute the approximation for \Sag(t) here, to get

\langle \timefunction_r(t) \stiffnessmatrix(t) \timefunction_c(t) \rangle = \langle \timefunction_r(t)  \timefunction_c(t) \rangle \stiffnessmatrix_0 + \sum\limits_{k=1}^N \langle \timefunction_r(t) \timefunction_k(t) \timefunction_c(t) \rangle \Sag_k.

In other words, we now longer have a time-average of a matrix – we have the time average of the product of three trigonometric functions \langle \timefunction_r \timefunction_k \timefunction_c \rangle. And this can be expressed in closed form quite easily.

And even more better, there are only a limited number of terms k for which the expression is non-zero. And these values of k can also be determined analytically. Thus, it is usually sufficient to compute something like N = 4P of the air-gap matrix Fourier coefficients.

And this is obviously a much easier task to computing a total of (2P+1)^2 of the matrix block integrals, like we would have been forced to do before.


Need help with electric motor design or design software? Let's get in touch - satisfaction guaranteed!
Making time-stepping FEM faster, part 4

3 thoughts on “Making time-stepping FEM faster, part 4

Leave a Reply

Your email address will not be published. Required fields are marked *