A novel modal superposition method with response dependent nonlinear modes for periodic vibration analysis of large MDOF nonlinear systems

Mechanical Systems and Signal Processing 135 (2020) 106388

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

A novel modal superposition method with response dependent nonlinear modes for periodic vibration analysis of large MDOF nonlinear systems Erhan Ferhatoglu 1, Ender Cigeroglu ⇑, H. Nevzat Özgüven Middle East Technical University, Department of Mechanical Engineering, 06800 Ankara, Turkey

a r t i c l e

i n f o

Article history: Received 28 February 2019 Received in revised form 15 September 2019 Accepted 21 September 2019

Keywords: Modal superposition of nonlinear systems Response dependent nonlinear modes Nonlinear vibrations Describing function method Nonlinear normal modes

a b s t r a c t Design of complex mechanical structures requires to predict nonlinearities that affect the dynamic behavior considerably. However, finding the forced response of nonlinear structures is computationally expensive, especially for large ordered realistic finite element models. In this paper, a novel approach is proposed to reduce computational time significantly utilizing Response Dependent Nonlinear Mode (RDNM) concept in determining the steady state periodic response of nonlinear structures. The method is applicable to all type of nonlinearities. It is based on the use of RDNM which is defined as a varying modal vector with changing vibration amplitude. At steady-state, due to periodic motion, it is possible to define equivalent stiffness due to nonlinear elements as a function of response level which enables one to create new linear systems at each response level by modifying original stiffness matrix of the underlying linear system. In this method, a new linear system is defined at each response level corresponding to each excitation frequency step, and modal information of these equivalent linear systems is used to construct RDNMs which forms a very efficient basis for the nonlinear response space. The response of the nonlinear system is then written in terms of these RDNMs instead of the modes of the underlying linear system. This reduces the number of modes that should be retained in modal superposition method for accurate representation of solution of the nonlinear system, which decreases the number of nonlinear equations, hence the computational effort, significantly. Dual Modal Space method is employed to decrease the computational effort in the calculation of RDNMs for realistic finite element models, i.e. for large MDOF systems. In the solution, nonlinear differential equations of motion are converted into a set of nonlinear algebraic equations by using Describing Function Method, and the numerical solution is obtained by using Newton’s method with arc-length continuation. The method is demonstrated on two different systems. Accuracy and computational time comparisons are performed by applying different case studies which include several different nonlinear elements such as gap, cubic spring and dry friction. Results show that the proposed method is very effective in determining periodic response of nonlinear structures accurately reducing the computational time considerably compared to classical modal superposition method that uses the modes of the underlying linear system. It is also observed that the variation of natural frequency with energy level in a nonlinear system can be approximately obtained by using RDNM concept. Ó 2019 Elsevier Ltd. All rights reserved.

⇑ Corresponding author. 1

E-mail address: [email protected] (E. Cigeroglu). Presently at the Department of Mechanical and Aerospace Engineering of Politecnico di Torino, Italy.

https://doi.org/10.1016/j.ymssp.2019.106388 0888-3270/Ó 2019 Elsevier Ltd. All rights reserved.

2

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

1. Introduction Dealing with nonlinearity in structures has always been a challenge in the industrial applications. Although the performance of computers has considerably increased with developing technology, models are also becoming complex and detailed, yielding large number of degrees of freedom (DOFs). Although, in the past three decades, researchers have developed methods to predict the dynamic behavior of nonlinear structures pretty accurately within reasonable computational time, solution techniques of nonlinear dynamic equations still need to be improved since the demand is also increasing. Dynamic response analysis of nonlinear multi degree of freedom (MDOF) systems started in the late 1980s by using step by step numerical integration techniques to solve second order differential equation of motion in time domain [1,2]. Although these methods are quite successful to determine vibration response of structures, they are not suitable for very large MDOF systems due to very high computational time required. Hence, in most of the studies frequency domain solution techniques were proposed. The earliest studies [3–10] focused on the steady state harmonic vibration analysis of MDOF nonlinear structures. Reduced order models were obtained either by partitioning the receptance matrices and making iteration only for nonlinear DOFs [7,11,12], or by using modal superposition approach with linear normal modes (MSM-LM) [8]. In literature, various solution approaches were proposed for periodically excited nonlinear systems possessing different types of nonlinearities. In Describing Function Method (DFM) [7], internal nonlinear forces are expressed as a multiplication of the so-called nonlinearity matrix (which is a function of unknown response amplitude vector) with displacement vector. In this way, nonlinear internal forces are included into the analysis in terms of an equivalent stiffness matrix which made it possible to extend several methods, such as structural modification and coupling methods, developed for linear systems to nonlinear systems. Later, the very-well known Harmonic Balance Method (HBM) was applied to MDOF systems by representing both the periodic solution and the internal nonlinear forces using Fourier series, and balancing similar terms. DFM and HBM are mathematically equivalent to each other. Internal nonlinear forcing vector is written as a vector in HBM; whereas, it is written as a matrix–vector multiplication, which again results in a vector, in DFM. There exist many studies which focus on to increase the usability and robustness of HBM by calculating internal nonlinear forces in more efficient ways (e.g., see [13–15]). Further different concepts and strategies were also developed to improve the solution methods used in the estimation of nonlinear vibration response of structures. A considerable effort has been spent on the construction of reduced order models (ROMs) [16–19]. ROMs provide basis functions that transform the equations of motion into a new domain in which the size of the nonlinear system to be solved is much smaller. Thus, the computational time required for the solution decreases directly [20,21]. Furthermore, when the number of nonlinear algebraic equations is reduced, stability of the nonlinear solver also improves [22]. Nonlinear Normal Mode (NNM) concept was first introduced by Rosenberg [23]. NNM concept has received more attention in recent years, and in various studies on dynamic analysis of MDOF systems this concept has been used [24–30]. NNM directly depends on the total energy in the system; therefore, mode shapes and resonance frequencies vary as the response level changes. In [28], it is shown that the use of just one NNM in the neighborhood of the resonance region captures the dynamic behavior accurately. However, calculation of NNMs for large MDOF systems is not practical and it is time consuming [20], particularly when the nonlinearities are distributed across the degrees of freedom. For further studies and more comprehensive reviews on NNMs, the reader can refer to [21,31]. In this paper, a new concept called Response Dependent Nonlinear Mode (RDNM) is introduced, and a novel method is proposed to obtain RDNMs of large ordered nonlinear structures by making use of stiffness variation of nonlinear elements as the response level and pattern changes, and then to use them for the computation of steady state response of the nonlinear system to periodic forcing. Equivalent stiffness matrix representing nonlinear elements depends on the vibration amplitudes; hence, it changes with varying displacement amplitude levels and patterns. In a recent study of authors [21], hybrid mode shape concept is developed and it is shown that using less number of hybrid mode shapes is sufficient to capture the response of the nonlinear system, which decreases the computational time, significantly. In [21], different linear systems are defined by using equivalent stiffness matrices corresponding to specific vibration amplitudes where the nonlinear force saturates. Then, hybrid mode shapes (HM) are constructed by using a linear combination of the mode shapes of these linear systems. In this study, this approach is extended further by introducing RDNMs which results in a much better basis for the nonlinear response space than the previously introduced hybrid mode shapes. It should be noted that since the equivalent stiffness matrix of the nonlinear system changes throughout the frequency range considered for a periodic excitation, so does the RDNM. In this study, RDNMs are used in nonlinear response computations, which reduces the number of nonlinear equations to be solved, and hence, the computational effort. In order to compute RDNMs for large finite element models (FEMs), the change in the stiffness matrix due to the variation in vibration amplitude is considered as a structural modification, and then by using Dual Modal Space approach, the computational effort is significantly reduced. The application and the performance of the proposed method are demonstrated firstly on a lumped parameter model and on a finite element model of a structure by performing several case studies. It is observed from these case studies that RDNMs obtained using modal information of successive equivalent linear systems provide excellent results by keeping the number of modes used in the solution at a minimum.

3

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

2. Theory Consider the general expression for the differential equation of motion of a nonlinear structure under periodic excitation modelled as an n DOF discrete system

€ðt Þ þ C
x_ ðt Þ þ iH
xðt Þ þ K
xðtÞ þ f N ðxðt Þ; x_ ðt Þ; :::Þ ¼ fðtÞ; M
x

ð1Þ

where M, C, H and K are mass, viscous damping, structural damping, and stiffness matrices representing the linear part of the system, respectively. f N ðxðt Þ; x_ ðt Þ; :::Þ and fðtÞ represent internal nonlinear forcing and external periodic excitation vectors, respectively. xðt Þ is the vector of the displacement response of the system, dot denotes differentiation with respect to time and i is the unit imaginary number. It should be noted that structural damping matrix, H, is only applicable to harmonic motion. Since periodic vibrations are of interest, excitation force and response of the nonlinear system can be written in complex form as follows

!

nh X

f ðtÞ ¼ f 0 þ Im

f h eihxt ;

ð2Þ

h¼1

xðt Þ ¼ x0 þ Im

nh X

! xh e

ihxt

;

ð3Þ

h¼1 th

where f h and xh are complex vectors of external forcing and response corresponding to the h harmonic, respectively. f 0 and x0 are the bias components of the external forcing and nonlinear modal response vectors, which are real. nh represents the total number of harmonics in the above expressions which is considered to be the same for both external forcing and response. Similarly, for periodic motion, internal nonlinear forcing vector can be written as

f N ðxðt Þ; xÞ ¼ f N;0 þ Im

nh X

!

f N;h e

ihxt

;

ð4Þ

h¼1

where f N;0 is the real bias amplitude vector and f N;h is the complex vector of the internal nonlinear forcing. It should be noted that f N;0 and f N;h are functions of xh ðh ¼ 0; 1;


; nh Þ and x. Substituting Eqs. (2), (3) and (4) into Eq. (1), the following set of nonlinear algebraic equations in frequency domain are obtained

K
x0 þ f N;0 ¼ f 0 ;

ð5Þ

2

½K  ðhxÞ M þ ihxC þ iH
xh þ f N;h ¼ f h for h ¼ 1; 2;


; nh :

ð6Þ

Nonlinear internal forcing vector f N;h ðh ¼ 0; 1;


; nh Þ can be calculated by using Describing Function Method (DFM) developed by Tanrıkulu et al. [7]. As the theory of the DFM used in this study is given in detail in [7], only the main equations are summarized here for the sake of completeness. DFM enables us to write internal nonlinear forcing vector as a multiplication of a response dependent complex nonlinearity matrix and complex displacement vector as follows

h i im ^Þ ¼ Dh ðxÞ
xh ¼ Dre f N;h ðx h ðxÞ þ iDh ðxÞ
xh ;

ð7Þ

h iT th ^ ¼ xT0 xT1


xTn where x and xh is the complex displacement vector corresponding to the h harmonic which is real for h th

h ¼ 0 and complex for h  1. Similarly, Dh ðxÞ is the displacement dependent nonlinearity matrix of the h harmonic which is real for h ¼ 0 and complex for h  1. Superscripts re and im indicate the real and imaginary parts, respectively. For single harmonic motion, elements of the complex nonlinearity matrix can be obtained as [7]

Dk;k ¼ mkk þ

Ne X

mkj and Dk;j ¼ mkj ;

ð8Þ

j¼1 j–k

where mkj is the harmonic input describing function of a nonlinear element in the system and can be described as equivalent kj

th

th

th

complex stiffness of the internal nonlinear force, f N , acting between the k and the j DOFs. The k component of f N can be kj fN

th

obtained by summing all for j ¼ 1 . . . N e . N e represents the number of nonlinear elements attached to the k DOF. Describing function value, mkj , for single harmonic motion can be calculated as follows

mkj ¼

Z 2p

i

pY

kj

0

 kj  f N ykj ; h eih dh ;

ð9Þ

4

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

where h ¼ xt and ykj ¼ Asinðh þ wÞ. Here A and w represent the sine coefficient and phase of the relative motion between the th

k

th

and the j

th

DOFs, respectively. Y kj is the complex amplitude of the relative motion between the k

th

and the j

th

DOFs. th

Similarly, for multi harmonic motion, the h harmonic describing function of the nonlinear element between the k and th

the j

DOFs can be approximated as follows

mkjh ¼

  kj kj kj F kj N;h X 0 ; X 1 ; . . . ; X nh Y kj h

¼

i

Z 2p

Y kj h

p

0

 kj  f N ykj ; h eihh dh ;

ð10Þ

where

ykj ¼ xk  xj ¼

nh   X kj ykj s;h sinðhhÞ þ yc;h cosðhhÞ ;

ð11Þ

h¼1

( Y kj h

¼

X kh  X jh

for k–j

X kh

for k ¼ j

:

ð12Þ th

th

kj j k Here, ykj s;h and yc;h are the sine and cosine coefficients of the relative motion between the k and the j DOFs. X h and X h are th

th

th

harmonic, respectively. Y kj is the  h th th th kj kj kj kj complex amplitude of the h harmonic of the relative motion between the k and the j DOFs. F N;h X 0 ; X 1 ; . . . ; X nh repre-

the complex displacement amplitude of the k

and the j

DOFs corresponding to the h

th

th

sents the complex amplitude of the internal nonlinear forcing acting between the k and the j DOFs corresponding to the th

h harmonic. It should be noted that Eq. (10) enables DFM to become completely correlated with HBM, in which Fourier   kj kj kj coefficient of the internal nonlinear force, F kj N;h X 0 ; X 1 ; . . . ; X nh , can be computed in advance by utilizing any methodology such as analytical integration, numerical integration or Alternating Frequency/Time approach as in HBM, from which the th

describing function values can then be obtained by employing Eq. (10). Therefore, nonlinearity matrix for the h harmonic can be obtained as follows

Dhk;k ¼ mkk h þ

Ne X

mkjh and Dhk;j ¼ mkjh ;

j¼1 j–k

ð13Þ

Substituting Eq. (7) into Eqs. (5) and (6), the following set of nonlinear algebraic equations are obtained

  K þ Dre 0
x0 ¼ f 0 ;

ð14Þ

h i 2 im K  ðhxÞ M þ ihxC þ iH þ Dre
xh ¼ f h ðh ¼ 1; 2;


; nh Þ: h þ iDh

ð15Þ

The number of nonlinear algebraic equations defined by Eqs. (14) and (15) is ð2nh þ 1Þn, where n is the number of total DOFs in the system. It should be noted that Eqs. (14) and (15) are coupled through complex nonlinearity matrices. The details of the DFM, given for fundamental harmonic only, can be found in [7]. In classical modal superposition method, utilizing expansion theorem, response of the system is expressed as follows

xh ¼ U
gh ðh ¼ 0; 1;


; nh Þ;

ð16Þ

where U is the ðn  mÞ mass normalized modal matrix of the underlying linear system, gh is the ðm  1Þ vector of modal th

coordinates of the h harmonic which is real for h ¼ 0 and complex for h  1. m is the number of modes retained in the modal superposition method and may take different values for different modal transformations used in this study. It should be noted that the total DOFs, n, is quite large for most engineering applications, and m is chosen as small as possible while preventing any loss of accuracy in the solution. Substituting Eq. (16) into Eqs. (14) and (15), and pre-multiplying both sides by UT , following equations are obtained T ½X þ UT Dre 0 U
g0 ¼ U f 0 ;

h

  i 2 T im X þ UT Dre
gh ¼ UT f h ðh ¼ 1; 2;


; nh Þ: h U  ðhxÞ I þ i Hd þ hxCd þ U Dh U

ð17Þ ð18Þ

Here I is identity matrix, Cd and Hd are modal viscous and structural damping matrices, respectively, which are diagonal if classical damping is assumed, and X is a diagonal matrix composed of squares of natural frequencies. The number of nonlinear algebraic equations defined by Eqs. (17) and (18) is ð2nh þ 1Þm. Since the number of modes used in modal superposition method is significantly less than the number of total DOFs (m  n), the number of nonlinear algebraic equations is

5

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

decreased drastically. The number of linear system modes, m, that should be retained in modal expansion depends on the excitation forcing vector, f and the nonlinear forcing vector, f N . It should be noted that f N changes with the excitation frequency. Due to the nonlinear internal forcing in the system, mode shapes of the linear system do not form a good basis to express the nonlinear solution, and therefore the number of linear modal vectors that should be used in the modal superposition method increases, especially, if the effect of the nonlinearity increases. Choosing a better set of basis vectors, the number of nonlinear equations can be kept at a minimum which decreases the computational effort drastically. Using DFM, the nonlinear forcing vector is written as a displacement dependent complex matrix multiplied by complex displacement vector results in the set of nonlinear equations defined in Eqs. (14) and (15). For the first harmonic, Eq. (15) can be written as follows

h i  im 2 K þ Dre
x1 ¼ f 1 : 1  x M þ ixC þ iH þ iD1

ð19Þ

It can be seen from Eq. (19) that the stiffness matrix is modified by the real part of the nonlinearity matrix. Considering this modification in the stiffness matrix, a new eigenvalue problem (EVP) can be formed at each response level as follows



 2 K þ Dre 1 u ¼ xN M
u:

ð20Þ

It should be noted that Dre 1 is a function of the displacement amplitude vector and if it is known, standard eigenvalue solvers can be used to obtain the natural frequencies, xN , and modal vectors, u, corresponding to the displacement amplitude vector considered. Then the response of the system can be written as

xh ¼ UN
qh ðh ¼ 0; 1;


; nh Þ;

ð21Þ th

where qh is the vector of the new modal coordinates for the h harmonic. UN ¼ ½ u1 u2


um  is the mass normalized nonlinear modal matrix, where ui are the RDNMs calculated by solving the EVP given in Eq. (20). Moreover each RDNM is defined for the response level considered in the determination of Eq. (20). They are recalculated at each frequency step (solution point) once and kept unvaried during the iterations in the solution of nonlinear algebraic equations. Therefore, they are a function of displacement vector, and change with the total energy level of the system similar to NNMs. Hence, the new nonlinear modal domain defined by the transformation given by Eq. (21) provides a better basis for the nonlinear response space which decreases the number of modes used in the modal superposition method significantly. It is important to note that RDNMs or UN obtained from Eq. (20) fully uncouples the undamped nonlinear system defined by Eq. (19). Therefore, it can be hereby concluded that the number of modes used in the MSM can be evidently minimized similar to linear systems. Moreover, if the nonlinear modes of the system are well separated, a single RDNM would be sufficient to capture the nonlinear response around the corresponding resonance frequency as in the case of linear systems. It should also be noted that since the focus of this study is on the determination of steady-state periodic vibration response of nonlinear structures, no emphasis is given to identify the stability of the solutions. One may refer to [22,33] for comprehensive stability analysis of periodic solutions of nonlinear structures. Utilizing Eq. (21), in Eqs. (5) and (6), and pre-multiplying both sides by UTN yields the following result

UTN KUN
q0 þ UTN f N;0 ¼ UTN f 0 ;

ð22Þ

h i 2 UTN K  ðhxÞ M þ ihxC þ iH UTN
qh þ UTN f N;h ¼ UTN f h ðh ¼ 1; 2;


; nh Þ:

ð23Þ

Eqs. (22) and (23) define a system of nonlinear algebraic equations whose size is proportional to the number of modes used, m, which is significantly smaller than the number of total DOFs. Although this approach decreases the computational effort drastically, the set of nonlinear algebraic equations defined by Eqs. (22) and (23) requires an EVP of size n  n to be solved at every solution point which normally increases the computational cost. In order to overcome this problem, Dual Modal Space method [32] is used. Let us consider Eqs. (17) and (18) which are in modal domain. Now, a second modal transformation is defined as follows 

gh ¼ U
qh ðh ¼ 0; 1;


; nh Þ;

ð24Þ T

and it is substituted into Eqs. (17) and (18). Note that both sides of these equations are multiplied by U and since f N;h ¼ Dh
xh as shown in Eq. (7), the nonlinear forcing vector defined by Dh
xh utilizing DFM can be replaced with f N;h as done in Eqs. (22) and (23) which makes equations easier to follow. It should be noted in this form HBM can be easily utilized in the solution of the nonlinear algebraic equation set, as well. Then, the following equations are obtained 



UT X U
q0 þ UTN f N;0 ¼ UTN f 0 ;

ð25Þ

h   i 2 UT X  ðhxÞ I þ iHd þ ihxCd U
qh þ UTN f N;h ¼ UTN f h ðh ¼ 1; 2;


; nh Þ:

ð26Þ

6

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388 



Here, UN ¼ U U. It is the mass normalized nonlinear modal matrix or RDNM matrix. U is the modal matrix obtained from the solution of the following eigenvalue problem 





2 ½X þ UT Dre 1 U
u ¼ x I
u :

ð27Þ

Eq. (27) can be obtained from Eq. (20) by using the Dual Modal Space method defined by transformations given by Eqs. h i      (16) and (24) where U ¼ u1 u2


um and xN is replaced with x. It should be noted that Eqs. (25) and (26) define the final nonlinear algebraic equation set to be solved for the vibration response of the structure and it does not contain any system matrices due to the utilization of Dual Modal Space method. th

Here, qh is the vector of unknown modal coefficients defined by the transformation given in Eq. (24) representing the h harmonic. The total number of real nonlinear algebraic equations becomes ð2nh þ 1Þm. It should be noted that, different from the frequency response analysis of linear systems, here the solution requires iterations since the nonlinear forcing vector (or the nonlinearity matrix) is response level dependent. In the approach proposed in this study, the modal matrix of the linear system, U, is obtained once by performing modal analysis to the underlying linear structure before starting nonlinear frequency response analysis, for which a commercial finite element software can be used when realistic finite element models are considered. Utilizing Eq. (27), the EVP that should be solved in following steps will have drastically reduced order. If all the modes of the linear system were used in Eq. (27), the computational effort would become identical to the EVP defined by Eq. (20). However, in practical applications, only p number of modes of the linear system will be used so that the size of the EVP defined in Eq. (27) will be p  p, where p  n. Therefore, the computational expense required for the solution of the EVP in following steps will be decreased significantly. For the solution of the resulting nonlinear algebraic equations given by Eqs. (25) and (26), Newton’s Method with Arclength Continuation is used in this study. This set of nonlinear algebraic equations can be written as a real valued residual vector function in the following form

b ¼ ½qre qim T :

Rðb; xÞ ¼ 0;

ð28Þ

Here qre and qim represent the real and imaginary parts of the unknown modal coordinate amplitudes corresponding to all harmonics, respectively. Iterative formula for the Newton’s method with arc-length continuation can be written as follows

" @Rðb;xÞ ajþ1 k

¼

ajk

@b @hðb;xÞ @b



@Rðb;xÞ @x @hðb;xÞ @x

#1     

(  j

j

bk ;xk

j

Rðbk ; xjk Þ j

hðbk ; xjk Þ

) ;

ð29Þ

where

( ajk

¼

j

bk

xjk

) ;

j

ð30Þ T

h ðbk ; xjk Þ ¼ ðajk  ak1 Þ
ðajk  ak1 Þ  s2 ¼ 0 :

ð31Þ

Here, s is the arc length parameter, which defines the radius of the n-dimensional hypothetical sphere in which the next solution is sought, k is the current solution point and j is the iteration number. More information on Newton’s Method with Arc-length Continuation can be found in [33–35]. 3. Case Studies In this section, the application of the method developed is demonstrated on two different systems. The first one is a 20DOF lumped parameter model on which the validation of the proposed approach is presented. The second system is a finite element model of a realistic structure with 5400 DOF, where the performance of the method can be highlighted clearly. Since the method is applicable to systems with all types of nonlinearities, different nonlinear elements are utilized using different case studies. 3.1. Lumped Parameter Model The 20-DOF underlying linear lumped parameter system used in the following case studies is shown in Fig. 1. Parameters of the model are k ¼ 30000N=m and m ¼ 0:1kg. Proportional viscous damping is utilized in the system. Viscous damping values, forcing amplitudes and locations of the external forcing are varied for each case study. Dynamic behavior of the system is strongly affected by the force location; hence, it is used to change the participation of the modes to the total solution. The purpose is to investigate the performance of the proposed method under different conditions. Two case studies are presented by using this lumped parameter system.

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

7

Fig. 1. 20-DOF Underlying Linear System.

3.1.1. Case Study 1: Gap Nonlinearity In the first case study, a unilateral elastic stop element with stiffness kg is considered. Internal nonlinear force of the nonlinear element can be expressed as follows

_ Þ¼ f N ðxðtÞ; xðtÞ



0

xðtÞ < d

kg ðxðtÞ  dÞ xðtÞ  d

;

ð32Þ

where d is breakaway amplitude. The main aim of this case study is to validate the proposed method on an example where the nonlinear effects are clearly seen. For this purpose, a moderately nonlinear element is used, whose properties are specified as kg ¼ 1500N=m; d ¼ 0:02m. The nonlinear element is attached between the 10th DOF and the ground in order to increase the effect of nonlinearity. It should be noted that due to symmetry of the underlying linear system, maximum deflection at the first vibration mode occurs at the 10th and the 11th DOFs. The nonlinear system is investigated around the first natural frequency of the linear system, i.e. around the first resonance frequency of the nonlinear system, where the effect of the gap nonlinearity is the strongest. Three different external periodic forces are applied to the 10th DOF with 5 N, 10 N and 20 N amplitudes. A relatively small proportional viscous damping value, i.e. 1%, is applied to the system and 7 harmonics are utilized in the response calculations. Fig. 2 shows the displacement amplitude of the 10th DOF, which is defined as half of the peak-to-peak response in one period. The nonlinear response given in the figure is obtained by employing MSM-LM and taking all modes of the system. It is used as the basis of comparison for the results obtained by the proposed method, i.e. modal superposition method with RDNM. It is clearly seen that nonlinearity strongly affects the response at the resonance region where hardening effects and jump phenomena are observed. The first resonance frequency shifts approximately 6% for 20 N excitation level. Participation of the zeroth and higher harmonics has also important effects on the response (which are not given here for brevity), although single resonance peak is observed. It is seen from Fig. 2 that using even 1 RDNM calculated by the proposed method captures the nonlinear response of the system very accurately for all excitation levels. This demonstrates that RDNMs form a very strong basis for the nonlinear response space since it evolves as the energy level of the system changes throughout the frequency range. One can also increase the number of RDNMs used in order to increase the response accuracy; however, it is not needed for this particular case since using just one RDNM is sufficient to determine the nonlinear response of the system accurately. Results obtained are also compared with the values calculated by time integration which validate the method developed. It should be noted that, since similar results are obtained for the other DOFs, only the response of the 10th DOF is presented here. 3.1.2. Case Study 2: Dry Friction and Cubic Stiffness Together In this case study, dry friction and cubic stiffness nonlinear elements are used simultaneously. Here, it is specifically intended to increase the nonlinear effect in the system by coupling different types of nonlinearities. Moreover, the model is excited by two different external forces in order to increase the contributions of higher modes. External forces are applied to the 6th DOF and the 12th DOF which are defined as 15sinðxt Þ and 5sinðxt þ p=2Þ N, respectively. 0.5% of proportional viscous damping is used in the system. In this case study, only the fundamental harmonic is used in DFM for simplicity, since for dry friction and cubic stiffness nonlinearities majority of the contribution to the response comes from the fundamental harmonic. It should also be noted that the effect of higher harmonics on system response is out of the scope for this study. Readers may refer to [11,15] to investigate the effects of higher harmonics on system response. For dry friction, 1D model with constant normal force is used. Expression of the internal nonlinear force is as follows

_ Þ¼ f N ðxðtÞ; xðtÞ



_ lN sign ðwðtÞÞ

slip

kd ðxðtÞ  wðtÞÞ

stick

;

ð33Þ

where xðtÞ, wðtÞ, lN and kd represent the relative motion between the two ends of the element, slip motion, slip force and tangential contact stiffness, respectively. Internal nonlinear forcing for cubic stiffness nonlinearity having stiffness of kc can be given as

_ Þ ¼ kc xðtÞ3 : f N ðxðtÞ; xðtÞ

ð34Þ

8

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

Fig. 2. Normalized displacement amplitudes of the 10th DOF vs. frequency.

Dry friction element is attached between the 6th DOF and the 10th DOF; whereas, cubic stiffness element is inserted between the 13th DOF and the ground. Parameters of the nonlinear elements are kd ¼ 1000N=m, lN ¼ 5N and kc ¼ 5  107 N=m3 for cubic spring element. Fig. 3 shows the displacement amplitude of the 5th DOF, at which the effect of nonlinearity also around the second resonance region, in addition to the first one, can be observed clearly. It can be seen from Fig. 3 that resonance frequency shifts more than 50% around the first resonance region as a result of strong nonlinear effects. In Fig. 4, responses of the 5th DOF around the first resonance region are given by using all modes of the linear system (exact nonlinear response), using 1 RDNM and 10 LMs. It is observed from the results obtained that using just 1 RDNM gives very accurate results throughout the frequency range considered. It also gives better results than using 10 LMs. It is expected that LMs do not form a strong basis for this specific case, since the nonlinearity is quite strong. However, as observed from the results obtained, even a single RDNM captures almost the exact dynamics. Similar results are also obtained at the second resonance region as can be seen from Fig. 5. Similar to the first resonance region, using just 1 RDNM gives quite satisfactory results. Although the vibration amplitude is relatively low compared to the first resonance and a very small resonance frequency shift occurs, a similar accuracy in the response is hardly obtained by using even 10 LMs. It should be noted that in the response computations, since Figs. 4 and 5 depict the response of the system around the first and the second resonance regions, only the first RDNM is utilized in the former case and only the second RDNM is used in the latter case. Moreover, the number of RDNMs used may be increased in order to enhance the response accuracy if required. Since the computational time spent in modal superposition method is directly a function of the number of modes used in the response calculations, using just 1 RDNM has drastic reduction in the computational time spent compared to the use of 10 LMs or all linear modes. In the analyses, it is observed that using 1 RDNM reduces the computational time 73% and 87% compared to using 10 LMs and all linear modes, respectively. Reduction in computational effort becomes very important especially for very large systems which are considered in the third case study. 

It should be noted that in this study, the EVP given by Eq. (27) is solved once, i.e. U is computed once, at the beginning of every frequency step (solution point) of the nonlinear solution process by using the displacement pattern obtained in the previous solution point. However, the computational expense of the EVP solution can be reduced if it is not calculated at every solution point. Instead, the computational cost can be decreased further by solving the EVP at every r number of solution points, i.e. by freezing the RDNMs used in the solution process for r solution points. Selection of r affects the accuracy of

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

Fig. 3. Displacement amplitudes of the 5th DOF with respect to frequency.

Fig. 4. Displacement amplitudes of the 5th DOF around the first resonance.

9

10

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

Fig. 5. Displacement amplitudes of the 5th DOF around the second resonance.

the solution as well. Hence, an optimum r value for each problem at hand can be found by obtaining a balance between the computational cost and the accuracy obtained. The effect of r (number of solution points RDNMs are frozen) is studied and the results obtained are given in Table 1 for the response of the 5th DOF around the first resonance region. In the analyses, arc-length parameter is selected so small that there is no precision loss in response calculations. As a result, the response figures that are obtained for each analysis exactly overlaps with each other in the frequency range considered. They are not presented here as a separate figure for brevity. As shown in Table 1, computational time is reduced from 100% to almost 60% by freezing the RDNM used as observed from the second column of Table 1. Moreover, convergence problems arise if the RDNM used is frozen for 10 solution points in the resonance region in this particular case study. This problem can be overcome by implementing an adaptive arc-length parameter; however, for the purpose of time comparisons, constant arc-length parameter is used in this case study. In the above computations, in each step, even though the number of RDNMs (UN ) used in the response computation is only 1, the EVP defined in Eq. (27) has the same size with the original system order. Hence, the number of modes, m, utilized 

in U and U are both 20. However, the computational cost of the proposed method obviously depends also on the size of the EVP solved in each step. Therefore, next, the effect of the size of the reduced EVP on the computational effort by decreasing 

the number of modes, m, in U and U is studied. While reduction in the size of the EVP reduces the computational time, it may also reduce the accuracy of the results. Therefore, before investigating the change in computational effort, firstly, the effect of the size of the reduced EVP on the accuracy of response is studied. The results are shown in Fig. 6. It is observed that using 1 RDNM obtained by solving a reduced size EVP decreases the solution accuracy, not only due to using truncation in modal summation, but also due to omitting some terms in the nonlinear algebraic equations given by Eqs. (25) and (26). The decrease in computational time by reducing the sizes of the EVP is given Table 2. It is observed that the computational time in this specific example is reduced, approximately 8%, 12.5% and 16% compared to using full size EVP, by reducing the size of the EVP by 25%, 50% and 75%, respectively. The trade-off here is again between decreased response accuracy and enhanced computational time. As shown in Case Study 1 and 2, RDNM calculated with the proposed method forms a very strong basis for the periodic response analysis of nonlinear systems. Consequently, by using just 1 RDNM, nonlinear response of the system can be captured very accurately in systems with well-separated modes. However, when there is modal coupling, the method proposed can still be employed by using more than one RDNM. In order to demonstrate the use of the method proposed in such a case,

11

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388 Table 1 Computational Time Comparison for Different Number of EVP Solution. EVP Points

Total Number of Solution Points that EVP is solved

Comp. Time with respect to Reference Solution (%)

Every Solution Point (r ¼ 0) Once in every 2 Solution Points (r ¼ 1) Once in every 3 Solution Points (r ¼ 2) Once in every 5 Solution Points (r ¼ 4) Once in every 10 Solution Points (r ¼ 9)

2233 1114 742 445 222

100 (Ref.) 79.7 70.3 64.0 60.9

Fig. 6. Displacement amplitudes of the 5th DOF around the first resonance for different sizes of the EVP solved.

Table 2 Comparison of Computational Time Spent for Different Sizes of the EVP. EVP Size

Comp. Time with respect to Reference Solution (%)

1 RDNM with 20  20 1 RDNM with 15  15 1 RDNM with 10  10 1 RDNM with 5  5

100 (Ref.) 92.2 87.5 84.4

the system is modified by adding a stiffness of k ¼ 30000N=m between the 10th DOF and ground, and as a result, the first and the second natural frequencies of the linear system become close to each other. As the 10th DOF is almost stationary at the second mode and is not affected by the addition of the new spring, the response at the 5th DOF is studied again for this modified system. Fig. 7 shows the displacement amplitude of the 5th DOF for the system with close modes. The first mode shifts considerably towards the second one; whereas, the change in the second resonance frequency is not significant. However, compared to the result of the previous case, a significant nonlinear behavior is observed around the second resonance region, which clearly shows strong modal coupling. The reference nonlinear response of the 5th DOF is obtained by including all modal contributions, and then compared with the responses calculated by using 1 RDNM, 2 RDNMs and 10 LMs around the first and the second resonance regions in Figs. 8

12

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

Fig. 7. Displacement amplitudes of the 5th DOF in stiff system.

and 9, respectively. It is seen from the figures that using 1 RDNM does not capture the nonlinear behavior accurately, especially at out of resonance regions due to strong modal coupling between the first two modes. Hence, one should increase the number of RDNMs used in the response calculations. When even only 2 RDNMs are used, however, the response is predicted very accurately in the frequency range considered. Using 2 RDNMs gives better results than using even 10 LMs. It should be noted that, using only 2 RDNMs reduces the number of nonlinear equations solved considerably, which has a significant effect on the computation time as well, in addition to improving accuracy. 3.2. Finite Element Model In this section, the applicability and the performance of the method proposed is demonstrated on a large-scale finite element model. The cantilever beam analyzed here (Fig. 10) is the same as the one used in [21]. It has extensions in all three perpendicular axes, and it is modeled with a commercial finite element software using 5400 DOFs. The first 30 modes of the finite element model are used to calculate the system response in the frequency range of interest for comparison purposes. The model is excited at the two nodes shown in Fig. 10. The components of the first excitation force in X, Y and Z directions are 5 N, 30 N and 30 N, respectively. The second excitation force has two nonzero components with amplitudes of 6 N and 3 N in Y and Z directions, respectively. It should be noted that with this type of forcing more modes of the system are excited which will be a good example to show the performance of the proposed modal superposition method with RDNMs. 1% proportional structural damping is assumed in the system. Only the fundamental harmonic is considered in the response calculations. Case Study 3: Piecewise Linear Stiffness, Dry Friction and Cubic Stiffness Together 9 nonlinear elements are attached to the structure. These elements are placed at different locations of the system so that capturing the nonlinear dynamics of the system requires high number of modes in the modal superposition approach. The types of the nonlinearities and their properties are given in Table 3. The piecewise linear stiffness nonlinearity used in the system is defined as follows:

_ Þ¼ f N ðxðtÞ; xðtÞ

8 > <

kp1 xðtÞ

kp2 xðtÞ þ ðkp1  kp2 Þd > : kp2 xðtÞ  ðkp1  kp2 Þd

xðtÞ < d & xðtÞ > d xðtÞ  d xðtÞ d

ð35Þ

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

13

Fig. 8. Displacement amplitude of the 5th DOF around the first resonance.

In this case study, it is aimed to demonstrate the superior performance of the proposed modal superposition method with RDNMs in the prediction of nonlinear periodic forced response of very large ordered systems; hence, the results obtained are compared with those of two recent studies. In a very recent study [36], several modal reduction methods for the periodic response analysis of nonlinear systems are compared. One of the outcomes of this study is that when the nonlinearity is not local, modal superposition method using linear modes (MSM-LM) [8] and modal superposition method using hybrid modes (MSM-HM), which is a method proposed by the authors of the present work [21], is computationally the most efficient ones, while maintaining the high accuracy. Furthermore, it is concluded that when a single harmonic is dominant, MSM-HM speeds up the computations providing the best method to apply. Therefore, the results obtained in this case study are compared with those obtained by MSM-LM and MSM-HM. In order to evaluate the performance of the proposed method, three different error criteria are defined to quantify errors in the results. The first one is integral error, which is defined as follows

Z Integral Error ¼

xf xi

!1=2 2

ðx  xex Þ dx

;

ð36Þ

where x, xex , xi and xf represent the displacement amplitude by using limited number of modal vectors, exact nonlinear displacement amplitude, initial and final frequencies of the frequency range of interest, respectively. Secondly, amplitude error, which is the error between the exact maximum displacement amplitude and the maximum displacement amplitude obtained by using limited number of modal vectors in the frequency range of interest. Mathematically it is expressed as

 ex  x  xmax  Amplitude Error ¼  max ex   100: x

ð37Þ

max

The last one is the frequency error, which is the error between the frequency of the maximum displacement amplitudes, i.e. resonance frequencies. It is expressed as follows

 ex  x  xmax  Frequency Error ¼  max ex   100:

xmax

ð38Þ

14

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

Fig. 9. Displacement amplitude of the 5th DOF around the second resonance.

Fig. 10. Finite Element Model.

Table 3 Nonlinear Element Types Used in the Third Case Study. Nonlinearity Type

Quantity

Properties

Piecewise Linear Stiffness Dry Friction Cubic Stiffness

3 3 3

kp1 ¼ 0; kp2 ¼ 500N=m; d ¼ 0:005m; kd ¼ 500N=m; lN ¼ 5N; kc ¼ 5  106 N=m3

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

15

Fig. 11. X component of the displacement amplitude of response node vs. frequency.

Fig. 11 shows the displacement amplitude of X component of the response node (shown in Fig. 10) for the underlying linear system and that of the nonlinear system using 30 modes of the linear system. It should be noted that since a relatively large number of DOF is considered in this case study, it makes sense to investigate the case where the nonlinearity affects the response at several modes. Comparison is performed around the first four resonance regions where the nonlinearity is the most dominant. Cubic stiffness nonlinearity is highly effective at the first resonance region as can be seen from Fig. 11. A shift of about 53% in the resonance frequency is also an indication of strong nonlinearity. The frequency shifts around the second and the third resonance regions are approximately 10% and 4%, respectively. The responses of the response node shown in Fig. 10 are calculated by using the modal superposition method with RDNMs and also by using MSM-LM, and are presented for the first four resonance regions, in Figs. 12–15, respectively. It is seen from these results that using even just 1 RDNM gives very accurate results. MSM-HM is also used for the calculation of nonlinear response; however, they are not presented in the figures given for brevity. Tables 4–7 give the error values as well as the computational time required by each method for different number of modes used. As can be seen from Tables 4– 7, similar accuracy can be obtained by using more number of LM or HM. Since MSM-HM provides more efficient basis than MSM-LM, two or three HMs are sufficient to be able to obtain accurate results compared to the use of 8 LMs. However, it can be observed from the tables that using a single RDNM gives the closest results to the exact response. It should also be noted that using RDNM decreases the computational time required by the solution significantly. The analyses are performed on a computer with Intel(R) Core(TM) i7 CPU 950 @3.07 GHz processor, 8 GB RAM and 64-bit operating system. Analyses are repeated for five times and an average of these five calculations is recorded as the computational time. For a fair comparison, total number of frequency steps is chosen as 200 ± 5 in the range of frequency interval analyzed in each case. It should be noted that due to the use of arc-length continuation, it is not possible to exactly satisfy the number of frequency steps. It can be seen that the computational time increases as the number of mode shapes used in the modal superposition approach increases, as expected. It can be seen from Tables 4–7 that modal superposition method with RDNM and MSM-LM utilizing 1 RDNM and 1 LM have similar computation times. Therefore, it can be concluded that computational time required by the EVP solution at each step using structural modification with Dual Modal Space method becomes insignificant compared to the nonlinear solution itself. However, using the proposed modal superposition method with RDNMs decreases the computational time required approximately 75% to 80% compared to using MSM-LM with 8 LMs, while maintaining the same level of accuracy in the response. Similarly, the total reduction obtained in computation time by the use of RDNM is approximately 95%, compared to using MSM-LM with 30 LMs. It should be noted that the method proposed is also significantly faster than MSM-HM, which has recently been shown to be the fastest nonlinear solution method for single harmonic analysis [36].

16

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

Fig. 12. X component of the displacement amplitude of response node vs. frequency around the first resonance region.

Fig. 13. X component of the displacement amplitude of response node vs. frequency around the second resonance region.

4. Conclusion In this study, a modal superposition method using response dependent nonlinear modes (RDNM) is proposed. The method is applicable to large ordered nonlinear MDOF systems, even with global nonlinearity. What makes this method novel is both the concept of RDNM, which is in a way similar to the well-known nonlinear normal modes (NNM), as well as the use of a structural modification with Dual Modal Space method, which further reduce the computational time of

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

17

Fig. 14. X component of the displacement amplitude of response node vs. frequency around the third resonance region.

Fig. 15. X component of the displacement amplitude of response node vs. frequency around the fourth resonance region.

RDNMs that need to be calculated at every frequency step in determining FRF curves of a nonlinear system. In this method, nonlinear internal forces are represented as a matrix, called nonlinearity matrix, multiplied by the displacement vector, using DFM. The nonlinearity matrix is an equivalent complex stiffness matrix representing the nonlinearity at a given displacement amplitude level. RDNMs of the system are calculated by solving the EVP corresponding to the displacement level and pattern for which the nonlinearity matrix is calculated. As the excitation frequency changes for a given harmonic forcing, the steady state response level and pattern will change resulting in a different equivalent stiffness matrix and therefore a different EVP. Structural modification with Dual Modal Space method is employed in order to decrease the computational effort in the evaluation of RDNM, which has upmost importance for realistic large DOF models of nonlinear structures. RDNMs are then used to predict steady state periodic forced response of the nonlinear structure. The method developed

18

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

is applicable to systems with any type of nonlinearity, making it very attractive for the determination of periodic forced response of real large nonlinear structures. Global nonlinearity is not a limitation for the method. Since RDNM forms a very efficient basis for the response space of a nonlinear system, the method developed uses minimum possible number of RDNMs, as a result of which the computational effort is reduced significantly. It is also observed in this study that the change of natural frequency with changing energy level, which is calculated using NNMs, can also be predicted by using RDNMs, although at some energy levels the results may not be so accurate. The preliminary studies on the relation between RDNM and NNM showed that it may be interesting to investigate the relation in a future study which may lead to the computation of NNMs by using RDNMs that can be calculated even for large ordered systems. In order to demonstrate the applicability and performance of the modal superposition method using RDNMs, case studies with four different nonlinear elements (gap, piecewise linear stiffness, dry friction and cubic stiffness nonlinearities) are presented. The case study with a 20 DOF lumped parameter system shows that by using even 1 RDNM, nonlinear response of the system can be captured very accurately when there is not modal interaction. Furthermore, it is seen that when there is strong modal interaction, the method proposed can still be employed by using more than one RDNM. It is observed in this case study that using only 2 RDNMs may be sufficient to predict the response very accurately, more accurately than using MSM-LM with even 10 LMs, which means significantly reduced computation time. It is also observed in further case studies using simple discrete systems that the variation of natural frequency with energy level in a nonlinear system can be approximately obtained by using RDNM concept. It is noted that it may be interesting to investigate the relation between RDNMs and widely used nonlinear normal modes (NNMs), which require a further study. Such a relation would make it much easier and more practical to calculate NNMs by employing RDNM concept, especially for large ordered systems. The case study with the finite element model having 5400 DOFs demonstrates the applicability of the method to large ordered systems with several different nonlinearities. Using even 1 or 2 RDNMs is sufficient to obtain very accurate results. This case study shows that the number of LMs in MSM-LM should be increased at least 4 to 8 times of the number of RDNMs used, depending on the nonlinear system considered, type of nonlinearities and their locations. It should be noted that, since

Table 4 Error analysis for the response node around the first resonance. Number of Modes

Integral Error

Amplitude Error (%)

Frequency Error (%)

Computational Time [s]

Exact Solution (30 LM) 1 LM 4 LM 8 LM 1 RDNM 3 HM 4 HM

– 47.9 29.5 7.76 6.87 14.9 10.4

– 41.6 0.96 0.12 0.05 0.95 0.43

– 5.59 0.81 0.05 0.03 0.21 0.08

212 9.2 24.7 47.7 11.5 23.7 28.8

Table 5 Error analysis for the response node around the second resonance. Number of Modes

Integral Error

Amplitude Error (%)

Frequency Error (%)

Computational Time [s]

Exact Solution (30 LM) 1 LM 4 LM 8 LM 1 RDNM 1 HM 2 HM

– 16.3 14.3 6.65 1.84 13.2 6.83

– 1.51 2.19 0.70 0.03 5.20 0.50

– 0.66 0.41 0.09 0.001 0.41 0.10

195 8.4 23.8 45.2 9.8 9.1 14.9

Table 6 Error analysis for the response node around the third resonance. Number of Modes

Integral Error

Amplitude Error (%)

Frequency Error (%)

Computational Time [s]

Exact Solution (30 LM) 1 LM 3 LM 6 LM 1 RDNM 2 HM 3 HM

– 3.22 0.73 0.47 0.75 1.05 0.59

– 9.14 1.02 0.40 0.003 4.42 0.49

– 0.23 0.06 0.02 0.004 0.03 0.03

192 8.7 18.0 32.9 9.6 13.9 18.6

19

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388 Table 7 Error analysis for the response node around the fourth resonance. Number of Modes

Integral Error

Amplitude Error (%)

Frequency Error (%)

Computational Time [s]

Exact Solution (30 LM) 1 LM 2 LM 4 LM 1 RDNM 1 HM 2 HM

– 0.52 0.07 0.12 0.37 0.57 0.07

– 2.02 0.24 0.56 0.0009 3.62 0.01

– 0.015 0.014 0.002 0.015 0.006 0.001

184 8.3 12.6 21.9 9.1 8.9 13.1

the computation of RDNM is performed using structural modification with Dual Modal Space method, the computational time required to solve EVP at each step becomes insignificant compared to the nonlinear solution, which makes the method developed practical for periodic forced vibration analysis of large order nonlinear FEM of real structures.

References [1] M.A. Dokainish, K. Subbaraj, A survey of direct time-integration methods in computational structural dynamics: I, explicit methods, Comput. Struct. 32 (6) (1989) 1371–1386. [2] M.A. Dokainish, K. Subbaraj, A survey of direct time-integration methods in computational structural dynamics: II, implicit methods, Comput. Struct. 32 (6) (1989) 1387–1401. [3] C.-H. Menq, J.H. Griffin, A comparison of transient and steady state finite element analyses of the forced response of a frictionally damped beam, J. Vib., Acoust., Stress, and Reliab. 107 (1) (1985) 19–25. [4] C.-H. Menq, J.H. Griffin, J. Bielak, The forced response of shrouded fan stages, J. Vib., Acoust., Stress, and Reliab. 108 (1) (1986) 50–55. [5] S. Setio, H.D. Setio, L. Jezequel, Modal analysis of nonlinear multi-degree-of-freedom structures, Int. J. Anal. Exp. Modal Analysis 7 (2) (1992) 75–93. [6] E. Budak, H.N. Özgüven, Iterative receptance method for determining harmonic response of structures with symmetrical non-linearities, Mech. Syst. Signal Process. 7 (1) (1993) 75–87. [7] O. Tanrikulu, B. Kuran, H.N. Özgüven, M. Imregun, Forced harmonic response analysis of nonlinear structures using describing functions, AIAA J. 31 (7) (1993) 1313–1320. [8] B. Kuran, H.N. Özgüven, A modal superposition method for non-linear structures, J Sound Vib. 189 (3) (1996) 315–339. [9] J.V. Ferreira, D.J. Ewins, Algebraic nonlinear impedance equation using multi-harmonic describing function, 15th Int Modal Anal. Conf. (1997) 1595. [10] J.J. Chen, C.H. Menq, Prediction of the resonant response of frictionally constrained blade systems using constrained mode shapes, Int. Gas Turbine Aeroengine Congress Exhibition 5 (1998) 1–10. [11] E. Cigeroglu, N. An, C.H. Menq, Forced response prediction of constrained and unconstrained structures coupled through frictional contacts, J. Eng. Gas Turbines Power 131 (2) (2009) 022505. [12] E.P. Petrov, A high-accuracy model reduction for analysis of nonlinear vibrations in structures with contact interfaces, J. Eng. Gas Turbines Power 133 (10) (2011) 102503. [13] L. Xu, M.W. Lu, Q. Cao, Nonlinear vibrations of dynamical systems with a general form of piecewise-linear viscous damping by incremental harmonic balance method, Phys. Letters A 301 (1) (2002) 65–73. [14] M. Guskov, J.J. Sinou, F. Thouverez, Multi-dimensional harmonic balance applied to rotor dynamics, Mech. Research Communications 35 (8) (2008) 537–545. [15] A. Grolet, F. Thouverez, On a new harmonic selection technique for harmonic balance method, Mech. Syst. Signal Process. 30 (2012) 43–60. [16] J.J. Hollkamp, R.W. Gordon, Reduced-order models for nonlinear response prediction: Implicit condensation and expansion, J Sound Vib. 318 (4–5) (2008) 1139–1153. [17] M. Krack, L. Panning-von Scheidt, J. Wallaschek, C. Siewert, A. Hartung, Reduced order modeling based on complex nonlinear modal analysis and its application to bladed disks with shroud contact, J. Eng. Gas Turbines Power 135 (10) (2013) 102502. [18] R.J. Kuether, B.J. Deaner, J.J. Hollkamp, M.S. Allen, Evaluation of geometrically nonlinear reduced-order models with nonlinear normal modes, AIAA J 53 (11) (2015) 3273–3285. [19] C. Gastaldi, S. Zucca, B.I. Epureanu, Jacobian projection reduced-order models for dynamic systems with contact nonlinearities, Mech. Syst. Signal Process. 100 (2018) 550–569. [20] C.S.M. Sombroek, P. Tiso, L. Renson, G. Kerschen, Numerical computation of nonlinear normal modes in a modal derivative subspace, Comput Struct 195 (2018) 34–46. [21] E. Ferhatoglu, E. Cigeroglu, H.N. Özgüven, A new modal superposition method for nonlinear vibration analysis of structures using hybrid mode shapes, Mech. Syst. Signal Process. 107 (2018) 317–342. [22] L. Peletan, S. Baguet, M. Torkhani, G. Jacquet-Richardet, A comparison of stability computational methods for periodic solution of nonlinear problems with application to rotordynamics, Nonlinear Dyn. 72 (3) (2013) 671–682. [23] R.M. Rosenberg, On nonlinear vibrations of systems with many degrees of freedom, Advances Appl. Mech. 9 (1966) 155–242. [24] S.W. Shaw, C. Pierre, Non-linear normal modes and invariant manifolds, Journal of Sound and Vibration 150 (1991) 170–173. [25] A.F. Vakakis, Non-linear normal modes and their applications in vibration theory: an overview, Mech. Syst. Signal Process. 11 (1997) 3–22. [26] G. Kerschen, M. Peeters, J.C. Golinval, C. Stéphan, Nonlinear modal analysis of a full-scale aircraft, J Aircraft 50 (5) (2013) 1409–1419. [27] R.J. Kuether, M.S. Allen, A numerical approach to directly compute nonlinear normal modes of geometrically nonlinear finite element models, Mech. Syst. Signal Process. 46 (2014) 1–15. [28] C. Touzé, Normal form theory and nonlinear normal modes: Theoretical settings and applications, Modal Anal. Nonlinear Mech. Syst. 555 (2014) 75– 160. [29] L. Renson, J.P. Noël, G. Kerschen, Complex dynamics of a nonlinear aerospace structure: numerical continuation and normal modes, Nonlinear Dyn. 79 (2) (2015) 1293–1309. [30] G. Haller, S. Ponsioen, Nonlinear normal modes and spectral submanifolds: existence, uniqueness and use in model reduction, Nonlinear Dyn. 86 (3) (2016) 1–42. [31] G. Kerschen, Computation of nonlinear normal modes through shooting and pseudo-arclength computation, Modal Analysis of Nonlinear Mechanical Systems, CISM International Centre for Mechanical Sciences, Springer, Vienna 555 (2014) 215–250. [32] Y.W. Luk, L.D. Mitchell, Implementation of the dual modal space structural modification method, 2nd Int, Modal Anal. Conf. 2 (1984) 930–936. [33] G. Von Groll, D.J. Ewins, The harmonic balance method with arc-length continuation in rotor/stator contact problems, J Sound Vib. 241 (2) (2001) 223– 233.

20

E. Ferhatoglu et al. / Mechanical Systems and Signal Processing 135 (2020) 106388

[34] P. Sundararajan, S.T. Noah, Dynamics of forced nonlinear systems using shooting/arc-jength continuation method-application to rotor systems, J. Vib. Acoust. 119 (1) (1997) 9–20. [35] E. Cigeroglu, H. Samandari, Nonlinear free vibration of double walled carbon nanotubes by using describing function method with multiple trial functions, Phys. E 46 (2012) 160–173. [36] E. Ferhatoglu, T. Dreher, E. Cigeroglu, M. Krack, H. Nevzat Özgüven, Comparison of linear and nonlinear modal reduction approaches, Nonlinear Struct. and Syst. 1 (2020) 229–233.