An equivalent model of a nonlinear bolted flange joint

Abstract The dynamic response of individual components in an assembled structure shows high accuracy compared to experimental measurements of the system response. However, when it comes to assemblies, the conventional linear approaches fail to deliver good accuracy, due to the uncertain linear and nonlinear mechanisms in the contact interface of the joints. Therefore, the inherent dynamics of the contact interfaces needs to be considered in modeling assembled structures. In this paper the prediction of the nonlinear dynamic response in a bolted flange joint was obtained in two ways. First, a 3D detailed finite element model capable of representing the micro-slip mechanism was made using a quasi-static time stepping analysis. The linear characteristics and nonlinear mechanisms developing in the contact interface of a bolted joint are investigated by using the 3D detailed model. Moreover, the natural frequencies of the assembled structure (representing the linear response) and the micro-slip behavior in terms of hysteresis loops (representing the nonlinear response) are obtained using the detailed model. Second, an equivalent model composed of beam elements and an appropriate joint model is then constructed for the assembled structure. An identification approach is proposed, and the parameters of the joint model are identified using both linear and nonlinear characteristics, i.e. natural frequencies and hysteresis loops. Comparing the hysteresis loops obtained from the detailed and equivalent models verifies the accuracy of the joint model used to represent the contact interface and the identification approach proposed for parameter quantification.


Introduction
Most mechanical structures are composed of various substructures assembled by mechanical joints. Amongst the various types of mechanical joints, bolted joints are widely used due to their simplicity in assembling mechanical structures and their ability to transfer forces and moments between substructures. The main advantages of bolted joints over other types of connections such as welded joints are their acceptable strength, safety and durability. However, despite their simplicity, their inherent dynamics is too complex to be easily modeled and analyzed. In contrast to welded joints, bolted joints do not show reliable performance under dynamic loads and vibration, including the loosening phenomenon which is often avoidable by using locking mechanisms.
The damping, stiffness, and dynamic response of mechanical assemblies are significantly affected by the physics of the contact interfaces of joints and is one of the main sources of nonlinearity and uncertainty in structural dynamics. In spite of having received considerable attention in the literature, the constitutive behavior within the contact interface of such https://doi.org/10.1016/j.ymssp.2020.107507 0888-3270/Ó 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). joints is still a source of uncertainty. The behavior of the contact interface of bolted lap and flange joints is displacement dependent. At low excitation force amplitudes, the behavior is linear because the relative displacement at the contact interface is not enough to initiate nonlinearities. Bolted lap and flange joints show specific nonlinear behavior particularly at higher excitation force amplitudes resulting from nonlinear slip and slap mechanisms. Usually these mechanisms manifest themselves as nonlinear stiffness [1] and damping [2] effects and also as the transfer of mechanical energy from lower to higher frequencies [3]. Due to the friction between the mating surfaces, micro/macro slip can develop in the tangential direction at the contact interfaces resulting in hysteretic behavior when a cyclic load is applied. The area enclosed in the hysteresis loop represents the amount of dissipated energy in the contact interface in one cycle [4] when a harmonic force is applied. The specific behavior of bolted joints originates from the two nonlinear slip and slap mechanisms occurring simultaneously in the contact interface. This is mainly due to the special configuration of bolted joints. These mechanisms have an interchange effect on each other such that when slap (or micro-vibro-impacts) happens in the normal direction, the distribution of normal forces in the contact interface changes, which alters the state of slip behavior in the tangential direction. The investigation of a proper mathematical formulation for the dynamics of mechanical joints started more than a half century ago [5,6,7]. However, this is an ongoing effort and has received increasing attention recently. Many models have been proposed that consider the nonlinearities associated with hysteretic and are capable of reproducing the properties of the contact interfaces in terms of energy dissipation and nonlinear stiffness. Moreover, the model parameters must be identified from joint level experimental data or an accurate finite element modeling of joint interface and the integration of the joint model into structural-level models must be practical [3].
The first attempt to model hysteretic behavior was undertaken by Timoshenko [8]. Later, Iwan [6,9] extended the approach suggested by Timoshenko which assumes that a general hysteretic system is composed of a large number of elastoplastic elements with different yield levels. Iwan [9] introduced one of the earliest models of mechanical joint interfaces consisting of a series of Jenkin's elements in parallel or series combinations. This model has been the basis of many other joint models developed thenceforth [10,11,12,13]. Gaul et al. [14] experimentally studied the nonlinear damping and response function characteristics of bolted joints. Lenz and Gaul [15] and Gaul and Lenz [16] studied the behavior of a bolted joint when it is subjected to longitudinal and torsional forces. They used the Valanis model, which is well known in plasticity [17], as a reduced-order lumped model, and showed that this model is capable of describing the nonlinear dynamic behavior of a bolted joint contact interface. Segalman [18] used the parallel series Iwan model in series with a soft element to represent the contact patch and the rest of the joint respectively and showed that this model is a good candidate to capture the nonlinear physics in the contact interface. Considering the softening effect at the contact interface of a bolted lap-joint, Ahmadian and Jalali [19] proposed a lumped joint model consisting of linear and cubic springs and linear damping elements. Yuan et al. [20] used the adjusted Iwan model to simulate the behavior of a bolted joint. Bogard et al. [21] provided a review on various method of modeling the behavior of mechanical joints in assembled structures. Jalali et al. [22] identified the parameters of a Valanis model governing the nonlinear restoring force raised by the micro-slip mechanism at the contact interface of the clamped end of a cantilever beam by using the describing function inversion method. Song et al. [23] developed a general joint interface element to predict the dynamics of structures at which both friction and impact at the contact interfaces cause energy dissipation. In this model, they considered various behaviors of stick, slip, and slap in the contact. Ahmadian and Mohammadali [24] proposed a zero thickness interface model to predict the slip and slap dynamic behaviors of joint contact interfaces in the normal and shear directions of the interface, respectively. A comparison was made by Lacayo et al. [25] between the performance of a time-domain whole-joint approach and a frequency-domain node-to-node approach. They introduced a numerical benchmark in order to assess the efficiency of the two joint modeling approaches in predicting and simulating the dynamics of a mechanical joint. Lingfei et al. [26] introduced an order-reduction modeling approach to estimate the dynamic response of a bolted joint structure. This model takes into account the effects of nonlinearity in the contact interface, damping due to various phenomena and geometrical characteristics of joint structures. Schwingshackl et al. [27] investigated the nonlinear dynamic response of bolted flange joints. They proposed a model with three dimensional elements and verified the capability of the model in regenerating both the nonlinear frequency response and the underlying nonlinear mechanism of the flange joint. The results of their study showed that the flange bolt connection is the main region at which energy is dissipated in a bolted flange joint. Meisami et al. [28] introduced a detailed model of the dynamic behavior of bolted flange joints. This model consists of cantilever beams and springs, respectively representing the joint laps and bolts, and showed that bilinear stiffness is required in modeling the flange joints. The dynamic behavior of joint laps can be modelled using Euler-Bernoulli beam theory. Beaudin and Behdinan [29] suggested a model with analytical parameters related to the flange geometry to consider the nonlinearity due to friction and partial clearance at the contact interface. By introducing this model, they avoided the complexity of previous models, that require experimental testing for model updating or numerous nonlinear parameters.
Many types of mechanical joints, such as bolted lap joints, under platform dampers, blade roots, shroud contacts and bolted flange joints have been investigated by many researchers [29][30][31][32][33][34][35][36][37][38]. However, bolted flange joints are widely used in assembling mechanical structures such as aero-engine casings. An accurate model is required to understand the uncertainties involved in the contact interfaces of a bolted joint. However, many researchers showed that this type of models are numerically expensive and therefore this motivates the search for equivalent models. The development of a reducedorder model based on a detailed model have received less attention by investigators in the past. The numerical model of the joint should be computationally efficient so that the costs of the numerical simulation are minimized. In this paper, the micro-slip behavior at the contact interface between two flange substructures is analyzed using two models. In the first N. Jamia, H. Jalali, J. Taghipour et al.
Mechanical Systems and Signal Processing 153 (2021) 107507 model, the substructures and the contact interface are considered in detail by using 3D solid and point-to-point contact elements and a detailed 3D FE model is constructed in ANSYS. A quasi-static analysis is used to extract the hysteretic behavior. By using this model different behavior in the tangential direction of the contact interface, such as sticking, micro-slip and macro-slip, can be simulated. However, the dynamic analysis of the detailed model of a joint is computationally expensive, hence there is an urgent need for development of equivalent reduced-order models. To this end, a reduced-order model is proposed for the bolted structure in which beam elements and an appropriate joint model consisting of Jenkins elements are used. An identification scheme is proposed and the joint model parameters in the reduced-order model are characterized by using the linear and nonlinear behavior of the contact interface in the detailed model. Good agreement is achieved between the detailed model and the identified reduced-order model.

3D detailed FE model
In this section a detailed three-dimensional finite element model of two bolted rectangular flanges is developed. The FEM analysis is performed using the commercial software ANSYS Workbench.
A bolted structure consisting of three M10 through bolts, three nuts, three washers and two rectangular flanges connected by the bolts as shown in Fig. 1 is considered. Structural steel with the following parameters was assigned as the material for whole model, i.e. bolts, nuts, washers and substructures; E = 200 GPa, m = 0.3 and q = 7850 kg/m 3 , where E is the Young's modulus, q is the mass density and m is the Poisson's ratio.
The dimensions of the bolted structure are given in Fig. 1. For simplicity, the threads are not modeled in this paper. The two L-shaped substructures have the same geometric and material properties.
To provide high accuracy in the contact analysis, particular care was taken to ensure that the pairs of nodes on the two contact surfaces using nonlinear elements are coincident to ensure an accurate interface mesh. The FE model is constructed using a fine mesh in the contact area, bolts, nuts and washers, as shown in Fig. 2(a). A very fine mesh was added near the flange holes to capture the high stress gradients as shown in Fig. 2(b). To model the bolts, a solid bolt simulation type was considered. This type requires the use of solid elements to model the bolt stud, head and nut. Therefore, a minimum of 3 elements in the thickness direction of the flanges, bolt heads, nuts and washers was used to provide accurate results. The circular edges of the bolt's holes and the external edges of the lower and upper surface of the contact interface were divided into a high number of elements to refine the mesh in the contact interface and around the bolt holes. A coarse mesh was used far from the contact interface. Tetrahedral elements were used in the mesh as shown in Fig. 2. The assembly has 17,440 elements and 91,864 nodes.
Three types of contact elements embedded in ANSYS are used in the model; bonded, frictionless and frictional. Frictionless behavior allows the parts to slide relative to one another without any resistance. This type of behavior was assigned to the contact between the bolt-shank and flange holes and between the bolt-shank and the washer. The frictional contact provides shearing forces between the parts in contact and was assigned between the contact bolt head to the top surface of the upper flange, the nut interface to the bottom surface of the lower flange and the interface between the upper flange and the lower flange. Finally a bonded contact was assigned to the contact nut interface to the bolt.  The contact interfaces of the bolted joint are modelled in ANSYS based on the type of the contact between parts. In the normal direction, ANSYS prevents the contacting bodies from interpenetration by enforcing contact compatibility using contact algorithms [39]. For the case of nonlinear solid body contact, the Augmented Lagrange and Pure Penalty formulations are chosen for frictional or frictionless and bonded contacts respectively. These methods combine robustness and high accuracy of the nonlinear contact solution. These two formulations are penalty-based contact formulations. Friction behavior is considered following Coulomb's law F tangential lF normal where m is the coefficient of static friction. Once the tangential force, F tangential , exceeds lF normal sliding will occur. If m is zero then the contact behaves as frictionless. In this model, a coefficient of friction of m=0.2 for steel-on-steel contact was assigned to the contacting area between the two flanges. It is worth mentioning that, in reality, oblique contact can be formed between asperities and causes nonlinear dependence of the friction force on sliding speed. This phenomenon is known as the Stribeck effect. However, in this paper, and for the sake of simplicity in the modeling, Coulomb friction in the tangential direction at each point in the contact interface is considered.
Furthermore, since the goal of this paper is to model a full structure containing the contact interface, it is not possible to model the detailed contact of asperities because the asperities are very small in size and there are some aspect ratio limitations in FE when modeling a detailed part in a structure. Modeling the actual contact of asperities, i.e. oblique contact, is only possible when a small portion of the contact interface is modeled in a macro scale.
Regarding loads and boundary conditions, the lower side of the bottom flange was constrained and a bolt axial load was applied to the three bolts. A pretension in the bolts was generated at the mid-plan of the bolt to initially displace the upper  and lower segment of the bolt stud towards each other and therefore create a preload effect. To specify the pretension in the bolt, a local coordinate system was defined with the Z-axis along the bolt length as shown in Fig. 3(a). Due to the presence of the bolt pretension, the simulation was performed over two load steps. During the first load step, the different assembly interferences are calculated to provide the prescribed preload and reaction forces into the bolts and mimic the assembly of the two flanges. Once the bolt pretension load is achieved, then the calculated reaction force agrees with the applied preload in the bolts. During the second load step, the external load was applied to the bolted joint as shown in Fig. 3(b). All the loading steps are quasi-static where a time integration period of one second is assigned to step 1 and a period of 3 s is assigned to the second load step. This second step is ramped with time. The force convergence is monitored, and convergence is achieved before moving to the next load step. The time response is monitored in all runs to guarantee that the steady state is obtained.

Detailed model simulation results
In this section the effects of bolt preload and external applied static force on the status of the contact interface (i.e. sticking or sliding) and the pressure distribution is studied. The different hysteresis loops obtained from the detailed model are then analyzed. The contact status given in Fig. 4(a) shows the status of the contact interface due to only pretension. The status present in a contact can take four forms: Sticking, where the two contact surfaces cannot move tangentially, Sliding, when there is a tangential movement between the two contact surfaces, Near, when the contact surface has a normal separation within a given radius, and Far, when the contact surface has considerable separation.
The boundary between the near field and far field open status is defined by the pinball region size defined in ANSYS (i.e. pinball radius). In this model, the pinball region radius is equal to 1.0383 mm.
A presence of well distributed sticking behavior is shown around the flange holes due to the closing effect of the preload in the bolts. After an external load is applied in two opposite directions, the stiffness in the joint interface is asymmetric as shown in Fig. 4(b) and (c). We notice the presence of sliding behavior in the contacting area between the two flanges, as shown in Fig. 4(b) and (c). Sliding is caused by the elements that slide on the surface of the contact; however sticking is caused by the contact elements that cannot move. As expected, a symmetric pressure distribution is obtained about the vertical centerlines of the bolts. The resulting contact pressure at the contact interface given in Fig. 5(a) shows a nearly uniform pressure distribution localized around the edges of the flange holes. A maximum of pressure occurs close to flange holes and decreases concentrically to an area of lower pressure between the holes. Due to the opening gap caused by the pretension of the bolts, we observe a zero pressure at the ends of the contact interface. Fig. 5(b) and (c) show the pressure distribution due to the external load applied in the two opposite directions, respectively. The total contact area and the shape of the pressure contact zone is clearly altered due to the presence of external force. A significant concentration around the bolt holes is still observed but a non-uniform distribution is due to the effect of the external force on the stiffness of the contact interface. In the following, the modal properties of the assembled structure are presented. To determine the natural frequencies and mode shapes of the system, an eigenfrequency study was performed using ANSYS Workbench. A fixed constraint was applied to the base of the lower flange. The first 6 natural frequencies and their corresponding mode shapes are given in Table 1 and Fig. 6, respectively.
The interest in this paper is the bending modes, and their corresponding natural frequencies will be compared to the equivalent model given in the next section. Therefore, only modes 1, 2 and 5 will be taken in consideration in this study. In addition, the focus in this paper is only on the slip mechanism and its nonlinearities in the joint interface. The slap mechanism is not considered and therefore no large deflection nonlinearities occur. Thus the scope of this paper is within the Table 1 Simulated natural frequencies (Hz).   linear region of the joint, apart from the interface friction, and the large deflections shown in the different modes in Fig. 6 is only for visualization purpose. Finally, the behavior of the contact interface under static and dynamic loadings is examined. A parametric study was performed where the amplitude of the static load shown in Fig. 3(b) was varied and the relative displacement at a point of the contact interface between the two flanges was monitored. Fig. 7 shows the behavior of the relative displacement at the contact interface between the two flanges under the action of the static external load. A linear behavior occurs for forces up to 1.5 kN and beyond that range a nonlinear behavior occurs. This static study gave an approximation of the value of static force starting from which the structure starts to behave nonlinearly.
A dynamic analysis was also performed where a harmonic distributed force is applied to the bolted joint in the transverse direction of the contact surface over the highlighted area shown in Fig. 8(a) (Note that the arrow in Fig. 8(a) is just an indicator of the force direction). This force is defined as Initially, a first load case was considered where one force was applied to the lower face of the upper flange, close to the contact interface between the two flanges, as shown in Fig. 8(a). This will help to excite the nonlinear behavior of the joint interface and creates a sliding behavior in the contact interface.
A time integration simulation was performed for a period of 4 s with a time step equal to 0.0025 s and the relative displacement of a point from the contact area of the bolted flanges is calculated for each time step. The displacement in the xdirection was collected from a set of points from the upper and lower surface of the contact interface. For the sake of brevity, only two points, P1 and P2 shown in Fig. 8(b), are considered to present the corresponding results in this study. Three different excitation frequencies, namely x = 2 Hz, x = 7 Hz and x = 15 Hz, were considered. Fig. 9 shows the external load versus the relative displacement at the interface contact between the two flanges for a number of loading-unloading cycles for different excitation frequencies. It can be observed that even for low frequency the results are not well converged and the steady state is not obtained; hence a noisy data is shown in Fig. 9(a). The noise in the hysteresis loops is mainly caused by the load applied to the contact interface. Increasing the excitation frequency will increase the opening and closing behavior (i.e. micro-vibro-impacts) at the contact interface and this will lead to more noisy data. In addition we can observe the absence of symmetry of the hysteresis loop at 15 Hz in Fig. 9(c) were the curve is shifted far from the center. The 'drifting' of hysteresis loops is the result of micro-vibro-impacts happening at the contact interface due to opening and closing. This phenomenon was first observed and reported by Ma et al. [1]. In order to eliminate this 'drifting' phenomenon in hysteresis loops the loading should be changed such that the level of micro-vibro-impact decreases. This is achieved by considering the load case shown in Fig. 10.
To improve the convergence of the results, a more refined mesh and a longer period of simulation is needed. However, this will increase the cost of the model in terms of computational time, which was not possible for this paper. In addition, to achieve a robust equivalent model of the detailed model, which will be presented in the next section, the detailed model needs to provide well-converged results. Therefore a second load case was considered in the following part of the paper, where a second harmonic distributed force, F 2 (t), that is equal to the first force, F 1 (t), but with opposite sign was applied to the lower face of the lower flange, as shown in Fig. 10. The choice of this type of load is explained by the need to amplify   the friction effect in the contact interface, which leads to the identification of the micro slip mechanism in the contact interface, which is the main interest in this paper. Hysteresis loop at two points, P1 and P2, are given in Figs. 11 and 12, respectively, for the three different excitation frequencies. It can be observed that the curve shows a hysteresis loop which has a parallelogram shape where two stages can be identified. The first stage occurs when the variation of the relative displacement is very slow relative to the variation of the applied load. This behavior corresponds to sticking. Increasing the applied load leads to gross slip and a microslip occurs during the transition between the stick and gross slip. The behavior of the contact interface is linear in the sticking and macroslip stages and non-linear in the micro-slip stage. Although for macro-slip, the whole contact interface is slipping, the hysteresis loops show a linear stiffness, which is due to the bending resistance of the bolt beam. This loop is used to quantify the energy dissipated in the joint contact due to the damping of the joint. The hysteresis loop in Fig. 12 shows that the detailed model is able to identify the behavior of microslip transition between the elastic contact deformation zone and the full macro-slip zone.
To verify the reliability of the detailed model, a parametric study was performed were the effect of bolt preload, force amplitude and friction coefficient on energy dissipated at the contact interface between the two flanges was depicted in Fig. 13. Fig. 13(a) shows the hysteresis loop obtained from the detailed model at an excitation frequency of x ¼ 2 Hz and a bolt preload of P ¼ 8 kN for different amplitudes of periodic force applied to the flanges. At low force amplitude, the hysteresis loop is almost linear and the contact surface stays in partial slip status. By increasing the force amplitude, the hysteresis loop shape becomes a parallelogram and the macro-slip behavior occurs and therefore the energy dissipated increases. In addition, the slope of the different loops remains constant since the pretension is the same in this case. Fig. 13(b) shows the hysteresis loops at an excitation frequency of x ¼ 2 Hz and a force amplitude of F ¼ 4 kN for different bolt preloads. At low pretension, the friction forces are small and therefore the contact is sliding and the load remains constant. Increasing the preload leads to an increase in the friction forces and therefore a less slipping occurs between the contact surfaces. Hence, the higher the pretension in the bolt, the lower the energy is dissipated in the contact interface. Fig. 13(c) shows the hysteresis loop at an excitation frequency of x ¼ 2 Hz, a bolt preload of P ¼ 5 kN and a force amplitude of F ¼ 8 kN for different friction coefficients in the contact surface between the two flanges. The energy dissipated decreases with an increase of the friction coefficient as expected since increasing friction coefficient prevents slippage at the contact interface and hence less energy will be dissipated. Finally, Fig. 13(d) shows the hysteresis loop at excitation frequencies of ¼ 7 Hz , x ¼ 15 Hz, x ¼ 20 Hz and x ¼ 25 Hz for P ¼ 5 kN and F ¼ 4 kN. The effect of the variation of the frequency of excitation on the energy dissipated is very small. In addition, it is shown in Fig. 13 that the maximum transmitted force is the same for all hysteresis loops since a constant external force is applied in all cases.
In general, by increasing the friction coefficient, the contact interface undergoes micro and macro slip conditions at a higher excitation level and also the amount of energy dissipated in the contact interface decreases. The effect of the friction coefficient would be to expand or shrink the hysteresis loops, as Fig. 13c shows, and not on the shape of the hysteresis loops. Therefore, the equivalent model and parameter identification approach proposed in this paper in the next sections can be effectively used for different friction coefficients. The results obtained from the detailed model under load case 2 suggests this model is a good test case to model the microslip and macroslip behavior with an FE approach. The main drawback of this type of model is the high computational cost where a run of the detailed model takes up to 14 h for a time step equal to 0.0025 s and a low excitation frequency of x ¼ 2 Hz. Therefore, an equivalent model of the 3D detailed model is formulated in the next section, which reduces the run time to less than 2 mins. 4. 1D equivalent model

Equivalent model formulations
The detailed model presented in the previous section is not suitable for dynamic response analysis under general loading conditions because of its large number of DOFs and the complexity of the contact elements used. Although, the search algorithms used in contact elements allow for an update of the contacting surfaces pairing, they increase the solution runtime drastically. Therefore, constructing equivalent reduced-order models for structures containing frictional contact interfaces has long been desirable. The challenge of this task is to find appropriate reduced-order models for contact surfaces and formulate reliable identification schemes to characterize the models used. In the following a reduced-order model is constructed for the structure whose detailed model was described in the previous section.
A representation of the reduced-order model of the bolted flange structure is shown in Fig. 14. The bolted flange structure is composed of two L-shaped plates which are connected using three bolts. Since the bending behavior of this structure is of interest, the L-shaped plates are reduced in the equivalent model to Timoshenko beam elements [40,41]. Timoshenko beam elements are placed in the mid-plane of the plate sections, and rigid constraints (or MPCs) are used to relate the DOFs in each L-shaped beam section. Before considering the modeling of the contact interface, the accuracy of the L-shaped sub-structures in the reduced-order model is evaluated by comparing their natural frequencies with those of their solid model counterparts and a good agreement is obtained. Finally, appropriate joint models have to be considered to model the behavior of the contact interface in normal (i.e. z) and tangential (i.e. x) directions. The joint model used in this paper is introduced in the following. There are two well-known mechanisms that develop in a contact interface: micro/macro slip in the tangential direction and micro-vibro-impacts in the normal direction [42]. The model used to consider the slip mechanism in the tangential direction is described first. From a microscopic point of view, a contact interface is composed of a number of discontinuous contact points, as shown in Fig. 15(a). The behavior of each of these contact points is linear under small loading conditions due to local microscopic elastic deformations. When the applied load is increased beyond a certain level, such that the portion of the load carried by a contact point exceeds the critical slip force, the contact point starts to slip. Therefore the behavior of a single contact point can initially be modeled as a combination of a linear spring k 1 in series with a frictional slider element as shown in Fig. 15(b). It is worth mentioning that there could be microscopic plastic deformation and wear in real contact interfaces prior to the macro slip, which are disregarded for the sake of simplicity in the analysis in this paper. The force-displacement behavior of a contact point remains linear even after slippage occurs at that point. This is mainly due to the linear behavior of the other contacting points which are in the pre-slippage stage. This linear behavior is also physically justifiable after all contact points start to slip, i.e. the macro-slip stage of the contact interface, in the case of a bolted joint.
In bolt jointed structure this post slip elastic behavior is due to the presence of the bolt [43]. When all points at the contact interface start to slip -i.e. the state of macro-slip is reached -and the relative displacement of the upper and lower surfaces exceeds the bolt hole and bolt beam clearance, the bolt beam resists further relative movement. Therefore, a linear force-displacement behavior is observed due to linear bending of the bolt beam. This linear post slip behavior of a contact point can be modeled as a linear spring element k 2 parallel to the combination of the linear spring and frictional slider elements to form the well-known Jenkin's element shown in Fig. 15(b) [44]. The force-displacement behavior of the Jenkin's  element is shown in Fig. 15(c). Therefore, a number of Jenkin's elements representing the behavior of contact points can be used to model the tangential behavior in a joint interface. These Jenkin's elements must have different critical slip forces, i.e. F s . One of the challenges in modeling contact interfaces is the identification of the distribution of the critical slip forces over the contact interface. In the Iwan model, the number of sliders is considered infinite and a probability distribution function is identified for the critical slip forces [45,12]. In this paper, a finite number is considered and a method is proposed to identify the distribution of the critical slip forces based on data from the detailed model. The force displacement relationship for the Jenkin's element can be expressed as, where u 1 and u 2 in Fig. 15(b) and in Eq. (1) represent respectively the tangential movement of the upper and lower plates at the contact interface as shown in Fig. 16(a). F 1 and F 2 are forces in spring elements k 1 and k 2 and, If the slider is not sliding but has slid before, then where u s is the amount of sliding. The behavior of the contact interface in the normal direction is considered linear in this paper. In other words, it is assumed that, for a specific range of frequency excitation and applied force, the nonlinear mechanism is such that microvibro-impacts do not occur in the contact interface in the normal direction. Fig. 16(a) shows a schematic of the contact interface modeling. The linear behavior in the normal direction is modeled using a linear spring k n , as shown in Fig. 16(b).
The equation governing the dynamic response of the reduced-order beam model introduced in Fig. 14 can be expressed as, where M ½ and K ½ are the mass and stiffness matrices corresponding to the beam sections of the structure modeled using beam elements, d f g and F t ð Þ f gare the vectors of global structural displacement and force fields, K À k 1 ; k 2 ; k n ð Þ h i is the contribution of the contact interface to the global stiffness matrix and F À F si ð Þ n o is a nonlinear force vector which is a function of the distribution of the critical slip forces. In this paper, it is assumed that all the Jenkin's elements have the same stiffness properties in the normal and tangential directions, but they differ in their critical slip forces. It is worth mentioning that the critical slip force, F si , has its maximum close to the bolt and it reduces in regions far from the bolt. This is due to the distribution of the normal pressure in the contact interface. This fact will be used in the identification of a distribution function governing the critical slip forces in next section. In Eq. (5) the parameters associated with contact interface, i.e. k 1 , k 2 , k n and the distribution of F si , must be identified. In the next section an identification scheme is introduced for the contact interface parameters.

Parameter identification approach
When a tangential force is applied to the contact interface, as is shown in Fig. 16(a), the force-displacement curve shows a softening effect which is common in bolted structures (Fig. 16(c)). In fact, at low levels of applied force, where all of the contact points are in the pre-slippage stage, the behavior of the contact interface is linear. As the level of applied force, F, is increased, the contact points start to slide, and the stiffness of the contact interface gradually decreases resulting in a nonlinear behavior. This stage is known as micro-slip. At macro-slip stage and when the whole contact interface starts to slide, the force-displacement curve is again linear due to the elastic shear deformation of the bolt. The linear parameters of the contact interface can be identified from the linear parts of the force-displacement curve. The non-linear portion of the force-displacement curve is used to estimate the distribution of the critical slip forces, as described in the following. It is worth mentioning that similar to the force-displacement curve, the hysteresis loops contain two linear parts and one nonlinear part which can be used equivalently for parameter identification.
The hysteresis loops of the bolted flanges are obtained from the detailed model and m Jenkin's element are used to model the contact interface, as shown in Fig. 16(b). The variable m was chosen based on a convergence study by increasing the number of elements until there was no significant change in the hysteresis loops. The slopes of the force-displacement curve in the two linear portions are used to determine k 1 and k 2 . Since at low level forces (i.e. before any slippage takes place) both of the springs k 1 and k 2 contribute to the tangential stiffness of the contact interface, the total tangential stiffness in the reduced-order model is m k 1 þ k 2 ð Þ . Similarly, at high levels of force and after all of the slider elements reach their critical slip forces, only the k 2 springs are active and therefore the total tangential stiffness of the reduced-order model is mk 2 . By comparing these stiffnesses, i.e. m k 1 þ k 2 ð Þand mk 2 , with the slopes of the two linear portions in force-displacement curve in Fig. 16 one obtains k 1 ¼ K 1 m and k 2 ¼ K 2 m : After identifying the linear tangential stiffness coefficients of the contact interface, there is one more linear stiffness in the normal direction left to be identified. This is done by comparing the natural frequencies obtained from the detailed model and the ones obtained from the reduced-order model. Since k 1 and k 2 are known, the equation governing the linear free vibration behavior of the reduced-order model is, By considering a harmonic response of the form d t ð Þ f g¼ X f gsin x n t ð Þin Eq. (6), the natural frequencies are obtained as a function of k n , the contact interface normal stiffness. k n can be determined by minimizing the difference between the actual and predicted natural frequencies for the first bending mode using a sensitivity-based identification scheme [46].
Next, the identification of the distribution of the critical slip forces F si is considered. First the minimum and maximum values of F si are determined. As described earlier in this section, at the end of the first linear part in force-displacement curve the first slider has the minimum critical slip force and slides. Since slider i slides when F si ¼ k 1 Du, therefore the minimum critical slip force in the contact interface is obtained as F s ð Þ min ¼ k 1 Du 1 . The same logic can be used to obtain the maximum critical slip force in the contact interface. The last slider with maximum critical slip force slides at the beginning of the second linear part in the force-displacement curve. Therefore, the maximum critical slip force is obtained as F s ð Þ max ¼ k 1 Du 2 . Du 1 and Du 2 are the relative displacement in the contact interface corresponding respectively to the end of the first linear part and beginning of the second linear part, as shown in Fig. 17(a).
Knowing the maximum and minimum values of the critical slip force, the distribution of the critical slip forces needs to be identified. As mentioned before, the distribution of the critical slip force in the contact interface is a function of the contact pressure in the normal direction. The normal contact pressure is high close to the bolt and gradually decreases to reach its minimum value at the furthest point from the bolt. Therefore it is assumed that the maximum and minimum critical slip forces correspond respectively to the right and left corners of the contact interface, as shown in Fig. 17(b). In addition, it is assumed that the critical slip force changes exponentially over the contact interface as (Fig. 17(b)), where z is a distance parameter and a and b are obtained by substituting the boundary values into Eq. (7) (i.e. @z ¼ 0; F s z ð Þ ¼ ðF s Þ min and @z ¼ L; F s z ð Þ ¼ ðF s Þ max ) and are defined as follows, where L is the length of the contact interface and c is an unknown parameter to be identified. c and k are tuned in the next section so that a good agreement is obtained between the hysteresis loops from the detailed and equivalent models. The form of the function in Eq. (7) is chosen such that it provides a smooth transition from minimum to maximum values of Fs over the contact interface, and also contains enough unknown parameters, which allows a good fit between hysteresis loops in the identification process. The equivalent model hysteresis loops are obtained by solving Eq. (5). The identification results are provided in the next section. It is worth mentioning that the distribution of critical slip forces considered in Eq. (7) is directly related to the contact pressure distribution in bolted joints. Many researchers have investigated the contact pressure distribution in the past by using experimental or numerical methods. The contact pressure distribution is a function of different parameters such as material properties, bolt size and hole clearances. Different function such as exponential [47] and polynomial [48] functions have been used to represent the distribution of contact pressure in bolted joints.

Parameter identification results
The identification of the equivalent model parameters is performed using the hysteresis loops shown in Fig. 13. Based on these hysteresis loops, it is observed that the stiffness of the contact interface in the initial linear behavior, i.e. K 1 þ K 2 , and in the final linear behavior, i.e. K 2 , are the same for all bolt pre-tensions, force excitation amplitudes and excitation frequencies.
As explained in the previous section, in the final linear behavior, the contact interface is in the macro-slip condition. It is worth mentioning that reaching the macro-slip state in the contact interface depends on several factors such as bolt pretension, applied force amplitude and excitation frequency. For specific combinations of bolt pretension, applied force amplitude and excitation frequency, the macro-slip at the contact interface does not occur. This can be seen for the case where F ¼ 4kN, P ¼ 8kN and x ¼ 2Hz, as shown in Fig. 13(b). Fig. 13(b) also shows that the effect of bolt pre-tension is more significant in the micro-slip stage, and at high bolt pre-tension, the contact interface starts to slip at higher relative displacement. This means that different distributions of the critical slip forces in the equivalent model are expected to be obtained at different bolt-pretensions. The same approach is applied to the hysteresis loops at different excitation frequencies (Fig. 13  (d)). The hysteresis loops in Fig. 13 are first used to determine the linear tangential stiffness of the contact interface in the equivalent model, i.e. k 1 and k 2 , by following the method described in the previous section. The identified stiffness coefficients for the tangential direction of the contact interface are presented in Table 2.
The contact interface in the equivalent model was modeled using m Jenkin's elements, where m ¼ 51. Knowing the stiffness coefficients in the tangential direction, the normal contact stiffness of the equivalent model, i.e. k n , is identified by updating the linear part of the equivalent model using the first three bending natural frequencies as described in previous section. The comparison between detailed and equivalent model natural frequencies is given in Table 3.
The results shown in Table 3 indicate that the equivalent model can predict the natural frequencies of the detailed model with an acceptable accuracy. Regarding the identification of the distribution of the critical slip forces in the contact interface     of the equivalent model, the backbone or force-displacement curves shown in Fig. 13(a) and (b) was used by employing the identification scheme proposed in the previous section. The identified equivalent model is then used to predict the hysteresis loops. The identified distribution of the critical slip forces and the comparison between the simulated and predicted hysteresis loops are now presented.  19(a) show that the identified distributions for critical slip forces at a constant bolt pre-tension or a constant force amplitude follow a similar pattern. Figs. 18(b) and 19(b) show that the identified distributions of the critical slip forces for the equivalent model can predict the simulated hysteresis loops obtained by the detailed model with a high accuracy.
In order to evaluate the predictability of the equivalent model, the prediction of the hysteresis loops at a different excitation frequency are considered. Fig. 20 shows the identification of the critical slip force distribution for different excitation frequencies at a constant force amplitude of F = 4 kN and bolt pre-tension of P = 5 kN. The similarity of the different curves  shows that, if the excitation amplitude and the bolt pre-tension are constant, then the same force distribution is identified for the slip forces for different excitation frequencies. In addition, the identified distribution of slip force at excitation frequency of x ¼ 2Hz may be used to predict the hysteresis loops for other excitation frequencies, as shown in Fig. 21.
The results presented in Figs. [18][19][20][21] show that the proposed method for modeling the contact interface of bolted flanges effectively captures the physics and predicts the dynamic response of the assembled structure. It is worth mentioning that since the equivalent model is 2D and also has many fewer DOFs than the detailed model, especially at the contact interface, a detailed comparison of the contact interface behavior for the two models is not achievable.
The detailed model presented in this paper successfully simulates the dynamic response of an assembled structure containing frictional contact interfaces. Therefore, the equivalent model, which represents the detailed model with acceptable accuracy, can be used in real applications including a frictional contact interface. Two main advantages of the equivalent model are that it avoids expensive computational FE analysis for the assembled structure, and it can be used in different loading conditions when its parameters are appropriately identified in one loading condition. The observed nonlinear behavior reported in the paper is inherent in the frictional contact interface. This means that for any assembled structure with a frictional contact interface a similar dynamic behavior should be observed. Therefore, the equivalent modeling method and its parameter identification approach presented in this paper can be effectively used in constructing reduced-order models for assembled structures in real application.

Conclusion
In this paper, a detailed model of bolted flanges was constructed to capture the friction behavior in the contact interface due to the micro-slip phenomena. Static and dynamic loads were applied to the joint and the force-displacement curves and hysteresis loops were obtained for the contact interface. Since detailed modeling is often not possible in real engineering applications, an equivalent model was constructed by using beam elements and a proposed joint model consisting of Jenkin's elements. Based on an identification scheme, the linear and nonlinear parameters of the equivalent model were identified using the results from the detailed modeling approach. The equivalent model showed a good ability to predict the dynamic response of the assembled bolted flange. It was shown that the equivalent model identified at an excitation frequency of x ¼ 2Hz, can predict the hysteresis loop for other excitation frequencies. This predictability characteristic of the equivalent model makes it applicable in real loading conditions where the excitation force is usually a combination of different frequencies. Also, by using the equivalent modeling approach, it is possible to establish a relationship between the joint geometry and the equivalent model parameters which is a great help in using bolted joint models in a real design application. In this paper the contact mechanism in the normal direction was considered linear and the equivalent model parameters were identified using the simulation results. Considering the nonlinear slap mechanism in modeling the contact interface in the normal direction and identifying the model parameters by using experimental results will be considered in future work.

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