Nonlinear structure-extended cavity interaction simulation using a new version of harmonic balance method

This study addresses the nonlinear structure-extended cavity interaction simulation using a new version of the multilevel residue harmonic balance method. This method has only been adopted once to solve a nonlinear beam problem. This is the first study to use this method to solve a nonlinear structural acoustic problem. This study has two focuses: 1) the new version of the multilevel residue harmonic balance method can generate the higher-level nonlinear solutions ignored in the previous version and 2) the effect of the extended cavity, which has not been considered in previous studies, is examined. The cavity length of a panel-cavity system is sometimes longer than the panel length. However, many studies have adopted a model in which the cavity length is equal to the panel length. The effects of excitation magnitude, cavity depth, damping and number of structural modes on sound and vibration responses are investigated for various panel cases. In the simulations, the present harmonic balance solutions agree reasonably well with those obtained from the classical harmonic balance method. There are two important findings. First, the nonlinearity of a structural acoustic system highly depends on the cavity size. If the cavity size is smaller, the nonlinearity is higher. A large cavity volume implies a low stiffness or small acoustic pressure transmitted from the source panel to the nonlinear panel. In other words, the additional volume in an extended cavity affects the nonlinearity, sound and vibration responses of a structural acoustic system. Second, if an acoustic resonance couples with a structural resonance, nonlinearity is amplified and thus the insertion loss is adversely affected.


Introduction
In recent decades, many researchers have tackled various vibro-acoustic and fluid-structure interaction problems (e.g., [1][2][3][4][5][6][7][8]). Most of these studies that have included sound and vibration analyses of a panel backed by a cavity have concerned the transmission losses of various enclosure panels and have adopted the assumption that the cavity size is equal to the panel size. For example, Pan et al. [9], and Nehete et al. [10] considered a structure coupled with a cavity and subjected to external excitation. The logarithm difference between the external and internal sound pressures is defined as the transmission loss. In fact, the insertion loss of a a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 panel is a very practical indicator and is defined as the logarithm difference between the sound pressures with and without the panel. This can truly reflect the sound reduction ability. In practice, the cavity length of a panel-cavity system is sometimes longer than the panel length, which may significantly affect sound reduction, but has not been considered in the aforementioned work. Although very few studies have examined the nonlinear structural acoustic problem considered in this study, numerous researchers have studied linear structural acoustic or nonlinear structure problems (e.g., [11][12][13][14][15][16][17][18]). The only research work directly related to this study of a nonlinear panel coupled with an extended cavity was done by Lee [19], who investigated the free vibration of a panel backed by an extended cavity only. It did not contain any new solution method and considered the forced vibration and sound responses. Nonlinear vibration (or large amplitude vibration) must be considered in structural acoustic problems for two reasons. First, when a thin panel is used, it vibrates nonlinearly because its weak structural stiffness can cause large amplitude vibrations. Second, a panel mounted very close to a sound source may also vibrate nonlinearly due to the high excitation level imposed by the sound source.
Many studies have applied various solution methods (e.g., the perturbation method, multiple scales method and elliptic integral method) to various nonlinear vibration/oscillation problems or differential equations (e.g., [20][21][22][23][24][25][26]). The total classical harmonic balance method and the incremental harmonic balance method have also been commonly used to solve nonlinear problems (e.g., [27][28][29][30][31]). Because the incremental harmonic balance method eliminates all nonlinear terms during the variational process, some nonlinear behaviour types are missing. The classical harmonic balance method retains all of the nonlinear terms to produce the multiple solutions possible in a set of nonlinear algebraic equations. The nonlinear solutions are assumed by Fourier series expansion to produce a set of nonlinear algebraic equations. From the results of Srirangarajan [27], the results from the classical harmonic balance method closely agreed with those of the closed-form solution. However, it is extremely time consuming to use the classical harmonic balance method to construct higher-order analytical approximations. Hence, in this study, the other harmonic balance method (i.e., the multilevel residue harmonic balance method), which was developed by Leung and Guo [32] and by Hansan et al. [33], is modified and used to solve the nonlinear structural acoustic problem. The main advantage of this method is that the higher-level solutions to any desired accuracy can be obtained easily by solving only one nonlinear algebraic equation and a set of linear algebraic equations in each solution level. Table 1 shows the numbers of algebraic equations required by the three harmonic balance methods to solve the governing equation of a single-mode cubic nonlinear undamped panel vibration. In the previous multilevel residue harmonic balance method, the  [34]. Fig 1 shows that the computational effort required for the new method is less than that for the classical harmonic balance method, which generates more coupled nonlinear algebraic equations. Tables 2 and 3 shows the comparison between the numbers of nonlinear terms generated in the two harmonic balance methods. There are much more nonlinear terms in the classical harmonic balance method. In fact, in the total time of the entire solution process should include the time spent on the formulation derivation as well as the time on the computation. The more nonlinear terms are generated, the more time is spent on the formulation derivation. Unlike finding the roots of the nonlinear algebraic equations (which is done by computer), the formulation derivation and setup of initial guess value for nonlinear algebraic equations are done manually. In this study, only two or three nonlinear algebraic equations are considered for each case. Nowadays, personal computers are very fast so that the computational time is not the bottleneck in the entire solution process and much less than the times spent on the manual steps, which depend on the complexity of the algebraic equations.  Two modes 2 6 10 Using the classical approach and modal analysis, the nonlinear modal governing equations are developed to represent the large-amplitude structural vibration of an enclosure panel coupled with a cavity. Then, the proposed harmonic balance method is used to solve the nonlinear differential equations. Generally, the classical approach is used for those structures and acoustic cavities with regular geometries (e.g. rectangle and circle). In fact, there were some finite element methods developed recently for structural/acoustic problems (e.g. [35][36][37]). They are more suitable for handling models with complicated geometries and boundary conditions. If one of these finite elements is used for the problem in this study, a set of nonlinear differential equations will also be derived from the finite element process. The proposed harmonic balance method is still suitable for solving them.  where k h is the wave number; p Q,h (x,y,z) is the acoustic pressure at the position of (x,y,z); and Q and h represent the panel mode number and harmonic component order, respectively. The total acoustic pressure field is assumed to be equal to the sum of the modal pressure fields:

Theory
where P QJ,h is the acoustic pressure amplitude of the J th acoustic mode and " J is the number of acoustic modes used. φ J is the J th acoustic mode function, which is expressed as follows: where l x , l y and l z are even non-negative integers and a, b and c are the cavity dimensions in the x, y and z directions, respectively. Using the technique of integration by parts in [38], Eq (1) can be expressed in the following form: where The derivative term on the right side of Eq (4) is the pressure gradient on the source panel or the nonlinear panel, which are expressed as follows: At z = 0, At z = c, for Δa +a' > x > Δa and Δb +b' > y > Δb, Otherwise, where ρ o is the air density; Δa = a-a' and Δb = b-b'; and v o (x,y) and v Q (x,y) are the normal velocities of the source panel and the nonlinear panel, respectively. w o (x,y) and w Q (x,y) are the displacements of the source panel and the nonlinear panel, respectively, which are expressed as follows: where A o is the vibration amplitude of the source panel;ϕ o (x,y) is the mode shape of the source panel and assumed as a double sine function (i.e. sin (πx/a') sin (πy/b'); ϕ Q (x,y) is the Q th mode shape of the panel at z = c and also assumed as a double sine function (i.e. sin (mπ(x-Δa)/a')sin (nπ(x-Δb)/b')); and m and n are the panel mode numbers. A Q,h is the h th harmonic component of the panel vibration amplitude and equal to the sum of the h th harmonic components of the vibration amplitudes at different solution levels: By substituting Eqs (5-11) into Eq (4), the acoustic pressure amplitude can be expressed as follows: where Then, the modal contribution of the acoustic pressure acting on the panel is expressed as follows: where Hence, the modal governing equation of the nonlinear panel, which is directly excited by the interior acoustic pressure force and expressed as follows [2,19]: where A Q and ω Q are the modal amplitude and resonant frequency of the Q th panel mode, respectively; m and n are the panel mode numbers; a' is the panel length; γ is the aspect ratio; E is Young's modulus; ν is Poisson's ratio; τ is the panel thickness; ρ is the panel surface density; F Q is the modal force magnitude; and ω is the excitation frequency, and β Q is the nonlinear stiffness coefficient which is given by The proposed solution method [34] was newly presented for nonlinear beam problem in 2017, modified from the previous harmonic method in [32,33], is applied to the nonlinear forced vibration of a flexible panel. The solution form of A Q is expressed as follows: where ε is an embedding parameter, the terms associated with ε 0 , ε 1 and ε 2 represent the zero, 1 st and 2 nd levels terms, respectively; A0 Q , A1 Q and A2 Q are the zero, 1 st and 2 nd level panel responses, respectively; A0 Q,1 , A1 Q,1 and A2 Q,1 are the zero, 1 st and 2 nd level amplitudes of the sin(ωt) responses, respectively; A1 Q, 3 and A2 Q,3 are the 1 st and 2 nd level amplitudes of sin(3ωt) responses, respectively; A2 Q,5 is the 2 nd level amplitude of the sin(5 ωt) response; and A Q,1 , A Q, 3 and A Q,5 are the sin(ωt), sin(3 ω t) and sin(5 ω t) components of the modal amplitude.
Substitute Eq (17) into Eq (15) and collect the terms associated with ε 0 to set up Eq (24). Then, consider the harmonic balance of sin(ωt) to set up Eq (26) and find the zero level solution to A0 Q,1 where R0 Q is the sum of all terms associated with ε 0 Rewrite Eq (26) as the following cubic algebraic equation form in terms of one unknown only where Z0 1 and Z0 3 are the coefficients that precede the linear and nonlinear terms in Eq (26) and Z0 0 is the constant term. If damped cases are considered, rewrite Eq (27) to obtain the following equation: where C0 and D0 are the real and imaginary parts of A0 Q,1 , respectively, and ω po is the resonant frequency of the zero-level resonance. The zero-level resonance occurs when the magnitude of A1 Q,3 is maximum and ω = ω po . Thus, the expression Z0 1 + Z0 3 (C0 2 −D0 2 ), which is function of ω, is equal to 0.
Then, consider the harmonic balance of sin(ωt) to set up the following equation: Note that Eq (35) is a linear equation that contains two unknowns. A0 Q,1 has been found at the zero level.
When substituting Eq (17) into Eq (15), collect those terms which contain the harmonic component of sin(3ωt), to set up the following 1 st level equation (excluding those terms with A2 Q or higher level terms): where R1 Q,3 is the sum of the terms which contain the harmonic component of sin(3ωt). Then, consider the harmonic balance of sin(3ωt) to set up the following equation: By substituting Eq (35) into Eq (38), it can be expressed in the following cubic algebraic equation form, in terms of one unknown only: where Z1 1 , Z1 2 , and Z1 3 are the coefficients that precede the linear and nonlinear terms in Eq (38) and Z1 0 is the zero-order term. If damped cases are considered, rewrite eq (39) to obtain the following equation where C1 and D1 are the real and imaginary parts of A1 Q,3, respectively. The 1 st level resonance occurs (i.e., the magnitude of A1 Q,3 is maximum and ω = ω p1 ) in the frequency range at which the zero-level resonance does not occur. Thus, the damping term in Z1 0 in Eq (41) can be omitted and the following expression is equal to 0 The right side of Eq (41) is purely imaginary (; C1 = 0): Hence, the resonant frequency, ω p1 can be obtained by solving Eqs (43-44) Again, substitute Eq (17) into Eq (15) and collect the terms associated with ε 2 to set up the following 2nd-level equations: where R2 Q is the sum of the terms associated with ε 2 . Then, consider the harmonic balances of sin(ωt) and sin(3ωt) Note that the 1 st nonlinear term in Eq (45) has been used for the harmonic balance of sin (3ωt) in the 1 st level. Thus, it is not considered for the harmonic balance of sin(3ωt) in the 2 nd level. Eqs (47) and (49) are linear equations that contain three unknowns, as A0 Q and A1 Q have been found in the zero and 1 st levels.
When substituting Eq (17) into Eq (15), collect all terms with A0 Q , A1 Q and A2 Q , which contain the harmonic component of sin(5ωt), to set up the following 2nd-level equation: where R2 Q,S is the sum of the terms with A0 Q, , A1 Q, and A2 Q, , which contain the harmonic component of sin(5ωt): Then, consider the harmonic balance of sin(5ωt) to set up the following equation: By substituting Eqs (47) and (49) into Eq (52), it can be expressed as a cubic algebraic equation in terms of one unknown only where Z2 1 , Z2 2, and Z2 3 are the coefficients that precede the linear and nonlinear terms. Z2 0 is the constant term. If damped cases are considered, rewrite Eq (15d) to obtain the following equation: where C2 and D2 are the real and imaginary parts of A2 Q,5 , respectively, and ω p2 is the resonant frequency of the 2nd-level resonance. The 2nd-level resonance occurs (i.e., the magnitude of A1 Q,3 is maximum and ω = ω p1 ) in the frequency range at which the zero level and 1st-level resonances do not occur. Thus, the damping term in Z2 0 in Eq (55) can be omitted and the following expression is equal to 0.
Note that the right side of Eq (55) is therefore purely imaginary (; C2 = 0): Hence, the resonant frequency ω p2 can be obtained by solving Eqs (57-58). Similarly, in the higher-level solution procedures, the harmonic components from the 3rd-or higher-level responses (e.g. A3 Q and A4 Q ) are used to set up the harmonic balance equations. One of these harmonic balance equations is a cubic algebraic equation and the others are linear. By solving these algebraic equations, the harmonic components of higher-level vibration amplitudes can be computed. For example, in the 3rd-level solution procedure, the three algebraic equations from the harmonic balances of sin(ωt), sin(3ωt), and sin(5ωt) are linear, whereas the one from the harmonic balance of sin(7ωt) is cubic. The unknowns in the four algebraic equations are A3 Q,1 , A3 Q,3 , A3 Q,5 , and A3 Q,7. Using the three algebraic equations, A3 Q,1 , A3 Q,3 , and A3 Q,5 can be expressed in terms of A3 Q,7. . Then substitute the expressions into the cubic algebraic equation to find A3 Q,7. Once A3 Q,7 is known, A3 Q,1, A3 Q,3 and A3 Q,5 can be determined using the three algebraic equations. Finally, the overall vibration amplitude of the Q th mode is expressed as follows: where " h ¼ highest solution level or harmonic order: After the nonlinear panel vibration amplitudes are obtained, the modal sound radiation and radiation efficiency can be computed using the modified Rayleigh's integral method [29]: where r is the distance between the panel corner and the observer point; ϕQ is the Q th panel mode shape; θ 1 and θ 2 are the angles between the observer vector and y-axis and between the observer vector and x-axis, respectively (see [29] for details); k h is the wave number; and C a is the speed of sound. Finally, the insertion loss is defined as the logarithm difference between the sound energy radiated from the nonlinear panel and the sound energy radiated from the source panel (i.e., the difference between the radiated sound levels in the cases with/without the nonlinear panel).
where " Q ¼ number of panel modes used σ o = radiation efficiency of the source panel mode.

Results and discussion
In this section, the material properties in the numerical cases considered are as follows: Young's modulus = 7.1 × 10 10 N/m 2 , Poisson's ratio = 0.3 and mass density = 2,700 kg/m 3 . In Tables 4-6, the case is the nonlinear vibration of a single undamped panel subject to uniform excitation (i.e., no cavity).
The modal external excitation term in Eq (2), F Q , is expressed as α Q /α QQ κρg, where κ is the dimensionless excitation parameter and g is 9.81 ms -2 . The panel dimensions are 1 m × 1 m × 2 mm. The first two structural modes are used. Tables 4-6 show the convergence studies of normalised panel vibration amplitudes for various excitation magnitudes and frequencies.
The 2nd-level solutions are normalised as 100. In the small excitation (κ = 0.1), the zero level solutions can achieve an error rate of less than 2% for the three different excitation frequencies.
The 1st-and 2nd-level solutions to two-decimal place accuracy are almost identical. In the other two excitation cases (i.e., κ = 0.5 and κ = 1.2), the maximum difference between the zero-level and 2nd-level solutions is 4.35%, whereas the maximum difference between the 1 st and 2 nd level solutions is only 0.79%. It can be seen that the 1st-level solutions are good enough for three-digit accuracy. Tables 7 and 8 show the modal contributions of the nonlinear panel coupled with an extended cavity for various excitation magnitudes and frequencies. The panel Nonlinear structure-extended cavity interaction simulation dimensions are 1 m × 1 m × 3 mm. The cavity length and width are equal to 1.5 × the panel length. In each case, the sum of the three modal contributions is 100. Note that the 2nd and 4th modes are asymmetrical and are not considered for this symmetric load case. In the lowfrequency excitations, the single-mode approach is sufficient because the contributions from the 3rd and 5th modes are minimal. When the excitation frequency is higher and closer to the resonant frequency of the 3rd mode (see the case of ω = 4 ω 3 in Table 7), the 3rd-mode contribution is significantly higher and the 5th-mode contribution is still minimal. Because the frequency range considered in this study is 0 to 6 ω 1 , the two-mode approach is appropriate. Figs 3 and 4 present the comparisons between the 1st-mode vibration amplitudes of a single undamped panel subject to uniform excitation (i.e., no cavity) obtained with the proposed classical harmonic balance method [31] and the previous multilevel residue harmonic balance method [32,33]. The dimensionless excitation parameter value κ = 0.1. In Fig 3, the results obtained from the two methods are generally in good agreement for both zero-and 1st-level nonlinear solutions. In Fig 4, the 1st-level solution of the previous multilevel residue harmonic balance method is linear and much different from those from the proposed method and the classical harmonic balance method. Only linear equations are set up in the 1st level of the previous multilevel residue harmonic balance method. Fig 5 shows the contribution of sin(3ωt)    approximately ω/ω l = 3. The higher the excitation level, the higher the degree of nonlinearity and the higher the peak frequency that can be seen at the vibration peak. The three solution lines converge closely around the zero-level nonlinear resonance. The slope of the solution line of 3τ around the zero-level nonlinear resonance is the shallowest, whilst the slope of the Nonlinear structure-extended cavity interaction simulation solution line of 0.2τ varies sharply (i.e., from shallow to deep). Generally, the 1st-mode vibration peaks are much higher than the 3rd-mode vibration peaks. The 3rd-mode vibration peak of 0.2τ is very small and more symmetric than that of high excitation. Furthermore, there is no jump phenomenon at the 3rd-mode vibration peak and no 1st-level nonlinear resonance. Unlike that in Fig 7, the higher the excitation level or degree of nonlinearity, the higher the insertion loss dip value that can be seen in Fig 8. In the non-resonant frequency ranges from ω/ω l = 2.5 to ω/ω l = 4.5, the three curves in Fig 8 are very close. Figs 7 and 8 show that the resonant bandwidths are wider for higher excitation magnitudes (note that it is a negative effect on the insertion loss of an acoustic panel). Thus, the linear deign of an acoustic panel may incorrectly estimate the noise reduction performance under high excitation. Figs 9-12 show the vibration amplitudes and insertion losses plotted against the excitation frequency for various cavity lengths and depths. The zero-and 1st-level nonlinear peak amplitudes/insertion loss dip values and frequencies decrease with the cavity length and depth. Unlike those in Fig 7 around the zero-level resonance, the three curves in Figs 9 and 11 are very close but intercept at approximately ω/ω l = 1.5. When ω/ω l < 1.5, the curves of a = a' and c = 0.5a' are the highest in Figs 9 and 11, respectively. When ω/ω l > 1.5, the curves of a = a' and c = 0.5a' are the lowest. The 1st-level nonlinear resonance is not significantly affected by changing the cavity dimensions. It is obvious that the nonlinearity of the structural acoustic system depends greatly upon the cavity size. If the cavity size is smaller, the nonlinearity is higher. A large cavity volume implies a low stiffness or small acoustic pressure transmitted from the source panel to the nonlinear panel. Unlike those in Fig 8, in the non-resonant frequency ranges from ω/ω 1 = 2.5 to ω/ω l = 4.5, the three curves in Figs 10 and 12 are clearly separate. The 2nd structural resonance occurs from ω/ω l = 5.0 to ω/ω 1 = 5.5 in the three cases in Figs 9-10 and the case of c = a' in Figs 11 and 12, whereas the 2nd structural resonance occurs around ω/ω l = 5.8 to ω/ω 1 = 6.2 for the cases of c = 0.5a' and c = 2.5a' in Figs 11 and 12. As mentioned, in the case of c = 0.5a' in Figs 11 and 12, the smaller cavity results in a higher Nonlinear structure-extended cavity interaction simulation Nonlinear structure-extended cavity interaction simulation nonlinearity so that the resonant peak/insertion loss dip and the corresponding peak and dip frequencies are higher. In the case of c = 2.5a' in Figs 11 and 12 (note that the cavity depth is very long and it looks like a tube), although the cavity is much bigger, the nonlinear phenomenon (i.e., jump phenomenon) around the 2nd structural resonance is obvious. The 1st acoustic resonance around ω/ω 1 = 4.7, which is close to and strongly coupled with the 2nd structural resonance, amplifies the nonlinearity. Note that in the case of c = 2.5a' in Fig 5A and 5B, there is no solution found from ω/ω l = 4.7 to 4.8 and thus the solution line is discontinuous there. In Fig 12, there is an anti-resonant peak at approximately ω/ω l = 2.7 in the case of c = 2.5a'. The acoustic pressure forces of the zero-frequency cavity mode and the 1st non-zero frequency cavity mode (their mode numbers are l x = l y = l z = 0 and l x = l y = 0; l z = 1) acting on the nonlinear panel are opposite. Fig 13 shows the two normalised acoustic modal force magnitudes (i.e., | F a,0 | and -|F a,1 |) against the normalised excitation frequency. The two curves intercept at approximately ω/ω l = 2.7. Figs 14 and 15 show the insertion loss dip frequencies and values plotted against the excitation magnitude for various damping ratios. It can be seen that the insertion loss dip values and the corresponding dip frequencies of the three cases increase monotonically with the excitation magnitude. When the excitation magnitude is low, the insertion loss dip frequencies in the three cases converge and the slopes of the three curves deepen. A low excitation magnitude results in linear panel vibrations in the system; the resonant frequency is not significantly affected by the damping. When the excitation magnitude is high, the three curves are almost linear and far from each other. Unlike the resonant frequency, the insertion loss dip value always highly depends on the damping ratio in the system. Thus, the three dip value curves are separate for the entire range of excitation magnitude. Figs 16 and 17 show the insertion loss dip frequencies and values plotted against the cavity depth for various cavity lengths. The insertion loss dip frequencies and values decrease and increase with the cavity depth, Nonlinear structure-extended cavity interaction simulation  ' = b' = c = 1m, a = b = 1.5a'). https://doi.org/10.1371/journal.pone.0199159.g014 Nonlinear structure-extended cavity interaction simulation