Stochastic assessment of electric powertrain whining noise under early-stage design uncertainties

Abstract

Despite the advantage of being quieter than traditional internal combustion engine vehicles, electric vehicles are often distinguished by high-frequency tonal components, which can be perceived as unpleasant to the occupants and the environment. To ensure optimal acoustic comfort in electric vehicles, it is important to analyze the NVH behavior of e-powertrains during the early stages of the design process which poses inherent uncertainties, such as varying operating conditions, partial knowledge of design parameters, dispersion in measurement-based data, etc. To effectively address these uncertainties, it is necessary to use fast and comprehensive stochastic models during the design phase. In this work, a probabilistic framework is presented to estimate the electric powertrain’s interior whining noises considering the structure-borne contribution. Hence, two different stochastic metamodels are developed for efficient quantification and propagation of uncertainties from electric motor stage to powertrain mounting system. Multivariate Bayesian regression models help to incorporate prior knowledge on the uncertain parameters and generate the respective posterior distributions using Markov chains Monte Carlo (MCMC) techniques. For this particular application, the data is generated through weakly-coupled multi-physical domains estimated using semi-analytical approaches and combined with measured vehicle transfer functions. Importantly, the validation of each domain is conducted separately to ensure accurate representation. The results obtained from the developed probabilistic framework will aid in the early design stages by guiding engineers in making informed decisions to optimize NVH performance.

Introduction

In the mobility sector, Battery Electric Vehicles (BEVs) have emerged as a promising solution in transitioning toward sustainable future utilizing clean energy, eliminating the need of fossil fuels and reducing greenhouse gas emissions. Moreover, the ever increasing imposition of stringent policies and global pressure to go net zero-emission have propelled vehicle manufacturers to expand their EV fleets at a rapid pace. ¹ As a result, automakers are pivoting their focus on enhancing various attributes of EVs to fulfill customer’s needs and comfort.

Among various characteristics, Noise, Vibration, and Harshness (NVH) is a crucial aspect that significantly influences the customer’s perception of quality and overall image of the vehicle. The absence of an Internal Combustion Engine (ICE) in EVs results in reduced noise levels, but introduces a distinct acoustic signature. The interior noise contribution is characterized by high-frequency components, commonly referred to as whining noise, which are often perceived as intrusive and causes discomfort to the occupants. ² Additionally, the noises that were previously masked by the ICE are now far more audible. From a broader perspective, the major sources of noise in EVs (see Figure 1) can be classified under three categories, namely tire-road interaction noise, aerodynamic wind noise, and e-powertrain noise. The first two sources form the background (or masking) noise and the primary source of tonal noise is from electrified powertrains. Further classification is possible for noises originating from electric powertrains: noises of mechanical origins, such as gear meshing and shaft misalignment, aerodynamic origins from air flow through and over the motor, and electromagnetic origins from electromagnetic interaction, refer to Gieras et al. ³

Figure 1.

Global scheme summarizing the perceived interior cabin noise due to three major sources in a typical BEV.

During vehicle design phase, computer simulations are often used to estimate the vibration response and acoustic levels using detailed full-vehicle structural-acoustic computational (refined 3D finite element) models, which are generally time-consuming. However, during early stages of the design process, when quick assessments of different design alternatives is desired, such detailed 3D design-based simulations are not feasible due to their time-intensive nature. Some other challenges during early-stage BEV design process include the limited (or partial) knowledge of unknown design parameter space, lack of simulated/measured data, and variable driving conditions. The outcome of this is a substantial lack of credibility in the deterministically estimated dynamic responses due to uncertainties originating from multiple sub-systems. Therefore, in order to have more informative predictions, it becomes necessary to quantify such uncertainties early in the design phase and assess the influence of the key design parameters on the resulting outputs. This is achieved using surrogate models (or metamodels) which are simpler approximations of a complex model depicting the functional relationship between the inputs and the corresponding outputs. ⁴

In automotive NVH domain, different metamodels have been employed to achieve minimal engine noise, optimal vehicle mass, learning an aerodynamic wind-noise model, etc. ⁵ Nowadays, such surrogates are commonly needed for performing optimization studies. With respect to the noise generated from electrical machines, neural network based surrogate was built to predict the natural frequencies of the stator. ⁶ Likewise, Mohammadi et al. ⁷ used three neural networks to predict the average torque, the torque ripple, and the sound pressure level (SPL) for multi-objective optimization. In another study by Ibrahim et al., ⁸ multiple surrogate models were compared to predict the acoustic noise. As the electromagnetic excitations in e-machines can be sensitive to even slight variations of the geometric and control parameters of the active magnetic parts, it is important to consider such variations in the output predictions. In this regard, Jeannerot et al. ⁹ conducted a time-consuming FEM-based probabilistic robust optimization of an e-machine to reduce SPL, taking into account the variability of random parameters. Beltrán-Pulido et al., ¹⁰ on the other hand, developed a Gaussian process surrogate model of an e-machine using nonlinear FEM to account for uncertainty in torque, flux linkage, and core loss. Almost all prior studies have only evaluated SPL from e-machine without taking into account the various transfer paths contributing to interior cabin noise. As a result, a “global” perspective of metamodels with uncertain parameters in BEVs global acoustic response was missing. This paper addresses this gap.

Probabilistic modeling provides flexibility to evaluate design alternatives by yielding a probability density function (PDF) of the output responses. To accomplish this, Monte Carlo (MC) simulation techniques have been primarily used by researchers for the probabilistic quantification of uncertainties. In automotive context, for instance, Durand et al. ¹¹ built a non-parametric model to capture the variability in the booming noise prediction through random matrices, ¹² quantified the variability in booming noise and body in-white, ¹³ and recently, Brogna et al.^14,15 used Bayesian approaches with Gibbs sampling to model global vibro-acoustic behaviors. When partial objective knowledge/data about the system is available, the Bayesian approach is particularly useful as it allows the incorporation of prior-knowledge (coming from domain expertise/literature/measurement database, etc) in the form of PDF on the uncertain parameters, making it a reasonable choice to be investigated. The work presented here builds upon the ideas and research established by Brogna et al. ¹⁵

In the Bayesian framework, the learning process involves modifying the initial probability statements/assumptions about the parameters before observing the data to updated or posterior probabilities that combine both the prior-knowledge and the data available. However, the posterior distribution can not always be represented in a closed form and is difficult to compute due to intractability issues. Therefore, sampling methods such as MC are used to sample from the unscaled posterior distribution which approximates the true distribution when the sample size is large enough, refer to Bolstad ¹⁶ and Gelman et al. ¹⁷ Also, it is important to build models where the parameters have physical sense and are interpretable, which in turn allows more control and flexibility over the output responses and allows design loops (from an optimization-problem point of view). Therefore, a parametric surrogate modeling approach is investigated here, characterized by a fixed number of physically pertinent key parameters that do not grow with the size of the input data set. Such parametric models make stronger assumptions on the nature of the data distribution and are generally faster than the non-parametric models which are more flexible but often computationally intractable, refer to Murphy. ¹⁸

The SPL perceived inside a BEV is, indeed, dependent on many parameters, for instance the operating speed of the e-machine, thereby defining the vehicle’s wheel speed, environmental factors (e.g. road conditions), EV body type (e.g. sedan, hatchback), etc. As can be seen from the global scheme in Figure 1, all three major categories of sources of noise are dependent on the operating conditions and when combined, different key performance indicators can be estimated, like the prominence ratio, which results in the prominent tones in the SPL spectra. In a recent article, Prakash et al. ¹⁹ developed Bayesian surrogates to consider such background noises (dashed blue box in Figure 1) using measurement databases. Hence, the focus of this article is on developing surrogate models for the electromagnetic whining noise originating from e-powertrains, as depicted in the dashed red box in Figure 1. As can be seen, whining noise is a combination of airborne and structure-borne contributions. In terms of SPL, it is given by their quadratic average as

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{L}}_p(\omega )=10\mathrm{lo}{\mathrm{g}}_{10}}\\ \multicolumn{1}{c}{\frac{〈{|[{\mathit{H}}^{\mathrm{NTF}}(\omega ){\mathit{Q}}^{\mathrm{S}}(\omega )]\oplus [{\mathit{H}}^{\mathrm{VTF}}(\omega ){\mathit{F}}^{\mathrm{B}}(\omega )]|}^2〉}{p_{\mathrm{ref}}^2},}\end{array}$$

(1)

where ⊕ denotes the complex addition of airborne and structure-borne contributions, ${\mathit{L}}_p$ is the SPL in dB(A) inside the cabin as a complex function of frequency $\omega$ in Hz, ${\mathit{H}}^{\mathrm{NTF}}$ corresponds to the measured noise transfer function (NTF) in Pa/(m³/s²), ${\mathit{Q}}^{\mathrm{S}}$ is the volume acceleration at the source location in m³/s ² , ${\mathit{H}}^{\mathrm{VTF}}$ the measured vibration transfer function (VTF) in Pa/N, ${\mathit{F}}^{\mathrm{B}}$ the excitation force applied on the car-body side in N, and $p_{\mathrm{ref}}$ is the reference sound pressure equal to 20  $\mu$ Pa.

Due to the complexity of global noise assessment involving the interaction between different weakly-coupled physical mechanisms and design parameters, a certain level of assumptions/simplifications is necessary. They are as follows:

1. The work considers an interior PMSM (IPMSM) in which the permanent magnets are embedded in the rotor core. IPMSMs have been widely used in EVs for traction application due to their high torque to inertia and volume ratio and high efficiency. ²⁰ Skewing and the magnetic saturation effects are not taken into account for this early stage study.

2. The mechanical (e.g. gear whine) and aerodynamic noises in the powertrain are neglected and the focus of this article is only on the noise due to electromagnetic origins.

3. Assessment of airborne noise contribution involves the computation of the radiated acoustic power, the equations for which are readily available in the literature. ³ Therefore, this study focuses on assessing the interior SPL considering only the structure-borne path that is, on estimating $\mathit{F}(\omega )$ in equation (1) and its dispersion. Nevertheless, the developed methodology can be applied to airborne noise contribution as well.

The subsequent sections of this paper are structured as follows. Section “Physical problem description and modeling” describes the physical problem and the modeling approach chosen for deterministic evaluation of output responses. The methodology encompassing input parameter sampling, simulated data generation, and the Bayesian-driven metamodeling procedure is explained in Section “Metamodeling for electromagnetic domain.” Within Section “Metamodeling for electric powertrain suspension,” the metamodeling strategy for the electric powertrain suspension is explained, wherein the Bayesian framework harmonizes with a dimensionality reduction technique. The uncertainties from electric motor stage and eletric powertrain suspension stage are coupled together to produce uncertain interior SPL, as described in Section “Interior sound pressure level through coupled stochastic metamodels.” The paper concludes with the validation of each physical block.

Physical problem description and modeling

In this paper, a lower-case character denotes a scalar variable, a lower-case bold character represents a vector, and a matrix is expressed as an upper-case bold character. Later in the paper, random variables are introduced which are expressed as upper-case characters.

Input parameters

The noise perceived inside the cabin due to electric powertrain is assumed to be largely dependent on three sets of parameters:

1. parameters related to the operating conditions (OC) belonging to the set ${\mathcal{P}}_{\mathrm{OC}}\subseteq {\mathbb{R}}^2$ which are typical client usage profiles collected during real driving conditions (from existing similar vehicles). They are represented by the pair $\{\mathrm{\Omega },\tau \}$ , where $\mathrm{\Omega }$ is the speed in RPM and $\tau$ the torque generated by the e-motor in Nm.

2. geometrical parameters of IPMSM such as stator length, rotor outer diameter, etc. belonging to the set ${\mathcal{P}}_{\mathrm{geo}}\subseteq {\mathbb{R}}^{n_{\mathrm{geo}}}$ , where $n_{\mathrm{geo}}$ is the total number of geometrical parameters considered.

3. control parameters that define the information related to the three-phase sinusoidal current excitation, ${\mathcal{P}}_{\mathrm{control}}\subseteq {\mathbb{R}}^2$ . They consist of the pair $\{I,\varphi \}$ where, $I$ is the root-mean-square (RMS) amplitude of the current in A and $\varphi$ is its phase angle in electrical degrees.

Multiphysical nature of electric powertrain noise assessment

In IPMSMs, the input current harmonics, combined with the motor’s geometry and the winding pattern, result in the harmonic distribution of the magnetic flux density in the airgap, thereby exerting electromagnetic (EM) forces on the stator core and frame of the machine. As a consequence, the stator and the frame vibrates at the corresponding excitation frequencies, generating magnetic noise and vibration.

The noise assessment of electric powertrains involves an interplay among different weakly-coupled physical mechanisms, as shown in Figure 2. Generally, the methods used for noise assessment are categorized under numerical, semi-analytical, and analytical methods. The computation of magnetic flux densities in the airgap is performed in the EM domain. In the structural domain, the excitations coming from the EM-domain are first transformed from time-domain to frequency domain and are then mapped onto the structural mesh. This can easily be achieved using commercially available FE solvers as was done in many previous studies.^9,21 Despite being more accurate than analytical or semi-analytical methods, such purely-numerical methods are time consuming and prediction on wide-speed range becomes a challenge. ²⁰ Therefore, a common technique is to resort to semi-analytical methods by calculating the EM force through a simplified 2D-FE model and then obtain vibration and acoustic predictions using analytical approaches.^22,23 Indeed, the magnetic flux densities in the airgap can also be estimated analytically using a magnetomotive force function and a permeance function. ³ However, the precision of such purely-analytical models are usually limited due to their simplification. Therefore, in this study, a more informative nominal model (using FEM) is preferred for computing the EM-domain.

Figure 2.

Overview of the multi-physical mechanism involved in the generation of whining noise in BEVs. NTF and VTF are the measured noise transfer function and vibration transfer function, respectively.

Semi-analytical nominal model for e-motor stage

The vector potential approach is utilized to solve the magnetic problem. The use of the magnetic vector potential in 2D-FEM has a well-established theoretical foundation. ²⁴ A 2D analysis of the IPMSM is generally preferred since only the magnetic flux densities need to be estimated, neglecting thermal and time-dependent effects (which are generally observed in 3D analyses). Also, under magnetostatic assumption, eddy current and temperature dependent effects are not considered. The radial and tangential Maxwell pressure, also referred to as the Airgap Surface Force (AGSF) in N/m² applying on the stator core along the airgap $\delta$ , are given by^25,26:

$$P_{\mathrm{rad}}(R_{\delta },\theta ,t)=-\frac{1}{2{\mu }_0}(B_{\mathrm{rad}}^2(R_{\delta },\theta ,t)-B_{\tan }^2(R_{\delta },\theta ,t)),$$

(2)

$$P_{\tan }(R_{\delta },\theta ,t)=-\frac{1}{2{\mu }_0}\left(B_{\mathrm{rad}}(R_{\delta },\theta ,t)B_{\tan }(R_{\delta },\theta ,t)\right),$$

(3)

where $R_{\delta }$ is the radius at which the pressure components are computed, usually corresponding to the middle of the airgap, $B$ is the magnetic flux density in T, ${\mu }_0$ is the magnetic permeability of vacuum, $\theta \in [0,2\pi )$ is the angular position, and $t$ is time. Taking the Fourier transform in both time and space, one gets ²⁵ :

$$P_{\mathrm{rad}}\left(R_{\delta },\theta ,t\right)=\sum _{r=0,s=-\infty }^{+\infty }{\widehat{P}}_{r,s}^{\mathrm{rad}}\exp \left(r\theta \pm s{\omega }_{e}t+{\varphi }_{r,s}^{\mathrm{rad}}\right),$$

(4)

$$P_{\tan }\left(R_{\delta },\theta ,t\right)=\sum _{r=0,s=-\infty }^{+\infty }{\widehat{P}}_{r,s}^{\tan }\exp \left(r\theta \pm s{\omega }_{e}t+{\varphi }_{r,s}^{\tan }\right).$$

(5)

where $r$ is the spatial order of the force (also referred to as the angular “wavenumber”) with respect to the angular position $\theta$ which determines the periodic shape of the force distribution, $s$ is the temporal order with respect to the mechanical frequency ${\omega }_e$ , and ${\varphi }_{r,s}$ is the phase information associated with each harmonic. The angular velocity of rotor in electrical degree is given by ${\omega }_e=2\pi f_e$ , where $f_e$ is the fundamental stator winding electrical frequency. With $\mathrm{Npp}$ being the number of pole-pairs and $\mathrm{\Omega }$ the rotational speed in RPM, the electrical frequency is given by

$$f_e=\frac{\mathrm{\Omega }}{60}{N}_{\mathrm{pp}},$$

(6)

where $f_{\mathrm{mech}}=\frac{\mathrm{\Omega }}{60}$ is the mechanical frequency of rotation.

In this article, a progressive wave of spatial order $r$ and frequency $f$ is denoted by a pair $(r,f)$ , where $r$ is an integral multiple of the greatest common divisor of $(2N_{\mathrm{pp}},Z_s)$ with $Z_s$ being the number of stator slots, and can be given by^3,27,28

$$r=-2\nu N_{\mathrm{pp}}\pm \nu Z_s.$$

(7)

It corresponds to the slotting effect in PMSMs and mainly occurs at

$$f=s{f}_e,s=2\nu ,$$

(8)

with $\nu \in \{0,1,2,...\}$ . It is also worth noting that the mechanical order, given by $k_{\mathrm{mech}}$ , and the electrical orders are related by the equation

$$k_{\mathrm{mech}}=s{N}_{\mathrm{pp}}.$$

(9)

In the structural domain, the dynamic response of the stator is calculated by the modal frequency response. To estimate the natural frequencies analytically, the stator system is simplified by considering the stator core with teeth, winding, and the frame which are modeled separately. Then, the natural frequencies of the stator system can be approximated as ³

$$f_{\mathrm{mn}}^{\mathrm{stat}}\approx \frac{1}{2\pi }\sqrt{\frac{K_m^{\mathrm{core}}+K_{\mathrm{mn}}^{\mathrm{frame}}+K_m^{\mathrm{winding}}}{M^{\mathrm{core}}+M^{\mathrm{frame}}+M^{\mathrm{winding}}}},$$

(10)

where $m\in \{0,1,2,...,N_m\}$ and $n\in \{1,2,3..,N_n\}$ denote the circumferential and axial nodes of the simplified subsystems, respectively, and $K,M$ represent the calculated stiffness in N/m and mass in kg of each subsystem, respectively. As explained by Le Besnerais, ²⁹ the radial displacements of the stator-system depends on the ratio of the mean yoke radius and mean thickness as well as on the spatial order $r$ . The vibration amplitude is inversely proportional to the fourth power of the spatial order and therefore, low spatial orders that is, $r\in \{0,\pm 8\}$ are of interest and used for further investigation in this study.

Metamodeling for electromagnetic domain

Problem statement

Metamodeling is a process that involves creating a simplified representation of a complex system. In principle, it involves three main steps: first step is the Design of Experiments (DOE), where the most important input features are selected and sampled as per prior-knowledge; second step is to collect the data using measurements or simulations using the sampled input features from the first step; and the final step is the metamodel development by choosing an appropriate functional mapping that relates input features and the output responses using statistical techniques, followed by metamodel validation and exploitation. ³⁰ In the context of this work, these steps are explained in the subsequent subsections.

Sampling of input parameters

The parameters defined by ${\mathcal{P}}_{\mathrm{OC}}$ , ${\mathcal{P}}_{\mathrm{geo}}$ , and ${\mathcal{P}}_{\mathrm{control}}$ are uncertain and their uncertainty can be taken into account by modeling them as random variables. In this study, real-life client driving profiles are made available a priori which is represented as a joint-PDF of $\{\mathrm{\Omega },\tau \}$ . This prior-knowledge is used to consider uncertainty in the OCs by drawing samples from the available joint-PDF. This can be achieved using conditional probability laws and since the problem is not high-dimensional, the required number of samples can be drawn with a simple-inverse transform sampling technique, refer to the Section “Sampling of operating conditions.” Note that the AGSF is primarily determined by the magnetic field distribution within the airgap, which, in turn, is determined by the geometry of the motor, the magnetic properties of the rotor and stator materials, and the current flowing through the stator windings. These factors are not directly dependent on the rotational speed $\mathrm{\Omega }$ , and therefore the magnitude of the AGSF should remain constant under steady-state operating conditions. However, several other factors do depend on $\mathrm{\Omega }$ such as the electrical frequency of the stator current and thereby influencing the frequency at which AGSF acts. Therefore, in this article, a total of $N_{\mathrm{\Omega },\tau }=300$ samples are drawn from the available client-profiles and the kernel density estimated joint-PDF of the drawn samples can be seen in Figure 3.

Figure 3.

Joint distribution of the sampled OCs $\{\mathrm{\Omega },\tau \}$ . The two shaded curves at the edges represent the marginal distributions of the respective variables.

Here, the number of geometrical parameters considered is $n_{\mathrm{geo}}=9$ , which corresponds to the macro-parameters of the stator, the magnets and the frame of the IPMSM under consideration. To incorporate randomness in ${\mathcal{P}}_{\mathrm{geo}}$ , each of geometrical parameters is assumed to follow a reference-distribution. Note that the hyper-parameters (shape and scale) of the reference distribution is chosen in such a way that it reflects some prior knowledge and is broad enough to accommodate various possible designs. The required samples $N_{\mathrm{sp}}$ are then drawn and used as inputs to perform EM computation. In the present work, the reference distribution for each geometrical parameter is taken to be inverse-gamma with support $(0,\infty )$ , which justifies its choice for physical parameters lying in ${\mathbb{R}}^{+}$ . In addition, the right skewed and heavier tailed (as compared to Gaussian distribution) nature of the distribution that is, long tail on the positive side of the distribution, allows to consider the possibility of including extreme positive samples. Other choices are well possible within the same framework, if the user has prior knowledge of the upper bounds. The same argument is applicable for the parameters in ${\mathcal{P}}_{\mathrm{control}}$ ; this can be collectively written as

$${\mathit{x}}_i~\mathrm{InvGamma}(a_i,b_i),\forall \mathit{x}\in {\mathcal{P}}_{\mathrm{geo}}\cup {\mathcal{P}}_{\mathrm{control}},$$

(11)

where, $a$ and b are the hyper-parameters controlling the shape and scale of the distribution. The hyper-parameters are often known based on some prior-knowledge of the analyst (for instance, the allowable tolerance during machining process) depicting the level of uncertainty for a particular design variable or through empirical trial and error methods. The mean $\left({\mu }_i\right)$ and the variance $\left({\sigma }_i^2\right)$ for an inverse-gamma distribution allows the analyst to control the sampling range of random variables and is given by:

$${\mu }_i=\frac{b}{a-1},\forall a>1;{\sigma }_i^2=\frac{b^2}{{(a-1)}^2(a-2)},\forall a>2.$$

A nominal IPMSM is considered with the same parameters as used for the study. ²⁷ The architectural details along with the macro-design parameters are presented in Tables 1 and 2, respectively. The mean of the distribution for each parameter in ${\mathcal{P}}_{\mathrm{geo}}$ and ${\mathcal{P}}_{\mathrm{control}}$ is set to the nominal value and the variance can be assigned as per analyst’s knowledge. For instance, the input parameters can be drawn from a distribution with large variance, if precise information is not available about the uncertainty of the input parameters. As mentioned, under steady-state condition, the AGSF can be computed at a specific $\mathrm{\Omega }$ using the simplified 2D FE model (described in Section “Multiphysical nature of electric powertrain noise assessment”), as AGSF remains constant for all rotational speed samples. For instance, in Figure 4, the radial and tangential components of the AGSF for the nominal design are plotted along the dominant spatial and mechanical orders for $\mathrm{\Omega }=463$  RPM, which gives $f_e=30.87$  Hz and $f_{\mathrm{mech}}=7.72$  Hz. Therefore, as per equation (8), the resulting frequency follows the integral multiples of $f=62$  Hz. In this study, $N_{\mathrm{sp}}=500$ is considered and the EM computations are performed using an open-source electromagnetic solver FEMM ³¹ coupled with Python-based open-source package Pyleecan. ³²

Table 1.

IPMSM architectural details.

Parameters	Value
Number of pairs of poles $\left(N_{\mathrm{pp}}\right)$	4
Number of slots $Z_{\mathrm{s}}$	48
Stator phase number	3
$I$	250 A
$\varphi$	140 °

Table 2.

IPMSM macro-design parameters for the e-machine considered here as an example.

Machine component	Parameter	Value (mm)
Stator	Outer radius $\left(R_{\mathrm{os}}\right)$	134.62
	Inner radius $\left(R_{\mathrm{is}}\right)$	80.95
	Stack length $\left(L_{\mathrm{s}}\right)$	83.82
Frame	Outer radius $\left(R_{\mathrm{of}}\right)$	144
	Inner radius $\left(R_{\mathrm{if}}\right)$	136
	Length $\left(L_{\mathrm{f}}\right)$	150
Magnet	Height $\left(H_{\mathrm{m}}\right)$	6.5
	Length $\left(L_{\mathrm{m}}\right)$	18.9
	Distance $\left(D_{\mathrm{m}}\right)$	14

Figure 4.

Air-gap surface force (N/m²) along dominant spatial and mechanical orders, computed for $\mathrm{\Omega }=463$ RPM: (a) shows the radial component of AGSF $|{\hat{P}}_{r,s}^{\mathrm{rad}}|$ and (b) shows the tangential component of AGSF $|{\hat{P}}_{r,s}^{\tan }|$ .

Simulated data generation for EM forces

For the set of input parameters ${\mathcal{P}}_{\mathrm{geo}}\subseteq {\mathbb{R}}^{N_{\mathrm{sp}}\times n_{\mathrm{geo}}}$ and ${\mathcal{P}}_{\mathrm{control}}\subseteq {\mathbb{R}}^{N_{\mathrm{sp}}\times 2}$ at the first sampled speed, $\mathrm{\Omega }=463$  RPM, the resulting AGSF at the dominant spatial frequencies is given by $\mathcal{Y}\in {\mathbb{R}}^{N_{\mathrm{sp}}\times N_f}$ , where $N_f$ is the number of frequency bins. Since the combined effect of both radial and tangential pressure components is considered, at each spatial order $r\in \{-8,0,+8\}$ , the magnitude of AGSF is considered as

$$\mathrm{AGS}{\mathrm{F}}_{r,s}=\sqrt{|{\hat{P}}_{r,s}^{\mathrm{rad}}{|}^2+|{\hat{P}}_{r,s}^{\tan }{|}^2}.$$

(12)

Figure 5 depicts the AGSF at different spatial orders along the first 10 frequencies. Each frequency bin corresponds to a specific mechanical or electrical order. For simplicity, only the data at $r=0$ (i.e. the data depicted in Figure 5(b)) is considered first for further analysis, but the same approach holds true for other spatial orders as well.

Figure 5.

Air-gap surface force (N/m²) samples generated for different spatial orders at $\mathrm{\Omega }=463$  RPM for three different spatial orders: (a) shows AGSF for $r=-8$ , (b) for $r=0$ , and (c) for $r=+8$ .

Multivariate multi-target parametric metamodel

In the case of stochastic setting, let $X_i,i=1,...,N_p$ denote the random input parameters (from $\{{\mathcal{P}}_{\mathrm{OC}}\cup {\mathcal{P}}_{\mathrm{geo}}\cup {\mathcal{P}}_{\mathrm{control}}\}$ ) and $\mathit{X}={(X_1,X_2,..,X_{N_p})}^{\mathrm{T}}$ be the random vector corresponding to such uncertain inputs defined on the probability space $(\varpi ,\mathcal{A},\mathbb{P})$ with $\varpi$ the underlying sample space, $\mathcal{A}$ the $\sigma$ -algebra, and $\mathbb{P}:\mathcal{A}\mapsto [0,1]$ the probability measure, such that $\mathit{X}:\varpi \mapsto {\mathbb{R}}^{N_p}$ . A realization of $\mathit{X}$ is denoted by $\mathit{x}:=\mathit{X}(\xi )\in {\mathbb{R}}^{N_p},\mathrm{for}\xi \in \varpi$ . It is assumed that the input parameters are a priori statistically independent, guided by two key considerations in design engineering domain. First, due to partial knowledge during early-stages of design, it is highly recommended that the models are parametrized without dependencies. Second, it is easier in practice to model each input parameter separately by determining its marginal distribution $\{p_{X_i},i=1,...,N_p\}$ based on the prior engineering knowledge. ³³ The joint-PDF of the random input vector is then given by

$$p_{\mathit{X}}(\mathit{x})=\underset{i=1}{\overset{N_p}{\Pi }}{p}_{X_i}(x).$$

(13)

With the random input vector $\mathit{X}$ , the simulation model returns random output response $\mathit{Y}={(Y_1,Y_2,...,Y_{N_f})}^{\mathrm{T}}$ and can be written as

$$\mathcal{M}(\mathit{X}):\varpi \mapsto {\mathbb{R}}_{+}^{N_f}.$$

(14)

A realization of $\mathcal{M}(\mathit{X})$ is denoted by $\mathit{Y}:=\mathcal{M}(\mathit{X}(\xi ))$ resulting in a sample of the form $(\mathit{x},\mathit{y})$ . Let $\mathcal{D}\stackrel{\Delta }{=}{\left\{({\mathit{x}}_i,{\mathit{y}}_i)\right\}}_{i=1}^{N_{\mathrm{sp}}}$ denote the data consisting of $N_{\mathrm{sp}}$ samples with the input parameters ${\mathit{x}}_i\in {\mathbb{R}}^{N_p}$ and the corresponding output responses ${\mathit{y}}_i\in {\mathbb{R}}_{+}^{N_f}$ generated from the randomized simulation model as stated in equation (14) . Since $\mathcal{M}$ is computationally expensive, the goal here is to find the unknown forward functional mapping or metamodel

$$\begin{array}{c}\multicolumn{1}{c}{\tilde{\mathcal{M}}:{\mathbb{R}}^{N_p}\mapsto {\mathbb{R}}_{+}^{N_f},}\\ \multicolumn{1}{c}{\mathit{x}\mapsto \mathit{y},}\end{array}$$

(15)

which is a simplified parametric approximation of the simulation model described in equation (14), that predicts $\mathit{y}$ for any future value of $\mathit{x}$ . In this particular case, it becomes a multi-target multivariate regression problem, where the vector valued input predictors are mapped to vector valued output responses. Multi-target regression problem considers each target (or output response, for instance in this work, AGSF value at a particular frequency point corresponds to a single target) as a separate regression problem. This is analogous to obtaining a typical response from multiple-input multiple-output (MIMO) system in structural dynamics. Several different approaches have been proposed in the past for parameterizing $\tilde{\mathcal{M}}(\mathit{X})$ , refer Murphy, ¹⁸ Sudret et al., ³³ and Wang et al. ³⁴ Specifically, in the context of multi-target regression, a broad classification of methods is presented by Borchani et al. ³⁵ From a wide plethora of methods available in the literature, problem transformation method based on single-target method is exploited in this work, where a multi-target model $\tilde{\mathcal{M}}$ is transformed into several single-target problems. Since the problem is not high-dimensional, polynomial basis functions can be used to represent the observed output responses.

For the given $\mathcal{D}$ , the multi-target model $\tilde{\mathcal{M}}$ comprises of $N_f$ single-target models given by

$${\tilde{\mathcal{M}}}_k:{\mathbb{R}}^{N_p}\mapsto {\mathbb{R}}_{+},k\in \{1,2,..,N_f\}$$

(16)

where each model ${\tilde{\mathcal{M}}}_k$ is trained on ${\mathcal{D}}_k\stackrel{\Delta }{=}\{({\mathit{x}}_i,y_{i(k)}){\}}_{i=1}^{N_{\mathrm{sp}}}$ to predict the value of the single target variable $Y_k$ . Considering polynomial basis functions, each output variable can then be written in a generic form as ³⁶

$$Y_k={\tilde{\mathcal{M}}}_k(\mathit{X},\mathit{q})+{ϵ}_k,$$

(17)

$${\tilde{\mathcal{M}}}_k(\mathit{X},\mathit{\theta })={\theta }_{0,k}+\sum _{j=1}^{N_p}\sum _{l=1}^{N_l}{X}_j^l{\theta }_{j,l,k},$$

(18)

where $N_l$ is the degree of the polynomial, ${\theta }_{0,k}$ the intercept term, ${\theta }_{j,l,k}$ represent the unknown coefficients, and ${ϵ}_k$ the fitting error consisting of modeling error and is assumed to be a zero mean Gaussian noise with variance ${\sigma }_y^2$ . In this work, polynomial with degree $N_l=1$ is considered for further analysis as discussed in subsequent sections.

Deterministic check for multi-target regression model

Before proceeding with the stochastic analysis using Bayesian approach, a fundamental practice is to perform the deterministic model evaluation to check the generalization error of the model developed. As discussed Prakash et al., ¹⁹ $k$ -fold cross-validation (CV) technique can be used to validate the model with sufficient accuracy. Note that, the number of folds depends on the amount of data available, so that each training set is representative of the entire available dataset. The CV process is repeated multiple times to consider the shuffling of the original dataset which allows the model to learn from a more representative data sample in the CV-training batch which leads to better performance and generalization property of the model. Considering the available simulated dataset $\mathcal{D}$ and formulating a linear regression problem, as shown in equation (14), the mean value of the coefficient of determination -score is close to $94\mathrm{\%}$ and the root-mean-squared-error (RMSE) $\approx 0.6$ . Some other metrics, for instance, the relative-RMSE, can also be evaluated. ³⁵

Once the deterministic model is validated with sufficient accuracy, stochasticity in the model can be induced using Bayesian modeling approach, which is described in the next subsection.

Bayesian hierarchical model for forward analysis

Bayesian networks are probabilistic directed acyclic graphical models used to represent the relationships between random variables. It consists of nodes representing random variables and directed edges between them representing the probabilistic dependencies between those variables. ¹⁸ The main goal of performing the Bayesian analysis is to combine the prior knowledge with the available data to quantify the uncertainties in the output responses. In principle, the wide-range of input parameters used to generate the simulated data through model $\mathcal{M}$ (i.e. the simplified 2D-FE model) for the training process (using reference distribution) can be different from the ones that would potentially be generated by the analyst intending to use the metamodel for answering particular questions. In this regard, a two-stage Bayesian network is developed in which the first-stage represents the classical metamodel training stage with reference distributions along with Bayesian model evaluation checks, and the second-stage deals in generating the posterior samples from the learned uncertain parameters, thereby propagating the uncertainties to other domains.

Metamodel training stage

The available simulated data $\mathcal{D}$ is first divided into training and validation datasets. Consider ${\mathcal{D}}^{†}$ to be the training data with ${\mathcal{X}}^{†}\in {\mathbb{R}}^{N_{\mathrm{sp}}^{†}\times N_p}$ the input training matrix of predictor variables, ${\mathcal{Y}}^{†}\in {\mathbb{R}}^{N_{\mathrm{sp}}^{†}\times N_f}$ the corresponding output training matrix, and $N_{\mathrm{sp}}^{†}$ the training samples. Following equation (18), the multi-target model for the output training matrix can be written as

$${\mathcal{Y}}^{†}=\tilde{\mathcal{M}}\left({\tilde{\mathcal{X}}}^{†},\mathrm{\Theta }\right)+E,$$

(19)

$$\tilde{\mathcal{M}}({\mathcal{X}}^{†},\mathrm{\Theta })=\mathbb{X}\mathrm{\Theta },$$

(20)

where, $\mathbb{X}=\left[1{\mathcal{X}}^{†}{\mathcal{X}}^{†\odot 2}\dots {\mathcal{X}}^{†\odot N_l}\right]\in {ℝ}^{N_{\mathrm{sp}}\times (N_l{N}_p+1)}$ is the augmented predictor matrix, $\mathit{Q}=[{\mathit{\theta }}_0{\mathit{\theta }}_1{\mathit{\theta }}_2...{\mathit{\theta }}_{{\mathbf{N}}_{\mathbf{l}}}]\in {\mathbb{R}}^{(N_l{N}_p+1)\times N_f}$ is the matrix containing the unknown coefficients, and $E\in {\mathbb{R}}^{N_{\mathrm{sp}}^{†}\times N_f}$ is the matrix of errors.

Let the “prior” PDF of the parameters $\mathrm{\Theta }$ be given by $p(\mathrm{\Theta })$ , $f({\mathcal{Y}}^{†}|\mathrm{\Theta },{\mathcal{X}}^{†})$ denote the “likelihood” of observing the data given the parameters, and $p(\mathrm{\Theta }|{\mathcal{D}}^{†})$ denote the “full posterior distribution” of the parameters conditional on the observed (simulated) training data ${\mathcal{D}}^{†}$ . Using Bayes’ theorem, the posterior distribution $p(\mathrm{\Theta }|{\mathcal{D}}^{†})$ can be approximated from the unscaled distribution given by^16,17

$$p(\mathrm{\Theta }|{\mathcal{D}}^{†})\propto f({\mathcal{Y}}^{†}|\mathrm{\Theta },{\mathcal{X}}^{†})p(\mathrm{\Theta }).$$

(21)

Markov Chain Monte Carlo (MCMC) methods are used to sample from this unscaled distribution which approximates the true distribution provided that the sample size is large enough.

During the training stage, the following steps are followed which are also described as a hierarchical model in Figure 7.

• 1. Firstly, the training data is log-transformed in order to reduce the skewness and then it is standardized as shown in Figure 6. It is assumed that this transformed training data is distributed as per the Normal distribution with mean given by the model shown in equation (20), and variance ${\sigma }_y^2$ . Indeed, even after transformation, the data shows slight skewness (see Figure 6(b)) and could be modeled using a skew-Normal distribution, which is a generalization of the Normal distribution that allows for non-zero skewness. However, for the sake of simplicity and without the loss of generality, the normality assumption is considered in this work. Furthermore, several normality checks exist in the literature that can be evaluated and necessary transformations can be applied. ³⁷

• 2. Formulate the function $\tilde{\mathcal{M}}({\mathcal{X}}^{†},\mathrm{\Theta })$ as a linear model with $N_l=1$ in equation (20). It is the vector-valued mean of the model output given the parameters and the training data ${\mathcal{X}}^{†}$ . Each unknown coefficient is modeled as a random variable following a prior-PDF characterized by its own hyper-parameters controlling the shape of the distribution. In this study, the unknown model coefficients (the intercept term is ignored as the data are centered before performing the analysis) are assumed to follow a Normal distribution, where the hyper-parameters are denoted by $\{a,b\}$ , representing the mean and variance of each distribution. Hyper-priors are also considered, for instance, ${\mu }_{\alpha }$ in Figure 7, depicting the level of uncertainty on the mean of the hyper-parameters. Likelihood and the prior-distributions are written as

$${\mathcal{Y}}^{†}|\mathrm{\Theta },{\mathcal{X}}^{†}\sim \mathcal{N}\left(\widetilde{ℳ}({\mathcal{X}}^{†},\mathrm{\Theta }),{\sigma }_y^2\right),$$

(22)

$${\theta }_{j,k}~\mathcal{N}({\mu }_{\theta },b_{\theta }),$$

(23)

$$\begin{array}{c}{\mu }_{\theta }~\mathcal{N}(a_{{\mu }_{\theta }},b_{{\mu }_{\theta }}),\hfill \\ \forall j\in \{1,..,N_p\},k\in \{1,..,N_f\}\hfill \end{array}$$

(24)

• 3. The model considers heteroscedasticity (Hetero-scedasticity refers to the property of a dataset where the dispersion of the dependent variable (output response) around its mean is not constant along the space of independent variables.) in the observations and therefore, the variance is distributed accordingly along the input frequency range and follows the inverse-gamma distribution:

$${\sigma }_{y,k}^2~\mathrm{InvGamma}(a_{\sigma },b_{\sigma }),\forall k\in \{1,..,N_f\}$$

(25)

• 4. Compute the posterior distribution of each random parameter using Hamiltonian Monte Carlo approach, specifically NUTS algorithm, ³⁸ modeled using an open-source library PyMC3. ³⁹ A total of four different MCMC chains were simulated resulting in 10,000 samples from which first 500 samples from each chain were discarded (commonly referred to as the “burn-in” process) to let the chain converge to its stationary distribution. The marginalized posterior distribution of each random parameter in the model is shown in Figure 8. Due to the high-dimensionality and random nature of the coefficient matrix $\mathit{\theta }\mathbf{\in }{\mathbb{R}}^{\mathbf{(}{\mathbf{N}}_{\mathbf{p}}\mathbf{\times }{\mathbf{N}}_{\mathbf{f}}\mathbf{)}}$ , where each variable represents a random parameter, presenting all PDFs on a single plot may lead to potential confusion. Consequently, in order to enhance comprehension, the box-plots of the PDFs associated with the predictor variables are considered at each frequency point, as illustrated in Figure 8(b). This also provides an estimate on the degree of influence each parameter has on the output at each frequency point.

Figure 6.

Data transformation to better represent the assumption of Bayesian modeling approach: (a) shows the occurrences of the original output training data (AGSF) and(b) the occurrences of the output training data after their standardized transformation.

Figure 7.

Bayesian hierarchical model. Key: rectangular boxes for deterministic data, random variables are shown in circles, diamonds represent the hyper-parameters of prior-distributions, and rounded boxes represent the replication of variables following the plate notation.

Figure 8.

Posterior distributions of random variables in the Bayesian hierarchical model: (a) depicts the posterior distribution of 1000 posterior samples of the random variables ${\sigma }_y^2$ and ${\mu }_{\theta }$ and (b) shows the box-and-whiskers plot of posterior samples of all predictor variables at each frequency point (only first eight dominant targets are shown) contained in the random variable $\mathit{\theta }$ .

Bayesian metamodel validation

The trained Bayesian model must be validated on the basis of the approximated posterior distribution of each unknown random parameter and its predictive capability on the validation dataset given by ${\mathcal{D}}^{\mathrm{‡}}\triangleq \left({\mathcal{X}}^{\mathrm{‡}},{\mathcal{Y}}^{\mathrm{‡}}\right)$ such that $\mathcal{D}={\mathcal{D}}^{†}\cup {\mathcal{D}}^{†}$ , with ${\mathcal{X}}^{†},{\mathcal{Y}}^{†}$ being the input/output validation matrices, respectively. There are several numerical and visual tools available in the literature to diagnose the MCMC sampling process (refer to Gelman et al. ¹⁷ ). For the current developed model, the Gelman-Rubin statistic ¹⁷ $\hat{R}=1.0$ , indicating that all chains are sampling from the same posterior distribution with the effective sample size ⁴⁰ of more than 5000 (out of 10,000 samples) for each random variable, thereby indicating low level of autocorrelation among the posterior samples. Model validation by considering the uncertainty about the parameter values as well as the data generating process is also performed using posterior predictive checks, where the idea is to draw the samples from the posterior distribution of parameters, generate the replicated data (with ${\mathcal{X}}^{†}$ ) and compare it with the observed data. The posterior predictive distribution of the replicated data $\left({\mathcal{Y}}^{\mathrm{rep}}\right)$ is given by

$$p({\mathcal{Y}}^{\mathrm{rep}}|{\mathcal{D}}^{†})=\int p({\mathcal{Y}}^{\mathrm{rep}}|\mathrm{\Theta },{\mathcal{X}}^{†})p(\mathrm{\Theta }|{\mathcal{D}}^{†})\mathrm{d}\mathrm{\Theta }.$$

(26)

The distribution of replicated data is compared against the distribution of original data, both for the overall marginalized distribution and also for the marginalized distribution at each frequency bin (or target), as shown in Figure 9. It is noticed that the model is able to capture the pattern in the training data. However, when the data distribution exhibits bimodality, for instance at $0$ and $185$  Hz, such a simple model may not be able to capture different modes. In such cases, mixture models ⁴¹ are more suitable. Nevertheless, the assumed model is sufficiently accurate for the application in this work. Indeed, Bayesian $p$ -values can also be used by specifying a test statistic (e.g. difference between true and predicted values) to check if the replicated data is more extreme than the observed data. ¹⁷

Figure 9.

Posterior predictive check: (a) shows the marginalized distributions of the overall training data and the overall replicated data and (b) shows the marginalized distributions of the training data and the replicated data at each frequency bin.

Exploitation of Bayesian metamodel

Conditioning on the observed training data ${\mathcal{D}}^{†}$ , Bayesian approach results in the posterior distribution of unknown random parameters in the model $p(\mathrm{\Theta }|{\mathcal{D}}^{†})$ , which is then used to determine the posterior predictive distribution (a reformulation of equation (26)):

$$p(\tilde{\mathcal{Y}}|\tilde{\mathcal{X}},{\mathcal{D}}^{†})=\int p(\tilde{\mathcal{Y}}|\mathrm{\Theta },\tilde{\mathcal{X}})p(\mathrm{\Theta }|{\mathcal{D}}^{†})\mathrm{d}\mathrm{\Theta },$$

(27)

where, $\tilde{\mathcal{X}}$ is the newly generated input parameters (depending on the analyst’s prior-knowledge about the physical parameters) and $\tilde{\mathcal{Y}}$ are the corresponding predicted samples. This exploitation stage of the metamodel is shown as a rather “light” Bayesian scheme in Figure 10, where the hyper-parameters describing the uncertainty on the physical parameters are denoted by ${\tilde{a}}_{\gamma },{\tilde{b}}_{\gamma },\forall \gamma \in \{1,2,...,N_p\}$ . As an example, considering the validation data set ${\mathcal{D}}^{†}$ , and letting, then ${\tilde{N}}_{\mathrm{sp}}$ samples can be drawn from the posterior predictive distribution as shown in Figure 11. The uncertainty estimates using box-and-whiskers are plotted and it is clear that the inter-quartile range varies across the spatial frequencies and the prediction follow the pattern, for instance high AGSF at $0$  Hz for $r=0$ , as observed in the labeled validation data ${\mathcal{Y}}^{†}$ (shown as blue dots in Figure 11).

Figure 10.

Scheme depicting the exploitation phase of Bayesian metamodel (key to identify different elements and relationships remains same as used in Figure 7).

Figure 11.

Posterior predictive samples for the validation dataset ${\mathcal{D}}^{†}$ at $r=0$ . Blue dots are the AGSF values retained in ${\mathcal{Y}}^{†}$ . Box and whiskers represent the extent of samples with quantiles. The axes shown on top of the figure depict the dominant orders (both mechanical and electrical) as a pair $(r,s{f}_e)$ or $(r,k_{\mathrm{mech}}{f}_{\mathrm{mech}})$ at the corresponding frequency.

On the basis of these posterior predictive samples $\tilde{\mathcal{Y}}$ and the uncertain input parameters generated in $\tilde{\mathcal{X}}$ , the uncertainties can be propagated to the dynamic displacement of the machine under consideration, which is discussed in the following subsection.

Propagation of uncertainties in electric motor stage

The dynamic vibration displacement of the IPMSM surface is directly dependent on the harmonic characteristics of magnetic excitation and the natural frequencies and mode shapes of the stator-system (refer to equation (10)). The uncertainties in the magnetic excitation is captured through the Bayesian metamodel developed in the previous subsection. The uncertain geometrical parameters ${\mathcal{P}}_{\mathrm{geo}}$ directly influences the natural frequency of the system. Since a semi-analytical approach is exploited here, a simplified homogenized stator-system is assumed as a finite length cylindrical shell clamped at both ends and then the root-mean-square (RMS) vibration displacement at frequency $\omega =2\pi f$ can be estimated analytically as

$${\mathcal{U}}_k^{\mathrm{RMS}}(\omega )=\sqrt{\frac{1}{3}\sum _i(|{\mathcal{U}}_{i,k}(\omega {)|}^2)},\forall k\in N_f,$$

(28)

where ${\mathcal{U}}_k^{\mathrm{RMS}}\in {\mathbb{R}}_{+}$ is the RMS vibration displacement at each mechanical order indexed by $k$ , $N_f$ is the number of mechanical orders considered, and ${\mathcal{U}}_{i,k}(\omega )$ is the vibration displacement in axial, circumferential, and radial directions, indexed by $i\in \{u,v,w\}$ , respectively.^3,42 For $N_{\mathrm{\Omega }}$ sampled speeds, ${\mathcal{U}}^{\mathrm{RMS}}(\omega )$ is stacked in a matrix ${\mathbf{U}}^{\mathrm{RMS}}(\omega )\in {\mathbb{R}}_{+}^{N_{\mathrm{\Omega }}\times N_f}$ . These vibration displacements are transmitted to the car-body through powertrain mounting system as described in the next section.

Metamodeling for electric powertrain suspension

Problem statement

In order to mitigate the vibrational energy transmitting from the e-powertrain to the vehicle body side, powertrain mounts are generally installed between them. The primary function of the engine mounts is to reduce the transmission of dynamic forces at lower frequencies and provide a static support to the powertrain weight. In the literature, the analysis of vibration characteristics of powertrain mount systems focused mostly on internal combustion engines (ICE). ⁴³ However, the vibrations transmitted from e-powertrain is quite different from the classical ICE powered drivetrains as e-motors have lower inertia, vibrates at higher frequencies and lower amplitudes along with the reverse torque. Nevertheless, globally for e-powertrains, some design considerations regarding the mounting architecture remains the same with certain modifications to adapt the electrification. Different models have been realized in the literature to study the dynamic behavior of e-powertrain mount resonances with its dynamic mass participation and with different mounting architectures.^44–46

As the scope of this paper is rather global, a few degrees of freedom (dof) linearized model is sufficient to model the forces transmitting from the e-powertrain to the mounts on the vehicle body. In this context, an internal Stellantis tool has been used to estimate these forces by considering a three-point suspension architecture, consisting of three mounts each with three translational dofs ( $x,y,z$ directions), thereby accounting for a total of nine input dofs, hereafter denoted as $N_d$ . As shown in Figure 12, for this typical three-point architecture, the vibrational energy from the surface of the e-motor is first transmitted to the powertrain mounts $\{M_i:i\in \{1,2,3\}\}$ through the respective lumped stiffnesses given by $\{{\overline{K}}_i:i\in \{1,2,3\}\}$ and then to the car-body interface through the mount stiffnesses $\{K_i,i\in \{1,2,3\}\}$ . Let $\{{\overline{K}}_i=1:\forall i\in \{1,2,3\}\}$ to consider zero-loss transmission from e-motor to its respective mounts and assume a unit excitation being imposed on each dof. The model acts as a transfer function for each dof, whose output is to be multiplied by the true vibration displacements computed (semi-analytically) using equation (28). Note that the simplified three-point suspension architecture considers independent transfer functions at each dof. In other words, the forces acting on the car-body through a particular dof, say for $M_1$ in $x$ -direction, is dependent only on the displacements imposed on this dof. This simplification of uncoupled mount responses is justified as the same RMS displacements are acting on each dof, which in turn requires only to learn the transfer function metamodel for each dof. The validation of this simplified mechanism is presented later in Section “Validation” with respect to the measured data. If needed, consideration of more advanced coupled systems are well possible within the proposed framework without altering its principles.

Figure 12.

A simplified schematic of a typical three-point e-powertrain suspension architecture showing three mounts ( $M_1$ : Left mount, $M_2$ : Right mount, and $M_3$ : Cradle point).

Let ${\mathrm{H}}_{d,k}^{\mathrm{MB}}(\omega )\in {ℂ}^{N_{\mathrm{\Omega }}\times 1},\forall d\in \left\{1,2,\dots ,N_d\right\},\forall k\in \left\{1,2,\dots ,N_f\right\}$ be a frequency dependent transfer function that maps the RMS vibration displacements given by ${\mathit{U}}_k^{\mathrm{RMS}}(\omega ):={\mathit{U}}_{[:,k]}^{\mathrm{RMS}}\in {\mathbb{R}}_{+}^{N_{\mathrm{\Omega }}\times 1}$ to the output car-body forces ${\mathit{F}}_{d,k}^{\mathrm{B}}(\omega )\in {\mathbb{C}}^{N_{\mathrm{\Omega }}\times 1}$ for each dof and each spatial order such that

$${\mathit{F}}_{d,k}^{\mathrm{B}}(\omega )=|{\mathit{H}}_{d,k}^{\mathrm{MB}}(\omega )|\odot {\mathit{U}}_k^{\mathrm{RMS}}(\omega ),$$

(29)

where ⊙ represents the Hadamard product for element- wise multiplication.

To account for the uncertainties in the mounting behavior through the transfer function model ${\mathit{H}}_{d,k}^{\mathrm{MB}}(\omega )$ , the spatial coordinates and stiffness parameters of the rubber mounts $M_i$ are assumed random and distributed as per their respective reference distributions as

$${\mathrm{c}}_d~\mathcal{N}(g_d^{\mathrm{c}},h_d^{\mathrm{c}}),$$

(30)

$${\mathrm{k}}_d~\mathrm{InvGamma}(g_d^{\mathrm{k}},h_d^{\mathrm{k}}),$$

(31)

$$\forall d\in \{1,...,N_d\},$$

(32)

where “c” and “k” represents the position coordinate and the stiffness parameter of the rubber mount, respectively characterized with their hyper-parameters. Then, Latin hypercube sampling, ⁴⁷ a widely used stratified sampling technique, is employed to sample from these distributions once the values for the hyperparameters are decided.

A similar Bayesian framework is applied as the one used for training and validating the electromagnetic forces in e-motor stage. However, the output space for EM-domain depends on the number of mechanical orders $N_f$ , which is usually low while dealing with electromagnetic NVH issues. In the case of e-powertrain suspension, the dimensions of the output space is directly proportional to the number of sampled speed profiles $N_{\mathrm{\Omega }}$ in ${\mathcal{P}}_{\mathrm{OC}}$ , which in turn defines the frequency resolution of the transfer function spectra. Hence, to efficiently compute the training phase, data reduction techniques are needed to be applied to the response; then Bayesian training will be performed on the transformed components.

Bayesian principal component regression with data reduction

Bayesian principal component regression (BPCR) essentially involves two steps: First, the principal component analysis (PCA) ⁴⁸ is performed on the output training matrix to obtain the orthogonal principal components (PCs). Then, a small subset of PCs containing most of the “information” is chosen (resulting in dimensionality reduction) to perform Bayesian regression to predict the response variable. To determine the number of PCs to retain, multiple methods have been proposed in the literature and grouped under subjective (e.g. scree plots), distribution-based (e.g. Barlett’s test), and computational procedures (e.g. cross-validation). ⁴⁹ In this article, percentage of total variance explained by a subset $\mathcal{S}$ of PCs is used as a measure to retain the number of components.

Let the complex-valued output training data (transfer function) be defined as ${\mathbb{Y}}_d^{†},\forall d\in \{1,2,..,N_d\}$ . Applying PCA and retaining $q=|\mathcal{S}|<N_{\mathrm{\Omega }}$ PCs such that

$${\mathbb{Y}}_d^{†}\in {\mathbb{C}}^{N_{\mathrm{sp}}^{†}\times N_{\mathrm{\Omega }}}\mapsto {\overline{\mathbb{Y}}}_d^{†}\in {\mathbb{C}}^{N_{\mathrm{sp}}^{†}\times q},$$

(33)

where ${\overline{\mathbb{Y}}}_d^{†}$ denotes the column-stacking of PCs and $N_{\mathrm{sp}}^{†}$ the total number of training samples. For a given PC, the proportion of total variance it accounts for is given by:

$${\mathcal{C}}_j=\frac{{\lambda }_j}{\sum _{i=1}^q{\lambda }_i},{\lambda }_1\le {\lambda }_2\le {\lambda }_3\le \dots \le {\lambda }_{N_{\mathrm{\Omega }}},$$

(34)

where ${\lambda }_j$ is the corresponding eigenvalue. The percentage of total variance explained by the PCs in $\mathcal{S}$ is given by, $\sum _{j\in \mathcal{S}}{\mathcal{C}}_j\times 100\%$ and $q$ is selected as per

$$q=\underset{|\mathcal{S}|}{\arg min}|\sum _{j\in \mathcal{S}}{\mathcal{C}}_j-V_{\mathrm{thresh}}|,$$

(35)

where $V_{\mathrm{thresh}}$ is the specified threshold value for the total variance explained by the subset.

The same approach is followed for the metamodel training as presented in Section “Metamodel training stage” with the difference that, it is now trained on PCs. The Bayesian model for each dof reads

$${\overline{\mathbb{Y}}}_d^{†}|\overline{\mathrm{\Theta }},{\overline{\mathbb{X}}}^{†}~\mathcal{N}(\tilde{\mathbb{M}}({\overline{\mathbb{X}}}^{†},\overline{\mathrm{\Theta }}),{\overline{\mathit{\sigma }}}_y^2),$$

(36)

$${\overline{\theta }}_{0,k}~\mathcal{N}({\overline{a}}_0,{\overline{b}}_0),$$

(37)

$${\overline{\theta }}_{i,\nu ,\kappa }~\mathcal{N}({\overline{a}}_i,{\overline{b}}_i),$$

(38)

$${\overline{\sigma }}_{y,\kappa }^2~\mathrm{InvGamma}({\overline{a}}_{\sigma },{\overline{b}}_{\sigma }),$$

(39)

$$\forall i\in \{1,..,N_l\},\nu \in \{1,..,2N_d\},\kappa \in \{1,..,q\},$$

where ${\overline{\mathbb{X}}}^{†}\in {\mathbb{R}}^{N_{\mathrm{sp}}^{†}\times 2N_d}$ is the predictor matrix (different from e-motor stage metamodel), denotes that the parameters involved in this training stage is related to the PCs and not to the real training data, and the surrogate function is formulated using the same basis expansion as used in Section “Metamodel training stage,”

$$\tilde{\mathbb{M}}({\overline{\mathbb{X}}}^{†},\overline{\mathrm{\Theta }})={\widehat{\overline{\mathbb{X}}}}^{†}\overline{\mathrm{\Theta }},$$

(40)

where the augmented predictor matrix is given by ${\widehat{\overline{\mathbb{X}}}}^{†}\in {\mathbb{R}}^{N_{\mathrm{sp}}^{†}\times (2N_l{N}_d+1)}$ and $\overline{\mathrm{\Theta }}=\left[{\overline{\theta }}_0{\overline{\theta }}_1{\overline{\theta }}_2\dots {\overline{\theta }}_{N_l}\right]\in {ℝ}^{(2N_l{N}_d+1)\times q}$ is the matrix containing the unknown coefficients.

The predictive distribution of the transformed posterior samples is as per equation (26). However, these samples are needed to be transformed back to obtain the true predictions. Let $\mathbb{V}\in {\mathbb{C}}^{N_{\mathrm{\Omega }}\times q}$ be the matrix of $q$ PCs, then the reduced transformed data are given by

$${\overline{\mathbb{Y}}}_d^{†}={\mathbb{Y}}_d^{†}\mathbb{V}.$$

(41)

Consider $N_{\mathrm{post}}$ samples are drawn from the posterior distribution $p\left({\overline{\mathbb{Y}}}_d^{†}|\overline{\mathrm{\Theta }},{\overline{\mathbb{X}}}^{†}\right)$ such that a particular sample is given by ${\overline{Y}}_d\in {\mathbb{C}}^q$ , then the true prediction for a particular dof is obtained using

$${\tilde{\mathbb{Y}}}_d={\overline{Y}}_d{\mathbb{V}}^{\mathrm{T}}.$$

(42)

Application of BPCR on a typical e-powertrain architecture

In this article, $V_{\mathrm{thresh}}=97.5\%$ and a Scree plot for a particular dof is shown in Figure 13, taking this as a criterion to select the number of PCs to retain.

Figure 13.

Scree plot depicting the total variance explained by each PC. The highlighted star-marked point indicates that $97.5\%$ of the total variation is explained by the first 6 PCs.

The magnitude of the posterior samples drawn from the posterior predictive distribution through the left mount $(M_1)$ to the vehicle-body at mechanical order $k=8$ , under unit displacement imposition, can be seen in Figure 14. This metamodel helps in understanding the dynamic behavior of the forces acting on the car-body. It can be observed that the model is able to predict the peaks due to resonance along with the uncertainty bands, as confirmed by the pattern observed in validation data. It is worth noting that this transformation introduces an inherent reconstruction error. However, considering that the selected PCs effectively account for 97.5% of the cumulative data variance, the inference drawn is that the reconstruction quality is deemed acceptable. Note that the results are demonstrated only for the transfer function through the left mount. Similar results can be plotted for the posterior samples through the right mount and the cradle point considered for the three-point mounting scheme. The results can then be combined to obtain the global acoustic level through structure-borne contribution.

Figure 14.

Magnitude of posterior predicted samples of the transfer function ${\mathit{H}}_{d,k}^{\mathrm{MB}}$ from the left mount to the vehicle body in all three directions for the mechanical order $k=8$ :(a) $x$ -direction, (b) $y$ -direction, and (c) $z$ -direction.

Interior sound pressure level through coupled stochastic metamodels

In this section, an example use-case is considered where low prior-knowledge is assumed on the physical parameters. Let $\tilde{\mathcal{X}},\tilde{\mathbb{X}}$ be the random input realizations from their respective input PDFs in e-motor and vehicle stage respectively, then the two stochastic metamodels developed in this work are

$$\tilde{\mathcal{M}}:\tilde{\mathcal{X}}\in {\mathbb{R}}^{N_p}\mapsto \tilde{\mathcal{Y}}\in {\mathbb{R}}^{N_f},$$

(43)

$$\widetilde{ℳ}:\tilde{\mathbb{X}}\in {ℝ}^{2N_d}\mapsto \tilde{\mathbb{Y}}\in {ℂ}^{N_{\mathrm{\Omega }}},$$

(44)

where $\tilde{\mathcal{Y}},\tilde{\mathbb{Y}}$ are the corresponding random responses. The model in equation (43) results in random EM-forces which in turn produces uncertain vibration displacements (refer to Section “Metamodeling for electromagnetic domain”) and the model in equation (44) captures the uncertainties in the transfer function from powertrain mounting system to the vehicle-body through BPCR framework including the inverse PCA map. Using equation (29), the posterior samples from each of the metamodels are combined to obtain the uncertain forces acting on the vehicle-body, which when multiplied with the measured transfer function data gives the SPL inside the cabin due to structure-borne contribution, as shown in Figure 15. The frequency range for the measurement was $[20,2048]$  Hz and the response was recorded for multiple locations inside the cabin (two front left seat, two back left seat, etc.). The mean of front left seat measurements is considered here.

Figure 15.

Interior SPL through left mount $\left(M_1\right)$ for a particular mechanical order $k=8$ : (a) interior SPL contribution in $x$ direction, (b) for $y$ direction, (c) for $z$ direction, and (d) shows the total contribution considering all three directions.

At each dof, the e-powertrain suspension metamodel has been trained and the acoustic levels are predicted at each specific mechanical order. The frequency resolution depends on the number of sampled speed profiles and therefore, $N_{\mathrm{\Omega }}=300$ is chosen which provides sufficient resolution and the BPCR framework can train models efficiently, which would have been time-consuming otherwise. The total contribution from a particular mount can be computed by taking the quadratic average of contributions in all three directions as a global NVH estimator.

Validation

The validation of the developed stochastic scheme for the structure-borne whining noise can be done by employing two distinct approaches: first, direct validation by measuring the SPL inside the cabin and then overlaying the measured data on the estimated dispersed data; second approach is, a block-wise validation scheme considering the electromagnetic, structural, and coupling domains separately, for the different domains shown in Figure 2. Despite serving as a crucial benchmark for method evaluation, with direct validation approach, the measured SPL data must be readily available for specific e-machine configurations, which is usually limited or nonexistent in the early stages of design. On the other hand, block-wise validation approach can be used in the absence of measured data and in identifying any limitations or discrepancies in each domain, allowing for targeted improvements or adjustments to the method developed.

Validation of EM-domain

The EM-domain is validated by comparing the average EM torque (as briefly described in Section “Validation of the simulation of EM-domain” and also refer to Pile et al. ²⁶ for more details on EM torque computation using Maxwell Stress Tensor method) from the simulated data using Pyleecan and the measured data used in the reference article by Yang et al., ²⁷ in which the measurement was done on the same e-machine configuration as used in this study (refer to the nominal parameters in Table 2). Figure 16 shows that the simulated EM torque values and thereby the AGSF agree well with the measured data at different OCs.

Figure 16.

Comparison of average electromagnetic torque from simulated data using the open-source tool and the reference data (measured in Yang et al. ²⁷ ) at $\mathrm{\Omega }=1000$  RPM for 13 OCs.

Validation of structural domain

In structural domain, the analytical estimation of natural frequencies $f_{\mathrm{mn}}^{\mathrm{stat}}$ in equation (10) and mode shapes of the simplified stator-system follow the classical thin-cylindrical shell theory equations as used in several references.^3,42 Figure 17 shows good correlation, especially for $n=1$ , of the estimated natural frequencies with the referenced data available from Gieras et al. ³ Some deviation from the reference data is expected as the e-machine is modeled using the open-source tool, where the micro-geometrical parameters are assumed to be fixed at a certain value. The estimation of RMS vibration displacement through modal expansion follows the classical approach of using the modal functions for a finite length cylindrical shell clamped at both ends. In the literature, such analytical estimations have been validated by using 2D/3D-finite element analyses. ⁵⁰

Figure 17.

Comparison of natural frequencies of the stator system (design parameters same as the one used in Gieras et al., ²² p. 117) between the referenced data and the estimated data for different $(m,n)$ pairs.

Validation of the in-house developed e-powertrain suspension model

For the validation of the in-house developed model for transmitting vibrations from e-motor to the mounts on the car-body side, the calculated vibration accelerations (using the internal Stellantis tool) at each dof is compared against the measured data for each spatial order. The measurement is done for a speed ramp-up $[600-13000]$  RPM, considering an in-house IPMSM (under mass-production) with 48 slots and eight poles but with different dimensions than the one used in Section “Metamodeling for electromagnetic domain.” The excitations acting on each mount are in the form of vibration displacements measured at the motor surface for different mechanical orders. The acceleration responses are simulated using the simplified linearized three-point suspension model depicted in Figure 12. In Figure 18, the vibration acceleration plots for an order of 8.86 are presented. It is evident from the analysis that the simulated (assumed deterministic) responses align closely with the measured data for each dof. Similar plots have been studied for different mechanical and spatial orders.

Figure 18.

Comparison of vibration accelerations on three mounts at a specific mechanical order 8.86: (a) shows the vibration acceleration on mount $M_1$ , (b) for mount $M_2$ , and(c) for mount $M_3$ .

Conclusions and outlook

In this article, stochastic predictive metamodels are developed within the Bayesian paradigm to efficiently incorporate a priori domain knowledge along with the available data. The developed framework considers Bayesian model training, validation, and exploitation through multivariate-regression along with the objective engineering based assessment of the hyper-parameters. Multiple metamodels allow the quantification and propagation of uncertainties from the EM-domain to the structural domain and then through the e-powertrain mounting system to the vehicle-body. The data used for modeling are generated using a semi-analytical approach. Due to the time-efficient sampling from the posterior distribution using MCMC process, further accelerated by BPCR framework, NVH risks can be assessed during the early stages of the design process. The knowledge gained could facilitate in prioritizing suitable noise mitigation techniques that can be implemented, as well as state of the art methodologies for addressing technical requirements.

For perspectives, uncertainties coming from the available measured data could be considered along with efficient data reduction techniques improving further the performance of the metamodel training stage. A comprehensive “global” metamodel would then include not just the estimates from structure-borne contribution but also from airborne contribution. Moreover, such metamodels would allow, as a straight forward usage of the stochastic metamodel, performing efficient sensitivity analysis to infer information about the most influential design parameters.