# Study of the Response Prediction of a Parametric Vibration Based on Modulation Feedback

Authors: Dishan Huang
Journal: Advances in Acoustics and Vibration
 (2010)

---

## Abstract

The problem of response prediction is investigated for parametric vibration in terms of a new concept. The response solution is presented in the special form of Fourier series for signal degree freedom of parametric vibration based on modulation feedback. By applying harmonic balance and limitation operation, all coefficients of a harmonic component are fully determined with a set of series. Meanwhile, some important dynamic behaviors are exposed through mathematical deduction, and an instability phenomenon can be discussed through given frequency factors. The investigation result shows that the new approach has an advantage in the complete response expression, and it is very significant for the theoretical research and engineering application concerning parametric vibration.


---

## Body

## 1. Introduction

The problem of parametric vibration arises in many branches of physics and engineering, and investigation of stability and response prediction are the two most significant dynamic problems in the parametric vibration system. In the past, several methods have been used to study the stability of systems with periodic coefficients. These include Hill’s Method [1], the perturbation method [2], the averaging approach, Floquet theory with numerical integration [3], and Sinha’s numerical scheme with the shifted Chebyshev polynomial [4, 5].Hill’s infinite determinant was used to compute the transition curve or the instability boundary. However, it does not yield the response vector at an arbitrary point in the parametric space. The perturbation method supplies a method for obtaining the solution vector and instability boundary. The solution vector is expanded as a power series in terms of the small parametric and substituted into the original equation. The method is applicable only if the parametric multiplying the period term is small [6]. Authors [7–9] have tried to use the approximate approach to determine the stability and response from an approximate system of equations, which is usually obtained by replacing the elements of the periodic coefficients matrix with piecewise constants or linear functions. Sinha and Wu presented a numerical scheme for computing the Floquet transition matrix associated with a class of periodic systems. In the authors’ approach, the solution vector is expanded in terms of the shifted Chebyshev polynomials. This technique reduces the original differential system to a system of linear algebraic equation from which the solutions in the interval of one period can be easily obtained.In science and engineering practice, it will be very useful to express the response solution in the form of a Fourier series, which is convenient in an online vibration monitoring system. Especially in a mechanical fault diagnosis, such expression is important and not replaceable. One typical example is to diagnose the mechanical fault in a rotational rotor with a cross crack by analyzing a characteristic of the parametric system [10]. However, so far, no related response solution has directly been considered in all of the above approaches.In this paper, the concept of modulation feedback is introduced to set a response form that is similar to a Fourier series. Using harmonic balance in the equation and limitation operation, a full response expression is conducted. Though this concept comes from the content of a modulation system, it reflects the intrinsic physical character of a parametric system. Moreover, the presented solution conducts the response prediction that could be seen as an expansion of the Perturbation solution.


## 2. Modulation Feedback Conception

Consider a classical parametric vibration equation:(1)d2xdt2+ωn2(1+βcos⁡ωt)x=0,
and rewrite it as
(2)d2xdt2+ωn2x=-xωn2βcos⁡ωt.
Based on (2), the problem of parametric vibration can be schematically described as a system drawn in Figure 1, and it is a special feedback system that contains a second-order linear unit and a frequency modulation. The output of the system is the general solution of parametric vibration.Figure 1
Modulation feedback system.In this system, in addition to enhancing the input and feedback signal, the second-order linear unit will inherently generate the harmonic componentC01cosωnt after the moment t=0. However, the generated component C01cosωnt will be multiplied by ωn2βcosωt, and this operation is called the frequency modulation, which conducts frequency fission in the system.The phenomenon of frequency fission can be stated as the following: the componentωn is fissile into the combination of harmonic components ωn-ω and ωn+ω by a modulation operation within the first Δt. Then, the modulation result directly feeds back to the second-order linear unit as an input of the new harmonic components. Thus, within the second Δt, three components ωn,ωn-ω, and ωn+ω will be the output of the system and be involved in another fission operation, a sequential combination of five harmonic components ωn,ωn-ω,ωn-2ω, ωn+ω, and ωn+2ω will become the modulation results and will feed back to the second-order linear unit, and the iterate fission operation continues. The whole physical procedure of the frequency fission is schematically described in Figure 2.Figure 2
Frequency fission procedure (t≥0,Δt→0).Assuming the frequency fission procedure continues, except for the case of the instable state, the frequency fission will attain a state of balance. Therefore, there are many linear combinations of harmonic components related toωn and ω in the output of the system as the frequency fission balance, and they can be mathematically expressed as follows:
(3)x=∑k=-∞∞Ckcos⁡(ωn+kω)t.
Because the energy of the harmonic components is concentrated on the narrow band range in the frequency domain, the coefficient Ck→0 while k→∞.Therefore, the response solution problem of the parametric vibration equation is simplified to the determination of coefficient values ofCk in (3).


## 3. Fission Coordinate Equation

After substituting the solution (3) into (1), we obtain the following infinite set of linear, algebraic, homogeneous equations for the Ck by equating the coefficients of each harmonic on both sides of (1):(4)ωn2β2Ck-1+[ωn2-(ωn+kω)2]Ck+ωn2β2Ck+1=0(k=⋯-m,-(m-1),…,-3,-2,-1,0,1,2,3⋯m-1,m⋯).
Here, we introduce notations:(5)ϖk=ωn2-(ωn+kω)2(k=⋯-m,-(m-1),…,-3,-2,-1,0,1,2,3⋯m-1,m⋯),γ=ωn2β2.
First, we consider 2m+1 finite equations, and put all coefficient relations together and form a linear algebra equation called a fission coordinate equation (which is similar to Hill’s form).(6)[ϖ-mγγϖ-m+1γ⋅⋅⋅γϖ-3γγϖ-2γγϖ-1γγ0γγϖ1γγϖ2γγϖ3γ⋅⋅⋅γϖm-1γγϖm][C-mC-m+1⋅C-3C-2C-1C0C1C2C3⋅Cm-1Cm]=[-γC-(m+1)0⋅0000000⋅0-γCm+1].


## 4. Determination of CoefficientCk

To get a full response expression, the determination of coefficientCk is carried out by the following operation.(a) Consider the lower half of (6):(7)[ϖ1-ϖ-1γγϖ2γγϖ3γ⋅⋅⋅γϖm][C1C2C3⋅Cm]=[-γ(2C0+C-2)00⋅-γCm+1].
Let frequency factors be
(8)α1=ϖ1-ϖ-1,α2=ϖ2-γ2α1,⋮αm=ϖm-γ2α(m-1).
Then
(9)[α1γα2γα3γ⋅⋅⋅αm][C1C2C3⋅Cm]=[-γ(2C0+C-2)γ2(2C0+C-2)α1-γ3(2C0+C-2)α1α2⋅(-1)mγm(2C0+C-2)∏i=1m-1αi-γCm+1].
One of the solutions to (9):
(10)C1=(2C0+C-2)×(-γα1-γ3α12α2-γ5α12α22α3-⋯-γ2k-1(∏i=1k-1αi2)αk-⋯+(-1)m+2γm∏im-1αiCm+1).
As mentioned in the above section, Cm+1→0 while m→∞. Let k=n; thus the coefficient C1 can be expanded as
(11)C1=(2C0+C-2)S1,
where
(12)S1=-γα1-γ3α12α2-γ5α12α22α3-⋯-γ2n-1(∏i=1n-1αi2)αn-⋯.
For the same reason, C2,C3,…Ck can be determined. They are
(13)C2=(2C0+C-2)S2,(14)C3=(2C0+C-2)S3,⋮(15)Ck=(2C0+C-2)Sk,
where
(16)S2=γ2α1α2+γ4α1α22α3+γ6α1α22α32α4+⋯+γ2n(∏i=2nαi2)α1αn+1-⋯,S3=-γ3α1α2α3-γ5α1α2α32α4-γ7α1α2α32α42α5-⋯-γ2n+1(∏i=3n+1αi2)α1α2αn+2-⋯,⋮Sk=(-1)k∑n=1∞γ2n+k-2(∏j=1k-1αj)(∏i=kn+k-2αi2)αn+k-1.(b) Consider the upper half of (6):(17)[α-1γα-2γα-3γ⋅⋅⋅α-m][C-1C-2C-3⋅C-m]=[-γ(2C0+C2)γ2(2C0+C2)α1-γ3(2C0+C2)α1α2⋅(-1)mγm(2C0+C2)∏i=1m-1αi-γC-(m+1)],
where the frequency factors
(18)α-1=ϖ-1-ϖ1,α-2=ϖ-2-γ2α-1,⋮α-m=ϖ-m-γ2α-(m-1).
The solution of (17) is
(19)C-1=(2C0+C2)R1,(20)C-2=(2C0+C2)R2,(21)C-3=(2C0+C2)R3,⋮(22)C-k=(2C0+C2)Rk,
where
(23)R1=-γα-1-γ3α-12α-2-γ5α-12α-22α-3-⋯-γ2n-1(∏i=1n-1α-i2)α-n-⋯,(24)R2=γ2α-1α-2+γ4α-1α-22α-3+γ6α-1α-22α-32α-4+⋯+γ2n(∏i=2nα-i2)α-1α-(n+1)-⋯,(25)R3=-γ3α-1α-2α-3-γ5α-1α-2α-32α-4-γ7α-1α-2α-32α-42α-5-⋯-γ2n+1(∏i=3n+1α-i2)α-1α-2α-(n+2)-⋯,⋮(26)Rk=(-1)k∑n=1∞γ2n+k-2(∏j=1k-1α-j)(∏i=kn+k-2α-i2)α-(n+k-1).
From (11), (13), (19), and (20), we obtain(27)C-2=-2C0(1+R1S2R1+S1),C2=-2C0(1+S1R2S1+R1).
The final solution form of (4) is
(28)C-k=-2C0S1R2S1+R1Rk,(k=1,2,3,…),Ck=-2C0R1S2R1+S1Sk,(k=1,2,3,…).


## 5. Special Case

Whileω=2ωn/q (q is an integral number) and β>0, the response solution (3) is changed to
(29)x=∑k=-∞∞Ckcos⁡(q+2k)ωnqt.
If q is an even number, one gets
(30)x=∑k=-∞-(1+q/2)Ckcos⁡(q+2k)ωntq+∑k=-q/2∞Ckcos⁡(q+2k)ωntq=C-(q/2+1)cos⁡2ωntq+C-(q/2+2)cos⁡4ωntq+C-(q/2+3)cos⁡6ωntq+⋯+C-q/2+C-q/2+1cos⁡2ωntq+C-q/2+2cos⁡4ωntq+C-q/2+3cos⁡6ωntq+⋯=C-q/2+∑m=1∞(C-q/2-m+C-q/2+m)cos⁡2mωntq.
If q is an odd number, one gets
(31)x=∑k=-∞-(q+1)/2Ckcos⁡(q+2k)ωntq+∑k=-(q-1)/2∞Ckcos⁡(q+2k)ωntq=C-(q+1)/2cos⁡ωntq+C-(q+1)/2-1cos⁡3ωntq+C-(q+1)/2-2cos⁡5ωntq+⋯+C-(q-1)/2cos⁡ωntq+C-(q-1)/2+1cos⁡3ωntq+C-(q-1)/2+2cos⁡5ωntq+⋯=∑m=0∞(C-(q+1)/2-m+C-(q-1)/2+m)cos⁡(2m+1)ωntq.(a)
There is a direct componentC-q/2 in the response x while q is an even number in (29). This is a kind of particular characteristic exposed by the given approach in the response study of parametric vibration.(b)
Based on (30) and (31), the theoretical spectrum of the response of the parametric vibration can be easily predicted.(c)
Whilem=q/2 in (30) or m=(q-1)/2 in (31), if the coefficients
(32)C0+C-q≫C0,
the response will be the case of the instability state.


## 6. Complete Form of the Solution

### 6.1. Complete Form of the Solution

According to the theory of second-order differential equation, the strict form of a general solution for the given parametric vibration equation should be(33)x=∑k=-∞∞Ckcos⁡(ωn+kω)t+Dksin(ωn+kω)t.
Meanwhile, it needs to meet the sufficient and necessary condition that the Wronskian determinant is equal to the constant.In fact,(34)dw(t)dt=ddt|u1(t)u2(t)u̇1(t)u̇2(t)|=0.w(t) has no relation with time t, where
(35)u1(t)=∑k=-∞∞Ckcos⁡(ωn+kω)t,u2(t)=∑k=-∞∞Dksin(ωn+kω)t.
Because w(t) is a constant, (33) is the complete form of the general solution of the given parametric vibration equation, where C0 and D0 are the arbitrary constants in the solution in (33) and depend on the initial condition.


### 6.2. Determination of CoefficientDk

The coefficientDk is determined by using the harmonic balance as the determination of the coefficient value of Ck:
(36)D-k=-2D0S1R2S1+R1Rk,Dk=-2D0R1S2R1+S1Sk.
Thus, in the special case of ω=2ωn/q, if q is an even number, the complete form of the solution will be
(37)x=C-q/2+∑m=1∞(C-q/2-m+C-q/2+m)cos⁡2mωntq+D-q/2+∑m=1∞(-D-q/2-m+D-q/2+m)sin2mωntq.
If q is an odd number, the complete form of the solution will be
(38)x=∑m=0∞(C-(q+1)/2-m+C-(q-1)/2+m)cos⁡(2m+1)ωntq+∑m=0∞(-D-(q+1)/2-m+D-(q-1)/2+m)sin(2m+1)ωntq.


### 6.3. Initial Condition

Equation (33) is the parametric vibration solution after the system attains a state of balance through frequency fission. However, when t=0, equation (1) degenerates to a linear equation. Therefore, the coefficients C0 and D0 in (33) can be determined by
(39)C0=x|t=0,D0=1ωndxdt|t=0.


### 6.4. Comparison with Another Method

Ifk=-1,0,1 in (33), we obtain
(40)x=∑k=-11Ckcos⁡(ωn+kω)t+Dksin(ωn+kω)t.
This is a form of the parametric vibration solution conducted by singular Perturbation called the two variables expansion. Thus, the presented method conducts the response solution that is the expansion of the solution obtained by Perturbation.In fact, equation (6) is similar to Hill’s form, which was used in the study of the instability boundary of parametric vibration during that time [11]. However, this paper determines the coefficients of the harmonic series by using a special property of ω̅0=0 in (6).


## 6.1. Complete Form of the Solution

According to the theory of second-order differential equation, the strict form of a general solution for the given parametric vibration equation should be(33)x=∑k=-∞∞Ckcos⁡(ωn+kω)t+Dksin(ωn+kω)t.
Meanwhile, it needs to meet the sufficient and necessary condition that the Wronskian determinant is equal to the constant.In fact,(34)dw(t)dt=ddt|u1(t)u2(t)u̇1(t)u̇2(t)|=0.w(t) has no relation with time t, where
(35)u1(t)=∑k=-∞∞Ckcos⁡(ωn+kω)t,u2(t)=∑k=-∞∞Dksin(ωn+kω)t.
Because w(t) is a constant, (33) is the complete form of the general solution of the given parametric vibration equation, where C0 and D0 are the arbitrary constants in the solution in (33) and depend on the initial condition.


## 6.2. Determination of CoefficientDk

The coefficientDk is determined by using the harmonic balance as the determination of the coefficient value of Ck:
(36)D-k=-2D0S1R2S1+R1Rk,Dk=-2D0R1S2R1+S1Sk.
Thus, in the special case of ω=2ωn/q, if q is an even number, the complete form of the solution will be
(37)x=C-q/2+∑m=1∞(C-q/2-m+C-q/2+m)cos⁡2mωntq+D-q/2+∑m=1∞(-D-q/2-m+D-q/2+m)sin2mωntq.
If q is an odd number, the complete form of the solution will be
(38)x=∑m=0∞(C-(q+1)/2-m+C-(q-1)/2+m)cos⁡(2m+1)ωntq+∑m=0∞(-D-(q+1)/2-m+D-(q-1)/2+m)sin(2m+1)ωntq.


## 6.3. Initial Condition

Equation (33) is the parametric vibration solution after the system attains a state of balance through frequency fission. However, when t=0, equation (1) degenerates to a linear equation. Therefore, the coefficients C0 and D0 in (33) can be determined by
(39)C0=x|t=0,D0=1ωndxdt|t=0.


## 6.4. Comparison with Another Method

Ifk=-1,0,1 in (33), we obtain
(40)x=∑k=-11Ckcos⁡(ωn+kω)t+Dksin(ωn+kω)t.
This is a form of the parametric vibration solution conducted by singular Perturbation called the two variables expansion. Thus, the presented method conducts the response solution that is the expansion of the solution obtained by Perturbation.In fact, equation (6) is similar to Hill’s form, which was used in the study of the instability boundary of parametric vibration during that time [11]. However, this paper determines the coefficients of the harmonic series by using a special property of ω̅0=0 in (6).


## 7. Instability

There are infinite terms of frequency factorαk(k=…,-3,-2,-1,1,2,3,…) in the expression of the obtained solution, which causes instability in the system while αk=0. The frequency factors αk are related to ω and β, and the roots of αk=0 distribute around the instability areas, which are summarized as follows.


### 7.1. Principal Instability Area

Whileω=2ωn and k=-1 are selected, the frequency factor α-1=ωn2-(ωn-ω)2∣ω=2ωn=0. As a result, the values of γ/α-1 and γ3/α-12α-2 in (23) will tend to be infinite while β≠0. This means the harmonic response corresponding to C-1 tends to be infinite in the principal instability area.


### 7.2. High-Order Instability Boundary Curve

Based on the roots ofαk=ωn2-(ωn+kω)2-γ2/αk+1=0 (k<-1), the group of harmonic resonance (higher-order instability) boundary curves are drawn in the parameter plan of ω and β, as shown in Figures 3(a)–3(d) (generated by commercial software Maple 9.5), which have the following features.High-order instability boundary curve.
(a)
Parameter plan ofω and β while α-2=0(b)
Parameter plan ofω and β while α-3=0(c)
Parameter plan ofω and β while α-4=0(d)
Parameter plan ofω andβ while α-5=0Left Bias
Due to the termγ2/αk+1 having influence on frequency factor αk, each group of multiple harmonic resonance curves is under bias toward 2ωn/m (m=2,3,4,…, and k=-2,-3,-4,…). In other words, each harmonic resonance boundary curve is located at the left side of the point 2ωn/m that is the high-order instability point while β is close to 0.Instability Region Overlap
There are different order harmonic resonance boundary curves in each instability area around2ωn/m (m=1,2,3,…). For example, as shown in Figure 4, the multiple harmonic resonance regions exit near the third-order instability area, and they are based on α-3=0, α-4=0⋯; but for α-4=0, α-5=0⋯ almost overlap in the same region while 0<β<1. Wherein the curve α-3=0 is the left side of the instability region, and the curve α-k=0 (k→∞, in theory) is the right side of the instability region. Therefore, the harmonic resonance boundary curves (α-k=0) are composed of the boundaries of the instability region.Figure 4
Third instability region and boundary curve.Cross Direction
As known from (26), the harmonic resonance strength is different between αk→0 and αk+1→0. Thus, the harmonic resonance strength is sensitive to the ω cross direction (from left to right or from right to left) to the high-order instability area.


## 7.1. Principal Instability Area

Whileω=2ωn and k=-1 are selected, the frequency factor α-1=ωn2-(ωn-ω)2∣ω=2ωn=0. As a result, the values of γ/α-1 and γ3/α-12α-2 in (23) will tend to be infinite while β≠0. This means the harmonic response corresponding to C-1 tends to be infinite in the principal instability area.


## 7.2. High-Order Instability Boundary Curve

Based on the roots ofαk=ωn2-(ωn+kω)2-γ2/αk+1=0 (k<-1), the group of harmonic resonance (higher-order instability) boundary curves are drawn in the parameter plan of ω and β, as shown in Figures 3(a)–3(d) (generated by commercial software Maple 9.5), which have the following features.High-order instability boundary curve.
(a)
Parameter plan ofω and β while α-2=0(b)
Parameter plan ofω and β while α-3=0(c)
Parameter plan ofω and β while α-4=0(d)
Parameter plan ofω andβ while α-5=0Left Bias
Due to the termγ2/αk+1 having influence on frequency factor αk, each group of multiple harmonic resonance curves is under bias toward 2ωn/m (m=2,3,4,…, and k=-2,-3,-4,…). In other words, each harmonic resonance boundary curve is located at the left side of the point 2ωn/m that is the high-order instability point while β is close to 0.Instability Region Overlap
There are different order harmonic resonance boundary curves in each instability area around2ωn/m (m=1,2,3,…). For example, as shown in Figure 4, the multiple harmonic resonance regions exit near the third-order instability area, and they are based on α-3=0, α-4=0⋯; but for α-4=0, α-5=0⋯ almost overlap in the same region while 0<β<1. Wherein the curve α-3=0 is the left side of the instability region, and the curve α-k=0 (k→∞, in theory) is the right side of the instability region. Therefore, the harmonic resonance boundary curves (α-k=0) are composed of the boundaries of the instability region.Figure 4
Third instability region and boundary curve.Cross Direction
As known from (26), the harmonic resonance strength is different between αk→0 and αk+1→0. Thus, the harmonic resonance strength is sensitive to the ω cross direction (from left to right or from right to left) to the high-order instability area.


## 8. Example

Letωn=1, ω=0.5, and β=0.3 in (1), the initial conditions are x∣t=0=1 and dx/dt∣t=0=0.From (37), the response solution can be expressed as(41)x(t)=C-2+(C-3+C-1)cos⁡0.5t+(C-4+C0)cos⁡t+(C-5+C1)cos⁡1.5t+(C-6+C2)cos⁡2t+⋯.
On the initial condition, obtain C0=1.The calculation result of coefficient valuesCk and the corresponding frequency component ωn+kω are summarized in Table 1 in which the theoretical spectrum of the response solution of parametric vibration can be predicted. The spectrum is schematically shown in Figure 5.Table 1
Calculation result of the valueCk and corresponding frequency component ωn+kω (ωn=1,β=0.3, and ω=0.5).
KFrequencyωn+kωCk value⋯⋯⋯−8-3.0ωn-1.785e-7−7-2.5ωn-9.519e-6−6-2ωn-3.329e-4−5-1.5ωn−0.0067−4-ωn-0.0551−3-0.5ωn0.0067−200.0218−10.5ωn−0.15220ωn111.5ωn0.152222ωn0.007632.5ωn2.1798e-443.0ωn4.086e-6⋯⋯⋯Figure 5
Spectrum of the parametric vibration response solution (ωn=1, β=0.3, ω=0.5).The response solution is very sensitive toω while ω is closed to 0.5. If ω is a little less than 0.5, the high-order instability will be reflected in the solution (33). For example, in the case of ω=0.49845877 in the above example, we get C-4=2825.8. The result can be seen as a case of infinity, and it is governed by the left bias of higher-order instability. Simultaneously, the direct component C-2 cannot be neglected. The detailed results are listed in Table 2, and the theoretical response spectrum is drawn in Figure 6.Table 2
Calculation result of the valueCk and corresponding frequency component ωn+kω (ωn=1, β=0.3, and ω=0.49845877).
KFrequencyωn+kωCk value⋯⋯⋯−8≅-3ωn0.0096−7≅-2.5ωn0.5090−6≅-2ωn17.6233−5≅-1.5ωn347.62−4≅-ωn2825.8−3≅-0.5ωn−579.18−2≅087.8963−1≅0.5ωn−6.78870ωn11≅1.5ωn6.78872≅2ωn0.34133≅2.5ωn0.00984≅3ωn0.00018⋯⋯⋯Figure 6
Spectrum of the parametric vibration response solution in the case of higher-order harmonic resonance (ωn=1, β=0.3, and ω=0.49845877).


## 9. Conclusion

The response prediction of the single freedom of parametric vibration can be determined with a set of combinations of harmonic components related toωn and ω. Once the coefficients C0 and D0 are given, whether the general case or the special case of ω=2ωn/q, the harmonic components coefficients Ck and Dk can be fully determined by a set of series. On the coefficients Ck and Dk, the corresponding theoretical spectrum can be drawn; it will be very useful in identifying and predicting the characteristics in the mechanical system on a model of parametric vibration, especially in the case of ω=2ωn/q.The instability region can be determined when the frequency factorαk is equal to zero. The investigation result matches the case of Hill’s instability boundary curve. The harmonic resonance in the high-order instability region has three features, which are left bias, region overlap, and cross direction.The presented concept of modulation feedback reflects the physical inherency in the parametric vibration system, though it comes from the scope of a modulation system. This can become an available tool for us to analyze the response in the parametric vibration. Therefore, it is expected to obtain a response prediction of parametric vibration with an external force excitation and of a multiple freedom of parametric vibration, and so on.



---

*Source: 284187-2011-03-01.xml*