Boundary Integral Analysis for the Non-homogeneous 3D Stokes Equation

TL;DR

This paper introduces a boundary integral method for efficiently solving the non-homogeneous 3D Stokes equation by transforming volume integrals into boundary integrals, validated through numerical tests.

Contribution

It presents a novel boundary integral reformulation for the 3D Stokes equation that simplifies volume integrals into boundary integrals, enabling efficient numerical computation.

Findings

The method successfully reduces volume integrals to boundary integrals.

Numerical validation confirms the accuracy of the boundary integral approach.

The approach facilitates potential extensions to Navier-Stokes equations.

Abstract

A regular-grid volume-integration algorithm is developed for the non-homogeneous 3D Stokes equation. Based upon the observation that the Stokeslet ${\mathcal U}$ is the Laplacian of a function ${\mathcal H}$, the volume integral is reformulated as a simple boundary integral, plus a remainder domain integral. The modified source term in this remainder integral is everywhere zero on the boundary and can therefore be continuously extended as zero to a regular grid covering the domain. The volume integral can then be evaluated on the grid. Applying this method to the Navier-Stokes equations will require obtaining velocity gradients, and thus an efficient algorithm for post-processing these derivatives is also discussed. To validate the numerical implementation, test results employing a linear element Galerkin approximation are presented.

Tables

Table 1. Table 1: Mean square errors for boundary velocity and force for a homogenous Stokes problem posed on the unit sphere. The mixed boundary conditions are from two exterior point source locations, P 1 = ( − 2.0 , 0 , 0 ) subscript 𝑃 1 2.0 0 0 P_{1}=(-2.0,0,0) and P 2 = ( 1.5 , 0 , 0 ) subscript 𝑃 2 1.5 0 0 P_{2}=(1.5,0,0) .

		$N_E=376$			$N_E=1504$
		$x$	$y$	$z$	$x$	$y$	$z$
$P_1$	$𝐮$	1.094E-4	5.670E-5	3.634E-5	3.796E-5	1.409E-5	7.812E-6
	$\tau$	5.520E-4	2.042E-4	2.347E-4	2.214E-4	7.840E-5	9.376E-5
$P_2$	$𝐮$	2.487E-4	8.597E-5	1.320E-4	6.050E-5	1.985E-5	2.009E-5
	$\tau$	1.515E-3	1.327E-3	1.136E-3	4.502E-4	2.751E-4	2.504E-4

Table 2. Table 2: Mean square errors for boundary velocity and force for a homogenous Stokes problem posed exterior to the unit sphere. The mixed boundary conditions are from two exterior point source locations, P 3 = ( 0 , 0.7 , 0 ) subscript 𝑃 3 0 0.7 0 P_{3}=(0,0.7,0) and P 4 = ( 0 , 0 , 0.8 ) subscript 𝑃 4 0 0 0.8 P_{4}=(0,0,0.8) .

		$N_E=376$			$N_E=1504$
		$x$	$y$	$z$	$x$	$y$	$z$
$P_3$	$𝐮$	1.964E-3	1.749E-3	3.412E-3	4.472E-5	3.674E-5	4.184E-5
	$\tau$	1.468E-2	1.530E-2	1.152E-2	7.571E-4	6.693E-4	5.617E-4
$P_4$	$𝐮$	8.042E-4	7.977E-4	1.458E-3	2.281E-5	1.991E-5	4.091E-5
	$\tau$	4.614E-2	1.222E-2	2.792E-2	5.633E-3	2.772E-3	5.118E-3

Table 3. Table 3: Mean square errors for boundary velocity and force for a nonhomogenous Stokes problem posed on the unit sphere. The mixed boundary conditions are from the function ℋ ​ ( Q , P ) ℋ 𝑄 𝑃 {\mathcal{H}}(Q,P) with P 𝑃 P the exterior points P 1 = ( 2.0 , 0 , 0 ) subscript 𝑃 1 2.0 0 0 P_{1}=(2.0,0,0) , P 2 = ( 1.5 , 0 , 0 ) subscript 𝑃 2 1.5 0 0 P_{2}=(1.5,0,0) , P 3 = ( 0 , 1.3 , 0 ) subscript 𝑃 3 0 1.3 0 P_{3}=(0,1.3,0) , P 4 = ( 0 , 0 , 1.2 ) subscript 𝑃 4 0 0 1.2 P_{4}=(0,0,1.2) .

		$N_E=376$			$N_E=1504$
		$x$	$y$	$z$	$x$	$y$	$z$
$P_1$	$𝐮$	2.815E-4	3.495E-5	8.432E-5	7.382E-5	8.187E-6	1.908E-5
	$\tau$	6.332E-4	2.033E-4	3.076E-4	2.199E-4	7.934E-5	1.016E-4
$P_2$	$𝐮$	3.314E-4	3.644E-5	1.039E-4	8.111E-5	8.700E-6	2.392E-5
	$\tau$	6.721E-4	2.400E-4	3.604E-4	2.165E-4	8.181E-5	1.512E-4
$P_3$	$𝐮$	2.579E-4	2.300E-5	5.769E-4	6.686E-5	7.441E-6	1.336E-5
	$\tau$	5.881E-4	2.913E-4	2.790E-4	2.010E-4	9.461E-5	8.569E-5
$P_4$	$𝐮$	1.643E-4	1.322E-5	5.353E-5	5.299E-5	3.347E-6	1.315E-5
	$\tau$	6.246E-4	1.984E-4	3.690E-4	2.043E-4	7.412E-5	1.100E-4

Table 4. Table 4: Average square errors for boundary gradient components for (a) interior homogeneous problem with source P = ( 2 , 0 , 0 ) 𝑃 2 0 0 P=(2,0,0) ; (b) exterior homogenous with source P = ( 0 , 0.7 , 0 ) 𝑃 0 0.7 0 P=(0,0.7,0) ; and (c) non-homogeneous interior with P = ( 0 , 0 , 1.2 ) 𝑃 0 0 1.2 P=(0,0,1.2) for the unit sphere.

	$N_E=376$			$N_E=1504$
$m$	${𝐮}_{1,m}$	${𝐮}_{2,m}$	${𝐮}_{3,m}$	${𝐮}_{1,m}$	${𝐮}_{2,m}$	${𝐮}_{3,m}$
$1$	1.002E-4	6.311E-5	6.833E-5	1.817E-5	1.199E-5	1.400E-5
$2$	5.547E-5	5.718E-5	3.893E-5	9.240E-6	1.116E-5	7.820E-6
$3$	6.335E-5	4.120E-5	5.443E-5	9.940E-6	8.053E-6	1.168E-5
$1$	4.428E-4	6.620E-4	3.751E-4	3.957E-5	5.462E-5	3.926E-5
$2$	1.046E-3	5.296E-4	4.392E-4	4.657E-5	3.758E-5	3.472E-5
$3$	5.526E-4	4.901E-4	4.443E-4	2.818E-5	3.620E-5	2.649E-5
$1$	1.135E-4	6.056E-5	6.881E-5	2.187E-5	1.251E-5	1.435E-5
$2$	4.324E-5	5.383E-5	3.878E-5	9.990E-6	1.094E-5	8.033E-6
$3$	5.934E-5	3.758E-5	5.362E-5	1.198E-5	7.676E-6	1.114E-5

Equations

$${\int }_{\Omega }{\mathcal{U}}_{kj}(Q,P){\mathcal{F}}_j(Q)\mathrm{d}{\Omega }_Q,$$

$${\mathcal{F}}_j={\rho }_d{\mathbf{u}}_k\frac{\partial {\mathbf{u}}_j}{\partial x_k}.$$

$${\int }_{\mathcal{B}}{\mathcal{U}}_{kj}(Q,P)\left({\mathcal{F}}_j-{\mathcal{F}}_j^0\right)(Q)\mathrm{d}{\Omega }_Q.$$

$$\mathcal{E}(\mathcal{H})=\mathcal{G},$$

$$\mu {\nabla }^2\mathcal{H}$$

$$\nabla \stackrel{\bullet }{}\mathcal{H}$$

$$\mu {\nabla }^2\mathbf{u}-\nabla p$$

$$\nabla \stackrel{\bullet }{}\mathbf{u}$$

$$\displaystyle\int_{\Sigma}\bigg{[}{\mathcal{T}}_{ijk}(Q,P)\mathbf{u}_{i}(Q)\mathbf{n}_{j}(Q)-{\mathcal{U}}_{kj}(Q,P)\mathbf{\tau}_{j}(Q)\bigg{]}\,{\rm d}\Sigma_{Q}=\hphantom{---------}$$

$${\int }_{\Omega }{\mathcal{U}}_{kj}(Q,P){\mathcal{F}}_j(Q)\mathrm{d}{\Omega }_Q,$$

$$-\mathbf{u}(P)$$

$${\mathcal{U}}_{kj}(Q,P)$$

$${\mathcal{T}}_{ijk}(Q,P)$$

$${\mathcal{H}}_{kj}(Q,P)$$

$$\mu {\nabla }^2{\mathcal{H}}_{kj}$$

$$\nabla \stackrel{\bullet }{}\mathcal{H}$$

$${\nabla }^2{\mathcal{F}}^0=0$$

$$\mu {\int }_{\Omega }{\nabla }^2{\mathcal{H}}_{kj}(Q,P){\mathcal{F}}_j^0(Q)\mathrm{d}{\Omega }_Q+{\int }_{\Omega }{\mathcal{U}}_{kj}(Q,P)\left({\mathcal{F}}_j-{\mathcal{F}}_j^0\right)(Q)\mathrm{d}{\Omega }_Q.$$

$$\mu {\int }_{\Sigma }\left(\frac{\partial {\mathcal{H}}_{kj}}{\partial \mathbf{n}}(Q,P){\mathcal{F}}_j^0(Q)-{\mathcal{H}}_{kj}(Q,P)\frac{\partial F_j^0}{\partial \mathbf{n}}(Q)\right)\mathrm{d}{\Sigma }_Q.$$

$$\frac{\partial {\mathcal{H}}_{kj}}{\partial \mathbf{n}}=\frac{1}{32\pi {\mu }^{2}r}\left(3\mathbf{n}\stackrel{\bullet }{}R{\delta }_{kj}-\left({\mathbf{n}}_k{R}_j+{\mathbf{n}}_j{R}_k\right)+\mathbf{n}\stackrel{\bullet }{}R\frac{R_k{R}_j}{r^2}\right)$$

$${\int }_{\Sigma }{\psi }_k(P){\int }_{\Omega }{\mathcal{U}}_{kj}(Q,P)\left({\mathcal{F}}_j-{\mathcal{F}}_j^0\right)(Q)\mathrm{d}{\Omega }_Q\mathrm{d}{\Sigma }_P=$$

$${\int }_{\Omega }\left({\mathcal{F}}_j-{\mathcal{F}}_j^0\right)(Q){\int }_{\Sigma }{\psi }_k(P){\mathcal{U}}_{kj}(Q,P)\mathrm{d}{\Sigma }_P\mathrm{d}{\Omega }_Q.$$

$${\tilde{\mathcal{H}}}_{kj}={\delta }_{kj}{r}^2\left(a+b\log (r)\right)+\frac{{\partial }^2}{\partial q_k\partial q_j}\left\{r^4\left(c+d\log (r)\right)\right\}.$$

$${\tilde{\mathcal{U}}}_{kj}(Q,P)=-{\delta }_{kj}\log (r)+\frac{R_k{R}_j}{r^2},$$

$${\tilde{\mathcal{H}}}_{kj}$$

$${\left[\frac{1}{N}\sum _{k=1}^N{ϵ}_k^2\right]}^{1/2}$$

$$\mu {\nabla }^2{\mathcal{H}}_{kj}$$

$$\nabla \stackrel{\bullet }{}\mathcal{H}$$

$$u_{kj}=\frac{1}{32\pi {\mu }^2}r\left(3{\delta }_{kj}-\frac{R_k}{r}\frac{R_j}{r}\right).$$

$${\sigma }_{kjl}=\frac{1}{16\pi \mu }\left(\frac{R_k{R}_l{R}_j}{r^3}+\frac{R_l{\delta }_{kj}+R_k{\delta }_{jl}-R_j{\delta }_{kl}}{r}\right)$$

$$-{\mathbf{u}}_{k,m}(P)$$

$$0$$

$$\mathbf{u}_{k,m}(P)=Lim_{\leftrightarrow}\int_{\Sigma}\bigg{[}\,{\mathcal{T}}_{kjl,m}(Q,P)\mathbf{u}_{i}(Q)\mathbf{n}_{j}(Q)-{\mathcal{U}}_{kj,m}(Q,P)\mathbf{\tau}_{j}(Q)\,\bigg{]}\,{\rm d}\Sigma_{Q}\,\,,$$

$${\mathcal{U}}_{kj,m}(Q,P)$$

Full text

1 Introduction

Boundary integral analysis for fluids has primarily involved potential (inviscid irrotational) and slow viscous (linear Stokes) flow [1]. In particular, the break-up and coalescence of free surface flows using integral formulations of the Laplace [2, 3, 4, 5, 6] and Stokes [7, 8, 9, 10, 11] equations have a long and successful history. Other examples of Stokes flow applications include rigid particle suspensions [12], vesicle suspensions [13, 14, 15], surfactant flows [16], electrohydrodynamics [17], microswimmers [18, 19], biological cell modeling [20] including mitotic cell division [21, 22], viscous erosion [23, 24], and micro-electro-mechanical systems (MEMS)[25, 26].

With a more general fluid equation, the conversion to a boundary integral statement cannot be completed exactly. To take an important example, the Navier-Stokes equation for a viscous incompressible fluid [27] on the domain $\Omega$ becomes the standard Stokes boundary integral together with a volume term of the form

[TABLE]

where ${\mathcal{U}}$ is the Stokes Green’s function and (with the summation convention and $\rho_{d}$ the fluid density)

[TABLE]

is the non-linear convection term [28]. Similarly, the integral formulation for a non-Newtonian fluid must also contend with a ‘pseudo body force’ domain integral of the form in Eq. (1) [29]. For the moving boundary simulations mentioned above, constructing a finite element volume mesh at every time step would be computationally demanding, if not impossible considering the extreme geometries. To overcome this difficulty, we present here a ‘body-fitted’ volume discretization for the non-homogeneous 3D Stokes equation.

Methods that have been employed for volume integral evaluation include direct finite element mesh calculations [28, 29, 30, 31], and approximate conversion to a boundary integral using either dual reciprocity [32, 33] or interior line integration [34, 35]. The embedded boundary method in [36] and the fast Poisson solvers [37, 38] exploit an easily constructed regular grid covering the domain, in conjunction with a Fast Multipole Method.

In particular, the volume integration in [37] for the Poisson solver and the work herein both rely on the continuous extension of a function from $\Omega$ to the grid. In the former, the source ${\mathcal{F}}$ is extended by solving an exterior Dirichlet boundary integral equation, with boundary values given by ${\mathcal{F}}$. For the Stokes algorithm herein, as well as the previous Laplace and elasticity treatments discussed in [39, 40, 41, 42], a function ${\mathcal{F}}^{0}$ is defined on $\Omega$ by means of an interior boundary integral solution (again with boundary values from ${\mathcal{F}}$). With ${\mathcal{B}}$ denoting the covering box, and with the help of Green’s Theorem, the modified body force integral becomes

[TABLE]

Here, ${\mathcal{F}}-{\mathcal{F}}^{0}$ is defined as zero outside $\Omega$ and is therefore continuous on the entire grid. The evaluation of this integral can be carried out by decomposing ${\mathcal{B}}$ into (uniform) cuboid cells, and using simple linear interpolation over all cells. In other words, with the continuous zero extension the ‘partial cells’ that straddle the boundary are treated as any other cell and the location of the domain boundary can effectively be ignored. Moreover, as ${\mathcal{F}}_{j}-{\mathcal{F}}^{0}_{j}\equiv 0$ on the boundary, the Green’s function singularity does not present any problem. Thus, unlike for the Poisson algorithm in [37], it is not necessary to refine the cells near the boundary. The calculations for ${\mathcal{F}}^{0}$, including interior point evaluations at cell vertices, are however expensive; nevertheless, as in [36, 37, 38, 43], fast methods can be employed [41].

A second key aspect of the work in [39, 40, 41, 42] is the construction (in a simple analytic form) of a function ${\mathcal{H}}$ that satisfied

[TABLE]

where ${\mathcal{E}}$ is the partial differential equation and ${\mathcal{G}}$ the corresponding Green’s function. For the 3D Stokes equation, it will turn out that the function ${\mathcal{H}}$ will satisfy the Stokes equations with zero pressure. That is,

[TABLE]

where again ${\mathcal{U}}$ is the Stokeslet, the Laplacian is applied componentwise and the divergence is applied to the columns of ${\mathcal{H}}$. The numerical treatment of the Stokes volume integral therefore has much in common with the Laplace algorithm in [39], and this somewhat simplifies the numerical implementation.

It is expected that this domain integral algorithm can form the basis of an effective numerical solution of the two non-linear problems mentioned above, incompressible Navier-Stokes and non-Newtonian flows. For the former, the body force is given in Eq. (2), and consequently the computation of ${\mathcal{F}}^{0}$ will require obtaining velocity gradients on the boundary. The evaluation of these derivatives will therefore be considered in Section 4. Although a non-linear solution algorithm will not be considered herein, it is worth noting that previous papers [44, 45] have successfully dealt with non-linear equations in conjunction with the regular grid algorithm.

2 $3$D Stokes’ Equation

The equations for Stokes flow in three dimensions are [1]

[TABLE]

where the problem domain is $\Omega$ and $\{\mathbf{u},p,\mu\}$ are fluid velocity, pressure, and viscosity. The source function ${\mathcal{F}}$ can be a given body force and/or nonlinear terms stemming from, for example, the Navier-Stokes or non-Newtonian equations.

The Stokes boundary/volume integral statement is [1, 28, 36]

[TABLE]

where $\Sigma$ is the boundary of $\Omega$, $\mathbf{n}$ the exterior unit normal , $\mathbf{\tau}(Q)$ the surface force (normal component of the stress), ${\mathcal{U}}_{kj}(Q,P)$ the Stokeslet (Green’s function) and ${\mathcal{T}}_{ijk}(Q,P)$ the corresponding stresslet. In this equation, it is assumed that the singular integrals are defined as a limit as $P$ approaches the boundary from outside the domain (this only comes into play for the ${\mathcal{T}}$ kernel). For later use (see Section 4), there is also an interior limit equation, ([1], Eq. (2.3.11))

[TABLE]

the only difference being the explicit appearance of the ‘free term’ $\mathbf{u}(P)$.

With $Q=\{q_{\ell}\}$, $P=\{p_{\ell}\}$, $\mathbf{R}=\{R_{\ell}\}=Q-P$ and $r=\|\mathbf{R}\|$ the distance, the fundamental solutions are given by [1]

[TABLE]

where $\delta_{kj}$ is the Kronecker delta function.

As with the Laplace and elasticity volume algorithms [39, 41], we seek to represent the kernel ${\mathcal{U}}$ as a derivative in such a way that Green’s Theorem can be used to reformulate the volume integral in Eq. (2). Based upon the isotropic elasticity analysis in [41], and the fact that the Stokeslet is the elasticity Green’s function with Poisson ratio $\nu=1/2$ and shear modulus equal to viscosity, we define the function ${\mathcal{H}}$ as

[TABLE]

Direct calculation (see Appendix A) shows that

[TABLE]

which says that ${\mathcal{H}}$ satisfies the zero pressure non-homogeneous Stokes equations, with the Stokeslet as the pseudo body force. This is completely analogous to the Green’s functions representations for Laplace and elasticity that were employed in [39, 41]. Fom this point the volume integral treatment follows as in these previous articles. However, as Eq. (2) is componentwise, the Stokes formulation more closely resembles the scalar Laplace implementation than the vector elasticity.

Briefly then, defining the vector ${\mathcal{F}}^{0}$ by

[TABLE]

with boundary conditions {\mathcal{F}}^{0}\big{|}_{\Sigma}={\mathcal{F}}\big{|}_{\Sigma}, and using Eq. (2) we can write the volume integral as

[TABLE]

Invoking Green’s Theorem and Eq. (12), the first term becomes the boundary integral

[TABLE]

The ${\mathcal{F}}^{0}$ boundary flux is obtained from a boundary integral solution of Eq. (12). Note that unlike their counterparts in Eq. (9), ${\mathcal{H}}(Q,P)$ and its normal derivative

[TABLE]

are rather innocous at $Q=P$. ${\mathcal{H}}$ behaves as $r$, and is therefore continuous, while the normal derivative is merely discontinuous.

2.1 Numerics: Remainder Volume Integral

All boundary integrals in this work are approximated using a standard Galerkin procedure with linear triangular elements. As just noted, the kernel functions in Eq. (14) do not diverge at $Q=P$, and thus the treatment of Eq. (14) is somewhat simplified. In this case all integrals are computed numerically, except for the coincident integral involving the normal derivative of ${\mathcal{H}}$. Techniques for the singular and non-singular boundary integrals are well known, so this section will briefly discuss the remainder volume integral Eq. (16). Further details can be found in previous papers [39, 41].

For the volume term, the Galerkin form can be exploited by interchanging of the order of integration [31]

[TABLE]

The integrand in the boundary integral is known in analytic form and thus this integral can be handled partially analytically. In the regular grid approach, the points $Q$ will be the cell vertices, and the interchange is especially useful for the ‘near-singular’ case when $Q$ is close to the boundary.

The other key aspect of Eq. (16) is that by design ${\mathcal{F}}_{j}-{\mathcal{F}}^{0}_{j}\equiv 0$ on $\Sigma$. This allows a continuous extension of this function as zero outside the domain , and permits evaluation using a regular grid that covers the domain $\Omega$. Using simple linear interpolation, the volume integration becomes simply a sum over the vertex values. The required interior values of ${\mathcal{F}}^{0}$ can be computed from the interior point boundary integral equation once the boundary ${\mathcal{F}}^{0}$ flux has been obtained. Again, the reader is asked to consult [39, 41] for further details.

2.2 $2$D Stokes’ Equation

Although the discussion herein is for three dimensions, for completeness we include the expression for the 2D version of Eq. (10). The ${\mathcal{H}}$ function for 2D Stokes, denoted by $\tilde{{\mathcal{H}}}$, is related to that for elasticity and, as shown in [41], the general form for elasticity is (dropping the normalization constant)

[TABLE]

In terms of the Poisson ratio, $\nu$, the coefficients are $b=\nu-1$, $d=(3-4\nu)/(64(1-\nu))$ and $a$, $c$ are only required to satisfy $64c(1-\nu)+2a(3-4\nu)=8\nu^{2}-9\nu+2$. For elasticity there is therefore leeway in the choice of these last two coefficients.

The Stokes equation corresponds to $\nu=1/2$ and in this case $b=-1/2$, $d=1/32$, and $16c+a=-1/4$. However, to satisfy $\nabla\!\stackrel{{\scriptstyle\scriptscriptstyle\bullet}}{{{}}}\!\tilde{{\mathcal{H}}}=0$ and $\mu\nabla^{2}\tilde{{\mathcal{H}}}=\tilde{{\mathcal{U}}}$, where $\tilde{{\mathcal{U}}}_{kj}$ is the two dimensional the Stokelet [1] (Eq. $2.6.17$)

[TABLE]

it is necessary that $c=-1/16$ and hence $a=3/4$. Thus, in two dimensions

[TABLE]

A Maple code that verifies this result can be found in the Appendix.

3 Test Calculations

3.1 ${\mathcal{F}}=0$

Before testing the evaluation of the volume integral, we first provide evidence that the boundary integral in Eq. (2) (i.e., without the volume integral) has been implemented correctly. To this end, we solve the homogeneous Stokes equations interior and exterior to the unit sphere, with velocity boundary conditions on the upper half, $z\geq 0$, and force boundary data on the lower half. These boundary values, as well as the corresponding exact solution, are taken from a point source placed exterior to the problem domain, specifically, the first column of the Green’s functions in Eq. (9). Table 1 lists the average square nodal errors

[TABLE]

for each component of the computed velocity and force. Discretizations with $N_{E}=376$ and $N_{E}=1504$ elements and two exterior points $P_{1}=(-2.0,0,0)$ and $P_{2}=(1.5,0,0)$ were employed. As expected, the errors are smaller for the velocity, and there is reasonable decay in the errors for the finer mesh. Moreover, as $P_{2}$ is closer to the boundary, these errors are appropriately larger than for $P_{1}$.

Table 2 lists the corresponding errors for the exterior Stokes problem, the two point sources located at $P_{3}=(0,0.7,0)$ and $P_{4}=(0,0,0.8)$. The results again indicate that the Stokes boundary integral equation has been implemented correctly.

3.2 ${\mathcal{F}}\neq 0$

To verify the evaluation of the volume integral, Eq. (14) and Eq. (16), a problem similar to the above point source tests is employed. Note that Eq. (2),

[TABLE]

states that the columns of ${\mathcal{H}}$ satisfy the Stokes equations, with zero pressure and with forcing function ${\mathcal{F}}$ given by the corresponding column of ${\mathcal{U}}$. Thus, for an exterior point $P$, $\mathbf{R}=Q-P$, $r=\|\mathbf{R}\|$, and selected column $j$, the velocity field is

[TABLE]

The corresponding stress field for ${\mathcal{H}}$ is

[TABLE]

As with the point source test, the domain is the unit sphere, and velocity is specified for $z\geq 0$ and traction for $z<0$. For the remainder volume integral, the box that covers the sphere is $-1.1<\{x,y,z\}<1.1$, subdivided with a $40\times 40\times 40$ grid. Table 3 lists the errors in the computed velocity and traction for four exterior points located at $P_{1}=(2.0,0,0)$ $P_{2}=(1.5,0,0)$, $P_{3}=(0,1.3,0)$ and $P_{4}=(0,0,1.2)$. The results confirm that the volume terms, the $F_{0}$ boundary integral equation and the remainder volume integral, have been implemented correctly.

4 Velocity Gradients

An important potential application of the regular grid algorithm would be a ‘boundary mesh only’ solution of the Navier-Stokes equations for viscous incompressible flow. The pseudo body force ${\mathcal{F}}$, Eq. (2), is a function of velocity gradients [28, 32], and thus the boundary conditions for ${\mathcal{F}}^{0}$ require knowledge of these derivatives on $\Sigma$ (gradient values will also required at interior grid vertices, but herein we deal solely with the boundary). In this section, the Stokes implementation of the gradient method proposed in [46] will be discussed.

To start, it is a simple matter to differentiate the interior limit equation Eq. (8) to obtain an expression for the gradient components,

[TABLE]

where the subscript $,m$ indicates differentiation with respect to the coordinate $P_{m}$. Note that once Eq. (2) has been solved, everything on the right hand side of Eq. (24) is known. Nevertheless, a simple implementation of this equation is problematic: the computation is clearly very expensive and it involves the hypersingular kernel ${\mathcal{T}}_{kjl,m}(Q,P)$. Moreover, given the approximations involved, Eq. (24) would probably not produce a highly accurate result. As a result, a variety of alternative boundary integral gradient algorithms that have been considered, [47, 48, 49, 50, 51] being just a partial list. The papers [52, 53] and the references therein provide a more complete overview of this topic.

To address the computational cost, the key observation in [46] is that, unlike for the basic velocity and traction boundary integral equations, the exterior limit gradient equation

[TABLE]

and the interior limit are not the same equation. Thus, subtracting the two equations, interior minus exterior, yields a new and useful expression for the gradient. The only terms that do not cancel in the limit difference are those that are discontinuous crossing the boundary, and this immediately eliminates all non-singular boundary integrals. Moreover, in the volume integral, the Stokeslet singularity is sufficiently weak that this term is also continuous across the boundary, and we can therefore write

[TABLE]

where $Lim_{\leftrightarrow}$ indicates the limit difference.

A Galerkin implementation is employed to deal with the second issue, the hypersingular kernel ${\mathcal{T}}_{ijk,m}(Q,P)$. With Galerkin, Eq. (26) reduces to the coincident and adjacent edge singular integrations for the ${\mathcal{T}}$ integral, while for the weaker ${\mathcal{U}}_{,m}$ singularity it is solely the coincident integral that contributes. Although the integration work is minimal, there is an additional computational expense: note that the Galerkin implementation of Eq. (26) couples the nodal values of a gradient component. As a consequence, for each of the nine components $\mathbf{u}_{k,m}$, $1\leq k,m\leq 3$ on the boundary, the solution of an $N\times N$ linear system is required, $N$ the number of boundary nodes. However, the coefficient matrix is the same for each component, and moreover it is sparse, symmetric positive definite. Thus, for moderate sized problems, only one matrix factorization is required, while for large scale problems an efficient iterative solver can be employed.

Differentiating Eq. (9) with respect to the coordinates of $P$ the derivative kernels are

[TABLE]

The coincident and adjacent edge singular integrals are evaluated as in [46, 54], and the details will not be repeated here. However, it is worth noting that much of the calculation is exact: the boundary limit is handled analytically, and the singular integrals are computed partially analytically. Importantly, this allows the exact cancellation of a potentially divergent term (of the form $1/\epsilon$, $\epsilon\rightarrow 0$ the distance to the boundary) that arises in the coincident ${\mathcal{T}}_{,m}$ integration. Further details about the singular integration can be found in the cited references.

4.1 Gradient Tests

To verify the implementation of the algorithm, complete velocity gradients $\mathbf{u}_{k,m}$ have computed for several of the test problems discussed in Section 3. Table 4 lists the mean square gradient errors for (a) the interior homogeneous problem with point source at $P=(2,0,0)$; (b) the exterior homogeneous problem with point source at $P=(0,0.7,0)$; and (c) the non-homogeneous interior problem with $P=(0,0,1.2)$. For the last example the exact solution is the ${\mathcal{H}}$ derivatives

[TABLE]

with $j=1$ . The results demonstrate that the volume integral does in fact cancel out of the limit difference equation.

At first it might appear strange that the gradient errors are smaller than for the corresponding computed velocity solution. However, recall that these tests are mixed boundary value problems, with velocity boundary conditions specified on half the sphere. Naturally, the gradient calculation is strongly dependent on the input surface velocity, and thus in these examples the algorithm is working with exact data on half the boundary.

5 Conclusion

A regular grid volume integration algorithm for the non-homogeneous 3D Stokes equation has been presented. The key to modifying the original volume integral is to represent the Green’s function (Stokeslet) as the Laplacian of a function ${\mathcal{H}}$. This is analogous to previous volume integral treatments for the Laplace [39, 40] and elasticity equations [41, 42], as the Laplacian can be viewed the zero-pressure Stokes equations. With the function ${\mathcal{H}}$, the domain integral exactly transforms to a simple boundary integral, plus a volume term wherein the modified source function is everywhere zero on the boundary. The continuous zero extension of this source allows this volume integral to be computed on a regular grid of cells covering the domain.

An effective Stokes volume integral technique allows the possibility of treating nonlinear equations, e.g., Navier-Stokes and non-Newtonian equations. Nonlinear analyses with the regular grid approach have been previously carried out [44, 45], and as well there have been nonlinear integral equation solutions obtained with other volume methods [36, 28, 32, 29]. It is therefore reasonable to expect that this work will lead to an efficient integral equation algorithm for the nonlinear fluids noted above. Regarding the Navier-Stokes solver, the grid algorithm requires boundary values of the surface gradient, and it has been shown herein that this post-processing calculation can be executed efficiently.

It would be of interest to repeat the successful potential flow studies of coalescence [6] and Rayleigh-Taylor break-up [4, 5] with the more complex fluids mentioned above. The Laplace calculations employ cylindrical coordinates $\{\rho,\theta,z\}$, and with $\theta$ integrated out they become two-dimensional $\{\rho,z\}$ analyses. The axi-symmetry also requires that any body force is independent of the polar angle, and in this situation it is straightforward to implement the regular grid method [39]. For Stokes however, axi-symmetry means that the body force vector, rather than being independent of $\theta$, rotates properly with $\theta$ [1]. As an example, the surface normal for an axi-symmetric geometry is not independent of $\theta$, but rather of the form $\mathbf{n}=\left(n_{x}\cos(\theta),n_{y}\sin(\theta),n_{z}\right)$. With ${\mathcal{F}}$ of this form, the volume algorithm does not immediately carry over, and a modified algorithm must be developed. We hope to return to this issue in the future.

As a second application, we hope to combine the framework developed here with the ideas laid out in [55] to simulate general viscoelastic flows. In this context, a polymeric stress field develops in conjunction with the flow and serves as the body force in the non-homogeneous Stokes system. Given an initial polymeric stress, the boundary-integral framework could be used to determine the corresponding flow, which would then be used to update the stress field at the next time step. Such a method would thus capture the nonlinear feedback between the polymeric stress and the flow with accuracy and efficiency.

Acknowledgements

The authors are grateful to Prof. B. Quaife for kindly pointing out several important references. L. J. Gray gratefully acknowledges an MTS Visiting Professorship Grant at the University of Minnesota, and he would like to thank Profs. S. Mogilevskaya, J. Labuz and S. Crouch for the hospitality at the Department of Civil Engineering. The participation of J. Jakowski was facilitated by J. Williams and Dr. D. Pickel of the Mathematics Department, Oak Ridge High School. M. N. J. Moore acknowledges support from the Simons Foundation Collaboration Grants for Mathematicians, award ID 524259.

Appendix A

The Maple codes below can be used to confirm Eq. (2) for 3D and 2D. That is, they establish that the Stokes equations, with zero pressure, applied to ${\mathcal{H}}_{kj}$ and $\tilde{{\mathcal{H}}}_{kj}$ yield the corresponding Stokes Green’s function, and these functions satisfy conservation of mass.

${\mathcal{H}}$ for $3D$

q := array(1..3); p := array(1..3); r := array(1..3); lapH := array(1..3); divH := array(1..3); EQ := array(1..3,1..3); H := array(1..3,1..3); G := array(1..3,1..3); del := array(1..3,1..3):

for k from 1 to 3 do for j from 1 to 3 do del[k,j] := 0; od; del[k,k] := 1; od;

r[1] := q[1] - p[1]; r[2] := q[2] - p[2]; r[3] := q[3] - p[3]; rr := sqrt( r[1]*r[1] + r[2]*r[2] + r[3]*r[3] );

Stokeslet

for k from 1 to 3 do for j from 1 to 3 do G[k,j] := del[k,j]/rr + r[k]r[j]/rr^3; G[k,j] := G[k,j] / (8Pi*mu); od; od;

Elasticity function

a := 2 - 2nu; b := -(3-4nu)/(24*(1-nu)); nu := 1/2; for k from 1 to 3 do for j from 1 to 3 do H[k,j] := arrdel[k,j] + bdiff(rr^3,q[k],q[j]); H[k,j] := H[k,j] / (8Pi*mu^2); od; od;

loop over columns of H

for j from 1 to 3 do

lapH[1] := diff(H[1,j],q[1],q[1]) + diff(H[1,j],q[2],q[2]) + diff(H[1,j],q[3],q[3]); lapH[2] := diff(H[2,j],q[1],q[1]) + diff(H[2,j],q[2],q[2]) + diff(H[2,j],q[3],q[3]); lapH[3] := diff(H[3,j],q[1],q[1]) + diff(H[3,j],q[2],q[2]) + diff(H[3,j],q[3],q[3]);

lapH[1] := factor(normal(lapH[1])); lapH[2] := factor(normal(lapH[2])); lapH[3] := factor(normal(lapH[3]));

divH[j] := diff(H[1,j],q[1]) + diff(H[2,j],q[2]) + diff(H[3,j],q[3]); divH[j] := factor(normal(divH[j]));

Stokes’ equation mu*Lap(u) = grad(p) with p = 0

Stokes(H) = G

EQ[1,j] := factor(normal(expand( mulapH[1] - G[1,j] ))); EQ[2,j] := factor(normal(expand( mulapH[2] - G[2,j] ))); EQ[3,j] := factor(normal(expand( mu*lapH[3] - G[3,j] )));

od; # end column loop

$\tilde{{\mathcal{H}}}$ for $2D$

R := array(1..2);

del[1,1] := 1; del[2,2] := 1; del[1,2] := 0; del[2,1] := 0;

R[1]:=x[1]-p[1]: R[2]:=x[2]-p[2]: rsq := R[1]*R[1] + R[2]*R[2]: r := sqrt(rsq):

nu := 1/2; b := nu - 1; d := (3 - 4nu) / ( 64(1-nu) ); a := -1/4 -16*c; c := -1/16;

for k from 1 to 2 do for j from 1 to 2 do H[k,j]:= rsq*( a + bln(rsq)/2 ) * del[k,j] + diff( rsq^2( c + dln(rsq)/2 ),x[k],x[j]): LapH[k,j] := diff(H[k,j],x[1],x[1]) + diff(H[k,j],x[2],x[2]): LapH[k,j] := subs( (2x[1]-2p[1])^2 = 4(x[1]-p[1])^2, (2x[1]-2p[1])^4 = 16*(x[1]-p[1])^4, (2x[2]-2p[2])^2 = 4*(x[2]-p[2])^2, (2x[2]-2p[2])^4 = 16*(x[2]-p[2])^4,LapH[k,j]); G[k,j] := -del[k,j]*ln(rsq)/2 + R[k]*R[j]/rsq; chk[k,j] := normal(expand( LapH[k,j] - G[k,j] )); od; DivH[k] := diff(H[k,1],x[1]) + diff(H[k,2],x[2]): DivH[k] := normal(expand( DivH[k] )); od;