Transient and Steady-State Responses of an Asymmetric Nonlinear Oscillator

Abstract

We study the dynamical response of an asymmetric forced, damped Helmholtz-Duffing oscillator by using Jacobi elliptic functions, the method of elliptic balance, and Fourier series. By assuming that the modulus of the elliptic functions is slowly varying as a function of time and by considering the primary resonance response of the Helmholtz-Duffing oscillator, we derived an approximate solution that provides the time-dependent amplitude-frequency response curves. The accuracy of the derived approximate solution is evaluated by studying the evolution of the response curves of an asymmetric Duffing oscillator that describes the motion of a damped, forced system supported symmetrically by simple shear springs on a smooth inclined bearing surface. We also use the percentage overshoot value to study the influence of damping and nonlinearity on the transient and steady-state oscillatory amplitudes.

1. Introduction

Since most of the nonlinear differential equations that characterized the motion of several physical systems do not have closed-form solution, we have to use numerical or perturbation techniques to study the dynamical response of these systems. However, most of the perturbation techniques such as multiple scales, averaging, and harmonic balance, to say a few, focus on only the determination of steady-state approximate solutions because of the complexity involved in finding transient solutions of nonlinear differential equations [1, 2].

The aim of this paper is to investigate the influence on the dynamical behavior of the transient and steady-state solutions of the system 𝑥̈+2𝜈𝑥̇+𝐴𝑥+𝐵𝑥²+𝜀𝑥³=𝐹cos𝜔_𝑓𝑡 (1) close to the primary resonance region. Here, 𝑥 denotes the displacement of the system, 𝐴 is the natural frequency, 𝜈 is the damping coefficient, 𝜀 is a dimensionless nonlinear parameter, 𝜔_𝑓 is the driving frequency, 𝑡 is the running time, 𝐹 is the amplitude of the driving force, and 𝐵 is a system parameter. To solve the homogeneous Helmholtz-Duffing oscillator for which 𝐹=0 in (1), Hu used the harmonic balance method to calculate the first-order approximations to the periodic solutions of this equation [3]. Belhaq and Lakrad used the harmonic balance method involving Jacobian elliptic functions to obtain the approximate solution of (1) by taking 𝜈=0 and 𝜔_𝑓=0 [4]. Tamura [5] and Hu [6] developed the exact solution of a quadratic nonlinear oscillator that is part of (1) by using an elliptic function. Cao and coworkers investigated in [7] the various symmetry breaking phenomena associated with the Helmholtz-Duffing oscillator (1) in the case for which 𝐵=1−𝜀 and for different values of the so-called symmetric parameter 𝜀. They also used the second-order averaging method to investigate its local bifurcation behavior. By considering a rational form elliptic solution to (1) when 𝐹=0, Elías-Zúñiga derived its analytical solution which is similar in form to that of its exact solution when 𝜈=0 [8].

Recently, Kovacic and coworkers studied the primary resonance response of (1) by applying the harmonic balance method and derived nonlinear algebraic equations for the steady-state system response [9], while Jeyakumari et al. analyzed how the potential well of an asymmetric Duffing oscillator affects the vibrational resonance response [10].

Here the approximate solution of (1) is derived by taking into account the transient and the steady-state responses without the simplifications regarding undamped and unforced system included in previously developed solutions such as [4, 6]. The approximate solution is based on trigonometric and Jacobi elliptic functions with slowly varying parameters that will help us to obtain amplitude-frequency response curves that evolve with time. Then, the influence of the nonlinear transient responses is investigated since recent studies show that transient vibrations can not only provide additional information to fully predict the system stable behavior [11], but also can be used to predict the system overshoot value [12, 13]. The determination of this value is of practical interest in understanding the importance of time in controlling the dynamical behavior of oscillatory systems [14, 15] therefore, the percentage overshoot value will be computed by considering the influence of the system parameters such as nonlinear and damping effects.

In the next section, the approximate general solution of (1) is derived by using Jacobi elliptic functions.

2. Approximate Solution

In order to obtain the general solution of (1) in the region of the primary resonance, we will consider that this equation can be written in equivalent form as 𝑥̈+2𝜈𝑥̇+𝐴𝑥+𝐵𝑥²+𝜀𝑥³=𝐹cn𝜔_𝑓𝑡,𝑘²_𝑓, (2) where cn(𝜔_𝑓𝑡,𝑘²_𝑓) is the Jacobian elliptic function that has a period in 𝜔_𝑓𝑡 equal to 4𝐾(𝑘²_𝑓), and 𝐾(𝑘²_𝑓) is the complete elliptic integral of the first kind for the modulus 𝑘_𝑓 [16, 17]. Note that 0≤|𝑘²_𝑓|<1 must hold, and when the modulus of cn(𝜔_𝑓𝑡,𝑘²_𝑓) is zero, then cn(𝜔_𝑓𝑡,𝑘²_𝑓) and the trigonometric function cos(𝜔_𝑓𝑡) coincide and, thus, (2) reduces to (1). Now, we make the assumption that (2) has an approximate general solution of the form 𝑥=𝑏(𝑡)+𝑐(𝑡)cn𝜔₁𝑡+𝜙₁,𝑘²₁+𝐷₁(𝑡)cn𝜔_𝑓𝑡,𝑘²_𝑓+𝐷₂(𝑡)sn𝜔_𝑓𝑡,𝑘²_𝑓, (3) where 𝜙₁, 𝜔₁(𝑡), 𝑏(𝑡), 𝐷₁(𝑡), 𝐷₂(𝑡), 𝑘₁(𝑡), and 𝑘_𝑓(𝑡) need to be determined. For the sake of simplicity, it is assumed that 𝑏(𝑡), 𝜔₁(𝑡), 𝐷₁(𝑡), 𝐷₂(𝑡), 𝑘₁(𝑡), and 𝑘_𝑓(𝑡) are slowly varying as a function of time so that 𝑏̇, 𝜔̇₁, 𝐷̇₁, 𝐷̇₂, 𝑘̇₁, and 𝑘̇_𝑓 or higher order are small enough to be ignored [18]. Also, for simplicity, we will use the following notation cn_𝑓≡cn𝑢_𝑓=cn(𝜔_𝑓𝑡,𝑘²_𝑓), and cn₁≡cn𝑢₁=cn(𝜔₁𝑡+𝜙₁,𝑘²₁) with similar notation for the functions sn and dn. Thus, substitution of (3) into (2) gives 𝐴𝑏+𝑏²𝐵+𝑏³𝜀+cn₁𝐴𝑐+2𝑏𝐵𝑐+3𝑏²𝜀𝑐+cn_𝑓2𝐵𝑐𝐷₁+6𝑏𝜀𝑐𝐷₁+3𝜀𝑐𝐷²₁cn²_𝑓−𝜔²₁𝑐+2𝑘²₁𝜔²₁𝑐+sn_𝑓2𝐵𝑐𝐷₂+6𝑏𝜀𝑐𝐷₂+3𝜀𝑐𝐷²₂sn²_𝑓+2𝜈𝑐̇+𝑐̈+cn²₁𝐵𝑐²+3𝑏𝜀𝑐²+cn³₁𝜀𝑐³−2𝜔²₁𝑘²₁𝑐−2𝜔₁dn₁sn₁[𝜈𝑐+𝑐]̇+cn_𝑓dn_𝑓2𝜈𝜔_𝑓𝐷₂+cn_𝑓−𝐹+𝐴𝐷₁+2𝑏𝐵𝐷₁+3𝑏²𝜀𝐷₁+3𝜀𝑐²𝐷₁cn²₁+2𝜔²_𝑓𝑘²_𝑓𝐷₁−𝜔²_𝑓𝐷₁+sn_𝑓2𝐵𝐷₁𝐷₂+6𝑏𝜀𝐷₁𝐷₂+3𝜀𝐷₁𝐷²₂sn²_𝑓+cn²_𝑓𝐵𝐷²₁+3𝑏𝜀𝐷²₁+cn³_𝑓𝜀𝐷³₁−2𝜔²_𝑓𝑘²_𝑓𝐷₁−dn_𝑓sn_𝑓2𝜈𝜔_𝑓𝐷₁+cn₁cn_𝑓sn_𝑓6𝜀𝑐𝐷₁𝐷₂+sn_𝑓𝐴𝐷₂+2𝑏𝐵𝐷₂+3𝑏²𝜀𝐷₂+3𝜀𝑐²𝐷₂cn²₁+3𝜀𝐷²₁𝐷₂cn²_𝑓−𝑘²_𝑓𝜔²_𝑓𝐷₂−𝜔²_𝑓𝐷₂+sn²_𝑓𝐵𝐷²₂+3𝑏𝜀𝐷²₂+sn³_𝑓𝜀𝐷³₂+2𝜔²_𝑓𝑘²_𝑓𝐷₂=0 (4) in which the following identities sn²+cn²=1 and dn²+𝑘²_𝑗sn²=1 have been used [19].

Notice that (4) depends on elliptic functions with different modulus. To simplify (4), we first compute its average with respect to 4𝐾_𝑓 while keeping terms with period 4𝐾₁ constant; this yields the following equation: 𝐴𝑏+𝑏²𝐵+𝑏³𝜀+1−𝐼_𝑓𝐵𝐷²₂+3𝑏𝜀𝐷²₂+𝐼_𝑓𝐵𝐷²₁+3𝑏𝜀𝐷²₁+cn₁𝐴𝑐+2𝑏𝐵𝑐+3𝑏²𝜀𝑐+3𝜀𝑐𝐷²₁𝐼_𝑓−𝜔²₁𝑐+2𝑘²₁𝜔²₁𝑐+3𝜀𝑐𝐷²₂1−𝐼_𝑓+2𝜈𝑐̇+𝑐̈+cn²₁𝐵𝑐²+3𝑏𝜀𝑐²−2𝜔₁dn₁sn₁[𝜈𝑐+𝑐]̇+cn³₁𝜀𝑐³−2𝜔²₁𝑘²₁𝑐=0. (5) We now take the average of (5) with respect to 4𝐾₁, to get 𝐴𝑏+𝑏²𝐵+𝑏³𝜀+1−𝐼_𝑓𝐵𝐷²₂+3𝑏𝜀𝐷²₂+𝐼_𝑓𝐵𝐷²₁+3𝑏𝜀𝐷²₁+𝐼₁𝑐²(𝐵+3𝑏𝜀)=0, (6) where the 𝐼_𝑗 terms are defined in the appendix.

Next, the average of (4) is computed with respect to 4𝐾₁ while terms with period 4𝐾_𝑓 are held constant; this yields 𝐴𝑏+𝑏²𝐵+𝑏³𝜀+𝐼₁𝐵𝑐²+3𝑏𝜀𝑐²+cn_𝑓sn_𝑓2𝐵𝐷₁𝐷₂+6𝑏𝜀𝐷₁𝐷₂×cn_𝑓−𝐹+𝐴𝐷₁+2𝑏𝐵𝐷₁+3𝑏²𝜀𝐷₁+3𝜀𝑐²𝐷₁𝐼₁+2𝜔²_𝑓𝑘²_𝑓𝐷₁−𝜔²_𝑓𝐷₁+3𝜀𝐷₁𝐷²₂+cn²_𝑓𝐵𝐷²₁+3𝑏𝜀𝐷²₁+cn³_𝑓𝜀𝐷³₁−2𝜔²_𝑓𝑘²_𝑓𝐷₁−3𝜀𝐷₁𝐷²₂+cn_𝑓dn_𝑓2𝜈𝜔_𝑓𝐷₂−dn_𝑓sn_𝑓2𝜈𝜔_𝑓𝐷₁+sn²_𝑓𝐵𝐷²₂+3𝑏𝜀𝐷²₂+sn_𝑓𝐴𝐷₂+2𝑏𝐵𝐷₂+3𝑏²𝜀𝐷₂+3𝜀𝑐²𝐷₂𝐼₁+𝑘²_𝑓𝜔²_𝑓𝐷₂−𝜔²_𝑓𝐷₂+𝜀𝐷³₂−sn_𝑓cn²_𝑓𝜀𝐷³₂+2𝜔²_𝑓𝑘²_𝑓𝐷₂−3𝜀𝐷²₁𝐷₂=0. (7) Then, the average of (7) is computed with respect to 4𝐾_𝑓, to get the following expression: 𝐴𝑏+𝑏²𝐵+𝑏³𝜀+𝐼_𝑓𝐵𝐷²₁+3𝑏𝜀𝐷²₁−𝐵𝐷²₂−3𝑏𝜀𝐷²₂+𝐼₁𝐵𝑐²+3𝑏𝜀𝑐²+𝐵𝐷²₂+3𝑏𝜀𝐷²₂=0. (8)

To determine the unknown parameters 𝜙₁, 𝜔₁(𝑡), 𝑏(𝑡), 𝐷₁(𝑡), 𝐷₂(𝑡), 𝑘₁(𝑡), and 𝑘_𝑓(𝑡), Fourier series are applied to (7) to transform the terms sndn and sncn² into sn, cnsn and cndn into cn Jacobian elliptic functions, respectively [20–22]. Therefore, (7) simplifies to cn_𝑓−𝐹+𝐴𝐷₁+2𝑏𝐵𝐷₁+3𝑏²𝜀𝐷₁+3𝜀𝑐²𝐷₁𝐼₁+2𝜔²_𝑓𝑘²_𝑓𝐷₁−𝜔²_𝑓𝐷₁+3𝜀𝐷₁𝐷²₂+2𝜈𝜔_𝑓𝐷₂𝐻₂+cn²_𝑓𝐵𝐷²₁+3𝑏𝜀𝐷²₁+cn³_𝑓𝜀𝐷³₁−2𝜔²_𝑓𝑘²_𝑓𝐷₁−3𝜀𝐷₁𝐷²₂+sn²_𝑓𝐵𝐷²₂+3𝑏𝜀𝐷²₂+sn_𝑓𝐴𝐷₂+2𝑏𝐵𝐷₂+3𝑏²𝜀𝐷₂+3𝜀𝑐²𝐷₂𝐼₁+𝑘²_𝑓𝜔²_𝑓𝐷₂−𝜔²_𝑓𝐷₂+𝜀𝐷³₂−2𝜈𝜔_𝑓𝐷₁𝐻₁−𝜀𝐷³₂+2𝜔²_𝑓𝑘²_𝑓𝐷₂−3𝜀𝐷²₁𝐷₂𝐻₃+𝐴𝑏+𝑏²𝐵+𝑏³𝜀+𝐼₁𝐵𝑐²+3𝑏𝜀𝑐²+higherharmonics=0, (9) where 𝐻₁, 𝐻₂, and 𝐻₃ are defined in the appendix. If we now use the amplitude for Jacobi elliptic functions so that cos𝜑_𝑗=cn(𝜔_𝑗𝑡;𝑘²_𝑗)=cn_𝑗, and sin𝜑_𝑗=sn(𝜔_𝑗𝑡;𝑘²_𝑗)=sn_𝑗, then (5) and (9) become 𝐴𝑏+𝑏²𝐵+𝑏³𝜀+1−𝐼_𝑓𝐵𝐷²₂+3𝑏𝜀𝐷²₂+𝐼_𝑓𝐵𝐷²₁+3𝑏𝜀𝐷²₁+cos𝜑₁𝐴𝑐+2𝑏𝐵𝑐+3𝑏²𝜀𝑐+3𝜀𝑐𝐷²₁𝐼_𝑓−𝜔²₁𝑐+2𝑘²₁𝜔²₁𝑐+3𝜀𝑐𝐷²₂1−𝐼_𝑓+2𝜈𝑐̇+𝑐̈+cos²𝜑₁𝐵𝑐²+3𝑏𝜀𝑐²−2𝜔₁1−𝑘²₁sin𝜑²₁sin𝜑₁[𝜈𝑐+𝑐]̇+cos³𝜑₁𝜀𝑐³−2𝜔²₁𝑘²₁𝑐=0,cos𝜑_𝑓−𝐹+𝐴𝐷₁+2𝑏𝐵𝐷₁+3𝑏²𝜀𝐷₁+3𝜀𝑐²𝐷₁𝐼₁+2𝜔²_𝑓𝑘²_𝑓𝐷₁−𝜔²_𝑓𝐷₁+3𝜀𝐷₁𝐷²₂+2𝜈𝜔_𝑓𝐷₂𝐻₂+34𝜀𝐷³₁−2𝜔²_𝑓𝑘²_𝑓𝐷₁−3𝜀𝐷₁𝐷²₂+sin𝜑_𝑓𝐴𝐷₂+2𝑏𝐵𝐷₂+3𝑏²𝜀𝐷₂+3𝜀𝑐²𝐷₂𝐼₁+𝑘²_𝑓𝜔²_𝑓𝐷₂−𝜔²_𝑓𝐷₂+𝜀𝐷³₂−2𝜈𝜔_𝑓𝐷₁𝐻₁−𝜀𝐷³₂+2𝜔²_𝑓𝑘²_𝑓𝐷₂−3𝜀𝐷²₁𝐷₂𝐻₃+𝐴𝑏+𝑏²𝐵+𝑏³𝜀+𝐼₁𝐵𝑐²+3𝑏𝜀𝑐²+12𝐵𝐷²₁+3𝑏𝜀𝐷²₁+𝐵𝐷²₂+3𝑏𝜀𝐷²₂+higherharmonics=0. (10) Since our original system (1) is subjected to a driving force of sinusoidal type, then the modulus 𝑘_𝑓≡0 and the following identities hold: cn_𝑓=cos_𝑓, sn_𝑓=sin_𝑓, 𝐼_𝑓=1/2, 𝐻₁=1, and 𝐻₂=1. If we also assume that 0<|𝑘²₁|<1/2, then, from (A.1), 𝐼₁≃1/2 and (6) becomes exactly the same as (8). This provides the following simplifications for (10): cos𝜑₁3𝐴𝑐+4𝜀𝑐³−2𝜔²₁𝑘²₁𝑐+2𝑏𝐵𝑐+3𝑏²3𝜀𝑐+2𝜀𝑐𝐷²₁+𝐷²₂−𝜔²₁𝑐+2𝑘²₁𝜔²₁𝑐+2𝜈𝑐̇+𝑐̈−2𝜔₁sin𝜑₁[𝑐]̇+1𝜈𝑐+4cos3𝜑₁𝜀𝑐³−2𝜔²₁𝑘²₁𝑐+higherharmonics=0,cos_𝑓−𝐹+𝐴𝐷₁+2𝑏𝐵𝐷₁+3𝑏²𝜀𝐷₁+32𝜀𝑐²𝐷₁−𝜔²_𝑓𝐷₁+3𝜀𝐷₁𝐷²₂+2𝜈𝜔_𝑓𝐷₂+34𝜀𝐷³₁−3𝜀𝐷₁𝐷²₂+sin_𝑓𝐴𝐷₂+2𝑏𝐵𝐷₂+3𝑏²𝜀𝐷₂+32𝜀𝑐²𝐷₂−𝜔²_𝑓𝐷₂+𝜀𝐷³₂−2𝜈𝜔_𝑓𝐷₁−14𝜀𝐷³₂−3𝜀𝐷²₁𝐷₂+higherharmonics=0. (11)

We then follow the harmonic balance procedure and ignore higher harmonics terms in (11), to find the following expressions: 3𝐴𝑐+4𝜀𝑐³−2𝜔²₁𝑘²₁𝑐+2𝑏𝐵𝑐+3𝑏²+3𝜀𝑐2𝜀𝑐𝐷²₁+𝐷²₂−𝜔²₁𝑐+2𝑘²₁𝜔²₁𝑐+2𝜈𝑐̇+𝑐̈=0,−𝐹+𝐴𝐷₁+2𝑏𝐵𝐷₁+3𝑏²𝜀𝐷₁+32𝜀𝑐²𝐷₁−𝜔²_𝑓𝐷₁+34𝜀𝐷₁𝐷²₂+2𝜈𝜔_𝑓𝐷₂+34𝜀𝐷³₁=0,𝐴𝐷₂+2𝑏𝐵𝐷₂+3𝑏²𝜀𝐷₂+32𝜀𝑐²𝐷₂−𝜔²_𝑓𝐷₂+34𝜀𝐷³₂−2𝜈𝜔_𝑓𝐷₁+34𝜀𝐷²₁𝐷₂=0,𝜀𝑐³−2𝜔²₁𝑘²₁̇𝑐=0,𝜈𝑐+𝑐=0. (12)(13)(14)(15)(16) Notice that the variable 𝑐 may be determined by integration of (16); this yields 𝑐(𝑡)=𝐶exp(−𝜈𝑡), (17) where 𝐶 is a constant of integration. Substitution of (17) into (8), as well as (12)–(15), yields −𝐷²₂2(𝐵+3𝑏𝜀)=𝐴𝑏+𝑏²𝐵+𝑏³𝜀+12𝐶²()+𝐷exp−2𝜈𝑡)(𝐵+3𝑏𝜀₂₁2(),𝜔𝐵+3𝑏𝜀₂₁=𝐴+2𝑏𝐵+3𝑏²3𝜀+2𝜀𝐷²₁+𝐷²₂+𝜀𝐶²()−𝜈exp−2𝜈𝑡²4,−34𝜀𝐷²₁+𝐷²₂𝐹=−𝐷₁+𝐴+2𝑏𝐵+3𝑏²𝜀+2𝜈𝜔_𝑓𝐷₂𝐷₁+32𝜀𝐶²exp(−2𝜈𝑡)−𝜔²_𝑓,−34𝜀𝐷²₁+𝐷²₂=𝐴+2𝑏𝐵+3𝑏²𝜀+32𝜀𝐶²exp(−2𝜈𝑡)−𝜔²_𝑓−2𝜈𝜔_𝑓𝐷₁𝐷₂,𝑘²₁=𝜀𝐶²exp(−2𝜈𝑡)2𝜔²₁. (18)(19)(20)(21)(22) Then, for given parameters values of 𝜈, 𝜀, 𝐴, 𝐵, 𝐹, and 𝜔_𝑓, the unknown parameters 𝜙₁, 𝜔₁, 𝐶, 𝐷₁, 𝐷₂, 𝑏, and 𝑘₁ can be found from (18)–(22) and from the initial conditions (I. C.) 𝑥(0)=𝑥₁₀, 𝑥̇(0)=𝑥̇₁₀. If we assume that the value of the initial velocity 𝑥̇(0)=𝑥̇₁₀ is such that 𝜙₁=0 in (3), then our approximate solution for (2) can be written as 𝑥=𝑏(𝑡)+𝐶exp(−𝜈𝑡)cn𝜔₁𝑡,𝑘²₁+𝐷₁(𝑡)cos𝜔_𝑓𝑡+𝐷₂(𝑡)sin𝜔_𝑓𝑡. (23) Since 𝑥(0)=𝑥₁₀ at 𝑡=0, we get from (23) that 𝐶=𝑥₁₀−𝑏(0)−𝐷₁(0). (24) We next use (20) and (21) to find 𝑏(𝑡) and 𝐷₂(𝑡) as 𝐷₂(𝑡)=𝐹∓𝐹²−4𝜈𝜔_𝑓𝐷₁²4𝜈𝜔_𝑓,𝑏(𝑡)=1384𝐷₁𝜀𝜈²𝜔²_𝑓×−128𝐵𝐷₁𝜈²𝜔²_𝑓∓256𝐷₁𝜈²𝜔²_𝑓×64𝐵²𝐷₁𝜈²𝜔²_𝑓×+3𝜀32𝜈²𝜔²_𝑓𝐹±𝑅₁×−𝐷₁6𝜀𝐹²∓6𝜀𝐹𝑅₁+32𝜈²𝜔²_𝑓×2𝐴+3𝐶²exp(−2𝜈𝑡)𝜀−2𝜔²_𝑓^1/2,𝑅₁=𝐹−4𝐷₁𝜈𝜔_𝑓𝐹+4𝐷₁𝜈𝜔_𝑓, (25)(26)(27) where the sign in (25) and (26) must be appropriately chosen to ensure that the higher elliptic terms in (23) have small amplitudes relative to the leading ones [23]. Substituting (25) and (26) into (18), we can easily prove, after some algebraic computations, that the values of 𝐷₁(𝑡) are determined from a ninth-order polynomial equation. Once the values of 𝐷₁(𝑡) are known, we may use (19) to determine 𝜔₁(𝑡). Notice that when the running time 𝑡 increases, the term 𝑒^−𝜈𝑡 approaches to zero, and thus (18), (19), and (22)–(26) become time-independent. In this case, the values of 𝑏, 𝐷₁, 𝐷₂, 𝜔₁, and 𝑘₁ in (23) will remain constant. This condition provides the steady-state solution of (1).

We will next discuss the accuracy of our derived approximate solution by studying the dynamical response of an asymmetric Duffing oscillator by plotting the amplitude-frequency and the percentage overshoot response curves and show how the evolution of time influences the systems behavior.

3. Numerical Simulations

In order to verify the accuracy of our proposed solution described in (23), we study the dynamical response of a rigid body supported symmetrically by a simple shear spring that is sliding over a smooth inclined bearing surface [24]. The system under consideration contains the following critical conditions: (a) finite amplitude forced vibration and (b) damping, demonstrating the general nature of the proposed solution for the resonance response of an asymmetric Duffing oscillator. In accordance with Elías-Zúñiga and Beatty [24], the equation for the damped motion of the load is a nonlinear, ordinary differential equation of the forced Duffing type with a constant static shift term of the form 𝜎̈+2𝜈𝜎̇+𝜎+𝜀𝜎³=𝜎_𝑒1+𝜀𝜎²_𝑒+𝐹cos𝜔_𝑓𝑡, (28) where the dot denotes the derivative with respect to 𝑡, 𝜎 represents the amount of simple shear deformation, 𝜎_𝑒 denotes the amount of static shear deflection of the load, 𝜈 is the damping ratio, 𝐹 is a dimensionless driving force, and 𝑡 is the dimensionless running time. If we transform (28) relative to 𝜎_𝑒 by using 𝑥=𝜎−𝜎_𝑒, then (28) becomes similar to the Helmholtz-Duffing equation (1) where the parameters 𝐴 and 𝐵 are described by 𝐴≡1+3𝜀𝜎²_𝑒,𝐵≡3𝜀𝜎_𝑒. (29)

To find the amplitude-frequency response curves of (1), we select the following parameter values: 𝜀=0.02, 𝜎_𝑒=0.497, 𝐹=0.1, 𝜈=0.01, 𝑥₀=1, and 𝑥̇₀=0.022, and plot these curves at the values of 𝑡=1, 15, 30, 80, and 100. Figure 1 shows the evolution of these amplitude-frequency curves with time. In the curves showed in Figure 1, we have plotted 𝐴_𝑚(𝑡)=|𝑏(𝑡)+𝐷₁(𝑡)+𝐷₂(𝑡)| versus 𝜔_𝑓. These curves characterize the evolution of the general motion of the system for the selected set of system parameter values. In these curves, the blue dots represent unstable system behavior determined by following the bifurcation analysis described in [9, 24]. Notice from Figure 1, that when the running time 𝑡 increases, that is, 𝑡≈100, the unstable regions in the amplitude-frequency response curves become smaller since the influence of time in (18), (19), and (22)–(26) is almost negligible. In Figure 1, the brown dots represent the stable steady-state amplitude-frequency response curve that characterizes the damped, forced nonlinear motion of a body supported by viscohyperelastic shear mountings [24]. Figure 2 illustrates the amplitude versus time response curve by considering the parameters values of 𝜎_𝑒=0.497, 𝜀=0.02, 𝜈=0.01, 𝐹=0.1, and 𝜔_𝑓=1.1 and the initial conditions 𝑥₀=0 and 𝑥̇₀=0.0663 for which 𝜙₁=0. There, the solid lines represent the numerical integration solution of (28) while the dashed lines represent the approximate solution provided by (23). In this case, the values of 𝑏, 𝐶, 𝐷₁, 𝐷₂, 𝑘₁, and 𝜔₁ were obtained from (18), (19), (22), (24), (25), and (26). In fact, at 𝑡=0, we get that 𝑏(0)=−0.0086, 𝐶=0.5496, 𝐷₁(0)=−0.541, 𝐷₂(0)=0.0653, 𝑘₁(0)=0.05418, 𝑥̇(0)=0.0663, and 𝜔₁(0)=1.0144. As shown in Figures 2(a) and 2(c), the highest vibration amplitudes occur during the transient oscillations. This transient oscillatory behavior provides information to ensure that the system design constraints are not exceeded, even if the shear suspension system performs well during its steady-state behavior. Figure 2(b) shows the phase plane of the transient oscillatory motion on the time interval 0≤𝑡≤31 plotted by considering the following initial conditions of 𝑥₀=0 and 𝑥̇₀=0.0663. It is clear from Figure 2(b) that the transient oscillations experience a fast increase in their amplitude values followed by a decrease when 𝑡>31, as illustrated in Figure 2(c). When the system reaches its steady-state behavior, the system phase portrait, depicted in Figure 2(d), exhibits stable behavior about the steady-state equilibria condition.

Figure 1 

Evolutive amplitude-frequency response curves of a rigid body supported symmetrically by a simple shear spring and sliding over a smooth inclined bearing surface with moderate static shear deflection 𝜎𝑒=0.497 and parameter values of 𝜀=0.02, 𝐹=0.1, 𝜈=0.01, and initial conditions 𝑥0=1 and 𝑥̇0=−0.0091. At 𝑡=0, the initial values of 𝑏, 𝐶, 𝐷1, 𝐷2, 𝑘1, and 𝜔1 are given as 𝑏(0)=−0.0121, 𝐶=0.9951, 𝐷1(0)=0.0962, 𝐷2(0)=0, 𝑘1(0)=0.0902, and 𝜔1(0)=1.0154. The blue dot curves represent unstable system response.

(a)

(b)

(c)

(d)

Figure 2 

Phase portraits and amplitude-time response curves for parameter values of 𝜀=0.02, 𝜎𝑒=0.497, 𝜈=0.01, 𝐹=0.1, 𝜔𝑓=1.1, 𝐷1(0)=−0.541, 𝐷2(0)=0.0653, 𝑏(0)=−0.0086, 𝑘1(0)=0.05418, 𝜔1(0)=1.0144, 𝐶=0.5496, and initial conditions 𝑥(0)=0 and 𝑥̇(0)= 0.0663. (a) Amplitude-time transient curves on 𝑡[0,31.5]; (b) amplitude-time response curve on 𝑡[0,700]; (c) transient response phase portrait; (d) steady-state system response phase portrait on 𝑡[650,700]. The solid lines represent the numerical solution, and the dashed lines represent the proposed EBM solution.

We next investigate the influence of the transient solution on the dynamics response of (28) by introducing the percentage of overshoot which is defined as 𝑥%overshoot=^tran−𝑥_st𝑥_st×100, (30) where 𝑥_trans and 𝑥_st represent the peak transient and steady-state amplitudes of the system, respectively, that can be estimated by setting to zero the time derivative of (23). Figure 3 shows the calculated system overshoot plotted against the system angular frequency, 𝜔_𝑓, by considering the following system parameter values: 𝜎_𝑒=0.497, 𝐹=0.1, 𝜈=0.01, and 𝑥₀=0. As we can see from Figure 3, the influence of the transient solution on the system response is evident since the percentage of overshoot values is close to a 100% at the approximate frequency values of 0.1, 0.5, 1.5, and 2. Here the blue, the red, and the purple dashed lines were calculated by considering the nonlinear parameter values of 𝜀=0.02, 0.1, and 0.25, respectively. It is interesting to observe in Figure 3 the influence of the nonlinearity of the system in the overshoot qualitative and quantitative curves behavior. In fact, one can conclude from Figure 3 that the smallest percentage values of the overshoot occur with 𝜀=0.02 and near to the resonance frequency. Furthermore, the predicted overshoot curves oscillate about the value of 92% on the interval 1.5≤𝜔_𝑓≤2.7 that corresponds to the region depicted in Figure 4 at which the ratio of the transient and the steady-state amplitudes has almost the constant value of 52%.

Figure 3 

Percentage overshoot versus the driving frequency for parameter values of 𝜎𝑒=0.497, 𝐹=0.1, 𝜈=0.01, 𝑥0=0, and 𝜀=0.02, 0.1, and 0.25.

Figure 4 

Transient and steady-state amplitude ratio (𝑥tran/𝑥st) versus the driving frequency. Here the parameter values are 𝜎𝑒=0.497, 𝐹=0.1, 𝜈=0.01, 𝑥0=0, and 𝜀=0.02, 0.1, and 0.25.

To examine damping effects on the system overshoot, we next use the parameter values of 𝜎_𝑒=0.497, 𝐹=0.1, 𝜀=0.02, and 𝑥₀=0 and consider the damping ratio values of 𝜈=0.01, 0.05, and 0.1. The influence of the transient response on the peak overshoot values for the curve in which 𝜀=0.02, 𝜈=0.01, 𝜔_𝑓=0.5 and 1.5 is illustrated in Figure 5(a) and 5(b). Figure 5(c) shows the percentage overshoot curves plotted against the driving frequency 𝜔_𝑓. As me can see from Figures 5(c) and for a damping value of 𝜈=0.1, the overshoot values do not exceed 73%. It is also observed that for increasing damping values, the percentage of the system overshoot becomes smaller. Figures 5(d) and 5(e) illustrate, by using the Continuous Wavelet Transform, the influence of the transient and the driving frequencies on the nonlinear behavior of the system. It is evident from Figures 5(d) and 5(e) that at the beginning of the system motion, the transient frequency 𝜔₁ influences the fast increase of the system oscillations. Then, its effects on the system non-linear motion decrease when the running time exceeds the value of 𝑡≈200. Based on these findings, it is concluded that transient amplitude values are bigger than those of the steady-state ones as showed in Figures 3 and 5. Therefore, if one wants to reduce the impact of the transient oscillations amplitudes, we need, when possible, a fast tuning of the system to its operating bandwidth frequency region which can be achieved by using suitable control algorithms.

Figure 5 

System nonlinear dynamic behavior for parameter values of 𝜎𝑒=0.497, 𝐹=0.1, 𝜀=0.02, 𝑥0=0. Here, (a) and (c) show the amplitude versus time response curves for which the solid black and the dashed black lines represent the numerical and the EBM solutions, respectively; (b) illustrates the percentage overshoot versus the driving frequency for damping values of 𝜈=0.01, 0.05, and 0.1; (d) and (e) exhibit the influence of the transient and driving frequencies on the nonlinear system motion by using the Continuous Wavelet Transform.

4. Conclusions

In this paper, we have used elliptic functions to obtain the transient and the steady-state solutions of an asymmetric Duffing oscillator and studied its influence on the region of the primary resonance response. The theoretical predictions of our derived solution given by (23) are based on the assumption that the physical parameters 𝜔₁(𝑡), 𝐷₁(𝑡), 𝐷₂(𝑡), 𝑘₁(𝑡), and 𝑏(𝑡) are slowly varying as a function of time. To obtain computational tractable expressions to determine the aforementioned parameters, we have computed the average of the corresponding equations with respect to the complete elliptic integral of the first kind for the moduli 𝑘₁ and 𝑘_𝑓. Then, we have used Fourier series and followed the harmonic balance method to balance the lowest harmonic terms. These steps plus the assumption that 0<|𝑘²₁|<1/2 provide us with the time-dependent equations to determine the unknown parameters. The effect of the running time on the shapes of the amplitude-frequency response curves becomes evident during the study of the dynamical oscillator response with no linear stiffness term and with hardening characteristics. We have also shown that when the running time increases, the frequency range at which unstable motion occurs becomes smaller. Furthermore, we have obtained the percentage overshoot charts to address the influence of the nonlinear and damping terms on the transient and steady-state system responses when the system is suddenly switched on. We have found that when the nonlinearity is small, the percentage overshoot values becomes smaller at the overshoot bandwidth regions closer to the resonance frequency with increasing values of the steady-state and the transient amplitudes ratio. The usage of Continuous Wavelet Transform helps to identify the influence of the transient response signal on the system behaviour as time evolve. We also found that by adding damping to the system, the percentage of overshoot in the system does not exceed, in the case of 𝜈=0.1, 52% if the system is tuned into the overshoot bandwidth region of 0.5≤𝜔_𝑓≤1.5.

Based on these findings, we can conclude that if one really wants to have a better understanding of the system dynamical behavior, transient oscillations amplitudes must be estimated to avoid undesirable system effects such as peak transient amplitude values that could violate design constraints or unstable system behavior when the excitation is suddenly switched on. In this case, the percentage overshoot charts could be used to identify the frequency bandwidth region at which the system can be tuned to have a reliable dynamical response.

Appendix

The terms 𝐼_𝑗 in (6) can be computed from 𝐼_𝑗≡14𝐾_𝑗^{4𝐾𝑗+𝜙𝑗}_𝜙𝑗cn²_𝑗𝑑Ψ_𝑗1=1−𝑘²_𝑗𝐸1−^𝑗𝐾_𝑗;Ψ_𝑗≡𝜔_𝑗𝑡+𝜙_𝑗, (A.1) where 𝐸_𝑗 represents the complete elliptic integral of the second kind for the modulus 𝑘_𝑗 [20, 21].

The terms 𝐻₁, 𝐻₂, and 𝐻₃ in (9) are defined as 𝐻₁=1𝜋₀^2𝜋1sn𝑢dn𝑢sin𝜑𝑑𝜑=𝜋₀^4𝐾sn𝑢²dn𝑢²=4d𝑢3𝑘²_𝑓𝜋2𝑘²_𝑓−1𝐸_𝑓+𝐾_𝑓1−𝑘²_𝑓,𝐻₂=1𝜋₀^2𝜋1cn𝑢dn𝑢cos𝜑𝑑𝜑=𝜋₀^4𝐾cn𝑢²dn𝑢²=4d𝑢3𝑘²_𝑓𝜋1+𝑘²_𝑓𝐸_𝑓−𝐾_𝑓1−𝑘²_𝑓,𝐻₃=1𝜋₀^2𝜋sn𝑢cn𝑢²=1sin𝜑d𝜑𝜋₀^4𝐾sn𝑢²cn𝑢²1dn𝑢d𝑢=4. (A.2) Here 𝜑_𝑗=am(𝜔_𝑗𝑡;𝑘²_𝑗) is called the amplitude for Jacobi elliptic functions.

Acknowledgment

This work was partially funded by the Tecnológico de Monterrey, Campus Monterrey, through the Research Chairs in Nanotechnology and Intelligent Machines. Additional support was provided by the project FOMIX NL-2010-C30-145429.