Fast Ewald summation for Stokes flow with arbitrary periodicity

2025-04-22 0 0 2.71MB 54 页 10玖币

侵权投诉

Fast Ewald summation for Stokes ﬂow with arbitrary periodicity

Joar Bagge∗, Anna-Karin Tornberg

KTH Mathematics, Linn´e FLOW Centre/Swedish e-Science Research Centre,

Royal Institute of Technology, SE-100 44 Stockholm, Sweden

Abstract

A fast and spectrally accurate Ewald summation method for the evaluation of stokeslet, stresslet and rotlet

potentials of three-dimensional Stokes ﬂow is presented. This work extends the previously developed Spectral

Ewald method for Stokes ﬂow to periodic boundary conditions in any number (three, two, one, or none) of the

spatial directions, in a uniﬁed framework. The periodic potential is split into a short-range and a long-range

part, where the latter is treated in Fourier space using the fast Fourier transform. A crucial component of the

method is the modiﬁed kernels used to treat singular integration. We derive new modiﬁed kernels, and new

improved truncation error estimates for the stokeslet and stresslet. An automated procedure for selecting

parameters based on a given error tolerance is designed and tested. Analytical formulas for validation in the

doubly and singly periodic cases are presented. We show that the computational time of the method scales

like O(Nlog N) for Nsources and targets, and investigate how the time depends on the error tolerance

and window function, i.e. the function used to smoothly spread irregular point data to a uniform grid. The

method is fastest in the fully periodic case, while the run time in the free-space case is around three times

as large. Furthermore, the highest eﬃciency is reached when applying the method to a uniform source

distribution in a primary cell with low aspect ratio. The work presented in this paper enables eﬃcient and

accurate simulations of three-dimensional Stokes ﬂow with arbitrary periodicity using e.g. boundary integral

and potential methods.

Keywords: fast summation, Stokes potentials, creeping ﬂow, reduced periodicity, Fourier analysis,

boundary integral equations

1. Introduction

Stokes ﬂow, also known as creeping ﬂow or viscous ﬂow, is a model of ﬂuid ﬂow in which inertial forces

are assumed to be negligible in comparison to viscous forces (i.e., the Reynolds number is very small).

This is often a valid assumption for phenomena involving suspension ﬂows on the micro- and nanoscales

(often in combination with Brownian motion), such as swimming microorganisms [1], cell dynamics [2],

microﬂuidic devices [3], gels [4,5], dynamics of nanoparticles and nanoﬁbrils [6,7,8,9], electrolytes [10],

and antibodies [11]. For Stokes ﬂow, the Navier–Stokes equations reduce to the Stokes equations, which for

an incompressible Newtonian ﬂuid are given by

−∇p(x) + µ∇2u(x) + g(x) = 0,(1)

∇ · u(x) = 0.(2)

Here, pis the pressure, uis the ﬂuid velocity, gis the body force per unit volume acting on the ﬂuid,

and µis the viscosity of the ﬂuid. We will here consider the nondimensionalized Stokes equations, which is

equivalent to setting µ= 1.

∗Corresponding author

Email addresses: joarb@kth.se (Joar Bagge), akto@kth.se (Anna-Karin Tornberg)

Preprint submitted to Journal of Computational Physics October 5, 2022

arXiv:2210.01255v1 [math.NA] 3 Oct 2022

In boundary integral and potential methods for Stokes ﬂow, the fundamental solution (Green’s function)

of the Stokes equations appears, namely the stokeslet kernel. In the three-dimensional case, the stokeslet is

a 3 ×3 tensor given by

Sjl(r) = δjl

|r|+rjrl

|r|3,(3)

where δjl is the Kronecker delta. Also commonly used are the stresslet and rotlet kernels, which can be

seen as derivatives of the stokeslet and will be introduced in section 2. In this paper, we are interested in

problems with periodic boundary conditions, in which Dof the three spatial directions will be periodic, and

the remaining 3 −Ddirections will be free, in this context meaning that the domain extends to inﬁnity

and that no boundary conditions are enforced in these directions in the summation procedure. We will call

the problem triply, doubly and singly periodic if D= 3,2,1, respectively, and free-space if D= 0. We will

consider a system of Npoint sources of strengths f(xn) located at positions xn. The velocity ﬁeld generated

by this system is given as a periodic sum over the point sources, i.e.

uDP(x) =

n=1 X

p∈PDP

S(x−xn+p)f(xn),(4)

where the set PDPis a D-dimensional lattice containing all periodic images, to be properly deﬁned in

section 2. We assume that we want to evaluate (4) also in Ntarget points, which may or may not be the

same as the source points. The periodic sum can be computed using Ewald summation, which splits the sum

into a short-range part, to be summed directly, and a smooth long-range part, to be treated in Fourier space.

For the stokeslet, the split was derived by Hasimoto [12]. A decomposition parameter ξ, to be introduced

in section 2, controls the Ewald summation split. For ﬁxed ξ, computing the Ewald sums directly leads to a

method with time complexity O(N2). Adjusting ξproperly reduces the complexity of the direct summation

to O(N3/2), see e.g. [13,14]. Fast Ewald summation methods such as the Particle–Mesh–Ewald (PME)

method [13] and smooth Particle–Mesh–Ewald (SPME) method [15,16] reduce the complexity further to

O(Nlog N), a signiﬁcant improvement.

The Spectral Ewald (SE) method is a PME-type method with spectral accuracy. As in all PME methods,

the interactions between source points and target points are computed via a uniform grid and the fast Fourier

transform (FFT). A so-called window function is used to interpolate between source/target points and the

uniform grid, much like in the nonuniform fast Fourier transform (NUFFT) [17,18]. The spectral accuracy

of the SE method comes from the choice of window function, which was originally a truncated Gaussian. In

contrast, other methods such as e.g. SPME [15,16] use Cardinal B-splines for interpolation, which leads to

algebraic accuracy.

The SE method for Stokes ﬂow has been developed in a series of papers, ﬁrst for the triply periodic

stokeslet [19], doubly periodic stokeslet [20], triply periodic stresslet [21], and triply periodic rotlet [22].

The method was extended to the free-space case for all three kernels by [23]. The current paper serves to

complete this development, much in the same way as was recently done for electrostatics by [24], by adding

the missing pieces (singly periodic case for all three kernels, and doubly periodic case for the stresslet

and rotlet), and unifying all periodic cases within a single framework. This opens up the possibility to

perform eﬃcient simulations of three-dimensional Stokes ﬂow with arbitrary periodicity (D= 3,2,1,0),

using boundary integral and potential methods.

The SE method has also been adapted to related models, such as Brinkman ﬂow [25], and Stokesian

dynamics [26]. To facilitate the inclusion of Brownian motion, [27] proposed the Positively Split Ewald

(PSE) method for the Rotne–Prager–Yamakawa (RPY) tensor of Stokesian dynamics, thus ensuring that

both the short-range and long-range parts are symmetric positive deﬁnite. The PSE method has found

much use in Brownian dynamics [28,29,30,4,7,5,11,31,2,8,9].

In this paper, the aperiodic directions are assumed to be free and extend to inﬁnity. If they are instead

bounded, one possibility is to explicitly discretize the boundary, and enforce boundary conditions on it, as

done e.g. in [32], without modifying the underlying SE method. Another option is the general geometry

Ewald-like method (GGEM) [33,34], which uses the same split into a short-range and long-range part, but

treats the long-range part in real space using a mesh-based solution method. GGEM can in principle handle

nontrivial boundary conditions, but involves an expensive correction solve if high accuracy is desired. A

more recent solution is presented by [35], in which a Chebyshev method is used for the boundary value

problem in the aperiodic direction; this method has been demonstrated for electrostatics but is expected

to generalize to Stokes ﬂow. Related work on Ewald-type methods for reduced periodicity in electrostatics

includes [36,37,38]. For half-space Stokes problems, the methods by [39] or [40] can be used. The SE

method has also been implemented for two-dimensional Stokes ﬂow [41].

An alternative to Ewald-like summation methods is the fast multipole method (FMM) [42,43,44,45],

which in general achieves O(N) complexity. Unlike Ewald-type methods, which reach their highest eﬃciency

for fully periodic problems, the FMM is most eﬃcient and most natural to formulate in the fully aperiodic

(free-space) setting. Nevertheless, the FMM has also been generalized to arbitrary periodicity [46,47,48].

An advantage of the FMM is that it is spatially adaptive, while the SE method requires a uniform spatial

grid due to the FFT. Thus, the FMM will typically be faster for highly nonuniform point distributions

(especially in free space), while the SE method may be faster for uniform distributions, as shown e.g. by

[23,24]. Yet another alternative method is found in [49], in which the long-range interaction is represented

by auxiliary sources. The principle is the same in both two and three dimensions, but it has not been

demonstrated that this method would be competitive with Ewald-type methods in three dimensions.

The contribution of the current paper is, as mentioned above, to complete and unify the previous work

on the SE method for the kernels of Stokes ﬂow (stokeslet, stresslet, rotlet), by adding the singly and doubly

periodic cases, and treat all periodic cases within the same framework. An important part of the uniﬁed

SE method is the modiﬁed kernels that are used to treat the singular integration in cases with reduced

periodicity, based on an idea by [50]. In this paper, we improve the convergence of the modiﬁed kernels

used by [23] in the free-space case, and derive new modiﬁed kernels for the singly and doubly periodic

cases. We also derive a new improved truncation error estimate for the stokeslet and stresslet, valid in all

periodic cases, based on techniques from [23]. Analytical formulas useful for validation are derived in the

singly and doubly periodic cases, completing the formulas previously derived by [20] for the doubly periodic

stokeslet. The SE method presented in this paper furthermore uses the polynomial Kaiser–Bessel (PKB)

window introduced by [51] and [24], and the adaptive Fourier transform (AFT) introduced by [52] for the

singly and doubly periodic cases.

The paper is organized as follows. In section 2, the Ewald splits and Ewald sums are given for the three

kernels. In section 3, the fast method to compute the Fourier-space Ewald sums is presented, i.e. the SE

method; this section also describes the modiﬁed kernels, PKB window, and AFT. In section 4, we give error

estimates for the SE method and describe an automated procedure to select the parameters of the method

given an error tolerance; this procedure is also tested by numerical examples. In section 5, more numerical

results follow, regarding the pointwise error, computational time and complexity, and the window function.

Finally, conclusions are drawn in section 6. In the appendices, we derive analytical formulas for validation

in Appendix A and Appendix B, and the improved truncation error estimate for the stokeslet and stresslet

in Appendix C.

2. Ewald summation for Stokes ﬂow

We consider three fundamental solutions of Stokes ﬂow, namely the stokeslet S, rotlet Ω, and stresslet T,

which are tensorial kernels given by

Sjl(r) = δjl

|r|+rjrl

|r|3,Ωjl(r) = jlm

|r|3, Tjlm(r) = −6rjrlrm

|r|5,(5)

where j, l, m ∈ {1,2,3}. Here, δjl denotes the Kronecker delta and jlm denotes the Levi-Civita symbol;

Einstein’s summation convention is used, with repeated indices implicitly summed over {1,2,3}. Given a

point force 8πfacting on the ﬂuid at x0, the Stokesian velocity ﬁeld is given by uj(x) = Sjl(x−x0)fl, the

vorticity is given by ωj(x) = 2Ωjl(x−x0)fl, and the stress ﬁeld is given by σjl(x) = Tjlm(x−x0)fm[53,

ch. 2.2]. In boundary integral equations, the kernels S,Ωand Tare multiplied by diﬀerent source quantities

which all give rise to velocity ﬁelds, namely

j(x) = Sjl(x−x0)fS

l,(6)

uΩ

j(x) = Ωjl(x−x0)fΩ

l,(7)

j(x) = Tjlm(x−x0)fT

lm,(8)

where notably the stresslet source fT

lm has two indices. Commonly, the stresslet source is of the form

lm =qlνm, where qand νare vectors, and this will be assumed in this paper.

In the following, we will make frequent use of the fact that the fundamental solutions of Stokes ﬂow can

be related to the harmonic (Laplace) Green’s function H(r) = 1/|r|or the biharmonic Green’s function

B(r) = |r|, via [54, Appendix E, p. 113][23][55]

Sjl(r)=(δjl∇2− ∇j∇l)B(r) =: KS

jlB(r),(9)

Ωjl(r) = −jlm∇mH(r) =: KΩ

jlH(r),(10)

Tjlm(r) = [(δjl∇m+δmj ∇l+δlm∇j)∇2−2∇j∇l∇m]B(r) =: KT

jlmB(r).(11)

We have here introduced a linear diﬀerential operator Kfor each kernel. To be able to treat all kernels

together where appropriate, we now introduce some special notation following [23]; we will write

u(x) = G(x−x0)·f,(12)

where the kernel Gmay be the stokeslet S, rotlet Ωor stresslet T, and the notation G·fwill be understood

to mean one of

Sjlfl,Ωjlfl, Tjlmflm,(13)

depending on the actual kernel. (Note that G·fis always a vector quantity.)

We are interested in eﬃciently evaluating periodic potentials of the form

uDP(x) =

n=1 X

p∈PDP

G(x−xn+p)·f(xn),(14)

generated by Nsources of strengths f(xn) located at points xninside a box B= [0, L1)×[0, L2)×[0, L3),

called the primary cell. The set PDPof periodic images depends on the number of periodic directions D

and is given by

PDP=









{(¯p1L1,¯p2L2,¯p3L3) : ¯pi∈Z},if D= 3 (triply periodic),

{(¯p1L1,¯p2L2,0) : ¯pi∈Z},if D= 2 (doubly periodic),

{(¯p1L1,0,0) : ¯pi∈Z},if D= 1 (singly periodic),

{(0,0,0)},if D= 0 (free space).

(15)

Note that the ﬁrst Dcoordinate directions are the periodic ones; the remaining 3 −Ddirections are called

the free directions. (The notation DPintroduced here is to be read as “D-periodic”.) The sum (14) may

be evaluated either at one of the source locations xm(m= 1, . . . , N), in which case the term corresponding

to m=n,p=0is omitted, or at other arbitrary locations. In the former case, we write

uDP?(xm) =

n=1

p∈PDP

G(xm−xn+p)·f(xn),(16)

where the star above the second sum denotes precisely that the term p=0is omitted when n=m.

Since the kernels (5) decay slowly as |r|→∞(the stokeslet decays like 1/|r|, and the rotlet and stresslet

like 1/|r|2), the periodic sums (14) and (16) are absolutely convergent only in the case D= 0 (for all

kernels), and D= 1 for the rotlet and stresslet (but not the stokeslet). In all other cases, the sums are either

conditionally convergent (i.e. their values depend on the order of summation) or even divergent. Special

care must then be taken when interpreting (14) and (16); an overview of these considerations is given by

[56]. Ewald summation corresponds to a spherical order of summation.

Ewald summation is based on the idea to split the periodic sum into two parts, where the ﬁrst part

converges fast in real space (the “real-space part”), and the other part converges fast in Fourier space (the

“Fourier-space part”). Ewald decompositions for the kernels considered here are given in section 2.1. To aid

the reader, we then present the Ewald sums ﬁrst in the triply periodic setting (D= 3) in section 2.2, and

then for arbitrary periodicity (D= 3,2,1,0) in section 2.3.

2.1. Ewald decompositions

In general, there are two diﬀerent but equivalent ways to derive an Ewald decomposition, called screening

and splitting [23]. The real-space part is easier to derive using the splitting approach, while the Fourier-

space part is easier using the screening approach. In this paper, we use the screening approach since the

Fourier-space part is our main focus.

In the screening approach, the kernel is convolved with a screening function γ(r;ξ), which is a smooth

function with the property RR3γ(r;ξ) dr= 1. The positive parameter ξis called the Ewald decomposition

parameter, and controls the decay of the resulting decomposition

G=G∗(δ−γ) + G∗γ, (17)

where GR:= G∗(δ−γ) is the real-space part and GF:= G∗γis the Fourier-space part; here, δis the Dirac

delta distribution, and ∗denotes convolution. In our case where G=KA, with Kbeing a linear diﬀerential

operator and Aeither the harmonic Hor biharmonic B, cf. (9)–(11), the decomposition can be written as

G=GR+GFwith

GR=K[A∗(δ−γ)] =: KAR,(18)

GF=K[A∗γ] =: KAF.(19)

Thus, we start by writing down the Ewald decompositions for the harmonic and biharmonic kernels, and

may then apply the appropriate diﬀerential operator Kto get the decompositions for the stokeslet, rotlet,

and stresslet.

For the harmonic, we select the classical Ewald screening function

γE(r;ξ) = ξ3π−3/2e−ξ2|r|2

bγE(k;ξ)=e−|k|2/(2ξ)2,(20)

where the hat denotes the Fourier transform

f(k) = F{f}(k) = ZR3

f(r)e−ik·rdr.(21)

This leads to the decomposition originally derived by Ewald [57], namely

HR(r;ξ)=[H∗(δ−γE)](r;ξ) = erfc(ξ|r|)

|r|,(22)

HF(r;ξ)=[H∗γE](r;ξ) = erf(ξ|r|)

|r|,(23)

where erf(·) is the error function and erfc(·)=1−erf(·) is the complementary error function. Here, the

singular behaviour of 1/|r|is contained in HR, while the long-range behaviour is contained in HF. Note

that HRdecays fast as |r|→∞in real space, while HFis smooth and therefore its Fourier transform decays

fast. By the convolution theorem, the Fourier transform (in the distributional sense) of HFis

HF(k;ξ) = b

H(k)bγE(k;ξ) = 4π

|k|2e−|k|2/(2ξ)2,(24)

文档加载中……请稍候！
如果长时间未打开，您也可以点击刷新试试。

下载文档到电脑，查找使用更方便

10 玖币 0人已下载

立即下载

摘要：

FastEwaldsummationforStokesowwitharbitraryperiodicityJoarBagge,Anna-KarinTornbergKTHMathematics,LinneFLOWCentre/Swedishe-ScienceResearchCentre,RoyalInstituteofTechnology,SE-10044Stockholm,SwedenAbstractAfastandspectrallyaccurateEwaldsummationmethodfortheevaluationofstokeslet,stressletandrotletpote...

展开>> 收起<<

Fast Ewald summation for Stokes flow with arbitrary periodicity.pdf

共54页,预览5页

还剩页未读，继续阅读

声明：本站为文档C2C交易模式，即用户上传的文档直接被用户下载，本站只是中间服务平台，本站所有文档下载所得的收益归上传人(含作者)所有。玖贝云文库仅提供信息存储空间，仅对用户上传内容的表现方式做保护处理，对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私，请立即通知玖贝云文库，我们立即给予删除！

Fast Ewald summation for Stokes flow with arbitrary periodicity

相关推荐

开通VIP享超值会员特权

作者详情

相关内容

热门标签

举报选择: