Piecewise constrained optimization harmonic balance method for predicting the limit cycle oscillations of an airfoil with various nonlinear structures

Journal of Fluids and Structures 55 (2015) 324–346

Contents lists available at ScienceDirect

Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs

Piecewise constrained optimization harmonic balance method for predicting the limit cycle oscillations of an airfoil with various nonlinear structures Haitao Liao n Chinese Aeronautical Establishment, Beijing 100012, China

a r t i c l e in f o

abstract

Article history: Received 17 November 2014 Accepted 9 March 2015 Available online 21 April 2015

A method is proposed to calculate the periodic solutions of piecewise nonlinear systems. The method is based on analytical derivation of nonlinear multi-harmonic equations of motion. Since periodic variations of nonlinear forces are characterized by different states, the vibration cycle is broken into sequential transition intervals according to the instant sets of state transitions. Analytical formulations of the harmonic coefficients of the nonlinear forces and its derivatives with respect to the harmonic coefficients of displacements are developed. Sensitivities of the harmonic coefficients of periodic solutions are determined for constructing explicit expressions for vibration amplitude levels as a function of structural parameters. Numerical investigations of the limit cycle oscillations and its sensitivities of an airfoil with different piecewise nonlinearities have been performed. The results show that the developed method is capable of determining the periodic solutions and its sensitivities with respect to the structural parameters. In order to guarantee time continuity of the nonlinear force, for the hysteresis model it is not right to track the periodic solutions by using the preload or freeplay as the continuation parameters. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Piecewise nonlinearity Limit cycle Harmonic balance method Sensitivity

1. Introduction The influence of structural nonlinearities such as cubic, freeplay and hysteresis on the aeroelastic response has received significant attention in aerospace community (Lee et al., 2005). The two-dimensional wing section (Liu and Dowell, 2004; Liu et al., 2007), restrained by two springs in plunge and pitch motions, is a common structural model that is widely used for the aeroelastic modeling. Extensive investigations have been carried out on the aeroelastic response of cubic nonlinearity. Although the analysis of the dynamic behavior in aeroelastic systems with freeplay and hysteresis is hampered by the nonsmooth nonlinearity feature, research on these kinds of non-smooth dynamic systems has also been performed. For example, Zhang et al. (2013) studied the limit cycle oscillations (LCOs) of an airfoil with torsional freeplay and the effects of the freeplay value and contact stiffness ratio on the frequency, amplitude and phase of LCO are revealed. The averaging method and the Floquet theory are employed by Guo and Chen (2012) to investigate the LCOs of a 2-D of airfoil containing a combination of freeplay and cubic nonlinearities in supersonic flow. It has been demonstrated that the presence of freeplay

n

Tel.: þ86 10 15810533731. E-mail addresses: [email protected], [email protected]

http://dx.doi.org/10.1016/j.jfluidstructs.2015.03.008 0889-9746/& 2015 Elsevier Ltd. All rights reserved.

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

325

nonlinearity modifies the aeroelastic behaviors, resulting in interesting phenomena such as amplitude jump. In addition, it should be noted that there have been some studies for the analysis of aeroelastic dynamics under uncertainties, such as the work of Liao (2013). Several numerical methods have been developed to treat the piecewise nonlinearity airfoil system. In several works of Liu et al. (2002a, 2002b), a mathematical technique based on the point transformation method is presented to investigate the aeroelastic model with a freeplay nonlinearity as well as a hysteresis nonlinearity. Two formulations are developed to detect the periodic motions with harmonics, period doubling, chaotic motions and the coexistence of stable limit cycles. Moreover, the perturbation incremental method is introduced by Chung et al. (2007, 2009) to deal with structural nonlinearities like freeplay and hysteresis. With the help of the time discretization technique of the cyclic method, the incremental harmonic balance technique in conjunction with the Fourier expansion approach of hysteresis nonlinearity is used by Liu et al. (2012) to simulate the dynamics of an airfoil with a pitching hysteresis nonlinearity. In Chen and Liu (2014) and Cui et al. (2015), bifurcation and LCOs of the aeroelastic systems with freeplay nonlinearity are analyzed and a predictor–corrector algorithm based on the precise integration method is proposed to determine the switching points of the piecewise nonlinearity system. Chaotic response of a two-Dof airfoil in subsonic flow with structural nonlinearities was predicted by Li et al. (2012). The rational polynomials have been introduced to properly approximate both the freeplay and hysteresis nonlinearities. To overcome these computational limitations for approximation of the nonlinear forces in Li et al. (2012) and Liu et al. (2012), a new approach using piecewise integration is proposed in the present contribution. A number of investigators have contributed to the development of the harmonic balance method(HBM) or have used the technique to solve LCOs of wings and cascades. The HBM expands periodic response of a nonlinear system as truncated Fourier series, whose coefficients are determined by solving a set of nonlinear algebraic equations. Recently, with the aid of the symbolic calculation software Mathematica, the complex Fourier coefficients associated with geometric nonlinearity term is derived by Dai et al. (2014) and a fast HBM based on the optimal iterative algorithm and explicit Jacobian matrix has been developed to produce periodic solutions of an airfoil with cubic nonlinearity. Continuation of periodic orbits has now been recognized as the most efficient method for numerically computing the family of periodic orbits and their implementation is well documented within numerous standard methods for step continuation and adaptation (arclength, pseudo-arclength) (Sarrouy and Sinou, 2011). However, if some branches of solutions are disconnected (in the sense that they do not arise from bifurcation points of previous branches), then the continuation algorithm would fail to detect them and new tools would be needed to solve the algebraic system of equations induced by the HBM. A method named constrained optimization harmonic balance method (COHBM) has been developed by Liao (2013) and Liao and Sun (2013) for the numerical prediction of the periodic solutions of nonlinear mechanical systems. The method is essentially a nonlinear constrained optimization problem and consists of a two-step procedure: In the first step, the harmonic balance method and the stability Hill method are used to construct the nonlinear equality and inequality constraints. In the second step, the MultiStart optimization algorithm is used to find the optimization solution. The proposed methodology is applied to several nonlinear mechanical systems. Sensitivity analysis (Narimani et al., 2004) for the determination of the gradients of objective and constraints is the dominant process in the accuracy and computational time of many optimization problems. In addition, the sensitivity analysis provides the trends of variation of the objective and constraint functions versus system parameters efficiently. Therefore, the sensitivity analysis might be used as a useful design tool to evaluate the system response to changing parameters. Numerical methods such as finite difference technique have been used to determine the gradients of the objective and constraints function. However, finite difference evaluation of gradients may result in inaccurate optimal solution and is also computationally expensive. Therefore, an efficient analytical evaluation of sensitivity gradients is required. The main objective of this paper is to derive the explicit Jacobian matrix of HB algebraic system, which is used for sensitivity analysis. As applications to aeroelastic systems are targeted, three types of nonlinearities, cubic, freeplay and hysteresis, are considered. The rest of this paper is organized as follows: a detailed description of the airfoil model used for simulation is given in Section 2. The general formulation of the developed method for determining the periodic solutions is presented in Section 3. Demonstrations of the proposed method are then conducted in Section 4 where numerical examples are given. Finally, concluding remarks are presented and discussed in Section 5.

2. The airfoil model A two-degree-of-freedom pitch-and-plunge airfoil model is considered and the prototypical aeroelastic wing section is shown in Fig. 1. The plunge deflection is denoted by h, which is positive in the downward direction, and α is the pitch angle between the elastic axis and the positive nose up direction. The elastic axis is located at a distance ah b from the midchord, positive toward the trailing edge of the airfoil.

326

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

Fig. 1. Schematic diagram of the airfoil system.

The nonlinear aeroelastic equations of an airfoil motion can be found in the literature (Liu and Dowell, 2004; Lee et al., 2005). The aeroelastic motion of a 2-dof airfoil is governed by the following coupled bending–torsion equations:
 ω 1 ξ€ þ xα α€ þ 2ζ ξ ξ_ þ GðξÞ ¼  C L ðtÞ; πμ Un
2 xα € 2ζ α 1 2 _ € α þ MðαÞ ¼ C M ðtÞ; ð1Þ ξ þ α þ πμr 2α r 2α Un Un where t ¼ Ut=b (t is the real time), ξ ¼ h=b and U n ¼ U=ðbωα Þ with U as the flow speed are respectively the time, plunge displacement and velocity in the dimensionless form. The dot represents the derivative with respect to t. ωξ and ωα represent the natural frequencies of the uncoupled plunge and pitch oscillations respectively and ω ¼ ωξ =ωα . r α denotes the radius of gyration about the elastic axis. ζ ξ and ζ α are the damping ratios. In Eq. (1), the aerodynamic force C L ðtÞ is the lift force (positive up) and C M ðtÞ is the aerodynamic moment about the elastic axis (positive nose up). Based on the lab experiments, lift C L ðtÞ and pitch moment C M ðtÞ are given as follows (Liu and Dowell, 2004; Lee et al., 2005):   _ þ ð1=2  ah Þαð0Þ C L ðtÞ ¼ πðξ€  ah α€ þ αÞ _ þ2π αð0Þ þ ξð0Þ _ ϕðtÞ Z π   € þ ð1=2  ah ÞαðσÞ þ 2π ϕðt σÞ αðσÞ _ þ ξðσÞ € dσ; 0

   π π π _ þ ð1=2  ah Þαð0Þ C M ðtÞ ¼ πð1=2 þah Þ αð0Þ þ ξð0Þ _ ϕðtÞ þ ξ€  ah α€  α€ ð1=2  ah Þ α_ 2 16 2 Z π   € þð1=2  ah ÞαðσÞ _ þ ξðσÞ € πð1=2 þah Þ ϕðt  σÞ αðσÞ dσ;

ð2Þ

0

where μ is the airfoil–air mass ratio, and the Wagner function ϕðtÞ is given by Jone's approximation ϕðtÞ ¼ 1  ψ 1 e  ε1 t  ψ 2 e  ε2 t with the constants ψ 1 ¼0.165, ψ 2 ¼ 0.335, ε1 ¼ 0.0455, and ε2 ¼0.3. In order to eliminate the integral terms in the integro-differential Eq. (2), the following four new variables are introduced: Z t e  ε1 ðt  σÞ αðσÞdσ; w1 ¼ 0 Z t e  ε2 ðt  σÞ αðσÞdσ; w2 ¼ Z

0

Z

0

t

w3 ¼

t

w4 ¼ 0

e  ε1 ðt  σÞ ξðσÞdσ; e  ε2 ðt  σÞ ξðσÞdσ:

ð3Þ

Then, the plunging-and-pitching motion of the airfoil in Eq. (2) can be cast as follows: c0 ξ€ þc1 α€ þ c2 ξ_ þ c3 αþ _ c4 ξ þ c5 α þ c6 w1 þc7 w2 þ c8 w3 þ c9 w4 þ c10 GðξÞ ¼ f ðtÞ; € d2 ξ_ þ d3 αþ _ d4 ξ þ d5 αþ d6 w1 þ d7 w2 þd8 w3 þ d9 w4 þ d10 MðαÞ ¼ gðtÞ; d0 ξ€ þ d1 αþ _ 1 ¼ α ε1 w1 ; w _ 2 ¼ α ε2 w2 ; w _ 3 ¼ ξ  ε1 w3 ; w _ 4 ¼ ξ  ε2 w4 : w

ð4Þ

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

327

The coefficients ci and di which are the functions of system parameters are detailed in Appendix. f ðtÞ and gðtÞ are functions depending on the initial conditions, Wagner's function, and the forcing terms. When transients are damped out, and the steady solutions are obtained f ðtÞ ¼ gðtÞ ¼ 0. In the present study, structural nonlinearity is considered only in the pitching degree of freedom and GðξÞ ¼ ξ. Three kinds of structure nonlinear models are considered to describe the restore force M(α): freeplay, hysteresis and hybrid models. For these models, the functional curves of the moment M(α) versus displacement α are illustrated in Fig. 2.

2.1. Freeplay model The freeplay nonlinearity, as shown in Fig. 3(a), is determined by the following equation: 8 M þα  αf ; > < 0 MðαÞ ¼ M Fp ðαÞ ¼ M 0 þM f ðα  αf Þ; > : M þα  α þ δðM  1Þ; 0 f f

α oαf αf r αr αf þ δ ; α 4αf þ δ

where M0 means the preload; αf denotes the beginning of the freeplay; the freeplay δ and slope M f are constants.

Fig. 2. Schematic representation of the piecewise torque versus pitch angle. (a) Freeplay model, (b) hysteresis model, and (c) hybrid model.

ð5Þ

328

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

Fig. 3. The optimization results at U n =U nL ¼ 0.7 obtained by the proposed method. (a) Pupp and (b) Plow.

2.2. Hysteresis model The hysteresis function in Fig. 2(b) can be described by 8 _ 0 M 0 þ α αf ; α oαf ; α4 > > > > > _ 0 ; α r αr α þ δ; α4 M > 0 f f > > > < M 0 þα  αf  δ; _ 0 α4 α þ δ; α4 f ; MðαÞ ¼ M Hys ðαÞ ¼ þδ  M ; α oα þ δ; α _ o0 α þα > 0 f f > > > >  M0 ;  αf  δr αr αf ; α_ o0 > > > > : α4 αf ; α_ o0 α þ αf  M0 ;

ð6Þ

where the preload M0, the beginning of the freeplay αf , and the freeplay δ are constants.

2.3. Hybrid model The hybrid model in Fig. 2(c) introduced to the aeroelastic system is a hybrid expression of cubic and freeplay nonlinearities (Chen and Liu, 2014; Chen et al., 2014).  M Cub ðαÞ ¼ α þ γα3 8   αþ γα3 þ M 0 þα  αf ; > > <   αþ γα3 þ M 0 þM f ðα  αf Þ ; MðαÞ ¼ MHM ðαÞ ¼ M Cub ðαÞ þ M Fp ðαÞ ¼ >    > : αþ γα3 þ M 0 þα  αf þ δðM f  1Þ ;

α oαf αf r α rαf þ δ ;

ð7Þ

α 4αf þδ

where γ is the nonlinear coefficient of cubic term. In the following, three cases of the freeplay force M Fp ðαÞ in Eq. (7) are considered and Case 0 corresponding to purely cubic nonlinearity is used for comparison purpose: Case 0. MðαÞ ¼ M Cub ðαÞ

8 > < σα; MðαÞ ¼ M Cub ðαÞ þ l1 ðαÞ; l1 ðαÞ ¼ 0; > : σα;

Case 1.

ð8Þ

αo δ  δr α rδ : α4δ

ð9Þ

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

8 > < α þδ; MðαÞ ¼ M Cub ðαÞ þl2 ðαÞ; l2 ðαÞ ¼ 0; > : α δ;

Case 2.

Case 3. MðαÞ ¼ M

Cub

ðαÞ þl3 ðαÞ; l3 ðαÞ ¼

8 > <

329

α o δ  δ rα r δ:

ð10Þ

α 4δ

α þ δ;

αo  δ

M f ðα þδÞ;

 δr α rδ :

> : α þ δM ; f

ð11Þ

α 4δ

The freeplay nonlinear l1 ðαÞ is taken from Chen et al., (2014) and the coefficient σ is a representative of the level of nonlinearity of the dynamic system. The simulation case l2 ðαÞ is borrowed from Li et al. (2012) and l3 ðαÞ is simplified from Eq. (5). In a general case, time history of piecewise nonlinear forces is made of several parts. Each part corresponds to a cleardefined state. Therefore, the total period of the limit cycle can be broken into several time intervals. 3. The proposed method The piecewise constrained optimization harmonic balance method with analytical gradients is presented in this section. The analytical gradients of the objective function and the nonlinear equality constraints with respect to the optimization variables are derived. Based on the analytical gradients of the nonlinear equality constraints, the sensitivity of the Fourier coefficients with respect to influence parameters including the vibration frequency can also be derived by using the implicit function theory. 3.1. Analytical derivations of the nonlinear equality constraints 3.1.1. The evaluation of the nonlinear equality constraints In the following, the HBM is used to study the LCOs of the airfoil in Eq. (4). For the sake of brevity, only a brief discussion about harmonic balance approach is presented here; further details about this approach regarding the solution procedure can be found in Liu and Dowell (2004), Liu et al. (2007), and Liao and Sun (2013). In the HBM, αðtÞ is represented by a truncated Fourier series: αðtÞ ¼ TðtÞU; h where TðtÞ ¼ 1 cosðωtÞ

ð12Þ sinðωtÞ

⋯

cosðkωtÞ

sinðkωtÞ

⋯

cosðN H ωtÞ

sinðNH ωtÞ

i 1ð1 þ 2NH Þ

is a vector of dimension  U 0 U c1 U s1 ⋯ U ck

(1 þ2NH) where NH denotes the highest retained harmonics and k refers to the harmonic index. U ¼

U sk ⋯U cNH U sNH T with size (1 þ2NH) represents the Fourier coefficients. Upon introducing the new independent variable τ ¼ ωt, the dimensionless form of αðtÞ can be written as αðτÞ ¼ TðτÞU by replacing ωt with τ in TðtÞ. € Differentiating Eq. (12) with respect to t twice yields 3.1.1.1. Expressions of α_ and α. _ ¼ TðtÞ∇U; αðtÞ

€ ¼ TðtÞ∇2 U αðtÞ

ð13Þ

where  ∇ ¼ diag 0; ∇1 ; …; ∇k ; …; ∇NH ;

∇k ¼ ωk

 ∇2 ¼ ∇n∇ ¼ diag 0; Ξ1 ; …; Ξk ; …; ΞNH ;

0

1

1

0

 ð14Þ

Ξk ¼ ∇k n∇k ¼ ðωkÞ2

1

0

0

1

 ð15Þ

3.1.1.2. General formula of αk in the frequency domain. In the following, the general formula of αk is derived for the cubic term in Eq. (7). h iT c s c s c s Let βðtÞ ¼ TðtÞV with V ¼ V 0 V 1 V 1 ⋯ V k V k ⋯ V NH V NH . Using the following trigonometric identities cos ðiτÞ cos ðjτÞ ¼ 12 ½ cos ði  jÞτ þ cos ði þ jÞτ; sin ðiτÞcosðjτÞ ¼ 12 ½ sin ði þ jÞτ þ

sin ði jÞτ;

cos ðiτÞ sin ðjτÞ ¼ 12 ½ sin ði þjÞτ  sin ði  jÞτ; sin ðiτÞ sin ðjτÞ ¼ 12 ½ cos ði  jÞτ  cos ðiþ jÞτ:

ð16Þ

330

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

   The product of U ck cos ðiτÞ þU sk sin ðiτÞ V ck cos ðjτÞ þ V sk sin ðjτÞ is given as follows:  c    U k cos ðiτÞ þU sk sin ðiτÞ V ck cos ðjτÞ þ V sk sin ðjτÞ  c c  c c U V þ U sk V sk U V U sk V sk þ cos ði þ jÞτ k k ¼ cos ði jÞτ k k  s2 c  c 2s U k V k  U ck V sk U V þ U sk V ck þ sin ði þ jÞτ k k : þ sin ði jÞτ 2 2

ð17Þ

With the aid of Eq. (17), neglecting the higher order harmonics and equating the terms with the same harmonic index gives n o   αðtÞβðtÞ ¼ ½TðtÞUT ½TðtÞV ¼ UT ½TðtÞT ½TðtÞ V ¼ TðtÞEðUÞ V; ð18Þ where EðUÞ is called the operational 2 ϒ1 U0 6 H1 U0I þ Q 2 6 6 6 H2 L1 þ Q 3 6 6 ⋱ L2 þ Q 4 6 6 ⋱ Hk EðUÞ ¼ 6 6 6 ⋱ ⋱ 6 6 6 HNH  2 ⋱ 6 6H 4 N H  1 LNH  2 þ Q N H HNH with

" Hk ¼

U ck U sk

ϒ2

⋱

ϒk

⋱

ϒNH  2

ϒNH  1

N1 þ Q 3

N2 þQ 4

⋱

⋱

⋱

NNH  2 þ Q NH

U0I þ Q 4 L1 þ Q 5

N1 þQ 5 ⋱

⋱ ⋱

⋱ ⋱

Nk þ Q N H ⋱

⋱ Nk

⋱

⋱

U 0 Iþ Q NH

N1

N2

⋱

⋱

L2 þ Q NH

L1

U0 I

N1

N2

Lk þ Q NH

⋱

L2

L1

U0 I

N1

⋱

Lk

⋱

L2

L1

U0I

LNH  2

⋱

Lk

⋱

L2

L1

LNH  1

#

1h c ; ϒk ¼ U k 2

For instance, if NH ¼8, Eq. 2 ϒ1 U0 6H U Iþ Q2 0 6 1 6 6 H2 L1 þ Q 3 6 6H 6 3 L2 þ Q 4 6 EðUÞ ¼ 6 6 H4 L3 þ Q 5 6H 6 5 L4 þ Q 6 6 6 H6 L5 þ Q 7 6 6H 4 7 L6 þ Q 8 H8

matrix given as follows:

L7

" c i 1 Uk U sk ; Lk ¼ s 2 Uk

U sk U ck

#

" c 1 Uk ; Nk ¼ s 2  Uk

U sk U ck

#

" c 1 Uk ; Qk ¼ s 2 Uk

U sk  U ck

(19) reduces to ϒ2

ϒ3

ϒ4

ϒ5

ϒ6

ϒ7

N1 þQ 3

N2 þ Q 4

N3 þQ 5

N4 þ Q 6

N5 þ Q 7

N6 þ Q 8

U 0 I þQ 4 L1 þQ 5

N1 þ Q 5 U0I þ Q 6

N2 þQ 6 N1 þQ 7

N3 þ Q 7 N2 þ Q 8

N4 þ Q 8 N3

N5 N4

L2 þQ 6

L1 þ Q 7

U 0 I þQ 8

N1

N2

N3

L3 þQ 7

L2 þ Q 8

L1

U0 I

N1

N2

L4 þQ 8

L3

L2

L1

U0I

N1

L5

L4

L3

L2

L1

U0I

L6

L5

L4

L3

L2

L1

ϒ8

ϒNH

3

NNH  1 7 7 7 NNH  2 7 7 ⋱ 7 7 7 N k 7; 7 ⋱ 7 7 7 N2 7 7 N1 7 5 U0 I

ð19Þ

# :

3

N7 7 7 7 N6 7 7 N5 7 7 7 N4 7 7: N3 7 7 7 N2 7 7 N1 7 5 U0I

ð20Þ

Thus, the analytical expansion of α2 is n o α2 ¼ αα ¼ ½TðtÞUT ½TðtÞU ¼ UT ½TðtÞT TðtÞ U ¼ TðtÞEðUÞU

ð21Þ

The cubic term α3 is obtained by repeated application of EðUÞ: o   n α3 ¼ αα2 ¼ ½TðtÞUT TðtÞEðUÞU ¼ UT ½TðtÞT TðtÞ ½EðUÞU   ¼ TðtÞEðUÞ ½EðUÞU ¼ TðtÞ½EðUÞ2 U:

ð22Þ

Therefore, the polynomial nonlinear terms satisfy the recurrence relation n o αk ¼ ½TðtÞUT αk  1 ¼ TðtÞ ½EðUÞk  1 U :

ð23Þ

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

3.1.1.3. Expression for various piecewise nonlinearities. h i 8 > ; ðUÞ; SHys ðUÞ ¼ U þSHys_0 TðtÞSHys > 1 1 1 > > > > Hys Hys Hys_0 > > ; TðtÞS2 ðUÞ; S2 ðUÞ ¼ S2 > > h i > > > Hys Hys Hys_0 > ; ðUÞ; S ðUÞ ¼ U þS TðtÞS < 3 3 3 h i M Hys ðαÞ ¼ Hys Hys_0 > > ; TðtÞSHys > 4 ðUÞ; S4 ðUÞ ¼ U þS4 > > > > Hys Hys Hys_0 > > ; TðtÞS5 ðUÞ; S5 ðUÞ ¼ S5 > > h i > > Hys Hys_0 > : TðtÞSHys ; 6 ðUÞ; S6 ðUÞ ¼ U þS6

The Fourier expression of the hysteresis model M Hys ðαÞ is αo αf ; α_ 4 0 αf rα r αf þ δ; α_ 4 0 α4 αf þ δ; α4 _ 0 αo αf þ δ; αo _ 0

;

_ 0 α4 αf ; αo

6

Using Eq. (22), the Fourier expression of M Cub ðαÞ is given as follows:    M Cub ðαÞ ¼ α þ γα3 ¼ TðtÞSCub ðUÞ; SCub ðUÞ ¼ U þ γ½EðUÞ2 U : Then, SCub ðUÞ can be used to derive the Fourier expression of the hybrid model in Eq. (7) as follows: 8 n o Fp HM HM Cub > α oαf > > TðtÞS1 ðUÞ; S1 ðUÞ ¼ S ðUÞ þ S1 ðUÞ ; > > n o < Fp HM HM Cub M HM ðαÞ ¼ TðtÞS2 ðUÞ; S2 ðUÞ ¼ S ðUÞ þ S2 ðUÞ ; αf r αr αf þ δ ; > n o > > > Fp HM HM Cub > α 4 αf þ δ : TðtÞS3 ðUÞ; S3 ðUÞ ¼ S ðUÞ þ S3 ðUÞ ; where SFp j ðUÞ; j ¼ 1; 2; 3 are the Fourier coefficients related to the three states in Eq. (5): h iT h i Fp_0 ; SFp_0 ¼ M0  αf 0 0 ⋯ 0 0 ⋯ 0 0 ; SFp 1 ðUÞ ¼ U þS1 1 h iT h i Fp_0 SFp ¼ M 0 M f αf 0 0 ⋯ 0 0 ⋯ 0 0 ; ; SFp_0 2 ðUÞ ¼ M f Uþ S2 2 h iT h i SFp ðUÞ ¼ U þSFp_0 ; SFp_0 ¼ M0  αf þ δðMf  1Þ 0 0 ⋯ 0 0 ⋯ 0 0 ; 3

ð24Þ

 αf  δr α r αf ; α_ o0

where the elements of vectors SHys_0 ; j ¼ 1; 2; …; 6 are zeros except for the first element. j h iT Hys_0 ¼ M 0  αf 0 0 ⋯ 0 0 ⋯ 0 0 ; S1  T SHys_0 ¼ M0 0 0 ⋯ 0 0 ⋯ 0 0 ; 2 h iT SHys_0 ¼ M 0  αf  δ 0 0 ⋯ 0 0 ⋯ 0 0 ; 3 h iT SHys_0 ¼ αf þ δ M 0 0 0 ⋯ 0 0 ⋯ 0 0 ; 4  T SHys_0 ¼ M 0 0 0 ⋯ 0 0 ⋯ 0 0 ; 5 h iT SHys_0 ¼ αf  M 0 0 0 ⋯ 0 0 ⋯ 0 0 :

3

331

3

ð25Þ

ð26Þ

ð27Þ

ð28Þ

and except for the first element all the others of SFp_0 ; j ¼ 1; 2; 3 are 0. j 3.1.1.4. Expression of the nonlinear equality constraints. Substituting Eqs. (24), (27) and (28) into Eq. (4) and applying the Galerkin method yield Z 2π   TT ðτÞTðτÞdτ A1 U þB1 Uξ ¼ 0; 0 Z 2π Z 2π   TT ðτÞTðτÞdτ A2 U þB2 Uξ þ d10 TT ðτÞMðαÞdτ ¼ 0; ð29Þ 0

where Uξ ¼

0

B1 1 A1 U

denotes the Fourier coefficients related to plunge dof ξ and

A1 ¼ c1 ∇ þ c3 ∇ þc5 I þ c6 ð∇ þ ε1 IÞ  1 þc7 ð∇ þ ε2 IÞ  1 ; 2

B1 ¼ c0 ∇2 þc2 ∇ þ ðc4 þc10 ÞI þ c8 ð∇ þ ε1 IÞ  1 þ c9 ð∇ þ ε2 IÞ  1 ; A2 ¼ d1 ∇2 þ d3 ∇ þ d5 Iþ d6 ð∇ þ ε1 IÞ  1 þd7 ð∇ þ ε2 IÞ  1 ; B2 ¼ d0 ∇2 þd2 ∇ þ d4 I þ d8 ð∇ þ ε1 IÞ  1 þ d9 ð∇ þ ε2 IÞ  1 :

ð30Þ

The explicit expressions for the linear elements of the above matrices A1 , B1 , A2 and B2 are worked out the same as in Dai R 2π et al. (2014). In the following, efforts are made to deal with the nonlinearity term 0 TT ðτÞMðαÞdτ.

332

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

During one period of oscillation, the piecewise nonlinear function MðαÞ has several states. According to the condition at each state, the period [0,2π] is divided into a set of N subintervals [τi ; τi þ 1 ] where τi and τi þ 1 denote the time instants at which the non-smooth nonlinearities take place. Then, the Fourier representations Fi ðUÞ are obtained by the definition of each state. Finally, Fourier expansion of the piecewise nonlinear force can be obtained as a sum of contributions given by each part of the piecewise nonlinear force. MðαÞ ¼

Fi ðUÞ A

8 > > > > > < > > > > > :

N 1 X

T½τi ; τi þ 1 Fi ðUÞ;

i¼0

SFp j ðUÞ; j ¼ 1; 2; 3

for the freeplay model in Eq: ð5Þ

SHys j ðUÞ; j ¼ 1; 2; …; 6 SHM j ðUÞ; j ¼

for the hysteresis model in Eq: ð6Þ ;

1; 2; 3

ð31Þ

for the hybrid model in Eq: ð7Þ

where the harmonic functions between τi and τi þ 1 are represented by T½τi ; τi þ 1 . It should be noted that the nonlinear force transition conditions have a nonlinear dependence on the unknowns U. Transition conditions related to the continuity of physical displacement and force at the instances when the system transitions between different states should be enforced. Using the expression in Eq. (31), the nonlinear term can be converted into the following state-space matrix form: ( ) Z 2π Z 2π N 1 X T T T ðτÞMðαÞdτ ¼ T ðτÞ T½τi ; τi þ 1 Fi ðUÞ dτ 0

0

¼

N  1 Z τi þ 1 X τi

i¼0

i¼0



NX 1 n o ½τ ;τ  TT ½τi ; τi þ 1 T½τi ; τi þ 1 dτ Fi ðUÞ ¼ WMi i þ 1 Fi ðUÞ ;

ð32Þ

i¼0 ½τ ;τi þ 1 

where the elements of the matrix WMi Z τi þ 1 ½τ ;τ  TT ðτÞTðτÞdτ; WMi i þ 1 ¼

are given as follows:

τi

½τ ;τ  WMi i þ 1 ð1; 1Þ ¼ τ=2;

½τ ;τi þ 1 

WMi

ð1; 2iÞ ¼

sin ðiτÞ ; 2i

½τ ;τi þ 1 

WMi

ð1; 2iþ 1Þ ¼ 

cos ðiτÞ ; 2i

sin ðiτÞ cos ðiτÞ ½τ ;τ  ; WMi i þ 1 ð2i þ1; 1Þ ¼  ; i i τ sin ð2iτÞ 1 þ cos ð2iτÞ ½τ ;τ  ½τ ;τ  WMi i þ 1 ð2i; 2iÞ ¼ þ ; WMi i þ 1 ð2i; 2i þ 1Þ ¼  ; 2 4i 4i 1 þ cos ð2iτÞ τ sin ð2iτÞ ½τ ;τ  ½τ ;τ  WMi i þ 1 ð2iþ 1; 2iÞ ¼  ; WMi i þ 1 ð2i þ 1; 2iþ 1Þ ¼  ; 4i 2 4i sin ðði jÞτÞ sin ðði þ jÞτÞ cos ðði  jÞτÞ cos ðði þjÞτÞ ½τ ;τ  ½τ ;τ  WMi i þ 1 ð2i; 2jÞ ¼ þ ; WMi i þ 1 ð2i; 2j þ 1Þ ¼  ; 2ði  jÞ 2ðiþ jÞ 2ði jÞ 2ði þ jÞ cos ðði  jÞτÞ cos ðði þjÞτÞ sin ðði  jÞτÞ sin ðði þ jÞτÞ ½τ ;τ  ½τ ;τ  WMi i þ 1 ð2iþ 1; 2jÞ ¼  þ ; WMi i þ 1 ð2iþ 1; 2j þ 1Þ ¼  : 2ði jÞ 2ði þ jÞ 2ði  jÞ 2ðiþ jÞ ½τ ;τi þ 1 

WMi

ð2i; 1Þ ¼

ð33Þ

Substituting Eq. (32) into Eq. (29) yields N 1 n n o o X ½τ ;τ  WMi i þ 1 Fi ðUÞ ¼ 0: CE ðUÞ ¼ diagð2π; π; …; πÞ ðA2  B2 B1 1 A1 ÞU þ d10

ð34Þ

i¼0

Eq. (34) is a nonlinear equation that represents the form of governing equations in the harmonic balance approach. Due to the nonlinear nature of Eq. (34), the nonlinear algebraic equations in Eq. (34) are used to construct the nonlinear equality constraints of the optimization problem. 3.1.2. Gradients of the nonlinear equality constraints In the following, the gradients of the nonlinear equality constraints are estimated analytically through the direct differentiation of Eq. (34) with respect to each influence parameter. According to Eqs. (19) and (34), analytical formulation of the Jacobian matrix is possible since derivatives of Fourier expansion coefficients of nonlinear forces can be performed analytically. Based on Eq. (34), the gradient of CE ðUÞ with respect to U is calculated as follows:

N 1 n o X ∂CE ðUÞ ½τ ;τ  ∂F ðUÞ ¼ diagð2π; π; …; πÞ ðA2  B2 B1 1 A1 ÞU þd10 : ð35Þ WMi i þ 1 i J¼ ∂U ∂U i¼0

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346 i ðUÞ Based on Eq. (31), ∂F∂U in Eq. (35) is calculated as follows: 8 Fp ∂Sj ðUÞ > > for the freeplay model in Eq: ð5Þ > ∂U ; j ¼ 1; 2; 3 > < Hys ∂Fi ðUÞ ∂Sj ðUÞ : A ∂U ; j ¼ 1; 2; …; 6 for the hysteresis model in Eq: ð6Þ > ∂U > > HM > ∂S ðUÞ : j ; j ¼ 1; 2; 3 for the hybrid model in Eq: ð7Þ

333

ð36Þ

∂U

By virtue of Eqs. (24) and (28), the calculations of evaluation of

∂SHys ðUÞ j ∂U

∂SHM ðUÞ j ∂U

and

∂SFp ðUÞ j ∂U

are generally straight forward. On the contrary, the

∂SCub ðUÞ ∂U

is transformed into the calculation of according to Eq. (27). ∂f½EðUÞ2 U g ∂fU þ γ½EðUÞ2 U g ∂ EðUÞU ∂SCub ðUÞ in ∂U ¼ , the analytical form of the partial derivative f ∂U g is derived as In order to compute ∂U ∂U follows:   ∂ EðUÞU ∂EðUÞ U þ EðUÞ: ð37Þ ¼ ∂U ∂U Using the expression in Eq. (19), each column of the matrix ∂EðUÞ ∂U U can be derived rather than calculating the individual elements: " # ∂EðUÞ ∂EðUÞ ∂EðUÞ ∂EðUÞ ∂EðUÞ ∂EðUÞ ∂EðUÞ ∂EðUÞ U¼ U; U; U; …; U; U; …; U; U ; ð38Þ ∂U ∂U 0 ∂U ck ∂U sk ∂U ck ∂U sk ∂U cNH ∂U sNH where 2

0 6 6 6 6 6 6 6 6 6 6 ∂EðUÞ 6 6 ∂Hk ¼ 6 Χ 6 ∂U k ∂U Χ k 6 6 6 6 6 6 6 6 6 4 0 with the elements

 1 ∂Hk ; c ¼ ∂U k 0

21 0 ∂Hk ¼ ; ∂U sk 1 21

∂ϒk ∂U Χ k

0

∂Q k ∂U Χ k

7 7 7 7 7 7 7 7 7 7 7 ∂Nk 7 7 Χ 7 ∂U k 7 7 7 7 7 7 7 7 7 5 0

∂Nk ∂U Χ k ∂Nk ∂U Χ k

⋰ ∂Q k ∂U Χ k

⋱

∂Lk ∂U Χ k ∂Lk ∂U Χ k

⋱ ∂Lk ∂U Χ k





 0 ∂Nk 1 1 0 ∂Q k 1 1 ; ; c ¼ c ¼ 1 ∂U k 2 0 1 ∂U k 2 0  1 



  1 ∂Nk 1 0 1 ∂Q k 1 0 1 ; ; : s ¼ s ¼ 0 ∂U k 2  1 0 ∂U k 2 1 0   k ∂ ½EðUÞ U for k42: Therefore, the following recurrence relation allows to calculate ∂U n o n h io h i k  1 k i ∂ ½EðUÞk U ∂ EðUÞ ðEðUÞÞ U ∂ ðEðUÞÞ  1 U ∂EðUÞh k1 ðEðUÞÞ ¼ ¼ : U þ EðUÞ ∂U ∂U ∂ðUÞ ∂U ∂ϒk 1 ¼ 1 ∂U ck 2 ∂ϒk 1 ¼ 0 ∂U sk 2

∂Lk 1 1 c ¼ ∂U k 2 0

 ∂Lk 1 0 1 12 ; s ¼ ∂U k 2 1

3

0



12

;

ð39Þ

0

ð40Þ

Following an analogous methodology, the gradients of CE ðUÞ with respect to other influence parameters can be obtained. 3.1.3. General scheme for calculation of vibration response sensitivity The sensitivity of the Fourier coefficients with respect to the influence parameters can easily be determined by analytical gradients. Using the implicit function theory (or the chain rule of differentiation), the total derivative of CE ðUÞ with respect to the kth structural parameter bk are given as dCE ðUÞ ∂CE ðUÞ ∂CE ðUÞ ∂U ∂CE ðUÞ ∂U ¼ þ ¼ þJ ¼ 0: dbk ∂bk ∂U ∂bk ∂bk ∂bk

ð41Þ

Taking the inverse of the J matrix in Eq. (41), the sensitivity can be expressed as a first order explicit form: ∂U ∂CE ðUÞ ¼ J1 : ∂bk ∂bk

ð42Þ

To summarize, using the developed analytical formulation, the gradients of the nonlinear equality constraints can be computed. A system sensitivity analysis based on the developed analytical gradients has been presented.

334

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

3.2. Stability analysis for the piecewise nonlinear systems Stability is a central issue in dynamical system. After the limit cycle is constructed, the stability analysis of such limit cycle should be conducted. According to the Floquet theory (Brown et al., 2012; Bittanti and Colaneri, 2009; Wang and Hale, 2001), the stability analysis can be performed by studying the complex eigenvalues (Floquet multipliers) of the monodromy matrix associated with the limit cycle of Eq. (4). In order to compute the monodromy matrix, Eq. (4) can be rewritten in the state space form as follows: z_ ðtÞ ¼ hðt; zðtÞÞ;   _ w1 ; w2 ; w3 ; w4 T denotes the state space vector. where z ¼ α; α; _ ξ; ξ; Dividing both sides of Eq. (43) by z0 which is the initial condition of zðtÞ and some rearrangements yield: 

 d ∂zðtÞ ∂hðt; zðtÞÞ ∂zðtÞ ¼  ∂z ; dt ∂z0 ∂z 0 zðtÞ

ð43Þ

ð44Þ

with ∂zð0Þ ∂z0 ¼ I where I represents the identity matrix. With the application of Floquet theory, Eq. (44) is integrated over one period of the LCOs and the eigenvalues (related to ρ ¼ ½ρ1 ; ρ2 ; …; ρ8 ) of the matrix ∂zðtÞ ∂z0 at t ¼ 2π=ω are used to determine the stability of LCOs. When the absolute values of all the eigenvalues are smaller than 1, the system is stable. In contrast to this condition, if at least one of the absolute values of ρ is larger than 1, the system is unstable. Therefore, the necessary and sufficient conditions for solution stability require jρn jo 1 for all n and the following stability criteria can be inferred: g S ðU; ωÞ ¼ maxðjρjÞ 1 o 0:

ð45Þ

3.3. Analytical gradients of the objective function In this section, the gradients of the maximum displacement are given. _ max Þ ¼ 0. Since the objective Assume that αðτÞ attains its maximum over the interval [0,2π] at the point τmax where αðτ function has a known expression, the gradient of the objective function with respect to U can be obtained as  ∂½TðτÞU τ ¼ τmax ∂αðτmax Þ ¼ ¼ TðτÞ .Obviously, the gradients of αðτmax Þ with respect to other influence parameters are zeros. ∂U

∂U

τ ¼ τmax

3.4. Optimization problem formulation Within the framework of the previous proposed method(Liao, 2013; Liao and Sun, 2013), the nonlinear algebraic equations in Eq. (34) and the stability criterion in Eq. (45) are treated as the generally nonlinear constraints. Therefore, the following nonlinear optimization problem can be formulated: f ðxÞ ¼ f ðU; ωÞ ¼ max αðτÞ ( gðxÞ ¼ CE ðUÞ ¼ 0   s:t ; g ðxÞ ¼ maxðρÞ  1 r 0 s

ð46Þ

where x ¼ fU; ωgT . gðxÞ and g s ðxÞ represent the nonlinear equality and inequality constraints, respectively. Once the Fourier coefficients and the vibration frequency are updated during each optimization iteration step, the displacement and velocity are calculated using Eqs. (12) and (13) respectively. Based on the state transition conditions in Eqs. (5), (6), (8)–(11), τi and τi þ 1 can be obtained and the period [0, 2*π] can be divided into a set of N subintervals according to the states defined in Eqs. (5), (6), (8)–(11). Finally, the harmonic balance equations are constructed and solved. In this way, the nonlinear function can be obtained accurately compared with the methods in Li et al. (2012) and Liu et al. (2012). 4. Numerical results In this section, numerical studies of an airfoil system with freeplay, hysteresis and hybrid piecewise nonlinearities are carried out to demonstrate the effectiveness of the proposed method. The results are compared to direct numerical simulation results. To align with the computational studies in Chung et al. (2007), Liu et al. (2012), and Chen et al. (2014), the following structural parameters were chosen: μ¼100, rα ¼0.5, ah ¼ 0.5, ϖ ¼ 0.2, xα ¼0.25, and ςα ¼ ςξ ¼ 0. 4.1. Numerical results for the freeplay model The freeplay model used is shown in Eq. (5) where the values of four parameters are M 0 ¼0, M f ¼0, δ¼0.5 and af ¼0.25. The system parameters are chosen from Chung et al. (2007) and a critical flutter point for linear system is found at U nL ¼6.2851. Similar to the value of NH in Liu et al. (2012), the highest harmonic retained NH is 10. In the following, the proposed method is used to find the upper and lower response bounds at a given flow velocity U n =U nL ¼ 0:7. While looking

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

335

for multiple solutions at a given flow velocity, the unknown variables that have to be determined are the unknown Fourier coefficients U and the vibration frequency ω of the periodic solution. It should be noted that the objective function in Eq. (46) is changed as the minimization of the vibration amplitude for searching the lower bound solution. After the optimization algorithm stops successfully, numerical results are obtained. The optimization results are reported in Fig. 3. It can be seen from Fig. 3 that the presence of several harmonic components can be detected in the frequency spectrums and there is one main harmonic component in the system response. Observing the Fourier coefficients in Fig. 3, the constant harmonic terms are identified as the most prominent after the fundamental. Compared with solution Plow, the solution Pupp shows high magnitude oscillation not only at the odd harmonic components but it also exhibits harmonic response at the even harmonic components. It should be noted that for calculating the solutions Pupp and Plow, the proposed method can be viewed as a root finding algorithm. In order to fully validate the proposed approach, direct numerical integrations have been performed. Numerical simulations are computed by the fourth-order Runge–Kutta scheme with fixed step-size of 0.001 here. The initial conditions can be readily supplied by the results of the presented method. In order to eliminate the transient part of the responses, the initial responses are discarded from the stored responses and results are plotted only for simulation time between 1000 and 1050 s. The time series of the airfoil system at these two solutions using the time integration method are presented in Fig. 4 while the comparison of the phase portraits between the two methods is depicted in Fig. 5 where the notations ‘RK’ and ‘PCOHBM’ represent the solutions obtained by the time integration method and the proposed method, respectively.

Fig. 4. Time series of the LCOs at U n =U nL ¼ 0.7 using the time integration method. (a) Pupp and (b) Plow.

Fig. 5. The comparison of the phase portraits of the two LCOs between the proposed method and time integration method. (a) Pupp and (b) Plow.

336

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

As shown in Fig. 5, the results from the proposed method and from the time domain integration method are difficult to distinguish when they are plotted on the same graph. The limit cycles from both methods agree satisfactorily with each other. Moreover, for solution Pupp, a maximum value of 1.5183 was produced by the airfoil system at frequency of 0.076632 while a maximum value of 1.2949 was experienced by the solution Plow at the frequency of 0.087209. It is obvious that the frequency corresponding to solution Pupp is lower than that of Plow. Figs. 6 and 7 show the time histories of nonlinear forces and the corresponding nonlinear force–displacement relationship at these two solutions. As illustrated in Figs. 6 and 7, the results clearly indicate that the system is acting with the freeplay model in Fig. 2(a) and the piecewise vibration movements of the three stages are identified. Using the continuation algorithm, LCO solutions can be tracked by the proposed method starting from U n =U nL ¼0.7 backward or forward. The effects of U n on the vibration amplitudes and frequencies of LCOs are depicted in Fig. 8 where the solid line and the dashed line denote the stable and unstable solutions, respectively. It can be seen from Fig. 8 that the response amplitude increases with the increase in U n and approaches to infinity at the upper limit of U n =U nL whereas the vibration frequency decreases. The results shown in Fig. 8 demonstrate that two independent branches of the amplitude and frequency curves are found and the two branches do not have any connection with each other. Two LCOs may exist simultaneously for certain values of U n =U nL . To show the curves and branches in detail, the magnification of plots is also presented in Fig. 8. Interestingly, the two LCOs are analogous to each other in shape and the amplitude of the upper branch is slightly greater than the counterpart of the lower branch. The frequencies of the LCOs in the lower branch are greater than that of the upper one. The stabilities of the LCOs obtained by the proposed method are judged by the Floquet theory in Section 3.2. Stability analysis in Fig. 8 shows that dynamic instability occurs for both the curves and the solutions located at the middle segments of the two branches are unstable. As can be seen from Fig. 8, the results are in agreement with those of Chung et al. (2007).

Fig. 6. The nonlinear force and force–displacement relationship curve for Pupp.

Fig. 7. The nonlinear force and force–displacement relationship curve for Plow.

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

Fig. 8. The LCO amplitude and frequency curves for the freeplay model.

Fig. 9. The sensitivity results of the two LCOs at U n =U nL ¼0.7. (a) Pupp and (b) Plow.

Fig. 10. The sensitivity results related to the curves in Fig. 8. (a) The lower branch and (b) the upper branch.

337

338

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

For the solutions in Fig. 3, the sensitivity coefficients of the harmonic contents for LCOs to variation of the structural parameters are calculated and plotted in Fig. 9. As shown in Fig. 9, the dominant contribution of the sensitivity is the first and second harmonic contents with negligible contributions from other higher terms. It can be noted that the solutions Pupp and Plow have great difference between the sensitivity. The solution Pupp has large sensitivity values while the sensitivities for solution Plow are relatively small. For all LCOs studied the Fourier coefficients appear to be very sensitive to the variation of frequency ω. The sensitivity results show that the sensitivities of U with respect to the structural parameters U* and δ have the lowest sensitivity. For the curves in Fig. 8, sensitivities of the Fourier coefficients with respect to the influence parameters are displayed in Fig. 10. Looking at the graphics in Fig. 10(a) which corresponds to the sensitivity results of the lower branch in Fig. 8, it is observed that the sensitivity values of U with respect to U n and δ increase monotonically with increasing U n . Conversely, the sensitivity coefficients for ∂U=∂U n increase sharply at flow velocities which are characterized by a change of stability (the bifurcation points). For ∂U=∂ω, the sensitivity magnitude of ∂U=∂ω decreases to a minimum value at U n =U nL ¼0.72 and then begins to increase. In Fig. 10(b), sensitivity coefficients with respect to U n and M0 are higher than those with respect to ω and δ. 4.2. Numerical results for the hysteresis model After the success of the proposed method in the above freeplay model, the proposed method is now applied to a hysteresis model. A hysteresis nonlinearity which is modeled by Eq. (6) is introduced to the governing equations in Eq. (4). Numerical examples which have been taken from recent publication (Liu et al., 2012) are considered with the following structural parameters: M0 ¼ 0:5π=180, δ ¼ π=180 and af ¼0.

Fig. 11. The optimization results at U n =U nL ¼ 0.9 for the hysteresis model.

Fig. 12. The time history and phase portraits for the LCO at U n =U nL ¼ 0.9.

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

Fig. 13. The nonlinear force and hysteresis curve for the LCO at U n =U nL ¼ 0.9.

Fig. 14. The LCO amplitude and frequency curves for the hysteresis model.

Fig. 15. Sensitivity results for the hysteresis model. (a) Sensitivity of LCO at U n =U nL ¼ 0.9 and (b) sensitivity results for Fig. 14.

339

340

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

Fig. 16. The optimization results obtained by the proposed method.

Fig. 17. Time histories of the LCOs at U n =U nL ¼ 1.5 for the four cases. (a) Case 0, (b) Case 1, (c) Case 2 and (d) Case 3.

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

341

Fig. 18. The comparison of the phase portraits of the LCO at U n =U nL ¼1.5 for the four cases between the two approaches. (a) Case 0, (b)Case 1, (c)Case 2 and (d)Case 3.

To perform a validation example, the LCO at U n =U nL ¼0.9 is first considered. The optimization results at U n =U nL ¼0.9 are shown in Fig. 11. As shown in Fig. 11, the frequency spectrum is primarily composed of the first harmonic term, the influence of other harmonic components being insignificant. Resembling the previous numerical analysis in Section 4.1, the time history obtained by the time integration method and the phase portrait comparison between the proposed method and the time integration method are plotted in Fig. 12. The comparison shows the trajectories on the phase plane between the two methods coinciding with each other. When compared with the results in Fig. 6 of Liu et al. (2012), the results in Fig. 11 are consistent with those of Liu et al. (2012). The agreement of comparison between the two approaches and the published results given in Fig. 6 of Liu et al. (2012) indicate that the present analysis is accurate. The nonlinear force and hysteresis curve associated with the optimized solution in Fig. 11 are illustrated in Fig. 13. As illustrated in Fig. 13, it can be seen that the periodic solution is made up of five distinct segments and classical hysteresis curve of the nonlinear force–displacement relationship can be observed. In a single period, a nonlinear hysteretic loop forms between the nonlinear force and the vibration displacement, which demonstrates its nonlinear characteristics. By taking U n as the continuation parameter, the influence of U n on the vibration amplitude and frequency of LCO is depicted in Fig. 14. In Fig. 14 the vibration amplitude is found to increase as the flow velocity is increased and reaches a maximum value at the upper limit of U n =U nL . Furthermore, the vibration frequency increases as the flow velocity increases. In addition, when the flow velocity is smaller than 0.8174, the airfoil system shows unstable periodic motions. Consequently, they are shown by the dashed lines.

342

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

Fig. 19. Time history of nonlinear restoring force and force–displacement relationship curve for Case 0 when U n =U nL ¼1.5.

In Fig. 15, the sensitivities of the Fourier coefficients with respect to the structural parameters are depicted. It is observed in Fig. 15(a) that the two parameters M0 and δ have a significant impact on the sensitivity results. The first order sensitivity results associated with the curves in Fig. 14 are shown in Fig. 15(b). It can be seen that, for all the cases considered, some sudden drops in these sensitivities are noted and the harmonic coefficients are the most sensitive to variation of all parameters at 0.8491. The sensitivity results show that the sensitivities of U is the lowest with respect to ω and sensitivity coefficients with respect to M0 are greater than those with respect to δ and U n . For the hysteresis model in Liu et al. (2012), the LCOs are followed by taking M0 or δ as the continuation parameter. In Figs. 13 and 14 of Liu et al. (2012), it is seen that the LCO amplitudes and frequencies are descending as the preload M0 increases up to a certain critical level. Moreover, studies in Figs. 16 and 17 of Liu et al. (2012) show that the LCO amplitudes and frequencies increase with the freeplay δ. However, these vibration behaviors in Figs. 13, 14, 16 and 17 of Liu et al. (2012) cannot be predicted by the continuation method since the continuity conditions are violated. During the vibration cycle, the system undergoes transitions between different states. At the moment of transition, compatibility conditions in terms of displacement and nonlinear restoring force should be satisfied. As one of system parameters M0 or δ varies while all the other parameters are kept constant, the continuity of piecewise nonlinear force at the transition boundaries in Eq. (6) is not satisfied. 4.3. Numerical results for the hybrid model To gain more insight into the combined effects of the cubic and freeplay nonlinearities on the system behavior, simulations are carried out for the hybrid model given in Eqs. (8)–(11). For the structural coefficients in these equations, the structural parameter values used are ϖ ¼0.25, σ ¼ 0:08; Mf ¼ 0:08; δ ¼ 0:01 and these parameters are taken from Chen et al. (2014). The flutter velocity of corresponding linear system is U nL ¼ 6.0385. Using the proposed method, the numerical results of the LCO at U n =U nL ¼1.5 for the four cases are given in Fig. 16. As can be seen from Fig. 16, the LCOs of the airfoil system are dominated by the first harmonic. The numerical results of the developed method are summarized in Table 1 where the Floquet multipliers ρj of the stability analysis are also listed. For all attained LCOs, there is always one multiplier 1, which complies with the Floquet theory. A pair of complex conjugate Floquet multipliers exists and five of these Floquet multipliers are purely real parts. Except for the solution corresponding to Case 1, all Floquet multipliers stay inside the unit circle and therefore all LCOs are stable. In order to verify the theoretical predictions, time numerical simulations are conducted using the fourth-order Runge– Kutta algorithm. The time histories of the LCOs for the four cases are shown in Fig. 17. For comparison purpose, the portraits predicted by the proposed method and time integration method are reported in Fig. 18. It is worth noting that these phase orbits are commonly of regular shape, but their amplitudes are different. For Cases 2 and 3 the vibration amplitudes are reduced compared to the results related to Cases 0 and 1. The computed displacement amplitudes for Cases 0 and 1 are about twice higher than that for Cases 2 and 3. However, it can be seen from Fig. 18 that the frequencies of LCOs change slightly, no matter which models are considered. Fig. 18 shows that the simulation results from the proposed method match well with numerical integration results. Therefore, the validity of the proposed piecewise harmonic balance method is verified. Using the force–displacement function given in Eqs. (8)–(11), the nonlinear forces as a function of non-dimensional time and the force–displacement curves associated with the solutions in Fig. 16 are depicted in Figs. 19–22 for the four cases. Fig. 19 shows a typical time history of the cubic nonlinear force and the amplitude has a notably larger amplitude. As

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

343

Fig. 20. Time history of nonlinear restoring force and force–displacement relationship curve for Case 1 when U n =U nL ¼ 1.5. (a)The total nonlinear force, (b) Zoom-in part showing the freeplay piecewise force.

Fig. 21. Time history of nonlinear restoring force and force–displacement relationship curve for Case 2 when U n =U nL ¼ 1.5

344

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

illustrated in Figs. 20–22, the total behavior of the nonlinear force can be obtained with the superposition of cubic nonlinear force and piecewise freeplay force. The time responses of the nonlinear forces contain three distinctive states and three stages are predicted by all the models. The time history of nonlinear force and force–displacement relationship curve for Case 1 (Fig. 20) is characterized by inconsistency transitions of different states. Fig. 20(b) shows zoomed-in sections of the freeplay signals presented in Fig. 20(a), where the discontinuity transition occurs. It is indicated that mistakes occur if the model in Eq. (9) is employed to obtain the LCOs. For the solutions corresponding to Cases 2 and 3, the nonlinear piecewise forces l2 ðαÞ and l3 ðαÞ computed by means of Eqs. (10) and(11) are continuous and the continuity conditions of state transition are met. It is observed in Figs. 21 and 22 that the motion patterns of the airfoil system for l2 ðαÞ and l3 ðαÞ are slightly different from each other. Moreover, the cubic nonlinear components are mostly dominant and their amplitudes are higher than those of the piecewise components. Sensitivity analysis is employed to determine how variations of independent design variables impact dependent responses. By using Eq. (42), sensitivity results for the solutions in Fig. 16 are plotted and shown in Fig. 23 for ∂U=∂ω and in Fig. 24 for ∂U=∂U n . In Figs. 23 and 24, Fig. (23) show enlargements for Cases 1–3 in the preceding Fig. (23). Note that the numerical results are the amplitudes of the first and third harmonics and the maximum amplitudes of the higher harmonics are far less than that for the first harmonic. From Figs. 23 and 24 and Table 1 it can be noted that the sensitivities for the cubic and freeplay nonlinearities hybrid systems(Cases 1–3) are relatively lower compared to the counterparts of the purely cubic nonlinearity system(Case 0) and the sensitivities of the cubic nonlinear systems are greatly suppressed in the presence of freeplay nonlinearity.

Fig. 22. Time history of nonlinear restoring force and force–displacement relationship curve for Case 3 when U n =U nL ¼1.5.

Fig. 23. Sensitivity results for ∂U=∂ω. (a) The results for Cases 0–3 and (b) zoom-in part showing Cases 1–3.

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

345

Fig. 24. Sensitivity results of the Fourier coefficients to variations in the flow velocity. (a) The results for Cases 0–3 and (b) zoom-in part showing Cases 1–3.

5. Conclusions The piecewise constrained optimization harmonic balance method has been developed to obtain the accurate periodic solutions of piecewise nonlinear systems. Piecewise expressions of the nonlinear forces in the frequency domain are presented. Exact derivatives of the Fourier coefficients of the nonlinear forces with respect to the counterparts of structural displacement are also derived analytically. The analytical derivation allows accurate and fast evaluation of the sensitivity coefficients. The sensitivity characteristics are calculated in frequency domain over the flow velocity range. The methodology developed is illustrated on an airfoil with different piecewise nonlinearities. The obtained solutions are in good agreement with the numerical integration solutions. The sensitivities of the harmonic contents for LCOs to variation of the structural parameters are revealed. It is found that for the hysteresis model the continuation method is not a proper tool for the analysis of the nonlinear dynamics of an airfoil when one of the structural parameters such as the preload or freeplay parameter is varied.

Acknowledgments The author is grateful to the editor and anonymous referees for their careful reading and comments. This study has been financially supported by National Natural Science Foundation of China (Project no. 10904178). Appendix 1 a 2 c0 ¼ 1 þ ; c1 ¼ xα  h ; c2 ¼ xα  ; μ μ μ   1 2 c3 ¼ 1 þ ð1  2ah Þð1 ψ 1  ψ 2 Þ ; c4 ¼ ε1 ψ 1 þ ε2 ψ 2 ; μ μ  2 c5 ¼ 1  ψ 1 ψ 2 þ ð1=2  ah Þðε1 ψ 1 þ ε2 ψ 2 Þ ; μ     2 2 c6 ¼ ε1 ψ 1 1 ð1=2  ah Þε1 ; c7 ¼ ε2 ψ 2 1  ð1=2 ah Þε2 ; μ μ
2 2 2 ω c8 ¼  ε21 ψ 1 ; c9 ¼  ε22 ψ 2 ; c10 ¼ μ μ Un 1 þ 8a2h x α ah 1 þ2ah  2 ; d1 ¼ 1 þ ; d2 ¼  ð1 ψ 1  ψ 2 Þ; 2 r α μr α 8μr 2α μr 2α  1  2ah ð1 4a2h Þð1 ψ 1  ψ 2 Þ 2ζ α 1 þ 2ah   þ n ; d4 ¼  ε1 ψ 1 þ ε2 ψ 2 ; d3 ¼ 2μr 2α 2μr 2α μr 2α U d0 ¼

d5 ¼ 

1  4a2h 1 þ2ah ð1 ψ 1  ψ 2 Þ  ðε1 ψ 1 þ ε2 ψ 2 Þ; 2 μr α 2μr 2α

346

H. Liao / Journal of Fluids and Structures 55 (2015) 324–346

  ð1 þ 2ah Þε1 ψ 1  ð1 þ 2ah Þε2 ψ 2  1  ð1=2  ah Þε1 ; d7 ¼  1  ð1=2  ah Þε2 ; μr 2α μr 2α
2 ð1 þ 2ah Þε21 ψ 1 ð1 þ 2ah Þε22 ψ 2 1 ; d9 ¼ ; d10 ¼ d8 ¼ 2 2 μr α μr α Un d6 ¼ 

References Bittanti, S., Colaneri, P., 2009. Periodic Systems: Filtering and Control. Springer, London. Brown, B.M., Eastham, M.S.P., Schmidt, K.M., 2012. Periodic Differential Operators. Springer, Basel, Switzerland. Chen, Y.M., Liu, J.K., 2014. Nonlinear aeroelastic analysis of an airfoil-store system with a freeplay by precise integration method. Journal of Fluids and Structures 46, 149–164. Chen, Y.M., Liu, Q.X., Liu, J.K., 2014. Limit cycle analysis of nonsmooth aeroelastic system of an airfoil by extended variational iteration method. Acta Mechanica 225 (7), 2151–2159. Chung, K.W., Chan, C.L., Lee, B.H.K., 2007. Bifurcation analysis of a two-degree-of-freedom aeroelastic system with freeplay structural nonlinearity by a perturbation-incremental method. Journal of Sound and Vibration 299 (3), 520–539. Chung, K.W., He, Y.B., Lee, B.H.K., 2009. Bifurcation analysis of a two-degree-of-freedom aeroelastic system with hysteresis structural nonlinearity by a perturbation-incremental method. Journal of Sound and Vibration 320 (1), 163–183. Cui, C.C., Liu, J.K., Chen, Y.M., 2015. Simulating nonlinear aeroelastic responses of an airfoil with freeplay based on precise integration method. Communications in Nonlinear Science and Numerical Simulation 22 (1), 933–942. Dai, H.H., Yue, X.K., Yuan, J.P., Xie, D., 2014. A fast harmonic balance technique for periodic oscillations of an aeroelastic airfoil. Journal of Fluids and Structures 50, 231–252. Guo, H.L., Chen, Y.S., 2012. Dynamic analysis of two-degree-of-freedom airfoil with freeplay and cubic nonlinearities in supersonic flow. Applied Mathematics and Mechanics 33 (1), 1–14. Lee, H.K., Liu, L., Chung, K.W., 2005. Airfoil motion in subsonic flow with strong cubic nonlinear restoring forces. Journal of Sound and Vibration 281 (3), 699–717. Li, D., Guo, S., Xiang, J., 2012. Study of the conditions that cause chaotic motion in a two-dimensional airfoil with structural nonlinearities in subsonic flow. Journal of Fluids and Structures 33, 109–126. Liao, H.T., 2013. Uncertainty quantification and bifurcation analysis of an airfoil with multiple nonlinearities. Mathematical Problems in Engineering 2013. (Article ID 570947). Liao, H.T., Sun, W., 2013. A new method for predicting the maximum vibration amplitude of periodic solution of non-linear system. Nonlinear Dynamics 71 (3), 569–582. Liu, J.K., Chen, F.X., Chen, Y.M., 2012. Bifurcation analysis of aeroelastic systems with hysteresis by incremental harmonic balance method. Applied Mathematics and Computation 219 (5), 2398–2411. Liu, L., Wong, Y.S., Lee, B.H.K., 2002a. Nonlinear aeroelastic analysis using the point transformation method. Part1: Freeplay model. Journal of Sound and Vibration 253 (2), 447–469. Liu, L., Wong, Y.S., Lee, B.H.K., 2002b. Non-linear aeroelastic analysis using the point transformation method. Part 2: Hysteresis model. Journal of Sound and Vibration 253 (2), 471–483. Liu, L.P., Dowell, E.H., 2004. The secondary bifurcation of an aeroelastic airfoil motion: effect of high harmonics. Nonlinear Dynamics 37 (1), 31–49. Liu, L.P., Dowell, E.H., Thomas, J.P., 2007. A high dimensional harmonic balance approach for an aeroelastic airfoil with cubic restoring forces. Journal of Fluids and Structures 23 (7), 351–363. Narimani, A., Golnaraghi, M.F., Jazar, G.N., 2004. Sensitivity analysis of the frequency response of a piecewise linear system in a frequency island. Journal of Vibration and Control 10, 175–198. Sarrouy, E., Sinou, J.J., 2011. Nonlinear periodic and quasi-periodic vibrations in mechanical systems – on the use of the harmonic balance methods. In: Ebrahimi, Farzad (Ed.), Advances in Vibration Analysis Research, InTech Press. Wang, X., Hale, J.K., 2001. On monodromy matrix computation. Computer Methods in Applied Mechanics and Engineering 190 (18), 2263–2275. Zhang, S.J., Wen, G.L., Peng, F., Liu, Z.Q., 2013. Analysis of limit cycle oscillations of a typical airfoil section with freeplay. Acta Mechanica Sinica 29 (4), 583–592.