On an improved adaptive reduced-order model for the computation of steady-state vibrations in large-scale non-conservative systems with friction joints

Joints are commonly used in many large-scale engineering systems to ease assembly, and ensure structural integrity and effective load transmission. Most joints are designed around friction interfaces, which can transmit large static forces, but tend to introduce stick-slip transition during vibrations, leading to a nonlinear dynamic system. Tools for the complex numerical prediction of such nonlinear systems are available today, but their use for large-scale applications is regularly prevented by high computational cost. To address this issue, a novel adaptive reduced-order model (ROM) has recently been developed, significantly decreasing the computational time for such high fidelity simulations. Although highly effective, significant improvements to the proposed approach is presented and demonstrated in this paper, further increasing the efficiency of the ROM. An energy-based error estimator was developed and integrated into the nonlinear spectral analysis, leading to a significantly higher computational speed by removing insignificant static modes from the stuck contact nodes in the original reduced basis, and improving the computational accuracy by eliminating numerical noise. The effectiveness of the new approach was shown on an industrial-scale fan blades system with a dovetail joints, showing that the improved adaptive method can be 2–3 times more computationally efficient than the original adaptive method especially at high excitation levels but also effectively improve the accuracy of the original method.


Introduction
The assembly of single components into a more complex structure always leads to the presence of mechanical joints. Most of mechanical joints are realized through friction interfaces, which are often regarded as the main source of non-conservative nonlinearities in an assembled [6]. Other types of nonlinearities can exist such as geometric nonlinearities, material nonlinearities and aerodynamic forces, but they still remain less important than friction joints, especially for the applications to turbo-machinery [1]. The contact friction can result in a significant reduction in the global stiffness of the whole structure through the relative motion on the joint interface, which also introduces amplitude-dependent damping [11,16,19,28]. For an accurate prediction of the dynamic behavior of the assembly, the dynamics of the joint of interest must be captured accurately [6]. In spite of the advance in com-putational capability, these inherent nonlinearities are often ignored or linearized in the simulation because of the large computational expense and numerical instability problems [35]. One feasible approach to incorporate the contact friction nonlinearities represented by physically constitutive contact friction models into the large-scale jointed structures is to employ ROM techniques. In fact, they have been effectively proved that it can reduce the size of associated high fidelity models by ten even hundred orders with high accuracy [5].
Component mode synthesis-based ROMs are traditionally used for model-order reduction for the assembly with contact friction interfaces [3,6,15,20,35]. Its main idea is to consider the components in an assembly as distinct sub-structural elements. The reduced basis for such an assembly can be then constructed with the linear normal modes and the static modes associated with the retained interface nodes from the sub-structural elements [38]. With this set of reduced basis, Galerkin Projection is then performed for modelorder reduction [23,35]. The subsequent reduction can effectively eliminate most internal DOFs for each component while retaining the whole DOFs on the joint interfaces [39]. For example, Rubin method and Craig-Bampton (CB) method are two of most efficient methods in a class of the CMS technique [8,25]. The benefit of these CMS methods is that it is convenient to integrate nonlinear contact constitutive model on the interface. It however makes the size of final ROMs highly depend on the size of contact interface, especially for a large-scale system when contact interfaces are intensive or meshed in a high density in order to recover the stress level [31,35]. This would greatly limit the computational efficient of classic CMS methods.
As a result, the concept of interface reduction was put forward to condense static impulse modes [13]. The idea is to replace full size of static modes with a much smaller subset of eigenvectors from a second modal analysis of the system projected by interface static modes. Becker and Gaul applied a common interface reduction technique to jointed structures in order to optimize the classical CB constrain modes [4]. Witteveen proposed another interface reduction approach which is referred as joint interface mode (JIM) method [31]. JIM is obtained also from an additional modal analysis but from a statically condensed system through respecting the Newton's third law on the joint interface. A comparative study of two method shows JIM method has superior convergence rate over the com-mon interface reduction method [4]. This is mainly because JIM method allows more local flexibility on joint interface leading to a better description of the relative motions on the joint interface. In split of a smaller size of interface modes, these methods all require a second eigenvalue analysis that would further increase the off-line cost in model-order reduction. Witteveen also proposes another approach called trial vector derivative (TVD) method to improve the convergent rate in ROM [22,32]. The idea is to enrich the reduced subspace by accounting the effects of nonlinearities on the linear CMS subspace. This TVD method is essentially based on the modal derivatives which is obtained by the firstorder Tayler expansion to the linear CMS subspace. The study shows the TVD method can achieve a better convergence rate than the JIM method. This modal derivative technique is also suitable for ROMs in geometrical nonlinear problems [30,33]. However, for the contact friction problems, it would further increase the size of the reduced system and the associated modal derivatives can be very much pre-loading dependent [22].
To overcome some of the above challenges, the adaptive CMS method was developed [36,37] that allowing the set of static modes update adaptively according to the varying contact conditions on the interface during the on-line computation. This adaptivity can lead to a considerable computational saving for jointed structures with micro-slip motions. The essence of this method consists in a re-written equation of motion describing the assembled structure with an underling linearized system and an adaptive internal variable that accounts the nonlinear effects from contact interface. Such a formulation would allow the subsequent ROM to remove a number of redundant static modes associated with sticking nodes adaptively. This adaptive ROM can also be conveniently integrated with the popular harmonic balance method (HBM) to obtain the most interesting steady-state response when the system is subjected to the periodic excitation. The study show this adaptive CMS method can make the simulation 50-100 times faster than the classic CMS-based method. It provides a very promising and appealing approach for the application to largescale joint assemblies when the contact interface is enormous or meshed in a high density. Although the original adaptive method works very well, there is still room for improvement. The original adaptive method [36] includes all the static modes associated with the stick-slip nodes on the contact interface. Many of these stick-slip nodes are actually in a sticking dominated mode. It means they contribute very tiny amount of energy dissipation which does not affects the dynamical behavior. These related static modes can be further removed from adaptive reduced-order model with a proper error estimator. Alternatively, some of the nonzero small nonlinear penalties in the original adaptive method might be not actually related to the partially sliding nodes but due to the errors from automatic size updating algorithm. Including these penalties however would compromise the accuracy of this adaptive CMS method.
Given that, this work proposes an improved adaptive ROM technique to enhance the computational efficiency and accuracy of the original adaptive method. An energy-based error estimator is developed and introduced into automatic size update algorithm to filter out those sliding nodes with little effects on the dynamic behavior. We consider an industrial-scale problem as the test case. A large-scale finite element fan blade system with dovetail joints is used. The classical CMSbased Rubin method and original adaptive method are used to benchmark the results obtained from the proposed method. The paper is organized as follows: firstly, the formulation of the adaptive ROM and the proposed error estimator is presented; it is followed by the detailed description of automatic size update algorithm with error estimator; after that, the fan blade FE model as well as its modal characteristics and the general contact mechanics in its joint will be described; the comparison of forced responses and associated computational cost between the proposed and benchmark methods is then presented and discussed followed by the conclusion. Figure 1 illustrates a general assembled structure consist of two deformable solids coupled by a friction joint. Joint contact interface between two solids is in red where the localized contact friction nonlinearity is introduced into the whole dynamical system. The equation of motion between these two deformable bodies with the friction joint can be discretized with finite element methods as:

Equation of motion
where M 1 , M 2 , K 1 , K 2 is the mass and stiffness matrix of these two deformable bodies; F 1 ex , F 2 ex is the harmonic external forced excitation on these two components, respectively; B 1 , B 2 is the signed Boolean matrix acting on the contact interface associated with the contact nodes on joint interface; F nl is the contact friction force vector on the joint interface, which is a function of the joint nodes displacement u 1 b , u 2 b .
2.2 3D node-to-node contact model The contact interfaces are discretized into dense mesh grids by using 3D node-to-node contact elements shown in Fig. 2 to have a realistic physical representation of contact friction behaviors. Each 3D friction element is made of two coupled Jenkins elements to model a two-dimensional in-plane motion, which are coupled with a spring in the normal direction. Normal contact conditions are defined as two states: in-contact mode and separation mode. If predicted normal force F N (t) becomes negative, the separation occurs between two contact points and the normal contact force becomes zero. Whereas if F N (t) is positive, the contact pair is in contact and thus there are two types of tangential contact states: one is the sticking mode and the other is the slipping mode. The two surfaces would remain stick if the total contact force (F T (t) and F N (t)) stays inside the friction cone. However, when the predicted friction force reaches the maximum static friction μF N (t), two surfaces would start to slip, and the energy is dissipated through the contact interface. Six parameters are used to characterize the properties of this element: the element area A, the friction coefficient μ, the normal static load N , the normal static gap g 0 and the tangential k t and normal k n contact stiffness (also known as 'penalty' coefficients). The detailed formulation of this 3D node-to-node contact elements can be found in [21,26]. There are also many friction rigs that are used to experimentally identify related contact parameters [9,10,18].

Modified equation of motion
The adaptive ROM approach is formulated based on a linearized system. The general EOM in Eq. (2) is linearized with contact stiffness k t and k n assuming the whole contact interface is in a sticking mode. The linearized EOM of the system with joint contact stiffness matrix can be expressed as: where K Joint is a full transformed global stiffness matrix of the local stiffness matrix K nl,nl Joint where 'nl' is associated with joint interface DOFs. The linearized stiffness matrix in Eq. (2) is denoted as K L .
The prediction from this linearized system is accurate if the whole contact interface is in a purely sticking condition. However, if any of the contact pairs is in a slipping or separation condition, the prediction might be not accurate. As such, an internal variable Δp is introduced into the Eq. (2). The full dimension of Δp is the half number of the total dovetail DOFs (3× N c ). N c is the number of contact nodes. The expression of the internal variable for i th contact pair in a local coordinate system can be formulated as: where Δu b is the assembly of relative displacement of all the contact pairs in joint interfaces. Δp i becomes a zero vector when the contact pair is in a full sticking condition. By adding the internal variable into linearized system, the final modified EOM for the adaptive approach can be written as: where K 1,nl Joint is the transformed joint matrix relating to the DOFs of Solid 1 and the internal variable Δp i . Similarly, K 2,nl Joint is the one associated to the Solid 2. K nl,nl Joint is the squared joint stiffness matrix. K nl,1 Joint , K nl,2 Joint is the transposed matrix of K 1,nl Joint , K 2,nl Joint . M G , K G are used to denote as the new assembled global mass and stiffness matrix.

Adaptive transformational matrix
The adaptive reduced-order model can be then constructed from the new reformulated dynamic system in Eq. (4). The reduced basis for the adaptive approach includes two types of reduced basis: one is the linear normal modes and the other is constrain modes related to the interface DOFs. The ROM can be therefore constructed with ψ and φ as: where φ are the modes of the linearized system from eigenanalysis of Eq. (2); ψ are the constrain modes obtained by imposing unit displacement on the DOFs associated with internal variable Δp; η is the modal participation factors of the selected dynamic modes; Δp R is the nonzero part of Δp; B is the Boolean matrix that to capture the nonzero part of Δp; is the overall transformation matrix for the adaptive reduced-order model. The size of depends on the nonzero part of Δp namely the contact conditions of the contact nodes. This would allow a significantly second model-order reduction to save computational cost. The evaluation of φ and ψ can be found in [36]. The adaptive ROM system is then obtained using the modal projection with the : where M R , K R is the reduced global mass and stiffness matrix.

Error estimator
As described in Sect. 2.4, the second reduction (from Δp to Δp R in Eq. (5)) in original adaptive method is achieved by the adaptive removal of zero terms in internal variable Δp. This process can be actually further improved on the following three aspects: -For the nonlinear vibrational computation of a large-scale system, numerical errors are always existing due to the use of iterative solvers to achieve pre-setting tolerance, especially for a large-scale system. Δp therefore would tend to include a lot of small nonzero terms due to such numerical errors. Using the original adaptive method, these unavoidable numerical noise would unnecessarily increase the size of nonzero part in the adaptive ROM, and that would result in a lower computational speed. Furthermore, it is difficult to set an explicit tolerance level beforehand within the original updating method to eliminate such numerical errors as it is pretty much case dependent. -The algorithm for real-time automatic size update [36] may still introduce some errors itself. This is because the contact conditions to update the reduced system is predicted from previous solution. Although acceptable, this could lead to some unavoidable updating errors. Using the original updating method, nonzero terms in the internal variable due to these errors would be regarded as the nonlinear penalties from the contact nodes. This would compromise the accuracy and computational speed of this adaptive ROM method. -The second reduction in original method (shown in Eq. (5)) retains all the sliding contact nodes in Δp R . Many of these sliding nodes are actually in a partially sliding (namely stick-slip) condition. Some of these sliding nodes seem to be very small energy dissipation on the interface leading to very little contribution to the overall dynamic performance. Therefore, removing the static mode associated with these node gives an additional opportunity to reduce the model further.
The function of error estimator in classical CMS method is to quantify to what degree each interface static mode influences the numerical error in the reduced solution [13]. It would allow for automatic control of the error in the reduced model and thus optimize suitable dimensions of the subspace associated with static modes. To overcome some of the above limitations in original adaptive ROM, a novel real time error estimator is proposed herein on the improvement of the second reduction from the Δp to Δp R in Eq. (5). The error estimator is based on the energy dissipation from the contact nodes in the joint interface. For i th in-contact joint node, the corresponding part in error estimator e i can be expressed from the integration of tangential friction force as: (7) where T is a period of vibrational cycle. The amount of energy dissipation e i indicates to what extent the associated static modes influence the vibrational behavior.
Each contact node has a corresponding part in error estimator that would be updated at each excitation frequency point. e i is supposed to be very close to zero of the ith contact pair is in a pure sticking mode. The associated static mode can be then removed from the reduced subspace without having much impacts on the dynamics. In this way, this error estimator is used as a physical indicator to remove or add static modes related to all the contact nodes adaptively. To be noted, in case of separation in a contact pair i, a large value is assigned to e i . The associated static modes for this pair would be retained. The Δp in Eq. (3) can be then updated using this error estimator expressed as: where ε is the energy dissipation level tolerance (unit: J) for each contact node. The energy dissipation level e i at the ith contact node less than the ε would be removed from the reduced basis . ε needs to be defined beforehand based on the allowable case-independent ratio of the energy dissipation level to the external energy level, e.g., 10 −6 . In this way, the reduced internal variable Δp R shown in Eq. (5) can be then further reduced to eliminate all the nodes that fall under this threshold.

Harmonic balance method
Multi-harmonic balance method (HBM) with alternating frequency time scheme (AFT) is used to evaluate the steady-state dynamic response of the system [14]. The basic idea of HBM is to decompose nonlinear dynamic response using truncated Fourier series with N H harmonic series: where q(t) is time-dependent dynamical response of the reduced system;Q n is the complex harmonic coefficients of n harmonic. AFT technique is used to switch the response between frequency and time domain in order to calculate the contact friction forces in time domain where the contact law is defined [15]. The continuation technique that include prediction and correction stages is employed to track the nonlinear frequency response with any soften or harden behavior. More details with respect to HBM, AFT and continuation techniques can be found in [15,27,35]. In the following case study, the Secant method is used for predictor while Arc-length method for the corrector.

Automatic size update with error estimator
The algorithm to achieve the automatic size update with error estimator is presented in Algorithm 1. The basic steps of this implementation are the same as that in original adaptive method [36] which would be briefly described in this section. The main difference from the original adaptive method in [36] is that the proposed error estimator is used in contact condition reevaluation process allowing to update the reduced system more efficiently. This part is highlighted in blue in Algorithm 1. The computation starts from the initial converged solution (Q o , ω o ) where the initial B o shown in Eq. (5) is full including all the joint DOFs. After the initialization, a while loop begins to compute the forced response for the excitation frequency between ω start and ω end .Q R is used to note the harmonic coefficient vector after the second adaptive reduction. For each new frequency, the classical continuation technique (include the prediction and iterative correction stage) is performed firstly to obtain the converged solution (Q R , ω c ). The converged solution is then expanded to ensure that all the solutions have the same dimension and saved.
It is then followed by the contact condition reevaluation to update B for the proposed adaptive ROM. The nonlinear force function coupled with AFT techniques is called to re-calculate F nl using the converged solution. It is worth noting that not only the predicted sliding or separation nodes, but all the contact nodes are evaluated. Different from the original algorithm, both the internal variable Δp and error estimator e are evaluated with the formulation in Eqs. (3) and 7 separately. It is followed by a further update of the Δp with e using Eq. (8) to remove the nodes with extremely low energy dissipation. After that, B new can be obtained with the updated Δp by retaining the nonzero terms. It is worth noting that such a contact condition re-evaluation is only performed once at a frequency. The additional computing cost is ignorable comparing to the overall computational cost.
After the re-evaluation, the reduced system needs to be updated for the next frequency point. However, if the new B new is the same as the B in last frequency, the reduced system remains the same and the computation jumps to the evaluation of the next frequency point. Otherwise, the reduced system M R , K R , needs to be updated according to the B new in Eq. (5). This iterative process continues until the tracing frequency reaches the end of defined excitation frequency range.
One of big advantages in this automatic updating algorithm is the classic continuation part in the HBM framework remains unchanged. The error estimator can successfully detect contact nodes that never slip but which were erroneously retained in the original adaptive reduced system. These identified nodes can then also be removed due to their very limited energy dissipation which has a minimal effect on the overall dynamic performance. This can further increase the computational efficiency despite the evaluation of error estimator would take some additional computation cost. Furthermore, the error estimator can potentially help to reduce the errors caused by updating errors. The following case study is going to evaluate how much improvement the new approach really is.

Case study: fan blade system
Fan blade system is the largest bladed disc assembly in a modern turbofan engine that provides the over 60% thrust to propel the aircraft [2,34,38]. This fan assembly consists of a number of blades and a disc on the rotating shaft. They are commonly assembled through a disk via curved or straight dovetail roots. Such a design ensures easy assembly and safe load distribution, and also provides essential damping to the system. The state-of-the-art linear vibration analysis often leads to an over-design of the components by ignoring the nonlinearities from the dovetail joint due to the complex and strong nonlinear dynamic nature in the friction joint. To further improve the fan blade root design, it is therefore necessary to take these nonlinearities into account for a better dynamic design [29]. However, performing nonlinear vibrational analysis is a well-known academic and industrial challenge partly due to the enormous computational expense [7]. This large-scale fan blade system is therefore used as the test case in this paper. Figure 3a shows a full scale fan blade system in turbofan aero-engine that includes more than two dozens of blades (in blue) and one disc (in red) connected via curved dovetails joints. The fan blade system is cyclic symmetric, and one sector of the assembly is considered for this study. Figure 3b shows a sector of the bladed disc assembly where the blade and disc are entirely made of Titanium considered as a homogeneous and isotropic material. The blade has a low slenderness aspect ratio, i.e., 4, and an increasing twist from root to tip. This high fidelity finite element model was built with the quadratic hexahedral element and each node has 3 DOFs. The model of each sector consists of a blade of 27,707 nodes and a cyclic symmetric part of the disc of 13,826 nodes. Figure 3c show the geometry of zoomed dovetail joint that connects the blade  The first four mode shapes of the fan bladed-disk and disc. The matching mesh is used in the joint contact interface between the blade and disc in order to employ the Jenkin element for accurate prediction of contact friction forces. The description of the contact element is detailed in Sect. 2.2. Figure 3d shows the root contact interfaces (in blue) between the blade and disc. For each dovetail joint, there are 208 contact node pairs on the joint interface in this model leading to 208 × 6 contact friction DOFs in total. Each of the contact pair has its own local coordinate system due to the curved contact interface. Figure 4 shows the first four linearized modes of one fan bladed disc sector where the whole joint surface is in a fully sticking condition. The relative spacing of the first four frequencies is 1:1.3:2:4. Due to the twist nature of fan blade, the type of the first four mode shape are very mixed. The first four mode shapes are dominated by out-of-plane flapping, edgewise and torsion motions. In-plane bending, and axial stretching do not appear in the first six modes. In this paper, the nonlinear effects of dovetail joint on the dynamics related to the first flapping mode is further studied for the assessment of proposed ROM approach. Figure 5 shows the contact loads and external forces on the dovetail joints in a fan blade system. The fan is mainly loaded by the centrifugal force T generated by the rotation of the shaft around the engine axis. The vibration motion on the blade is mainly derived from Fig. 6 Comparison of forced frequency response between adaptive ROM-and CMS-based Rubin method at different excitation levels the unsteady aerodynamic loading due to the rotation. The vibration is then transmitted to the disc through a closed-form connection via the dovetail joint that would produce high local contact normal and tangential stress fields P and Q on the contact surface [24]. In spite of the dovetail joint with high centrifugal loading, micro-slip motion is often present at the edges of the interface due to the vibrational load. This stick-slip zone in the joint allows for local slipping at the extremities of the contact zone. On the positive side, this microslip motion can dissipate the vibrational energy through the friction mechanism that can contribute to the reduction of maximum dynamic stress in the blades [12,17]. Alternatively, this would give rise to fretting fatigue problems leading to early cracking initialization as a result of a combination of low cycle centrifugal loading during taking off and landing and high cycle vibrational loading. This study is focused on the assessment of the effects of such a micro-slip motion in dovetail joints on the dynamic response with high fidelity FE model.

Results and discussions
This section is to demonstrate how well the improved ROM technique works with the fan blade system comparing to the original method and show whether the adaptive ROM method can be used for an industrial scale problem. The results using the original adaptive method is going to be present and analyzed firstly. It is because that the original method has not been applied to the large-scale structures yet. It would be therefore sensible to see firstly how the original adaptive method works for large models, before moving on to assess the improved method with error estimator in comparison. Figure 6 shows the FRFs of the fan bladed disc system using the original adaptive method and Rubin method at four excitation levels, namely 1 N, 3 N, 5 N and 8 N. The point force is applied to the node at the bottom of the fan blade close to the dovetail joint in the out-of-plane direction. The dynamic response is measured from the node in the middle of the fan blade system. The CMSbased Rubin method is used so that the effectiveness of the original adaptive method can be demonstrated. 3 harmonics are used in HBM and maintained for all the simulations. The pre-loading pressure distribution at the contact interface under the centrifugal loading is assumed to be uniform. It should be calculated from the nonlinear static analysis in a real engineering design. As a classical CMS method, the formulation of the Rubin method can be found in [25,35,36]. The resonance amplitude of forced response decreases gradually with the increase in the excitation levels. This is due to the increase in the energy dissipation from the dovetail joint. Comparing to the Rubin method, the original adaptive method can obtain very similar nonlinear dynamic response of fan blade disc at all the four excitation levels. However, there are still some discrepancy in the forced response between the original adaptive Fig. 7 The contact condition of interface nodes at different excitation levels (Red circle: stick-slip node) method and Rubin method. The maximum difference occurs in the frequency close to resonance at high excitation level of 8 N, which is around 2.2% lower than that from Rubin method. It means the original adaptive method overestimates the energy dissipation on the joint interface. This is due to the errors introduced from automatic size updating algorithm as discussed previously in Sect. 3. The adaptive method utilizes the predictions from the current convergent solution to update the ROM for the evaluation at the next frequency point. Although this method has approved accurate enough for most small scale systems [36], it may still lead to that some elements are eliminated that may want to slip next time round for a large-scale system. Figure 7 show the tangential contact condition of 208 contact pairs on the joint interface at each frequency point. The frequency point is corresponding to the frequency in forced response in Figure 6. The red circle indicates the node is in a stick-slip condition. It shows the number of stick-slip nodes increase significantly when the frequency approaches the resonance. As expected, the region where the stick-slip node locates expands on both sides of the resonance frequency with the increase in excitation levels. The increase in stick-slip nodes contributes to a higher energy dissipation through the contact friction in the dovetail. This is why the peak of forced response decreases in spite of the increase in the excitation level. Figure 8 compares the distribution of the stick-slip node (in red) and purely sticking node (in blue) on the contact surface in the resonance frequency between different excitation levels. The sliding contact node starts to appear from the edge of the inner surface and the middle of the outboard surface. This is because the relative displacement in these areas is relatively large when the blade is under the first flapping motion. As the excitation level increases, the region with the stick-slip node spreads towards to the middle of the inner surface and the edge of the outboard surface. The more these stickslip nodes, the higher the energy dissipation. This trend of the increased number of stick-slip nodes is consistent with the reduced resonance amplitude in Fig. 6. Figure 9 compares the size of the adaptive reducedorder model at each frequency point for these four excitation levels. One can clearly observe that the size of reduced system is updated consistently according to the number of the sliding contact nodes as shown in Fig.  7. At a low excitation level, for example 1 N, all the contact nodes are in a sticking condition. The size of system remains at a minimum size of 700 for all fre- Fig. 8 The distribution of stick-slip nodes (red) and stick nodes (blue) in the projected contact surface at the resonance a F = 8 N b Fig. 9 The size of adaptive reduced-order model at different excitation levels Fig. 10 Computing speed-up factor of adaptive ROM against the CMS-based Rubin method quency points, which is 100 × (2N H +1). 100 is the total number of vibrational modes from the blade and disc. N H is the number of harmonics used in harmonic balance method. Consistent with contact condition, as the excitation level increases, the size of adaptive reduced system increases and expands toward to the further sides of resonance frequency point. Figure 10 shows the computational speedup between the original adaptive ROM-and CMS-based Rubin method at different excitation levels. The computing speedup varies from 36 to 120 as the excitation level decreases from 8 to 1 N. The advantages are somewhat reduced as full slip is reached, since in that case all nonlinear modes need to maintain in the model. This demonstrates that the original adaptive ROM method can achieve impressive computing speedups for such large-scale industrial models as it is similar to the small scale systems studied in [36]. Now that it can be confirmed, that the previous original method can be used for large-scale systems, the effectiveness of the improved adaptive ROM on the further improvement of computational efficiency will be then evaluated. Figure 11 shows the comparison of the FRFs between the improved adaptive ROM-and CMS-based Rubin method at the same four excitation levels. The improved adaptive ROM with error estimator can more accurately capture the forced responses obtained from the Rubin method at all the excitation levels. In comparison with Fig. 6, one can clearly observe the improved adaptive method reduces most of the difference between the original adaptive method and Rubin method. This improvement is more obvious at the excitation levels when the sliding takes place on the joint interface. The proposed error estimator can help to reduce 2.2% difference at the excitation level of 8 N. This suggests that the improved adaptive method can not only retain the accuracy of the forced response but also help reduce to the errors from the auto updating algorithm. More importantly, it also indicates that the numerical errors due to the updating algorithm may actually contribute to the tiny amounts of "dissipating energy" captured in the contact interface that should have been removed. Otherwise, the accuracy of the forced response would be compromised. Table 1 compares the computational time between improved and original adaptive methods at the different excitation levels. The computing speedup of the improved adaptive method against the original method is highlighted in the Table. It shows the improved adaptive method can speed up the simulations by 2.38 and 2.07 times at 8 and 5 N, respectively, when compared to the original method that has already achieved the impressive speed-ups shown in Fig. 10. As expected, the speedup of the improved adaptive method gradually decreases to one as the excitation level is reduced. It is because the number of sliding nodes reduces to zero if the whole contact interface is in a sticking condition, and hence the novel approach cannot remove any more  elements that have not already been removed from the original approach. Figure 12 compares the retained slipping nodes between the original adaptive method and the improved one with error estimator at the excitation level of 8 N. The contact nodes in a stick-slip condition are plotted in a colored circle. Different colors represent the different levels of energy dissipation as it is shown in the colorbar. One can clearly observe that the number of retained stick-slip nodes is reduced significantly using the improved adaptive method. It is worth noting here that a very conservative value of ε, 10 −10 (J), is used in error estimator. It means a sliding contact node needs to reach the dissipated energy of 10 −10 (J) before it is retained in the reduced set. This suggests that there is a large number of retained sliding nodes in the original method having actually very tiny contributions to the energy dissipation. This tiny energy dissipation may also be due to the errors in the automatic size update and/or numerical errors from the iterative solver. As shown in Fig. 12b, removing these sliding nodes from the reduced set does not actually affect the capture of  the sliding nodes with significant energy dissipation during the automatic size updating. Figure 13 compares the cumulated total energy dissipation of all the contact nodes between these two methods at the frequency points where the sliding occurs. The improved adaptive ROM can retain the almost same dissipated energy level as the original adaptive method for each frequency point. Figure 14 compares the distribution of included stick-slip nodes between the original and improved adaptive methods in the resonance frequency. The dissipated energy of each retained stick-slip contact node is also displayed. As it is similar in Fig. 8, due to the twist in the fan blade and asymmetry of the contact surface, the distribution of the stick-slip is rather asymmetric. With the help of the error estimator, 54 stick-slip nodes with extremely low energy dissipation levels are successfully removed from the original adaptive method that accounts about 40% of total number of contact nodes. As expected, the sliding nodes with the energy dissipation over than ε are retained in the improved reduced-order model. Figure 15 compares the size of the reduced model between the improved and original adaptive ROM, respectively. One can clearly see that, compared to improved adaptive ROM, the size of original method is reduced tremendously at the frequency close to the resonance. When the whole contact surface is in a sticking condition, the size of reduced model for both methods are kept same. In this case, the improved adaptive method does not improve the computational efficiency further. This is why, at the excitation level of 1N, the improved adaptive method does not offer a further speedup against the original method.

Improved adaptive ROM method
Overall, it confirms that the proposed error estimator indeed plays a very effective role in removing the low energy stick-slip nodes from the original adaptive ROM contributing to a further computing speedup in compared with the original adaptive method without the compromise of the accuracy in vibration predictions.

Conclusions
This paper described an improved adaptive ROM technique to enhance the computational efficiency and accuracy of nonlinear vibrations of the large-scale assembly. An energy-based adaptive error estimator was proposed and integrated into the original adaptive formulation for an improved automatic size update. It was effectively used as a monitoring indicator to update the set of reduced basis through adding or removing associated static modes during the on-line computation. The integration of this error estimator into the automatic size update algorithm within the HBM framework was also presented.
The proposed adaptive ROM method was studied with an industrial scale problem in jet engines. A fan blade system with curved dovetail joints was used as a test case. The study confirms that the adaptive ROM method can work well for the large-scale problem. It achieves an extremely high speed up when compared to the classical CMS technique. In comparison with the Rubin method, the original adaptive ROM can reduce the computational time by up to 100 times. Even at a high excitation level when most of the contact is in a sliding condition, the original adaptive method is still 50 times fast.
It also proves that the improved adaptive method is able to further reduce the computing time when compared to the original method thanks to the proposed error estimator. The error estimator is able to help remove those static modes associated with the stick-slip nodes with extremely low energy dissipation. It leads to a considerable computational boost by up to 2.3 times especially at a high excitation level. It was found, at an excitation level of 8N, the number of retained sliding nodes is reduced by 40% when compared to the original adaptive method. In terms of the computational accuracy, the FRFs obtained from the improved adaptive ROM are actually more accurate than that from the original method. The case study showed up to 2.2% difference between the original method and Rubin method is completely eliminated on the use of improved adaptive method. It indicates that the tiny amounts of so-called "dissipating energy" captured on the contact interface may also be due to the automatic updating errors or the numerical errors from the iterative solver. This evidences that the proposed error estimator can also improve the computational accuracy significantly as the error estimator can exclude the update error and/or numerical errors from the original automatic size update algorithm.
In conclusion, the study confirms that the proposed improved adaptive method is an efficient and accurate ROM technique for nonlinear vibration analysis of large-scale assembly with contact friction interfaces, particularly when the joint is in a micro-slip motion.