Boundary regularized integral equation formulation (BRIEF) of Stokes   flow

Q. Sun; E. Klaseboer; B. C. Khoo; D. Y. C. Chan

arXiv:1902.02886·physics.comp-ph·February 11, 2019

Boundary regularized integral equation formulation (BRIEF) of Stokes flow

Q. Sun, E. Klaseboer, B. C. Khoo, D. Y. C. Chan

PDF

TL;DR

This paper introduces a boundary regularized integral equation method for Stokes flow that eliminates singularities, leading to more accurate and easier-to-implement solutions, especially in complex boundary contact scenarios.

Contribution

It presents an exact, singularity-free boundary integral formulation for Stokes flow using auxiliary fields, simplifying implementation and enhancing numerical accuracy.

Findings

01

Eliminates mathematical singularities in boundary integral equations.

02

Provides significant coding and computational savings.

03

Maintains high accuracy even with nearly contacting boundaries.

Abstract

Single-phase Stokes flow problems with prescribed boundary conditions can be formulated in terms of a boundary regularized integral equation that is completely free of singularities that exist in the traditional formulation. The usual mathematical singularities that arise from using the fundamental solution in the conventional boundary integral method are removed by subtracting a related auxiliary flow field, $w$ , that can be constructed from one of many known fundamental solutions of the Stokes equation. This approach is exact and does not require the introduction of additional cutoff parameters. The numerical implementation of this boundary regularized integral equation formulation affords considerable savings in coding effort with improved numerical accuracy. The high accuracy of this formulation is retained even in problems where parts of the boundaries may almost be in…

Tables2

Table 1. Table 1: The auxiliary flow field w i = u i 0 + 1 μ M j k Q i j k subscript 𝑤 𝑖 superscript subscript 𝑢 𝑖 0 1 𝜇 subscript 𝑀 𝑗 𝑘 subscript 𝑄 𝑖 𝑗 𝑘 w_{i}=u_{i}^{0}+\frac{1}{\mu}{M}_{jk}{Q}_{ijk} and the corresponding pressure, 𝒫 𝒫 \mathcal{P} , constructed using fundamental solutions, Q i j k subscript 𝑄 𝑖 𝑗 𝑘 {Q}_{ijk} , of the Stokes equation that are third order tensors. The domain of applicability are for interior problems bounded by a closed surface, exterior problems in an infinite domain or both. All other symbols are defined in the text or in Fig. 1.

$Q_{i j k} = \hat{x} {}_{k}δ_{i j}$
Solution name	Pressure, $𝒫$	$x_{i}^{d}$	Domain
(I) Simple linear solution(Klaseboer, Sun, and Chan, 2012)	Constant	No	Both
$M_{j k} = f_{j}^{0} n_{k}^{0} - \frac{1}{4} f_{l}^{0} n_{l}^{0} [δ_{j k} + n_{j}^{0} n_{k}^{0}]$
(II) Stresson(Chwang and Wu, 1975; Pozrikidis, 2011)	Constant	No	Interior
$Q_{i j k} = 4 \hat{x} {}_{j}δ_{k i} - \hat{x} {}_{k}δ_{i j} - \hat{x} {}_{i}δ_{j k}$
$M_{j k} = \frac{1}{6} (f_{j}^{0} n_{k}^{0} + n_{j}^{0} f_{k}^{0}) + \frac{1}{3} f_{l}^{0} n_{l}^{0} [δ_{j k} - 4 n_{j}^{0} n_{k}^{0}]$
(III) Stokeson dipole(Chwang and Wu, 1975; Pozrikidis, 2011)	Constant	No	Interior
$Q_{i j k} = 4 \hat{x} {}_{j}δ_{k i} - \hat{x} {}_{k}δ_{i j} - \hat{x} {}_{i}δ_{j k}$
$M_{j k} = \frac{1}{3} (n_{j}^{0} f_{k}^{0}) + \frac{1}{48} f_{l}^{0} n_{l}^{0} [δ_{j k} - 19 n_{j}^{0} n_{k}^{0}]$
(IV) Source potential(Pozrikidis, 2011)	Constant	Yes	Both
$Q_{i j k} = (\frac{\hat{x}_{k}^{d}}{r_{d}^{3}} - \frac{\hat{a}_{k}^{d}}{a^{3}}) δ_{i j}$
$M_{j k} = - \frac{a^{3}}{3 b} (f_{j}^{0} \hat{a} {}_{k}^{d}- \hat{a} {}_{j}^{d}f_{k}^{0}) + \frac{a^{3}}{2} [\frac{a^{2} (f_{l}^{0} n_{l}^{0}) - 2 b (f_{l}^{0} \hat{a} {}_{l}^{d})}{a^{2} + b^{2}}] [\frac{2}{3 b} (n_{j}^{0} \hat{a} {}_{k}^{d}- \hat{a} {}_{j}^{d}n_{k}^{0}) + δ_{j k}]$
(V) Stresslet(Chwang and Wu, 1975; Pozrikidis, 2011)	$2 M_{j k} (3 \frac{\hat{x} {}_{j}^{d}{\hat{x}}_{k}^{d}}{r_{d}^{5}} - \frac{δ_{j k}}{r_{d}^{3}})$	Yes	Both
$Q_{i j k} = (3 \frac{\hat{x} {}_{i}^{d}{\hat{x}} {}_{j}^{d}{\hat{x}}_{k}^{d}}{r_{d}^{5}} - \frac{\hat{x} {}_{i}^{d}δ_{j k}}{r_{d}^{3}}) - (3 \frac{\hat{a} {}_{i}^{d}{\hat{a}} {}_{j}^{d}{\hat{a}}_{k}^{d}}{a^{5}} - \frac{\hat{a} {}_{i}^{d}δ_{j k}}{a^{3}})$
$M_{j k} = \frac{a}{12 b^{2}} \hat{a} {}_{j}^{d}{\hat{a}} {}_{k}^{d}{[- 6 b (f_{l}^{0} \hat{a} {}_{l}^{d}) + a^{2} (f_{l}^{0} n_{l}^{0})]} + \frac{a^{3}}{6 b} (\hat{a} {}_{j}^{d}f_{k}^{0} + f_{j}^{0} \hat{a} {}_{k}^{d})$
(VI) Stokes doublet(Chwang and Wu, 1975; Pozrikidis, 2011)	$2 M_{j k} (3 \frac{\hat{x} {}_{j}^{d}{\hat{x}}_{k}^{d}}{r_{d}^{5}} - \frac{δ_{j k}}{r_{d}^{3}})$	Yes	Both
$Q_{i j k} = [(3 \frac{\hat{x} {}_{i}^{d}{\hat{x}} {}_{j}^{d}{\hat{x}}_{k}^{d}}{r_{d}^{5}} - \frac{\hat{x} {}_{i}^{d}δ_{j k}}{r_{d}^{3}}) - (3 \frac{\hat{a} {}_{i}^{d}{\hat{a}} {}_{j}^{d}{\hat{a}}_{k}^{d}}{a^{5}} - \frac{\hat{a} {}_{i}^{d}δ_{j k}}{a^{3}})] + [(\frac{\hat{x} {}_{k}^{d}δ_{i j} - \hat{x} {}_{j}^{d}δ_{i k}}{r_{d}^{3}}) - (\frac{\hat{a} {}_{k}^{d}δ_{i j} - \hat{a} {}_{j}^{d}δ_{i k}}{a^{3}})]$
$M_{j k} = \frac{a^{5}}{6 b^{2}} (n_{j}^{0} f_{k}^{0}) + \frac{a^{3}}{b} \hat{a} {}_{j}^{d}{\hat{a}} {}_{k}^{d}{(\frac{f_{l}^{0} n_{l}^{0}}{6 b} - \frac{5}{12 a^{2}} f_{l}^{0} \hat{a} {}_{l}^{d})}$

Table 2. Table 2: The auxiliary flow field w i = u i 0 + 1 μ D j S i j subscript 𝑤 𝑖 superscript subscript 𝑢 𝑖 0 1 𝜇 subscript 𝐷 𝑗 subscript 𝑆 𝑖 𝑗 w_{i}=u_{i}^{0}+\frac{1}{\mu}D_{j}{S}_{ij} and the corresponding pressure, 𝒫 𝒫 \mathcal{P} , constructed using fundamental solutions, S i j subscript 𝑆 𝑖 𝑗 {S}_{ij} , of the Stokes equation that are second order tensors. All other symbols are defined in the text or in Fig. 1.

$S_{i j} = (2 r_{d}^{2} δ_{i j} - \hat{x} {}_{i}^{d}{\hat{x}} {}_{j}^{d}) - (2 a^{2} δ_{i j} - \hat{a} {}_{i}^{d}{\hat{a}} {}_{j}^{d})$
Solution name	Pressure, $𝒫$	$x_{i}^{d}$	Domain
(I) Stokeson(Chwang and Wu, 1975; Pozrikidis, 2011)	$10 D_{j} \hat{x}_{j}^{d}$	Yes	Interior
$D_{j} = [\frac{2 f_{l}^{0} \hat{a}_{l}^{d}}{b (9 b^{2} - 6 a^{2})} - \frac{f_{l}^{0} n_{l}^{0}}{6 b^{2} - 4 a^{2}}] [- 2 b n_{j}^{0} + \hat{a} {}_{j}^{d}] + \frac{1}{3 b} [f_{j}^{0} - (f_{l}^{0} n_{l}^{0}) n_{j}^{0}]$
(II) Source potential doublet(Chwang and Wu, 1975; Pozrikidis, 2011)	Constant	Yes	Both
$S_{i j} = [\frac{δ_{i j}}{r_{d}^{3}} - \frac{3 \hat{x} {}_{i}^{d}{\hat{x}}_{j}^{d}}{r_{d}^{5}}] - [\frac{δ_{i j}}{a^{3}} - \frac{3 \hat{a} {}_{i}^{d}{\hat{a}}_{j}^{d}}{a^{5}}]$
$D_{j} = - \frac{a^{5} f_{j}^{0}}{6 b} + \frac{a^{5} n_{j}^{0}}{6 b} [\frac{a^{2} (f_{l}^{0} n_{l}^{0}) - 2 b (f_{l}^{0} \hat{a} {}_{l}^{d})}{a^{2} + b^{2}}] + \frac{\hat{a}_{j}^{d}}{6 b} [\frac{- 2 a^{5} b (f_{l}^{0} n_{l}^{0}) + (a^{5} + 5 a^{3} b^{2}) (f_{l}^{0} \hat{a} {}_{l}^{d})}{a^{2} + b^{2}}]$

Equations92

- \frac{\partial p}{\partial x _{i}} + μ \frac{\partial ^{2} u _{i}}{\partial x _{k} \partial x _{k}} = 0

- \frac{\partial p}{\partial x _{i}} + μ \frac{\partial ^{2} u _{i}}{\partial x _{k} \partial x _{k}} = 0

\frac{\partial u _{i}}{\partial x _{i}} = 0.

\frac{\partial u _{i}}{\partial x _{i}} = 0.

σ_{ik} = - p δ_{ik} + μ [\frac{\partial u _{i}}{\partial x _{k}} + \frac{\partial u _{k}}{\partial x _{i}}],

σ_{ik} = - p δ_{ik} + μ [\frac{\partial u _{i}}{\partial x _{k}} + \frac{\partial u _{k}}{\partial x _{i}}],

c_{0} u_{j}^{0} + \int_{S} u_{i} T_{ij k} n_{k} d S (x) = \frac{1}{μ} \int_{S} σ_{ik} n_{k} U_{ij} d S (x) \equiv \frac{1}{μ} \int_{S} f_{i} U_{ij} d S (x) .

c_{0} u_{j}^{0} + \int_{S} u_{i} T_{ij k} n_{k} d S (x) = \frac{1}{μ} \int_{S} σ_{ik} n_{k} U_{ij} d S (x) \equiv \frac{1}{μ} \int_{S} f_{i} U_{ij} d S (x) .

U_{ij} (x, x_{0}) = \frac{δ _{ij}}{r} + \frac{x ^ _{i} x ^ _{j}}{r ^{3}},

U_{ij} (x, x_{0}) = \frac{δ _{ij}}{r} + \frac{x ^ _{i} x ^ _{j}}{r ^{3}},

T_{ij k} (x, x_{0}) = - 6 \frac{x ^ _{i} x ^ _{j} x ^ _{k}}{r ^{5}},

T_{ij k} (x, x_{0}) = - 6 \frac{x ^ _{i} x ^ _{j} x ^ _{k}}{r ^{5}},

- \frac{\partial P}{\partial x _{i}} + μ \frac{\partial ^{2} w _{i}}{\partial x _{k} \partial x _{k}} = 0

- \frac{\partial P}{\partial x _{i}} + μ \frac{\partial ^{2} w _{i}}{\partial x _{k} \partial x _{k}} = 0

\frac{\partial w _{i}}{\partial x _{i}} = 0

\frac{\partial w _{i}}{\partial x _{i}} = 0

Σ_{ik} = - P δ_{ik} + μ [\frac{\partial w _{i}}{\partial x _{k}} + \frac{\partial w _{k}}{\partial x _{i}}] .

Σ_{ik} = - P δ_{ik} + μ [\frac{\partial w _{i}}{\partial x _{k}} + \frac{\partial w _{k}}{\partial x _{i}}] .

w_{i} (x_{0}) \equiv w_{i}^{0} = u_{i}^{0} \equiv u_{i} (x_{0}),

w_{i} (x_{0}) \equiv w_{i}^{0} = u_{i}^{0} \equiv u_{i} (x_{0}),

F_{i} (x_{0}) \equiv F_{i}^{0} = f_{i}^{0} \equiv f_{i} (x_{0}),

c_{0} w_{j}^{0} + \int_{S} w_{i} T_{ij k} n_{k} d S = \frac{1}{μ} \int_{S} F_{i} U_{ij} d S .

c_{0} w_{j}^{0} + \int_{S} w_{i} T_{ij k} n_{k} d S = \frac{1}{μ} \int_{S} F_{i} U_{ij} d S .

\int_{S} (u_{i} - w_{i}) T_{ij k} n_{k} d S = \frac{1}{μ} \int_{S} (f_{i} - F_{i}) U_{ij} d S .

\int_{S} (u_{i} - w_{i}) T_{ij k} n_{k} d S = \frac{1}{μ} \int_{S} (f_{i} - F_{i}) U_{ij} d S .

w = u_{0} + \frac{1}{μ} M \cdot (x - x_{0}) \equiv u_{0} + \frac{1}{μ} M \cdot \hat{x},

w = u_{0} + \frac{1}{μ} M \cdot (x - x_{0}) \equiv u_{0} + \frac{1}{μ} M \cdot \hat{x},

w_{i} = u_{i}^{0} + \frac{1}{μ} M_{ik} \overset{x}{^}_{k} = u_{i}^{0} + \frac{1}{μ} M_{j k} \overset{x}{^}_{k} δ_{ij} .

w_{i} = u_{i}^{0} + \frac{1}{μ} M_{ik} \overset{x}{^}_{k} = u_{i}^{0} + \frac{1}{μ} M_{j k} \overset{x}{^}_{k} δ_{ij} .

M_{l l} \equiv Tr (M) = 0.

M_{l l} \equiv Tr (M) = 0.

Σ_{i l} = μ [\frac{\partial w _{i}}{\partial x _{l}} + \frac{\partial w _{l}}{\partial x _{i}}] = M_{j k} δ_{k l} δ_{ij} + M_{j k} δ_{k i} δ_{l j} = M_{i l} + M_{l i},

Σ_{i l} = μ [\frac{\partial w _{i}}{\partial x _{l}} + \frac{\partial w _{l}}{\partial x _{i}}] = M_{j k} δ_{k l} δ_{ij} + M_{j k} δ_{k i} δ_{l j} = M_{i l} + M_{l i},

F_{i} \equiv Σ_{i l} n_{l} = (M_{i l} + M_{l i}) n_{l} .

F_{i}^{0} \equiv [M_{i l} + M_{l i}] n_{l}^{0} = f_{i}^{0} .

F_{i}^{0} \equiv [M_{i l} + M_{l i}] n_{l}^{0} = f_{i}^{0} .

M_{j k} = f_{j}^{0} n_{k}^{0} + f_{l}^{0} n_{l}^{0} [c_{1} δ_{j k} + c_{2} n_{j}^{0} n_{k}^{0}],

M_{j k} = f_{j}^{0} n_{k}^{0} + f_{l}^{0} n_{l}^{0} [c_{1} δ_{j k} + c_{2} n_{j}^{0} n_{k}^{0}],

2 c_{1} + 2 c_{2} + 1

2 c_{1} + 2 c_{2} + 1

1 + 3 c_{1} + c_{2}

M_{j k} = f_{j}^{0} n_{k}^{0} - \frac{1}{4} f_{l}^{0} n_{l}^{0} [δ_{j k} + n_{j}^{0} n_{k}^{0}],

M_{j k} = f_{j}^{0} n_{k}^{0} - \frac{1}{4} f_{l}^{0} n_{l}^{0} [δ_{j k} + n_{j}^{0} n_{k}^{0}],

w_{i} = u_{i}^{0} + \frac{1}{μ} M_{j k} Q_{ij k},

w_{i} = u_{i}^{0} + \frac{1}{μ} M_{j k} Q_{ij k},

w_{i} = u_{i}^{0} + \frac{1}{μ} M_{j k} [(3 \frac{x ^ _{i}^{d} x ^ _{j}^{d} x ^ _{k}^{d}}{r _{d}^{5}} - \frac{x ^ _{i}^{d} δ _{j k}}{r _{d}^{3}}) - (3 \frac{a ^ _{i}^{d} a ^ _{j}^{d} a ^ _{k}^{d}}{a ^{5}} - \frac{a ^ _{i}^{d} δ _{j k}}{a ^{3}})],

w_{i} = u_{i}^{0} + \frac{1}{μ} M_{j k} [(3 \frac{x ^ _{i}^{d} x ^ _{j}^{d} x ^ _{k}^{d}}{r _{d}^{5}} - \frac{x ^ _{i}^{d} δ _{j k}}{r _{d}^{3}}) - (3 \frac{a ^ _{i}^{d} a ^ _{j}^{d} a ^ _{k}^{d}}{a ^{5}} - \frac{a ^ _{i}^{d} δ _{j k}}{a ^{3}})],

P = 2 M_{j k} (3 \frac{x ^ _{j}^{d} x ^ _{k}^{d}}{r _{d}^{5}} - \frac{δ _{j k}}{r _{d}^{3}}) .

P = 2 M_{j k} (3 \frac{x ^ _{j}^{d} x ^ _{k}^{d}}{r _{d}^{5}} - \frac{δ _{j k}}{r _{d}^{3}}) .

F_{i} = \frac{6}{r _{d}^{5}} [M_{l l} (\overset{x}{^}_{l}^{d} n_{l}) \overset{x}{^}_{i}^{d} + (\overset{x}{^}_{j}^{d} M_{j k} n_{k}) \overset{x}{^}_{i}^{d} + \overset{x}{^}_{j}^{d} M_{j i} (\overset{x}{^}_{l}^{d} n_{l}) - \frac{5}{r _{d}^{2}} (\overset{x}{^}_{j}^{d} M_{j k} \overset{x}{^}_{k}^{d}) (\overset{x}{^}_{l}^{d} n_{l}) \overset{x}{^}_{i}^{d}] .

F_{i} = \frac{6}{r _{d}^{5}} [M_{l l} (\overset{x}{^}_{l}^{d} n_{l}) \overset{x}{^}_{i}^{d} + (\overset{x}{^}_{j}^{d} M_{j k} n_{k}) \overset{x}{^}_{i}^{d} + \overset{x}{^}_{j}^{d} M_{j i} (\overset{x}{^}_{l}^{d} n_{l}) - \frac{5}{r _{d}^{2}} (\overset{x}{^}_{j}^{d} M_{j k} \overset{x}{^}_{k}^{d}) (\overset{x}{^}_{l}^{d} n_{l}) \overset{x}{^}_{i}^{d}] .

f_{i}^{0} = \frac{6}{a ^{5}}

f_{i}^{0} = \frac{6}{a ^{5}}

M_{j k} = \frac{a}{12 b ^{2}} \overset{a}{^}_{j}^{d} \overset{a}{^}_{k}^{d} [- 6 b (f_{l}^{0} \overset{a}{^}_{l}^{d}) + a^{2} (f_{l}^{0} n_{l}^{0})] + \frac{a ^{3}}{6 b} (\overset{a}{^}_{j}^{d} f_{k}^{0} + f_{j}^{0} \overset{a}{^}_{k}^{d}),

M_{j k} = \frac{a}{12 b ^{2}} \overset{a}{^}_{j}^{d} \overset{a}{^}_{k}^{d} [- 6 b (f_{l}^{0} \overset{a}{^}_{l}^{d}) + a^{2} (f_{l}^{0} n_{l}^{0})] + \frac{a ^{3}}{6 b} (\overset{a}{^}_{j}^{d} f_{k}^{0} + f_{j}^{0} \overset{a}{^}_{k}^{d}),

w_{i} = u_{i}^{0} + \frac{1}{μ} D_{j} S_{ij}

w_{i} = u_{i}^{0} + \frac{1}{μ} D_{j} S_{ij}

w_{i} = u_{i}^{0} + \frac{1}{μ} D_{j} {[\frac{δ _{ij}}{r _{d}^{3}} - \frac{3 x ^ _{i}^{d} x ^ _{j}^{d}}{r _{d}^{5}}] - [\frac{δ _{ij}}{a ^{3}} - \frac{3 a ^ _{i}^{d} a ^ _{j}^{d}}{a ^{5}}]}

w_{i} = u_{i}^{0} + \frac{1}{μ} D_{j} {[\frac{δ _{ij}}{r _{d}^{3}} - \frac{3 x ^ _{i}^{d} x ^ _{j}^{d}}{r _{d}^{5}}] - [\frac{δ _{ij}}{a ^{3}} - \frac{3 a ^ _{i}^{d} a ^ _{j}^{d}}{a ^{5}}]}

F_{i} = - \frac{6}{r _{d}^{5}} [D_{i} (\overset{x}{^}_{l}^{d} n_{l}) + \overset{x}{^}_{i}^{d} (n_{l} D_{l}) + n_{i} (D_{l} \overset{x}{^}_{l}^{d}) - \frac{5}{r _{d}^{2}} \overset{x}{^}_{i}^{d} (D_{l} \overset{x}{^}_{l}^{d}) (\overset{x}{^}_{k}^{d} n_{k})] .

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Boundary regularized integral equation formulation (BRIEF) of Stokes flow

Q. Sun

Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, 119260, Singapore

E. Klaseboer

Institute of High Performance Computing, 1 Fusionopolis Way, 138632, Singapore

B. C. Khoo

Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, 119260, Singapore

D. Y. C. Chan

[email protected].

Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, 119260, Singapore

Institute of High Performance Computing, 1 Fusionopolis Way, 138632, Singapore

Department of Mathematics and Statistics, The University of Melbourne, Parkville 3010 VIC Australia

Department of Chemistry and Biotechnology, Swinburne University of Technology, Hawthorn 3122 VIC Australia

Abstract

Single-phase Stokes flow problems with prescribed boundary conditions can be formulated in terms of a boundary regularized integral equation that is completely free of singularities that exist in the traditional formulation. The usual mathematical singularities that arise from using the fundamental solution in the conventional boundary integral method are removed by subtracting a related auxiliary flow field, $\bm{w}$ , that can be constructed from one of many known fundamental solutions of the Stokes equation. This approach is exact and does not require the introduction of additional cutoff parameters. The numerical implementation of this boundary regularized integral equation formulation affords considerable savings in coding effort with improved numerical accuracy. The high accuracy of this formulation is retained even in problems where parts of the boundaries may almost be in contact.

††preprint: Phys. Fluids

I Introduction

A common feature in the modelling of microfluidic design,Bazhlekov, Anderson, and Meijer (2004) the motion of particles with complex mixed stick-slip boundary conditions,Sun et al. (2013) the motility of biological cells Li and Spagnolie (2014) and the dynamics of cell-cell and cell-substrate interactions Bouffanais and Yue (2010); Subramaniam, Gee, and King (2013) is the need to determine low-Reynolds-number Stokes flow in domains with boundaries having arbitrary shapes. In many applications, such boundaries may also deform as the quasi-static flow progresses, for instance, in modelling mobile deformable droplets or biological cells.Chan, Klaseboer, and Manica (2011) Furthermore, such problems may require precision in resolution over very different characteristic length scales, for example, from micrometers for cell or droplet dimensions to nanometers for the deformation and spacing between interacting surfaces.Chan, Klaseboer, and Manica (2011)

It is evident that modelling of the above problems using numerical methods based on the discretization of the 3D spatial domain will encounter a number of challenges. These include the need to insert and delete grid points when the boundaries deform and multi-grid methods might be needed to give the desired resolution and precision over different length scales. When high accuracy in locating boundaries is required, the 3D domain based discretisation algorithms can become impractical.

On the other hand, an approach based on the boundary integral formulation has a number of advantages for such problems. The most obvious is the reduction in dimension from 3D to 2D as the focus is entirely on the boundaries. Thus precision tracking of their deformation or motion becomes easier. Although the boundary integral formulation will give rise to a full coefficient matrix, the advent of the fast multipole method Liu (2009) has reduced the computational cost to a very competitive $O(N\log N)$ level.

However, the conventional boundary integral formulation of Stokes flow problems involves singularities in the kernels that originate from the use of fundamental solutions of the Stokes equation.Becker (1992) Although such singularities are integrable, they do require careful analysis and additional coding effort in numerical implementation.Telles (1987); Gaul, Kögl, and Wagner (2003) Furthermore, when two boundaries are close together, the singular behaviour on one surface can adversely affect the integral taken over the other nearby surface even in cases in which the field quantities are expected to be bounded. Also, the singular behavior in the flow domain near boundaries is often more difficult to deal with than the singularities on the boundaries.

Recently we have developed a reformulation of the boundary integral equations for the potential problem,Sun et al. (2014) the Helmholtz equation, the Stokes and linear elasticity problems in which the traditional singularities can be eliminated analytically.Klaseboer, Sun, and Chan (2012) This means that none of the above mentioned issues associated with the integrable singularities in the traditional boundary integral approach will arise and the integrals can be evaluated using any convenient quadrature method. Furthermore, the term involving the solid angle in the conventional boundary integral equation has also been eliminated, thus avoiding the need to calculate the solid angle at each node that is a complex function of local geometry.Mantic̆ (1993)

In this paper, we show that corresponding to each of the many known fundamental solutions of the single-phase Stokes equation, a different regularized boundary integral equation for the Stokes problem can be derived. The regularization process results in numerically robust equations and this affords a choice of different numerically equivalent approaches for different applications.

In the next section, we give the general formulation of the boundary regularized integral equation formulation (BRIEF) for Stokes flow that involves finding an auxiliary flow field $\bm{w}$ that is used to remove all singularities analytically. In Sec. III, we give a new derivation of an earlier result using a linear flow fieldKlaseboer, Sun, and Chan (2012) that is a special case of a third order tensor fundamental solution to serve as a template on how to derive boundary regularized integral equations. In Sec. IV, we show how to construct $\bm{w}$ for other fundamental solutions that are also third order tensors and in Sec. V, we consider the construction of $\bm{w}$ from fundamental solutions that are second order tensors. For easy reference, all key results are summarized in Tables 1 and 2. Numerical examples are given in Sec. VI to demonstrate the advantages of the BRIEF of Stokes flow. The structure of the algebraic system that arises from our approach is discussed in the Appendix.

II Boundary regularized integral equation formulation (BRIEF) of Stokes flow

The Stokes equation for low-Reynolds-number flow in a Newtonian fluid with dynamic shear viscosity, $\mu$ : $-\nabla p+\mu\nabla^{2}\bm{u}=0$ can be written in Cartesian tensor notation with the summation convention of repeated indices as

[TABLE]

and the incompressibility condition: $\nabla\cdot\bm{u}=0$ is

[TABLE]

At the field point, $\bm{x}$ (with components $x_{i}$ ) where the pressure is $p(\bm{x})$ and the velocity field is $\bm{u}(\bm{x})$ (with components $u_{i}$ ), the stress tensor $\sigma_{ik}$ is given by

[TABLE]

in which $\delta_{ik}$ is the Kronecker delta function.

In the classic boundary integral formulation, the Lorentz reciprocal theorem (Lorentz, 1907) is used to give the following integral equation evaluated on the bounding surface(s), $S$ , of the flow domain relating the velocity and the stress tensor on the boundary (Pozrikidis, 1992)

[TABLE]

In Eq. (4), $c_{0}$ is related to the solid angle at $\bm{x}_{0}$ that is located on the surface $S$ . The component of the velocity at $\bm{x}_{0}$ in the $j$ -th direction is $u_{j}^{0}\equiv u_{j}(\bm{x}_{0})$ , the component of the velocity at $\bm{x}$ in the $i$ -th direction is $u_{i}\equiv u_{i}(\bm{x})$ , and $n_{k}$ is the $k$ -th component of the unit normal vector at position $\bm{x}$ on the surface, pointing out of the fluid domain. The $i$ -th component of the traction vector $\bm{f}$ , is defined as $f_{i}\equiv{\sigma}_{ik}n_{k}$ . The kernels ${U}_{ij}$ and ${T}_{ijk}$ are the fundamental solutions for the 3D Stokes equation under a point force $\bm{g}$ : $-\nabla p+\mu\nabla^{2}\bm{u}+\bm{g}\delta(\bm{x}-\bm{x}_{0})=\bm{0}$ , (Pozrikidis, 1992)

[TABLE]

where ${\hat{x\mkern 3.0mu}\mkern-3.0mu}{}_{i}$ is the $i$ -th component of ${\hat{\bm{x}\mkern 3.0mu}\mkern-3.0mu}{}\equiv\bm{x}-\bm{x}_{0}$ , $r\equiv|{\hat{\bm{x}\mkern 3.0mu}\mkern-3.0mu}{}|\equiv|\bm{x}-\bm{x}_{0}|$ (see Fig. 1), with corresponding velocity and traction fields $u_{i}=(1/8\pi\mu)U_{ij}g_{j}$ and $\sigma_{ik}=(1/8\pi)T_{ijk}g_{j}$ . The kernels ${U}_{ij}$ and ${T}_{ijk}$ diverge as $1/r$ and $1/r^{2}$ respectively as $\bm{x}\rightarrow\bm{x}_{0}$ , but the integrals over these singularities in Eq. (4) are finite even though their numerical evaluation requires careful treatment.Gaul, Kögl, and Wagner (2003) One method to deal with the singularities is to use a nonsingular contour-integral representation of the surface integrals, but the resulting computational performance is about an order of magnitude slower than the conventional approach.Bazhlekov, Anderson, and Meijer (2004) Another way to circumvent the divergence as $\bm{x}\rightarrow\bm{x}_{0}$ is to replace the Dirac $\delta$ -function by spreading the applied force over a small ball of radius $\epsilon$ centered at $\bm{x}_{0}$ .Cortez (2001) The new length scale, $\epsilon$ has to be chosen carefully to give convergence without affecting the fidelity of solution since the fundamental solutions will also be modified. A “near-singularity” subtraction method has also been proposed to handle the singularities in Eq. (4), but the method cannot completely eliminate the unbounded behavior of the double-layer integrand.Zinchenko, Rother, and Davis (1997)

In earlier work,Klaseboer, Sun, and Chan (2012) we have shown that all the singularities associated with the kernels can be removed analytically, thus obviating the need to alter the nature of the fundamental solutions or to develop special integration algorithms to deal with the singularities. The way to achieve this is to construct an auxiliary known flow field, $\bm{w}(\bm{x})$ (with components $w_{i}$ ), for a given $\bm{x}_{0}$ in Eq. (4), that also satisfies the governing equations for the Stokes flow in Eqs. (1) and (2). Thus the field $\bm{w}$ satisfies

[TABLE]

and the compressibility condition

[TABLE]

in which $\mathcal{P}$ is the pressure corresponding to flow field $w_{i}$ . The associated stress tensor ${\Sigma}_{ik}$ of this flow field is

[TABLE]

We further require $w_{i}$ to obey the following two constraints:

[TABLE]

where $\mathcal{F}_{i}\equiv{\Sigma}_{ik}n_{k}$ is the traction vector of the flow field $w_{i}$ .

Since $w_{i}$ satisfies the governing equations for incompressible Stokes flow, Eqs. (7) and (8), it also satisfies the boundary integral equation similar to Eq. (4)

[TABLE]

Subtracting Eq. (11) from Eq. (4), we get

[TABLE]

This is the key general result of the boundary regularized integral equation formulation (BRIEF) for Stokes flow that relates the traction, $f_{i}$ , and the surface velocity, $u_{i}$ . The major advantages are that the term containing the solid angle $c_{0}$ is no longer present and that the integrands in Eq. (12) are regular over the entire integration surface $S$ .Klaseboer, Sun, and Chan (2012) As a consequence, the surface integrals can be evaluated using any convenient quadrature method. The remaining task is to construct the flow field $w_{i}$ that will satisfy the two conditions given by Eq. (10). Then Eq. (12) can be solved when either the traction, $f_{i}$ , or the surface velocity, $u_{i}$ , or a relation between $f_{i}$ and $u_{i}$ is specified from the prescribed boundary conditions of the problem.

III BRIEF of Stokes flow using a linear flow field

The simplest flow field $w_{i}$ that can be constructed to satisfy Eq. (10) is the linear solution

[TABLE]

with an appropriate choice for the constant matrix ${M}$ . The explicit form of ${M}$ for this case has been given earlier without proof.Klaseboer, Sun, and Chan (2012) Here we give a general derivation of this result that will provide guidance on how to construct other forms of $w_{i}$ using different fundamental solutions of the Stokes equation.

The linear flow field in Eq. (13) can be written as that due to a third order tensor

[TABLE]

We can see immediately that $w_{i}$ now satisfies Eq. (7) with $\mathcal{P}=0$ , and Eq. (8) gives

[TABLE]

The stress tensor, ${\Sigma}_{il}$ , and traction, $\mathcal{F}_{i}$ , that correspond to $w_{i}$ are

[TABLE]

Now ${M}$ must be chosen to satisfy: $\text{Tr}({M})=0$ , and the traction, $\mathcal{F}_{i}$ , must, according to Eq. (10b), be equal to $f^{0}_{i}$ at $\bm{x}=\bm{x}_{0}$ , that is,

[TABLE]

Noting that the only vectors at our disposal are $\bm{f}_{0}$ and $\bm{n}_{0}$ , then owing to the linear nature of the problem, ${M}_{jk}$ can only be a linear combination of the tensors: ${\delta}_{jk}$ , $f_{j}^{0}n_{k}^{0}$ , $n_{j}^{0}f_{k}^{0}$ and $n_{j}^{0}n_{k}^{0}$ . One can see immediately that ${M}_{jk}$ should contain $f^{0}_{j}n^{0}_{k}$ , since $M_{jk}{\delta}_{kl}{\delta}_{ij}n^{0}_{l}=(f^{0}_{j}n^{0}_{k}){\delta}_{kl}{\delta}_{ij}n^{0}_{l}=f_{i}^{0}$ . However, terms with $f^{0}_{l}n^{0}_{l}$ will also appear due to the term $M_{jk}{\delta}_{ki}{\delta}_{lj}$ that gives $M_{jk}{\delta}_{ki}{\delta}_{lj}n^{0}_{l}=(f^{0}_{j}n^{0}_{k}){\delta}_{ki}{\delta}_{lj}n^{0}_{l}=f^{0}_{j}n^{0}_{i}n^{0}_{j}=(f^{0}_{l}n^{0}_{l})n^{0}_{i}$ . Thus the most general form of ${M}_{jk}$ is

[TABLE]

where $c_{1}$ and $c_{2}$ are constants and $c_{1}{\delta}_{jk}+c_{2}n^{0}_{j}n^{0}_{k}$ is the most general tensor, independent of $\bm{f}_{0}$ since it already appears in $f^{0}_{l}n^{0}_{l}$ , that can be constructed from $n^{0}_{j}n^{0}_{k}$ (see Landau and Lifshitz (1966)). Combining Eqs. (15), (17) and (18), we have

[TABLE]

with solution $c_{1}=c_{2}=-\frac{1}{4}$ . Thus the desired form for ${M}_{jk}$ is

[TABLE]

a result that was given earlier.Klaseboer, Sun, and Chan (2012) The linear solution for $\bm{w}$ can be used for both interior problems or exterior problems in an infinite domain, though not for semi-infinite domains.

IV BRIEF of Stokes flow using a third order tensor fundamental solution

We have seen in the preceding section that the linear flow field $w_{i}$ considered in Eq. (14) can be written in a more general form as a third order tensor

[TABLE]

in which ${M}_{jk}$ is a constant matrix and ${Q}_{ijk}$ is a third order tensor function of $\bm{x}$ that can be one of the fundamental solutions of Stokes flow.Blake and Chwang (1974); Chwang and Wu (1975); Pozrikidis (2011) The various known forms of ${Q}_{ijk}$ together with the corresponding pressure field $\mathcal{P}$ and coefficient matrix ${M}_{jk}$ are summarised in Table 1.

We now show how the constant matrix ${M}_{jk}$ can be determined by ensuring that the conditions in Eq. (10) are satisfied. Consider as an example, the Stokes stresslet fundamental solution (Chwang and Wu, 1975; Pozrikidis, 2011) given by item (V) in Table 1 for which $w_{i}$ is

[TABLE]

with the corresponding pressure field

[TABLE]

In Eqs. (22) and (23), $\bm{x}_{d}\equiv x_{i}^{d}$ is the source position of the stresslet that is located outside the flow domain (see Fig. 1), ${\hat{\bm{x}\mkern 3.0mu}\mkern-3.0mu}{}_{d}\equiv{\hat{x\mkern 3.0mu}\mkern-3.0mu}{}^{d}_{i}=x_{i}-x_{i}^{d}$ , $r_{d}=|{\hat{\bm{x}\mkern 3.0mu}\mkern-3.0mu}{}_{d}|$ , ${\hat{\bm{a}\mkern 3.0mu}\mkern-3.0mu}{}_{d}\equiv{\hat{a\mkern 3.0mu}\mkern-3.0mu}{}^{d}_{i}=x^{0}_{i}-x_{i}^{d}$ , $a=|{\hat{\bm{a}\mkern 3.0mu}\mkern-3.0mu}{}_{d}|$ . Obviously, the flow field given in Eqs. (22) and (23) satisfies Eqs. (7) and (8) as well as the first requirement of Eq. (10a): $w^{0}_{i}=u^{0}_{i}$ , as $\bm{x}\rightarrow\bm{x}_{0}$ . Noting that a stresslet is symmetric, the matrix ${M}_{jk}$ in Eq. (22) must also be symmetric. Thus, we can write the traction that corresponds to $w_{i}$ as

[TABLE]

The condition in Eq. (10b) implies

[TABLE]

There are a few different choices for ${M}_{jk}$ that can satisfy the above constraint. The vectors at our disposal are $\bm{f}_{0}$ , $\bm{n}_{0}$ and ${\hat{\bm{a}\mkern 3.0mu}\mkern-3.0mu}{}_{d}$ , so to preserve the linear nature of ${M}_{jk}$ with respect to $f_{i}^{0}$ , it should be a linear combination of ${\delta}_{jk}$ , $f_{j}^{0}{\hat{a\mkern 3.0mu}\mkern-3.0mu}{}^{d}_{k}$ , ${\hat{a\mkern 3.0mu}\mkern-3.0mu}{}^{d}_{j}f_{k}^{0}$ , $f_{j}^{0}n_{k}^{0}$ , $n_{j}^{0}f_{k}^{0}$ , $n_{j}^{0}{\hat{a\mkern 3.0mu}\mkern-3.0mu}{}^{d}_{k}$ , ${\hat{a\mkern 3.0mu}\mkern-3.0mu}{}^{d}_{j}n_{k}^{0}$ , ${\hat{a\mkern 3.0mu}\mkern-3.0mu}{}^{d}_{j}{\hat{a\mkern 3.0mu}\mkern-3.0mu}{}^{d}_{k}$ and $n_{j}^{0}n_{k}^{0}$ . One possible option for such a linear combination is

[TABLE]

in which $b={\hat{a\mkern 3.0mu}\mkern-3.0mu}{}^{d}_{l}n_{l}^{0}\neq 0$ . This is the result given in item (V) in Table 1.

Different formulations of the regular form of the boundary integral equation (12) using other fundamental solutions, ${Q}_{ijk}$ , of the Stokes equation that are third order tensors together with the corresponding pressure field, $\mathcal{P}$ , and coefficient matrix, ${M}_{jk}$ , are given in Table 1. They can be readily derived by following the steps outlined above.

V BRIEF of Stokes flow using a second order tensor fundamental solution

Fundamental solutions for Stokes flow that are second order tensor functions, ${S}_{ij}$ , can also be used to construct the flow field, $w_{i}$ of the general form

[TABLE]

where the constant vector, $D_{j}$ , is found by imposing the two conditions in Eq. (10).

One example of this class of solutions is the source potential doublet,Chwang and Wu (1975); Pozrikidis (2011) see item (II) of Table 2, that has zero pressure, $\mathcal{P}=0$ . The flow field $w_{i}$ that can be constructed is

[TABLE]

that satisfies Eq. (10a). The corresponding traction is

[TABLE]

The requirement of Eq. (10b) on the traction then implies

[TABLE]

in which $b={\hat{a\mkern 3.0mu}\mkern-3.0mu}{}^{d}_{l}n_{l}^{0}\neq 0$ . This gives

[TABLE]

Table 2 also contains the result for $w_{i}$ obtained using a Stokeson.

VI Examples

We now furnish examples to illustrate the implementation of the boundary regularized integral equation formulation (BRIEF) of single-phase Stokes flow problems and to demonstrate the many unique advantages of this approach. The key result of the BRIEF given by Eq. (12) relates the surface velocity to the surface traction. Therefore we first demonstrate the solution of two problems that are defined either by specifying the velocity field or the surface traction. We also use this example to quantify the precision that can be gained by using quadratic instead of linear surface elements while keeping the number of unknowns or degrees of freedom consistent, and to show the numerical equivalence of using different fundamental solutions to regularize the integral equation.

A second example will be used to demonstrate that it is straightforward to calculate the velocity field accurately not only far from but also close to boundaries using the BRIEF because of the complete absence of singular behavior. This is a very distinct advantage over the conventional boundary integral method in which singularities of the kernel can adversely affect the accurate evaluation of the velocities at field points close to the surface.

Another advantage of the BRIEF is that extreme geometric aspect ratios in the boundaries do not degrade the numerical precision of the solution. We illustrate this by considering the case of two nearly touching spheres in an external flow field that also serves the purpose of illustrating how forces and torques can be computed accurately by the BRIEF. In this and subsequent examples, we also consider cases in which neither the surface velocity nor the surface traction are specified, but rather a relation between these two quantities obeying the Navier slip boundary condition is imposed.

We also consider examples of lubrication flow between closely spaced surfaces where the field quantities become unbounded in the limit of zero separation. These examples illustrate the fact that in absence of mathematical singularities, the BRIEF is better positioned to handle the unavoidable physical divergences.

Finally, we consider the calculation of forces and torques experienced by bodies of varying geometric aspect ratios with varying Navier slip boundary conditions imposed on the surfaces.

Where appropriate, we compare results with that obtained from the conventional boundary integral method (CBIM) in which the rigid body solution has been subtracted:Becker (1992)

[TABLE]

There is a general point to note in selecting the desired flow field $\bm{w}$ given in Tables 1 and 2 to construct the BRIEF. For external problems, it is also necessary to consider the integral over the surface at infinity, $S_{\infty}$ , that grows at $r^{2}$ as $r\rightarrow\infty$ . For the cases of $\bm{w}$ that have been identified as applicable to both internal and external problems under the Domain category in Tables 1 and 2, the terms that involve $r_{d}$ decay faster than $r^{-2}$ and so such terms do not contribute to the integral over $S_{\infty}$ . The integral on the right hand side of Eq. (12) also vanishes on $S_{\infty}$ because $U_{ij}$ vanishes as $r^{-1}$ , $f_{i}$ as $r^{-2}$ and $\mathcal{F}_{i}$ as $r^{-3}$ as $r\rightarrow\infty$ . Finally, the limiting form of the left hand side of Eq. (12) on $S_{\infty}$ will depend on the choice of the flow field $w_{i}$ constructed from the fundamental solutions. As an example, for the case of a stresslet (entry V in Table I), the left hand side of Eq. (12) on $S_{\infty}$ is

[TABLE]

A notable exception is the linear solution for which the contribution on $S_{\infty}$ also vanishes due to cancellations owing to the symmetry of the angular integration when the domain is infinite. If this is not the case, for instance, in a semi-infinite domain, the linear solution cannot be used.

VI.1 Velocity or traction boundary conditions

We consider two examples in which we determine the velocity on a boundary that is either a sphere with free slip or zero tangential stress boundary conditions placed in an external flow field, $U\bm{k}$ , that is uniform at infinity, or on a boundary that is a sphere centered at the origin in a quiescent fluid and subjected to a surface traction $\bm{f}=(3\mu U/R^{2})z\bm{n}$ . In the former case, the velocity components parallel and perpendicular to the surface are given by Lamb (1932): $u_{\bot}=U\cos\theta\{1-(R/r_{c})\}$ , $u_{||}=U\sin\theta\{-1+\tfrac{1}{2}(R/r_{c})\}$ , where $r_{c}$ is the radial distance from the center of the sphere and $\theta$ is the angle between the vector $\bm{k}$ and the radial direction.

In these two examples, we compare in Fig. 2, the variation of the maximum absolute error in the surface velocity as a function of the number of degrees of freedom (or the number of unknowns) using linear and quadratic elements based on the same triangular surface mesh. With linear elements, we see that results from the conventional boundary integral method (CBIM) and from the boundary regularized integral equation formulation (BRIEF) using a linear function, a source potential or a stresslet with $\bm{x}_{d}$ at the center of the sphere to remove the singularities, gave practically identical results. However, if quadratic elements are used, the magnitude of the maximum absolute error decreases by almost an order of magnitude, or equivalently, if quadratic elements are used instead of linear elements, the same precision can be attained be reducing the degree of freedom by about a factor of 10. Thus the combination of the BRIEF with quadratic elements offers significant advantages in terms of precision and computational effort.

VI.2 Velocity field near boundaries

The boundary regularized integral equation formulation (BRIEF) of Stokes flow can also be used to calculate accurately the velocity at field points close to a boundary by alleviating the loss of precision in the conventional boundary integral method (CBIM) due to the near singular behavior of the integrals. We achieve this by using a simple extension of the method developed for the solution of the Laplace equationSun et al. (2014) to give a robust way to calculate the flow velocity by BRIEF anywhere within the flow domain using:

[TABLE]

Here $\bm{x}_{p}$ is a point in the flow domain for which the fluid velocity $u_{j}(\bm{x}_{p})$ is to be calculated, $w_{j}(\bm{x}_{p})$ is value of the constructed flow field at the same position, $U^{p}_{ij}=U_{ij}(\bm{x},\bm{x}_{p})$ , $T^{p}_{ijk}=T_{ijk}(\bm{x},\bm{x}_{p})$ , and $\bm{x}_{0}$ in $U_{ij}(\bm{x},\bm{x}_{0})$ , $T_{ijk}(\bm{x},\bm{x}_{0})$ is taken to be the node on the boundary, $S$ . When $\bm{x}_{p}$ is close to the surface, we choose $\bm{x}_{0}$ to be related to $\bm{x}_{p}$ by: $\bm{x}_{p}=\bm{x}_{0}+\alpha\bm{n}_{0}$ , where $\alpha$ is a small constant. With this choice, the near singular behavior of the term $({T}^{p}_{ijk}-{T}_{ijk})n_{k}$ is alleviated.

To illustrate the accuracy that can be obtained by using the BRIEF, we consider a sphere of radius, $R$ , with the free slip (zero tangential stress) boundary condition in a uniform flow field, $U\bm{k}$ , at infinity and compare the velocity in the equatorial plane $(z=0)$ as a function of position with the known analytical resultLamb (1932)

[TABLE]

From Fig. 3, we see that using the BRIEF plus the source potential (item IV in Table 1 with $\bm{x}_{d}$ at the center of the sphere) with 2352 linear trangular elements and 1178 nodes, the relative error is less than 1% at all positions. In contrast, the relative error using the CBIM with the same surface mesh can exceed 20% close to the sphere surface. This indicates that the BRIEF also alleviates the near singular behavior in the flow domain close to boundaries that is often more difficult to deal with than the singular behavior on the boundaries in CBIM.

VI.3 Nearly touching surfaces

Another significant advantage of the BRIEF of Stokes flow is in cases in which boundaries are very close together whereby the problem has disparate but important characteristic length scales. In the conventional boundary integral formulation (CBIM), the singular behavior of the kernel on one boundary will invariably have an adverse effect on the precision of integrals evaluated on the other nearby boundary even in instances in which the physical problem does not have any singular behavior. On the other hand, the complete absence of singular terms in the BRIEF for Stokes flow means that such numerical problems do not arise.

We illustrate this with the example of two nearly touching spheres of radius $R$ in a uniform flow field at infinity oriented either parallel or perpendicular to the line of centers. The sphere surfaces either have a free slip (zero tangential stress) condition or a no slip (zero velocity) condition. Results are given for the drag force, $\bm{F}$ , on each sphere calculated by integrating the traction on the surface, $\bar{S}$ , of that sphere

[TABLE]

and similarly for the torque, $\bm{N}$ , about the center of the sphere at $\bm{x}_{c}$ ,

[TABLE]

The force and torque are calculated as functions of the distance of closest approach, $h/R$ , between the spheres. The results are obtained using the source potential (item IV in Table 1 with $\bm{x}_{d}$ at the center of the sphere on which $\bm{x}_{0}$ is located) solution in the BRIEF with 2352 linear elements and 1178 nodes on each sphere.

For two spheres aligned with their line of centers along the uniform flow field, we show in Fig. 4, the drag force as a function of sphere separation for the combinations of free slip boundary condition on both spheres, no slip boundary conditions on both spheres, no slip boundary condition on the upstream sphere with free slip boundary condition on the downstream sphere as well as for two spheres that obey the Navier slip boundary condition with a prescribed slip length, $s$

[TABLE]

Here $u_{n}$ is normal component and $u_{t1}$ and $u_{t2}$ are the tangential components of the fluid velocity on the surface of the sphere. The tangential components of the traction acting on the sphere are denoted by $f^{S}_{t1}$ and $f^{S}_{t2}$ . The limit of the slip length $s=0$ , reduces to the familiar no-slip condition with $\bm{u}=0$ on the surface. The limit $s\rightarrow\infty$ is the perfect free-slip or zero tangential stress case with the boundary conditions

[TABLE]

In cases for which exact analytic solutions are available Reed and Morrison (1974), the agreement with BRIEF is excellent (see Fig 4).

In Fig. 5, we show corresponding results for the surface velocity vector field for the free slip/free slip and free slip/no slip combinations at a small distance of closest approach between the spheres: $h/R=0.01$ to illustrate the capabilities of the BRIEF. In Fig. 6, the traction field for all three combinations of free slip/free slip, free slip/no slip and no slip/no slip at the same separation are shown.

In Figs. 7-9, we show similar results for the case when the uniform external flow field, $U\bm{i}$ , is oriented perpendicular to the line of centers between the spheres. The drag force on each sphere in the direction of the external flow and the magnitude of the torque about the $y$ -axis are shown in Fig. 7 as a function of sphere separation, $h/R$ . Again the agreement with results from available exact analytic solutionsO’Neill and Majumdar (1970) is excellent. The surface velocity vector field for the free slip/free slip and free slip/no slip combinations and the traction field for all three combinations at $h/R=0.01$ are shown in Figs. 8 and 9.

A further test of the numerical precision of our implementation is to compute components of the forces and torques in the above examples that are expected to be zero from symmetry considerations. We found that the magnitudes of these components are of the order $10^{-10}$ times smaller than the non-zero component. These results illustrate the flexibility and accuracy for all types of boundary conditions.

VI.4 Lubrication flow

Problems in which lubrication effects within thin gaps between closely spaced boundaries are important provide additional challenges to the boundary integral methods. In addition to the mathematical singularities in the conventional boundary integral method, the physics of the problem can give rise to divergences such as an unbounded pressure field as the gap width approaches zero. Although the BRIEF removed the mathematical singularities, the divergent behavior of the relevant physical quantities must remain.

To illustrate such situations, we consider two examples: two nearly touching spheres with no-slip boundary condition moving with different velocities, $U$ and $2U$ , perpendicular to their line of centers (Fig. 10a), and two nearly touching spheres rotating with the same angular velocity, $\omega$ , about axes perpendicular to their line of centers (Fig. 11a). Due to the lubrication effect, the pressure within the gap becomes very high as the sphere separation decreases, and can become challenging numerically. To obtain satisfactory results corresponding to a minimum sphere separation, scaled by the radius, $h/R$ down to $10^{-3}$ , we employ a finer uniform mesh with 4800 linear triangular elements connected by 2402 nodes on each sphere - about twice as dense as previous examples. The results obtained using BRIEF with the source potential (Item IV in Table I) are in excellent agreement with exact bipolar coordinate system solution,O’Neill and Majumdar (1970) as shown in Figs. 10 and 11.

VI.5 Forces and torques on spheroids

As a final example, we consider forces experienced by a solid prolate spheroid translating with constant velocity, $U$ , along its axes as functions of the aspect ratio and varying Navier slip boundary conditions (Fig. 12). We also investigate the torque experienced by a prolate spheroid rotating about the major axes at constant angular velocity, $\omega$ (Fig. 13). The surface of the spheroid, with semi-major and semi-minor axes $a$ and $b$ , is defined by

[TABLE]

and is represented by a linear triangular mesh with 2352 elements connected by 1178 nodes. Analytical expressions for the forces and torques of such particles with the no slip boundary condition are given by Chwang and Wu.Chwang and Wu (1974, 1975) For comparison we also give results for spheroids with the Navier slip boundary condition for different slip lengths, $s$ , in which, the analytical solutions are given by LambLamb (1932) when $a=b$ (sphere).

Results for the forces on the prolate spheroid in external flow fields of different orientations are given in Fig. 12 and for the torques on rotating prolate spheroids are given in Fig. 13. In those cases for which analytical results are available, the agreement with the numerical results obtained using BRIEF is excellent.

VII Conclusions

In this paper, we have shown how the Stokes equation for single-phase flow with prescribed boundary conditions can be cast as a boundary regularized integral equation (12) with integrands that are completely free of singularities that normally exist in the traditional form of the boundary integral equation. The singularities are removed by subtracting an auxiliary known flow field $\bm{w}$ that has the same singular behaviour as in the traditional formulation of the problem. As a result, the singularities as well as the solid angle that appear in the conventional boundary integral formulation are removed analytically without the need to introduce additional regularization parameters.

The present formulation is therefore fundamentally different to the regularized Stokeslet approach Cortez (2001) that removes the singularities in the fundamental solutions by spreading a $\delta$ -function force over a small ball of radius $\epsilon$ , where the optimal value of $\epsilon$ that minimises error needs to be determined. The numerical effort in solving the linear system in BRIEF is comparable to that in CBIM but no special provisions are needed to compute the matrix coefficients. In this sense, it is more efficient that the contour integral formulation Bazhlekov, Anderson, and Meijer (2004).

We show that different auxiliary flow fields, $\bm{w}$ , can be constructed using different known fundamental solutions of the Stokes equation as summarized in Eqs. (21), (27) and Tables 1 and 2. The absence of singular terms in the integrals means that even problems with surfaces that are very close together will not suffer any loss of numerical precision due to the adverse influence of the singularity of one surface upon a nearby surface. The regular nature of the BRIEF also provides a numerically robust way to evaluate field quantities near boundaries that is often more difficult to deal with than the singular behavior on the boundaries of the conventional boundary integral method.

By reformulating the boundary integral equation to remove the traditional singularities analytically, rather than developing integration algorithms to handle the integration over the element that contains the singular behavior, the same triangular mesh can be readily used to represent linear or quadratic elements. In practical implementations, the absence of any singular terms allows a significant reduction in the amount of computer code required with a corresponding reduction in the possibility of coding errors. And as we have seen, the use of quadratic elements offers a substantial gain in numerical precision. This flexibility therefore provides a convenient way to check the accuracy of a calculation by comparing results obtained from using linear and quadratic elements.

In all our examples, there is no practical difference in the numerical results obtained from using the different forms of the auxiliary function $\bm{w}$ given in Tables I and II. This is an expected result and indicates the analytical stability of BRIEF in removing the mathematical singularities in the conventional formulation of the boundary integral equations. Any convenient choice of $\bm{w}$ will result in a numerically robust scheme for solving the boundary integrals equation.

With all the above advantages, the boundary regularized integral equation formulation (BRIEF) of the Stokes problem presented here should always be used in preference to the traditional approach. Since the original physical fluid problem has no singularities, it is also intuitively satisfying that the corresponding numerical scheme can indeed be free of any singular behavior of a purely mathematical origin. Although we have only discussed application of the BRIEF to single-phase flow, an obvious development is to adopt this approach to multi-phase flow problems.

Acknowledgements.

This work is supported in part by a Discovery Project Grant from the Australian Research Council to DYCC who is is a Visiting Scientist at the Institute of High Performance Computing and an Adjunct Professor at the National University of Singapore.

Appendix A Numerical implementation of BRIEF

We now give the details of the numerical implementation of BRIEF of Eq. (12). Take the BRIEF with simply linear solution (Item I in Table I) for example, after introducing Eqs. (14) and (III) into Eq. (12) by using the expression for $M_{jk}$ in Eq. (20) that is listed at Item I of Table I, the boundary regularized integral formulation equation can be rearranged so that all velocities appear on the left-hand side and all tractions on the right-hand side

[TABLE]

When the surface $S$ is partitioned into elements defined by $N$ nodes, the integrals in Eq. (A) over each element can be evaluated by quadrature and give rise to the matrix equation

[TABLE]

for the components of $\bm{u}$ and $\bm{f}$ at each node, see Fig 14. Here, $\mathsfbi{H}$ is the $3N\times 3N$ influence matrix corresponding to the velocity on the left-hand side of Eq. (A), and $\mathsfbi{G}$ is the $3N\times 3N$ influence matrix corresponding to the traction on the right-hand side of Eq. (A). The vectors $\widetilde{\bm{u}}$ and $\widetilde{\bm{f}}$ , both of size $3N$ , contain every component of $\bm{u}$ and $\bm{f}$ at all the nodes. The coefficients corresponding to terms $u_{i}^{0}$ and $f_{i}^{0}$ appear in the “x” items of the $3\times 3$ sub-matrices along the diagonals of the matrices $\mathsfbi{H}$ and $\mathsfbi{G}$ , respectively, and the coefficients corresponding to terms $u_{i}$ and $f_{i}$ appear in the “o” items in the matrices $\mathsfbi{H}$ and $\mathsfbi{G}$ , respectively.

In the traditional BIM formulation, the elements of the $3\times 3$ sub-matrices along the diagonals of the matrices $\mathsfbi{H}$ and $\mathsfbi{G}$ , indicated with an “x” within the square boxes in Fig. 14, involved integrals over the singularities of the fundamental solutions. Numerous schemes have been developed to ensure the accurate numerical evaluation of such integrals Gaul, Kögl, and Wagner (2003) and is one of the more arduous tasks in the numerical implementation of the BIM. But since the integrands in Eq. (A) are now free of singularities, all matrix entries can be obtained using, for example, the Gauss-Quadrature scheme to evaluate the integrals over each surface element.

To implement BRIEF with the linear solution for the auxiliary function $\bm{w}$ , the number of lines of code needed to compute the matrix elements in squares that contain “x” along the diagonals of the matrices $\mathsfbi{H}$ and $\mathsfbi{G}$ is less than 40% of that needed in the CBIM, where the local coordinate system transformationGaul, Kögl, and Wagner (2003) is applied to deal with the singularities.

Bibliography26

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Bazhlekov, Anderson, and Meijer (2004) I. B. Bazhlekov, P. D. Anderson, and H. E. H. Meijer, “Nonsingular boundary integral method for deformable drops in viscous flows,” Physics of Fluids 16 , 1064 (2004).
2Sun et al. (2013) Q. Sun, E. Klaseboer, B. C. Khoo, and D. Y. C. Chan, “Stokesian dynamics of pill-shaped Janus particles with stick and slip boundary conditions,” Physical Review E 87 , 043009 (2013).
3Li and Spagnolie (2014) L. Li and S. E. Spagnolie, “Swimming and pumping of rigid helical bodies in viscous fluids,” Physics of Fluids 26 , 041901 (2014).
4Bouffanais and Yue (2010) R. Bouffanais and D. K. P. Yue, “Hydrodynamics of cell-cell mechanical signaling in the initial stages of aggregation,” Physical Review E 81 , 041920 (2010).
5Subramaniam, Gee, and King (2013) D. R. Subramaniam, D. J. Gee, and M. R. King, “Deformable cell-cell and cell-substrate interactions in semi-infinite domain,” Journal of Biomechanics 46 , 1067–1074 (2013).
6Chan, Klaseboer, and Manica (2011) D. Y. C. Chan, E. Klaseboer, and R. Manica, “Film drainage and coalescence between deformable drops and bubbles,” Soft Matter 7 , 2235–2264 (2011).
7Liu (2009) Y. J. Liu, Fast Multiple Boundary Element Method: Theory and Applications in Engineering (Cambridge University Press, 2009).
8Becker (1992) A. A. Becker, The Boundary Element Method in Engineering: A complete Course (Mc Graw-Hill International (UK) Limited, 1992).