Nonlinear substructuring in the modal domain: numerical validation and experimental verification in presence of localized nonlinearities

In many systems of interest, most of the structure is well approximated as linear but some parts must be treated as nonlinear to get accurate response predictions: significant nonlinear effects are due to the connections between coupled subsystems, such as in automotive or aerospace structures. The present work aims at predicting the nonlinear behavior of coupled systems using a substructuring technique in the modal domain. This study focuses on the effects of nonlinear connections on the dynamics of an assembly in which the coupled subsystems can be considered as linear. Each connection is instead considered as a quasi-linear substructure with stiffness that is function of amplitude or energy. The iterative procedure used here is enhanced with respect to previous works by enforcing a better control of the total energy at each iteration allowing to obtain the solution for a prescribed set of energy levels. Also, the initial guess and the convergence criterion have been modified to speed up the procedure. This technique is applied to a system made of two continuous linear subsystems coupled by nonlinear connections. The numerical results of the coupling are first compared to the ones obtained by using the Harmonic Balance technique on the model of the complete assembly to evaluate its effectiveness and understand the effects of modal truncation. Besides, a nonlinear connecting element, specifically designed in order to have a nearly cubic hardening behavior, is used in an experimental setup. Substructuring results are compared to experimental results measured on the assembled system, in order to evaluate the correlation between mode shapes and the accuracy in the resonance frequency at several excitation levels.


Introduction
During the last decades, the use of Finite Element Methods spread in many engineering fields, allowing to achieve very accurate results. However, the main draw-back of these methods is that both the computational time and the reliability of the results are directly related to the spatial discretization chosen for the system [1]. Moreover, if transient phenomena are under evaluation and time integration of the equations is required, the correct time step needs to be carefully chosen in order to achieve the needed information [2]. This is true, for instance, for vibrating mechanical systems that must satisfy many types of requirements. In the early phases of the design process, quick results are needed to evaluate the performances of the considered system, thus FE methods might not be the proper choice to make.
Dynamic Substructuring techniques represent a valid alternative to analyze large engineering systems, such as aircrafts or cars, providing accurate data in a short time. They were introduced during the 1960s by Hurty [3], Craig and Bampton [4] and Rubin [5] and have become widely used because they allow to analyze large numerical linear systems. On the other hand, also experimental approaches to the substructuring problem using frequency response functions were developed [6,7]. In 2008 De Klerk et al. gave an extensive review of the substructuring techniques developed so far [8], and proposed a general framework to formulate substructuring problems either in physical, modal and frequency domains. These methods have become very popular since their range of application is very broad: it is possible to obtain the dynamic behavior of complex systems starting from the known dynamic behavior of its component substructures (coupling) [9,10]; on the other hand the behavior of one substructure can be achieved from the known behavior of the complete structure and that of the residual substructures (decoupling) [11][12][13]. It is the case, for example, of the decoupling technique developed by Saeed et al. [14] used to identify the linear characteristics of a joint in an assembled structure. Besides, coupling and decoupling can be performed by combining experimental models of some substructures and numerical models of the other substructures [15][16][17], as done when using the Transmission Simulator approach [18].
Although these techniques have been constantly improved to solve increasingly challenging problems, they are all based on the assumption that the underlying systems are linear. This is a strong limitation, especially if it is known that the considered components do not behave linearly in the working range in which their behavior is analyzed. For this reason, the presence of nonlinearities has to be accounted for if more accurate results are needed.
Some attempts have been made during the last years to allow dynamic substructuring techniques to deal also with nonlinear problems. Mahdiabadi et al. [19] and Kuether et al. [20] propose a non intrusive method to deal with geometrical nonlinearities of nonlinear systems to obtain a Nonlinear Reduced Order Model (NLROM), using the Implicit Condensation and Expansion (ICE) technique proposed by Hollkamp [21]. This method allows to express the nonlinear terms as polynomials whose coefficients are estimated using lest squares on force-displacement data obtained from a full finite element model. Two different approaches, named "construction" and "extraction" are instead suggested by Chong and Imregun [22] to model the behavior of nonlinear substructures. Both of them are used to obtain the modal parameters describing the system to be used for the coupling in the modal domain. The construction technique is used to obtain nonlinear modal parameter variations from a physical nonlinear stiffness or damping laws while the extraction technique allows to obtain the nonlinear modal parameters via a nonlinear modal analysis of measured response data [23]. The whole procedure using the extraction technique is then applied to deal with an experimental setup composed of a beam with cubic nonlinearity [24]. Another way to account for the presence of nonlinearities in substructuring analyses is by using Describing Functions (DFs) [25]. Tanrikulu et al. [26] first used DFs developing an approach to account for the effects of different types of nonlinearities on the dynamics of a nonlinear system. Wei and Zheng [27] used an approach based on DFs to reduce the nonlinear problem into a set of complex algebraic equations to be iteratively solved such that it is possible to tackle those cases in which the nonlinearity is localized [28]. The use of the DFs is fundamental in the works of Kalaycioglu and Özgüven, who developed a frequency-based approach to account for the presence of nonlinearities in coupling and decoupling procedures. In particular, in [29,30] the coupling problem between multiple substructures with localized nonlinearities is addressed, while in [31] the decoupling one is investigated. In the former case the structural modification technique developed by Özgüven [32] is modified and the presence of the nonlinearity is accounted for in a nonlinear matrix using DFs that is updated at each iteration. In the latter case different formulations are used depending on where the nonlinear element is placed: if in the known subsystem, a modal identification technique [33] is used; if it is in the unknown subsystem the nonlinearity matrix [26] is used; if it is at the connection between them, it can be reduced to one of the previous cases, depending on whether the properties of the nonlinear element are known or not.
Wang et al. [34] focused the attention on the decoupling between components in which the nonlinearity could be addressed as a nonlinear element between two coordinates. In particular, they managed to obtain the Frequency Response Function (FRF) of the nonlinear substructure assuming to know the system-level FRFs of the coupled nonlinear system. It can be very helpful to perform the identification of nonlinear components in a complex nonlinear system. Inspired by the nonlinear coupling method in the frequency domain, Kuether and Allen [35] developed an iterative modal substructuring technique to account for the presence of nonlinearities, modeling each subsystem using a set of Nonlinear Normal Modes (NNMs). Nonlinear Normal Modes were first introduced by Rosenberg during 1960s [36][37][38], studied by Vakakis [39] and later by Kerschen et al. [40], who defined a nonlinear normal mode as the non-necessarily synchronous periodic motion of a conservative nonlinear system which is characterized by energy dependent frequency and amplitudes. Over the years, NNMs have been widely used to better describe the dynamic behavior of systems whenever the linear approximation was not sufficient [41,42]. As a result, many methods have been established to calculate the NNMs of a structure: the shooting and pseudo-arclength continuation [43], the Harmonic Balance (HB) [44,45] and the Modal Derivatives (MDs) [46,47]. Peeters et al. in [48] propose some ways to experimentally measure NNMs of nonlinear structures.
The present work aims at predicting the nonlinear behavior of coupled systems using substructuring techniques. In particular, the focus of this study is to analyze the effects that nonlinear connections have on the dynamics of an assembly in which the coupled subsystems can be considered as linear. This assumption is valid in many engineering situations in which the connecting element introduces a nonlinearity much more relevant than those of the substructures to be coupled, as it happens with bolted joints [49], wire rope isolators [50] or for the connections of turbine blades to a rotor [51]. In this context it is then possible to isolate the nonlinearity and consider it as a substructure itself to be included in the coupling process. The method proposed here is inspired by the work of Kuether and Allen [35] who define an energetic approach to perform the nonlinear coupling in the modal domain to obtain the NNMs of the complete assembly. Due to the assumption of considering only the connecting elements as nonlinear, in the present case the NNMs have to be evaluated only for the nonlinear connecting elements, whereas the linear substructures are described through a possibly truncated set of Linear Normal Modes (LNMs). Note that modal reduction can be also performed to reduce the number of DoFs included in LNMs. These assumptions lead to a reduction of the computational burden, especially if large systems are involved. The NNMs obtained from the coupling are then compared to the ones calculated applying the HB method on the whole system. This method was already used to couple lumped parameter models connected by single or multiple nonlinear springs, presenting hardening and softening connections [52,53].
In the present work the method has been significantly enhanced with respect to the previous work by enforcing a better control of the total energy at each step of the procedure through an adjustment of the energy distribution among substructures. Also, the initial guess and the convergence criterion have been modified in order to speed up the procedure. Furthermore, the method is here applied to a more realistic system made up of two continuous linear subsystems coupled by nonlinear connections. The numerical results of the coupling are first compared to the ones obtained by using the HB on the model of the complete assembly to evaluate its effectiveness and understand the effects of modal truncation. Besides, a nonlinear connecting element has been designed and manufactured in order to have a nearly cubic hardening behavior. An experimental setup, composed by two beams joined by two of the aforementioned connecting elements, was used. Tests have been performed on both the beams and the nonlinear connecting elements to identify the parameters of the corresponding numerical models, to be used in the coupling procedure. Substructuring results are compared to experimental results measured on the assembled system, in order to evaluate the correlation between mode shapes and the accuracy in the resonance frequency at several excitation levels.
The present method is based on the same theoretical assumption made by Kalaycioglu and Özgüven in [30,31] of considering only the fundamental frequency of the response of each mode. This assumption does not allow to account for higher harmonics effects and internal resonances. However, it allows to investigate the dynamics of many engineering systems at the transition between linear and nonlinear behavior. The proposed method is developed using a modal approach whereas the method in [30] is set in the frequency domain. Using the frequency-based approach, experimental data can be easily used to perform the analysis and retrieve the variation of the resonance frequency at different excitation levels. However, a dense discretization in frequency is necessary to correctly estimate the resonance frequency, thus possibly leading to an increase in the computational time, especially if a system with many DoFs is considered. Using the modal approach, instead, experimental data might not be straightforward to obtain [54], but the method itself does not require heavy computation and the convergence is reached with reduced computational burden.
The advantages of using nonlinear substructuring techniques, such as the one presented in this paper, with respect to other methods that require the model of the entire system, such as the HB method, shooting and pseudo-arclength continuation, FEM, are implicit in the substructuring approach: combining experimental and numerical models, considering simpler subsystems instead of a unique complex one, performing indirect analyses and allowing a reduction of the computational burden.
The work is organized as follows: Section 2 highlights the main characteristic of substructuring and the steps of the enhanced nonlinear coupling procedure are presented. Section 3 introduces the models used in the two case studies: the fully-clamped system and the free-free system, the latter being inspired to the experimental setup. Section 4 presents the results of the simulation on the fully-clamped system, highlighting the effects of modal truncation. Section 5 shows the comparison between the numerical results and the preliminary experimental results on the free-free system.

Theoretical background
The present approach extends the classical linear substructuring techniques to include the effects of nonlinear phenomena according to NNMs theory. The linear substructuring technique in the modal domain and the nonlinear coupling procedure are introduced in the following section.

Linear coupling procedure in the modal domain
Substructuring methods are used to analyze complex structures composed of multiple parts either to evaluate the behavior of an unknown subsystem from the one of the complete assembly (decoupling) [12] or to obtain the behavior of the assembly knowing the single subcomponent ones (coupling) [8], as it is done here. The dynamic equation for an undamped substructure is: where M s) is the mass matrix, K (s) is the stiffness matrix, f (s) is the vector of external forces and g (s) is the vector of connecting forces with the adjacent substructures. The equations of motion of each subsystem can be collected in a block diagonal form as: where x = {x (1)T · · · x (s)T · · · } T , so that M and K are, respectively, the mass and stiffness block diagonal matrices, f collects the external forces and g collects the connecting forces. The coupling between the substructures is performed by imposing the compatibility and equilibrium conditions. The first one implies that any pair of matching DoFs, e.g. DoF a on subsystem s and DoF b on subsystem t must have the same displacement, that is x This condition is in general expressed as: where B is a signed boolean matrix. The equilibrium condition is referred to internal constraint forces and implies that the sum of connecting forces at a pair of matching DoFs must be zero, e.g. g (s) a + g (t) b = 0: this holds for any pair of matching DoFs. Furthermore, if DoF l on subsystem r is not a connecting DoF, it must be g (r ) l = 0: this is valid for any non-connecting DoF. Overall, the above conditions can be expressed as: The compatibility condition is implicitly satisfied if a unique set of interface DoFs q is defined as: where L is the same matrix introduced earlier. Combining Eqs. (3) and (5): then BL = 0, i.e. L = null(B). Considering a large number of degrees of freedom can result in an increase of the computational cost for the method. Thus, modal reduction is in general performed by considering only a truncated set of the free interface modes of each substructure. However, if no truncation takes place, the number of degrees of freedom does not change and it becomes a simple transformation into modal coordinates that might be convenient to solve some cases. Modal reduction is performed by building the block diagonal matrix R, where each block contains the set of modes of the corresponding substructure. Thus, it is possible to write: then, by substituting Eqs. (7) and premultiplying by R T , Eq. (2) becomes: where The compatibility condition can be expressed in modal coordinates as: According to Eq. (10) the compatibility is enforced between the modes of the various substructures and it is imposed exactly if no modal truncation is performed, otherwise there is some constraint relaxation.
A unique set of generalized modal coordinates z satisfying the compatibility condition can be defined using the same approach followed for the physical coordinates: Satisfying the compatibility condition implies in both cases that: and so: Hence, the equilibrium condition can be written as: By substituting Eq. (12) and premultiplying by L T m Eq. (8) becomes: with: As it will be outlined in 2.2, Eq. (16) will be used to couple nonlinear substructures too.

Nonlinear coupling procedure in the modal domain
Nonlinear coupling [30,31] and nonlinear modal substructuring [35,55] account for the effects of the nonlinear internal forces considering them as harmonic at the same frequency of the excitation, thus higher harmonics effects are not captured. The same assumption is valid for the procedure discussed in this section, that is meant to couple linear substructures connected through nonlinear connecting elements.
The procedure is applied to analyze the dynamics of systems like the one shown in Fig. 1, where each linear substructure L i has n i DoFs and each nonlinear substructure has 2 connecting DoFs. It follows that the The coupling in the modal domain, as seen in Sect. 2.1, requires the modal bases for the considered substructures. In this case, LNMs are used for the linear subsystems and NNMs for the nonlinear ones. The advantage of the Nonlinear Coupling Procedure relies on the fact that instead of computing the NNMs of the whole nonlinear structure, thus requiring a lot of time, it suffices to compute the NNMs only for the simpler nonlinear substructure, for example using the HB method. The coupling procedure is then used to propagate the local nonlinear behavior on the dynamics of the complete assembly, thus allowing to obtain its NNMs. Significant improvements have been made with respect to the method used in [52,53] in order to control the energy levels at which the NNM is evaluated. These lead also to a relevant enhancement in the convergence speed of the numerical procedure. Since the NNMs are energy dependent [40], the coupling is performed by considering increasing energy levels collected in the array where e k indicates the k-th energy level at which the coupling is performed. The idea is to reconstruct the NNM curve for each mode of the coupled system by using an iterative continuation technique, in which the solutions at the two previous energy levels are used to derive the initial guess for the following one, as shown in Fig. 2 and detailed in Sect. 2.2.3. In order to start the procedure, two initial points of the NNM of the complete structure are needed, as shown in Fig. 3. These points can be obtained by choosing the first two energy levels of the array e to be very low. In this way the first two points on the NNM are  computed by using a linear approximation for the initial guess (detailed in Sect. 2.2.2). Since all the subsystems are considered as undamped, the energy associated with each substructure is constant during the free oscillation and equal to the sum of the elastic and kinetic energy. The energy at each iteration is computed by considering that the eigenvector coordinates represent the amplitude of the displacement oscillation of each degree of freedom and that the velocity can be derived since an harmonic oscillation is being considered.
The complete iterative algorithm of the Nonlinear Coupling Procedure is shown schematically in Fig. 4, where: the blue area indicates the operations of the nonlinear coupling procedure (Sect. 2.2.1); the yellow area indicates how the initial guess is obtained for the first two energy levels based on the linear approximation (Sect. 2.2.2); the green area indicates the predic-

Nonlinear coupling in modal domain
The coupling in the modal domain requires the definition of the matrix R introduced in Eq. (7). Since nonlinear substructures are involved, the reduction matrix includes both LNMs and NNMs. This part of the algorithm aims at iteratively obtaining the energy distribution of the substructures in the coupled configuration, i.e. the correct energy associated with the NNMs included in the matrix R, and at comparing it to the distribution found at the previous iteration, until convergence is achieved.
Once the distribution of the energy e k among the substructures is estimated, it is possible to perform the coupling in the modal domain. The equations of motion describing the nonlinear system are and f nl (x) is the vector representing the nonlinear term. It is now possible to perform the coupling procedure according to the following steps: Step 1: normalization of the nonlinear modes according to the energy associated with the correspondent substructure; Step 2: coupling of the linear and nonlinear substructures; Step 3: convergence check based on an energy criterion that will be discussed in the sequel of this Section.
The reduction matrix introduced in Eq. (7) becomes energy dependent and can be written for the k-th energy level for the assembly and iteration i using the matrices of linear eigenvectors A and B and mode shapes C E of the nonlinear modes computed at energy E C E k,i obtained from Eq. (42) for the case where all the nonlinearity is in the connection: The eigenvectors of linear substructures are mass normalized and the t-th nonlinear mode ψ t in C E is normalized accordingly (Step 1): being χ t (E C E k,i ) the mode shape taken as the initial displacement of the substructure for the periodic solution at a given energy level E C E k,i . Note that in the present approach the nonlinear connection is considered as a substructure, thus the mass matrix must be necessarily defined and the NNMs can be always computed. This is not a limitation since if the connecting element is a mass-less nonlinear spring, end masses can be added by subtracting them from the linear substructures [52]. The denominator in Eq. (22) is a normalization coefficient and it will be addressed as c s t since it is computed for the t-th mode of the s-th substructure, and all can be gathered to form c k : that is computed only at the first iteration of the kth energy level and it will be used during the whole iterative procedure to perform the transformation from modal to physical domain. Similarly the matrix of the corresponding "eigenvalues" becomes: where A and B are the eigenvalues of the linear substructures, whereas C E (E C E k,i ) is the matrix of the squared angular frequencies for the NNMs of the nonlinear substructure at a given energy level E C E k,i , computed using the HB method.
The matrix R k,i in Eq. (21) can be used to define the nonlinear modal transformation: that holds approximately in case of modal truncation. By forcing the assumption of orthogonality of the NNMs with respect to the mass and stiffness matrices [35], the equation of motion in Eq. (19) becomes: where R k,i accounts for the nonlinear force present in Eq. (19). Hence Eq. (12) becomes: whereL k,i depends on the energy level. It is possible to perform the coupling (Step 2) following the same approach outlined in Sect. 2.1 from Eq. (8) to Eq. (16): where the matricesM k,i andK k,i are updated at each iteration according to Eq. (17). The system in (28) is a linear approximation of the nonlinear system at the energy e k , so it is possible to compute its eigenvalues and eigenvectors. A given eigenvalue-eigenvector pair represents the natural angular frequency ω r and the eigenvector ξ r expressed in generalized modal coordinates z of the r-th mode for the coupled system. This result can be brought back in the physical domain by performing a coordinate transformation according to Eq. (25) and (27): wherec is the participation factor obtained as being () + the pseudo-inverse.
Since the system is nonlinear, the predicted vector χ k,i is associated with an energy level that may be different from e k . Thus, it is necessary to scale the predicted vector χ k,i by introducing a scaling coefficient α, obtaininĝ such that, by defining x k,i (t) = χ k,i e jω k,i t , the energy level associated withχ k,i is forced to be equal to e k introduced in Eq. (18), i.e.
Although Eq. (32) holds at any instant, it is useful to evaluate it when the amplitude of oscillation is maximum, i.e. for t = 0: For a cubic spring, f nl (αχ k,i ) = k nl α 3 δ 3 k,i , where δ represents the deformation of the nonlinear spring. Therefore, Eq. (32) becomes: This represents one of the main differences with the approach outlined in [35,52], because in that case the result from Eq. (29) would be directly used to compute the energy of the substructures in the coupled configuration. This new approach, instead, guarantees that the total energy of the assembly is kept constant throughout the iterative procedure. At this point it is possible to evaluate the energyÊ k,i at iteration i of the substructures in the coupled configuration. The iterative algorithm stops when a stopping criterion is satisfied. While the convergence crite-rion in [35,52] is based on the comparison between the energy distribution of the substructures in the coupled and uncoupled configuration, here it is defined as: where the norm of the energy increment at iteration i is related to the norm of the energy at iteration i − 1.
In this way, when the difference between the energy distribution obtained at two consecutive iterations is below a certain tolerance, the algorithm stops. Note that if i = 1,Ê k,0 is assumed to be initial guess vector for the energy distribution E k,1 . If Eq. (35) is not satisfied, the initial energy for Eqs. (21)-(30) is updated as: where λ is a weighting factor ranging from 0 to 0.5 to improve the numerical stability of the iterative procedure. This updating expression is different from the one outlined in [35,52], where the convergence criterion is developed as a combination of the energy distributions of the substructures in the uncoupled and coupled configuration. In Eq. (36), instead, the energy distributions are always referred to the coupled configuration. In fact, since the initial guess is specifically chosen to be close to the real solution, in this correction the initial estimate of the energy distribution E k,i is given more importance with a weighting factor greater than 0.5, and the other termÊ k,i is used to perturb it for the next iteration. If Eq. (35) is satisfied, the iterative procedure is stopped, the solution for the energy level k is stored and the next energy level k + 1 can be discussed. Also, in order to enforce that the energy is always equal to e k , the following normalization is performed to obtain the input energy distributionẼ k for iteration i + 1: This represents also an improvement with respect to the method in [35,52] because it enforces a better control of the energy at each iteration. It is important remarking that the passages from Eq.(31) to Eq. (37) have been modified with respect to [35,52] such that the energy at each iteration is kept constant and equal to the desired value e k . Also, the convergence is verified by comparing the energy distribution of the substructures between two consecutive iterations.

Estimation of the initial guess, k ≤ 2
This procedure is used to initialize the computation for the r-th mode and it is performed at the first iteration, only for the first two energy levels. It is convenient to set their correspondent energy value to be sufficiently low such that the LNM and NNM of the complete system almost coincide in frequency and mode shape. In this way the coupled system can be considered as linear and its eigenvectors are easily computed. Referring to Fig. 1, where K tot and M tot are the stiffness and mass matrices of the linear assembled system defined on the unique set of DoFs q, and V lin and lin are the eigenvectors and eigenvalues matrices. The considered r-th mode is associated with an energy E r equal to Since the linear eigenvector v r needs to approximate the nonlinear NNM at a given low energy level, it is necessary to normalize it such that its energy is equal The index r is now dropped since the solution is found for one mode of the coupled system by varying the energy in all of the component modes [35]. Instead, the index referring to the iteration is introduced, which for the first one is i = 1. Using Eq. (5) it is possible to expand the eigenvector v in the primal formulation to the eigenvectorū in the dual domain, thus obtaining in the case of two linear subsystems connected through one nonlinear connection: whereū A ,ū B andū C E represent the portions of the eigenvectorū corresponding to the substructures A, B and CE, respectively. With this partition it is possible to compute how the energy e k is distributed among the substructures: where the energy e k is the sum of the components E j k,1 of the vector and the matrices K A , K B represent the linear stiffness of each substructure and K C E represents the tangent stiffness matrix evaluated at the static equilibrium position. In doing so, the initial guess describing how the energy e k for the energy levels k = 1, 2 is distributed among the substructures is obtained. This approach can easily be extended to the case with multiple nonlinear connecting substructures.

Estimation of the initial guess, k > 2
One of the novelties introduced in the procedure with respect to the one proposed by Kuether and Allen [35] is in the initial guess of the energy distribution for the iterative algorithm. In fact, it is estimated by using the energy distribution found at the two previous energy levels instead of considering the energy as equally distributed among the substructures. The initial energy distribution associated with the first two energy levels is found using the linear eigenvectors for all the substructures. For higher energy levels (k > 2) the initial energy distributionÊ k and the mode shape of the NNM of the complete structureχ k can be extrapolated using a linear projection of the two previous energy levels. They are used as the predictions for the first iteration because it is reasonable to suppose that the solutions between two consecutive energy levels do not differ much from each other if considering that the increase of energy is low. This approach allows to account for the trend of variation of each component of either the initial energy distribution and the mode shape of the NNM. Thus, considering for example the h-th element ofχ k,1 and -th element of E k,1 : In doing so, it is not necessary to perform the coupling for the first iteration, thus it is possible to directly go to Eq. (35). This projection is performed only for the first iteration when k > 2, while for the following iterations the procedure shown in Sect. 2.2.1 holds.

Models
In a previous work the Nonlinear Coupling Procedure (NLCP) was applied to lumped parameter models [52].
Here, instead, it is used to couple two continuous systems through nonlinear connecting elements. Two different cases are considered, the first considering a fully clamped system and the second being a free-free system. In particular, the second model reproduces an experimental system that was built in the laboratory of the Department of Engineering Physics at University of Wisconsin-Madison.

Fully clamped system
The beams are modeled using 20 Euler-Bernoulli finite elements of equal length, as shown in Fig. 5, considering a 2D problem. The axial DoF of the beams can be neglected, therefore each node has two degrees of freedom, the one related to the vertical displacement and the one accounting for the rotation. The beams are connected together through two cubic springs, as shown in Fig. 5, and their dimensions and material properties are summarized in Table 1. In order to have two different displacement fields for the two beams at a given frequency, the lower beam's thickness is twice the upper one's. In this way there is no occurrence in which the nonlinear springs are not subjected to deformation. Besides, the nonlinear elements are not symmetric with respect to the center where k l = 1000 N/m and k nl = 5 × 10 6 N/m 3 that are chosen such that the stiffness of the springs has the same order of magnitude of the static stiffness of the beams.

Free-free experimental system
It is composed of two beams of 910 mm x 30 mm x 6,35 mm made up of steel, connected together through two Nonlinear Connecting Elements (NLCEs), as shown in Fig. 6. Since is hung from the frame through elastic cords, the system can be considered as free-free. The crucial element in this setup is the NLCE, shown in Fig. 7a. The beams' thicknesses are such that they would yield before exhibiting nonlinearity whilst the NLCE has been designed such that it exhibits significant nonlinearity before yielding. The NLCE is composed of a thin plate of spring steel connected at its ends to a C-shaped element and, at the middle, to a T-shaped element, as shown in Fig. 7.
The NLCE dimensions are listed in Table 2.
The configuration of the nonlinear element is inspired to the cubic behavior of a beam clamped at both ends and subjected to large displacement at mid-span. In this case, the plate plays the role of the clamped beam. Every component of the assembled structure is tested to obtain the necessary information to build a numerical model that fit the experimental data. Each one of  Fig. 8. Every node has two degrees of freedom, the one related to the vertical displacement and the one accounting for the rotation. The NLCEs are modeled as a nonlinear cubic springs (orange lines) and are located in accordance with their position in the real system.

Nonlinear coupling using complete modal basis
The NLCP and the HB method with one harmonic are performed to compute the first five NNMs of the coupled system. The results can be compared in a Frequency-Energy Plot (FEP), as shown in Fig. 9. The modal basis used to reduce each beam contains all the modes derived from the discretization (i.e. 38 modes per beam, considering that both beams are clamped at both ends). The NLCP curves manage to retrace the HB main backbone in a wide energy range,  even in the areas where there is a significant increase in the frequency, as it is for mode 5. As it is shown in Fig. 10, the deviation in frequency associated with the results is below 0.2% and its average is about 0.05%. Note that the NLCP results well retrace the main backbone of the NNMs computed using the HB method with 3 harmonics, as shown in Fig. 11.

Enhancements of the NLCP
The NLCP technique outlined in Sect. 2.2 has been significantly enhanced with respect to the one described in [52]. It is possible to evaluate the effects of these improvements by applying the two methods to the fully-clamped system and compare the results. They share the same theoretical approach, thus the NNMs obtained using either one of them are the same. However, as already explained, major improvements have been made to have a better control on the total energy of the coupled system and to increase the rate of convergence. For both the procedures it is necessary to define the array of energy levels at which the coupling is performed. The chosen array is composed of 400 logarithmically spaced points between 10 −4 and 10 4 (roughly 50 points per order of magnitude). Figure 12 shows the energy of the five NNMs of the coupled structure computed using the coupling procedure with both the methods.
The results of the old coupling procedure show that the energy of the obtained NNMs is always different from the energy levels given as input. Furthermore, the set of energy levels at which each NNM is defined is different from the others. Figure 13 shows, for example, the FEP for the fifth mode obtained using both the methods.
Even if the first energy level given as input is equal to 10 −4 J , the old NLCP converges at ≈ 10 −1 J , thus the solution for lower energy levels cannot be achieved. Also, there are wide energy intervals in which the solution is not available. In engineering applications both circumstances are a limit, especially if the lack of information affects those energy levels at which the system works. Using the enhanced NLCP, instead, it is possible to get the solution at the required energy levels by just inserting them in the energy array e defined in Eq. (18). Conversely, the enhanced NLCP provides the NNMs at the same energy levels used as input. Figure 14 shows for the second mode, as an example, the cumulative  The overall number of iterations performed using the enhanced NLCP (solid line) is much lower than those required for the old coupling procedure (dashed line) using the same number of energy levels. Consequently, the computational time employed for the analysis of each mode decreases if using the enhanced NLCP, as summarized in Fig. 15. The enhanced NLCP performs the analysis in less than 3 minutes, while the other requires about 30 minutes to compute the five NNMs. This comparison between the two methods highlights the importance of the performed improvements: on one hand, the computational time decreases with respect to the previous case, and on the other hand, it is possible to obtain the results at the desired energy levels. The latter statement might even be considered of major importance, because the two consecutive solutions are at pre-set energy levels, thus there is the possibility to obtain the NNM for each required energy value. Using the enhanced NLCP, the NNM can be obtained at any given energy level, whereas the use of the previous procedure would most likely require the interpolation of the results to get the NNM at a specific energy value.

Modal truncation effects on nonlinear coupling accuracy
It is interesting to evaluate how modal truncation affects the accuracy of the substructuring procedure. Six truncation frequencies are chosen and the coupling is performed considering the proper number of modes for each beam with respect to each truncation frequency, as summarized in Table 3. The curves describing how the error in frequency varies for each mode at each truncation frequency are shown in Fig. 16. As it can be seen from the figures, the error depends either on the truncation frequency and also on the energy level at which the curve is computed. The error decreases as the truncation frequency increases, and the opposite happens as the energy grows. Considering for example the first mode, whose linear frequency is 12 Hz, it is possible to see that using a truncation frequency of 100 Hz (yellow line) the error is about 0.5% at 0.001 J, but it reaches almost 5% above 1000 J. In this case, the truncation frequency is almost ten times the linear frequency, and yet the error grows of one order of magnitude. This holds also for the second mode, when a truncation frequency of 50 Hz is considered (orange curve), resulting in a error that starts below 0.5% and overcomes 100%. The curves in Fig. 16 highlight that this is a general trend for all the modes, especially if the used truncation frequency is low with respect to the frequency of the linear mode under evaluation. However, it is not possible to simply relate the characteristic error for a mode with the choice of the truncation frequency because it is important to take also into account the energy at which the dynamics of the system needs to be evaluated. These curves show how the modal truncation affects the error in frequency, but they do not provide any information regarding the effects on the nonlinear mode shapes of the system. The correlation between the mode shapes can be estimated by selecting an energy level and a truncation frequency and computing the Modal Assurance Criterion (MAC) [56] of the modes obtained with the NLCP (shortened as CP in the following formula) considering the modes computed with HB as reference: The MAC is evaluated at low and high energy levels for the first five modes. Referring to Fig. 9, the low energy level for each mode corresponds to a point on the curve at which the hardening effect is low (left end of the curve). The high energy level for each mode, instead, corresponds to a point at which there is a significant hardening effect (right end of the curve). The selected energy values for each mode are listed in Table 4. Figure 17 shows different evaluations of the Modal Assurance Criterion for the five modes, computed for each truncation frequency at low and high energy.
In particular, blue bars in the graphs refer to the MAC of the assembly and show how the decrease in the truncation frequency leads to a reduction of the MAC number, especially for modes 1 and 2 at both low and high energy. For the considered five nonlinear modes, 200 Hz (i.e. ten retained modes) seems to be the lowest possible value for the truncation frequency such that the corresponding mode shapes are not affected much by the modal truncation in the considered energy range.
Also, the partial Modal Assurance Criterion (pMAC), defined in [57], is computed for the upper and lower beams and for the nonlinear connections in order to correlate parts of the modal vectors (orange, yellow and purple bars, respectively, in the histograms of Fig. 17). A low value of the pMAC indicates the portion of modal vector (substructure) in which there is poor correlation between the mode shapes. As it happens for the MAC, also the pMACs are affected by the number of modes retained during the reduction process. It is interesting to highlight that the MAC can be lower than the pMACs: for instance, in the present case it happens for mode 1 at high energy, when it is computed using a modal basis corresponding to a truncation frequency of 25 and 50 Hz. This means that the poor correlation indicated         by the MAC does not depend only on the partial correlation (i.e. local mode shapes) between the involved parts, but also on the energy distribution among the substructures (i.e. the local amplitudes of the mode shape). Figures 18a and b show the mode shapes for the second and fifth NNMs of the whole structure evaluated at the corresponding high energy values for the considered truncation frequencies. It is possible to relate the graphical representation of the mode shapes with the results of the pMACs. For example, considering the shape of the fifth mode computed for the truncation frequency at 100 Hz (green curve in Fig. 18b), the local shape of the upper beam is similar to the shape computed with the HB but with a different amplitude, while the local shape of the lower beam is not. This agrees with the corresponding results of the pMACs shown in Fig. 17. Besides, for the second mode shown in Fig. 18a, at low values of the truncation frequency, (25 and 50 Hz) there is a significant difference in the local shape of the upper beam, highlighted also by the corresponding pMAC. On the other hand, for the lower beam there is a good correlation in shape but a significant difference   in amplitude that is not detected by the corresponding pMAC but only by the MAC. The energy distribution defined in Eq. (42) provides a useful insight about the error in the amplitude of the local shapes for each portion of the assembly. Figure 18c and d shows that as the truncation frequency decreases, the energy associated with each substructure differs from that computed using the HB (dashed lines), thus explaining the difference in amplitude for the local mode shapes.
These results indicate that, in case of modal truncation, the number of modes in the modal basis affects the dynamics of oscillation (frequency and amplitude) depending on the considered energy level. Also, the results show a substantial convergence of the results for an increasing number of modes included in the modal basis. Thus, the CP can be considered as reliable to obtain meaningful information regarding the dynamic behavior of the system.

Nonlinear coupling procedure using experimental data
In this Section, the Coupling Procedure is used to predict the behavior of the experimental system introduced in Sect. 3.2. The nonlinear coupling procedure, and more in general the experimental substructuring techniques, requires data that can be difficult to measure. For instance, rotational coupled DoFs are very tricky to obtain experimentally, both in terms of displacements (rotation) and forces (momentum). Thus, it might be useful to deal with numerical models of the substructures, properly adjusted by using the experimental data of the real subsystems. In the present case, the assembly is composed of two beams and two NLCEs and in the two following sections, the identification process is outlined for each of these components. The last section, instead, provides a comparison between the nonlinear  Fig. 19 Comparison between the numerical and experimental FRF at the driving point of one beam experimental response of the free-free system and the results of the nonlinear coupling procedure performed between the updated beams and the identified NLCEs.

Experimental tuning of beam models
A linear modal analysis test is performed on both beams separately to estimate their resonance frequencies. The comparison between the numerical and experimental FRF of one beam is shown in Fig. 19.
Four different frequencies are distinguishable in the range from 0 to 500 Hz. These results are used to identify the correct value for the Young Modulus such that the Euler-Bernoulli beam model fits the experimental data.

Experimental identification of the nonlinear connecting element
Since each NLCE is regarded as a substructure during the substructuring procedure, it is necessary to determine their NNMs. Since the NLCEs are realized in laboratory, their NNMs can be obtained experimentally. After having measured the NNMs of interest, they are used to fit the coefficients of the cubic law that better approximates the NLCEs' behavior. Here the whole procedure is briefly explained for the NLCE addressed as "left spring" (referring to the coupled configuration in Fig. 6); the same procedure is valid also for the right element. A burst random analysis is performed using a laser scanner vibrometer to obtain a frequency spectrum used to estimate the linearized frequency corresponding to the mode of interest. The set of measurement nodes is shown in Fig. 20, where nodes 1-20 are referred to the T-shaped element, nodes 21-30 to the thin plate and nodes 31-40 to the C-shaped element.
Since the NNMs are frequency-energy dependent, the measurement of the NNM is performed by setting the amplitude and by adjusting the frequency of the sinusoidal excitation such that the resonance condition is attained: this occurs when the velocity measured by the scanner is in-phase with respect to the input signal. This procedure is performed for each energy level by step-wise increasing the amplitude of the excitation, for a total of seven energy levels for each NLCE. For each amplitude-frequency pair of the NNM the velocity of oscillation at the resonance condition of each measurement DoF is recorded. The velocity signal is integrated to obtain the displacement time series related to the nodes of the C and T shaped element, after removing the low frequency contents due to rigid body oscillation in the free-free configuration. The measured NNM is that of a translational spring, thus, in order to represent it, the amplitude of the relative displacement between the C and T shaped elements needs to be associated with the corresponding resonance frequency. For each node belonging to the C-shaped element, there are two nodes on the T-shaped one sharing the same value of abscissa: referring to Fig. 21  (b) NNM of the right NLCE. Step 1: Compute the average instantaneous displacement of the nodes belonging to the T-shaped element that share the same abscissa (for example, (u 1 + u 11 )/2); Step 2: Compute the absolute value of the difference between the displacement of the node belonging to the C-shaped element and the result obtained from Step 1; Step 3: Compute the amplitude of oscillation of the result obtained from Step 2 averaged over the recorded periods.
Each of these steps is performed for every group of 3 nodes and for each of the seven levels of excitation. As a result there are ten curves, each one containing seven points corresponding to the measured data. The NNM is obtained as a the average of this curves, as shown in Figs. 22a and b for the left and right NLCE, respectively.
As it is evident from the two curves, both the NLCEs experience an hardening behavior. These sets of data are used as reference to identify the parameters k l and k nl for the two cubic springs performing a minimization through a least-square error method. Each NLCE is considered as a 2 DoFs oscillator, as shown in Fig. 23, with masses m c = 0.297 kg and m t = 0.232 kg. Thus, it is possible to use the numerical law describing the frequency of oscillation as a function of the deformation of the spring [54,58] in the minimization. In this way the NNM computed through the HB described by that law is the same as the one measured. The coeffi-  cients resulting from the minimization are gathered in Table 5 and in Fig. 24a and b the comparison between the experimental NNM and the numerical one is shown for both the left and right NLCE.

Experimental nonlinear response of the free-free system
The dynamic response of the free-free system, described in Sect. 3.2, is measured in order to find its experimental NNMs. The comparison with the experimental data is performed with the results obtained using the numerical models defined in the previous sections. Finally, a new set of parameters for the numerical models of the NLCEs is introduced to better retrace the experimental results of the assembly.  The scanning laser vibrometer is used to obtain the experimental response of the system on a finite set of points of both the beams and the NLCEs with respect to an excitation provided by the shaker acting on a fixed point (Fig. 6). The measurement point locations on the bottom beam and on the NLCEs of the assembled structure are shown in Fig. 25. As for the NLCE, the measurement of the NNMs for the free-free system is performed by setting the amplitude and by adjusting the frequency of the sinusoidal excitation to identify the resonance conditions. After removing the low frequency contents due to rigid body oscillation in the free-free configuration, the velocity signal is integrated to achieve the output in terms of displacement for each node.
The NNMs of the system are also obtained numerically by using the coupling procedure on the updated substructures. The experimental identification of NNMs provides the resonance frequency and the amplitude of oscillation of measurement points, for each excitation level. Thus, the comparison between numerical and experimental NNMs cannot be performed using the FEP diagram since it is not possible to express the measured data as function of the energy. In the present case, since the right NLCE is significantly more deformed than the left NLCE, it is the main source of nonlinearity for the considered mode in the assembled system. Thus, it is possible to express the resonance frequency of the NNM as function of the deformation of the right NLCE providing a meaningful representation of the nonlinear behavior of the assembly. Figure 26a shows the experimental results (x-marked blue line) compared to the numerical results (solid green line) obtained performing the Coupling Procedure. The full modal basis of the numerical model of the beams (Sect. 5.1) is used, and the NNMs of the NLCEs are obtained using the coefficients reported in Table 5. The use of experimentally identified parameters leads to a general overestimation of the resonance frequency expressed as function of the deformation of the right nonlinear connection, with an average error of ≈ 9%.
However, as shown in Fig. 26b, the MAC between the experimental mode shapes and the numerical ones (green bars) is almost 1, indicating a perfect correlation, up to a deformation of the right spring of 2.51 mm and it decreases to 0.91 at 4.35 mm. It should be noted that the numerical and experimental mode shapes are consistent for values of deformation of the right spring inside the range for which its parameters (k l and k nl ) are experimentally evaluated (Fig. 24b).
In order to investigate the causes of the error in frequency highlighted in Fig. 26a, an optimization of the parameters for both the NLCEs is carried out. In fact, although the substructures to be coupled (i.e. beams and NLCEs) are singularly well identified, some difference between the numerical coupled model and the real one can arise. It could be due either to the fact that the real connections between the beams and the NLCEs cannot be assumed as perfectly rigid or to the clamping that could introduce some pre-stress on the NLCEs, changing their effective linear/nonlinear stiffness. Hence the optimization aims at modifying the linear and nonlinear terms of the NLCEs' law to account for the loosening effects of the real interface. The differences between the numerical frequencies and the experimental ones have been simultaneously minimized through a multi-objective optimization technique. In Fig. 26a the red line represents the Frequency-Deformation Plot obtained using the optimized coefficients for the cubic law, as reported in Table 6.
For both the NLCEs the optimized coefficients are lower than the ones experimentally identified, especially for the right one. This may be ascribed to the additional flexibility of the connection between the beams Table 6 Coefficients of the cubic law for the two NLCEs after optimization with percentage difference with respect to the ones in Table Fig. 27.
In order to highlight the differences between the mode shapes, the scale along the y axis is magnified. The similarity in the mode shapes shown in Fig. 27 agrees with the MAC values shown in Fig. 26b. Note that the NLCP manages to retrace the experimental mode shapes with high confidence (MAC is always ≥ 0.9). The experimental results agree with the hypothesis of considering only the response due to one harmonic in the NLCP.

Concluding remarks
The substructuring technique in the modal domain was successfully employed to study the dynamics of linear subsystems connected through nonlinear connecting elements. The adopted method has been significantly enhanced by enforcing a better control of the total energy at each step of procedure through an adjustment of the energy distribution among substructures. Also, the choice of the initial guess based on the results obtained at the two previous energy levels and the definition of a stopping criterion that compares the energy distribution between two consecutive iterations lead to a more rapid convergence of the procedure. The dynamic behavior of the coupled system is evaluated by computing its Nonlinear Normal Modes. The coupling procedure is validated by comparing its results with those of the HB using a single harmonic, thus it can be used to evaluate the main dynamic behavior of the system already in the early phases of the design. The effects of modal truncation have been studied by comparing the resonance frequencies and by analyzing the correlation between the mode shapes, using several indicators. The results show that the error in frequency depends both on the choice of the truncated modal basis and on the energy level at which the NNM is computed: generally, the errors increase for an increment of the energy level and for a decrease in the number of modes retained in the reduction matrix.
The specifically designed nonlinear connecting elements have been employed in the experimental setup to connect two linear beams. The hardening behavior of the connecting elements have been evaluated experimentally by measuring its NNM. Substructuring results compared to the experimental results show that there is a good correlation with the mode shapes measured on the assembled system. Furthermore, a good accuracy in the resonance frequency at different excitation levels requires a careful modeling of the NLCEs taking into account the effects of the connecting interfaces. In fact, there seems to be some uncertainty in the properties of the NLCEs, leading to some errors in the resonance frequency. The optimized parameters of the physical models of the NLCEs highlight an overestimation of the linear and nonlinear stiffness properties that could be due to loosening of the connecting interfaces or to the pre-stress effects introduced by the clamping. Note that, although the NLCP does not include the effects of higher harmonics, in the present case there is a good correlation between the numerical and experimental mode shapes. In fact, no significant effects of internal resonances and higher harmonics have been experimentally observed although the nonlinear elements were subjected to quite large deformations.
The results show that the substructuring procedure is altogether suitable to predict the dynamic behavior of linear systems coupled through nonlinear connections. Further studies could be addressed to take into account both material damping and damping localized at interfaces. The present approach seems to be suitable to be applied also to those cases in which the connections between the substructures need a more complex description, for example by considering multi-DoFs connection or hard to model interactions.