# A Reduced Order Model for Joint Assemblies by Hyper-Reduction and Model-
Driven Sampling

Ahmed Amr Morsy111Corresponding author.
E-mail address<EMAIL_ADDRESS>,† Institute for Mechanical Systems, ETH
Zürich, Leonhardstrasse 21, 8092 Zurich, Switzerland  Mariella Kast222Both
authors contributed equally to this work. Institute for Mechanical Systems,
ETH Zürich, Leonhardstrasse 21, 8092 Zurich, Switzerland  Paolo Tiso
Institute for Mechanical Systems, ETH Zürich, Leonhardstrasse 21, 8092 Zurich,
Switzerland

Abstract: The dynamic behavior of jointed assemblies exhibiting friction
nonlinearities features amplitude-dependent dissipation and stiffness. To
develop numerical simulations for predictive and design purposes, macro-scale
High Fidelity Models (HFMs) of the contact interfaces are required. However,
the high computational cost of such HFMs impedes the feasibility of the
simulations. To this end, we propose a model-driven method for constructing
hyper-reduced order models of such assemblies. Focusing on steady-state
analysis, we use the Multi-Harmonic Balance Method (MHBM) to formulate the
equations of motion in frequency domain. The reduction basis is constructed
through solving a set of vibration problems corresponding to fictitious
interface conditions. Subsequently, a Galerkin projection reduces the order of
the model. Nonetheless, the necessary fine discretization of the interfaces
represents a bottleneck for achieving high speedups. For this reason, we
implement an adapted Energy Conserving Weighing and Sampling (ECSW) technique
for Hyper Reduction (HR), thereby allowing significant speedups for meshes of
arbitrary fineness. This feature is particularly advantageous since analysts
typically encounter a trade-off between accuracy and computational cost when
deciding on the mesh size, whose estimation is particularly challenging for
problems of this type. To assess the accuracy of our method without resorting
to the HF solution, we propose an error indicator with thresholds that have
proven reliable in our analyses. Finally, the accuracy and efficiency of the
method are demonstrated by two case studies.
Keywords: jointed structures, friction, model order reduction, harmonic
balance method, hyper-reduction

## 1\. Introduction

Friction occurs on contact interfaces of joints, which are commonly found in
mechanical, civil and aerospace engineering applications. The dynamic response
of jointed assemblies can be heavily influenced by the nonlinearities of
frictional contact resulting in amplitude-dependent dissipation and stiffness
[1]. As demands increase on optimization for performance and cost, including
minimizing experimental testing for calibration of models during design
phases, predictive numerical simulations of joints become increasingly
important. Any attempt to model friction inevitably depends on the length and
time scales considered, as well as the purpose of the analysis. For instance,
regarding the time scale, Fantetti et al. [2] have recently presented
insightful experimental evidence demonstrating the time evolution of both the
coefficient of friction and the tangential stiffness. A significant increase
in both parameter values during an initial number of fretting cycles, followed
by a steady-state stabilization of the parameters has been observed. Yuan et
al. [3] have built upon the experimental results reported in [2], using the
statistical distribution of the stabilized values from the different tests, to
perform a sensitivity analysis with respect to the contact parameters for
steady-state experimental and numerical analysis of an underplatform damper
system. The variability of the friction coefficient was shown to significantly
affect the amplitude of the responses which featured considerable slipping
along the interfaces, while the impact of the variability of the tangential
stiffness was predominate for responses that featured low amount of slipping.
The time scale we consider for our analysis corresponds to the steady-state
behavior of the system, which is assumed to emerge after an initial number of
fretting cycles.

As for the length scale, our work focuses on the accurate prediction of
friction-induced dissipation, which demands in turn the capacity to simulate
the microslip behavior in joints [4]. Microslip refers to the isolated,
partial slipping of the interfaces relative to each other, in contrast to
gross slip, which is the complete slipping of the interfaces. Different
modelling strategies exist to reproduce/predict the microslip behavior. A
recent review of damping models for friction has been presented in [5], where
connections between different damping models were highlighted. For predictive
analyses, Macro-scale High-Fidelity Models (HFMs) which involve dense meshes
of nonlinear elements, are widely used. In this context, modelling of
frictional contact differs depending on assumptions as to whether, the
tangential motion on the interface is one-dimensional or two-dimensional,
coupling is assumed for tangential planar motion, and if contact pressure
temporally evolves across the interface. Firrone and Zucca compared features
of the responses obtained by employing different formulations in [6]. Brake et
al. [7] showed experimentally that a correct simulation of joint behavior must
take into account time varying contact pressure across the interface. An
interfacial mesh whose elements can reproduce this effect, equipped with an
elastic or rigid Coulomb law, makes the depiction of global microslip behavior
possible [8]. Yet, despite the accuracy of such HFMs, a major drawback is
their high computational cost, which stands in the way of efficient,
predictive analyses.

Projection-based Reduced Order Models (ROMs) reduce the size of dynamical
systems by projecting them onto suitable low-dimensional subspaces, thus
providing accurate and efficient solutions. An efficient family of methods
used extensively in substructuring is Component Mode Synthesis (CMS) methods
[9], which include the widely common Hurty/Craig-Bampton (HCB) method [10]. In
fact, many ROM techniques aimed for frictional contact problems implement a
CMS first step, by which the degrees of freedom (DOFs) are reduced into fixed-
interface modes and physically-retained boundary DOFs, which include
interfacial DOFs. A Galerkin projection of the system onto the subspace
defined by the transformation matrix is then carried out. HCB provides an
example of a model-driven ROM, as opposed to a data-driven ROM, since the
construction of the reduction basis is based solely on the model and no data
from HF simulations is required. Due to the high number of retained
interfacial DOFs in the case of HFM of joints, further reductions are
necessary. Therefore, attempts have been made to supplement CMS methods using
the concept of modal derivatives[11, 12]. The concept, which has been
implemented for problems featuring geometric nonlinearity [13], has been
adapted for friction problems by Pichler et al. [14]. In their work, they
formulate Trial Vector Derivatives (TVDs) for time domain analysis and use
them for test cases where the normal contact, in particular, is accorded
predominant importance. A comparison of the efficiency of different contact
and friction models has also been presented in [15]. Hughes et al. [16]
proposed an interface reduction method where they use System-level
Characteristic Constraint (SCC) modes and Properly Orthogonal Interface Modal
Derivatives (POIMDs). Their results show that global quantities can be
reliably assessed using SCC modes with impressive speedups reaching to the
order of thousands, while a high accuracy on local quantities such as slipping
areas is possible with the addition of POIMDs, with overall speedups ranging
between 10-20.

In parallel, ROM techniques for frequency domain analysis focusing on steady-
state analysis have been developing. Krack presented the concept of extended
periodic motion [17], which is essential to his extension of nonlinear modes
to friction-damped systems [18]. Using this method, the dynamics of the system
are assumed to take place on a manifold that is constituted of a union of
periodic motions coming from the computed nonlinear mode. The external forcing
and linear damping are subsequently superimposed to the ROM. The method proved
accurate and efficient for the test cases in [18], where reasonable
assumptions regarding linear damping and isolated resonance are met. Mitra et
al. [19] proposed the Amplitude Microslip Projection (AMP) method, where the
construction of the reduction basis is based on eigensolutions at linearized
configurations which are obtained by artificially probing the contact
interface. Gastaldi et al. [20] presented the Jacobian Projection (JP)
technique, which extends AMP to the multi-harmonic context, accounting for
static-dynamic coupling. Along this line, we present an Augmented JP method,
whose augmented vectors have proven necessary for accurate detection of local
quantities such as stresses on the interface, as well as global amplitude-
dependent dissipation and frequency particularly in a case where multiple
contact interfaces are present. Yet, in spite of the reduced size of a ROM, a
typical problem in case of HF modelling of joints is the large number of
nonlinear elements, which in turn demands expensive computations in full-space
to compute the nonlinear forces. Therefore, attention is directed to Hyper
Reduction (HR) techniques that target reducing the computational cost of the
nonlinear forces.

Balaji et al. [21] have presented a ROM strategy through HR of the interfaces.
In their work, HR refers to a re-presentation of the contact interface such
that the computation of nonlinear forces is based solely on the reduced DOFs.
They present two strategies: a novel whole-joint formulation that is
representative of the HF model, where the interface is divided into patches
based on level sets of a field objective such as the contact pressure after a
static step, and another strategy by which coarser re-meshing of the interface
with contours of a chosen field objective are used to place new nodes on the
interface. Their results captured modal characteristics using a system of
significantly reduced size. Still in the spirit of reducing nonlinear
computations, Yuan et al. [22] present an adapted component mode synthesis
method where the idea is to reduce the interface DOFs by adapting the
reduction basis based on the nodes that are expected to actually exhibit
nonlinear behavior, i.e. not sticking. Their novel approach is conceptually in
line with continuation methods from the standpoint of optimizing the choice of
initial conditions as the selection of (nonlinear) reduced interface DOFs is
updated at each new frequency point. Speedups ranging from 50-100 have been
reported considering the cost per iteration for a test case with relatively
restrictive assumptions regarding contact conditions. A further improvement
was presented in [23] that excludes numerical noise, as well as nodes that are
dominated by sticking behavior, from influencing the choice of reduced
interface DOFs. To contrast our work with existing methods, we regard that
relative to [21], we explore employing a HR technique that retains the
original interfacial mesh, and in relation to [22, 23], we test considering
only one reduced set of elements for the nonlinear computations at all
frequency points.

(Discrete) Empirical Interpolation Method (D)EIM [24] and Energy Conserving
Sampling and Weighting (ECSW) [25] are popular techniques to approximate
nonlinear terms which require long computation times. In particular, the ECSW
method approximates the nonlinear forces via energy-based quantities, where
the respective weights are computed from HF snapshots of the forces [25]. If
the ROM is data-driven, the HF force snapshots required for the generation of
training vectors will present no additional costs. On the other hand, for a
model-driven ROM, the idea is to avoid costly HF solutions. Therefore, for
instance, in the context of geometrically nonlinear problems, it was proposed
to perform a lifting of the linear modal subspace onto a quadratic manifold,
thereby eliminating the need for HF snapshots [26]. To cover a wider class of
geometric nonlinearities, the training generation proposed later in [27] comes
from solutions to nonlinear static systems where the applied forcing belongs
to Krylov subspace and is stochastically produced. In our work, we present an
adaptation of ECSW aimed at friction nonlinearities, where the training
generation presents no additional costs with respect to the ROM construction.
Moreover, it is well-known that the determination of the mesh size for
frictional contact problems is a daunting task and convergence studies with
respect to interfacial mesh size are rarely carried out [8]. To this end, we
will show that such studies become completely within reach using our method.

The paper is organized as follows: in Section 2.1, we present the mathematical
formulation of the problem, along with the adopted contact law for frictional
contact, and the computation of contact displacement is illustrated. An
overview of frequency-domain analysis using MHBM is presented in Section 2.2.
Next, we express two typical preliminary computations in Section 3.1.
Subsequently, we describe the Augmented JP ROM strategy in 3.2, then review
ECSW HR in Section 3.3 and adapt it to frictional contact in Section 3.4.
Figures 2 and 3 depict a complete flowchart of the method. In Section 3.5, we
present an error indicator that depends only on the results of HR-ROM. The
accuracy and efficiency of the method is tested on two case studies in
sections 4 and 5. Finally, we present the conclusions in Section 6.

## 2\. Joint modelling framework for vibrations

### 2.1. Problem Formulation

The momentum balance of a structural continuum consisting of an assembly of
jointed structural components, discretized using Finite Element Method (FEM),
can be expressed as

$\boldsymbol{M}\boldsymbol{\mathrm{\ddot{u}}}\text{(t)}+\boldsymbol{C}\boldsymbol{\mathrm{\dot{u}}}\text{(t)}+\boldsymbol{K}\boldsymbol{\mathrm{u}}\text{(t)}+\boldsymbol{\mathrm{f}}(\boldsymbol{\mathrm{u}}\text{(t)},\boldsymbol{\mathrm{\dot{u}}}\text{(t)})=\boldsymbol{\mathrm{p}}_{\text{ext,}i}\text{(t)}$
(1)

with

$\boldsymbol{\mathrm{p}}_{\text{ext,}i}\text{(t)}=\boldsymbol{\mathrm{p}}_{\text{0}}+\alpha_{i}\boldsymbol{\mathrm{p}}_{\text{E}}^{\text{}}(\Omega\text{t})$
(2)

where $\boldsymbol{\mathrm{u}}\text{(t)}\in\mathbb{R}^{n}$ is the vector of
nodal displacements, with $n$ being the number of Degrees of Freedom (DOFs).
The time dependence of the kinematic variables in (1) will be dropped from now
on for clarity. We assume linear elasticity and small displacements. Matrices
$\boldsymbol{K},\boldsymbol{C},\boldsymbol{M}\in\mathbb{R}^{n\times n}$ are
the ”free” stiffness, damping, and mass matrix of the mechanical system, which
means that no coupling effects are considered in the assembly of the matrices
due to the different components. The vector of external forces
$\boldsymbol{\mathrm{p}}_{\text{ext,i}}\in\mathbb{R}^{n}$ is split into a
vector of static forces $\boldsymbol{\mathrm{p}}_{\text{0}}$, including
possibly bolt preclamping forces, and a periodic forcing of unit value in
$\boldsymbol{\mathrm{p}}_{\text{E}}\ (\Omega\text{t})$, with an excitation
frequency $\Omega$, amplified with a set of $a$ amplitudes of excitation
$\alpha_{i}$ where $i=1,...,a$.

The vector
$\boldsymbol{\mathrm{f}}(\boldsymbol{\mathrm{u}},\boldsymbol{\mathrm{\dot{u}}})\in\mathbb{R}^{n}$
represents the modelled frictional contact forces. Since we are studying
joints in the microslip regime, we carry over the assumption of small
displacements to contact kinematics, therefore contact is assumed to occur
between fixed nodal pairs on the interface during vibrations. Accordingly, we
adopt the node-to-node Jenkins element formulation with a unilateral spring in
the normal direction [28]. The normal contact force
$\text{f}_{\text{n}}^{\text{}}$ is given by

$\text{f}_{\text{n}}^{\text{}}=\begin{cases}\text{k}_{\text{n}}^{\text{}}\Delta\text{u}_{\text{n}}^{\text{}}&\text{for
penetration},\\\ 0&\text{for separation},\end{cases}$ (3)

where $\Delta\text{u}_{\text{n}}^{\text{}}$ is the relative normal
displacement of the contact nodal pair, $\text{k}_{\text{n}}^{\text{}}$ is
normal contact stiffness. As for the tangential force,

$\text{f}_{\text{t}}^{\text{}}=\begin{cases}\text{k}_{\text{t}}^{\text{}}(\Delta\text{u}_{\text{t}}^{\text{}}-\Delta\text{u}_{\text{t}}^{\text{st}})+\text{f}_{\text{t}}^{\text{st}}&\text{for
sticking,}\\\ \mu\text{f}_{\text{n}}^{\text{}}\zeta&\text{for slipping},\\\
0&\text{for separation},\end{cases}$ (4)

where $\text{f}_{\text{t}}^{\text{}}$ is the tangential force,
$\text{k}_{\text{t}}^{\text{}}$ is the tangent contact stiffness,
$\Delta\text{u}_{\text{t}}^{\text{st}}$ and $\text{f}_{\text{t}}^{\text{st}}$
refer to the relative displacement and the friction force between two
interface nodes at the start of the current stick phase. Regarding the
slipping case, the Coulomb friction coefficient is denoted by $\mu$, and
$\zeta=\text{sgn}(\text{f}_{\text{t}}^{\text{slip}})$ represents the sign of
$\text{f}_{\text{t}}^{\text{slip}}$, which is the tangential force at the end
of the previous stick phase [28]. The transition between sticking and slipping
is thus defined by a slip limit force, which depends on the normal contact
force $|\mu\text{f}_{\text{n}}^{\text{}}|$. Note that despite the non-smooth
evolution of contact forces in time due to the possibility of separation and
slipping behaviors, the contact forces remain both piece-wise linear and
$C^{0}$ continuous in time. Figure 1 illustrates the computation of relative
displacements, and represents the contact law in case of contact.

Physically, this model allows for three contact behaviors on a local level:
sticking, slipping and lift-off. For small relative tangential displacements
between the node pairs, the friction force is elastic, whereas slippage
between the nodes, such as depicted by the original Coulomb type friction, is
assumed to occur once a threshold of displacement is reached depending on the
normal contact force. The modelled friction force satisfies the Masing
hypothesis leading to a closed hysteresis curve for periodic motions[29, 30].
The full contact force vector
$\boldsymbol{\mathrm{f}}(\boldsymbol{\mathrm{u}})$ is constructed by taking
into account the behavior of all Jenkins elements across the contact
interfaces, which leads to a nonlinear, non-smooth, spatial variation of
forces.

Figure 1: Tangential and normal forces in case of contact

### 2.2. Overview of Multi-Harmonic Balance Method

The dynamic response of mechanical systems can be analyzed in time or
frequency domains. For time domain analysis, time integration schemes are
employed to transform the system of ODEs in (1) into a set of nonlinear
algebraic equations to be solved at each time step. If one is interested in
investigating a periodic response under steady-state conditions due to one
excitation frequency, a large number of time steps is required until the
transient response decays, especially when structural damping is low.
Alternatively, one can efficiently study the effect of multiple excitation
amplitudes across a wide range of frequencies by resorting to a frequency
domain approach. In particular, we implement MHBM, which we review in this
section.

To compute a periodic solution of (1), the steady-state displacement solution
$\boldsymbol{\mathrm{u}}\text{(t)}$ is approximated using Fourier base
functions as

$\boldsymbol{\mathrm{u}}_{\text{h}}(\text{t})=\boldsymbol{\mathrm{U}}_{0}+\sum_{\text{j}=1}^{\text{H}}\left({\boldsymbol{\mathrm{U}}_{\text{j}}^{\text{c}}\cos(\text{j}\Omega\text{t})+\boldsymbol{\mathrm{U}}_{\text{j}}^{\text{s}}\sin(\text{j}\Omega\text{t})}\right),$
(5)

where the vector $\boldsymbol{\mathrm{U}}_{0}\in\mathbb{R}^{n}$ contains the
static coefficients, H is the chosen number of harmonics, and the vectors
$\boldsymbol{\mathrm{U}}_{\text{c}}^{\text{j}},\boldsymbol{\mathrm{U}}_{\text{s}}^{\text{j}}\in\mathbb{R}^{n}$
contain the Fourier coefficients of the cosine and sine components of the j-th
harmonic, respectively. These components are collected in a column vector
$\boldsymbol{\mathrm{U}}$ as
$\boldsymbol{\mathrm{U}}=[\boldsymbol{\mathrm{U}}_{\text{0}}^{\text{}}\
\boldsymbol{\mathrm{U}}_{\text{1}}^{\text{c}}\
\boldsymbol{\mathrm{U}}_{\text{1}}^{\text{s}}\ ...\
\boldsymbol{\mathrm{U}}_{\text{H}}^{\text{c}}\
\boldsymbol{\mathrm{U}}_{\text{H}}^{\text{s}}]^{\text{T}}$. Similarly, the
nonlinear forces are approximated as

$\boldsymbol{\mathrm{f}}_{\text{h}}^{\text{}}(\text{t})=\boldsymbol{\mathrm{F}}_{\text{0}}^{\text{}}+\sum_{\text{j}=1}^{\text{H}}\left({\boldsymbol{\mathrm{F}}_{\text{j}}^{\text{c}}\cos({\text{j}\Omega\text{t}})+\boldsymbol{\mathrm{F}}_{\text{j}}^{\text{s}}\sin{(\text{j}\Omega\text{t}})}\right),$
(6)

where

$\begin{split}\boldsymbol{\mathrm{F}}_{\text{0}}^{\text{}}&=\dfrac{1}{\text{T}}\int_{0}^{\text{T}}\boldsymbol{\mathrm{f}}\boldsymbol{\mathrm{(u)}}\mathrm{d}\text{t},\\\
\boldsymbol{\mathrm{F}}_{\text{j}}^{\text{c}}&=\dfrac{2}{\text{T}}\int_{0}^{\text{T}}\boldsymbol{\mathrm{f}}\boldsymbol{\mathrm{(u)}}\cos({\text{j}\Omega\text{t}})\mathrm{d}\text{t},\\\
\boldsymbol{\mathrm{F}}_{\text{j}}^{\text{s}}&=\dfrac{2}{\text{T}}\int_{0}^{\text{T}}\boldsymbol{\mathrm{f}}\boldsymbol{\mathrm{(u)}}\sin({\text{j}\Omega\text{t}})\mathrm{d}\text{t}.\end{split}$
(7)

However, the closed-form evaluation of the integrals in (7) is complicated due
to the coupling of the different harmonics, and the diverse phase transitions
that an element potentially undergoes during one time period. Therefore, we
employ the Alternating-Frequency Time (AFT) scheme [31, 32], which can be
summarized as follows: given the displacement vector $\boldsymbol{\mathrm{U}}$
in the frequency-domain, one reconstructs
$\boldsymbol{\mathrm{u}}_{\text{h}}^{\text{}}(\text{t})$ across a sampled time
period using (5). Subsequently, the nonlinear forces are evaluated at the
sampled time instants using the adopted contact law. One can then use a FFT to
obtain the Fourier coefficients of the contact force vector in (7).

MHBM refers to a weighted residual formulation of the equations of motion that
are now represented using the Ansatz, where the residual forces are imposed to
be orthogonal to the weight functions (Fourier base functions) in a weak sense
(see [33], for instance, for more details). Consequently, one obtains

$\boldsymbol{Z}\boldsymbol{\mathrm{U}}+\boldsymbol{\mathrm{F}}(\boldsymbol{\mathrm{U}})-\boldsymbol{\mathrm{P}}_{\text{ext,}i}=\boldsymbol{\mathrm{0}},$
(8)

where
$\boldsymbol{\mathrm{F}},\boldsymbol{\mathrm{P}}_{\text{ext,}i}\in\mathbb{R}^{n(2H+1)}$
are vectors collecting all Fourier coefficients of the nonlinear forces and
the external forces as

$\boldsymbol{\mathrm{F}}=[\boldsymbol{\mathrm{F}}_{\text{0}}^{\text{}}\
\boldsymbol{\mathrm{F}}_{\text{1}}^{\text{c}}\
\boldsymbol{\mathrm{F}}_{\text{1}}^{\text{s}}\ ...\
\boldsymbol{\mathrm{F}}_{\text{H}}^{\text{c}}\
\boldsymbol{\mathrm{F}}_{\text{H}}^{\text{s}}]^{\text{T}},\\\ $ (9)
$\boldsymbol{\mathrm{P}}_{\text{ext,}i}=\boldsymbol{\mathrm{P}}_{\text{0}}^{\text{}}\
+\alpha_{i}\boldsymbol{\mathrm{P}}_{\text{E}}^{\text{}}$ (10)

where

$\begin{split}\boldsymbol{\mathrm{P}}_{\text{0}}^{\text{}}=[\boldsymbol{\mathrm{P}}_{\text{0}}^{\text{}}\
\boldsymbol{\mathrm{0}}_{\text{}}^{\text{}}\
\boldsymbol{\mathrm{0}}_{\text{}}^{\text{}}\ ...\
\boldsymbol{\mathrm{0}}_{\text{}}^{\text{}}\
\boldsymbol{\mathrm{0}}_{\text{}}^{\text{}}]^{\text{T}},\\\
\boldsymbol{\mathrm{P}}_{\text{E}}^{\text{}}=[\boldsymbol{\mathrm{0}}_{\text{}}^{\text{}}\
\boldsymbol{\mathrm{P}}_{\text{E,1}}^{\text{c}}\
\boldsymbol{\mathrm{P}}_{\text{E,1}}^{\text{s}}\
\boldsymbol{\mathrm{0}}_{\text{}}^{\text{}}\ ...\
\boldsymbol{\mathrm{0}}_{\text{}}^{\text{}}\
\boldsymbol{\mathrm{0}}_{\text{}}^{\text{}}]^{\text{T}},\\\ \end{split}$ (11)

and $\boldsymbol{Z}$ is a blockdiagonal matrix, with block matrices
$\boldsymbol{Z}_{\text{j}}^{\text{}}\in\mathbb{R}^{n\times n}$:

$\boldsymbol{Z}=\begin{bmatrix}\boldsymbol{Z}_{\text{0}}^{\text{}}&&&&\\\
&\boldsymbol{Z}_{\text{1}}^{\text{}}&&&\\\ &&...&&&\\\
&&&\boldsymbol{Z}_{\text{j}}^{\text{}}&&\\\ &&&&..&\\\
&&&&&\boldsymbol{Z}_{\text{h}}^{\text{}}\end{bmatrix},$ (12)

where $\boldsymbol{Z}_{\text{0}}^{\text{}}=\boldsymbol{K}$, and

$\boldsymbol{Z}_{\text{j}}^{\text{}}=\boldsymbol{Z}_{\text{MC,j}}^{\text{}}(\Omega)+\boldsymbol{Z}_{\text{K}}^{\text{}}=\begin{bmatrix}-(\text{j}\Omega)^{2}\boldsymbol{M}&\text{j}\Omega\boldsymbol{C}\\\
-\text{j}\Omega\boldsymbol{C}&-(\text{j}\Omega)^{2}\boldsymbol{M}\end{bmatrix}+\begin{bmatrix}\boldsymbol{K}&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{K}\end{bmatrix}$ (13)

We refer to the solution of (8) as the High Fidelity (HF) solution. Note that
due to the nonlinearity of (8), the solution is sought using an iterative
Newton-Raphson (N-R) scheme:

$\begin{split}\boldsymbol{\mathrm{R}}(\boldsymbol{\mathrm{U}}^{(r)})=\boldsymbol{Z}\boldsymbol{\mathrm{U}}^{(r)}&+\boldsymbol{\mathrm{F}}(\boldsymbol{\mathrm{U}}^{(r)})-\boldsymbol{\mathrm{P}}_{\text{ext,}i}\\\
\
\boldsymbol{\mathrm{R}}(\boldsymbol{\mathrm{U}}^{(r+1)})\approx\boldsymbol{\mathrm{R}}(\boldsymbol{\mathrm{U}}^{(r)})&+\left.\dfrac{\partial\boldsymbol{\mathrm{R}}(\boldsymbol{\mathrm{U}})}{\partial\boldsymbol{\mathrm{U}}}\right|_{\boldsymbol{\mathrm{U}}^{(r)}}\left({\boldsymbol{\mathrm{U}}^{(r+1)}-\boldsymbol{\mathrm{U}}^{(r)}}\right)=\boldsymbol{\mathrm{0}},\\\
\end{split}$ (14)

where superscript (r) denotes the r-th iteration, and

$\dfrac{\partial\boldsymbol{\mathrm{R}}(\boldsymbol{\mathrm{U}})}{\partial\boldsymbol{\mathrm{U}}}=\boldsymbol{Z}+\dfrac{\partial\boldsymbol{\mathrm{F}}(\boldsymbol{\mathrm{U}})}{\partial\boldsymbol{\mathrm{U}}},$

where the evaluation of the matrix
$\dfrac{\partial\boldsymbol{\mathrm{F}}(\boldsymbol{\mathrm{U}})}{\partial\boldsymbol{\mathrm{U}}}$
is carried out using the AFT procedure, with the derivatives analytically
computed in the time domain, then transformed to the frequency domain using
the FFT [32, 34]. For the first iteration of (14), we use as an initial guess
the converged solution of the previous frequency point in the frequency sweep.
As for the very first iteration, we use a linear dynamic solution as an
initial guess, as we discuss in the next section.

## 3\. Hyper Reduced ROM for Joints

### 3.1. Preliminary computations

  * •

_Modal analysis of preclamped configuration_
In order to determine the approximate frequency of the mode shape of interest
of the jointed structure, we first solve the nonlinear static step:

$\boldsymbol{K}\boldsymbol{\mathrm{u}}_{\text{0}}^{\text{}}+\boldsymbol{\mathrm{f}}(\boldsymbol{\mathrm{u}}_{\text{0}}^{\text{}})=\boldsymbol{\mathrm{p}}_{\text{0}}^{\text{}},$
(15)

where $\boldsymbol{\mathrm{u_{\text{0}}}}\in\mathbb{R}^{n}$ is the solution to
the displacements under the effect of preclamping forces
$\boldsymbol{\mathrm{p}}_{\text{0}}^{\text{}}$. Having computed the static
configuration $\boldsymbol{\mathrm{u_{\text{0}}}}$, we define a stiffness
matrix $\boldsymbol{K}_{\text{\emph{l}}}^{\text{}}$ that is linearized around
$\boldsymbol{\mathrm{u}}_{\text{0}}^{\text{}}$, which now takes into account
the contact states at the interface, as

$\boldsymbol{K}_{\text{\emph{l}}}^{\text{}}=\boldsymbol{K}+\dfrac{\partial\boldsymbol{\mathrm{f}}(\boldsymbol{\mathrm{u}}_{\text{}}^{\text{}})}{\partial\boldsymbol{\mathrm{u}}}\bigg{\rvert}_{\text{u}_{\text{0}}^{\text{}}}$
(16)

Next, we perform a modal analysis of the linearized, preclamped structure by
the solving the eigenvalue problem of size $n$:

$(\boldsymbol{K}_{\text{\emph{l}}}^{\text{}}-\omega^{2}_{\text{i}}\boldsymbol{M})\boldsymbol{\mathrm{\phi_{\text{i}}}}=\boldsymbol{\mathrm{0}}.$
(17)

From the set of eigenfrequencies $\omega_{\text{i}}$, we choose that
corresponding to the vibration mode of interest, and denote it by
$\overline{\omega}$.

  * •

_Linearized dynamic solution_
We compute the dynamic part of the response of the linearized structure using

$\boldsymbol{Z}_{\text{\emph{l}}}^{\text{}}\boldsymbol{\mathrm{U}}_{\text{\emph{l}}}^{\text{}}=\boldsymbol{\mathrm{P}}_{\text{E}}^{\text{1-H}},\\\
$ (18)

where 1-H refers to retaining only the dynamic components of the vector
$\boldsymbol{\mathrm{P}}_{\text{E}}^{\text{}}$ such that
$\boldsymbol{\mathrm{P}}_{\text{E}}^{\text{1-H}}\in\mathbb{R}^{n(2H)}$. Vector
$\boldsymbol{\mathrm{U}}_{\text{\emph{l}}}^{\text{}}\in\mathbb{R}^{n(2H)}$ is
the displacement due to the unitary force vector
$\boldsymbol{\mathrm{P}}_{\text{E}}^{\text{1-H}}\in\mathbb{R}^{n(2H)}$, while
the matrix
$\boldsymbol{Z}_{\text{\emph{l}}}^{\text{}}\in\mathbb{R}^{n(2H)\times n(2H)}$
now consists of diagonal blocks, where the j-th block written as

$\boldsymbol{Z}_{\text{\emph{l},j}}^{\text{}}=\boldsymbol{Z}_{\text{MC,j}}^{\text{}}(\Omega)+\begin{bmatrix}\boldsymbol{K}_{\text{\emph{l}}}^{\text{}}&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{K}_{\text{\emph{l}}}^{\text{}}\end{bmatrix}.$ (19)

Note that regardless of the number of harmonics H in the MHBM Ansatz, the
dynamic response $\boldsymbol{\mathrm{U}}_{\text{\emph{l}}}^{\text{}}$ only
features one harmonic, due to the block diagonal structure of
$\boldsymbol{Z}_{\text{\emph{l}}}^{\text{}}$ (19). Expressing (18) and (19) is
therefore chosen for convenience and clarity later on. Appending the static
component $\boldsymbol{\mathrm{u}}_{\text{0}}^{\text{}}$ (15), and accounting
for an amplitude of excitation $\alpha_{i}$, the total linear response becomes

$\boldsymbol{\mathrm{\theta}}_{\text{}}^{\text{}}(\alpha_{i},\Omega)=[\boldsymbol{\mathrm{u}}_{\text{0}}^{\text{T}}\
\alpha_{i}\boldsymbol{\mathrm{U}}_{\text{\emph{l}}}^{\text{}}(\Omega)^{\text{T}}].$
(20)

where
$\boldsymbol{\mathrm{\theta}}_{\text{}}^{\text{}}=\boldsymbol{\mathrm{\theta}}_{\text{}}^{\text{}}(\alpha_{i},\Omega_{1})$
is used to compute the initial guess for the first frequency $\Omega_{1}$ for
a frequency sweep corresponding to the i-th amplitude. Moreover, (20) is used
with
$\boldsymbol{\mathrm{\theta}}_{\text{}}^{\text{}}=\boldsymbol{\mathrm{\theta}}_{\text{}}^{\text{}}(\alpha_{i},\overline{\omega})$
in the construction of the reduction basis as discussed in the next section.

### 3.2. Augmented Jacobian Projection Method

Although jointed assemblies exhibit non-smooth dynamics, the adopted contact
law remains piecewise linear. This feature is leveraged in the JP method [20],
since the contribution of each nonlinear element to the jacobian stiffness
matrix becomes straightforward to determine, for given interface condition.
Conceptually, the JP method is based on solutions of free vibration problems
of a set of simulated linear systems, whose eigenvectors, in aggregate, are
assumed to well approximate the response of the nonlinear system in steady-
state condition.

One starts by assuming a trial vector that reflects an expected baseline
configuration of the interface. Next, this trial vector is amplified to
simulate different boundary conditions on the interface, leading to systems
with diverse contact states. Gastaldi et al. [20] suggest using amplifications
of the stuck mode, which refers to a state where the interfaces are fully in
contact, and are in sticking condition333In case the constructed basis results
in a ROM containing high errors, they suggest using a computed displacement
solution as a new trial vector to construct a new basis.. For our purposes, we
construct $m$ linear solutions of (20), corresponding to $m$ amplitudes
$\alpha_{k=1,...m}$, computed at resonance, and denote each by
$\boldsymbol{\mathrm{\theta}}^{k}$, such that

$\boldsymbol{\mathrm{\theta}}^{k}=\boldsymbol{\mathrm{\theta}}_{\text{}}^{\text{}}(\alpha_{k},\overline{\omega}).$
(21)

where $k=1,...m$. Even though our default strategy is to set $m=a$ , where $a$
is the number of excitation amplitudes in (2), an enrichment of the basis with
$m>a$ may be required in some cases as will be shown in Section 5. Next, and
in line with [20], once we choose the amplified displacement vector, it is
then used as an input to the AFT algorithm to compute a nonlinear force vector
$\boldsymbol{\mathrm{F}}^{k}$, and the nonlinear stiffness contribution
$\dfrac{\partial\boldsymbol{\mathrm{F}}^{k}(\boldsymbol{\mathrm{U}})}{\partial\boldsymbol{\mathrm{U}}}$,
for each contact state associated to $\boldsymbol{\mathrm{\theta}}^{k}$. Note
that while the input features dynamics consisting only of one harmonic
component, the outputs will feature higher harmonic contributions due to the
anticipated slipping and separation stimulated by the amplified trial vectors
$\boldsymbol{\mathrm{\theta}}^{k}$. Next, an eigenvalue problem describing the
free vibrations of each intermediary, linear, autonomous system is now posed
as

$\displaystyle(\boldsymbol{J}^{k}-\lambda_{\text{i}}^{k}\
\overline{\boldsymbol{M}})\boldsymbol{\mathrm{\gamma}}_{\text{i}}^{k}=\boldsymbol{\mathrm{0}},$
(22) $\displaystyle\|\boldsymbol{\mathrm{\gamma}}_{\text{i}}^{k}\|=1,$

where

$\boldsymbol{J}^{k}=\boldsymbol{K}{}{}+\dfrac{\partial\boldsymbol{\mathrm{F}}(\boldsymbol{\mathrm{U}}^{k})}{\partial\boldsymbol{\mathrm{U}}},$
(23)

and

$\overline{\boldsymbol{M}}=\begin{bmatrix}\boldsymbol{0}&&&&&\\\
&\boldsymbol{M}_{\text{}}^{\text{}}&&&&&&&&\\\
&&\boldsymbol{M}_{\text{}}^{\text{}}&&&&&&&\\\ &&&\ddots&&&&&&\\\
&&&&\text{j}_{\text{}}^{\text{2}}\boldsymbol{M}_{\text{}}^{\text{}}&&&&&\\\
&&&&&\text{j}_{\text{}}^{\text{2}}\boldsymbol{M}_{\text{}}^{\text{}}&&&&\\\
&&&&&&\ddots&&&\\\ &&&&&&&H^{2}\boldsymbol{M}_{\text{}}^{\text{}}&&\\\
&&&&&&&&H^{2}\boldsymbol{M}_{\text{}}^{\text{}}&\end{bmatrix}=\begin{bmatrix}\boldsymbol{0}&\boldsymbol{0}\\\
\boldsymbol{0}&\overline{\boldsymbol{M}}_{\text{0,H}}\end{bmatrix},$ (24)

with $\boldsymbol{J}^{k},\overline{\boldsymbol{M}}\in\mathbb{R}^{n(2H+1)\times
n(2H+1)}$. To avoid the ill-conditioning of (22) that is associated to the
first zero block diagonal term of $\overline{\boldsymbol{M}}$, a static
condensation step is performed (superscript $k$ is dropped for clarity):

$\begin{bmatrix}\boldsymbol{J}_{\text{0,0}}&\boldsymbol{J}_{\text{0,H}}\\\
\boldsymbol{J}_{\text{H,0}}&\boldsymbol{J}_{\text{H,H}}\end{bmatrix}\begin{bmatrix}\boldsymbol{\mathrm{\gamma}}_{\text{i,0}}^{\text{}}\\\
\boldsymbol{\mathrm{\gamma}}_{\text{i,H}}^{\text{}}\end{bmatrix}-\lambda_{i}\begin{bmatrix}\boldsymbol{0}&\boldsymbol{0}\\\
\boldsymbol{0}&\overline{\boldsymbol{M}}_{\text{H,H}}\end{bmatrix}\begin{bmatrix}\boldsymbol{\mathrm{\gamma}}_{\text{i,0}}^{\text{}}\\\
\boldsymbol{\mathrm{\gamma}}_{\text{i,H}}^{\text{}}\end{bmatrix}=\boldsymbol{\mathrm{0}}$
(25)

$\boldsymbol{J}_{\text{0,0}}\
\boldsymbol{\mathrm{\gamma}}_{\text{i,0}}^{\text{}}+\boldsymbol{J}_{\text{0,H}}\
\boldsymbol{\mathrm{\gamma}}_{\text{i,H}}^{\text{}}=\boldsymbol{\mathrm{0}}\Rightarrow\boldsymbol{\mathrm{\gamma}}_{\text{i,0}}^{\text{}}=-(\boldsymbol{J}_{\text{0,0}})^{-1}\
\boldsymbol{J}_{\text{0,H}}\
\boldsymbol{\mathrm{\gamma}}_{\text{i,H}}^{\text{}}\\\ $ (26)
$\left(\boldsymbol{J}_{\text{H,H}}-\left(\boldsymbol{J}_{\text{H,0}}\
\boldsymbol{J}^{-1}_{\text{0,0}}\ \boldsymbol{J}_{\text{0,H}}\right)\right)\
\boldsymbol{\mathrm{\gamma}}_{\text{i,H}}^{\text{}}-\ \lambda_{i}\
\overline{\boldsymbol{M}}_{\text{H,H}}\
\boldsymbol{\mathrm{\gamma}}_{\text{i,H}}^{\text{}}=\boldsymbol{\mathrm{0}}.$
(27)

Since the approximate frequency $\overline{\omega}$ of the mode shape of
interest is already available from (17), we select from (27) the pair of
eigenvectors $\boldsymbol{\mathrm{\gamma}}_{\text{1,H}}^{\text{}}$,
$\boldsymbol{\mathrm{\gamma}}_{\text{2,H}}^{\text{}}$ with frequencies
$\omega_{1}$, $\omega_{2}$ closest to $\overline{\omega}$. Note that we select
two eigenvectors because we are interested in a response that is dominated by
one harmonic (e.g fundamental harmonic), and the number of modes is 2$n$,
owing to the cosine/sine representation. In case of complex eigenvectors, as
could be already anticipated by the absence of symmetry in
$\boldsymbol{J}^{k}$ in case slipping and/or separation occur, we only need
the real and imaginary parts of one complex mode to construct the
corresponding two-dimensional invariant subspace in $\mathbb{R}^{n(2H)}$ for
the $k$-th linear system. Subsequently, we group the two modes as

$\boldsymbol{\Gamma}_{\text{H}}^{k}=[\boldsymbol{\mathrm{\gamma}}_{\text{1,H}}^{k}\
,\boldsymbol{\mathrm{\gamma}}_{\text{2,H}}^{k}\ ],$ (28)

where $\boldsymbol{\Gamma}_{\text{H}}^{k}\in\mathbb{R}^{n(2H)\times 2}$. As
for the static part of the eigenvectors, it should be highlighted that the JP
method [20] takes into account the coupling between static and dynamic
components in the off-diagonal blocks $\boldsymbol{J}_{\text{0,H}}^{\text{}}$,
$\boldsymbol{J}_{\text{H,0}}^{\text{}}$. Therefore, we compute the static
components of the eigenvectors using (26), and group them as

$\boldsymbol{\Gamma}_{0}^{k}=[\boldsymbol{\mathrm{\gamma}}_{1,0}^{k}\
,\boldsymbol{\mathrm{\gamma}}_{2,0}^{k}\ ],$ (29)

where $\boldsymbol{\Gamma}_{0}^{k}\in\mathbb{R}^{n\times 2}$. Assembling the
components in one matrix $\boldsymbol{\Gamma}^{k}$, we write

$\boldsymbol{\Gamma}^{k}=\begin{bmatrix}\boldsymbol{\Gamma}_{0}^{k}\\\
\boldsymbol{\Gamma}_{\text{H}}^{k}\end{bmatrix},$ (30)

where $\boldsymbol{\Gamma}^{k}\in\mathbb{R}^{n(2H+1)\times 2}$. In comparison
to the original method [20], we augment the JP basis $\boldsymbol{\Gamma}^{k}$
by computing the displacement solutions of the forced systems

$\displaystyle\boldsymbol{Z}^{k}\boldsymbol{\mathrm{\overline{\gamma}}}_{\text{i}}^{k}=\boldsymbol{\mathrm{P}}_{\text{ext,}i},$
(31)

where $i=1,...,a$, and
$\boldsymbol{\mathrm{P}}_{\text{ext,}i}\in\mathbb{R}^{n(2H+1)}$ (see (10)).
Matrix $\boldsymbol{Z}^{k}\in\mathbb{R}^{n(2H+1)\times n(2H+1)}$ involves
cross-harmonic coupling in the stiffness term, unlike (19), and consists of:

$\boldsymbol{Z}^{k}=\begin{bmatrix}\boldsymbol{0}&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{Z}_{\text{MC}}^{\text{}}\end{bmatrix}+\boldsymbol{J}^{k}.$
(32)

The vectors are normalized as:

$\|\boldsymbol{\mathrm{\overline{\gamma}}}_{\text{i}}^{k}\|=1.$ (33)

Note that the LU decomposition of the dense matrix $\boldsymbol{Z}^{k}$ in the
linear system (31) should be carried out only once for each k-system. We then
augment the matrix $\boldsymbol{\Gamma}^{k}$ with
$\boldsymbol{\mathrm{\overline{\gamma}}}_{\text{i}}^{k}$. In fact, these
augmented vectors have proven necessary for an enhanced accuracy on the forces
such that that an error metric on the residual forces, such as what we present
in Section 3.5, could be used on the test cases. In addition, it has proven
essential as well for an increased accuracy on the displacements for the more
complicated test cases in Section 5.

Next, and similar to [20], we apply this procedure for $k=1,...,m$
amplifications to construct $m$ matrices $\boldsymbol{\mathrm{\Gamma}}^{k}$
that are then assembled to form $\boldsymbol{\Gamma}_{\text{}}^{\text{}}$ (see
Algorithm 1). Then, $\boldsymbol{\Gamma}$ is partitioned according to its
harmonic components. Finally, a Singular Value Decomposition (SVD) procedure
is performed on each basis to avoid ill-conditioning, thereby the reduction
basis for $\boldsymbol{V}_{\text{j}}^{\text{c/s}}$ for each harmonic component
is constructed.

As a result, each harmonic component can be represented as

$\boldsymbol{\mathrm{U}}_{\text{j}}^{\text{c/s}}\approx\boldsymbol{V}_{\text{j}}^{\text{c/s}}\boldsymbol{\mathrm{Q}}_{\text{j}}^{\text{c/s}},$
(34)

or, for the whole system, as

$\boldsymbol{\mathrm{U}}\approx\boldsymbol{W}\boldsymbol{\mathrm{Q}},$ (35)

where

$\boldsymbol{W}=\begin{bmatrix}\boldsymbol{V}_{\text{0}}^{\text{}}&&&&&\\\
&\boldsymbol{V}_{\text{1}}^{\text{c}}&&&&\\\
&&\boldsymbol{V}_{\text{1}}^{\text{s}}&&&\\\ &&&\dots&&\\\
&&&&\boldsymbol{V}_{\text{H}}^{\text{c}}&\\\
&&&&&\boldsymbol{V}_{\text{H}}^{\text{s}}\end{bmatrix}.$

Matrix $\boldsymbol{W}\in\mathbb{R}^{n(2H+1)\times p}$ has $p$ total number of
columns in $\boldsymbol{W}$. Note that the bases
$\boldsymbol{V}_{\text{0}}^{\text{}}\in\mathbb{R}^{n\times m_{0}}$,
$\boldsymbol{V}_{\text{1}}^{\text{c}}\in\mathbb{R}^{n\times m_{1c}}$, etc.,
have in general different number of columns. 444We note that an alternative
strategy would be to construct $\boldsymbol{W}$ as a full matrix, in which
there is no longer an association between a set of reduced DOFs and a distinct
harmonic component. In our experience, this has resulted in a low accuracy of
the associated ROM. A possible explanation is that the SVD procedure reduces
then the contributions of higher harmonics inaccurately, since they typically
have small magnitudes. The reduced representation of the displacements in time
domain can now be expressed as

$\begin{split}\boldsymbol{\mathrm{u}}_{\text{h}}^{\text{}}(\text{t})&=\boldsymbol{\mathrm{U}}_{\text{0}}^{\text{}}+\sum_{\text{j}=1}^{\text{H}}\left({\boldsymbol{\mathrm{U}}_{\text{j}}^{\text{c}}\cos(\text{j}\Omega\text{t})+\boldsymbol{\mathrm{U}}_{\text{j}}^{\text{s}}\sin{(\text{j}\Omega\text{t}})}\right)\approx\boldsymbol{\mathrm{V}}_{\text{0}}^{\text{}}\boldsymbol{\mathrm{Q}}_{\text{0}}^{\text{}}+\sum_{\text{j}=1}^{\text{H}}\left(\boldsymbol{\mathrm{V}}_{\text{j}}^{\text{c}}{\boldsymbol{\mathrm{Q}}_{\text{j}}^{\text{c}}\cos({\text{j}\Omega\text{t}})+\boldsymbol{\mathrm{V}}_{\text{j}}^{\text{s}}\boldsymbol{\mathrm{Q}}_{\text{j}}^{\text{s}}\sin{(\text{j}\Omega\text{t}})}\right).\end{split}$
(36)

Next, the subspace of $\boldsymbol{W}{}{}$ is used to construct the ROM via
the Galerkin projection:

$\tilde{\boldsymbol{Z}}\boldsymbol{\mathrm{Q}}+\tilde{\boldsymbol{\mathrm{F}}}(\boldsymbol{\mathrm{WQ}})=\tilde{\boldsymbol{\mathrm{P}}}_{\text{ext}},$
(37)

where

$\displaystyle\tilde{\boldsymbol{Z}}=\boldsymbol{W}_{\text{}}^{\text{T}}\boldsymbol{Z}\boldsymbol{W},$
$\displaystyle\tilde{\boldsymbol{\mathrm{F}}}=\boldsymbol{W}_{\text{}}^{\text{T}}\boldsymbol{\mathrm{F}}(\boldsymbol{\mathrm{WQ}}),$
$\displaystyle\tilde{\boldsymbol{\mathrm{P}}}_{\text{ext}}=\boldsymbol{W}_{\text{}}^{\text{T}}\boldsymbol{\mathrm{P}}_{\text{ext}}^{\text{}}.$

Algorithm 1 Basis construction

0: $a$ excitation amplitudes (2), external force vector
$\boldsymbol{\mathrm{P}}_{\text{ext,}i}$ (10), static solution
$\boldsymbol{\mathrm{u}}_{\text{0}}^{\text{}}$ (15), modal frequency of
interest $\overline{\omega}$ (17), $\boldsymbol{Z}_{\text{MC}}^{\text{}}$
(19), linear dynamic solution
$\boldsymbol{\mathrm{U}}_{\text{\emph{l}}}^{\text{}}$ (18), $m$ amplification
factors (21), multi-harmonic mass matrix $\overline{\boldsymbol{M}}$ (24)

0: reduction basis $\boldsymbol{W}\in\mathbb{R}^{n(2H+1)\times p}$_construct
basis vectors_

1: for $k=1:m$ do

2:
$\boldsymbol{\mathrm{\theta}}^{k}\leftarrow[\boldsymbol{\mathrm{u}}_{\text{0}}^{\text{T}},\
\alpha_{\text{k}}\boldsymbol{\mathrm{U}}_{\text{\emph{l}}}^{\text{}}(\overline{\omega})^{T}]$

3: $\boldsymbol{J}^{k},\boldsymbol{\mathrm{F}}^{k}\leftarrow$ AFT1
$(\boldsymbol{\mathrm{\theta}}^{\text{k}},\overline{\omega})$

4:
$\boldsymbol{\mathrm{\gamma}}_{1}^{k},\boldsymbol{\mathrm{\gamma}}_{2}^{k}\leftarrow$
EigenCond2 $(\overline{\boldsymbol{M}},\boldsymbol{J}^{k},\overline{\omega})$

5:
$\boldsymbol{\Gamma}_{\text{H}}^{k}\leftarrow[\text{Real}\\{{\boldsymbol{\mathrm{\gamma}}_{\text{1,H}}^{k}}\\},\text{Im}\\{{\boldsymbol{\mathrm{\gamma}}_{\text{1,H}}^{k}}\\},\text{Real}\\{{\boldsymbol{\mathrm{\gamma}}_{\text{2,H}}^{k}}\\},\text{Im}\\{{\boldsymbol{\mathrm{\gamma}}_{\text{2,H}}^{k}}\\}]$

6:
$\boldsymbol{\Gamma}_{0}^{k}\leftarrow\texttt{StaticRet}^{3}(\boldsymbol{\Gamma}_{\text{H}}^{k},\boldsymbol{J}^{k})$

7:
$\boldsymbol{\Gamma}^{k}\leftarrow\begin{bmatrix}\boldsymbol{\Gamma}_{0}^{k}\\\
\boldsymbol{\Gamma}_{\text{H}}^{k}\end{bmatrix}$

8: for $i=1:a$ do

9: $\boldsymbol{\mathrm{\overline{\gamma}}}_{\text{i}}^{k}\leftarrow$
AugForced4
$[\boldsymbol{Z}_{\text{MC}}^{\text{}},\boldsymbol{J}^{k},\boldsymbol{\mathrm{P}}_{\text{ext,}i}]$

10:
$\boldsymbol{\Gamma}^{k}\leftarrow[\boldsymbol{\Gamma}^{k},\boldsymbol{\mathrm{\overline{\gamma}}}_{\text{i}}^{k}]$

11: end for _assign reduction basis of each harmonic_

12:
$\boldsymbol{V}_{\text{0}}^{\text{}}\leftarrow[\boldsymbol{V}_{\text{0}}^{\text{}},\boldsymbol{\Gamma}_{0}^{k}]$$\boldsymbol{V}_{\text{j}}^{\text{c}}\leftarrow[\boldsymbol{V}_{\text{j}}^{\text{c}},\boldsymbol{\Gamma}_{\text{j}}^{k\text{,c}}]$$\boldsymbol{V}_{\text{j}}^{\text{s}}\leftarrow[\boldsymbol{V}_{\text{j}}^{\text{c}},\boldsymbol{\Gamma}_{\text{j}}^{k\text{,s}}]$

13: end for

14: $\boldsymbol{V}_{\text{0}}^{\text{}}\leftarrow$ SVD5
$(\boldsymbol{V}_{\text{0}}^{\text{}})$$\boldsymbol{V}_{\text{j}}^{\text{c}}\leftarrow$
SVD
$(\boldsymbol{V}_{\text{j}}^{\text{c}}$)$\boldsymbol{V}_{\text{j}}^{\text{s}}\leftarrow$
SVD $(\boldsymbol{V}_{\text{j}}^{\text{s}})$

1 AFT uses $\boldsymbol{\mathrm{\theta}}^{k}$ as an input displacement to the
AFT algorithm to compute the nonlinear forces $\boldsymbol{\mathrm{F}}^{k}$
and the Jacobian of nonlinear forces $\boldsymbol{J}^{k}$.
2 EigenCond solves the eigenvalue problem (22) after applying the static
condensation procedure using (26) and (27).
3 StaticRet retrieves the static components of the eigenvectors (26).
4 AugForced computed the forced response of the $k$-th linear system according
to (31).
5 SVD performs the singular value decomposition on the basis of each harmonic
component.

### 3.3. Hyper Reduction by Energy Conserving, Sampling, and Weighing (HR
ECSW)

Even though the size of the Augmented JP ROM in (37) is significantly smaller
than that of the original system, the speedups remain limited due to the high
density of the mesh of nonlinear elements across the interface and the
computational cost associated to the high number of sampling points required
in the AFT scheme. Hyper Reduction (HR) refers to techniques developed to
overcome the computational burden associated with the evaluation of the vector
of nonlinear forces. In this section, we briefly review the Energy Conserving
Mesh Sampling and Weighting (ECSW) method [25] method, and explain why we find
it suitable for jointed structures. For these purposes, let us consider a
data-driven ROM constructed in the time-domain. Using HF simulations, a
reduction basis $\boldsymbol{S}\in\mathbb{R}^{n\times r}$, which approximates
the displacements as

$\boldsymbol{\mathrm{u}}\text{(t)}\approx\boldsymbol{S}\boldsymbol{\mathrm{q}}\text{(t)},$
(38)

can be constructed using, for instance, Proper Orthogonal Decomposition (POD),
where $\boldsymbol{\mathrm{q}}\text{(t)}$ is the vector of reduced
coordinates. Subsequently, a Galerkin projection of the system onto
$\boldsymbol{S}$ can be carried out. Since the dynamics of joints are known to
be remarkably dissipative, it would be desirable to still preserve the
dissipation of the system upon the approximation of the nonlinear forces. The
dissipation over one time period due to frictional contact is

$\int_{\text{t}}^{\text{t+T}}\boldsymbol{\mathrm{f(t)}}_{\text{}}^{\text{T}}\boldsymbol{\mathrm{\dot{u}(t)}}\
\text{dt}\approx\int_{\text{t}}^{\text{t+T}}\boldsymbol{\mathrm{f(t)}}_{\text{}}^{\text{T}}\boldsymbol{S}\boldsymbol{\mathrm{\dot{q}(\text{t})}}\
\text{dt}=\int_{\text{t}}^{\text{t+T}}\boldsymbol{\mathrm{\dot{q}(\text{t})}}_{\text{}}^{\text{T}}\boldsymbol{S}_{\text{}}^{\text{T}}\boldsymbol{\mathrm{f(t)}}\
\text{dt}.$ (39)

Accurately approximating the projection of the nonlinear forces
$\boldsymbol{S}_{\text{}}^{\text{T}}\boldsymbol{\mathrm{f(t)}}$ would
therefore mean preserving the dissipation computed by the ROM. Indeed, this is
what ECSW aims at, as it seeks an energetically equivalent approximation of
the projection of the nonlinear forces
$\tilde{\boldsymbol{f}}\text{(t)}=\boldsymbol{S}_{\text{}}^{\text{T}}\boldsymbol{f}\text{(t)}$,
such that

$\tilde{\boldsymbol{f}}=\sum_{e=1}^{\text{n}_{\text{e}}^{\text{}}}\boldsymbol{S}_{\text{e}}^{\text{T}}\boldsymbol{\mathrm{f}}_{\text{e}}^{\text{}}(\boldsymbol{S}_{\text{e}}^{\text{}}\boldsymbol{\mathrm{q}}_{\text{}}^{\text{}})\approx\sum_{e\in
E}\xi_{\text{e}}\boldsymbol{S}_{\text{e}}^{\text{T}}\boldsymbol{\mathrm{f}}_{\text{e}}^{\text{}}(\boldsymbol{S}_{\text{e}}^{\text{}}\boldsymbol{\mathrm{q}}_{\text{}}^{\text{}})=\tilde{\boldsymbol{f}}_{\text{HR}},$
(40)

where $\text{n}_{\text{e}}^{\text{}}$ is the total number of nonlinear
elements, and E represents only a small subset of elements for which the
weights $\xi_{\text{e}}$ satisfy the condition $\xi_{e\in E}>0$. The training
required to compute the vector
$\boldsymbol{\mathrm{\xi}}\in\mathbb{R}^{n_{e}}$, which consists of the
weights of all nonlinear elements, is based on $n_{s}$ force snapshots
$\begin{bmatrix}\boldsymbol{\mathrm{f}}^{(i)}\end{bmatrix}_{i=1}^{n_{s}}$.
Referring to the original method, these are snapshots of the nonlinear forces
that are extracted from a HF time integration scheme. To compute the weights,
a Non-Negative Least Squares (NNLS) optimization problem is set as

$\boldsymbol{\mathrm{\xi}}^{*}=\text{argmin}_{\boldsymbol{\mathrm{\xi}}\in\mathbb{R}^{n_{e}},\boldsymbol{\mathrm{\xi}}\geq
0}\|\boldsymbol{G}\boldsymbol{\mathrm{\xi}}-\boldsymbol{\mathrm{b}}\|_{2}^{2},$
(41)

where

$\boldsymbol{G}=\left.\begin{bmatrix}{\boldsymbol{\mathrm{g}}}_{1,1}&\dots&\boldsymbol{\mathrm{g}}_{1,n_{e}}\\\
\vdots&\ddots&\vdots\\\
\boldsymbol{\mathrm{g}}_{n_{s},1}&\dots&\boldsymbol{\mathrm{g}}_{n_{s},\text{n}_{\text{e}}^{\text{}}}\\\
\end{bmatrix}\right\\}\text{$n_{s}$ snapshots,}$ (42)

with

$\begin{split}\boldsymbol{g}_{\text{i,e}}^{\text{}}=\boldsymbol{S}_{\text{e}}^{\text{T}}\boldsymbol{\mathrm{f}}_{\text{e}}^{\text{}}(\boldsymbol{S}_{\text{e}}^{\text{}}\boldsymbol{\mathrm{q}}^{\text{i}}),\\\
\boldsymbol{\mathrm{b}}_{\text{i}}^{\text{}}=\sum_{\text{e}}^{\text{n}_{\text{e}}}{\boldsymbol{g}_{\text{i,e}}^{\text{}}},\end{split}$
(43)

where $\boldsymbol{G}\in\mathbb{R}^{r(n_{s})\times n_{e}}$,
$\boldsymbol{\mathrm{b}}\in\mathbb{R}^{r(n_{s})\times 1}$,
$\boldsymbol{\mathrm{g}}_{\text{i,e}}^{\text{}}\in\mathbb{R}^{r\times 1}$, and
$\boldsymbol{\mathrm{b}}_{\text{i}}^{\text{}}\in\mathbb{R}^{r\times 1}$.
However, the optimal solution to (41) is NP-hard. Therefore, a sub-optimal,
yet feasible, algorithm that is known as the sparse NNLS [35], which
terminates the search for a solution once the objective function returns a
value below a certain threshold, is used. This threshold is controlled by a
user-set tolerance $\tau$, and is given by

$\|\boldsymbol{G}\boldsymbol{\mathrm{\xi}}-\boldsymbol{\mathrm{b}}\|_{2}\leq\tau\|\boldsymbol{\mathrm{b}}\|_{2}.$
(44)

### 3.4. Training for Frictional Contact in Frequency Domain

We propose HR training that does not rely on HF simulations, and is in line
with frequency domain analysis in MHBM. Intuitively, the information conveyed
in a set of multiple time snapshots of the solutions can also be contained in
the harmonic components of the nonlinear force vector. Therefore, for a given
amplitude of excitation, we use the harmonic components of
$\boldsymbol{\mathrm{F}}^{k}$, which is an output of the the AFT in step 2 of
Algorithm 1, as our 2H+1 force snapshots. This means that this training set
does not incur any additional costs after the construction of the reduction
basis. Considering the reduction basis $\boldsymbol{W}$ and the Augmented JP
ROM (37), the approximation we seek becomes

$\tilde{\boldsymbol{F}}=\sum_{e=1}^{\text{n}_{\text{e}}^{\text{}}}\boldsymbol{W}_{\text{e}}^{\text{T}}\boldsymbol{\mathrm{F}}_{\text{e}}^{\text{}}(\boldsymbol{W}_{\text{e}}^{\text{}}\boldsymbol{\mathrm{Q}}_{\text{}}^{\text{}})\approx\sum_{e\in
E}\xi_{\text{e}}\boldsymbol{W}_{\text{e}}^{\text{T}}\boldsymbol{\mathrm{F}}_{\text{e}}^{\text{}}(\boldsymbol{W}_{\text{e}}^{\text{}}\boldsymbol{\mathrm{Q}}_{\text{}}^{\text{}})=\tilde{\boldsymbol{F}}_{\text{HR}}.$
(45)

The same sparse NNLS algorithm [35] is used to determine
$\boldsymbol{\mathrm{\xi}}$. However, we define an alternative matrix
$\boldsymbol{\tilde{G}}$ as

$\boldsymbol{\tilde{G}}=\left.\begin{bmatrix}\tilde{\boldsymbol{\mathrm{g}}}_{1,1}&\dots&\tilde{\boldsymbol{\mathrm{g}}}_{1,\text{n}_{\text{e}}^{\text{}}}\\\
\vdots&\ddots&\vdots\\\
\tilde{\boldsymbol{\mathrm{g}}}_{\text{(2H+1),1}}&\dots&\tilde{\boldsymbol{\mathrm{g}}}_{\text{(2H+1)},\text{n}_{\text{e}}^{\text{}}}\\\
\end{bmatrix}\right\\}(\text{$2H+1$) snapshots,}$ (46)

where

$\begin{split}\tilde{\boldsymbol{g}}_{\text{1,e}}&=[\boldsymbol{V}_{\text{0,e}}^{\text{}}]^{\text{T}}[\boldsymbol{\mathrm{F}}_{\text{e}}^{\text{}}]_{0},\\\
\tilde{\boldsymbol{g}}_{\text{2,e}}&=[\boldsymbol{V}_{\text{1,e}}^{\text{c}}]^{\text{T}}[\boldsymbol{\mathrm{F}}_{\text{e}}^{\text{}}]_{1}^{c},\\\
\tilde{\boldsymbol{g}}_{\text{3,e}}&=[\boldsymbol{V}_{\text{1,e}}^{\text{s}}]^{\text{T}}[\boldsymbol{\mathrm{F}}_{\text{e}}^{\text{}}]_{1}^{s},\end{split}$
(47)

and $\boldsymbol{\mathrm{\tilde{b}}}_{\text{i}}^{\text{}}$ then becomes

$\boldsymbol{\mathrm{\tilde{b}}}_{\text{i}}=\sum_{\text{e}=1}^{\text{n}_{\text{e}}^{\text{}}}{\boldsymbol{\tilde{g}}_{\text{i,e}}}.$
(48)

Regarding the value of $\tau$, a value of 0.01 has been found to lead to
sparse hyper meshes capable of producing accurate results. Figure 3 shows a
flowchart of the hyper reduction step.

Figure 2: Flowchart of the method 1/2 - Basis Construction at
$\Omega=\overline{\omega}$. Figure 3: Flowchart of the method 2/2 - Hyper
Reduction determined at $\Omega=\overline{\omega}$. The same hyper mesh and
weights are used then across the frequency sweep.

### 3.5. Error Metric

Since we aim to construct reliable ROMs which no longer require HF
simulations, we propose an error metric on the ROM solution that does not rely
on a HF reference solution. Since approximating the displacements using (5)
and (36) leads to residual forces $\boldsymbol{\mathrm{r}}(\text{t})$, where

$\boldsymbol{\mathrm{r}}(\text{t})=\boldsymbol{M}\boldsymbol{\mathrm{\ddot{u}_{h}}}(\text{t})+\boldsymbol{C}\boldsymbol{\mathrm{\dot{u}_{h}}}(\text{t})+\boldsymbol{K}\boldsymbol{\mathrm{u}}_{\text{h}}^{\text{}}(\text{t})+\boldsymbol{\mathrm{f}}(\boldsymbol{\mathrm{u}}_{\text{h}}^{\text{}}(\text{t}),\boldsymbol{\mathrm{\dot{u}}}_{\text{h}}^{\text{}}(\text{t}))-\boldsymbol{\mathrm{P}}_{\text{ext,}i}(\text{t}),$
(49)

we therefore evaluate the error $E$ on the residual forces at the DOF of
excitation (subscript y) using a root mean square (RMS) measure as

$E=\frac{\sqrt{\frac{1}{N}\sum\limits_{\text{i}}^{N}[\text{r}_{\text{y}}^{\text{}}(\text{t}_{\text{i}}^{\text{}})]^{2}}}{\sqrt{\frac{1}{N}\sum\limits_{\text{i}}^{N}[\text{p}_{\text{E,y}}^{\text{}}(\text{t}_{\text{i}}^{\text{}})]^{2}}},$
(50)

where $N$ is the number of sampling points in the time domain. It will be
shown in Section 4 and 5, that for the test cases we studied, errors below 15%
indicate a HR-ROM solution with an accuracy capable of capturing the global
amplitude-dependent stiffness and damping of the system.

## 4\. Case Study: Jointed Beam

### 4.1. FE Model

Figure 4 shows a prototype structure of the Brake-Reuß beam [1]. The structure
is composed of two beams joined together by the preclamping action of 3 bolts.
Each beam has dimensions 42x2.5x2.5 cm. The material is steel with Young’s
Modulus E = 189 GPa and density $\rho$ =7820 Kg/$\text{m}^{3}$. The beams are
modelled using planar Euler-Bernoulli beam elements. The contact interface is
12 cm long, and the frictional contact is modelled using an array of node-to-
node Jenkins elements with unilateral springs. The contact parameters are
assumed as $\text{k}_{\text{t}}=3\times 10^{11}\ \text{N/m}^{2}$,
$\text{k}_{\text{n}}=4\times 10^{12}\ \text{N/m}^{2}$, and $\mu$=0.4. Each
bolt exerts a force of 1.25 KN, which is modelled as concentrated forces
applied at nodes lying within the tributary area of each bolt. The material
damping is modelled as Rayleigh damping with a modal damping ratio
$\zeta=3\times 10^{-3}$. The excitation force acts in the y-direction, with
amplitudes $\text{A}_{\text{}}^{\text{}}=0.1,2,5,10$ N, frequencies $\Omega$
ranging from 205 to 230 Hz, and is positioned at a distance 22 cm from the
right end. We assume 5 harmonics in the Ansatz. The analysis is carried out
using an in-house MATLAB code and implements sequential continuation555An
attempt to apply a numerical continuation method (e.g. arc-length
parameterization) for the HR-ROM did not lead to converging results. Unlikely
to be a limitation, since the same procedure is problem-free for the next test
case we will present, we report this experience as an anomaly. It should also
be kept in mind that a failure in convergence is a flag to the analyst, as
opposed to an inaccurate result, which if present in a converged solution,
would be detected by the error metric. Lastly, the application of arc-length
continuation for the HF simulation of the beam test case did not show any new
insights compared to sequential continuation, therefore the presented results
have been concluded valid.. We monitor the ratio between the maximum
displacement at the DOF of excitation and the amplitude of excitation.

We perform our analysis using two different meshes, namely:

  * •

Mesh (1): 121 Jenkins elements for the contact interface and a total of 300
beam elements.

  * •

Mesh (2): 241 Jenkins elements for the contact interface and a total of 540
beam elements.

Figure 4: Jointed beam - Test case 1

### 4.2. Results

Figure 6 shows the FRFs for Mesh (1) around the first bending mode using the
HFM (9900 DOFs), and the HR-ROM (75 DOFs). As shown, the results of HR-ROM
closely match those of the HFM, and the online speedups range from 5.5 to 7.0,
as shown in Table 1. The offline costs reported in Table 3. Note that the
basis is constructed only once regardless of the number of analysis performed,
and that the total time for the training for hyper reduction is less than 1
second. As for the hyper-mesh, Table 3 shows that the ratio of hyper-elements
to total number of Jenkins elements ranges from 12% to 24%. Using the error
metric, the maximum error around resonance for A=10 N is less than 9%. Since
the error metric considers residual forces, and it is known that the
convergence on the forces is slow compared to the displacements [8], one
should in general expect the presence of errors even in the case of
overlapping FRFs. Nonetheless, the ROM has the capacity for reproducing the
evolution of friction forces with good accuracy. For instance, the
distributions of friction forces at different time instances for A= 10N and
$\Omega$= 216.2 Hz, according to the HFM and HR-ROM, respectively, are plotted
in Figure 8. This case corresponds to an error of 3.25% for the HR-ROM;
nevertheless, the agreement is good.

Figure 5: FRF for Mesh (1) using HF and HR-ROM

Figure 6: FRFs for Mesh (1) and Mesh (2) using HR-ROM

Figure 7: Errors around resonance for A= 10N using HR-ROM for Mesh (1) and
Mesh (2). The red square locates an error of 3.25%, for which the
corresponding response along the contact interface is plotted in Figure 8.

To test the convergence of the model with respect to the number of nonlinear
elements, we analyze the structure using Mesh (2). The total number of reduced
DOFs increases only marginally, from 75 to 77. Figure 6 indicates that Mesh
(1) is already convergent, whereas Table 1 shows speedups that now range from
9.5 to 42.8. Noting the significant increase in speedups for the finer mesh,
we highlight another advantage of the method: the computational burden of
simulations aimed to determine an appropriate mesh size for a structure, whose
response is unforeseen, is largely alleviated. Indeed, although the nonlinear
elements in Mesh (2) are twice as much as in Mesh (1), the hyper mesh only
increases by around 10% for A= 10N, for example, as shown in Table 3. In
practical situations, it is not trivial to determine apriori the right mesh
size for a desired accuracy, and often tedious mesh convergence studies need
to be carried out. The example discussed here shows that doubling the number
of nonlinear elements results in a tenfold more expensive HFM. The HR-ROM
counterpart, however, features a relatively marginal increase in solution time
from 95.1s to 129.7s. This is mainly due to hyper-reduction, which selects
just enough elements to represent the nonlinearity, irrespective of the
underlying mesh density.

Figure 8: Snapshots of friction forces along the interface for Mesh (1) across
the time period T due to $\text{A}=\ 10\text{N}$, $\Omega=\ 216.2\text{Hz}$

.

| Mesh (1) | Mesh (2) |
---|---|---|---
Force [N] | HR [s] | HF [s] | Speedups | HR [s] | HF[s] | Speedups | Ratio of Speedups
0.1 | 29.5 | 206.4 | 7.0 | 57.7 | 550.7 | 9.5 | 1.4
2 | 68.7 | 404.2 | 5.9 | 100.9 | 1864 | 18.5 | 3.1
5 | 85.7 | 473.5 | 5.5 | 115.7 | 3217.6 | 27.8 | 5.0
10 | 95.1 | 520.4 | 5.5 | 129.7 | 5547.8 | 42.8 | 7.8
Table 1: HF and ROM run times for Mesh (1) and Mesh (2)

Offline cost | Mesh (1) | Mesh (2)
---|---|---
Total training time [s] | 0.35 | 0.86
Basis time [s] | 28.4 | 97.5
Table 2: ROM offline costs for Mesh (1) and Mesh (2)

F [N] | Hyper Mesh (1) | Hyper Mesh (2)
---|---|---
0.1 | 14 | 19
2 | 23 | 26
5 | 27 | 29
10 | 29 | 32
Table 3: Hyper mesh for Mesh (1) and Mesh (2)

## 5\. Case Study: Frame

### 5.1. FE Model

Figure 10 shows a frame that is composed of a vertical column and a beam, both
of length 2.80m, and made of a steel HEA 300 cross-section, with an area of
$112\ \text{cm}^{2}$, and a second moment of area of $18260\ \text{cm}^{4}$.
The column and the beam are joined by a bolted steel bracket having a T-cross
section, whose area is $37.5\ \text{cm}^{2}$ and second moment of area of
$3652\ \text{cm}^{4}$. The structure is modeled using planar Euler-Bernoulli
beam elements. Each contact interface is 40 cm long, is subject to the
preclamp action of 5 bolts, and discretized using 121 Jenkins elements. The FE
discretization results in 1743 spatial DOFs. Each bolt exerts a force of 6.75
tons, which is modelled in terms of concentrated forces applied at the nodes
existing within the tributary area of each bolt. The contact parameters are
assumed as $\text{k}_{\text{t}}^{\text{}}=3e11\ \text{N/m}^{2}$,
$\text{k}_{\text{n}}^{\text{}}=4e12\ \text{N/m}^{2}$, and $\mu=0.3$. The
excitation force is applied at a distance 1.4m from the fixed end, with
amplitudes ranging from 100N to 1500N. The excitation frequency targets the
second mode shape of the linearized structure, which is displayed in Figure
10. The modal damping ratio is 3e-3. Numerical continuation is used to solve
the MHBM system of equations, using NLvib [33]. The continuation uses tangent
predictor and arc-length corrector steps. Concerning the number of harmonics,
we carried out a convergence study which showed that an Ansatz of 3 harmonics
delivers satisfactory accuracy. The total number of DOFs for the MHBM
equations becomes 12201. The HR tolerance $\tau$ is set to $10e-3$. Two
reduction basis are considered for this study, namely, Basis (1) with 4
amplifications (i.e., $\alpha_{\text{k}=1,..,4}$ ), and Basis (2) with 6
amplifications. These are denoted in the tables and figures by B(1) and B(2).

Figure 9: Frame structure- Test case 2

Figure 10: $2^{\text{nd}}$ mode shape

### 5.2. Results

Figure 12 shows the FRFs of the HFM and HR ROM using Basis (1). The total
number of reduced DOFs is 90 in this case. Even though the agreement is good
for low amplitudes, notable deviations of the responses, in particular, at
A=1500N are shown. Since the ROM should be operable without resorting to HFM,
one should be able to detect such inaccuracies through the error metric. The
errors, shown in Figure 14, reach 21% around resonance. Furthermore, one
should expect converging results upon the enrichment of the basis of the ROM.
Consequently, two additional amplifications are included to construct the
second reduction basis B(2). The total number of reduced DOFs increases from
90 to 110. As shown in Figure 12, the FRFs show significant improvement.
Moreover, the overall errors are shown in Figure 14 to decrease with the
maximum error around resonance falling below 15%. In fact, an error threshold
of 15% has repeatedly proved indicative of a reliable depiction of the
amplitude-dependent dissipation and stiffness for the test cases we studied.

As discussed in Test Case (1), the forces are expected to exhibit slower
convergence. Accordingly, a lower threshold for an accurate spatial and
temporal evolution of frictional contact forces was shown to be around 5%. For
instance, the time snapshots of friction forces for $\Omega\approx 267.8$ Hz
are plotted in Figure 13 along both interfaces. This is to say, that depending
on the purpose of the analysis, the user can enrich and accordingly assess the
ROM. The online and offline costs are reported in Tables 4 & 6. The speedups
using Basis (2) range from 11.9 to 18.5. Lastly, the number of hyper elements
for each case is shown in Table 6.

Figure 11: FRFs for HF vs HR - Basis (1)

Figure 12: FRFs for HF vs HR - Basis (2)

(a) Column interface

(b) Beam interface

(c)

Figure 13: Time snapshots of friction forces along the contact interface of the column (top) and of the beam (bottom) for A= 1500N at $\Omega$= 267.83 Hz for HR-ROM (red rectangle in Figure 14), and $\Omega$=267.88 Hz for HFM Figure 14: Errors computed around resonance for A= 1500N using HR-ROM for Basis (1) and Basis (2). The black circle locates the error corresponding to $\Omega$=267.83 Hz whose corresponding response along the contact interface is plotted in Figure 13. |  | Basis 1 (90 DOFs) | Basis 2 (110 DOFs)
---|---|---|---
Force [N] | HF [s] | HR [s] | Speedups | HR [s] | Speedups
100 | 2761 | 156 | 17.7 | 111 | 18.5
500 | 6030 | 407 | 14.8 | 501 | 13.2
1000 | 11329 | 937 | 12.1 | 973 | 12.1
1500 | 16078 | 1269 | 12.7 | 1335 | 11.9
Table 4: HF and ROM run times for Basis (1) and Basis (2)

Offline cost [s] | Basis (1) | Basis (2)
---|---|---
Total training time | 1.38 | 1.76
Basis time | 41.63 | 65.1
Table 5: ROM offline costs for Basis (1) and Basis (2)

F [N] | Basis (1) | Basis (2)
---|---|---
100 | 43 | 47
500 | 52 | 57
1000 | 68 | 66
1500 | 68 | 73
Table 6: Hyper mesh for Basis (1) and Basis (2)

## 6\. Conclusions

We presented a hyper-reduced order model for steady-state analysis of jointed
structures. The reduction basis is constructed using the Augmented Jacobian
Projection method. This involves solutions to a set of autonomous and forced
linear problems, corresponding to different contact states covering the
expected range of nonlinearity. After the projection of the dynamical system,
we adapt the energy conserving, sampling and weighting method for Hyper
Reduction (HR) in the frequency domain and use training vectors that are
already available from the construction of the Reduced Order Model (ROM). An
error indicator that is based on the HR-ROM solution was presented. The main
contributions of this work are:

  * •

A HR-ROM that is capable of capturing not only global features such as overall
dissipation and stiffness changes with high accuracy, but also local
quantities, whose convergence is known to be more challenging, such as the
evolution of friction forces along the interface.

  * •

An error metric that depends solely on the HR-ROM solution, with reported
thresholds that have proved reliable for the different analyses we performed.

  * •

Allowing convergence studies with respect to interface mesh size for steady-
state analysis. This is possible by virtue of the sparse selection of
nonlinear element, using the proposed HR scheme.

Regarding the limitations of the method, since the construction of the basis
involves the selection of a mode shape of interest, the ROM remains valid
across a limited range of frequencies around that corresponding eigenfrequency
of the linearized structure. An extension of this range can be sought using
linear mode shapes of the off-resonant modes, and further, by considering
basis vectors that represent other modes of interest in the reduction basis.

#### Declaration of Competing Interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

#### Acknowledgments

The authors acknowledge the funding of the Swiss National Science Foundation
project “Meso-scale modeling of Friction in reduced non-linear interface
Dynamics: MesoFriDy”.

## References

  * [1] Matthew R W Brake.  The Mechanics of Jointed Structures.  2018\.
  * [2] A. Fantetti, L. R. Tamatam, M. Volvert, I. Lawal, L. Liu, L. Salles, M. R.W. Brake, C. W. Schwingshackl, and D. Nowell.  The impact of fretting wear on structural dynamics: Experiment and Simulation.  Tribology International, 138(February):111–124, 2019.
  * [3] Jie Yuan, Alfredo Fantetti, Enora Denimal, Shubham Bhatnagar, Luca Pesaresi, Christoph Schwingshackl, and Loic Salles.  Propagation of friction parameter uncertainties in the nonlinear dynamic response of turbine blades with underplatform dampers.  Mechanical Systems and Signal Processing, 156:107673, 2021.
  * [4] L. Gaul and J. Lenz.  Nonlinear dynamics of structures assembled by bolted joints.  Acta Mechanica, 125(1-4):169–181, 1997.
  * [5] Allen T. Mathis, Nidish N. Balaji, Robert J. Kuether, Adam R. Brink, Matthew R.W. Brake, and D. Dane Quinn.  A Review of Damping Models for Structures with Mechanical Joints.  Applied Mechanics Reviews, 72(4):1–29, 2020.
  * [6] Christian Maria and Stefano Zucca.  Modelling Friction Contacts in Structural Dynamics and its Application to Turbine Bladed Disks.  Numerical Analysis - Theory and Application, (September), 2011.
  * [7] M. R.W. Brake, J. G. Stark, S. A. Smith, D. P.T. Lancereau, T. W. Jerome, and T. Dossogne.  In situ measurements of contact pressure for jointed interfaces during dynamic loading experiments.  Conference Proceedings of the Society for Experimental Mechanics Series, 4:133–141, 2017.
  * [8] Malte Krack, Loic Salles, and Fabrice Thouverez.  Vibration Prediction of Bladed Disks Coupled by Friction Joints.  Archives of Computational Methods in Engineering, 24(3):589–636, 2017.
  * [9] Walter C. Hurty.  Dynamic analysis of structural systems using component modes.  AIAA Journal, 3(4):678–685, 1965.
  * [10] Roy R. Craig and Mervyn C.C. Bampton.  Coupling of substructures for dynamic analyses.  AIAA Journal, 6(7):1313–1319, 1968.
  * [11] Sergio R. Idelsohn and Alberto Cardona.  A load-dependent basis for reduced nonlinear structural dynamics.  Computers and Structures, 20(1-3):203–210, 1985.
  * [12] P. M.A. Slaats, J. de Jongh, and A. A.H.J. Sauren.  Model reduction tools for nonlinear structural dynamics.  Computers and Structures, 54(6):1155–1171, 1995.
  * [13] Long Wu and Paolo Tiso.  Nonlinear model order reduction for flexible multibody dynamics: a modal derivatives approach.  Multibody System Dynamics, 36(4):405–425, 2016.
  * [14] Florian Pichler, Wolfgang Witteveen, and Peter 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(5):1–12, 2017\.
  * [15] Florian Pichler, Wolfgang Witteveen, and Peter 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.
  * [16] Patrick J. Hughes and Robert J. Kuether.  Nonlinear Interface Reduction for Time-Domain Analysis of Hurty/Craig-Bampton Superelements with Frictional Contact.  Journal of Sound and Vibration, 507:116154, 2021.
  * [17] Malte Krack.  Nonlinear modal analysis of nonconservative systems: Extension of the periodic motion concept.  Computers and Structures, 154:59–71, 2015.
  * [18] Malte Krack.  The Concept of Nonlinear Modes Applied to Friction-Damped Systems.  2014\.
  * [19] Mainak Mitra, Stefano Zucca, and Bogdan I. Epureanu.  Adaptive Microslip Projection for Reduction of Frictional and Contact Nonlinearities in Shrouded Blisks.  Journal of Computational and Nonlinear Dynamics, 11(4):1–15, 2016\.
  * [20] Chiara Gastaldi, Stefano Zucca, and Bogdan I. Epureanu.  Jacobian projection reduced-order models for dynamic systems with contact nonlinearities.  Mechanical Systems and Signal Processing, 100:550–569, 2018.
  * [21] Nidish Narayanaa Balaji, Tobias Dreher, Malte Krack, and Matthew R.W. Brake.  Reduced order modeling for the dynamics of jointed structures through hyper-reduced interface representation.  Mechanical Systems and Signal Processing, 149:107249, 2021.
  * [22] Jie Yuan, Loic Salles, Fadi El Haddad, and Chian Wong.  An adaptive component mode synthesis method for dynamic analysis of jointed structure with contact friction interfaces.  Computers and Structures, 229:106177, 2020.
  * [23] Jie Yuan, Christoph Schwingshackl, Chian Wong, and Loïc Salles.  On an improved adaptive reduced-order model for the computation of steady-state vibrations in large-scale non-conservative systems with friction joints.  Nonlinear Dynamics, 103(4):3283–3300, 2021.
  * [24] Saifon Chaturantabut and Danny C. Sorensen.  Nonlinear model reduction via discrete empirical interpolation.  SIAM Journal on Scientific Computing, 32(5):2737–2764, 2010.
  * [25] Charbel Farhat, Philip Avery, Todd Chapman, and Julien Cortial.  Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy-based mesh sampling and weighting for computational efficiency.  International Journal for Numerical Methods in Engineering, 98(9):625–662, 2014.
  * [26] Shobhit Jain and Paolo Tiso.  Simulation-free hyper-reduction for geometrically nonlinear structural dynamics: A quadratic manifold lifting approach.  Journal of Computational and Nonlinear Dynamics, 13(7):1–12, 2018\.
  * [27] J. B. Rutzmoser and D. J. Rixen.  A lean and efficient snapshot generation technique for the Hyper-Reduction of nonlinear structural dynamics.  Computer Methods in Applied Mechanics and Engineering, 325:330–349, 2017.
  * [28] E. P. Petrov and D. J. Ewins.  Analytical formulation of friction interface elements for analysis of nonlinear multi-harmonic vibrations of bladed disks.  Journal of Turbomachinery, 125(2):364–371, 2003.
  * [29] Georg Masing.  Zur Heyn’schen Theorie der Verfestigung der Metalle durch verborgen elastische Spannungen.  Wissenschaftliche Veröffentlichungen aus dem Siemens-Konzern, pages 231–239, 1923.
  * [30] P Jayakumar.  Modeling and Identification in Structural Dynamics, 1987.
  * [31] T. M. Cameron and J. H. Griffin.  An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems.  Journal of Applied Mechanics, Transactions ASME, 56(1):149–154, 1989.
  * [32] Alberto Cardona, Thierry Coune, Albert Lerusse, and Michel Geradin.  A multiharmonic method for non‐linear vibration analysis.  International Journal for Numerical Methods in Engineering, 37(9):1593–1608, 1994.
  * [33] Malte Krack and Johann Gross.  Harmonic Balance for Nonlinear Vibration Problems.  2019\.
  * [34] Dominik Süß and Kai Willner.  Investigation of a jointed friction oscillator using the Multiharmonic Balance Method.  Mechanical Systems and Signal Processing, 52-53(1):73–87, 2015\.
  * [35] Robert Peharz and Franz Pernkopf.  Sparse nonnegative matrix factorization with $\emph{l}^{0}$ -constraints.  Neurocomputing, 80:38–46, 2012.


$(g,m)\cdot(h,n)=(gh,n),\qquad\forall m,n\in M,g,h\in G,$

whenever $m=\Phi(h)(n)$.

### References

  * [1] Angulo, C., A cohomology theory for Lie $2$-algebras and Lie $2$-groups, Ph. D. Thesis, S$\mathrm{\tilde{a}}$o Paulo: Instituto de Mathemática e Estatística, Universidade de S$\mathrm{\tilde{a}}$o Paulo, 2018.
  * [2] Baez, J. and Crans, A., Higher-dimensional algebra VI: Lie $2$-algebras, Theory Appl. Categ., 2014, 12: 492-538.
  * [3] Baez, J., Hoffnung, A. and Rogers, C., Categorified symplectic geometry and the classical string, Comm. Math. Phys., 2010, 293(3): 701-725.
  * [4] Baez, J. and Lauda, A., Higher-dimensional algebra V: $2$-groups, Theory Appl. Categ., 2014, 12: 423-491.
  * [5] Baez, J. and Schreiber, U., Higher gauge theory, In: Categories in Algebra, Geometry and Mathematical Physics, Contemp. Math., Vol. 431, Providence, RI: AMS, 2007, 7-30.
  * [6] Baez, J., Stevenson, D., Crans, A. and Schreiber, U., From loop groups to $2$-groups, Homology Homotopy Appl., 2007, 9(2): 101-135.
  * [7] Bai, C., Sheng, Y. and Zhu, C., Lie $2$-bialgebras, Comm. Math. Phys., 2013, 320(1): 149-172.
  * [8] Berwick-Evans, D. and Lerman, E., Lie $2$-algebras of vector fields, arXiv: 1609.03944.
  * [9] Bonechi, F., Ciccoli, N., Laurent-Gengoux, C. and Xu, P., Shifted Poisson structures on differentiable stacks, arXiv: 1803.06685.
  * [10] Cai, X., Liu, Z. and Xiang, M., Cohomology of hemistrict Lie $2$-algebras, Comm. Algebra, 2020, 48(8): 3315-3341.
  * [11] Callies, M., Frégier, Y., Rogers, C. and Zambon, M., Homotopy moment maps, Adv. Math., 2016, 303: 954-1043.
  * [12] Cattaneo, A. and Felder, G., Relative formality theorem and quantisation of coisotropic submanifolds, Adv. Math., 2007, 208(2): 521-548.
  * [13] Chen, S., Sheng, Y. and Zheng, Z., Non-abelian extensions of Lie 2-algebras, Sci. China Math., 2012, 55(8): 1655-1668.
  * [14] Chen, Z., Stiénon, M. and Xu. P., Weak Lie $2$-bialgebras, J. Geom. Phys., 2013, 68: 59-68.
  * [15] Chen, Z., Stiénon, M. and Xu. P., Poisson 2-groups, J. Differential Geom., 2013, 94(2): 209-240.
  * [16] Courant, T., Dirac manifolds, Trans. Amer. Math. Soc., 1990, 319(2): 631-661.
  * [17] Crainic, M. and Fernandes, R., Integrability of Lie brackets, Ann. of Math. (2), 2003, 157(2): 575-620.
  * [18] Dorfman, I., Dirac Structures and Integrability of Nonlinear Evolution Equations, Nonlinear Science: Theory and Applications, Chichester: John Wiley-Sons, 1993.
  * [19] Doubek, M. and Lada, T., Homotopy derivations, J. Homotopy Relat. Struct., 2016, 11(3): 599-630.
  * [20] Frégier, Y., Laurent-Gengoux, C. and Zambon, M., A cohomological framework for homotopy moment maps, J. Geom. Phys., 2015, 97: 119-132.
  * [21] Gerstenhaber, M., A uniform cohomology theory for algebras, Proc. Nat. Acad. Sci. U.S.A., 1964, 51: 626-629.
  * [22] Getzler, E., Lie theory for nilpotent $L_{\infty}$-algebras, Ann. of Math. (2), 2009, 170(1): 271-301.
  * [23] Gualtieri, M., Generalized complex geometry, Ann. of Math. (2), 2011, 174(1): 75-123.
  * [24] Henriques, A., Integrating $L_{\infty}$-algebras, Compos. Math., 2008, 144(4): 1017-1045.
  * [25] Iglesias-Ponte, D., Laurent-Gengoux, C. and Xu, P., Universal lifting theorem and quasi-Poisson groupoids, J. Eur. Math. Soc. (JEMS), 2012, 14(3): 681-731.
  * [26] Kelly, G. and Street, R., Review of the elements of $2$-categories, In: Category Seminar (Proc. Sem., Sydney, 1972/1973), Lecture Notes in Math., Vol. 420, Berlin: Springer-Verlag, 1974, 75-103.
  * [27] Khudaverdian, H. and Voronov, T., Higher Poisson brackets and differential forms, In: Geometric Methods in Physics, AIP Conf. Proc., Vol. 1079, Melville, NY: Amer. Inst. Phys., 2008, 203-215.
  * [28] Kinyon, M. and Weinstein, A., Leibniz algebras, Courant algebroids and multiplications on reductive homogeneous spaces, Amer. J. Math., 2001, 123(3): 525-550.
  * [29] Kravchenko, O., Strongly homotopy Lie bialgebras and Lie quasi-bialgebras, Lett. Math. Phys., 2007, 81(1): 19-40.
  * [30] Lada, T. and Markl, M., Strongly homotopy Lie algebras, Comm. Algebra, 1995, 23(6): 2147-2161.
  * [31] Lada, T. and Stasheff, J., Introduction to SH Lie algebras for physicists, Internat. J. Theoret. Phys., 1993, 32(7): 1087-1103.
  * [32] Lang, H. and Liu, Z., Crossed modules for Lie $2$-algebras, Appl. Categ. Structures, 2016, 24(1): 53-78.
  * [33] Lang, H., Liu, Z. and Sheng, Y., Integrations of derivations for Lie $2$-algebras, Transform. Groups, 2016, 21(1): 129-152.
  * [34] Lang, H., Liu, Z. and Sheng, Y., Affine structures on Lie groupoids, Pacific J. Math., 2020, 307: 353-382.
  * [35] Lang, H., Sheng, Y. and Xu, X., Strongly homotopy Lie algebras, homotopy Poisson manifolds and Courant algebroids, Lett. Math. Phys., 2017, 107(5): 861-885.
  * [36] Liu, Z., Sheng, Y. and Xu, X., The Pontryagin class for pre-Courant algebroids, J. Geom. Phys., 2016, 104: 148-162.
  * [37] Liu, Z., Sheng, Y. and Zhang, T., Deformations of Lie $2$-algebras, J. Geom. Phys., 2014, 86: 66-80.
  * [38] Liu, Z., Weinstein, A. and Xu, P., Manin triples for Lie bialgebroids, J. Differential Geom., 1997, 45(3): 547-574.
  * [39] Lu, J.H., Multiplicative and affine Poisson structures on Lie groups, Ph. D. Thesis, Berkeley, CA: University of California, Berkeley. 1990.
  * [40] Mackenzie, K., General Theory of Lie Groupoids and Lie Algebroids, London Mathematical Society Lecture Note Series, Vol. 213, Cambridge: Cambridge University Press, 2005.
  * [41] Mammadova, L. and Zambon, M., Lie $2$-algebra moment maps in multisymplectic geometry, Differential Geom. Appl., 2020, 70: 101631, 22 pp.
  * [42] Markl, M., A cohomology theory for $A(m)$-algebras and applications, J. Pure Appl. Algebra, 1992, 83(2): 141-175.
  * [43] Mill${\rm\grave{e}}$s, J., André-Quillen cohomology of algebras over an operad, Adv. Math., 2011, 226(6): 5120-5164.
  * [44] Noohi, B., Integrating morphisms of Lie $2$-algebras, Compos. Math., 2013, 149(2): 264-294.
  * [45] Ortiz, C. and Waldron, J., On the Lie 2-algebra of sections of an LA-groupoid, J. Geom. Phys., 2019, 145: 103474, 34 pp.
  * [46] Penkava, M., $L_{\infty}$-algebras and their cohomology, arXiv: q-alg/9512014.
  * [47] Rogers, C., $L_{\infty}$-algebras from multisymplectic geometry, Lett. Math. Phys. 2012, 100(1): 29-50.
  * [48] Rogers, C., $2$-plectic geometry, Courant algebroids, and categorified prequantization, J. Symplectic Geom., 2013, 11(1): 53-91.
  * [49] Rotenberg, D., Courant algebroids, derived brackets and even symplectic supermanifolds, Ph.D. Thesis, Berkeley, CA: University of California, Berkeley. 1999.
  * [50] Rotenberg, D., On the structure of graded symplectic manifolds and Courant algebroids, In: Quantization, Poisson Brackets and Beyond (Manchester, 2001), Contemp. Math., Vol. 315, Providence, RI: AMS, 2002, 169-185.
  * [51] Rotenberg, D., On weak Lie $2$-algebras, In: XXVI Workshop on Geometrical Methods in Physics, AIP Conf. Proc., Vol. 956, Melville, NY: Amer. Inst. Phys., 2007, 180-198.
  * [52] Rotenberg, D. and Weinstein, A., Courant algebroids and strongly homotopy Lie algebras, Lett. Math. Phys., 1998, 46(1): 81-93.
  * [53] Schlessinger, M. and Stasheff, J., The Lie algebra structure of tangent cohomology and deformation theory, J. Pure Appl. Algebra, 1985, 38(2-3): 313-322.
  * [54] $\rm{\check{S}}$evera, P. and Weinstein, A., Poisson geometry with a $3$-form background, Progr. Theoret. Phys. Suppl., 2001, 144: 145-154.
  * [55] Sheng, Y. and Liu, Z., From Leibniz algebras to Lie $2$-algebras, Algebr. Represent. Theory, 2016, 19(1): 1-5.
  * [56] Sheng, Y. and Zhu, C., Semidirect products of representations up to homotopy, Pacific J. Math., 2011, 249(1): 211-236.
  * [57] Sheng, Y. and Zhu, C., Integration of semidirect product Lie $2$-algebras, Int. J. Geom. Methods Mod. Phys., 2012, 9(5): 1250043, 31 pp.
  * [58] Sheng, Y. and Zhu, C., Integration of Lie $2$-algebras and their morphisms, Lett. Math. Phys., 2012, 102(2): 223-244.
  * [59] Voronov, T., Graded manifolds and Drinfeld doubles for Lie bialgebroids, In: Quantization, Poisson Brackets and Beyond (Manchester, 2001), Contemp. Math., Vol. 315, Providence, RI: AMS, 2002, 131-168.
  * [60] Wagemann, F., On Lie algebra crossed modules, Comm. Algebra, 2006, 34(5): 1699-1722.
  * [61] Waldorf, K., A global perspective to connections on a principal $2$-bundles, Forum. Math., 2018, 30(4): 809-843.
  * [62] Weinstein, A., Affine Poisson structure, Internat. J. Math., 1990, 1(3): 343-360.
  * [63] Weinstein, A., Omni-Lie algebras, S${\bar{u}}$rikaisekikenky${\bar{u}}$sho K${\bar{u}}$ky${\bar{u}}$roku, 2000, 1176: 95-102.