Self-adaptive macroslip array for friction force prediction in contact interfaces with non-conforming meshes

The steady-state nonlinear forced response (NFR) of finite element (FE) models with friction joints is usually computed in the frequency domain through the combination of node-to-node contact elements and the Harmonic Balance Method (HBM). In the current state of the art, rare are the cases where the friction forces are estimated for contact interfaces with non-conforming mesh grids. This need is nowadays taking place due to the improving capability of commercial FE software to manage any kind of boundary condition (i.e., either coupling or contact), without requiring coincident pairs of nodes at the joint interfaces. Such an advantage becomes a drawback when the analysts are requested to investigate the NFR of the assembly by using build-in codes, where the contact forces prediction usually requires node-to-node contact elements whose parameters (i.e., the contact stiffnesses and friction coefficients) can be easily identified by means of experiments. This paper addresses the mentioned limitation, and proposes a novel self-adaptive macroslip array (SAMA) model for the estimation of the nonlinear friction forces on FE contact interfaces with non-conforming meshes. The SAMA model consists on a set of node-to-node contact elements ordered in parallel, whose contact parameters and normal preloads are identified through a step-by-step self-adaptive weighting algorithm that depends on the topology of the meshes in contact. The goodness of the proposed model is assessed on the calculation of the NFR of a bladed disk with shroud contacts, under the hypotheses of cyclic symmetry and HBM. The nonlinear dynamic behavior of the bladed disk is evaluated in two different cases. First, in the case of lack of node-to-node congruence at the contact interface for the structure being in its undeformed configuration, and second, in the case of a relevant static misalignment of the contact interfaces due to the application of large static loads.


Introduction
A dependable mechanical design against high cycle fatigue (HCF) requires understanding the causes of static and dynamic loads for an accurate estimation of the fatigue limit. Static and dynamic loads have different nature depending on the field of application. In lowpressure and high-pressure turbine bladed disks for aircraft engines, static loads are mainly represented by the centrifugal force, temperature field, and steady pressure distribution on the blades airfoils, whereas dynamic loads can be identified with the unsteady component of the gas flow resulting in traveling wave engine order (EO) excitations, which are the major responsible for forced vibrations in blade arrays [1][2][3].
In order to reduce the occurrence of HCF damages, the design of a bladed disk assemblies is driven by the aim of mitigating the blades vibration amplitude. One of the commonest solutions for this task is the design optimization of the joints involved in the mechanical assembly, for the maximization of the friction damping taking places at the contact interfaces due to the occurrence of macro and micro-slip phenomena [4][5][6][7]. Typical applications where such an approach is successfully applied are the blade root and shroud joints for rotating bladed disks [8,9], and the interlocking and hook joints for stator vane segments [10].
Plenty of scientific papers have been published on this topic in the recent past years, and several numerical methods have been developed in order to efficiently predict the FR of structures in the presence of friction contacts. Most of these have been purposely developed for industrial applications, mainly in the turbomachinery field [11,12], but also concerning mechanical assemblies where bolted flanges and lap joints are used [7,13,14]. The key idea behind the NFR prediction is the estimation of the contact forces by means of appropriate contact models [15][16][17][18], and the solution of a reduced set of nonlinear equation of motion (EQM) in the frequency domain by using the harmonic balance method (HBM) [15] or the multi-harmonic balance method (MHBM) [19,20]. These methods allow to overcome the expensive computational costs typical of a direct time integration (DTI), which require to find the steady-state response of the structure for each excitation frequency. For large FE models (FEMs) with extended contact interfaces, DTI in commercial FE software might barely get the convergence, and sometimes make the dynamic analyses unfeasible due to the huge number of nonlinear DOFs involved in the numerical simulations.
Nowadays the use of commercial FE software is unavoidable when dealing with complex engineering systems, and lots of efforts are needed for the creation of (FEMs) able to capture the structural behavior that the analysts want to investigate. In this regard, an accurate meshing of the contact interfaces can make the difference for an easier and more precise evaluation of the contact forces that neighboring components mutually exert on each other during vibrations. A widely accepted and used practice requires to adopt conforming meshes at the contact interface. The term conformity here denotes the existence of well-defined pairs of nodes sharing the same geometric location on the contact plane. This condition allows the application of node-to-node contact elements whose parameters (i.e., the contact stiffnesses and the friction coefficient) can be experimentally identified [21,22]. However, the mentioned approach holds unless the mesh conformity remains preserved despite to the variable, but negligible mesh misalignments caused by vibrations.
FEMs of mechanical systems of engineering relevance are often characterized by the lack of mesh conformity at the contact interfaces. This might happen when the FEMs of the components belonging to the same structure are created and meshed independently by different engineering teams. In other circumstances, despite the initial mesh conformity, the application of large static loads might force the contact interfaces to achieve a relative positioning where conformity of the meshes is no more guaranteed. In both cases, node-tonode contact elements cannot be applied as usual, but new solution strategies are needed to let the contact interfaces interact without changing the contact element type, and still preserve the reliability of the predicted NFR of the system. This paper presents a novel self-adaptive macroslip array (SAMA) model for the contact force calculation on FE surfaces with non-conforming meshes and provides an analytical formulation of the Jacobian matrix computed at each iteration step of the Newton-Raphson method (NRM) for the solution of the EQM in the frequency domain. The proposed method is applied for the NFR prediction of a bladed disk assembly with shroud contacts.
In Sect. 2, the reduced nonlinear EQM of a rotationally periodic structure with friction contact are obtained assuming the hypothesis of cyclic symmetry (CS) and the MHBM. Section 3 recalls the usage of the commonest node-to-node contact elements in the technical literature and gives an insight on how the EQM are solved in the frequency domain using the NRM. Section 4 describes the algorithm behind the novelty of the paper, i.e., the SAMA model. Sections 5 and 6 are focused on the application of the SAMA model for the contact forces calculation at the shroud of a bladed-disk reduced order model (ROM). A set of benchmark NFRs on a bladed disk ROM with conforming meshes at the shroud (i.e., the ROM1) and node-to-node contact elements is computed to achieve two different outcomes. The first is the validation of the SAMA model employed on a ROM with non-conforming shroud meshes (i.e., ROM2), by comparing its NFRs versus the benchmark ones. The second is to prove the relevance of considering a misalignment of the shroud due to the application of a set of static loads, and its effects on the NFR of the structure. In this regard, a nonlinear FE static analysis on the model with conforming meshes is performed to find the misalignment and normal preload at the shroud, which represent the equilibrium contact conditions for the nonlinear forced vibrations. The SAMA model is employed in the NFR to overcome the loss of mesh conformity due to the shroud misalignment. The resulting NFRs are compared to the benchmark ones in order to quantify the differences in terms of amplitude, resonance frequency and Q-factor.

Equation of motion of a cyclic symmetric sector
The nonlinear EQM of a structure with friction contacts can be written as: where M, K and C are the mass, stiffness, and viscous damping matrices of the system, q is the generalized coordinates vector, while f e and f c are the vectors of the excitation and nonlinear contact forces, respectively. The dot above the vector q denotes the derivative with respect to the time variable t. For large FEMs with a high density of nodes, reduction algorithms in the class of component mode synthesis (CMS) [23,24] can be used to write the Eq. 1 in a reduced form, meaning that the system's dynamics can be condensed on a limited set of master DOFs suitably defined. For a rotationally periodic structure with n s identical sectors the nonlinear EQM can be obtained in a more compact form as described in the following. The EQM of the sth sector can be written as: where the matrices M s , K s and C s are the mass, stiffness and viscous damping matrices of the single sector, q s is the displacement vector of the sector expressed in a cylindrical reference frame fixed to the rotor (Fig. 1), f e s is the corresponding excitation force and f c s (q s ,q s , t) is the vector of the nonlinear forces resulting from the interaction of the contact interfaces of the sth sector.
In regular operating conditions, turbine blade arrays experience a wide spectrum of traveling forces whose Fig. 1 Cyclic symmetric shrouded bladed disks. The cyclic symmetry constraints are enforced to the DOFs at the right and left sector frontiers spatial distribution consists of an integer number of wavelengths having excitation frequencies equal to an integer multiple of the engine's rotation speed [3]. Such time-varying forces are caused by the non-uniform gas flow interacting with the blades. Because of the periodicity of traveling wave excitation, the excitation force for the full bladed disk can be expressed as a truncated Fourier series of n h harmonics: Assuming a matrix partitioning such that the n s DOFs of the sth sector are followed by the n s DOFs of the sector s + 1, and so on, f e s can be written as: where s = 1, . . . , n s , φ m = 2π m/n s is the phase shift between the force components on two adjacent sectors, τ = m r t the non-dimensional time variable, with m the spatial periodicity of the traveling wave force (i.e., the engine order) and r the rotation speed.
The prediction of the bladed disk's mode shapes for a prescribed number of nodal diameters h requires to enforce CS constraints to the DOFs at the right and left frontiers of the bladed disk sector [19,[25][26][27].
With reference to Fig. 1, the CS constraints can be written as: where: q o s is the partition with the physical DOFs non belonging to the disk frontiers and the modal coordinates of the CMS method eventually used to create the single sector ROM; -I is the identity matrices with dimensions consistent either to the vector x r s or q o s ; -ϕ h = 2π h/n s is the inter-blade phase angle; -T h is the CS constraint matrix; -q h s is the vector of generalized coordinates for the fundamental sector in CS conditions. The application of CS constraints requires to substitute Eq. 5 into Eq. 2 and the pre-multiplication of each of its term by the conjugate transpose of T h . The mentioned transformation coincides to a Galerkin's projection, which leads to the following EQM in CS conditions: where: Note that the transformation applied to the viscous damping matrix holds if C s is defined under the hypothesis of proportional damping, i.e., C s = αM s + βK s . A possible alternative for the definition of C h s would be applying the inverse modal transformation to the modal damping matrix C h defined as: where h s is the modal matrix with the n m eigenvectors of the eigenproblem associated with Eq. 8, ω n j and ζ j the eigenvalue and the modal damping ratio of the jth CS mode shape, respectively, [3]. Furthermore, due to the orthogonality of the trigonometric functions an excitation force f e with a pure spatial harmonic distribution of order k can excite a mode shape having periodicity h if the following condition is met: For z = 0 the engine order of the excitation force coincides to the number of nodal diameters h of the mode shape (i.e., k = h). For z = 0 aliasing occurs and the force projection onto the mode shape is not null even if k = h [28]. This means that the harmonics components of the force f h e s that actually excites the CS structure are those for which k = h · n satisfies Eq. 9.
Since the EQM are formulated for a single sector in CS conditions, in the following the superscript h and the subscript s will be omitted for brevity.
The steady state response of the system can be obtained solving the EQM in the frequency domain using the harmonic balance method (HBM) [15,19,20]. Due to the periodicity of f e , the DOFs vector q and the nonlinear contact forces f c can be expressed as a truncated Fourier series: where Q (n) and F (n) c are the n th -order complex amplitudes of q and f c , respectively. Substituting the Fourier series of Eqs. 4 and 10 into Eq. 6, the following system of nonlinear complex algebraic equations is obtained: with: where the coefficients F (n) c , which are typically obtained from the time domain contact forces, depend on the complex amplitudes Q (n) and their derivatives for n = 1, . . . , n h .

Contact model and equation of motion in the frequency domain
The problem of modeling periodical contact forces and their implementation in numerical solvers for the NFR prediction of structures in the presence of friction contacts has been addressed by several authors. In the past years the research has led to several node-to-node contact models for the calculation of contact forces given the relative displacement of a node pair [18]. The choice of the suitable contact element depends on the kinematics of the contact pair. For a 3-D periodic relative displacement, a contact element able to capture the normal load variation on the contact surface has to be preferred [29,30]. In this regard, two different modeling techniques can be considered: -to use two independent 1-D penalty contact elements with normal load variation [4,31], whose tangential directions (i.e., the x and y direction on the contact plane) are considered orthogonal to each other, in order to capture the in-plane 2-D trajectory of the relative displacement on the contact surface ( Fig. 2a); -to use a 2-D penalty contact element with normal load variation, which accounts for the coupling between the two orthogonal in-plane components of the tangential relative displacements [5] ( Fig. 2b).
According to the schemes reported in Fig. 2, the contact parameters for both the models are the tangential and normal stiffness k t and k n respectively, the friction coefficient μ and the static normal preload n 0 . u(τ ), v(τ ) and g(τ ) are the input relative tangential and normal displacements, and the slider's displacement at the time instant τ , while f t (τ ) and f n (τ ) represent the tangential and normal contact forces, respectively. The subscripts x and y denote in-plane direction along which the displacements and forces are projected. The projection of a 2-D in-plane trajectory onto two straight and orthogonal directions leads to an underestimation of the friction damping. In fact, at the time instant τ , the actual tangential friction force f t (τ ) (Fig. 2), defined as the vector sum of f t x (t) and f t y (τ ), in general is not oriented as x or y, and has a magnitude larger than |f t x (τ )| or |f t y (τ )|. This means that there might occur cases where |f t (τ )| gets the Coulomb limit μ|f n (τ )| (and sliding occurs), even if the amplitudes of f t x (τ ) and f t y (τ ) are so that |f t x (τ )| < μ|f n (τ )| and |f t y (τ )| < μ|f n (τ )|. Therefore, although less precise, the first modeling technique is more conservative than the second one from the point of view of the dynamic design, since it predicts larger vibration amplitudes as described in [18]. Furthermore, in the technical literature, there are several studies where the assumption of independent tangential component of the displacement is used [32,33]. For these reasons in Sect. 4 the SAMA model will be presented with reference to the 1-D contact elements with normal load variation, but no restrictions occur when the 2-D contact element with normal load variation is used.
The 1-D node-to-node contact element with normal load variation is used to compute the periodic contact forces for a given periodic relative displacement, by taking into account possible separation of the nodes in contact. The contact element allows to model three different contact states: stick, slip and separation. A schematic view of this contact model is provided in Fig. 3, where the two-dimensional relative displacement is decomposed into two perpendicular directions: two in-plane tangential component denoted by the u(τ ) and g(τ ) components, i.e., the nodes' tangential relative displacement and the slider's displacement respectively, and one out-of-plane normal component v(τ ). The contact model parameters are the same described for Fig. 2.
At every time instant the normal contact force f n (τ ) is defined as: If n 0 is positive, the bodies are in contact before vibration starts, while if n 0 is negative an initial gap g 0 = − n 0 k n exists between the two nodes. In the tangential direction, the contact force is defined as: The stick, slip and separation states alternate each other during the vibration period according to the transition criteria reported in [4]. The contact models work only in the time domain, but the Fourier coefficients of the contact forces F (n) c participate to a dynamic equilibrium that is written in the frequency domain (i.e., Eq. 12). The alternating frequency time (AFT) method [34] is therefore used to switch between the frequency and time domain by using the fast Fourier transform (FFT) and the inverse fast Fourier transform (iFFT) algorithms as shown in Fig. 4.
Since the nonlinear contact force only depends on the relative displacement of the contact DOFs, the size of the EQM can be further reduced by keeping only the nonlinear DOFs in the system. This is done by partitioning the vector Q (n) and Eq. 12 into their linear and nonlinear (i.e., those related to the contacts) components: being the linear partition of the nonlinear contact forces equal to 0 L , since no contact forces are usually applied to the linear DOFs. As reported by several authors (e.g., [33]), the linear DOFs depend on the nonlinear ones. It is therefore possible to solve the nonlinear partition of Eq. 15 for the nonlinear entries only: with: while the linear DOFs are evaluated from the nonlinear ones by using the following relationship: The set of Eq. 16 can be solved by applying an iterative solution scheme such as the NRM [35], in order to minimize the norm of the residual vector r: where: r = (r (0) ) T , . . . , (r (n h ) ) T T , where the nth-order residual vector r (n) is associated with the nth-order Eq. 16; The NRM generates an approximate solution that at each iteration step approaches to the roots of the system of Eq. 16. The iterative step for Eq. 19 can be expressed as: where Q N| i is the response vector at the ith iteration, r| i-1 is the residual vector at the iteration i-1 and J| i-1 is the Jacobian matrix at the iteration i-1 that reads: It must be pointed out that although the equilibrium equation are written using complex arithmetic, the vectors and matrices of Eqs. 19 and 20 allow for real entries only, meaning that the vectors in Eq. 19 have to be separated into their real and imaginary components and thē D matrix rearranged accordingly. In particular, if n N is the total number of nonlinear DOFs, Q N is turned into its real counterpart Q N as follows: and the real counterpart ofD becomes: where: The core of the Jacobian matrix J is represented by the matrix J a , which contains the partial derivatives of the Fourier coefficients of the nonlinear contact forces with respect to the Fourier coefficients of the absolute displacements of the nonlinear DOFs. The definition of J a can be carried out from the matrix J r , which contains the partial derivatives of the Fourier coefficients of the nonlinear contact forces with respect to the Fourier coefficients of the relative displacements of the nonlinear DOFs. For a single node-to-node contact with 3 DOFs (x, y and z), J r can be expressed as: where: -X , Y and Z are the complex amplitudes of the physical relative displacement of the node pair in the x, y, and z direction, respectively. -F tx , F ty are the complex amplitudes of the tangential contact forces in the local x and y directions, while F nz denotes the complex amplitudes of the normal contact force in the local z direction.
The finite difference approximation of J r is highly expensive from a numerical point of view. Therefore, by following the approach proposed by Cardona et al. [36], for the efficient iterative minimization of the residual vector norm, the Fourier coefficients of the contact forces as well as the partial derivatives of Eq. 25 can be computed using the AFT algorithm. Further details on the analytical formulation of the partial derivatives in Eq. 25 can be found in [19]. J a can be obtained by applying the following relative-to-absolute coordinate transformation to J a : where: The subscripts used for I denote the size of the identity matrices (6 is the number of DOFs of two nodes in a three-dimensional space), while the symbol ⊗ defines the Kronecker product.

Self-adaptive macroslip array model
In the technical literature, most of the contact problems in structural dynamics are solved assuming the conformity of the FE meshes at the contact interfaces. From a practical point of view, this assumption comes from the capability of the meshing tools in commercial FE software that easily allow to replicate the mesh topology of one contact interface on its counterpart. Such an approach is extremely convenient since it allows any node-to-node boundary condition by just exploiting automatic routines that are able to match the node pairs given a certain tolerance on their mutual distance. The contact surfaces so obtained are therefore suitable for the application of node-to-node contact elements (i.e that described in Sect. 2) and allow to consider valid all the considerations carried out for the solution of the nonlinear EQM in the frequency domain. The great advantage of using node-to-node elements is due to the possibility of an easy experimental estimation the contact parameters μ, k t and k n . However, that are applications where the non-conformity of the contact meshes does not guarantee the identification of well-defined pairs of nodes, and new solution strategies are needed to let the usage of node-to-node elements possible. The novel self-adaptive macroslip array (SAMA) model is here formulated to overcome the mentioned limitation. In the following, the SAMA method is presented with reference to two 1-D contact interfaces (i.e., the contact and the target) discretized by 1-D elements, whose nodes represent the locations where the SAMA model applies. This choice aims for an easier interpretation of the following schemes, and the algorithm behind the SAMA model is presented for surface elements with an arbitrary number of nodes. The main hypothesis behind the usage of node-tonode contact elements is to have two nodes at the same geometric location within a certain tolerance. Although this hypothesis is not so strict when the contact interfaces behave as rigid bodies during in-plane vibratory motion, it is necessary when the contact interfaces are flexible or particularly wide. In the latter cases, the displacement of a certain node cannot be considered representative of the whole displacement field of the corresponding contact interface. This idea can be better visualized in Fig. 5 where the motion of two contact interfaces is shown. The red and blue dashed segments represent the undeformed contact and target elements respectively, while the continuous lines denote the same deformed interfaces at a certain time instant τ during vibration. It can be noted that no relative displacement can be computed at the location of the contact node c, since no node counterpart exists on the target surface at the same geometric location. Furthermore, it is also worth to mention that a non-reliable estimation of the contact forces would be made if the relative displacements between the contact node c and the closest target t are fed into the penalty contact models introduced in Sect. 3.
The problem of studying the contact interaction between the non-conforming contact and target interfaces is solved by the following step-by-step algorithm: 1. given a contact node c on the contact side, this is projected onto the target side, so that a virtual node t v with coordinates x t v and y t v in the local reference frame of the contact can be identified (Fig. 6). 2. for each t v all the corresponding target nodes must be identified: these are in the number of n T and define the 2-D area on the target surface (i.e., the target element) containing the virtual node t v (Fig. 6). Note that in the 1-D representation of Fig. 6 n T = 2.
For a 2-D case, depending on the order of the elements used for the surface mesh, n T ≥ 3. 3. for each node pair t v -t j the distance d j between the virtual node t v and the jth target is found as: where x t j and y t j are the x and y coordinates of the jth target node. 4. given all the distances d j as in Eq. 28, for each target node the weights w j are computed as: so that: Moreover: -if d j = 0 and w j = 1, the nodes c and t j share the same geometric location and all the other weights are equal to zero (i.e., w i = 0 for i = j). From a physical point of view, the nodes pair c-t j is defined as for conforming meshes; -w j = 0 if at least one weight w i with i = j is equal to 1; -if neither of the previous two conditions occur it happens that 0 < w j < 1. The smaller d j the larger the weight w j ; 5. once all the weights w j are computed, any of the penalty contact elements introduced in Sect. 3 can be applied for each pair c − t j , assuming the following weighted contact parameters: friction coefficient k t j = w j k t tangential contact stiffness k n j = w j k n normal contact stiffness n 0 j = w j n 0 static normal preload (31) Fig. 6 The contact node c is projected onto the target interface, so that the virtual node t v is found and the corresponding target nodes t 1 and t 2 are identified Note that for the condition of Eq. 30 it holds that: meaning that the global contact stiffness (in the tangential and normal direction) and preload are preserved regardless the non-conformity of the meshes; 6. the contact forces are finally computed by feeding the tangential and normal relative displacements u c − u t j and v c − v t j to the contact model with weighted contact parameters.
The application of the weighted parameters of Eq. 31 results in weighted contact forces at the target nodes t j . It is worth to mention that the application of the contact element in Fig. 3 does not involve strong approximations of the contact forces if the target mesh is sufficiently refined. Figure 7a shows a schematic representation of the SAMA model which resembles the well-known macroslip array [11,15] (Fig. 7b).
Both contact models consists of a set of parallel node-to-node contact elements (Fig. 3). The macroslip array aims to divide a single contact element with contact parameters μ, k t and k n , in an array of N contact elements, whose contact stiffnesses k t i , k n i with i = 1, . . . , N , are computed by following the procedure described in [5,15]. Each component of the macroslip array is subjected to the same relative tangential and normal displacement, but the different contact parameters lead to different hysteresis loops. In the SAMA model, the different contact parameters (Eq. 31) and relative displacements self-adapt to the topology Fig. 7 Self-adaptive macroslip array model (a) and macroslip array model [11,15] (b) of contact and target meshes. Note that the weighting algorithm previously described is not only valid for the node-to-node nonlinear contact elements, but it can also be used for the assembly of two interfaces with non-conforming meshes by means of linear spring elements having local tangential and normal stiffness k t j and k n j , respectively.
For each contact pair c − t j the matrix J r can be written exactly as in Eq. 25, being X , Y and Z the Fourier coefficients of the relative displacement u c −u t j and F t x , F ty and F nz the Fourier coefficient of the contact forces resulting from the contact model with weighted contact parameters (Eq. 31). After the coordinate transformation of Eq. 28, the global Jacobian matrix is obtained by assembling all the elementary Jacobian matrices J a . Given a single SAMA element, due to the multiple correspondence (c − t j ), the assembly of the global Jacobian matrix requires to sum up all partial derivatives corresponding to the contact DOFs, while no repetitions occurs for the partial derivative corresponding to the jth target node t j .
It must be pointed out that the contact forces at the location of the contact node c could have been computed by evaluating the relative displacement between c and a virtual target node t v , whose displacement could be evaluated from all the u t j through the element shape functions (Fig 8).
Although this process would have been convenient for the force evaluation, a further redistribution of the contact forces at the target nodes t j would have been performed, since the DOFs vector just contains the displacements of the target non-virtual nodes. However, such a force redistribution is not mathematically possi- Fig. 8 Displacement of the virtual node t v on the target surface ble, because the shape functions cannot be used to find the nodal forces at the target nodes t j from the known contact force inside the target element at the location of the virtual node t v .

Alternative procedure for cyclic symmetry boundary conditions
The dynamic behavior of the CS structures in turbomachinery is significantly influenced by the presence of the mechanical joints between neighboring sectors. Typical example are the turbine bladed disks where the blades interact through the shroud contact (Fig. 9a) and the stators where the contact between the vane segments is guaranteed by the preload at the interlocking (Fig. 9b). The numerical investigation of the friction damping mechanism in such joints during vibration requires the estimation of the relative displacement at the contact interface for the evaluation of the friction forces. If the disk dynamics is studied in CS conditions as described in Sect. 2, the relative displacement at the contact interfaces of a joint located at the sector frontiers (e.g., the shroud or interlocking) can be found by exploiting the same CS equations applied to the frontiers of the sector where the continuity of the material is guaranteed (i.e., the disk and the casing). According to the scheme of Fig. 10, the relative displacement u rel for a contact node pair at the shroud is computed as: (33) where the left contact interface displacements u l of the sth sector are shifted of an inter-blade phase angle ϕ h in order to obtain the homologous displacement of the sector s + 1. u rel is finally fed into the contact model and processed with the AFT method for the prediction The industrial practices for the mechanical assessment of a structure are based on the results of severe analyses performed on FE commercial software. In this regard, FE contact analyses are nowadays taking hold to better predict the role that the joints play on the static For instance, Ansys APDL provides well-established routines for contact analyses which require the application of layers of "contact elements" on the surface of the solid elements defining the contact interfaces. In most of the cases such routines work properly when the contact interfaces are close to each other so that just a small gap or interference occurs between them. For a contact FE analyses in CS conditions aimed at studying the blades interaction through the shrouds, the physical distance between the left and right contact interface does not allow to apply the contact elements as the practices require. This problem finds a practical solution in a method largely employed in the industry for the easy application of any boundary condition at the shroud of a blade sector in CS conditions. The method is schematized in Fig. 11 where the elements associated with one of the contact interfaces (e.g., the left one) are removed from their original position and translated in a cylindrical reference frame of a sector angle towards the other, non-modified contact interface (e.g., the right one).
This operations leads to two different results. First, the CS contact pair is turned into a standard contact pair, where the term standard is used to denote a pair with contact surfaces that are in proximity. Second, CS constraints originally applied to the disk frontiers must be extended to the DOFs lying of the surfaces along which the element cut has occurred. According to the nomenclature adopted in Fig. 1, the frontiers DOFs vectors x r s and x l s of Eq. 5 become: where the subscripts d and c denote the disk and cut frontiers DOFs.

Numerical example
In this Section the SAMA model is used to calculate the contact forces at the shroud of a bladed disk, for the prediction of the NFR of the system, for different values of the static preload n 0 .

Model order reduction
The bladed disk has 100 identical sectors consisting of a disk sector and a shrouded blade assembled together at the blade root joint (Fig. 12). Two different FEMs of the fundamental sector were created in Ansys by using 10-node-tetrahedral elements and the standard material properties of steel (Young's modulus E = 210 GPa, Poisson's ratio ν = 0.33 and density ρ = 7800 kg/m 3 ). The first FEM (FEM1) has conforming meshes at the shroud and consists of 51763 elements with 90042 nodes. The second FEM (FEM2) has non-conforming meshes and consists of 50548 elements with 88100 nodes. For both FEMs the method described in Sect. 5 was used in order to have the contact and target interfaces sharing the same contact plane, i.e., overlapped one on the other, as shown in Fig. 13. Figure 14 shows the locations of the contact and target nodes at the shroud interface, for FEM1 and FEM2 respectively.
A preliminary nonlinear static analysis in CS conditions for the harmonic index h = 0 was set on both FEMs by applying a radial temperature profile, a steady-pressure distribution at the blade's airfoil, a rotational speed of 250 rad/s and a mechanical interference of 2 mm at the shroud. A contact pair based on the surface-to-surface Ansys contact elements, i.e., the CONTA174 and TARGE170 in the augmented Lagrangian formulation [37], was applied at the shroud to achieve two different outcomes. The first was to prestress both FEMs as happens for real structures due to the application of static loads. The second was to evaluate the pressure distribution at the shroud (i.e., the normal preload), and the static misalignment of the contact interfaces that represents the equilibrium configuration from where the forced vibrations are computed as described in Sect. 3. Figures 15 and 16 show the interface misalignment for both FEMs as well as the preload distribution at the contact side.
Two reduction techniques were used to condense the disks dynamics on few master DOFs. First, a prestressed Craig-Bampton (CB-CMS) reduction was carried out in Ansys to into account the stiffening effect due to steady loads in the static analysis. This reduction led to a huge compression of the FEMs DOFs, since just the contact, the CS frontiers and few accessory DOFs were retained as a master. The details on the CB-CMS reductions can be find in Table 1.  Second, after having imposed the CS constraints for a harmonic index 10, the Gram-Schmidt interface (GSI) reduction method [38] was used to reduce the frontiers DOFs (i.e., the disk sector and cutted shrouded frontiers together) with 100 GSI modes. An alternative to the GSI method is the use of the Tran technique [39]. At the end of the reduction process two ROMs, i.e., the ROM1 and ROM2 for FEM1 and FEM2, respectively,  Table 2.
For both the FEMs and ROMs a modal analysis for h = 10 was performed without enforcing any boundary condition at the shroud, meaning that the interaction of the contact interfaces was completely prevented (free condition). The Modal Assurance Criterion (MAC) was used to compare the ROMs and FEMs modal displacements at the physical master DOFs (i.e., accessory and contact) (Fig. 17), while the natural frequencies were assessed by computing the relative percentage difference of the ROMs quantities with respect to the FEMs ones. The large values of the diagonal MAC matrix (∼ 1) as well as the small maximum percentage difference on the natural frequencies (< 0.1%) confirm the goodness of the reduction process employed for both FEMs. The NFR analyses were performed for a harmonic excitation applied at the mid-leading edge of the blade (Fig. 17). The force amplitude was tuned in order to get realistic displacements at resonance for the two mode shapes of interest, i.e., the 2nd bending (M1) and torsional (M2) modes shown in Fig. 18. Two sets of forced responses (i.e., one per mode shape) were computed on ROM1 without considering the interface misalignment resulting from the nonlinear static analysis, and assuming the parameters in Table 3 for the penalty contact elements in Fig. 2a.
Due to the mesh conformity at the shroud of ROM1, these sets were obtained by using the classic approach based on the application of the elements as shown in Fig. 2a, and represent the benchmark analyses used to achieve two different outcomes. The first is the validation of the SAMA model employed on the ROM2, Fig. 18 Mode shapes of interest: 2-nd bending (a) and torsional (b) mode shape for a fully stuck contact condition at the shroud Fig. 19 Nonlinear FRFs of the node where the force is applied for the mode shape M1, for different n 0 at the shroud contact interface. The vibration amplitudes are normalized with respect to the highest amplitude FRF obtained for the second mode shape, i.e., that corresponding to n 0 = 300 N (Fig. 20) where no conformity of the mesh at the shroud occurs, while the second is to prove the relevance of considering a misalignment of the shroud due to the application of static loads, and their effects on the NFR of the bladed disk.
The benchmark FRs are shown in Figs. 19 and 20 where the typical stick-slip phenomenon occurs for both mode shapes. For a fixed value of the external load an increase of the contact preload has the effect to shift the resonance towards higher frequencies, moving from a contact condition where sliding occurs to the so called fully-stuck condition, i.e., the condition for which no sliding occur since the contact elements behave as springs.
The same forced response analyses were then performed on the ROM2, where the shroud contact forces were computed using the SAMA model. The analyses results for n 0 = 113 N (i.e., the reference static preload) are shown in Figs. 21 and 22, where the comparison with respect to the benchmark FRFs is provided.
The comparison shows negligible differences in terms of vibration amplitude and resonance frequency that are probably due to numerical approximations.
A full comparison of the forced responses obtained from the two ROMs is reported in Figs. 23, 24 and 25, where the resonance frequencies, the response amplitudes at resonance and the Q-factor are plotted as a function of the static preloads reported in Table 3.
The overlapped plots on each diagram denote the total consistency of the forced responses and confirm the goodness of the SAMA method for the prediction of  the contact forces on surfaces where no node-to-node contact pairs can be identified.
A further set of NFRs was computed on the ROM1 considering the interface misalignment obtained as a result of the nonlinear static analysis (Fig. 15a). Due to the non-conformity of the meshes caused by the interfaces misalignment, the contact forces prediction required the use of the SAMA method. The comparison versus the benchmark FRFs is here provided for n 0 = 113 N (Figs. 26 and 27).
For the mode shape M1 (Fig. 26) it can be noted a decrease of the resonance frequency of the 2.12% (M1 -conforming taken as reference) and an increase of the vibration amplitude of the 1.13%. The reason behind such differences is due to the resulting smallest contact area at the shroud caused by the misalignment, which makes the joint's looseness higher and decreases its equivalent stiffness during vibration. For the mode shape M2 the misaligned contact condition does not have so much effect on the resonance fre- quency (0.14% difference), but leads to a decrease in the vibration amplitude larger than the 15%.
In this case the kinematics of the misaligned shroud is such that slipping is privileged, which causes an increase of the friction damping and a decrease of the vibration amplitude at resonance.
The full comparison of the FRFs is shown in Figs. 28, 29 and 30, where the performance plots in terms of vibration amplitude at resonance, resonance frequency and Q-factor are provided.
In general the trends still confirm the occurrence of differences on the FRFs between the cases of conforming and misaligned shroud interfaces. In particular, for the mode M1 no remarkable changing occur for the vibration amplitudes and Q-factors, while differences close to the 2% still occur in terms of resonance frequencies for all the preloads. For the mode M2, no important variations of the resonance frequency take place, while an increasing difference in terms of  . 26 Comparison of the FRFs obtained from ROM1 for conforming and misaligned meshes, for a preload value of n 0 = 113 N and for the mode shape M1 Fig. 27 Comparison of the FRFs obtained from ROM1 for conforming and misaligned meshes, for a preload value of n 0 = 113 N and for the mode shape M2 vibration amplitude and Q-factor occurs as the preload increases. Note that the smaller amplitudes for the case with misaligned shroud (red plot in Fig. 28b) corresponds to smaller Q-factors in Fig. 30b, meaning that the amplitude variation is due to an increase of friction damping.
The presented results prove the relevance of capturing to consider the effects of the static relative position of the contact interface in un-coupled dynamic calculations. If the static and dynamic equilibrium equations are solved together (coupled approach [18]), the stepby-step algorithm described in Sect. 4 still holds. This must be repeated at each iteration step of the NRM, so that the static misalignment is continuously updated, and the contact forces can be evaluated for different relative positions of the contact meshes until the numerical converge is found.

Conclusions
This paper presents a novel self-adaptive macroslip array (SAMA) model for the calculation of friction contact forces in FE contact interfaces with nonconforming meshes. The SAMA model consists of a set of parallel node-to-node contact elements whose parameters (i.e., the contact stiffnesses and the normal preload) and relative displacements are identified by means of a self-adaptive weighting algorithm depending on the topology of the meshes in contact. In this paper the SAMA model has been used to predict the NFR of a ROM of a bladed disk in CS conditions, in the presence of friction contact at the shroud. Two different cases have been presented and developed. First, the contact interfaces at the shroud are considered overlapped to each other, but with a different mesh grids that do not allow the identification of well-defined node pairs. This condition might happen when the FEMs of two or more components in contact are created independently. In such a case the SAMA model has been successfully validated on a set of benchmark FRFs obtained on the same bladed disk, but with conforming meshes at the shroud. The FRFs computed in the lack of node-to-node coincidence have been found identical to benchmark ones. Second, starting from a condition where the contact meshes are considered conforming, due to the application of a set of static loads, the shroud interfaces are led to a misaligned relative positioning where the initial conformity of the meshes is no more guaranteed. In this case that SAMA model has been used to prove the effectiveness of the contact interface misalignment in modifying the nonlinear benchmark FRFs. It has been observed that depending on the mode shape, a non-negligible changing in the resonance frequency or in the vibration amplitude might occur. In particular, for the bending mode shape the misalignment has led to a 2.12% decrease of the resonance frequency and to a 1.13% increase of the vibration amplitude. For the torsional mode shape (M2), although no meaningful changing in the resonance frequency have been observed, the mesh misalignment has caused a decrease of the vibration amplitude that has been quantified in more than a 15%. These results suggest the importance of capturing the nonlinear dynamic behavior of friction joint in working conditions far from the nominal one. The SAMA model has been here applied with reference to the 1-D node-to-node contact element with normal load variation. However, the self-adaptive weighting algorithm used for the identification of the contact parameters can be applied without any restriction to even more advanced node-to-node contact elements.
Funding Open access funding provided by Politecnico di Torino within the CRUI-CARE Agreement.
Code availability The code is written according to the proposed model.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.