Finite element formulations for acoustical radiation

Journal of Sound and Vibration (1983) 88(l), 47-64

FINITE

ELEMENT

FORMULATIONS

ACOUSTICAL

FOR

RADIATION

R. J. ASTLEY AND W. EVERSMAN University of Missouri-Rolla, Rolla, Missouri 65401, U.S.A. (Received 11 March 1982, and in revised form 20 July 1982)

Finite and infinite element techniques are applied to linear acoustical problems involving infinite anechoic boundaries. Theory is presented for a simple one dimensional model based on Webster’s horn equation. Results are then presented both for the one dimensional model and for two axisymmetric test cases. Comparisons with exact solutions indicate that both the infinite element and wave envelope schemes are effective in correctly predicting the near field. The wave envelope scheme is also shown to be capable of resolving the far field radiation pattern.

1. INTRODUCTION The treatment of radiation boundary conditions at the “open” boundaries of wave propagation problems presents significant difficulties when numerical solutions are sought. “Sommerfeld” type boundary condition, which The imposition of a straightforward implies an assumption of locally plane, outward travelling waves, is appropriate only at large distances from the generating or refracting mechanism. This constraint frequently requires an unacceptably large computational domain and is particularly troublesome if the length scale of disturbances is already much smaller than the overall geometrical length scale of the problem. An important aeroacoustic problem which exhibits these characteristics is that of the radiation of fan noise from the nacelle inlet of a turbofan aircraft engine where the directivity and intensity of the acoustical far field are known to be sensitive to the geometry and lining configuration of the nacelle, the dimensions of which are generally many times larger than the acoustical wavelengths of the major energy carrying frequencies [l]. In modelling this problem for realistic sound levels the acoustical field may be regarded as linear and for a given frequency the acoustical pressure may then be expressed as the solution of a modified Helmholtz equation, The problem so formed is similar in general character to many other problems involving linear wave phenomena and although the investigation reported in this paper is directed primarily towards the development of numerical models for turbofan inlets, the methods described are clearly applicable to a broader range of physical problems involving wave radiation at infinite boundaries. The imposition of radiation boundary conditions other than by the application of an approximate Sommerfeld condition to a finite numerical solution [2, 31 may be achieved by a variety of alternative methods. Closed form analytic solutions are available for a limited number of simple geometries. Such solutions range from the variables separable and source distribution methods of Rayleigh and others [4, 51 to more recent studies in which the Wiener-Hopf technique [6-81 has been used. None of these solutions offers any extension to arbitrary nacelle geometries. 47

0022-460X/83/090047 + 18 $03.00/O

@ 1983 Academic

Press Inc. (London)

Limited

48

R. .I. ASTLEY

AND

W.

EVERSMAN

An alternative approach is to be found in a number of formulations in which the problem is solved in purely numerical form only within a relatively small inner domain. This encloses any awkwardly shaped boundaries and non-homogeneities (e.g., density or mean flow variations in acoustical problems, depth changes in analogous water wave problems). At the boundary of this region the numerical solution is matched to an analytical far field solution with the correct radiation properties. The matching may be achieved either through an appropriate eigenfunction expansion at the matching surface [9] or via source distributions or related boundary integral equation techniques [ 10-121. Although this type of approach has proved successful for a number of problems the requirement of an analytic solution in the outer domain is an unwelcome constraint on the initial field equations. Such methods also tend to destroy the sparsity of the stiffness matrix when a finite element (FE) discretization is used for the numerical portion of the solution. An approach which does not suffer from these disadvantages and which has proved effective in modelling the refraction of water waves about arbitrarily shaped objects [13] is the application of “infinite” elements. This modification of “finite” element methodology involves the extension of the computational domain to infinity and the inclusion of an exponentially decaying, outward travelling, wavelike variation in the shape functions for the dependent variables. The assumption of exponential decay ensures that the integrals involved in the calculation of the stiffness terms are, indeed, finite over the infinite domain. The same assumption has the disadvantage, however, of violating the known asymptotic behavior of the radiated component of the solution at large distances from the refracting or generating mechanism. The infinite element approach may therefore be expected to capture the solution adequately within the inner numerical domain, but will give little reliable information regarding far field radiation. In this paper a new approach is proposed whereby a conventional FE solution in the inner region is compatibly matched to a “wave envelope” FE solution in a large but finite outer region. The term “wave envelope” is used since this approach closely resembles the FE and finite difference wave envelope schemes used successfully for internal duct acoustics problems [14,15]. The correct asymptotic behavior of the solution at large distances is preserved by incorporating a wavelike variation with the correct asymptotic decay in the shape functions for an outer region. The Sommerfeld radiation condition is then applied in the usual way at a distant, but finite boundary. The formulation of the problem in the outer domain and the evaluation of the appropriate stiffness terms are simplified if the complex conjugates of the shape functions rather than the shape functions themselves are used as weighting functions in a modified Galerkin process. The method is initially presented for a simple one dimensional model based on the solution of Webster’s horn equation. FE solutions are detailed for the wave envelope, conventional and infinite element schemes and computed results compared with exact solutions. Details of the extension of the method to two and three dimensions are omitted. but results are presented for two axisymmetric test cases. 2. ONE 2.1.

GEOMETRY

AND

GOVERNING

DIMENSIONAL

MODEL

EQUATIONS

A one dimensional model which incorporates the necessary features for the presentation and testing of the alternative FE formulations is provided by the solution of Webster’s horn equation for an asymptotically conical duct. The equation to be solved is an integrated low frequency approximation for acoustical transmission in a duct of variable section A(x) with essentially plane wave propagation at all points. If the acoustical

RADIATION

FINITE

ELEMENT

49

FORMULATIONS

pressure is assumed to vary with time as eiot the complex pressure amplitude p(x) is a solution of (1) where the reduced frequency k is given by k = w/c (c denotes the local sound speed). If it is assumed that A(x) -x2 as x becomes large, equation (1) becomes a spherically symmetric form of Helmholtz equation with x representing a radial co-ordinate. Although the numerical formulations are presented for a general function A(x) the results which follow are calculated for the specific case of a uniform cylindrical section joined to a conical expansion. For this particular configuration, shown in Figure l(a), A(x) is given by

A(x)=

i

“7

p,ston

I)- ----

A07 O~X~X()

Aoblxo)*,

x

I

(2)

>xn ’

Rigid bcWda'Y /IA(~)

-__~,~___~-_l+__~___

_(;~xj------,

((]I x

II

Nodej

-------?

(b)

L Figure 1. One dimensional discretization (infinite element

model. (a) Geometry; (b) FE discretization scheme); (d) FE discretization (wave envelope

(conventional scheme).

scheme);

(c) FE

p(x) = (A’ eplkx +B’ eikx)/kx, x >x,,,

(3)

and exact solutions are readily found of the form p(x) =A epikx +B eikx, Ocx 4x0,

where the constants A, B, A’, and B’ are determined by boundary conditions at x = 0 and x = CO. Appropriate matching of pressure and its axial derivative at x =x0, the junction of the cylindrical and conical sections, is also required. At the left-hand end of the duct (see Figure l(a)) a vibrating piston of velocity amplitude a0 is assumed to be present, giving a boundary condition

dp=

-ika

dx

0

atx =O.

50

R. J. ASTLEY

AND

W.

EVERSMAN

The boundary condition at the open end must simulate the existence at large distances from the piston of outwardly travelling waves only. For the analytic solution this is achieved by setting B’ = 0. An approximate Sommerfeld type boundary condition for this termination is given by (5) dp/dx+ikp=O atx =x1, where x1 is a distant boundary. For specific values of k and x0 the imposition of either of the above boundary conditions (at x = co or x = xi) fully determines the values of the constants A, B, A’, and B’. The resulting analytic solution comprises a standing wave pattern in the region x x0. If the approximate Sommerfeld condition is used (instead of the exact anechoic condition) a small negatively propagating component is generated in the region x > x0. A secondary standing wave pattern is then superimposed on the outwardly moving component. If x >>xothe effect is small and vanishes as x L-+ 00. 2.2. FINITE ELEMENT DISCRETIZATION AND SOLUTION Solutions are presented in this section for conventional, infinite element and wave envelope FE formulations of the one dimensional model. All three approaches may be regarded as variants of a general Galerkin scheme and commence with the initial assumption of a trial expansion p”(x) of the form F(X)= i

UJVi(X).

16)

i-1

The unknown coefficients ai constitute the degrees of freedom for the idealization. The basis functions N,(x) may then be chosen in a variety of ways. For an FE discretization the basis functions are identified with appropriate global shape functions and the coefficients ai then correspond to nodal values of p. Equation (6) is assumed valid throughout the solution region 0
A Galerkin orthogonalization W,(i=l,...,n)thenyields

(7)

of R with respect to a complete set of weighting functions

I,” W,[&$(A(x)g)+k’jT]Adx=O,

i=l,...,n,

and integration by parts gives I,“A(~)[k~W~~-~~]dx+~‘[WiA~x~~]=0, *o

i = 1,. . . , II.

The incorporation of boundary equations (4) and (5) is then effected by the replacement of dp’/dx in the second term of the above equations by -ikao and -ikp’ at x =x0 and x = xi, respectively. This yields a final form of the weighted equations which includes both field and boundary residuals,

i = 1,2, . . , n.

(8)

RADIATION

FINITE

ELEMENT

51

FORMULATIONS

The substitution of the trial expansion for p”then completes the solution and yields a set of n linear equations, which may be written in matrix form as

Kl @I= PI, k2WiNj-Tz]

(9)

dx-ikA(x)WiNjlx=,,,

F; = -%a0 WiA(x)jIzo,

i,i=1,2

,...,

i = 1, 2, . * . ) n,

n,

(10) (11)

{a} being a vector whose components are the coefficients Ui of the original trial expansion. Appropriate choices of Wi and Ni now yield the various formulations.

The conventional formulation is obtained by selecting the basis functions Ni and weight functions Wi to be standard global FE shape functions. For the solutions presented in this paper two noded linear elements are used and the resulting global shape function Njvi,associated with a typical node i, has the form shown in Figure l(b): i.e., Ni and Wj are given for the conventional formulation by Nj = Wj = Nj = Wj = (x

-Xj+l)/(Xj

-Xj+l),

(X -Xi-l)/(Xj Xi SX

-Xj-l)r

S Xj+l,

xi-1

S x 6 Xi,

Nj= Wj=O,~>xj+lOrx
(12)

where xi-1, xj and xi+1 are the co-ordinates of adjacent nodes. An expansion of the trial solution in piecewise linear functions of this type will clearly require many terms to represent adequately a simple wave-like spatial variation. A relatively fine grid is therefore required throughout the solution region. The terms in a typical element stiffness matrix [k]’ for this formulation are given by 1; kLp =

I x;

A(x) (N:N;k2

- NLXN& ) dx,

ff, P = 132,

(13)

where N’,, (Y= 1,2, are conventional element shape functions for a linear element with nodes 1 and 2 at xi and x5. A further contribution to [k]” occurs if the element extends to the outer boundary (i.e., if x; =x1) and is given by kZ2 = -ikA(x,). The components of [k]’ are assembled in the usual way to form the system matrix [K] of equation (9). The infinite element formulation proceeds along similar lines, but requires modification to the definition of Ni within an outer region, x >xj say. For x
52

R. J. ASTLEY

AND

W.

EVERSMAN

The weighting functions are again chosen to be the shape functions themselves and the Galerkin process proceeds as before. When the expressions for Nj and Nj+r of equations (14) are substituted into equation (1Oj an infinite element stiffness matrix is obtained, of the form cc + y(N&N; + N;xN:) - NzxN&] dx, k:, = A(x) e -v’2x-x’--xB ‘[(k2 - $)N:N; I x; fx, p = 1, 2, (,15) where y = ik + l/L

and the functions, Nz, (Y= 1,2, are element shape functions for a

conventional linear element with nodes at XT and x; given by xj and xj+r (all contributions

from the boundary terms at x = x1 vanish as x1 + CO).Some practical difficulties arise in the calculation of these stiffness contributions since they include integrals of the form co e -vxx “A (x) dx, m = 0, 1, 2. (16) I 1; Gauss-Legendre numerical integration is no longer appropriate and if numerical integration is required a more advanced scheme must be used. A modified Newton-Coates formula has been proposed [13] requiring five integration points for the current problem, with weights which must be calculated as functions of L and k. For the specific test case considered in this paper the existence of an analytic expression for A(x) permits :he evaluation of such integrals without recourse to numerical integration. The wave envelope formulation requires two distinct modifications to the preceding schemes. The shape functions must again be redefined in an outer, but finite, region and a different choice of weighting functions must be made. The inner region will again be denoted by x
= (Xk/X)

ePik(xP*S)N~(x),

(17)

where N; is the conventional global shape function for that node (given explicitly by the expression for Nj of equation (12) with the subscript i replaced by k). Once again the property Ni(xi) = Sij is preserved and the parameters ai retain their significance as nodal values of p. The trial expansion so formed can moreover represent exactly the known asymptotic form of the pressure field for large values of x (i.e., (l/x) eeikx). The global shape function Nj for node i at the interface between the conventional and wave envelope regions is shown in Figure 1 (d). The formulation is completed by the selection of the complex conjugates of the shape functions to form the weighting functions Wi. This choice has the attractive property of removing harmonic spatial variations from the integrals involved in the calculation of the stiffness matrix. The resulting expressions for the components of [k]” are then given by

RADIATION

FINITE

ELEMENT

FORMULATIONS

53

where NZ, LY= 1,2, is the shape function for a conventional linear element with nodes at XI and x;. Simple Gauss quadrature may be used to evaluate these stiffness contributions irrespective of the size of each element and the number of spatial wavelengths of the solution over which the element extends. The contributions to [k]’ from boundary terms at x =x1 no longer vanish as x1 becomes large and a finite but distant outer boundary must therefore be retained. An additional element contribution, k& = -ikA(xl),

then occurs if x; corresponds 2.3.

NUMERICAL

to the outer boundary at x =_Y~.

RESULTS

One dimensional results were calculated for the duct shown in Figure l(a) for a range of values of duct length and non-dimensional frequency. The general features of these solutions were similar in all cases and for ease of presentation the results which follow were calculated for a “typical” frequency defined by kxo = 27r, giving a single axial wavelength variation in the solution between the piston and the start of the conical section. For the conventional and wave envelope formulations the finite outer boundary is defined by x = xi, and is specified in the majority of cases by xi = 5x0. This choice is sufficiently large for the Sommerfeld radiation condition to introduce only minor spurious reflections. Results are presented as axial plots of the real part and absolute value of the computed nodal pressures. The real part gives an indication of the extent to which the various schemes are resolving the locally wave-like behavior of the pressure field whereas the absolute value is a measure of the large scale decay of the solution and is usually the quantity of most practical interest, relating as it does to measured sound pressure levels. Exact solutions calculated from equations (3) are presented in all plots for comparison with the computed results. Results for a conventional formulation with 10 elements per wavelength are shown in Figures 2(a) and (b). The outer boundary is placed at a distance 5x0 from the piston. Except for a small phase lag good correspondence exists between the exact and computed solutions. Slight undulations in the absolute value of the computed solution for x )x0 (Figure 2(b)) represent the effect of reflection due to the approximate radiation condition imposed at x =x1. Clearly for x1 = 5x0 this effect is small. Conventional results for a coarser mesh with only five elements per wavelength are shown in Figures 3(a) and (b). The phase error in the computed solution is noticeably larger and the computed absolute value of pressure shows significant deterioration. In order to compare like with like in the assessment which follows of the various schemes the resolution of Figures 2(a) and (b) (i.e., 10 elements per wavelength) will be considered as an acceptable “standard” case for regions where conventional elements are used. Before considering the solutions obtained from the infinite element and wave envelope schemes it is interesting to note the effect of reducing the dimensionality of the conventional solution by applying the approximate Sommerfeld condition at a boundary closer to the generating piston. A plot of absolute pressure for a conventional formulation with x1 = 2x0 is shown in Figure 4. The resolution used is the same as in Figures 2(a) and (b) (i.e., 10 elements per wavelength) but the total number of degrees of freedom is reduced from 51 to 21. It can be seen that the cost of a reduction of the dimensionality of the problem in this way is the introduction of spurious standing waves superimposed on the outward travelling wave (indicated by the “wiggles” in the computed solution compared with the exact anechoic solution). Analytic solutions are shown, in this case only, both for an anechoic termination and for an approximate Sommerfeld type condition at .Y= x’1

54

R. J. ASTLEY

AND

W. EVERSMAN ,

I

,

I

/ i

Anal

co-ordtnale,.x/x0

Figure 2. Comparison of exact and computed axial pressures (a) real part and (b) absolute value; conventional scheme, xt = 5x0, fine mesh, 10 elements/wavelength. - - -, Exact solution; -0-, F.E. solution.

1.5

(bl

Awol

co-ordlnote,

x/x,,

Figure 3. Comparison of exact and computed axial pressures (a) real part and (b) absolute value; conventional scheme, x1 = 5x0, coarse mesh, 5 elements/wavelength. - - -, Exact solution; -0-, F.E. solution.

RADIATION

FINITE

ELEMENT

55

FORMULATIONS

J Figure 4. Comparison of exact and computed axial pressures (absolute value); conventional scheme, x1 = 2x0. fine mesh. - - -, Exact solution, a-, F.E. solution (node values), --, exact solution (finite length, RO-C termination).

(termed an RO-C termination in Figure 4). As one would expect the numerical solution converges well to the latter, but diverges significantly from the true solution. Infinite element solutions for the problem are presented in Figures 5(a)-(d). In applying the infinite element scheme a choice of both the decay length scale L and the position of the outermost node (node i + 1 in Figure l(c)) must be made. The arbitrary choice of these parameters may justifiably be regarded as a disadvantage of the method, although reported results for water wave problems have indicated that the solutions produced tend to be insensitive to major adjustments in either [13]. As the asymptotic attenuation of the solution is known to be of the form l/x as x + co it is clearly desirable to choose L so that the decay term ePX’L approximates this behavior in some sense. Expanding l/x and ePXjL as Taylor series in the vicinity of x =f, say, gives, f~,!~=_~/x=l/(X+Ax~=(l/X)[l-Ax/n+~~~]~f(l)[l-Ax/x+~~~], and g(x) = e - e-(n+dx”L =e -““[l -Ax/L +. 9a]= g(_f)[ 1 -Ax/L +. . .I. The correct value of L in the vicinity of i would therefore appear to be 2 itself. For the infinite element formulation the region in which the exponential term is intended to model the actual decay extends from xi to infinity. A single value of L will clearly not serve to model accurately the whole region, but a choice of L = xi or L = xi+1should give some representation of the solution in the region near the interface with the conventional elements and hence, one might hope, a valid solution within the inner region. The choice of xi+1 itself is arbitrary subject only to the constraint xi+1 > xP In most of the results that follow the value of xj+l is chosen to correspond to the values of x1 used in the other two formulations. A typical infinite element solution with L = xi = 2x0 and Xj+i = 5x0 is shown in Figures 5(a) and (b). As with the conventional FE formulation real and absolute values of pressure are plotted separately and compared with exact solutions. It can be seen that the FE solution matches the analytic solution well within the conventional inner domain, x < 2x0 and in the vicinity of its-boundary with the infinite element. As the axial distance increases, however, the exponential decay built into the FE formulation exaggerates the actual amplitude decay of the solution and results in an underestimate of the pressure amplitude at x = xi”1 (because of the sparsity of nodal values in the infinite and wave envelope regions internodal values of pressure are plotted at regular intervals in the outer region of Figure 5(b) and in all figures which follow). The numerical damping of the solution exhibited by the results of Figures 5(a) and (b) is only slightly modified by alterations in the magnitude of the decay length and the location of the outermost node. This is illustrated in Figures 5(c) and (d) showing results for two different combinations of L

56

R. J. ASTLEY

AND

W. EVERSMAN

Figure 5. Comparison of exact and computed axial pressures; infinite element scheme, xi = 2x0, x,+ I = 5x0. L = 2~; (a) real part and (b) absolute value; (c) as (b) but with L = 5x0; (d) as (b) but with relocated outermost node, xltl = 3~. - - -, Exact solution; 0, F.E. solution (node values); A, F.E. solution (internodal values).

and xj+l. In Figure 5(c) the value of L is increased to 5x0 (instead of 2x0 as in Figure 5(b)). In Figure 5(d) the decay length retains its original value, but the outermost node is moved closer to the finite/infinite element boundary (Xj+l = 3x0, instead of 5x0 as for Figures 5(b) and (c)). It can be seen that these fairly severe modifications to L and to

RADIATION

FINITE

ELEMENT

57

FORMULATIONS

the location of the outermost node have little perceptible effect on the solution in the vicinity of the conventional mesh. They do, however, modify the attenuation of the solution in the outer region. It is worth observing at this point that the infinite element solutions extend to infinity and have only been plotted in the preceding figures for x < 5x0 or x <.x~+~, whichever is greater. If the plots are extended to larger values of x the disparity between the exact and computed solutions increases rapidly whatever values are selected for L and xj+l. This is a consequence of the assumption for exponential decay in the infinite element shape functions, which must inevitably overestimate the attenuation of the solution if x is sufficiently large. Wave envelope results for the same problem are presented in Figures 6(a)-(c). Ten elements per wavelength are used in the conventional inner region and results are presented for an outer boundary at 5x0 with an inner region bounded by x = 2x0. The region 2x0
0

, 0

1

/

I

2

3

4

Am

co-ordinate,

5

x/x0

Figure 6. Comparison of exact and computed axial pressures; wave envelope scheme, one outer element, x, = 2x0, xi+1 = 5x0; (a) real part; (b) absolute value; (c) as (b) but with xi+, = 3.5~~ and x,+2 = 5x0 giving two outer elements.

58

R. J. ASTLEY

AND

W. EVERSMAN

for this formulation is then the same as that for the infinite element scheme of Figures 5(a)-(d) and one more than for the conventional scheme of Figure 5(a). A comparison of the three sets of results serves as a direct measure of the accuracy attainable by each method for the same computational effort. In the near field (x < 2x0) the wave envelope and infinite element schemes both give good results without producing the artificial reflections generated by the conventional scheme of equal dimensionality. The correspondence between exact and computed solutions extends to the outer region, however, only when the wave envelope scheme is used. The accuracy of the wave envelope scheme may in fact be further increased at little extra cost in matrix dimensionality if more than one element is placed in the outer region. An axial plot of absolute value of pressure for a mesh with two elements in the outer region is shown in Figure 6(c). A close comparison with Figure 6(b) reveals a decrease in the (small) error between the exact and computed solutions. At the outer boundary, for example, the value of the error decreases from 5% to 1%. A similar, but more pronounced effect will be noted in the axisymmetric results which are presented in the following section.

3. AXISYMMETRIC

TESTS CASES: THEORY

AND RESULTS

The extension of the present theory to two and three dimensional problems presents no conceptual difficulties. The Sommerfeld condition, where appropriate, is applied at a distant cylindrical or spherical boundary and the wave-like and decay components in the shape functions for the infinite element and wave envelope schemes are appropriately defined in two or three dimensions. A two dimensional infinite element scheme for water waves is described in some detail in reference [13] and the analogous wave envelope scheme may easily be deduced. Specific details of the extension of the current analysis to two and three dimensions are not included in this paper. Results will, however, be presented for two test cases. The first of these comprises the calculation of acoustical pressures generated by a vibrating circular piston located in the centre of an infinite rigid baffle. The exact solution of this problem may be expressed as a uniform monopole distribution over the area of the piston [2] and numerically evaluated with little difficulty. FE solutions for the problem are presented in Figures 7(a)-(d), for which axisymmetric equivalents of the conventional, infinite element and wave envelope schemes of the preceding sections have been used. The elements used are annular, being axially symmetric about an axis perpendicular to the plane of the baffle and passing through the centre of the piston. A plane of revolution is divided into distorted rectangles each of which has sides parallel to lines of constant r and 8, where r and 8 denote radial and latitudinal co-ordinates in a spherical polar system. The geometry of the problem and a discretization of the region adjacent to the piston into a regular mesh of 576 conventional elements is shown in Figure 7(a). Within each element linear shape functions in radial and latitudinal directions are assumed. The elements so formed become direct axisymmetric analogues of those used in the one dimensional model and the discretization along the axis of symmetry is in fact identical in both cases. This type of subdivision is not of course suitable for more general geometries, but serves to demonstrate the method for this specific problem and may be generated relatively easily from the program used to calculate the one dimensional results of the preceding section. The parameters r, and rl in Figures 7(a)-(d) denote the radii of the piston and of the bounding outer spherical surface (where applicable). For the infinite element scheme rl denotes the location of the outermost nodes. The parameter ro denotes the boundary between the conventional

(absolute

Inflnlte baffle

boundary)

variation at centrelme)

p(x) (pressure

ot outer

7 =rp

T

i

P(x)

p(s)

PlShll

vorlatlon boundaryi

voriatm at centrehe)

pressure

Piston

InfInIte baffle

(pressure varmtlon at outer bwndary)

varmt1on 01 centrehe)

p(x) pressure

at o&r

p(s) (pressure

r = rp

I = r.

7. Axisymmetric test case; comparison of exact and computed pressures; (a) conventional scheme, r 1= 6r,, rP = A; (b) infinite element scheme, rl = 2.5r,, rP = A ; scheme, one outer element, ro = 2<5r,, rl = 6r,, r, = A ; (d) wave envelope scheme, two outer elements, rO= 2.5r,, rl= 6r,, rP = A. -, Exact solution value of pressure); 0, F.E. solution (absolute value of pressure); - --, exact solution ireal value of pressure); A, F.E. solution (real value of pressure).

(c) wave envelope

Figure

(pressure varmtlon at centrelml

Piston

boundary!

p(s) (pressure wmatm

p(x)

01 outer

p(s) (pressure variation

60

R. J. ASTLEY

AND

W.

EVERSMAN

inner region and the outer infinite element or wave envelope regions; r. and rl are therefore analogous to x0 and xl in the one dimensional schemes. The acoustical field was calculated for a number of values of the natural nondimensional frequency parameter of the problem, k = wr,,/co (w is the actual frequency, co the local sound speed). The general characteristics of the resulting solutions are similar in all cases and for ease of presentation a single value, k A 27r, will be considered for the results which follow. The resolution used in the conventional regions of the various schemes is similar to that applied to the one dimensional problem, although coarsened slightly to give a clearer representation of the meshes themselves in the output plots. Radially the mesh comprises eight elements per wavelength, (“wavelength” being defined as the wavelength of a plane wave at a large distance from the source), compared with the “standard” 10 elements per wavelength of the previous section. The choice of element spacings in the latitudinal direction is then dictated by an attempt to provide at least equal resolution in the latitudinal and radial directions in the region r
RADIATION

FINITE

ELEMENT

FORMULATIONS

61

interface of the conventional and wave envelope regions. The pressure pattern at the outer boundary, however, is in excellent agreement with the exact solution. The accuracy of the inner solution may be greatly improved at little computational cost by the inclusion of one or more additional layers of wave envelope elements, as illustrated by the results of Figure 7(d), which show a wave envelope solution for the same problem with two layers of elements in the outer region. The far field pattern is again well represented, but in this case the near field also shows good correspondence with the exact solution, the reflections at r = r. having been almost entirely eliminated. The wave envelope results may be further improved by the addition of more elements if desired. This behavior is clearly an axisymmetric manifestation of similar characteristics exhibited by the one dimensional solutions of Figures 6(b) and (c). The ability to vary the mesh resolution in the outer region in this way is a useful feature of the method. It is worth observing that an equivalent refinement in the infinite element solutions would involve changing the number of nodes in the outer elements and hence modifying the shape functions themselves.

Figure 8. Hyperboloid duct geometry.

The second test case, shown in Figure 8, is one which more closely resembles the situation encountered in real turbofan inlets. The geometry is axisymmetric and consists of a cylindrical duct of radius a which is smoothly joined to a hyperboloidal expansion of asymptotic cone angle Bo. This particular geometry permits an exact solution for a given set of incident modes in the cylindrical duct [5]. Numerical simulation of this solution involves substantial modification to the programs used to obtain the results already presented. A modal boundary condition must be applied at the inlet plane (AB in Figure 8) to permit the generation of reflected waves. An angular dependence e’“’ (4 is an azimuthal co-ordinate as indicated in Figure 8) must also be included throughout the solution to model the presence of spinning nodes. In order to simulate realistic frequencies and spinning mode numbers it also becomes necessary at this stage to employ more efficient higher order elements in place of the simple linear elements of the previous analyses. The full details at these modifications lie outside the scope of the present paper but results are presented which were obtained by using such a scheme. Conceptually the conventional, wave envelope and infinite element schemes retain their significance. A sequence of computed conventional, wave envelope and infinite element sound pressure level contours for the case e. = 70°, ka = 14.5 (k denotes the reduced frequency as before, and a the duct radius as indicated in Figure 8) and a spinning mode number of

62

R. J. ASTLEY

AND

W. EVERSMAN

_ I(b)

7

I

_I

1

!I

wovelenqth

0

3 I

, /

$2

i’ /’

1i 6

,’ // /

I .I

4

2

,

/I

/

,

/

,

,’

-

Plane

wavelength

Figure 9. Hyperboloid duct. Contours of computed absolute acoustical pressure, @a= 70”, &a = 14.5, m = 8; first mode incident; (a) conventional scheme (rmax= 4~); (b) wave envelope scheme (rmax= 4a, r, = 1.6a); (c) infinite element scheme (rr = 1.6~).

eight are shown in Figures 9(a)-(c). The pressure levels represented by the contours are the same in each case and represent uniform pressure increments. Their obvious agreement with each other is striking and their correspondence with the exact solution of Y. C. Cho is demonstrated in Figure 10 which shows polar plots of the “far field” directivity function f(f3) for the first two incident modes (f(e) is defined by ~(r, /3)-f(0)/r as r + co where p(r, f3) e’“’ is the acoustical pressure amplitude normalized with respect to the

ot/” 0

I

0.1

I

I

,

0.2

0.3

0.4

I

0.6

lfW(

Figure 10. Pressure directivity pattern for hyperboloid duct, B0= 70”, ka = 14.5, m = 8. -‘-, Conventional scheme, rmax= 1.6a ; -, conventional scheme, I,, = 4~; - - -, wave envelope scheme, r,,, = 4~2,rr = 1.6a ; 0 0, exact asymptotic solution (Y. C. Cho).

RADIATION

FINITE

ELEMENT

FORMULATIONS

63

amplitude of the incident mode). For the FE solution f(e) is calculated at the outer boundary of the computational region (denoted by r = rmax in Figure 8). Results obtained from the conventional scheme are shown in Figure 10 for rmax = 4a (corresponding to the solution of Figure 9(a)) and rmax= 1.6a. The latter is included to demonstrate the distortion of the directivity pattern due to near field effects. Wave envelope results are shown for rmax = 4a with conventional elements in the region r < rl = 1.6~. These results correspond to the solution represented by the pressure contours of Figure 9(b). The correspondence between the computed and the exact asymptotic directivity patterns is seen to be excellent for both the conventional and wave envelope schemes. The infinite element scheme of Figure 9(c) does not predict the far field in a meaningful way and is not represented in Figure 10. Although these values of frequency and spinning mode number are still relatively low for real inlets (values of ka in excess of 50 and spinning mode numbers of 30 or more would not be unrealistic) the solutions presented do demonstrate the utility and accuracy of the numerical schemes for a configuration which at least approaches that of the real situation.

4. DISCUSSION

AND CONCLUSIONS

The results presented in this paper for simple one dimensional and axisymmetric test cases enable some tentative conclusions to be drawn regarding the general application of infinite and wave envelope elements to acoustical problems. Within obvious limitations, both methods appear to be effective in dealing with the radiation boundary condition. The infinite element approach is inappropriate in itself for the calculation of far field radiation, but is capable of resolving the near field and hence of adequately modelling the radiative effect of the far field. In the context of the turbofan inlet problem infinite elements may well prove a useful technique for situations where only the reflection conditions at the nacelle inlet are required. A more general class of acoustical problems to which such methods might be applied are those of structuralacoustical coupling where a representation of the near field interaction is frequently all that is required. The wave envelope approach appears to incorporate most of the advantages of the infinite element method with the added useful characteristic of providing accurate information on the far field. It has in addition the practical advantage of not requiring a choice of artificial parameters (such as a decay length). A significant computational advantage of the wave envelope scheme lies in the absence from the expressions for element stiffness of infinite limits of integration and of harmonic spatial terms within the integrand. Both of these features ease the problem of numerically integrating such contributions for problems involving more complex field equations. Both methods are clearly a great deal more satisfactory than the straightforward conventional schemes and both appear to have more general application than the boundary element and boundary integral equation methods, with which they may obviously be compared, in that they permit an outer field governed by equations for which simple source solutions are not available (e.g., acoustical problems with mean flow gradients). The above conclusions are tentative in that the test cases presented are not representative of aeroacoustic problems of real practical interest. The current studies do, however, justify further investigation of these methods with more realistic geometries and more complex field equations. Such studies are now in hand.

R. J. ASTLEY

64

AND W. EVERSMAN

ACKNOWLEDGMENTS The research program reported here has been supported by the NASA Langley Research Center under grant NAG-1-198. The authors would also like to thank Dr Y. C. Cho for his generous assistance in providing comparison results.

REFERENCES and C. RAYL 1979 American Institute of Aeronautics and Astronautics Paper 79-0581. The influence of the inlet contour on forward radiated fan noise. 2. A. CRAGGS 1976 Journal of Sound and Vibration 48, 377-392. A finite element method for

1. D. SLOAN, B. W. FARQUHAR

damped acoustic systems: an application to evaluate the performance of reactive mufflers. 3. Y. KAGAWA, T. YAMABUCHI and T. MORI 1977 Journal of Sound Vibration 53, 357-374. Finite element simulation of an axisymmetric acoustic transmission system with a sound absorbing wall. 4. LORD RAYLEIGH 1877 Theory of Sound, Volume 2 (reissued Dover Publications, 1945). 5. Y. C. CHO 1979 American Institute of Aeronautics and Astronautics Paper 79-0677. Sound

radiation from a hyperboloidal inlet duct. 6. H. LEVINE and J. SCHWINGER 1948 Physical Review 73, 383-406. On the radiation of sound from an unflanged circular pipe. 7. L. A. WEINSTEIN 1969 Theory of Diffraction and the Factorization Method, (Generalized Weiner-Hopf Technique). Boulder, Colorado: Golem Press. and J. A. LORDI 1975 Journal of Sound and Vibration 41, 283-290.

8. G. F. HOMICZ

A note on the directivity of duct acoustic modes. 9. H. S. CHEN and C. C. MEI 1975 Proceedings of the Modelling Techniques Conference, (Modelling 1975), San Francisco, 3-5 September 1975, 1, 63-81. Hybrid element for water waves. 10. J. C. W. BERKHOFF 1975 in Finite Elements in Fluids (editor R. H. Gallagher). London: Wiley. Linear wave propagation problems and the finite element method. 11. Y. KAGAWA, T. YAMABUCHI, T. YOSHIKAWA, S. OOIE and N. KYOUNO 1980 Journal of Sound and Vibration 69, 207-228. Finite element approach to acoustic transmissionradiation systems and application to horn and silencer design. 12. S. J. HOROWITZ, R. K. SIGMAN and B. T. ZINN 1981 American Institute of Aeronautics and Astronautics Paper 81- 1981. An iterative finite element technique for predicting sound radiation from turbofan inlets. 13. P. BETTESS and 0. C. ZIENKIEWICZ 1977 International Journal of Numerical Methods in Engineering 11,1271-1290. Diffraction and refraction of surface waves using finite and infinite elements. 14. K. J. BAUMEISTER 1977 NASA Technical Paper 1001. Finite difference theory for sound propagation in a lined duct with uniform flow using the wave envelope technique. 15. R. J. ASTLEY and W. EVERSMAN 1981 Journal of Sound and Vibration 76, 595-601. A note on the utility of a wave envelope approach in finite element duct transmission studies.