A Novel Approach to Predict the Structural Dynamics of E-Bike Drive Units by Innovative Integration of Elastic Multi-Body-Dynamics

Abstract

This paper presents a novel approach to address noise, vibration, and harshness (NVH) issues in electrically assisted bicycles (e-bikes) caused by the drive unit. By investigating and optimising the structural dynamics during early product development, NVH can decisively be improved and valuable resources can be saved, emphasising its significance for enhancing riding performance. The paper offers a comprehensive analysis of the e-bike drive unit’s mechanical interactions among relevant components, culminating—to the best of our knowledge—in the development of the first high-fidelity model of an entire e-bike drive unit. The proposed model uses the principles of elastic multi body dynamics (eMBD) to elucidate the structural dynamics in dynamic-transient calculations. Comparing power spectra between measured and simulated motion variables validates the chosen model assumptions. The measurements of physical samples utilise accelerometers, contactless laser Doppler vibrometry (LDV) and various test arrangements, which are replicated in simulations and provide accessibility to measure vibrations onto rotating shafts and stationary structures. In summary, this integrated system-level approach can serve as a viable starting point for comprehending and managing the NVH behaviour of e-bikes.

1. Introduction

Transportation alternatives are becoming increasingly important. Electric bicycles (e-bikes) have shown their potential in terms of sustainable and healthy mobility. Higher product expectations accompany the customer market growth [1], such as having a relatively silent ride without disturbing Noise, Vibration and Harshness (NVH).

Regarding NVH, the e-bike drive unit is mainly responsible for system excitation and the e-bike frame for sound radiation (see Figure 1) [2,3]. In detail, the drive unit’s rotor dynamic forces act along a transfer path towards the interface section of the frame connection. The resulting interface forces excite the frame structure, and the frame’s orthogonal surface velocities stimulate the surrounding air. Since a high excitation level is crucial for the overall NVH behaviour, this study focuses on the structural dynamics of the e-bike drive unit. These dynamics comprise the rotor dynamic interactions and excitations. The interface forces towards the bicycle frame define the system boundary.

An identical drive unit raises up to 10 dB higher sound pressure levels of single frequencies when mounted to a bicycle frame [4]. Depending on the excited mode shapes, the frame’s resonances amplify the vibrations [5]. Moreover, various similar manufactured drive units can vary in the level of sound radiation, assuming tolerance deviations and temperature impact [6,7,8]. Since the housing covers and seals all inside laying components, the investigation does not incorporate the direct air-borne sound emission. We assume that the housing reduces considerably direct sound emission through a high acoustic impedance, according to [5].

In order to have the most attainable influence on potential product optimisations, it is essential to understand system interrelationships and to recognise critical dynamic system states at an early stage. Conducting measurements to identify these relationships is often complex, time consuming or impossible due to accessibility issues or non-available samples. Therefore, a simulation approach is reasonable, which can predict the structural dynamics. Currently, simulation software is available for such structural dynamic calculations, but there are no approaches for applying them to an e-bike drive unit system. Hence, the primary motivation of this study is to identify the system interactions and to develop a structural dynamic and predictive simulation model of the e-bike drive unit. The target values are the vibration characteristics resulting from both the system motion and the excited mode shapes during operation.

The state of the art of NVH simulations is the coupling of computer-aided engineering (CAE) methods across multi-physical domains. The principal simulation fields considering electric drive trains are electromagnetics, structural dynamics and acoustics [9,10,11,12]. Alternative methodologies are evident in distinct domains. As a representative example, Pasch et al. (2020) [13] consider the fluid dynamics instead of electro-magnetic excitation.

Various electromagnetic simulation approaches study the structural dynamic impact on the motor NVH performance [14,15,16,17,18,19,20]. These calculations are significantly computation intensive and refer only to the rotor and stator geometry. Jaeger et al. (2020) [10] show a viable approach to retroactively couple electrical and structural dynamic domains, forming a NVH system model using co-simulations.

Given the complexity and extensive research in the domain of electric motor modelling, this study aims to concentrate exclusively on the structural dynamics. To address the electromagnetic excitation, a simplified approach is employed. The emphasis on the accurate modelling of these dynamics is essential to make any acoustic assessment. Consequently, acoustic aspects have not been taken into consideration so far. In their work, Kamper and Wegerhoff et al. [2,3] demonstrate a hybrid approach for acoustic evaluation. They incorporate the excitation of an e-bike drive unit as in situ blocked forces (interface forces to the e-bike frame), but this approach deviates from the objective of developing a predictive drive unit model.

The structural dynamics are usually modelled applying finite element methods (FEMs) or elastic multi-body dynamics (eMBD) combining rigid body motions with flexible bodies. The integration of flexible bodies allows small and linear body deformations relative to a local reference frame, while the local reference frame can perform large and non-linear global motions [21]. The eMBD simulation is faster than a complex FEM simulation due to the reduced number of motion equations. The flexible bodies are pre-processed by modal superposition based on Craig–Bampton [22] and mode shape orthonormalisation. Through the system size of the e-bike drive unit and the expectation of the non-linear system behaviour, e.g., by clearance and geometry-based non-linear stiffness, an eMBD approach is chosen.

To ensure precise modelling, a prior understanding of the system’s behaviour is mandatory to effectively map relevant components and fit their parameters accurately [23]. Therefore, this research deals with two necessary objectives: on the one hand, system identification, and, on the other hand, modelling and simulation. Given the absence of a comprehensive analysis of the structural dynamics from mid-mounted e-bike drive units, this research establishes a fundamental basis for achieving a deeper knowledge of the system’s dynamics. The presented methodology utilising eMBD enables building an integrated system model of the e-bike drive unit, while revealing the main interactions and valid simplifications. Based on this, various investigations and optimisations can be carried out, contributing to the mitigation of NVH and consequently enhancing the overall riding experience. Section 2 continues to describe the drive unit system and derived methodology.

2. Materials and Methods

Within this section, the aim of the research is specified in a hypothesis, first. This is followed by given insights of the e-bike drive unit’s functionality (see Section 2.2, Section 2.3 and Section 2.4), and the proposed methodology for developing the eMBD model (see Section 2.5, Section 2.6, Section 2.7, Section 2.8 and Section 2.9).

2.1. Hypothesis

Through the utilisation of the eMBD methodology, it is expected that facilitated insights of the e-bike drive unit enable a deeper system comprehension. This, in turn, establishes the foundation for various investigations and optimisations leading to reduction in NVH, which is anticipated to enhance the overall riding experience for e-bike users.

2.2. E-Bike Drive Unit

The referred drive unit is a mid-mounted engine system utilised in 90% of all e-bikes and, thus, the most common system [24]. Figure 2 illustrates such a system.

The mid-mounted engine system amplifies the rider’s torque directly when pedalling. The bicycle chain then transmits the overall torque to the rear wheel, likewise a standard bicycle without a supportive engine. In general, the recently commercialised mid-mounted engine systems, alias e-bike drive units, represent a new use case area. The function of the individual components is principally known, but the configuration and interaction are new. The components of the drive unit and their function in the system are presented next. Figure 3 illustrates the schematic system layout, exemplarily for a two-gear stage.

Table 1. Description of the e-bike drive unit component’s function in the system and its NVH significance.

Designation	Component Function and System Compliance	NVH Significance
Engine freewheel	Decouples the torque of the rotor uni-directional from the pedalling torque of the rider. The decoupling allows the resistance torque created by the gear transmission to avoid being dragged along when pedalling without engine support.	No noises are known due to the clamping mode of operation.
Rider freewheel	Decouples the rider’s pedalling torque uni-directional from the torque of the chain wheel, resulting from the bicycle chain force. It allows the chain to spin while the pedal remains in a fixed position.	There is a clicking sound when the freewheel opens due to the form-fitting mode of action.
Electrical motor (rotor and stator)	Generates a torque to support the rider’s pedalling torque. Requires special controlling of the system, as the timing when the rider is pedalling is indeterminate.	The motor vibrations lead to direct acoustic sound emissions, including the pulse width modulation, which controls the average power delivered by the electrical signal.
Gear stage (shafts, gears and bearings) and chain	Transfers and increases the rotor torque by the gear ratio. Decreases the revolution speed.	The tooth meshing of the gears induces rolling and sliding vibrations. Over-roll frequencies of the bearings are negligible to perceive. The chain braces the system and indicates vibration forces.
Housing	The housing absorbs the rotor dynamic forces transmitted by the bearings. It must resist the rider’s weight, the overall system torque and emerging dynamical forces.	The housing is mounted directly to the e-bike frame without the possibility of further acoustic sealing. Moreover, the stator is assembled directly inside the housing, whose vibrations during operation result in structure-borne acoustic emission.

explains the related functions.

The transfer path of the aforementioned components is presented next, to comprehend the system interactions.

2.3. Mechanical Transfer Path

Figure 4 shows the transfer path in the force direction. It starts at the inner parts of the rotor, shafts, gears, freewheel(s) and bearings, continues at the stator and housing and ends at the e-bike frame section.

The rotor and stator interact through their electromagnetic forces. These forces act on the housing and along the shaft, which ties the rotor. The transmission of the gears transfers the resulting forces of the rotorshaft, the further shafts, and the freewheel(s). Thereby, the bearings transfer the forces towards the housing. The housing reaction forces retroactively impact the inner drive train.

The connection plates transfer the housing forces to the e-bike frame, whereas dynamical reaction forces act retroactively on the housing. The rotor-based torque accumulates with the rider’s pedalling torque, which the chain transfers to the rear wheel. The resulting angular acceleration at the rear wheel, after eliminating all resistance moments and inertia of the system, then leads to an acceleration of the bicycle. Section 2.4 explains the resulting excitation effects, which occur out of the system interrelationship.

2.4. Excitation Effects

The excitation of the system originates from the rotor dynamics. The causative effects are explained in the following. Figure 5 reveals the aforementioned acting forces in detail, applied for a two-gear stage drive unit.

The torque of the electric motor $T_{Rotor,res}$ results from the dynamic air gap forces between the rotor and stator (e.g., $F_{Rotor,nodalforces}$). A torque ripple and radial pulsing forces emerge through the energisation and control, as well as through the non-equidistant mechanical position of the rotor relative to the stator [15,25]. Axial acting forces may occur when the rotor does have an axial offset [26].

The gear forces excite the system due to the tooth meshing and its transmission error, resulting in elastic deformation, misalignment, and dynamic gear mesh stiffness [27]. Caused by the present hexagonal tooth mesh, additional axial forces occur, which bend and axial squeeze the combination of shafts, bearings and freewheel. There can be additional non-linear effects at each coupling point within the system, e.g., due to clearance, displacement-variable stiffness, damping or even friction-based excitation [28].

The engine freewheel’s torsional rigidity depends on the deflection angle related to the transferred torque $T_{Freewheel}$ [29]. This load-dependent stiffness characteristic varies with larger deflection angles, leading to a parametric excitation in the system (comparable to [30]). The overall gear impact of bending and tilting onto the freewheel has yet to be investigated. It is likely that these effects further impair freewheeling and hence impact the system behaviour. The axial acting forces negatively affect the functionality of the freewheel, according to [31]. How these effects affect the drive unit system will be further investigated. The rider freewheel is not further analysed since this freewheel opens when the rider applies no torque, in which case the electric motor also does not provide support. The rider applies forces at both pedals, acting as reaction forces at the crankshaft ($F_{P.,left/right}$) [32]. Torques result around the crankshaft ($T_{P.,left/right}$), in combination with the force application point and the effective pedal length. These pedalling forces may additionally excite the system.

The bicycle chain applies additional forces ($F_{Chain,res}$) to the drive unit system. On the one hand, there is a reaction force due to the chain link and emerging system stiffness to the output drive [33]. On the other hand, there are further acting forces caused by the oscillation behaviour of the chain [34]. Rolling of the chain around the chainring additionally excites the system, as a multiple of the number of chainring teeth. The chain forces are transformed into a torsional torque ($T_{Chain,res}$) through the effective radius of the chainring. Figure 5 does not depict the bearings. The cross frequencies of the rolling elements do not significantly excite the system in the case of a functioning bearing. However, the bearing’s clearance, misalignment and damping cause a variable stiffness and, hence, lead to a parametric excitation in the system comparable to the freewheel.

2.5. Model Development by Deduction and Induction

Our principal approach is derived from the V-model [35], containing surrogate models which are integrated into a holistic system level model (see Figure 6). The applied method allows to investigate sub-systems of the model and to extract and identify the influence of one single parameter variable.

Due to the interaction of all drive unit components with each other (see Section 2.3), the evolved method is a mixture of an induction and a deduction approach. Induction considers a minimal model led by a hypothesis. This model is simplified and contains only utterly necessary information. Its strategy is to conclude towards a general calculation approach, which is valid for the used application. Contrary to this, deduction means starting with an oversized model and reducing it towards a simplified one by evaluating the system parameters and analysing its sensitivities [21,28].

Indeed, the mechanical transfer path and excitation effects were elaborated on before, but the performance of the dynamic behaviour is to be based on the holistic system. Therefore, it is mandatory to model the whole system with all assumptions known, even accepting erroneous approaches in the first induction step. Subsequently, further deduction steps are applied. In each step, specific parameter values are modified or added. The gained results are compared to the previous step. In this way, the impact of the new parameters on the dynamic system behaviour can be identified. In the case that the assumption is correct and verified, the model step becomes the basis for the next loop. This proceeding must be repeated for all relevant parameters. Consequently, a certain number of iterations is essential before the model reaches a sufficient level of accuracy. Though, if the results match the target at any point, the model is designated as validated, and no further loops are required. The final model state may be reached theoretically after the very first induction step, but that is less assumable for real applications.

2.6. Overall Approach

Figure 7 illustrates and specifies the overall approach to develop and validate the simulation model derived from [28]. It mainly contains the process of modelling and simulation, physical testing, parameter comparison and model optimisation. A measurement series alias measuring the campaign or design of the experiment is principally possible but not considered as explained in the following. The illustrated loop indicates the approach mentioned in the previous Section 2.5.

Generally, the drive unit’s design data provide the basis for the modelling and simulation. For this paper, we choose and modify physical samples to correspond to the model state and measure both the sample and the simulation model equally. The resulting values are used to derive qualitative corrections for the simulation model, to correctly map resonance frequencies and orders, representing multiples of the fundamental frequency. A quantitative comparison is necessary to determine deviations due to manufacturing tolerances and temperature impacts, statistically. Nevertheless, this is not considered, as it requires a more extensive measurement campaign with many samples, which are not available. Additionally, the needed capacity and effort would exceed the scope of this contribution. Hence, the model shall predict the phenomenological effects and vibration characteristics without claiming that the values are correct in their amplitude.

2.7. Modelling and Simulation

The modelling and simulation process of the previous Figure 7 is elaborated on in the following Figure 8. The eMBD is necessary to represent the frequency-dependent elastic properties of the system. Additionally, the contacts in the model are required to capture non-linearities and transient processes.

The geometry information is obtained from CAD data. Thereby, only components which are mechanically functioning are selected (compare Section 2.3). The non-bearing structure is neglected since their contribution on the structural dynamics is expected to be insignificant in terms of stiffness. Their mass impact needs to be further investigated but exceeds the scope of this work. The principles of elastic MBD are applied by substituting all rigid parts with the elastic characteristics of substructures and connecting them at the interface nodes via joint elements or stiffness properties. The detailed process of the FEA modelling is given in Section 3.2 and the eMBD modelling in Section 3.3. The gained differential algebraic equation (DAE) of motions is solved by the Hilber–Hughes–Taylor (HHT) solver in MSC Adams. It allows to efficiently calculate fast movements and shocks by explicit methods and slow and stiffer movements by implicit methods [36,37]. It is chosen due to the contact in the gears, where both fast and slow motions can occur [38]. For validation, the simulation results are compared to the measurement results.

2.8. eMBD Model Development and System Identification

Using the following significant simplifications enables step-wise development and verification of the model (see Figure 9), assuming a holistic system model. Thereby, the output data of one step are iteratively compared to those of the next step, through which the impact of one single parameter change can be seen. For example, the stiffness of the freewheel is load-dependent. A constant stiffness parameter is compared to the interpolated stiffness function for the same load points under equal initial conditions.

The step-wise modelling procedure is described below:

(a)

First, the system is simplified as rigid but with gear stiffness.The bearings are implemented as revolute joints, which just allows a rotational degree of freedom (DOF). By this, the associated tooth meshing can be adjusted and optimised for better convergence.

(b)

After that, the engine freewheel’s torsional rigidity is added. The angular system deflection measured at the rotor should correspond to a physical measurement of the system’s twist. The rider’s freewheel is considered rigid since backpedalling complies with a dynamical state without torque transmission. This state is yet out of our scope of interest.

(c)

The flexibility of the shafts is added. As a result, the shafts’ torsional twisting can be verified with an analytical analysis. Through the possible bending of the shafts, the bearing constellation is mapped as an inline and orientation joint, each allowing three rotational DOFs and one bearing an additional axial translational DOF. Thereby, the bending of the rotor-dynamic effects can be analysed without interference by the bearing stiffness.

(d)

Next, the stiffness of the bearings is included. Through the hexagonal tooth meshing, the bearing stiffness may affect the total system behaviour [39]. The influence of the rotor dynamics can be analysed by comparing it to motion variables of the previous step (c).

(e)

Since the previous steps cannot be built within a real test rig, a particular test rig structure is added in both the simulation and the measurement. By this, a first qualitative comparison from the motion variables of the rotor dynamics is possible throughout the measurements (see Section 4.2).

(f)

Finally, the fundamental housing structure substitutes the test rig. By identifying the dynamical rotor dynamic effects within the previous steps, it can be ensured and reversibly determined which resonance comes from the housing structure or the rotor dynamic effects. With the determination of the prior dynamic behaviour, this is possible. This step is also used for validation by comparing the measurement and simulation motion variables (see Section 4.3).

2.9. Verification of the Method

To verify the methodological approach, two major requirements need to be discussed: The first requirement pertains to conducting dynamic-transient calculations in the time domain. The second requirement contains the accuracy through simplifications in the approximation of the relevant physical effects.

Dynamic-transient effects occur through contact in a gear stage, like shocks due to premature tooth meshing caused by load-induced deformation [27]. Moreover, parametric excitations through time-variable stiffness result. This is the case of the tooth meshing (compare [27]), the torsional rigidity of the freewheel (compare [29]) and the load dependency of the bearings. The electromagnetic motor forces are contingent upon the current position of the shaft, which, in its turn, relies on all the preceding effects. These aforementioned effects influence each other, resulting in non-linear system behaviour. Hence, a separate consideration of all features and the linear superposition in the frequency domain calculation is not adequate. For accurate modelling and subsequently achieving the aim of NVH prediction, a holistic system approach requiring dynamic-transient calculation is needed. Calculating these effects within a reasonable computation time and sufficient accuracy is the strength of elastic multi-body-dynamic calculations [40] and, therefore, the tool of choice.

On the component base, it is mandatory to reach a necessary degree of precision with adequate simplification for the shafts, housing, bearings, gears, freewheel(s) and chain. Due to the acting motor forces and transferred gear forces, the shafts twist and bend. This must be taken into account by the elastic shaft properties and cannot be covered by rigid bodies (compare Section 5.1) since occurring misalignment and offset significantly influence the rotor dynamics. The same applies to the housing since the stiffness at the bearing seats and at the clamping intersection have a significant impact on the system behaviour. This approach is the state of the art within FEM and, hence, is seen as given. The used bearing simplification only represents the total bearing stiffness and therefore no bearing orders, which is, however, discussed as sufficient since no condition monitoring is operated, and in a faultless bearing, these orders play no role (compare [6]). The used gear simplification concerns the tooth stiffness, owing to the profile modifications, and contact, which is a decisive base for the transient and non-linear effects. Both the bearings and gears do use numerical viscous damping optimised for convergence and are pre-processed with the internal software of MSC Adams. It is a quasi-static approach providing sufficient precision at the system level, seen as given, as it is state of the art. The engine freewheel’s rigidity is similar to that in the literature (compare [29]) and thus defined as a plausible approach. In contrast, the rider’s freewheel is locked via form fitting when motor torque is applied, which does not lead to additional rigidity, and hence neglect is assumed to be valid. As an initial approach, inserting the chain’s stiffness while disregarding its dynamic behaviour is assumed to be valid (compare [34]).

3. Modelling

This section first analyses the principal rotor dynamics within a simplified torsional system model (see Section 3.1), first. Thereafter, it presents the needed model assumptions to create the holistic eMBD model (see Section 3.3).

3.1. Torsional Oscillator System

For essential investigation of the excitation effects from Section 2.4, the shown rotor dynamic system of Figure 5 is described as a nonlinear oscillator system with time-variable coefficients, derived from [28,41]. The approach is simplified by considering the torsional rotation only. Figure 10 presents the accompanying free-cut system of the rotor dynamic parts from the e-bike drive unit. The non-rotating surrounding structure is neglected.

$\ddot{\varphi }$ is the angular acceleration, $\dot{\varphi }$ the angular velocity, and $\varphi$ the angular deflection. $J_{}$ is the mass moment of inertia from the rotor, gear body or chain wheel. $c_T$ is the torsional shaft rigidity, $c_{FW}$ the freewheel stiffness, and $c_{1,2}$ the gear stiffness in the radial direction. r is the effective gear radius, and $s_g$ the gear clearance. $T_{in/out}$ represents the amplified torque for the input and output. The system is presented as a forced and non-linear oscillator system. With M as the mass matrix, D as the damping matrix, C as the stiffness matrix, and f as the force vector, the arising equation in matrix notation is as follows:

$$M\ddot{\mathrm{\Phi }}+D\dot{\mathrm{\Phi }}+C\mathrm{\Phi }=f$$

(1)

with the motion variables

$${\mathrm{\Phi }}^{\top }=[{\varphi }_1,\cdots ,{\varphi }_8];{\dot{\mathrm{\Phi }}}^{\top }=[\dot{{\varphi }_1},\cdots ,\dot{{\varphi }_8}];{\ddot{\mathrm{\Phi }}}^{\top }=[\ddot{{\varphi }_1},\cdots ,\ddot{{\varphi }_8}].$$

(2)

The mass matrix contains the moments of inertia

$$M=\left[\begin{array}{ccc}{J}_1& & \\ & \ddots & \\ & & J_8\end{array}\right]$$

(3)

Applying viscous damping according to [42] in the gear contact ($d_1,d_2$), the freewheel ($d_3$) and the chain ($d_4$), the damping matrix yields

$$D=\left[\begin{array}{cccccccc}\hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill -d_1{r}_2& \hfill d_1{r}_3& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill d_1{r}_2& \hfill -d_1{r}_3& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill -d_2{r}_4& \hfill d_2{r}_5& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill d_2{r}_4& \hfill -d_2{r}_5-d_3& \hfill d_3& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill d_3& \hfill -d_3& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill -d_4{r}_7& \hfill d_4{r}_8\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill d_4{r}_7& \hfill -d_4{r}_8\end{array}\right]$$

(4)

Assuming zero clearance ($s_{1,2}=0$), the stiffness matrix is

$$C=\left[\begin{array}{cccccccc}\hfill c_{T1}& \hfill -c_{T1}& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill -c_{T1}& \hfill c_{T1}+{\left(r_2\right)}^2·c_1& \hfill -r_2·r_3·c_1& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill -r_3·r_2·c_1& \hfill c_{T2}·{\left(r_3\right)}^2·c_1& \hfill -c_{T2}& \hfill 0& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill -c_{T2}& \hfill c_{T2}+{\left(r_4\right)}^2·c_2& \hfill -r_4·r_5·c_2& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill -r_5·r_4·c_2& \hfill c_{FW}+{\left(r_5\right)}^2·c_2& \hfill -c_{FW}& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill -c_{FW}& \hfill c_{FW}+c_{T3}& \hfill -c_{T3}& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill -c_{T3}& \hfill c_{T3}-{\left(r_7\right)}^2·c_{Ch}& \hfill r_7{r}_8·c_{Ch}\\ \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill 0& \hfill {\left(r_7\right)}^2·c_{Ch}& \hfill -r_7{r}_8·c_{Ch}\end{array}\right]$$

(5)

Assuming no pedal torque as depicted above, the input and output torques are

$$f^{\top }=\left[T_{in},0,0,0,0,0,(T_{Pedal}=0),T_{out}\right]$$

(6)

Presuming a certain constant stiffness $C_0$ allows to determine the undamped and unloaded eigenvalues $\omega$ and eigenmodes by building the characteristic equation according to [42]

$$det(C_0-{\omega }^{2}M)=0$$

(7)

and creating the modal mass matrix ${\mu }_i$ and the modal stiffness matrix ${\gamma }_i$ (${\omega }^2=\frac{\gamma }{\mu }$) with the help of eigenvectors $v_i$ [28]:

$${\mu }_i=v_i^{T}M{v}_i=\sum _{k=1}^n{J}_k{v}_{ki}^2{\gamma }_i=v_i^T{C}_0{v}_i=\sum _{k=2}^n{c}_k{(v_{ki}-v_{k+1,i})}^2$$

(8)

Out of ${\mu }_i$ and ${\gamma }_i$, it is possible to derive the sensitivity coefficients under partial derivation of the inertia and the stiffness [28].

$${\mu }_{ik}=\frac{\partial {\mu }_i}{\partial J_k}=\sum _{k=1}^n{v}_{ki}^2{\gamma }_i=\frac{\partial {\gamma }_i}{\partial c_k}=\sum _{k=2}^n{(v_{ki}-v_{k+1,i})}^2$$

(9)

This allows to analyse the sensitivity of the eigenfrequency resulting from a parameter change, which can be applied for a qualitative comparison as presented in Figure 7.

3.2. FEM Pre-Processing

The numerical methodology to create the pre-processed modal neutral files (MNF) as input for the eMBD model is explained in the following (see Figure 11), revealing information about the mesh, boundary conditions, and computational steps and cost.

For interface modelling, the reference points on which the structural parts are linked in the eMBD model are defined first. These points are the centre positions of the bearings seats, the gears, the freewheel(s), the chain, the rotor and the stator. They are linked to the surface on which the force is acting, with the help of a coupling constraint (see Figure 12). Special cases are those of the rotor and stator, at whose circumferential segments the reference points are placed, on which the electromagnetic force is initiated within the eMBD model.

The meshing is performed by tetrahedral elements allowing to mesh the complex geometries. The following Figure 13 exemplarily illustrates such meshed components, and

Table 2. Mesh information and calculation time from linear perturbation step of processed FEM structures with tetrahedral elements of type C3D10.

Designation	Global Mesh Size	Number of Nodes	Number of Elements	Total CPU Time	Wallclock Time 1
Rotor	≤1.7	57,110	39,029	2.23 × 10${}^2$ s	≈1 min
Rotorshaft	≤0.75	171,936	118,871	1.32 × 10${}^3$ s	≈6 min
Gearwheel	≤1.5	54,842	35,724	8.30 × 10${}^0$ s	≈1 min
Driveshaft	≤2.0	687,359	437,595	7.61 × 10${}^3$ s	≈28 min
Test rig	≤4.0	420,354	285,390	1.82 × 10${}^3$ s	≈8 min
Housing	≤3.9	426,763	249,725	2.29 × 10${}^3$ s	≈16 min

¹ Reference System: CPU @3.6 GHz, 4 cores, 8 logical processors, 32 GB RAM.

reveals the mesh information about all inserted components striking the best balance between accuracy and performance. A convergence study is performed according to [43] by reducing the mesh size by half for each simulation run, while comparing the accuracy and simulation duration.

Table 3. Convergence studies by analysing the mesh size impact towards the solution accuracy and the computational time exemplarily on the rotorshaft. There are no more modes until 20 kHz.

Run	Global Mesh Size	1st Mode	Relative Deviation 1	Total CPU Time	Wallclock Time
1	1.5 s	9539.6 Hz	0.2%	8.5 × 10${}^0$	≈1 min
2	0.75 s	9520.5 Hz	0.07%	1.32 × 10${}^3$	≈ 6 min
3	0.375	9513.3 Hz	-	6.50 × 10${}^4$	≈167 min

${}^1\mathrm{Relative}\mathrm{deviation}(\%)=\left|\frac{\mathrm{Value}\mathrm{run}x-\mathrm{Value}\mathrm{run}x+1}{\mathrm{Value}\mathrm{run}x+1}\right|\times 100\%.$

exemplarily presents this analysis at one component.

The next step is to perform a natural frequency calculation alias modal analysis receiving the modal shapes of each component. This is performed by the Lanczos solver, which provides the most accurate solutions compared to other available solvers within Abaqus [44]. All eigenvalues until 10 kHz are chosen. No boundary conditions or loads are applied, leaving the component’s free vibration. The final step is to generate a substructure. This reduces the problem size through eliminating all degrees of freedom but the retained ones of the reference points [44]. The format of the substructure is chosen to be compatible with its integration into MSC Adams and transferred as the MNF structure.

3.3. Holistic eMBD System

The modelling focus is on the structural dynamics of the e-bike drive unit, including the rotor dynamics. Figure 14 exemplarily shows such a holistic model.

The evolved methodical process to establish the elastic MBD model is explained in the following and illustrated in the respective Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19. The used software for modal reduction of the flexible bodies is SIMULIA Abaqus and MSC Adams for the transient eMBD-simulation of the holistic model. Gear and bearing stiffness are pre-processed with the help of MSC Adams.

The final model assembly contains all mechanically relevant components of the e-bike drive unit, which can be seen as forces acting on the oscillator system similar to Section 3.1: there are rotor and stator forces, gear forces, bearing forces, freewheel forces, and reaction forces caused by shafts and housing.

Rotor and stator forces are calculated in advance by electromagnetic simulations. Thus, they are considered as given. Within the eMBD simulation model, there is a look-up table, which gives the electromagnetic forces for each local node, depending on the relative position from the rotor to the stator. Values in between the discretisation of the table values are interpolated with the help of a spline function. These forces are coupled to a surface section of the structure from the rotor and stator.

The bearings and gears are pre-processed models in which a surrogate stiffness is calculated out of a finite element (FE) model. This stiffness is then applied in the eMBD simulation model and continuously adjusted according to the relative positions. For the bearing, the misalignment and deflection are determined as an additional function of the clearance, the number of rolling elements and the bearing seat, as well as the preload.

In the case of the gear, a contact algorithm is used to determine the distance between a pair of gears by using node points from the pre-processed FE-meshing. These node points lie on the surface of the gear, and their relative position is interpolated. Depending on the pressure angle and the relative position of the two gears, the stiffness results and gear force are output.

Instead of creating a simulation model of the freewheel, the quasi-static torsional stiffness characteristic is measured and pre-processed to access it within the final eMBD model. This characteristic is then interpolated, depending on the difference in the torsional–angular position of the freewheeling components. The resulting torque is the output.

The structural parts, like the shafts and the housing, are pre-processed by FEM. They are modal reduced within a frequency calculation, transferred into sub-structures, and exported/imported as modal neutral files (MNFs) to the eMBD model containing the surrogate stiffness. The momentum of inertia is also considered for all components. Structural damping is applied for the shafts and the housing. Additionally, viscous damping is applied for the relative velocity of the tooth meshing, the bearing motion and the freewheel twisting. The damping terms are fitted for the numerically stable simulation.

The following parameter values are used as a base for simulations (see

Table 4. Used parameter values as reference for input within the eMBD simulation.

Parameter Name	Parameter Value	Related Parts
Solver HHT 1	error $=1.0\times {10}^{-6}$, max step $=5.0\times {10}^{-5}$, $\alpha =-0.3$	Holistic model
Steps size output	$d{t}_{out}=1/14.000$
Duration (1 s Ref)	Calculation time ≤ 2 h
Structural damping 2	2%	Shafts, Gear bodies, Housing
Gear damping	variable (damping $\approx 1.0\times {10}^{-2}$ Ns/mm, $d_{exp}=2$)	Gear contact
Gear stiffness	variable	Gear tooth
Bearing damping	variable (max. $\approx 1.0\times {10}^{-2}$ Ns/mm)	Bearing
Bearing stiffness	variable (max. $2.0\times {10}^{-5}$ N/mm)	Bearing
Freewheel damping	<1.0 Nmms/∘	Engine freewheel
Freewheel stiffness	variable	Engine freewheel
Torque Rotor 3	<5000 Nmm	Holistic model
Speed Rotor (RPM) 3	<3000 1/min	Holistic model

¹ Optimal solver settings after convergence studies. ² Applied to flexible structures. ³ Bosch internal.

), giving a concise overview.

4. Measurement Arrangements

To identify the transient system dynamics, the motion variables of the different components must be measured. Therefore, two setups are chosen, which allow an alignment of the simulation results. The principle approach is explained first.

4.1. Principle Validation Approach

Firstly, the validation criteria are defined (see Figure 20). For the e-bike drive unit, the motion variables, e.g., acceleration and velocity, are chosen. Measured in the time domain, the signal analysis in the time–frequency domain allows to identify transient (short-term occurring) and non-linear effects (e.g., parametric excitations).

Secondly, different measurement setups are chosen, which are replicated by the simulation. The first arrangement is a self-developed test rig assembly, which allows decisively measuring the rotor dynamics inside the drive unit under negligible impact through the housing (see Section 4.2). The second arrangement considers a complete e-bike drive unit to identify the holistic system’s dynamics (see Section 4.3). Since the vibration measurement is performed on the fast rotating components in the first arrangement, contactless measurement techniques are required. State-of-the-art laser vibrometry is used, providing velocity information at one shaft allocation at a time, which then can be examined in the simulation model. Tri-axial accelerometers are employed on the outer structure, such as the housing, facilitating conclusions about the rotor dynamics inside the structure.

Thirdly, to align the data from both the measurement and simulation, the following metrics are used (see

Table 5. Metrics to align measurement and simulation results.

Designation	Initial Conditions	Analysis Method
Steady state	Constant RPM, different load	Power Spectral Density (PSD)
Ramp up	Increasing RPM, different load	PSD Spectrogram
		Order spectrum

):

The frequency peaks from the spectral analysis of the steady-state run or from the order spectrum needs to match at similar frequencies. Missing orders of knowingly modelled effects, e.g., gear mesh, indicate that something is fundamentally wrong in the simulation. In this context, the rpm-dependent effects exert a significant influence, whereas the spectral analysis of the ramp up additionally visualises excited and even non-linear mode shape resonances as horizontal bands. Non-existent or shifted resonances in the simulation results are an indicator that necessary degrees of freedom are missing, or that stiffness and mass as well as damping do not have the same parameter size as the measurement. Generally, the power spectral density (PSD) allows to relate the energy that each frequency has [45] and hence is used for comparison.

At last, based on the alignment metrics, the simulation model is iteratively refined to minimise discrepancies between the measured and simulated motion variables. This iterative process involves, for example, gear mesh adjusting, bearing stiffness adjusting, adapting the numerical damping, adjusting the DOFs on particular parts like the freewheel, replication of more accurate elastic structures in FEM, or capturing and adapting more accurate force interactions, to optimise the model’s performance.

4.2. Measurement of Rotor Dynamics

Since the excitation significantly influences the total vibration behaviour and originates from all rotating components, these are to be examined first. Therefore, an arrangement is chosen which allows an approximated blocked force approach at the bearings seat (see Figure 21). By this procedure, it can be excluded that excitation or resonance comes from the enclosed structure.

The thick plates substitute the drive unit’s housing, no pedalling loads are implemented, and no chain is used. The plates are very stiff compared to the initial drive unit’s housing. As a result, the total stiffness at the bearing seat is significantly affected by the bearing properties themselves, and the surrounding structure’s impact can be neglected. Furthermore, the deflection of the stator during load is prevented chiefly. These conditions allow measuring the dynamical shaft vibrations without significantly affecting the results by structural resonances of the housing. The test rig resonances must be monitored to value and prove the assumption.

Accelerometers capture the test rig vibrations at the bearing seats, a laser Doppler vibrometer (LDV) captures the translational shaft vibrations, and a laser rotation vibrometer captures the rotational vibrations. The measurement position is directly on the shaft’s surface by preparing them with a reflector film. Figure 21c shows these positions, whereas LDV indicates the lasers’ measurement positions, P denotes the accelerometer sensor position, and $x,y,z$ represents the measured direction.

4.3. Measurement of Drive Unit Dynamics

The load case scenario, the chain as an additional drive train component, and the surrounding housing structure, affect the rotor dynamics. Hence, the measurement retains an entire drive unit and uses tri-axial accelerometers, which are closely placed to the outer bearing seat (see Figure 22). One accelerometer is engaged at only one bearing of each shaft since the other positions are unreachable. Through transient effects, the motion variables at each measurement position are almost similar to all other positions. A test setup without an e-bike frame is chosen to eliminate the frame’s structural dynamic impact on the measurement results.

5. Results

First, the results of the simulations derived from the developed model steps from Section 2.8 are presented (see Section 5.1), on the one hand, to exemplarily show the identification and verification process and, on the other hand, to prove that both the measurement and the comparison of the simulated and measured motion variables works. Second, the model assumptions are validated through appropriate cross checking of the measurements and simulations in Section 5.2.

5.1. Verification of Model States

The model variants in

Table 6. Simulations of different states used for verification of the model.

Simulation, Variant	Gears	Shafts	Bearings	Freewheel	Torque
A	flex tooth	rigid	revolute joints	none	const.
B	flex tooth	rigid	revolute joints	rigidity 1	const.
C	flex tooth	flex	joints 2	rigidity	const.
D	flex tooth	flex	flex	rigidity	const.
E	flex tooth	flex	flex	rigidity	EMAG
(F) 3	flex tooth	flex	flex	rigidity	EMAG
(G) 4	flex tooth	flex	flex	rigidity	EMAG

¹ Torsional. ² Inline and orientation. ³ containing the flexible test rig structure. ⁴ containing the flexible housing structure.

are used to point out the impact of increasing the model complexity to approximate the natural vibration behaviour. The chosen sampling resolution is (at least) 14,000 Hz in the simulation and 20,480 Hz in the measurement to capture all relevant frequencies. These resolutions ensure that each tooth mesh is scanned at least 20 times, even at maximum gear speed. The analyses are performed at a steady state since the occurring orders are a multiple of the rotor speed.

The motion variable used for comparison is the rotational speed of the rotorshaft. Due to the transient interrelationships, it should contain the frequency information of all components affected by the rotational motion. Figure 23 presents the different simulation states.

The exemplarily presented analysis is carried out solely up to 600 Hz since the tooth meshing and motor frequencies are lower and thus assessed in the evaluation. All rotational speed-dependent frequencies are present in at least their first order. If the motor force is applied (simulation variant E), the motor orders dominate, and the tooth meshing by the gear becomes less visible. The information about the frequency peaks can be obtained from

Table 7. Comparison of the rotational oscillation from both the simulation and the measurement at similar speed.

Designation	Frequency, Measurement [Hz]	Frequency, Simulation [Hz]	Deviation [%] 4
1. RPM 1	27.7	28.4	2.53
2. RPM	56.6	56.0	−1.06
3. RPM	84.1	84.1	0.00
4. 2 RPM	112.2	113.8	1.43
5. RPM	140.4	139.1	−0.93
1. 1st gear stage (GS)	133.4	134.1	0.53
6. RPM	168.0	166.3	−1.01
7. RPM	195.4	196.5	0.56
8. RPM	224.5	225.0	0.22
9. RPM	252.3	251.6	−0.27
10. RPM	280.5	280.9	0.01
SubH 3 1st GS (+20 Hz)	153.0	155.8	1.83
2. 1st GS	267.2	268.1	0.34
1. 2nd GS	392.4	393.8	0.36
1. Electro-magnetic (EMAG)	561.6	562.3	0.12
SubH 2 EMAG	533.4	536.2	0.52
SubH 3 EMAG	504.4	504.9	0.10
SubH 4 EMAG	477.1	479.5	0.50

¹ Revolutions per minute. ² Additionally SubH 1st GS (−20 Hz). ³ Sub-harmonics. 4 Deviation (%) = $\left|\frac{\mathrm{Simulated}\mathrm{frequency}-\mathrm{Measured}\mathrm{frequency}}{\mathrm{Measured}\mathrm{frequency}}\right|\times 100\%$.

.

The next step is the comparison with measured values. Therefore, the test rig arrangement from Section 4.2 is used under similar initial conditions. It is accepted that there are deviations in the operating parameters, e.g., bearing clearance, unknown gear mesh damping (even though the lubricant is left out for testing), and misalignment. At least the primary orders have to match the simulation and the testing. Figure 24 presents the results.

The spectra of the steady-state analysis of the simulation exhibit the same characteristics as the measurements, specifically about the frequencies at which peak amplitudes occur (see Figure 24). The evaluated frequency peaks are compared in Table 7.

The maximum frequency deviation of the simulation is 2.53% from the particular measured value within the analysed frequency range.

5.2. Validation of Model Assumptions

In order to prove that the model assumptions are correct, a holistic e-bike drive unit is physically measured and simulated. Since measuring inside the drive unit is impossible in the existing configuration, the exemplarily chosen sensor position is $P{1}_z$ (see Figure 22c). The translational velocity at this position is compared from a measurement with a locked rider freewheel replicating the simulation state and a similar simulation, including the flexible housing. The drive unit’s housing structure is connected by special metal sheets to fix it in the measurement, whereas these connections are rigid in the model. Figure 25 shows the Campbell plots of the measured motion variables during a linear ramp up from zero rpm to the maximum rpm at medium torque. The applied model state is the simulation variant G (see Table 6).

An approximate correspondence can be assumed by visually comparing both plots but it is insignificant. Moreover, the simulation has some numerical convergence issues at the beginning, visible as the broadband-like vertical line within the plot.

For more accurate analysis, the average order spectrum from various simulations is plotted and compared to the measurement in Figure 26. This order spectrum analysis enables the comparison of speed-dependent effects arising from tooth meshing, motor excitation, and rotor dynamic vibrations. It is used to prove that the simulation has the same trend as the measurement. In the simulation model, some parameters are varied to exemplarily show the variations influenced by the bearing clearance, freewheel rigidity and housing stiffness. These modifications do not affect the orders directly but the amplitude of the order, which can be seen.

It is noticeable that the order spectrum differs depending on those parameters mentioned above. Generally, the simulation results approximate the trend of the curve characteristics. It becomes clear that the most relevant orders exist in the measurement and the simulation. Consequently, this study concludes that the simulation model includes the principal correct assumption of the system definition and parameters. A design of experiment allowing a quantitative comparison would be needed to validate the model correctly.

6. Discussion

The comparison of the previously depicted spectra from simulation and measurement shows that the principal dynamical impact of most parameters, like electromagnetic forces, gear stiffness, bearing stiffness and freewheel stiffness, is obtained in the model (see Figure 24 and Figure 26). Optimising the stiffness and damping terms would enable better accuracy since the amplitude or sub-harmonics of the motion variables do not exactly match. The mode shapes of the rotor dynamic system in conjunction with the bearing stiffness exist but do not precisely match the measurement results. For example, there are resonance frequencies, visible as horizontal lines, shown in Figure 25, which result from the excitation of divergent mode shapes. Through the parametric excitation caused by the load-dependent stiffness of the gears, bearings, freewheel or even housing, these mode shapes oscillate at an integer multiple of their eigenfrequencies.

Furthermore, this research has yet not questioned using only one physical e-bike drive unit sample for each measurement. Considering deviations caused by manufacturing tolerances from different samples or the repeatability of the identical measurement performed at one sample would provide a statistically grounded statement. Therefore, the results show an acceptable accuracy. A more accurate measurement of the operating vibration modes is still required. However, accomplishing this task, in reality, is not straightforward, primarily due to the necessity of multiple laser Doppler vibrometers (LDVs) and the requirement for accessible locations on the shafts. Matching the operation vibration modes remains an open point for further improvement. The aim is that the amplitudes of the compared motion variables match in the frequency domain and contain the structural resonance effects. Achieving this task requires a design of experiments (DoEs) and sensitivity analysis to quantitatively investigate disturbances through, for example, manufacturing tolerances and temperature, to gain the correct parameters for the model validation.

Consequently, the model parameters represent the operating conditions of the e-bike drive unit with nominal dimensions in the current model state. Adjusting model parameters associated with mechanical properties could enable the indirect consideration of both temperature effects and tolerance influence for a specific and constant state. The motor model does not contain changes in the electromagnetic force by dynamically caused structural effects, which would have a further impact. A high-fidelity motor model and retroactive coupling would enable to map electromagnetic forces impacted by dynamic effects. However, the primary excitation orders of the motor exist in the model, enabling essential analyses of the holistic structural dynamics.

The holistic eMBD model is comparable to the drive-unit assembly of Section 4.3. It contains a housing structure modelled as a merged component in the contact area, which is a significantly simplified assumption. The contact area of the junction points and junction damping impacts this structure and should be further researched. Additionally, the connection of the holistic e-bike drive unit to the frame intersection should contain the flexibility of the connection plates. The test rig structure from Section 4.2 does work, but generally, it is necessary to apply a load. Additionally, through the different mounting situations of the bearing and the assembly, there are tolerance deviations of the bearing seats, which impact the results. Hence, this arrangement is only usable for a principal comparison. The mentioned aspects need to be regarded and prospectively analysed for final validation.

In summary, this research successfully identifies relevant components for developing the eMBD model. It achieves sufficient accuracy in prediction, even though non-load-bearing components are neglected. The caused fault by these assumptions still needs to be checked. The aim is achieved: to predict occurring transient and non-linear effects with the chosen simplifications. However, this simulation exclusively incorporates caused effects of elements that are explicitly included in the model. Hence, it is important to assess whether other effects hold significance under changed parameter values and boundary conditions. Overall, a foundation can be created which serves as a base for further research work, fulfilling the demand to enlighten the e-bike drive unit’s structural dynamics.

7. Conclusions

This paper provides a clear and systematic approach for creating an elastic multi-body-dynamic (eMBD) model, which can be applied to various e-bike drive units and potentially other dynamic systems. It offers valuable insights into the critical model assumptions and simplifications that are necessary for practical modelling, thereby creating an initial foundation for further model development. It makes a valuable contribution to a better understanding of the e-bike drive unit’s structural dynamics by elaborating on the holistic mechanical transfer path and excitations, allowing improved statements of noise, vibration and harshness (NVH) and future acoustic modelling. Through the presented approach of simultaneously identifying and building the model using deduction and induction methods, we can save time and enhance efficiency in modelling. The steps for model verification and validation during the modelling process are outlined, ensuring the model’s accuracy and reliability. An advantage of the methodology is that transient and non-linear dynamic effects can be simulated within a reasonable amount of time. Furthermore, the chosen model assumptions and simplifications do require few data, apart from the geometric data, enabling predictive usage in early product development, thereby making a substantial contribution. This can facilitate the optimisation of various NVH issues, such as avoiding critical resonances by a customised design of the layout, or resolving error patterns that occur, which leads to improved riding performance. One tricky aspect of this work is striking a balance between model simplifications for practicality and the pursuit of higher quantitative accuracy. We acknowledge the need for further quantitative investigations, evaluating disturbances like deviations through tolerances and temperature. While the model simplifications are valid, there is room for improvement in terms of achieving a higher level of accuracy, either through simulation or measurement. These limitations highlight the potential for future research in refining the model.