- Email: [email protected]
An efficient model for coupling structural vibrations with acoustic radiation
Journal of Sound and Vibration (1995) 182(5), 741–757
AN EFFICIENT MODEL FOR COUPLING STRUCTURAL VIBRATIONS WITH ACOUSTIC RADIATION A. F Analytical Services and Materials, Inc., Hampton, Virginia 23666, U.S.A.
L. M NASA Langley Research Center, Hampton, Virginia 23681-0001, U.S.A.
L. T Courant Institute of Mathematical Sciences, New York, New York 10012, U.S.A. (Received 6 December 1993, and in final form 4 April 1994) In this paper, the problem of coupling between panel vibration and near and far field acoustic radiation is studied. The panel vibration is governed by the non-linear plate equations while the loading on the panel, which is the pressure difference across the panel, depends on the reflected and transmitted waves. Two models are used to solve this structural–acoustic interaction problem. One solves the three-dimensional non-linear Euler equations for the acoustic field coupled with the non-linear plate equations (the fully coupled model). The second uses the linear wave equation for the acoustic field and expresses the load as a double integral involving the panel oscillation (the decoupled model). The panel oscillation governed by a system of integro-differential equations is solved numerically and the acoustic field is then defined by an explicit formula. Numerical results are obtained using the two models for linear and non-linear panel vibration regimes excited by incident waves having different sound pressure levels. The predictions given by these two models are in good agreement, but the computational time needed for the ‘‘fully coupled model’’ is 60 times longer than that for ‘‘the decoupled model’’.
1. INTRODUCTION
The physical problem under consideration is that of structural–acoustic interaction commonly encountered in engineering. A typical example of such problem is the transmission of an external acoustic noise through an airframe into the aircraft’s interior. The interaction occurs because the incident waves excites panel oscillations, which induce transmitted and reflected waves, which, in turn, contribute to the loading on the panel. As shown on Figure 1, the planar interface, the (x, z) plane, is a rigid surface except on the flexible thin panel, D, with panel thickness h much smaller than the size L of the panel, or the size of D, h/L W 1.
(1)
In typical engineering situations, the amplitude of the transverse displacement, h, is much smaller than the thickness of the panel, i.e., h/h W 1.
(2)
741 0022–460X/95/200741 + 17 $08.00/0
7 1995 Academic Press Limited
742
. .
Figure 1. The computational domain.
in this case the panel oscillation is linear. However, the structural–acoustic interaction problem can be non-linear in three different aspects. (1) The panel oscillation can be nonlinear when the transverse displacement is of the order of the panel thickness h: that is, h 1 h.
(3)
In this case, the non-linear plate equations are needed to describe the panel oscillations. (2) In the near field of length scale L, the acoustic field is linear and obeys the simple wave equation but in the far field of a length scale much larger than L, the acoustic field may become non-linear when the second order terms are needed to account for the gradual wave steepening. (3) If the acoustic field is non-uniform and/or the initial pressure variation is not much smaller than the ambient pressure, then the near field acoustic is non-linear. In this case, one needs to solve the Euler equations with appropriate initial and boundary conditions to simulate an incident wave. When trying to simulate the simple problem of panel oscillations with no mean flow on either side, the acoustic fields are governed by the simple wave equation with the same ambient pressure and density. The validity of the linear theory for the acoustic fields is confirmed by the recent investigations of Frendi et al. [1–3] and Maestrello et al. [4]. They showed that, even when the panel vibration is non-linear, the acoustic field can remain linear. This method can be extended to the problem of sound transmission through a moving airframe by simply using the convective wave equation for the external acoustic field and the simple wave equation for the interior one. In addition, the pressure and density are different between the interior and exterior. The pressure difference across the panel induces the panel oscillation which in turn excites the transmitted wave and an additional reflected wave. Thus, the solution of the panel oscillation is coupled with the solutions of the scattered and transmitted waves. For periodic solutions, the waves satisfy the far field radiation conditions for outgoing plane waves (see, for example, reference [5]). For initial value problems, the waves fulfill Huygen’s principle.
–
743
In general, the numerical solution for the scattering of an incident wave by an elastic scatterer, an interface or a panel requires the introduction of a finite computational domain for both acoustic fields. Higher order radiation conditions were derived [6, 7] so that they can be imposed on the boundary of the finite computational domain to give more accurate approximation to the solution in the unbounded domain. Since the size of the computational domain has to be much larger than the size of the panel or a scatterer, the numerical solution of this three-dimensional unsteady problem is very tedious, especially when the computation has to be continued for a long time relative to the period of oscillation of the panel. Furthermore, the accuracy of the solution depends not only on that of the numerical solution of the differential equations but also on the approximate boundary conditions. Refinements of the grid size and time step have to go with an enlargement of the computational domain. It is desirable to find exact boundary conditions so that the computational domain does not have to be much larger than the scatterer and is independent of the choice of grid size. Such exact conditions were presented in reference [8]. With the scatterer inside the computational domain, the integral representation of the solution of the simple wave equation is applicable to the region outside and on the boundary of the computational domain. In this paper, the scatterer is a panel embedded in a rigid plane. The integral representation of the acoustic field can be applied on the panel and the integrand involves only the normal velocity of the panel [9, 10]. Thus we have a system of integro-differential equations for the panel oscillation decoupled from the acoustic field. This system is referred to as the decoupled model. The integro-differential equations can be reduced to partial differential equations when the acoustic wavelength is much smaller than L [9–11]. This simplification is not used in this paper because the acoustic wavelength is of the same order as the plate length, L. Therefore, a set of integro-differential equations is integrated to obtain the panel oscillation and then using the integral representation the acoustic field is calculated. An extension of this work to the problem of noise transmission in a supersonic flight at high altitude is given in Frendi et al. [12]. Note that the solution of the decoupled panel oscillation is an unsteady two-dimensional problem in the finite domain D and is several orders of magnitude simpler than the solution of the fully coupled problem. To show the efficiency of the decoupled model and the accuracy of the solution even when the panel oscillation is non-linear, we compare the solution of this decoupled system with the solution of the fully coupled problem for which we solve the three-dimensional non-linear Euler equations for the acoustic field coupled with the non-linear equations for the panel oscillations. In the next section, we present the mathematical formulations of these two models. The various numerical techniques used to solve the problems involved are described in section 3. The results and discussion are given in section 4, and the conclusions are in section 5.
2. FORMULATION OF THE TWO MODELS
The physical problem being studied is that of linear and non-linear oscillations of a typical aircraft panel excited by harmonic plane waves at normal incidence. The panel is clamped on to a large rigid plate. Let the typical panel be represented by a rectangular domain D in the (x, z) plane, with length L, width W and one vertex located at (x0 , z0 ), i.e., D = {(x, z) = x0 Q x Q x0 + L,
z0 Q z Q z0 + W}.
(4)
. .
744
The complement of D in the (x, z) plane represents the rigid plate. The incident domain is the half-space y q 0 and the transmitted domain is the half-space y Q 0. Both domains have the same ambient pressure pa and temperature Ta, and hence the same speed of sound C and density ra . We denote the pressure and velocity potential of the acoustic field by p and f and use the superscripts + and − to denote the quantities on the sides y q 0 and y Q 0. On the side y q 0, the superscripts i and r are used to denote the quantities associated with the incident and reflected waves, respectively. These two waves are mirror images with respect to the (x, z) plane. When the incident waves are plane waves advancing in the direction opposite to the y-axis and hitting the (x, z) plane at t = 0, we have f (i)(t, x, y, z) = f (Ct + y)
f (r)(t, x, y, z) = f (Ct − y),
and
(5)
with f(j) = 0 for j E 0. Therefore f (r) = 0 when 0 e y Q Ct. Transverse oscillation, h, of the rectangular flexible panel, D, is excited by the pressure difference across the panel for t q 0 and the oscillation in turn induces scattered waves f s in the incident side and transmitted waves f t. Both waves have the homogeneous initial conditions at t = 0, f (s) = ft(s) = 0,
y q 0,
f (t) = ft(t) = 0,
y Q 0.
(6)
Under linear theory, we have f (r)(t, x, y, z) = f (i)(t, x, −y, z),
f (t)(t, x, −y, z) = −f (s)(t, x, y, z),
for y q 0.
(7)
The wavefronts of f (s) and f (t) are the envelopes of the sonic spheres with radius Ct and centered on the panel in y q 0, and y Q 0, respectively. The pressures in the incident and transmitted sides are p+ − pa = p (i) + p (r) + p (s),
y q 0,
p− − pa = p (t),
y Q 0,
(8)
and the pressure difference across the panel is Dp = p− = p+ = 2p (t) − [p (i) + p (r)],
y = 0,
(x, z) $ D.
(9)
In the next sections we formulate two mathematical models to compare with and complement the experimental studies. We present the system of non-linear equations governing the acoustic fluid and panel oscillation (the fully coupled model) in section 2.1 and formulate the system of integro-differential equations for non-linear panel oscillation, including the effect of a linear acoustic field (the decoupled model), in Section 2.2. 2.1. In this model, the structural–acoustic interaction is analyzed by solving the three-dimensional non-linear Euler equations together with the non-linear plate equations. The configuration of the computational domain is shown in Figure 1. In Cartesian co-ordinates, x, y and z, the compressible, non-linear Euler equations can be written in conservation form as Qt = Fx + Gy + Hz ,
(10)
where Q is the vector (r, ru, rv, rw, e), r is the density, ru, rv and rw are the x, y and z momenta respectively, and e is the total energy per unit volume, given by e = 12 r(u 2 + v 2 + w 2) + rcv T,
(11)
with cv being the specific heat at constant volume. In equation (10), the functions F, G and H are:
–
F J G ru G 2 ru + p G G F =G ruv G, G ruw G u(e + p)G G f j
F G rv J G G ruv G G =G rv 2 + p G, G rvw G v(e + p)G G f j
745
F G rw J G G ruw G H = G rvw G. 2 G rw + p G w(e + p)G G f j
(12)
In addition to equation (10), the equation of state of an ideal gas is used: p = rRT,
(13)
where p is the pressure, R the gas constant, and T the temperature. Since we assume that the incident wave hits the panel at t = 0, the flow on the transmitted side, y Q 0, is at rest for t E 0. The initial data at t = 0 for y Q 0 are u = v = w = 0,
p = pa
and
r = ra .
(14)
On the rigid plate, we have zero normal velocity, v(x, 0, z) = 0
for
(x, z) ( D.
(15)
The motion of the flexible panel, D, is described by a system of three, non-linear partial differential equations, given by reference [13]: D9 4h + rp hhtt + ght = Dp +
Eh [(u 0 + 1 h 2 )(h + nhzz ) 1 − n 2 x 2 x xx
+(wz0 + 12 hz2 )(hzz + nhxx ) + (1 − n)hxz (uz0 + wx0 + hx hz )],
(16)
0 0 uxx + d1 uzz0 + d2 wxx = −hx (hxx + d1 hzz ) − d2 hz hxz , 0 0 + d2 uxz = −hz (hzz + d1 hxx ) − d2 hx hxz , wzz0 + d1 wxx
where 9 4h = hxxxx + 2hxxzz + hzzzz , d1 = (1 − n)/2,
d2 = (1 + n)/2,
(17)
D = Eh 3/12(1 − n 2),
(18)
u and w are the in-plane displacements, and h is the transverse displacement. The physical constants of the panel appearing in equation (16) are the stiffness (D), the density (rp ), the thickness (h), the physical damping (g), the modulus of elasticity (E) and the Poisson ratio (n). The system of equations (16) is solved subject to the homogeneous initial conditions 0
0
at t = 0,
u 0 = w 0 = h = ht = 0,
(19)
and the clamped boundary conditions on C given by x = x0 , x0 + L,
u 0 = w 0 = h = hx = 0,
z = z0 , z0 + W,
u 0 = w 0 = h = hz = 0.
(20)
In equation (16), the load Dp defined by equation (9) contains the coupling with the acoustic field p t and the forcing term (p (i) + p (r) ), which represents the load on the panel if the panel were rigid. We assume the forcing term to be a harmonic wave of the form p (i) + p (r) = (e sin(vt) + e* sin(2vt))H(t)
on y = 0+,
(21)
where e, e* and v are the amplitudes and frequency of the wave and H(t) denotes the Heaviside unit step function. For small values of e, e* = 0 (linear), while for large values
. .
746
of e, e* = e/2 (non-linear). This input is used based on the knowledge of a typical experimental input. The load on the panel D given by equation (9) becomes Dp = [2p t − e sin(vt) − e* sin(2vt)]H(t).
(22)
Another condition coupling the acoustic fields and the panel oscillation is the kinematic condition, ht (t, x, z) = v2(t, x, 02, z),
(x, z) $ D.
(23)
2
It is imposed on y = 0 because of the small panel displacement, equation (1). The numerical scheme for the solution of the non-linear plate equation (16) and the Euler equations (10) for the transmitted waves with homogeneous initial data will be described in section 3. 2.2. Now we analyze the non-linear panel oscillation excited by a weak incident pressure wave, p (i), under the assumption that the pressure fluctuation remains much smaller than the ambient pressure pa , i.e., (p − pa )/pa W 1.
(24)
Consequently, the small disturbance theory is applicable to the acoustic fields in the incident and transmitted sides and the panel/acoustic interaction problem is described by equations (4)–(9). The velocity potential F2(t, x, y, z) is governed by the simple wave equation, 2 2 (C−21tt2 − 1xx − 1yy − 1zz2 )F2 = 0,
for 2y q 0,
(25)
and the acoustic pressure and velocity v2 are related to the potential by p2 = −pa 1t F2
v2 = 9F2.
and
(26)
In particular, we have v = 1y F, where v denotes the transverse velocity. The transverse displacement of the panel, h(t, x, z), is governed by the system of partial differential equations (16). To produce the special forcing term, (p (i) + p (r) ), acting on the panel specified by equation (21) in section 2.1, the incident potential should be f (i)(t, x, y, z) =
0 $ 0 1%
e y cos v t + 2ra v C
+2
$ 0 1%1 0 1
e* y cos 2v t + e C
H t+
y . (27) C
Note that the incident wave f (i) is a solution of equation (25) in the whole space. When the panel is rigid, the incident wave is reflected by the rigid (x, z) plane and the reflected wave is f (r)(t, x, y, z) = f (i)(t, x, −y, z)
for y e 0.
(28)
On the side y Q 0, the acoustic fluid remains at rest because there is no transmitted wave. For a flexible panel, the oscillation by the pressure difference across the panel −2p (i)(t, x, 0+, z) induces a scattered wave f (s) on the incident side and a transmitted wave f (t). Thus we write F = f (i) + f (r) + f (s),
y q 0,
(29)
and F = f (t),
y Q 0.
(30)
–
747
Note that f and f are governed by the simple wave equation (25) in y q 0 and y Q 0, respectively. The kinematic conditions on the (x, z) plane, equations (23) and (26), become (s)
(t)
1y f (s)(t, x, 0+, z) = 1t h(t, x, z)
1y f (t)(t, x, 0−, z) = 1t h(t, x, z),
and
(31)
where h denotes the extension of the transverse displacement of the panel, D, to the (x, z) plane, i.e., h = 0,
(x, z) ( D.
(32)
If the incident wavefront hits the panel at t = 0, we can impose the homogeneous initial conditions on f (s), f (t) and h, f (s) = 0,
1t f (s) = 0,
f (t) = 0,
1t f (t) = 0
and h = 0,
1t h = 0,
for t E 0.
(33)
Since the ambient fluid above and below the panel are the same, equation (31) implies that the velocity potential induced by the panel oscillation has to be antisymmetric in y, i.e., f (t)(t, x, y, z) = −f (s)(t, x, −y, z),
u q 0.
(34)
Then, the pressure difference across can be written as Dp = 2ra [1t f (i)(t, x, 0+, z) + 1t f (s)(t, x, 0+, z)],
(x, z) $ D.
(35)
The derivation of the system of equations (24)–(35) completes the formulation of the structural–acoustic interaction problem using linear theory for the acoustic fields. The panel oscillation, which can be non-linear, is governed by the system of equations (16)–(20). The velocity potential f (s) induced by the panel oscillation h is governed by the wave equation (25), the initial conditions (33) and the boundary condition (31). The solution is given by the Kirchhoff formula, 1 f (s)(t, x, y q 0, z) = − 2p
gg
{ht (t, x', z')} dx' dz', R G
(36)
where R = [(x − x')2 + y 2 + (z − z')2]1/2 denotes the distance from a point (x, y q 0, z) to a source at (x', 0, z') and {·} denotes the retarded value, i.e., {ht (t, x', 0, z')} = ht (t − R/C, x', 0, z').
(37)
The domain of dependence of f (s)(t, x, y, z) is the circular disc H in the (x', z') plane, i.e., H= R E Ct
or
r 2 = (x' − x)2 + (z' − x)2 E C 2t 2 − y 2.
(38)
The domain of integration in equation (36) is the intersection of H and the panel, i.e., G=H + D
(39)
Now we introduce the polar co-ordinates, r, u, centered at (x, z), i.e., x' − x = r cos u,
z' − z = r sin u,
(40)
and equation (36) becomes, 1 2p
f (s)(t, x, y q 0, z) = −
g$g 0
dugt t −
G
R , r, u C
1%
r dr , R
(41)
with g(t, r, u) = h(t, x + r cos u, z + r sin u).
(42)
748 . . Here, R = (r 2 + y 2)1/2 represents the slant height of a circular cone with vertex P(x, y, z), a vertical axis and a base circle of radius r in the (x, z) plane. On account of equation (32), we can extend the domain of integration of H and rewrite equation (41) as an iterated integral, 1 2p
f (s)(t, x, y q 0, z) = −
g $g Ct
0
2p
0
dugt t −
0
R , r, u C
1%
r dr . R
(43)
The integral in u represents the contribution of the sources on a base circle of radius r. In equation (35), we need to relate the unknown, ft(s) (x, 0+, z) on the panel, to h(t, x, y) for (x, z) $ D. This is obtained by differentiating the above equation with respect to t and using equation (32). This results in replacing f (s) on the left side of equation (43) and gt in the integrand by ft(s) and gtt , respectively. As y :0+, R :r and the circular disc H is bounded by the sonic circle r = Ct. Equation (43) yields ft(s) (t, x, 0+, z) = −
1 2p
gg 0 2p
Ct
gtt t −
0
0
1
r , r, u dr du. C
(44)
Note that in this form we remove the kernel 1/R in equation (36), which becomes singular as y :0 and r :0 and shows that the area of the domain of integration is bounded above by 2pCt. By using equations (35) and (44), equation (16) becomes a system of integro-differential equations for the panel oscillation, h. The initial and boundary conditions for h are equations (19) and (20). Thus, we complete the formulation of the decoupled model. The numerical solution of this system will be described in section 3. 3. NUMERICAL METHODS
For the fully coupled model, the non-linear Euler equations (10) are solved using an explicit finite difference scheme. The scheme, which is a generalization of MacCormack’s scheme obtained by Gottlieb and Turkel [14], is fourth order accurate in space and second order accurate in time. Further details on the implementation of the scheme can be found in Frendi et al. [1–3]. The physical boundary of the computational domain (the bottom boundary; see Figure 1) is composed of a flexible panel clamped between rigid plates. The boundary conditions employed on both the flexible and rigid plates for the Euler equations are v = 0,
T = Tw ,
(45)
over the rigid plates, and v = 1h/1t,
T = Tw ,
(46)
over the flexible panel. In equations (45) and (46), Tw is a specified wall temperature; in this paper it is taken to be the free-stream temperature, Ta . Over the rigid plates, the x and z components of the velocity (u and w, respectively) are obtained through linear extrapolation from the interior of the computational domain. The pressure boundary condition is obtained using the normal momentum equation by simply imposing the normal gradient of the sum of pressure and vertical momentum flux to be zero, i.e., 1y (p + rv 2) = 0. Over the flexible panel, the x and z velocity components and the pressure are extrapolated from the interior of the computational domain. The remaining numerical boundary conditions of the computational domain are derived using the method of characteristics [15]. One should mention that the characteristic boundary conditions and the extrapolation are only first order accurate, since the information on the boundary
–
749
depends on the point next to the boundary only, while the interior scheme is fourth order. These boundaries are believed to be a source for numerical error in the fully coupled model. The truncation error in the interior of the domain is of order (Dx)4, while near the boundary it is of order (Dx); this jump in truncation error is a major source for the spurious waves that propagate back into the domain. The numerical error can show up in the results in the form of additional non-physical frequencies. The dimension of the computational domain is 3.02 m in the x, y and z directions, respectively. The number of computational points used are 141, 181, and 141 in the respective directions. The non-linear plate equations (16) are solved using a finite element method developed by Robinson [16]. The panel is 30·5 cm long, 20·32 cm wide and 0·102 cm thick, and the number of elements used are six and eight in the length and width, respectively. Since the grid used for solving the plate equations is rectangular, it is easier to evaluate the integral in equation (44) in Cartesian co-ordinates: ft(2) (t, x, 0+, z) = −
1 2p
g g 0 Xmax
Xmin
Zmax
Zmin
htt t −
1
dx' dz' r , , x', 0, z' r C
(47)
where (Xmin , Zmin ) and (Xmax , Zmax ) lie within the borders of the smallest and largest circles on the panel, respectively. The radius of the smallest circle is the smaller of Dx or Dz, and that of the largest circle is R = Ct. The presence of the small circle results from the singularity at r = 0, which appears in the integral (47). The contribution to the pressure field of the point at r = 0, (Xp , Zp ), is found by expanding htt (t − r/C, x', 0, z') in a Taylor
Figure 2. (a) The time history and (b) power spectral density of the panel center displacement for an excitation amplitude of 130 dB (linear). ——, Fully coupled model; . . . . . . , decoupled model.
750
. .
Figure 3. (a) The time history and (b) power spectral density of the transmitted near field pressure, 2·54 cm from the panel center, excitation amplitude 130 dB (linear). ——, Fully coupled model; . . . . . . , decoupled model.
series. The expansion is truncated after the third order term. The contribution of the leading term in the expansion is htt (t, Xp , 0, Zp )R1 Du, where R1 is equal to the smaller of Dx or Dz and Du is equal to 2p for an interior point, p for a boundary point and p/4 for a corner point. The contribution of the various points on the panel that lie within the sonic circle (R E CT) is calculated by first integrating in x using a combination of Simpson and trapezoidal rules. The result is then integrated in z using Simpson’s rule of integration. Because of the presence of the retarded time in the integral, the acceleration of the plate (ht t) is stored at each point for several time-steps. The evaluation of the double integral can use a much larger time step than the Dt for the integration of the differential equations, yet having the same degree of accuracy. For a plate of given dimensions (L, W), and for a fixed time-step Dt, the maximum number of time-steps to be stored is N = zL 2 + W 2/(CDt). The number N has to be changed when calculating the radiated pressure away from the plate. The number N is related to the radius of the largest sonic circle in the domain that contains the flexible panel. For the cases presented in the next section, N = 1503 for a time-step Dt = 4 × 10−6 s. 4. RESULTS AND DISCUSSION
The numerical schemes for the two models presented in section 3 are used to predict the vibration of a flexible panel and the resulting acoustic radiation. The panel is forced
–
751
to vibrate by harmonic plane acoustic waves at normal incidence. The frequency of the source is 751 Hz for low amplitude excitations, and a combination of 751 Hz and 1502 Hz for high amplitude excitations. This method of input was chosen based on the knowledge of the experimental input, however, no direct comparisons with experiments are made in this paper. The frequency 751 Hz corresponds to a natural frequency of the flexible panel. Two different amplitudes of the incident waves are used. The properties of the panel, which are considered to be uniform, are density rp = 4450·15 kg/m3, modulus of elasticity E = 1·10316 × 1011 N/m2, Poisson ratio n = 0·33 and a damping ratio of 0·01 is used. The acoustic fluid properties are temperature Ta = 288·33 K, density ra = 1·23 kg/m3, pressure ra = 1·013 × 105 N/m2 and sound speed C = 340 m/s. The specific heat at constant volume is cv = 1·004 kJ/(kg K), the ratio of specific heats is g = cp /cv = 1·4. In the far field, the fluid is at rest. The variables plotted in Figures 2–5 are non-dimensional. The reference quantities are given by (x, y, z, h)ref = lref , rref = ra ,
tref = lref /C,
(u, v, w, 1h/1t)ref = C,
Tref = C 2/cv , (p, e)ref = ra C 2,
(48)
where the reference length is lref = 0·3048 m. The configuration of the computational domain, a square box with the lower side composed of a flexible panel clamped between rigid plates, is shown in Figure 1. For a low excitation amplitude, 130 dB or 6·8 × 10−4 atm, the panel response is linear as shown
Figure 4. (a) The time history and (b) power spectral density of the transmitted far field pressure, 203·2 cm from the panel center, excitation amplitude 130 dB (linear). ——, Fully coupled model; . . . . . . , decoupled model.
752
. .
Figure 5. (a) The time history and (b) power spectral density of the strain response of the panel for an excitation amplitude of 174 dB (non-linear). ——, Fully coupled model; . . . . . . , decoupled model.
in Figures 2(a) and (b). The time history of the normalized panel center displacement (h/h, where h is the panel thickness) shows a periodic behavior with one dominant frequency f (see Figure 2(a)). In this case, both models predict the same panel response as shown by the figure. This result is further confirmed by the power spectral density of the displacement, Figure 2(b), which shows a strong spike at f = 751 Hz. The near field transmitted pressure, 2·54 cm from the panel center or L/12, is shown in Figures 3(a) and (b). The time history of the transmitted near field pressure, Figure 3(a), shows a periodic behavior with one dominant frequency. The power spectrum of this time signal shows a strong spike at the fundamental frequency f = 751 Hz (see Figure 3(b)). Both models predict the same near field pressure. A small wiggle is shown by the power spectrum of the decoupled model near the fundamental, its magnitude is negligible compared to that of the fundamental and therefore has no effect on the results as is shown in Figure 3(a). This wiggle is believed to be caused by the difference in the length of the time samples used to calculate the power spectrum. In Figures 4(a) and (b) is shown the transmitted far field pressure, 203·2 cm from the panel center or 5L. Similar to the near field, the far field pressure is periodic with one dominant frequency f = 751 Hz as indicated by both the time history (Figure 4(a)) and the power spectrum (Figure 4(b)). The two figures also show that both models predict the same far field pressure. When the level of the excitation is increased to 174 dB, or 0·22 atm, the response of the panel becomes non-linear. The strain response of the panel at the middle of the long edge is shown in Figures 5(a) and (b). The time history of the strain, Figure 5(a), shows a
–
753
periodic response with a period 4T corresponding to the frequency f/4, and a strong harmonic 2f characterized by the break-up of the positive peaks. This observation is further confirmed by the power spectral density of the strain, Figure 5(b), which shows the presence of several frequencies. Most of the frequencies are either harmonics or subharmonics of the fundamental frequency f. This spectrum is characteristic of a non-linear panel response. The panel response obtained by the ‘‘decoupled model’’ is in fairly good agreement with that obtained by the ‘‘fully coupled model’’. In particular, there is a good agreement in the power spectra of the response predicted by the two models. The predicted panel center displacement is shown in Figures 6(a) and (b). The time history of the normalized panel center displacement is shown on Figure 6(a). Contrary to the strain response, the displacement response shows less non-linearity. In Figure 6(a) is shown a dominant fundamental frequency f with a weak subharmonic f/2. This is further confirmed by the power spectrum of the displacement, which shows a weak harmonic 2f in addition to f and f/2. The predictions of the two models are in fairly good agreement. Since the coupling between the acoustic fluid and structure is obtained through the out-of-plane velocity of the panel, a comparison between the time histories of the non-dimensional panel velocities at the panel center is shown in Figure 7. As expected, the centerline panel velocity shows a more pronounced effect of the harmonic 2f. This is particularly apparent in the bottom part of the time histories. The predictions of both models are in fairly good agreement.
Figure 6. (a) The time history and (b) power spectral density of the panel center displacement for an excitation amplitude of 174 dB (non-linear). ——, Fully coupled model; . . . . . . , decoupled model.
754
. .
Figure 7. The time history of the panel center out-of-plane velocity for an excitation amplitude of 174 dB. ——, fully coupled model; . . . . . . , decoupled model.
In Figures 8(a) and (b) is shown the transmitted near field pressure, 2·54 cm form the panel center or L/12. The time history, Figure 8(a), shows a periodic behavior with a dominant fundamental f and a strong harmonic 2f. The effect of the harmonic is noticeable
Figure 8. (a) The time history and (b) power spectral density of the transmitted near field pressure, 2·54 cm from the panel center, excitation amplitude 174 dB (non-linear). ——, Fully coupled model; . . . . . . , decoupled model.
–
755
Figure 9. (a) The time history and (b) power spectral density of the transmitted far field pressure, 203·2 cm from the panel center, excitation amplitude 174 dB (non-linear). ——, Fully coupled model; . . . . . . , decoupled model.
in the break-up of the positive peaks of the time history. The power spectrum of the near field pressure, Figure 8(b), shows a weak subharmonic f/2 with a strong fundamental f and harmonic 2f. Notice that the difference in level between f and 2f is about 10 dB. There is a fairly good agreement between the predictions of the two models. The transmitted far field pressure, 203·2 cm from the center of the panel, is shown in Figures 9(a) and (b). The time history of the far field pressure shows a periodic behavior with period T/2 corresponding to 2f. It is important to note that in the near field (Figure 8(a)), the fundamental frequency f is dominant; however, in the far field (Figure 9(a)) it is the harmonic 2f that is dominant. The power spectral density of the far field pressure, Figure 9(b), shows that the fundamental and the harmonic 2f have nearly the same level. This indicates that the two frequencies do not decay at the same rate into the far field. This behavior is due to the fact that the panel is not a compact source but rather a set of distributed sources. The predictions of the two models are in farily good agreement. The level of 2f predicted by the fully coupled model is 3 dB higher than that predicted by the decoupled model. This is believed to be due to a combination of numerical errors and weakly non-linear wave propagation effects. From a computational view point, it is important to compare the performance of the two models based on the CPU time required by each calculation. In the linear vibration regime, the ‘‘fully coupled model’’ used 36 000 s of CPU time on a Cray-ymp to advance the calculation by 10 000 time-steps, whereas the ‘‘decoupled model’’ used only 1000 s for
756
. .
the same calculation. In the non-linear vibration regime, grid refinements were needed to resolve the large gradients both on the panel and in the radiation field. Therefore, in order to advance the calculation by 10 000 time-steps, 72 000 s were used by the ‘‘fully coupled model’’ while the ‘‘decoupled model’’ used only 1200 s.
5. CONCLUSIONS
An efficient model for coupling the vibration of a panel to the on surface acoustic radiation is derived. The model uncouples the panel vibrations from the acoustic wave propagation problem. The results show that this model, referred to as ‘‘decoupled model’’, accurately predicts the panel response and acoustic radiation in the linear and non-linear vibration regimes as long as the pressure variation in the acoustic fluid remains much smaller than the ambient pressure. For the cases studied in this paper, the computational cost of the numerical integration of this model is 36 times cheaper in the linear regime and 60 times cheaper in the non-linear regime than the cost of the ‘‘fully coupled model’’.
ACKNOWLEDGMENTS
The first author acknowledges the support of NASA Langley Research Center under contract NAS1-19700. Partial support was also provided by NASA Langley to the third author under contracts NAS1-18605 and NAS1-19480 while in residence at ICASE.
REFERENCES 1. A. F and J. R 1993 American Institute of Aeronautics and Astronautics Journal 31, 1992–1997. On the effect of acoustic coupling on random and harmonic plate vibrations. 2. A. F, L. M and A. B 1994 Journal of Sound and Vibration 177, 207–226. Coupling between plate vibration and acoustic radiation. 3. A. F, L. M, J. R and A. B 1993. American Institute of Aeronautics and Astronautics Paper 93–4367. On acoustic radiation from a vibrating panel. 4. L. M, A. F and D. E. B 1992. American Institute of Aeronautics and Astronautics Journal 30, 2632–2638. Nonlinear vibration and radiation from a panel with transition to chaos induced by acoustic waves. 5. P. M. M and K. U. I 1968 Theoretical Acoustics. New York: McGraw-Hill. 6. G. A. K and C. M 1980 SIAM Journal for Scientific and Statistical Computing 1, 371–385. Solving the Helmholtz equation for exterior problems with variable index of refraction: I. 7. A. B and E. T 1980 Communications in Pure and Applied Mathematics 23, 707–725. Radiation boundary conditions for wave-like equations. 8. L. T and M. J. M 1986 Journal of the Acoustical Society of America 80, 1825–1827. Exact boundary conditions for scattering problems. 9. M. J. M and L. T 1989 Wave Motion 11, 545–557. Scattering of an incident wave from an interface separating two fluids. 10. M. J. M and L. T 1992 Northwestern University Report 9130. Scattering of acoustic wave with a flexible wall. 11. L. T 1993. Paper presented at a workshop on Perturbation Methods in Physical Mathematics, Renselear Polytechnic Institute, Troy, New York. On surface conditions for structural acoustic interactions in moving media. 12. A. F, L. M and L. T 1993. Institute for Computer Applications in Science and Engineering (ICASE) , Report No. 93–18. An efficient model for coupling structural vibrations with acoustic radiation. 13. C. Y. C 1990 Nonlinear Analysis of Plates. New York: McGraw-Hill. 14. D. G and E. T 1976 Mathematics of Computation 30, 703–723. Dissipative two-four methods for time dependent problems.
–
757
15. S. S. A, W. S. D, D. G, D. H. R and J. C. T 1990 Institute for Computer Applications in Science and Engineering (ICASE) , Report 90-16. Secondary frequencies in the wake of a circular cylinder with vortex shedding. 16. J. H. R 1990 Master’s Thesis, Mechanical Engineering Department, Old Dominion University, Norfolk, Virginia. Finite element formulation and numerical simulation of the random response of composite plate.
AN EFFICIENT MODEL FOR COUPLING STRUCTURAL VIBRATIONS WITH ACOUSTIC RADIATION A. F Analytical Services and Materials, Inc., Hampton, Virginia 23666, U.S.A.
L. M NASA Langley Research Center, Hampton, Virginia 23681-0001, U.S.A.
L. T Courant Institute of Mathematical Sciences, New York, New York 10012, U.S.A. (Received 6 December 1993, and in final form 4 April 1994) In this paper, the problem of coupling between panel vibration and near and far field acoustic radiation is studied. The panel vibration is governed by the non-linear plate equations while the loading on the panel, which is the pressure difference across the panel, depends on the reflected and transmitted waves. Two models are used to solve this structural–acoustic interaction problem. One solves the three-dimensional non-linear Euler equations for the acoustic field coupled with the non-linear plate equations (the fully coupled model). The second uses the linear wave equation for the acoustic field and expresses the load as a double integral involving the panel oscillation (the decoupled model). The panel oscillation governed by a system of integro-differential equations is solved numerically and the acoustic field is then defined by an explicit formula. Numerical results are obtained using the two models for linear and non-linear panel vibration regimes excited by incident waves having different sound pressure levels. The predictions given by these two models are in good agreement, but the computational time needed for the ‘‘fully coupled model’’ is 60 times longer than that for ‘‘the decoupled model’’.
1. INTRODUCTION
The physical problem under consideration is that of structural–acoustic interaction commonly encountered in engineering. A typical example of such problem is the transmission of an external acoustic noise through an airframe into the aircraft’s interior. The interaction occurs because the incident waves excites panel oscillations, which induce transmitted and reflected waves, which, in turn, contribute to the loading on the panel. As shown on Figure 1, the planar interface, the (x, z) plane, is a rigid surface except on the flexible thin panel, D, with panel thickness h much smaller than the size L of the panel, or the size of D, h/L W 1.
(1)
In typical engineering situations, the amplitude of the transverse displacement, h, is much smaller than the thickness of the panel, i.e., h/h W 1.
(2)
741 0022–460X/95/200741 + 17 $08.00/0
7 1995 Academic Press Limited
742
. .
Figure 1. The computational domain.
in this case the panel oscillation is linear. However, the structural–acoustic interaction problem can be non-linear in three different aspects. (1) The panel oscillation can be nonlinear when the transverse displacement is of the order of the panel thickness h: that is, h 1 h.
(3)
In this case, the non-linear plate equations are needed to describe the panel oscillations. (2) In the near field of length scale L, the acoustic field is linear and obeys the simple wave equation but in the far field of a length scale much larger than L, the acoustic field may become non-linear when the second order terms are needed to account for the gradual wave steepening. (3) If the acoustic field is non-uniform and/or the initial pressure variation is not much smaller than the ambient pressure, then the near field acoustic is non-linear. In this case, one needs to solve the Euler equations with appropriate initial and boundary conditions to simulate an incident wave. When trying to simulate the simple problem of panel oscillations with no mean flow on either side, the acoustic fields are governed by the simple wave equation with the same ambient pressure and density. The validity of the linear theory for the acoustic fields is confirmed by the recent investigations of Frendi et al. [1–3] and Maestrello et al. [4]. They showed that, even when the panel vibration is non-linear, the acoustic field can remain linear. This method can be extended to the problem of sound transmission through a moving airframe by simply using the convective wave equation for the external acoustic field and the simple wave equation for the interior one. In addition, the pressure and density are different between the interior and exterior. The pressure difference across the panel induces the panel oscillation which in turn excites the transmitted wave and an additional reflected wave. Thus, the solution of the panel oscillation is coupled with the solutions of the scattered and transmitted waves. For periodic solutions, the waves satisfy the far field radiation conditions for outgoing plane waves (see, for example, reference [5]). For initial value problems, the waves fulfill Huygen’s principle.
–
743
In general, the numerical solution for the scattering of an incident wave by an elastic scatterer, an interface or a panel requires the introduction of a finite computational domain for both acoustic fields. Higher order radiation conditions were derived [6, 7] so that they can be imposed on the boundary of the finite computational domain to give more accurate approximation to the solution in the unbounded domain. Since the size of the computational domain has to be much larger than the size of the panel or a scatterer, the numerical solution of this three-dimensional unsteady problem is very tedious, especially when the computation has to be continued for a long time relative to the period of oscillation of the panel. Furthermore, the accuracy of the solution depends not only on that of the numerical solution of the differential equations but also on the approximate boundary conditions. Refinements of the grid size and time step have to go with an enlargement of the computational domain. It is desirable to find exact boundary conditions so that the computational domain does not have to be much larger than the scatterer and is independent of the choice of grid size. Such exact conditions were presented in reference [8]. With the scatterer inside the computational domain, the integral representation of the solution of the simple wave equation is applicable to the region outside and on the boundary of the computational domain. In this paper, the scatterer is a panel embedded in a rigid plane. The integral representation of the acoustic field can be applied on the panel and the integrand involves only the normal velocity of the panel [9, 10]. Thus we have a system of integro-differential equations for the panel oscillation decoupled from the acoustic field. This system is referred to as the decoupled model. The integro-differential equations can be reduced to partial differential equations when the acoustic wavelength is much smaller than L [9–11]. This simplification is not used in this paper because the acoustic wavelength is of the same order as the plate length, L. Therefore, a set of integro-differential equations is integrated to obtain the panel oscillation and then using the integral representation the acoustic field is calculated. An extension of this work to the problem of noise transmission in a supersonic flight at high altitude is given in Frendi et al. [12]. Note that the solution of the decoupled panel oscillation is an unsteady two-dimensional problem in the finite domain D and is several orders of magnitude simpler than the solution of the fully coupled problem. To show the efficiency of the decoupled model and the accuracy of the solution even when the panel oscillation is non-linear, we compare the solution of this decoupled system with the solution of the fully coupled problem for which we solve the three-dimensional non-linear Euler equations for the acoustic field coupled with the non-linear equations for the panel oscillations. In the next section, we present the mathematical formulations of these two models. The various numerical techniques used to solve the problems involved are described in section 3. The results and discussion are given in section 4, and the conclusions are in section 5.
2. FORMULATION OF THE TWO MODELS
The physical problem being studied is that of linear and non-linear oscillations of a typical aircraft panel excited by harmonic plane waves at normal incidence. The panel is clamped on to a large rigid plate. Let the typical panel be represented by a rectangular domain D in the (x, z) plane, with length L, width W and one vertex located at (x0 , z0 ), i.e., D = {(x, z) = x0 Q x Q x0 + L,
z0 Q z Q z0 + W}.
(4)
. .
744
The complement of D in the (x, z) plane represents the rigid plate. The incident domain is the half-space y q 0 and the transmitted domain is the half-space y Q 0. Both domains have the same ambient pressure pa and temperature Ta, and hence the same speed of sound C and density ra . We denote the pressure and velocity potential of the acoustic field by p and f and use the superscripts + and − to denote the quantities on the sides y q 0 and y Q 0. On the side y q 0, the superscripts i and r are used to denote the quantities associated with the incident and reflected waves, respectively. These two waves are mirror images with respect to the (x, z) plane. When the incident waves are plane waves advancing in the direction opposite to the y-axis and hitting the (x, z) plane at t = 0, we have f (i)(t, x, y, z) = f (Ct + y)
f (r)(t, x, y, z) = f (Ct − y),
and
(5)
with f(j) = 0 for j E 0. Therefore f (r) = 0 when 0 e y Q Ct. Transverse oscillation, h, of the rectangular flexible panel, D, is excited by the pressure difference across the panel for t q 0 and the oscillation in turn induces scattered waves f s in the incident side and transmitted waves f t. Both waves have the homogeneous initial conditions at t = 0, f (s) = ft(s) = 0,
y q 0,
f (t) = ft(t) = 0,
y Q 0.
(6)
Under linear theory, we have f (r)(t, x, y, z) = f (i)(t, x, −y, z),
f (t)(t, x, −y, z) = −f (s)(t, x, y, z),
for y q 0.
(7)
The wavefronts of f (s) and f (t) are the envelopes of the sonic spheres with radius Ct and centered on the panel in y q 0, and y Q 0, respectively. The pressures in the incident and transmitted sides are p+ − pa = p (i) + p (r) + p (s),
y q 0,
p− − pa = p (t),
y Q 0,
(8)
and the pressure difference across the panel is Dp = p− = p+ = 2p (t) − [p (i) + p (r)],
y = 0,
(x, z) $ D.
(9)
In the next sections we formulate two mathematical models to compare with and complement the experimental studies. We present the system of non-linear equations governing the acoustic fluid and panel oscillation (the fully coupled model) in section 2.1 and formulate the system of integro-differential equations for non-linear panel oscillation, including the effect of a linear acoustic field (the decoupled model), in Section 2.2. 2.1. In this model, the structural–acoustic interaction is analyzed by solving the three-dimensional non-linear Euler equations together with the non-linear plate equations. The configuration of the computational domain is shown in Figure 1. In Cartesian co-ordinates, x, y and z, the compressible, non-linear Euler equations can be written in conservation form as Qt = Fx + Gy + Hz ,
(10)
where Q is the vector (r, ru, rv, rw, e), r is the density, ru, rv and rw are the x, y and z momenta respectively, and e is the total energy per unit volume, given by e = 12 r(u 2 + v 2 + w 2) + rcv T,
(11)
with cv being the specific heat at constant volume. In equation (10), the functions F, G and H are:
–
F J G ru G 2 ru + p G G F =G ruv G, G ruw G u(e + p)G G f j
F G rv J G G ruv G G =G rv 2 + p G, G rvw G v(e + p)G G f j
745
F G rw J G G ruw G H = G rvw G. 2 G rw + p G w(e + p)G G f j
(12)
In addition to equation (10), the equation of state of an ideal gas is used: p = rRT,
(13)
where p is the pressure, R the gas constant, and T the temperature. Since we assume that the incident wave hits the panel at t = 0, the flow on the transmitted side, y Q 0, is at rest for t E 0. The initial data at t = 0 for y Q 0 are u = v = w = 0,
p = pa
and
r = ra .
(14)
On the rigid plate, we have zero normal velocity, v(x, 0, z) = 0
for
(x, z) ( D.
(15)
The motion of the flexible panel, D, is described by a system of three, non-linear partial differential equations, given by reference [13]: D9 4h + rp hhtt + ght = Dp +
Eh [(u 0 + 1 h 2 )(h + nhzz ) 1 − n 2 x 2 x xx
+(wz0 + 12 hz2 )(hzz + nhxx ) + (1 − n)hxz (uz0 + wx0 + hx hz )],
(16)
0 0 uxx + d1 uzz0 + d2 wxx = −hx (hxx + d1 hzz ) − d2 hz hxz , 0 0 + d2 uxz = −hz (hzz + d1 hxx ) − d2 hx hxz , wzz0 + d1 wxx
where 9 4h = hxxxx + 2hxxzz + hzzzz , d1 = (1 − n)/2,
d2 = (1 + n)/2,
(17)
D = Eh 3/12(1 − n 2),
(18)
u and w are the in-plane displacements, and h is the transverse displacement. The physical constants of the panel appearing in equation (16) are the stiffness (D), the density (rp ), the thickness (h), the physical damping (g), the modulus of elasticity (E) and the Poisson ratio (n). The system of equations (16) is solved subject to the homogeneous initial conditions 0
0
at t = 0,
u 0 = w 0 = h = ht = 0,
(19)
and the clamped boundary conditions on C given by x = x0 , x0 + L,
u 0 = w 0 = h = hx = 0,
z = z0 , z0 + W,
u 0 = w 0 = h = hz = 0.
(20)
In equation (16), the load Dp defined by equation (9) contains the coupling with the acoustic field p t and the forcing term (p (i) + p (r) ), which represents the load on the panel if the panel were rigid. We assume the forcing term to be a harmonic wave of the form p (i) + p (r) = (e sin(vt) + e* sin(2vt))H(t)
on y = 0+,
(21)
where e, e* and v are the amplitudes and frequency of the wave and H(t) denotes the Heaviside unit step function. For small values of e, e* = 0 (linear), while for large values
. .
746
of e, e* = e/2 (non-linear). This input is used based on the knowledge of a typical experimental input. The load on the panel D given by equation (9) becomes Dp = [2p t − e sin(vt) − e* sin(2vt)]H(t).
(22)
Another condition coupling the acoustic fields and the panel oscillation is the kinematic condition, ht (t, x, z) = v2(t, x, 02, z),
(x, z) $ D.
(23)
2
It is imposed on y = 0 because of the small panel displacement, equation (1). The numerical scheme for the solution of the non-linear plate equation (16) and the Euler equations (10) for the transmitted waves with homogeneous initial data will be described in section 3. 2.2. Now we analyze the non-linear panel oscillation excited by a weak incident pressure wave, p (i), under the assumption that the pressure fluctuation remains much smaller than the ambient pressure pa , i.e., (p − pa )/pa W 1.
(24)
Consequently, the small disturbance theory is applicable to the acoustic fields in the incident and transmitted sides and the panel/acoustic interaction problem is described by equations (4)–(9). The velocity potential F2(t, x, y, z) is governed by the simple wave equation, 2 2 (C−21tt2 − 1xx − 1yy − 1zz2 )F2 = 0,
for 2y q 0,
(25)
and the acoustic pressure and velocity v2 are related to the potential by p2 = −pa 1t F2
v2 = 9F2.
and
(26)
In particular, we have v = 1y F, where v denotes the transverse velocity. The transverse displacement of the panel, h(t, x, z), is governed by the system of partial differential equations (16). To produce the special forcing term, (p (i) + p (r) ), acting on the panel specified by equation (21) in section 2.1, the incident potential should be f (i)(t, x, y, z) =
0 $ 0 1%
e y cos v t + 2ra v C
+2
$ 0 1%1 0 1
e* y cos 2v t + e C
H t+
y . (27) C
Note that the incident wave f (i) is a solution of equation (25) in the whole space. When the panel is rigid, the incident wave is reflected by the rigid (x, z) plane and the reflected wave is f (r)(t, x, y, z) = f (i)(t, x, −y, z)
for y e 0.
(28)
On the side y Q 0, the acoustic fluid remains at rest because there is no transmitted wave. For a flexible panel, the oscillation by the pressure difference across the panel −2p (i)(t, x, 0+, z) induces a scattered wave f (s) on the incident side and a transmitted wave f (t). Thus we write F = f (i) + f (r) + f (s),
y q 0,
(29)
and F = f (t),
y Q 0.
(30)
–
747
Note that f and f are governed by the simple wave equation (25) in y q 0 and y Q 0, respectively. The kinematic conditions on the (x, z) plane, equations (23) and (26), become (s)
(t)
1y f (s)(t, x, 0+, z) = 1t h(t, x, z)
1y f (t)(t, x, 0−, z) = 1t h(t, x, z),
and
(31)
where h denotes the extension of the transverse displacement of the panel, D, to the (x, z) plane, i.e., h = 0,
(x, z) ( D.
(32)
If the incident wavefront hits the panel at t = 0, we can impose the homogeneous initial conditions on f (s), f (t) and h, f (s) = 0,
1t f (s) = 0,
f (t) = 0,
1t f (t) = 0
and h = 0,
1t h = 0,
for t E 0.
(33)
Since the ambient fluid above and below the panel are the same, equation (31) implies that the velocity potential induced by the panel oscillation has to be antisymmetric in y, i.e., f (t)(t, x, y, z) = −f (s)(t, x, −y, z),
u q 0.
(34)
Then, the pressure difference across can be written as Dp = 2ra [1t f (i)(t, x, 0+, z) + 1t f (s)(t, x, 0+, z)],
(x, z) $ D.
(35)
The derivation of the system of equations (24)–(35) completes the formulation of the structural–acoustic interaction problem using linear theory for the acoustic fields. The panel oscillation, which can be non-linear, is governed by the system of equations (16)–(20). The velocity potential f (s) induced by the panel oscillation h is governed by the wave equation (25), the initial conditions (33) and the boundary condition (31). The solution is given by the Kirchhoff formula, 1 f (s)(t, x, y q 0, z) = − 2p
gg
{ht (t, x', z')} dx' dz', R G
(36)
where R = [(x − x')2 + y 2 + (z − z')2]1/2 denotes the distance from a point (x, y q 0, z) to a source at (x', 0, z') and {·} denotes the retarded value, i.e., {ht (t, x', 0, z')} = ht (t − R/C, x', 0, z').
(37)
The domain of dependence of f (s)(t, x, y, z) is the circular disc H in the (x', z') plane, i.e., H= R E Ct
or
r 2 = (x' − x)2 + (z' − x)2 E C 2t 2 − y 2.
(38)
The domain of integration in equation (36) is the intersection of H and the panel, i.e., G=H + D
(39)
Now we introduce the polar co-ordinates, r, u, centered at (x, z), i.e., x' − x = r cos u,
z' − z = r sin u,
(40)
and equation (36) becomes, 1 2p
f (s)(t, x, y q 0, z) = −
g$g 0
dugt t −
G
R , r, u C
1%
r dr , R
(41)
with g(t, r, u) = h(t, x + r cos u, z + r sin u).
(42)
748 . . Here, R = (r 2 + y 2)1/2 represents the slant height of a circular cone with vertex P(x, y, z), a vertical axis and a base circle of radius r in the (x, z) plane. On account of equation (32), we can extend the domain of integration of H and rewrite equation (41) as an iterated integral, 1 2p
f (s)(t, x, y q 0, z) = −
g $g Ct
0
2p
0
dugt t −
0
R , r, u C
1%
r dr . R
(43)
The integral in u represents the contribution of the sources on a base circle of radius r. In equation (35), we need to relate the unknown, ft(s) (x, 0+, z) on the panel, to h(t, x, y) for (x, z) $ D. This is obtained by differentiating the above equation with respect to t and using equation (32). This results in replacing f (s) on the left side of equation (43) and gt in the integrand by ft(s) and gtt , respectively. As y :0+, R :r and the circular disc H is bounded by the sonic circle r = Ct. Equation (43) yields ft(s) (t, x, 0+, z) = −
1 2p
gg 0 2p
Ct
gtt t −
0
0
1
r , r, u dr du. C
(44)
Note that in this form we remove the kernel 1/R in equation (36), which becomes singular as y :0 and r :0 and shows that the area of the domain of integration is bounded above by 2pCt. By using equations (35) and (44), equation (16) becomes a system of integro-differential equations for the panel oscillation, h. The initial and boundary conditions for h are equations (19) and (20). Thus, we complete the formulation of the decoupled model. The numerical solution of this system will be described in section 3. 3. NUMERICAL METHODS
For the fully coupled model, the non-linear Euler equations (10) are solved using an explicit finite difference scheme. The scheme, which is a generalization of MacCormack’s scheme obtained by Gottlieb and Turkel [14], is fourth order accurate in space and second order accurate in time. Further details on the implementation of the scheme can be found in Frendi et al. [1–3]. The physical boundary of the computational domain (the bottom boundary; see Figure 1) is composed of a flexible panel clamped between rigid plates. The boundary conditions employed on both the flexible and rigid plates for the Euler equations are v = 0,
T = Tw ,
(45)
over the rigid plates, and v = 1h/1t,
T = Tw ,
(46)
over the flexible panel. In equations (45) and (46), Tw is a specified wall temperature; in this paper it is taken to be the free-stream temperature, Ta . Over the rigid plates, the x and z components of the velocity (u and w, respectively) are obtained through linear extrapolation from the interior of the computational domain. The pressure boundary condition is obtained using the normal momentum equation by simply imposing the normal gradient of the sum of pressure and vertical momentum flux to be zero, i.e., 1y (p + rv 2) = 0. Over the flexible panel, the x and z velocity components and the pressure are extrapolated from the interior of the computational domain. The remaining numerical boundary conditions of the computational domain are derived using the method of characteristics [15]. One should mention that the characteristic boundary conditions and the extrapolation are only first order accurate, since the information on the boundary
–
749
depends on the point next to the boundary only, while the interior scheme is fourth order. These boundaries are believed to be a source for numerical error in the fully coupled model. The truncation error in the interior of the domain is of order (Dx)4, while near the boundary it is of order (Dx); this jump in truncation error is a major source for the spurious waves that propagate back into the domain. The numerical error can show up in the results in the form of additional non-physical frequencies. The dimension of the computational domain is 3.02 m in the x, y and z directions, respectively. The number of computational points used are 141, 181, and 141 in the respective directions. The non-linear plate equations (16) are solved using a finite element method developed by Robinson [16]. The panel is 30·5 cm long, 20·32 cm wide and 0·102 cm thick, and the number of elements used are six and eight in the length and width, respectively. Since the grid used for solving the plate equations is rectangular, it is easier to evaluate the integral in equation (44) in Cartesian co-ordinates: ft(2) (t, x, 0+, z) = −
1 2p
g g 0 Xmax
Xmin
Zmax
Zmin
htt t −
1
dx' dz' r , , x', 0, z' r C
(47)
where (Xmin , Zmin ) and (Xmax , Zmax ) lie within the borders of the smallest and largest circles on the panel, respectively. The radius of the smallest circle is the smaller of Dx or Dz, and that of the largest circle is R = Ct. The presence of the small circle results from the singularity at r = 0, which appears in the integral (47). The contribution to the pressure field of the point at r = 0, (Xp , Zp ), is found by expanding htt (t − r/C, x', 0, z') in a Taylor
Figure 2. (a) The time history and (b) power spectral density of the panel center displacement for an excitation amplitude of 130 dB (linear). ——, Fully coupled model; . . . . . . , decoupled model.
750
. .
Figure 3. (a) The time history and (b) power spectral density of the transmitted near field pressure, 2·54 cm from the panel center, excitation amplitude 130 dB (linear). ——, Fully coupled model; . . . . . . , decoupled model.
series. The expansion is truncated after the third order term. The contribution of the leading term in the expansion is htt (t, Xp , 0, Zp )R1 Du, where R1 is equal to the smaller of Dx or Dz and Du is equal to 2p for an interior point, p for a boundary point and p/4 for a corner point. The contribution of the various points on the panel that lie within the sonic circle (R E CT) is calculated by first integrating in x using a combination of Simpson and trapezoidal rules. The result is then integrated in z using Simpson’s rule of integration. Because of the presence of the retarded time in the integral, the acceleration of the plate (ht t) is stored at each point for several time-steps. The evaluation of the double integral can use a much larger time step than the Dt for the integration of the differential equations, yet having the same degree of accuracy. For a plate of given dimensions (L, W), and for a fixed time-step Dt, the maximum number of time-steps to be stored is N = zL 2 + W 2/(CDt). The number N has to be changed when calculating the radiated pressure away from the plate. The number N is related to the radius of the largest sonic circle in the domain that contains the flexible panel. For the cases presented in the next section, N = 1503 for a time-step Dt = 4 × 10−6 s. 4. RESULTS AND DISCUSSION
The numerical schemes for the two models presented in section 3 are used to predict the vibration of a flexible panel and the resulting acoustic radiation. The panel is forced
–
751
to vibrate by harmonic plane acoustic waves at normal incidence. The frequency of the source is 751 Hz for low amplitude excitations, and a combination of 751 Hz and 1502 Hz for high amplitude excitations. This method of input was chosen based on the knowledge of the experimental input, however, no direct comparisons with experiments are made in this paper. The frequency 751 Hz corresponds to a natural frequency of the flexible panel. Two different amplitudes of the incident waves are used. The properties of the panel, which are considered to be uniform, are density rp = 4450·15 kg/m3, modulus of elasticity E = 1·10316 × 1011 N/m2, Poisson ratio n = 0·33 and a damping ratio of 0·01 is used. The acoustic fluid properties are temperature Ta = 288·33 K, density ra = 1·23 kg/m3, pressure ra = 1·013 × 105 N/m2 and sound speed C = 340 m/s. The specific heat at constant volume is cv = 1·004 kJ/(kg K), the ratio of specific heats is g = cp /cv = 1·4. In the far field, the fluid is at rest. The variables plotted in Figures 2–5 are non-dimensional. The reference quantities are given by (x, y, z, h)ref = lref , rref = ra ,
tref = lref /C,
(u, v, w, 1h/1t)ref = C,
Tref = C 2/cv , (p, e)ref = ra C 2,
(48)
where the reference length is lref = 0·3048 m. The configuration of the computational domain, a square box with the lower side composed of a flexible panel clamped between rigid plates, is shown in Figure 1. For a low excitation amplitude, 130 dB or 6·8 × 10−4 atm, the panel response is linear as shown
Figure 4. (a) The time history and (b) power spectral density of the transmitted far field pressure, 203·2 cm from the panel center, excitation amplitude 130 dB (linear). ——, Fully coupled model; . . . . . . , decoupled model.
752
. .
Figure 5. (a) The time history and (b) power spectral density of the strain response of the panel for an excitation amplitude of 174 dB (non-linear). ——, Fully coupled model; . . . . . . , decoupled model.
in Figures 2(a) and (b). The time history of the normalized panel center displacement (h/h, where h is the panel thickness) shows a periodic behavior with one dominant frequency f (see Figure 2(a)). In this case, both models predict the same panel response as shown by the figure. This result is further confirmed by the power spectral density of the displacement, Figure 2(b), which shows a strong spike at f = 751 Hz. The near field transmitted pressure, 2·54 cm from the panel center or L/12, is shown in Figures 3(a) and (b). The time history of the transmitted near field pressure, Figure 3(a), shows a periodic behavior with one dominant frequency. The power spectrum of this time signal shows a strong spike at the fundamental frequency f = 751 Hz (see Figure 3(b)). Both models predict the same near field pressure. A small wiggle is shown by the power spectrum of the decoupled model near the fundamental, its magnitude is negligible compared to that of the fundamental and therefore has no effect on the results as is shown in Figure 3(a). This wiggle is believed to be caused by the difference in the length of the time samples used to calculate the power spectrum. In Figures 4(a) and (b) is shown the transmitted far field pressure, 203·2 cm from the panel center or 5L. Similar to the near field, the far field pressure is periodic with one dominant frequency f = 751 Hz as indicated by both the time history (Figure 4(a)) and the power spectrum (Figure 4(b)). The two figures also show that both models predict the same far field pressure. When the level of the excitation is increased to 174 dB, or 0·22 atm, the response of the panel becomes non-linear. The strain response of the panel at the middle of the long edge is shown in Figures 5(a) and (b). The time history of the strain, Figure 5(a), shows a
–
753
periodic response with a period 4T corresponding to the frequency f/4, and a strong harmonic 2f characterized by the break-up of the positive peaks. This observation is further confirmed by the power spectral density of the strain, Figure 5(b), which shows the presence of several frequencies. Most of the frequencies are either harmonics or subharmonics of the fundamental frequency f. This spectrum is characteristic of a non-linear panel response. The panel response obtained by the ‘‘decoupled model’’ is in fairly good agreement with that obtained by the ‘‘fully coupled model’’. In particular, there is a good agreement in the power spectra of the response predicted by the two models. The predicted panel center displacement is shown in Figures 6(a) and (b). The time history of the normalized panel center displacement is shown on Figure 6(a). Contrary to the strain response, the displacement response shows less non-linearity. In Figure 6(a) is shown a dominant fundamental frequency f with a weak subharmonic f/2. This is further confirmed by the power spectrum of the displacement, which shows a weak harmonic 2f in addition to f and f/2. The predictions of the two models are in fairly good agreement. Since the coupling between the acoustic fluid and structure is obtained through the out-of-plane velocity of the panel, a comparison between the time histories of the non-dimensional panel velocities at the panel center is shown in Figure 7. As expected, the centerline panel velocity shows a more pronounced effect of the harmonic 2f. This is particularly apparent in the bottom part of the time histories. The predictions of both models are in fairly good agreement.
Figure 6. (a) The time history and (b) power spectral density of the panel center displacement for an excitation amplitude of 174 dB (non-linear). ——, Fully coupled model; . . . . . . , decoupled model.
754
. .
Figure 7. The time history of the panel center out-of-plane velocity for an excitation amplitude of 174 dB. ——, fully coupled model; . . . . . . , decoupled model.
In Figures 8(a) and (b) is shown the transmitted near field pressure, 2·54 cm form the panel center or L/12. The time history, Figure 8(a), shows a periodic behavior with a dominant fundamental f and a strong harmonic 2f. The effect of the harmonic is noticeable
Figure 8. (a) The time history and (b) power spectral density of the transmitted near field pressure, 2·54 cm from the panel center, excitation amplitude 174 dB (non-linear). ——, Fully coupled model; . . . . . . , decoupled model.
–
755
Figure 9. (a) The time history and (b) power spectral density of the transmitted far field pressure, 203·2 cm from the panel center, excitation amplitude 174 dB (non-linear). ——, Fully coupled model; . . . . . . , decoupled model.
in the break-up of the positive peaks of the time history. The power spectrum of the near field pressure, Figure 8(b), shows a weak subharmonic f/2 with a strong fundamental f and harmonic 2f. Notice that the difference in level between f and 2f is about 10 dB. There is a fairly good agreement between the predictions of the two models. The transmitted far field pressure, 203·2 cm from the center of the panel, is shown in Figures 9(a) and (b). The time history of the far field pressure shows a periodic behavior with period T/2 corresponding to 2f. It is important to note that in the near field (Figure 8(a)), the fundamental frequency f is dominant; however, in the far field (Figure 9(a)) it is the harmonic 2f that is dominant. The power spectral density of the far field pressure, Figure 9(b), shows that the fundamental and the harmonic 2f have nearly the same level. This indicates that the two frequencies do not decay at the same rate into the far field. This behavior is due to the fact that the panel is not a compact source but rather a set of distributed sources. The predictions of the two models are in farily good agreement. The level of 2f predicted by the fully coupled model is 3 dB higher than that predicted by the decoupled model. This is believed to be due to a combination of numerical errors and weakly non-linear wave propagation effects. From a computational view point, it is important to compare the performance of the two models based on the CPU time required by each calculation. In the linear vibration regime, the ‘‘fully coupled model’’ used 36 000 s of CPU time on a Cray-ymp to advance the calculation by 10 000 time-steps, whereas the ‘‘decoupled model’’ used only 1000 s for
756
. .
the same calculation. In the non-linear vibration regime, grid refinements were needed to resolve the large gradients both on the panel and in the radiation field. Therefore, in order to advance the calculation by 10 000 time-steps, 72 000 s were used by the ‘‘fully coupled model’’ while the ‘‘decoupled model’’ used only 1200 s.
5. CONCLUSIONS
An efficient model for coupling the vibration of a panel to the on surface acoustic radiation is derived. The model uncouples the panel vibrations from the acoustic wave propagation problem. The results show that this model, referred to as ‘‘decoupled model’’, accurately predicts the panel response and acoustic radiation in the linear and non-linear vibration regimes as long as the pressure variation in the acoustic fluid remains much smaller than the ambient pressure. For the cases studied in this paper, the computational cost of the numerical integration of this model is 36 times cheaper in the linear regime and 60 times cheaper in the non-linear regime than the cost of the ‘‘fully coupled model’’.
ACKNOWLEDGMENTS
The first author acknowledges the support of NASA Langley Research Center under contract NAS1-19700. Partial support was also provided by NASA Langley to the third author under contracts NAS1-18605 and NAS1-19480 while in residence at ICASE.
REFERENCES 1. A. F and J. R 1993 American Institute of Aeronautics and Astronautics Journal 31, 1992–1997. On the effect of acoustic coupling on random and harmonic plate vibrations. 2. A. F, L. M and A. B 1994 Journal of Sound and Vibration 177, 207–226. Coupling between plate vibration and acoustic radiation. 3. A. F, L. M, J. R and A. B 1993. American Institute of Aeronautics and Astronautics Paper 93–4367. On acoustic radiation from a vibrating panel. 4. L. M, A. F and D. E. B 1992. American Institute of Aeronautics and Astronautics Journal 30, 2632–2638. Nonlinear vibration and radiation from a panel with transition to chaos induced by acoustic waves. 5. P. M. M and K. U. I 1968 Theoretical Acoustics. New York: McGraw-Hill. 6. G. A. K and C. M 1980 SIAM Journal for Scientific and Statistical Computing 1, 371–385. Solving the Helmholtz equation for exterior problems with variable index of refraction: I. 7. A. B and E. T 1980 Communications in Pure and Applied Mathematics 23, 707–725. Radiation boundary conditions for wave-like equations. 8. L. T and M. J. M 1986 Journal of the Acoustical Society of America 80, 1825–1827. Exact boundary conditions for scattering problems. 9. M. J. M and L. T 1989 Wave Motion 11, 545–557. Scattering of an incident wave from an interface separating two fluids. 10. M. J. M and L. T 1992 Northwestern University Report 9130. Scattering of acoustic wave with a flexible wall. 11. L. T 1993. Paper presented at a workshop on Perturbation Methods in Physical Mathematics, Renselear Polytechnic Institute, Troy, New York. On surface conditions for structural acoustic interactions in moving media. 12. A. F, L. M and L. T 1993. Institute for Computer Applications in Science and Engineering (ICASE) , Report No. 93–18. An efficient model for coupling structural vibrations with acoustic radiation. 13. C. Y. C 1990 Nonlinear Analysis of Plates. New York: McGraw-Hill. 14. D. G and E. T 1976 Mathematics of Computation 30, 703–723. Dissipative two-four methods for time dependent problems.
–
757
15. S. S. A, W. S. D, D. G, D. H. R and J. C. T 1990 Institute for Computer Applications in Science and Engineering (ICASE) , Report 90-16. Secondary frequencies in the wake of a circular cylinder with vortex shedding. 16. J. H. R 1990 Master’s Thesis, Mechanical Engineering Department, Old Dominion University, Norfolk, Virginia. Finite element formulation and numerical simulation of the random response of composite plate.