Iterative techniques for analyzing nonlinear vibrating dynamical systems in the frequency domain

Abstract

The nonlinear response of prototypical structures experiencing harmonic excitation is studied using novel techniques called iterative harmonic analysis (IHA) and iterative modal analysis (IMA). First, a simple damped oscillator with a cubic hardening stiffness nonlinearity is studied, and IHA is used in this single-degree-of-freedom system. In this first section, a high-order harmonic balance is applied, and IHA is employed in order to find the amplitude coefficients for different harmonics and their codependence. Additionally, a set of nested sums are identified that describe the harmonic coupling explicitly. Secondly, a pinned–pinned nonlinear beam of rectangular cross-section is studied, and IMA is applied to find the amplitude coefficients for different modes and their codependence. The nonlinearity is introduced through the membrane effect, where axial strain due to transverse deflection becomes a significant contributor to the system behavior. Typical frequency-domain methods cannot be easily applied to these systems as the solutions of the differential equations lead to intricate coupling between coefficients of the solution, and no analytical expression exists for those coefficients. Based upon these examples, other nonlinear systems may also be considered in future work using either a modal-based or finite element model. Finally, the advantages of the new method (reduced computational cost) as well as the limitations (effectively the same as those of an nth-order harmonic balance) are emphasized.

Appendices

Appendix A: Formulation for base excitation

The partial differential equation changes form when excited by base excitation, as the time-dependent input is now in the boundary conditions. This may still be converted into an ordinary differential equation for each modal response.

$$\begin{aligned}&EI\frac{\partial ^4w}{\partial x^4}-N_x\frac{\partial ^2w}{\partial x^2}+\rho A\frac{\partial ^2w}{\partial t^2}=0 \end{aligned}$$

(68)

$$\begin{aligned}&N_x=\frac{EA}{2L}\int _{0}^{l_x}\left( \frac{\partial w}{\partial x}\right) ^2\mathrm{d}x \end{aligned}$$

(69)

subject to

$$\begin{aligned} \begin{aligned}&w\left( t,0\right) =w\left( t,l_x\right) =H\sin \omega t \\&\frac{\partial ^2w\left( t,0\right) }{\partial x^2}=\frac{\partial ^2w\left( t,l_x\right) }{\partial x^2} = 0 \end{aligned} \end{aligned}$$

(70)

Now, projecting onto the adjoint eigenfunction yields boundary terms. This projection

$$\begin{aligned} \left\langle EI\frac{\partial ^4w}{\partial x^4}-N_x\frac{\partial ^2w}{\partial x^2}+\rho A\frac{\partial ^2w}{\partial t^2}=0,\sin \frac{n\pi x}{l_x} \right\rangle \end{aligned}$$

(71)

generates the following equality

$$\begin{aligned}&\left\langle EI\frac{\partial ^4w}{\partial x^4},\sin \frac{n\pi x}{l_x}\right\rangle =-EIH\left[ 1-\left( -1\right) ^n\right] \left( \frac{n\pi }{l}\right) ^3\nonumber \\&\quad \sin \omega t+EI\left( \frac{n\pi }{l}\right) ^4\left| \phi _n^2\right| a_n \end{aligned}$$

(72)

Including the nonlinear analysis,

$$\begin{aligned}&\rho A\ddot{a}_n+b\dot{a}_n+EI\left( \frac{n\pi }{l_x}\right) ^4a_n\left[ 1+\frac{1}{4n^2}\frac{A}{I}\sum _{p}p^2a_{p}^2\right] \nonumber \\&\quad =\frac{2EIH}{l_x}\left[ 1-\left( -1\right) ^n\right] \left( \frac{n\pi }{l}\right) ^3\sin \omega t \end{aligned}$$

(73)

or

$$\begin{aligned}&M_n\left( \ddot{a}_n+2\left( \zeta \omega \right) _n\dot{a}_n+\omega _n^2a_n\left[ 1+\alpha _n\sum _{p}p^2a_{p}^2\right] \right) \nonumber \\&\quad =F_n\left[ 1-\left( -1\right) ^n\right] \sin \omega t \end{aligned}$$

(74)

where

$$\begin{aligned} F_n\equiv EIH\left( \frac{n\pi }{l_x}\right) ^3 \end{aligned}$$

(75)

Fig. 20

RMS response of the linear beam model at query point $x_F$ for increasing forcing levels ranging from 1 to 100 N. $\zeta _1=0.02$

Appendix B: Iterative modal analysis for broadband signals

For a bounded, broadband signal,

$$\begin{aligned} \Phi _F=\left\{ \begin{array}{ll} 0 &{}\quad 0<\omega<\omega _\mathrm{min} \\ \frac{F_0^2}{2\Delta \omega } &{}\quad \omega _\mathrm{min}<\omega<\omega _\mathrm{max} \\ 0 &{}\quad \omega _\mathrm{max}<\omega <\infty \end{array}\right. \end{aligned}$$

(76)

Dowell and Kubota showed that the mean-square response of a mode in a system excited by such a signal is [4]

$$\begin{aligned} \bar{A}_n^2=\int _{0}^{\infty }\left| H_{A_n}\right| ^2\Phi _F\mathrm{d}\omega =\frac{\pi \psi _n^2\left( x_F\right) F_0^2}{8M_n^2\Delta \omega \left( \zeta \omega \right) _n\omega _n^2}. \end{aligned}$$

(77)

However, outside of the high-frequency regime, the bounds of the integral cannot be assumed to be large for a bounded broadband signal. The value of the mean-square response of a given mode is found in Appendix C, showing that

$$\begin{aligned} \bar{A}_n^2= & {} \frac{\psi _{nF}^2F_0^2}{8M_n^2\Delta \omega }\frac{1}{\left( \zeta \omega \right) _n\omega _n^2\sqrt{\zeta _n^2-1}}\nonumber \\&\times \left[ \arctan \left( \frac{\omega _1}{\omega _n}s_-\right) s_--\arctan \left( \frac{\omega _1}{\omega _n}s_+\right) s_+\right. \nonumber \\&\left. -\arctan \left( \frac{\omega _2}{\omega _n}s_+\right) s_++\arctan \left( \frac{\omega _2}{\omega _n}s_-\right) s_-\right] \nonumber \\ \end{aligned}$$

(78)

where

$$\begin{aligned} s_+= & {} \sqrt{2\zeta _n^2-1+2\zeta _n\sqrt{\zeta _n^2-1}} = e^{-i\frac{1}{2}\arctan \frac{2\zeta _n\sqrt{1-\zeta _n^2}}{1-2\zeta _n^2}} \nonumber \\ s_-= & {} \sqrt{2\zeta _n^2-1-2\zeta _n\sqrt{\zeta _n^2-1}} = e^{i\frac{1}{2}\arctan \frac{2\zeta _n\sqrt{1-\zeta _n^2}}{1-2\zeta _n^2}} \nonumber \\ s_+= & {} \frac{1}{s_-} \end{aligned}$$

(79)

Its approximate peak value is therefore

$$\begin{aligned} A_n^2= & {} \frac{\psi _{nF}^2F_0^2}{4M_n^2\Delta \omega }\frac{1}{\left( \zeta \omega \right) _n\omega _n^2\sqrt{\zeta _n^2-1}}\nonumber \\&\times \left[ \arctan \left( \frac{\omega _1}{\omega _n}s_-\right) s_--\arctan \left( \frac{\omega _1}{\omega _n}s_+\right) s_+\right. \nonumber \\&\left. -\arctan \left( \frac{\omega _2}{\omega _n}s_+\right) s_++\arctan \left( \frac{\omega _2}{\omega _n}s_-\right) s_-\right] \nonumber \\ \end{aligned}$$

(80)

where $\omega _m$ is now the natural frequencies between $\omega _\mathrm{min}$ and $\omega _\mathrm{max}$. Using the same update relation for $\omega _{n,j}$

$$\begin{aligned} \omega _{n,j}^2= & {} \omega _{n,j-1}^2+\left( 1-\gamma \right) \nonumber \\&\times \left( \omega _{n,0}^2\left[ 1+\alpha _n\sum _pp^2A_{p,j-1}^2\right] -\omega _{n,j-1}^2\right) \nonumber \\ \end{aligned}$$

(81)

The vector of these updated natural frequencies are truncated to include only the subset that lies within $\omega _\mathrm{min}<\omega _{n,j}<\omega _\mathrm{max}$, and the method proceeds. Time marching will be set up in the same way as in the work of Culver and Dowell [2], using a linearly swept frequency input from $\omega _\mathrm{min}$ to $\omega _\mathrm{max}$.

1.1 Results

The behavior of the system subject to a bounded broadband excitation is reported in the following figures. The independent axis now represents the center frequency of a 1/3 octave band normalized by the value of the fundamental natural frequency of the system.

Fig. 21

RMS response of the nonlinear beam model to a 1/3 octave input at query point $x_F$ for increasing forcing levels ranging from 1 to 100 N. $\zeta _1=0.02$

Fig. 22

RMS response of the fundamental mode of the nonlinear beam model to a 1/3 octave input at query point $x_F$ for increasing forcing levels ranging from 1 to 100 N. $\zeta _1=0.02$

First, we compare the estimates of IMA and time marching for the linear beam to verify the computational test foundations. Figure 20 shows a comparison of time marching to IMA for the linear case ($\alpha _n=0$). This reduces to classical modal analysis. The only disparities are either a consequence of the transient homogeneous response partially remaining even after a 95% settling time (in the very low response magnitude regions) or the modeling of a broadband signal with a frequency sweep (small errors in the peak response regions).

Fig. 23

Runtimes for time marching and IMA for various force levels and excitation frequencies of the broadband forcing case

Now we may compare the results for the beam including the nonlinearity. Figures 21 and 22 show the comparison in the nonlinear regime. Figure 22 illustrates very interesting behavior of the system even for time marching as it exhibits chaotic behavior in the highest forcing levels as the system chooses one of two stable equilibria. The harmonic balance method and hence IMA is not able to capture this chaotic behavior of course.

The disagreement at the higher forcing levels, once again, is caused in part by a breakdown of the harmonic balance trial solution assumption. Figure 22 highlights regions in the responses to excitation bands near the fundamental, second, and third frequency of the beam. The agreement at these forcing levels is not as close as that of the system experiencing single-harmonic forcing, and suggests that IMA may not be suitable for low-frequency (below the fourth natural frequency) broadband analysis of beams. However, at higher frequencies, the method converges with time marching and is readily applicable.

As seen in the main body of this work, one of the main advantages of IMA is improved computational cost. The runtimes for the broadband forced nonlinear case are shown in Fig. 23. For the very high forcing levels’ peak responses, IMA could not converge on a solution for $\gamma =90\%$. For almost all converged solutions, IMA is less than 10% of the computational cost of time marching.

Appendix C: Broadband transfer function integral

The power spectrum for a bounded broadband input is

$$\begin{aligned} \Phi _F=\frac{F_0^2}{2\Delta \omega }H\left( \omega -\omega _1\right) \left[ 1-H\left( \omega -\omega _2\right) \right] \end{aligned}$$

(82)

So the integral for the mean-square response of a modal amplitude is

$$\begin{aligned} \bar{A}_n^2=\int _{0}^{\infty }\frac{\psi _{nF}^2F_0^2}{2M_n^2\Delta \omega }\frac{H\left( \omega -\omega _1\right) \left[ 1-H\left( \omega -\omega _2\right) \right] }{\left( \omega _n^2-\omega ^2\right) ^2+4\left( \zeta \omega \right) _n^2\omega ^2}\mathrm{d}\omega \nonumber \\ \end{aligned}$$

(83)

or

$$\begin{aligned} \bar{A}_n^2=\frac{\psi _{nF}^2F_0^2}{2M_n^2\Delta \omega }\int _{\omega _1}^{\omega _2}\frac{1}{\left( \omega _n^2-\omega ^2\right) ^2+4\left( \zeta \omega \right) _n^2\omega ^2}\mathrm{d}\omega \nonumber \\ \end{aligned}$$

(84)

Fig. 24

a (Upper left) the response history (time is normalized by the forcing period and response is normalized by the beam thickness) indicates linear behavior at approximately 8% of the beam thickness. b (Upper right) the linear Poincare map converges to a point. c (Lower left) a simple helical time-phase history is indicative of linear behavior. d (Lower right) the velocity response (normalized by $\frac{2h\omega _f}{\pi }$) also show simple harmonic behavior

The value of Eq. 83 is no longer capable of resolving Eq. 84 may be found by partial fraction expansion into the sum of two expressions whose form is that of the arctangent integrand Therefore, the value of this integral is

$$\begin{aligned}&\frac{1}{4}\frac{1}{\sqrt{\omega _n^4\zeta _n^2\left( \zeta _n^2-1\right) }}\frac{1}{\sqrt{2\zeta _n^2\omega _n^2-\omega _n^2+2\sqrt{\omega _n^4\zeta _n^2\left( \zeta _n^2-1\right) }}}\nonumber \\&\qquad \frac{1}{\sqrt{2\zeta _n^2\omega _n^2-\omega _n^2-2\sqrt{\omega _n^4\zeta _n^2\left( \zeta _n^2-1\right) }}}\ldots \nonumber \\&\quad \left[ \arctan \left( \frac{\omega _1}{\sqrt{2\zeta _n^2\omega _n^2-\omega _n^2+2\sqrt{\omega _n^4\zeta _n^2\left( \zeta _n^2-1\right) }}}\right) \right. \nonumber \\&\qquad \sqrt{2\zeta _n^2\omega _n^2-\omega _n^2-2\sqrt{\omega _n^4\zeta _n^2\left( \zeta _n^2-1\right) }}\nonumber \\&\quad -\arctan \left( \frac{\omega _1}{\sqrt{2\zeta _n^2\omega _n^2-\omega _n^2-2\sqrt{\omega _n^4\zeta _n^2\left( \zeta _n^2-1\right) }}}\right) \nonumber \\&\quad \sqrt{2\zeta _n^2\omega _n^2-\omega _n^2+2\sqrt{\omega _n^4\zeta _n^2\left( \zeta _n^2-1\right) }}\\&\quad -\arctan \left( \frac{\omega _2}{\sqrt{2\zeta _n^2\omega _n^2-\omega _n^2+2\sqrt{\omega _n^4\zeta _n^2\left( \zeta _n^2-1\right) }}}\right) \nonumber \\&\quad \sqrt{2\zeta _n^2\omega _n^2-\omega _n^2-2\sqrt{\omega _n^4\zeta _n^2\left( \zeta _n^2-1\right) }} \nonumber \\&\quad +\arctan \left( \frac{\omega _2}{\sqrt{2\zeta _n^2\omega _n^2-\omega _n^2-2\sqrt{\omega _n^4\zeta _n^2\left( \zeta _n^2-1\right) }}}\right) \nonumber \\&\quad \left. \sqrt{2\zeta _n^2\omega _n^2-\omega _n^2+2\sqrt{\omega _n^4\zeta _n^2\left( \zeta _n^2-1\right) }}\right] \nonumber \end{aligned}$$

(85)

This can be simplified recognizing that the variables are, themselves, real.

$$\begin{aligned} \begin{aligned} \bar{A}_n^2&=\frac{\psi _{nF}^2F_0^2}{8M_n^2\Delta \omega }\frac{1}{\left( \zeta \omega \right) _n\omega _n^2\sqrt{\zeta _n^2-1}}\\&\quad \left[ \arctan \left( \frac{\omega _1}{\omega _n}s_-\right) s_--\arctan \left( \frac{\omega _1}{\omega _n}s_+\right) s_+\right. \\&\qquad \left. -\arctan \left( \frac{\omega _2}{\omega _n}s_+\right) s_++\arctan \left( \frac{\omega _2}{\omega _n}s_-\right) s_-\right] \end{aligned} \end{aligned}$$

(86)

where

$$\begin{aligned} s_+&= \sqrt{2\zeta _n^2-1+2\zeta _n\sqrt{\zeta _n^2-1}} = e^{-i\frac{1}{2}\arctan \frac{2\zeta _n\sqrt{1-\zeta _n^2}}{1-2\zeta _n^2}}\nonumber \\ s_-&= \sqrt{2\zeta _n^2-1-2\zeta _n\sqrt{\zeta _n^2-1}} = e^{i\frac{1}{2}\arctan \frac{2\zeta _n\sqrt{1-\zeta _n^2}}{1-2\zeta _n^2}}\nonumber \\ s_+&= \frac{1}{s_-} \end{aligned}$$

(87)

Fig. 25

a (Upper left) the response history indicates quasi-linear behavior at approximately 35% of the beam thickness. b (Upper right) the Poincare map shows an egg-shaped attractor. c (Lower left) the time-phase history is no longer a simple helix in the last ten forcing periods of the simulation. d (Lower right) the velocity response (normalized by $\frac{2h\omega _f}{\pi }$) also shows quasi-linear behavior

Appendix D: Onset of chaos

Two assumptions have been made in the development of IMA. The first is that the beam slope $\frac{\partial w}{\partial x}$ is very small relative to 1. The second requires the modal amplitude time-dependent functions to be harmonic with the same frequency as the forcing in the steady solution. The validity of the former may be estimated with response amplitude and wavelength. The second is more complicated. Consider the following figures of response attributes for the prototypical beam excited by a 700 radian per second point forcing of specified magnitude.

Figure 25 shows the system behavior in an approximately linear limit. The response and response velocity histories are simple harmonic with the same frequency as the excitation, the Poincare map converges to a point, and the time-phase history is a simple helix (Fig. 24).

As the excitation amplitude increases, so does the response. As the response approaches one-sixth of the beam thickness (see Fig. 25), the beam modal amplitude responses begin to shift, but harmonic balance is still quite accurate. However, as the response approaches one-third of the beam thickness, the system’s Poincare map may have attractors with varying shapes.

Fig. 26

a (Upper left) Quasiperiodicity. b (Upper right) The Poincare map shows an egg-shaped attractor. c (Lower left) complicated helical time-phase response. d (Lower right) Quasiperiodicity in the velocity response

Fig. 27

a (Upper left) period-doubled simple periodicity. b (Upper right) The Poincare map shows two point attractors. c (Lower left) alternating elliptical helix behavior in the time-phase history. d (Lower right) period-doubling in the velocity response

Fig. 28

a (Upper left) period-doubled simple periodicity. b (Upper right) The Poincare map shows a single point attractor. c (Lower left) alternating elliptical helix behavior in the time-phase history. d (Lower right) period-doubling in the velocity response

Fig. 29

a (Upper left) Quasiperiodicity. b (Upper right) The Poincare map shows an strange attractor. c (Lower left) complicated helical time–phase response. d (Lower right) Quasiperiodicity in the velocity response

Fig. 30

a (Upper Left) Quasiperiodicity. b (Upper right) The Poincare map shows an egg-shaped attractor. c (Lower left) complicated helical time–phase response. d (Lower right) Quasiperiodicity in the velocity response

Responses whose amplitude ranges from 30% of the thickness to 60% of the thickness exhibit a number of characteristics. Some are simple periodic, others show period-doubling, and others still show quasiperiodicity. See Figs. 26, 27, 28, 29, 30 and 31.

Fig. 31

a (Upper left) the response history returns to simple periodic behavior at 60% of the beam thickness. b (Upper right) The Poincare map, again, converges to a point. c (Lower left) the time–phase history is a simple helix. d (Lower right) the velocity response is also simple harmonic

For increasing forcing amplitudes, the system response amplitude varies non-monotonically between 60 and 80% of the beam thickness for excitation forces between 500 and 650 N. In this range, simple harmonic, quasiperiodic, and even chaotic behavior appears. Consider Fig. 32a, d. Figure 31a shows simple periodic behavior 1600 forcing periods after the linear estimate for transient settling. However, Fig. 31d shows chaos 200 forcing periods after this settling time has passed. The chaos now disappears (see Fig. 32) well after the linear settling time.

Fig. 32

a (Upper left) the response history returns to simple periodic behavior at 60% of the beam thickness. b (Upper right) The Poincare map exhibits chaos, but eventually converges to a point. c (Lower left) the time-phase history at the end of the simulation is a simple helix. d (Lower right) the velocity response early in the estimated steady solution is chaotic

Near 650 N, the response amplitude suddenly jumps to 100% of the beam thickness, and the response is again periodic (though it exhibits higher harmonic behavior). See figure (Fig. 33)

Fig. 33

a (Upper left) the response history is periodic behavior with an amplitude that is 100% of the beam thickness. b ( Upper right) The Poincare map converges to a point. c (Lower left) the time-phase history at the end of the simulation is helical, with small return coils indicative of higher harmonic influences on the response. d (Lower right) the velocity response is also periodic, showing higher harmonic behavior

This periodicity appears to continue as the forcing amplitude increases, even up to 900 N where, again, chaos is only evident early in the expected steady response.

Fig. 34

a (Upper left) the response history is periodic behavior with an amplitude that is over 100% of the beam thickness. b (Upper right) The Poincare map converges to a point after capturing chaotic behavior early in the predicted steady response. c (Lower left) the time-phase history at the end of the simulation is helical, with a small return coils indicative of higher harmonic influences on the response. d (Lower right) the velocity response early in the expected steady response shows chaos

The onset of chaos early in the predicted steady response in Fig. 32 suggests that harmonic balance may no longer be valid when vibrations amplitude exceed 60% of the beam thickness (Fig. 34).