<figure><img src="image_1.png" /></figure>

<figure><img src="image_2.png" /></figure>

# New Method for Accurate and Efficient Transient Simulation of Finite Element Structures with Interference Fits

W. Witteveen<sup>*1*</sup>, D. Müller<sup>2</sup>, S. Dietz<sup>2</sup>, L. Koller<sup>3</sup>, and F. Pichler<sup>3</sup>

## Abstract

An interference fit is a common joining technique used to connect a shaft and a hub. In the presence of dynamic loads and vibrations, characteristic variables such as contact pressure and slippage are load and state-dependent quantities. Such effects have either not been investigated using previous simulation methods or have only been addressed in a simplified manner. The reason for this is the nonlinear contact between the shaft and the hub, which makes a Finite Element simulation with fine meshing, while taking all dynamic effects into account, very demanding and significantly increases the computational effort. This work also offers a new and alternative view of contact modes. This perspective is particularly advantageous for structures with initial stresses that occur in the presence of an interference fit. In this paper, so-called contact modes are applied to interference fits with some modifications. This closes the previously mentioned gap in the simulation landscape because it allows nonlinear, accurate, and fast numerical time integration of finely meshed Finite Element models with interference fits, without the need for simplifications regarding contact, friction, and dynamics due to vibrations and nonlinear rigid body motion. Local plasticity and temperature fields were not taken into consideration.

**Keywords:** Interference Fit, Dynamics, Contact Modes, Numerical Time Integration, Tribomechodynamics

Received on May 7, 2025, Accepted on January 5, 2026, Published on January 7, 2026

## 1 Introduction

Interference fits (also known as press fits or shrink fits) are widely used in engineering applications, especially in the assembly of mechanical components where tight tolerances and secure connections are crucial. Interference fits involve the insertion of one component into another such that the inner component has a slightly larger diameter than the outer one, resulting in a compressive force that holds the components together without the need for fasteners. The numerical time integration of Finite Element (FE) structures with interference fits is highly demanding due to the nonlinear contact problem between the two components involved. If plastic deformations or temperature fields are taken into account, the computation becomes even more complex. In this work, we restrict ourselves to an accurate characterization of the contact problem, neglecting plastic material behavior and temperature fields. This is a condition that exists in many press fits. Of course, this kind of problem also exists in similar form in other connections (e.g., bolted connections), so some of the findings of this work can also be applied to them, but we limit ourselves here to interference fits alone.

Due to the computational effort involved, publications can be found that deal with quasi-static loading in which the shrinkage process and additional loads are considered while the inertia forces are neglected. In [1], a ring is

<footnote n="1"><sup>1</sup>Upper Austria University of Applied Sciences, Department of Mechanical Engineering, Roseggerstraße 15, 4600 Wels, Austria</footnote>

<footnote n="2"><sup>2</sup>Dassault Systems Deutschland GmbH, Friedrichshafener Straße 1, 82205 Gilching, Germany</footnote>

<footnote n="3"><sup>3</sup>FH OOE Research and Development Ltd., Roseggerstraße 15, 4600 Wels, Austria</footnote>

<footnote>*wolfgang.witteveen@fh-wels.at</footnote>
shrunk onto a shaft, but no additional load case is considered. It turns out that the Finite Element Method (FEM) provides much more realistic results than analytical expressions, such as Lamés equation. If the geometry in or near the contact area becomes even more complex (e.g., holes or grooves), it can be assumed that only FEM can provide useful results. In [2], the shaft is subjected to bending after shrinking has been applied. The relevant conclusions are that the FE mesh needs to be fine, and that the most common contact formulations (penalty contact, augmented Lagrange, and pure Lagrange) can provide accurate results. In [3], the contact pressure at the edge of an oil groove is investigated because damage occurs at this location due to high stress concentrations. One of the findings from this study is that the FE mesh has to be extremely fine to accurately capture such phenomena. The final static investigations to be discussed here are taken from [4]. In this publication, a shrinking process with various interferences is examined. Under smaller interferences, the stresses in the contact area remain in the linear elastic range, whereas local plastification occurs with larger interferences. After shrinking, the hub is subjected to tensile loading until global sliding occurs and the interference fit fails. In this work, experiments were compared with FEM and analytical formulas. The conclusions relevant to this work are: (1) on the test bench and in the FEM simulation, failure of the interference fit occurs under similar tensile forces. This means that, in principle, FEM is suitable for the simulation of such processes, and (2) analytical formulas provide poor estimates and do not represent a reliable alternative to FEM.

The following conclusions can be drawn from the cited publications relating to static loading: (1) FEM can be used to compute interference fits. (2) FEM provides much better results than analytical formulas. (3) Relatively fine FE meshes are required. If local stress concentrations are to be resolved, very fine meshes are necessary.

There is no reason to assume that these conclusions should be different in dynamics. The application of FEM to transient processes requires numerical time integration, taking into account the nonlinear contact and friction forces, as well as a fine FE mesh. The immense computational effort behind such simulations is likely the reason why the authors of this paper have not been able to find a publication reporting on such simulations. However, the availability of such simulation technology would be very relevant, as the contact pressures, the friction stresses, and thus the reliability of an interference fit are certainly influenced by vibrations and varying loads. This paper proposes a strategy that closes this gap in the simulation landscape and makes this type of computation possible. The following citations are examples of studies that deal with dynamics but avoid time integration of FE models, likely for the same reasons as mentioned previously. The influence of an interference fit in a transfer path is investigated numerically and experimentally in [5]. In their experiment, the transmission behavior was measured under different interferences, and an influence was observed. No FE model was used for the simulation; rather, this quasi-static simulation is based on a numerical solution of analytical expressions. A hysteresis is calculated by unloading and loading as a function of the interference. The experimental investigations reveal an influence of the interference on the dynamics and thus emphasize the importance of a realistic consideration of an interference joint in a transient simulation. In [6], the interference fit between the inner ring of a rolling bearing and a shaft is investigated. The change in interference due to centrifugal forces and temperature fields was taken into account. A nonlinear contact force was modeled, and inertia forces were taken into account; however, the simulations were not carried out using a 3D FE model, rather self-derived analytical equations were solved numerically. The conclusions regarding the global vibration behavior were verified by experiment. Therefore, it is not expected that this approach can provide an accurate contact- and friction-force distribution in the contact zone. In [7], a similar problem was solved using comparable approaches (without FEM). The conclusions are also similar, as is the fact that the stress distributions on the contact surfaces are not sufficiently accurate.

The few publications on the dynamics of interference fits can be summarized as follows: (1) Interference fits influence the dynamics of a system. (2) The solution strategies are very problem-specific. (3) Numerical time integration using 3D FEM models is not commonly performed likely due to the computational effort required, although this is not explicitly stated.

This literature review revealed that, in cases where static and dynamic investigations of interference fits are concerned, a simulation strategy with the following characteristics would be reasonable: (1) the use of FE models with fine meshes and a high degree of detail; (2) correct nonlinear consideration of contact and friction forces; (3) consideration of inertia forces due to vibrations and nonlinear rigid body motion; and (4) consideration of local plastification and temperature fields. This paper presents an accurate and efficient simulation strategy that fulfils points (1) to (3) and concludes with brief comments on point (4).

The paper is organized as follows: In the theory chapter, after a few introductory definitions, it is shown that contact modes can be computed not only with modal derivatives (as in previous studies) but also with a Taylor series expansion of the nonlinear contact force. This approach is particularly useful for prestressed situations, such as interference fits. This section concludes with the presentation of a special modeling technique that makes it possible to model the contact in the interference fit as a self-contact. This is of significant importance for the efficiency of the
method. The interference fit of a shaft with a pulley was chosen as a numerical example. The contact area was finely meshed, and static load cases were used to show that the solution converges to the FEM results with an increasing number of contact modes. A dynamic comparison with FEM is not possible. However, it can be observed that the results converge with an increasing number of contact modes. Finally, it is demonstrated that a modal basis, computed for a particular interference, can also be used for other interferences.

## 2 Theory

The use of so-called contact modes is proposed in this work so that realistic deformations and stresses can be computed in the contact area of an interference fit. This section consists of four subsections. In the beginning, some general relations and definitions are introduced. The second subsection recapitulates the concept of contact modes as has been published in the literature. The next subsection presents a new, equivalent but alternative view of contact modes. The advantage of this view lies particularly in its application to structures with initial stress, as is the case with interference fits. The last section explains how two FE bodies (shaft and hub) can be treated as one body. This is crucial for the efficiency of the method.

The following two subsections provide a very compact summary of contact modes based on modal derivatives. More information, derivations, and examples can be found in publications of the author and his co-workers, see [8], [9] and [10]. To avoid misunderstandings, it is briefly explained what contact modes are used for in this work — and what they are not used for. The objective is not the separate computation of the properties of two substructures, A and B, to subsequently couple them somehow via compatibility conditions. Contact modes pursue a different goal: the efficient and accurate consideration of nonlinear contact at the interface between A and B. The modes of A and B are not computed separately; instead, A and B always form a combined structure, even during mode computation. The connection between A and B is achieved through nonlinear contact and friction forces.

### 2.1 Introduction

Efficient dynamics of FE structures requires model order reduction because nodal degrees of freedom would lead to far too many equations. Model order reduction via projection approximates a ($N \times 1$) vector $x$ of nodal displacements (or rotations) in the form of

$$x = \Phi q \qquad (1)$$

where the ($N \times Q$) matrix $\Phi$ holds $Q$ trial vectors in its columns. The $i$-th column of such a matrix is denoted hereafter on as $\Phi_i$. These trial vectors are scaled by $Q$ scaling factors that are collected in the ($Q \times 1$) vector $q$. The $i$-th entry of this vector will be denoted as $q_i$. The trial vectors are commonly called 'modes' and the scaling factors 'modal coordinates'. Assuming a given modal basis $\Phi$, the undamped differential equation of a FE model with an interference fit

$$M\ddot{x} + Kx = f + f_{nl} \qquad (2)$$

can be transformed into a reduced system

$$M_r\ddot{q} + K_rq = f_r + \Phi^T f_{nl} \qquad (3)$$

by inserting (1) into (2). Here, $M$ and $K$ are the ($N \times N$) mass and stiffness matrices of the linear FE model, $f$ is the ($N \times 1$) vector of external forces, and $f_{nl}$ is the ($N \times 1$) vector resulting from contact and friction forces in the interference fit. The double dots over $x$ indicate second time derivatives. In the reduced system, $M_r(= \Phi^T M \Phi)$ and $K_r(= \Phi^T K \Phi)$ represent the ($Q \times Q$) mass and stiffness matrices of the reduced system, and $f_r(= \Phi^T f)$ is the ($Q \times 1$) vector of reduced external forces. The goal is to find a reduction matrix $\Phi$ such that, when the solution of (3) is inserted into (1), the difference from the result of the full (non-reduced) system (2) is small. If additionally $Q \ll N$, then the reduced equation can be solved much faster than the full system, while still yielding sufficiently accurate results. More information on model order reduction via projection can be found in sections 3 and 4 of the doctoral thesis of Rutzmoser [11].

Two groups of modes must be considered in order to obtain useful solutions for contact problems. Therefore, the matrix $\Phi$ is subdivided into

$$\Phi = [ \begin{matrix} \Phi_G & \Phi_C \end{matrix} ] \qquad (4)$$
where the ($N \times G$) submatrix $\Phi_G$ holds common modes, such as the well-known ones of Craig and Bampton [12]. A reduction using only $\Phi_G$ leads to internal degrees of freedom that correlate with vibration modes, and to interface degrees of freedom on which external forces can act. These interface degrees of freedom are usually FE nodal degrees of freedom that are preserved, for example, where the structure connects to other structures. Of course, it is theoretically conceivable that all degrees of freedom in the contact region of the interference fit are treated as the aforementioned interface degrees of freedom. For a finely meshed FE structure, this is not a practically feasible approach since the number of reduced degrees of freedom $Q$ remains large, so that the system cannot be integrated efficiently in time. In the literature, several approaches have been proposed to reduce the number of interface degrees of freedom in such cases. One example is the work by Craig and Chang [13]. In their study, three solution strategies are proposed, all of which are based on a static condensation of the involved interface degrees of freedom. First, an eigenvalue problem is solved using the mass and stiffness matrices reduced to the interface. The deformation field obtained in this way can then be used as trial vectors to approximate the interface displacements. A similar idea is used in [14]. Second, Craig and Chang propose a Guyan reduction, in which the displacement field is reconstructed based on a few remaining displacement degrees of freedom at the interface. Third, they suggest that arbitrary Ritz trial vectors may also be employed. In a certain sense, the present work can be regarded as a realization of this latter approach. However, the Ritz trial vectors here are not computed on the basis of a static condensation of all interface degrees of freedom. Consequently, the matrix $\Phi_G$ is computed as if the contact region does not exist. The extension of the modal basis with the contact modes in $\Phi_C$ has the aim that the deformations within the contact region can be represented accurately enough, even though the number of modes in this submatrix is much smaller than the number of nodal degrees of freedom in the contact region. The computation of these contact modes is addressed later. For the moment, they can be taken as given and the following assumptions are now made for all further considerations: (1) Both contact surfaces belong to the same body. They are therefore in self-contact, and the deformations of both surfaces are described by corresponding lines of $\Phi$.

(2) There is small sliding contact between the surfaces.

(3) For simplicity, it is assumed that both contact surfaces are congruently meshed with $P$ node-to-node contact pairs.

(4) Also, for the sake of simplicity, it is assumed that the normal and tangential directions of the contact surface are represented as degrees of freedom in the displacement vector. This means that a planar contact surface is normal to one of the three axes of the reference coordinate system of the FE model.

Thus, the ($P \times 1$) vectors $x_n^1$ and $x_n^2$ can be constructed from $x$, containing all the normal displacements of the nodes on contact surface 1 and contact surface 2, respectively. In our case, contact surface 1 would belong to the hub and contact surface 2 to the shaft— or vice versa. In a reduced model, these vectors are also approximated using (1) so that $x_n^1 = \Phi_n^1 q$ and $x_n^2 = \Phi_n^2 q$. The ($P \times Q$) matrices $\Phi_n^1$ and $\Phi_n^2$ contain the corresponding rows from $\Phi$ that are associated with the degrees of freedom in $x_n^1$ and $x_n^2$. Based on these, the relative normal displacements of the two contact surfaces can be computed by a single matrix vector multiplication in the form

$$g = x_n^2 - x_n^1 = \Phi_n q \qquad (5)$$

where the ($P \times 1$) vector $g$ holds the normal gap at all $P$ node-to-node contact pairs. The $i$-th entry of this vector is denoted hereafter on as $g_i$. More details on the computation of the ($P \times Q$) matrix $\Phi_n = (\Phi_n^2 - \Phi_n^1)$ can be found in Pichler et al. [10], especially in cases where the previously mentioned simplifications do not apply.

In a completely analogous manner, the two ($P \times 1$) tangential displacement vectors $t_1$ and $t_2$ can also be computed according to

$$t_1 = \Phi_{t_1} q \qquad t_2 = \Phi_{t_2} q \qquad (6)$$

In order to better illustrate some of the considerations, an exemplary penalty contact model is useful. This is given by

$$p_i = p_i(g_i) = \begin{cases} 0 & g_i > 0 \\ k_p g_i^\alpha & g_i \le 0 \end{cases} \qquad (7)$$

The value $p_i$ stands for the contact pressure at the $i$-th node pair. The quantities $k_p$ and $\alpha$ are parameters of the contact model. All contact pressures are collected in the ($P \times 1$) vector $p$.

As mentioned, the matrix $\Phi$ (including the contact modes) is used for model reduction. Numerical time integration is carried out with the reduced model, in which the contact and friction forces are considered nonlinearly, similar to the FEM. As this is described in great detail in the literature, see [10] and [9], it will not be discussed here, with the exception of a few comments in Section 3.2.2. The focus in the following theory section is on the computation of contact modes as presented in the literature and on an alternative approach for doing this.
<figure><img src="image_3.png" /><figcaption>Fig. 1: Computational flow</figcaption></figure>

## 2.2 Some Comments on the Computational Flow

Fig. 1 shows the procedure for computing a mode basis including contact modes.

(1) The process starts with a given matrix $\Phi_G$, as described in the previous section.

(2) Based on modal derivatives or a Taylor series expansion, $D(P \times 1)$ pressure load vectors $p_1$ to $p_D$ are computed. Their exact computation will be explained later.

(3) Since the space spanned by $p_1$ to $p_D$ is highly redundant, a principal component analysis (PCA) is applied in the next step. The first C result vectors of the PCA (those with the highest eigenvalues) are denoted as $P_1$ to $P_C$.

(4) The pressure distributions $P_1$ to $P_C$ are applied as load cases in a linear FE analysis. The result is a set of C ($P \times 1$) deformation vectors $\Phi_{C,1}$ to $\Phi_{C,C}$.

(5) The final mode basis is constructed by extending $\Phi_G$ with the vectors $\Phi_{C,1}$ to $\Phi_{C,C}$.

To better understand the overall workflow, some additional details are provided next:

* It should be noted that there are also alternative PCA approaches. If the displacements according to the former pressures are used and weighted with the stiffness matrix, the strain energy is approximated in the best possible way. This was examined in detail in [9]. It turns out that using strain energy actually requires fewer contact modes for comparable result quality. However, the deformations must be computed in that case. This requires an additional FE computation with many load cases. To avoid this extra effort, PCA is applied to the loads, and it is accepted that more contact modes are required.

* Before PCA is applied, the vectors $p_1$ to $p_D$ are scaled so that their Euclidean norm is one.

* For this work, the space spanned by $\Phi$ is of central importance. For the actual time integration, $\Phi$ is further orthogonalized with respect to mass and stiffness, so that the reduced mass and stiffness matrices become diagonal. This is a well-documented process and not of particular interest for this work.

* As shown in equation (4), the final mode basis $\Phi$ is formed by combining $\Phi_G$ and $\Phi_C$. In this context, the modes in $\Phi_C$ can also be interpreted as correction vectors. These correction vectors reduce the displacement error in the contact zone caused by the omission of high-frequency modes and are referred to as residual vectors in many FE programs, including the software used in this work.

## 2.3 Contact Modes Based on Modal Derivatives: A Brief Review

As can be seen in detail in [9], the pressures of step (2) in Fig. 1 can be computed along

$$p_{ij} = \left[ \frac{\partial K_n}{\partial q_i} \right] \Phi_{n,j} = K_{i,n} \Phi_{n,j} \qquad (8)$$

where the ($P \times 1$) vector $p_{ij}$ contains the contact pressure of the modal derivative according to the modes with the numbers i and j. The ($P \times 1$) vector $\Phi_{n,j}$ contains the j-th column of the matrix $\Phi_n$ as it was defined in Eq. 5. The ($P \times P$) matrix $K_{i,n}$ contains the change in stiffness within the contact area due to mode i. In the undeformed reference state, $K_{i,n}$ is the zero matrix, as both sides of the contact area do not interact with each other. The change in the stiffness matrix with respect to the undeformed reference configuration is the contact stiffness itself. If (7) is derived with respect to the gap, the contact stiffness of node pair k ($1 \le k \le P$) is $2k_p g_k$ for $\alpha = 2$. A small change to $q_i$ (all other modes remain undeformed) leads to penetration in some areas of the contact region and to gaping in others. If there is penetration in the pair of nodes k, this leads to an entry $2k_p \Phi_{n,i,k}$ at the k-th position of the main diagonal of
$K_{i,n}$. The scalar $\Phi_{n,i,k}$ holds the k-th entry of $\Phi_{n,i}$. The multiplication of this matrix with the j-th mode $\Phi_{n,j}$ finally leads to the entry

$$p_{i,j,k} = \begin{cases} 0 & \Phi_{n,i,k} > 0 \\ 2k_p \Phi_{n,i,k} \Phi_{n,j,k} & \Phi_{n,i,k} \le 0 \end{cases} \qquad (9)$$

at position k of $p_{ij}$. Of course, the perturbation of the reference state is also conceivable with a negative $q_i$. In this case, all the areas of the contact region that were just closed would open up and vice versa. Due to this asymmetry, the modal derivatives in the negative direction must also be taken into account. Thus, in addition to $p_{ij}$, $p_{-ij}$ needs to be considered as well with the entry

$$p_{-ij,k} = \begin{cases} 2k_p \Phi_{n,i,k} \Phi_{n,j,k} & \Phi_{n,i,k} > 0 \\ 0 & \Phi_{n,i,k} \le 0 \end{cases} \qquad (10)$$

at the k-th position.

The pressures $p_{11}$, $p_{12}$ to $p_{GG}$ and $p_{-11}$ to $p_{-GG}$ obtained in this way correspond to the pressures $p_1$ to $p_D$ of step (2) in Fig. 1.

Numerous computations have already been successfully carried out by the author and his co-workers with these contact modes. One example is [15], where it was shown that this approach is capable of catching measured nonlinear effects, such as amplitude-dependent frequency and damping of modes. In none of those examples were modal derivatives computed around a preloaded state (e.g., bolt preload), although improved convergence could be expected in such cases. In Section 3.1 of [16], it is explained in more detail why this approach is not straightforward. To still obtain good results, the mode basis was extended in almost all examples by an additional mode representing the deformation at the operating point. This trial vector is a result of a nonlinear FE computation. The Taylor-series-based computation of contact modes, which is presented in the following section, does not have this drawback, as it can be easily applied to both unloaded and pre-loaded structures.

## 2.4 An Alternative View on Contact Modes: Taylor Series of Contact Pressure

Now it is assumed that the contact pressure $p$ is expanded around an operating point $x_0$ with the help of a Taylor series. This can be the undeformed reference configuration or a preloaded configuration (e.g., due to bolt loading or shrinking). The corresponding variables at the operating point are labeled $p_0$, $q_0$ and $g_0$. Formally, this can be given as

$$p = p_0 + \sum_{i=1}^{G} \left. \frac{\partial p}{\partial q_i} \right|_{p_0} \Delta q_i + \frac{1}{2} \sum_{j=1}^{G} \sum_{i=1}^{G} \left. \frac{\partial p}{\partial q_i \partial q_j} \right|_{p_0} \Delta q_i \Delta q_j + \dots \qquad (11)$$

The series expansion is terminated after the quadratic terms. For the linear terms, the contact law defined in (7) results in the entry

$$\left( \left. \frac{\partial p}{\partial q_i} \right|_{p_0} \right)_k = p_{i,k} = \begin{cases} 0 & g_{0,k} > 0 \\ k_p \alpha g_{0,k}^{\alpha-1} \Phi_{n,i,k} & g_{0,k} \le 0 \end{cases} \qquad (12)$$

at the k-th position of the vector $p_i$. The quantity $g_{0,k}$ holds the penetration of node pair k at the operating point. For the quadratic terms

$$\left( \left. \frac{\partial p}{\partial q_j \partial q_i} \right|_{p_0} \right)_k = p_{ji,k} = \begin{cases} 0 & g_{0,k} > 0 \\ k_p \alpha (\alpha - 1) g_{0,k}^{\alpha-2} \Phi_{n,i,k} \Phi_{n,j,k} & g_{0,k} \le 0 \end{cases} \qquad (13)$$

can be given as entry at the k-th position of the vector $p_{ji}$. Note that all vectors $p_i$ and $p_{ji}$ corresponds to the pressures $p_1$ to $p_D$ of step (2) in Fig. 1.

### 2.4.1 Undeformed Reference Configuration

This configuration is characterized by $p_0 = 0$, $q_0 = 0$ and $g_0 = 0$. The application of $g_0 = 0$ to equation (12) reveals that the first-order derivatives vanish for $\alpha > 2$. It is interesting to observe, that the second-order derivatives in (13) take on a form very similar to the modal derivatives (9) when $\alpha = 2$ is used. The only difference is the limit of when the
<figure><img src="image_4.png" /><figcaption>Fig. 2: Bending of a shaft-hub connection</figcaption></figure>

pressure is zero and when it is not. Since the gap is zero in the undeformed state anyway, these limits make no sense. The expression $\partial p/\partial q_i$ should then be better interpreted as the change in contact pressure due to perturbation by mode $i$. With this reasoning, the detection of the contact region changes from $g_{0,k} \le 0$ into $\Phi_{n,i,k} \le 0$. This establishes formal equality between (9) and (13). The argument for the necessity of the negative perturbation is similar to that for the modal derivatives. This results in exactly the same pressure vectors as for the modal derivatives. Thus, it could be shown that the previously used modal derivatives can also be interpreted as a Taylor series expansion of the contact pressure around the undeformed reference configuration. Although the two approaches are formally equivalent at this point, the advantage of the Taylor series approach is clear in the case of contact modes around preloaded situations.

## 2.4.2 Deformed Configuration

Examples of deformed configurations are press fits or preloaded bolt connections. In this case, there are already deformations and contact pressures in the contact area without a dynamic load acting ($p_0 \ne 0, q_0 \ne 0$ and $g_0 \ne 0$). In this case, first-order and second-order derivatives do not vanish. The relations (12) and (13) only need to be evaluated with a certain contact law. Note that now there is a clearly defined contact area. Therefore, the negative derivatives ($p_{-i}$ and $p_{-ji}$) can be omitted.

It should be noted that $x_0$ must of course be part of the mode basis. This requires a nonlinear static FE computation. It is recommended that the deformation $x_0$ is already taken into account as a mode in $\Phi_G$.

Now, there are two ways of how the final contact modes can be computed:

Variant 1: The first- and second-order derivatives $p_i$ and $p_{ji}$ are plugged together into a single PCA. The first $C$ result vectors obtained in this way (the ones with the highest eigenvalues) serve as loads for the subsequent linear FE computation.

Variant 2: The first-order derivatives are fed into a separate PCA, and the second-order derivatives are also fed into a separate PCA. The $C$ load vectors for the final FE computation are determined by first using the result vectors of the first-order PCA and then those of the second-order PCA.

Which of the two variants converges better is investigated numerically in Section 3. It turns out that variant 2 converges better.

With respect to literature, it can be stated that the idea of using the derivative of the contact pressure with respect to certain displacements as additional basis vectors is not new. In [17] and [18], mode bases are computed using a comparable idea. Both publications focus on the multi-harmonic balance method (MHBM) around a preselected frequency of interest. Therefore, there are no modal coordinates as responses, but rather amplitudes of different harmonics with increasing frequency. Consequently, the change in the contact force is computed with respect to the amplitude of the harmonics involved, and not, as in the present case, with respect to modal coordinates. It can be assumed that this mode basis performs better in the case of a response analysis using the MHBM, since it is more specifically tailored to it. However, if a mode basis is sought for a time integration with arbitrary and varying excitations, then the approach cannot be directly used, and the proposed method seems more promising.

## 2.5 Relative Rigid Body Modes Between Shaft and Hub

In this work, the Floating Frame of Reference Formulation (FFRF) is used to represent the flexible body in the multibody simulation, see Shabana [19] for more information on that well-established approach. The main idea of the FFRF is to split the overall motion of a body into a large nonlinear rigid body movement and small linear deformations. These small deformations are computed based on trial vectors (=‘modes’) which stem from the FEM and are stored column-wise in the matrix $\Phi$.
One conceivable approach would be to treat the shaft and the hub as separate bodies. In that case it is not possible to compute the deformations in the contact area by a highly efficient matrix-vector multiplication because both bodies have their own nonlinear rigid body motion which must be considered. This is explained in detail in Section 2 of [20]. This holds even though the two bodies are connected by an interference fit. The use of two separate bodies can therefore be considered as less interesting in terms of efficiency. It is, of course, conceivable to somehow couple the two rigid body motions with a kind of constraint. This would also represent a general approach; however, it has the disadvantage that enforcing a constraint between the global rotation angles of two flexible bodies requires some implementation effort in existing software. In contrast, the approach presented below, is nonintrusive. Another approach would be to model both bodies together as one FE model without any interaction. This would lead to 12 rigid body modes which can be deactivated in the multibody simulation. In this case, only the modes corresponding to the elastic deformations remain. For further considerations, it must be considered that these modes occur independently of each other, as there is no interaction between the two bodies. A small thought experiment can be taken from Fig. 2 to illustrate the problems of this approach. It is assumed that both the shaft (black) and the hub (blue) have been modeled using the approach just outlined and the deformation shown in the figure occurs. In this particular deformation state, the hub would make a small rigid body movement with respect to the shaft. Since this is not possible, a solver would try to use elastic modes for the missing rigid body motion. It can be expected that the result will be poor. This is demonstrated in Section 3.3.4.

Rigid body motion between shaft and hub must therefore be possible. Due to the properties of an interference fit and the assumption of small elastic deformations within the framework of the FFRF, it is to be expected that these relative rigid body movements will remain small. Therefore, a modeling technique is required that allows relative rigid body motion between shaft and hub but treats both bodies as one. This can be achieved with a kind of trick. Basically, the two bodies are modeled in one model, but independently of each other, without interaction. The shaft and hub are then joined together with very soft springs. The springs should be attached in places where they do not have any influence. This modeling technique results in three groups of modes. The first group contains six pure rigid body modes with frequencies close to zero. The second group holds six rigid body modes of the two parts relative to each other. The frequencies are very low, but higher than the pure rigid body modes. A good choice of spring stiffness would be frequencies around 1Hz for this group of modes. Then, with significant frequency distance, comes the third group of elastic modes. It is crucial that the second group remains active in the multibody simulation. Since the frequencies are different from the first group, its selection is simple and because the frequencies are significantly below the first elastic mode, they have no effect on its behavior. One could argue that a large frequency difference to the first elastic mode may not be possible. However, this seems difficult to imagine for metallic shaft-hub connections. The first elastic frequencies of shaft or hub are certainly well above 1Hz. Fig. 5 in Section 3.2.1 contains images of all three mode groups for the example under consideration.

# 3 Numerical Example

## 3.1 Finite Element Model

The FE model has been created using ABAQUS® [21]. Images of the model can be found in Fig. 3. On the left-hand side, a wireframe representation of the entire model can be seen, while on the right-hand side, two solid representations of the shaft and the pulley are depicted. The blue lines in the image on the left-hand side symbolize the soft springs that lead to the required relative rigid body modes (RRBM) between the shaft and the pulley. The model consists of 218,226 elements and 320,384 nodes. Ten-node quadratic tetrahedral volume elements (C3D10) were used to discretize the structure. The interference fit was modeled with an oversize of 0.025mm. In the contact area, 9792 coincident nodes have been modeled, see Fig. 3 for the quality of the mesh.

Three distributed couplings were modeled for the bearing and load application. The position of these couplings is shown in Fig. 4. Two couplings link the nodes on the surface of the shaft to a middle node (A and B). These are mainly used to support the structure. A further coupling connects the nodes on the circumference of the pulley to a middle node C. This is mainly used for load application. As already mentioned in the theory section, it is necessary to couple the shaft and the pulley with soft springs. These lead to RRBMs that enable small relative rigid body movements in the numerical simulation. As mentioned above, these are shown in Fig. 3 as blue lines. For this purpose, four small distributed couplings around the circumference were modeled in the belt pulley. Both the center node A and the center node B have been connected to the four center nodes of these distributed couplings. The eight springs were modeled with an element that connects all six degrees of freedom of two nodes with linear springs. The translational spring stiffness was set to 0.01N/mm and the rotational spring stiffness was set to 0.01Nmm/rad.
<figure><img src="image_5.png" /><figcaption>Fig. 3: Images of the Finite Element model</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4: Distributed couplings for mounts and load application</figcaption></figure>
In the case of the computation of the contact modes around the shrink state, a nonlinear static computation must be performed. The result of this corresponds to state $x_0$ in the theory section and is taken into account as an additional mode that guarantees a very accurate representation of the shrinkage. To ensure a fair comparison, the same contact and friction models were used in both the FE code and the numerical time integration. The corresponding settings are as follows:

* Small sliding contact

* Node to surface contact

* Nonlinear contact pressure model with a penalty stiffness of 1.0e+5N/mm3 in case of closed contact and zero stiffness in case of open contact.

* A Coulomb-like friction model has been used with a friction coefficient of 0.2 and a maximal allowable relative motion of 0.01mm before sliding occurs.

## 3.2 Multibody Simulation

This work is a result of a research project using the multibody simulation software SIMPACK® [22]. The computation of the contact modes and the consideration of the nonlinear contact forces in the numerical time integration were made available in preliminary test versions. Following successful tests, it is planned to release these features for customers.

### 3.2.1 Extension of the Craig-Bampton Mode Basis with Contact Modes

In a first step, a Craig-Bampton [12] mode basis is computed. For this, the number of fixed-boundary normal modes and the interface degrees of freedom must be specified. For this work, 20 vibration modes and all six degrees of freedom of the nodes A, B, and C are specified for the Craig-Bampton reduction. The deformation state due to shrinkage (see Section 3.1) is used as an additional trial vector. This results in a total of 39 trial vectors (= modes).

Once this mode basis (still without contact modes) has been imported into the multibody simulation, the effect of the soft springs can already be evaluated. The first six modes should be rigid body modes (RBM) of the overall system with frequencies close to zero. This should be followed by a group of six relative rigid body modes (RRBM) characterized by relative movements between shaft and pulley. The frequencies should be as low as possible, but high enough so that this group can be clearly distinguished from the first group. The first elastic mode should then follow with the largest possible frequency spacing. The first group of RBMs is deactivated in the multibody simulation. According to the FFRF, the multibody simulation provides the degrees of freedom for the rigid body motion. The next group of RRBM must not be deactivated under any circumstances. It enables small rigid body motion between shaft and pulley and is essential for the quality of the results. This will be demonstrated in one of the later sections. For the present example, exemplary modes are given in Fig. 5. The upper row shows three of the six rigid body modes of the first group. The frequency ranges from 0.004Hz to 0.006Hz. It can be seen that the modes represent rigid body motion of the entire structure. Since the rigid body modes are not computed with the global FE mass and stiffness matrices, but on the basis of an already reduced model, they are not exactly zero. The lower row shows two of the six relative rigid body modes ranging from 0.36Hz to 1.48Hz. The pictures indicate relative motion between the pulley and the shaft. Finally, the first elastic mode is depicted having an eigenfrequency of 1,777Hz.

The contact mode forces are now computed according to the theory outlined earlier. As mentioned, there are two variants that are examined in this paper:
Variant 1: The first- and second-order contact pressure derivatives are fed together into a single PCA. The vectors with the highest singular values are then used as contact modes.

Variant 2: The first-order and second-order contact pressure derivatives are each fed to separate PCAs. The mode basis is first extended by vectors according to the first-order PCA and then by the vectors of the second-order PCA.

After the contact mode forces have been exported, a new FE computation is invoked. The contact mode forces are added to the former input file. The result of the FE computation is a mode basis with the contact modes, which can be re-imported into the multibody simulation. In addition to the two previously mentioned groups of rigid body modes, the Craig-Bampton modes with a more global deformation character and the contact modes with a more local deformation character can now be found in the mode basis. This is illustrated in Fig. 6 using randomly selected modes. The upper row contains the modes that can be assigned to the Craig-Bampton method, whereas the bottom row shows examples of contact modes. It is good to see that contact modes add local compliance in the potential contact region.
<figure><img src="image_7.png" /><figcaption>Fig. 5: Rigid body modes, relative rigid body modes and first elastic mode</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6: Craig-Bampton modes and contact modes</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 7: Contact pressure due to shrinking and bending</figcaption></figure>

### 3.2.2 Nonlinear Contact and Friction Forces

The computation of the contact and friction forces during the numerical time integration is well documented, see Pichler [10] and in great detail his doctoral thesis [9]. For this reason, only the most important key points are mentioned here:

* A small sliding contact is assumed.

* The contact kinematics is computed according to a node-to-surface formulation.

* The relative nodal displacements in the contact area are calculated from the modal coordinates. The corresponding contact and friction force laws are then applied. Finally, the contact and friction forces are projected back into the modal space.

The specific contact and friction models as well as the selected parameters are the same as in the finite element model, see Section 3.1. Note that the contact and friction model implementation is not exactly identical in the FE and multibody simulation software. However, the results presented later indicate that the contact model is almost identical and the friction model is very similar.

All multibody simulations follow the same pattern. First, the shrinkage process is simulated in a static computation. The resulting deformation is applied to the model as an initial displacement for all subsequent dynamic simulations.

## 3.3 Static Bending

In order to be able to compare the results of the proposed method directly with the FEM, static computations were carried out. A comparison with the Finite Element Method (FEM) in the case of dynamics has not been considered because of the huge computational burden.

### 3.3.1 Introduction

For these computations, nodes A and B (see Fig. 4) were fixed to ground. A force acts in the y-direction at node C. In both FEM and multibody simulation, the shrinkage process was computed first and then the load was applied. Two computations were carried out, one for 40kN load and one for 100kN. Fig. 7 shows how the contact pressure changes according to the load.

All computations were carried out with $C = 25, 50, 75, 100, 150$ and $200$ contact modes. In this way, it is possible to determine how many contact modes are necessary for a sufficiently accurate solution.

For quantifying the method's accuracy, an error $p_i$ in each of the $P = 9792$ contact pairs is computed along

$$p_i = \frac{|p_{FEM,i} - p_{MBS,i}|}{|p_{FEM,i}|} 100\% \qquad (14)$$

The value $p_{FEM,i}$ stands for the contact pressure according to the FEM at the i-th node pair ($1 \le i \le P$). The value $p_{MBS,i}$ corresponds to the same quantity computed with contact modes in the multibody simulation. Note that the error
<figure><img src="image_10.png" /><figcaption>Fig. 8: Error of contact pressure due to shrinking only</figcaption></figure>

$p_i$ is given in percent. Finally, the mean value

$$p_{\text{mean}} = \frac{1}{P} \sum_{i=1}^{P} p_i \qquad (15)$$

the standard deviation,

$$p_{\text{sd}} = \sqrt{\frac{1}{P} \sum_{i=1}^{P} (p_i - p_{\text{mean}})^2} \qquad (16)$$

and the maximum error

$$p_{\text{max}} = \max(p_i) \qquad (17)$$

are used as quality criteria.

### 3.3.2 Results Due to Shrinking Only

Fig. 8 shows the error due to shrinkage without bending. As expected, the error is very small because the deformation due to shrinkage is part of the mode basis. The fact that the error is nevertheless not exactly zero is due to the interplay between the different solvers and the slightly different formulations for the contact and friction forces.

### 3.3.3 Results Due to Shrinking and Bending

Fig. 9 shows the error measures when the bending load is applied after shrinkage. The left column contains the results for 40kN load and the right column for 100kN. The following conclusions can be drawn from the depicted curves:

* The result converges to the FE solution with an increasing number of contact modes.

* Variant 2 (separate PCAs for the first- and second-order derivatives) converges faster. With this variant, fewer contact modes are required for an acceptable solution quality.
<figure><img src="image_11.png" /><figcaption>Fig. 9: Error of contact pressure due to shrinking and bending</figcaption></figure>

* The errors for the 100kN load case are larger. This is also to be expected, as the contact pressure is further away from the development point of the Taylor series $p_0$.

* All in all, 75 or 100 contact modes seem to provide an acceptable solution quality when the contact modes due to Variant 2 are used. Both variants lead to an accurate solution when 150 contact modes are considered.

Now it is checked whether the converged solution also produces good results with regard to the tangential direction in which the frictional forces act. Unfortunately, it is not possible to directly compare friction forces with the standard output options of the FE and multibody simulation software. However, this is possible with tangential displacements. If these are similar, then there is no reason why the same friction model should not deliver the same result. Fig. 10 contains a comparison of the two tangential displacements in the contact area. For the mode-based simulation, 100 contact modes of Variant 2 have been used. It can be seen that there is good agreement.

For the 100kN bending load case 180 seconds wall-clock time are required when 100 contact modes are used. This simulation has been performed using SIMPACK® with one CPU. For the FE computations with ABAQUS®, 9210 seconds are required with one CPU and 2987 seconds with four CPUs. All computations have been performed on the same computer (Windows PC, 32 GB RAM, Intel 3820 CPU).

### 3.3.4 Importance of Relative Rigid Body Modes (RRBM)

In the theory section, it was stated that small rigid body movements between the parts must be possible in order to achieve good results. To make this possible without losing the efficiency of a single flexible body, the soft springs were introduced. Fig. 11 shows numerical evidence for this strategy. All three stress distributions result from the 100kN load case. The left picture corresponds to the result of the FEM. The center image is a result with 100 contact modes of Variant 2 including the RRBM. The right image is the result of a computation with 100 contact modes and deactivated RRBM. It can be seen that the RRBM have a significant effect on the quality of the results.

## 3.4 Static Torsion

This section analyses and compares the transition from local sticking to local sliding. Local sticking means that the local frictional stress is less than the local contact pressure multiplied by the static friction coefficient. In the other case, there is local sliding. The transition from sticking to sliding therefore depends on both the contact and the friction force. Correct reproduction of this transition indicates proper contact and friction forces.

For this purpose, a moment around the z-axis is applied at node C. Nodes A and B are fixed to ground. In a quasi-static computation, the torque is increased in 500Nm steps from 0 to 3,500Nm. At the beginning, directly after shrinkage, there is only sticking in the entire interference fit. As the torque increases, areas within the interference fit begin to slide. As the load continues to increase, these areas become larger and larger until there are no more
<figure><img src="image_12.png" /><figcaption>Fig. 10: Slip in both directions: FEM and multibody simulation with 100 contact modes</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11: Contact pressure due to shrinkage and 100kN load. Importance of relative rigid body modes (RRBM)</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 12: Transition from sticking to sliding with increasing torque</figcaption></figure>

areas of sticking. In this case, the interference fit would fail, and a rigid body rotation between shaft and pulley occurs. The comparison with the FEM is performed with 100 contact modes, as the former section indicates that this number delivers accurate results.

Fig. 12 shows how the sliding zones expand with increasing torsion load. The top row contains the FEM results and the bottom row the results of the contact-mode-based multibody simulation. Up to 1,000Nm there is no slipping in either case. At 3,500Nm there are no more sticking areas in both cases. In between, the results are very similar, although the sliding zones are slightly overestimated in the contact-mode-based multibody simulation. In view of the sensitivity of this parameter, however, this can be considered as a good result.

## 3.5 Dynamics

The latter comparison with the FEM was limited to static loads for the reason of acceptable computational cost. A dynamic computation is many times more complex and therefore it was previously considered as practically impossible to numerically integrate an interference fit with a fine FE mesh, nonlinear contact and friction forces and taking into account the full dynamics.

This section demonstrates that the presented method makes efficient dynamic computations with fine meshes possible. For this purpose, a space-fixed 50kN bending load is applied at node C (see Fig. 4) after the interference fit has been activated. At nodes A and B, revolute joints were applied between the shaft and ground. The shaft rotates at 1,000rpm. After approximately two revolutions, the force is abruptly set to zero and the rotating shaft begins to vibrate. The contact pressure at a selected node (with label 20315, see Fig. 13) is evaluated. The computation is carried out with different numbers of contact modes (25, 50, 75, 100, and 150). All modal damping values are set to 10% of critical damping. The SODASRT2 integrator has been used with an absolute tolerance of 1e-5 and a relative tolerance of 1e-7. SODASRT2 is a BDF integrator (backward differentiation formula) that was originally based on DASRT [23] but includes a number of significant (and non-published) improvements that make highly efficient multibody simulation possible. The simulation covers 0.2sec and the output frequency is 50,000Hz. The computations have been performed using SIMPACK® 2025 on a Windows PC with 32 GB RAM having an Intel 3820 CPU.

Fig. 14 shows the contact pressure curves of the selected node when 0, 25, 50, 75, 100 and 150 contact modes according to Variant 2 (separate PCAs for 1st and 2nd order derivatives) are used. It can be seen that the space-fixed load is deactivated after approximately 2 revolutions. The lower image contains a detailed view in the time from 0.149sec to 1.555sec. The following conclusions can be drawn:

* The solution converges with an increasing number of contact modes.

* For this load case, 75 contact modes deliver an acceptable result quality. It is interesting to see that this is a similar number to that evaluated for the static bending.
<figure><img src="image_15.png" /><figcaption>Fig. 13: Contact pressure in N/mm2 due to shrinkage and 50kN bending load. Location of FE node 20315</figcaption></figure>

<table><caption>**Table 1:** Error measures for 100 contact modes (computed at an interference 0.025mm) and different applied interference values</caption><thead><tr><th rowspan="2"></th><th colspan="3">Shrinking Only</th><th colspan="3">Shrinking and 40KN Bending</th><th colspan="3">Shrinking and 100KN Bending</th></tr><tr><th>p<sub>mean</sub>[%]</th><th>p<sub>sd</sub>[%]</th><th>p<sub>max</sub>[%]</th><th>p<sub>mean</sub>[%]</th><th>p<sub>sd</sub>[%]</th><th>p<sub>max</sub>[%]</th><th>p<sub>mean</sub>[%]</th><th>p<sub>sd</sub>[%]</th><th>p<sub>max</sub>[%]</th></tr></thead><tbody><tr><td>15µm</td><td>0.3</td><td>0.2</td><td>1.6</td><td>0.7</td><td>0.6</td><td>4.3</td><td>1.7</td><td>1.9</td><td>12.6</td></tr><tr><td>20µm</td><td>0.2</td><td>0.1</td><td>0.9</td><td>0.5</td><td>0.4</td><td>3.2</td><td>1.2</td><td>1.1</td><td>8.7</td></tr><tr><td>25µm</td><td>0.1</td><td>0.1</td><td>0.4</td><td>0.4</td><td>0.3</td><td>2.4</td><td>1.1</td><td>0.9</td><td>6.7</td></tr><tr><td>30µm</td><td>0.2</td><td>0.1</td><td>0.6</td><td>0.4</td><td>0.3</td><td>2.5</td><td>1.0</td><td>0.8</td><td>5.7</td></tr><tr><td>35µm</td><td>0.2</td><td>0.2</td><td>1.1</td><td>0.5</td><td>0.3</td><td>2.7</td><td>1.0</td><td>0.8</td><td>6.0</td></tr></tbody></table>

* The wall-clock time for the simulations with 100 contact modes was 733 seconds. Considering that such simulations have not been possible before, this is a very short simulation time.

## 3.6 Usability of a Mode Basis for Different Interference

In order to obtain an accurate result with as few contact modes as possible, (a) the deformation due to shrinkage is included in the mode basis and (b) the derivatives are formed around this state. The consideration of the deformation due to shrinkage requires at least one nonlinear FE computation. The question arises whether a new mode basis is necessary if existing contact modes and another interference are taken into account. In other words: How sensitive is the quality of the result if the actual interference does not correspond to the interference that was considered for the calculation of the mode basis? A mode basis with 100 contact modes (Variant 2: separate PCAs for 1st and 2nd order derivatives) is used for this investigation. The additional mode due to shrinkage was computed at an interference of 0.025mm. The modal derivatives were also formed around this state. The bending load case from Section 3.3 was computed with 40kN and 100kN for an interference of 0.015mm, 0.02mm, 0.025mm, 0.3mm and 0.035mm. The errors defined in Section 3.3 are given in Table 1. Note that the first column contains the result when shrinking is applied without bending.

A closer look at the results outlined in Table 1 reveals that the error increases as the bending load increases. This is to be expected. The error also tends to increase when the difference between the actual interference and the interference in the computation of the mode basis increases. This is also not surprising.

Although the interference was varied significantly, the errors remain relatively small. It therefore appears to be possible to vary the interference without the necessity of computing a new mode basis (and thus also a nonlinear FEM computation). This brief investigation naturally raises interesting questions for future research. One such question, for example, would be: Can a conservative estimate for a maximum permissible deviation be given?
<figure><img src="image_16.png" /><figcaption>Fig. 14: Dynamic load: pressure at node 20315</figcaption></figure>
# 4 Conclusion

In this work, a new strategy for the efficient time integration of transient processes of FE models with interference fits is presented. The contact and friction forces are considered nonlinearly, the FE mesh can be fine, and no simplification of the inertia forces due to vibrations and nonlinear rigid body motion is assumed. The key technology herein is the use of contact modes, which can be easily computed and enable the accurate representation of the deformations inside a contact zone in a reduced model. For the applicability of contact modes to interference fits, two issues were discussed in the theory section:

* It has been shown that contact modes based on modal derivatives are equivalent to contact modes based on a Taylor series development. The second one enables a simpler computation of contact modes when the operating point is not the undeformed configuration, as is the case with interference fits.

* A highly efficient computation requires self-contact. Therefore, a strategy is proposed that results in the shaft and hub being one body, but still allows small, relative rigid body motion between these two parts.

The conclusions and findings from the theoretical foundation were then numerically tested on a shaft-pulley connection. The following was observed and verified:

* Comparisons were made with the FEM (reference solution) for static bending and torsional loads. In all cases, the contact quantities (contact pressure, frictional stress, status) correlate well with the FEM solution when an increasing number of contact modes is used. The number of contact modes is much smaller than the number of node pairs in the contact area.

* A dynamic load case was investigated by continuously increasing the number of contact modes. Convergence was observed and the number of necessary contact modes corresponds well with the static load case.

* The numerical comparison shows that fewer contact modes are needed if two separate PCAs are performed, one for the first-order and one for the second-order derivatives.

* It was also shown that the results remain accurate, even if the interference is varied without recomputing the modes and, thus, also the nonlinear shrinkage.

Considering that the dynamics of finely meshed interference fits have not been previously computed economically, a simulation time in the range of minutes represents a significant step forward for such investigations.

Two important topics were not addressed in this publication. One is the consideration of temperature fields, the other the consideration of local plastification. Both are important future issues. The consideration of temperature fields with additional trial vectors in the mode basis seems plausible, as studies have already been published, for example [24] and [25]. As far as we know, there is nothing in literature on the consideration of local plasticizing with trial vectors. Nevertheless, this does not seem to be a completely erroneous approach, because plasticizing inside interference fits should be relatively smoothly distributed over a certain subarea of the contact zone.

# Authors' Contributions

The theory is based on ideas of Witteveen, Koller and Pichler. The first author is responsible for the structure of the work, for the numerical simulations, for the evaluation of the results and for the writing of this paper. Müller and Dietz are responsible for the implementation of the contact mode computation and the numerical time integration in ABAQUS® and SIMPACK®. Koller supported the project with consulting activities and practical work related to the FE model.

# References

[1] Y. Zhang, B. McClain, and X. Fang. Design of interference fits via finite element method. *International Journal of Mechanical Sciences*, 42(9):1835–1850, 2000. ISSN 0020-7403. doi:https://doi.org/10.1016/S0020-7403(99)00072-7.

[2] F. Lanoue, A. Vadean, and B. Sanschagren. Finite element analysis and contact modeling considerations of interference fits for fretting fatigue strength calculations. *Simul. Model. Pract. Theory*, 17:1587–1602, 2009. doi:https://doi.org/10.1016/j.simpat.2009.06.017.
[3] B. Xinxiao, W. Chen, and L. Yuanrong. Research on the fretting contact fatigue strength of interference fit in high speed rail traction motor shafts. *Engineering Failure Analysis*, 157:107875, 2024. ISSN 1350-6307. doi:https://doi.org/10.1016/j.engfailanal.2023.107875.

[4] J. Madej and M. Śliwka. Analysis of interference-fit joints. *Applied Sciences*, 11(23), 2021. doi:10.3390/app112311428.

[5] N. Ajavakom and A. Sillapapinij. The effect of interference fit on vibration transmission from stator coil to base of a spindle motor in a hard disk drive. *Journal of Research and Applications in Mechanical Engineering*, 1(1):25-29, 2018, URL https://ph01.tci-thaijo.org/index.php/jrame/article/view/150707.

[6] Z. Wang, Z. Wang, X. Bai, X. Zhang, and Y. Wang. Effect of interference fit on dynamic characteristics of spindle rotor system. *Journal of the Brazilian Society of Mechanical Sciences and Engineering*, 44(316), 2022. doi:https://doi.org/10.1007/s40430-022-03545-4.

[7] G. Liu, J. Hong, W. Wu, and Y. Sun. Investigation on the influence of interference fit on the static and dynamic characteristics of spindle system. *Int J Adv Manuf Technol*, 99:1953-1966, 2018. doi:https://doi.org/10.1007/s00170-018-2567-8.

[8] W. Witteveen and F. Pichler. Efficient model order reduction for the dynamics of nonlinear multilayer sheet structures with trial vector derivatives. *Shock and Vibration*, 2014, 2014. doi:https://doi.org/10.1155/2014/913136.

[9] F. Pichler. *Multibody dynamics of jointed flexible structures*. Ph.D. Thesis, TU Graz, Institut fuer Fahrzeugtechnik, 2018.

[10] F. Pichler, W. Witteveen, and P. Fischer. A complete strategy for efficient and accurate multibody dynamics of flexible structures with large lap joints considering contact and friction. *Multibody System Dynamics*, 40(4):407-436, 2017. doi:10.1007/s11044-016-9555-2.

[11] J. B. Rutzmoser. *Model Order Reduction for Nonlinear Structural Dynamics Simulation-free Approaches*. Ph.D. Thesis, Technical University Munich, Chair of Applied Mechanics, 2017.

[12] R. R. Craig and M. C. C. Bampton. Coupling of Substructures Using Component Modes. *AIAA Journal*, 3(4):678-685, 1968. doi:10.2514/3.4741.

[13] R. R. Craig and C.-J. Chang. Substructure coupling for dynamic analysis and testing. NASA-CR-2781, page 89, 1977, URL https://ntrs.nasa.gov/api/citations/19770010568/downloads/19770010568.pdf.

[14] P. Castanier, M., C. Pierre, and Y.-C. Tan. Characteristic constraint modes for component mode synthesis. *AIAA Journal*, 39: 1182-1187, 2001. doi:https://doi.org/10.2514/2.1433.

[15] W. Witteveen, M. Kuts, and L. Koller. Can transient simulation efficiently reproduce well known nonlinear effects of jointed structures? *Mechanical Systems and Signal Processing*, 190, 2023. doi:10.1016/j.ymssp.2023.110111.

[16] F. Pichler, W. Wolfgang, and P. Fischer. Reduced-order modeling of preloaded bolted structures in multibody systems by the use of trial vector derivatives. *Journal of Computational and Nonlinear Dynamics*, 12:051032, 2017. doi:https://doi.org/10.1115/1.4036989.

[17] A. A. Morsy, M. Kast, and P. Tiso. A frequency-domain reduced order model for joints by hyper-reduction and model-driven sampling. *Mechanical Systems and Signal Processing*, 185:109744, 2023. ISSN 0888-3270. doi:https://doi.org/10.1016/j.ymssp.2022.109744.

[18] C. Gastaldi, S. Zucca, and B. I. Epureanu. Jacobian projection reduced-order models for dynamic systems with contact nonlinearities. *Mechanical Systems and Signal Processing*, 100:550-569, 2018. ISSN 0888-3270. doi:https://doi.org/10.1016/j.ymssp.2017.07.049.

[19] A. A. Shabana. *Dynamics of Multibody Systems, 5th Edition*. Cambridge University Press, 2020. ISBN 1108485642.

[20] W. Witteveen and H. Irschik. Efficient mode-based computational approach for jointed structures: Joint interface modes. *AIAA Journal*, 47(1):252-263, 2009. doi:https://doi.org/10.2514/1.38436.

[21] URL https://www.3ds.com/products/simulia/abaqus.

[22] URL https://www.3ds.com/de/produkte-und-services/simulia/produkte/simpack/.

[23] K. E. Brenan, S. L. Campbell, and L. R. Petzold. *Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations*. Society for Industrial and Applied Mathematics, 1995. doi:10.1137/1.9781611971224.

[24] S. Jain and P. Tiso. Model order reduction for temperature-dependent nonlinear mechanical systems: A multiple scales approach. *Journal of Sound and Vibration*, 465:115022, 2020. ISSN 0022-460X. doi:https://doi.org/10.1016/j.jsv.2019.115022.

[25] H. Yamashita, R. Arora, H. Kanazawa, and H. Sugiyama. Reduced-order thermomechanical modeling of multibody systems using floating frame of reference formulation. *Proceedings of the Institution of Mechanical Engineers, Part K*, 233(3):617-630, 2019. doi:https://doi.org/10.1177/1464419318810886.