A robust and non-singular formulation of the boundary integral method for the potential problem

TL;DR

This paper introduces a non-singular boundary integral method for the Laplace equation that eliminates singularities analytically, improving accuracy and efficiency in multi-scale potential problems.

Contribution

It presents a new formulation of the boundary integral method that removes singularities, enabling standard quadrature and robust near-boundary potential calculations.

Findings

60% reduction in coding effort

Enhanced accuracy near boundaries

Applicable to 2D, 3D, and axisymmetric problems

Abstract

A non-singular formulation of the boundary integral method (BIM) is presented for the Laplace equation whereby the well-known singularities that arise from the fundamental solution are eliminated analytically. A key advantage of this approach is that numerical errors that arise due to the proximity of nodes located on osculating boundaries are suppressed. This is particularly relevant in multi-scale problems where high accuracy is required without undue increase in computational cost when the spacing between boundaries become much smaller than their characteristic dimensions. The elimination of the singularities means that standard quadrature can be used to evaluate the surface integrals and this results in about 60% savings in coding effort. The new formulation also affords a numerically robust way to calculate the potential close to the boundaries. Detailed implementations of this approach are illustrated with problems involving osculating boundaries, 2D domains with corners and a wave drag problem in a 3D semi-infinite domain. The explicit formulation of problems with axial symmetry is also given.

Figures

Tables

Table 1. Table 1: Absolute error of the 2D non-singular boundary integral mehtod for solving Cases I, II and III given by Eqs. (19)-(21) for the domain in Fig.  3 for different corner angle β 𝛽 \beta .

Domain angle	Case I	Case II	Case III
$\beta ={45}^{\text{o}}$	0.0046%	0.29%	0.87%
$\beta ={90}^{\text{o}}$	0.012%	0.019%	0.77%

Equations

$${\nabla }^2\varphi =0.$$

$$c_0\varphi ({\mathbf{x}}_0)+{\int }_S\varphi (\mathbf{x})\frac{\partial G}{\partial n}\text{d}S(\mathbf{x})={\int }_S\frac{\partial \varphi }{\partial n}G(\mathbf{x},{\mathbf{x}}_0)\text{d}S(\mathbf{x}),$$

$$\psi (\mathbf{x})\equiv \varphi ({\mathbf{x}}_0)+{\left(\frac{\partial \varphi }{\partial n}\right)}_{0}f(\mathbf{x})$$

$${\nabla }^{2}f(\mathbf{x})=0,f({\mathbf{x}}_0)=0,\nabla f({\mathbf{x}}_0)\cdot {\mathbf{n}}_0=1.$$

$${\int }_S\left[\varphi (\mathbf{x})-\varphi ({\mathbf{x}}_0)-{\left(\frac{\partial \varphi }{\partial n}\right)}_{0}f(\mathbf{x})\right]\frac{\partial G(\mathbf{x},{\mathbf{x}}_0)}{\partial n}\text{d}S(\mathbf{x})$$

$$={\int }_S\left[\frac{\partial \varphi }{\partial n}-{\left(\frac{\partial \varphi }{\partial n}\right)}_0\nabla f(\mathbf{x})\cdot \mathbf{n}\right]G(\mathbf{x},{\mathbf{x}}_0)\text{d}S(\mathbf{x}).$$

$$\psi (\mathbf{x})$$

$$\psi (r,z)$$

$$\frac{\partial \psi (r,z)}{\partial n}$$

$$\int \frac{\chi (r,z)rE(m)}{(1-m){\bar{R}}^3}\left[(r-r_0)n_r+(z-z_0)n_z-2r_0{n}_r(1-m)/m\right]\text{d}\Gamma$$

$$=\int \frac{2r{r}_0{n}_r\chi (r,z)K(m)}{m{\bar{R}}^3}\text{d}\Gamma +\int \frac{\partial \chi (r,z)}{\partial n}\frac{rK(m)}{\bar{R}}\text{d}\Gamma$$

$$m\equiv \frac{4r{r}_o}{(r+r_0{)}^2+(z-z_0{)}^2}\equiv \frac{4r{r}_o}{{\bar{R}}^2},0\le m\le 1.$$

$$4\pi \left\{\varphi ({\mathbf{x}}_p)-\left[\varphi ({\mathbf{x}}_0)+{\left(\frac{\partial \varphi }{\partial n}\right)}_{0}f({\mathbf{x}}_p)\right]\right\}$$

$$+\int \left\{\varphi (\mathbf{x})-\left[\varphi ({\mathbf{x}}_0)+{\left(\frac{\partial \varphi }{\partial n}\right)}_{0}f(\mathbf{x})\right]\right\}\frac{\partial G(\mathbf{x},{\mathbf{x}}_p)}{\partial n}\text{d}S(\mathbf{x})$$

$$=\int \left\{\frac{\partial \varphi }{\partial n}-{\left(\frac{\partial \varphi }{\partial n}\right)}_0\nabla f(\mathbf{x})\cdot \mathbf{n}(\mathbf{x})\right\}G(\mathbf{x},{\mathbf{x}}_p)\text{d}S(\mathbf{x}).$$

$$\varphi ({\mathbf{x}}_p)=\varphi ({\mathbf{x}}_0)+{\left(\frac{\partial \varphi }{\partial n}\right)}_{0}f({\mathbf{x}}_p)$$

$$-\frac{1}{4\pi }\int \left\{\varphi (\mathbf{x})-\left[\varphi ({\mathbf{x}}_0)+{\left(\frac{\partial \varphi }{\partial n}\right)}_{0}f(\mathbf{x})\right]\right\}\left\{\frac{\partial G(\mathbf{x},{\mathbf{x}}_p)}{\partial n}-\frac{\partial G(\mathbf{x},{\mathbf{x}}_0)}{\partial n}\right\}\text{d}S(\mathbf{x})$$

$$+\frac{1}{4\pi }\int \left\{\frac{\partial \varphi }{\partial n}-{\left(\frac{\partial \varphi }{\partial n}\right)}_0\nabla f(\mathbf{x})\cdot \mathbf{n}(\mathbf{x})\right\}\left\{G(\mathbf{x},{\mathbf{x}}_p)-G(\mathbf{x},{\mathbf{x}}_0)\right\}\text{d}S(\mathbf{x}).$$

$$\psi (\mathbf{x})=\varphi ({\mathbf{x}}_0)+{\left(\frac{\partial \varphi }{\partial n}\right)}_0^L{f}^L(\mathbf{x})+{\left(\frac{\partial \varphi }{\partial n}\right)}_0^R{f}^R(\mathbf{x})$$

$${\nabla }^2{f}^L(\mathbf{x})=0,f^L({\mathbf{x}}_0)=0,\nabla f^L({\mathbf{x}}_0)\cdot {\mathbf{n}}_0^L=1,\nabla f^L({\mathbf{x}}_0)\cdot {\mathbf{n}}_0^R=0.$$

$${\nabla }^2{f}^R(\mathbf{x})=0,f^R({\mathbf{x}}_0)=0,\nabla f^R({\mathbf{x}}_0)\cdot {\mathbf{n}}_0^R=1,\nabla f^R({\mathbf{x}}_0)\cdot {\mathbf{n}}_0^L=0.$$

$$\{\begin{array}{rl}& \nabla \varphi ({\mathbf{x}}_0)={\left(\frac{\partial \varphi }{\partial n}\right)}_0^L{\mathbf{n}}_0^L+{\left(\frac{\partial \varphi }{\partial t}\right)}_0^L{\mathbf{t}}_0^L,\\ & \nabla \varphi ({\mathbf{x}}_0)={\left(\frac{\partial \varphi }{\partial n}\right)}_0^R{\mathbf{n}}_0^R+{\left(\frac{\partial \varphi }{\partial t}\right)}_0^R{\mathbf{t}}_0^R;\end{array}$$

$$f^L(\mathbf{x})=\frac{-({\mathbf{n}}_0^L\cdot {\mathbf{n}}_0^R){\mathbf{n}}_0^R+{\mathbf{n}}_0^L}{1-({\mathbf{n}}_0^L\cdot {\mathbf{n}}_0^R{)}^2}\cdot (\mathbf{x}-{\mathbf{x}}_0),$$

$$f^R(\mathbf{x})=\frac{-({\mathbf{n}}_0^R\cdot {\mathbf{n}}_0^L){\mathbf{n}}_0^L+{\mathbf{n}}_0^R}{1-({\mathbf{n}}_0^R\cdot {\mathbf{n}}_0^L{)}^2}\cdot (\mathbf{x}-{\mathbf{x}}_0).$$

$$\varphi =1-x,\text{Case I};$$

$$\varphi =1-xy,\text{Case II};$$

$$\varphi =\frac{\sinh (\pi y)}{\sinh (\pi )}\sin (\pi x)+1,\text{Case III}.$$

$$f^L(\mathbf{x})=\frac{|{\mathbf{x}}_0-{\mathbf{x}}_D^L{|}^2}{{\mathbf{n}}_0^L\cdot ({\mathbf{x}}_0-{\mathbf{x}}_D^L)}\ln \left(\frac{|\mathbf{x}-{\mathbf{x}}_D^L|}{|{\mathbf{x}}_0-{\mathbf{x}}_D^L|}\right)$$

$$f^R(\mathbf{x})=\frac{|{\mathbf{x}}_0-{\mathbf{x}}_D^R{|}^2}{{\mathbf{n}}_0^R\cdot ({\mathbf{x}}_0-{\mathbf{x}}_D^R)}\ln \left(\frac{|\mathbf{x}-{\mathbf{x}}_D^R|}{|{\mathbf{x}}_0-{\mathbf{x}}_D^R|}\right)$$

$$\{\begin{array}{r}{\mathbf{n}}_0^L\cdot ({\mathbf{x}}_0-{\mathbf{x}}_D^L)\ne 0,{\mathbf{n}}_0^R\cdot ({\mathbf{x}}_0-{\mathbf{x}}_D^L)=0;\\ {\mathbf{n}}_0^R\cdot ({\mathbf{x}}_0-{\mathbf{x}}_D^R)\ne 0,{\mathbf{n}}_0^L\cdot ({\mathbf{x}}_0-{\mathbf{x}}_D^R)=0.\end{array}$$

Full text

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 [1, 2], the description of waves [3, 4], cavitation or supercavitation phenomena [5, 6], in civil, marine and ocean engineering and oscillating bubble dynamics in sonophysics and sonochemistry [7].

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 [8] can be used to reduce both CPU time and memory requirement from $O(N^{2})$ to $O(N\log N)$. Thus in spite of being a simplification of the full Navier-Stokes 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” [9]. Since the physical problem itself is perfectly well behaved on the boundaries, such singularities are numerical inconveniences generated by the mathematical formulation [10]. Traditionally, the singular behaviour is dealt with by a local change of variables in the evaluation of the surface integrals [11] 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 [10, 12]. Another method to remove the singularities requires finding additional parameters that have to be determined on a fictitious ‘nearby’ boundary [13, 14, 15].

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 well-known 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 semi-infinite 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 Non-singular formulation of the boundary integral method

To develop a non-singular formulation of the boundary integral method, consider the internal problem in a 3D domain that is bounded by the closed surface $S$, as shown in Fig. 1. The potential, $\phi$, is governed by the Laplace equation

[TABLE]

By using the 3D free space Green’s function $G(\boldsymbol{x},\boldsymbol{x}_{0})=1/|\boldsymbol{x}-\boldsymbol{x}_{0}|$ and with the help of Green’s second identity, the solution of Eq. (1) can be found by solving the conventional boundary integral equation

[TABLE]

where $c_{0}$ is the solid angle at the observation point $\boldsymbol{x}_{0}$ on $S$ and the surface element $\text{d}S(\boldsymbol{x})$ is at $\boldsymbol{x}$. The normal derivatives are defined by ${\partial{\phi}}/{\partial{n}}\equiv\boldsymbol{\nabla}\phi(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{n(\boldsymbol{x})}$ and ${\partial{G}}/{\partial{n}}\equiv\boldsymbol{\nabla}G(\boldsymbol{x},\boldsymbol{x}_{0})\boldsymbol{\cdot}\boldsymbol{n(\boldsymbol{x})}$, where $\boldsymbol{n}\equiv\boldsymbol{n(\boldsymbol{x})}$ is the unit normal vector pointing out of the internal domain at $\boldsymbol{x}$ [9]. If either $\phi$ or ${\partial{\phi}}/{\partial{n}}$ or a mixed boundary condition is given on the whole boundary $S$, the corresponding ${\partial{\phi}}/{\partial{n}}$ or/and $\phi$ can be obtained from Eq. (2). Although numerical integration over the singularities in $G(\boldsymbol{x},\boldsymbol{x}_{0})$ and ${\partial{G}}/{\partial{n}}$ can be effected by established methods [11], our aim is to remove such mathematical singularities analytically at the outset.

Corresponding to a given point $\boldsymbol{x}_{0}$ on the boundary, we construct a function $\psi(\boldsymbol{x})$ that also satisfies the Laplace equation and hence Eq. (2), with the properties that as $\boldsymbol{x}\rightarrow\boldsymbol{x}_{0}$, $\psi(\boldsymbol{x})\rightarrow\phi(\boldsymbol{x}_{0})\equiv\phi_{0}$ and ${\partial{\psi}}/{\partial{n}}\rightarrow({\partial{\phi}}/{\partial{n}})_{0}\equiv\boldsymbol{\nabla}\phi(\boldsymbol{x}_{0})\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x}_{0})$, where $\boldsymbol{n}_{0}\equiv\boldsymbol{n}(\boldsymbol{x}_{0})$ is the outward unit normal vector at $\boldsymbol{x}_{0}$. We choose $\psi(\boldsymbol{x})$ to be of the form

[TABLE]

so the function $f(\boldsymbol{x})$ must satisfy

[TABLE]

Taking the difference between the conventional boundary integral equations for $\phi(\boldsymbol{x})$ and for $\psi(\boldsymbol{x})$ we obtain an integral equation relating $\phi(\boldsymbol{x})$ and $\partial{\phi}/\partial{n}$ on $S$ that replaces the conventional boundary integral equation in Eq. (2):

[TABLE]

The key point is that both integrands in Eq. (2) are now non-singular and thus any convenient quadrature method can be used to evaluate the integrals. Note also that the solid angle $c_{0}$ no longer appears. Implicit in the derivation of Eq. (2) is that the outward unit normal $\boldsymbol{n}(\boldsymbol{x}_{0})$ is uniquely defined at $\boldsymbol{x}_{0}$. 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 de-singularisation method using the linear function: $f(\boldsymbol{x})=\boldsymbol{n}_{0}\boldsymbol{\cdot}(\boldsymbol{x}-\boldsymbol{x}_{0})$ has been given elsewhere [16] 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 [17]. This approach has also been extended to de-singularise boundary integral equations that arise in Stokes flow, in solving the Helmholtz equation and the equations of linear elasticity using a linear function for $f(\boldsymbol{x})$ [16].

However, for potential problems involving domains of infinite extent (external problems), a different choice of $f(\boldsymbol{x})$ is required since the linear function is unbounded at infinity. One possible choice for $f(\boldsymbol{x})$ is

[TABLE]

The constant vector $\boldsymbol{x}_{D}$ is the position of any convenient point that is located outside the solution domain and satisfies $\boldsymbol{n}_{0}\boldsymbol{\cdot}(\boldsymbol{x}_{0}-\boldsymbol{x}_{D})\neq 0$. In the next section, we will use this form of $\psi(\boldsymbol{x})$ to formulate a non-singular BEM for axisymmetric problems and in Section 6, we demonstrate how Eq. (6) can be used to formulate a non-singular boundary integral problem in a semi-infinite domain to solve a 3D time-dependent 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 $f(\boldsymbol{x})$ to construct non-singular versions of the boundary integral equation. In fact, any form of $f(\boldsymbol{x})$ 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 $G$ and ${\partial{G}}/{\partial{n}}$ at $\boldsymbol{x}=\boldsymbol{x}_{0}$.

3 Non-singular axisymmetric boundary integral equation

For problems that possess axial symmetry whereby in the cylindrical variables: $r,\theta$ and $z$, we have $\phi(\boldsymbol{x})=\phi(r,z)$ and $\partial{\phi(\boldsymbol{x})}/\partial{n}=\partial{\phi(r,z)}/\partial{n}$, the integration over the azimuthal angle, $\theta$, can be evaluated analytically [9]. Since a point $(r,z)$ on the axisymmetric boundary surface is specified by some given equation $S(r,z)=0$, the surface integrals for the axisymmetric case can be reduced to 1D integrals.

The $\psi(\boldsymbol{x})$ given by Eq. (6) will be an axisymmetric function if we choose (in Cartesian coordinates) $\boldsymbol{x}_{D}=(0,0,z_{D})$ to lie on the $z$-axis of symmetry with $z_{D}$ located outside the solution domain. We measure the azimuthal angle, $\theta$, relative to $\boldsymbol{x}_{0}=(r_{0},0,z_{0})$ so that $\boldsymbol{x}=(r\cos\theta,r\sin\theta,z)$ and the surface normals are given by $\boldsymbol{n}_{0}=(n_{r0},0,n_{z0})$ and $\boldsymbol{n}=(n_{r}\cos\theta,n_{r}\sin\theta,n_{z})$. With these definitions, the axisymmetric function $\psi(\boldsymbol{x})$ in Eq. (6) and its normal derivative are given explicitly by

[TABLE]

where $\rho\equiv\sqrt{r^{2}+(z-z_{D})^{2}}$, $\rho_{0}\equiv\sqrt{r_{0}^{2}+(z_{0}-z_{D})^{2}}$, $s\equiv rn_{r}+(z-z_{D})n_{z}$ and $s_{0}\equiv r_{0}n_{r0}+(z_{0}-z_{D})n_{z0}$.

The axisymmetric version of the non-singular boundary integral equation obtained after performing the $\theta$-integration in Eq. (2) using Eq. (6) is (see also [9])

[TABLE]

where $\chi(r,z)\equiv\phi(r,z)-\psi(r,z)$. The arc length element $\text{d}{\Gamma}$ may be expressed as $\sqrt{1+(\text{d}r/\text{d}z)^{2}}\text{ d}z$ using the equation: $S(r,z)=0$ that defines the axisymmetric surface. The quantities $K(m)$ and $E(m)$ are the complete elliptic integrals of the first and second kind of the parameter, $m$ (see [9, 18]), with

[TABLE]

The axisymmetric boundary integral equation given by Eq. (3) is non-singular because the $\log(1-m)$ divergence in $K(m)$ as $m\rightarrow 1$ in the limit $\boldsymbol{x}\rightarrow\boldsymbol{x}_{0}$ is now suppressed by the terms containing the difference between $\phi$ and $\psi$ and between ${\partial{\phi}}/{\partial{n}}$ and ${\partial{\psi}}/{\partial{n}}$ that vanish as $(1-m)$ as $\boldsymbol{x}\rightarrow\boldsymbol{x}_{0}$. Thus the line integral in Eq. (3) can be evaluated using any quadrature method. This is a consequence of using an axisymmetric form of $\psi(r,z)$ given by Eq. (6) that simplified to Eqs. (7) and (8).

An additional bonus of the non-singular 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 [9].

4 Robust method to calculate the potential near a boundary

The method of de-singularising 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, $\phi(\boldsymbol{x}_{p})$, at an observation point $\boldsymbol{x}_{p}$ 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 $[\phi(\boldsymbol{x}_{p})-\psi(\boldsymbol{x}_{p})]$, with $\psi(\boldsymbol{x})$ given by Eq. (3) and $c_{0}=4\pi$ for points inside the domain,

[TABLE]

The point $\boldsymbol{x}_{0}$ is located on the boundary and its relation to $\boldsymbol{x}_{p}$ is specified below. The nearly singular behaviour of the integrands when $\boldsymbol{x}_{p}$ is close to the boundary can be eliminated by subtracting Eq. (2) from Eq. (4) to give

[TABLE]

When the observation point $\boldsymbol{x}_{p}$ is near the boundary, a suitable choice for the boundary point, $\boldsymbol{x}_{0}$, is to project $\boldsymbol{x}_{p}$ onto the boundary using the relation $\boldsymbol{x}_{p}=\boldsymbol{x}_{0}-\varepsilon\boldsymbol{n}_{0}$, with $\varepsilon>0$.

Another important benefit arising from the numerical robustness of the present de-singualarised formulation of the BIM occurs when different parts of the domain boundary become close together as often happens in multi-scale moving boundary problems. A simple example of this is when the boundary comprises of two nearly touching spheres and the point $\boldsymbol{x}_{0}$ 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 $\boldsymbol{x}_{0}$. A numerical demonstration of this is given in Section 6.

5 Corner nodes in 2D

The above formulation of the non-singular form of the BEM for the potential problem assumes that the unit normal $n_{0}\equiv n(\boldsymbol{x}_{0})$ is uniquely defined at $\boldsymbol{x}_{0}$. When $\boldsymbol{x}_{0}$ is a corner node on a 2D boundary, we choose $\psi(\boldsymbol{x})$ to have the form

[TABLE]

where $\phi_{0}\equiv\phi(\boldsymbol{x}_{0})$ is the potential at $\boldsymbol{x}_{0}$, and the constants $(\partial{\phi}/\partial{n})^{L}_{0}$ and $(\partial{\phi}/\partial{n})^{R}_{0}$ are the values of normal derivatives as one approaches $\boldsymbol{x}_{0}$ from the ‘left’ ($L$) and from the ‘right’ ($R$). The functions $f^{L}(\boldsymbol{x})$ and $f^{R}(\boldsymbol{x})$ satisfy the following conditions

[TABLE]

[TABLE]

The value of $\nabla\phi(\boldsymbol{x}_{0})$ at the corner $\boldsymbol{x}_{0}$ can be expressed in two equivalent forms

[TABLE]

in which, $(\partial{\phi}/\partial{t})^{L}_{0}$ and $(\partial{\phi}/\partial{t})^{R}_{0}$ are the two tangential derivatives of the potential on side $L$ and $R$ at $\boldsymbol{x}_{0}$, respectively, and $\boldsymbol{t}^{L}_{0}$ and $\boldsymbol{t}^{R}_{0}$ are the unit tangential vectors along side $L$ and $R$ at $\boldsymbol{x}_{0}$, respectively. This compatibility condition in Eq. (16) provides an additional relation between $(\partial{\phi}/\partial{n})^{L}_{0}$ and $(\partial{\phi}/\partial{n})^{R}_{0}$  [19] in the formulation of the BIM whereas the tangential derivatives can be constructed from the values of $\phi$ at neighbouring nodes.

A choice for $f^{L}(\boldsymbol{x})$ and $f^{R}(\boldsymbol{x})$ for a 2D corner problem can be constructed in terms of the surface normals $\boldsymbol{n}^{L}_{0}$ and $\boldsymbol{n}^{R}_{0}$ on either side of the corner at $\boldsymbol{x}_{0}$

[TABLE]

and

[TABLE]

6 Numerical demonstrations and examples

6.1 Two nearly touching spheres

To demonstrate the utility and robustness of our non-singular 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 $R$ and $3R/2$ and are at a separation of $10^{-4}R$ at the point of closest approach (Fig. 2a). We compare the velocity potential obtained using the standard axisymmetric boundary integral method [20], with our non-singular version given by Eq. (3) with $\boldsymbol{x}_{D}$ 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 non-singular 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 non-singular formulation, such effects do not arise.

6.2 Corner problem in 2D

To illustrate the implementation of our non-singular 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, $\beta$. At the corner nodes, the function $\psi(\boldsymbol{x})$ given by Eq. (13) used in the non-singular boundary integral equation (2) and the double node technique [19] 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 $\phi$ along the boundaries. For the remaining nodes on the edges, the form of $\psi(\boldsymbol{x})$ in Eq. (3) with $f(\boldsymbol{x})=\boldsymbol{n}_{0}\boldsymbol{\cdot}(\boldsymbol{x}-\boldsymbol{x}_{0})$ is used in Eq. (2).

We tested three cases corresponding to analytical solutions:

[TABLE]

[TABLE]

[TABLE]

We solved the above corner problems with Dirichlet boundary conditions. The comparisons between the results for the normal derivatives obtained by our non-singular BIM with 21 nodes on each edge and the analytical solutions are shown in Fig. 4 when $\beta=45^{\text{o}}$ and $\beta=90^{\text{o}}$. The absolute errors for these cases can be found in Table 1.

Other pairs of $f^{L}(\boldsymbol{x})$ and $f^{R}(\boldsymbol{x})$ that satisfy Eqs. (14) and (15) can also be used, for example

[TABLE]

[TABLE]

where $\boldsymbol{x}^{L}_{D}$ and $\boldsymbol{x}^{R}_{D}$ are located outside of the calculation domain, and

[TABLE]

This pair of $f^{L}(\boldsymbol{x})$ and $f^{R}(\boldsymbol{x})$ are also used to solve the above three cases with Dirichlet boundary conditions when $\beta=90^{\text{o}}$. Once again, 21 nodes are employed on each edge. The absolute errors between the results for the normal derivatives obtained by our non-singular 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 semi-infinite deformable boundary

To illustrate the generality of our non-singular formulation of the BEM, we examine the following wave drag problem that has sufficient complexity to be interesting [21]. Consider two identical spheres of radius, $R$, separated by a constant distance $D$ between their centres and moving with constant speed $\boldsymbol{u}_{0}=U_{0}\boldsymbol{i}$ at the same depth $H$ 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 $z=0$) with the surface elevation $z$ 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 time-dependent 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, $\boldsymbol{u}=\boldsymbol{\nabla}\phi$, generated by the moving spheres is found by solving the Laplace equation for the velocity potential, $\phi$ at each time step. On the spheres, the velocity is given and the potential is calculated. Initially, $\phi=0$ 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 $\partial{\phi}/\partial{n}$ at the free surface that is then used to predict its shape at the next time step. The value of $\phi$ on the free surface at the next time step is found from the unsteady Bernoulli equation evaluated on the surface: $\rho\text{D}\phi/\text{D}t=\rho u^{2}/2-\rho gz$ where $\text{D}/\text{D}t=\partial{}/\partial{t}+\boldsymbol{u}\cdot\boldsymbol{\nabla}$ is the material derivative, $\rho$ the fluid density and $g$ the gravitational acceleration. Such time-stepping then gives the spatio-temporal 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 $p$ on the sphere surfaces.

The deformable interface extends to infinity in the $xy$-direction where it asymptotes to a flat surface. Therefore far from the spheres, the potential and its derivative vanish asymptotically like $1/|\boldsymbol{x}|$ and $1/|\boldsymbol{x}|^{2}$ as $|\boldsymbol{x}|\rightarrow\infty$. We therefore assume that beyond a radius of 40$R$ 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 [22, 23]. Since Eq. (2) is non-singular, 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 [23]. Finally, the integral over the half spherical surface at infinity will contribute two terms $2\pi\phi(\boldsymbol{x}_{0})$ and $2\pi(\partial{\phi}/\partial{n})_{0}|\boldsymbol{x}_{0}-\boldsymbol{x}_{D}|^{2}/[\boldsymbol{n}_{0}\boldsymbol{\cdot}(\boldsymbol{x}_{0}-\boldsymbol{x}_{D})]$ to the left-hand 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 $\boldsymbol{x}_{0}$ is on a sphere, $\boldsymbol{x}_{D}$ in Eq. (6) is taken to be the centre of the sphere. When $\boldsymbol{x}_{0}$ is on the free surface, we choose the Cartesian components of $\boldsymbol{x}_{D}=(x_{D},y_{D},z_{D})$ to be $z_{D}=3R$ and any convenient $x_{D}$ and $y_{D}$ that ensure $\boldsymbol{n}_{0}\boldsymbol{\cdot}(\boldsymbol{x}_{0}-\boldsymbol{x}_{D})\neq 0$.

With a constant time step of $0.02R/U_{0}$, the position of each surface node is updated with $\text{D}\boldsymbol{x}/\text{D}t=\boldsymbol{u}\equiv\boldsymbol{\nabla}\phi(\boldsymbol{x})$ and the spheres are translated using $\boldsymbol{u}_{0}$. The calculation continues until the spheres have travelled a distance of $30R$. The elastic mesh technique [25] is applied on the free surface to ensure a uniform mesh even after many time-steps.

The drag force acting on each sphere is calculated by $\boldsymbol{F}=\int_{S}p\text{ }\boldsymbol{n}\text{ d}S$ where $S$ is the surface of the sphere. In Fig. 6a we present the time variation of the force $\boldsymbol{F}=F_{x}\boldsymbol{i}$ on each of the two spheres travelling with constant velocity $\boldsymbol{u}_{0}=U_{0}\boldsymbol{i}$. The three non-dimensional parameters that govern this problem are: $H/R$, $D/R$ and the Froude number, $Fr=U_{0}/\sqrt{gH}$ (Fig. 5). For our illustrative example, we chose a sphere separation of $D/R=2.4$, a submerged depth of $H/R=2$, and $Fr=1$. Also shown are results for $H/R=\infty$ that corresponds to the absence of the deformable interface and the force on an isolated single sphere that corresponds to $D/R=\infty$. The leading sphere experiences a retarding force in the $x$-direction whereas the trailing sphere experiences a force in the direction of travel. At $H/R=\infty$, 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 non-zero wave drag on the pair of moving spheres [21]. 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 non-singular 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, $D$ at $H/D=\infty$. There is excellent agreement between our non-singular boundary integral method and the analytic result of [24].

Snapshots of the spatio-temporal 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 [21].

7 Conclusions

In this paper, we have provided a fully non-singular 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 [9], the amount of computer code needed to implement the non-singular 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 [9]. 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.