Nonlinear hybrid modal synthesis based on branch modes for dynamic analysis of assembled structure

Mechanical Systems and Signal Processing 99 (2018) 624–646

Contents lists available at ScienceDirect

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

Review

Nonlinear hybrid modal synthesis based on branch modes for dynamic analysis of assembled structure Xing-Rong Huang a,⇑, Louis Jézéquel a, Sébastien Besset a, Lin Li b, Olivier Sauvage c a

Laboratoire de Tribologie et Dynamique des Systèmes, Ecole Centrale de Lyon - 36, Avenue Guy de Collongues, 69134 Ecully, France School of Energy and Power Engineering, Beihang University, 37 Xueyuan Road, 100191 Beijing, China c PSA Group, Scientific and Future Technologies Department/StelLab, F-78943 Vélizy Villacoublay Cedex, France b

a r t i c l e

i n f o

Article history: Received 9 July 2016 Received in revised form 23 June 2017 Accepted 1 July 2017

Keywords: Nonlinear mode Branch mode Reduced model Steady state response Nonlinear interface Assembled structures

a b s t r a c t This paper describes a simple and fast numerical procedure to study the steady state responses of assembled structures with nonlinearities along continuous interfaces. The proposed strategy is based on a generalized nonlinear modal superposition approach supplemented by a double modal synthesis strategy. The reduced nonlinear modes are derived by combining a single nonlinear mode method with reduction techniques relying on branch modes. The modal parameters containing essential nonlinear information are determined and then employed to calculate the stationary responses of the nonlinear system subjected to various types of excitation. The advantages of the proposed nonlinear modal synthesis are mainly derived in three ways: (1) computational costs are considerably reduced, when analyzing large assembled systems with weak nonlinearities, through the use of reduced nonlinear modes; (2) based on the interpolation models of nonlinear modal parameters, the nonlinear modes introduced during the first step can be employed to analyze the same system under various external loads without having to reanalyze the entire system; and (3) the nonlinear effects can be investigated from a modal point of view by analyzing these nonlinear modal parameters. The proposed strategy is applied to an assembled system composed of plates and nonlinear rubber interfaces. Simulation results have proven the efficiency of this hybrid nonlinear modal synthesis, and the computation time has also been significantly reduced. Ó 2017 Elsevier Ltd. All rights reserved.

Contents 1. 2.

3.

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Theoretical aspects of reduced nonlinear modal synthesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1. Nonlinear modes based on the normal form theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2. Nonlinear modal parameters derived from nonlinear modal analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3. Mode synthesis of forced responses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4. Strategy guiding the reduced nonlinear modal synthesis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation and results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1. Analysis of the truncation effects of higher-order nonlinear modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2. Analysis of truncation effects on forced responses of internal modes, branch modes in the reduced basis . . . . . . . . . . .

⇑ Corresponding author. E-mail address: [email protected] (X.-R. Huang). http://dx.doi.org/10.1016/j.ymssp.2017.07.002 0888-3270/Ó 2017 Elsevier Ltd. All rights reserved.

625 626 626 628 628 629 631 633 636

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

4.

3.3. Analysis of nonlinear phenomena Conclusion . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgments . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . .

from a ...... ...... ......

modal point of view . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ............................................................... ............................................................... ...............................................................

625

642 645 645 646

1. Introduction In the context of engineering structural dynamics, the systems are generally large-sized and of nonlinear nature. A typical nonlinear industrial case consists of assembled structures containing nonlinearities located along interfaces between substructures. To analyze the nonlinear effects and calculate the steady state responses of these systems under a constrained computation time, fast and efficient methods are needed. To obtain a reduced-order model of industrial structures, the linear mode synthesis (LMS) theory has been developed into a powerful tool by various authors for analyzing mechanical systems. A brief review of dynamic substructuring can be found in Klerk et al. [26]. For systems with flexible physical interfaces, a clarification about component mode synthesis methods has been provided in Ohayon and Soize [33], namely: (1) attachment modes were first introduced to describe substructures by Hurty [18]; (2) a method combining static response modes (constraint modes) and free vibration modes of the substructure clamped on its interface (internal modes) was proposed by Craig and Bampton [12]; (3) branch modes, which are the first boundary modes obtained by condensing the modes of the structure on its interface thanks to the global structure’s static response modes, were investigated by Jézéquel et al. [21,9,16]; and (4) model reduction techniques for structural dynamics were compared in Besselink et al. [7]. These techniques are all limited to linear systems. Whereas the nonlinear effect on global structural dynamics cannot be neglected in most cases and linear assumptions in these cases may lead to considerable discrepancies. In order to solve the nonlinear problem, numerical time integration methods, such as the Newmark method and RungeKutta method, are often used [3]. These methods yield quite accurate results though appear to be time-consuming, especially when analyzing large nonlinear systems containing many degrees of freedom (DOFs) and when the steady state responses are required. To overcome this drawback, Rosenberg proposed the nonlinear normal mode (NNMs) concept to study nonlinear systems [37]. His premise was that resonance occurs in the neighborhood of the normal mode vibration regardless of whether the system is linear or nonlinear; moreover, the response around the main resonance can be represented by this main resonant mode response. The NNMs concept has been further developed by various authors based on theoretical aspects [32,5], numerical aspects [17,24,36] and experimental aspects [39,13,14,25,34,11]. A thorough review of theoretical developments of NNMs based on the determination of modal lines in configuration space and the computation of invariant manifolds of motion is conducted in Mikhlin and Avramov [29]. A certain summary of NNM applications is furthermore given by Avramov et al. [4]. A single nonlinear mode approach has been legitimized by Szemplinska-Stupnicka [41] to solve steady state responses of nonlinear differential equations by using approximate methods, e:g. the harmonic balance method, averaging method and asymptotic method [30]. In relying on this single nonlinear mode approach, Jézéquel and Lamarque have normal form theory to describe the forced response of harmonically excited systems for a two-DOF system [20]. Following this, Neild and Wagg [2] have then applied the normal form method to forced nonlinear vibration problems with the system equations expressed in the second-order form. Setio et al. [40] extended this normal form method to express general transient responses as the algebraic addition of nonlinear modal responses of general systems. These generalized nonlinear modal analysis methods are efficient in reducing the computation time when analyzing nonlinear systems, yet prove to be still time-consuming since the iterations are integrated into the numerical approaches. Reduction techniques are therefore required. Reduced-order modeling techniques relying on nonlinear modes are presented in Touzé and Amabili [42], Lülf et al. [28]. Kuether et al. have developed a reduction technique for geometrically nonlinear models [27]. A reduction method, based on the enrichment of reduction basis constituted of vibration modes with modal derivatives, is presented in Wu and Tiso [43] for flexible multibody systems. In general, these methods are devoted to analyzing complex nonlinear phenomena. While set up for weak, large nonlinear systems, the nonlinear problem may be simplified under appropriate assumptions and lead to a faster algorithm. Nonlinear models can subsequently be reduced by integrating reduction techniques that are analogous to those employed in the LMS method. Various reduction techniques in component mode synthesis can be used to extend the normal form theory [26]. The chosen reduction technique should comply with the concrete nonlinear problem. A discussion on all reduction techniques available has not be included herein, and special attention has been paid to Craig & Bampton reduction and double modal synthesis based on branch modes. Since the reduction techniques developed in LMS theory are often applied to analyzing engineering systems and since the NNMs obtained from free vibrations can be employed to calculate the forced responses, we propose herein to integrate reduction techniques with the NNMs concept in studying the dynamic performance of assembled systems. Before calculating nonlinear modes, the reduction technique is first applied to the nonlinear model in the case of free vibrations: physical displacements of the system are projected onto the generalized modal coordinates. A second reduction proceeds by selecting the dominant nonlinear modes obtained by solving the reduced order nonlinear problem with numerical iteration methods. The nonlinear modal parameters, i:e. mode shapes, natural frequencies and damping ratios, are deduced based on the

626

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

reduced nonlinear modes and then interpolated as a function of their corresponding modal amplitudes. Analyzing these modal parameters makes it possible to describe the system from a modal point of view; such an overview proves to be very useful when evaluating the dynamic performance of the system or identifying the responsible modes for nonlinear phenomena. Furthermore, forced responses can be computed by using the interpolation models of these modal parameters. This work is the extension of the double modal synthesis presented in Besset and Jézéquel [8], and it is devoted to the study of large assembled nonlinear mechanical systems. Our work program is outlined as follows. First, Section 2 will present the theoretical aspects of the reduced nonlinear modal synthesis: The calculation of nonlinear modes, based on normal form theory for slightly damped nonlinear systems, is presented; Nonlinear modal parameters are extracted from the nonlinear modes and then used to compute forced responses; Expressions for generating the steady state forced responses are given; Descriptions on how reduction techniques are integrated into the nonlinear modal analysis approach are provided. The proposed strategy will then be applied to a model composed of plates and nonlinear rubber material interfaces for validation purposes in Section 3. Lastly, some general conclusions and an outlook will be shared in Section 4. 2. Theoretical aspects of reduced nonlinear modal synthesis The systems under consideration are large assembled structures, with slight damping and weak nonlinearity continuously located in connections between substructures, as shown in Fig. 1. The governing motion equation is expressed as [30]:

€ þ Du_ þ Ku þ ~fðuÞ ¼ F cosðxt þ uÞ Mu

ð1Þ

where M is the mass matrix, D the damping matrix, K the stiffness matrix, and u the vector of unknown physical displacements. F cosðxt þ uÞ is the harmonic force applied on the system, with x being the excitation angular frequency and u the phase delay between forced response and excitation effort. ~fðuÞ is the nonlinear restoring force, which depends on the unknown physical displacements and renders the system nonlinear. The flow diagram illustrating the architecture of the proposed reduction technique is depicted in Fig. 2 to give a comprehensible global picture of the reduction strategy. The nonlinear modal analysis approach (NLMA) devoted to calculating the nonlinear modes is depicted in Section 2.1; the expressions of modal parameters containing essential nonlinear information of the system are given in Section 2.2; single nonlinear mode responses are computed based on the interpolation models of the modal parameters with respect to modal amplitude in Section 2.3; the strategy used for model reduction is presented in Section 2.4. The whole simulation process can be divided into two stages: offline simulation which permits obtaining the reduced-order model before conducting the nonlinear modal synthesis; online simulation which permits obtaining the nonlinear modes, nonlinear modal parameters and forced responses, based on the reduced-order model generated during the offline stage. 2.1. Nonlinear modes based on the normal form theory Assuming that the vibration of the nonlinear system in resonance condition can be approximated by the corresponding resonant mode and that all the non-resonant coordinates can be neglected except the single resonant coordinate [41], the nonlinear modes depending on the movement amplitude are used to minimize the coupling between mode shapes in modal ~ j is approached by: base. The stationary solution of the system in the resonance condition at x ¼ x

~ j cosðx ~ j tÞ uj ðtÞ  qj U

ð2Þ

Fig. 1. Substructures with nonlinear interfaces.

627

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

Fig. 2. Flow diagram of reduced nonlinear modal synthesis method.

where uj ðtÞ is the stationary solution in resonance condition. Notations analogous to those used in linear mode synthesis ~ j and the eigenvector for nonlinear mode-j is denoted (LMS) have been introduced herein: the eigenfrequency is denoted x ~ j; x ~ j are the unknown parameters depending on the modal amplitude qj . ~ j and U U Inserting Eq. (2) into Eq. (1) and taking the first harmonic component of the nonlinear restoring force ~fðuÞ, the nonlinear second order differential equation for mode-j is written by:

h

i ~ j Þ þ ~fðqj U ~ jÞ ¼ 0 ~ 2j M þ KÞðqj U ðx

ð3Þ

~ j Þ. Given that the nonDistinct from the linear modal analysis, Eq. (3) is a nonlinear problem with the nonlinear vector ~fðqj U ~ j are of a ~ j Þ is involved in the motion equation, both the natural frequencies x ~ j and modal shapes U linear restoring force ~fðqj U nonlinear nature and depend on the modal amplitude qj . For the sake of simplicity, this dependence on qj has neither been denoted herein nor in the following discussion. The eigensolutions of this nonlinear problem can only be solved by numerical procedures. The nonlinear frequency and nonlinear mode are obtained as a function of qj by progressively increasing the modal amplitude. The equivalent linearization method of Kryloff and Bogoliubov [10] is employed to obtain these unknown parameters. The idea of the equivalent linearization approach is to replace the nonlinear differential equation by a linear equation expressed ~ is the stiffness matrix of the equivalent linear system. in the following form, where K

h

i ~ ~ ~ 2j M þ KÞðq ðx j Uj Þ ¼ 0

ð4Þ

where

~U ~ j qj ¼ KU ~ j qj þ ~fðqj U ~ jÞ K

ð5Þ

The unknown parameters are determined by minimizing the difference between these two systems. The residue function between the equivalent linear equation and the nonlinear differential equation is expressed as:

ðqj Þ ¼

h

i h i ~ ~ ~ ~ ~ ~ 2j M þ KÞðq ~2 ðx j Uj Þ  ðxj M þ KÞðqj Uj Þ þ fðqj Uj Þ

ð6Þ

The nonlinear modes are solved by minimizing the residue function ðqj Þ according to the Newton-Raphson procedure

~ j ). The linear modes (xj ; Uj ) are set as initial conditions ~ j; U outlined below. The unknown parameters of this equation are: (x ~ k ) corresponding to qk are set as initial ~ k; U for the nonlinear problem corresponding to q0 ¼ 0. The solution of iteration-k (x j

j

j

j

~ kþ1 ) corresponding to qkþ1 . Assuming the model ~ jkþ1 ; U conditions so as to calculate the solution of the next iteration-k þ 1 (x j j

628

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

dimension to be of size N, Eq. (6) is a N degree equation. While N þ 1 unknowns are involved in the numerical approach to solve the nonlinear problem: 1 eigenvalue and N eigenvector components. A normalization condition with respect to the mass matrix is imposed to solve this problem:

~ i Þ MU ~ j ¼ 1 ði ¼ jÞ ðU T

ð7Þ

2.2. Nonlinear modal parameters derived from nonlinear modal analysis Once the nonlinear modes are obtained, they should be interpolated as a function of the corresponding modal amplitude for further calculations. In order to reduce the unknown parameters to be retained, it is more convenient to express the non~jk . It should ~ j as a linear combination of N l associated linear modes and N l nonlinear participation coefficients b linear mode U be noted that an implicit reduction is taking place here, given that the impacts of higher order linear normal modes (superior to N l ) on this nonlinear mode are being neglected when using this representation.

~j ¼ U

Nl X ~jk Uk b

ð8Þ

k¼1

~ j is expressed as: The participation of linear normal mode Uk in the nonlinear mode U

~jk ¼ ðUk ÞT MU ~j b

ðj – kÞ

ð9Þ

The linear normal modes in U should satisfy the following relationship:

ðUi ÞT MUj ¼ dij ~ j is 1: As a consequence, the participation of linear normal mode Uj in the nonlinear mode U

~jj ¼ 1 ðk ¼ jÞ b The nonlinear perturbation effect has been integrated into the nonlinear natural frequency and nonlinear modes, thus resulting in a nonlinear stiffness matrix:

~U ~j ~ jÞ K ~ 2j ¼ ðU x T

ð10Þ

The nonlinear modal damping is given by:

~ ¼ ðU ~ j ÞT DU ~j d j

ð11Þ

The nonlinear modal damping factor is thus:

~ T ~ ~nj ¼ 1 ðUj Þ DUj ~ 2j 2 x

ð12Þ

Rather than constraining all components of nonlinear mode shapes, we propose conserving the participation coefficients ~jk and nonlinear natural frequency x ~ j as a function of modal amplitude qj . When q approaches 0, the nonlinear effort conb ~jk turns out to be 0 ~ j converge to the linear normal mode shapes Uj . Consequently, b verges to 0 and nonlinear mode shapes U ~jk with respect to modal amplitude qj are subsequently used ~ j and b for all k – j and 1 for k ¼ j. The interpolation models of x to calculate the forced responses in the following section. 2.3. Mode synthesis of forced responses When the modes are slightly coupled, the nonlinear forced responses can be obtained from a superposition approach of the forced responses induced by single nonlinear modes. According to Ritz’s approximation method and left multiplying Eq. ~ j ÞT , we have: (1) with ðU

~ jÞ ðU

T

  ~ U ~ j qj ¼ ðU ~ j ÞT F x2 M þ iD þ K

ð13Þ

The modal amplitude corresponding to the single nonlinear mode response of mode-j can thus be expressed as:

qj ¼

~f ðjq jÞ j j ~ ~ 2j ðjqj jÞ x þ idj ðjqj jÞ þ x 2

ð14Þ

where ~f j is the nonlinear modal force of mode j:

~f ðjq jÞ ¼ ðU ~ j ðjqj jÞÞT F j j

ð15Þ

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

629

~ ;x ~ ~ j; d The modal amplitude qj is a nonlinear function of modal parameters (l j ~ j ; f j ). These modal parameters are already interpolated with respect to qj in the last step. The determination of qj turns out to be a nonlinear scalar problem pertaining to excitation frequency x in Eq. (14). By using a superposition procedure [40], the steady state responses of the nonlinear system can be approximated as a ~ j depending on their modal amplitudes qj : linear combination of N r nonlinear modes U

uðtÞ 

Nr Nr X X ~ j cosðxtÞ uj ðtÞ  qj U j¼1

ð16Þ

j¼1

where N r is the number of selected nonlinear modes used to approximate the response of the nonlinear system. This approach is based on the assumption that the coupling between the nonlinear modes is neglected and the system is slightly damped. 2.4. Strategy guiding the reduced nonlinear modal synthesis The above analysis has been based on the full order model containing N DOFs. Both the determination of nonlinear modes in Section 2.1 and the forced responses in Section 2.3 pertain to iterative procedures. If fewer variables were involved in the iterations, then the numerical computation process would be accelerated, thus leading to a significant savings in computational resources. Transformation matrices, such as those with a reduced basis in the component mode synthesis, are employed herein to project physical coordinates onto the generalized modal coordinates prior to the calculation of nonlinear modes. The single nonlinear modal analysis approach described in Section 2.1 was further applied to the reduced model, upon which the modal coordinates of reduced model’s forced responses were obtained by superimposing the single nonlinear mode responses, as indicated in Section 2.3. Lastly, the physical coordinates of the forced responses are obtained by projecting the modal space onto the physical space by use of the reduced basis. To facilitate the numerical protocol in u; M; K; D and F, the nodes have been partitioned into the master DOFs and the slave DOFs. Master DOFs are those on the interfaces containing nonlinearities and conserved as physical DOFs; slave DOFs are those in the substructures that would be projected onto the modal space. The numbers of master DOFs and slave DOFs are respectively denoted N J and N S . In the following analysis, a first reduction is carried out by transforming the interior DOF displacements into generalized modal coordinates, while boundary DOF displacements remain the same in physical space. This step is conducted utilizing the first free vibration modes of the substructure clamped on its interface (denoted by internal modes) and the boundary node functions (denoted by constraint modes), referred to as Reduced Nonlinear Modal Synthesis based on Constraint modes and Internal modes (RNLMS-CI). A second reduction is applied by condensing the modes of the structure on its interface thanks to static response modes of the global structure (denoted by branch modes), referred to as Reduced Nonlinear Modal Synthesis based on Branch modes and Internal modes (RNLMS-BI). The RNLMS-CI is based on the Craig & Bampton method [15]. The idea is to transform the FE model from a set of physical coordinates to a hybrid set of physical displacements on the boundary and modal coordinates at the interior DOFs. The physical displacements of the interior DOFs are composed of two parts: the rigid body displacements due to the boundary DOFs (WSJ uJ ) and the displacements relative to the fixed interface base (UII qI ).

uS ¼ WSJ uJ þ UII qI

ð17Þ

Matrix WSJ represents the rigid body displacements of the interior DOFs due to unit motion at one of the boundary DOFs. The matrix is determined by fixing all boundary DOFs and limiting consideration to the static problem. The governing equation in this case reduces to:

KSJ uJ þ KSS uS ¼ 0

ð18Þ

The set of internal displacements u1S , due to unit displacement of boundary DOFs when releasing the boundary DOFs one by one, may be solved as:

u1S ¼ K1 SS KSJ uJ ¼ WSJ uJ

ð19Þ

WSJ ¼ K1 SS KSJ

ð20Þ

with

The internal modes matrix UII is formed by eigenvectors UIj to the following equation, when the substructures are clamped on their interfaces with no forces acting on the interior points:

ðx2Ij MSS þ KSS ÞUIj ¼ 0 The set of internal displacements

ð21Þ u2S

relative to the fixed interface base is thus:

630

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

u2S ¼ UII qI

ð22Þ

The RNLMS-CI method entails selecting the primary internal modes: the physical displacements of slave DOFs are expressed as the addition of a matrix of static response modes (WSJ ) multiplied by the physical coordinates of master DOFs (uJ ), plus a reduced fixed interface modes matrix (UII ) multiplied by the modal coordinates of slave DOFs (qI ). The selected fixed interface mode number is denoted N I . In applying the reduced basis (TCI ), the physical displacements of the system is formed as:

 u¼

uJ



 ¼ TCI

uS

uJ

 ð23Þ

qI

where

 TCI ¼ ½ TC

TI  TC ¼

IJJ



WSJ

 TI ¼

0 UII

 ð24Þ

The above reduction base serves to reduce internal DOFs, while boundary DOFs remain the same. Boundary DOFs can also be reduced by using branch mode analysis [22,16]: whereby branch modes are introduced to condense the modes of the global structure on its interfaces. The matrix of branch modes is the solution of the global structure eigenvalue problem projected on constraint modes (TC ). The master DOFs involving nonlinearities are thus also projected onto the modal space. The selected branch mode number is denoted N B . Branch modes matrix XB is formed by eigenvectors XBj to Eq. (25).

ðx2Bj MB þ KB ÞXBj ¼ 0

ð25Þ

where

 MB ¼ ðTC ÞT MTC ;

KB ¼ ðTC ÞT KTC ;

TC ¼

IJJ



WSJ

WSJ ¼ K1 SS KSJ

The transformation matrix TBI is applied as follows:



u¼

uJ uS





¼ TBI

qB qI



ð26Þ

where

 TBI ¼ ½ TB

TI  TB ¼

XB WSJ XB



 TI ¼

0



UII

In the following analysis, the two transformation matrices TCI and TBI are both denoted Tr . Both ½qC ; qI  and ½qB ; qI  are denoted qr . The reduced motion equation is expressed as:

~ 2r Mr þ iDr þ Kr Þqr þ ðTr ÞT ~fðTr qr Þ ¼ 0 ðx

ð27Þ

where the reduced mass matrix, damping matrix and stiffness matrix are:

Mr ¼ ðTr ÞT MTr ;

Dr ¼ ðTr ÞT DTr ;

Kr ¼ ðTr ÞT KTr

By applying Ritz’s approximation based on a single nonlinear mode approach, the modal coordinates of forced responses are obtained by the superposition of the first N r nonlinear modes:

qr 

Nr Nr X X ~ rj qrj  qj U j¼1

ð28Þ

j¼1

~ rj Þ, which are calculated using the same procedure as ~ rj ; U In this way, we’ve introduced the reduced nonlinear modes ðx described in Section 2.1. The modal coordinates of forced responses are derived by the same rule presented in Section 2.3. ~S . ~ J and U The reduced nonlinear mode shape can be separated into two parts: U

" ~ rj ¼ U

~J U rj

rj

#

rj

~S U rj

~ J allows studying the contributions of linear constraint modes or linear branch modes in Analyzing the components in U rj ~ S yields the contributions of internal modes in the nonlinthe nonlinear mode-j; moreover, evaluating the components in U rj

ear mode-j. It should be noted that even the nonlinearities are only located on the interface, the components representing internal modes contributions also vary with qj . This finding can be explained by the way in which the reduced basis has been integrated into the nonlinear problem: since WSJ in RNLMS-CI or WSJ XB in RNLMS-BI is nonzero in the reduced basis, nonlin~ S through the coupling term in the reduced mass matrix Mr . earity would be transmitted to the last N I components in U rj

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

631

Let’s note that only N J þ N I variables are involved in the iterations when deducing the eigenvectors of the reduced nonlinear problem with RNLMS-CI, instead of computing N J þ N S variables using the NLMA method. Furthermore, N B þ N I variables are part of the iterative process when employing the RNLMS-BI, rather than introducing N J þ N S variables into the nonlinear problem by applying the RNLMS-CI method. When N B ¼ N J , the RNLMS-BI method is equivalent to the RNLMSCI method. From the full model described in Section 2, when N I ¼ N S , the RNLMS-CI method is equivalent to the NLMA method applied on the full model. The algorithm for calculating forced responses of the nonlinear system by the proposed reduced nonlinear modal synthesis approach is outlined as follows: Step 1: Construct the Finite Element (FE) model of the nonlinear system.

€ þ Du_ þ Ku þ ~fðuÞ ¼ F cosðxt þ uÞ Mu Step 2: Perform reductions on the full-order FE model, e:g. Craig & Bampton reduction technique, Branch mode analysis [16] or another reduction technique.

u ¼ Tr qr ~ rj Þ of the reduced model by implementing Kryloff and Bogoliubov’s ~ 2rj ; U Step 3: Compute the reduced nonlinear modes ðx equivalent linearization method with the Matlab solver.

~ rj þ ðTr ÞT ~fðTr qj U ~ rj Þ ¼ 0 ~ 2rj Mr þ iDr þ Kr Þqj U ðx ~rj ; x ~ rj , and ~f rj from the reduced nonlinear modes. Step 4: Extract nonlinear modal parameters d

~ ¼ ðU ~ rj ÞT Dr U ~ rj d rj

~f ¼ ðU ~ rj ÞT F rj

Step 5: Interpolate the extracted nonlinear modal parameters as a function of modal amplitude. Step 6: Solve the modal amplitude qj of the single nonlinear mode response by using the interpolated functions in Step 5 along with the Matlab solver and compute the generalized modal coordinates qrj .

~ rj qrj ¼ qj U

qj ¼

~f ðjq jÞ j j ~ ~ 2j ðjqj jÞ x þ idj ðjqj jÞ þ x 2

Step 7: Project the generalized modal coordinates onto the physical space with the reduced basis used in Step 2.

uj ¼ Tr qrj Step 8: Superimpose the physical coordinates of steady state responses of the selected nonlinear modes; the steady state responses of the nonlinear system are then obtained.

u

Nr X uj j¼1

Fig. 3. Illustration of the assembled system composed of Kirchhoff plates and nonlinear rubber layer interface.

632

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

Fig. 4. Rubber layer model.

Fig. 5. Finite element model of the assembled system composed of Kirchhoff plates and nonlinear rubber layer interface.

3. Simulation and results For purpose of description and validation of the proposed strategy, a model containing two Kirchhoff plates and one rubber layer interface has been investigated, as shown in Fig. 3. Plate 1 on the left is denoted by S1 : its length a1 is 0.66 m, the width b1 ¼ 0:6 m and the thickness e is 0.002 m; Plate 2 on the right is denoted by S2 : its length a2 is 0.44 m, the width b2 ¼ 0:6 m and thickness e ¼ 0:002 m; The rubber layer interface is denoted by J, with a thickness (h) of 0.002 m and a width (l) of 0.003 m: the rubber material is padded on both the top and bottom of the plates so as to ensure the connections between S1 and S2 . The thickness and width variation of the rubber layer acts upon the change in its stiffness value. Fig. 4 shows the details of this rubber layer model. The model is excited with a tire balance type loading [1] F ¼ mx2 R cos xt at a frequency x, with m ¼ 6 g, and R ¼ 4 cm. This loading corresponds to the real case when applying excitation on rotating machinery with increasing speed and has a lot of applications. The excitation is located at x ¼ 0:44 m and y ¼ 0:2 m, hence lying on S1 . The investigated excitation frequency band is ½3; 50 Hz. Structural damping has been integrated into the model, with a damping ratio of 0.02 in the nonlinear interface J, and 0.01 in S1 and S2 . The FE model has been built with a Structural Dynamics Toolbox (SDTools) [6]. To compute the forced responses of the FE model, Runge-Kutta’s time integration method takes a considerable amount of time to obtain the steady state response. As a compromise between the accuracy of the reduced model and the computation time, 4 elements per wavelength were utilized for the mesh size. The entire model contains 72 DOFs and 15 elements, of which 9 are in plate 1 on the left (i:e. S1 : a ¼ 0:66 m; b ¼ 0:6 m; e ¼ 0:002 m), 6 in plate 2 on the right (i:e. S2 : a ¼ 0:44 m; b ¼ 0:6 m; e ¼ 0:002 m), and 4 elements in the joints on the interface (J : e ¼ 0:002 m; l ¼ 0:003 m), as shown in Fig. 5. Rubber bushings are extensively used to link parts in a vehicle chassis, thus making it possible to filter noise and vibration [35]. For the sake of simplicity, scalar springs have been introduced to model the rubber layer, which is considered as reasonable for describing the stress-strain relationship of the rubber specimen. Since the rubber layer mass is negligible compared to the plate mass and since the rubber layer stiffness is far less than the plate stiffness. The rubber layer material is

633

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

assumed to be orthotropic, in considering the Young’s modulus value (E) along connection axis x and shear modulus (G) along both the y and z axes. This rubber material can be deemed as elastic for small deformations. To ensure the quasiincompressibility of the rubber material in volume term, a Poisson’s ratio value of 0.49 has been adopted. Young’s modulus varies with carbon content in the rubber. The Young’s modulus value adopted in this paper agrees with that of a real rubber material in order of magnitude [19]. The equivalent linear element stiffness (i:e. kx ; ky , and kz ) are given by:

8 dy > < kx ¼ E h; ky ¼ G dy h; > : kz ¼ G dy h:

ð29Þ

where dy is the mesh width. Nonlinearities are continuously distributed along the interface. The nonlinear restoring forces are assumed to lie in the translation direction of z. Cubic nonlinearity is one of the most popular means to describe the nonlinear effect on the dynamic behavior of the system, as revealed in many works [38,31,23]. The nonlinearities included herein are Duffing oscillators in the form of a cubic nonlinear restoring force that depends on the relative displacement between the nodes on the interface. The cubic nonlinearity ratio is denoted a, and the nonlinear restoring force is given by:

~fðuÞ ¼ aDu3

ð30Þ

Further details on the geometric dimension and material properties are provided in Table 1. 3.1. Analysis of the truncation effects of higher-order nonlinear modes By applying this nonlinear modal analysis on the full FE model, steady state responses of the model within the targeted frequency band (i:e. ½3; 50 Hz) can be compared to those output by Runge-Kutta’s classical time integration method. Given the spatial limitations, responses in z direction of the excitation point will be analyzed herein. The reduction coefficient, denoted by r c , is a common criterion in the mode synthesis method. It is defined as follows: if the maximum frequency of the targeted frequency band is 50 Hz, then the modes less than r c  50 Hz are all to be selected when describe the model performance. A convergence study of forced responses to the reduction coefficient has been conducted for the linear reduced-order model. A sensitivity analysis of the forced responses to the nonlinear mode reduction coefficient has been performed here for rc ¼ 1; rc ¼ 1:5; rc ¼ 2 and r c ¼ 4 (see Fig. 6). The forced response curves corresponding to the excitation frequency sweeping across the targeted frequency band from 3 Hz to 50 Hz are marked in solid lines, while those as the excitation frequency sweeps across the targeted band from 50 down to 3 Hz in the reverse direction are plotted in dotted lines. It can be observed that r c ¼ 2 appears to suffice for the study of dynamic behavior up to 50 Hz, since the response curves obtained by the NLMA method with r c ¼ 2 show only small differences with respect to the response curves given by r c ¼ 4. A reduction coefficient of rc ¼ 2 for nonlinear modes has thus been assigned to construct the forced responses. Fig. 7 compares the response curves obtained by summing the steady state responses of the first 18 single nonlinear modes with those obtained by the classical Runge-Kutta’s method. Two curves were generated for each method: as the excitation frequency sweeps across the targeted frequency band in ascending direction from 3 Hz up to 50 Hz, and as this excitation frequency sweeps across the targeted band in descending direction from 50 down to 3 Hz in the opposite direction. The amplitude of steady state responses is plotted on the upper subfigure, and the phase on the lower subfigure. These results demonstrate the effectiveness and accuracy of the NLMA method. The absolute and relative errors for the peak amplitude and corresponding resonance frequency of modes 5, 7 and 9 are outlined in Table 2 when the excitation frequency sweeps from 50 Hz to 3 Hz, and in Table 3 when the excitation frequency sweeps from 3 Hz to 50 Hz. By comparing Tables 2 and 3, it is observed that the errors in both frequency and peak amplitude

Table 1 Outlines of model geometry and meterial properties. SubStructure

Physical quantity

Plates

Material Thickness Poisson’s ratio Density Young’s modulus Shear modulus

Joint

Thickness Length Nonlinear coefficient Poisson’s ratio Young’s modulus Shear modulus

Unit mm kg/m3 Pa Pa mm mm

Pa Pa

Value Steel 2 0.285 7800 2.10E+11 8.17E+11 2 3 0.9 0.49 1.0E+08 3.4E+07

634

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646 −3

3.5

x 10

r =1 Raise

Physical Displaecment (m)

c

rc=1 Decrese

3

r =1,5 Raise c

r =1,5 Decrese c

2.5

r =2 Raise c

r =2 Decrese c

r =4 Raise

2

c

r =4 Decrese c

1.5 1 0.5 0

0

5

10

15

20

25

30

35

40

45

50

Frequency (Hz) 4

Phase

3 2 1 0

0

10

20

30

40

50

60

70

80

Frequency (Hz) Fig. 6. Comparison of the forced responses obtained by NLMA method applied on the full FE model, corresponding to different r c . r c corresponds to the reduction coefficient on nonlinear modes. Results corresponding to rc ¼ 1; r c ¼ 1:5; r c ¼ 2 and r c ¼ 4 are compared herein.

Physical Displaecment (m)

NLMA Raise

NLMA Decrease

4

Runge−Kutta Raise

Runge−Kutta Decrease

Excitation DDL32

−3

x 10

3 2 1 0

0

5

10

15

20

25

30

35

40

45

50

35

40

45

50

Frequency (Hz) 4

Phase

3

2

1

0

0

5

10

15

20

25

30

Frequency (Hz) Fig. 7. Comparison of the forced responses obtained by nonlinear modal analysis corresponding to r c ¼ 2 and those obtained by Runge-Kutta method.

635

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

Table 2 Comparison of natural frequency and forced response at the peak amplitude of modes 5, 7 and 9 using nonlinear modal synthesis and Runge–Kutta method when the excitation frequency sweeps from 50 Hz down to 3 Hz. 50 Hz ! 3 Hz

Mode

5

7

9

Natural frequency

NLMA Frequency (Hz) Runge–Kutta Frequency (Hz) Absolute Error (Hz) Relative Error (%)

27.36 27.33 0.03 0.11

35.21 35.16 0.05 0.14

47.54 47.33 0.21 0.44

Peak amplitude

NLMA Amplitude (mm) Runge-Kutta Amplitude (mm) Absolute Error (mm) Relative Error (%)

0.74 0.73 0.01 1.4

0.76 0.84 0.08 9.5

0.89 0.77 0.12 13

Table 3 Comparison of natural frequency and forced response at the peak amplitude of modes 5, 7 and 9 using nonlinear modal synthesis and Runge-Kutta method when the excitation frequency sweeps from 3 Hz up to 50 Hz. 3 Hz ! 50 Hz

Mode

5

7

9

Natural frequency

NLMA Frequency (Hz) Runge-Kutta Frequency (Hz) Absolute Error (Hz) Relative Error (%)

27.74 27.52 0.22 0.8

37.26 36.83 0.43 1.2

48.29 47.59 0.7 1.4

Peak amplitude

NLMA Amplitude (mm) Runge-Kutta Amplitude (mm) Absolute Error (mm) Relative Error (%)

1.02 0.92 0.1 10

2.4 2.2 0.2 8

1.10 1.09 0.01 0.9

Table 4 Comparison of natural frequency and forced response of mode 7 around 35 Hz with nonlinear modal synthesis and Runge-Kutta method for different excitation amplitudes, when the excitation frequency sweeps from 50 Hz down to 3 Hz. 50 Hz ! 3 Hz

m (g)

3

6

8

Natural frequency

NLMA Frequency (Hz) Runge-Kutta Frequency (Hz) Absolute Error (Hz) Relative Error (%)

34.73 34.68 0.04 0.12

35.21 35.16 0.05 0.14

35.49 35.49 0 0

Peak amplitude

NLMA Amplitude (mm) Runge-Kutta Amplitude (mm) Absolute Error (mm) Relative Error (%)

0.48 0.52 0.04 7.7

0.76 0.84 0.08 9.5

0.97 1.1 0.13 12

Table 5 Comparison of natural frequency and forced response of mode 7 around 35 Hz with nonlinear modal synthesis and Runge-Kutta method for different excitation amplitudes, when the excitation frequency sweeps from 3 Hz up to 50 Hz. 3 Hz ! 50 Hz

m (g)

3

6

8

Natural frequency

NLMA Frequency (Hz) Runge-Kutta Frequency (Hz) Absolute Error (Hz) Relative Error (%)

35.54 35.21 0.33 2.8

37.26 36.83 0.43 1.2

37.38 37.31 0.07 0.2

Peak amplitude

NLMA Amplitude (mm) Runge-Kutta Amplitude (mm) Absolute Error (mm) Relative Error (%)

0.76 0.72 0.04 5

2.4 2.2 0.2 8

3.6 3.7 0.1 3

are larger when the excitation sweeps in the ascending direction in comparison to those in the descending direction. The relative error in natural frequency is smaller than that in peak amplitude and discrepancy grows with the frequency value in general. For all these three modes, the relative errors in frequency are below 0:5% when the excitation sweeps in the descending direction and under 2% in the ascending direction. This demonstrates that the instable region can be quite accurately predicted using the proposed method. Moreover, since the largest error in the natural frequency and peak amplitude, predicted by the reduced model and the reference solution, is found to be round 35 Hz. Therefore, a quantitative analysis for several excitation amplitudes around

636

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

Table 6 Comparison of normal modes frequency of the full FE and the reduced model. Mode number

1

2

3

4

5

6

7

8

9

NLMA (Hz) RNLMS-CI (Hz) RNLMS-BI (Hz)

3.97 3.97 3.97

12.91 12.91 12.91

15.06 15.06 15.06

24.42 24.43 24.44

26.73 26.75 26.76

31.50 31.59 31.64

34.04 34.14 34.16

43.36 43.47 43.64

46.24 46.45 46.51

this mode is provided. Results corresponding to three excitation amplitudes are investigated: m ¼ 3 g; 6 g and 9 g in the tire balance loading F ¼ mx2 R cos xt. The relative errors are negligible in frequency, while there is an increase in peak amplitude when excitation amplitudes increase, as shown in Table 4 when the excitation sweeps in the descending direction. The relative errors decrease in frequency when excitation amplitudes increase, as seen in Table 5 when the excitation sweeps in the ascending case. 3.2. Analysis of truncation effects on forced responses of internal modes, branch modes in the reduced basis Before applying the reduced nonlinear modal synthesis (RNLMS) to the nonlinear model, a convergence study of forced responses to the reduction coefficient corresponding to component modes has been conducted on the linear model. Setting rc ¼ 2 as the reduction coefficient for all types of component modes, i:e. constraint modes, branch modes and internal modes, a quality analysis of the reduced model is evaluated as follows: The natural frequencies of the reduced model are compared with those of the full FE model. The natural frequency error is negligible up to 50 Hz, as displayed in Table 6. The modal assurance criterion (MAC) is also calculated in order to compare the reduced mode shapes and the exact mode shapes, with the MAC value being greater than 0.95 up to 50 Hz. Based on the above simulation results, a reduced nonlinear modal synthesis has been performed with rc ¼ 2 for all mode types involved in the transformation matrix T r used in the following analysis. Note that a higher r c value may be selected in case that a more accurate reduced model is needed. The reduction effect of internal modes and branch modes on the forced responses of the nonlinear system will be presented in the next section. By using the RNLMS-CI method, both the constraint modes and internal modes are involved in the transformation matrix. The reduction step is performed by condensing internal DOFs. Fig. 8 depicts the forced responses of the excitation point in z translation direction corresponding to various reduction coefficients: this Fig. 8 indicates good convergence of these responses around resonant frequencies when the selected internal mode number varies between 1 and 6, which has also revealed a reduction effect of the internal modes. It appears that rc ¼ 2 suffices for studying the dynamic performance of a nonlinear system up to 50 Hz. No higher-order internal modes need to be selected since the forced responses remain the same when more internal modes are retained with a reduction coefficient greater than 2. More precisely, the relative

−3

Physical Displaecment (m)

3.5

x 10

r =1 Raise c

rc=1 Decrese

3

rc=2 Raise

2.5

r =2 Decrese c

r =4 Raise c

2

r =4 Decrese c

r =6 Raise c

1.5

rc=6 Decrese

1 0.5 0 0

5

10

15

20

25

30

35

40

45

50

35

40

45

50

Frequency (Hz) Phase

4 3 2 1 0 0

5

10

15

20

25

30

Frequency (Hz) Fig. 8. Comparison of the forced responses obtained by RNLMS-CI method corresponding to different r c . r c corresponds to the reduction coefficient on internal modes. Results corresponding to r c ¼ 1; r c ¼ 2; r c ¼ 4 and r c ¼ 6 are shown herein.

637

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

Physical Displaecment (m)

RNLMS−CI Raise

RNLMS−CI Decrese

4

Runge−Kutta Raise

Runge−Kutta Decrese

Excitation − DDL32

−3

x 10

3 2 1 0

0

5

10

15

20

25

30

35

40

45

50

35

40

45

50

Frequency (Hz) 4

Phase

3

2

1

0

0

5

10

15

20

25

30

Frequency (Hz) Fig. 9. Comparison of the forced responses on excitation point obtained by Constraint NLMS method with r c ¼ 2 and by Runge-Kutta method. The point is located at ðx; yÞ ¼ ð0:22 m; 0:4 mÞ and its z-translation is shown in the figure.

Physical Displaecment (m)

RNLMS−CI Raise

RNLMS−CI Decrese

8

Runge−Kutta Raise

Runge−Kutta Decrese

Interface − DDL13

−3

x 10

6 4 2 0

0

5

10

15

20

25

30

35

40

45

50

35

40

45

50

Frequency (Hz) 4

Phase

2

0

−2

−4

0

5

10

15

20

25

30

Frequency (Hz) Fig. 10. Comparison of the forced responses on an arbitrary point of J obtained by Constraint NLMS method with r c ¼ 2 and by Runge-Kutta method. The point is located at ðx; yÞ ¼ ð0:66 m; 0 mÞ and its z-translation is shown in the figure.

638

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

Physical Displaecment (m)

RNLMS−CI Raise

RNLMS−CI Decrese

Runge−Kutta Raise

Runge−Kutta Decrese

Plate 1 − DDL48 0.02 0.015 0.01 0.005 0

0

5

10

15

20

25

30

35

40

45

50

35

40

45

50

Frequency (Hz) 0

Phase

−1

−2

−3

−4

0

5

10

15

20

25

30

Frequency (Hz) Fig. 11. Comparison of the forced responses on an arbitrary point of S1 obtained by Constraint NLMS method with r c ¼ 2 and by Runge-Kutta method. The point is located at ðx; yÞ ¼ ð0:22 m; 0:6 mÞ and its y-rotation is shown in the figure.

Physical Displaecment (m)

RNLMS−CI Raise

RNLMS−CI Decrese

Runge−Kutta Raise

Runge−Kutta Decrese

Plate 2 − DDL66 0.02 0.015 0.01 0.005 0

0

5

10

15

20

25

30

35

40

45

50

35

40

45

50

Frequency (Hz) 4

Phase

2

0

−2

−4

0

5

10

15

20

25

30

Frequency (Hz) Fig. 12. Comparison of the forced responses on an arbitrary point of S2 obtained by Constraint NLMS method with r c ¼ 2 and by Runge-Kutta method. The point is located at ðx; yÞ ¼ ð0:88 m; 0:4 mÞ and its x-rotation is shown in the figure.

639

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646 −3

x 10

Physical Displaecment (m)

4.5

rc=1 Raise

4

r =1 Decrese c

r =2 Raise

3.5

c

rc=2 Decrese

3

r =4 Raise c

2.5

r =4 Decrese c

all Raise all Decrese

2 1.5 1 0.5 0

0

5

10

15

20

25

30

35

40

45

50

35

40

45

50

Frequency (Hz)

Phase

4 3 2 1 0

0

5

10

15

20

25

30

Frequency (Hz) Fig. 13. Comparison of the forced responses obtained by RNLMS-BI method corresponding to different r c . r c corresponds to the reduction coefficient on branch modes. Results corresponding to rc ¼ 1; r c ¼ 2; r c ¼ 4 and no truncation are shown herein.

−3

Physical Displaecment (m)

3.5

x 10

r =2 Raise c

3

r =2 Decrese c

r =4 Raise c

2.5

r =4 Decrese c

r =6 Raise

2

c

r =6 Decrese c

1.5 1 0.5 0

0

5

10

15

20

25

30

35

40

45

50

35

40

45

50

Frequency (Hz)

Phase

4 3 2 1 0

0

5

10

15

20

25

30

Frequency (Hz) Fig. 14. Comparison of the forced responses obtained by RNLMS-BI method corresponding to different r c . r c corresponds to the reduction coefficient on internal modes. Results corresponding to rc ¼ 2; rc ¼ 4, and r c ¼ 6 are shown herein.

error in frequency between r c ¼ 1 and r c ¼ 2 is 3% around 27 Hz, 10% around 37 Hz; while no significant disparity appears between the forced responses obtained with r c ¼ 2 and r c ¼ 4. By employing the reduced nonlinear modal synthesis with constraint modes and internal modes (RNLMS-CI), the forced responses at four points of the model are calculated and compared to results obtained by the Runge-Kutta method, i:e.: an excitation point in Fig. 9, a randomly selected point on J in Fig. 10, a randomly selected point on S1 in Fig. 11, and a randomly selected point on S2 in Fig. 12. The accuracy of the method appears to be well within the purpose of practical engineering.

640

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

Physical Displaecment (m)

RNLMS−BI Raise

RNLMS−BI Decrese

4

Runge−Kutta Raise

Runge−Kutta Decrese

Excitation − DDL32

−3

x 10

3 2 1 0

0

5

10

15

20

25

30

35

40

45

50

35

40

45

50

Frequency (Hz) 4

Phase

3

2

1

0

0

5

10

15

20

25

30

Frequency (Hz) Fig. 15. Comparison of the forced responses on excitation point obtained by Branch NLMS method with r c ¼ 2 for internal modes and with r c ¼ 2 for branch modes and by Runge-Kutta method. The point is located at ðx; yÞ ¼ ð0:22 m; 0:4 mÞ and its z-translation is shown in the figure.

Physical Displaecment (m)

RNLMS−BI Raise

RNLMS−BI Decrese

6

Runge−Kutta Raise

Runge−Kutta Decrese

Interface − DDL13

−3

x 10

5 4 3 2 1 0

0

5

10

15

20

25

30

35

40

45

50

35

40

45

50

Frequency (Hz) 4

Phase

2

0

−2

−4

0

5

10

15

20

25

30

Frequency (Hz) Fig. 16. Comparison of the forced responses on an arbitrary point of J obtained by Branch NLMS method with r c ¼ 2 for internal modes and with r c ¼ 2 for branch modes and by Runge-Kutta method. The point is located at ðx; yÞ ¼ ð0:66 m; 0 mÞ and its z-translation is shown in the figure.

641

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

Physical Displaecment (m)

RNLMS−BI Raise

RNLMS−BI Decrese

Runge−Kutta Raise

Runge−Kutta Decrese

Plate 1 − DDL48 0.02 0.015 0.01 0.005 0

0

5

10

15

20

25

30

35

40

45

50

35

40

45

50

Frequency (Hz) 0

Phase

−1

−2

−3

−4

0

5

10

15

20

25

30

Frequency (Hz) Fig. 17. Comparison of the forced responses on an arbitrary point of S1 obtained by Branch NLMS method with r c ¼ 2 for internal modes and with r c ¼ 2 for branch modes and by Runge-Kutta method. The point is located at ðx; yÞ ¼ ð0:22 m; 0:6 mÞ and its y-rotation is shown in the figure.

Physical Displaecment (m)

RNLMS−BI Raise

RNLMS−BI Decrese

Runge−Kutta Raise

Runge−Kutta Decrese

Plate 2 − DDL66 0.02 0.015 0.01 0.005 0 0

5

10

15

20

25

30

35

40

45

50

35

40

45

50

Frequency (Hz) 4

Phase

2

0

−2

−4 0

5

10

15

20

25

30

Frequency (Hz) Fig. 18. Comparison of the forced responses on an arbitrary point of S2 obtained by Branch NLMS method with r c ¼ 2 for internal modes and with r c ¼ 2 for branch modes and by Runge-Kutta method. The point is located at ðx; yÞ ¼ ð0:88 m; 0:4 mÞ and its x-rotation is shown in the figure.

642

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

A further reduction is performed on the boundary DOFs using the reduced basis combined with branch modes and internal modes in the nonlinear modal synthesis (RNLMS-BI). To study the truncation effect of branch modes on the forced responses, stationary responses at the excitation point are plotted in Fig. 13. The reduced model investigated herein corresponds to r c ¼ 1; r c ¼ 2; r c ¼ 4 and no truncation for branch modes and r c ¼ 2 for internal modes. When observing the response curves around 27 Hz, it appears that a reduced model with rc ¼ 4 for branch modes provides sufficient convergence to study the stationary responses. To be more exact, the relative error in frequency between r c ¼ 2 and rc ¼ 4 is 1:4% around 27 Hz; while no significant disparity appears between the forced responses obtained with rc ¼ 4 and r c ¼ 6. The truncation effect of internal modes on the forced responses of the reduced model has also been analyzed when setting rc ¼ 4 for branch modes. Stationary responses at the excitation point corresponding to rc ¼ 2; rc ¼ 4 and rc ¼ 6 for internal modes are compared in Fig. 14. It appears that the reduced model with rc ¼ 2 for internal modes is adequate to study the stationary responses.The forced responses are observed to be more sensitive to the selected branch mode number than to the internal mode number, which can be explained by the fact that branch modes represent the overall model performance. The amplitude and phase of physical displacement for the 4 points obtained by RNLMS-BI with rc ¼ 4 for branch modes and r c ¼ 2 for internal modes are compared to those obtained using Runge-Kutta method: the excitation point (Fig. 15), a randomly selected point on J (Fig. 16), a randomly selected point on S1 (Fig. 17), and a randomly selected point on S2 (Fig. 18). These figures all suggest fairly good match between response curves obtained with the RNLMS-BI method and the temporal integration reference method, which proves the efficiency of the RNLMS-BI method. Whenever the amplitude reaches a peak on the upper subfigure, there is violent phase change occurring on the lower subfigure. In another finding, the model shows an obvious nonlinear performance when resonance is encountered, i:e. around 27 Hz, 36 Hz and 48 Hz. The response curve inflection leads to multivalued amplitudes and hence a ‘‘jump” phenomenon, which is actually a consequence of nonlinearities located along the joints.

3.3. Analysis of nonlinear phenomena from a modal point of view Fig. 19 displays the forced responses of the first 9 nonlinear modes. It is observed that mode 1 around 4 Hz and mode 2 around 13 Hz are nearly linear, while mode 7 around 36 Hz and mode 9 around 48 Hz show a relatively strong nonlinear performance.

mode−1

mode−2

−50

mode−3

−50

0

Raise Decrease

−50 −100

−100 −100

Physical Displacement (10lg(m))

−150

20

40

60

−150

mode−4

20

40

60

−150

mode−5

−100

−50

−150

−100

20

40

60

mode−6 −50 −100 −150

−200

20

40

60

−150

mode−7

20

40

60

−200

mode−8

20

40

60

mode−9

−100

−50

−60 −100

−80 −150

−100

−150

−120 20

40

60

−200

20

40

60

−200

Frequency (Hz) Fig. 19. Forced responses of the nine first nonlinear modes.

20

40

60

643

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

natural frequency variation

mode 1

mode 2

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

-0.1

0

0.02

0.04

-0.1

0

mode 7 0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0 0

0.02

0.04

mode 9

0.4

-0.1

0.02

0.04

-0.1

0

0.02

0.04

modal amplitude (m) Fig. 20. Nonlinear mode frequency variation in function of the modal amplitude for nonlinear modes-1; 2; 7; 9.

modal damping variation

mode 1

mode 2

0

0

-0.1

-0.1

-0.2

-0.2

-0.3

0

0.02

0.04

-0.3

0

mode 7

0

-0.1

-0.1

-0.2

-0.2

0

0.02

0.04

mode 9

0

-0.3

0.02

0.04

-0.3

0

0.02

0.04

modal amplitude (m) Fig. 21. Damping ratio variation in function of the modal amplitude for nonlinear modes-1; 2; 7; 9.

644

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

modal participations

mode 1

mode 2

1

1

0.5

0.5

0

0

-0.5

0

0.0 2

0.04

-0.5

0

0.02

mode 7

0.04

mode 9

1

1

0.5

0.5

0

0

-0.5

-0.5

0

0.02

0.04 0 modal amplitude (m)

0.02

0.04

Fig. 22. Modal participations of linear normal modes in function of the modal amplitude for nonlinear modes-1; 2; 7; 9.

Mode 1

Mode 2

8

8

6

6

4

4

2

2

0

0

−2

−2

−4

−4

−6

−6

Modal Participations

constraint modes

−8 0

0.005

0.01

0.015

0.02

−8 0

internal modes

0.005

Mode 7 8

6

6

4

4

2

2

0

0

−2

−2

−4

−4

−6

−6 0.005

0.01

0.015

0.02

0.015

0.02

Mode 9

8

−8 0

0.01

0.015

0.02

−8 0

0.005

0.01

Modal Amplitude (m) Fig. 23. Variation of modal components versus their corresponding modal amplitude for nonlinear modes-1; 2; 7; 9 obtained with reduced constraint NLMS strategy.

645

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646 Table 7 Comparison of retained modes number and CPU time with different methods. Method

Branch Modes number

S1 Internal Modes number

S2 Internal Modes number

Total modes number

CPU time (s)

NLMA RNLMS-CI RNLMS-CB

24 24 10

30 7 7

18 5 5

72 36 24

6394 4098 2636

Fig. 20 reveals how nonlinear natural frequency varies with its modal amplitude. It is observed that nonlinear natural frequencies all stem from the natural frequencies of the linear model without Duffing oscillators, and ultimately approximate the models’ linear natural frequencies when the nodes on the interface connecting each substructure are clamped to each other. The variation in nonlinear modal damping ratios as a function of modal amplitude is also plotted in Fig. 21. This variation tends to be nearly reciprocal to that of the nonlinear natural frequency, which can be explained by the mathematical expression of the modal damping ratio in Eq. (12). ~k1 ; b ~k2 ; b ~k7 ; b ~k8 where k ¼ 1; . . . ; 18) is The participation of linear normal modes in reduced nonlinear modes 1, 2, 7 and 9 (b captured in Fig. 22. The participation of the corresponding linear mode in the nonlinear mode always equals 1, which is in ~jj ¼ 1. agreement with the relation that b Figs. 20–22 explain why mode 1 around 4 Hz and mode 2 around 13 Hz are nearly linear in Fig. 7, while mode 7 around 36 Hz and mode 9 around 48 Hz exhibit a relatively strong nonlinear performance. The natural frequency, damping factor and mode participation parameters almost do not vary with modal amplitude for modes 1 and 2, while those for modes 7 and 9 vary substantially. Furthermore, with a transformation matrix composed of constraint modes and fixed interface modes, another modal overview can be derived on the nonlinear phenomena of the system. Similar to the above analysis, modes 1, 2, 7 and 9 have been investigated. Fig. 23 shows the variation trend of the components in the reduced nonlinear modes vs. the modal amplitude. The components pertaining to internal modes are marked by red lines, while those representing constraint modes are marked in blue in Fig. 23. Nonlinear modes vary with constraint mode components, yet they do not vary considerably with internal modes components. It can be concluded that constraint modes are responsible for the model’s nonlinear performance. To summarize, the RNLMS is proven to be accurate and efficient in analyzing assembled systems with nonlinearities continuously located at interfaces between substructures. For this assembled system composed of Kirchhoff plates and a nonlinear rubber layer, 50% of DOFs are truncated by applying the RNLMS-CI method and, furthermore, 67% of DOFs are truncated by applying the RNLMS-BI method, as indicated in Table 7. Simulations have been conducted on a server containing 32 Xeon (R) processors running at 2.9 GHz. Considering the 10,000 linearly-spaced modal amplitudes from 0 to 1, the CPU time consumed in computing one nonlinear mode is 6394 s when employing the NLMA method on the full FE model, 4098 s when applying RNLMS-CI, and 2636 s when using reduced RNLMS-BI. All reduction techniques serve to lower computation time, with the branch mode analysis providing even more savings since branch modes are also truncated compared to Craig & Bampton method. This finding may be especially valuable for analyzing large assembled systems with RNLMS-BI method, since DOFs on the interfaces represent the majority of DOFs and the RNLMS-BI method demonstrates considerable computational timesaving. 4. Conclusion The proposed hybrid modal synthesis approach for analyzing assembled nonlinear systems has combined the nonlinear mode concept with component mode synthesis based on branch modes. The dependences of nonlinear modes on modal amplitude were established. The nonlinear modal parameters extracted from the nonlinear modes were highlighted, which have been investigated to describe the nonlinear behavior of the system. The steady state responses were approximated by the proposed nonlinear mode approach and simulation results have proven the efficiency and accuracy of this derived strategy. Reduced nonlinear modal synthesis methods employing constraint modes or branch modes have been applied herein to an assembled system. This model has been reduced to a smaller size while maintaining sufficient information to describe model performance over the targeted frequency band. The truncation effects of internal modes and branch modes on forced responses have been investigated as well. A significant advantage of the reduced nonlinear modal synthesis lies in its computational timesaving in analyzing large assembled systems. The hybrid nonlinear modal synthesis developed in this work is applicable to optimizing weakly nonlinear systems during the pre-design process. Acknowledgments The authors gratefully acknowledge the financial support of the China Scholarship Council. This work program is ongoing in the context of OpenLab PSA ‘‘VAT@Lyon”.

646

X.-R. Huang et al. / Mechanical Systems and Signal Processing 99 (2018) 624–646

References [1] M. Agnieszka, Rotordynamics, CRC Press, 2005. [2] S.A. Neild, D.J. Wagg, Applying the method of normal forms to second-order nonlinear vibration problems, in: Proceedings of the Royal Society. No. October 2010, 2011, pp. 1141–1163. [3] U. Ascher, L. Petzold, Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, 1998. [4] K. Avramov, Y. Mikhlin, K.V. Avramov, Review of applications of nonlinear normal modes for vibrating mechanical systems, ASME Appl. Mech. Rev. (2013) 0–20. [5] J. Awrejcewicz, I. Andrianov, L. Manevitch, Asymptotic Approaches in Nonlinear Dynamics: New Trends and Applications, Springer, Berlin/Heidelberg, 2012. [6] E. Balmès, Structural Dynamics Toolbox. Tech. Rep., 2009. [7] B. Besselink, U. Tabak, A. Lutowska, N. Van De Wouw, H. Nijmeijer, D.J. Rixen, M.E. Hochstenbach, W.H.A. Schilders, A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control, 2013. [8] S. Besset, L. Jézéquel, Dynamic substructuring based on a double modal analysis, J. Vib. Acoust. 130 (1) (2008) 011008. [9] S. Besset, L. Jézéquel, Vibroacoustical analysis based on a multimodal strategy: triple modal synthesis, J. Vib. Acoust. 130 (3) (2008) 031009. [10] T.K. Caughey, Equivalent linearization techniques, J. Acoust. Soc. Am. 34 (12) (1962) 2001. [11] J. Ciambella, F. Vestroni, The use of modal curvatures for damage localizationin beam-type structures, J. Sound Vib. 340 (2015) 126–137. [12] R.-R. Craig, M.C.C. Bampton, Coupling of substructures for dynamics analyses, AIAA J. 6 (7) (1968) 1313. [13] D. Ewins, Modal Testing: Theory, Practice and Application, John Wiley & Sons, 2000. [14] C. Gibert, Analyse modale nonlineaire experimentale Ph.D. thesis, Ecole centrale de Lyon, 2001. [15] W.B. Haile, Primer on the craig-bampton method. Tech. rep., 2000. [16] X. Huang, L. Jézéquel, S. Besset, L. Li, Optimization of the dynamic behavior of vehicle structures by means of passive interface controls, J. Vib. Control I26 (2016). [17] E.S. Hugo Ramon, Non-Linear Modal Analysis Methods for Engineering Ph.D. thesis, Imperial College, London, 2004. [18] W.C. Hurty, J.D. Collins, G.C. Hart, Dynamic analysis of large structures by modal synthesis techniques, Special Issue Struct. Dyn. 1 (4) (1971) 535–563. [19] H.M. James, E. Guth, Theory of the elastic properties of rubber, J. Chem. Phys. 11 (10) (1943) 455. [20] L. Jézéquel, C. Lamarque, Analysis of non-linear dynamical systems by the normal form theory, J. Sound Vib. 149 (3) (1991) 429–459. [21] L. Jézéquel, H. Setio, Component modal synthesis methods based on hybrid models, Part I: Theory of hybrid models and modal truncation methods, ASME J. Appl. Mech. 61 (1) (1994) 100–108. [22] L. Jézéquel, H. Setio, Component modal synthesis methods based on hybrid models, Part I: Theory of hybrid models and modal truncation methods, ASME J. Appl. Mech. 61 (1) (1994). [23] D. Jiang, C. Pierre, S.W. Shaw, Nonlinear normal modes for vibratory systems under harmonic excitation, J. Sound Vib. 288 (4–5) (2005) 791–812. [24] G. Kerschen, M. Peeters, J.C. Golinval, a.F. Vakakis, Nonlinear normal modes, Part I: A useful framework for the structural dynamicist, Mech. Syst. Signal Process. 23 (2009) 170–194. [25] G. Kerschen, K. Worden, A.F. Vakakis, J.-C. Golinval, Past, present and future of nonlinear system identification in structural dynamics, Mech. Syst. Signal Process. 20 (3) (2006) 505–592. [26] D.D. Klerk, D.J. Rixen, S.N. Voormeeren, General framework for dynamic substructuring: history, review, and classification of techniques, AIAA J. 46 (5) (2008). [27] R.J. Kuether, M.S. Allen, Modal substructuring of geometrically nonlinear finite element models with interface reduction, AIAA J. 54 (2) (2015) 1–30. [28] F.A. Lülf, D.-M. Tran, R. Ohayon, Reduced bases for nonlinear structural dynamic systems: a comparative study, J. Sound Vib. 332 (15) (2013) 3897– 3921. [29] Y. Mikhlin, K. Avramov, Nonlinears normal modes for vibrating mechanical systems. review of theoretical developments, ASME Appl. Mech. Rev. (2011) 0–21. [30] A. Nayfeh, D. Mook, Nonlinear Oscillations, John Wiley & Sons, 1979. [31] A. Nayfeh, J. Nayfeh, On methods for continuous systems with quadratic and cubic nonlinearities, Nonlinear Dyn. 3 (1992) 145–162. [32] A.H. Nayfeh, B. Balachandran, Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods, John Wiley & Sons, 2007. [33] R. Ohayon, C. Soize, Clarification about component mode synthesis methods for substructures with physical flexible interfaces, Int. J. Aeronaut. Space Sci. 15 (2) (2014) 113–122. [34] M. Peeters, R. Viguié, G. Sérandour, G. Kerschen, J.C. Golinval, Nonlinear normal modes, Part II: Toward a practical computation using numerical continuation techniques, Mech. Syst. Signal Process. 23 (2009) 195–216. [35] G. Puel, B. Bourgeteau, D. Aubry, Parameter identification of nonlinear time-dependent rubber bushings models towards their integration in multibody simulations of a vehicle chassis, Mech. Syst. Signal Process. 36 (2) (2013) 354–369. [36] L. Renson, G. Kerschen, B. Cochelin, Numerical computation of nonlinear normal modes in mechanical engineering, J. Sound Vib. 364 (2016) 177–206. [37] R. Rosenberg, The normal modes of nonlinear n-degree-of-freedom systems, ASME 29 (1) (1962) 7–14. [38] S. Sangriyadi, Nonlinear Modal Analysis Ph.D. thesis, Ecole Centrale de Lyon, 1990. [39] H.D. Setio, S. Setio, L. Jezequel, A method of non-linear modal identification from frequency response tests, J. Sound Vib. 182 (1995) 336–341. [40] S. Setio, H.D. Setio, L. Jezequel, Modal analysis of nonlinear multi-degree-of-freedom structure, Int. J. Anal. Exp. Modal Anal. 7 (0) (1992) 75–93. [41] W. Szemplinska-Stupnicka, The resonant vibration of homogeneous non-linear systems, Int. J. Non-Linear Mech. 15 (4–5) (1980) 407–415. [42] C. Touzé, M. Amabili, Nonlinear normal modes for damped geometrically nonlinear systems: application to reduced-order modelling of harmonically forced structures, J. Sound Vib. 298 (2006) 958–981. [43] L. Wu, P. Tiso, Nonlinear model order reduction for flexible multibody dynamics: a modal derivatives approach, Multibody Syst. Dyn. (2016) 405–425.