Non-singular boundary integral methods for fluid mechanics applications

TL;DR

This paper introduces a boundary integral method that analytically removes singularities, enabling more accurate and efficient solutions for various fluid mechanics problems like potential flow, Helmholtz, Stokes, and elasticity equations.

Contribution

It develops a non-singular boundary integral formulation applicable to multiple PDEs in fluid and continuum mechanics, improving computational robustness.

Findings

Successfully removes singularities analytically from integral equations.

Applicable to Laplace, Helmholtz, Stokes, and elasticity problems.

Demonstrates broad applicability in practical fluid mechanics scenarios.

Abstract

A formulation of the boundary integral method for solving partial differential equations has been developed whereby the usual weakly singular integral and the Cauchy principal value integral can be removed analytically. The broad applicability of the approach is illustrated with a number of problems of practical interest to fluid and continuum mechanics including the solution of the Laplace equation for potential flow, the Helmholtz equation as well as the equations for Stokes flow and linear elasticity.

Equations

$$4\pi \varphi ({\mathbf{x}}_0)+{\int }_S\varphi (\mathbf{x})\nabla G(\mathbf{x},{\mathbf{x}}_0)\cdot \mathbf{n}\text{ d}S(\mathbf{x})={\int }_{S}G(\mathbf{x},{\mathbf{x}}_0)\nabla \varphi (\mathbf{x})\cdot \mathbf{n}\text{ d}S(\mathbf{x}).$$

$$(4\pi -c)\varphi ({\mathbf{x}}_0)+{\int }_{S,\text{PV}}\varphi (\mathbf{x})\nabla G(\mathbf{x},{\mathbf{x}}_0)\cdot \mathbf{n}\text{ d}S(\mathbf{x})={\int }_{S}G(\mathbf{x},{\mathbf{x}}_0)\nabla \varphi (\mathbf{x})\cdot \mathbf{n}\text{ d}S(\mathbf{x}),$$

$$\psi (\mathbf{x})=\varphi ({\mathbf{x}}_0)+\mathbf{a}({\mathbf{x}}_0)\cdot (\mathbf{x}-{\mathbf{x}}_0),$$

$$4\pi \psi ({\mathbf{x}}_0)+{\int }_S\psi (\mathbf{x})\nabla G(\mathbf{x},{\mathbf{x}}_0)\cdot \mathbf{n}\text{ d}S(\mathbf{x})={\int }_{S}G(\mathbf{x},{\mathbf{x}}_0)\nabla \psi (\mathbf{x})\cdot \mathbf{n}\text{ d}S(\mathbf{x}),$$

$${\int }_S[\varphi (\mathbf{x})-\psi (\mathbf{x})]\nabla G(\mathbf{x},{\mathbf{x}}_0)\cdot \mathbf{n}\text{ d}S(\mathbf{x})={\int }_{S}G(\mathbf{x},{\mathbf{x}}_0)\nabla [\varphi (\mathbf{x})-\psi (\mathbf{x})]\cdot \mathbf{n}\text{ d}S(\mathbf{x}).$$

$$\mathbf{a}({\mathbf{x}}_0)=[\nabla \varphi ({\mathbf{x}}_0)\cdot {\mathbf{n}}_0]{\mathbf{n}}_0\equiv {\left(\frac{\partial \varphi }{\partial n}\right)}_0{\mathbf{n}}_0,$$

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

$$=$$

$${\int }_S[\varphi (\mathbf{x})-\psi (\mathbf{x})]\nabla H(\mathbf{x},{\mathbf{x}}_0)\cdot \mathbf{n}\text{ d}S(\mathbf{x})={\int }_{S}H(\mathbf{x},{\mathbf{x}}_0)\nabla [\varphi (\mathbf{x})-\psi (\mathbf{x})]\cdot \mathbf{n}\text{ d}S(\mathbf{x}),$$

$$\psi (\mathbf{x})=\varphi ({\mathbf{x}}_0)\cos [k{\mathbf{n}}_0\cdot (\mathbf{x}-{\mathbf{x}}_0)]+\frac{b({\mathbf{x}}_0)}{k}\sin [k{\mathbf{n}}_0\cdot (\mathbf{x}-{\mathbf{x}}_0)].$$

$$b({\mathbf{x}}_0)=\nabla \varphi ({\mathbf{x}}_0)\cdot {\mathbf{n}}_0\equiv {\left(\frac{\partial \varphi }{\partial n}\right)}_0.$$

$${\int }_S\left\{\varphi (\mathbf{x})-\varphi ({\mathbf{x}}_0)\cos \chi -\frac{1}{k}{\left(\frac{\partial \varphi }{\partial n}\right)}_0\sin \chi \right\}\nabla H(\mathbf{x},{\mathbf{x}}_0)\cdot \mathbf{n}\text{d}S(\mathbf{x})$$

$$=$$

$$-\nabla p+\mu {\nabla }^2\mathbf{u}=0\text{and}\nabla \cdot \mathbf{u}=0.$$

$$\text{mathsfbi}{\sigma }_{ik}=-p\text{mathsfbi}{\delta }_{ik}+\mu \left[\frac{\partial u_i}{\partial x_k}+\frac{\partial u_k}{\partial x_i}\right],$$

$$8\pi u_j^0+{\int }_S{u}_i\text{mathsfbi}{T}_{ijk}{n}_k\text{ d}S$$

$$\text{mathsfbi}{U}_{ij}(\mathbf{x},{\mathbf{x}}_0)=\frac{\text{mathsfbi}{\delta }_{ij}}{r}+\frac{\hat{x}{}_i\hat{x}{}_j}{r^3},$$

$$\text{mathsfbi}{T}_{ijk}(\mathbf{x},{\mathbf{x}}_0)=-6\frac{\hat{x}{}_i\hat{x}{}_j\hat{x}{}_k}{r^5},$$

$$\mathbf{w}(\mathbf{x})=\mathbf{u}({\mathbf{x}}_0)+\frac{1}{\mu }\text{mathsfbi}\mathbf{M}({\mathbf{x}}_0)\cdot (\mathbf{x}-{\mathbf{x}}_0),$$

$$\mathbf{\Sigma }({\mathbf{x}}_0)=\mathbf{M}({\mathbf{x}}_0)+{\mathbf{M}}^T({\mathbf{x}}_0)\text{ or }\text{mathsfbi}{\Sigma }_{ik}^0=\text{mathsfbi}{M}_{ik}^0+\text{mathsfbi}{M}_{ki}^0.$$

$$\text{Tr}[\text{mathsfbi}\mathbf{M}({\mathbf{x}}_0)]=\text{Tr}[\text{mathsfbi}\mathbf{\Sigma }({\mathbf{x}}_0)]/2=0$$

$${\int }_S\left[u_i-u_i^0-\frac{1}{\mu }\text{mathsfbi}{M}_{il}^0(x_l-x_l^0)\right]\text{mathsfbi}{T}_{ijk}{n}_k\text{ d}S=\frac{1}{\mu }{\int }_S\left(f_i-\text{mathsfbi}{\Sigma }_{il}^0{n}_l\right)\text{mathsfbi}{U}_{ij}\text{ d}S.$$

$$\mathbf{M}({\mathbf{x}}_0)=\mathbf{f}({\mathbf{x}}_0)\mathbf{n}({\mathbf{x}}_0)-\frac{1}{4}\mathbf{f}({\mathbf{x}}_0)\cdot \mathbf{n}({\mathbf{x}}_0)\left[\mathbf{I}+\mathbf{n}({\mathbf{x}}_0)\mathbf{n}({\mathbf{x}}_0)\right],$$

$$\text{mathsfbi}{M}_{il}^0=f_i^0{n}_l^0-\frac{1}{4}(f_k^0{n}_k^0)(\text{mathsfbi}{\delta }_{il}+n_i^0{n}_l^0),$$

$$\mathbf{\Sigma }({\mathbf{x}}_0)=\mathbf{f}({\mathbf{x}}_0)\mathbf{n}({\mathbf{x}}_0)+\mathbf{n}({\mathbf{x}}_0)\mathbf{f}({\mathbf{x}}_0)-\frac{1}{2}\mathbf{f}({\mathbf{x}}_0)\cdot \mathbf{n}({\mathbf{x}}_0)\left[\mathbf{I}+\mathbf{n}({\mathbf{x}}_0)\mathbf{n}({\mathbf{x}}_0)\right],$$

$$\text{mathsfbi}{\Sigma }_{il}^0=\text{mathsfbi}{M}_{il}^0+\text{mathsfbi}{M}_{li}^0=(f_i^0{n}_l^0+f_l^0{n}_i^0)-\frac{1}{2}(f_k^0{n}_k^0)(\text{mathsfbi}{\delta }_{il}+n_i^0{n}_l^0).$$

$$\text{mathsfbi}{\epsilon }_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right).$$

$$\text{mathsfbi}{\sigma }_{ij}=\frac{2\text{mathsfbi}\mu \nu }{1-2\nu }\text{mathsfbi}{\delta }_{ij}\text{mathsfbi}{\epsilon }_{ll}+2\mu \text{mathsfbi}{\epsilon }_{ij},$$

$$\text{mathsfbi}{U}_{ij}=\frac{1}{2(1-\nu )}\left[(3-4\nu )\frac{\text{mathsfbi}{\delta }_{ij}}{r}+\frac{\hat{x}{}_i\hat{x}{}_j}{r^3}\right],$$

$$\text{mathsfbi}{T}_{ijk}=-\frac{1}{1-\nu }\left[3\frac{\hat{x}{}_i\hat{x}{}_j\hat{x}{}_k}{r^5}+\frac{1-2\nu }{r^3}\left(-{\delta }_{ij}{\hat{x}}_k+{\delta }_{jk}{\hat{x}}_i+{\delta }_{ki}{\hat{x}}_j\right)\right].$$

$$S=\zeta +\frac{1}{2}{a}_s{\xi }^2+\frac{1}{2}{b}_s{\eta }^2=0$$

Full text

††volume: 696

1 Introduction

The boundary integral formulation is an efficient method of representing the solutions of certain linear partial differential equations by reducing the dimensionality of the problem by one. The solution in a volume or area domain is represented in terms of an integral over surface(s) or line(s) that enclose the domain. Solutions to fluid dynamics problems that can be modelled by the Laplace equation for potential flow or the Stokes equation for low-Reynolds-number flow as well as continuum mechanics problems such as the Helmholtz equation in scattering problems or problems in linear elasticity can all be represented in terms of such boundary integrals (Becker, 1992). Symm (1963) provided a practical way to solve the integral equations by treating the surfaces or lines as discrete elements. Since the 1970s the boundary integral method (BIM) has gained increasing prominence (see Cheng & Cheng, 2005, for a historical overview). The advantage of the BIM is self-evident. The reduction in the dimensionality of the problem from a volume (surface) mesh to a surface (line) mesh provides a substantial gain in computational efficiency. However, this gain is offset by the fact that the numerical implementation of the BIM is not straightforward because the approach is plagued by ‘a mathematical monster that leaps out of every page’ (Becker, 1992). In essence, the boundary element formulation uses the Green’s function that has a $1/r$ divergence and a $1/r^{2}$ divergence in its derivative around the source point. The integral over the $1/r$ divergence gives rise to a weak singularity that can be evaluated using semi-analytical techniques. The term from the $1/r^{2}$ divergence gives rise to a Cauchy principal value (PV) integral that requires careful numerical treatment.

In this communication, we develop a general non-singular boundary integral formulation that is applicable to the Laplace equation for the potential problem, the Helmholtz equation, and equations associated with Stokes flow and linear elastic deformations. The approach is based on removing the singularities in the BIM formulation by subtracting the solution of a special related problem. We demonstrate the details of our approach using the potential problem from which it is easy to see how the method can be extended to the more complicated cases of Stokes flow and linearly elastic deformations. Validation of the approach is obtained by comparing numerical results for problems in Stokes flow for which analytic results are known. This non-singular boundary integral formulation simplifies numerical solutions based on this popular technique.

2 Potential problem — non-singular boundary integral formulation

In fluid dynamics, the potential problem arises in incompressible, inviscid or high-Reynolds-number flows. The velocity field $\boldsymbol{u}=\boldsymbol{\nabla}\phi$ can be expressed in terms of a scalar potential that satisfies the Laplace equation $\boldsymbol{\nabla}^{2}\phi=0$. The corresponding free space Green’s function $G(\boldsymbol{x},\boldsymbol{x}_{0})=1/r$, $r=|\boldsymbol{x}-\boldsymbol{x}_{0}|$ satisfies $\boldsymbol{\nabla}^{2}G=-4\pi\delta(\boldsymbol{x}-\boldsymbol{x}_{0})$, where $\delta(\boldsymbol{x}-\boldsymbol{x}_{0})$ is the Dirac $\delta$-function. Consider the solution in the fluid domain enclosed by the surface $S$. With the help of Green’s second identity (Becker, 1992), the solution at $\boldsymbol{x}_{0}$ inside the domain can be written as the surface integral (see Figure 1a)

[TABLE]

This integral relates the potential $\phi$ at $\boldsymbol{x}_{0}$ inside the domain to integrals over $\phi$ and its normal derivative $\boldsymbol{\nabla}\phi\boldsymbol{\cdot}\boldsymbol{n}=\partial{\phi}/\partial{n}$, on the surface $S$ with $\boldsymbol{n}$ being the outward unit normal. The vector $\boldsymbol{x}$ points to the integration position on the surface $S$. By letting $\boldsymbol{x}_{0}$ onto the surface $S$, we have an equation that can be solved for $\phi$ (or $\partial{\phi}/\partial{n}$) on the surface if $\partial{\phi}/\partial{n}$ (or $\phi$) is specified. This corresponds to the Dirichlet (or Neumann) problem. This is the boundary integral formulation. However, with $\boldsymbol{x}_{0}$ on the surface, then as $\boldsymbol{x}\rightarrow\boldsymbol{x}_{0}$, the integral involving $G$ with a $1/r$ singularity (the single layer term) has a weak singularity that can be handled numerically by changing to local polar coordinates on the surface whereas the integral involving $\boldsymbol{\nabla}G$ with a $1/r^{2}$ singularity (the double layer term) gives rise to a Cauchy principal value (PV) integral and a Dirac $\delta$-function contribution. Thus with $\boldsymbol{x}_{0}$ on the surface $S$ as illustrated in Figure 1b, the boundary integral equation that needs to be solved is

[TABLE]

where $c$ is the solid angle subtended at $\boldsymbol{x}_{0}$ with $c=2\pi$ if the surface has a defined curvature at $\boldsymbol{x}_{0}$. The numerical evaluation of the weak singularity associated with $G$ and the Cauchy principal value integral associated with $\boldsymbol{\nabla}G$ in Eq. (2) requires special considerations as ordinary integration methods such as Gaussian quadrature can no longer be used (Becker 1992). Different numerical methods have been developed to handle these singularities (see for example Lean & Wexler (1985), Bazhlekov et al. (2004)). Thus, if either the potential or the normal velocity, $\partial{\phi}/\partial{n}$, is known on the surface $S$, the other unknown quantity can be calculated (see for example Gonzalez-Avila et al. (2011), Fong et al. (2009), Blake et al. (1986), Wrobel (2002), Wang (1998) and Zhang et al. (2001)). The aim therefore is to avoid the numerical effort needed when having to deal with these singularities.

We recapitulate the earlier work of Klaseboer et al. (2009) and show that both singular terms associated with $G$ and $\boldsymbol{\nabla}G$ in Eq. (2) can be removed by considering the linear potential function

[TABLE]

where the vector $\boldsymbol{a}(\boldsymbol{x}_{0})$ will be chosen to eliminate the singularities in Eq. (2). Clearly $\psi(\boldsymbol{x})$ satisfies $\boldsymbol{\nabla}^{2}\psi=0$ and the Green’s identity, Eq. (1):

[TABLE]

with $\boldsymbol{x}_{0}$ inside the domain. Thus subtracting Eq. (4) from Eq. (1) and using Eq. (3) gives

[TABLE]

When $\boldsymbol{x}_{0}$ is located on the surface $S$, the singularities in Eq. (5) can be eliminated with the following choice of the vector

[TABLE]

that depends on the point $\boldsymbol{x}_{0}$, where the outward unit normal is $\boldsymbol{n}_{0}\equiv\boldsymbol{n}(\boldsymbol{x}_{0})$ (see Figure 1b). Finally the required result of the non-singular formulation of the boundary integral equation, with $\boldsymbol{x}_{0}$ now located on the surface S, takes the form

[TABLE]

This result supersedes the traditional form of the boundary integral formulation given in Eq. (2) because all singularities have now been removed. In particular, as $\boldsymbol{x}\rightarrow\boldsymbol{x}_{0}$ which is now on the surface $S$, the weak singularity associated with the integral over $G$ has been eliminated because $[(\partial\phi/\partial n)-(\partial\phi/\partial n)_{0}\boldsymbol{n}_{0}\boldsymbol{\cdot}\boldsymbol{n}]\rightarrow 0$. In the integral over $\boldsymbol{\nabla}G$, there will no longer be a Cauchy principal value integral or Dirac $\delta$-function contribution as in Eq. (2). The approach in Eq. (2) that we present here will give a relationship between $\phi$ and ${\partial{\phi}}/{\partial{n}}$ as in the original problem. The numerical algorithm for solving Eq. (2) using the BIM is given by Klaseboer et al. (2009). A proof of the convergence of the integrals in Eq. (2) is given in the Appendix. Liu & Rudolphi (1999) have suggested a similar method of removing the singularities via a Taylor series expansion. Unfortunately, when their approach is implemented, it would give a relationship between $\phi$, ${\partial{\phi}}/{\partial{n}}$ and $\boldsymbol{\nabla}\phi(\boldsymbol{x}_{0})\boldsymbol{\cdot}\boldsymbol{n}$, which requires knowledge of $\boldsymbol{\nabla}\phi(\boldsymbol{x}_{0})\boldsymbol{\cdot}\boldsymbol{n}$.

The surface integral in Eq. (2) is taken over all surfaces that enclose the domain. In particular, for problems in an infinite domain outside the surface $S$, one must also take into account the ‘surface at infinity’ which will give an additional term $4\pi\phi(\boldsymbol{x}_{0})$ on the left hand-side of Eq. (2), see Klaseboer et al. (2009), Liu & Rudolphi (1991) and Liu & Rudolphi (1999).

3 Helmholtz problem — non-singular boundary integral formulation

For the solution of the Helmholtz equation $\boldsymbol{\nabla}^{2}\phi+k^{2}\phi=0$, in which $k$ is the wave number, we have the free space Green’s function $H=\cos(kr)/r$, $r=|\boldsymbol{x}-\boldsymbol{x}_{0}|$ that satisfies $\boldsymbol{\nabla}^{2}H+k^{2}H=-4\pi\delta(\boldsymbol{x}-\boldsymbol{x}_{0})$. In the same way as we obtained Eq. (5) for the potential problem, we find, for $\boldsymbol{x}_{0}$ in the domain

[TABLE]

where

[TABLE]

Upon putting $\boldsymbol{x}_{0}$ onto the surface $S$ in Eq. (8), we can eliminate all singular terms with the choice

[TABLE]

Thus the non-singular formulation of boundary integral equation for the Helmholtz problem when $\boldsymbol{x}_{0}$ is now located on the surface S, with $\chi\equiv k\boldsymbol{n}_{0}\boldsymbol{\cdot}(\boldsymbol{x}-\boldsymbol{x}_{0})$, is

[TABLE]

This is the key result in which all singularities associated with the boundary integral formulation of the Helmholtz problem have been removed. As $\boldsymbol{x}\rightarrow\boldsymbol{x}_{0}$ , the analytic structure of the integrands Eq. (3) is essentially the same as that in Eq. (2) where the integrands do not diverge. In the limit $k\rightarrow 0$, this reduces to result Eq.(2) for the potential problem.

4 Stokes problem — non-singular boundary integral formulation

The governing equations for the pressure, $p$, and velocity field, $\boldsymbol{u}$ for incompressible Stokes flow in a Newtonian fluid of dynamic viscosity $\mu$ are

[TABLE]

The stress tensor $\sigma_{ik}$ is given by

[TABLE]

where $\delta_{ik}$ is the Kronecker delta function. Using the Lorentz reciprocal theorem (Lorentz, 1907), the velocity component in the $j$-th direction, $u_{j}^{0}\equiv u_{j}(\boldsymbol{x}_{0})$, at position $\boldsymbol{x}_{0}$ in the fluid domain can be written as a boundary integral over the enclosing surface $S$ (Pozrikidis, 1992) as

[TABLE]

The $i^{th}$ component of the traction vector $\boldsymbol{f}$, is defined as $f_{i}=\mathsfbi{\sigma}_{ik}n_{k}$. Equation (14) for the Stokes problem is the analogue of (1) for the potential problem. The fundamental solutions for Stokes flow $\mathsfbi{U}_{ij}$ and $\mathsfbi{T}_{ijk}$ are given by (Pozrikidis, 1992):

[TABLE]

[TABLE]

where ${\hat{x\mkern 3.0mu}\mkern-3.0mu}{}_{i}$ etc, are the components of ${\hat{\boldsymbol{x}\mkern 3.0mu}\mkern-3.0mu}{}=\boldsymbol{x}-\boldsymbol{x}_{0}$, $r=|{\hat{\boldsymbol{x}\mkern 3.0mu}\mkern-3.0mu}{}|$ and $n_{k}$ is the $k$-th component of the unit normal of the surface pointing out of the flow domain. The functions $\mathsfbi{U}_{ij}$ and $\mathsfbi{T}_{ijk}$ diverge as $1/r$ and $1/r^{2}$, respectively, with the same behaviour as $G$ and $\boldsymbol{\nabla}G$ for the potential problem and give rise to singular behaviour in Eq. (14) when $\boldsymbol{x}_{0}$ is on the surface $S$. The traditional boundary integral formulation is obtained by putting $\boldsymbol{x}_{0}$ onto the surface $S$ in Eq. (14), and as in Eq. (2), this will give rise to Cauchy principal value integrals and Dirac $\delta$-function contributions on the left hand-side and weakly singular integrands on the right hand-side of Eq. (2). To remove such singularities, consider a zero-pressure linear velocity field

[TABLE]

where the matrix $\mathsfbi{\boldsymbol{M}}(\boldsymbol{x}_{0})$ will be chosen to cancel the arising singularities in Eq. (14) when $\boldsymbol{x}_{0}$ is on the surface $S$. The symmetric stress tensor corresponding to this linear flow field is

[TABLE]

For the velocity field $\boldsymbol{w}$ to meet the incompressibility condition: $\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{w}(\boldsymbol{x})=0$,

[TABLE]

must hold. Since $\boldsymbol{w}$ satisfies the equations for Stokes flow, the difference $(\boldsymbol{u-w})$ also satisfies Eq. (14). In component form the difference becomes ($x^{0}_{j}=j^{th}$ component of $\boldsymbol{x}_{0}$)

[TABLE]

This integral equation, with $\boldsymbol{x}_{0}$ located on the surface S, will have no singular behaviour if we choose (adopting the convention of implicit summation over repeated indices)

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

The relation between the stress tensor $\boldsymbol{\Sigma}(\boldsymbol{x}_{0})$ and the matrix $\boldsymbol{M}(\boldsymbol{x}_{0})$ in terms of the traction $\boldsymbol{f}(\boldsymbol{x}_{0})$ and the surface normal $\boldsymbol{n}(\boldsymbol{x}_{0})$ needed to ensure Eq.(20) is non-singular is in fact not unique. It is easy to verify that the expressions in Eqs. (21) to (24) obey those constraints in Eq. (19) and as a result, the integrands in Eq.(20) are not singular as $\boldsymbol{x}\rightarrow\boldsymbol{x}_{0}$ on the surface (see the proof in Appendix). Thus, Eqs. (20) to (24) form the non-singular boundary integral formulation of the Stokes problem. Analogous to the potential problem, the elements of the matrix $\boldsymbol{M}(\boldsymbol{x}_{0})$ vary for each $\boldsymbol{x}_{0}$ on the surface $S$.

5 Linear elasticity problem — non-singular boundary integral formulation

The regularisation method described thus far is quite general. The non-singular boundary integral formulation of the Stokes problem can be adapted to the linear elastic problem in solid mechanics as follows. The strain tensor $\mathsfbi{\varepsilon}$ in the elastic problem is defined in terms of the position vector field $\boldsymbol{u}$

[TABLE]

For a linear elastic material in equilibrium and in the absence of body forces, the stress tensor satisfies $\partial{{\sigma}_{ij}}/\partial{x_{j}}=0$ and is given by

[TABLE]

where the Poisson ratio $\nu$, and the shear modulus $\mu$ are related to the Young’s modulus, $E=2\mu(1+\nu)$. The displacement field $\boldsymbol{u}(\boldsymbol{x}_{0})$ at an interior point of the elastic material can be expressed in terms of integrals over the displacement field and the surface traction on the enclosing surface $S$ by the same equation as Eq. (14), except the fundamental solutions $\mathsfbi{U}_{ij}$ and $\mathsfbi{T}_{ijk}$ for the linear elastic problem are now given by Becker (1992):

[TABLE]

and

[TABLE]

With these replacements, Eqs. (20) to (24) are also the non-singular boundary integral formulation of the linear elastic problem for the displacement field $\boldsymbol{u}$ with the traction vector $\boldsymbol{f}$ defined by $f_{i}=\mathsfbi{\sigma}_{ik}n_{k}$. For an incompressible material: $\nu=1/2$, then $\mathsfbi{U}_{ij}$ and $\mathsfbi{T}_{ijk}$ for the linear elastic problem and the Stokes problem become identical.

6 Discussion and implementation

In this communication we have developed a non-singular boundary integral formulation for solving four common and related problems in hydrodynamics and solid mechanics. The common theme in the formulation is the removal of the singularities associated with the traditional boundary integral formulation by subtracting a simpler solution of a related problem with an appropriate choice of the free parameter in the solution.

The numerical implementation of our non-singular formulation for the potential problem, Eq. (2), has been described in Klaseboer et al. (2009). The Stokes problem is very similar, except that the matrix elements appear in blocks of sub-matrices. When the surface is discretised in nodes and elements, the usual Gaussian-quadrature integration procedure can be applied for all elements, including the singular ones. This will result in a system of equations relating the potential and its normal derivative through two influence matrices, which is the discretized equivalent of Eq. (2). The previously singular contributions can be found on the diagonals of the influence matrices, one corresponding to $G(\boldsymbol{x},\boldsymbol{x}_{0})$ and one to its normal derivative. All of the terms corresponding to $\boldsymbol{x}_{0}$ now correspond to those contributions and can be obtained by simple summation. Several examples for which analytical solutions exist were tested (see Klaseboer et al. (2009) for more details).

For the Stokes problem, the singularities appear in blocks of $3\times 3$ around the diagonals of the influence matrices. A procedure very similar to that followed for the potential flow can be followed to get those values by summation once more, based on Eq. (20).

Two examples are provided for the Stokes flow implementation. Both use flat three-noded linear elements in which the surface representation and the shape functions are linear. The first example is that of a sphere moving with a constant velocity $\boldsymbol{U}$. The mesh is similar to that used in Klaseboer et al. (2009). Thus the velocity vectors, $\boldsymbol{u}$, are given at all nodes as $\boldsymbol{u}=\boldsymbol{U}$. The traction, $\boldsymbol{f}$, is then calculated and, within the expected discretisation error, agrees with the analytic solution $\boldsymbol{f}=3\mu\boldsymbol{U}/(2R)$, where $R$ is the radius of the sphere.

In the second test case, the traction $\boldsymbol{f}$ is given instead and the velocity $\boldsymbol{u}$ is calculated. We use the same test case as presented in Pigeonneau & Sellier (2011), and take $\boldsymbol{f}=-\rho gz\boldsymbol{n}$ that corresponds to a spherical bubble rising under buoyancy force, where $g$ is the magnitude of gravity, and $\rho$ is the fluid density. The exact solution in cylindrical coordinates is given by Eqs. (34) and (35) in Pigeonneau & Sellier (2011): $u_{r}=U\sin(2\theta)/4$ and $u_{z}=U(1-\sin^{2}\theta/2)$, where $U=\rho gR^{2}/(3\mu)$, and $\theta$ is the angle between the unit vector in the $z$-direction and the radial direction. The results are shown in Figure 2. Even for a mesh consisting of only $252$ nodes ($500$ elements) the accuracy is within $2$%.

The present non-singular formulation therefore offers all the advantages associated with the reduction of dimension afforded by the boundary integral technique without the extra numerical effort needed to handle the singularities that arise with the traditional boundary integrals formulation. Although the size of the numerical problem remains the same, the absence of singularities means that there will be a significant reduction in coding effort which will minimise the opportunity for coding error.

We believe the present contribution is a novel advance that will have both pedagogical and practical implications.

Acknowledgements.

EK would like to thank A. Prosperetti for stimulating discussions and B. C. Khoo for maintaining his interest in practical applications of boundary element methods. DYCC is a Visiting Scientist at the IHPC 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.

\oneappendix

7 Non-singular proof

We show that the integrands in Eq. (2) for the potential flow problems and Eq. (20) for the Stokes flow problems are non-singular by analysing the analytic behaviour of the integrand in the neighbourhood of $\boldsymbol{x}_{0}$. Define a Cartesian system $(\xi,\eta,\zeta)$ with $\boldsymbol{x}_{0}=(0,0,0)$ as the origin and $\boldsymbol{n}_{0}=(0,0,1)$. In the neighbourhood of $\boldsymbol{x}_{0}$, a point $\boldsymbol{x}=(\xi,\eta,\zeta)$ that lies on the surface $S$ with a suitable choice of the local coordinates, $\xi$, $\eta$ and $\zeta$, must satisfy

[TABLE]

where higher order terms of $O(\xi^{3},\eta^{3})$ have been omitted. The constants $a_{s}$ and $b_{s}$ are related to the principal curvatures of $S$ at $\boldsymbol{x}_{0}$ and $\zeta$ is quadratic in $\xi$ and $\eta$. The unit normal vector at $\boldsymbol{x}$ is $\boldsymbol{n}=\boldsymbol{\nabla}{S}/|\boldsymbol{\nabla}{S}|=\cos{\gamma}\text{ }(a_{s}\xi,b_{s}\eta,1)$, where $\cos{\gamma}=[1+(a_{s}\xi)^{2}+(b_{s}\eta)^{2}]^{-1/2}$ is the direction cosine. The differential surface is $\text{d}S(\boldsymbol{x})=\text{d}\xi\text{d}\eta/\cos{\gamma}$. The Green function has the form: $G(\boldsymbol{x},\boldsymbol{x}_{0})=(\xi^{2}+\eta^{2}+\zeta^{2})^{-1/2}$.

With these preliminary results, the integral on the left hand-side of Eq. (2), has the following limiting form obtained by using a Taylor expansion about $\boldsymbol{x}_{0}$ in terms of the local coordinates

[TABLE]

We see that both the numerator and the denominator of the integrand are of $O(\xi^{3},\eta^{3})$, thus the integrand remains finite, as $\xi\rightarrow 0$ and $\eta\rightarrow 0$, that is $\boldsymbol{x}\rightarrow\boldsymbol{x}_{0}$. Furthermore, if the surface around $\boldsymbol{x}_{0}$ is a planar element, the constants $a_{s}$ and $b_{s}$ will be zero and the integrand vanishes. For non-planar elements, the point-wise discontinuity at $(\xi,\eta)=(0,0)$ is a set of measure zero and therefore does not contribute to the value of the integral.

Similarly, as $\boldsymbol{x}\rightarrow\boldsymbol{x}_{0}$, the integral on the right hand-side of Eq. (2), has the form

[TABLE]

We see that both the numerator and the denominator of the integrand vanish linearly with $\xi$ and $\eta$, as $\xi\rightarrow 0$ and $\eta\rightarrow 0$. Thus the integrand has no divergences. Again, if the surface around $\boldsymbol{x}_{0}$ is a planar element and the normal derivative $(\partial{\phi}/\partial{\zeta})$ is constant over that element, then the integrand vanishes. For non-planar elements, the point-wise discontinuity at $(\xi,\eta)=(0,0)$ is a set of measure zero and therefore does not contribute to the value of the integral. This has been demonstrated for quadratic elements (Klaseboer et al., 2009).

This completes the proof that Eq. (2) is non-singular. In the same way, Eq. (3) for the Helmholtz problem can also be shown to be non-singular.

Using the same local coordinate system, we can also show that Eq. (20) for the Stokes problem is not singular. First we consider the integral on the left hand-side of Eq. (20). As $\boldsymbol{x}\rightarrow\boldsymbol{x}_{0}$, it is straightforward to show using a Taylor expansion that the two terms in the integrand have the limiting form using the notation ${\hat{\boldsymbol{x}\mkern 3.0mu}\mkern-3.0mu}{}\equiv\boldsymbol{x}-\boldsymbol{x}_{0}=(\xi,\eta,\zeta)$, $\boldsymbol{f}(\boldsymbol{x})=(f_{1},f_{2},f_{3})$ and $\boldsymbol{f}(\boldsymbol{x}_{0})=(f^{0}_{1},f^{0}_{2},f^{0}_{3})$

[TABLE]

The numerator and the denominator of both terms are of $O(\xi^{5},\eta^{5})$ and, thus, they approach constant values as $\boldsymbol{x}\rightarrow\boldsymbol{x}_{0}$, and as a consequence, the integral on the left hand-side of Eq. (20) is non-singular.

Turning now to the integral on the right hand-side of Eq. (20) where in the limit $\boldsymbol{x}\rightarrow\boldsymbol{x}_{0}$, the integrand has the limiting form

[TABLE]

where $L_{i}\equiv\xi(\partial{f_{i}}/\partial{\xi})_{0}+\eta(\partial{f_{i}}/\partial{\eta})_{0}+\zeta(\partial{f_{i}}/\partial{\zeta})_{0}$. The numerator and the denominator of the first term are of $O(\xi,\eta)$ and those of the second and third terms are of $O(\xi^{3},\eta^{3})$ and thus all terms approach constant values as $\boldsymbol{x}\rightarrow\boldsymbol{x}_{0}$, and so it follows that the integral on the right hand-side of Eq. (20) is non-singular.

For the case in which the surface $S$ is a planar element for which the curvatures $a_{s}$ and $b_{s}$ are zero, $\cos{\gamma}=1$ and the traction $\boldsymbol{f}(\boldsymbol{x}_{0})$ is constant within the plane, the integrands of both integrals in Eq. (20) will vanish. For non-planar elements, the point-wise discontinuity at $(\xi,\eta)=(0,0)$ is a set of measure zero and therefore does not contribute to the value of the integral.

This completes the proof that Eq. (20) is non-singular. In the same way, the linear elasticity problem can also shown to be non-singular.