1 Introduction
Solution of the Laplace equation for the potential problem underpins many applications in electrostatics and heat conduction. It is also central to modelling moving or deformable boundaries in the high Reynolds number regime in fluid mechanics. There, viscous and boundary layer effects are not dominant whereby a description based on potential flow can therefore be used as a first approximation or as a base case for further refinement. Numerous examples can be found in hydrofoil dynamics (Faltinsen and Semenov, 2008; Xu and Wu, 2013), the description of waves (Xue et al., 2001; Liu et al., 2001), cavitation or supercavitation phenomena (Blake et al., 1986, 1987), in civil, marine and ocean engineering and oscillating bubble dynamics in sonophysics and sonochemistry (Leighton, 1994).
Many of the above applications in multiphase fluid mechanics require the accurate tracking of moving interfaces that can be cumbersome and expensive to implement using grid based methods, especially in 3D. Consequently the use of the boundary element method (BEM) is advantageous because computational effort can focus on modelling all the important interfaces with the additional benefit of reducing the dimension of the problem by one, thus obviating the need to compute solutions in the whole flow domain. To track moving and deforming interfaces with precision, the Laplace equation is solved at each time step and the boundaries evolve according to the unsteady Bernoulli equation. This approach can readily be adapted to handle large or infinite domains or boundaries. The BEM is especially appealing for infinite fluid domains since the behaviour at infinity can be accounted for analytically. Although the BEM generates dense matrix equations, the Fast Multipole Boundary Element Method (Liu, 2009) can be used to reduce both CPU time and memory requirement from to . Thus in spite of being a simplification of the full NavierStokes description, the theory of potential flow together with the Bernoulli equation to describe unsteady problems occupies an important role in multiphase fluid dynamics.
The inherent use of the fundamental solution in the formulation of the BEM means that the integral equation contains singular kernels. This characteristic feature has been described as “a mathematical monster that leaps out of every page” due to “very unfamiliar and complex mathematics” (Becker, 1992). Since the physical problem itself is perfectly well behaved on the boundaries, such singularities are numerical inconveniences generated by the mathematical formulation (Liu and Rudolphi, 1999). Traditionally, the singular behaviour is dealt with by a local change of variables in the evaluation of the surface integrals (Telles, 1987) that comes at the expense of additional coding effort. Previous attempts to remove such singularities analytically required the introduction of additional unknowns such as derivatives tangential to the surface that have to be found by developing and solving extra integral equations (Liu and Rudolphi, 1999; Chen et al., 2005). Another method to remove the singularities requires finding additional parameters that have to be determined on a fictitious ‘nearby’ boundary (Cao et al., 1991; Zhang et al., 1999; Chen et al., 2009).
Moreover, when two different boundaries or two parts of the same boundary become close to each other, the traditional implementation of the BIM does not prevent the deleterious influence of singularities that originate from nodes of one boundary on the other. In particular, for moving boundary potential problems, it is highly desirable to eliminate all singular terms that arise in traditional formulations of the BEM as this avoids the need to track the spatial separation of different parts of the boundary and to determine when remedial action may be required.
Here we show that the wellknown mathematical singularities that arise in the BEM for solving the Laplace equation can be removed analytically without generating additional unknowns or equations to be solved. The desingularised formulation is given for both general 3D and axisymmetric cases. The approach can also be applied to evaluate the potential at points near boundaries in a numerically robust way. The implementation and resulting improvement in accuracy are illustrated with a number of examples: a problem with osculating boundaries, a 3D problem with a semiinfinite domain that arises in the study of wave drag near a deformable surface (with movies in the electronic supplement) and problems involving domains with corners in 2D.
2 Nonsingular formulation of the boundary integral method
To develop a nonsingular formulation of the boundary integral method, consider the internal problem in a 3D domain that is bounded by the closed surface , as shown in Fig. 1. The potential, , is governed by the Laplace equation
(1) 
By using the 3D free space Green’s function and with the help of Green’s second identity, the solution of Eq. (1) can be found by solving the conventional boundary integral equation
(2) 
where is the solid angle at the observation point on and the surface element is at . The normal derivatives are defined by and , where
is the unit normal vector pointing out of the internal domain at
(Becker, 1992). If either or or a mixed boundary condition is given on the whole boundary , the corresponding or/and can be obtained from Eq. (2). Although numerical integration over the singularities in and can be effected by established methods (Telles, 1987), our aim is to remove such mathematical singularities analytically at the outset.Corresponding to a given point on the boundary, we construct a function that also satisfies the Laplace equation and hence Eq. (2), with the properties that as , and , where is the outward unit normal vector at . We choose to be of the form
(3) 
so the function must satisfy
(4) 
Taking the difference between the conventional boundary integral equations for and for we obtain an integral equation relating and on that replaces the conventional boundary integral equation in Eq. (2):
(5) 
The key point is that both integrands in Eq. (2) are now nonsingular and thus any convenient quadrature method can be used to evaluate the integrals. Note also that the solid angle no longer appears. Implicit in the derivation of Eq. (2) is that the outward unit normal is uniquely defined at . The implementation of Eq. (2) at nodes where the normal is not defined, e.g. at a corner, is considered in Section 5. A detailed proof of this desingularisation method using the linear function: has been given elsewhere Klaseboer et al. (2012) and establishes the theoretical basis of the numerical scheme constructed earlier to regularise the system of linear equations that arise from the standard implementation of the BEM (Klaseboer et al., 2009). This approach has also been extended to desingularise boundary integral equations that arise in Stokes flow, in solving the Helmholtz equation and the equations of linear elasticity using a linear function for (Klaseboer et al., 2012).
However, for potential problems involving domains of infinite extent (external problems), a different choice of is required since the linear function is unbounded at infinity. One possible choice for is
(6)  
The constant vector is the position of any convenient point that is located outside the solution domain and satisfies . In the next section, we will use this form of to formulate a nonsingular BEM for axisymmetric problems and in Section 6, we demonstrate how Eq. (6) can be used to formulate a nonsingular boundary integral problem in a semiinfinite domain to solve a 3D timedependent potential flow problem. Before proceeding, we note that we are not restricted to using a linear function or the form given by Eq. (6) for to construct nonsingular versions of the boundary integral equation. In fact, any form of that satisfies the conditions given by Eq. (4) can be considered for use in Eq. (2) to remove the singular behaviour due to the presence of and at .
3 Nonsingular axisymmetric boundary integral equation
For problems that possess axial symmetry whereby in the cylindrical variables: and , we have and , the integration over the azimuthal angle, , can be evaluated analytically (Becker, 1992). Since a point on the axisymmetric boundary surface is specified by some given equation , the surface integrals for the axisymmetric case can be reduced to 1D integrals.
The given by Eq. (6) will be an axisymmetric function if we choose (in Cartesian coordinates) to lie on the axis of symmetry with located outside the solution domain. We measure the azimuthal angle, , relative to so that and the surface normals are given by and . With these definitions, the axisymmetric function in Eq. (6) and its normal derivative are given explicitly by
(8) 
where , , and .
The axisymmetric version of the nonsingular boundary integral equation obtained after performing the integration in Eq. (2) using Eq. (6) is (see also (Becker, 1992))
(9) 
where . The arc length element may be expressed as using the equation: that defines the axisymmetric surface. The quantities and are the complete elliptic integrals of the first and second kind of the parameter, (see (Becker, 1992; Abramowitz and Stegun, 1972)), with
(10) 
The axisymmetric boundary integral equation given by Eq. (3) is nonsingular because the divergence in as in the limit is now suppressed by the terms containing the difference between and and between and that vanish as as . Thus the line integral in Eq. (3) can be evaluated using any quadrature method. This is a consequence of using an axisymmetric form of given by Eq. (6) that simplified to Eqs. (3) and (8).
An additional bonus of the nonsingular axisymmetric formulation in Eq. (3) is that no special effort is required to handle node points that are located on the axis of symmetry that has been a technical inconvenience of the conventional axisymmetric form of the boundary integral equation (Becker, 1992).
4 Robust method to calculate the potential near a boundary
The method of desingularising the boundary integral equation by subtracting out the singular behaviour can also be used to give a numerically robust method to calculate the value of the potential at points near boundaries. To find the potential, , at an observation point that is located inside the domain but may be close to the boundary, we begin with the conventional boundary equation, Eq. (2), for the function , with given by Eq. (3) and for points inside the domain,
(11) 
The point is located on the boundary and its relation to is specified below. The nearly singular behaviour of the integrands when is close to the boundary can be eliminated by subtracting Eq. (2) from Eq. (4) to give
(12) 
When the observation point is near the boundary, a suitable choice for the boundary point, , is to project onto the boundary using the relation , with .
Another important benefit arising from the numerical robustness of the present desingualarised formulation of the BIM occurs when different parts of the domain boundary become close together as often happens in multiscale moving boundary problems. A simple example of this is when the boundary comprises of two nearly touching spheres and the point is located near where the spherical boundaries are close together. The present formulation eliminates the near singular behaviour of the kernel at nodes that are on a different part of boundary that is close spatially to the observation point . A numerical demonstration of this is given in Section 6.
5 Corner nodes in 2D
The above formulation of the nonsingular form of the BEM for the potential problem assumes that the unit normal is uniquely defined at . When is a corner node on a 2D boundary, we choose to have the form
(13) 
where is the potential at , and the constants and are the values of normal derivatives as one approaches from the ‘left’ () and from the ‘right’ (). The functions and satisfy the following conditions
(14) 
(15) 
The value of at the corner can be expressed in two equivalent forms
(16) 
in which, and are the two tangential derivatives of the potential on side and at , respectively, and and are the unit tangential vectors along side and at , respectively. This compatibility condition in Eq. (16) provides an additional relation between and (Grilli and Svendsen, 1990) in the formulation of the BIM whereas the tangential derivatives can be constructed from the values of at neighbouring nodes.
A choice for and for a 2D corner problem can be constructed in terms of the surface normals and on either side of the corner at
(17) 
and
(18) 
6 Numerical demonstrations and examples
6.1 Two nearly touching spheres
To demonstrate the utility and robustness of our nonsingular boundary integral method for problems with near osculating boundaries, we consider the potential problem associated with two nearly touching spheres translating along their line of centre at identical constant speed. The spheres have radii and and are at a separation of at the point of closest approach (Fig. 2a). We compare the velocity potential obtained using the standard axisymmetric boundary integral method (Wang et al., 1996), with our nonsingular version given by Eq. (3) with set to the centre of each sphere. In Fig. 2b, we see that the error in the standard axisymmetric boundary integral method is very large in the region where the spheres are close together. The standard 3D boundary integral method gives errors similar to the standard axisymmetric method whereas the nonsingular 3D version given by Eq. (2) has the same accuracy as the axisymmetric version. This large error arises in the standard version of the boundary integral method because of the influence of the observation point on one sphere from the singular kernel centred at nearby nodes located on the other sphere. As we have shown in Section 4, with our nonsingular formulation, such effects do not arise.
6.2 Corner problem in 2D
To illustrate the implementation of our nonsingular formulation of the BIM, we consider three different interior potential problems in rectangular domains of different shapes, as shown in Fig. 3. The length of all four edges are set to be 1 unit. Uniform linear elements are employed on the boundary of the parallelogram specified by the angle, . At the corner nodes, the function given by Eq. (13) used in the nonsingular boundary integral equation (2) and the double node technique Grilli and Svendsen (1990) has been applied. The tangential derivatives of the potential in the compatibility condition in Eq. (16) are obtained by a fourth order finite difference scheme using values of along the boundaries. For the remaining nodes on the edges, the form of in Eq. (3) with is used in Eq. (2).
We tested three cases corresponding to analytical solutions:
(19) 
(20) 
(21) 
We solved the above corner problems with Dirichlet boundary conditions. The comparisons between the results for the normal derivatives obtained by our nonsingular BIM with 21 nodes on each edge and the analytical solutions are shown in Fig. 4 when and . The absolute errors for these cases can be found in Table 1.
Domain angle  Case I  Case II  Case III 

0.0046%  0.29%  0.87%  
0.012%  0.019%  0.77% 
Other pairs of and that satisfy Eqs. (14) and (15) can also be used, for example
(22) 
(23) 
where and are located outside of the calculation domain, and
(24) 
This pair of and are also used to solve the above three cases with Dirichlet boundary conditions when . Once again, 21 nodes are employed on each edge. The absolute errors between the results for the normal derivatives obtained by our nonsingular BIM and the analytical solutions for these cases are the same as Table 1.
This method of treating corner nodes can be extended to handle edges and vertices in 3D domains, but will not be considered here.
6.3 Wave drag at a semiinfinite deformable boundary
To illustrate the generality of our nonsingular formulation of the BEM, we examine the following wave drag problem that has sufficient complexity to be interesting (Landau and Lifshitz, 1966). Consider two identical spheres of radius, , separated by a constant distance between their centres and moving with constant speed at the same depth below an infinite deformable free surface in a gravity field, as shown in Fig. 5. The origin of the global reference frame is set to coincide with the initially undisturbed free surface (at ) with the surface elevation pointing upwards (Fig. 5). For simplicity, the upper phase (air) is assumed to have negligible mass density and interfacial tension effects have been omitted, although it is easy to dispense with such simplifications at the expense of introducing more physical parameters. This is a timedependent potential flow problem as surface waves will be generated on the deformable free surface while the pair of spheres travelling beneath it.
The potential flow velocity field, , generated by the moving spheres is found by solving the Laplace equation for the velocity potential, at each time step. On the spheres, the velocity is given and the potential is calculated. Initially, on the undisturbed free surface. Solving the mixed boundary value problem defined by the free surface and the two spheres, provides values of the normal velocity at the free surface that is then used to predict its shape at the next time step. The value of on the free surface at the next time step is found from the unsteady Bernoulli equation evaluated on the surface: where is the material derivative, the fluid density and the gravitational acceleration. Such timestepping then gives the spatiotemporal evolution of the surface waves. The wave drag force experienced by the two spheres that are moved at constant velocity in this 3D problem is calculated by integrating the pressure on the sphere surfaces.
The deformable interface extends to infinity in the direction where it asymptotes to a flat surface. Therefore far from the spheres, the potential and its derivative vanish asymptotically like and as . We therefore assume that beyond a radius of 40 from the spheres, we can take the interface to be flat and the potential and its derivative also take on their asymptotic forms so that we can evaluate analytically the contribution to the surface integrals from the far field, see for example (Gao and Davies, 1998; Ribeiro and Paiva, 2009). Since Eq. (2) is nonsingular, there is no need to map the far field elements to two triangular elements or introduce any artificial points at infinity as was done in Ribeiro and Paiva (2009). Finally, the integral over the half spherical surface at infinity will contribute two terms and to the lefthand side of Eq. (2).
On each sphere, 1280 linear triangular elements with 642 nodes are employed. On the free surface, 15000 elements with 7651 nodes are used. When is on a sphere, in Eq. (6) is taken to be the centre of the sphere. When is on the free surface, we choose the Cartesian components of to be and any convenient and that ensure .
With a constant time step of , the position of each surface node is updated with and the spheres are translated using . The calculation continues until the spheres have travelled a distance of . The elastic mesh technique (Wang et al., 2003) is applied on the free surface to ensure a uniform mesh even after many timesteps.
The drag force acting on each sphere is calculated by where is the surface of the sphere. In Fig. 6a we present the time variation of the force on each of the two spheres travelling with constant velocity . The three nondimensional parameters that govern this problem are: , and the Froude number, (Fig. 5). For our illustrative example, we chose a sphere separation of , a submerged depth of , and . Also shown are results for that corresponds to the absence of the deformable interface and the force on an isolated single sphere that corresponds to . The leading sphere experiences a retarding force in the direction whereas the trailing sphere experiences a force in the direction of travel. At , these two forces are equal and opposite as expected from the d’Alembert Paradox of potential flow. The proximity of the moving spheres to the deformable interface provides a net nonzero wave drag on the pair of moving spheres (Landau and Lifshitz, 1966). We point out that the treatment of the initial condition of this problem has been simplified by omitting the detailed effects of acceleration from rest that will give rise to an added mass term. The inclusion of this effect complicates the equation of motion, but does not affect the demonstration of the utility of the present nonsingular formulation of the BIM.
In Fig. 6b we show the variation of the drag force on the trailing sphere as a function of the sphere separation, at . There is excellent agreement between our nonsingular boundary integral method and the analytic result of Miloh (1977).
Snapshots of the spatiotemporal variations of the surface waves generated by the pair of moving submerged spheres are shown in Fig. 7. Variations of the surface wave amplitude along the direction of travel at two different times are shown in Fig. 7 (a) and (b). Three dimensional representations of the corresponding surface waves in Fig. 7 (c) and (d) show the lateral extent of surface disturbance that obviously depends on the sphere spacing and the depth of immersion. Movies of the wave amplitude and the surface wave corresponding to these figures are available as online supplementary material. Although potential flow is conservative with no energy dissipation or damping mechanism, the drag waves on the free surface appear damped because of kinetic energy of the moving sphere being distributed into the infinite fluid domain (Landau and Lifshitz, 1966).
7 Conclusions
In this paper, we have provided a fully nonsingular formulation of the boundary integral method for potential problems. The usual singularities in the kernels that arise from the fundamental solutions are removed analytically without introducing additional unknowns or extra equations to be solved. Apart from excising the “mathematical monsters” in the conventional boundary integral formulation (Becker, 1992), the amount of computer code needed to implement the nonsingular boundary integral equation has been reduced by about 60%. For the special case of axisymmetric problems, our formulation also removed the technical inconvenience associated with nodes on the axis of symmetry (Becker, 1992). The robustness and generality of our approach has been demonstrated with examples involving osculating boundaries, domains with corners and a wave drag due to an infinite deformable boundary.
Acknowledgement
DYCC is a Visiting Scientist at the Institute of High Performance Computing and an Adjunct Professor at the National University of Singapore. This work is supported in part by the Australian Research Council Discovery Project Grant Scheme.
References
 Faltinsen and Semenov (2008) O. M. Faltinsen, Y. A. Semenov, The effect of gravity and cavitation on a hydrofoil near the free surface, Journal of Fluid Mechanics 597 (2008) 371–394.
 Xu and Wu (2013) G. D. Xu, G. X. Wu, Boundary element simulation of inviscid flow around an oscillatory foil with vortex sheet, Engineering Analysis with Boundary Elements 37 (2013) 825–835.
 Xue et al. (2001) M. Xue, H. Xü, Y. Liu, D. K. P. Yue, Computations of fully nonlinear threedimensional wavewave and wavebody interactions. Part 1. Dynamics of steep threedimensional waves, Journal of Fluid Mechanics 438 (2001) 11–39.
 Liu et al. (2001) Y. Liu, M. Xue, D. K. P. Yue, Computations of fully nonlinear threedimensional wavewave and wavebody interactions. Part 2. Nonlinear waves and forces on a body, Journal of Fluid Mechanics 438 (2001) 41–66.
 Blake et al. (1986) J. R. Blake, B. B. Taib, G. Doherty, Transient cavities near boundaries. Part 1. Rigid boundary, Journal of Fluid Mechanics 170 (1986) 479–497.
 Blake et al. (1987) J. R. Blake, B. B. Taib, G. Doherty, Transient cavities near boundaries. Part 2. Free surface, Journal of Fluid Mechanics 181 (1987) 197–212.
 Leighton (1994) T. G. Leighton, The Acoustic Bubble, Academic Press, 1994.
 Liu (2009) Y. J. Liu, Fast Multiple Boundary Element Method: Theory and Applications in Engineering, Cambridge University Press, 2009.
 Becker (1992) A. A. Becker, The Boundary Element Method in Engineering: A complete Course, McGrawHill International (UK) Limited, 1992.
 Liu and Rudolphi (1999) Y. J. Liu, T. J. Rudolphi, New identities for fundamental solutions and their applications to nonsingular boundary element formulations, Computational Mechanics 24 (1999) 286–292.
 Telles (1987) J. C. F. Telles, A selfadaptative coordinate transformation for efficient numerical evaluation of general boundary element integrals, International Journal for Numerical Methods in Engineering 24 (1987) 959–973.
 Chen et al. (2005) H. B. Chen, J. F. Jin, P. Q. Zhang, P. Lü, Multivariable nonsingular BEM for 2D potential problems, Tsinghua Science and Technology 10 (2005) 43–50.
 Cao et al. (1991) Y. Cao, W. W. Schultz, R. F. Beck, Threedimensional desingularized boundary integral methods for potential problems, International Journal for Numerical Methods in Fluids 12 (1991) 785–803.
 Zhang et al. (1999) Y. L. Zhang, K. S. Yeo, B. C. Khoo, W. K. Chong, Simulation of threedimensional bubbles using desingularised boundary integral method, International Journal for Numerical Methods in Fluids 31 (1999) 1311–1320.

Chen et al. (2009)
W. Chen, Z. Fu, X. Wei,
Potential problems by singular boundary method satisfying moment condition,
CMES: Computer Modeling in Engineering & Sciences 54 (2009) 65–85.  Klaseboer et al. (2012) E. Klaseboer, Q. Sun, D. Y. C. Chan, Nonsingular boundary integral methods for fluid mechanics applications, Journal of Fluid Mechanics 696 (2012) 468–478.
 Klaseboer et al. (2009) E. Klaseboer, C. RosalesFernandez, B. C. Khoo, A note on true desingularization of boundary element methods for threedimensional potential problems, Engineering Analysis with Boundary Elements 33 (2009) 796–801.
 Abramowitz and Stegun (1972) M. Abramowitz, I. A. Stegun (Eds.), Handbook of Mathematical Functions With Formulas, Graphs, and Mathematical Tables, number 55 in Applied Mathematics Series, 10th ed., National Bureau of Standards, 1972.
 Grilli and Svendsen (1990) S. Grilli, I. Svendsen, Corner problems and global accuracy in the boundary element solution of nonlinear wave flows, Engineering Analysis with Boundary Elements 7 (1990) 178–195.
 Wang et al. (1996) Q. X. Wang, K. S. Yeo, B. C. Khoo, K. Y. Lam, Strong interaction between a buoyancy bubble and a free surface, Theoretical and Computational Fluid Dynamics 8 (1996) 73–88.
 Landau and Lifshitz (1966) L. D. Landau, E. M. Lifshitz, Volume 6 of Course of Theoretical Physics, Fluid Mechanics, 3rd ed., Pergamon Press, 1966.
 Gao and Davies (1998) X. W. Gao, T. G. Davies, 3D infinite boundary elements for halfspace problems, Engineering Analysis with Boundary Elements 21 (1998) 207–213.
 Ribeiro and Paiva (2009) D. B. Ribeiro, J. B. Paiva, A new infinite boundary element formulation applied to threedimensional domains, in: Proceedings of the World Congress on Engineering 2009, London, U.K., 2009, pp. 1468–1473.
 Miloh (1977) T. Miloh, Hydrodynamics of deformable contiguous spherical shapes in an incompressible inviscid fluid, Journal of Engineering Mathematics 11 (1977) 349–372.
 Wang et al. (2003) C. Wang, B. C. Khoo, K. S. Yeo, Elastic mesh technique for 3D BIM simulation with an application to underwater explosion bubble dynamics, Computers & Fluids 32 (2003) 1195–1212.
Comments
There are no comments yet.