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 M0 and M0t operators have regular kernels, hence the aim is to find expressions for the following:

{ L0 ~
e
 
}DG(p) =
ó
õ
DG 
 G0(p,qdSq  ,
(57)

{ N0 ~
e
 
}DG(p; vp) =
vp

ó
õ
DG 
  G0
nq
(p,qdSq   ,
(58)
where DG is a conical element, p Î DG (though not on an edge of the element).

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 45o right circular cone is added to each flat side of the element. The integrals over the two 45o 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 45o 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