- Email: [email protected]
Spectral density of oscillator with bilinear stiffness and white noise excitation
Probabilistic Engineering Mechanics 18 (2003) 215–222 www.elsevier.com/locate/probengmech
Spectral density of oscillator with bilinear stiffness and white noise excitation Finn Ru¨dinger*, Steen Krenk Department of Mechanical Engineering, Technical University of Denmark, Building 101E, DK-2800 Kgs. Lyngby, Denmark Received 9 October 2002; revised 28 February 2003; accepted 1 April 2003
Abstract The power spectral density of an oscillator with bilinear stiffness excited by Gaussian white noise is considered. A method originally proposed by Krenk and Roberts [J Appl Mech 66 (1999) 225] relying on slowly changing energy for lightly damped systems is applied. In this method an approximate solution for the power spectral density at a given energy level is obtained by considering local similarity with the free undamped response. The total spectrum is obtained by integrating over all energy levels weighting each with the stationary probability density of the energy. The accuracy of the approximate analytical solution is demonstrated by comparing with results obtained by stochastic simulation. It is shown how the method successfully captures the broadening of the resonance peak and the presence of higher harmonics in the power spectral density of strongly non-linear systems. q 2003 Elsevier Ltd. All rights reserved. Keywords: Stochastic oscillator; Bilinear stiffness; Spectral density
1. Introduction In structural and mechanical engineering, problems involving unpredictable or stochastic variables or processes are frequently encountered, and in these cases a probabilistic analysis may be the most rational way of approaching the problem. In many problems involving the dynamical behaviour of mechanical systems, the dominating source of uncertainty or unpredictability is the excitation. Examples of such cases are wave and wind load on offshore structures, wind load on bridges and high rise buildings and earthquake excitation. If the excitation is given in terms of a stochastic process, the response of the mechanical system is also a stochastic process. If the excitation process is Gaussian and the system is linear, a fairly complete theory exists for the evaluation of the statistical properties of the response. If, on the other hand, the excitation is non-Gaussian or the system is non-linear, it is generally difficult to find exact solutions for the response statistics. In non-linear random vibration much effort has been directed towards the establishment of solutions in terms of * Corresponding author. Tel.: þ 45-45-25-19-67; fax: þ 45-45-88-43-25. E-mail address: [email protected] (F. Ru¨dinger). 0266-8920/03/$ - see front matter q 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0266-8920(03)00015-8
the probability density of the response. If the excitation is a Gaussian white noise, the state space vector is a Markov process, and the joint probability density of the state space variables can be obtained as the solution to the corresponding Fokker – Planck equation. Analytical solutions to this partial differential equation are generally difficult to find. A number of exact solutions are reported by e.g. Caughey [2] or Lin and Cai [3]. Most of these solutions focus on the stationary probability density of single-degree-of-freedom systems. Solutions in terms of the power spectral density of the response has attracted less attention, though this function yields important information concerning the dynamics of the system. The function is defined as the Fourier transform of the covariance function of the stationary response, and represents the distribution of energy on the frequency spectrum. For single-degree-of-freedom systems the width of the resonance peak of the power spectral density is a measure of the damping. Non-linearities in the stiffness will appear in the form of resonance peaks at the higher harmonics, and will furthermore tend to broaden the resonance peak, since the eigenfrequency of the system in this case is a function of energy level. Analytical solutions exist for linear systems with external excitation. For systems
216
F. Ru¨dinger, S. Krenk / Probabilistic Engineering Mechanics 18 (2003) 215–222
not belonging to this class only very few exact solutions can be obtained, [4,5]. Statistical linearization offers a classical way of obtaining the power spectral density of the response. In this method the non-linear equations of motion are replaced by a linear set of equations such that the variance of the difference between the two systems is minimized [6]. This procedure gives an accurate prediction of the mean frequency of the response, but only for small levels of non-linearity will the spectrum be accurately represented over a large part of the frequency range. For systems with a strongly non-linear behaviour, this method will generally give a poor representation of the spectrum, [7,8]. An early approach involving eigenfunction expansions and variational principles was proposed by Atkinson [9]. For systems with non-linear restoring force, the accuracy of this method seems to be decreasing with decreasing damping. Cai and Lin [10] have proposed a cumulant-neglect closure scheme, which shows good agreement with results obtained from digital simulation. However, the method relies on the convergence of the closure scheme. In a method proposed by Roy and Spanos [11] a perturbation method is used to approximate the spectrum of a Duffing oscillator. The method relies on Wiener-Hermite expansions of the non-linear terms in combination with a linearization procedure. In a later paper, [12], the power spectral density is expressed in terms of a power series expansion. The coefficients of the series are given by use of a Pade´ approximation and the accuracy of the method is shown to depend on the closure of the Pade´ approximation. The method is demonstrated for the singleand double-well Duffing oscillator and the Van der Pol oscillator and seems to give accurate results, though spurious peaks tend to appear in the spectrum for low levels of damping. An extension of the statistical linearization procedure was proposed by Miles [13]. In this method the natural frequency of the equivalent system is given as a function of the amplitude, which is a stochastic process. The spectrum is thereby evaluated as the linear spectrum for a given amplitude, and the total spectrum is evaluated by integrating over all amplitude levels weighting each by the probability density at that amplitude. The technique was demonstrated for a system with cubic stiffness (Duffing oscillator) and in a following paper for an oscillator with bilinear stiffness [14]. In both cases the system was assumed to be lightly damped, and the excitation was assumed to be Gaussian white noise. Methods of this type have been further investigated by Fogli et al. [15] and Fogli and Bressolette [8] considering an oscillator with bilinear stiffness. Soize [16] proposed a similar method, where the non-linear system is replaced by an equivalent linear system with stochastic stiffness and damping parameters. A joint probabilistic density of the stiffness and damping parameters is determined by assuming a parametric representation of the probability density, and requiring that the first and second moments of
the state space variables of the original system and the equivalent linear system are identical. The spectrum for given values of stiffness and damping is obtained as the linear spectrum. The total spectrum is then obtained by a weighted average using the joint probability density of the stiffness and damping parameters and the technique is demonstrated by considering the power spectral density of the Duffing oscillator. For a system with non-linear stiffness the eigenfrequency is a function of the amplitude. Since the stochastic response takes place at various levels of the amplitude, a shift of frequency will be observed in the power spectral density, i.e. a broadening of the resonance peak. The techniques mentioned above, Miles [13,14], Soize [16], Fogli et al. [15], Fogli and Bressolette [8], generally capture the broadening of the peak at the fundamental frequency quite well, but are not capable of capturing the resonance peaks at higher harmonics. A more accurate approximation to the power spectral density is obtained by a method proposed by Bouc [7], where the higher harmonics are accounted for via a perturbation technique. The method was demonstrated for the Duffing oscillator and in a later paper, [17], for an oscillator with bilinear asymmetric stiffness. In the present case a method proposed by Krenk and Roberts [1], originally applied to a system with linear-cubic damping, is investigated for a system with bilinear stiffness, i.e. a system with discontinuity in the tangential stiffness. Initially an exact solution for the stationary probability density of the mechanical energy is obtained. The free undamped vibration at a given energy level is expressed in terms of a set of modified phase plane variables conserving polar symmetry, which are expanded in Fourier series. It is shown how the free undamped response is approximated very accurately with only a few terms even in the case of systems with strong non-linearities. An approximation to the covariance function at a given energy level is obtained by assuming local similarity with the free undamped response. This approximate expression is also given in terms of a Fourier series expansion, and the spectrum at a given energy level is evaluated by applying the Fourier transform. The total spectrum is obtained by an integration over all energy levels weighting each with the probability density of the energy. This procedure is shown to give an accurate prediction of the power spectral density where inclusion of additional terms in the series expansion corresponds to including the resonance peaks at the higher harmonics.
2. Theory The theoretical background for the method is given in detail by Krenk and Roberts [1]. In this section a brief summary of the theory is presented. An oscillator with linear damping and non-linear stiffness is governed by an equation
F. Ru¨dinger, S. Krenk / Probabilistic Engineering Mechanics 18 (2003) 215–222
of motion of the following type,
can be expanded in sine and cosine series as
x€ þ gx_ þ gðxÞ ¼ WðtÞ
ð1Þ
where xðtÞ is the displacement and a dot indicates the derivative with respect to time. WðtÞ is the external excitation and is assumed to be a Gaussian white noise. g is the damping coefficient and gðxÞ is the stiffness function. The mechanical energy of the system is given by ðx Eðx; x_ Þ ¼ 12 x_ 2 þ GðxÞ; GðxÞ ¼ gðhÞdh ð2Þ 0
2
where ð1=2Þ_x is the kinetic energy and GðxÞ is the potential energy obtained by integration of the stiffness function.
The stationary joint probability density of the phase plane variables ðx; x_ Þ is given by, Caughey [2], 1 0 g GðxÞ þ 12 x_ 2 A px;_x ðx; x_ Þ ¼ C exp@2 ð3Þ pS0 where S0 is the intensity of the white noise process WðtÞ: C is a normalizing constant. It is seen, that the phase plane variables only enter the expression through the mechanical energy. The probability density of the mechanical energy is obtained from Eq. (3) as, gE pE ðEÞ ¼ TE px;_x ðEÞ ¼ CTE exp 2 ð4Þ pS0 where TE is the period of free undamped vibration at energy level E: TE is evaluated by the following integral, ðxmax xmin
dx pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ðE 2 GðxÞÞ
1 X z ðtÞ 2pjt p1ffiffiffiffi ¼ sin wðtÞ ¼ sj ðEÞsin TE 2E j¼1;3;… 1 X z ðtÞ 2pjt p2ffiffiffiffi ¼ cos wðtÞ ¼ cj ðEÞcos TE 2E j¼1;3;…
ð5Þ
xmin and xmax are the minimum and maximum displacement values for free undamped response at energy level E:
1 X k_x2 l ¼ 2kcos2 wl ¼ cj ðEÞ2 ¼ cðEÞ E j¼1;3;…
2E ¼ 2GðxÞ þ x_ 2 ¼ z21 þ z22 ¼ 2E sin2 w þ 2E cos2 w
ð6Þ
Free pffiffiffiffi undamped vibration thus describes a circle with radius 2E in the modified phase plane. z1 and z2 are given by pffiffiffiffiffiffiffi pffiffiffiffi z1 ¼ signðxÞ 2GðxÞ ¼ 2E sin w ð7Þ pffiffiffiffi z2 ¼ x_ ¼ 2E cos w with wt¼0 ¼ 0: The roles of sine and cosine are interchanged relative to the formulation by Krenk and Roberts [1]. In the case of free undamped vibration at energy level E; z1 and z2
ð9Þ
where k l indicates the mean value over a period of oscillation. It is observed that sðEÞ ¼ cðEÞ ¼ 2: In the case of a linear system s1 ðEÞ ¼ c1 ðEÞ ¼ 1; and all other coefficients vanish. 2.3. Spectral density The spectral density is defined as the Fourier transform of the covariance function. In the present case the one-sided spectral density is considered. It is defined as Sx ðvÞ ¼
1 ð1 E½xðtÞxðt þ tÞcosðvtÞdt p 21
ð10Þ
where E½ is the mean value. An approximation to the spectral density at the energy level l can be obtained as the sum of the spectral densities of each of the harmonics as
gcðEÞE p
Sx ðvlEÞ ¼
1 X j¼1;3;…
2cj ðEÞ2 ; ðj2 v2E 2 v2 Þ2 þ ðgcðEÞvÞ2
2.2. Modified state space variables In the case of a free undamped response, the energy level is constant. A set of modified phase plane variables ðz1 ; z2 Þ are introduced as
ð8Þ
since both variables are periodic with period TE : The coefficients of the sine and cosine expansions are related to the potential and kinetic energy. Taking the mean value of the square of the expressions (8) the following relations are obtained, 1 X k2GðxÞl ¼ 2ksin2 wl ¼ sj ðEÞ2 ¼ sðEÞ E j¼1;3;…
2.1. Probability density
TE ¼ 2
217
vE ¼
ð11Þ
2p TE
where cj ðEÞ are the Fourier coefficients introduced in Eq. (8). A more detailed derivation of Eq. (11) is given by Krenk and Roberts [1] by means of an approximate expression for the covariance function. The total spectrum is obtained by integrating with respect to the energy, weighting the spectrum at each energy level by the probability density of the energy, Sx ðvÞ ¼
ð1 0
Sx ðvlEÞpE ðEÞdE
ð12Þ
The method has also been extended to include systems with parametric excitation, Krenk et al. [18].
F. Ru¨dinger, S. Krenk / Probabilistic Engineering Mechanics 18 (2003) 215–222
218
3. Oscillator with bilinear stiffness A system with bilinear stiffness is considered. The stiffness function gðxÞ is shown in Fig. 1. Mathematically the function is described the following way, ( gðxÞ ¼
v21 x; v22 x
þ
aðv21
2
v22 ÞsignðxÞ;
lxl # a
ð13Þ
lxl . a
The initial stiffness is given by v21 : When the absolute value of the displacement exceeds a; the stiffness changes to v22 : It is assumed that both v1 ; v2 and a are larger than zero. For v1 ¼ v2 ; a ! 0 or a ! 1 the linear system is retrieved. The potential energy is obtained by integration of Eq. (13) as
GðxÞ ¼
8 < 12 v21 x2 ;
lxl # a
: 1 v2 x2 þaðv2 2 v2 Þlxlþ 1 a2 ðv2 2 v2 Þ; lxl . a 1 2 2 1 2 2 2
The potential energy is thus described by two parabolas with continuity of the zeroth and first derivative at the intersection. Since gðxÞ is an odd function, GðxÞ is even. 3.1. Probability density For free undamped vibration at energy level E; xmax is obtained by Gðxmax Þ ¼ E as E # 12 v21 a2 xmax ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > a > 2E 2 > 2 1 2 2 > : ð12 V Þþ V V 21þ a2 v2 ; E . 2 v1 a
ð15Þ
1
where the parameter V ¼ v1 =v2 will be referred to as the frequency ratio. Due to the symmetry in the potential
ð16Þ
1
The probability density of the energy is now determined directly from Eq. (4), and is conveniently investigated considering a non-dimensional formulation. Non-dimensional variables are introduced as
l¼ ð14Þ
8 pffiffiffiffi > 2E > > ; > < v1 a
energy, the natural period is determined as 4 ðxmax dx pffiffiffiffiffiffiffiffiffiffiffi TE ¼ pffiffi E2GðxÞ 2 0 8 2p > > ; E# 12 v21 a2 > > v > 1 > > > > > v1 a > < 4 sin21 p ffiffiffiffi 2E ¼ v1 > > 0 1 > > > > > V C > þ 4 cos21 B > @ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A; E. 12 v21 a2 > > v : 2 V2 21þ2E=a2 v2
E gE gv2 a2 ; E0 ¼ 12 v21 a2 ; b ¼ 0 ¼ 1 E0 pS0 2pS0
ð17Þ
E0 is thus a reference energy level and corresponds to the transition from the region with stiffness v21 to the region with stiffness v22 : For free undamped vibration the condition l #1 yields a harmonic vibration. b is a non-dimensional measure of the reference energy level relative to the input intensity. A non-dimensional natural period is defined as
vT T~ l ¼ 1 E 2p 8 1; l #1 > < ! ¼ 2 21 1 2V 21 V > pffiffi þ pffiffiffiffiffiffiffiffiffiffiffiffiffi ; l .1 cos : sin 2 p p l V 21þ l ð18Þ The probability density of the non-dimensional energy variable l is obtained from Eq. (4) as pl ðlÞ¼E0 pE ðE=E0 Þ¼C T~ l expð2blÞ
ð19Þ
It is thus seen that the probability distribution of l only depends on the damping parameter b and on the frequency ratio V: In order to verify the validity of the results comparison is made with records obtained by stochastic simulation. In Fig. 2 the probability density of l is shown for V ¼ 0:2 and b ¼ 0:2: The solid line represents the theoretical solution. The dots correspond to histogram points obtained by processing a simulated record. The results are seen to agree very well. 3.2. Free undamped vibration
Fig. 1. Bilinear stiffness function.
Free undamped vibration at energy level E is now considered. Only the range 0 # t # ð1=4ÞTE needs to be investigated due to symmetry. For E # ð1=2Þv21 a2 the system is linear and the phase plane variables ðx; x_ Þ are
F. Ru¨dinger, S. Krenk / Probabilistic Engineering Mechanics 18 (2003) 215–222
219
where t2 is the remaining part of the quarter period after the crossing of the level x ¼ a: t2 is determined from Eq. (16). Evaluation of A2 from Eq. (24a) gives the following value, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð25Þ A2 ¼ V a2 V2 2 a2 þ A21 The same result is obtained if continuity in x_ ðtÞ at time t ¼ t1 is used to determine this constant. The modified phase plane variables are determined from Eq. (6) as rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 2 2 2 z1 ðtÞ ¼ 2E 2 v2 A2 sin v2 4 TE 2 t ; ð26Þ 1 z2 ðtÞ ¼ v2 A2 sin v2 4 TE 2 t
Fig. 2. Probability density of non-dimensional energy, V ¼ 0:2; b ¼ 0:2:
simply given by xðtÞ ¼ A1 sinðv1 tÞ; pffiffiffiffi 2E A1 ¼ v1
x_ ðtÞ ¼ v1 A1 cosðv1 tÞ; ð20Þ
where A1 is the amplitude. For E . ð1=2Þv21 a2 the solution above is still valid until the time t1 where xðtÞ crosses into the range where the stiffness changes from v21 to v22 ; xðtÞ ¼ A1 sinðv1 tÞ;
x_ ðtÞ ¼ v1 A1 cosðv1 tÞ; 1 v1 a sin21 pffiffiffiffi 0 # t # t1 ¼ v1 2E
ð21Þ
For t . t1 the solution in terms of xðtÞ can be obtained by solving the following linear differential equation, x€ þ v22 x ¼ aðv22 2 v21 Þ
ð22Þ
The equation is non-homogeneous due to the non-zero equilibrium condition in the case where v1 – v2 : The initial conditions are given by Eq. (21) at t ¼ t1 since both xðtÞ and x_ ðtÞ are continuous. The solution is given by xðtÞ ¼ að1 2 V2 Þ þ A2 cos v2 14 TE 2 t ; ð23Þ x_ ðtÞ ¼ v2 A2 sin v2 14 TE 2 t ; t1 , t # 14 TE where only one arbitrary constant A2 appears, since the solution must have a form where the velocity is zero at t ¼ ð1=4ÞTE : The constant A2 is the amplitude of the harmonic part of the motion. Continuity in xðtÞ at t ¼ t1 yields aV2 ¼ A2 cosðv2 t2 Þ;
0 1 1 1 V B C t2 ¼ TE 2 t1 ¼ cos21 @ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A 4 v2 2 2 2 V 2 1þ 2E=a v 1
ð24Þ
Again, the solution is most easily investigated by introducing a non-dimensional formulation. The time t and z1 and z2 are made non-dimensional by the following rescaling, z ðtÞ 1 ~t ¼ v1 t; ~t1 ¼ v1 t1 ¼ sin21 pffiffi ; z~ j ðtÞ ¼ pjffiffiffiffi ð27Þ 2E l where the non-dimensional energy variable l is introduced in Eq. (17). z~ 1 ð~tÞ and z~2 ð~tÞ are expressed as 8 l #1_0# ~t # ~t1 sinð~tÞ; > > < vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! 1 ~ 2 z~1 ð~tÞ¼ u u > t12 V 21þ l sin2 2 pTl 2 ~t ; otherwise > : l V ð28Þ 8 cosð~tÞ; > > < sffiffiffiffiffiffiffiffiffiffiffiffiffi z~2 ð~tÞ¼ V2 21þ l > > sin : l
l #1_0# ~t # ~t1 1 ~ ~ 2 pTl 2 t
V
! ;
otherwise
Finally, the modified phase plane variables are expanded in a Fourier series. The Fourier coefficients sj ðlÞ and cj ðlÞ introduced in Eq. (8) can be evaluated by the following non-dimensional integrals, 4 ðð1=2ÞpT~ l j~t sj ðlÞ¼ z~1 ð~tÞsin d~t; ð29Þ ~ ~ p Tl 0 Tl 4 ðð1=2ÞpT~ l j~t z~2 ð~tÞcos cj ðlÞ¼ d~t ~ ~ p Tl 0 Tl where the non-dimensional period T~ l is defined in Eq. (18). Due to symmetry the integration can be reduced to the first quarter of the period and the result multiplied by 4. The above formulation shows that the coefficients sj ðlÞ and cj ðlÞ only depend on l; V and j: In the following examples the coefficients will be obtained by numerical integration. The Fourier series expansions are now compared to the exact solutions in order to estimate the number of terms necessary for an accurate prediction of the response. In Figs. 3 and 4 the non-dimensional phase plane variables z~1 and z~2 are shown as function of the non-dimensional time 4t=TE : The solid line is the analytical solution. The first
220
F. Ru¨dinger, S. Krenk / Probabilistic Engineering Mechanics 18 (2003) 215–222
conveniently expressed in the following non-dimensional formulation, Sx ðrÞv31 g~ ð1 ¼ cðlÞlpl ðlÞ E0 p 0 1 X 2cj ðlÞ2 £ dl; 2 2 2 ~ 22 2 2 2 j¼1;3;… ðj Tl 2r Þ þ g~ cðlÞ r g v g~ ¼ ; r ¼ ð30Þ v1 v1
Fig. 3. Modified state space variables, l ¼ 10; V ¼ 10; (—) exact, ( – –) 1 term, ( – · –) 3 terms.
modified phase plane variable z~1 starts at 0, while the second z~2 starts at one. In Fig. 3 the parameters are l ¼ 10 and V ¼ 10: The solution for an expansion with one term is shown by the dashed line and the solution for an expansion with three terms (j ¼ 1; 3; 5) is shown by the dashed/dotted line. Inclusion of only one term in the series expansion gives a rather inaccurate representation, while inclusion of three terms leads to an accurate representation of the exact function. In Fig. 4 the functions are shown for the parameters l ¼ 10 and V ¼ 0:2 and expansions with one term and three terms have been included in this case as well. The accuracy of the truncated series expansions is of the same order of magnitude as in the previous case. 3.3. Spectral density The spectral density is obtained from the expression (11) and (12). In the present case the spectral density is most
Fig. 4. Modified state space variables, l ¼ 10; V ¼ 0:2; (—) exact, (– –) 1 term, ( – · –) 3 terms.
In the following numerical examples the integration in this equation is performed numerically, so the equation takes the form of a double summation. From the above expression and from the earlier expressions for pl ðlÞ; cj ðlÞ and T~ l it is seen that the spectral density in this non-dimensional form is a function of V; b and g~ only. In Figs. 5 – 9 a number of numerical examples are shown for different combinations of the parameters V; b and g~: The dots indicate solutions obtained by stochastic simulation. In order to demonstrate the effect of the truncation of the series in (40) both the solution with one term and the solution with three terms (j ¼ 1; 3; 5) are included, indicated by a dashed line and a solid line, respectively. In Fig. 9 the solution with two terms (j ¼ 1; 3) is also included indicated by the dashed/dotted line. In Fig. 5 the values are chosen as V ¼ 0:6; b ¼ 1 and g~ ¼ 0:05: The degree of non-linearity in this case is relatively small, as seen by the magnitude of the peak of the higher harmonic. It is clearly seen that inclusion of only one term in the solution fails to capture the second peak, and that the additional terms have no influence (are negligibly small) near the fundamental frequency. The analytical results are seen to agree very well with the results obtained by stochastic simulation, except for the fact that the theory predicts a peak at v ¼ v1 : This corresponds to oscillations taking place at small energy levels where the response does not enter into the region of changed stiffness v22 : However, this peak is not as pronounced in the simulated
Fig. 5. Spectral density, V ¼ 0:6; b ¼ 1; g~ ¼ 0:05: (X) Stochastic simulation, ( – –) theory (1 term), (—) theory (3 terms).
F. Ru¨dinger, S. Krenk / Probabilistic Engineering Mechanics 18 (2003) 215–222
221
Fig. 6. Spectral density, V ¼ 0:2; b ¼ 1; g~ ¼ 0:05: (X) Stochastic simulation, (– –) theory (1 term), (—) theory (3 terms).
Fig. 8. Spectral density, V ¼ 5; b ¼ 0:5; g~ ¼ 0:05: (X) Stochastic simulation, ( – –) theory (1 term), (—) theory (3 terms).
response. This observation was also made by Fogli and Bressolette [8]. In Fig. 6 the level of non-linearity is increased by decreasing the frequency ratio to V ¼ 0:3; and letting the other parameters remain unchanged. The bandwidth of the fundamental frequency is now larger (due to the frequency sweep), and both the magnitude and bandwidth of the higher harmonics are increased. Again, the pronounced peak predicted by the theoretical result does not show in the spectrum of the simulated record. Furthermore, a small deviation is seen at zero frequency, a feature which is also observable in results given by Fogli and Bressolette [8], Krenk and Roberts [1]. A discrepancy reported by Krenk and Roberts [1] at high frequencies could be due to truncation after the first two terms. In the above examples the system is characterized by hardening stiffness, i.e. v2 . v1 or V , 1: In Figs. 7 and 8 systems with softening stiffness are considered. In Fig. 7 the values are chosen as V ¼ 2; b ¼ 0:5 and g~ ¼ 0:05:
As opposed to the two previous cases, the peak is now located below v ¼ v1 ; due to the fact that v2 , v1 : The non-linearity is relatively small, as seen by the magnitude of the higher harmonic. Again a small deviation is seen at zero frequency. In Fig. 8 the non-linearity is increased by increasing the value of the frequency ratio to V ¼ 5; while the other parameters remain unchanged. This is seen to broaden both the fundamental frequency peak and the peak of the first higher harmonic. The discrepancy at zero frequency is larger in this case compared to the previous, and the accuracy a bit lower, though still satisfactory. The peak at v ¼ v1 observed in Figs. 5 and 6 in the theoretical solution is less pronounced for the results shown in Figs. 7 and 8. Finally the parameter combination V ¼ 0:3; b ¼ 1 and g~ ¼ 0:01 is investigated. This relatively strong non-linearity in combination with a decreasing value of the damping parameter g~ results in a spectrum, where it is possible to distinguish both the second and third higher harmonic.
Fig. 7. Spectral density, V ¼ 2; b ¼ 0:5; g~ ¼ 0:05: (X) Stochastic simulation, (– –) theory (1 term), (—) theory (3 terms).
Fig. 9. Spectral density, V ¼ 0:3; b ¼ 1; g~ ¼ 0:01: (X) Stochastic simulation, ( – –) theory (1 term), ( – · –) theory (2 terms), (—) theory (3 terms).
222
F. Ru¨dinger, S. Krenk / Probabilistic Engineering Mechanics 18 (2003) 215–222
In Fig. 9 both the solutions with two terms and three terms have been included, and the effect of including the third term on the third harmonic is quite evident. From this plot it seems that the deviations observed by Krenk and Roberts [1] at high frequencies are due to neglecting the third term.
4. Conclusion An exact solution for the probability density of the energy and an approximate analytical expression for the power spectral density of an oscillator with bilinear stiffness and Gaussian white noise excitation has been derived following a method proposed by Krenk and Roberts [1]. The result has been presented in a convenient non-dimensional form, and the accuracy of the theory has been investigated by comparing with results obtained from simulated records. The approximation of the free undamped response by a series expansion is shown to be very accurate. Even in the case of a strong non-linearity in the stiffness only a few terms are needed to represent the response. This accuracy is reflected in the approximation of the power spectral density, where each term in the series expansion corresponds to including one additional resonance peak. The examples clearly show the presence of higher harmonics and the broadening of the resonance peak and demonstrate the ability of the theory to capture these characteristics.
Acknowledgements The project is supported by the Danish Technical Research Council. Part of the work was carried out while the first author was a Marie Curie grant holder at the European Laboratory for Structural Assessment, Joint Research Centre of the European Commission at Ispra, Italy. This support is gratefully acknowledged.
References [1] Krenk S, Roberts JB. Local similarity in non-linear random vibration. J Appl Mech 1999;66:225 –35. [2] Caughey TK. Nonlinear theory of random vibrations. Advances in applied mechanics. New York: Academic Press; 1971. p. 209–53. [3] Lin YK, Cai GQ. Probabilistic structural dynamics. New York: McGraw-Hill; 1995. [4] Dimentberg MF, Hou Z, Noori M. Spectral density of a non-linear single-degree-of-freedom system’s response to a white-noise random excitation: a unique case of an exact solution. Int J Non-linear Mech 1995;30:673–6. [5] Dimentberg MF, Lin YK. An exact solution for response spectral density of a single-degree-of-freedom system under both parametric and additive white noise excitation. J Appl Mech 2002;69:399–400. [6] Roberts JB, Spanos PD. Random vibration and statistical linearization. Chichester: Wiley; 1990. [7] Bouc R. The power spectral density of response for a strongly non-linear random oscillator. J Sound Vibration 1994;175:317– 31. [8] Fogli M, Bressolette PH. Spectral response of a stochastic oscillator under impacts. Meccanica 1997;32:1–12. [9] Atkinson JD. Eigenfunction expansions for randomly excited nonlinear systems. J Sound Vibration 1973;30:153 –72. [10] Cai GQ, Lin YK. Response spectral densities of strongly nonlinear systems under random excitation. Probabilistic Engng Mech 1997;12: 41– 7. [11] Roy RV, Spanos PD. Wiener-hermite functional representation of nonlinear stochastic systems. Struct Safety 1989;6:187–202. [12] Roy RV, Spanos PD. Power spectral density of nonlinear system response: the recursion method. J Appl Mech 1993;60:358– 65. [13] Miles RN. An approximate solution for the spectral response of Duffing’s oscillator with random input. J Sound Vibration 1989;132: 43– 9. [14] Miles RN. Spectral response of a bilinear oscillator. J Sound Vibration 1993;163:319–26. [15] Fogli M, Bressolette P, Bernard P. The dynamics of a stochastic oscillator with impacts. Eur J Mech A: Solids 1996;15:213–41. [16] Soize C. Stochastic linearization method with random parameters for SDOF nonlinear dynamical systems: prediction and identification procedures. Probabilistic Engng Mech 1995;10:143–52. [17] Bellizzi S, Bouc R. Spectral response of asymmetrical random oscillators. Probabilistic Engng Mech 1996;11:51–9. [18] Krenk S, Lin YK, Ru¨dinger F. Effective system properties and spectral density in random vibration. J Appl Mech 2002;69:161 –70.
Spectral density of oscillator with bilinear stiffness and white noise excitation Finn Ru¨dinger*, Steen Krenk Department of Mechanical Engineering, Technical University of Denmark, Building 101E, DK-2800 Kgs. Lyngby, Denmark Received 9 October 2002; revised 28 February 2003; accepted 1 April 2003
Abstract The power spectral density of an oscillator with bilinear stiffness excited by Gaussian white noise is considered. A method originally proposed by Krenk and Roberts [J Appl Mech 66 (1999) 225] relying on slowly changing energy for lightly damped systems is applied. In this method an approximate solution for the power spectral density at a given energy level is obtained by considering local similarity with the free undamped response. The total spectrum is obtained by integrating over all energy levels weighting each with the stationary probability density of the energy. The accuracy of the approximate analytical solution is demonstrated by comparing with results obtained by stochastic simulation. It is shown how the method successfully captures the broadening of the resonance peak and the presence of higher harmonics in the power spectral density of strongly non-linear systems. q 2003 Elsevier Ltd. All rights reserved. Keywords: Stochastic oscillator; Bilinear stiffness; Spectral density
1. Introduction In structural and mechanical engineering, problems involving unpredictable or stochastic variables or processes are frequently encountered, and in these cases a probabilistic analysis may be the most rational way of approaching the problem. In many problems involving the dynamical behaviour of mechanical systems, the dominating source of uncertainty or unpredictability is the excitation. Examples of such cases are wave and wind load on offshore structures, wind load on bridges and high rise buildings and earthquake excitation. If the excitation is given in terms of a stochastic process, the response of the mechanical system is also a stochastic process. If the excitation process is Gaussian and the system is linear, a fairly complete theory exists for the evaluation of the statistical properties of the response. If, on the other hand, the excitation is non-Gaussian or the system is non-linear, it is generally difficult to find exact solutions for the response statistics. In non-linear random vibration much effort has been directed towards the establishment of solutions in terms of * Corresponding author. Tel.: þ 45-45-25-19-67; fax: þ 45-45-88-43-25. E-mail address: [email protected] (F. Ru¨dinger). 0266-8920/03/$ - see front matter q 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0266-8920(03)00015-8
the probability density of the response. If the excitation is a Gaussian white noise, the state space vector is a Markov process, and the joint probability density of the state space variables can be obtained as the solution to the corresponding Fokker – Planck equation. Analytical solutions to this partial differential equation are generally difficult to find. A number of exact solutions are reported by e.g. Caughey [2] or Lin and Cai [3]. Most of these solutions focus on the stationary probability density of single-degree-of-freedom systems. Solutions in terms of the power spectral density of the response has attracted less attention, though this function yields important information concerning the dynamics of the system. The function is defined as the Fourier transform of the covariance function of the stationary response, and represents the distribution of energy on the frequency spectrum. For single-degree-of-freedom systems the width of the resonance peak of the power spectral density is a measure of the damping. Non-linearities in the stiffness will appear in the form of resonance peaks at the higher harmonics, and will furthermore tend to broaden the resonance peak, since the eigenfrequency of the system in this case is a function of energy level. Analytical solutions exist for linear systems with external excitation. For systems
216
F. Ru¨dinger, S. Krenk / Probabilistic Engineering Mechanics 18 (2003) 215–222
not belonging to this class only very few exact solutions can be obtained, [4,5]. Statistical linearization offers a classical way of obtaining the power spectral density of the response. In this method the non-linear equations of motion are replaced by a linear set of equations such that the variance of the difference between the two systems is minimized [6]. This procedure gives an accurate prediction of the mean frequency of the response, but only for small levels of non-linearity will the spectrum be accurately represented over a large part of the frequency range. For systems with a strongly non-linear behaviour, this method will generally give a poor representation of the spectrum, [7,8]. An early approach involving eigenfunction expansions and variational principles was proposed by Atkinson [9]. For systems with non-linear restoring force, the accuracy of this method seems to be decreasing with decreasing damping. Cai and Lin [10] have proposed a cumulant-neglect closure scheme, which shows good agreement with results obtained from digital simulation. However, the method relies on the convergence of the closure scheme. In a method proposed by Roy and Spanos [11] a perturbation method is used to approximate the spectrum of a Duffing oscillator. The method relies on Wiener-Hermite expansions of the non-linear terms in combination with a linearization procedure. In a later paper, [12], the power spectral density is expressed in terms of a power series expansion. The coefficients of the series are given by use of a Pade´ approximation and the accuracy of the method is shown to depend on the closure of the Pade´ approximation. The method is demonstrated for the singleand double-well Duffing oscillator and the Van der Pol oscillator and seems to give accurate results, though spurious peaks tend to appear in the spectrum for low levels of damping. An extension of the statistical linearization procedure was proposed by Miles [13]. In this method the natural frequency of the equivalent system is given as a function of the amplitude, which is a stochastic process. The spectrum is thereby evaluated as the linear spectrum for a given amplitude, and the total spectrum is evaluated by integrating over all amplitude levels weighting each by the probability density at that amplitude. The technique was demonstrated for a system with cubic stiffness (Duffing oscillator) and in a following paper for an oscillator with bilinear stiffness [14]. In both cases the system was assumed to be lightly damped, and the excitation was assumed to be Gaussian white noise. Methods of this type have been further investigated by Fogli et al. [15] and Fogli and Bressolette [8] considering an oscillator with bilinear stiffness. Soize [16] proposed a similar method, where the non-linear system is replaced by an equivalent linear system with stochastic stiffness and damping parameters. A joint probabilistic density of the stiffness and damping parameters is determined by assuming a parametric representation of the probability density, and requiring that the first and second moments of
the state space variables of the original system and the equivalent linear system are identical. The spectrum for given values of stiffness and damping is obtained as the linear spectrum. The total spectrum is then obtained by a weighted average using the joint probability density of the stiffness and damping parameters and the technique is demonstrated by considering the power spectral density of the Duffing oscillator. For a system with non-linear stiffness the eigenfrequency is a function of the amplitude. Since the stochastic response takes place at various levels of the amplitude, a shift of frequency will be observed in the power spectral density, i.e. a broadening of the resonance peak. The techniques mentioned above, Miles [13,14], Soize [16], Fogli et al. [15], Fogli and Bressolette [8], generally capture the broadening of the peak at the fundamental frequency quite well, but are not capable of capturing the resonance peaks at higher harmonics. A more accurate approximation to the power spectral density is obtained by a method proposed by Bouc [7], where the higher harmonics are accounted for via a perturbation technique. The method was demonstrated for the Duffing oscillator and in a later paper, [17], for an oscillator with bilinear asymmetric stiffness. In the present case a method proposed by Krenk and Roberts [1], originally applied to a system with linear-cubic damping, is investigated for a system with bilinear stiffness, i.e. a system with discontinuity in the tangential stiffness. Initially an exact solution for the stationary probability density of the mechanical energy is obtained. The free undamped vibration at a given energy level is expressed in terms of a set of modified phase plane variables conserving polar symmetry, which are expanded in Fourier series. It is shown how the free undamped response is approximated very accurately with only a few terms even in the case of systems with strong non-linearities. An approximation to the covariance function at a given energy level is obtained by assuming local similarity with the free undamped response. This approximate expression is also given in terms of a Fourier series expansion, and the spectrum at a given energy level is evaluated by applying the Fourier transform. The total spectrum is obtained by an integration over all energy levels weighting each with the probability density of the energy. This procedure is shown to give an accurate prediction of the power spectral density where inclusion of additional terms in the series expansion corresponds to including the resonance peaks at the higher harmonics.
2. Theory The theoretical background for the method is given in detail by Krenk and Roberts [1]. In this section a brief summary of the theory is presented. An oscillator with linear damping and non-linear stiffness is governed by an equation
F. Ru¨dinger, S. Krenk / Probabilistic Engineering Mechanics 18 (2003) 215–222
of motion of the following type,
can be expanded in sine and cosine series as
x€ þ gx_ þ gðxÞ ¼ WðtÞ
ð1Þ
where xðtÞ is the displacement and a dot indicates the derivative with respect to time. WðtÞ is the external excitation and is assumed to be a Gaussian white noise. g is the damping coefficient and gðxÞ is the stiffness function. The mechanical energy of the system is given by ðx Eðx; x_ Þ ¼ 12 x_ 2 þ GðxÞ; GðxÞ ¼ gðhÞdh ð2Þ 0
2
where ð1=2Þ_x is the kinetic energy and GðxÞ is the potential energy obtained by integration of the stiffness function.
The stationary joint probability density of the phase plane variables ðx; x_ Þ is given by, Caughey [2], 1 0 g GðxÞ þ 12 x_ 2 A px;_x ðx; x_ Þ ¼ C exp@2 ð3Þ pS0 where S0 is the intensity of the white noise process WðtÞ: C is a normalizing constant. It is seen, that the phase plane variables only enter the expression through the mechanical energy. The probability density of the mechanical energy is obtained from Eq. (3) as, gE pE ðEÞ ¼ TE px;_x ðEÞ ¼ CTE exp 2 ð4Þ pS0 where TE is the period of free undamped vibration at energy level E: TE is evaluated by the following integral, ðxmax xmin
dx pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ðE 2 GðxÞÞ
1 X z ðtÞ 2pjt p1ffiffiffiffi ¼ sin wðtÞ ¼ sj ðEÞsin TE 2E j¼1;3;… 1 X z ðtÞ 2pjt p2ffiffiffiffi ¼ cos wðtÞ ¼ cj ðEÞcos TE 2E j¼1;3;…
ð5Þ
xmin and xmax are the minimum and maximum displacement values for free undamped response at energy level E:
1 X k_x2 l ¼ 2kcos2 wl ¼ cj ðEÞ2 ¼ cðEÞ E j¼1;3;…
2E ¼ 2GðxÞ þ x_ 2 ¼ z21 þ z22 ¼ 2E sin2 w þ 2E cos2 w
ð6Þ
Free pffiffiffiffi undamped vibration thus describes a circle with radius 2E in the modified phase plane. z1 and z2 are given by pffiffiffiffiffiffiffi pffiffiffiffi z1 ¼ signðxÞ 2GðxÞ ¼ 2E sin w ð7Þ pffiffiffiffi z2 ¼ x_ ¼ 2E cos w with wt¼0 ¼ 0: The roles of sine and cosine are interchanged relative to the formulation by Krenk and Roberts [1]. In the case of free undamped vibration at energy level E; z1 and z2
ð9Þ
where k l indicates the mean value over a period of oscillation. It is observed that sðEÞ ¼ cðEÞ ¼ 2: In the case of a linear system s1 ðEÞ ¼ c1 ðEÞ ¼ 1; and all other coefficients vanish. 2.3. Spectral density The spectral density is defined as the Fourier transform of the covariance function. In the present case the one-sided spectral density is considered. It is defined as Sx ðvÞ ¼
1 ð1 E½xðtÞxðt þ tÞcosðvtÞdt p 21
ð10Þ
where E½ is the mean value. An approximation to the spectral density at the energy level l can be obtained as the sum of the spectral densities of each of the harmonics as
gcðEÞE p
Sx ðvlEÞ ¼
1 X j¼1;3;…
2cj ðEÞ2 ; ðj2 v2E 2 v2 Þ2 þ ðgcðEÞvÞ2
2.2. Modified state space variables In the case of a free undamped response, the energy level is constant. A set of modified phase plane variables ðz1 ; z2 Þ are introduced as
ð8Þ
since both variables are periodic with period TE : The coefficients of the sine and cosine expansions are related to the potential and kinetic energy. Taking the mean value of the square of the expressions (8) the following relations are obtained, 1 X k2GðxÞl ¼ 2ksin2 wl ¼ sj ðEÞ2 ¼ sðEÞ E j¼1;3;…
2.1. Probability density
TE ¼ 2
217
vE ¼
ð11Þ
2p TE
where cj ðEÞ are the Fourier coefficients introduced in Eq. (8). A more detailed derivation of Eq. (11) is given by Krenk and Roberts [1] by means of an approximate expression for the covariance function. The total spectrum is obtained by integrating with respect to the energy, weighting the spectrum at each energy level by the probability density of the energy, Sx ðvÞ ¼
ð1 0
Sx ðvlEÞpE ðEÞdE
ð12Þ
The method has also been extended to include systems with parametric excitation, Krenk et al. [18].
F. Ru¨dinger, S. Krenk / Probabilistic Engineering Mechanics 18 (2003) 215–222
218
3. Oscillator with bilinear stiffness A system with bilinear stiffness is considered. The stiffness function gðxÞ is shown in Fig. 1. Mathematically the function is described the following way, ( gðxÞ ¼
v21 x; v22 x
þ
aðv21
2
v22 ÞsignðxÞ;
lxl # a
ð13Þ
lxl . a
The initial stiffness is given by v21 : When the absolute value of the displacement exceeds a; the stiffness changes to v22 : It is assumed that both v1 ; v2 and a are larger than zero. For v1 ¼ v2 ; a ! 0 or a ! 1 the linear system is retrieved. The potential energy is obtained by integration of Eq. (13) as
GðxÞ ¼
8 < 12 v21 x2 ;
lxl # a
: 1 v2 x2 þaðv2 2 v2 Þlxlþ 1 a2 ðv2 2 v2 Þ; lxl . a 1 2 2 1 2 2 2
The potential energy is thus described by two parabolas with continuity of the zeroth and first derivative at the intersection. Since gðxÞ is an odd function, GðxÞ is even. 3.1. Probability density For free undamped vibration at energy level E; xmax is obtained by Gðxmax Þ ¼ E as E # 12 v21 a2 xmax ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > a > 2E 2 > 2 1 2 2 > : ð12 V Þþ V V 21þ a2 v2 ; E . 2 v1 a
ð15Þ
1
where the parameter V ¼ v1 =v2 will be referred to as the frequency ratio. Due to the symmetry in the potential
ð16Þ
1
The probability density of the energy is now determined directly from Eq. (4), and is conveniently investigated considering a non-dimensional formulation. Non-dimensional variables are introduced as
l¼ ð14Þ
8 pffiffiffiffi > 2E > > ; > < v1 a
energy, the natural period is determined as 4 ðxmax dx pffiffiffiffiffiffiffiffiffiffiffi TE ¼ pffiffi E2GðxÞ 2 0 8 2p > > ; E# 12 v21 a2 > > v > 1 > > > > > v1 a > < 4 sin21 p ffiffiffiffi 2E ¼ v1 > > 0 1 > > > > > V C > þ 4 cos21 B > @ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A; E. 12 v21 a2 > > v : 2 V2 21þ2E=a2 v2
E gE gv2 a2 ; E0 ¼ 12 v21 a2 ; b ¼ 0 ¼ 1 E0 pS0 2pS0
ð17Þ
E0 is thus a reference energy level and corresponds to the transition from the region with stiffness v21 to the region with stiffness v22 : For free undamped vibration the condition l #1 yields a harmonic vibration. b is a non-dimensional measure of the reference energy level relative to the input intensity. A non-dimensional natural period is defined as
vT T~ l ¼ 1 E 2p 8 1; l #1 > < ! ¼ 2 21 1 2V 21 V > pffiffi þ pffiffiffiffiffiffiffiffiffiffiffiffiffi ; l .1 cos : sin 2 p p l V 21þ l ð18Þ The probability density of the non-dimensional energy variable l is obtained from Eq. (4) as pl ðlÞ¼E0 pE ðE=E0 Þ¼C T~ l expð2blÞ
ð19Þ
It is thus seen that the probability distribution of l only depends on the damping parameter b and on the frequency ratio V: In order to verify the validity of the results comparison is made with records obtained by stochastic simulation. In Fig. 2 the probability density of l is shown for V ¼ 0:2 and b ¼ 0:2: The solid line represents the theoretical solution. The dots correspond to histogram points obtained by processing a simulated record. The results are seen to agree very well. 3.2. Free undamped vibration
Fig. 1. Bilinear stiffness function.
Free undamped vibration at energy level E is now considered. Only the range 0 # t # ð1=4ÞTE needs to be investigated due to symmetry. For E # ð1=2Þv21 a2 the system is linear and the phase plane variables ðx; x_ Þ are
F. Ru¨dinger, S. Krenk / Probabilistic Engineering Mechanics 18 (2003) 215–222
219
where t2 is the remaining part of the quarter period after the crossing of the level x ¼ a: t2 is determined from Eq. (16). Evaluation of A2 from Eq. (24a) gives the following value, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð25Þ A2 ¼ V a2 V2 2 a2 þ A21 The same result is obtained if continuity in x_ ðtÞ at time t ¼ t1 is used to determine this constant. The modified phase plane variables are determined from Eq. (6) as rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 2 2 2 z1 ðtÞ ¼ 2E 2 v2 A2 sin v2 4 TE 2 t ; ð26Þ 1 z2 ðtÞ ¼ v2 A2 sin v2 4 TE 2 t
Fig. 2. Probability density of non-dimensional energy, V ¼ 0:2; b ¼ 0:2:
simply given by xðtÞ ¼ A1 sinðv1 tÞ; pffiffiffiffi 2E A1 ¼ v1
x_ ðtÞ ¼ v1 A1 cosðv1 tÞ; ð20Þ
where A1 is the amplitude. For E . ð1=2Þv21 a2 the solution above is still valid until the time t1 where xðtÞ crosses into the range where the stiffness changes from v21 to v22 ; xðtÞ ¼ A1 sinðv1 tÞ;
x_ ðtÞ ¼ v1 A1 cosðv1 tÞ; 1 v1 a sin21 pffiffiffiffi 0 # t # t1 ¼ v1 2E
ð21Þ
For t . t1 the solution in terms of xðtÞ can be obtained by solving the following linear differential equation, x€ þ v22 x ¼ aðv22 2 v21 Þ
ð22Þ
The equation is non-homogeneous due to the non-zero equilibrium condition in the case where v1 – v2 : The initial conditions are given by Eq. (21) at t ¼ t1 since both xðtÞ and x_ ðtÞ are continuous. The solution is given by xðtÞ ¼ að1 2 V2 Þ þ A2 cos v2 14 TE 2 t ; ð23Þ x_ ðtÞ ¼ v2 A2 sin v2 14 TE 2 t ; t1 , t # 14 TE where only one arbitrary constant A2 appears, since the solution must have a form where the velocity is zero at t ¼ ð1=4ÞTE : The constant A2 is the amplitude of the harmonic part of the motion. Continuity in xðtÞ at t ¼ t1 yields aV2 ¼ A2 cosðv2 t2 Þ;
0 1 1 1 V B C t2 ¼ TE 2 t1 ¼ cos21 @ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi A 4 v2 2 2 2 V 2 1þ 2E=a v 1
ð24Þ
Again, the solution is most easily investigated by introducing a non-dimensional formulation. The time t and z1 and z2 are made non-dimensional by the following rescaling, z ðtÞ 1 ~t ¼ v1 t; ~t1 ¼ v1 t1 ¼ sin21 pffiffi ; z~ j ðtÞ ¼ pjffiffiffiffi ð27Þ 2E l where the non-dimensional energy variable l is introduced in Eq. (17). z~ 1 ð~tÞ and z~2 ð~tÞ are expressed as 8 l #1_0# ~t # ~t1 sinð~tÞ; > > < vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! 1 ~ 2 z~1 ð~tÞ¼ u u > t12 V 21þ l sin2 2 pTl 2 ~t ; otherwise > : l V ð28Þ 8 cosð~tÞ; > > < sffiffiffiffiffiffiffiffiffiffiffiffiffi z~2 ð~tÞ¼ V2 21þ l > > sin : l
l #1_0# ~t # ~t1 1 ~ ~ 2 pTl 2 t
V
! ;
otherwise
Finally, the modified phase plane variables are expanded in a Fourier series. The Fourier coefficients sj ðlÞ and cj ðlÞ introduced in Eq. (8) can be evaluated by the following non-dimensional integrals, 4 ðð1=2ÞpT~ l j~t sj ðlÞ¼ z~1 ð~tÞsin d~t; ð29Þ ~ ~ p Tl 0 Tl 4 ðð1=2ÞpT~ l j~t z~2 ð~tÞcos cj ðlÞ¼ d~t ~ ~ p Tl 0 Tl where the non-dimensional period T~ l is defined in Eq. (18). Due to symmetry the integration can be reduced to the first quarter of the period and the result multiplied by 4. The above formulation shows that the coefficients sj ðlÞ and cj ðlÞ only depend on l; V and j: In the following examples the coefficients will be obtained by numerical integration. The Fourier series expansions are now compared to the exact solutions in order to estimate the number of terms necessary for an accurate prediction of the response. In Figs. 3 and 4 the non-dimensional phase plane variables z~1 and z~2 are shown as function of the non-dimensional time 4t=TE : The solid line is the analytical solution. The first
220
F. Ru¨dinger, S. Krenk / Probabilistic Engineering Mechanics 18 (2003) 215–222
conveniently expressed in the following non-dimensional formulation, Sx ðrÞv31 g~ ð1 ¼ cðlÞlpl ðlÞ E0 p 0 1 X 2cj ðlÞ2 £ dl; 2 2 2 ~ 22 2 2 2 j¼1;3;… ðj Tl 2r Þ þ g~ cðlÞ r g v g~ ¼ ; r ¼ ð30Þ v1 v1
Fig. 3. Modified state space variables, l ¼ 10; V ¼ 10; (—) exact, ( – –) 1 term, ( – · –) 3 terms.
modified phase plane variable z~1 starts at 0, while the second z~2 starts at one. In Fig. 3 the parameters are l ¼ 10 and V ¼ 10: The solution for an expansion with one term is shown by the dashed line and the solution for an expansion with three terms (j ¼ 1; 3; 5) is shown by the dashed/dotted line. Inclusion of only one term in the series expansion gives a rather inaccurate representation, while inclusion of three terms leads to an accurate representation of the exact function. In Fig. 4 the functions are shown for the parameters l ¼ 10 and V ¼ 0:2 and expansions with one term and three terms have been included in this case as well. The accuracy of the truncated series expansions is of the same order of magnitude as in the previous case. 3.3. Spectral density The spectral density is obtained from the expression (11) and (12). In the present case the spectral density is most
Fig. 4. Modified state space variables, l ¼ 10; V ¼ 0:2; (—) exact, (– –) 1 term, ( – · –) 3 terms.
In the following numerical examples the integration in this equation is performed numerically, so the equation takes the form of a double summation. From the above expression and from the earlier expressions for pl ðlÞ; cj ðlÞ and T~ l it is seen that the spectral density in this non-dimensional form is a function of V; b and g~ only. In Figs. 5 – 9 a number of numerical examples are shown for different combinations of the parameters V; b and g~: The dots indicate solutions obtained by stochastic simulation. In order to demonstrate the effect of the truncation of the series in (40) both the solution with one term and the solution with three terms (j ¼ 1; 3; 5) are included, indicated by a dashed line and a solid line, respectively. In Fig. 9 the solution with two terms (j ¼ 1; 3) is also included indicated by the dashed/dotted line. In Fig. 5 the values are chosen as V ¼ 0:6; b ¼ 1 and g~ ¼ 0:05: The degree of non-linearity in this case is relatively small, as seen by the magnitude of the peak of the higher harmonic. It is clearly seen that inclusion of only one term in the solution fails to capture the second peak, and that the additional terms have no influence (are negligibly small) near the fundamental frequency. The analytical results are seen to agree very well with the results obtained by stochastic simulation, except for the fact that the theory predicts a peak at v ¼ v1 : This corresponds to oscillations taking place at small energy levels where the response does not enter into the region of changed stiffness v22 : However, this peak is not as pronounced in the simulated
Fig. 5. Spectral density, V ¼ 0:6; b ¼ 1; g~ ¼ 0:05: (X) Stochastic simulation, ( – –) theory (1 term), (—) theory (3 terms).
F. Ru¨dinger, S. Krenk / Probabilistic Engineering Mechanics 18 (2003) 215–222
221
Fig. 6. Spectral density, V ¼ 0:2; b ¼ 1; g~ ¼ 0:05: (X) Stochastic simulation, (– –) theory (1 term), (—) theory (3 terms).
Fig. 8. Spectral density, V ¼ 5; b ¼ 0:5; g~ ¼ 0:05: (X) Stochastic simulation, ( – –) theory (1 term), (—) theory (3 terms).
response. This observation was also made by Fogli and Bressolette [8]. In Fig. 6 the level of non-linearity is increased by decreasing the frequency ratio to V ¼ 0:3; and letting the other parameters remain unchanged. The bandwidth of the fundamental frequency is now larger (due to the frequency sweep), and both the magnitude and bandwidth of the higher harmonics are increased. Again, the pronounced peak predicted by the theoretical result does not show in the spectrum of the simulated record. Furthermore, a small deviation is seen at zero frequency, a feature which is also observable in results given by Fogli and Bressolette [8], Krenk and Roberts [1]. A discrepancy reported by Krenk and Roberts [1] at high frequencies could be due to truncation after the first two terms. In the above examples the system is characterized by hardening stiffness, i.e. v2 . v1 or V , 1: In Figs. 7 and 8 systems with softening stiffness are considered. In Fig. 7 the values are chosen as V ¼ 2; b ¼ 0:5 and g~ ¼ 0:05:
As opposed to the two previous cases, the peak is now located below v ¼ v1 ; due to the fact that v2 , v1 : The non-linearity is relatively small, as seen by the magnitude of the higher harmonic. Again a small deviation is seen at zero frequency. In Fig. 8 the non-linearity is increased by increasing the value of the frequency ratio to V ¼ 5; while the other parameters remain unchanged. This is seen to broaden both the fundamental frequency peak and the peak of the first higher harmonic. The discrepancy at zero frequency is larger in this case compared to the previous, and the accuracy a bit lower, though still satisfactory. The peak at v ¼ v1 observed in Figs. 5 and 6 in the theoretical solution is less pronounced for the results shown in Figs. 7 and 8. Finally the parameter combination V ¼ 0:3; b ¼ 1 and g~ ¼ 0:01 is investigated. This relatively strong non-linearity in combination with a decreasing value of the damping parameter g~ results in a spectrum, where it is possible to distinguish both the second and third higher harmonic.
Fig. 7. Spectral density, V ¼ 2; b ¼ 0:5; g~ ¼ 0:05: (X) Stochastic simulation, (– –) theory (1 term), (—) theory (3 terms).
Fig. 9. Spectral density, V ¼ 0:3; b ¼ 1; g~ ¼ 0:01: (X) Stochastic simulation, ( – –) theory (1 term), ( – · –) theory (2 terms), (—) theory (3 terms).
222
F. Ru¨dinger, S. Krenk / Probabilistic Engineering Mechanics 18 (2003) 215–222
In Fig. 9 both the solutions with two terms and three terms have been included, and the effect of including the third term on the third harmonic is quite evident. From this plot it seems that the deviations observed by Krenk and Roberts [1] at high frequencies are due to neglecting the third term.
4. Conclusion An exact solution for the probability density of the energy and an approximate analytical expression for the power spectral density of an oscillator with bilinear stiffness and Gaussian white noise excitation has been derived following a method proposed by Krenk and Roberts [1]. The result has been presented in a convenient non-dimensional form, and the accuracy of the theory has been investigated by comparing with results obtained from simulated records. The approximation of the free undamped response by a series expansion is shown to be very accurate. Even in the case of a strong non-linearity in the stiffness only a few terms are needed to represent the response. This accuracy is reflected in the approximation of the power spectral density, where each term in the series expansion corresponds to including one additional resonance peak. The examples clearly show the presence of higher harmonics and the broadening of the resonance peak and demonstrate the ability of the theory to capture these characteristics.
Acknowledgements The project is supported by the Danish Technical Research Council. Part of the work was carried out while the first author was a Marie Curie grant holder at the European Laboratory for Structural Assessment, Joint Research Centre of the European Commission at Ispra, Italy. This support is gratefully acknowledged.
References [1] Krenk S, Roberts JB. Local similarity in non-linear random vibration. J Appl Mech 1999;66:225 –35. [2] Caughey TK. Nonlinear theory of random vibrations. Advances in applied mechanics. New York: Academic Press; 1971. p. 209–53. [3] Lin YK, Cai GQ. Probabilistic structural dynamics. New York: McGraw-Hill; 1995. [4] Dimentberg MF, Hou Z, Noori M. Spectral density of a non-linear single-degree-of-freedom system’s response to a white-noise random excitation: a unique case of an exact solution. Int J Non-linear Mech 1995;30:673–6. [5] Dimentberg MF, Lin YK. An exact solution for response spectral density of a single-degree-of-freedom system under both parametric and additive white noise excitation. J Appl Mech 2002;69:399–400. [6] Roberts JB, Spanos PD. Random vibration and statistical linearization. Chichester: Wiley; 1990. [7] Bouc R. The power spectral density of response for a strongly non-linear random oscillator. J Sound Vibration 1994;175:317– 31. [8] Fogli M, Bressolette PH. Spectral response of a stochastic oscillator under impacts. Meccanica 1997;32:1–12. [9] Atkinson JD. Eigenfunction expansions for randomly excited nonlinear systems. J Sound Vibration 1973;30:153 –72. [10] Cai GQ, Lin YK. Response spectral densities of strongly nonlinear systems under random excitation. Probabilistic Engng Mech 1997;12: 41– 7. [11] Roy RV, Spanos PD. Wiener-hermite functional representation of nonlinear stochastic systems. Struct Safety 1989;6:187–202. [12] Roy RV, Spanos PD. Power spectral density of nonlinear system response: the recursion method. J Appl Mech 1993;60:358– 65. [13] Miles RN. An approximate solution for the spectral response of Duffing’s oscillator with random input. J Sound Vibration 1989;132: 43– 9. [14] Miles RN. Spectral response of a bilinear oscillator. J Sound Vibration 1993;163:319–26. [15] Fogli M, Bressolette P, Bernard P. The dynamics of a stochastic oscillator with impacts. Eur J Mech A: Solids 1996;15:213–41. [16] Soize C. Stochastic linearization method with random parameters for SDOF nonlinear dynamical systems: prediction and identification procedures. Probabilistic Engng Mech 1995;10:143–52. [17] Bellizzi S, Bouc R. Spectral response of asymmetrical random oscillators. Probabilistic Engng Mech 1996;11:51–9. [18] Krenk S, Lin YK, Ru¨dinger F. Effective system properties and spectral density in random vibration. J Appl Mech 2002;69:161 –70.