Introduction

In the automotive industry, rubber bushings are commonly used passive damping devices to suppress vibrations between components [29]. One such example is the assembly of a rear subframe to the body-in-white (BIW), which are connected through rubber bushings, e.g. see Bylin et al. [9] for such an application. The finite element (FE) method [6] is usually used to analyse such assembled structures, for which accurate models with proper model parameter settings are required. Furthermore, for connected structures, a good damping estimation can be important in the assembly’s predictive capabilities. The damping is usually not modelled based on first principles due to its complexity. Rather, it can be estimated from test data in calibration of FE models [15]. It is therefore of interest to update FE models towards experimentally obtained data from vibration tests [14], so that the FE model parameters can be identified, but also to obtain the damping characteristics of the structure. However, in the automotive industry, many nominally identical cars are built. Hence, using one parameter setting to analyse a whole population of cars can lead to faulty predictions if the spread between the nominally identical components is large. Then, it becomes important to understand which component properties are causing the spread and how large the variability is. Bushings being partly constructed from rubber materials will, to a large extent, share rubber’s properties, known to exhibit potentially large variation between samples [29]. It is therefore natural to assume that the largest variation in the rear subframe stem from the rubber in the bushings. Hence, a method to identify the rubber stiffness is of particular interest.

Various methods exist for estimation of dynamic stiffness of rubber bushings, only some are mentioned here. Hydraulic testing is many times used to estimate the dynamic stiffness in translational degrees-of-freedom (DOFs). However, stiffness around rotational DOFs are usually not measured, but estimated, e.g. see [28, 46]. In previous research, vibration testing has been used to estimate rubber’s dynamic stiffness, e.g. see Kari [31] and Ooi and Ripin [42]. Meggitt et al. [39] proposed an in situ method for estimating the dynamic stiffness of resilient elements using a substructuring approach. This approach was further investigated by Haeussler et al. in [26]. The substructuring approaches allow for estimating the dynamic stiffness of rotational DOFs, which can be especially important for higher frequency analysis. However, the use of rotational DOFs in coupling has been found to be important in experimental analytical dynamic synthesis also in the low-frequency range, e.g. see Bylin et al. [9].

In this paper, an FE model updating approach is taken to the bushing parameter estimation problem. This approach can also estimate the rotational stiffness components. Three nominally identical rear subframes of a Volvo XC90 (2015) are considered. The rear subframe is important from a vibration perspective of the assembled car as it connects the suspension link arms to the BIW. Therefore, there is interest in investigating the effect of the subframe variability on the system response in an assembly of the subframe and the BIW. Before component variability is propagated through the assembly a better understanding of the components can be necessary, validating the system assembly for a few parameter configurations. To this end, quantifying the overall variability between the three subframes, and also the variability between the bushing stiffness of the three components is a first step in such an analysis. In this approach, two variants of the FE model are updated. The subframe is first updated towards vibration data of each component, to minimise the subframe deviation to test data. Then, in a second step, the bushing stiffness parameters are updated towards test data, assuming that the fixture, in this case the main body of the subframe, is modelled correctly. To increase the bushing parameter identifiability in test data, the bushings are mass loaded, effectively bringing many modes with strong bushing deformation down in frequency where they can be clearly identified and easier analysed. Hence, two experimental configurations are considered for each component, and two variants of the FE model created, with and without mass loaded bushings.

Traditionally, bushings are modelled with generalised spring elements (known as a CBUSH element in Nastran). It is of interest to investigate how well the frequency response is predicted with a generalised spring element model compared to a more geometrically realistic solid model of the bushings. Such a study is performed here, where the model fit is compared to test data for the two types of bushing models. However, to simplify the analysis of the solid model, an isotropic linear elastic material model is assumed for the rubber in the bushing which Jones [29] states is sufficient for small strains, narrow frequency ranges and when no temperature variation is present. The model assumptions are further explained in “Solid Bushing Model”. Note that the generalised spring element allows for frequency dependency to be incorporated easily in the model, however this has not been investigated. Furthermore, damping is also easily introduced in the generalised spring element, but is not modelled in the FE models, mainly due to limitations in the calibration procedure explained below, but also due to its complexity for the remaining structure modelled with shell and solid elements. Instead, for both bushing models, mapped modal damping identified from experiments is used, assumed to be sufficient for this application.

FE model updating is a vast field, with many different approaches, e.g. see Friswell and Mottershead [15] for an overview of methods. Recently, stochastic methods have been proposed to the model updating problem, e.g. see Mares et al. [35], Mottershead et al. [41] and Govers and Link [23]. Beck and Katafygiotis [7] proposed a Bayesian approach to the model updating problem. In addition, perturbation methods have also been developed, e.g. see Mottershead et al. [40]. Non-probabilistic methods can be useful when little information is available, e.g. a fuzzy-interval method is presented by Haag et al. [25] and an info-gap approach by Ben-Haim [8]. A review of probabilistic and non-probabilistic uncertainty quantification methods for model updating can be found in Simoen et al. [43]. In this paper, an FE model updating procedure proposed by Abrahamsson and Kammer [3] is used. Successful results have been obtained in previous studies using this method on various structures, e.g. see [1, 2, 13, 16, 21, 33]. In the calibration procedure, a parametrised surrogate model of the FE model around the nominal parameter setting is first created. The calibration procedure then seeks to obtain a calibrated parameter setting minimising the deviation between the FE model’s frequency response functions (FRFs) and the experimentally identified model’s FRFs. Normalised damping is imposed on both models to circumvent the mode pairing problem. It was found in Abrahamsson and Kammer [3] that a relatively high damping yields a smooth deviation metric. However, a consequence of the introduction of equalised damping is that no damping parameters can be updated. Using FRFs with equalised damping has several advantages over modal based approaches. Modally dense models with mode crossings can be handled without any need for mode matching. Also, applying mode matching metrics to the possibly noise and error corrupted test modes is avoided. Consequently, model updating in the mid-frequency range is possible, see [3] for more details. To increase the confidence in the estimated parameters an uncertainty quantification (UQ) procedure has been proposed for the above FE model updating method by Vakilzadeh et al. [44]. That method, however, is not always feasible for large test data sets. Therefore, another method is proposed here, based on creating a linear-in-parameters surrogate model, at the calibrated parameter setting, which is used together with bootstrap data sampling to quantify the influence of measurement noise on the identified parameters. The FRF based approach also has advantages for the UQ step over modal approaches. The parameter uncertainties can be quantified with respect to the noise in the raw FRFs. Therefore, no intermediate step is necessary to quantify the modal parameter uncertainties first.

A brief review of the model updating theory is presented in “Finite Element Model Updating Theory”, with a section describing the linear-in-parameters UQ method. In “Finite Element Models”, the FE models are presented along with the parameter selection. Section “Experimental Modal Analysis” describes the vibration testing procedure and system identification. The results are presented in “Results” and in Conclusion” the paper is concluded.

Finite Element Model Updating Theory

The deterministic model updating procedure used here has been developed by Abrahamsson and Kammer [3], and will be briefly restated. A stochastic variant was proposed by Vakilzadeh et al. [44]. However, it has been observed that for very large data sets, not uncommon for test data consisting of FRFs, the stochastic updating procedure is not always feasible. Therefore, a new method is proposed here, using a linear-in-parameters surrogate model, described in “Model Parameter Uncertainty Quantification”.

The model updating problem can be formulated as a minimisation problem in which a parameter setting \(\hat {\boldsymbol {p}}\) minimising the deviation between the experimental FRFs \(\boldsymbol {H}^{{\mathrm {X}}} (\omega _{{{k}}}) \in \mathbb {C}^{{n_{\mathrm {y}}} \times {n_{\mathrm {u}}}}\) and FE FRFs \(\boldsymbol {H}^{{\text {FE}}} (\omega _{{{k}}}) \in \mathbb {C}^{{n_{\mathrm {y}}} \times {n_{\mathrm {u}}}}\) is sought, for \({{k}} = 1, \dots , {{K}}\) discrete frequencies \(\omega _{{{k}}} \in [\underline {\omega }, \overline {\omega }]\). Here \(\underline {\omega }\) and \(\overline {\omega }\) denote the lower and upper frequency of interest, respectively, and ny and nu denote the number of outputs and inputs, respectively. Below, the notation \(\boldsymbol {H}_{{{k}}} \triangleq \boldsymbol {H}(\omega _{{{k}}})\) will be used for brevity. For the UQ problem, an estimate of the probability distribution, or commonly only the first and second moments, of the model parameters is sought, with respect to experimental data. Following the expansion in, e.g. [32, 44], the experimental FRF matrix can be expressed for each frequency as

$$ \boldsymbol{H}_{{{k}}}^{{\mathrm{X}}} = \boldsymbol{H}_{{{k}}}^{{\mathrm{R}}} + \boldsymbol{N}_{{{k}}}^{{\mathrm{O}}} = \boldsymbol{H}_{{{k}}}^{{\text{FE}}} + \boldsymbol{N}_{{{k}}}^{{\mathrm{M}}} + \boldsymbol{N}_{{{k}}}^{{\mathrm{O}}} \triangleq \boldsymbol{H}_{{{k}}}^{{\text{FE}}} + \boldsymbol{N}_{{{k}}}^{{\mathrm{G}}} $$
(1)

where \(\boldsymbol {H}_{{{k}}}^{{\mathrm {R}}}\) denote the true frequency response of the real structure and \(\boldsymbol {N}_{{{k}}}^{{\mathrm {O}}}\) denote the measurement noise. The true system can be expressed in terms of a possibly inaccurate FE model \(\boldsymbol {H}_{{{k}}}^{{\text {FE}}}\) and a model prediction error \(\boldsymbol {N}_{{{k}}}^{{\mathrm {M}}}\). The model error and noise can be grouped into a term called the observed prediction error \(\boldsymbol {N}_{{{k}}}^{{\mathrm {G}}}\). It is assumed that \(N_{ij}^{{\mathrm {G}}} (\omega _{{{k}}})\), for \(i = 1, \dots , {n_{\mathrm {y}}}\) and \(j = 1, \dots , {n_{\mathrm {u}}}\), can be modelled as independent, zero mean, circularly complex normally distributed random variables with variance \(\sigma _{ij}^{2}\). Validating that assumption through very many tests is beyond the scope of this paper.

The calibration procedure can be divided in roughly four equally important steps. Collecting vibration test data, system identification, deterministic calibration to find the calibrated parameter setting p and a UQ step to asses the parameter uncertainty with respect to measurement noise and errors. The vibration test data acquisition and system identification steps are described in detail in “Experimental Modal Analysis”.

Deterministic Model Updating Procedure

The third step in the calibration procedure consists of a deterministic calibration procedure to find p from the nominal FE model parameter setting p0. Here, the identified model’s FRFs with equalised damping \(\tilde {\boldsymbol {H}}_{{{k}}}^{{\text {ID}}}\) are used as test data. With equalised damping between the FE model and test data, it is possible to circumvent the issue of mode pairing during the calibration process. Using the identified model, the frequency response can be evaluated at any frequency, and so the calibration can be made more computationally efficient. In [3] the authors propose choosing a few frequencies per half-band-width of the eigenmodes. The calibration procedure is defined as a constrained, to the parameter set \(\mathbb {P}\), minimisation of a log least squares function

$$ \begin{array}{@{}rcl@{}} \boldsymbol{p}^{*} & =& \underset{\boldsymbol{p} \in \mathbb{P}}{\arg\min} \frac{1}{{N}} \sum\limits_{{{k}} = 1}^{{{K}}} \boldsymbol{\epsilon}_{{{k}}}^{\mathrm{H}} (\boldsymbol{p}) \boldsymbol{\epsilon}_{{{k}}} (\boldsymbol{p}) \end{array} $$
(2a)
$$ \begin{array}{@{}rcl@{}} \boldsymbol{\epsilon}_{{{k}}}(\boldsymbol{p}) & = & \log_{10} \left( \text{vec}\left( \boldsymbol{H}_{{{k}}}^{{\text{FE}}}(\boldsymbol{p}) \right) \right) - \log_{10} \left( \text{vec}\left( \tilde{\boldsymbol{H}}_{{{k}}}^{{\text{ID}}} \right) \right) \end{array} $$
(2b)

The deviation metric 𝜖k(p) is set-up to be smooth and not discriminate against deviations at frequencies where the structural response is small. The motivation behind the particular deviation metric is further explained in Abrahamsson and Kammer [3]. Here HFE(p) represents the FE model, at parameter setting p, the superscript \(\cdot ^{^{\mathrm {H}}}\) denote the conjugate transpose and vec(⋅) stands for the vectorisation operation in which a matrix is transformed into a vector by stacking its columns. The deviation metric is scaled with the number of total data points N = nynuK. To minimise the non-linear least squares function in Eq. 2, the Levenberg-Marquardt optimisation algorithm [4, 34, 36] is used. Many start parameter settings are used, generated by the Latin hypercube sampling technique [37], from which p is selected from the best calibration outcome. After the calibration, a mode pairing algorithm [47] is used to map the experimental damping to the correct FE model modes.

The calibration procedure utilises the mass and stiffness matrices formed from an FE representation of the structure of interest. The equations of motion of the considered mechanically vibrating systems can be written as [11]

$$ \boldsymbol{M}\ddot{\boldsymbol{q}}(t) + \boldsymbol{V}\dot{\boldsymbol{q}}(t) + \boldsymbol{K}\boldsymbol{q}(t) = \boldsymbol{f}(t) $$
(3)

where the dot notation is used for time differentiation and \(\boldsymbol {M} \in \mathbb {R}^{m \times m}\) represent the positive definite mass matrix and \(\boldsymbol {V} \in \mathbb {R}^{m \times m}\) and \(\boldsymbol {K} \in \mathbb {R}^{{{m}} \times {{m}}}\) represent the positive semi-definite viscous damping and stiffness matrices, respectively. The general displacement vector is denoted by q(t) and the external force vector by f(t). Explicit time dependence is dropped from here on. A first-order equivalent system for acceleration outputs y of Eq. 3 can be obtained by forming the state vector \(\boldsymbol {x} = [\boldsymbol {q}^{\mathrm {T}}, \dot {\boldsymbol {q}}^{\mathrm {T}}]^{\mathrm {T}}\) such that

$$ \begin{array}{@{}rcl@{}} \dot{\boldsymbol{x}} & =& \boldsymbol{A}\boldsymbol{x} + \boldsymbol{B}\boldsymbol{u} \end{array} $$
(4a)
$$ \begin{array}{@{}rcl@{}} \boldsymbol{y} & =& \boldsymbol{C}\boldsymbol{x} + \boldsymbol{D}\boldsymbol{u} \end{array} $$
(4b)

with \(\boldsymbol {A} \in \mathbb {R}^{{n} \times {n}}\) the system matrix, \(\boldsymbol {B} \in \mathbb {R}^{{n} \times {n_{\mathrm {u}}}}\) the input matrix, \(\boldsymbol {C} \in \mathbb {R}^{{n_{\mathrm {y}}} \times {n}}\) the output matrix and \(\boldsymbol {D} \in \mathbb {R}^{{n_{\mathrm {y}}} \times {n_{\mathrm {u}}}}\) the direct throughput matrix. Here, y represents the system acceleration outputs and the load vector f in Eq. 3 is related to the system input vector u with the Boolean transformation matrix U as f = Uu. The system matrices are then

$$ \begin{array}{@{}rcl@{}} \boldsymbol{A} & =& \left[\begin{array}{cc} \boldsymbol{0} & \boldsymbol{I} \\ -\boldsymbol{M}^{-1}\boldsymbol{K} & -\boldsymbol{M}^{-1}\boldsymbol{V} \end{array}\right] \end{array} $$
(5a)
$$ \begin{array}{@{}rcl@{}} \boldsymbol{B} & =& \left[\begin{array}{cc} \boldsymbol{0} \\ \boldsymbol{M}^{-1} \boldsymbol{U} \end{array}\right] \end{array} $$
(5b)

where C and D are formed appropriately so that linear combinations of the system states x and inputs u form the system acceleration outputs y. Here, I is the identity matrix of appropriate dimension. The quadruple

$$ \boldsymbol{\Sigma}^{{\text{FE}}} = \left\{ \boldsymbol{A}, \boldsymbol{B}, \boldsymbol{C}, \boldsymbol{D} \right\} $$
(6)

denotes the state-space system built from the FE model. The transfer function matrix \(\boldsymbol {H}_{{{k}}}^{{\text {FE}}}\) at frequency ωk of the FE model can be computed as

$$ \boldsymbol{H}_{{{k}}}^{{\text{FE}}} = \boldsymbol{C} (\mathrm{i} \omega_{{{k}}} \boldsymbol{I} - \boldsymbol{A})^{-1} \boldsymbol{B} + \boldsymbol{D} $$
(7)

with the imaginary number i2 = − 1.

To avoid mode pairing and to regularise the deviation metric Abrahamsson and Kammer [3] proposed to impose a damping equalisation on the experimentally identified system ΣID, and the FE system ΣFE. To this end, the experimentally identified system ΣID is diagonalised by a similarity transformation [22] using the eigenvector matrix X of the eigenvalue problem

$$ \boldsymbol{A}^{{\text{ID}}} \boldsymbol{X} = \boldsymbol{X} \boldsymbol{\Lambda}^{{\text{ID}}} = \boldsymbol{X} \text{diag}\left( \lambda_{1}^{{\text{ID}}},\dots,\lambda_{r}^{{\text{ID}}}\dots,\lambda_{{n}}^{{\text{ID}}} \right) $$
(8)

where ΛID is the diagonal eigenvalue matrix, with \(\lambda _{r}^{{\text {ID}}}\) the r:th complex eigenvalue. For small damping values the relative damping \(\xi _{r}^{{\text {ID}}}\) is

$$ \xi_{r}^{{\text{ID}}} = - \frac{ \operatorname{Re}\left( \lambda_{r}^{{\text{ID}}} \right) }{ \left\vert \operatorname{Im}\left( \lambda_{r}^{{\text{ID}}} \right) \right\vert } $$
(9)

with \(\left \vert \cdot \right \vert \) the absolute value and \(\operatorname {Re}\left (\cdot \right )\) and \(\operatorname {Im}\left (\cdot \right )\) the real and imaginary part of a complex number, respectively. In the damping equalisation step, the modal damping is set to a fixed value for all modes, i.e \(\xi _{r}^{{\text {ID}}} = \xi _{0} \forall r\). A new state-space system with equalised damping is denoted with the quadruple set

$$ \tilde{\boldsymbol{\Sigma}}^{{\text{ID}}} = \left\{ \tilde{\boldsymbol{\Lambda}}^{{\text{ID}}}, \boldsymbol{X}^{-1}\boldsymbol{B}^{{\text{ID}}}, \boldsymbol{C}^{{\text{ID}}}\boldsymbol{X}, \boldsymbol{D}^{{\text{ID}}} \right\} $$
(10)

with the associated FRF matrix \(\tilde {\boldsymbol {H}}_{{{k}}}^{{\text {ID}}}\). The eigenvalues \(\tilde {\lambda }_{r}^{{\text {ID}}}\) in the state matrix \(\tilde {\boldsymbol {\Lambda }}^{\text {ID}}\) for \(r = 1,\dots ,{{n}}\) are

$$ \tilde{\lambda}_{r}^{{\text{ID}}} = \left\{\begin{array}{ll} \operatorname{Im}\left( \lambda_{r}^{{\text{ID}}}\right)\left( -\xi_{0} + i\right) &\forall \operatorname{Im}\left( \lambda_{r}^{{\text{ID}}}\right) > 0 \\ \operatorname{Im}\left( \lambda_{r}^{{\text{ID}}}\right)\left( \xi_{0} + i\right) &\forall \operatorname{Im}\left( \lambda_{r}^{{\text{ID}}}\right) < 0 \end{array}\right. $$
(11)

For the FE model with given mass M and stiffness K matrices it is possible to impose the same level of relative damping by forming the modal viscous damping as described in [11]

$$ \boldsymbol{V} = 2\xi_{0} \boldsymbol{M} \boldsymbol{T} \text{diag}\left( {\Omega}_{1}, \dots, {\Omega}_{{{m}}} \right) \boldsymbol{T}^{\mathrm{T}} \boldsymbol{M} $$
(12)

where \(\text {diag}\left ({\Omega }_{1},\dots ,{\Omega }_{{{m}}}\right )\) is the diagonal matrix of eigenfrequencies and T the mass normalised eigenvector matrix from the undamped eigenvalue problem of the system in Eq. 3

$$ \boldsymbol{K}\boldsymbol{T} = \boldsymbol{M}\boldsymbol{T}\text{diag}\left( {{\Omega}^{2}_{1}},\dots,{\Omega}^{2}_{{{m}}}\right) $$
(13)

Usually full FE models are too expensive to solve at each parameter setting in an optimisation procedure. Therefore, parametrically reduced-order models must be constructed. The parametrised mass and stiffness matrices can be expressed as \(\boldsymbol {M} = \boldsymbol {M}\left (\boldsymbol {p}\right )\) and K = K(p). In the implementation, the parameters are normalised. The physical parameter vector \(\boldsymbol {p} \in \mathbb {R}^{{{P}} \times 1}\) is related to a normalised parameter vector \(\bar {\boldsymbol {p}}\) and some fixed non-zero nominal parameter p0 so that \(\boldsymbol {p} = \boldsymbol {p}_{0} \circ \left (1 + \bar {\boldsymbol {p}} \right )\) holds [3], with ∘ the Hadamard product. Here P denotes the number of parameters. The eigenvector matrix T0 at the nominal parameter setting, computed as in Eq. 13, is used as a reduction basis and is kept constant during the calibration procedure. Using the coordinate transformation q = T0z and pre-multiplying Eq. 3 by \(\boldsymbol {T}_{0}^{\mathrm {T}}\) the reduced mass and stiffness matrices at any parameter setting p are given by

$$ \begin{array}{@{}rcl@{}} \bar{\boldsymbol{M}}(\boldsymbol{p}) & =& \boldsymbol{T}_{0}^{\mathrm{T}} \boldsymbol{M}(\boldsymbol{p}) \boldsymbol{T}_{0} \end{array} $$
(14a)
$$ \begin{array}{@{}rcl@{}} \bar{\boldsymbol{K}}(\boldsymbol{p}) & =& \boldsymbol{T}_{0}^{\mathrm{T}} \boldsymbol{K}(\boldsymbol{p}) \boldsymbol{T}_{0} \end{array} $$
(14b)

At the nominal parameter setting p0, the reduced mass and stiffness matrices are represented by M0 and K0. Gradients of the reduced mass and stiffness matrices are computed as

$$ \begin{array}{@{}rcl@{}} \bar{\boldsymbol{M}}_{,s} & =& \boldsymbol{T}_{0}^{\mathrm{T}} \left( \frac{\partial \boldsymbol{M}(\boldsymbol{p})}{\partial p_{s}} |_{\boldsymbol{p}=\boldsymbol{p}_{0}} \right) \boldsymbol{T}_{0} \end{array} $$
(15a)
$$ \begin{array}{@{}rcl@{}} \bar{\boldsymbol{K}}_{,s} & =& \boldsymbol{T}_{0}^{\mathrm{T}} \left( \frac{{\partial \boldsymbol{K}(\boldsymbol{p})}}{{\partial p_{s}}} |_{\boldsymbol{p}=\boldsymbol{p}_{0}} \right) \boldsymbol{T}_{0} \end{array} $$
(15b)

for the s:th calibration parameter. From the first-order expansions of the Taylor series of the reduced mass and stiffness matrices about po it is possible to form a surrogate model that is linear in the parameters

$$ \begin{array}{@{}rcl@{}} \tilde{\boldsymbol{M}}(\boldsymbol{p}) & =& \bar{\boldsymbol{M}}_{0} + \sum\limits^{{{P}}}_{s=1} \left( p_{s} - p_{s,0}\right) \bar{\boldsymbol{M}}_{,s} \end{array} $$
(16a)
$$ \begin{array}{@{}rcl@{}} \tilde{\boldsymbol{K}}(\boldsymbol{p}) & =& \bar{\boldsymbol{K}}_{0} + \sum\limits^{{{P}}}_{s=1} \left( p_{s} - p_{s,0}\right) \bar{\boldsymbol{K}}_{,s} \end{array} $$
(16b)

Here, ps,0 is the s:th parameter at the nominal parameter setting. The parametrised state-space matrices with equalised damping, forming the system \(\tilde {\boldsymbol {\Sigma }}^{{\text {FE}}}(\boldsymbol {p})\), can be computed analogously to Eq. 5b with the parametrised equalised viscous damping matrix computed as in Eq. 12. Then, the parametrised FE FRFs \(\tilde {\boldsymbol {H}}_{{{k}}}^{{\text {FE}}}(\boldsymbol {p})\) at some frequency ωk, used in the calibration procedure, can be computed as in Eq. 7 from the state-space matrix quadruple \(\tilde {\boldsymbol {\Sigma }}^{{\text {FE}}}(\boldsymbol {p})\).

Model Parameter Uncertainty Quantification

From the deterministic calibration procedure, a best parameter configuration p is obtained, with respect to the identified and equalised damping model \(\tilde {\boldsymbol {\Sigma }}^{{\text {ID}}}\). From a decision-making perspective it is important to asses the uncertainties in the obtained parameters so that confidence in the predictions can be obtained. Therefore, the last step considers the parameter uncertainties, in which a new surrogate model is computed around the calibrated parameter configuration p. In [44] an undamped Gauss-Newton optimisation procedure [4] was proposed together with bootstrapping [27] to estimate the parameter uncertainties. However, for large data sets that procedure can be computationally infeasible. Therefore, a less computationally intensive method is proposed here, based on a first-order Taylor expansion. As such, it is assumed to work well for small parameter variations.

Consider the following deviation metric computed at the parameter configuration p, for the K frequencies, as

$$ \boldsymbol{\gamma}(\boldsymbol{p}) = \left[\begin{array}{cc} \text{vec}\left( \boldsymbol{H}_{1}^{{\text{FE}}}(\boldsymbol{p}) - \boldsymbol{H}_{1}^{{\mathrm{X}}} \right) \\ \vdots\\ \text{vec}\left( \boldsymbol{H}_{{{K}}}^{{\text{FE}}}(\boldsymbol{p}) - \boldsymbol{H}_{{{K}}}^{{\mathrm{X}}} \right) \end{array}\right] $$
(17)

for which a linear-in-parameters model takes the form

$$ \boldsymbol{\gamma}(\boldsymbol{p}) = \boldsymbol{E}(\boldsymbol{p} - \boldsymbol{p}^{*}) + \boldsymbol{\gamma}(\boldsymbol{p}^{*}) $$
(18)

where the data matrix is

$$ \boldsymbol{E} = \left[ \frac{{\partial \boldsymbol{\gamma}(\boldsymbol{p})}}{{\partial p_{1}}} |_{\boldsymbol{p}=\boldsymbol{p}^{*}}, \dots, \frac{{\partial \boldsymbol{\gamma}(\boldsymbol{p})}}{{\partial p_{s}}} |_{\boldsymbol{p}=\boldsymbol{p}^{*}}, \dots, \frac{{\partial \boldsymbol{\gamma}(\boldsymbol{p})}}{{\partial p_{{{P}}}}} |_{\boldsymbol{p}=\boldsymbol{p}^{*}} \right] $$
(19)

Bootstrapping is used to estimate the parameter uncertainties, which is cheap given the linear-in-parameters model.

A bootstrap parameter setting for \(b = 1,\dots ,{n_{\mathrm {b}}}\) is found as

$$ \boldsymbol{p}_{b} = \boldsymbol{E}_{b}^{+} \boldsymbol{\gamma}_{b} $$
(20)

where ⋅+ denote the Moore-Penrose inverse [22]. The procedure works by drawing nb random datasets γb and Eb with replacement from the original data sets γ and E, respectively. See Hastie et al. [27] for a detailed explanation of the re-sampling procedure. Thus, nb vectors of parameter values pb are obtained, from which the expected value can be computed as

$$ \hat{\boldsymbol{p}} = \operatorname{E} \left[\boldsymbol{p}_{b}\right] = \frac{1}{{n_{\mathrm{b}}}} \sum\limits_{b = 1}^{{n_{\mathrm{b}}}} \boldsymbol{p}_{b} $$
(21)

and the variance as

$$ \operatorname{Var} \left[\boldsymbol{p}_{b}\right] = \frac{1}{{n_{\mathrm{b}}} - 1} \sum\limits_{b = 1}^{{n_{\mathrm{b}}}} \left( \boldsymbol{p}_{b} - \operatorname{E} \left[\boldsymbol{p}_{b}\right] \right)^{2} $$
(22)

Finite Element Models

The rear subframe FE model with generalised spring element and solid bushing models is described here. Also, the FE model with mass loaded bushings is described. Both FE model configurations have free-free boundary conditions, to represent the experimental data as much as possible. In both configurations, dissipation in the updated models is modelled with modal damping identified from experimental data, assumed sufficient for this application.

Rear Subframe

The rear subframe FE model, seen in Fig. 1, is built up of mainly 3-noded and 4-noded Mindlin-Reissner plate elements and solid elements. The mesh size is approximately 5 mm, resulting in a total of 420324 DOFs for the subframe model with generalised spring elements. In the same figure, the sensor locations are also shown with circles, and excitation locations with a grey square, further described in “Rear Subframe”. In Fig. 2, the subframe parametrisation is indicated. The first seven parameters are thickness parameters corresponding to subcomponents of the subframe. The last three parameters are Young’s modulus parameters of the marked subcomponents. Note that, thickness parameters 1 to 4 are mainly used because it was assumed that the bushing models did not correctly represent the bushing mass. Therefore, these parameters are used as surrogate parameters to control the mass distribution at the extremities of the model. Parameter 7 (thickness) and 8 (Young’s modulus) are related to the same subcomponent, which actively contributed to the dynamics of many modes. The nominal parameter settings can be found in Table 1.

Fig. 1
figure 1

FE model of the rear subframe (a) top view and (b) bottom view. Circular markings indicate positions of accelerometers measuring normal to the surface and the rectangular marking indicate the force input location (normal to the surface) with a direct accelerometer configuration

Fig. 2
figure 2

Representative subframe model. To the left (a) thickness parameters and to the right (b) stiffness parameters

Table 1 Parameter values of the rear subframe FE model with generalised spring element bushings, see Fig. 2. COV [%] in parenthesis

Mass Loaded Bushings

The subframe model with generalised spring elements and mass loaded bushings consist of 475554 DOFs. The only addition in this model compared to the subframe model is the extra quasi-rigid bodies attached to the four bushings. In Fig. 3, the FE model with generalised spring elements and mass loaded bushings is shown. The sensor locations are also shown. The four excitation positions are marked with grey squares (located at the very corners of the quasi-rigid bodies). Three orthogonal inputs per excitation location were used. For each output, marked with white circles, a triaxial accelerometer was used. This is further described in “Mass Loaded Bushings”. Note that the two front bushings are nominally identical, as are the two rear bushings. Therefore, the mass loading is different for two bushings of the same type, to as much as possible avoid mode multiplicity.

Fig. 3
figure 3

FE model of the rear subframe with mass loaded bushings (a) front view and (b) rear view. Circular markings indicate positions of accelerometers measuring normal to the corner surfaces and rectangular markings indicate the force input locations (normal to the corner surfaces)

In Fig. 4a, the generalised spring elements used to model the rubber bushings are numbered. Note that, one generalised spring element is indicated with each marking in that figure. Every generalised spring element has six dynamic stiffness parameters, three translational and three rotational. Hence, 24 dynamic stiffness parameters can be updated in this model.

Fig. 4
figure 4

Mass loaded rear subframe models. To the left (a) the model with generalised spring element bushings and to the right (b) the model with solid bushing models. Notice that the front and rear solid bushings are not geometrically identical in (b), also shown in Fig. 5

Solid Bushing Model

A subframe with bushings modelled with solid elements from the bushing geometry is considered as well. This subframe FE model has 622578 DOFs, and the mass loaded variant has 681474 DOFs. This is considerably more DOFs compared to the generalised bushing model. In Fig. 4b, a cut-through view of the solid bushing models is shown. A more detailed cross-section view of the solid bushings can be seen in Fig. 5. The area in dark grey indicates the rubber material which is parametrised. The light grey areas are plastic parts, one of which is enclosed in the rubber, and the other outside. In mid grey (in the middle of the bushings) is a hollow aluminium cylinder. The plastic and aluminium parts are not parametrised. Rubber can be characterised as a viscoelastic material, see Jones [29]. It is known that rubber’s Young’s modulus is frequency and temperature-dependent, again see Jones [29], which is also shown in Cunningham and Ivey [12] for some specific rubber materials. Here, small strains are assumed and non-linear effects are considered negligible. Temperature effects are also disregarded. It is furthermore assumed that the frequency range of interest is small enough for the Young’s modulus to be considered frequency independent. Therefore, the rubber is modelled with solid elements and an isotropic linear material model, which according to Jones [29] is sufficient under the above-mentioned assumptions. Note that, only the Young’s modulus of the rubber is parametrised in the four bushings, dark grey area in Fig. 5. Furthermore, while the rubber is modelled as an isotropic material, the geometry of the whole bushing will allow the bushing to have different stiffness in different directions. Hence, four stiffness parameters can be varied in this model. This is significantly fewer parameters compared to the model with generalised spring element bushings, which allowed for 24 stiffness parameters.

Fig. 5
figure 5

Cross-section of the front (a) and rear (b) FE models of the bushings. Plastic parts in light grey, aluminium part in mid grey and rubber part in dark grey

Experimental Modal Analysis

Vibration tests have been performed on three rear subframes of the Volvo XC90 (2015) with and without mass loaded bushings. Identical procedures were used for all three components. Calibration of sensors, data acquisition and FRF estimation were performed using the open-source toolbox abraDAQ [17]. The rear subframe without mass loading is shown in Fig. 6a and with mass loaded bushings in Fig. 6b. The components were hung in long thin high-strength lines attached to helical springs on a supporting steel structure. Five support modes were all below 2 Hz, while a support bouncing mode was just under 5 Hz, i.e. the quasi rigid body modes were well under the first flexible mode of the component in both configurations (approximately 75 and 40 Hz without and with mass loaded bushings, respectively).

Fig. 6
figure 6

To the left (a) experimental set up of the rear subframe and to the right (b) experimental set up of subframe with mass loaded bushings

Rear Subframe

The rear subframe FE model was used in pre-test planning to find good sensor locations. A total of 170 candidate sensor locations were selected, from which 20 uniaxial sensor locations were found for the 20 first flexible modes using the expansion effective independence method [30] with an added gramian rejection step for rejection of locations with similar modal information [19]. In addition, six more accelerometer locations were added for visualisation purposes. Also, an additional 10 triaxial sensors were placed on the 6 bushings. However, these sensors were not used in system identification or calibration. The sensor locations can be seen in Fig. 1. The input location, see Fig. 1, was selected based on shaker attachment convenience. An accelerometer was placed on the opposite side (from the shaker) of the sheet metal to obtain a direct accelerance.

Two types of accelerometers were used, 10 triaxial PCB Piezotronic type 356A03 weighing 1 gramme each and 26 uniaxial PCB Piezotronic type PCB 352C22/NC weighing 0.5 grammes each. The uniaxial accelerometers were attached with synthetic wax while the triaxial accelerometers were glued. The accelerometer masses were included in the FE models. The shaker used in the vibration test was of Ling Dynamic Systems make and of type V201, with a metallic stinger approximately 5 mm in length, as recommended by Ewins [14]. The excitation force was measured with a Brüel & Kjær force sensor of type 8203 with an IEPE converter 2647B, attached to the component through a stinger attachment plate which was glued and had a mass of around 0.2 grammes. The force sensor and stinger attachment plate masses were not modelled in the FE model.

Two excitation methods were used, and in both tests, a sampling rate of 20 kHz was used. Periodic chirp tests at various amplitudes were performed to assess the linearity of the system to obtain proper excitation loads. A frequency range from 20 to 800 Hz was used. The calibration data was collected with a stepped sine excitation signal with 10 simultaneous sinusoidal components. For such an excitation signal, test data is collected reasonably fast and the measurement noise is kept minimal to reduce the influence on the estimated parameters in the model updating procedure. After the system reached stationary, the frequency responses were estimated from data collected in 25 consecutive instances. A frequency range from 60 to 500 Hz was selected with 3000 frequency lines. The frequency lines were selected based on the improved frequency selection method [45] for system identification with N4SID [38].

An identified system is necessary for the calibration procedure. A recently proposed state-space system identification method using residual states is used for this purpose [18, 20]. For single-input multi-output systems, as considered here, only one re-estimation of the system input and output matrices is necessary, according to [24]. The 26 uniaxial accelerometers were used for system identification, calibration and validation. System identification was performed on data from 60 to 400 Hz with 15 in-band modes, and 2 higher frequency residual modes. This produced a system that accurately predicted both resonance and antiresonance behaviour. In Fig. 1 calibration channels are represented by white circles, and validation channels with dark circles. Note that the direct accelerance channel 14 was used for validation. Approximately two-thirds of the channels were used for calibration.

Mass Loaded Bushings

Experiments on the mass loaded subframe were conducted as well. Here 12 triaxial accelerometers were used, placed according to Fig. 3. FE modes were analysed and accelerometers placed so that all modes could be observed. Input locations are also shown in Fig. 3, where three orthogonal inputs per location were used. Hence, a total of 12 inputs on the four different bushings were obtained. As obvious from the results, this was enough to control local modes in each bushing below 300 Hz.

The accelerometers used were of Brüel & Kjær make and type 4524-B, weighing 4.4 grammes, and attached with synthetic wax. A modal hammer of PCB make and type 086C01 was used for excitation with a hard plastic tip. Soft silicon buds were placed on the structure to mechanically filter out high-frequency force components, and also allow for better hit repeatability. A sample rate of 10 kHz was used, with a fade time of 10 s and a cut-off frequency at 500 Hz, which resulted in approximately 0.1 Hz frequency stepping. However, to reduce the computational time for system identification, the frequency content was filtered so that approximately 0.5 Hz frequency stepping was obtained. The FRFs were estimated from the average response from 3 repeated hits. The mass of the accelerometers was not included in the FE model.

The same system identification procedure was used here as for the rear subframe [18, 20]. The data is multi-input multi-output and therefore a few (typically four) iterations of the system input and output matrices were necessary to reach convergence of the least squares loss function. One-third of the data was used for validation, all outputs for inputs in the y-direction (see Fig. 3), resulting in one input per bushing. Roughly two-thirds of the channels were used for system identification and calibration, inputs in the x and z directions, with all outputs except for outputs in the y-direction for triaxial accelerometers at locations 5, 9, 11 and 12, see Fig. 3. System identification was performed on data from 40 to 340 Hz with 28 in-band modes, and 16 high-frequency residual modes. This produced an identified system that accurately predicted both resonance and antiresonance behaviour.

Results

Here, the calibration results are presented for the rear subframe without and with mass loaded boundaries for the case with bushings modelled with generalised spring elements. A calibration with solid bushing elements for the mass loaded case is also presented. The experimental model is denoted EX while the finite element model is denoted FE, and the three components with numbers 1, 2 and 3. The nominal model is indicated with a superscript N, the three updated FE models with a superscript C, and the uncertainty quantification models with superscript UQ. Note that often both a subscript and superscript are used indicating a comparison between an FE model and one experimental model.

Multiple Rear Subframes

The subframe model was updated from frequency data ranging from 60 to 265 Hz. In that range, 8 flexible modes were present for all three subframes. A 0.5% equalised damping was used with 4 frequency lines per half-band-width. The order of the reduced linear-in-parameters model was based on 50 modes. 1000 Latin hypercube samples were used to obtain good starting points for calibration minimisations. From these samples, 11 random starts, including the nominal parameter setting, were used in the optimisation.

The nominal and updated rear subframe parameter values are shown in Table 1. See Fig. 2 for parametrisation. The parameters were allowed to vary up to 40% from the nominal setting, apart from parameters 2 to 4 which were allowed to take any value. Note that parameter 1 is fixed, mainly due to high dependence on parameter 2. It was also manually updated from 1.4 mm to 2.5 mm. This was done to compensate for the generalised spring element bushing model, which had less mass than the solid bushing model. However, parameters 2 to 4 were updated, mainly because it was found that the mass distribution had a large influence on modes 3 and 4, and 5 and 6. Note that the Young’s modulus values are mainly used as surrogate parameters as the more likely cause of stiffness variation is due to uncertain geometry variation that is more natural but harder to parametrise.

In Fig. 7, the relative parameter change in relation to the nominal parameters is shown for the calibration and UQ parameter configurations. It can be seen that both calibration and UQ outcomes show similar change, and also that the Young’s modulus parameters show small relative change. The thickness parameters on the other hand show a large change compared to the nominal setting. They also indicate a larger variation. It is also noted in Table 1 that the coefficient of variation (COV) is larger for the thickness parameters.

Fig. 7
figure 7

Relative difference for calibrated and mean of bootstrapping parameter settings (relative the nominal parameter setting) for the rear subframe with generalised spring element bushings

In Table 2, the masses of the real subframes as measured by a scale and the masses of the nominal and updated FE models are shown. It can be seen that the mass was too low in the nominal model, at most 2.5% difference. However, in the calibrated models it is too high, but less than 1.4% difference for all three components.

Table 2 Rear subframe mass [kg] (FE model with generalised spring element bushings)

In Table 3, the eigenfrequencies for the first eight flexible modes are shown. Figure 8 shows the absolute differences in eigenfrequency between the nominal FE model and experimentally identified models, and the updated FE models and experimentally identified models. Note that modes 2, 6, 7 and 8 show the most variation. For example, mode 7 of the nominal FE model shows an error exceeding 10 Hz, for all three subframes. The variation in the updated models is significantly lowered, indicating that the calibrated models match their experimental counterpart better.

Table 3 Rear subframe eigenfrequencies [Hz] (FE model with generalised spring element bushings)
Fig. 8
figure 8

Absolute difference in eigenfrequency for the nominal and calibrated rear subframe with generalised spring element bushings (relative the experimentally identified eigenfrequencies)

Figure 9 shows the modal assurance criterion (MAC) matrix [5]. The top row shows the MAC matrix between the nominal and experimental models, while the bottom rows shows the updated counterpart. Note that all three updated models manage to correct the two mode shifts present in the nominal model.

Fig. 9
figure 9

MAC for first eight modes of the nominal (top) and calibrated (bottom) rear subframe with generalised spring element bushings. Subframe 1 to the left, 2 middle and 3 right. FE model along abscissa and experimental model along ordinate axes

In Table 4, the deviation metric from Eq. 2b is shown for the nominal and calibrated models with the same level of equalised damping. It can be observed that a significant reduction is achieved for both calibration and validation data, although slightly less for the validation data. Note that the deviation is also larger for the nominal model for validation data.

Table 4 Deviation metric (defined in Eq. 2b) for the rear subframe model with generalised spring element bushings

Figure 10 shows one example of an FRF, in this case for output 10 (see Fig. 1). For this FRF the nominal model’s deviation metric (defined in Eq. 2b) is 1.23 and the calibrated model’s deviation metric is 0.29, with mapped damping.

Fig. 10
figure 10

FRF for output 10 (see Fig. 1) of rear subframe 2 with generalised spring element bushings

Multiple Rear Subframes with Mass Loaded Bushings

The subframe model with mass loaded bushings was updated using frequency data from 40 to 265 Hz. In that range, 26 flexible modes were present for all three subframes and a reduced linear-in-parameters model of order 70 was used. The rotational stiffness parameters of the bushings were first found in an calibration with 3% equalised damping, 6 frequency lines per half-band-width and 100 random starts. In the final calibration, all parameters were free, apart from parameters 13, 14 and 19 that were fixed due to low identifiability. For the final calibration, a 2% equalised damping was used with 6 frequency lines per half-band-width. The nominal parameter setting together with 5 selected samples, out of 2000 Latin hypercube samples, were used in the optimisation procedure.

The nominal and updated rear subframe parameter values are shown in Table 5. In that table, the first group of 6 parameters represent bushing 1 as seen in Fig. 4a. The next group of 6 parameters represent bushing 2, the third group bushing 3 and fourth group bushing 4. The six parameters for each bushing are organised with stiffness in x, y and z first and then with rotational stiffness around x, y and z. The parameters related to the bushings were allowed to vary up to 150% from the nominal setting to reflect the uncertainty level of the bushing properties. The nominal values of the bushing rubber material are estimated from a separate test at 100 Hz and large strain. Therefore, for the relatively small strains in this calibration test, and large frequency range covered, it is reasonable to suspect clear variation from the nominal values.

Table 5 Stiffness [N/mm] parameter values of the generalised spring element models for the mass loaded subframe model, see Fig. 4. COV [%] in parenthesis

Figure 11 shows the relative parameter change in relation to the nominal parameters for the calibration and UQ parameters configurations. Note that results are mostly consistent for each parameter for the three subframe individuals. The rotational stiffness around the z-axis (parameters 6, 12, 18 and 24) show a decrease in stiffness up to 45% while all other parameters show an increase in stiffness. Note that the dynamic stiffness in x and y directions (parameters 1, 2, 7, 8, 13, 14, 19 and 20) show the largest variation between bushings and subframes. For those parameters the reported COV was also largest, see Table 5.

Fig. 11
figure 11

Relative difference for calibrated and mean of bootstrapping parameter settings (relative the nominal parameter setting) for the mass loaded rear subframe with generalised spring element bushings

Table 6 shows the eigenfrequencies for the first 26 flexible modes. Figure 12 shows the absolute difference in eigenfrequencies between the nominal and updated FE models and the experimentally identified models. The differences are significantly lowered, indicating that the calibrated models match their experimental counterpart better.

Table 6 Mass loaded subframe eigenfrequencies [Hz] with generalised spring element bushings
Fig. 12
figure 12

Absolute difference in eigenfrequency for the nominal and calibrated mass loaded rear subframe with generalised spring element bushings (relative the experimentally identified eigenfrequencies)

In Fig. 13, a set of MAC matrices are shown. The top row shows the MAC matrices between the nominal and experimental models, while the bottom rows shows the updated counterpart. The nominal model show very poor MAC correlation towards all three subframes. The calibrated models essentially capture all 26 modes well. The mode switch for models FE1 and FE3 are mainly due to the very close modes in test data, i.e. for FE1 it is at 128.83 Hz and at 128.96 Hz, and for FE3 it is at 171.50 Hz and at 171.91 Hz. It was found that changing eigenvalue solver could cause the mode switch.

Fig. 13
figure 13

MAC for the nominal (top) and calibrated (bottom) mass loaded rear subframe with generalised spring element bushings. Subframe 1 to the left, 2 middle and 3 right. FE model along abscissa and experimental model along ordinate axes

Table 7 shows the deviation metric from Eq. 2b for the nominal and calibrated models with the same level of equalised damping. It can be observed that a significant reduction is achieved for both calibration and validation data.

Table 7 Deviation metric (defined as in Eq. 2b) for the mass loaded rear subframe model with generalised spring element bushings

Figure 14 shows an example of an FRF, in this case for input at location 1 in the z direction and output at location 5 in z direction (see Fig. 3). The nominal deviation metric (defined in Eq. 2b) for this channel is 1.72, and for the calibrated model with mapped damping it is 0.27.

Fig. 14
figure 14

FRF of mass loaded rear subframe 2 with generalised spring element bushings for input at location 1 in the z-direction and output at location 5 in z-direction (see Fig. 3)

Comparison with Solid Bushing Model

The subframe model with solid bushing elements was calibrated using the same calibration settings as the subframe with generalised spring elements. However, only one calibration step was sufficient as only four parameters were calibrated. Here, a Poisson’s ration of 0.49 was used, which is judged to be a reasonably realistic value for the rubber material in the bushing, see Cardarelli [10]. The calibrated result from only one model is shown, as the solid model cannot capture the local bushing dynamics as well as the generalised spring element model. In Table 8, the parameter settings are shown. Note that the two rear bushings, corresponding to parameters 1 and 2, obtain calibrated values that are relatively close, which is reasonable. Also, the two front bushings, corresponding to parameters 3 and 4, show similar values. The mean value of the UQ is here close to the calibrated counterpart, and the COV is low. The nominal and calibrated Young’s modulus values are reasonable for rubber materials, see Cardarelli [10]. In Table 9, the eigenfrequencies are shown. Figure 15 shows the relative difference in eigenfrequencies for the nominal and calibrated models (relative the experimentally identified eigenfrequencies). Note that modes 9, 10, 11, 15, 17 and 18 of the updated model have a larger error compared to the nominal model. The MAC values are shown in Fig. 16. It can be seen that a slight improvement is achieved for the first few modes only. Table 10 shows the deviation metric from Eq. 2b for the nominal and calibrated models with the same level of equalised damping.

Table 8 Young’s modulus [MPa] parameter values of the solid bushing model for the mass loaded rear subframe model, see Fig. 4. COV [%] in parenthesis
Table 9 Mass loaded subframe eigenfrequencies [Hz] with solid bushing model
Fig. 15
figure 15

Absolute difference in eigenfrequency for the nominal and calibrated mass loaded rear subframe with solid bushing model (relative the experimentally identified eigenfrequencies)

Fig. 16
figure 16

MAC for nominal (left) and calibrated (right) mass loaded subframe with solid bushing model. FE model along abscissa and experimental model along ordinate axes

Table 10 Deviation metric (defined as in Eq. 2b) for the mass loaded rear subframe model with solid bushing model

Conclusion

A model updating procedure using damping equalisation has been successfully used to update the overall model of a rear subframe and find optimal bushing parameters, for three nominally identical subframes. An uncertainty quantification procedure using a linear-in-parameters surrogate model was introduced to assess the parameter uncertainties with respect to measurement noise and errors. Two bushing models have been compared, the generalised spring element and a more geometrically realistic solid model. The spring element model has six individual dynamic stiffness and therefore larger freedom to adapt to test data. That larger freedom is found necessary to obtain a good fit. The solid bushing model, for which the rubber was modelled with a linear elastic material model, shows poorer correlation with experimental data, and could not capture the many local bushing modes. It is possible that the material model is too simple or that the rubber material properties are not uniform. The overall subframe dynamics showed smaller variation in eigenfrequencies, compared to the bushings, which is expected. The largest variation was found in the translational stiffness. However, these parameters were also most uncertain, indicating that test data is not informative enough for these parameters to be updated accurately. The updated parameters for the three components were consistent and physically realistic. In conclusion, the simple mass loading considered here seem most effective for estimating rotational dynamic stiffness parameters, which are usually not measured by separate procedures. They are, therefore, often considered to be the most uncertain parameters.