A Pilot Fortran Software Library for the Solution of Laplace’s Equation by the Boundary Element Method

2.1 Boundary integral equation formulation of the interior Laplace problem

Laplace’s equation (2) governs the interior domain D enclosed by a boundary S. The solution must also satisfy a boundary condition, and it is important in terms of maintaining the generality of the method that this is in a general (Robin) form:

In the direct BEM, Laplace’s equation is replaced by an equivalent integral equation of the form:

The terminology ∂∗∂nqrepresents the partial derivative of the function* with respect to the unit outward normal at point q on the boundary. The function G is known as Green’s function. Physically, G(p, q) represents the effect observed at point p of a unit source at point q. For the Laplace equation, Green’s function is denoted by G and is defined as

for two-dimensional Laplace problems, where r=∣q−p∣.

Integral Eqs. (4) and (5) can be derived from the Laplace equation by applying Green’s second theorem. The power of the formulation lies in the fact that Eq. (4) relates the potential φ and its derivative on the boundary alone; no reference is made to φ at points in the domain in this particular boundary integral equation. In a typical boundary value problem, we may be given φ(q), ∂φq∂nqor a combination of such data on S. The boundary integral equation is a means of determining the unknown boundary function(s), followed by the domain solution from the given boundary data.

2.2 Operator notation

Operator notation is a useful shorthand in writing integral equations. Moreover, it will be shown that it is a very powerful notation in that it clearly demonstrates the connection between the integral equation and the linear system of equations that results from its discretisation.

Integral equations can always be written in terms of integral operators. For example, if ζ is a function defined on a (closed or open) boundary Г, then applying the following operation to ζ for all points p on Г

gives a function μ. This may be viewed as the application of an operator to the function ζ to return the function μ. More simply we may write

In Eq. (8)L represents the integral operator, and the subscript (Г) refers to the domain of integration. Г is used as a variable, representing either a whole boundary or a part of the boundary. The other three important Laplace integral operators are defined as follows:

where vpis any unit vector. In operator notation of the previous subsection, the integral equation formulation (3) can be written in the following form:

where vq=∂φq∂nq.

2.3 Direct boundary element method

For the direct boundary element method solution of the interior Laplace problem, that is, developed in this section, the initial stage involves solving boundary integral Eq. (4), returning (approximations to) both φ and ∂φ/∂n on S. The second stage of the BEM involved finding the solution at any chosen points in the domain D. The most straightforward method for solving integral equations like Eq. (4) is that of collocation. Collocation may be applied in a remarkably elementary form, which is termed C⁻¹ collocation in this text since it is derived by approximating the boundary functions by a constant on each panel. In this subsection the C⁻¹ collocation method is briefly outlined.

To begin with, the boundary S is assumed to be expressed as a set of panels:

Usually the panels have a characteristic form and cannot represent a given boundary exactly. For example, a two-dimensional boundary can be approximated by a set of straight lines. In order to complete the discretisation of the integral equations, the boundary functions also need to be approximated on each panel. In this work, it is the characteristics of the panel and the representation of the boundary function on the panel that together define the element in the boundary element method. By representing the boundary functions by a characteristic form on each panel, the boundary integral equations can be written as a linear system of equations of the form introduced earlier.

The term element refers not only to the form of ΔSjbut also to the method of representing the boundary functions on ΔSj. The C⁻¹ collocation method involves representing the boundary function by a constant on each panel:

The substitution of representations of this form for the boundary functions in the integral equation reduces it to discrete form. The simplifications allow us to rewrite Eq. (11) as the approximation:

where e is the unit function (e≡1). LeΔSj˜p, for example, for a specific point p, are the numerical values of definite integrals that together can be interpreted as the discrete form of the L integral operator.

The constant approximation is taken to be the value of the boundary functions at the representative central point (the collocation point) on each panel. By finding the discrete forms of the relevant integral operators for all the collocation points, a system of the form

for i = 1, 2 and n is obtained by putting p=pSiin the previous approximation. Note that because of the approximation of the boundary functions (and also the boundary approximation, if applicable), the discrete equivalent of Eq. (12) is an approximation relating the exact values of the boundary functions at the collocation points.

This system of approximations can now be written in the matrix-vector form:

with the matrix components defined by LSSij=LeΔSj˜pSi,MSSij=MeΔSj˜pSi. The vectors φ̂¯Sand v̂¯Sare representative or approximate values of φ and vat the collocation points. In the first stage of the boundary element, the system (18) is solved alongside the discrete form of the boundary condition (3):

The discrete forms are definite integrals that need to be computed usually by numerical integration. For the solution of Eqs. (18) and (19), the (approximation to) boundary data is known at the collocation points.

Once the (approximations to) functions on the boundary are known, after completing the initial stage of the boundary element method, the domain solution can be found. In the case of the interior Laplace problem, Eq. (13) will yield the domain solution. Similarly, the discrete equivalent of Eq. (11) may be derived:

for each point pDiin the domain D∼. Let the solution be sought at m domain points pDifor i=1,2,…m, and then the equation above, for all the domain points, is written as

where LDSij=LeΔSj˜pDi,MDSij=MeΔSj˜pDi.

2.4 LIBEM2 and L2LC modules, test problem and results

This section includes an outline of the Laplace interior BEM 2D (LIBEM2) and Laplace 2D linear boundary approximation constant element (L2LC) modules, and they are demonstrated by means of a test problem. The LIBEM2 module solves Laplace’s equation in an interior two-dimensional domain. L2LC is the most important component module.

L2LC. The L2LC module computes the discrete Laplace operators for two-dimensional problems. In the notation of this article, the routine computes LeΔSj˜p, MeΔSj˜p,MteΔSj˜pand NeΔSj˜p, where ΔSj∼is the panel that is the domain of integration and pis any point. The call to the subroutine has the following form:

SUBROUTINE L2LC(P,VECP,QA,QB,LPONEL, LVALID,EGEOM,LFAIL,

  * NEEDL,NEEDM,NEEDMT,NEEDN,L0,M0,M0T,N0),

where P is point p; VECP is a unit directional vector that passes through p; QA and QB are the points, either side of the panel and hence defining the panel; LPONEL is a logical switch that declares whether pis on the panel; and NEEDL is a logical switch that states whether the discrete Loperator is required and similar to the other operators. The computed values for the integrals are output in L0, M0, M0T and N0. For the straightforward direct BEM, developed in the previous section, only Land Moperators are required.

In general L2LC simply implements a Gaussian quadrature rule in order to determine the integral, using a higher-order rule when point pis close to the panel ΔSj∼. However, when point pis on the panel, then an exact integration is used [21, 22].

LIBEM2. The LIBEM2 module solves the interior Laplace problem and has the following form:

LIBEM2(MAXNODES,NNODE,NODES,MAXPANELS,NPANEL,PANELS,

  * MAXPOINTS,NPOINT,POINTS,

  * SALPHA,SBETA,SF,SINPHI,PINPHI,

  * LSOL,LVALID,TOLGEOM,

  * SPHI,SVEL,PPHI,

  * L_SS,M_SSPMHALFI,L_PS,M_PS,

  * PERM,XORY,C,workspace)

The boundary is set up through listing a set of nodal coordinates, and each panel is determined through the two nodal indices for the endpoints of the panel. The nodal coordinates are input through NODES and the panel information through PANELS. The nodes are oriented clockwise on each panel for an outer boundary and anticlockwise for any inner boundary. Usually, a solution in the domain is sought, and for this a set of (interior) domain points are set in POINTS. The boundary condition is set with the parameters SALPHA, SBETA and SF, setting αi, βiand fivalues in Eq. (19) for each panel.

Test problem. The test problem is that of solving Laplace’s equation on a unit square with the boundary conditions defined as shown in Figure 1. The solution is sought at the five interior points (0.25, 0.25), (0.75, 0.25), (0.25, 0.75), (0.75, 0.75) and (0.5, 0.5), and these are also illustrated in the figure.

The test problem is set up in the file LIBEM2_T.FOR. The boundary is defined by 32 nodes and panels. The nodes are indexed, starting with 1: (0.0, 0.0), 2: (0.0, 0.125), 3: (0.0, 0.25) and continue clockwise around the boundary until the final node 32: (0.125, 0.0). The panels are similarly set up in the clockwise sense with panel 1:1–2 (panel 1 links node 1 with node 2) and 2:2–3 until the final panel 32:32–1, linking the final node with the first node to complete the boundary. The boundary conditions shown in Figure 1 are then applied.

The exact solution is φp=10+10x; this is clearly a solution of Laplace’s equation and satisfies the boundary conditions. The exact solution at the interior points is therefore φ=12.5at the two points on the left, φ=17.5at the two points on the right and φ=15.0at the central point. The exact and computed results are shown in Table 1.

A set of nodal coordinates and each panel is determined through the two nodal indices for the endpoints of the panel.

The results from this test problem are also intuitively correct. With the left and right sides of the square at different potentials, it is common sense to expect the potential in the middle to be halfway between etc. The potentials can—most simply—be interpreted as temperatures in a steady-state heat conduction problem.

3.1 LBEMA and L3ALC modules, test problems and results

Let us start on the introduction of three-dimensional problems with the axisymmetric codes. These codes are used in a very similar manner. As with the two-dimensional problem, the component module L3ALC computes the integrals over the panels and is called as follows:

SUBROUTINE L3ALC(P,VECP,QA,QB,LPONEL,LVALID,EGEOM,LFAIL,

  * NEEDL,NEEDM,NEEDMT,NEEDN,DISL,DISM,DISMT,DISN).

For axisymmetric problems, the surface is defined by conical panels, which are defined by piecewise straight lines along the generator. The parameters follow a similar pattern as L2LC, except the points and vectors are in cylindrical rzcoordinates. QA and QB are the two points either side of the panel on the generator.

The LBEMA subroutine computes the solution of the Laplace equation by the direct boundary element method and has the following form:

LBEMA(MAXNODES,NNODE,NODES,MAXPANELS,NPANEL,PANELS,

  * LINTERIOR,MAXPOINTS,NPOINT,POINTS,

  * SALPHA,SBETA,SF,SINPHI,PINPHI,

  * LSOL,LVALID,TOLGEOM,

  * SPHI,SVEL,PPHI,

  * L_SS,M_SSPMHALFI,L_PS,M_PS,

  * PERM,XORY,C,workspace)

In LBEMA, NODES lists the rzcoordinates of the nodes on the generator of the surface, and PANELS states the two nodes that together define each panel. LINTERIOR is a logical input, which is set to TRUE if an interior problem is to be solved and FALSE for an interior problem.

The interior test problem is in file LBEMA_IT. The test problem is the unit sphere with the exact solution:

which is easily shown to be a solution of Laplace’s equation by writing r2=x2+y2.A Dirichlet boundary condition is applied; the solution is sought at four interior points, and the results for 18 elements are given in Table 2.

The exterior test problem is in file LBEMA_ET. The test problem is the unit sphere (approximated by 18 elements) with the exact solution:

where ris the distance from the origin. φis a solution of Laplace’s equation as it is a simple multiplication of Green’s function (22). A Dirichlet boundary condition is applied to the upper hemisphere, and a Neumann boundary condition is applied on the lower hemisphere. The solution is sought at four interior points, and the results are given in Table 3.

3.2 LBEM3 and L3LC modules, test problems and results

The LBEM3 and L3LC subroutines implement the boundary element method for general three-dimensional problems. As with the two-dimensional and axisymmetric codes, the component module L3LC computes the integrals over the panels. The L3LC subroutine is called as follows:

SUBROUTINE L3LC(P,VECP,QA,QB,QC,LPONEL,LVALID,EGEOM,LFAIL,

 * NEEDL,NEEDM,NEEDMT,NEEDN,DISL,DISM,DISMT,DISN)

The parameters follow a similar purpose as did in the L2LC, except that the points and vectors have three values. QA, QB and QC are the coordinates of the vertices of the triangular panel.

The LBEM3 module solves Laplace’s equation in a general interior or exterior three-dimensional domain and is called as follows:

LBEM3(MAXNODES,NNODE,NODES,MAXPANELS,NPANEL,PANELS,

 * LINTERIOR,MAXPOINTS,NPOINT,POINTS,

 * SALPHA,SBETA,SF,SINPHI,PINPHI,

 * LSOL,LVALID,TOLGEOM,

 * SPHI,SVEL,PPHI,

 * L_SS,M_SSPMHALFI,L_PS,M_PS,

 * PERM,XORY,C,WKSPC1,WKSPC2,WKSPC3)

As with LIBEM2 and LBEMA, NODES and PANELS define the boundary. However, in this case, NODES lists the three coordinates of each surface node, and PANELS lists the three nodal indices that make up each triangular panel.

The interior test problem is that of a unit sphere approximated by 36 triangular panels. The exact solution that is applied as a Dirichlet boundary condition is

The results at four interior points are given in Table 4.

The exterior test problem is that of a unit sphere approximated by 36 triangular panels, as in the previous test. The exact solution that is applied as a Dirichlet boundary condition is

where ris the distance from 000.5.The results at four exterior points are given in Table 5.

4.1 Integral equations and boundary element equations for thin shells

Following the work of Warham [27], the first step is to designate an ‘upper’ and ‘lower’ surface of a shell Γand denote them by ‘+’ and ‘−’. We then introduce the quantities of difference and average of the potential and its normal derivative across the surface:

The integral equation formulations for the Laplace equation in the exterior domain can now be written using the operator notation introduced earlier:

The boundary condition may be expressed in the following form:

The discrete equivalents of Eq. (21) are as follows:

4.2 LSEMA module, test problem and results

The LSEMA subroutine computes the solution of Laplace’s equation surrounding thin shells or discontinuities. As with the LBEMA, the subroutine relies on L3ALC to compute the matrix components in the systems (38)–(40). In this subsection, the LSEMA routine is demonstrated through solving a test problem.

The module LSEMA has the form:

LSEMA(MAXNODES,NNODE,NODES,MAXPANELS,NPANEL,PANELS,

 * MAXPOINTS,NPOINT,POINTS,

 * HA,HB,HF,HAA,HBB,HFF,

 * HIPHI,HIVEL,PINPHI,

 * LSOL,LVALID,TOLGEOM,

 * PHIDIF,PHIAV,VELDIF,VELAV,PPHI,

 * AMAT,BMAT,L_EH, M_EH,

 * PERM,XORY,C,WKSPC1,WKSPC2,WKSPC3).

The LSEMA parameters are similar to the LBEMA ones. However the expressions of the boundary condition and the boundary function are different.

HA stores the values of αon the shell panels, similarly HB, β; HAA, A; and HBB, B.The main output from the subroutine is PHIDIF that corresponds to δ̂¯Γ; PHIAV, Φ̂¯Γ; VELDIF, ν̂¯Γ; VELAV, V̂¯Γ; and PPHI, φ̂¯E.

The test problem is in file LSEMA_T. It consists of two circular coaxial parallel plates in the r,θplane, of radius 1.0 and a distance of 0.1 apart in the planes where z=0.0and z=0.1.A Dirichlet boundary condition is applied to both plates. On the plate at z=0.0, the potential of 0.0 is applied, and a potential (δ = 0, Φ=0) of 1.0 is applied on the other plate (δ = 0, Φ=1). A complete analytic solution is not available. However in the central region between the plates, a simple gradient of potential is intuitive, as discussed. The results from the test problem are listed in Table 6.

4.3 LSEM3 module, test problem and results

The LSEM3 module solves Laplace’s equation exterior to a thin shell in three dimensions. The subroutine call has the following form:

LSEM3(MAXNODES,NNODE,NODES,MAXPANELS,NPANEL,PANELS,

 * MAXPOINTS,NPOINT,POINTS,

 * HA,HB,HF,HAA,HBB,HFF,

 * HINPHI,HINVEL,PINPHI,

 * LSOL,LVALID,TOLGEOM,

 * PHIDIF,PHIAV,VELDIF,VELAV,PPHI,

 * AMAT,BMAT,L_EH,M_EH,

 * PERM,XORY,C,WKSPC1,WKSPC2,WKSPC3)

The definition of the important parameters can be found from the previous notes on LBEM3 and LSEMA. The test problem is in the file LSEM3_T, and it is similar to the test problem for LSEMA. This time, the open boundaries are two unit square plates of in x−yplanes. The two squares are 0.1 apart: one is at a potential of zero and the other is at a potential of one. The squares are each divided into 32 panels. The results at points between the squares, along a central axis, are shown in Table 7.