Abstract

The method development for Helmholtz problems can be considered to be an extension of the methods for Laplace problems. The BEMLAP on-line companion manual for the solution of Laplace is also available.

1  Introduction

The Fortran subroutines described here are useful in the implementation of integral equation methods for the solution of the general two-dimensional, the general three-dimensional and the axisymmetric three-dimensional Helmholtz equation,

The Helmholtz integral operators are defined as follows:

For the special case when k = 0 the Helmholtz equation (1) is the Laplace equation. In this particular case the chosen Green's functions are

For each particular case of boundary division, the discrete form of the operators is computed using the subroutines H2LC (two-dimensional), H3LC (three-dimensional) and H3ALC (axisymmetric three-dimensional). The subroutines in this paper are thus useful for the solution of the interior or exterior Helmholtz or Laplace equation via integral equation methods. Examples of situations where the subroutines in this paper could be useful are the problems addressed in references [32], [9], [29], [34], [31], [1], [5], [6], [13], [14], [35], [2], [12], [16], [17], [18], [19], [22], [24]. A full description of the solution of acoustic problems in interior, exterior domains and of the solution to modal problems is given in Kirkup [25]. The objective of this paper is to describe the underlying methods employed in computing the discrete form of the integral operators (2)-(5), to outline the three Fortran subroutines and examine their computational cost, explain how the subroutines may be utilised and to demonstrate the subroutines on the Burton and Miller equation [7]. The subroutines have been written to the Fortran77 standard and employ double-precision arithmetic.

2  Some Properties of the Kernel Functions

The results given in this section are extracted mainly from references [8], [9], and [14]. In the following r = p - q and r = |r|, G_k = G_k(p,q), G₀ = G₀(p,q).

2.1  Derivatives of G₀ with respect to r

In two dimensions we have

2.2  Derivatives of G_k (k ¹ 0) with respect to r

In two dimensions we have

In three dimensions we have

2.3  Expressions for the normal derivatives of G_k

The following expressions hold in both two and three dimensions and for all k:

2.4  Expressions for the normal derivative of r

The derivatives of r with respect to v_p and n_q may be written as follows:

2.5  Expressions for ¶² G₀/ ¶v_p ¶n_q

The following results can be derived from the substitution of (25) and (12),(13) or (14),(15) into (22) with k = 0:

2.6  Asymptotic Properties

In the following results, p, q Î G where G is a surface and G is smooth at p:

3  Discretization of the Integral Operators

The computational method is developed in a similar way to that described in reference [20]. In order to derive the discrete forms of the integral operators (2), (3), (4) and (5), G is approximated by a set of n elements G = åⁿ_{j = 1} DG_j . The boundary function m is replaced by its equivalent on the approximate boundary G. The function is then replaced by a constant on each element. Thus for the L_k integral operator:

The discrete forms are thus defined as follows:

The derivative operator in (35) can always be carried inside the integral. The same is true for the operator in (36) when p does not lie on the element DG_j. Thus we may write:

The special techniques applied in this paper involve `subtracting out' the singularity and evaluating the singular part and remaining regular part separately. The following results are immediate from the asymptotic properties of the kernel functions (28) and (31):

In summary, the evaluation of the integral operators requires a summation of a set of integrand values multiplied by quadrature weights. In the case when p Î DG_j the evaluation of the subtracted out part is also required for the L_k and N_k operators.

4  Evaluation of the Discrete Forms

In this section a list of operations for computing the discrete integral operators is given. This particular programme of computation is given prominence because it needs to be executed for each quadrature point and hence it is the key to determining the computational cost of the overall method. The following computational programme is intended to be optimal. Note that it is assumed that p and v_p are already set.

A: Set q the point on the element

B: Set n_q the unit outward normal to the element at q

C: Compute v_p.n_q

D: Compute r

E: Compute r, r²

F: Compute r³ (in the 3D case)

G: Compute ¶r/ ¶n_q via (23)

H: Compute ¶r/ ¶v_p via (24)

I: Compute

¶r/ ¶n_q *¶r/ ¶v_p

J: Compute ¶² r/ ¶v_p ¶n_q via (25)

K: Compute kr [k*r] and ikr [i*kr]

L: Compute skr [kr*kr] (in the 3D cases)

M: Compute H [Hankel function H₀⁽¹⁾ (kr), H₁⁽¹⁾ (kr) for 2D problems] or E [exp(ikr) for 3D problems].

N: Compute Green's function via (6) or (10) for 2D or via (7) or (11) for 3D.

O: Multiply L_k kernel by weight and add to sum

P: Compute ¶G_k/ ¶r via (16) or (12) for 2D or via (18) or (14) for 3D

Q: Compute value of quadrature weight multiplied by the result of operation N.

R: Compute M_k kernel multiplied by weight and add to sum

S: Compute M_k^t kernel multiplied by weight and add to sum

T: Compute ¶² G_k/ ¶r² via (17) or (13) for 2D or via (19) or (15) for 3D

U: Multiply N_k kernel by weight and add to sum

Apart from operations E and M, each operation can be directly costed in terms of floating-point operations. Operation E computes a square root and operation M computes a Hankel function in two-dimensional problems and a complex exponential in three-dimensional problems. The square root and the complex exponential functions are available in most programming languages, though in some cases it could be beneficial (in terms of computational cost) not to use the standard language functions.

In the two-dimensional case, operation M requires the computation of the spherical Hankel functions or log functions when k is zero. An external module is required for the evaluation of the Hankel function.

In the subroutines H2LC, H3LC and H3ALC the square root function and the exponential and/or Hankel function or log functions that are evaluated at each quadrature point need to be provided as external functions with the identifiers FNSQRT, FNEXP, FNHANK and FNLOG. The freedom to define the functions externally and to choose the quadrature rule allows the user to take full control of the efficiency of the subroutines.

5  General Introduction to the Subroutines

The subroutines have the identifiers H2LC, H3LC and H3ALC, for computing the discrete Helmholtz integral operators for the two-dimensional, three-dimensional and three-dimensional axisymmetric cases. The subroutines' parameter list have the following general form:

SUBROUTINE H{2 or 3 or 3A } LC(

complex wavenumber,

point (p and the unit vector v_p, if necessary),

geometry of the element (vertices which define element),

quadrature rule (weights and abscissae for the standard element),

validation and control parameters,

discrete Helmholtz integral operators (output) ).

5.1  General control of the subroutines

The parameter complex K passes the wavenumber k to the subroutines. More efficient methods are applied within the subroutines if the value of K is classed as zero, real or imaginary. These special cases of the wavenumber are discerned through testing whether the size of the real part or imaginary part is less than or greater than the parameter EK.

The logical variable LVALID allows the user to switch on ( .TRUE.) or off ( .FALSE.) the validation of the input parameters. Setting LVALID= .TRUE. effects on series of tests on the input data. The use of this facility has a small computational cost but could be useful when initially constructing a program that calls the subroutines or to facilitate the diagnosis of an error that appears in the calling program or in the subroutine.

The logical variables LLK, LMK, LMKT and LNK allows the user to choose the specific operators of L_k, M_k, M_k^t or N_k that are required. The LLK, LMK, LMKT and LNK that have .TRUE. values results in DISLK, DISMK, DISMKT and DISNK taking the values of the discrete operators on exit from the subroutine.

5.2  Geometrical Information

Real one-dimensional Fortran arrays define the point p and the points that define the geometry of the element. The arrays have two components in H2LC and H3ALC and three components in H3LC. In H2LC and H3LC, the components of the arrays define the Cartesian coordinates of the points and of the directional vector v_p. In H3ALC, the components define the cylindrical polar coordinates of the points and v_p.

In H2LC and H3ALC two points QA and QB, each one-dimensional real arrays in Fortran with two components, define the element. In H2LC, QA and QB are the edges of the straight line element. In H3ALC, QA and QB define the outer rims of the truncated conical element, the two points lying on the edges of the generator of the element. In H3LC, three points QA, QB and QC, each one-dimensional real arrays in Fortran with three components, are the vertices that define the triangular element.

The normal to the element is computed within the subroutines. In H2LC and H3ALC, the normal is defined to be [ QA(2)- QB(2), QB(1)- QA(1)] normalised. Note that in H3ALC the normal is the normal on the generator of the element. Hence the normal is to the right on the line QA- QB in H2LC and H3ALC. In H3LC the normal is defined to be [ QB- QA]×[ QC- QA] normalised, where × denotes the vector cross product. Hence the normal is presumed to be outward when QA- QB- QC are arranged anti-clockwise around the triangle.

The integrals that define L_k and N_k require special treatment when the point p lies on the element. The logical variable LPONEL is used to inform the subroutine as to whether the point p lies on the element ( .TRUE.) or not ( .FALSE.). In the case LPONEL= .TRUE., methods based on the expressions for the integral operators in equations (39) and (40) are used.

For the computation of the discrete M_k^t and N_k operators the vector v_p needs to be defined. In the subroutines it is denoted VECP. In the subroutines VECP must represent a unit vector. Results for a non-unit vector v_p can be obtained through scaling the vector v_p before calling the subroutine (so that it is a unit vector) and suitably scaling the results DISMKT and DISNK afterwards. If the point p ( P) lies on the element then v_p ( VECP) must be the unit normal to the element, as defined above.

5.3  Quadrature Rule

In subroutine H2LC a quadrature rule (a set of weights and abscissae) over the unit interval must be input. A quadrature rule over the standard triangle is required in H3LC and a quadrature rule over the unit square (consisting of standard quadrature rules applied in both the generator and angular directions) must be input to H3ALC. The input abscissae are mapped onto the element over which the integration is carried out. In general, the greater the number of quadrature points and the more efficient the choice of quadrature rule, the more accurate the resulting discrete integral operators will be. However, the computational cost of each subroutine is roughly proportional to the number of quadrature points.

In a practical computation, one of the subroutines will be called several times. The subroutines allow different quadrature rules to be input for each call. For example when p lies on the element a special quadrature rule is desirable, as discussed above. Also if one element is relatively large or relatively close to the point p then a relatively more accurate quadrature rule is desirable.

In subroutines H2LC and H3LC the parameter MAXNQ sets the limit on the number of quadrature points in the quadrature rule defined in the generator and in the angular direction. The value of this parameter must not be changed between calls of the subroutines. The actual number of quadrature points, denoted NQ, must be less than or equal to MAXNQ but can be altered between calls. For the H3ALC subroutine the number of quadrature points consists of NGQ along the generator and NTQ in the angular direction. The values of NGQ and NTQ may be changed between calls of H3ALC, but they must always be less than MAXNGQ and MAXNTQ respectively. The total number of quadrature points used in H3ALC is NGQ× NTQ.

When the point p does not lie on the element, the integrand is continuously differentiable and hence standard quadrature rules are satisfactory. However, when p lies on the element, the integrands are continuous (after the subtractions described by equations (39) and (40)), but they are generally not differentiable at the point p. If the input quadrature rule does not take this property into account then it could lead to an inaccurate evaluation of the discrete forms.

In subroutine H2LC, the parameters AQ and WQ store the abscissae and the weights of the quadrature rule. The parameters AGQ, WGQ and ATQ, WTQ store similar quadrature rules for the generator and angular directions in H3ALC. Note that in H3ALC the element is defined as the result of rotating the generator QA- QB about the z axis. In practice, the size of the element in the angular direction will vary greatly between calls of H3ALC. Hence, for an efficient overall method, the quadrature rule applied in the angular direction should depend on the radius of the conical element.

In subroutines H2LC and H3ALC, generally the most efficient quadrature rules to use are standard Gauss-Legendre rules. The weights and abscissae for such rules are listed in Stroud and Secrest (1966) and can also be generated using NAG routine D01BBF. When the point p lies on the element then an efficient quadrature rule is obtainable by dividing the domain at p and using suitably scaled Gauss-Legendre rules on the separate regions. Note that such division is only necessary in the generator direction in H3ALC.

In subroutine H3LC, the parameters XQ, YQ, WQ store the x coordinates, the y coordinates and the weights of the quadrature rule over the standard triangle. Quadrature rules of the Gauss-Legendre type are listed in Laursen and Gellert (1978). In the case when p lies on the element then a suitable quadrature rule can be obtained by dividing the triangle into three regions by connecting the point p to the vertices of the element and map suitably scaled standard quadrature rules onto the three separate triangular regions.

5.4  Validation of the input

If LVALID is .TRUE. then a file for storing and error messages or warnings must be open prior to calling the subroutine. For subroutine H2LC the file H2LC.ERR must be opened using a Fortran statement like OPEN(UNIT=10,FILE='H2LC.ERR',STATUS='UNKNOWN'). Similarly, for subroutines H3LC and H3ALC the files H3LC.ERR and H3ALC.ERR respectively should be opened.

If LVALID is .TRUE. then the parameters to the subroutines EGEOM and EQRULE set the maximum absolute tolerance allowed in the input geometrical and quadrature rule data respectively. For example if the corresponding components of QA and QB are all within EGEOM of each other (the points coincide) or the sum of the weights in the quadrature rules is greater than NQ× EQRULE (allowing for an error of EQRULE in each component) then an error message will be output to the error file. In such cases the subroutine will be aborted and the parameter LFAIL will register .TRUE..

6  Test problem - the Burton and Miller Equation

A direct or indirect integral equation formulation of the Helmholtz equation will usually consist of some or all of the Helmholtz integral operators introduced in section 2. The solution of the Helmholtz equation (1) in the exterior E to a closed boundary S with the appropriate radiation condition can be reformulated as the Burton and Miller integral equation [7],

The solution of the Helmholtz problem consists of two stages. The primary stage involves the solution of the integral equation (41). This yields the solution f(p) for points on S. The secondary stage involves computing the solution f(p) for points in E using (42). Methods of this type, known as the Boundary Element Method (see, for example Banerjee and Butterfield [4]), for solving the exterior Helmholtz equation have been considered in references [32], [8], [26], [9], [29], [36], [34], [31], [10], [1], [28], [14], [2], [17], [18], [19], [22] [24].

Applying collocation to (41) generally requires that the closed boundary S is replaced by an approximate boundary S made up of a set of n elements DS_j (j = 1,..,n) in the way described in section 3. Let the points p_j (j = 1,..,n) with p_j Î DS_j be the collocation points. In this paper we consider only the approximation of the boundary functions by a constant on each elment. It is helpful to adopt the following notation:

The application of collocation to the Burton and Miller equation give the following linear systems of approximations

The secondary stage of the boundary element method requires the calculation of the approximation to f(p) where p is a point in the approximate exterior domain E. For this the discrete forms are substituted into (42) to give

Test problems can easily be devised by setting f(p) = G_k(p, p^*) where p Î E ÈS, p^* Î D (D is the interior to S) with ¶f/ ¶v_p (p) = ¶G_k/ ¶v_p (p,p^*) for p Î S.

7  Subroutine H2LC

In this section the Fortran subroutine H2LC is described. The subroutine computes the discrete form of the two-dimensional Helmholtz integral operators. Details of the methods employed in the subroutine are given. The subroutine is applied to the Burton and Miller integral equation for the problem where the boundary is a circle and the results are given.

7.1  Background Methods

The regular integrals that arise are approximated by a standard quadrature rule such as a Gauss-Legendre rule which is specified in the parameter list to the subroutines. Tables of Gauss-Legendre rules are given in Stroud and Secrest (1966) and can also generated from the NAG library [30]. The non-regular integrals that arise in the formulae (39) and (40) are computed via the following methods. See Jaswon and Symm [11] and Kirkup [14] for the background to these methods.

The M₀ and M₀^t operators have regular kernels, hence the aim is to find expressions for:

7.2  Test Problem and Results

The test problem is in the file H2LC.FOR . The test problem consists of a boundary of a circle of unit radius, centred at the origin. The exterior solution and Neumann boundary condition are determined by a point source at (0,0.5) and a point sink at (0,-0.5). The circle is approximated by 32 straight line elements of equal length, the vertices of the regular polygon lying on the original circle. Eight point Gauss-Legendre quadrature rules are used for all but the computation of the diagonal components of the matrices where 2 ×8 point composite Gaussian-Legendre quadrature rules are used.

The functions FNSQRT and FNLOG simply call the corresponding standard Fortran functions. The subroutine FNHANK was constructed through the use of NAG routines S17AEF, S17AFF, S17ACF and S17ADF and the subroutine is only suitable for real values of k. Note that routines for computing the Hankel functions for complex arguments are now available in Mark 14 of the NAG Fortran library [30]. The subroutine FNHANK is listed in this section.

C SUBROUTINE TO COMPUTE HANKEL FUNCTIONS
C FOR THE MOMENT ONLY Z REAL IS CATERED FOR

      SUBROUTINE FNHANK(Z,H)
      COMPLEX*16 Z,H(0:1)
      COMPLEX*16 J(0:1),Y(0:1)
      COMPLEX*16 CIMAG
      REAL*8 ZERO,ONE
      REAL*8 X
      REAL*8 S17AEF,S17AFF,S17ACF,S17ADF
      INTEGER IFAIL
      IF (ABS(AIMAG(Z)).GT.1.0E-6) STOP
      IFAIL=0
      ZERO=0.0D0
      ONE=1.0D0
      CIMAG=CMPLX(ZERO,ONE)
      X=REAL(Z)
      J(0)=S17AEF(X,IFAIL)
      J(1)=S17AFF(X,IFAIL)
      Y(0)=S17ACF(X,IFAIL)
      Y(1)=S17ADF(X,IFAIL)
      H(0)=J(0)+CIMAG*Y(0)
      H(1)=J(1)+CIMAG*Y(1)
      END

Table 1A: Solution for the circle with k = 0.0, a = 1, b = 0
Point	Exact Solution	Computed Solution
1	0.173161	0.173408
2	0.160666	0.160913
3	0.139776	0.140010
4	0.114977	0.115174
5	0.089028	0.089176
6	0.063136	0.063236
7	0.037647	0.037704
8	0.012506	0.012524
(0,2)	0.081300	0.081008

Table 1B: Solution for the circle with k = 1.0, a = 1, b = i
Point	Exact Solution	Computed Solution
1	0.204820 + i0.106145	0.210022 + i0.106022
2	0.190535 + i0.102052	0.195224 + i0.101934
3	0.166397 + i0.094027	0.170345 + i0.093919
4	0.137368 + i0.082387	0.140541 + i0.082294
5	0.106643 + i0.067589	0.109068 + i0.067513
6	0.075745 + i0.050206	0.077453 + i0.050150
7	0.045200 + i0.030909	0.046214 + i0.030875
8	0.015018 + i0.010435	0.015355 + i0.010424
(0,2)	0.028905 + i0.140053	0.028868 + i0.141126

8  Subroutine H3LC

In this section the Fortran subroutine H3LC is described. The subroutine computes the discrete form of the three-dimensional Helmholtz integral operators. Details of the background methods employed by the subroutines are given. The subroutine is applied to the Burton and Miller equation for the test problem where the boundary is a cube and results are given.

8.1  Background Methods

The regular integrals that arise are approximated by a quadrature rule defined on a triangle. Laursen and Gellert (1978) contains a selection of Gauss-Legendre quadrature rules for the triangle. The non-regular integrals that arise in the formulae (41) and (42) are computed by the following methods. See Jaswon and Symm [11], Terai [34], Banerjee and Butterfield [4] and Kirkup [14] for the background to these methods.

The M₀ and M₀^t operators have regular kernels, hence the aim is to find expressions for:

The integrals (54) and (55) may be written in the form:

8.2  Test Program and Results

The test problem is in the file H3LC.FOR . The test problem consists of a boundary of a cube with vertices (1,1,1), (1,1,-1), (1,-1,1), (1,-1,-1), (-1,1,1), (-1,1,-1), (-1,-1,1) and (-1,-1,-1). The exterior solution and the Neumann boundary condition are determined by a point source at the centre of the cube. The cube is divided into 24 uniform triangular elements. Seven point Gaussian quadrature rules are applied to all but the computation of the diagonal components of the matrices, where a 3×7 point rule is used. Solutions are sought on the boundary and at the point (0,0,2) at wavenumbers k = 0, 1, i, 1+i. The values of a and b in the Burton and Miller equation (41) are set as a = 1, b = 1 for all four values of k.

The functions FNSQRT and FNEXP simply call the corresponding standard Fortran functions. Since the solution is similar at each collocation point, a comparison between one exact and one computed solution is given at one surface point and the exterior point. The results are given in table 2.

Table 2: Solution for the cube.
k	a and b	Point	Exact Solution	Computed Solution
0.0	a = 1, b = i	Coll. Pt.	0.832050	0.829128 + i0.119394
0.0	a = 1, b = i	(0,0,2)	0.500000	0.535075 - i0.000003
1.0	a = 1, b = i	Coll. Pt.	0.300064 + i0.776060	0.163472 + i0.840514
1.0	a = 1, b = i	(0,0,2)	-0.208073 + i0.454649	-0.242757 + i0.516164
1.0i	a = 1, b = i	Coll. Pt.	0.250145	0.252107 + i0.048411
1.0i	a = 1, b = i	(0,0,2)	0.067668	0.068024 + i0.002300
1.0 + 1.0i	a = 1, b = i	Coll. Pt.	0.090211 + i0.233313	0.049252 + i0.246489
1.0 + 1.0i	a = i, b = i	(0,0,2)	-0.002816 + i0.061530	-0.003732 + i0.061234

9  Subroutine H3ALC

In this section, the Fortran subroutine H3ALC is described. The subroutine computes the discrete form of the axisymmetric three-dimensional Helmholtz integral operators. Details of the methods employed by the subroutine are given. The subroutine is applied to Burton and Miller equation for the test problem where the boundary is a sphere with a Neumann boundary condition and results are given.

9.1  Background Methods

The regular integrals that arise are approximated by a two-dimensional quadrature rule defined on a rectangle which is specified in the parameter list to the subroutine. These integrals can be approximated using a Gauss-Legendre rule in the generator and q directions. The non-regular integrals that arise in the formula are computed by the following methods.

The M₀ and M₀^t operators have regular kernels, hence the aim is to find expressions for the following:

The integral in (56) is evaluated through dividing the integral with respect to the generator direction into two parts at p and transforming the integral through changing the power of the variable in line with a method described in references [3] and [14]. The resulting regular integral on both parts is computed via the quadrature rule supplied to the routine.

The integral in (57) is evaluated by using the result that if the surface of integration in (57) is extended to enclose a three-dimensional volume then the integral vanishes (see [14]). As each element is a truncated right circular cone, a 45^o right circular cone is added to each flat side of the element. The integrals over the two 45^o cones are regular and are computed by a composite rule based on the quadrature rule supplied to the subroutine. The solution is thus equal to minus the sum of the integrals over the two 45^o cones.

9.2  Subroutine Introduction

9.3  Test Problem and Results

The test problem is in the file H3ALC.FOR . The test problem consists of a boundary of a sphere of unit radius, centred at the origin. The exterior solution and Neumann boundary condition are determined by a point source at (0,0.5) and a point sink at (0,-0.5) in (R,z) coordinates.

The sphere is divided into 16 conical elements with the vertices of each element forming an equal angle about the origin. Eight point Gauss-Legendre rules are applied in the generator direction to all but the computation of the diagonal components of the matrices where a 2×8 point Gauss-Legendre rule is used. A composite 8 point Gauss-Legendre rule is applied in the q direction so that the density of points is approximately equal to the density of points in the generator direction.

Solutions are sought on the boundary and at the point (0,0,2) at wavenumbers k = 0, 1, i, 1+i. The values of a and b in the Burton and Miller equation (41) are set as a = 1, b = 0, a = 1, b = 1, a = 1, b = 0, a = 1, b = 0 for each respective value of k.

The functions FNSQRT and FNEXP simply call the corresponding standard Fortran functions. Since the solution is antisymmetric about the R axis, a comparison of exact and computed solution is given at 8 collocation points. The computed and exact solutions are also compared at the point (0,0,2).

Table 3A: Solution for the sphere with k = 0.0, a = 1, b = 0
Point	Exact Solution	Computed Solution
1	1.313632	1.319734
2	1.174095	1.178789
3	0.963395	0.967408
4	0.744850	0.748010
5	0.546051	0.548242
6	0.371109	0.372513
7	0.215024	0.215787
8	0.070406	0.070644
(0,0,2)	0.266667	0.265787

Table 3B: Solution for the sphere with k = 1.0, a = 1, b = i
Point	Exact Solution	Computed Solution
1	1.685653 + i0.292436	1.770379 + i0.296744
2	1.525811 + i0.281171	1.597400 + i0.285144
3	1.278466 + i0.259084	1.334422 + i0.262705
4	1.012108 + i0.227037	1.053972 + i0.230196
5	0.758596 + i0.186279	0.788669 + i0.188865
6	0.524952 + i0.138386	0.545144 + i0.140305
7	0.308010 + i0.085203	0.319625 + i0.086384
8	0.101502 + i0.028767	0.105292 + i0.029165
(0,0,2)	0.367616 + i0.425608	0.371234 + i0.431319

Table 3C: Solution for the sphere with k = 1.0i, a = 1, b = 0
Point	Exact Solution	Computed Solution
1	1.046648	1.051918
2	0.922641	0.926596
3	0.739520	0.742876
4	0.556228	0.558839
5	0.396960	0.398729
6	0.263713	0.264822
7	0.150322	0.150913
8	0.048804	0.048985
(0,0,2)	0.115919	0.115353

Table 3D: Solution for the sphere with k = 1.0 + 1.0i, a = 1, b = 0
Point	Exact Solution	Computed Solution
1	1.035865 + i0.429558	1.041259 + i0.431229
2	0.908338 + i0.402124	0.912356 + i0.403586
3	0.720630 + i0.354254	0.724026 + i0.355550
4	0.534374 + i0.294661	0.537001 + i0.295743
5	0.375229 + i0.229906	0.376987 + i0.230730
6	0.245410 + i0.163746	0.246498 + i0.164320
7	0.138172 + i0.097833	0.138743 + i0.098168
8	0.044555 + i0.032521	0.044728 + i0.032630
(0,0,2)	0.036827 + i0.128731	0.036316 + i0.128244

10  Concluding Discussion

Integral equation methods such as the boundary element method are becoming increasingly popular as methods for the numerical solution of linear elliptic partial differential equations such as the Helmholtz equation. The application of (discrete) collocation to the integral equation formulation of the Helmholtz equation requires the computation of the discrete operators. In this paper Fortran subroutines for the evaluation of the discrete Helmholtz integral operators resulting from the use of constant elements and the most simple boundary approximation to two-dimensional, three-dimensional and axisymmetric problems have been described and demonstrated.

The subroutines have been designed to be easy-to-use, flexible, reliable and efficient. It is the intention that the subroutines are to be used as a `black box' which can be utilised either for further analysis of integral equation methods or in software for the solution of practical physical problems which are governed by the Helmholtz or Laplace equations. The subroutines have already been applied to various problems in references [15], [19], [22], [24]. A full description of the application of the subroutines to interior, exterior and modal analysis problems is given in Kirkup [25].

References

[1]

S. Amini and D. T. Wilton (1986). An Investigation of the Boundary Element Method for the Exterior Acoustic Problem, Computer Methods in Applied Mechanics and Engineering, 54, 49-65.

[2]

S. Amini, P. J. Harris and D. T. Wilton (1992). Coupled Boundary and Finite Element Methods for the Solution of the Dynamic Fluid-Structure Interaction Problem, Lecture Notes in Engineering, C. A. Brebbia and S. A. Orszag (eds), 77, Springer-Verlag, (1992).

[3]

K. E. Atkinson (1989). An Introduction to Numerical Analysis, John Wiley and Sons.

[4]

P. K. Banerjee and R. Butterfield (1981). Boundary Element Methods in Engineering Science, McGraw-Hill.

[5]

J. Ben Mariem and M. A. Hamdi (1987). A new Boundary Element Method for Fluid-Structure Interaction Problems, International Journal of Numerical Methods in Engineering, 24, 1251-1267.

[6]

R. J. Bernhard, B. K. Gardner, C. G. Mollo and C. R. Kipp (1987). Prediction of Sound Fields in Cavities Using Boundary-Element Methods, AIAA Journal, 25, 1176-1183.

[7]

A. J. Burton and G. F. Miller (1971). The Application of Integral Equation Methods to the Numerical Solution of some Exterior Boundary Value Problems, Proc. Royal Society, London, A323, 201-210.

[8]

A. J. Burton (1973). The Solution of Helmholtz Equation in Exterior Domains using Integral Equations. NPL Report NAC30, National Physical Laboratory, Teddington, Middlesex, UK.

[9]

A. J. Burton (1976). Numerical Solution of Acoustic Radiation Problems, NPL Report OC5/535, National Physical Laboratory, Teddington, Middlesex, UK.

[10]

D. Colton and R. Kress (1983), Integral Equation Methods in Scattering Theory, Wiley/Interscience, New York.

[11]

M. A. Jaswon and G. T. Symm (1977). Integral Equation Methods in Potential Theory and Elastostatics, Academic Press.

[12]

R. A. Jeans and I. C. Mathews (1990). Solution of Fluid-Structure Interaction Problems using a Coupled Finite Element and Variational Boundary Element Technique, Journal of the Acoustical Society of America, 88(5), 2459-2466.

[13]

C. R. Kipp and R. J. Bernhard (1987). Prediction of Acoustical Behavior in Cavities Using an Indirect Boundary Element Method, ASME Journal of Vibration and Acoustics, 109, 22-28.

[14]

S. M. Kirkup (1989). Solution of Exterior Acoustic Problems by the Boundary Element Method, PhD thesis, Brighton Polytechnic, Brighton, UK.

[15]

S. M. Kirkup and S. Amini (1991). Modal Analysis of Acoustically-loaded Structures via Integral Equation Methods, Computers and Structures, 40(5), 1279-1285.

[16]

S. M. Kirkup (1991). The Computational Modelling of Acoustic Shields by the Boundary and Shell Element Method, Computers and Structures, 40(5), 1177-1183.

[17]

S. M. Kirkup and D. J. Henwood (1992). Computational Solution of Acoustic Radiation Problems by Kussmaul's Boundary Element Method, Journal of Sound and Vibration, 152(2), 388-402.

[18]

S. M. Kirkup and D. J. Henwood (1992). Methods for Speeding-Up the Boundary Element Solution of Acoustic Radiation Problems, Transaction of the ASME Journal of Vibration and Acoustics, 114, 374-380.

[19]

S. M. Kirkup and S. Amini (1993). Solution of the Helmholtz Eigenvalue Problem via the Boundary Element Method, International Journal for Numerical Methods in Engineering, 36(2), 321-330.

[20]

S. M. Kirkup and D. J. Henwood (1994). An Empirical Error Analysis of the Boundary Element Method with application to Laplace's Equation, Applied Mathematical Modelling, 18, 32-38.

[21]

S. M. Kirkup (1994) Computational Solution of the Acoustic Field Surrounding a Baffled Panel by the Rayleigh Integral Method, Applied Mathematical Modelling, 18, 403-407. See also RIM .

[22]

S. M. Kirkup and M. A. Jones (1996). Computational methods for the Acoustic Modal Analysis of an Enclosed Fluid with application to a Loudspeaker Cabinet, Applied Acoustics, 48(4), 295-299 (1996).

[23]

S. M. Kirkup (1996). Fortran Codes For Computing the Helmholtz Integral Operators: User Guide. Report MCS-96-06, Department of Computing and Mathematical Sciences, University of Salford, Salfrod M5 4WT, England.

[24]

S. M. Kirkup (1997). Solution of Discontinuous Interior Helmholtz Problems by the Boundary and Shell Element Method. Computer Methods in Applied Mechanics and Engineering, 48(4), 275-299, (1997).

[25]

S. M. Kirkup (1998). The Boundary Element Method in Acoustics, Integrated Sound Software.

[26]

R. E. Kleinmann and G. F. Roach (1974). Boundary Integral Equations for the three-dimensional Helmholtz Equation, SIAM Review, 16 (2), 79-90.

[27]

M. E. Laursen and M. Gellert (1978). Some Criteria for Numerically Integrated Matrices and Quadrature Formulas for Triangles, International Journal for Numerical Methods in Fluids, 12, 67-76.

[28]

I. C. Mathews (1986). Numerical Techniques for three-dimensional steady state Fluid-Structure Interaction, Journal of the Acoustical Society of America, 79(5), 1317-1325.

[29]

W. L. Meyer, W. A. Bell, B. T. Zinn and M. P. Stallybrass (1978). Boundary Integral Solutions of three dimensional Acoustic Radiation Problems, Journal of Sound and Vibration, 59(2), 245-262.

[30]

NAG Library, The Numerical Algorithms Group, Oxford, UK.

[31]

M. N. Sayhi, Y. Ousset and G. Verchery (1981). Solution of Radiation Problems by Collocation of Integral Formulations in terms of Single and Double Layer Potentials, Journal of Sound and Vibration, 74(2), 187-204.

[32]

H. A. Schenck (1968). Improved Integral Formulation for Acoustic Radiation Problems, Journal of the Acoustical Society of America, 44(1), 41-68.

[33]

A. H. Stroud and D. Secrest (1966). Gaussian Quadrature Formulas, Prentice Hall.

[34]

T. Terai (1980). On the Calculation of Sound Fields around three-dimensional Objects by Integral Equation Methods, Journal of Sound and Vibration, 69(1), 71-100.

[35]

A. G. P. Warham (1988). The Helmholtz Integral Equation for a Thin Shell, NPL Report DITC 129/88, National Physical Laboratory, Teddington, Middlesex.

[36]

D. T. Wilton (1978). Acoustic Radiation and Scattering from Elastic Structures, International Journal for Numerical Methods in Engineering, 13, 123-128.