Framework for efﬁcient dynamic analysis applied to a tubular generator for suspension applications

This paper considers a slotless three-phase tubular permanent magnet generator located in an automotive suspension system for the application of vibration energy harvesting. A two-dimensional ﬁnite element method model of the harvester is produced and an experimental setup that contains the generator is constructed. Signal decomposition methods are applied to measured suspension displacement data and the resulting signal components are used as input for the model. Two approaches for signal decomposition are discussed, namely, the discrete Fourier transform and the continuous wavelet transform. The individual emf responses of the model are reconstructed to a single output, while a sideband prediction algorithm accounts for the non-linearities in the system. The simulation results are compared with the reference measurements conducted on the setup to determine the accuracy of each of the signal decomposition methods.


INTRODUCTION
In the automotive industry, active electromechanical suspension systems are gaining interest [1,2]. These devices can apply active forces on the wheels and the vehicle body to oppose road irregularities, which improves handling and comfort [3,4]. Furthermore, the energy that is normally dissipated in passive suspension solutions can be partly recuperated by such electromechanical machines, making the latter particularly attractive from an efficiency point of view [5][6][7][8].
Finite element analysis (FEA) is a mandatory step to accurately model and design electromechanical devices [9][10][11]. Moreover, to fully analyze the energy harvesting capabilities of a vibration harvester utilized in vehicle suspensions, the corresponding finite element method (FEM) model should be assembled, solved, evaluated, and intensively iterated for different road conditions. Unfortunately, road profile data is random in nature, meaning that fine time-stepping is required to account for all the frequency content and transient phenomena, which translates to a massive computational burden. Therefore, analysis utilizing signal component separation and superposition is preferred to speed up FEA simulations, as the number of required solver iterations to account for all the frequency content is far less when simulations are performed in the frequency domain. However, transforming the linear translation movement into sinusoidal emf components is not possible without leading to some inaccuracies due to the cosine term in the emf equation for linear machines: the long-stroke movement of the original profile covers the whole cosine, while the short-stroke movement of the individual components enables the linearization of this term [12]. Additionally, the non-linear material characteristics, as well as the induced eddy currents in the harvester, add to the discrepancy.
The slotless tubular permanent magnet machine, as presented in [13], is selected as a candidate for vibration energy harvesting in automotive suspension applications. Besides the non-linear soft-magnetic material characteristics, the tubular permanent magnet generator (TPMG) has a slit in the back iron to account for the connection of the phase wires to the power electronics unit. This removes the axisymmetry, which will lead to additional inaccuracies when a 2D model is considered. As expected, the slit blocks the eddy current path, influencing the magnetic loading of the device. Moreover, the machine performance is enhanced due to the minimization of the eddy current flow in the circumferential direction.
This paper is an extended version of the paper entitled "Dynamic Analysis of a Tubular Generator for Automotive Suspension Applications," which was presented at the 12th International Symposium on Linear Drives for Industry Applications (LDIA) and, subsequently, published in the conference proceedings [14]. A framework for efficient dynamic analysis by utilizing decomposition and reconstruction of signal components is explained. The presented framework is applied to the aforementioned non-linear problem. Measured suspension displacement data is considered as input and is either decomposed into its harmonics using the discrete Fourier transform (DFT), as proposed in [14], or decomposed into daughter wavelets utilizing continuous wavelet transforms (CWT). Each of the resulting signal components is used as input for a transient 2D non-linear FEM model [15] of the TPMG, which provides an induced emf response per signal component. Reconstruction of the resulting emf elements into a single waveform provides meaningful output data that directly relates to the original input profile. Furthermore, a sideband prediction algorithm (SPA) deals with the emergence of higher harmonic content due to non-linearities in the system. From the resulting emf waveform, the harvested energy is calculated. The described approach requires less computational effort compared to a single simulation where the complete reference profile is considered. For validation, the simulation results of the approximation methods are compared with the measurement results that are obtained by applying the original reference signal to the experimental setup. These validation measurements demonstrate the potential of the framework and prove the functionality of dedicated solvers in the transformed domain.

MODELLING FRAMEWORK
Harmonic balance (HB) [16,17] and wavelet balance (WB) [18] are dedicated solvers that take full advantage of the transformed domain, meaning that the computational effort is mainly scaled by the number of coefficients or basis functions, which could greatly reduce the computational burden for FEM models. Unfortunately, such solvers do not exist yet in FEA software. The focus of this paper, therefore, lies with the demonstration of the benefit the harmonic or wavelet analysis can bring to electromagnetic simulations, as compared to the time-stepping FEM. Transient time-stepping FEM simulations are used in this paper to obtain the system response, while accounting for the non-linear material characteristics and the effects of magnetic saturation. The reduction in the number of time-steps is realized by signal decomposition based on the DFT and the CWT. The discrepancies of these approximation methods are proven to be small enough, demonstrating the use of the aforementioned dedicated solvers in FEA domains. Instead of solving the system with one large transient FEM simulation, the proposed modelling framework decomposes the random input signal in the time representations of a finite number of basis functions. These signal components are applied to multiple smaller transient FEM simulations, which once combined, require less computational effort than the single complete simulation. The individual simulation outputs are used for the reconstruction of the system response, while a SPA deals  with the emergence of higher frequency content due to nonlinear phenomena. The considered signal decomposition methods are the harmonic approximation using the DFT and the wavelet approximation utilizing the CWT. A schematic of the framework is provided in Figure 1, where, N 1 ≫ N 2 ≫ N 3 . Figure 2a illustrates a quarter car situation, where the TPMG, provided in Figure 3, is placed in parallel with the suspension system. Here, M s,q represents the quarter sprung mass, k a and d a , are the spring and damping coefficients of the suspension system, z is the suspension displacement due to the road vibrations, M u is the unsprung mass and, k t is the spring constant of the tire. The influence of the generator dynamics on the suspension displacement is neglected. Additionally, the power electronics (PE) of the harvester are modelled as a resistive circuit, which is illustrated in Figure 2b. Finally, a transient axisymmetric 2D FEM model of the TPMG is selected over a transient 3D FEM model, thus, excluding the slit, which limits the mesh size and, consequently, the computational effort. Figure 4 provides a single period of the TPMG, while the full stage schematic is illustrated in Figure 5.

Input profile
The suspension displacement data measured in [3] by driving a benchmark BMW 530i vehicle around the campus of Eindhoven University of Technology, is considered as input to the system. The measurement noise is filtered by limiting the frequency range of the input profile to the operating range of the vehicle. The signal is applied to a band-pass filter, where the cut-off frequencies are defined by the frequencies between the  Suspension displacement profile as measured in [3] after filtering sprung mass resonant frequency and the unsprung mass resonant frequency together with an offset of 20 %. The mechanical parameters correspond to the values described in [3,13]. The filtered suspension displacement profile is shown in Figure 6, from which it is observed that most of the signal is within one geometrical period, as the pole pitch of the TPMG, p , equals 8.33 mm.

HARMONIC APPROXIMATION
The first method that is considered is the harmonic approximation (HA), which decomposes the input signal into a sum of sinusoids, or harmonics, using the DFT algorithm. To make the graphs more visible, a window of the filtered input data is selected. Additionally, the windowed signal has equal magnitude for the first and the last time-step, indicating a periodic signal, which guarantees better accuracy of the Fourier transform [19].

Decomposition
Using the DFT, the windowed signal, as shown in Figure 7a, is transformed to the frequency domain. Each frequency component is extracted from the resulting frequency spectrum, given by Figure 7b, and transformed back to the time domain using the inverse discrete Fourier transform (IDFT). To visualize this, Figure 8a shows in red the individual frequency components that belong to the four harmonics with the most signal energy, Each harmonic has the same sampling as the original signal, which means that the number of time-steps of each sinusoid equals the number of time-steps of the original signal. However, the harmonics can be truncated to a single period, lowering the number of time-steps. Additionally, the number of time-steps can be reduced further by down-sampling these single period sinusoids according to the Nyquist-Shannon sampling theorem. The theorem states that any perfect band limited signal can be represented by at least a number of samples that is equal to twice the highest significant frequency component In this application, the highest significant frequency component of the suspension displacement equals the unsprung mass resonant frequency, provided by (2). Unlike the individual sinusoids, the measured suspension displacement is not perfectly band limited, hence, this theorem cannot be applied to it directly. Figure 9a illustrates the truncation of the sinusoids to their sin- of 1.9 % is induced, where, u re f is the reference signal, N is the number of time-steps of the signals and, u recon is the reconstructed signal utilizing the harmonic decomposition.

Reconstruction
Transient simulations in FEA software are conducted using each of the 29 suspension displacement harmonics, resulting in 29 induced emf waveforms in the phase windings of the harvester. These 29 emf responses are then superposed, which yields a single emf response that results from the original suspension displacement. However, additional steps are required to properly account for the higher frequency content that is induced by the non-linearities in the system. The single-phase emf response of a linear machine is given by where, k m is the motor constant, v is the speed that is linked to the suspension displacement z, and p is the pole pitch of the harvester. For each of the 29 suspension displacement harmonics, both z and v are single frequency sinusoids with a signal amplitude smaller than p , which is shown in Figure 9b. Thus, only small perturbations are applied to the generator. Consequently, the frequency spectrum of the resulting emf waveform will contain a fundamental frequency component, as well as its odd harmonics, where the latter ones are negligible in magnitude. However, in reality, the complete signal for the windowed suspension displacement of Figure 7a is applied to the system, so both z and v are signals that contain multiple frequencies, meaning that higher harmonic content, also known as sidebands, will emerge due to the frequency interactions. This higher harmonic content is not represented by the simulations and the linear superposition of the individual harmonics alone and, thus, it has to be predicted and added to the linear superposition of the emf components to yield proper reconstruction results. To realize such a SPA, the cosine term in (5) has to be simplified, for which the Taylor series can be used. From Figure 6 it is known that most of the signal is within one geometrical period. Hence, a Taylor approximation of order 4 of the cosine function is sufficient. Even if z would be larger than p , the Taylor approximation still holds by introducing more terms. Indeed, the convergence radius for the Taylor series of a trigonometric function is infinite. To demonstrate the approach for the SPA, a sum of two sinusoids is considered as input suspension displacement. The singlephase emf response, depending on frequencies x and y, can be expressed as e(x, y) = (a cos(x) + b cos(y)) cos(c sin(x) + d sin(y)). (7) Here, the substitutes are as in A Taylor approximation of order 4 of the cosine function is applied to the right-hand cosine term of (7), where The three resulting Taylor terms are represented by T 3 = (a cos(x) + b cos(y)) (c sin(x) + d sin(y)) 4 4! .
By rewriting the Taylor terms, the higher harmonic content due to the speed-displacement interaction is obtained. For example, the second term, T 2 , is written as From (16) it is observed that, according to the second Taylor term, four sidebands will emerge at frequencies, x − 2y, 2x − y, 2x + y, and, x + 2y, next to the two fundamentals at frequencies x and, y, and the two odd harmonics at frequencies 3x and 3y. Additionally, (16) provides information on how the signal amplitudes are affected by the frequency interaction of the input signal as well as the effects of saturation and the non-linear material characteristics, since k m in (9) is recalculated for each suspension displacement harmonic. Furthermore, (5) could be extended to allow for the prediction of even more harmonic content, which could be used to represent the effect of slotting in iron-cored machines. If the aforementioned approach is extended to all Taylor terms for the 29 components of the suspension displacement of Figure 7a, the SPA is realized. The SPA ensures that the higher harmonic content, as well as the deviations in the signal amplitudes, resulting from the non-linearities in the system can be incorporated into the emf reconstruction, which resolves the discrepancy induced by the linear superposition. This will be numerically and experimentally validated in Section 5.3. Moreover, the use of the SPA is useful to represent the cross-terms in an HB solver.

WAVELET APPROXIMATION
The second method that is considered is the wavelet approximation (WA). The WA decomposes the input signal utilizing the CWT algorithm, which is a sequence of projections of the input signal onto rescaled and translated versions of an analyzing function or wavelet [20]. The advantage of the CWT over the DFT is the fact that the latter only provides frequency localization, while wavelets offer time-frequency localization. Moreover, due to the finite support of wavelets, their application is not limited to periodic functions, which makes it more beneficial for the considered application.

Decomposition
The analyzed signal is decomposed in a scaled version of a basis function. In MATLAB three wavelet families are implemented for the CWT, from which the Morse wavelet is considered to be a superfamily of wavelets [21]. However, due to the limited implementation of the Morse wavelet in MATLAB, another basis, the Morlet wavelet is selected instead. Figure 10 shows the shape of the standard Morlet wavelet in blue. Additionally, due to the number of redundant time-steps of the Morlet wavelet, the possibility exists to down-sample this wavelet by a factor 2, which is illustrated by the red waveform in the figure. The down-sampled Morlet now consists of 56 time-steps and is used for further signal decomposition. The CWT algorithm, which is a convolution of the input signal with multiple scaled versions of the Morlet wavelet, is applied to the windowed signal, as shown in Figure 7a. It results in 81 output waveforms, each corresponding to the convolution of the input signal with one of the 81 scaled Morlets. Adding the 81 waveforms results in the original signal with a difference of less than 5.0×10 −3 %, but only seven convolution outputs and their Morlets are needed to approximate the original signal with a discrepancy of 3.7 %. Figure 11a provides four of these scaled Morlets as an example, while Figure 11b shows the corresponding convolutions. By applying this method the suspension displacement profile, consisting of 1890 time-steps, can now be simulated by utilizing seven Morlets, that once combined consist of 392 time-steps. A reduction in the number of time-steps of factor 4.8 is realized, which is more than double the factor achieved by the HA.

Reconstruction
Each of the seven suspension displacement Morlets is applied to the FEM model, which results in seven induced emf Morlets. These Morlets are again convolved with the original suspension displacement to provide emf transform outputs. Figure 12a provides four of the emf Morlets resulting from the FEM simula- tions, while Figure 12b shows the corresponding transform outputs. Summing the results provides the complete emf response that results from the suspension displacement profile. Additionally, due to the small amplitude of the suspension displacement Morlets of Figure 11a, the FEM simulations are within smallsignal approximation region, meaning that the higher harmonic content is not accounted for. Hence, the SPA, as derived before is added to the summation of the transforms to yield proper reconstruction results. Moreover, the SPA can be utilized to represent the cross-terms in a WB solver.

VALIDATION
To prove the accuracy of the two considered approximation methods for efficient FEA of systems with random input signals, a transient 2D FEM model of the harvester is created and an experimental setup that contains the generator is constructed.  (2) The TPMG together with the measurement equipment (3) The driver unit

FEM model
Utilizing the materials, parameters, and dimensions of the TPMG, as provided in [13], a transient 2D FEM model of the harvester is produced. A single period of the TPMG, together with the full device, is previously shown in Figures 4 and 5.
The generator consists of six periods and the number of pole pairs equals 8.5. Furthermore, the non-linear material characteristics and the end effects are included in the FEM model. On the other hand, the eddy currents are excluded from the model, as [13] indicates that the influence of the slit on the magnetic loading is negligible. The profile that is used for validation measurements in [13] is also applied to the FEM model, from which the accuracy of the model is observed. Utilizing (4), a discrepancy of 1.6 % is calculated between the induced emf waveforms, thus validating the FEM model.

Horizontal setup
The horizontal setup, provided in Figure 13, is constructed to obtain experimental data that act as reference, since both the slit and all other non-linearities are present in the setup. The setup consists of a tubular actuator that can mimic suspension displacement profiles using software coupled to a driver unit. The driver is connected to a PC from which the suspension displacement profiles can be uploaded. Consequently, the commutation of the actuator is controlled by the driver unit to allow for the simulation of the profiles with a tracking error of less than 1 mm. Due to the mechanical connection between the actuator and the TPMG, emf is induced in the phase windings of the harvester, which in its turn is measured using a dSPACE system. Moreover, an optical encoder measures the speed and the displacement of the generator, while a load cell that couples the moving parts allows for the measurement of the forces. Figure 14 shows the induced emf waveform of the generator as measured by the setup, as well as the emf waveform that is calculated by the FEM model when the displacement profile  of Figure 7a is considered. A discrepancy of 5 % is observed between the two waveforms. The harvested energy can be derived from the emf waveform, where, R phase , is the measured phase resistance of the generator. Load matching is applied to maximize the harvested energy, meaning that the load resistance and the phase resistance in Figure 2b are of equal magnitude (5.9 Ω) [12]. The emf waveform and the corresponding harvested energy measured with the setup acts as reference for the decomposition algorithms.

Approximation accuracies
With the reference magnitudes and waveforms for the emf and the harvested energy known from the setup, the accuracy of the HA, as well as the WA, can be verified. Table 1 provides the magnitudes for both the harvested energy levels, as well as the induced discrepancies obtained by the different methods. Furthermore, Figure 15a compares the reference emf measurement from the setup with the result of the superposition of the 29 individual harmonics, utilizing the FEM model, while Figure 15b provides the results when the SPA is added. Additionally, Figure 16 shows the related frequency spectrums. Finally, Figures 17 and 18 provide comparisons when the WA is considered. From the table and figures, it is observed that only the linear superposition induces large errors for both approaches. This is caused by the aforementioned frequency interaction due to the non-linearity of the reference input profile, which is not represented by the linear superposition of the individual signal components. Therefore, when the SPA is implemented, the results become much more accurate. Furthermore, the WA performs better than the HA as the reduction factor in the number of time-steps is twice as high, while the reconstruction discrepancy is lower.
Moreover, by applying both approximations to the complete signal in Figure 6, it is observed that the WA really outperforms the HA. The number of time-steps for the complete measured suspension displacement profile equals 60,000, and if the HA is considered for the complete profile, 750 harmonics, that once combined, consist of 21,750 time-steps, would be needed to simulate the emf waveform. This is a reduction in the number of time-steps of factor 2.8, which is of the same order as the windowed signal. Alternatively, only eight Morlets, which once combined consist of 446 time-steps, are needed to simulate the emf response: a reduction of factor 134 is reached. Additionally, Table 2 summarizes the harvested energies and discrepancies for the different methods when the full suspension displacement signal is considered. It is observed that the  WA with the SPA has lower discrepancies compared to the HA with the SPA, while the reduction in the number of time-steps is an order of magnitude higher. Furthermore, Figure 19 shows the convergence of the two methods for both the windowed signal and the complete suspension displacement: the reduction in emf discrepancy is illustrated when more basis functions and, thus, more time-steps are considered, which clearly shows the limit of the HA in terms of scalability.

CONCLUSION
This paper has investigated the numerical error induced by two approaches to speed up transient FEA simulations for a nonlinear system. Harmonic approximation based on the DFT, as well as wavelet approximation based on the CWT, has been applied to measured suspension displacement data. The corresponding signal components are used as input for a transient 2D non-linear FEM model of a vibration energy harvester in automotive suspension solutions. The outputs for both methods are reconstructed into a single waveform that defines the system response, while a sidebands prediction algorithm resolves the discrepancy due to non-linear phenomena in the system. The final results are compared with an experimental setup to quan-tify the accuracy of the approximations. It should be noted that for the reference input profile a deviation of 5% was observed in the induced emf waveform between the reference FEM model and the measurements. For the HA, a discrepancy of 7.5% is observed in the induced emf waveforms, along with a 5.1% difference in harvested energies, while the simulation time is reduced with a factor of 2.3 utilizing the same transient solver if the windowed signal is considered. The reduction factor of the HA remains the same order when the algorithm is applied to the complete suspension displacement profile, but the discrepancy in emf increases to 6.9% and the discrepancy in harvested energy increases to 9.2%. Applying the WA, a difference of 4.9% is observed in the induced emf waveforms, together with a 2.7% discrepancy in harvested energies, while simulation time is reduced with a factor of 4.8 for the windowed profile. This is increased to a factor 134 for the longer input signal, with a discrepancy of 5.8% in terms of emf and 4.5% discrepancy in terms of harvested energy. These numbers show the limitations of the HA and, on the other hand, the benefits of the WA. Moreover, this research paves the way for the use of dedicated solvers in the transformed domain: to simulate the full suspension displacement, 750 coefficients are needed for the harmonic balance, which is impractical. On the other hand, only eight coefficients are needed for the wavelet balance in the considered application, thereby reducing the computational effort greatly for random input signals.