#### Abstract

This study addresses the chaotic phenomena and nonlinear responses in a vibroacoustic system. It is the first study about the chaotic phenomena in a vibroacoustic system, which is formed by a flexible panel coupled with a cavity. A multimode formulation is developed from the acoustic governing equation and nonlinear structural governing equation. The chaotic and various nonlinear responses are computed from the multimode formulation using a numerical integration method. The results obtained from the proposed method and classical harmonic balance method are generally consistent. A set of modal convergence studies is performed to check the proposed method. The effects of various parameters on triggering the nonchaotic responses to chaotic responses in a vibroacoustic system are studied in detail.

#### 1. Introduction

Over the past decades, numerous researchers have been working on research topics related to chaos science, nonlinear vibration, and vibroacoustics (e.g., [1–6]). In practice, linear designs for chaotic and nonlinear machines or structures are inappropriate and result in unsafety. In fact, there are many physical machines and structures, which would undergo chaotic and nonlinear vibrations. For example, Tian et al. [7] studied the nonlinear aeroelastic characteristics of a trapezoidal wing in hypersonic flow. In their numerical results, it was found that the geometrical parameters of trapezoidal wing imposed significant effects on the nonlinear aeroelastic behaviors of wing structure; and the evolution processes of chaos exhibited remarkable difference for the wing configurations considered in the study. Rao et al. [8] presented a work about the dynamics of a cracked rotor system with oil-film force in parameter space. The “eye” of chaos was found in the cracked rotor system, emerging as the accumulation limit of forward and reverse period-doubling bifurcation cascades. Asemani and Vatankhah [9] proposed a new control system to stabilize the unstable periodic orbit of chaotic spinning disks with incomplete state information. The proposed control structure was developed according to the T-S fuzzy systems and its design procedure fulfilled the constraint in the T-S fuzzy dynamic output feedback control signal. Akbarimajd and Yousefi [10] proposed a new control strategy based on Takagi-Sugeno fuzzy model for deceasing the power system oscillation. In the control system, the stability of the whole closed-loop model was enhanced using a general Lyapunov-Krasovskii functional. The proposed strategy was applied to a 16-machine/68-bus power system. In the nonlinear time domain simulations, the effectiveness of the proposed method was checked.

Moreover, so far, there are few research works about both nonlinear vibration and structural acoustics [11–15], although numerous studies about linear vibroacoustic (e.g., [16–19]) and nonlinear vibration (e.g., [20–24]) have been carried out. In [11, 12], Lee et al. did the research works about the sound radiation and absorption of a curved panel, which underwent nonlinear vibrations. In the simulation results, there were no chaotic phenomena observed. Lee et al. [13] studied the sound radiation of a chaotically vibrating curved beam/panel. In the theoretical model, there was no structural acoustic coupling term considered. In other words, the problem in [13] is not structural-acoustic and is different from the one in this study. In most of the linear structural-acoustic studies, various panel absorbers and panel-cavity systems were investigated. It was assumed that the structural vibrations in the systems were small and mainly focused on the sound absorption and sound reduction. For example, Lee et al. [16] studied the acoustic absorption of a finite flexible microperforated panel backed by an air cavity. The absorption formula for the microperforated absorber was based on the modal analysis solution of the classical plate equation coupled with the acoustic wave equation. Choy et al. [19] proposed a compact flow-through plate silencer using reinforced composite plates. The light-weight and high stiffness property was a crucial element in the silencer design. The other concept in the design was that the sound reflection from the plate of the silencer caused a desirable noise reduction from low to medium frequency with wide broadband. In practice, the structural parts in these panel absorbers and panel-cavity systems might be very thin and undergoing nonlinear vibration. In the studies of nonlinear vibration, various structures and systems (e.g., beam, plate, shell, and spring-mass) were investigated. Some of them focused on the solution methods. For example, Fan et al. [20] studied the steady-state periodic and quasi-periodic responses of van der Pol-Mathieu system. They proposed combining the method of multiple scales and double perturbation technique to obtain the special periodic and quasi-periodic solutions. Huang and Zhu [21] investigated the nonlinear dynamic responses of an Euler-Bernoulli beam attached to a rotating rigid hub with a constant angular velocity. They used Lagrange’s equations based on discretized expressions of kinetic and potential energies of the system to develop the spatially discretized governing equations. Then, they used the incremental harmonic balance method to solve the governing equations and obtain the results of periodic responses and period-doubling bifurcations.

To the best of the author’s knowledge, this study is the first one about the chaotic phenomena in a vibroacoustic system which considers the structural acoustic coupling. The multimode formulation is developed and solved by the numerical integration method. The effects of various parameters on triggering the nonchaotic responses to chaotic responses in a vibroacoustic system are studied in detail.

#### 2. Theory

Figure 1 shows a vibroacoustic system that is formed by a curve structure backed by a cavity. The governing equation of the acoustic pressure within the cavity is the well-known wave equation [26, 27]where is the acoustic pressure, is the sound speed, and represents the structural mode number.

Consider the modal decomposition approach and then is expressed in terms of where is the th acoustic mode shape and is the corresponding modal response; is the number of acoustic modes used. According to [26, 27], the mode shape function is taken to be those in an enclosure with rigid boundaries (i.e., , where , , and are the acoustic mode numbers).

Multiply to the right side of (1) and take integration over the cavity volume, :Consider the technique of integration by parts: where represents the normal direction of the cavity boundary surface and is the surface area.

Put (4) into (3): where ; ; at ; at ; is the air density; , , , , , and are the acoustic modes numbers and cavity lengths in the , , and directions, respectively; and are the excitation source displacement and curved panel displacement, respectively, which are expressed in the following forms: where and are the harmonic excitation source amplitude and curved panel amplitude; and are the 1st and th vibration mode shape functions of the source panel and curved panel, respectively. The boundary condition is simply supported (i.e., ).

Put (7a) and (7b) into (6): where ; is the wave number. The acceleration of the source panel is assumed as ; is the excitation frequency; is the dimensionless excitation parameter and is the gravity of 9.81 ms^{−2}.

From (2), the acoustic pressure is expressed in terms of . The acoustic pressure acting on the panel surface is given by According to (8), can be rewritten in the following form: In [11, 12], the “beam-like” curved panel was adopted and the flexural modes along the direction were ignored. It was experimentally found that the flexural modes along the direction were not very important in the nonlinear phenomena. According to [11–13], the governing equation for the nonlinear curved panel is expressed in the following form: where is the transverse displacement of the panel; is the initial deflection; is Young’s modulus; *ρ* is the material density; is the damping coefficient; is the width () and is the thickness; and is the acoustic pressure acting on the panel surface.

Consider the following modal decomposition: Put (12) into (11); multiply by each term on the right side and take integration over the surface. In this study, the first three structural modes are considered (i.e., ). Note that, in the convergence study in the next section, it is proven that the contribution of the 4th mode is very minimal. Therefore, the three modal equations are developed and given in the following: where ; ; ; ; = modal damping coefficient; , , and are the resonant frequencies of the 1st to 3rd modes; and are the structural mode numbers. According to [11, 26, 27], the acceleration terms on the right sides of (14a)–(14c) can be rewritten as .

The above coupled modal mode differential equations can be solved using the Runge-Kutta time domain numerical integration [11–13]. The overall root-mean square amplitude, positive amplitude, and negative amplitude of the displacement responses at the steady state are defined bywhere , , and are the modal displacement responses at the steady state at th time step; is the number of time steps used; and are the maximum and minimum values within the steady state, respectively.

#### 3. Results and Discussion

In this section, the material properties in the numerical cases considered are as follows: Young’s modulus = N/m^{2}, Poisson’s ratio = 0.3, and mass density = 2,700 kg/m^{3}. The panel dimensions are 0.5 m 0.4 m 2 mm. The air density is 1.2 kg/m^{3}. The sound speed is 340 m/s. In Tables 1(a)–1(c), the case is the nonlinear forced vibrations of a curved panel backed by a cavity. The first two symmetrical and antisymmetrical structural modes are used. The damping ratio . The excitation pressure is evenly distributed over the panel surface. In Table 1(a), the nonlinear vibration responses are nonchaotic. The initial centre deflection is 2 mm. The modal convergence study shows that the contributions of the two antisymmetrical modes are zero for various excitation frequencies. The approach of two symmetric modes is good enough. In Table 1(b), the nonlinear vibration responses are chaotic for all values. The excitation frequency is 0.805. The contribution of the first antisymmetrical mode is higher than that of the 2nd symmetric mode and cannot be neglected. It is noted that, for chaotic cases, the approach of two symmetric modes and one antisymmetric mode is necessary. In Table 1(c), the nonlinear vibration responses are also chaotic. The modal contributions are not very sensitive to the number of acoustic modes in the chaotic case. It is because the excitation frequency is far from the first nonzero cavity resonant frequency. Figures 2(a)–2(d) present the comparisons between the backbone curves and time histories of the nonlinear free panel vibrations obtained from the numerical integration method and classical harmonic balance method [25]. Figures 2(a)-2(b) show the 1st and 2nd symmetrical mode backbone curves, where in each of which the amplitude is plotted against the corresponding resonant frequency. The 1st four structural modes and 1st sixteen acoustic modes are used. and are the 1st mode linear and nonlinear resonant frequencies, respectively. is the initial centre deflection. Figures 2(c)-2(d) show the time histories for the vibration amplitude set as and 1.4. In the solution procedures of the classical harmonic balance method, the two harmonic terms (i.e., and ) are considered. The backbone curve and time history results obtained from the two methods are generally in good agreement.

**(a)**Comparison of the backbone curve results from the proposed and classical harmonic balance methods [25] (1st symmetric mode, mm, m, m, mm, no cavity)

**(b)**Comparison of the backbone curve results from the proposed and classical harmonic balance methods [25] (2nd symmetric mode, mm, m, m, mm, no cavity)

**(c)**Comparison of the time history free vibration results from the proposed and classical harmonic balance methods [25] (1st mode, mm, m, m, mm, no cavity, initial amplitude = )

**(d)**Comparison of the time history free vibration results from the proposed and classical harmonic balance methods [25] (1st mode, mm, m, m, mm, no cavity, initial amplitude = 1.4)Figure 3(a) shows the RMS vibration amplitude plotted against the excitation frequency for various excitation magnitudes. The vibration responses in these cases are nonchaotic as the excitation magnitude is not large enough and the curvature is not deep enough. The peak frequency and RMS vibration amplitude increase with the excitation magnitude. Due to the zero frequency cavity mode, the RMS vibration amplitude is higher when the excitation frequency is set closer to zero.

**(a)**Vibration amplitude versus excitation frequency for various source excitation magnitudes ( mm, m, m, m, , mm)

**(b)**Vibration amplitude versus excitation frequency for various curvatures ( mm, m, m, m, , )

**(c)**Vibration amplitude versus excitation frequency for various cavity depths ( mm, m, m, , mm, )In the case of small excitation (i.e., ), the peak frequency is around 2.8 and looks more linear than the other two. The peak frequency is higher than one because of the cavity stiffness. When the excitation magnitude is set higher (e.g., and 4), the peaks are inclined to the right side. It is called “hardening effect,” which implies that the structural stiffness is stronger due to the nonlinearity. In the case of small excitation, the solution line at the low frequency range is smoother than those in the two other cases, which contain some other peaks due to the nonlinearity. For an example, in the case of large excitation (i.e., ), there is an obvious peak around 1.5 which does not appear in the case of small excitation. The simple harmonic solution line is found in each of the three cases for the excitation frequency higher than 3. Figures 4(a)-4(b) show the steady-state time histories of the superharmonic and simple harmonic cases. Note that as only the vibration response at the steady state is shown, the time does not start at zero in each of these figures. In Figure 4(a), there are many peaks as the higher harmonic components are very significant. In Figure 4(b), the time history shows a set of simple sine waves. The numbers of sine waves and harmonic cycles are equal. That is why it is considered as “simple harmonic.”

**(a)**Time history for the superharmonic case ( mm, m, m, m, , mm, , )

**(b)**Time history for the simple harmonic case ( mm, m, m, m, , mm, , )

**(c)**Time history for the chaotic case ( mm, m, m, m, , mm, , )Figure 3(b) shows the RMS vibration amplitude plotted against the excitation frequency for various curvatures. In the case of 2.5 mm curvature, the “softening effect” is seen (i.e., the peaks are inclined to the left side). It is implied that the structural stiffness is weaker due to the nonlinearity. The “jump down” phenomena are observed for the excitation frequency decreasing from 4 to 3.5 and from 2.8 to 1.5, respectively. The peak value around 3 is smaller than the one in the case of 1 mm curvature; and the peak value around 1.5 is higher. There is an abrupt jump down around . At the very low frequency range (), the chaotic responses are observed in the case of 2.5 mm curvature (see Figure 4(c)); and the solution line is not smooth. In the time history, the equilibrium position, vibration amplitude, and vibration period are varying abruptly. That is why it is considered as “chaotic.” The deeper the panel curvature, the less smooth the solution line. It is because the curved panel would chaotically vibrate at the low frequency. Figure 3(c) shows the RMS vibration amplitude plotted against the excitation frequency for various cavity depths. The first peaks and second peaks on the solution lines are observed at the excitation frequencies around 1.35, 1.45, 1.65, 3, 3.5, and 4, respectively. The longer the cavity depth, the lower the 1st and 2nd peak frequencies, lower 2nd peak amplitude, and higher 1st peak amplitude. At the very low frequency range (), the three solution lines almost overlap with each other. If the cavity depth is set longer, the simple harmonic and superharmonic solution lines are shifted to the left side.

Figure 5(a) shows the positive and negative vibration amplitudes plotted against the dimensionless excitation magnitude for various panel curvatures. Note that the positive and negative vibration amplitudes are not equal for curved panel. Generally, the vibration amplitudes are monotonically increasing with the excitation magnitude. In the cases of 2.5 mm and 4 mm curvatures, the vibration amplitudes abruptly increase around the critical values (i.e., and 5). It is implied that the vibration response changes from nonchaotic to chaotic at this excitation magnitude. The solution lines are not smooth because the amplitudes of the chaotic vibrations are so sensitive. If the excitation magnitude parameter is smaller than critical value, the vibration response is nonchaotic. In this situation, the curved panel vibrates without “snap-through motion.” Figures 6(a)-6(b) and 7(a)-7(b) show the time histories and phase plots for the various dimensionless excitation magnitudes. From these time histories and phase plots, the simple harmonic and quasi-chaotic vibration responses can be seen. In Figures 6(a) and 7(a), the excitation magnitude is far from the critical value. The vibration response is purely simple harmonic. In Figure 6(b), it is close to the critical value. The vibration response starts to change from nonchaotic to chaotic. That is why the phase plot in Figure 7(b) shows both chaotic and periodic properties. If the excitation magnitude is higher than and close to the critical value, the vibration response is sometimes “snap-through” and sometimes is not. It means that the vibration equilibrium position is varying (see Figures 6(c) and 7(c)). The phase plot looks like several ellipses of different sizes and centres overlapping with each other. The ellipses are filled with the solution lines. If the excitation magnitude is much higher than the critical value, the vibration response is clearly considered as “chaotic” and “snap-through” (see Figures 6(d) and 7(d)). The phase plot looks like two equal sized ellipses touching each other. The two ellipses are filled with the solution lines.

**(a)**Vibration amplitude versus excitation magnitude for various curvatures ( mm, m, m, m, , )

**(b)**Vibration amplitude versus excitation magnitude for various excitation frequencies ( mm, m, m, m, , mm)

**(c)**Vibration amplitude versus excitation magnitude for various damping ratios ( mm, m, m, m, mm, )

**(a)**Time history for the simple harmonic case ( mm, m, m, m, , mm, , )

**(b)**Time history for the transition case ( mm, m, m, m, , mm, , )

**(c)**Time history for the chaotic case ( mm, m, m, m, , mm, , )

**(d)**Time history for the snap-through case ( mm, m, m, m, , mm, , )

**(e)**Time history for the simple harmonic case ( mm, m, m, m, , mm, , )

**(f)**Time history for the superharmonic case ( mm, m, m, m, , mm, , )

**(a)**Phase plot for the simple harmonic case ( mm, m, m, m, , mm, , )

**(b)**Phase plot for the transition case ( mm, m, m, m, , mm, , )

**(c)**Phase plot for the chaotic case ( mm, m, m, m, , mm, , )

**(d)**Phase plot for the snap-through case ( mm, m, m, m, , mm, , )

**(e)**Phase plot for the simple harmonic case ( mm, m, m, m, , mm, , )

**(f)**Phase plot for the superharmonic case ( mm, m, m, m, , mm, , )Figure 8(a) shows the modal contributions plotted against the dimensionless excitation magnitude for mm. It is seen that, for the excitation magnitude parameter much higher than the critical value, the contributions of the 1st symmetric, 1st antisymmetric, and 2nd symmetric modes are quite constant (i.e., ≈81%, 13%, and 6%); for the excitation magnitude parameter much lower than the critical value, the modal contributions are ≈89%, 0%, and 11%, respectively (it is implied that the contribution of the 1st antisymmetric mode is zero in the simple harmonic case); and for the excitation magnitude parameter around the critical value, the contribution of the 1st antisymmetric mode abruptly jumps up to ≈38%. Figure 8(b) shows the modal contributions plotted against the dimensionless excitation magnitude for 0 mm curvature (i.e., flat panel). It is seen that the contribution of the 1st antisymmetric mode is always zero. The contribution of the 1st symmetric mode ranges from ≈94.5% to 97.5%.

**(a)**Modal contribution versus excitation magnitude ( mm, m, m, m, , = mm, )

**(b)**Modal contribution versus excitation magnitude ( mm, m, m,