Estimation Modal Parameter Variation with Respect to Internal Energy Variation Based on the Iwan Model

: Typical factors that cause nonlinear behavior in structures are geometric nonlinearity, force and displacement boundary condition nonlinearities, and material nonlinearity. The nonlinearity caused by an increase of the internal energy in built-up structures is mostly due to the displacement boundary condition induced by the contact interface region. This study proposes an experimental mode analysis technique that predicts changes in natural frequencies and damping ratios when the external excitation force increases in a structure’s contact surfaces. Speciﬁcally, the nonlinearity of the dynamic characteristics induced by the contact region is described by the constitutive Iwan model. Next, an estimation method was developed for two parameters among the four of the Iwan model. This study used a modal analysis method. As an extension of a previous study, the approximate form of the harmonic excitation-induced force was determined in closed form. The conﬁguration of the numerical model for the full structure was introduced from this resultant form. By using these numerical results, responses in the full structure, according to the harmonic excitation, have been represented in mode summation form. This research proposes an estimation method for two parameters among the four of the constitutive model. The proposed method was veriﬁed by simulations conducted with the lumped model and by experiments conducted on a partially connected double beam.


Introduction
Generally, experimental modal analysis was performed to ensure structural stability under operating conditions. The structure's stability under dynamic conditions was evaluated based upon the identified modal parameters, such as the natural frequency, mode shape, and damping ratio. However, most mechanical structures do not operate under steady-state conditions; moreover, the operating frequency varies with the magnitude of the excitation force [1,2]. In vibration analysis, an increase in the external force is interpreted as an increase in the internal energy, which changes the dynamic characteristics [3]. Finally, a nonlinear analysis of the structure is required. Nonlinearities are typically divided into geometric nonlinearities [4], force and displacement boundary conditions [5], and material nonlinearities [6]. The nonlinearity due to an increase in the internal energy of the general built-up structure is mostly due to a displacement boundary condition induced by the contact region inside the structure [7,8]. The study proposes a methodology for predicting changes in the natural frequency and damping ratio when the external force increases in a built-up structure with contact surfaces. Specifically, the frequency response function (FRF) in the modal domain is derived based on the Iwan model, which is a representative constitutive model describing the nonlinearity of the dynamic characteristics due to the contact surfaces. Based on the derived FRF, the natural frequencies and damping ratios at a certain magnitude of excitation force or internal energy are predicted.  Figure 1 shows the parallel-series Iwan model that has an equivalent stiffness k. The mathematical description of this model is as follows: In Equation (1), F denotes the applied force. φ denotes the sliding distance of each element and ρ(φ) denotes the population density of Jenkins element of strength φ. u denotes the imposed distance and x denotes the slider displacement. Each quantity has the following dimensions. F has the dimension of force (N). φ, u and x have the dimension of length (m). The slider displacement x(t,φ) evolves from the imposed displacement u(t). .
x and . u represent the slider velocity and the imposed velocity, respectively. If u is subject to frequency ω and amplitude A, it means Then, Equation (1) can be represented as follows: The second term on the right-hand side of Equation (4) can be obtained because x(t, φ) = 0 for all φ > A. If one re-arranges Equation (4), where K T = ∞ 0 ρ(φ)dφ, ∆F S (t) = A 0 ρ(φ)x(t, φ)dφ and u(t) = A sin(ωt − θ). In Equation (5), K T represents the added stiffness by joints in built-up structures. F S represents the variation force for the occurring slip. We adopt the representation of the density function ρ in Reference [12] to investigate in more detail the result in this equation; this representation is composed of the 4 parameters and is shown below: In Equation (6), R, χ, φ max and S are four parameters that describe the density function. Here, R is the parameter related to the magnitude of the density function. χ is a constant, which has a value between −1 and 0. φ max is a constant representing the border line that determines the micro-slip and macro-slip in the contact surfaces of two objects. S is a constant to describe the dynamic characteristics when macro-slip occurs. In addition, H and δ in this equation represent the Heaviside function and the Dirac delta function, respectively. These four parameters {R, χ, φ max and S} can be transformed into the following four parameters {χ, K T , F S and β}. In these newly introduced parameters, χ is the existing Appl. Sci. 2019, 9, 4290 4 of 22 parameter in the original set and K T represents the introduced added stiffness by joint in Equation (4). The other remaining parameters have the following relations [12]: In the results of Equation. (7), F S represents the critical force to generate macro-slip. β is the constant related to the curvature in the backbone curve. If one assumes that the constitutive model in Equation (3) behaves in the micro-slip region, then the variation force in Equation. (5) can be represented as the following: In Equation (8), the slider displacement x is related to the imposed velocity . u as follows [9,13]: Substituting Equation (9) in Equation (8), then, In Equation (10), the superscripts + and -refer to the moment when . u has positive and negative values, respectively. If one substitutes Equation (3) into the results of Equation (10), then, where . If one attempts to represent the variation force for one period, as shown in Equation (11), which is divided into two regions in the time domain by a Fourier series, the following is obtained: where Equation (12) shows the variation force, which is represented in the Fourier series. Equation (13) shows the forms of every coefficient of this equation. In Equation (13), T denotes one period of time (i.e., T = 2π/ω). If one applies Equation (11) to Equation (12) and represents the variation force using only the fundamental frequency component (i.e., the first harmonic component), then the following result can be obtained (a detailed description on how to obtain this result is in Appendix A): Configurations of κ 1 and κ 2 in Equation (14) refer to Equation (A15). If one considers the more practical case of the imposed displacement compared to the case in Equation (3), The variation force with respect to the case in Equation (15) can be obtained as follows if one follows a similar procedure to obtain the result in Equation (14).
Re-representing Equation (16) after substituting Equation (15) into this equation, produces the following: It can be confirmed from Equation (A15) that κ 1 and κ 2 will have negative and positive values, respectively. Therefore, one can re-represent Equation (17) as follows: The above equation approximately represents the induced force when harmonic dynamic behavior occurs in the constitutive model.  Figure 1 shows the parallel-series Iwan model that has an equivalent stiffness k. The mathematical description of this model is as follows: In Equation (1), F denotes the applied force. denotes the sliding distance of each element and ( ) denotes the population density of Jenkins element of strength . u denotes the imposed distance and

Modal Formulation of the Full Structure
In this section, the configuration of the variation force, which is generated in the constitutive model during dynamic behavior, is examined through a 2 degrees-of-freedom (DOF) system. In Appl. Sci. 2019, 9, 4290 6 of 22 addition, this result is extended to a multi-degree-of-freedom system, which involves the constitutive model. The configuration of the FRF for this system with respect to the harmonic excitation is obtained. The methodology, which is intended in this study to estimate certain parameters out of the parameters of the constitutive model, is introduced through these results. Figure 2 represents a 2 DOF system, which is connected by the stiffness k ij and the Iwan model. Here, u i and u j represent the displacement at each node. The relative displacement between these two nodes when an external force is applied can be represented as follows: In Equation (19), A ij denotes the amplitude of the relative displacement between the i th and j th nodes. By using this equation, the restoring force in the stiffness and the generating force in the Iwan model can be represented in vector-matrix form using Equations (5) and (18): where Figure 2. Description of the 2 DOF system connected by the stiffness k and the Iwan model between the i th and j th nodes. Figure 2 represents a 2 DOF system, which is connected by the stiffness k ij and the Iwan model. Here, u i and u j represent the displacement at each node. The relative displacement between these two nodes when an external force is applied can be represented as follows: In Equation (19), denotes the amplitude of the relative displacement between the i th and j th nodes. By using this equation, the restoring force in the stiffness and the generating force in the Iwan model can be represented in vector-matrix form using Equations (5)  (22). In this equation, the matrix K ij T denotes a matrix that is composed of combinations of the added stiffness (K T ) by joints, which is one of the parameters defining the Iwan model. The matrix A ij denotes a matrix that is composed of the combinations of relative displacement between nodes. Actually, the results in these equations can be extended without a loss of generality for general built-up structures. This means that the results of Equations (20)-(22) can be regarded as the local matrix for the interface between the jointed surfaces of general built-up structures. If one represents the equilibrium of force in general built-up structures by a vector summation form according to the Newton's second law, then, In Equation.
(23), the subscripts I, e, c, r and J denote the force by inertia, the external force, the force by damping, the restoring force by stiffness and the force occurring in the joint component, respectively.
If the equation of motion for the full structure is represented with this equation, then the following is obtained: On the left side of Equation (24), M, C and K denote the mass, damping and stiffness matrices for the full system, respectively. Note that, in general, there are multiple numbers of joint components in the built-up structures. Therefore, as mentioned in [13], "if one wants to estimate information on a certain constitutive model for a certain joint, the other joint components should be excepted from the full built-up structures apart from the joint component of interest". Hence, the numerical model for the full structure, which includes the joint component shown in Equation (24), indicates that only the joint component of interest is connected in the full system. In this study, the full system stands for the same meaning of this statement. On the right side of Equation (24), every term apart from the external forcing vector f e represents forces generated from the joint in a global coordinate system. In particular, matrices K T and A are all of the symmetric matrices and are composed of real values. Re-arranging Equation (24) after transferring all of the terms in the right hand side of the equation apart from the term in the external forcing vector to the left side produces the following: where Appl. Sci. 2019, 9, x FOR PEER REVIEW 6 of 23 The above equation approximately represents the induced force when harmonic dynamic behavior occurs in the constitutive model.

Modal Formulation of the Full Structure
In this section, the configuration of the variation force, which is generated in the constitutive model during dynamic behavior, is examined through a 2 degrees-of-freedom (DOF) system. In addition, this result is extended to a multi-degree-of-freedom system, which involves the constitutive model. The configuration of the FRF for this system with respect to the harmonic excitation is obtained. The methodology, which is intended in this study to estimate certain parameters out of the parameters of the constitutive model, is introduced through these results.  Figure 2 represents a 2 DOF system, which is connected by the stiffness k ij and the Iwan model. Here, u i and u j represent the displacement at each node. The relative displacement between these two nodes when an external force is applied can be represented as follows: In Equation (19), denotes the amplitude of the relative displacement between the i th and j th nodes. By using this equation, the restoring force in the stiffness and the generating force in the Iwan model can be represented in vector-matrix form using Equations (5)  K J denotes the stiffness matrix in a static state and both |κ 1 | ω A and κ 2 A denote the effects in the joint component by the dynamic behavior of the system. It can be confirmed through the signs of these two matrices that each term of these two effects of dynamic behavior will cause a decrease of stiffness and an increase of damping in the full system. These confirmations are consistent with previous research [13,17]. In addition, assuming an external force with a harmonic component, then the velocity vector can be expressed by the multiplication of jω with the displacement vector. Here, j stands for the imaginary value. From this assumption, From Equation (26), it is also confirmed that the increased damping due to dynamic behavior shows the form of structural damping (hysteric damping). It is a natural consequence because the hysteric behavior in elastoplastic elements is described in bi-linear form in Reference [9]. In conclusion, the physical meaning of Equations (25) and (26) are representative of the full system; if the joint component is connected, then the original stiffness K is increased by K J . In addition, if dynamic behavior is initiated by external forces, then slip occurs at the interface of the jointed structures and results in an increase in damping and decrease in stiffness. The type of damping variation during dynamic behavior shows the form of structural damping. Hereafter, the frequency response function will be obtained for the full system. To this end, the eigenvalue problem of Equation (26) is solved. For the case of no external force, there are no dynamic effects (κ 1 ,κ 2 ), and therefore, this equation can be transformed into a simpler form. First, for the un-damped case, the n th eigenvalue and eigenvector can be represented by ω Jn and ϕ n , respectively. These two results are composed of all real values.
Eigenvalue : ω Jn Eigenvector : ϕ n where ϕ n = ϕ 1 n ϕ 2 n · · · ϕ N n T . Equation (27) shows the configuration of the n th eigenvalue and eigenvector. Here, the eigenvector represents the mass normalized eigenvector. The superscript T stands for the transpose. If assuming proportional damping (e.g., C = aM + bK J ; here, a and b are real constants), then the eigenvector will have the equivalent to the eigenvector in Equation (27), and the eigenvalue will have following form; If the effect of stiffness in the static state is much larger than the effect of the variation stiffness by dynamic behavior, it is suggested that the following conditions are satisfied, one can obtain the following results: where, l = 1, N.
From the results in Equation (30), Equation (26) can be represented by the model summation form.
In Equation (31), the term in the bracket represents the frequency response function of a built-up structure that includes the joint component. In the term, the values of ∆ω 2 n and ∆c n in the denominator represent a real and imaginary part, respectively, and indicate the amount of the decreased natural frequency and increased damping by slip at the interface of jointed structures during dynamic behavior. The amount of the damping variation during dynamic behavior will be investigated in more detail. The amount of this quantity for a certain mode is In Equation (32), P represents a quadratic polynomial with variables of the l th eigenvector. M denotes the number of the Iwan model. A represents the magnitude of the relative displacement between the connected nodes by the constitutive model. If the magnitude of the initially applied force was assumed to be f e 0 and the experiment was performed by increasing the applied force z times, then the relative displacement A also increased z times. By using this assumption, the variation of damping in each experiment could be expressed as: where B = M m=1 P m (ϕ l )A χ+1 m 0 . Equation (33) shows the variation in the magnitude of damping values between each sequential experiment when the magnitude of the external force was increased z times based on f e 0 . The total number of experiments shown in Equation (33) was ν+1. However, the variation of damping values in each experiment shown in Equation (33) was difficult to identify practically. Therefore, the relative variation of damping values between sequential experiments was investigated.
Taking the logarithm of the terms on both sides of Equation (34) and using z as the logarithm base, the following was obtained: Equation (35) shows the difference in the damping values of sequential experiments in which the external force was increased z times in each trial. The results of this equation were obtained experimentally. If the results are plotted after taking the logarithmic function, which has z as the logarithm base, then the slope of this plot is χ + 1. Here, the abscissa indicated a negative one order for the performed experiments, and the ordinate indicated the measured values. This concept was similar to the conventional method used to estimate the value of χ by converting the dissipated energy with respect to the external harmonic force to a logarithmic scale. To estimate χ using Equation (35), linear curve fitting was performed on the measured damping values converted to a logarithmic scale. Here, the y intercept was regarded as a 0 . Since z could be measured through experiments or as the magnitude ratio of the applied force in successive experiments, |κ 1 |B was estimated through a 0 in Equation (35). This was the same as ∆c l0 in Equation (33). As a result, every ∆c lv in Equation (33) could be estimated.
By combining Equations. (33), (34), and (A15), the following relations were obtained: In Equation (36), all of the variables except the constant B could be estimated by experiments. If a numerical model for the full system exists, then the stiffness K in Equations (24) and (25) can be estimated, which was the case for the absence of joint components. If an experiment was performed by applying a very small force to the built-up structure, which includes the joint component, then the added stiffness K T at the interface could be estimated. Eventually, a static stiffness K J could be obtained. From this result, the mode shape for an undamped system could be estimated. In addition, if the external forces used in the experiments were applied to the numerical model to obtain the results in Equation (33), then the relative displacement at the contact surface could be estimated. Using these results, Equation (36) was re-arranged: The superscript tilde in Equation (37) indicates estimated values from the numerical model. This equation indicates that if B was estimated using the mode shape and the relative displacements at the contact surface, which were obtained from the numerical model, then the parameter R could be estimated.

Degree-of-Freedom (DOF) Lumped Model
In this section, a simulation is conducted to verify the proposed method from the previous section. Figure 3 shows the configuration of the 10 DOF lumped model that is used in the simulation. Here, m and k represent the lumped mass and stiffness, respectively. In addition, u represents the displacement of each mass. As shown in the figure, two numbers of the Iwan model are assigned (one is connected between m 1 and m 6 and another is connected between m 2 and m 7 ). These two models show their dynamic characteristics according to the relative displacement of the connected nodes in the u direction. The frequency response function for this system is obtained by applying a force f at the lumped mass m 5 , and we define the step sine as the applied force. In addition, a total of twelve experiments were conducted in this simulation, and in each experiment, the magnitude of the applied force follows the step sine, which increases consistently. Detailed introductions and specifications on the system and the performed experiments are as follows: The four parameters on Iwan model include the following: Appl. Sci. 2019, 9, x FOR PEER REVIEW 10 of 23 estimated through a0 in Equation (35). This was the same as in Equation (33). As a result, every in Equation (33) could be estimated. By combining Equations. (33), (34), and (A15), the following relations were obtained: In Equation (36), all of the variables except the constant B could be estimated by experiments. If a numerical model for the full system exists, then the stiffness in Equations (24) and (25) can be estimated, which was the case for the absence of joint components. If an experiment was performed by applying a very small force to the built-up structure, which includes the joint component, then the added stiffness at the interface could be estimated. Eventually, a static stiffness could be obtained. From this result, the mode shape for an undamped system could be estimated. In addition, if the external forces used in the experiments were applied to the numerical model to obtain the results in Equation (33), then the relative displacement at the contact surface could be estimated. Using these results, Equation (36) was re-arranged: where = ( ) The superscript tilde in Equation (37) indicates estimated values from the numerical model. This equation indicates that if was estimated using the mode shape and the relative displacements at the contact surface, which were obtained from the numerical model, then the parameter could be estimated. In this section, a simulation is conducted to verify the proposed method from the previous section. Figure 3 shows the configuration of the 10 DOF lumped model that is used in the simulation. Here, m and k represent the lumped mass and stiffness, respectively. In addition, u represents the displacement of each mass. As shown in the figure, two numbers of the Iwan model are assigned (one is connected between m1 and m6 and another is connected between m2 and m7). These two models show their dynamic characteristics according to the relative displacement of the connected nodes in the u direction. The frequency response function for this system is obtained by applying a Moreover, damping is assigned to increase the accuracy of the measurements made by the simulation. Proportional damping of the mass and stiffness matrices is assumed. As mentioned in the previous section, the stiffness matrix involved in this proportional damping is related to the stiffness K J introduced in Equation (25):

Degree-of-Freedom (DOF) Lumped Model
In Equation (40), a and b are constants, and in this study, values of 5 × e −1 and 7 × e −5 are assigned as these constants, respectively. In addition, the frequency resolution in the step sine experiment is 0.5 Hz and the time integral in the simulation is performed by the Newmark method [18]. Figure 4 shows the obtained FRFs from performed simulations. As shown in this figure, it can be confirmed that each mode has a different effect with respect to the magnitude of the external force. For instance, some modes (i.e., 2 nd , 4 th , 6 th and 8 th modes) show that their natural frequencies are shifted and that their peak values are decreased according to the increased magnitude of the force, whereas the other modes (i.e., 1 st , 3 rd , 5 th , 7 th and 9 th modes) hardly change their natural frequencies or their peak values. This is because the effect of the Iwan model is activated when relative displacement occurs in the region of interest. However, the modes that show no variations in the FRF results with respect to the changes of the external force all have in-phase motion in the interface between the masses of the Iwan model. The results shown in this figure are consistent to the statements in Ref. [13]. In addition, this figure magnifies the results for the 2 nd , 4 th and 6 th modes where variations in the frequency results occur with respect to the magnitude of the external force. From these magnified plots, changes of damping and the natural frequency can be identified in each experiment. To identify the parameter χ with the proposed method in this study, the results of Equation (35) for these three modes have been identified. In this process, the modal parameters for these modes have to be estimated. Recently, polyreference least-squares complex frequency-domain method which is so-called PolyMAX [19] is introduced as a modal parameter estimation method conducted in the frequency domain. It is one of the experimental modal parameter estimation methods and it guarantees stable identification of the system poles and participation factors. PolyMAX is applied to these three modes.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 11 of 23 force f at the lumped mass m5, and we define the step sine as the applied force. In addition, a total of twelve experiments were conducted in this simulation, and in each experiment, the magnitude of the applied force follows the step sine, which increases consistently. Detailed introductions and specifications Moreover, damping is assigned to increase the accuracy of the measurements made by the simulation. Proportional damping of the mass and stiffness matrices is assumed. As mentioned in the previous section, the stiffness matrix involved in this proportional damping is related to the stiffness KJ introduced in Equation (25): In Equation (40), a and b are constants, and in this study, values of 5 × e -1 and 7 × e -5 are assigned as these constants, respectively. In addition, the frequency resolution in the step sine experiment is 0.5 Hz and the time integral in the simulation is performed by the Newmark method [18].  Figure 4 shows the obtained FRFs from performed simulations. As shown in this figure, it can be confirmed that each mode has a different effect with respect to the magnitude of the external force. For instance, some modes (i.e., 2 nd , 4 th , 6 th and 8 th modes) show that their natural frequencies are shifted and that their peak values are decreased according to the increased magnitude of the force, whereas the other modes (i.e., 1 st , 3 rd , 5 th , 7 th and 9 th modes) hardly change their natural frequencies or their peak values. This is because the effect of the Iwan model is activated when relative displacement occurs in the region of interest. However, the modes that show no variations in the FRF results with respect to the changes of the external force all have in-phase motion in the interface between the masses of the Iwan model. The results shown in this figure are consistent to the   Figure 5b,c, there are no values defined in the range of y ( ) from a value of zero (i.e., ν − 1 = 0) to the other end of the function domain. This is because the variation of damping, which is obtained from the first and second experiments, is too small to sense its value. Therefore, the variation cannot be represented by a logarithmic scale. In addition, in Figure 5c, this causes the similarity of values between 1 and 2 in the function domain. This is because the vibration magnitude at the higher modes (i.e., 4 th , 6 th modes) decreased compared to the lower mode (i.e., 2 nd mode). We performed other experiments to investigate the results in Figure 5 in more detail. PolyMAX was used to estimate the modal parameters of the twelve performed experiments. We have harmonically excited the system to the estimated natural frequency in each experiment with the same amplitude of force that was used in the twelve performed experiments at the same positions. Relative displacements are confirmed with respect to these excitations at the nodes where the Iwan model is connected (i.e., the relative displacement between m 1 and m 6 , and the relative displacement between m 2 and m 7 ). These results are listed in Table 1.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 12 of 23 statements in Ref. [13]. In addition, this figure magnifies the results for the 2 nd , 4 th and 6 th modes where variations in the frequency results occur with respect to the magnitude of the external force. From these magnified plots, changes of damping and the natural frequency can be identified in each experiment. To identify the parameter χ with the proposed method in this study, the results of Equation (35) for these three modes have been identified. In this process, the modal parameters for these modes have to be estimated. Recently, polyreference least-squares complex frequency-domain method which is so-called PolyMAX [19] is introduced as a modal parameter estimation method conducted in the frequency domain. It is one of the experimental modal parameter estimation methods and it guarantees stable identification of the system poles and participation factors.
PolyMAX is applied to these three modes.     Table 1 lists half of the peak to peak values of the measurements, which are measured in the time domain one second after a long period of time (i.e., 9 seconds was used in this study). The results in the table can be associated with the parameter ϕ max of the Iwan model. ϕ max can be obtained by the following relation [12] with the defined parameters in Equation (39): The obtained ϕ max from Equation (41) for the defined Iwan model in Equation (39) is 0.0042. From the results in the table, it can be confirmed that for the second mode in the second Iwan model of the seventh experiment (i.e., the relative displacement between m 2 and m 7 ), a larger value is obtained compared to the value of ϕ max (0.0042) for the first experiment. This suggests that the macro-slip occurs in the seventh experiment for the first time. This also indicates that the system remained in the micro-slip region through the sixth experiment. The sixth experiment is associated with the 4 (ν-1) in the function domain in Figure 5. Therefore, in Figure 5a, the constant increase continues until function domain 4, but it shows abrupt changes at function domain 5. The same trend is confirmed in results of other modes. In the previous section, it was assumed that if one increases the amplitude of the external force as z times, then the relative displacement will linearly increase by z times. As confirmed in Table 1, however, it is observed that the results are less than z times (in this study, 2 times). The results of the estimated slope of the curve shown in Figure 5 show lower values than the true value. Eventually, this results in a larger amplitude (0.52~0.53) of the estimated χ than the true value. Next, R is estimated using the estimated χ and Equation (37), as shown in Figure 5. Similar to ϕ max , R can be obtained with the defined parameters in Equation (39) by the following relation which is introduced in Reference [12].
The obtained R by Equation (42) is 1.69 × e 5 . In this study, the constant B in Equation (37) was obtained in two different ways. First, relative displacements that are necessary to estimate this constant are obtained by the same way that is introduced in the previous section. This indicates that the system matrices M, C, and K J are obtained and estimate the relative displacements at nodes associated with the Iwan model by solving the inverse problem with these matrices. Then, R can be estimated by using the estimated B and Equation (37). The errors between the estimated R and the obtained true value of R by Equation (42) are listed in Table 2 in percent format (%). In Table 2, the inverse matrix shows the obtained results using the method proposed in this study with relative displacement, whereas the Newmark methods results are given in Table 1. If one examines the characteristics observed in this table, the same region corresponding to the linear region where the curve fitting process was performed in Figure 5 can confirm that the results in the table show more closed values than the true value of the defined R. In addition, the results, which were estimated using the results in Table 1 (Newmark method), are more accurate than the estimated results using the proposed method (Inverse matrix). Furthermore, as the mode number increases, the error of the estimated results of the proposed method also increases. This result is related to the mode shape because the estimated mode shape ϕ n from the lumped model is different to the true mode shape in dynamic behavior. Based on the results in the table, the difference between these two mode shapes (i.e., the lumped model vs. practical behavior) increases as the mode number increases.

Experiments
Experiments were performed on a partially connected double beam in order to verify the proposed method in Section 2. Figure 6 shows the double layer beam, which was composed of two identical aluminum cantilever beams (350 mm × 50 mm × 1 mm) with one end clamped and the other end fixed with glue. The glue was evenly distributed along the longitudinal direction of the beam for 20 mm. As shown in the figure, the exciter was located on the clamped side and the force transducer was located between the exciter and the double beam system. The vibration induced by the exciter was measured with an accelerometer, which was located in same position as the accelerator (Acc. 1).
In order to verify the proposed method, several experiments were conducted, as introduced in the theory section. Experiments were performed while increasing the voltage applied to the exciter twice. In these experiments, a linear relationship between the force and the applied voltage was assumed. The linear relationship was confirmed by attaching accelerometers to the area where the beams were glued. In the figure, the magnitude of the relative acceleration measured by Acc. 2 and Acc. 3 doubled in successive experiments. Thus, there was a linear relationship between the applied force and the voltage. The experiment was performed ten times from 5 mV to 2.56 V. The excitation frequency varied from 50 to 450 Hz based on the reliability of the exciter and the step sine used 2 Hz intervals. In order to verify the proposed method, several experiments were conducted, as introduced in the theory section. Experiments were performed while increasing the voltage applied to the exciter twice. In these experiments, a linear relationship between the force and the applied voltage was assumed. The linear relationship was confirmed by attaching accelerometers to the area where the beams were glued. In the figure, the magnitude of the relative acceleration measured by Acc. 2 and Acc. 3 doubled in successive experiments. Thus, there was a linear relationship between the applied force and the voltage. The experiment was performed ten times from 5 mV to 2.56 V. The excitation frequency varied from 50 to 450 Hz based on the reliability of the exciter and the step sine used 2 Hz intervals.   Figure 7 shows the FRFs measured in the 10 experiments performed. The second mode in the figure shows that the peak values at the resonance frequency decreased as the applied force increased. That is, the mode shape at the contact region was predominantly governed by the slip phenomenon by out-of-mode shape. Other modes were independent to the magnitude of the external force. Typically, this tendency is confirmed by the peak value in the fifth mode. In order to estimate the variation of natural frequencies and damping ratios in response to a change in external force, the two parameters of the Iwan model were estimated. The proposed method was applied to the measurements from the first seven experiments. The natural frequencies and damping ratios of the remaining three experiments were predicted using the estimated two parameters. The modal parameters, natural frequency, and damping ratio, predicted using the estimated two parameters from the remaining three experiments, were compared with the estimated modal parameters from the real measurements of the same experiments. First, Equation (35) was applied to the second mode in the figure.
The resultant y v values showed a large deviation at a small force level and a nearly linear trend where the force was increased. This was because slip rarely occurs when the external force is small, and in this case, the estimated results by the proposed method were more sensitive to the noise in the experiment. This tendency was also confirmed in the simulation results. The objective χ was estimated using the damping ratios calculated from the first seven experiments. Figure 8 represents the difference in energy dissipation in successive experiments; thus, χ was estimated by curve fitting the first six results by a 1 st order straight line (-). The estimated result was −0.37366. Next, the parameter R of the Iwan model was estimated using Equation (37), but first the mode shape of the partially connected double beam shown in Figure 6 had to be determined. The mode shape of a partially connected double beam was studied analytically and numerically in [20]. By applying Equation (37) to the mode shape vector, R was estimated as 3.5431 × e +5 . Utilizing the estimated χ and R, the natural frequencies and damping ratios were estimated using the measurements in the 8 th , 9 th , and 10 th experiments. The variation of natural frequency, ∆ω 2 n , when increasing the external force is represented in Equation (30). It can also be represented as follows, similar to Equation (33): . . .
where B = M m=1 P m (ϕ l )A χ+1 m 0 . Using the y-intercept in Figure 8 and Equation (35), |κ 1 |B was found to be 29.431. By substituting the estimated χ and R into Equation (A15), κ 1 was obtained. Utilizing the estimated value and |κ 1 |B, the following result was determined: B = 2.787 × e −4 . Additionally, κ 2 was found using Equations. (A12) and (A15). ∆ω 2 n , which represents the change in natural frequency, could then be obtained for the 8 th through 10 th experiments. Based on obtained ∆ω 2 n , the natural frequency ω 2 n was estimated. In order to verify the proposed method, several experiments were conducted, as introduced in the theory section. Experiments were performed while increasing the voltage applied to the exciter twice. In these experiments, a linear relationship between the force and the applied voltage was assumed. The linear relationship was confirmed by attaching accelerometers to the area where the beams were glued. In the figure, the magnitude of the relative acceleration measured by Acc. 2 and Acc. 3 doubled in successive experiments. Thus, there was a linear relationship between the applied force and the voltage. The experiment was performed ten times from 5 mV to 2.56 V. The excitation frequency varied from 50 to 450 Hz based on the reliability of the exciter and the step sine used 2 Hz intervals.  two parameters of the Iwan model were estimated. The proposed method was applied to the measurements from the first seven experiments. The natural frequencies and damping ratios of the remaining three experiments were predicted using the estimated two parameters. The modal parameters, natural frequency, and damping ratio, predicted using the estimated two parameters from the remaining three experiments, were compared with the estimated modal parameters from the real measurements of the same experiments. First, Equation (35) was applied to the second mode in the figure. The resultant yv values showed a large deviation at a small force level and a nearly linear trend where the force was increased. This was because slip rarely occurs when the external force is small, and in this case, the estimated results by the proposed method were more sensitive to the noise in the experiment. This tendency was also confirmed in the simulation results. The objective was estimated using the damping ratios calculated from the first seven experiments. Figure 8 represents the difference in energy dissipation in successive experiments; thus, was estimated by curve fitting the first six results by a 1 st order straight line (-). The estimated result was −0.37366. Next, the parameter R of the Iwan model was estimated using Equation (37), but first the mode shape of the partially connected double beam shown in Figure 6 had to be determined. The mode shape of a partially connected double beam was studied analytically and numerically in [20]. By applying Equation (37) to the mode shape vector, R was estimated as 3.5431 × e +5 . Utilizing the estimated and R, the natural frequencies and damping ratios were estimated using the measurements in the 8 th , 9 th , and 10 th experiments. The variation of natural frequency, , when increasing the external force is represented in Equation (30). It can also be represented as follows, similar to Equation (33): The second column of Table 3 is the variation in the natural frequencies ∆ω 2 n between successive experiments. The third column shows the variation of the natural frequencies estimated by applying the proposed method, Equation (43), to the experimental results from the first to the seventh trials. The fourth column shows the error between the experimental and estimated values. On average, the error was less than 10%. Table 4 shows the natural frequencies ω 2 n for the individual experiments, as obtained through Table 3. The second column is the measured ω 2 n in the experiment and the third column is the result estimated from Table 4. The fourth column shows the error between the experimental and the estimated values. Since the natural frequency is large, even if there was an error of 10% in the variation of the natural frequency listed in Table 3, the estimated natural frequency was very accurate at an error of less than 0.1%. ∆ω 2 n in the denominator of Equation (31) was directly related to the variation in natural frequencies measured experimentally. However, the relation between ∆c n and the damping ratio ζ J measured through the experiment should be confirmed. The denominator of the FRF for the individual modes in the modal domain, as derived from the equations of motion for the linear structure is: In Equation (44), ω J and ζ J denote the natural frequencies and damping ratios measured, respectively, by applying the experimental modal analysis method to the experiments. Through the imaginary part of Equations (31) and (44), ∆c n and ζ J have the following relation.
In Equation (45), ω n and ζ n are the natural frequency and damping ratio of the structure, respectively, when there is no Iwan model in the structure or when the structure vibrates with very small forces.
Assuming that the ratio of ω J and ω n is close to one, the damping ratio measured from the experiment is expressed as the right-hand side of Equation (45). From Equation (33), ∆c l0 = |κ 1 |B; thus, every ∆c lv in Equation (33) can be determined. ω J and ζ J were measured through experiments and ∆c n was predicted through the proposed method.
Equation (46) was derived from Equation (45) and shows the damping ratio change when the external force increases. Subscript J is the damping ratio measured in the experiment, and subscript v is the index of the experiment performed with increasing force level. As a result, the damping ratio measured from the 8 th through 10 th experiments could be estimated from Equation (46). The second column of Table 5 is the variation of damping ratio ∆ζ Jv between consecutive experiments and the third column is the variation of the estimated damping ratio when applying the proposed method to the experimental results of the first seven experiments. The fourth column shows the error between the experimental and estimated values. The overall average error was less than 7%. Table 6 shows the damping ratios for individual experiments estimated based on Table 5. The second column was measured in the experiment and the third column is from Table 5. The fourth column shows the error between the experimental value and the estimated value. Since the damping ratio was larger than the variation in the damping ratio, the estimated damping ratio guaranteed a very accurate estimate with an error of less than 1%. This held even if the variation of the damping ratio estimated in Table 5 was 7%. Through the experiments performed, the estimation of variations in the natural frequency and the damping ratio had an error of 10% compared with the actual measurements. The natural frequencies and damping ratios; however, estimated based on the estimated variations were confirmed to be very accurate with less than 1% error. The reliability and effectiveness of the proposed method were, thus, verified.  Table 6. List of ζ Jv .

Order of Expt. Measured (%) Estimated (%) Error (%)
Funding: This work was supported by the Dong-A University research fund.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The coefficient a 0 for the dc term in Equations (12) and (13) will be investigated. Applying Equation (11) to Equation (13) The constant γ 2 in Equation (11) is located within the integral region [−T/2, T/2], and therefore, the terms cancel each other out. By accounting for the configuration of the sine function in the first and the third term on the right hand side, the following result can be obtained: From Equation (A2), it can be confirmed that there is no dc term in the variation force that occurs during harmonic behavior. Next, the coefficient a 1 for the first harmonic component will be investigated. Similar to the case of a 0 , applying Equation (11) to Equation (13) produces the following: Similar to the case of Equation (A3), if taking the sine and cosine function into account, the terms in the first and the third from the top of the right hand side of Equation (A3) can be transformed: The first term on the right hand side of Equation (A5) can be obtained by a change of variable. If transforming sin ωt to the variable z, then the integral region will be transformed to [−1, 1]. Eventually, the result of the first term will be as follows: The second term in Equation (A5) can be obtained without difficulty. Therefore, the coefficient a 1 is ∴ a 1 = 4 π γ 1 χ + 3 2 χ+2 − γ 2 (A7) The coefficient b 1 of the second harmonic component of the base component will be investigated. Similar to the previous case of applying Equation (11) to Equation (13), by taking the configurations of