- Email: [email protected]
A comparison of numerical collocation and variational procedures to the hypersingular acoustic integral operator
Computer Methods in Applied Mechanics and Engineering North-HoIland
101 (1992) 5-26
CMA 298
A comparison of numerical collocation and variational procedures to the hypersingular acoustic integral operator R. Jeans and I.C. Mathews Imperial College
of Tec~nolog~l,
Science and medicine, department Prime Consort Road. London SW7 2RY. UK
of Aeronautics,
Received 14 October 1991
A major difficulty in implementing many of the more innovative acoustic formulations for the solution of submerged three-dimensional thin shells of arbitrary geometry in an unbounded fluid medium has been the difficulty of numerically approximating a hypersingular integral operator. A variational approach given by Hamdi reduces this hypersingularity/Cauchy singularity to a weak singularity allowing for its accurate numerical approximation. Recently, Mariem and Hamdi [A new boundary finite element method for ~uid-st~cture interaction, Internat. J. Numer. Methods Engrg. 24 (1987) 1251-12673 presented an acoustic formuIation for thin shells, where the pressure difference across the shell is reiated to the shell’s normal velocity solely in terms of this integral operator. Despite its advantage of generating a symmetric fluid matrix, a double surface integration is required to form this matrix, some researchers have endeavoured to develop a collocation method of approximating this highly singular operator. In this paper, a novel solution to the numerical approximation of this operator is given. In the implementation of both the variational and collocation methodologies, quadratic quadrilateral elements are used to interpolate for the acoustic and geometric variables. Results are then presented for thin shell scattering and radiation problems for both closed and non-closed thin shells demonstrating the validity and accuracy of both methods.
1. Introduction
In this paper a collocation method is described for the solution of radiation and scattering problems for a submerged body of arbitrary geometry. A new way of numerically integrating the hypersingular kernels that occur in the boundary integral formulation is presented, and the method is shown to be independent of the interpolation used for the fluid and geometrical variables. This method is then compared with the author’s implementation of a variational method for the hypersingular operator Nk given by Mariem and Hamdi [l]. The variational method not only provides an elegant solution to the problem of the hypersingular kernel, its implementation results in the formation of a symmetric matrix equation. The numerical implementation is also independent of the type of numerical interpolation used. Correspondence to: R. Jeans, Department Consort Road, London SW7 2BY. UK. 004%7825/92/$05.00
0
of Aeronautics,
Imperial College of Science and Technology,
1992 Elsevier Science Publishers B.V. All rights reserved
Prince
6
R. Jeans, f.C. Mathews, A comparison of numerical coi~ocation and variational procedures
A recent acoustic formulation applicable for open thin shells developed by Terai [2] was used to examine the collocation and variational procedures for providing numerical approximations to the hypersingular operator JV~.
2. Acoustic problem Solutions are sought to the exterior pressure field generated by a submerged shell S with interior volume D surrounded by a fluid of density p in the exterior region E. The exterior pressure field is assumed to have a harmonic time dependence of e-I@“. The governing equations are defined by reduced wave equation ,
(V2 + k2)4 = 0 )
dS = 0 , radiation condition ,
!&/[z-ik$)’ 84 = an
(1)
Neumann
iwpv, ,
(2)
boundary condition.
(3)
In these equations, # is defined as the exterior pressure field and u, is the normal velocity of the shell. The surface and the differentiated surface Helmholtz exterior integral equations for the exterior domain are
(4) (5) where the integral operators
.L$, Ak and JV~are defined by
I lib(q)G,(~>4) dS,= =%h#W), I 4(q)aGAp 4) dSq= J4M(P) 4) dSq=ul\T,[cbl(P) > I 4(q) a2GAp, q s
s
s
an
an Y
7
an
P
with the free-space
Green’s function Gk( p, q) given by
G,(P,~)=&‘~‘, r=lp--qj.
(7)
3. Thin shell formulation A prime motivation for this work was the analysis of thin shell acoustic problems. In this paper, a thin shell is defined as a shell for which the through shell displacement field is
7
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variationalprocedures
assumed to be constant. both sides of the shell.
Acoustically,
this implies that the surface normal velocity is equal on
3.1. Boundary integral formulation The geometry of the thin shell is illustrated in Fig. 1. On such a shell three closely associated points may be defined. The points are p, p- and p+, where p represents a point midway through the thickness of the shell, p+ is a point on one surface, and p- is the corresponding point on the other surface. The normal n, is defined to be in the direction from p- to b’. The Green’s function and the normal derivative of Green’s function at these points will have the following simple relationships:
ad+?+)= ww) an,
an9
’
Using these relationships, the surface Helmholtz exterior integral equation formulation thin shell problem may be written as
4(P) = ##‘) +
Is+ ($(q+)
act;
q’)
-
+
li s_
4w)
9+
aG&T 4-l _ aw-)
an _
an,-
G,(P, q+)) dS,+
“;:+I
9+
for the
‘TAP, 4t))
dS,- -
4
By using the relationships in (8) and by setting a(p) = ~#(p+) - &( p-), the surface and the differentiated surface Helmholtz exterior integral equation formulations for the thin shell become
4w = hPl(p) W,(P) -=jlr,lQi](P)++,
an
+ #@> , a 4(P)
PEE,
(10)
PEE.
(11)
The surface domain S is now the surface of one side on the thin shell. Since S has been E
Fig. 1. The thin shell geometry.
8
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variational procedures
redefined, the domains D and E also need to be redefined. The domain D becomes the domain whose interface with S contains the points p-, and E becomes the domain whose interface with S contains the points p+. Both D and E are sub-domains of the exterior domain surrounding S. Taking the limit of P to p in (lo), (ll), (12) W(P)
~ an
W( P>
=
JCPWP)+ +---
PES.
>
(13)
P
3.2. Edge conditions When considering thin shells which do not enclose an interior volume, it is important to consider the behaviour of the integral equations at the edge of this shell. It has been proposed in previous work [3] that by taking the limit of P to p in from D and E, results in two equations that give a condition on @. However, if the limits are taken correctly then
1imJuk[@l(P) = 4-(P)
=
Jkl,c@w)+ (1- c+tP)>($+(P)- #-m
7
(14)
lim JukPwP) = h-(P)
=
JukPw? - (1 -
7
(15)
P-tpf
P-+p
-
C-.(0(4+(P)
- 4+(P))
which leads to the identity (C+(P) + C-(P)
- W(P)
= 0*
(16)
This equation simply states the fact that c_(p) = 1 - c+(p) and gives no condition on @. A more sophisticated argument is needed to define the edge conditions for the open plate problem. When the thickness of the plate is taken to zero, the edge around the plate becomes an additional boundary. Consequently there needs to be a supplementary boundary condition specified on this boundary. In his recent paper Martin [4] discusses the behaviour of one-dimensional hypersingular integral equations over finite intervals and this idea is discussed in more detail. This edge boundary condition is arbitrary in the mathematical sense, but for this case is governed by the original physical problem. Continuity of the pressure difference across the plate means that @ must be zero at the edge of the shell; cP(p)=O,
p on theedge.
(17)
In past work [2,3], this edge boundary condition has been satisfied by the choice of appropriate fluid basis functions. In the numerical work described later, this boundary condition is implicitly imposed upon the numerical formulation. 4. Numerical implemen~tion
of the collocation procedure
It is the Nk operator that needs the greatest amount of consideration when implementing a numerical boundary element method since there are high order singularities that need
9
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variational procedures
integrating. Based on a proof shown by Maue [5], it is possible to convert the kernel of this operator into tangential derivatives using the basic vector identity,
hp *V,)(n, *VJ = (np*n,)(V, By integrating
by parts, the transformed
4Ml(~) =
*v>- (nqx vq>* (n, x Vq). expression for the .Nk operator
(18) is given by
qM(q) - (n, x Y@(q))*(n, x Y$W? j-s Kn, *qJk2G,(P,
q)]dS,
The second term on the right-hand side of (19) can be transformed by Stokes’ theorem into an integral around the edge of S. Therefore, for closed surfaces it is equal to zero and can be neglected. The transformed expression for the Nk operator is finally given by
Jwl(~> =
j-sn,)k’G,(E Kn, *
qMq) - (nq x V&4 4))- hp
x V,)G,(J’~dl ds, . (20)
The validity of this expression has been shown in greater detail by Stallybrass [6] and Mitzner [7]. The JV~operator may be written in the following way:
where the first operator on the right-hand side contains terms of order O(r-‘) and the second terms of order 0(F2). Using the element shape functions, (21) can be approximated by the following expression:
4kW)
=
5{Isb,
j=l
-
-
. n,)~2G,(f’, q)W’}’ (nq x VqUOt)*(y,
x V,W,U’~ s>- G,(C s)>ldS,}b#d (22)
The O(r-$) singularity in the second term of this expression is ignored for the time being and (22) is evaluated at each collocation point to assemble the following matrix equation:
This represents the matrix equation for the JY~operator, and no further modification is made to integrate the singularity occurring in the [NJ term. It will be shown that there is cancellation of the inaccuracies due to the integration of these terms.
10
R. Jeans, I. C. Mathews, A comparison of numerical collocation and variational procedures
Consider the singular integration over an element for the NC,operator. The singularity is separated from the rest of the element by a disk segment of small radius 8. The singular integration over the element can be separated into a non-singular integration over the element minus the disk segment, Ii, and a singular integration over the disk segment, 1:: (nq x I’@( 4)). (n,, x V,G,(P, q)) dS, = Z; +
Ii =
Z; .
(24)
The integration 1: is assumed to be accurate within the limits of the numerical integration scheme. The inaccuracies for the element integration are assumed to be contained within the integration Z:. The small disk segment is assumed to be flat and the vector i is assumed to be in the plane of the disk, perpendicular to the normal. The gradients within (24) can be expressed as
The integration
in the second component
I.=j6;j 2 I
0
' a#(r,e)
0
ar
=,(r,
@>
of (24) can be rewritten rj.7/drde
as
,
dr
where the Jacobian I.!] can be taken as constant. Within the disk segment, it is assumed that (12 has a Iinear dependency on the global position vector and so in (26), the variation of 4 is independent of r and given by
w (4 dr
(27)
=cosO($)~=o+sin8($)Y=o. r = 0
Using this relationship,
Zj can be represented
by
For all elements j that contain the singularity, there will be similar errors in ‘r due to the inadequate singular integration. However, since the summation of the 1: terms is equal to zero, these errors cancel each other out. Therefore it is safe to evaluate the singular terms in the [No] matrix using the Lachat and Watson inverse distance singularity scheme [S]. In a recent paper by Wu et al. [9], a very similar collocation method has been described. This shows that Maue’s equation needs C’ continuity at the collocation points and to achieve this on Co continuous elements, the collocation points are put inside the elements to form an over-determined set of equations. To achieve the integration of the Cauchy principal value integral, they use additional regularization. This study uses a much more simple method. By integrating (22) for singular elements directly, it is recognized that there will be errors in the matrix approximation of the No operator. However, this work has shown that there will be cancellation of these errors when the acoustic problem is solved and so the sophisticated treatment of (22) in the work by Wu et al. is unnecessary.
11
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variationalprocedures
5. Variational boundary integral formulation The boundary
integral equation
problems
is defined by
ac;bx P>
G(P)
____ an
for acoustic thin shell scattering
=
JUWP) + ---&-P ,
where IQ,is the pressure difference across the shell. Following the method shown above, the Nk operator can be transformed into tangential derivatives:
For the variational formulation, p over the shell surface S:
(30) is multiplied by 6@(p) and integrated
S@(p)Nk[@](p)
dSP +
a44 P) S@(P) c dSP . P
The expression for the Nk operator given in (30) can be further transformed by parts to obtain the expression
-(npx’PpfwP))*(
n4
with respect to
x Y$%NlG,(p~d
(31)
by integrating
dSqdSP-
(32)
Equation (32) contains singularities of O(r-‘) and consequently it is possible to integrate this equation accurately using the modified integration scheme of Lachat and Watson. In arriving at (32), integrations around the edge of the thin shell are discarded. The justification for closed thin shells is that there are no edge boundaries and consequently these integrations must be equal to zero. For non-closed thin shells this is not the case but the pressure difference around the edge of the shell must be equal to zero. Therefore, for these shells the edge integration must also be equal to zero.
6. Numerical implementation of the variational formulation The discretization of the variational formulation is implemented using the same order of interpolation as that for the collocation method. The XL operator can be numerically approximated by
-
(np x ${N;})-
(n, x V,{~“,~‘W,(p, 91ds,
dSp){@? . (33)
12
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variational procedures
In (33), the gradient of the shape function is expressed in terms of the local coordinate axis system. In order to evaluate these gradients, it is necessary to use the Jacobian matrix relationship relating the local and curvilinear coordinate systems. This Jacobian matrix enables the gradient to be expressed in derivatives of the curvilinear system. For the integration in (33), the point p ranges over S, and the point q ranges over Tj: The integration with respect to p is performed using simple Gauss integration. The singularltles in the above numerical approximation will occur at these Gauss integration points when i = j. It is these auto-influence elements that need to be integrated using the Lachat and Watson singular integration scheme. For the integration with respect to q, the element, will be divided as shown in Fig. 2. The expression in (32) can be assembled into a matrix formulation of the variational problem,
(34) where the matrix [A] represents
the numerical approximation,
(35) In practice, matrix.
this banded
symmetric
matrix may be further
approximated
by the diagonal
a;; = .~Nf’( q) dS, . i
The matrix equation in (34) is symmetric and is factorized into an LL’ form using the standard Choleski factorization technique. The pressure differences are then evaluated by means of forward and backward substitution. The order of Gauss integration used to calculate the numerical results in this paper are summarized in Table 1. The pressure distribution exterior to the thin shell is given by
4(P) = JuJ@](f’)
3
(37)
I’ E E >
y
Fig. 2. Subelement
Singulari@ at Gauss point
division.
R. Jeans, I. C. Mathews, A comparison of numerical collocation and variational procedures
13
Table 1 Order of Gauss integration Element
p integration
q integration
Singular Non-singular
3x3 2x2
4X4X4” 2x2
a Integration
within subelement.
and the far field pressure distribution
7. Edge boundary
is approximated
by the vector expression
conditions
For non-closed thin shells, it is necessary to enforce an edge physical reality of this is that the pressure difference across the shell Since this boundary condition is not implicit in the collocation approximation of the .Nk operator, then it must be imposed on the set. This matrix equation set can be represented by
boundary condition. The must be zero at the edges. or variational numerical resulting matrix equation
where [H] is the matrix approximation of the collocation or variational formulation the appropriate right-hand side of this matrix equation:
Defining {GE} and {@r} as those elements of the vector that correspond nodes, (39) can be rewritten as
and { y} is
to edge or interior
(41) The interior pressure differences can now be written in terms of the known edge pressure differences and the known vector {y,},
If @ is equal to zero on the edge of the plate, this equation
reduces to
(43)
R. Jeans, Z.C. Mathews, A comparison of numerical collocation and variational procedures
14
Imposing the edge boundary condition in this way is equivalent to making the shape function for an edge node equal to zero so that there is no contribution to the interpolated pressure difference within the element from this node. Consequently the way the pressure difference goes to zero depends on the other shape functions within the edge element. This problem has not been examined extensively by earlier authors [l-3]. They have used constant value elements and for these elements there are no edge nodes and the edge boundary condition is satisfied when distributing out the pressure differences from the element constant value to the nodal values. The numerical results in this paper show that the method described above satisfies the edge boundary condition adequately. However, it was found that there are larger numerical errors at these edge points due to the imposition of edge boundary condition. One possible improvement not implemented in this study would be the modification of the shape functions in the edge elements, so that the edge boundary condition is included in the numerical formulation implicitly and thus more efficiently.
8. Numerical results Both the collocation and variational methods were implemented within a UNIX environment using FORTRAN. The two methods share a lot of common code concerning the element geometry and far field calculations. The discretization used for this study was based on 9-noded quadrilateral isoparametric elements. The extra integration for the variational method places an extra burden on the computational problem, which is significant for the small problems considered in this study, Despite this, both the collocation and variational methods have an O(n2) time dependency. The matrix factorizations, however, have time dependency of O(n’). A comparison between the collocation and variational methods applied to the spherical radiation problem with 6 elements and constant normal velocity is shown in Table 2. The accuracies and timings are tabulated for a range of integration schemes. The times for matrix assembly and factorization are tabulated in Table 3 and plotted in Fig. 3. These results clearly demonstrate the time dependencies of factorization and assembly, and show the increased burden of the extra integration for the variational method. However, the difference between the variational and collocation assembly times is not as significant as might be expected and for large problems with n >> 100, the time limitation will be the matrix factorization, for which the variational method has the advantage. 8.1.
Spheroids
The thin shell formulation is used to calculate the backscattered form function for different spheroidal geometries with an end on incident plane wave. The form function is given by
For these
closed geometries
there
will be problems
of uniqueness
at interior
resonant
R. Jeans, I.C. ~athew~~ A co~~ari~o~ of numerical collocation and variational procedures
Table 2 Comparison
of accuracy with different orders of integration
Method
Non-singular
Singular
Collocation
2x2
2x2 3x3 4x4 2x2 3x3 4x4 2x2 3x3 4x4
219 376 579 365 509 723 589 730 937
1.55 0.74 1.05 0.68 0.08 0.23 1.04 0.30 0.60
2x2 3x3 4x4 2x2 3x3 4x4 2x2 3x3 4x4
333 509 732 1157 1550 2204 3111 3847 4942
1.38 0.65 0.23 0.65 0.81 0.43 0.32 0.43 0.56
3x3
4x4
Variational
15
2x2
3x3
4x4
Timing (s x lo-“)
Accuracy (%)
frequencies [lo] and no method is implemented to remove these problems. This is likely to account on the breakdown of the numerical results with respect to frequency since the range of wavenumbers over which the numerics problem is ill-conditioned at a critical wavenumber increases with frequency until there is overlap. Loss of accuracy at low frequencies is unlikely to be evident unless the wavenumber is very close to the critical wavenumber.
Table 3 Comparison
of matrix assembly and matrix factorization
times CPU timing (s X 10F2)
Method
Elements
Nodes
Assembly
Collocation
3 6 12 24 48 96
19 33 61 113 217 417
172 367 1039 3071 9680 36708
1 3 7 49 457 5598
173 370 1046 3120 10137 42306
Variational
3 6 12 24 48 96
19 33 61 113 237 417
205 518 1437 4348 15258 58364
0 1 10 24 179 1205
205 519 1447 4372 15437 59569
Factorization
Total
16
R. Jeans, I. C. Mathews, A comparison of numerical collocation and variational procedures
CPU
Seconds x10-2
*o10
Number of Elements
Fig. 3. Graph showing matrix assembly and matrix factorization methods.
lco
times for the variational
(v) and collocation
(c}
Figure 4 shows the results calculated using both the variational and collocation method compared to the analytical solution. The 6 element results show high accuracy at low frequencies with a breakdown of both numerical solutions starting at ka = 2.5. There is no significant difference in accuracy between the two methods. The 24 element results show no breakdown of the solutions below ka = 5 but there is slightly better accuracy for the variational results at the higher frequencies. Figures 5 and 6 show the two methods compared to the results generated by Wu’s axisymmetric method [ll] using 20 high order elements applied to prolate and oblate spheroids, respectively. In the case of the prolate spheroid, the breakdown of both numerical results occurs at the same order of frequencies, and there is no significant difference in accuracy between them. For the oblate spheroid, there is no distinct breakdown of the numerical methods for both mesh densities. The 24 element case shows high accuracy for all frequencies. For the 6 element case, the accuracy of both methods degenerates at about ku = 3.0 with the variational method showing slightly higher accuracy. 8.2. Flat disk Figures 8-11 show the numerical calculation of the dimensionless radial surface pressure on a circular disk radiating with constant normal velocity without a baffle. The results are calculated using the mesh geometries shown in Fig. 7 using both the variational and collocation methods. For a flat continuous surface with no incident wave, the surface pressure on the disk is given by
M%(P)+ k(P)) = 4cPw4
=0
’
(45)
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variational procedures
0.8
f 0.8
a/b=1 6 Elements
0.0
1
0
2
3
ka
4
ka
5
1.0 -
0.8 -
f
,
a/b=1 24 Elements
0.0
0 Analytical
__
I
I
I
I
1
2
3
4
Variational
--*-**--- ,
Collocation
----
.
Fig. 4. The far field backscattered form function for a sphere discretized into 6 elements and 24 elements. numerical results are calculated using the variational and collocation thin shell formulations.
Therefore
the dimensionless P=
ik@ 2(&#&n),
’
The
surface pressure is given by
(46)
where (~#/~~)” represents the constant normal surface velocity. The numerical results are compared to results extracted from work by Weiner [12], which were calculated from diffraction data published by Leitner [13]. Similar results have also been published by Wu et al. [14] who use an axisymmetric variational procedure.
18
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variational procedures
0.6
a/b=O.S 6 Elements
V.”
0
1
2
3
ka
4
ka
4
1.2
0.8 -
.
f
0.6 -
0.4 -
a/b=OS 24 Elements
0.2 0.0 /.,
,
0 s.w.wu
1 -.
,
*
2
Variational
,
3
---------.
,:
Collocation
----
.
Fig. 5. The far field backscattered form function for a prolate spheroid discretized into 6 elements and 24 elements. The numerical results are calculated using the variational and collocation thin shell formulations with an end on incident plane wave.
All results show convergence between 20 and 80 elements and there is little that separates the accuracy between the collocation and variational methods. Overall, there is a high degree of accuracy. It is however, clear from these graphs that the accuracy is significantly lower at the edges of the disk. This edge inaccuracy seems not to decrease as the density of elements is increased. Figure 12 shows the radiation reactance and resistance of the disk calculated using both the variational and collocation methods. The radiated power and rate of radiation of kinetic
R. Jeans, Z.C. ~ut~ew~, A comparison of numerical collocation and variational procedures
19
a/b=2 6 Elements
s.w.wu
-.
Variational
-*--..*.*
.
Collocation
----
.
Fig. 6. The far field backscattered form function for an oblate spheroid discretized into 6 elements elements. The numerical results are calculated using the variational and collocation thin shell formulations end on incident plane wave,
and 24 with an
energy for the disk are given by
(47) The radiation
impedance
“=a;--icr,=-,
is defined as II K
(48)
20
R. Jeans,
I.C. curlews,
A comparison of numerical collocation and variational procedures
(a)
20 element
(c) 80
element
disk
(b) 28 element disk
disk
(d) 96 element
(e,J 24 element Terai problem
Fig. 7. The non-closed
plate
p&e
thin shell mesh geometries.
where a, is the positive radiation resistance, representing radiation damping, and q is the radiation reactance that is usually positive; representing added fluid mass rather than fluid damping. The numerical results in Fig. 12 compare well to the results of Weiner and Bouwkamp [El, with a slightly higher degree of accuracy with the variational method. The resistance results from the work by Bouwkamp were calculated for the complimentary aperture. The final results for the flat disk shown in Fig. 13 show the convergence of the radial pressure amplitude for the three mesh discretizations. The 28 element mesh is distinguished by the high density of elements near the edge, however all the mesh geometries show similar inaccuracies at the edge of the plate. To improve the accuracy at the edge of the disk, it seems that a more sophisticated way of imposing the edge boundary condition needs to be developed. 8.3.
Flat square plate
Figure 14 shows the radiation resistance and reactance calculated for the square plate Of side length L. The numerical results are calculated using both the variational and collocation methods and are compared to the results extracted from work by Warham [3].
R. Jeans, I.C. Mathews, A comparison
of numerical
collocation and variational procedures
21
0’1
8'0 9.0
VO 2'0 0’0
olu!paz!laJx!pysipaql30 JalJienb auo ylf~ pale~n~~~?~ aJaM Sl[nSaJ UO!lWOl[OZ7 aq& -9B3Jcq E lIIOql!M kl!30~% ~EUIJOU ~UI?lSUO3Yl!M %J!le!peJyup JeIn3J!3E uo apnlgdruv aJnssaJdpz!peJssa~uo!suaur~p ago, .II %A
VJ
0’ 1
olu!pazyarwp ys!payi30 JalJenbauo ql!~ palepu1w aJlaMslInsaJ uoy?~o[~osaqL .au3Jeq e lnoql!MdlI30Iaa ~t?UIJOU lWlSUO2 ql!M %I!ltXpeJyS!p JeInDJ!D1! UO apnlydwe aJnssaJd[e!peJssa~uo!suaur~p aye_ '01 'fh.3
I.
.._.
...”
..
.
.
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variational procedures
d II
m Y
23
24
R. /cans, 1.C. Mathews, A comparison of numerical collocation and variational procedures
7
8
R. Jeans, l.c.
Mathews, A comparison of numerical collocation and variationalprocedures
25
There is good agreement using all methods but the graph shows that there is higher accuracy with the variational method. The results from Warham correspond to his results for an asymptotically fine mesh. Comparison with the results for his ‘practical’ mesh shows that there is a significantly higher degree of accuracy and rate of convergence for both the collocation and variational methods. This would be expected since Warham’s method uses a piecewise constant element discretization. 8.4. Terai’s problem Terai’s paper [2] provides a comparison between his numerical results and an experimental result. Figure 15 shows the comparison between the collocation and variation results and the measured result extracted from Terai’s work. The results are the pressure gain in dB around a rectangular plate due to a point source. The results are calculated at a radius of 0.31 from the origin, in the X-Z plane, and the source is at a distance of 0.5 from the origin along the 2 axis. The rectangular plate in the X-Y plane has dimensions 0.3 x 0.2. Both the collocation results and the variational results agree well with the measured values. Also they show a significantly higher degree of accuracy than Terai’s results.
9. Conclusions For the numerical acoustic analysis of arbitrary three-dimensional structures, the boundary element method is becoming established as the preferred solution technique. In this study, a collocation and variational BEM have been presented. Both methods make use of isoparametric elements but unlike many previous BEM they are independent of the order of interpolation used. The dif~culty in implementing many acoustic formulations has been the accurate evaluation of the hypersinguIar integral operator. Both the collocation and variational formulations presented make use of the conversion to tangential derivatives first presented by Maue [5]. After conversion to tangential derivatives, there still remains a Cauchy type O(r-‘) singularity. In order to implement a collocation form of the hypersi~gular integral operator, an argument for essentially ignoring the degree of this singularity has been presented. It was shown that for a collocation point on the surface with strict C, continuity, there will be cancellation of errors in the assembled matrix equation. It is assumed that with just C, continuity, there is still sufficient continuity of first order derivatives to ensure numerical accuracy. The resulting collocation .approximation of the fik operator is used in a thin shell acoustic formulation. The thin shell acoustic problem is also formulated using the variational procedure of Mariem and Hamdi [l]. This formulation has the advantage that the hypersingular operator is transformed into an expression containing only weak singularities and the resulting matrix problem is symmetric. The disadvantage of the method is the increased computational burden of the extra integration. The collocation and variational methods are compared by applying to closed and open thin shell problems. Both methods show good agreement with established results, using relatively coarse mesh discretizations. The variational method shows slightly better accuracy than the collocation method, however the difference is small. The results for
26
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variational procedures
the plate problems indicate that the treatment of the edge boundary condition for these isoparametric boundary element methods is a source of error. The comparison of the collocation and variational methods, and consideration of their complimentary computational strengths and weaknesses, leads us to believe that a combination of the two procedures is an area of possible further research. One approach currently being considered by the authors is the formation of a variational matrix formulation through application of the collocation method at the element Gauss points, where by definition there is C, continuity. Such a formulation would combine the speed of the collocation procedure with the symmetry and robustness associated with variational methods.
References [l] J.B. Mariem and M.A. Hamdi, A new boundary finite element method for fluid-structure interaction problems, Internat. J. Numer. Methods Engrg. 24 (1987) 1251-1267. [2] T. Terai, On calculation of sound fields around three dimensional objects by integral equation methods, J. Sound Vibration 69 (1980) 71-100. [3] A.G.P. Warham, The Helmholtz integral equation for a thin shell, NPL Report DITC 129188, 1988. [4] P.A. Martin, End-point behaviour of solutions to hypersingular integral equations, Proc. Roy. Sot. London Ser. A 432 (1991) 301-320. [5] R.D. Maue, Zur Formulierung eines allgemeinen Beugungsproblems durch eine Integralgleichung, Z. Phys. 126 (1949) 601-618. [6] M.P. Stallybrass, On a pointwise variational principle for the approximate solution of linear boundary value problems, J. Math. Mech. 16 (1967) 1247-1286. [7] K.M. Mitzner, Acoustic scattering from an interface between media of greatly different density, J. Math. Phys. 7 (1966) 2053-2060. [S] J.C. Lachat and J.O. Watson, Effective numerical treatment of boundary integral equations: A formulation for three-dimensional elastostatics, Internat. J. Numer. Methods Engrg. 10 (1976) 991-1005. [9] T.W. Wu, A.F. Seybert and G.C. Wan, On the numerical implementation of a Cauchy principal value integral to insure a unique solution for acoustic radiation and scattering, J. Acoust. Sot. Am. 90 (1991). [lo] A.J. Burton and G.F. Miller, The application of integral equation methods to the numerical solution of some exterior boundary value problems, Proc. Roy. Sot. London Ser. A 323 (1971) 201-210. [ll] S.W. Wu, A fast, robust and accurate proceadure for radiation and scattering analyses of submerged elastic axisymmetric bodies, Ph.D. Thesis, Department of Aeronautics, Imperial College, London, 1990. [12] F.M. Wiener, On the relationship between the sound fields radiated and diffracted by plane obstacles, J. Acoust. Sot. Amer. 23 (1951) 697-700. [13] A. Leitner, Diffraction of sound by a circular disk, J. Acoust. Sot. Amer. 21 (1949) 331-334. [14] X.F. Wu. A.D. Pierce and J.H. Ginsberg, Variational method for computing surface acoustic pressure on vibrating bodies, applied to transversely oscillating disks, IEEE J. Oceanic Engrg. 12 (1987). [15] C.J. Bouwkamp, Theoretische en numerieke behandeling van de buiging door ronde opening, Thesis, University of Groningen. 1941.
101 (1992) 5-26
CMA 298
A comparison of numerical collocation and variational procedures to the hypersingular acoustic integral operator R. Jeans and I.C. Mathews Imperial College
of Tec~nolog~l,
Science and medicine, department Prime Consort Road. London SW7 2RY. UK
of Aeronautics,
Received 14 October 1991
A major difficulty in implementing many of the more innovative acoustic formulations for the solution of submerged three-dimensional thin shells of arbitrary geometry in an unbounded fluid medium has been the difficulty of numerically approximating a hypersingular integral operator. A variational approach given by Hamdi reduces this hypersingularity/Cauchy singularity to a weak singularity allowing for its accurate numerical approximation. Recently, Mariem and Hamdi [A new boundary finite element method for ~uid-st~cture interaction, Internat. J. Numer. Methods Engrg. 24 (1987) 1251-12673 presented an acoustic formuIation for thin shells, where the pressure difference across the shell is reiated to the shell’s normal velocity solely in terms of this integral operator. Despite its advantage of generating a symmetric fluid matrix, a double surface integration is required to form this matrix, some researchers have endeavoured to develop a collocation method of approximating this highly singular operator. In this paper, a novel solution to the numerical approximation of this operator is given. In the implementation of both the variational and collocation methodologies, quadratic quadrilateral elements are used to interpolate for the acoustic and geometric variables. Results are then presented for thin shell scattering and radiation problems for both closed and non-closed thin shells demonstrating the validity and accuracy of both methods.
1. Introduction
In this paper a collocation method is described for the solution of radiation and scattering problems for a submerged body of arbitrary geometry. A new way of numerically integrating the hypersingular kernels that occur in the boundary integral formulation is presented, and the method is shown to be independent of the interpolation used for the fluid and geometrical variables. This method is then compared with the author’s implementation of a variational method for the hypersingular operator Nk given by Mariem and Hamdi [l]. The variational method not only provides an elegant solution to the problem of the hypersingular kernel, its implementation results in the formation of a symmetric matrix equation. The numerical implementation is also independent of the type of numerical interpolation used. Correspondence to: R. Jeans, Department Consort Road, London SW7 2BY. UK. 004%7825/92/$05.00
0
of Aeronautics,
Imperial College of Science and Technology,
1992 Elsevier Science Publishers B.V. All rights reserved
Prince
6
R. Jeans, f.C. Mathews, A comparison of numerical coi~ocation and variational procedures
A recent acoustic formulation applicable for open thin shells developed by Terai [2] was used to examine the collocation and variational procedures for providing numerical approximations to the hypersingular operator JV~.
2. Acoustic problem Solutions are sought to the exterior pressure field generated by a submerged shell S with interior volume D surrounded by a fluid of density p in the exterior region E. The exterior pressure field is assumed to have a harmonic time dependence of e-I@“. The governing equations are defined by reduced wave equation ,
(V2 + k2)4 = 0 )
dS = 0 , radiation condition ,
!&/[z-ik$)’ 84 = an
(1)
Neumann
iwpv, ,
(2)
boundary condition.
(3)
In these equations, # is defined as the exterior pressure field and u, is the normal velocity of the shell. The surface and the differentiated surface Helmholtz exterior integral equations for the exterior domain are
(4) (5) where the integral operators
.L$, Ak and JV~are defined by
I lib(q)G,(~>4) dS,= =%h#W), I 4(q)aGAp 4) dSq= J4M(P) 4) dSq=ul\T,[cbl(P) > I 4(q) a2GAp, q s
s
s
an
an Y
7
an
P
with the free-space
Green’s function Gk( p, q) given by
G,(P,~)=&‘~‘, r=lp--qj.
(7)
3. Thin shell formulation A prime motivation for this work was the analysis of thin shell acoustic problems. In this paper, a thin shell is defined as a shell for which the through shell displacement field is
7
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variationalprocedures
assumed to be constant. both sides of the shell.
Acoustically,
this implies that the surface normal velocity is equal on
3.1. Boundary integral formulation The geometry of the thin shell is illustrated in Fig. 1. On such a shell three closely associated points may be defined. The points are p, p- and p+, where p represents a point midway through the thickness of the shell, p+ is a point on one surface, and p- is the corresponding point on the other surface. The normal n, is defined to be in the direction from p- to b’. The Green’s function and the normal derivative of Green’s function at these points will have the following simple relationships:
ad+?+)= ww) an,
an9
’
Using these relationships, the surface Helmholtz exterior integral equation formulation thin shell problem may be written as
4(P) = ##‘) +
Is+ ($(q+)
act;
q’)
-
+
li s_
4w)
9+
aG&T 4-l _ aw-)
an _
an,-
G,(P, q+)) dS,+
“;:+I
9+
for the
‘TAP, 4t))
dS,- -
4
By using the relationships in (8) and by setting a(p) = ~#(p+) - &( p-), the surface and the differentiated surface Helmholtz exterior integral equation formulations for the thin shell become
4w = hPl(p) W,(P) -=jlr,lQi](P)++,
an
+ #@> , a 4(P)
PEE,
(10)
PEE.
(11)
The surface domain S is now the surface of one side on the thin shell. Since S has been E
Fig. 1. The thin shell geometry.
8
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variational procedures
redefined, the domains D and E also need to be redefined. The domain D becomes the domain whose interface with S contains the points p-, and E becomes the domain whose interface with S contains the points p+. Both D and E are sub-domains of the exterior domain surrounding S. Taking the limit of P to p in (lo), (ll), (12) W(P)
~ an
W( P>
=
JCPWP)+ +---
PES.
>
(13)
P
3.2. Edge conditions When considering thin shells which do not enclose an interior volume, it is important to consider the behaviour of the integral equations at the edge of this shell. It has been proposed in previous work [3] that by taking the limit of P to p in from D and E, results in two equations that give a condition on @. However, if the limits are taken correctly then
1imJuk[@l(P) = 4-(P)
=
Jkl,c@w)+ (1- c+tP)>($+(P)- #-m
7
(14)
lim JukPwP) = h-(P)
=
JukPw? - (1 -
7
(15)
P-tpf
P-+p
-
C-.(0(4+(P)
- 4+(P))
which leads to the identity (C+(P) + C-(P)
- W(P)
= 0*
(16)
This equation simply states the fact that c_(p) = 1 - c+(p) and gives no condition on @. A more sophisticated argument is needed to define the edge conditions for the open plate problem. When the thickness of the plate is taken to zero, the edge around the plate becomes an additional boundary. Consequently there needs to be a supplementary boundary condition specified on this boundary. In his recent paper Martin [4] discusses the behaviour of one-dimensional hypersingular integral equations over finite intervals and this idea is discussed in more detail. This edge boundary condition is arbitrary in the mathematical sense, but for this case is governed by the original physical problem. Continuity of the pressure difference across the plate means that @ must be zero at the edge of the shell; cP(p)=O,
p on theedge.
(17)
In past work [2,3], this edge boundary condition has been satisfied by the choice of appropriate fluid basis functions. In the numerical work described later, this boundary condition is implicitly imposed upon the numerical formulation. 4. Numerical implemen~tion
of the collocation procedure
It is the Nk operator that needs the greatest amount of consideration when implementing a numerical boundary element method since there are high order singularities that need
9
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variational procedures
integrating. Based on a proof shown by Maue [5], it is possible to convert the kernel of this operator into tangential derivatives using the basic vector identity,
hp *V,)(n, *VJ = (np*n,)(V, By integrating
by parts, the transformed
4Ml(~) =
*v>- (nqx vq>* (n, x Vq). expression for the .Nk operator
(18) is given by
qM(q) - (n, x Y@(q))*(n, x Y$W? j-s Kn, *qJk2G,(P,
q)]dS,
The second term on the right-hand side of (19) can be transformed by Stokes’ theorem into an integral around the edge of S. Therefore, for closed surfaces it is equal to zero and can be neglected. The transformed expression for the Nk operator is finally given by
Jwl(~> =
j-sn,)k’G,(E Kn, *
qMq) - (nq x V&4 4))- hp
x V,)G,(J’~dl ds, . (20)
The validity of this expression has been shown in greater detail by Stallybrass [6] and Mitzner [7]. The JV~operator may be written in the following way:
where the first operator on the right-hand side contains terms of order O(r-‘) and the second terms of order 0(F2). Using the element shape functions, (21) can be approximated by the following expression:
4kW)
=
5{Isb,
j=l
-
-
. n,)~2G,(f’, q)W’}’ (nq x VqUOt)*(y,
x V,W,U’~ s>- G,(C s)>ldS,}b#d (22)
The O(r-$) singularity in the second term of this expression is ignored for the time being and (22) is evaluated at each collocation point to assemble the following matrix equation:
This represents the matrix equation for the JY~operator, and no further modification is made to integrate the singularity occurring in the [NJ term. It will be shown that there is cancellation of the inaccuracies due to the integration of these terms.
10
R. Jeans, I. C. Mathews, A comparison of numerical collocation and variational procedures
Consider the singular integration over an element for the NC,operator. The singularity is separated from the rest of the element by a disk segment of small radius 8. The singular integration over the element can be separated into a non-singular integration over the element minus the disk segment, Ii, and a singular integration over the disk segment, 1:: (nq x I’@( 4)). (n,, x V,G,(P, q)) dS, = Z; +
Ii =
Z; .
(24)
The integration 1: is assumed to be accurate within the limits of the numerical integration scheme. The inaccuracies for the element integration are assumed to be contained within the integration Z:. The small disk segment is assumed to be flat and the vector i is assumed to be in the plane of the disk, perpendicular to the normal. The gradients within (24) can be expressed as
The integration
in the second component
I.=j6;j 2 I
0
' a#(r,e)
0
ar
=,(r,
@>
of (24) can be rewritten rj.7/drde
as
,
dr
where the Jacobian I.!] can be taken as constant. Within the disk segment, it is assumed that (12 has a Iinear dependency on the global position vector and so in (26), the variation of 4 is independent of r and given by
w (4 dr
(27)
=cosO($)~=o+sin8($)Y=o. r = 0
Using this relationship,
Zj can be represented
by
For all elements j that contain the singularity, there will be similar errors in ‘r due to the inadequate singular integration. However, since the summation of the 1: terms is equal to zero, these errors cancel each other out. Therefore it is safe to evaluate the singular terms in the [No] matrix using the Lachat and Watson inverse distance singularity scheme [S]. In a recent paper by Wu et al. [9], a very similar collocation method has been described. This shows that Maue’s equation needs C’ continuity at the collocation points and to achieve this on Co continuous elements, the collocation points are put inside the elements to form an over-determined set of equations. To achieve the integration of the Cauchy principal value integral, they use additional regularization. This study uses a much more simple method. By integrating (22) for singular elements directly, it is recognized that there will be errors in the matrix approximation of the No operator. However, this work has shown that there will be cancellation of these errors when the acoustic problem is solved and so the sophisticated treatment of (22) in the work by Wu et al. is unnecessary.
11
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variationalprocedures
5. Variational boundary integral formulation The boundary
integral equation
problems
is defined by
ac;bx P>
G(P)
____ an
for acoustic thin shell scattering
=
JUWP) + ---&-P ,
where IQ,is the pressure difference across the shell. Following the method shown above, the Nk operator can be transformed into tangential derivatives:
For the variational formulation, p over the shell surface S:
(30) is multiplied by 6@(p) and integrated
S@(p)Nk[@](p)
dSP +
a44 P) S@(P) c dSP . P
The expression for the Nk operator given in (30) can be further transformed by parts to obtain the expression
-(npx’PpfwP))*(
n4
with respect to
x Y$%NlG,(p~d
(31)
by integrating
dSqdSP-
(32)
Equation (32) contains singularities of O(r-‘) and consequently it is possible to integrate this equation accurately using the modified integration scheme of Lachat and Watson. In arriving at (32), integrations around the edge of the thin shell are discarded. The justification for closed thin shells is that there are no edge boundaries and consequently these integrations must be equal to zero. For non-closed thin shells this is not the case but the pressure difference around the edge of the shell must be equal to zero. Therefore, for these shells the edge integration must also be equal to zero.
6. Numerical implementation of the variational formulation The discretization of the variational formulation is implemented using the same order of interpolation as that for the collocation method. The XL operator can be numerically approximated by
-
(np x ${N;})-
(n, x V,{~“,~‘W,(p, 91ds,
dSp){@? . (33)
12
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variational procedures
In (33), the gradient of the shape function is expressed in terms of the local coordinate axis system. In order to evaluate these gradients, it is necessary to use the Jacobian matrix relationship relating the local and curvilinear coordinate systems. This Jacobian matrix enables the gradient to be expressed in derivatives of the curvilinear system. For the integration in (33), the point p ranges over S, and the point q ranges over Tj: The integration with respect to p is performed using simple Gauss integration. The singularltles in the above numerical approximation will occur at these Gauss integration points when i = j. It is these auto-influence elements that need to be integrated using the Lachat and Watson singular integration scheme. For the integration with respect to q, the element, will be divided as shown in Fig. 2. The expression in (32) can be assembled into a matrix formulation of the variational problem,
(34) where the matrix [A] represents
the numerical approximation,
(35) In practice, matrix.
this banded
symmetric
matrix may be further
approximated
by the diagonal
a;; = .~Nf’( q) dS, . i
The matrix equation in (34) is symmetric and is factorized into an LL’ form using the standard Choleski factorization technique. The pressure differences are then evaluated by means of forward and backward substitution. The order of Gauss integration used to calculate the numerical results in this paper are summarized in Table 1. The pressure distribution exterior to the thin shell is given by
4(P) = JuJ@](f’)
3
(37)
I’ E E >
y
Fig. 2. Subelement
Singulari@ at Gauss point
division.
R. Jeans, I. C. Mathews, A comparison of numerical collocation and variational procedures
13
Table 1 Order of Gauss integration Element
p integration
q integration
Singular Non-singular
3x3 2x2
4X4X4” 2x2
a Integration
within subelement.
and the far field pressure distribution
7. Edge boundary
is approximated
by the vector expression
conditions
For non-closed thin shells, it is necessary to enforce an edge physical reality of this is that the pressure difference across the shell Since this boundary condition is not implicit in the collocation approximation of the .Nk operator, then it must be imposed on the set. This matrix equation set can be represented by
boundary condition. The must be zero at the edges. or variational numerical resulting matrix equation
where [H] is the matrix approximation of the collocation or variational formulation the appropriate right-hand side of this matrix equation:
Defining {GE} and {@r} as those elements of the vector that correspond nodes, (39) can be rewritten as
and { y} is
to edge or interior
(41) The interior pressure differences can now be written in terms of the known edge pressure differences and the known vector {y,},
If @ is equal to zero on the edge of the plate, this equation
reduces to
(43)
R. Jeans, Z.C. Mathews, A comparison of numerical collocation and variational procedures
14
Imposing the edge boundary condition in this way is equivalent to making the shape function for an edge node equal to zero so that there is no contribution to the interpolated pressure difference within the element from this node. Consequently the way the pressure difference goes to zero depends on the other shape functions within the edge element. This problem has not been examined extensively by earlier authors [l-3]. They have used constant value elements and for these elements there are no edge nodes and the edge boundary condition is satisfied when distributing out the pressure differences from the element constant value to the nodal values. The numerical results in this paper show that the method described above satisfies the edge boundary condition adequately. However, it was found that there are larger numerical errors at these edge points due to the imposition of edge boundary condition. One possible improvement not implemented in this study would be the modification of the shape functions in the edge elements, so that the edge boundary condition is included in the numerical formulation implicitly and thus more efficiently.
8. Numerical results Both the collocation and variational methods were implemented within a UNIX environment using FORTRAN. The two methods share a lot of common code concerning the element geometry and far field calculations. The discretization used for this study was based on 9-noded quadrilateral isoparametric elements. The extra integration for the variational method places an extra burden on the computational problem, which is significant for the small problems considered in this study, Despite this, both the collocation and variational methods have an O(n2) time dependency. The matrix factorizations, however, have time dependency of O(n’). A comparison between the collocation and variational methods applied to the spherical radiation problem with 6 elements and constant normal velocity is shown in Table 2. The accuracies and timings are tabulated for a range of integration schemes. The times for matrix assembly and factorization are tabulated in Table 3 and plotted in Fig. 3. These results clearly demonstrate the time dependencies of factorization and assembly, and show the increased burden of the extra integration for the variational method. However, the difference between the variational and collocation assembly times is not as significant as might be expected and for large problems with n >> 100, the time limitation will be the matrix factorization, for which the variational method has the advantage. 8.1.
Spheroids
The thin shell formulation is used to calculate the backscattered form function for different spheroidal geometries with an end on incident plane wave. The form function is given by
For these
closed geometries
there
will be problems
of uniqueness
at interior
resonant
R. Jeans, I.C. ~athew~~ A co~~ari~o~ of numerical collocation and variational procedures
Table 2 Comparison
of accuracy with different orders of integration
Method
Non-singular
Singular
Collocation
2x2
2x2 3x3 4x4 2x2 3x3 4x4 2x2 3x3 4x4
219 376 579 365 509 723 589 730 937
1.55 0.74 1.05 0.68 0.08 0.23 1.04 0.30 0.60
2x2 3x3 4x4 2x2 3x3 4x4 2x2 3x3 4x4
333 509 732 1157 1550 2204 3111 3847 4942
1.38 0.65 0.23 0.65 0.81 0.43 0.32 0.43 0.56
3x3
4x4
Variational
15
2x2
3x3
4x4
Timing (s x lo-“)
Accuracy (%)
frequencies [lo] and no method is implemented to remove these problems. This is likely to account on the breakdown of the numerical results with respect to frequency since the range of wavenumbers over which the numerics problem is ill-conditioned at a critical wavenumber increases with frequency until there is overlap. Loss of accuracy at low frequencies is unlikely to be evident unless the wavenumber is very close to the critical wavenumber.
Table 3 Comparison
of matrix assembly and matrix factorization
times CPU timing (s X 10F2)
Method
Elements
Nodes
Assembly
Collocation
3 6 12 24 48 96
19 33 61 113 217 417
172 367 1039 3071 9680 36708
1 3 7 49 457 5598
173 370 1046 3120 10137 42306
Variational
3 6 12 24 48 96
19 33 61 113 237 417
205 518 1437 4348 15258 58364
0 1 10 24 179 1205
205 519 1447 4372 15437 59569
Factorization
Total
16
R. Jeans, I. C. Mathews, A comparison of numerical collocation and variational procedures
CPU
Seconds x10-2
*o10
Number of Elements
Fig. 3. Graph showing matrix assembly and matrix factorization methods.
lco
times for the variational
(v) and collocation
(c}
Figure 4 shows the results calculated using both the variational and collocation method compared to the analytical solution. The 6 element results show high accuracy at low frequencies with a breakdown of both numerical solutions starting at ka = 2.5. There is no significant difference in accuracy between the two methods. The 24 element results show no breakdown of the solutions below ka = 5 but there is slightly better accuracy for the variational results at the higher frequencies. Figures 5 and 6 show the two methods compared to the results generated by Wu’s axisymmetric method [ll] using 20 high order elements applied to prolate and oblate spheroids, respectively. In the case of the prolate spheroid, the breakdown of both numerical results occurs at the same order of frequencies, and there is no significant difference in accuracy between them. For the oblate spheroid, there is no distinct breakdown of the numerical methods for both mesh densities. The 24 element case shows high accuracy for all frequencies. For the 6 element case, the accuracy of both methods degenerates at about ku = 3.0 with the variational method showing slightly higher accuracy. 8.2. Flat disk Figures 8-11 show the numerical calculation of the dimensionless radial surface pressure on a circular disk radiating with constant normal velocity without a baffle. The results are calculated using the mesh geometries shown in Fig. 7 using both the variational and collocation methods. For a flat continuous surface with no incident wave, the surface pressure on the disk is given by
M%(P)+ k(P)) = 4cPw4
=0
’
(45)
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variational procedures
0.8
f 0.8
a/b=1 6 Elements
0.0
1
0
2
3
ka
4
ka
5
1.0 -
0.8 -
f
,
a/b=1 24 Elements
0.0
0 Analytical
__
I
I
I
I
1
2
3
4
Variational
--*-**--- ,
Collocation
----
.
Fig. 4. The far field backscattered form function for a sphere discretized into 6 elements and 24 elements. numerical results are calculated using the variational and collocation thin shell formulations.
Therefore
the dimensionless P=
ik@ 2(&#&n),
’
The
surface pressure is given by
(46)
where (~#/~~)” represents the constant normal surface velocity. The numerical results are compared to results extracted from work by Weiner [12], which were calculated from diffraction data published by Leitner [13]. Similar results have also been published by Wu et al. [14] who use an axisymmetric variational procedure.
18
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variational procedures
0.6
a/b=O.S 6 Elements
V.”
0
1
2
3
ka
4
ka
4
1.2
0.8 -
.
f
0.6 -
0.4 -
a/b=OS 24 Elements
0.2 0.0 /.,
,
0 s.w.wu
1 -.
,
*
2
Variational
,
3
---------.
,:
Collocation
----
.
Fig. 5. The far field backscattered form function for a prolate spheroid discretized into 6 elements and 24 elements. The numerical results are calculated using the variational and collocation thin shell formulations with an end on incident plane wave.
All results show convergence between 20 and 80 elements and there is little that separates the accuracy between the collocation and variational methods. Overall, there is a high degree of accuracy. It is however, clear from these graphs that the accuracy is significantly lower at the edges of the disk. This edge inaccuracy seems not to decrease as the density of elements is increased. Figure 12 shows the radiation reactance and resistance of the disk calculated using both the variational and collocation methods. The radiated power and rate of radiation of kinetic
R. Jeans, Z.C. ~ut~ew~, A comparison of numerical collocation and variational procedures
19
a/b=2 6 Elements
s.w.wu
-.
Variational
-*--..*.*
.
Collocation
----
.
Fig. 6. The far field backscattered form function for an oblate spheroid discretized into 6 elements elements. The numerical results are calculated using the variational and collocation thin shell formulations end on incident plane wave,
and 24 with an
energy for the disk are given by
(47) The radiation
impedance
“=a;--icr,=-,
is defined as II K
(48)
20
R. Jeans,
I.C. curlews,
A comparison of numerical collocation and variational procedures
(a)
20 element
(c) 80
element
disk
(b) 28 element disk
disk
(d) 96 element
(e,J 24 element Terai problem
Fig. 7. The non-closed
plate
p&e
thin shell mesh geometries.
where a, is the positive radiation resistance, representing radiation damping, and q is the radiation reactance that is usually positive; representing added fluid mass rather than fluid damping. The numerical results in Fig. 12 compare well to the results of Weiner and Bouwkamp [El, with a slightly higher degree of accuracy with the variational method. The resistance results from the work by Bouwkamp were calculated for the complimentary aperture. The final results for the flat disk shown in Fig. 13 show the convergence of the radial pressure amplitude for the three mesh discretizations. The 28 element mesh is distinguished by the high density of elements near the edge, however all the mesh geometries show similar inaccuracies at the edge of the plate. To improve the accuracy at the edge of the disk, it seems that a more sophisticated way of imposing the edge boundary condition needs to be developed. 8.3.
Flat square plate
Figure 14 shows the radiation resistance and reactance calculated for the square plate Of side length L. The numerical results are calculated using both the variational and collocation methods and are compared to the results extracted from work by Warham [3].
R. Jeans, I.C. Mathews, A comparison
of numerical
collocation and variational procedures
21
0’1
8'0 9.0
VO 2'0 0’0
olu!paz!laJx!pysipaql30 JalJienb auo ylf~ pale~n~~~?~ aJaM Sl[nSaJ UO!lWOl[OZ7 aq& -9B3Jcq E lIIOql!M kl!30~% ~EUIJOU ~UI?lSUO3Yl!M %J!le!peJyup JeIn3J!3E uo apnlgdruv aJnssaJdpz!peJssa~uo!suaur~p ago, .II %A
VJ
0’ 1
olu!pazyarwp ys!payi30 JalJenbauo ql!~ palepu1w aJlaMslInsaJ uoy?~o[~osaqL .au3Jeq e lnoql!MdlI30Iaa ~t?UIJOU lWlSUO2 ql!M %I!ltXpeJyS!p JeInDJ!D1! UO apnlydwe aJnssaJd[e!peJssa~uo!suaur~p aye_ '01 'fh.3
I.
.._.
...”
..
.
.
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variational procedures
d II
m Y
23
24
R. /cans, 1.C. Mathews, A comparison of numerical collocation and variational procedures
7
8
R. Jeans, l.c.
Mathews, A comparison of numerical collocation and variationalprocedures
25
There is good agreement using all methods but the graph shows that there is higher accuracy with the variational method. The results from Warham correspond to his results for an asymptotically fine mesh. Comparison with the results for his ‘practical’ mesh shows that there is a significantly higher degree of accuracy and rate of convergence for both the collocation and variational methods. This would be expected since Warham’s method uses a piecewise constant element discretization. 8.4. Terai’s problem Terai’s paper [2] provides a comparison between his numerical results and an experimental result. Figure 15 shows the comparison between the collocation and variation results and the measured result extracted from Terai’s work. The results are the pressure gain in dB around a rectangular plate due to a point source. The results are calculated at a radius of 0.31 from the origin, in the X-Z plane, and the source is at a distance of 0.5 from the origin along the 2 axis. The rectangular plate in the X-Y plane has dimensions 0.3 x 0.2. Both the collocation results and the variational results agree well with the measured values. Also they show a significantly higher degree of accuracy than Terai’s results.
9. Conclusions For the numerical acoustic analysis of arbitrary three-dimensional structures, the boundary element method is becoming established as the preferred solution technique. In this study, a collocation and variational BEM have been presented. Both methods make use of isoparametric elements but unlike many previous BEM they are independent of the order of interpolation used. The dif~culty in implementing many acoustic formulations has been the accurate evaluation of the hypersinguIar integral operator. Both the collocation and variational formulations presented make use of the conversion to tangential derivatives first presented by Maue [5]. After conversion to tangential derivatives, there still remains a Cauchy type O(r-‘) singularity. In order to implement a collocation form of the hypersi~gular integral operator, an argument for essentially ignoring the degree of this singularity has been presented. It was shown that for a collocation point on the surface with strict C, continuity, there will be cancellation of errors in the assembled matrix equation. It is assumed that with just C, continuity, there is still sufficient continuity of first order derivatives to ensure numerical accuracy. The resulting collocation .approximation of the fik operator is used in a thin shell acoustic formulation. The thin shell acoustic problem is also formulated using the variational procedure of Mariem and Hamdi [l]. This formulation has the advantage that the hypersingular operator is transformed into an expression containing only weak singularities and the resulting matrix problem is symmetric. The disadvantage of the method is the increased computational burden of the extra integration. The collocation and variational methods are compared by applying to closed and open thin shell problems. Both methods show good agreement with established results, using relatively coarse mesh discretizations. The variational method shows slightly better accuracy than the collocation method, however the difference is small. The results for
26
R. Jeans, I.C. Mathews, A comparison of numerical collocation and variational procedures
the plate problems indicate that the treatment of the edge boundary condition for these isoparametric boundary element methods is a source of error. The comparison of the collocation and variational methods, and consideration of their complimentary computational strengths and weaknesses, leads us to believe that a combination of the two procedures is an area of possible further research. One approach currently being considered by the authors is the formation of a variational matrix formulation through application of the collocation method at the element Gauss points, where by definition there is C, continuity. Such a formulation would combine the speed of the collocation procedure with the symmetry and robustness associated with variational methods.
References [l] J.B. Mariem and M.A. Hamdi, A new boundary finite element method for fluid-structure interaction problems, Internat. J. Numer. Methods Engrg. 24 (1987) 1251-1267. [2] T. Terai, On calculation of sound fields around three dimensional objects by integral equation methods, J. Sound Vibration 69 (1980) 71-100. [3] A.G.P. Warham, The Helmholtz integral equation for a thin shell, NPL Report DITC 129188, 1988. [4] P.A. Martin, End-point behaviour of solutions to hypersingular integral equations, Proc. Roy. Sot. London Ser. A 432 (1991) 301-320. [5] R.D. Maue, Zur Formulierung eines allgemeinen Beugungsproblems durch eine Integralgleichung, Z. Phys. 126 (1949) 601-618. [6] M.P. Stallybrass, On a pointwise variational principle for the approximate solution of linear boundary value problems, J. Math. Mech. 16 (1967) 1247-1286. [7] K.M. Mitzner, Acoustic scattering from an interface between media of greatly different density, J. Math. Phys. 7 (1966) 2053-2060. [S] J.C. Lachat and J.O. Watson, Effective numerical treatment of boundary integral equations: A formulation for three-dimensional elastostatics, Internat. J. Numer. Methods Engrg. 10 (1976) 991-1005. [9] T.W. Wu, A.F. Seybert and G.C. Wan, On the numerical implementation of a Cauchy principal value integral to insure a unique solution for acoustic radiation and scattering, J. Acoust. Sot. Am. 90 (1991). [lo] A.J. Burton and G.F. Miller, The application of integral equation methods to the numerical solution of some exterior boundary value problems, Proc. Roy. Sot. London Ser. A 323 (1971) 201-210. [ll] S.W. Wu, A fast, robust and accurate proceadure for radiation and scattering analyses of submerged elastic axisymmetric bodies, Ph.D. Thesis, Department of Aeronautics, Imperial College, London, 1990. [12] F.M. Wiener, On the relationship between the sound fields radiated and diffracted by plane obstacles, J. Acoust. Sot. Amer. 23 (1951) 697-700. [13] A. Leitner, Diffraction of sound by a circular disk, J. Acoust. Sot. Amer. 21 (1949) 331-334. [14] X.F. Wu. A.D. Pierce and J.H. Ginsberg, Variational method for computing surface acoustic pressure on vibrating bodies, applied to transversely oscillating disks, IEEE J. Oceanic Engrg. 12 (1987). [15] C.J. Bouwkamp, Theoretische en numerieke behandeling van de buiging door ronde opening, Thesis, University of Groningen. 1941.