- Email: [email protected]
Coupled fluid-structure interaction analysis
Finite Elements in Analysis and Design 7 (1991) 331-342
331
Elsevier
Coupled fluid-structure interaction analysis Kamran Izadpanah, Robert L. Harder, Raj Kansakar and Mike Reymond MacNeal-Schwendler Corporation, 815 Colorado Blvd., Los Angeles, CA 90041, USA Revised June 1990 Abstract. A general-purpose fluid-structure interaction capability which can be used in analyzing internal acoustic problems in automobiles and aircraft has been formulated. This capability has become an important issue in the last few years because of the emphasis that has been placed on noise control in all vehicles. A variety of output features that are needed by the analyst in identifying the causes of the acoustic problems have been studied and implemented. The fluid (air) is modeled as an inviscid, irrotational fluid with pressure variables as the degrees of freedom. An acoustic absorber element and an acoustic barrier element are also implemented.
Introduction An internal acoustic analysis capability [1] which is fully automated and requires minimum user input has become an important tool in analyzing internal noise problems in automobiles, trucks, submarines and airplanes. In the automotive industry for example, this type of analysis can help the understanding of vehicle noise which is usually referred to as booming noise or road noise occurring in a frequency range of up to 500 Hz. The two major causes of this noise are the engine and the various structural panels which are excited by the suspension system. With the methodology outlined in this paper and implemented in Version 67 of MSC/ NASTRAN, these problems can be identified and corrected before the actual structure is built and tested. The analyst can perform eigenvalue analysis, frequency response analysis (direct, modal and mixed, i.e., modal fluid, direct structure and vice versa), and, transient response analysis (direct, modal and mixed, i.e., modal fluid, direct structure and vice versa). The fluid-structure coupling terms are automatically calculated for both matching and non-matching interface meshes. In addition, two new elements have been formulated and added to the element library of MSC//NASTRAN which can act as acoustic absorbers and acoustic barriers and are useful when analyzing acoustic noise control devices.
Theoretical background The governing differential equations describing compressible fluid flow are the Euler equation
p[(v.V)v+e] = - r e , and the continuityequation
ap/Ot + PoV " v =
O,
where p is density, v is velocity, P is pressure, and S7 = (O/ax)i + ( O / ~ y ) j + (O/Oz)k. Elsevier Science Publishers B.V.
(1) (2)
K. Izadpanah et al. / Coupled fluid-structure interaction analysis
332
For acoustic analysis, three assumptions are made at this stage: • small motion theory, i.e., v---~; • the convective m o m e n t u m terms (first terms of the Euler equation) are negligible, i.e., p/i = - VP; • the P-p behavior is locally linear, i.e., ( d P / d p ) = YPo/Po where y is the ratio of specific heats and c = ~ is the speed of sound. Assuming an homogeneous fluid, the continuity equation becomes
(dp/dP)(dP/dt)
+ Po V . v = O.
(3)
Using the third assumption and integrating, this equation will yield P = -B V-u,
(4)
where B = c02po is called the "coefficient of compressibility". This relation is now differentiated twice with respect to time and X7 • ii is substituted from the continuity equation to yield the wave equation describing the fluid behavior which is expressed in the pressure variable only:
( l / B ) , ~ = ( 1 / p ) V . VP,
(5)
and the pressure and displacement variables are related by the modified Euler equation: On = - V P .
(6)
The finite element discretization of this equation can now be carried out using standard procedures. Assume: e = E P,N,,
(7)
where Pi and N~ are the pressure and shape functions associated with the ith node of an element. The weak form or the virtual work statement of cqn. (5) can be written as:
f v o l [ ( 1 / B ) P - (1/p)X7 " VP] d V S P = O .
(8)
F r o m the divergence theorem we have:
V " ( aV) = V a . V+ aV . V.
(9)
Applying this theorem to the second term of eqn. (8) yields the following (where a = 8P and V = VP): V • (SP WP) = V S P - V P + 8P V • V P ,
(10)
or
= -v
p. v e + v .
re).
(11)
However, VSP- VP
=
(12)
½5VP- V P ,
and therefore ~l'(V. re)
= - ½~(vP'vP)
+ V " (SJ' r e ) .
(13)
Also Pap
=
(14)
and fvol v • x dV = fsx. ds,
(15)
K. lzadpanah et al. / Coupled fluid-structure interaction analysis
333
where ds is an incremental surface vector (outward). Substituting these relationships in eqn. (8) yields:
LoIS[(1/2B)PZ+(1/2p)VP.VP)] dV- fs(1/p)vP.dsSP=O.
(16)
The first integral in eqn. (16) represents the virtual work statement of the differential equation (5). The second integral in eqn. (16) is the boundary condition associated with this problem. To obtain a valid solution, these two components have to be equal to zero for all admissible 8P. In order for the second integral to be zero for all admissible 8P, either 8P or v P " n has to be zero. Since 8P is an arbitrary admissible variation, it could be zero only when P is specified on the boundary. This is the essential boundary condition for this problem. If P = 0, the essential boundary condition corresponds to the conditions on a free surface. If P is not specified on the boundary then, V P " n or must be zero on the boundary. The is the natural boundary condition associated with this problem and it occurs at rigid walls. From eqn. (16), the last term in this relation can be written as:
OP/an
fs(ii, d s )
8P,
8Plan
(17)
or in matrix form as: [A](~},
(18)
Aij = fsNjNi ds.
(19)
where
Noting that the acceleration of the fluid is in a direction opposite to that of the structure and denoting the acceleration of the structure by (fi},, the problem is reduced to solving the following set of linear equations: [ M ] ( t i} + [ K ] { P }
+[A]{/i}s=0,
(20)
where
Mij = fvo(1/B)NdVj dr,
(21)
Kij = fvolVN/(]//p)vNj dV.
(22)
Fluid-structure interface
At a fluid-structure interface, the boundary condition that the pressure P must satisfy is:
aP/0n = - p / i ° ,
(23)
where n is the direction of the outward normal. This condition is reflected in the third term of
eqn. (20). On the other hand, the fluid affects the structure by applying forces over the structure surface area. These force vectors can be expressed as: (F} =-[A]T{P},
(24)
The A matrix defined in eqn. (19) is therefore the coupling matrix that needs to be computed in order to do coupled fluid-structure analysis. The calculation of this matrix for the general case
334
K. lzadpanah et al. / Coupledfluid-structure interaction analysis
is not a trivial matter. In fact a substantial amount of effort is necessary for general shapes and contact areas.
Coupled fluid-structure system The general equations of motion of the entire structure can now be written as: (25)
Mki + K x = F ,
where
0]
M--
,26,
Mf
and _A T
,27,
with the solution vector x as: x=
p
,
and the external load vector F=
~.
These relationships are used in a frequency response analysis. The eigenvalues (known to be real) of these unsymmetric real matrices may also be computed. Similarly, in a transient response analysis these matrices are used with the standard time integration procedures of MSC/NAST~N without any modifications. The inclusion of structural damping in this analysis does cause the eigenvalues of the system to be complex numbers. When the fluid is highly compressible, it is possible to symmetrize the K and M matrices by a suitable transformation of variables. This process involves inversion of the fluid matrices and thus is recommended only when these matrides are not large or when they have been replaced by their smaller modal representation. There are several approaches to this problem. The simplest approach is by substituting r for P where r=P-
[gfl-'[A]{u
}.
(30)
The resulting matrices for K and M are:
0
and F=
Fs ff"
(33)
K Izadpanah et al. / Coupled fluid-structure interaction analysis
335
Modal decomposition In very large problems, it is sometimes desirable to apply the modal decomposition method in order to reduce the size of the analysis. In this case the fluid, the structure, or both are possible candidates for this operation. Applying dynamic reduction to the structure, where the modes of the structure used for reduction are computed in vacuum results in:
u = .~s, [m.] = . T M . , ,
(34)
[k~] = . V g . s.
(36)
(35)
For instance, eqns. (35) and (36) substituted in eqns. (26), (27) and (29) result in
[m~ 0 ] M=A. ~ Mr'
(37)
K=[~'
(38)
-':AT]Kf '
and F - - [ "TF]
(39)
Lr,/• Applying dynamic reduction to the fluid (the modes of the fluid in a rigid container) results in: p = .f~f,
[mf ] =
(40)
"TM'f,
(41)
[kf ] = .fTK.f,
(42)
which when substituted in eqns. (26), (27) and (29) yields: M=
"~A
mf
and
-=[o AT° kf
(~)
J
with the same force vector as before (eqn. (29)). It is also possible to perform the dynamic reduction in both the fluid and the structure. The resulting stiffness and mass matrices are a combination of the ones derived above with the resulting modal mass and stiffness matrices:
M= ~ A . ,
mf'
K=[~ -'[AT'f]kf"
(46)
The equation of motion becomes
ms "?'~'s
0 mf
~'~ + ~f
--
kf
] t ~f J
. ['fTFf I
(47)
K. lzadpanah et a L / Coupled fluid-structure interaction analysis
336
Note that the modal analysis method is a useful method if the analyst is fully aware of the reduction process that is being carried out. It is possible to conceive of certain problems where the results obtained by an inexperienced analyst will be completely unacceptable. In that case it is recommended that a direct analysis be performed. On the other hand, the direct method of solution involves either sparse but unsymmetric matrices or dense but symmetric matrices. This means that for certain type of problems the computational effort involved might be significantly more than for a correctly applied modal reduction method.
Modal and boundary participation factors The modal solution method provides valuable information regarding the structural and modal participation in the acoustic response [5]. The acoustic modal participation at the interior points ( p~ } can be computed by expanding eqn. (40) as: L
{p} =
= E
(48)
i=I
where L is the total number of acoustic modes used in the analysis. Each vector (p~ } then represents the contribution to the acoustic response from the acoustic mode i or the modal acoustic participation factor. The structural modal participation factor can be computed for a harmonic solution of radian frequency o~ as follows. First the equation of motion is solved for the fluid modal vector: {if } =
[z2](~o2[a]{t2~}),
(49)
where [a] = [*f]T[A][~s]
(50)
is the modal coupling matrix and [z21 = ( - , : [ m f ]
+ kf)-'
(51)
is the acoustic modal frequency response function. Substituting eqn. (49) into eqn. (40) results in: SM
(p}=J[*t][z2][a](~s}=
•
(p:},
(52)
j=l where S M is the total number of structural modes. Note that each vector { p: } represents the contribution of the structural mode j to the the acoustic response or the structural modal participation factor. The boundary panel participationfor a boundary panel b is basically the collectivestructure modal participation of all the structural nodes on that panel. The [Ab] matrix which is the coupling matrix due to panel b therefore has to be computed. The global [A] matrix will be the assembly of all the panels in the structure. Once the [Ab] are known, eqn. (50) can be written as: N
(53)
[a] = [q~flT y , [ A e ] [ ~ , l , b=l
which when substituted in eqn. (52) becomes: N
{p}
=~,2[~fl[z21[*flT E
b=l
N
= Y'. ( p b } , b=l
(54)
K. lzadpanah et al. / Coupled fluid-structure interaction analysis
337
where N is the total number of boundary panels. Note that each { Pb } vector represents the contribution of boundary panel b to the acoustic response.
Simple acoustic source
A "simple" acoustic source is defined as an enclosed surface vibrating with arbitrary velocity distribution but is of such size that all dimensions are much smaller than the wavelength of the emitted sound [23]. It can be shown that the pressure field of any irregular simple source m a y be approximated by the pressure field produced by a pulsating sphere of the same source strength. The source strength Q, is defined as the rate of flow of the fluid away from the source given by: Q e j'°' =
fsh.n
dS,
(55)
where ~i is the velocity vector and n is the outward normal. For a sphere pulsating at a given frequency to, the total power can be calculated as:
1I = ½~rOoc(Q/?~ )2,
(56)
where 11 is the total power, c is the speed of sound in the fluid, and h = 2~rc/to is the wavelength. This expression can be inverted to express the source strength as a function of total power H:
1 / 8~rcH O = -d
(57)
p
This is the expression that is used in MSC//NASTRANto convert the power versus frequency curve supplied for a given source into the right hand vector (loading) for the pressure variables. Note that whenever a phase angle ~b is specified by the user, this load vector should be viewed as the magnitude of a complex loading vector hating an angle ~. This complex vector is decomposed into real and imaginary load vectors.
Acoustic barrier
An acoustic barrier is the component that is usually placed between two compartments such as the barrier between the engine compartment and the passenger compartment. The purpose of the barrier is to reduce the noise transmission between the two compartments. Even though there are m a n y different types of barriers, the most common construction is a sandwich consisting of a steel backing, a foamy middle layer, and a plastic type airtight cover usually called the septum. The transmission phenomenon is a very complicated physical phenomenon. In the adopted approach, this behavior is simplified by assuming that the incident wave is normal to the surface and that there is no in-plane (parallel to the sandwich) transmission of energy. In other words, the barrier model constructed here is a relatively crude model that gives reasonable results in one-dimensional cases and should be used with care in order to avoid erroneous conclusions. There is a standard (ASTM) test to measure the transmission loss (TL) of any barrier defined as T L = 20 log( Pi/Pt
),
(58)
K. hadpanah et al. / Coupled fluid-structure interaction analysis
338
where Pi is the incident pressure and Pt is the transmitted pressure. Note that some of the pressure wave hitting the receiving side of the barrier is reflected and that in the definition of TL only the incident pressure and not the total pressure is used. In a typical TL versus frequency curve, it is seen that the TL increases about 6 dB per octave in the low-frequency range and about 18 dB per octave in the high-frequency range. There is a transition region between these two ranges which is characterized by very low TL. In this region the sandwich construction is basically vibrating in resonance and the barrier is least effective. The "equivalent structural" system selected to represent the transmission phenomena consists of a spring sandwiched by two masses. The mass of the backing and of the septum is usually known. However, the stiffness of the spring is going to be calculated from the resonant frequency of the sandwich as observed in experiment. The theoretical resonant frequency of this model can be calculated as:
1 //¢
=
(59)
/Me'
where
M,M Me -- MI + 3,/2
(60)
and M 1, M 2 are the masses of the backing and the septum respectively. This means that the equivalent stiffness can be calculated by
M~M2(2~Ir.~) 2 K=
3,/1 + M 2
(61)
Once the effective stiffness t¢ is calculated, the stiffness and mass matrices of the barrier element can be calculated in the usual manner, i.e., for face 1 (steel): K U -- ~ f u , uj d a , ,
(62)
M . = M, f
(63)
d& ;
and for face 2 (septum): Kij =
f u,% dA2,
(64)
M,j = M2fm% dA2
(65)
The justification for using this equivalent structural system is that when it is subjected to the same loading and boundary conditions as the TL test, the results are close to experimental results. The pressure distribution in the T L test when the source is in the left chamber (x < 0) can be thought of as:
P =
Pi e - i ~ ( t - x / c ) + Pr e-i~(t+x/c),
x = 0-,
Pt e - i ~ ( t - x / c ) ,
x = 0 +,
(66)
where Pi is the incident pressure, P, is the reflected pressure, Pt is the transmitted pressure, c is the speed of sound, and i = ~-Z-i-. In all of the complex equations, only the real part is to be used.
K. lzadpanah et al. / Coupledfluid-structure interaction analysis
339
This pressure distribution is basically a right traveling incident wave, a left traveling reflected wave, and a right traveling transmitted wave. Note that because of the special rooms used in this test there are no left traveling waves on the back side of the sandwich. The boundary conditions on the left and right side of the sandwich are: pii + OP/Ox = 0,
(67)
which means:
fi( t i t
x0
/i(0-~)
x=0 +
u=
(68)
e -'°',
Now applying the equilibrium equations results in: i
( - ml., 2 + K) ~--/~ (P, - Pr) - K0.,c P, = Pi + Pr, i
--Kpa~c ( P i - Pr)
"
+ (3/2 ¢°2 + K ) ~ c Pt = - P t .
(69) (70)
These relations can be solved to obtain: Pt = - 2iKpwc Pi [ M1M2to4 - ( ( M 1 + M 2 ) K + (pc)2) oa2] + i [ ( M 1 + M2)Pcto 3 - 2Kptoc] ' (71) and, also: Pi+Pr Pt
ptoc(-Mlto2+K)+i[M1M2oj4-(MI+M2)Kto2] p~cK
(72)
With the proper choice of the resonant frequency, there is good agreement with test results.
Acoustic
absorber
An acoustic absorber is a component that is usually placed inside a compartment in order to minimize the reflection of sound. In an automobile, for example, the interior is padded with different types of absorbers depending on the available space and the acoustic objective in mind. An acoustic absorber is characterized by its specific acoustic impedance z defined as the ratio of the acoustic pressure in a medium to the associated particle speed, i.e., z = P / u . This is a complex number with the real part usually called the specific acoustic resistance and the imaginary part the specific acoustic reactance. These parameters are measured in a standing wave tube (ASTM standards) and have some general characteristics. As in the case of an acoustic barrier, the behavior of an acoustic absorber is modeled by structural components. Here a single mass in conjunction with a spring and damper in parallel is the "equivalent structural" model that is used. In the standing wave test, the absorbing material is secured to the end of the tube and a pressure wave of P sin a~t is applied. Therefore in the equivalent structural model, the spring and damper components are grounded and the equilibrium equation for this one-degree-offreedom system becomes: M g + CYc + K x = P sin ~0t.
(73)
340
K. Izadpanah et aL / Coupledfluid-structure interaction analysis
The acoustic impedance of this system is u
:~
C+
(
=t
~ 0 M - -~ 1
(74)
The test data are interpreted as follows: the resistance ~ ( z ) is averaged and used as the damper C. The reactance data .St(z) is used a least square fit to the curve
Ma~ - K/~o.
(75)
The resulting error function to be minimized for this least square fit is
(
~i Wi M°°i -- K __ Bi
;
(76)
where B i are the experimental data and W~ is a weighting function. The simultaneous equations that have to be solved for this problem are:
- E Wi
E "72Wii ] [ K J
- E "7'W, Bi
(77)
The summation can be carried out over an entire range of data or for the data corresponding to frequencies below a cutoff frequency. Solving for M and K results in: E o~Z2WiEc.oiWiB, - E WiEcaoZ1Wi "`
M
E-72 -
K
=
(78)
- (E
E +:w,E,+-'w,B, E E -,-- _ (E w,) +
2
(79)
The resulting M and K values can now be used in the structural model of the acoustic absorber.
Implementation
This methodology was implemented in the general-purpose finite element code MSC//NASTRAN Version 67 with the goal of having a very clean interface with existing capabilities. The fluid medium is described by the use of the standard GRID and element entries. The fluid grid points are distinguished from the structural points by a parameter on the GRID entry. The fluid elements that are available are in the shape of hexahedra, pentahedra or tetrahedra. Both linear (8-noded HEXA, etc.) and quadratic (20-noded HEXA, etc.) elements are supported. The topology (connectivity) of these elements is specified using the existing Bulk Data entries for these shapes, i.e., CHEXA, CPENTA, CTETRA. The property entry for the solid elements is the key in distinguishing between structural and pressure elements. By designating the fluid elements on this entry, the pressure-based fluid elements will be used. The bulk modulus and the mass density of the fluid have to be specified for these elements. It is sometimes necessary to find the contribution of a set of grid points to the noise level in an acoustic cavity. For example, in the automotive industry the entire roof assembly contributes significantly to the noise level inside the passenger cabin. Consequently it is necessary to compute and minimize the contribution of this "panel".
K. Iza@anah et aL / Coupled fluid-structure interaction analysis
341
The user is able to define panels by specifying the structural grid points in each panel using a new entry PANEL and convenient set definition tools. The acoustic absorber and barrier elements were implemented as new elements CHACAB and CrtACBR. These elements consist of 8 to 16 structural grid points, and define two surfaces. The property entry for the absorber elements points to three tables which contain the resistance (Y,), reactance (I12) and the confidence level (Y3) of the material as a function of frequency as measured in the standing wave tube test (ASTM standard). This data will be used to calculate the value of the mass, spring, and damper used in the structural model that will emulate the behavior of the absorber material. The property entry for the barrier contains the mass per unit area of the backing, mass per unit area of the septum, and the resonant frequency of the barrier. The two mass values are readily available and the resonant frequency of the sandwich is defined as the frequency with the lowest values of transmission loss (ASTM test). This frequency will be used to calculate the value of the spring used in the structural model that will emulate the behavior of the absorber material. The loading of the fluid is in most cases analogous to enforced displacements in structural mechanics analysis. There are several types of loading that have to be provided for the fluid elements: (1) Constant enforced pressure at the grid points; (2) Frequency or time-dependent enforced pressure at the grid points; (3) Acoustic source characterized by a volumetric flow rate corresponding to a power spectral density function. The large mass approach that is currently used to model enforce displacements is used to have enforced pressure loading. A simple source is modeled as a pulsating sphere in infinite space. The user can conveniently model this type of source by supplying the power versus frequency curve characterizing the acoustic source and the phase angle. The phase angle is useful whenever multiple sources are present. The single-point constraints of the fluid ( P = 0.0) can be enforced using the existing constraint entries. This type of boundary condition occurs at free surfaces. The identification of the fluid-structure interface region is accomplished by using a series of sorting operations. Ideally, the computations involving the fluid-structure interface should be completely invisible to the user. This is, however, a computation intensive process. In the current design the user is requested to "assist" the process in order to expedite the calculations. This user "assistance" is in the identification of the "wetted area" of the structure. Its only effect will be the reduction of the sorting operations carried out for the computations involved. To expedite the solution process, the user needs to supply: (1) Type of interface, i.e., matching or non-matching fluid-structure meshes; the default interface in the program is the non-matching interface; (2) Actual structural and fluid grid points that are in contact. The program always outputs a table containing the interface information. The capability is fully operational with the substructuring (superelement) capability of MSC//NASTRAN. In these superelement analyses only the following restrictions apply (SOL 107, 108, 109, 110, 111, and 112 of MSC/NASTRAN): • Each upstream superelement can contain some of the fluid or some of the structure but not both. • The fluid-structure interface grid points can only be in the residual superelement. This means that upstream superelements of different types will not contain the coupling terms (the A matrix). It is possible to perform modal reduction on the upstream fluid or structure superelements or
342
K. lzadpanah et al. / Coupled fluid- structure interaction analysis
both using the standard superelement case control capabilities in conjunction with multiple subcases. This is a standard MSC/IqASTRAN procedure which is outlined in detail in the Superelement H a n d b o o k [4]. In case of the residual superelement, the user can specify the eigenvalue extraction method that is to be used for the fluid points and the structural points. The modal reduction is carried out only on the component specified on this entry. The following data blocks are available for post processing: (1) Pressure values. This is the solution vector for the fluid points. N o additional user interface is needed. (2) Particle velocity. This vector is recovered in an intermediate step when the fluid displacement vector is computed. N o additional user interface is needed. (3) RMS sound pressure level. The user requests this data by entering RMS in the SPL entry of the ACOUT entry. When this data is requested, the RMS values ( = P e a k / ~ - ) are stored in a data block in addition to the usual pressure values (which are peak values). (4) Sound pressure level in dB. When the user provides the peak reference pressure in the PREFDB field of the ACOUT entry, an output data block containing this information is supplied. The dB level is defined as dB = 20 I o g ( P R M s / P o ) (for air P0 = 20 × 10 - 6 Pa). (5) Acoustic modal participation factor. The user requests this output data block by specifying the number of acoustic modes for which the modal participation factor is computed. (6) Structural modal participation factor. The user requests this output by specifying the number of structural modes for which the modal participation factor is computed. (7) Load modal participation factor. This matrix is printed whenever the structural modal participation factor is requested. The load modal participation factor is defined as the contribution to the fluid solution by the loads which are applied to the fluid. The fluid solution is the summation of the SMP and LPF. (8) Panel participation factor (PPF) as defined earlier is computed and printed for all the panels specified. The RMS calculations are valid in modal and direct frequency response analysis only. The modal participation factors as defined earlier are valid only for modal frequency response analysis. The User request for these calculations will be ignored (with a warning message displayed) for other solutions.
References [1] IZADPANAH, K., "'Modifications/consolidation of the HEXA element", Internal Memo, MacNeal-Schwendler
Corporation, Los Angeles, CA, May 20, 1987. [2] KINSLER,L.E., R.A. FREY,A.B. COPPENSand J.V. SANDERS,Fundamentals of Acoustics, Wiley, New York, 1982.
[3] MORSE,P.M. and K.U. INGARD,Theoretical Acoustics, Princeton University Press, Princeton, NJ, 1968. [4] GOCKEL,M.A. (ed.), MSC/NASTRAN Handbook of Superelement Analysis. [5] SUNG, S.H. and D.J. NEFSKE, "'A coupled structural-acoustic finite element model for vehicle interior noise analysis", Trans. A S M E 106 (4), 1984.
331
Elsevier
Coupled fluid-structure interaction analysis Kamran Izadpanah, Robert L. Harder, Raj Kansakar and Mike Reymond MacNeal-Schwendler Corporation, 815 Colorado Blvd., Los Angeles, CA 90041, USA Revised June 1990 Abstract. A general-purpose fluid-structure interaction capability which can be used in analyzing internal acoustic problems in automobiles and aircraft has been formulated. This capability has become an important issue in the last few years because of the emphasis that has been placed on noise control in all vehicles. A variety of output features that are needed by the analyst in identifying the causes of the acoustic problems have been studied and implemented. The fluid (air) is modeled as an inviscid, irrotational fluid with pressure variables as the degrees of freedom. An acoustic absorber element and an acoustic barrier element are also implemented.
Introduction An internal acoustic analysis capability [1] which is fully automated and requires minimum user input has become an important tool in analyzing internal noise problems in automobiles, trucks, submarines and airplanes. In the automotive industry for example, this type of analysis can help the understanding of vehicle noise which is usually referred to as booming noise or road noise occurring in a frequency range of up to 500 Hz. The two major causes of this noise are the engine and the various structural panels which are excited by the suspension system. With the methodology outlined in this paper and implemented in Version 67 of MSC/ NASTRAN, these problems can be identified and corrected before the actual structure is built and tested. The analyst can perform eigenvalue analysis, frequency response analysis (direct, modal and mixed, i.e., modal fluid, direct structure and vice versa), and, transient response analysis (direct, modal and mixed, i.e., modal fluid, direct structure and vice versa). The fluid-structure coupling terms are automatically calculated for both matching and non-matching interface meshes. In addition, two new elements have been formulated and added to the element library of MSC//NASTRAN which can act as acoustic absorbers and acoustic barriers and are useful when analyzing acoustic noise control devices.
Theoretical background The governing differential equations describing compressible fluid flow are the Euler equation
p[(v.V)v+e] = - r e , and the continuityequation
ap/Ot + PoV " v =
O,
where p is density, v is velocity, P is pressure, and S7 = (O/ax)i + ( O / ~ y ) j + (O/Oz)k. Elsevier Science Publishers B.V.
(1) (2)
K. Izadpanah et al. / Coupled fluid-structure interaction analysis
332
For acoustic analysis, three assumptions are made at this stage: • small motion theory, i.e., v---~; • the convective m o m e n t u m terms (first terms of the Euler equation) are negligible, i.e., p/i = - VP; • the P-p behavior is locally linear, i.e., ( d P / d p ) = YPo/Po where y is the ratio of specific heats and c = ~ is the speed of sound. Assuming an homogeneous fluid, the continuity equation becomes
(dp/dP)(dP/dt)
+ Po V . v = O.
(3)
Using the third assumption and integrating, this equation will yield P = -B V-u,
(4)
where B = c02po is called the "coefficient of compressibility". This relation is now differentiated twice with respect to time and X7 • ii is substituted from the continuity equation to yield the wave equation describing the fluid behavior which is expressed in the pressure variable only:
( l / B ) , ~ = ( 1 / p ) V . VP,
(5)
and the pressure and displacement variables are related by the modified Euler equation: On = - V P .
(6)
The finite element discretization of this equation can now be carried out using standard procedures. Assume: e = E P,N,,
(7)
where Pi and N~ are the pressure and shape functions associated with the ith node of an element. The weak form or the virtual work statement of cqn. (5) can be written as:
f v o l [ ( 1 / B ) P - (1/p)X7 " VP] d V S P = O .
(8)
F r o m the divergence theorem we have:
V " ( aV) = V a . V+ aV . V.
(9)
Applying this theorem to the second term of eqn. (8) yields the following (where a = 8P and V = VP): V • (SP WP) = V S P - V P + 8P V • V P ,
(10)
or
= -v
p. v e + v .
re).
(11)
However, VSP- VP
=
(12)
½5VP- V P ,
and therefore ~l'(V. re)
= - ½~(vP'vP)
+ V " (SJ' r e ) .
(13)
Also Pap
=
(14)
and fvol v • x dV = fsx. ds,
(15)
K. lzadpanah et al. / Coupled fluid-structure interaction analysis
333
where ds is an incremental surface vector (outward). Substituting these relationships in eqn. (8) yields:
LoIS[(1/2B)PZ+(1/2p)VP.VP)] dV- fs(1/p)vP.dsSP=O.
(16)
The first integral in eqn. (16) represents the virtual work statement of the differential equation (5). The second integral in eqn. (16) is the boundary condition associated with this problem. To obtain a valid solution, these two components have to be equal to zero for all admissible 8P. In order for the second integral to be zero for all admissible 8P, either 8P or v P " n has to be zero. Since 8P is an arbitrary admissible variation, it could be zero only when P is specified on the boundary. This is the essential boundary condition for this problem. If P = 0, the essential boundary condition corresponds to the conditions on a free surface. If P is not specified on the boundary then, V P " n or must be zero on the boundary. The is the natural boundary condition associated with this problem and it occurs at rigid walls. From eqn. (16), the last term in this relation can be written as:
OP/an
fs(ii, d s )
8P,
8Plan
(17)
or in matrix form as: [A](~},
(18)
Aij = fsNjNi ds.
(19)
where
Noting that the acceleration of the fluid is in a direction opposite to that of the structure and denoting the acceleration of the structure by (fi},, the problem is reduced to solving the following set of linear equations: [ M ] ( t i} + [ K ] { P }
+[A]{/i}s=0,
(20)
where
Mij = fvo(1/B)NdVj dr,
(21)
Kij = fvolVN/(]//p)vNj dV.
(22)
Fluid-structure interface
At a fluid-structure interface, the boundary condition that the pressure P must satisfy is:
aP/0n = - p / i ° ,
(23)
where n is the direction of the outward normal. This condition is reflected in the third term of
eqn. (20). On the other hand, the fluid affects the structure by applying forces over the structure surface area. These force vectors can be expressed as: (F} =-[A]T{P},
(24)
The A matrix defined in eqn. (19) is therefore the coupling matrix that needs to be computed in order to do coupled fluid-structure analysis. The calculation of this matrix for the general case
334
K. lzadpanah et al. / Coupledfluid-structure interaction analysis
is not a trivial matter. In fact a substantial amount of effort is necessary for general shapes and contact areas.
Coupled fluid-structure system The general equations of motion of the entire structure can now be written as: (25)
Mki + K x = F ,
where
0]
M--
,26,
Mf
and _A T
,27,
with the solution vector x as: x=
p
,
and the external load vector F=
~.
These relationships are used in a frequency response analysis. The eigenvalues (known to be real) of these unsymmetric real matrices may also be computed. Similarly, in a transient response analysis these matrices are used with the standard time integration procedures of MSC/NAST~N without any modifications. The inclusion of structural damping in this analysis does cause the eigenvalues of the system to be complex numbers. When the fluid is highly compressible, it is possible to symmetrize the K and M matrices by a suitable transformation of variables. This process involves inversion of the fluid matrices and thus is recommended only when these matrides are not large or when they have been replaced by their smaller modal representation. There are several approaches to this problem. The simplest approach is by substituting r for P where r=P-
[gfl-'[A]{u
}.
(30)
The resulting matrices for K and M are:
0
and F=
Fs ff"
(33)
K Izadpanah et al. / Coupled fluid-structure interaction analysis
335
Modal decomposition In very large problems, it is sometimes desirable to apply the modal decomposition method in order to reduce the size of the analysis. In this case the fluid, the structure, or both are possible candidates for this operation. Applying dynamic reduction to the structure, where the modes of the structure used for reduction are computed in vacuum results in:
u = .~s, [m.] = . T M . , ,
(34)
[k~] = . V g . s.
(36)
(35)
For instance, eqns. (35) and (36) substituted in eqns. (26), (27) and (29) result in
[m~ 0 ] M=A. ~ Mr'
(37)
K=[~'
(38)
-':AT]Kf '
and F - - [ "TF]
(39)
Lr,/• Applying dynamic reduction to the fluid (the modes of the fluid in a rigid container) results in: p = .f~f,
[mf ] =
(40)
"TM'f,
(41)
[kf ] = .fTK.f,
(42)
which when substituted in eqns. (26), (27) and (29) yields: M=
"~A
mf
and
-=[o AT° kf
(~)
J
with the same force vector as before (eqn. (29)). It is also possible to perform the dynamic reduction in both the fluid and the structure. The resulting stiffness and mass matrices are a combination of the ones derived above with the resulting modal mass and stiffness matrices:
M= ~ A . ,
mf'
K=[~ -'[AT'f]kf"
(46)
The equation of motion becomes
ms "?'~'s
0 mf
~'~ + ~f
--
kf
] t ~f J
. ['fTFf I
(47)
K. lzadpanah et a L / Coupled fluid-structure interaction analysis
336
Note that the modal analysis method is a useful method if the analyst is fully aware of the reduction process that is being carried out. It is possible to conceive of certain problems where the results obtained by an inexperienced analyst will be completely unacceptable. In that case it is recommended that a direct analysis be performed. On the other hand, the direct method of solution involves either sparse but unsymmetric matrices or dense but symmetric matrices. This means that for certain type of problems the computational effort involved might be significantly more than for a correctly applied modal reduction method.
Modal and boundary participation factors The modal solution method provides valuable information regarding the structural and modal participation in the acoustic response [5]. The acoustic modal participation at the interior points ( p~ } can be computed by expanding eqn. (40) as: L
{p} =
= E
(48)
i=I
where L is the total number of acoustic modes used in the analysis. Each vector (p~ } then represents the contribution to the acoustic response from the acoustic mode i or the modal acoustic participation factor. The structural modal participation factor can be computed for a harmonic solution of radian frequency o~ as follows. First the equation of motion is solved for the fluid modal vector: {if } =
[z2](~o2[a]{t2~}),
(49)
where [a] = [*f]T[A][~s]
(50)
is the modal coupling matrix and [z21 = ( - , : [ m f ]
+ kf)-'
(51)
is the acoustic modal frequency response function. Substituting eqn. (49) into eqn. (40) results in: SM
(p}=J[*t][z2][a](~s}=
•
(p:},
(52)
j=l where S M is the total number of structural modes. Note that each vector { p: } represents the contribution of the structural mode j to the the acoustic response or the structural modal participation factor. The boundary panel participationfor a boundary panel b is basically the collectivestructure modal participation of all the structural nodes on that panel. The [Ab] matrix which is the coupling matrix due to panel b therefore has to be computed. The global [A] matrix will be the assembly of all the panels in the structure. Once the [Ab] are known, eqn. (50) can be written as: N
(53)
[a] = [q~flT y , [ A e ] [ ~ , l , b=l
which when substituted in eqn. (52) becomes: N
{p}
=~,2[~fl[z21[*flT E
b=l
N
= Y'. ( p b } , b=l
(54)
K. lzadpanah et al. / Coupled fluid-structure interaction analysis
337
where N is the total number of boundary panels. Note that each { Pb } vector represents the contribution of boundary panel b to the acoustic response.
Simple acoustic source
A "simple" acoustic source is defined as an enclosed surface vibrating with arbitrary velocity distribution but is of such size that all dimensions are much smaller than the wavelength of the emitted sound [23]. It can be shown that the pressure field of any irregular simple source m a y be approximated by the pressure field produced by a pulsating sphere of the same source strength. The source strength Q, is defined as the rate of flow of the fluid away from the source given by: Q e j'°' =
fsh.n
dS,
(55)
where ~i is the velocity vector and n is the outward normal. For a sphere pulsating at a given frequency to, the total power can be calculated as:
1I = ½~rOoc(Q/?~ )2,
(56)
where 11 is the total power, c is the speed of sound in the fluid, and h = 2~rc/to is the wavelength. This expression can be inverted to express the source strength as a function of total power H:
1 / 8~rcH O = -d
(57)
p
This is the expression that is used in MSC//NASTRANto convert the power versus frequency curve supplied for a given source into the right hand vector (loading) for the pressure variables. Note that whenever a phase angle ~b is specified by the user, this load vector should be viewed as the magnitude of a complex loading vector hating an angle ~. This complex vector is decomposed into real and imaginary load vectors.
Acoustic barrier
An acoustic barrier is the component that is usually placed between two compartments such as the barrier between the engine compartment and the passenger compartment. The purpose of the barrier is to reduce the noise transmission between the two compartments. Even though there are m a n y different types of barriers, the most common construction is a sandwich consisting of a steel backing, a foamy middle layer, and a plastic type airtight cover usually called the septum. The transmission phenomenon is a very complicated physical phenomenon. In the adopted approach, this behavior is simplified by assuming that the incident wave is normal to the surface and that there is no in-plane (parallel to the sandwich) transmission of energy. In other words, the barrier model constructed here is a relatively crude model that gives reasonable results in one-dimensional cases and should be used with care in order to avoid erroneous conclusions. There is a standard (ASTM) test to measure the transmission loss (TL) of any barrier defined as T L = 20 log( Pi/Pt
),
(58)
K. hadpanah et al. / Coupled fluid-structure interaction analysis
338
where Pi is the incident pressure and Pt is the transmitted pressure. Note that some of the pressure wave hitting the receiving side of the barrier is reflected and that in the definition of TL only the incident pressure and not the total pressure is used. In a typical TL versus frequency curve, it is seen that the TL increases about 6 dB per octave in the low-frequency range and about 18 dB per octave in the high-frequency range. There is a transition region between these two ranges which is characterized by very low TL. In this region the sandwich construction is basically vibrating in resonance and the barrier is least effective. The "equivalent structural" system selected to represent the transmission phenomena consists of a spring sandwiched by two masses. The mass of the backing and of the septum is usually known. However, the stiffness of the spring is going to be calculated from the resonant frequency of the sandwich as observed in experiment. The theoretical resonant frequency of this model can be calculated as:
1 //¢
=
(59)
/Me'
where
M,M Me -- MI + 3,/2
(60)
and M 1, M 2 are the masses of the backing and the septum respectively. This means that the equivalent stiffness can be calculated by
M~M2(2~Ir.~) 2 K=
3,/1 + M 2
(61)
Once the effective stiffness t¢ is calculated, the stiffness and mass matrices of the barrier element can be calculated in the usual manner, i.e., for face 1 (steel): K U -- ~ f u , uj d a , ,
(62)
M . = M, f
(63)
d& ;
and for face 2 (septum): Kij =
f u,% dA2,
(64)
M,j = M2fm% dA2
(65)
The justification for using this equivalent structural system is that when it is subjected to the same loading and boundary conditions as the TL test, the results are close to experimental results. The pressure distribution in the T L test when the source is in the left chamber (x < 0) can be thought of as:
P =
Pi e - i ~ ( t - x / c ) + Pr e-i~(t+x/c),
x = 0-,
Pt e - i ~ ( t - x / c ) ,
x = 0 +,
(66)
where Pi is the incident pressure, P, is the reflected pressure, Pt is the transmitted pressure, c is the speed of sound, and i = ~-Z-i-. In all of the complex equations, only the real part is to be used.
K. lzadpanah et al. / Coupledfluid-structure interaction analysis
339
This pressure distribution is basically a right traveling incident wave, a left traveling reflected wave, and a right traveling transmitted wave. Note that because of the special rooms used in this test there are no left traveling waves on the back side of the sandwich. The boundary conditions on the left and right side of the sandwich are: pii + OP/Ox = 0,
(67)
which means:
fi( t i t
x0
/i(0-~)
x=0 +
u=
(68)
e -'°',
Now applying the equilibrium equations results in: i
( - ml., 2 + K) ~--/~ (P, - Pr) - K0.,c P, = Pi + Pr, i
--Kpa~c ( P i - Pr)
"
+ (3/2 ¢°2 + K ) ~ c Pt = - P t .
(69) (70)
These relations can be solved to obtain: Pt = - 2iKpwc Pi [ M1M2to4 - ( ( M 1 + M 2 ) K + (pc)2) oa2] + i [ ( M 1 + M2)Pcto 3 - 2Kptoc] ' (71) and, also: Pi+Pr Pt
ptoc(-Mlto2+K)+i[M1M2oj4-(MI+M2)Kto2] p~cK
(72)
With the proper choice of the resonant frequency, there is good agreement with test results.
Acoustic
absorber
An acoustic absorber is a component that is usually placed inside a compartment in order to minimize the reflection of sound. In an automobile, for example, the interior is padded with different types of absorbers depending on the available space and the acoustic objective in mind. An acoustic absorber is characterized by its specific acoustic impedance z defined as the ratio of the acoustic pressure in a medium to the associated particle speed, i.e., z = P / u . This is a complex number with the real part usually called the specific acoustic resistance and the imaginary part the specific acoustic reactance. These parameters are measured in a standing wave tube (ASTM standards) and have some general characteristics. As in the case of an acoustic barrier, the behavior of an acoustic absorber is modeled by structural components. Here a single mass in conjunction with a spring and damper in parallel is the "equivalent structural" model that is used. In the standing wave test, the absorbing material is secured to the end of the tube and a pressure wave of P sin a~t is applied. Therefore in the equivalent structural model, the spring and damper components are grounded and the equilibrium equation for this one-degree-offreedom system becomes: M g + CYc + K x = P sin ~0t.
(73)
340
K. Izadpanah et aL / Coupledfluid-structure interaction analysis
The acoustic impedance of this system is u
:~
C+
(
=t
~ 0 M - -~ 1
(74)
The test data are interpreted as follows: the resistance ~ ( z ) is averaged and used as the damper C. The reactance data .St(z) is used a least square fit to the curve
Ma~ - K/~o.
(75)
The resulting error function to be minimized for this least square fit is
(
~i Wi M°°i -- K __ Bi
;
(76)
where B i are the experimental data and W~ is a weighting function. The simultaneous equations that have to be solved for this problem are:
- E Wi
E "72Wii ] [ K J
- E "7'W, Bi
(77)
The summation can be carried out over an entire range of data or for the data corresponding to frequencies below a cutoff frequency. Solving for M and K results in: E o~Z2WiEc.oiWiB, - E WiEcaoZ1Wi "`
M
E-72 -
K
=
(78)
- (E
E +:w,E,+-'w,B, E E -,-- _ (E w,) +
2
(79)
The resulting M and K values can now be used in the structural model of the acoustic absorber.
Implementation
This methodology was implemented in the general-purpose finite element code MSC//NASTRAN Version 67 with the goal of having a very clean interface with existing capabilities. The fluid medium is described by the use of the standard GRID and element entries. The fluid grid points are distinguished from the structural points by a parameter on the GRID entry. The fluid elements that are available are in the shape of hexahedra, pentahedra or tetrahedra. Both linear (8-noded HEXA, etc.) and quadratic (20-noded HEXA, etc.) elements are supported. The topology (connectivity) of these elements is specified using the existing Bulk Data entries for these shapes, i.e., CHEXA, CPENTA, CTETRA. The property entry for the solid elements is the key in distinguishing between structural and pressure elements. By designating the fluid elements on this entry, the pressure-based fluid elements will be used. The bulk modulus and the mass density of the fluid have to be specified for these elements. It is sometimes necessary to find the contribution of a set of grid points to the noise level in an acoustic cavity. For example, in the automotive industry the entire roof assembly contributes significantly to the noise level inside the passenger cabin. Consequently it is necessary to compute and minimize the contribution of this "panel".
K. Iza@anah et aL / Coupled fluid-structure interaction analysis
341
The user is able to define panels by specifying the structural grid points in each panel using a new entry PANEL and convenient set definition tools. The acoustic absorber and barrier elements were implemented as new elements CHACAB and CrtACBR. These elements consist of 8 to 16 structural grid points, and define two surfaces. The property entry for the absorber elements points to three tables which contain the resistance (Y,), reactance (I12) and the confidence level (Y3) of the material as a function of frequency as measured in the standing wave tube test (ASTM standard). This data will be used to calculate the value of the mass, spring, and damper used in the structural model that will emulate the behavior of the absorber material. The property entry for the barrier contains the mass per unit area of the backing, mass per unit area of the septum, and the resonant frequency of the barrier. The two mass values are readily available and the resonant frequency of the sandwich is defined as the frequency with the lowest values of transmission loss (ASTM test). This frequency will be used to calculate the value of the spring used in the structural model that will emulate the behavior of the absorber material. The loading of the fluid is in most cases analogous to enforced displacements in structural mechanics analysis. There are several types of loading that have to be provided for the fluid elements: (1) Constant enforced pressure at the grid points; (2) Frequency or time-dependent enforced pressure at the grid points; (3) Acoustic source characterized by a volumetric flow rate corresponding to a power spectral density function. The large mass approach that is currently used to model enforce displacements is used to have enforced pressure loading. A simple source is modeled as a pulsating sphere in infinite space. The user can conveniently model this type of source by supplying the power versus frequency curve characterizing the acoustic source and the phase angle. The phase angle is useful whenever multiple sources are present. The single-point constraints of the fluid ( P = 0.0) can be enforced using the existing constraint entries. This type of boundary condition occurs at free surfaces. The identification of the fluid-structure interface region is accomplished by using a series of sorting operations. Ideally, the computations involving the fluid-structure interface should be completely invisible to the user. This is, however, a computation intensive process. In the current design the user is requested to "assist" the process in order to expedite the calculations. This user "assistance" is in the identification of the "wetted area" of the structure. Its only effect will be the reduction of the sorting operations carried out for the computations involved. To expedite the solution process, the user needs to supply: (1) Type of interface, i.e., matching or non-matching fluid-structure meshes; the default interface in the program is the non-matching interface; (2) Actual structural and fluid grid points that are in contact. The program always outputs a table containing the interface information. The capability is fully operational with the substructuring (superelement) capability of MSC//NASTRAN. In these superelement analyses only the following restrictions apply (SOL 107, 108, 109, 110, 111, and 112 of MSC/NASTRAN): • Each upstream superelement can contain some of the fluid or some of the structure but not both. • The fluid-structure interface grid points can only be in the residual superelement. This means that upstream superelements of different types will not contain the coupling terms (the A matrix). It is possible to perform modal reduction on the upstream fluid or structure superelements or
342
K. lzadpanah et al. / Coupled fluid- structure interaction analysis
both using the standard superelement case control capabilities in conjunction with multiple subcases. This is a standard MSC/IqASTRAN procedure which is outlined in detail in the Superelement H a n d b o o k [4]. In case of the residual superelement, the user can specify the eigenvalue extraction method that is to be used for the fluid points and the structural points. The modal reduction is carried out only on the component specified on this entry. The following data blocks are available for post processing: (1) Pressure values. This is the solution vector for the fluid points. N o additional user interface is needed. (2) Particle velocity. This vector is recovered in an intermediate step when the fluid displacement vector is computed. N o additional user interface is needed. (3) RMS sound pressure level. The user requests this data by entering RMS in the SPL entry of the ACOUT entry. When this data is requested, the RMS values ( = P e a k / ~ - ) are stored in a data block in addition to the usual pressure values (which are peak values). (4) Sound pressure level in dB. When the user provides the peak reference pressure in the PREFDB field of the ACOUT entry, an output data block containing this information is supplied. The dB level is defined as dB = 20 I o g ( P R M s / P o ) (for air P0 = 20 × 10 - 6 Pa). (5) Acoustic modal participation factor. The user requests this output data block by specifying the number of acoustic modes for which the modal participation factor is computed. (6) Structural modal participation factor. The user requests this output by specifying the number of structural modes for which the modal participation factor is computed. (7) Load modal participation factor. This matrix is printed whenever the structural modal participation factor is requested. The load modal participation factor is defined as the contribution to the fluid solution by the loads which are applied to the fluid. The fluid solution is the summation of the SMP and LPF. (8) Panel participation factor (PPF) as defined earlier is computed and printed for all the panels specified. The RMS calculations are valid in modal and direct frequency response analysis only. The modal participation factors as defined earlier are valid only for modal frequency response analysis. The User request for these calculations will be ignored (with a warning message displayed) for other solutions.
References [1] IZADPANAH, K., "'Modifications/consolidation of the HEXA element", Internal Memo, MacNeal-Schwendler
Corporation, Los Angeles, CA, May 20, 1987. [2] KINSLER,L.E., R.A. FREY,A.B. COPPENSand J.V. SANDERS,Fundamentals of Acoustics, Wiley, New York, 1982.
[3] MORSE,P.M. and K.U. INGARD,Theoretical Acoustics, Princeton University Press, Princeton, NJ, 1968. [4] GOCKEL,M.A. (ed.), MSC/NASTRAN Handbook of Superelement Analysis. [5] SUNG, S.H. and D.J. NEFSKE, "'A coupled structural-acoustic finite element model for vehicle interior noise analysis", Trans. A S M E 106 (4), 1984.