Making time-stepping FEM faster, part 3

Hi you, and welcome the the third instalment in my blog series about my ongoing work! In the first post, I explained what time-stepping finite element analysis is and why it should be made faster in the first place. The second post basically discussed how the time-dependent FEM solution could be presented as a Fourier series – if only the coefficients could be calculated without having to generate a lot of samples in advance. So, enter the harmonic balance method.

The harmonic WHAT?

Sorry about the name. I didn’t come up with it.

Anyway, in the previous post I concluded that the time-dependent solution vector $\avec(t)$ to typical FEM problems could be written as a Fourier series

$\avec(t) \approx \avec_0 + \sum\limits_{k=1}^{P} \avec_{\mathrm{cos},k} \cos(\omega k t) + \avec_{\mathrm{sin},k} \sin(\omega k t)$ .

By renaming the terms a bit, this can be more compactly written as

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

where the time-function $\timefunction$ is either $\sin$ or $\cos$ , and $\avec_k$ are the corresponding Fourier coefficients.

No matter the notation, the number of terms $P$ could be quite reasonable, practically a small fraction of the number of samples produced by time-stepping analysis. Of course, the coefficients $\avec_k$ would first have to be calculated.

And this is where the harmonic balance method steps in. Originally developed for analysing non-linear circuits, the method consists of writing the solution a Fourier series (non-surprisingly), and then solving the unknown coefficients. For this purpose, the power of the method of weighted residuals is harnessed.

Galerkin’s method

…everybody calls it Galerkin’s method, though. This is technically correct, although sligthtly inaccurate, I guess.

Remember how the FEM problem was defined by a system of equations of form

$\stiffnessmatrix(t) \avec(t) = \loadvector(t)$ ?

The Galerkin’s method consists of first defining the residual – the difference between the two sides of the equation

$\resvector(t) = \stiffnessmatrix(t) \avec(t) - \loadvector(t)$ .

The inner product of the residual and some suitable weight function is then set to zero.

Say what? More simply put, in this case, the residual is multiplied with a suitable function of time, e.g. $\cos$ . This yields us

$\cos(\omega t) \resvector(t) = \cos(\omega t) \stiffnessmatrix(t) \avec(t) - \cos(\omega t) \loadvector(t)$ .

Since this expression is still a function of time, the time-average is still calculated, and set to equal zero. In other words (or mathematical terms, to be exact), we get

$\frac{1}{T} \int\limits_{0}^T { \cos(\omega t) \stiffnessmatrix(t) \avec(t) - \cos(\omega t) \loadvector(t) } \mathrm{d}t = \zeromatrix$ .

For this, the more concise $\langle \rangle$ is often adopted, so we can equivalently write

$\langle \cos(\omega t) \stiffnessmatrix(t) \avec(t) \rangle = \langle \cos(\omega t) \loadvector(t) \rangle$ .

This is often called the Galerkin’s orthogonality criterion.

Matrix system

The last expression alone is of course meaningless, since we don’t know the solution $\avec(t)$ . Thus, we finally substitute the Fourier series expression for $\avec(t)$ , and so get

$\langle \cos(\omega t) \stiffnessmatrix(t) \sum\limits_{k=1}^{2P+1} \avec_k \timefunction_k(t)\rangle = \langle \cos(\omega t) \loadvector(t) \rangle$

Then, we take the $\timefunction_k(t)$ functions of different frequencies that we have in the series, and in turn use each of them as a weight function in the orthogonality criterion. This, obviously, gives us as many equation-systems that we have unknown vectors $\avec_k$ . These can be expressed as a big block matrix system

$\begin{bmatrix} \stiffnessmatrix_{1,1} & \stiffnessmatrix_{1,2} & \ldots & \stiffnessmatrix_{1,2P+1} \\ \stiffnessmatrix_{2,1} & \stiffnessmatrix_{2,2} & \ldots & \stiffnessmatrix_{2,2P+1} \\ \vdots & \vdots & \ddots & \vdots \\ \stiffnessmatrix_{2P+1,1} & \stiffnessmatrix_{2P+1,2} & \ldots & \stiffnessmatrix_{2P+1,2P+1} \end{bmatrix} \begin{bmatrix} \avec_1 \\ \avec_2 \\ \vdots \\ \avec_{2P+1} \end{bmatrix} \begin{bmatrix} \end{bmatrix} \begin{bmatrix} \loadvector_1 \\ \loadvector_2 \\ \vdots \\ \loadvector_{2P+1} \end{bmatrix}$ ,

with the blocks corresponding to

$\stiffnessmatrix_{r,c} &= \langle \timefunction_r(t) \stiffnessmatrix(t) \timefunction_c(t) \rangle$
$\loadvector = \langle \timefunction_r(t) \loadvector(t) \rangle$ .

By solving this system, we get the coefficients $\avec_k$ . And when we have the coefficients, we have an approximation for the time-dependent behaviour $\avec(t)$ as

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

Results

And this finally brings us to the figure I have shown you two times already.

Harmonic balance simulation. — Time-stepping results (brute-force) compared to the harmonic balance method (reduced).

The blue curve was obtained with nothing else than the harmonic balance method, with the number of terms equal to 15.

And as can be seen, the results are quite good indeed.

But, as often, they could still be improved. And I’m working on precisely that. But before that, let’s have one post about some practicalities!

Check out EMDtool - Electric Motor Design toolbox for Matlab.

Need help with electric motor design or design software? Let's get in touch - satisfaction guaranteed!

Making time-stepping FEM faster, part 3

Antti Lehikoinen

Making time-stepping FEM faster, part 3

Galerkin’s method

Matrix system

Results

Leave a Reply Cancel reply