Coupled Numerical Model of Vibration-Based Harvester

Herein, the authors publish the complex design of a numerical coupled model of a vibration-based harvester that transforms mechanical vibrations into electric energy. A numerical model is based on usage of the finite element method, connecting analysis of the damped mechanical oscillation, electromagnetic field and electrical circuit. The model was demonstrated on the design of a microgenerator (MG), and then experimentally tested. The numerical model allows us to execute optimization of the design with many degrees of freedom. The transformation of the wave spreading in the form of mechanical vibrations was solved in the area of resonance of the electromechanical system.

One of the prerequisites for the design of the electromechanical system is a correct way of grasping the physical principle for maximum description of phenomena and processes; it is fundamental in harvesting extraction, as cases in the theses and works [17][18][19]. When comparing the majority of experimental methods and approaches [9,20], including hybrid design methods [6][7][8] with numerically modeled coupled tasks, the realized design using finite element methods (FEM) [2] is difficult to process but it leads to results with significant parameters relative to other quick approaches, as reported in work [19], Table 1, column "Effective power density [W/m 3 ]".
In Figure 1, an example is shown of a principal design of a harvester model based on the principle of electromagnetic induction using the Faraday induction law in its full range [2,19]. The analysis and results evaluation were based on the numerical model solved by the FEM. The critical parameter can be, for example, the boundary sensitivity of the generator to the minimum vibration amplitude. For the studied design the acceleration value is G = 0.01 g − 0.05 g (g = 9.81 m/s −2 ). Minimal external dimensions of the microgenerator (MG) had to be found and simultaneously, the expected volume was in the range of VMG = (10 -50) ×0.10 −6 m 3 . Further, the range of expected output effective power (RMS) was Pout =10-100mW, expected output voltage range Uout =2-20V and excitation frequency range fs =15-35 Hz. The principle of resonant arrangement of the moving MG core was used to attain a high efficiency of the vibration transformation, Figure 1. For moving path selection as linear part with non-linear magnetic braking system (Figure 2a), this concept was fully numerically modeled using the associated FEM model (112) described below. It has been shown that the technical solution of the MG core conductor designed within our approach leads to an increase in the damping coefficient and it does not have sufficient sensitivity to lowest acceleration values in the range G = 0.01 g − 0.02 g. Therefore, the concept of beam version (BV) in the arrangement, Figure 2b,c was approached in the solution. Thus the designed and tested MG achieved the expected parameters in the range G = 0.01 g − 0.02 g.
A swinging arrangement (BV, Figures 1 and 2b,c) based on bearings with a minimum achievable damping rate lb has been proposed. The swinging arrangement was damped by a non-linear element -magnetic dampers at the extreme positions, Figure 2b,c. For maximum sensitivity, the electromechanical system was tuned to enter the resonance in the supposed frequency range fs of mechanical vibrations. The critical parameter can be, for example, the boundary sensitivity of the generator to the minimum vibration amplitude. For the studied design the acceleration value is G = 0.01 g − 0.05 g (g = 9.81 m/s −2 ). Minimal external dimensions of the microgenerator (MG) had to be found and simultaneously, the expected volume was in the range of V MG = (10 − 50) × 0.10 −6 m 3 . Further, the range of expected output effective power (RMS) was P out =10-100mW, expected output voltage range U out =2-20V and excitation frequency range f s =15-35 Hz. The principle of resonant arrangement of the moving MG core was used to attain a high efficiency of the vibration transformation, Figure 1. For moving path selection as linear part with non-linear magnetic braking system (Figure 2a), this concept was fully numerically modeled using the associated FEM model (112) described below. It has been shown that the technical solution of the MG core conductor designed within our approach leads to an increase in the damping coefficient and it does not have sufficient sensitivity to lowest acceleration values in the range G = 0.01 g − 0.02 g. Therefore, the concept of beam version (BV) in the arrangement, Figure 2b,c was approached in the solution. Thus the designed and tested MG achieved the expected parameters in the range G = 0.01 g − 0.02 g.
A swinging arrangement (BV, Figures 1 and 2b,c) based on bearings with a minimum achievable damping rate l b has been proposed. The swinging arrangement was damped by a non-linear element-magnetic dampers at the extreme positions, Figure 2b,c. For maximum sensitivity, the electro-mechanical system was tuned to enter the resonance in the supposed frequency range f s of mechanical vibrations. The critical parameter can be, for example, the boundary sensitivity of the generator to the minimum vibration amplitude. For the studied design the acceleration value is G = 0.01 g − 0.05 g (g = 9.81 m/s −2 ). Minimal external dimensions of the microgenerator (MG) had to be found and simultaneously, the expected volume was in the range of VMG = (10 -50) ×0.10 −6 m 3 . Further, the range of expected output effective power (RMS) was Pout =10-100mW, expected output voltage range Uout =2-20V and excitation frequency range fs =15-35 Hz. The principle of resonant arrangement of the moving MG core was used to attain a high efficiency of the vibration transformation, Figure 1. For moving path selection as linear part with non-linear magnetic braking system (Figure 2a), this concept was fully numerically modeled using the associated FEM model (112) described below. It has been shown that the technical solution of the MG core conductor designed within our approach leads to an increase in the damping coefficient and it does not have sufficient sensitivity to lowest acceleration values in the range G = 0.01 g − 0.02 g. Therefore, the concept of beam version (BV) in the arrangement, Figure 2b,c was approached in the solution. Thus the designed and tested MG achieved the expected parameters in the range G = 0.01 g − 0.02 g.
A swinging arrangement (BV, Figures 1 and 2b,c) based on bearings with a minimum achievable damping rate lb has been proposed. The swinging arrangement was damped by a non-linear element -magnetic dampers at the extreme positions, Figure 2b,c. For maximum sensitivity, the electromechanical system was tuned to enter the resonance in the supposed frequency range fs of mechanical vibrations.

FEM Numerical Model
As mentioned in the principles and basics of the analytical design of the model [2,[17][18][19], Figure  2, it is possible to describe electromagnetic field more generally with the FEM numerical model (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11). The electromechanical vibrating harvester (34) works in special modes using resonances; its model includes nonlinearities and material nonlinearities depending on the temperature etc. As the first variant of the MG, the 3D structure with the application of magnetic damping elements was modeled and simulated. This variant did not used the cantilever beam, Figure 2a. Later, the second variant with cantilever beam swivel arrangement (BV) was also modeled a simulated.
The design of the numerical model was based on the reduced set of Maxwell equations in Heaviside's notation for quasi-stationary cases of electromagnetic field. MG elements with a high degree of non-linearity and hysteresis (permanent magnet, ferromagnetic material of pole pieces, etc.) were used in the model and the MG worked in resonance mode, which is a therefore strongly nonlinear task [2]. The FEM model has taken into account all these criteria and the results of the analyzes were used to find sensitive parameters of the mathematical model. It was necessary to take into account the accuracy of the analysis and the parameters of the MG bond with the source of vibration.
In the numerical model of the discussed MG for quasi-stationary analysis, the effect of displacement currents (Equations (4a) and (6a)) was further retained. This displacement current effect was not used for the present MG model based on Faraday's law of induction. It was taken in to account within analysis of another type MG with different parameters (higher vibration frequency) based on piezo-electric phenomenon. However, this piezo-electric generator is subject of other research. Therefore, the displacement current effect will not be used in the analysis of MG, Figure  2b,c. Next, a brief derivation of the numerical model for the tetrahedral, pentahedral and hexahedral elements of FEM will be briefly introduced. This model is using the Galerkin method of functional

FEM Numerical Model
As mentioned in the principles and basics of the analytical design of the model [2,[17][18][19], Figure 2, it is possible to describe electromagnetic field more generally with the FEM numerical model (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11). The electromechanical vibrating harvester (34) works in special modes using resonances; its model includes nonlinearities and material nonlinearities depending on the temperature etc. As the first variant of the MG, the 3D structure with the application of magnetic damping elements was modeled and simulated. This variant did not used the cantilever beam, Figure 2a. Later, the second variant with cantilever beam swivel arrangement (BV) was also modeled a simulated.
The design of the numerical model was based on the reduced set of Maxwell equations in Heaviside's notation for quasi-stationary cases of electromagnetic field. MG elements with a high degree of non-linearity and hysteresis (permanent magnet, ferromagnetic material of pole pieces, etc.) were used in the model and the MG worked in resonance mode, which is a therefore strongly nonlinear task [2]. The FEM model has taken into account all these criteria and the results of the analyzes were used to find sensitive parameters of the mathematical model. It was necessary to take into account the accuracy of the analysis and the parameters of the MG bond with the source of vibration.
In the numerical model of the discussed MG for quasi-stationary analysis, the effect of displacement currents (Equations (4a) and (6a)) was further retained. This displacement current effect was not used for the present MG model based on Faraday's law of induction. It was taken in to account within analysis of another type MG with different parameters (higher vibration frequency) based on piezo-electric phenomenon. However, this piezo-electric generator is subject of other research. Therefore, the displacement current effect will not be used in the analysis of MG, Figure 2b,c. Next, a brief derivation of the numerical model for the tetrahedral, pentahedral and hexahedral elements of FEM will be briefly introduced. This model is using the Galerkin method of functional minimization with conversion to a mathematical problem. In the design of the basic model (Equation (12)) the relative motion of the magnetic field generated by the permanent magnet with magnetization M and the induction coil is considered as the source of excitation, Figure 1.
where H is the magnetic field intensity vector, B is the magnetic flux density vector, J is the current density vector.
where E is the electric field intensity vector, D is the electric flux density vector, ρ e is the electric charge density, which is equal to zero for the considered MG and the area Ω s , Ω v ρ eΩs, Ωv . Material relations are represented by the equations, whose respect the application of permanent magnets with magnetization M in both functional parts, the main part of the MG and their damping elements respectively.
where µ 0 is vacuum permeability, µ r is relative permeability of environment, µ = µ r µ 0 , M is magnetization of permanent magnet, γ is specific conductivity of environment, ε 0 is permittivity of vacuum, ε r is relative permittivity, ε = ε 0 ε r . The temporal changes of the electric and magnetic fields in the considered model of the MG, Figure 2, are negligible according to the expected parameters. That means, the relations (4b), (6b) are not respected in the proposed numerical model. Vector functions of electric and magnetic fields are expressed using scalar electric ϕ e and vector magnetic potential A The total current density from Equation (4a) J is superposed from the circuit's excitation current density J circ and the current density from the eddy currents J v . Movement is respected in the model by current density in both parts, electrically conductive parts and electrical windings of the MG respectively The electromagnetic field model is formulated from Equations (1) to (10). Based on (1) and (10) is For individual parts, Figure 1, of the model Ω holds Ω⊂Ω v ∪Ω s , where Ω v is the region with dominant eddy currents, according to Equation (6a) Ω s is the region with known current density distribution J s . In the model under consideration: Appl. Sci. 2020, 10, 2725

of 25
Then, Equation (11a) can be modified using formulas from (1) to (10) with respect to the source of the magnetic field or the damping elements based on permanent magnets that are represented by magnetization M, Figure 1: where v is the velocity of the motion area Ω, s is the displacement vector, M is the magnetization vector and according to (4a) holds div γ −gradφ e − ∂A ∂t = 0 in the area Ω.
From Equations (3a) and (3b), where the bond between the electric and magnetic field is captured and is expressed with the help of used potentials relation From the Equation (4b), where the distribution of the electric field is captured The boundary and initial conditions will be determined as: n · (γ i gradϕ e,i ) = 0 on the boundary Γ Ω , where i, j stands for interface indexes, n · γ i gradϕ e,i − γ j gradϕ e, j = 0 on the boundary Γ i,j , i j n · (ε i gradϕ e,i ) = K 0 on the boundary Γ Ω , n · ε i gradϕ e,i − ε j gradϕ e,j = K 1 on the bounday Γ i,j , i j Appl. Sci. 2020, 10, 2725 where n is a normal vector perpendicular to the boundary of the surface area Γ, Γ i,j is the interface between the area i and the area j, Γ Ω is the interface at the outer edge of the area; the indexes i, j denotes for the quantities their belonging to the areas Ω i Ω j . Then the initial conditions are The discretization of relations (12) to (14) can be accomplished by approximating the scalar electrical potential where ϕ e is the nodal value of the scalar electrical potential, W is the base function, N ϕ is the number of nodes of the discretization network, where S is the coordinate of the node, W is the base function, N s is the number of nodes of the discretization elements, where a is the node value of the vector magnetic potential, W is the base function, N A is the number of nodes of the discretization network. Applying the approximation (20) to (22) and the Galerkin method in relation (12) to (14) gives a semi-discrete solution for the region of Ω m model The Equation (24) is modified by application of 2nd Green's formula and Gauss's theorem to expression Respecting the boundary conditions of the problem according to the expressions (18), the relation (25) changes to Applying the approximation (20) to (22) and the Galerkin method in relation (17a) gives a semi-discrete solution for the region of Ω m model The expression (27) is modified by application of 2nd Green's formula and Gauss's theorem to expression Respecting the boundary conditions of the problem according to the expressions (18), the relation (25) changes to Substituting approximation functions, A, s, ϕ ε according to (20) to (22) into (26) and (28) gives a semi-discrete solution Relation (29) can be rewritten to form Appl. Sci. 2020, 10, 2725 8 of 25 The simplified notation of the system of Equations (30) is The expression is easy to rewrite to form To simplify the calculation algorithm in the numerical part of the solution of the system of equations and its acceleration, the system (33) is converted to the form The form of the coefficients is expressed in relations (30) and (31).
Appl. Sci. 2020, 10, 2725 9 of 25 where i, j, k are the base vectors of the used Cartesian coordinate space. The algorithm for constructing matrices of coefficients (35) to (44) is simplified when the system of Equations (31) is rewritten to form The coefficients of the system of the Equations (45) and (46) are expressed in the relations The system of Equations (45) and (46) is written in the form The procedure of quantification of the coefficients for tetrahedral, pentahedral and hexahedral element is described in detail [3].
The short description of the elements is shown in Figures 3-5. The tetrahedral element in Figure 3 has a base function W The system of Equations (45) and (46) is written in the form The procedure of quantification of the coefficients for tetrahedral, pentahedral and hexahedral element is described in detail [3].
The short description of the elements is shown in Figures 3-5. The tetrahedral element in Figure  3 has a base function W ( ) Make use of where the indexes i, j are node numbers of the element e. If a pentahedral element is selected according to Figure 4 and the function W is displayed on it, then their writing at the selected local Cartesian coordinate system o,ξ,η,ζ is as follows: and the gradient of the function W is expressed in the Cartesian coordinate system o, x, y, z After the transformation of the function W into the global Cartesian coordinate system o, x, y, z from the local coordinate system o,ξ,η,ζ from the expressions (30) and (33) are the coefficients of the system of Equations (57). In them are In the (77) relation, the derivative of the position vector Rp is represented by the following expression the coordinates x, y, z and their derivatives according to the local coordinates ξ,η,ζ are expressed  Integration in the relations (80) to (89) can be solved analytically. This method is difficult and time consuming. Avoiding the use of analytically expressed integrals in most cases will not significantly reduce the accuracy of the solution of the system of Equations (57). Another easier solution is to apply the Gaussian quadrature. Its form, usage and properties are known [3]. The coefficients of the system of Equations (57) are then calculated for a numerical solution where the indexes at x, y, z coordinates are local node numbers according to Figure 3. The coefficients of the function W can be expressed a i = x j y j z j x k y k z k x l y l z l i, j, k, l = 1, . . . , 4 where the indexes i, j, k, l change cyclically over a given interval. Index i is the natural number of the element function, indexes j, k, l are coordinate indexes of the element nodes. The coefficients of the model are for the tetrahedral element of Figure 3 k where the indexes i, j are node numbers of the element e. If a pentahedral element is selected according to Figure 4 and the function W is displayed on it, then their writing at the selected local Cartesian coordinate system o,ξ,η,ζ is as follows: and the gradient of the function W is expressed in the Cartesian coordinate system o, x, y, z After the transformation of the function W into the global Cartesian coordinate system o, x, y, z from the local coordinate system o,ξ,η,ζ from the expressions (30) and (33) are the coefficients of the system of Equations (57). In them are In the (77) relation, the derivative of the position vector R p is represented by the following expression ∂R p the coordinates x, y, z and their derivatives according to the local coordinates ξ,η,ζ are expressed When the hexahedral element from Figure 5 is selected, its base functions are written as For the hexahedral element, the sizes of the coefficients from the set of Equations (57) are calculated from Equations (80)-(89) using the functions (90).
Integration in the relations (80) to (89) can be solved analytically. This method is difficult and time consuming. Avoiding the use of analytically expressed integrals in most cases will not significantly reduce the accuracy of the solution of the system of Equations (57). Another easier solution is to apply the Gaussian quadrature. Its form, usage and properties are known [3]. The coefficients of the system of Equations (57) are then calculated for a numerical solution The number of integration points of the local coordinates for the tetrahedral element was chosen n g = m g = l g = 5. This number according to tests for accuracy and speed of quantification of quadrature relation fits the mathematical model. The expression of the auxiliary function L is identical to the function (58). By means of the L function, the coordinates of the integration points are determined according to the formulas (102) to (104) below.
Set values of the auxiliary function L and the weight function H for 5 selected points within the area Ω e are 1. point : The coordinates of the integration point are determined for the tetrahedral element Six points of integration n g = m g = l g = 6 were chosen for the pentahedral element. The number corresponds to the required accuracy and speed of the quadrature calculation. The coordinates of the points are determined ξ g1 = ±0.577350269189626, Similarly, for a hexahedral element, where n g = m g = l g = 8. The coordinates of the point for calculating the quadrature according to (n102) to (n104) is expressed by the auxiliary function L The electromagnetic part of the model of the task MG can be characterized by the expression (57). The description of the mechanical part of the model by means of the FEM is known for example from theses and works [1,2] and the system [21], then the electromechanically bounded equation is written in the form where f is the vector of specific force, λ 1 , λ 2 are auxiliary scalar functions. For the geometrical model, Figure 2, the external mechanical load was not considered. Such presumption will lead to f 0 = 0. But, to respect the position of the MG in the earth's gravitational field, with the acceleration g a static force, it was considered f 0 = M g V g ·g, where m g is the mass of the MG core, V g is the volume of the MG core. The specific magnetic force f m is expressed as There are valid boundary and initial conditions for the model, for example known [21]. Similar to the above, mentioned derivation of the electromagnetic part of the model, the coupled dynamical mechanical-electromagnetic model can be written in the form Appl. Sci. 2020, 10, 2725 20 of 25 and the corresponding semi-discrete model solution After deriving and quantifying the corresponding coefficients of the system of equations, we obtain the overall FEM model with the corresponding matrices of coefficients K, M. φ column matrix of unknown second derivatives according to time of scalar electric potentials in nodes of discretization network, F is column matrix of forces, p σ is column matrix of vectors of tangential pressure components, p n is column matrix of vectors of normal pressure components E Te is a column matrix of modulus vectors of specific deformation in temperature dependence of material.
Then the overall numerical FEM model of the electromechanical coupled task of the vibration MG is written as a system of equations This problem can be solved in ANSYS [21] using programming tools and the APDL internal language. The designed 3D model (112) was used to analyze and evaluate the transient excitation state of MG Figures 2a and 6. Then, a steady state analysis of resonance and the corresponding electrical load R L was performed [18,19]. Sensitive parts of the MG model design have emerged from 3D analyzes and experiments. Due to the magnitude and amount of cumulative computations, absolutely not all details are given, but only the key details (the magnetic field in the MG core region). In order to find a suitable (optimal) direction of the magnetic arrangement solution, an accelerated 2D FEM model and a magneto-static field design were used, which in turn lead to maximizing the induced voltage u(t) in the winding of the MG coil.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 23 of 27 and a magneto-static field design were used, which in turn lead to maximizing the induced voltage u(t) in the winding of the MG coil. The 3D model was designed as the primary version of MG, Figure 6, which was not optimally arranged in the electro-mechanical conceptual design. It was solved both for start-up as a transient task and for the so-called steady-state regime in the resonant state by the model (112). Given the previous experience with the complex analysis of the associated 3D model of the first version of MG, Figure 6, it became clear that critical part is the design of magnetic circuit of the MG (Figures 1 and 8). For the chosen concept of the beam version of MG, Figure 7, only 2D static nonlinear problem with hysteresis, Figure 8a was analyzed because of the higher speed of connected simulation. A conceptual solution of the magnetic circuit was proposed, using the results of performed analyzes, Figure 8b. Thus, a design with a maximum magnetic flux density B(t) and maximum amplitude of the electric voltage u(t) on the output of the MG's winding has been achieved The 3D model was designed as the primary version of MG, Figure 6, which was not optimally arranged in the electro-mechanical conceptual design. It was solved both for start-up as a transient task and for the so-called steady-state regime in the resonant state by the model (112).
Given the previous experience with the complex analysis of the associated 3D model of the first version of MG, Figure 6, it became clear that critical part is the design of magnetic circuit of the MG (Figure 1 and Figure 8). For the chosen concept of the beam version of MG, Figure 7, only 2D static non-linear problem with hysteresis, Figure 8a was analyzed because of the higher speed of connected simulation. A conceptual solution of the magnetic circuit was proposed, using the results of performed analyzes, Figure 8b. Thus, a design with a maximum magnetic flux density B(t) and maximum amplitude of the electric voltage u(t) on the output of the MG's winding has been achieved in a relatively rapid manner, Figure 7. According to these auxiliary FEM models and analyzes, the final concept of MG was designed, Figures 2b, 7 and 8.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 24 of 27 in a relatively rapid manner, Figure 7. According to these auxiliary FEM models and analyzes, the final concept of MG was designed, Figures 2b, 7 and 8.

Modeling the MG
A free-arm system (BV) and its magnetic damping (non-linear) [19], Figure 7, has been designed. Details of the MG evaluation of the analysis can be demonstrated for example in the magnetic field the distribution of the magnetic flux density module B [T] in Figure 8. The realized MG function sample, Figures 7 and 8, was tested for the frequency range f s = 30-45Hz. Typical frequency dependencies of output voltage and output power of the MG are shown in Figure 9. The graph shows the effective value of the output power P_RMS and the effective electric voltage V_RMS at the MG's winding output loaded by R L = 2.7KΩ at the vibration G = 0.5 g RMS .

Modeling the MG
A free-arm system (BV) and its magnetic damping (non-linear) [19], Figure 7, has been designed. Details of the MG evaluation of the analysis can be demonstrated for example in the magnetic field the distribution of the magnetic flux density module B[T] in Figure 8. The realized MG function sample, Figures 7 and 8, was tested for the frequency range fs = 30-45Hz. Typical frequency dependencies of output voltage and output power of the MG are shown in Figure 9. The graph shows the effective value of the output power P_RMS and the effective electric voltage V_RMS at the MG's winding output loaded by RL = 2.7KΩ at the vibration G = 0.5 g RMS.

Conclusions
The paper presents a detailed design of a numerical model of a FEM mini-generator for harvesting energy from vibrations within the range of G = 0.05 g − 0.08 g, for the frequency range of fs = 10-100 Hz, fmain = 36-41 Hz.
Based on the proposed beam version MG models and their analyses, evaluation of the parameters, qualitative conclusions were applied to design of the swinging arm of mini-generator, Figure 7. Some partial specific results from this systematic approach to generator design were published in theses and works [17][18][19]. This article describes a complete FEM numerical model that was used for vibration MG analysis. The results of the generator analysis were experimentally verified and subsequently utilized to build the functional MG. A specific harvesting system for mechanical energy extraction was proposed. It utilizes the principle of transformation of propagating mechanical wave to electrical energy. The designed numerical model was implemented in ANSYS simulation software [21]. The simulated effective transformed mass and volume power densities of the generator are peff = 266 Wm -3 and peffV = 0.31 Wkg −1 respectively.

Conclusions
The paper presents a detailed design of a numerical model of a FEM mini-generator for harvesting energy from vibrations within the range of G = 0.05 g − 0.08 g, for the frequency range of f s = 10-100 Hz, f main = 36-41 Hz.
Based on the proposed beam version MG models and their analyses, evaluation of the parameters, qualitative conclusions were applied to design of the swinging arm of mini-generator, Figure 7. Some partial specific results from this systematic approach to generator design were published in theses and works [17][18][19]. This article describes a complete FEM numerical model that was used for vibration MG analysis. The results of the generator analysis were experimentally verified and subsequently utilized to build the functional MG. A specific harvesting system for mechanical energy extraction was proposed. It utilizes the principle of transformation of propagating mechanical wave to electrical energy. The designed numerical model was implemented in ANSYS simulation software [21]. The simulated effective transformed mass and volume power densities of the generator are p eff = 266 Wm -3 and p effV = 0.31 Wkg −1 respectively.
The FEM model showed the key moment in the design of the first versions of MG concept and its design solutions. It also pointed out the key parameters of the model and directed the modifications of the MG's equipment to the optimal solution of the harvester. An associated model (112), was able to detect the sensitive points of the task. This was proved by measurement results from experiments with the MG sample. FEM models of lower geometrical levels (2D) and simplified models (non-linear with hysteresis, etc.) were used as a necessary tool for quick estimation of the changes of design or concept changes.
Author Contributions: P.F. contributed to the theoretical section, numerical model, and experiments; he also wrote the paper; Z.S. conceived and designed the experiments; J.D. contributed to the optimization procedures; P.E., J.Z. and R.P. evaluated the experiments and graphics. All authors have read and agreed to the published version of the manuscript.

Funding:
The research was financed by the BD 2020-2022, FEKT-S-20-6360, the infrastructure of the SIX Center was used.