1. Introduction
Jean Baptiste Fourier offered an incredibly useful tool to describe intricate phenomena in the dual frequency space. Fourier analysis is the basis of a huge list of experimental techniques, either based on spectral analyses (arising from either thermal or chaotic signals) or on the system response to an oscillatory force with monochromatic frequency. The latter approach allows tracking the phase delay between the external forcing
$F(t)$
and the system’s response
$V(t)$
whose relation generally involves a memory integral
$V(t) = \int \chi (t-t^{\prime }) F(t^{\prime }) {\rm d}t^{\prime }$
. The convolution theorem elegantly expresses this relation in the dual space using (complex) phasor quantities
$V(\omega) = \chi (\omega) F(\omega)$
to unveil the elastic (in-phase,
$\text{Re}[\chi]$
) and dissipative response (out-of-phase,
$\text{Im}[\chi]$
) encoded in the complex susceptibility
$\chi (\omega)$
. While this route is generally applicable (provided linear response), the susceptibility function gathers all possible (elastic and inelastic) mechanisms able to propagate forces to the system. This introduces a serious handicap in experimental analyses of soft matter (in general) and bio-matter (in particular), because these objects naturally exist in liquid environments. Thus, to unveil the intrinsic susceptibility of soft objects in liquid, one first necessarily needs to evaluate the fluid contribution and then to extract it from the total one
$\chi (\omega)$
measured in experiments. Yet, the analytical evaluation of the fluid contribution to
$\chi (\omega)$
under oscillatory regimes is a complicated hydrodynamic problem (particularly when boundaries are present), which has been studied by various generations of leading experts over at least the last 50 years (see Mazur, Bedeaux & Mazur Reference Mazur, Bedeaux and Mazur1974
b; Pozrikidis Reference Pozrikidis1989; Felderhof Reference Felderhof2005; Simha, Mo & Morrison Reference Simha, Mo and Morrison2018; Fouxon & Leshansky Reference Fouxon and Leshansky2018; Fouxon, Rubinstein & Leshansky Reference Fouxon, Rubinstein and Leshansky2023). This handicap is currently hampering and significantly blurring the interpretation of many techniques designed to study the mechanical response of soft biological samples. Some examples are quartz crystal microbalance (QCM) (Vázquez-Quesada et al. Reference Vázquez-Quesada, Schofield, Tsortos, Mateos-Gil, Milioni, Gizeli and Delgado-Buscalioni2020; Schofield & Delgado-Buscalioni Reference Schofield and Delgado-Buscalioni2021), vibrational spectra from atomic force microscopy (AFM) (de Beer et al. Reference de Beer, van den Ende and Mugele2008), AC magnetometry (HAC) and AC susceptibility (ACS) analysis (Zhong et al. Reference Zhong, Lucht, Hankiewicz, Schilling and Ludwig2019), electrochemical impedance spectroscopy (e.g. of DNA) in liquids, and other techniques where nonlinear oscillatory hydrodynamics create the driving forces, such as in ultrasound manipulation (Marmottant & Hilgenfeldt Reference Marmottant and Hilgenfeldt2003; Bruus Reference Bruus2012)and viscoelastic response of deformable objects (e.g. vesicles) under oscillatory electric fields (Dimova et al. Reference Dimova, Bezlyepkina, Jordö, Knorr, Riske, Staykova, Vlahovska, Yamamoto, Yang and Lipowsky2009; Vlahovska Reference Vlahovska2015), amongst others. The relevance of hydrodynamics in some of these experimental devices has attracted the attention of theoretical and simulation groups, which have proven the dominant contribution of the solute-induced perturbative flows in QCM (see Meléndez-Schofield et al. Reference Meléndez, Vázquez-Quesada and Delgado-Buscalioni2020; Fouxon et al. Reference Fouxon, Rubinstein and Leshansky2023; Delgado-Buscalioni Reference Delgado-Buscalioni2024), some of them being able to reproduce experimental results without fitting parameters (Vázquez-Quesada et al. Reference Vázquez-Quesada, Schofield, Tsortos, Mateos-Gil, Milioni, Gizeli and Delgado-Buscalioni2020; Delgado-Buscalioni Reference Delgado-Buscalioni2024). Other theoretical analyses have focused on liquid atomic force microscopy (AFM) (Clarke et al. Reference Clarke, Cox, Williams and Jensen2005; Kabarowski & Khair Reference Kabarowski and Khair2020), ultrasound (Bruus Reference Bruus2012; Balboa Usabiaga & Delgado-Buscalioni Reference Balboa Usabiaga and Delgado-Buscalioni2013) and optical spectroscopy of pinned elastic membranes (Janeš et al. Reference Janeš, Schmidt, Blackwell, Seifert and Smith2019), to mention a few. In this contribution, we tackle this problem by presenting a fast and accurate (spectral) scheme to solve arbitrary oscillatory flows, more precisely, to calculate the velocity field arising from arbitrary force field distributions in the fluid.
The Green function for oscillatory force propagation in the free space was analytically derived a long time ago (Pozrikidis Reference Pozrikidis1989; Kim & Karilla Reference Kim and Karilla1991) and Mazur et al. (Reference Mazur, Bedeaux and Mazur1974b
) derived exact analytical extensions for the self-mobility of an oscillating rigid sphere of radius
$a$
in incompressible fluids and in compressible fluids (Mazur et al. Reference Mazur, Bedeaux and Mazur1974a
). However, the boundary presence significantly complicates both analytical approaches (Pozrikidis Reference Pozrikidis1989; Fouxon et al. Reference Fouxon, Rubinstein and Leshansky2023) and numerical solutions, due to the formation of small boundary layers where the flow is governed by exponential corrections to the dominant potential flow outside (Fouxon et al. Reference Fouxon, Rubinstein and Leshansky2023). Felderhof (Reference Felderhof2005) (see also Felderhof Reference Felderhof2009,Reference Felderhof2012) offered important contributions to this problem by analytically deriving an exact expression for the in-plane Fourier transform of the semi-bounded Green (phasor) tensor. Felderhof used this result to derive the self-mobility of a point particle at some distance
$h$
from the wall, by including the flow reflected from the wall in the so-called reaction field tensor. The height-dependent particle mobility results in a diagonal tensor, with different in-plane and normal-to-plane components. Around a boundary oscillating at angular frequency
$\omega$
, momentum spreads exponentially up to a characteristic distance
$\delta =(2\nu /\omega)^{1/2}$
called the penetration length (here,
$\nu$
is the kinetic viscosity of the fluid). Simha et al. (Reference Simha, Mo and Morrison2018) discussed Felderhof’s point-particle approximation and warned about the existence of two different non-dimensional length scales
$\delta /a$
and
$h/a$
(here,
$a$
is the particle radius) whose smallness might induce the invalidation of the point-particle limit. Indeed, as a finite particle (large enough
$a/\delta$
) approaches the wall, the lack of symmetry in the normal-to-wall direction induces translational–rotational couplings which alter the self-mobility tensor of a spherical particle (Fouxon et al. Reference Fouxon, Rubinstein and Leshansky2023). These are absent in free-space (Mazur et al. Reference Mazur, Bedeaux and Mazur1974b
). Notably, the Green tensor of a slippery surface coincides with that of a no-slip wall if the particle is far enough from the boundary (Fouxon & Leshansky Reference Fouxon and Leshansky2018). However, close to the wall, the no-slip boundary condition introduces higher-order reflections, ultimately leading to a peculiar oscillatory lubrication regime to be further studied (see Fouxon et al. Reference Fouxon, Rubinstein and Leshansky2023; Delgado-Buscalioni Reference Delgado-Buscalioni2024). Nevertheless, the wall presence induces a
$z$
-dependent flow which changes in distances of order
$\delta$
, and thus, for small enough values of
$h/a$
, even the monopole approximation to the self-mobility (averaging the fluid velocity over the particle volume) differs from the point-particle limit, as we show here. As explained by Pozrikidis (Reference Pozrikidis1989), the Green function in real space cannot be solved analytically, which precludes a closed solution for the self-mobility and makes even more difficult the analysis of two-particle mobility tensor (which depends on both particles’ heights). Despite this limitation, several theoretical contributions have analysed finite-particles in oscillatory flows under certain asymptotic limits (small particle, long distance from wall, low frequency) (Fouxon & Leshansky Reference Fouxon and Leshansky2018; Fouxon et al. Reference Fouxon, Rubinstein and Leshansky2023; Zhang et al. Reference Zhang, Bertin, Essink, Zhang, Fares, Shen, Bickel, Salez and Maali2023). Recently, Leshansky’s group used the Lorentz reciprocal theorem to derive relevant analytical relations in several asymptotic regimes, and applied their analysis to hydrodynamic effects in the QCM response of suspended and adsorbed particles (Fouxon et al. Reference Fouxon, Rubinstein and Leshansky2023; Leshansky et al. Reference Leshansky, Rubinstein, Fouxon, Johannsmann, Sadowska and Adamczyk2024), which were also explored by some of us (Schofield & Delgado-Buscalioni Reference Schofield and Delgado-Buscalioni2021; Delgado-Buscalioni Reference Delgado-Buscalioni2024). In the temporal realm, and even in the free-space case, the difficulty to treat memory kernels makes intractable the use of the Green function formalism to solve the unsteady Stokes equation. Most standard numerical schemes have solved the oscillatory Stokes equations in the real space and time domains. QCM studies have been reported using lattice Boltzmann with fixed obstacles (Gillissen et al. Reference Gillissen, Jackman, Tabaei and Cho2018; Gopalakrishna et al. Reference Gopalakrishna, Langhoff, Brenner and Johannsmann2021) and using the immersed boundary method for flexible moving structures (Vázquez-Quesada et al. Reference Vázquez-Quesada, Schofield, Tsortos, Mateos-Gil, Milioni, Gizeli and Delgado-Buscalioni2020). Working in the time domain requires waiting long transient times
$L^2/\nu$
to solve vorticity diffusion over the box and reach the homogeneous oscillatory response. The idea of working in the frequency domain is then quite appropiate and has been used for QCM applications (Gopalakrishna et al. Reference Gopalakrishna, Langhoff, Brenner and Johannsmann2021). In the spatial domain though, standard schemes (based on finite boxes closed in vertical direction) have to pay a large computational cost arising from physical and practical reasons. First, one needs to solve the exponentially small boundary layer (size
$\delta$
) very close to the wall or the moving boundary (while adapted grids are useful (Zhang et al. Reference Zhang, Bertin, Essink, Zhang, Fares, Shen, Bickel, Salez and Maali2023), they restrict the generality of the scheme). Second, flow disturbances are propagated at long distances, both in the in-plane direction (as
$1/L_{\parallel }^2$
as we shall later see) and in the normal-to-wall direction (exponentially decaying as
$e^{-z/\delta }$
). These two facts conspire to complicate the calculations, imposing the use of small grids to resolve the flow near obstacles and long boxes to avoid finite-size effects. These facts significantly increase the computational time required to accurately solve unsteady Stokes dynamics near walls and obstacles, particularly when many parameters need to be explored (Reference Vázquez-Quesada, Schofield, Tsortos, Mateos-Gil, Milioni, Gizeli and Delgado-BuscalioniVázquez-Quesada et al. 2020). Practical applications (e.g. QCM) introduce an even worse problem: the small size
$a$
of the adsorbed analytes (proteins, viruses) studied in experiments introduce a severe limiting factor for standard computational methods based on closed boxes (see Gillissen et al. Reference Gillissen, Jackman, Tabaei and Cho2018; Gopalakrishna et al. Reference Gopalakrishna, Langhoff, Brenner and Johannsmann2021; Vázquez-Quesada et al. Reference Vázquez-Quesada, Schofield, Tsortos, Mateos-Gil, Milioni, Gizeli and Delgado-Buscalioni2020; Zhang et al. Reference Zhang, Bertin, Essink, Zhang, Fares, Shen, Bickel, Salez and Maali2023). To resolve the obstacle, one requires a small mesh size
$h_{{grid}} \sim 0.1 \, a$
, but also one needs tall boxes of at least
$L_z \gt 4\delta$
to avoid finite-size hydrodynamic effects. For proteins,
$a \sim 5\, \text{nm}$
while
$\delta \sim 100 \,\text{nm}$
so both conditions cannot not be rationally satisfied (
$L_z/h_{{grid}} \sim 10^3$
). Ideally, one would like to restrict the mesh to the interest region (a layer of height
$H \sim a$
adjacent to the wall) and treat this computational domain as an open portion of a semi-infinite domain. In the present work, we develop an alternative scheme which addresses both spatial and temporal complications. We solve the oscillatory Stokes equation using a spectral scheme for doubly periodic (DP) domains in free or semi-confined spaces. In the periodic
$(xy)$
plane, equations are discretised in the reciprocal space using the Fourier transform, while for the non-periodic
$z$
direction, we use Chebyshev collocation points. Borrowing ideas from the work by Hashemi et al. (Reference Hashemi, Peláez, Natesh, Sprinkle, Maxian, Gan and Donev2023) for DP solvers of the steady Stokes equation, the system is decomposed into an inner ‘open’ domain, where forces are applied, coupled to the outer domain, which is solved analytically using plane-wave expansions. Concerning the temporal question, we avoid using time-stepping by working in the frequency domain using complex phasor fields. While this approach is limited to linear response, it allows for a huge reduction in computational time compared with time-stepping. Extensions to deal with nonlinear contributions are suggested in § 7. Here, we focus on the description of the fluid solver and its validation against available exact analytical relations for the mobility and the Green tensor in the reciprocal space. Coupling hydrodynamic fields with the dynamics of viscoelastic matter will be the subject of another contribution. We will present the computational scheme in § 3, then in § 6, we test its accuracy by comparison with the exact Felderhof analytical relations for the Fourier transform of the Green tensor (Felderhof Reference Felderhof2005) (collected in Appendix H) and analyse the mutual mobility, the self mobility of spherical objects in free and bounded space, analysing the validity of the reaction-field point-particle approximation (Felderhof Reference Felderhof2005) and the scaling of hydrodynamic finite size effects due to image interactions in the DP system.
2. System set-up and oscillatory Stokes equations
The system geometry is illustrated in figure 1. It is periodic in
$x$
and
$y$
directions with lengths given by
$L_x$
and
$L_y$
, respectively (we use square boxes and denote
$L_{\parallel }=L_x=L_y$
). The domain of interest (inner domain) where numerical discretisation is applied (in the reciprocal space for
$xy$
and in the real space for
$z$
) is a layer of height
$2H$
which might be either open at
$z=\pm H$
(open slit) or present a rigid wall at
$z=-H$
.
Figure 1. Diagram of the problem geometry. A fluid-filled three-dimensional domain is defined with periodic boundary conditions in the
$x$
and
$y$
directions (of lengths
$L_x,L_y$
, respectively). The domain is open in the positive
$z$
direction with an open boundary at
$z=H$
and walled at
$z=-H$
. The wall oscillates at frequency
$\omega$
and with amplitude
${\textbf v}_{{wall}}$
. An arbitrary oscillatory fluid force density with phasor
$\boldsymbol {f}(\boldsymbol {r})$
and frequency
$\omega$
is contained inside the domain, i.e.
$\boldsymbol {f}(x, y, |z|\geqslant H) = 0$
.
We solve the unsteady Stokes equations, corresponding to the zero Reynolds limit of the momentum equation for an incompressible fluid,
\begin{eqnarray} \rho \partial _t \textbf {v} - \eta \,\nabla ^2 \textbf {v} = -\nabla p + \textbf {f}, \end{eqnarray}
\begin{eqnarray} \nabla \cdot \textbf {v} =0, \end{eqnarray}
where
$\eta$
is the fluid viscosity,
$\rho$
is the fluid density,
$p(\textbf {r},t)$
is the hydrodynamic pressure,
$\textbf {v}(\textbf {r},t)$
is the fluid velocity field and
$\boldsymbol {f}(\textbf {r},t)$
is a time dependent force distribution. Oscillatory flows are forced by a monochromatic oscillation of the sources and/or the wall’s velocity at frequency
$\omega$
. Introducing the Fourier transform in time,
$v(r,\omega)=\int _{-\infty }^{\infty } \exp [-\mathrm{i} \omega t] v(\textbf {r},t) {\rm d}t$
, one gets a set of equations for the (complex) phasor fields. To simplify the notation (unless otherwise indicated), we use
$\textbf {v}(\textbf {r})$
to indicate the phasor field and its dependence on
$\omega$
is not explicitly indicated. Fourier transform in time leads to
\begin{eqnarray} \eta [\alpha ^2 - \nabla ^2] \textbf {v} = -\nabla p + \textbf {f}, \end{eqnarray}
\begin{eqnarray} \nabla \cdot \textbf {v} =0, \end{eqnarray}
where
$\alpha \equiv (1 -\mathrm{i})\, \delta ^{-1}$
and
$ \delta \equiv [2 \eta /(\rho \omega)]^{1/2}$
is the penetration length determining the width of the fluid layer where momentum diffusion balances the oscillatory momentum.
The external (oscillatory) force distribution
$\textbf {f}(\textbf {r})$
is restricted to the inner domain
$z\in (-H,H)$
so that
$\textbf {f} (|z|\gt H)=0$
. The open slit configuration (or ‘free space’) corresponds to an open layer connected to the exterior domain, solved by plane-wave expansions. The semi-bounded domain corresponds to a no-slip rigid boundary at
$z=-H$
and open conditions at
$z=H$
. In the semi-bounded set-up, the bottom wall is allowed to move so that
$\textbf {v}(z=-H)={\textbf v}_{{wall}}$
and the flow vanishes far away from the boundary,
$\textbf {v} \rightarrow 0$
and
$p \rightarrow 0$
for
$|z|\rightarrow \infty$
.
3. Numerical method
Following Hashemi et al. (Reference Hashemi, Peláez, Natesh, Sprinkle, Maxian, Gan and Donev2023), we decompose the solution of (2.3) and (2.2) as
\begin{eqnarray} \textbf {v}= {\textbf v}_{dp} + {\textbf v}_{c}, \end{eqnarray}
\begin{eqnarray} p = {\textbf v}_{dp} + p_{c}, \end{eqnarray}
where
${\textbf v}_{dp}$
and
$p_{dp}$
are the velocity and pressure fields in free space using doubly periodic boundaries, and
${\textbf v}_{c}$
and
$p_{c}$
are correction fields imposed to satisfy the desired boundary conditions (BCs) at
$z=\pm H$
.
The free-space solution for the doubly periodic system satisfies
\begin{eqnarray} \eta \left [\alpha ^2 - \nabla ^2\right] {\textbf v}_{dp} = -\nabla p_{dp} +{\textbf f}\;\; \mathrm{at} \;z\in [-H,H], \end{eqnarray}
\begin{eqnarray} \eta \left [\alpha ^2 - \nabla ^2\right] {\textbf v}_{dp} = -\nabla p_{dp} \;\;\;\;\; \mathrm{at} \;z\notin [-H,H], \end{eqnarray}
\begin{eqnarray} \nabla \cdot {\textbf v}_{dp} =0, \end{eqnarray}
\begin{eqnarray} {\textbf v}_{dp}(|z|=\infty) = p_{dp}(|z|=\infty)=0. \end{eqnarray}
This sub-problem is solved for a given set of sources
${\textbf f}({\textbf r})$
. We use a spectral method which applies a Fourier transform in the
$xy$
plane, leading to independent equations for a discrete set of wavevectors
${\textbf k}=2\pi \,(n_x/L_x,n_y/L_y)$
(with
$n_x$
and
$n_y$
integers) and Chebyshev transforms in
$z\in [-H,H]$
. The solution for the outer domain
$|z|\gt H$
, expressed as plane waves, is matched with the inner solution, providing the boundary conditions at the borders of the inner domain
$z=\pm H$
. Once the free-space flow
${\textbf v}_{dp}$
and
$p_{dp}$
are obtained, the correction flow can be derived analytically. As all sources have already been imposed in the free-space flow, the correction flow satisfies a homogeneous equation,
\begin{eqnarray} \eta \left [\alpha ^2 - \nabla ^2\right] {\textbf v}_{c} = -\nabla p_{c}\;\; \mathrm{at} \;z\in [-H,H], \end{eqnarray}
\begin{eqnarray} \nabla \cdot {\textbf v}_{c} =0, \end{eqnarray}
and, in this case, we need to impose the boundary conditions for the velocity
${\textbf v}$
at
$z=\pm H$
. In the most general case, the present method allows the imposition of boundary conditions in the reciprocal space
$\boldsymbol {\mathcal{B}}[{\textbf v}]({\textbf k})={\textbf B}_{wall}({\textbf k})$
, where
$\boldsymbol {\mathcal{B}}$
is a linear operator. The flow decomposition in (3.1) leads to
\begin{equation} \boldsymbol {\mathcal{B}}_{\pm H}[{\textbf v}_{c}] = {\textbf B}_{wall} -\boldsymbol {\mathcal{B}}_{\pm H}[{\textbf v}_{dp}]. \end{equation}
Here, we consider
$\mathcal{\boldsymbol {B}}={\textbf 1}$
and
${\textbf B}_{wall}({\textbf k})={\textbf 0}$
for any
${\textbf k}\neq {\textbf 0}$
, while
${\textbf B}_{wall}(k=0)={\textbf v}_{wall}$
corresponds to an oscillating no-slip rigid wall. Each wavevector is solved independently due to linearity; however, we see that the
${\textbf k}={\textbf 0}$
mode leads to a different set of equations providing the overall current at each
$z$
due to the net force
$\int {\textbf f}({\textbf r})\, {\rm d}x {\rm d}y$
.
3.1. Free-space solver
We first deal with the free-space solver, so that
${\textbf v}_{c}=0$
and
$p_{c}=0$
.
3.1.1. Solution for the pressure
Taking the divergence on the momentum equation,
\begin{equation} \nabla ^2 p_{dp}=\nabla \cdot {{\textbf f}}. \end{equation}
As stated, the in-plane dependence is worked out in the reciprocal space. Decomposing the coordinates as
${\textbf r} = {\textbf s} + z\, {\hat {{\textbf z}}}$
, the Fourier transform in
$xy$
of the velocity field is
\begin{align} {{\textbf v}_{dp}}({\textbf k},z) = \int e^{\mathrm{i}{\textbf k}\cdot {\textbf s}} {\textbf v}_{dp}({\textbf s},z) {\rm d}^2 {\textbf s}. \end{align}
Introducing this transformation into the pressure equation leads to a Laplace equation inside the domain,
\begin{equation} \left (\partial _z^2 - k^2\right) {p_{dp}} = (\mathrm{i} {\textbf k}, \partial _z) \cdot {{\textbf f}}\quad \mathrm{for}\,z \in [-H,H]. \end{equation}
We denote
${\textbf k}=(k_x,k_y)$
and
$k=(k_x^2+k_y^2)^{1/2}$
with a discrete set of wavevectors in the doubly periodic box.
Outside the domain, source terms vanish and
\begin{equation} (\partial _z^2 - k^2) p_{dp} = 0 \quad \mathrm{for}\,z \notin [-H,H], \end{equation}
whose solution is of the form
$p_{dp}({\textbf k},z) = C_1 \exp [-k z]+C_2 \exp [k z]$
. Using
$p=0$
for
$|z|\rightarrow \infty$
,
\begin{eqnarray} p_{dp}({\textbf k},z) = C_2 \exp [k z]\quad \mathrm{for}\;z\leqslant -H,\end{eqnarray}
\begin{eqnarray} p_{dp}({\textbf k},z) = C_1 \exp [-k z] \quad \mathrm{for}\;z\geqslant H .\end{eqnarray}
At the borders of the domain, this implies
$(\partial _z \pm k) p({\textbf k},z=\pm H)=0$
, which should be satisfied for each wavevector
${\textbf k}$
. Hence, the Chebyshev solver for the pressure
$p_{dp}$
in the inner domain solves
\begin{eqnarray} \left (\partial _z^2 - k^2\right) p_{dp} = (\mathrm{i} {\textbf k}+\partial _z) \cdot {{\textbf f}}\quad \mathrm{for}\,z \in [-H,H], \end{eqnarray}
\begin{eqnarray} (\partial _z \pm k) p_{dp}({\textbf k},z=\pm H) =0. \end{eqnarray}
The boundary value problem (BVP) (3.16) subject to boundary conditions in (3.17) is solved for each wavenumber using a Chebyshev discretisation in
$z$
(a similar approach to that used to solve the Poisson equation (Maxian et al. Reference Maxian, Peláez, Greengard and Donev2021)). Once the BVP is solved, the constants
$C_1$
and
$C_2$
depend on the wavevector
${\textbf k}$
, and can be obtained by matching the outwards and inwards pressure fields at the borders. From (3.14), and (3.15),
\begin{eqnarray} C_1({\textbf k}) = p_{dp}({\textbf k},z=H) \, e^{k H}\quad \mathrm{at}\; z=H, \end{eqnarray}
\begin{eqnarray} C_2({\textbf k}) = p_{dp} ({\textbf k},z=-H) e^{-k H} \quad \mathrm{at}\; z=-H .\end{eqnarray}
3.1.2. Solution for the horizontal (parallel) velocity component
Once
$p$
is known, we can solve for the velocity. In the case of the oscillatory Stokes equation, we need to take into account the time derivative, which, in phasor quantities, equals to
$\eta \alpha ^2 {\textbf v}_{dp}$
. Inside the domain,
\begin{equation} \eta \left (\partial ^2_z-\tilde {k}^2\right) {\textbf v}_{dp} = (\mathrm{i} {\textbf k} + \hat {{\textbf z}} \partial _z)\, p_{dp} - {\textbf f}\quad \mathrm{for}\,z\in [-H,H],\end{equation}
where we have defined the complex wavenumber
$\tilde {{\textbf k}}$
whose square is
$\tilde {k}^2 \equiv k^2 + \alpha ^2$
. Outside the domain,
\begin{equation} \eta \left (\partial ^2_z-\tilde {k}^2\right) {\textbf v}_{dp} = (\mathrm{i} {\textbf k} + \hat {{\textbf z}} \partial _z)\, p_{dp} \quad \mathrm{for}\,z\notin [-H,H]. \end{equation}
Equation (3.21) needs to be solved for both
${\textbf v}_{dp}^{\parallel }$
(velocity parallel to the plane) and for
${\textbf v}_{dp}^{\perp }$
(velocity perpendicular to the plane). For the parallel velocity inside the domain,
\begin{equation} \eta \left (\partial ^2_z-\tilde {k}^2\right) {\textbf v}^{\parallel }_{dp} = \mathrm{i} {\textbf k} \, p_{dp} \quad \mathrm{for}\,z\notin [-H,H],\end{equation}
where the pressure field is given by (3.14) and (3.15). Hence,
\begin{eqnarray} \left (\partial ^2_z-\tilde {k}^2\right) {\textbf v}^{\parallel }_{dp} = \mathrm{i} \frac {{\textbf k} C_1}{\eta } e^{-kz} \quad \mathrm{for}\,z\geqslant H, \end{eqnarray}
\begin{eqnarray} \left (\partial ^2_z-\tilde {k}^2\right) {\textbf v}^{\parallel }_{dp} = \mathrm{i} \frac {{\textbf k} C_2}{\eta } e^{kz} \quad \mathrm{for}\,z\leqslant -H,\end{eqnarray}
whose solution has the form (see Appendix A)
\begin{eqnarray} {\textbf v}^{\parallel }_{dp} ={\textbf A}_{\parallel,1} e^{-\tilde {k} z} -\frac {\mathrm{i} {\textbf k} C_1}{\eta \alpha ^2} e^{-kz} \quad \mathrm{for}\,z\geqslant H,\end{eqnarray}
\begin{eqnarray} {\textbf v}^{\parallel }_{dp} = {\textbf A}_{\parallel,2} e^{\tilde {k} z} -\frac {\mathrm{i} {\textbf k} C_2}{\eta \alpha ^2} e^{kz} \quad \mathrm{for}\,z\leqslant -H. \end{eqnarray}
After some algebra, one derives the Dirichlet-to-Neumann maps providing the boundary conditions for the inner domain,
\begin{eqnarray} \partial _z {\textbf v}^{\parallel }_{dp} + \tilde {k} {\textbf v}^{\parallel }_{dp} = -\frac {\mathrm{i} {\textbf k} (\tilde {k}-k)}{\eta \alpha ^2} p_{dp} \quad \mathrm{for}\;z =H, \end{eqnarray}
\begin{eqnarray} \partial _z {\textbf v}^{\parallel }_{dp} - \tilde {k} {\textbf v}^{\parallel }_{dp} = \frac {\mathrm{i} {\textbf k} (\tilde {k}-k)}{\eta \alpha ^2} p_{dp} \quad \mathrm{for}\;z = -H, \end{eqnarray}
where we recall that
$\tilde {k}-k$
is a complex number.
3.1.3. Solution for the vertical (perpendicular) velocity component
Outside the domain, in the
$\hat {{\textbf z}}$
direction,
\begin{equation} \eta \left (\partial ^2_z-\tilde {k}^2\right) v^{\perp }_{dp} = \partial _z \, p_{dp} \quad \mathrm{for}\,z\notin [-H,H],\end{equation}
with the outward pressure fields given by (3.14) and (3.15), leading to
\begin{eqnarray} \left (\partial ^2_z-\tilde {k}^2\right) v^{\perp }_{dp} = -\frac {k C_1}{\eta } \exp [-k z] \quad \mathrm{for}\;z\geqslant H, \end{eqnarray}
\begin{eqnarray} \left (\partial ^2_z-\tilde {k}^2\right) v^{\perp }_{dp} = \frac {k C_2}{\eta } \exp [k z]\quad \mathrm{for}\;z\leqslant -H, \end{eqnarray}
whose solution (Appendix A),
\begin{eqnarray} v^{\perp }_{dp} = A_{\perp,1} e^{-\tilde {k} z} + \frac {k C_1}{\eta \alpha ^2} e^{-kz} \quad \mathrm{for}\,z\geqslant H,\end{eqnarray}
\begin{eqnarray} v^{\perp }_{dp} = A_{\perp,2} e^{\tilde {k} z} -\frac {k C_2}{\eta \alpha ^2} e^{kz} \quad \mathrm{for}\,z\leqslant -H, \end{eqnarray}
has to satisfy the following boundary conditions for the inner domain,
\begin{eqnarray} \partial _z v^{\perp }_{dp} +\tilde {k} v^{\perp }_{dp} = \frac {\left(\tilde {k} -k \right)k}{\eta \alpha ^2} p_{dp}\quad \mathrm{for}\,z=H, \end{eqnarray}
\begin{eqnarray} \partial _z v^{\perp }_{dp} -\tilde {k} v^{\perp }_{dp} = \frac {\left(\tilde {k} -k \right)k}{\eta \alpha ^2} p_{dp}\quad \mathrm{for}\,z=-H. \end{eqnarray}
The BVP problems (3.25), (3.26) with BCs (3.27) and (3.28), and (3.32), (3.33) with BCs (3.34) and (3.35) uniquely determine the free-space solution which is solved using a spectral Chebyshev solver (Maxian et al. Reference Maxian, Peláez, Greengard and Donev2021; Hashemi et al. Reference Hashemi, Peláez, Natesh, Sprinkle, Maxian, Gan and Donev2023) for each wavevector
${\textbf k}$
used in the in-plain Fourier decomposition in an embarrasingly parallel way.
3.2. Correction flow
At this stage, we have a numerical solution for
${\textbf v}_{dp}$
in
$z\in [-H,H]$
based on a Chebyshev expansion (i.e. we have the corresponding coefficients) and the solution of the pressure, from the Laplace equation, in
$z\in [-H,H]$
. The solution of
${\textbf v}_{dp}$
in the inner domain and, in particular, at
$z=\pm H$
is needed to calculate
$\boldsymbol {\mathcal{B}}_{\pm }[{\textbf v}_{dp}]$
and impose the BC for the correction flow
${\textbf v}_c$
. Each discrete wavevector allowed in the doubly periodic system leads to an independent boundary condition which determines the spectral solution of the confined flow, i.e.
\begin{equation} \boldsymbol {\mathcal{B}}_{\pm H}[{\textbf v}_{c}]({\textbf k})={\textbf B}_{wall,\pm H}({\textbf k}) - \boldsymbol {\mathcal{B}}_{\pm H}[{\textbf v}_{dp}]({\textbf k}). \end{equation}
3.3. Boundary condition at the bottom wall
According to (3.36), boundary conditions are imposed by adding a correction flow to the free-space solution. Here, we consider a wall at the bottom of the domain
$z=-H$
(figure 1) so the correction flow applies to the domain
$z\in [-H, \infty]$
and satisfies
\begin{eqnarray} \eta \left (\partial _z^2 - \tilde {k}^2\right) {\textbf v}_{c} = \left (\mathrm{i} {\textbf k} +\hat {{\textbf z}} \partial _z\right) \, p_{c}\quad \mathrm{for}\; z\in [-H,\infty], \end{eqnarray}
\begin{eqnarray} \left ( \mathrm{i} {\textbf k} + \hat {{\textbf z}} \partial _z \right) {\textbf v}_{c} =0, \end{eqnarray}
\begin{eqnarray} \mathcal{B}_{-H}[{\textbf v}_{c}]({\textbf k},-H) = {\textbf B}_{wall}({\textbf k}) - \mathcal{B}_{-H}[{\textbf v}_{dp}]({\textbf k}, -H). \end{eqnarray}
The BC in (3.39) is written in its most general form to indicate that the scheme permits introduction of a complex velocity pattern at the surface, which may involve a linear operator
$\boldsymbol {\mathcal{B}}$
acting on
${\textbf v}$
and some pattern in the Fourier space
${\textbf B}_{wall}({\textbf k})$
. This allows for a broad range of situations (e.g. partial slippage on patterned surfaces). Here, we will apply an homogeneous velocity distribution at the wall, corresponding to
${\textbf B}_{wall}({\textbf k})=0$
for
${\textbf k} \ne 0$
and
${\textbf B}_{wall}(k=0)={\textbf v}_{wall}$
. A steady no-slip wall (
$\mathcal{B}_{-H}=1$
and
${\textbf B}_{wall} =0$
) corresponds to
${\textbf v}_{c}({\textbf k}, -H) = -{\textbf v}_{dp}({\textbf k},-H)$
.
We first deal with the correction pressure which satisfies,
$(\partial _z^2 -k^2) p_c =0$
, so
$p_c = B e^{-kz}$
and
$\partial _z p_c =- k B e^{-kz}$
. For the correction velocity,
${\textbf v}_c={\textbf v}_c^{\parallel }+ v_c^{\perp } \,\hat {{\textbf z}}$
, we have
\begin{eqnarray} (\partial _z^2 -\tilde {k}^2) {\textbf v}_c^{\parallel } = \frac {\mathrm{i} {\textbf k} B}{\eta } e^{-kz}, \end{eqnarray}
\begin{eqnarray} \left(\partial _z^2 -\tilde {k}^2 \right) v_c^{\perp } = -\frac {k B}{\eta } e^{-kz}, \end{eqnarray}
with solution
\begin{eqnarray} {\textbf v}_c^{\parallel } = {\textbf D}_{\parallel } e^{-\tilde {k} z} - \frac {\mathrm{i} {\textbf k} B}{\eta \alpha ^2} e^{-kz}, \end{eqnarray}
\begin{eqnarray} v_c^{\perp } = D _{\perp } e^{-\tilde {k} z} + \frac {k B}{\eta \alpha ^2} e^{-kz} .\end{eqnarray}
To ensure continuity,
\begin{equation} \tilde {k} D_{\perp }=\mathrm{i} {\textbf k}\cdot {\textbf D}_{\parallel }, \end{equation}
so we can define the vector of coefficients
${\textbf D}$
in terms of
${\textbf D}_{\parallel }$
(which lives in the plane),
\begin{equation} {\textbf D} = \left [{\textbf I}_{\parallel } + \left (\frac {\mathrm{i} k}{\tilde {k}}\right) \hat {{\textbf z}} \otimes \hat {{\textbf k}}\right] \cdot {\textbf D}_{\parallel } \end{equation}
and also
\begin{equation} {\textbf B} \equiv (-\mathrm{i} \hat {{\textbf k}} + \hat {{\textbf z}}) B. \end{equation}
Hence, the correction pressure and velocity can be written as
\begin{eqnarray} p_c = B e^{-kz}, \end{eqnarray}
\begin{eqnarray} {\textbf v}_c = {\textbf D} e^{-\tilde {k} z} + {\textbf B} \frac {k}{\eta \alpha ^2} e^{-kz}. \end{eqnarray}
At the boundary, we need to impose (3.39), which for a steady homogeneous wall (
${\textbf v}_{{wall}}({\textbf k})=0$
for
$k\gt 0$
) implies
${\textbf v}_c(z=-H) = -{\textbf v}_{dp}(z=-H)$
. Note that at this point, we already derived
${\textbf v}_{dp}$
in the form of an expansion of Chebyshev polynomials. Using (3.48), at
$z=-H$
, this leads to the following system of equations:
\begin{equation} e^{\tilde {k} H} {\textbf D}_{\parallel }({\textbf k}) + \frac {k}{\eta \alpha ^2} e^{k H} \, {\textbf B}({\textbf k}) = {\textbf v}_B\equiv -{\textbf v}_{dp}({\textbf k},-H),\end{equation}
where
${\textbf v}_B \equiv {\textbf v}_{c}({\textbf k},-H)$
is the correction velocity at the boundary, and
${\textbf B}$
and
${\textbf D}$
satisfy (3.46) and (3.45) or (3.44). This leaves us with a solvable system of equations with three equations and three unknowns (
$\boldsymbol {D}^\parallel$
and
$B$
). In more explicit form,
\begin{align} \begin{bmatrix}\exp ((\tilde {k}-k)H) & 0 & \dfrac {-ik_x}{\eta \alpha ^2}\\[9pt] 0&\exp ((\tilde {k}-k)H)&\dfrac {-ik_y}{\eta \alpha ^2}\\[9pt] \frac {ik_x}{\tilde {k}}\exp ((\tilde {k}-k)H)&\dfrac {ik_y}{\tilde {k}}\exp ((\tilde {k}-k)H)&\dfrac {k}{\eta \alpha ^2} \end{bmatrix} \begin{bmatrix}D_x\\[9pt] D_y\\[5pt] B\end{bmatrix} = \begin{bmatrix}v_B^{x}\\[5pt] v_B^y\\[5pt] v_B^z\end{bmatrix} \exp (-kH) \end{align}
with solution
\begin{align} D_x=\frac {\exp (-\tilde {k}H)}{k(k-\tilde {k})}\left [v_B^x(k^2_y-k\tilde {k}) -k_xk_yv_B^y -i\tilde {k}k_xv_B^z\right], \end{align}
\begin{align} D_y=\frac {\exp (-\tilde {k}H)}{k(k-\tilde {k})}\left [-v_B^xk_xk_y+v_B^y\left (k_x^2-k\tilde {k}\right)-v_B^zi\tilde {k}k_y\right], \end{align}
\begin{align} B = \frac {\eta \alpha ^2\exp (-kH)}{k(k-\tilde {k})}\left (v_B^xik_x+v_B^yik_y-v_B^z\tilde {k}\right). \end{align}
4. The
$ k={\textbf 0}$
mode
We finally solve the overall current at each
$z$
arising from the net force imposed on the system
$z\in [-H,H]$
. Let us first deal with the pressure.
4.1. Pressure
For the bottom wall we consider here, taking
${\textbf k} =0$
in (3.12) leads to
\begin{equation} \partial _z^2 p({\textbf k}=0,z) = \partial _z f_z(k=0,z), \end{equation}
where we have introduced the perpendicular component of the net force at each
$z$
,
\begin{align} f_z(k=0,z)=\int _{x,y} {\textbf f}(x,y,z) \cdot \hat {{\textbf z}}\,{\rm d}x\,{\rm d}y. \end{align}
The pressure is constant for
$z\lt -H$
and integrating (4.1) from
$z=-(H+\epsilon)$
to
$z=-H$
leads to the force balance
$\partial _z p({\textbf k},-H) = f_z(-H)$
. Integrating (4.1) from
$-H$
to
$z$
and using the force balance at
$z=-H$
reveals the balance between pressure gradient and net-force balance plane-by-plane,
$\partial _z p({\textbf k}=0,z) = f_z(z)$
. Integrating once again,
\begin{align} p({\textbf k}=0,z) = \int _{-H}^z f_{z} ({\textbf k}=0,z) {\rm d}z +C. \end{align}
We choose
$p({\textbf k}=0,-H)=0$
so that
$C=0$
. Hence, the BVP for the pressure involves (4.1) and the BCs,
\begin{eqnarray} p({\textbf k}=0,-H) = 0, \end{eqnarray}
\begin{eqnarray} p({\textbf k}=0,H) = \int _{-H}^H f_{z} ({\textbf k}=0,z) {\rm d}z, \end{eqnarray}
which is computed from the Chebyshev coefficients of
$f_z$
(see Appendix C).
4.2. Perpendicular velocity
For the perpendicular velocity
$k=0$
field,
\begin{equation} \eta \left (\partial _z^2 v^{\perp } - \alpha ^2 v^{\perp} \right) = \partial _z p -f_z =0, \end{equation}
whose solution
$v^{\perp }= a e^{\alpha z} + b e^{-\alpha z}$
should vanish at
$z=\infty$
. Imposing a moving wall,
$v^{\perp }(z=-H)={\textbf v}_{{wall}} \cdot \hat {{\textbf z}}$
, leads to
\begin{equation} v^{\perp }({\textbf k}=0,z)= {\textbf v}_{{wall}} \cdot \hat {{\textbf z}}\, e^{-\alpha (z+H)} .\end{equation}
In the present study, we consider
${\textbf v}_{{wall}} \cdot \hat {{\textbf z}}=0$
, which implies
$v^{\perp }({\textbf k}=0,z)=0$
and thus no-net momentum in the vertical direction.
4.3. Parallel velocity
In the case of the net parallel velocity, we have the following set:
\begin{eqnarray} \eta \left (\partial _z^2 {\textbf v}^{\parallel } - \alpha ^2 {\textbf v}^{\parallel} \right) = -{\textbf f}^{\parallel }({\textbf k}=0,z), \end{eqnarray}
\begin{eqnarray} {\textbf v}^{\parallel }({\textbf k}=0,z=-H) ={\textbf v}_{{wall}}^{\parallel }, \end{eqnarray}
\begin{eqnarray} {\textbf v}^{\parallel }({\textbf k},z=\infty) =0. \end{eqnarray}
We note that the homogeneous part of (4.8) (
${\textbf f}=0$
) contributes to
${\textbf v}^{\parallel }$
with the Stokes basal flow
${\textbf v}_{wall}^{\parallel } e^{-\alpha (z+H)}$
, which guarantees the BC (4.9). The inhomogeneous part of (4.8) (solving for
${\textbf f}\neq 0$
and
$v_{wall}=0$
) is a BVP which can be solved using the Green function formalism, as explained in Appendix B. Adding the homogeneous and inhomogeneous parts of the solution leads to
\begin{equation} {\textbf v}^{\parallel }({\textbf k}=0,z) = -\frac {1}{2\eta \alpha } \int _{-H}^H \left (e^{-\alpha (z+z'+2H)} -e^{-\alpha |z-z'|)}\right) {\textbf f}^{\parallel }(z') dz' + {\textbf v}_{{wall}}^{\parallel }\,e^{- \alpha (H+z)}. \end{equation}
Equation (4.11) provides the velocity of the zero mode at any
$z$
, and it can be evaluated at the position of the Chebyshev grid, if needed. However, using fast Chebyshev transforms and solving the BVP (i.e. imposing the BC at
$z=H$
) is a faster computational route. For this reason, we now derive the BC for (4.8) at
$z=H$
. We have implemented a Neumann condition at
$z=H$
; however, it is in principle possible to impose a Dirichlet BC at
$z=H$
. Both are illustrated below.
4.3.1. Dirichlet condition for the inhomogeneous part of (4.8) at
$z=H$
The two boundary conditions for the BVP problem of the inhomogeneous part of (4.8) (i.e.
${\textbf f}\neq 0$
and
$v_{wall}=0$
) are
\begin{eqnarray} {\textbf v}^{\parallel }({\textbf k}=0,z=-H) = 0, \end{eqnarray}
\begin{eqnarray} {\textbf v}^{\parallel }({\textbf k}=0,z=H) = -\frac {e^{-\alpha H}}{2\eta \alpha } \int _{-H}^H \left (e^{-\alpha (z'+2H)} -e^{\alpha z'}\right) {\textbf f}^{\parallel }(z') {\rm d}z' .\end{eqnarray}
The second relation can be simplified to
\begin{equation} {\textbf v}^{\parallel }({\textbf k}=0,z=H) = \frac {e^{-2\alpha H}}{\eta \alpha } \int _{-H}^H \sinh [\alpha (z'+H)] {\textbf f}^{\parallel }(z') {\rm d}z'. \end{equation}
For
$|\alpha z |\lt \lt 1$
, one gets, to first order in
$\alpha$
,
${\textbf v}^{\parallel }= \alpha \int _{-H}^H (z'+H){\textbf f}(z') {\rm d}z'$
, whose real part coincides with the solution for the Stokes limit. Finally, adding the homogeneous solution of (4.8),
$v_{wall} e^{-\alpha (z+H)}$
to impose
${\textbf v}^{\parallel }({\textbf k}=0,z=-H) = v_{wall}$
, leads to the complete result,
\begin{equation} {\textbf v}^{\parallel }({\textbf k}=0,z=H) = \frac {e^{-2\alpha H}}{\eta \alpha } \int _{-H}^H \sinh [\alpha (z'+H)] {\textbf f}^{\parallel }(z') {\rm d}z' + {\textbf v}_{{wall}}^{\parallel }\,e^{-2\alpha H}. \end{equation}
4.3.2. Neumann condition for the inhomogeneous part of (4.8) at
$z=H$
.
A Dirichlet BC at
$z=-H$
and a Neumann-type BC at
$z=H$
is obtained by taking the derivative in relation (4.11)
\begin{eqnarray} {\textbf v}^{\parallel }({\textbf k}=0,z=-H) = {\textbf v}_{{wall}}^{\parallel }, \end{eqnarray}
\begin{eqnarray} \partial _z {\textbf v}^{\parallel }({\textbf k}=0,z=H) = -\frac {e^{-2 \alpha H}}{\eta } \int _{-H}^H \sinh [\alpha (z'+H)] {\textbf f}^{\parallel }(z') {\rm d}z' - \alpha v_{wall}\,e^{-2\alpha H} .\end{eqnarray}
Note that for
$\alpha \rightarrow 0$
, we get
$\partial _z {\textbf v}^{\parallel }(H)=0$
as in the steady Stokes limit (Hashemi et al. Reference Hashemi, Peláez, Natesh, Sprinkle, Maxian, Gan and Donev2023). In our scheme, we have implemented the
$k=0$
BCs (4.16) and (4.17). The costly part (in either (4.15) or (4.17)) resides in the evaluation of the integral
$\int _{-H}^H \sinh [\alpha (z'+H)]\, {\textbf f}^{\parallel }(z') {\rm d}z'$
, which is solved in the preprocessing step using a monomial expansion (see Appendix C).
5. Implementation
For the sake of clarity, we now recapitulate the different steps of the scheme.
-
(i) Given a force distribution
${\textbf f}({\textbf r})$
(normally resulting from the spreading of blob forces
${\textbf f}({\textbf r})= \sum _i S({\textbf r}-{\textbf q}_i) {\textbf F}_i$
, see Appendix D), we first take the Fourier–Chebyshev transform (FCT) to obtain
${\textbf f}({\textbf k},z)$
and use Chebyshev differentiation to obtain
$\partial _z {\textbf f}({\textbf k},z)$
. The scheme extensively uses fast Fourier–Chebyshev transforms (Trefethen Reference Trefethen2000) and spectral operations, which have been thoroughly described in several works (Maxian et al. Reference Maxian, Peláez, Greengard and Donev2021; Hashemi et al. Reference Hashemi, Peláez, Natesh, Sprinkle, Maxian, Gan and Donev2023; Pérez Peláez Reference Peláez-Pérez2022). -
(ii) Then, solve the BVP for the free space pressure
$p_{dp}$
using (3.16) with BC (3.17) and reconstruct the pressure field
$p_{dp}({\textbf k},\pm H)$
using (3.14), (3.15) and (3.18). Use this to evaluate the BVP for the velocity
${\textbf v}_{dp}$
in (3.25), (3.26) with BCs (3.27), (3.28) and (3.32), (3.33) with BCs (3.34), (3.35). At this point, we have the solution for the free space
$p_{dp}$
and
${\textbf v}_{dp}$
in (3.3), (3.4). -
(iii) To evaluate the correction flow, first evaluate
${\textbf v}_{dp}({\textbf k},-H)$
from (3.25), (3.26), (3.32) and (3.33), then use (3.47), (3.48) to obtain
$p_{c}$
and
${\textbf v}_c$
, where
${\textbf D}$
and
$B$
are calculated from the system in (3.49). -
(iv) Solve the
$k=0$
mode for
$p(0,z)$
using Chebyshev integration on (4.5). Then, solve the BVP given by (4.8) and BCs (4.16) and (4.17) for the
$k=0$
parallel velocity. The integral in (4.17) is solved using the method explained in Appendix C. Finally, set
${\textbf v}_{\perp }(0,z)=0$
. -
(v) For all
${\textbf k}$
, take a one-dimensional (1-D) fast Fourier transform (FFT) in the z direction to compute the Chebyshev coefficients of the correction and zero-mode solutions, and add them to doubly periodic solution to obtain the three-dimensional (3-D) spectral representation of
$p$
and
${\textbf v}$
. -
(vi) At this point, performing a back transform to real space (iFCT) leads to
$p({\textbf r})$
and
${\textbf v}({\textbf r})$
, as sketched in figure 2(a). Figure 2(b) illustrates the velocity fields obtained from a Gaussian force distribution oscillating in normal and tangential directions. Numerical discretisation is only carried out in the inner domain
$|z| \leqslant H$
, and the flow outside is connected by the plane wave analytic continuation. As indicated in figure 2(a), one can also interpolate blob velocities
${\textbf u}=\mathcal{J}_{{\textbf q}}[{\textbf v}]=\int {\textbf v}({\textbf r}) S({\textbf r}-{\textbf q}) {\rm d}{\textbf r}$
to evaluate mobility fields. Also, spectral quantities can be directly used to evaluate local stress fields
$\boldsymbol {\sigma } = -p{\textbf I} +\eta \nabla (\boldsymbol {v}+ \boldsymbol {v}^{T})$
and drag forces
$\mathcal{J}_{{\textbf q}}[\nabla \cdot \boldsymbol {\sigma }]$
on the blobs.
The CUDA implementation of the present solver makes extensive use of the UAMMD library (Peláez et al. 2025) and will be available at the code website https://uammd.readthedocs.io/en/latest/.
6. Validation and results
6.1. Validation tests
Figure 2. (a) Sketch of the method used to measure particle mobility. An immersed boundary (IB) kernel (typically Gaussian, as described in Appendix D) located is employed to spread a force
$\mathbf{F}$
onto the fluid grid (with cell size
$h_{\text{grid}}$
) resulting in a force distribution
${\textbf f}({\textbf r}^{\prime })$
. The resulting flow
${\textbf v}({\textbf r})$
is interpolated at the location
${\textbf q}$
using a similar kernel, with an associated hydrodynamic radius
$a$
. Both are related by the mobility tensor
${\textbf u}={\mathcal{\boldsymbol M}}_{{\textbf q},{\textbf q}^{\prime }} {\textbf F}$
. (b) Streamlines of oscillatory flows induced by normal (left) and tangential forces (right) created by a Gaussian distribution with radius
$a=0.1\delta$
at a height
$h=2.5\,a$
and penetration length
$\delta =[2 \eta /(\rho \omega)]^{-1/2}$
. The streamlines are taken from the real part of the velocity field, i.e. at zero phase lag. We illustrate the location of the inner domain
$z\in [-H,H]$
(where forces are confined and the velocity field is numerically solved) and the outer domain
$z\gt H$
, where the field is analytically continued by a plane wave expansion. In this example,
$H=0.5\,\delta$
. Colors in the streamline indicate the flow speed in units of
$F/(6 \pi \eta \delta)$
, with
$F$
the force amplitude.
We have performed two types of validations of the scheme, presented in the following subsections. The first test compares the Green function in the reciprocal space (in-plane Fourier transform) with analytical expressions of Felderhof (Reference Felderhof2005). The second calculates mutual and self-mobilities of spherical particles (actually immersed boundary ‘blobs’, defined next) and compare with Mazur–Bedeaux relation (for free space) and Felderhof’s point-particle approximation (for semibounded domains). In these tests, we use immersed boundary (IB) kernels to spread localised forces and to interpolate the velocity field at different locations, to evaluate ‘particle’ mobilities. Most of the results presented have been obtained with Gaussian kernels (Hashemi et al. Reference Hashemi, Peláez, Natesh, Sprinkle, Maxian, Gan and Donev2023; Pérez Peláez Reference Peláez-Pérez2022), but we also compare with the three-point Peskin (Reference Peskin2002) kernels (see Appendix D for details). As shown in previous works (Pérez Usabiaga et al. Reference Balboa Usabiaga, Delgado-Buscalioni, Griffith and Donev2014; Peláez Reference Peláez-Pérez2022), the IB kernels have compact support and an associated hydrodynamic radius
$a$
, which can be mapped to the radius of a ‘rigid’ spherical particle, providing minimal particle-models (Usabiaga et al. Reference Balboa Usabiaga, Delgado-Buscalioni, Griffith and Donev2014; Pérez Peláez Reference Peláez-Pérez2022). Despite their simplicity, IB kernels provide surprisingly accurate ‘finite particle’ effects, being able to even recover the nonlinear drag forces of hard-spheres up to large particle Reynolds number (Balboa et al. Reference Balboa, Bell, Delgado-Buscalioni, Donev, Fai, Griffith and Peskin2012), or the extent of (nonlinear) ultrasound forces in compressible oscillatory flows (Bruus Reference Bruus2012; Balboa Usabiaga & Delgado-Buscalioni Reference Balboa Usabiaga and Delgado-Buscalioni2013), also, they can be generalised to provide a solenoidal spreading field Bao et al. (Reference Bao, Donev, Griffith, McQueen and Peskin2017). Figure 2(a) illustrates the use of these kernels in mobility evaluations. A force distribution
${\textbf f}({\textbf r}^{\prime })= {\textbf F} S({\textbf r}^{\prime } - {\textbf q}_0)$
is prescribed by an IB kernel
$S({\textbf r})$
centred at
${\textbf q}_0$
and activated to propagate momentum to the fluid. The resulting velocity field
${\textbf v}({\textbf r})$
is interpolated at another location
${\textbf q}$
, providing a ‘particle’ velocity
${\textbf u} = \int S({\textbf r} -{\textbf q}) {\textbf v}({\textbf r}) {\rm d}{\textbf r}$
. The nomenclature used in the following sections is indicated in figure 2(a) (for instance,
$h$
is the vertical component of the kernel at
${\textbf q}$
). Note that interpolations are carried out in the discrete setting by summing over a discrete set of fluid cells (i.e. in the real space) (Usabiaga et al. Reference Balboa Usabiaga, Delgado-Buscalioni, Griffith and Donev2014).
Discretisation errors are presented in Appendix E against the mesh size in
$x$
and
$y$
directions,
$h_{{grid}}$
(recall that Chebyshev points in
$z$
are not regularly spaced, being more compact near
$z=\pm H$
). We advance the main outcome of the error analysis, based on the ratio
$N_{\delta }\equiv \delta /h_{{grid}}$
, which is roughly the number of mesh points used to resolve the fluid penetration length
$\delta =(2\nu /\omega)^{1/2}$
(again, recall that Chebyshev is not regularly spaced). Notably, relative errors become smaller than
$0.05$
for just
$N_{\delta }\gt 4$
, which reveals the excellent accuracy of spectral schemes (in a 3-D staggered regular-mesh, a similar accuracy requires approximately
$N_{\delta }\gt 12$
(Delgado-Buscalioni Reference Delgado-Buscalioni2024)). In Appendix F, we also show that the choice of
$H$
is arbitrary (it does not introduce errors) as long as forces and sources are confined in the inner domain (
$0\lt z\lt H$
).
6.2. Green tensor
The solution of (2.3) and (2.2) can be expressed using the Green tensor of the problem
$\boldsymbol {G}({\textbf r},{\textbf r}^{\prime })$
, as
\begin{equation} {\textbf v}({\textbf r}) = \int \boldsymbol {G}({\textbf r},{\textbf r}^{\prime }) {\textbf f}({\textbf r}^{\prime }) {\rm d}^3{\textbf r}^{\prime }. \end{equation}
Using the isotropy of the problem in the
$xy$
plane (and decomposing
${\textbf r} = {\textbf s} + z\, \hat {{\textbf z}}$
), Felderhof (Reference Felderhof2005) derived the in-plane Fourier transform of the Green tensor for the semi-bounded oscillatory Stokes problem,
\begin{equation} \boldsymbol {{\mathcal{G}}}({\textbf k},z,z^{\prime })= \int \boldsymbol {G}({\textbf s}-{\textbf s}^{\prime },z,z^{\prime }) e^{\mathrm{i} {\textbf k}\cdot ({\textbf s}-{\textbf s}^{\prime })} {\rm d}{\textbf s}^{\prime }, \end{equation}
which we shall deploy to validate the numerical scheme. The tensor in the reciprocal space
$\boldsymbol {{\mathcal{G}}}$
is however not analytically invertible to the real space (see Pozrikidis Reference Pozrikidis1989; Fouxon & Leshansky Reference Fouxon and Leshansky2018). Recently, Fouxon and Leshanski have derived sound approximations to
${\textbf G}$
to study small particles or arbitrarily large spheres far enough from the surface (Fouxon & Leshansky Reference Fouxon and Leshansky2018; Fouxon et al. Reference Fouxon, Rubinstein and Leshansky2023). In any case, to validate the present spectral scheme, Felderhof’s exact expressions
$\boldsymbol {{\mathcal{G}}}({\textbf k},z,z^{\prime })$
(recapitulated in Appendix H) are quite convenient. To that end, we indicate
${\textbf r}={\textbf s}+z \,\hat {{\textbf z}}$
and apply a force distribution
${\textbf f}({\textbf r}^{\prime })$
localised at
${\textbf q}_0={\textbf s}_0 + {h}_0 \hat {{\textbf z}}$
. We will use a smeared delta function with a spatial dependence separable as a product of kernels in orthogonal coordinates:
$S_\parallel ({\textbf s})$
indicates the in-plane
$xy$
kernel product and
$S_z(z)$
the kernel used in the vertical direction (see Appendix D for details),
\begin{align} S({\textbf r}-{\textbf q})= S_{\parallel }({\textbf s}-{\textbf s}_{0})\times S_z(z-{h}_0). \end{align}
We apply a force density
${\textbf F}$
at
${\textbf q}_0={\textbf s}_0 + {h}_0 \hat {{\textbf z}}$
, resulting in a velocity field,
\begin{align} v_{\alpha }({\textbf r};{h}_0) = \int G_{\alpha \beta }({\textbf s}-{\textbf s}^{\prime },z,z^{\prime }) F_{\beta }\; [S_{\parallel }({\textbf s}^{\prime }-{\textbf s}_{0})\times S_z(z^{\prime }-{h}_0)] {\rm d}{\textbf s}^{\prime } {\rm d}z^{\prime } \end{align}
(where the
${h}_0$
dependence in the velocity indicates that it depends on the point of application of the force). Its plane-Fourier transform,
\begin{align} \int e^{-\mathrm{i} {\textbf k} \cdot {\textbf s}} v_{\alpha }({\textbf q};{h}_0) {\rm d}{\textbf s} & = \int dz^{\prime } \int {\rm d}{\textbf s}^{\prime } \int {\rm d}{\textbf s} e^{-\mathrm{i} {\textbf k} \cdot {\textbf s}} G_{\alpha \beta }({\textbf s}-{\textbf s}^{\prime },z,z^{\prime }) F_{\beta }\; [S_{\parallel }({\textbf s}^{\prime }-{\textbf s}_0) \nonumber\\& \quad \times\, S_z(z^{\prime } -{h}_0)], \end{align}
can be simplified using the isotropy in the
$xy$
plane, to introduce the variable
$\xi ^{\prime } ={\textbf s}^{\prime }- {\textbf s}_0$
with
${\rm d}\xi ^{\prime } = {\rm d}{\textbf s}^{\prime }$
,
\begin{align} v_{\alpha }({\textbf k},z;\,{h}_0) = \!\int\! {\rm d}z^{\prime } \!\int\! {\rm d}{\textbf s}^{\prime } e^{-\mathrm{i} {\textbf k} \cdot {\textbf s}^{\prime }} \!\int\! {\rm d}\boldsymbol {\xi }^{\prime } e^{-\mathrm{i} {\textbf k} \cdot \boldsymbol {\xi }} G_{\alpha \beta }(\boldsymbol {\xi }^{\prime },z,z^{\prime }) F_{\beta }\; [S_{\parallel }( \boldsymbol {\xi }^{\prime }-{\textbf s}_0)\times S_z( z^{\prime }-{h}_0)],\end{align}
to get
\begin{eqnarray} v_{\alpha }({\textbf k},z\,;{h}_0) =e^{-\mathrm{i} {\textbf k} \cdot {\textbf s}_0} S_{\parallel }({\textbf k}) \bar {\mathcal{G}}_{\alpha,\beta }({\textbf k},z\,;{h}_0), \end{eqnarray}
\begin{eqnarray} \bar {\mathcal{G}}_{\alpha,\beta }({\textbf k},z\,;{h}_0) &\equiv & \int {\rm d}z^{\prime } S_z(z^{\prime }-{h}_0) \mathcal{G}_{\alpha \beta }({\textbf k},z,z^{\prime })\, F_{\beta }. \end{eqnarray}
Proving the velocity field excited by unit forces in different
$\hat {\boldsymbol {\gamma }}$
orthogonal directions (i.e. using
$F_{\beta } = \delta _{\beta,\gamma }$
), we reconstruct the Green tensor components as
$v_\alpha ^{(\gamma)} \propto \bar {\mathcal{G}}_{\alpha,\gamma }$
, which is analytically derivable from (6.8) using Gaussian kernels for
$S_{\parallel }({\textbf s})$
and
$S_z(z)$
and Felderhof’s expressions for
$\mathcal{G}_{\alpha \beta }$
(see Appendix H).
Figure 3. Comparison between the analytical result for the Gaussian kernel smeared Green function
$\boldsymbol {{\mathcal{G}}}_S({\textbf k},z,\mathrm{h})$
((6.2), lines), obtained by Felderhof (Reference Felderhof2005Reference Felderhof2009, Reference Felderhof2012) (Appendix H) and numerical results (symbols). Results correspond to the flow at different positions in the vertical coordinate
$z/\delta$
with an imposed unit force at
$\mathrm{h}=0.25\delta$
spread to the fluid with a Gaussian kernel of hydrodyamic radius
$a=0.1\delta$
. Numerical results were obtained using a grid of
$\{150, 150, 65\}$
left with a cell size
$h_{{grid}} = \delta / 30$
.
Such a comparison is illustrated in figure 3, where a perfect agreement is found between the numerical results (circles) and the analytical relations (lines) for all the nine tensor components. The case corresponds to a force located at
${\textbf s}_0=0$
and
$h_0=0.1\delta$
spread to the fluid using a Gaussian kernel of hydrodynamic radius
$a=0.1\delta$
in a box with
$L_x=L_y=5\delta$
, and
$H=0.7\delta$
discretised with
${\textbf N}= \{150, 150, 65\}$
points.
6.3. Two-particle mobility
Figure 4. (a) Real and (b) imaginary parts of the components
$\mathcal{M}^{\alpha \beta }_{12}$
of the mobility. Particle
$1$
is located at
${\textbf q}_1=(0,0,3a)$
and particle
$2$
at
${\textbf q}_2=(x,0,6a)$
. According to the reciprocal relations (6.11), dots and lines of the same colour should be equal (and they match to machine precision). The box size is
$L = [48a, 48a, 16a]$
and the penetration length is
$\delta = 10a$
.
The previous validation of the Green tensor allows to ensure that the two-particle mobility is correctly captured by the numerical scheme. We now evaluate the mutual mobility of two ‘blob particles’ defined by their kernel distributions
$S(r)$
. The relevant issue for our purpose here is to evaluate the mutual mobility, or equivalently, the smeared Green tensor. To that end, we apply a force density field
${\textbf f}({\textbf r}^{\prime }) = {\textbf F} S({\textbf r}^{\prime }-{\textbf q}_2)$
associated with a unit force
$|{\textbf F}|=1$
centred at
${\textbf q}_2$
and evaluate the resulting velocity of a blob located at
${\textbf q}_1$
. Denoting
${\textbf u}_1$
the blob’s velocity, this operation results in
\begin{equation} {\textbf u}_1 = \int\! \int S({\textbf r}-{\textbf q}_1) \boldsymbol {G}({\textbf r},{\textbf r}^{\prime }) S({\textbf r}^{\prime }-{\textbf q}_2) {\textbf F}_2 {\rm d}^3 {\textbf r} \, {\rm d}^3{\textbf r}^{\prime } = \boldsymbol {{\mathcal{M}}}_{12} {\textbf F}_2 .\end{equation}
The two-particle mobility tensor is thus simply
\begin{equation} \boldsymbol {{\mathcal{M}}}_{12} \equiv {\textbf G}_S({\textbf q}_1,{\textbf q}_2) \equiv \int\! \int S({\textbf r}-{\textbf q}_1) \boldsymbol {G}({\textbf r},{\textbf r}^{\prime }) S({\textbf r}^{\prime }-{\textbf q}_2) {\rm d}^3{\textbf r} \,{\rm d}^3{\textbf r}^{\prime }, \end{equation}
where we have introduced the smeared Green tensor (field)
${\textbf G}_S({\textbf q}_1,{\textbf q}_2)$
associated with the kernel
$S({\textbf r})$
. We present results for the two-particle mobility using Gaussian blobs. In doing so, we also validate the numerical scheme against the reciprocal relations (Felderhof Reference Felderhof2009; Fouxon & Leshansky Reference Fouxon and Leshansky2018) which arise from the Green tensor,
$G_{\alpha \beta } ({\textbf r},{\textbf r}_0) =G_{\beta \alpha } ({\textbf r}_0,{\textbf r})$
. In terms of mobilities,
\begin{equation} \mathcal{M}^{\alpha \beta }_{12}= \mathcal{M}^{\beta \alpha }_{21}. \end{equation}
We have verified that our numerical scheme satisfies relations (6.11) to machine precision. This test is shown in figure 4, where the mutual mobility
$\mathcal{M}_{12}^{\alpha \beta }(x,\;{\mathrm h}_1,{\mathrm h}_2)$
is displayed against the in-plane particle–particle distance (
$x\equiv |x_2-x_1|$
). One particle is located at
${\textbf q}_1=(0,0,3a)$
and the other at
${\textbf q}_2=(x,0,6a)$
, where
$a$
is the hydrodynamic radius of the Gaussian kernel and
$x$
varies from
$0$
to
$15a$
. Note that in this example, the components
$yx$
and
$yz$
vanish because both particles are located at
$y=0$
.
6.4. Self-mobility
Before presenting results for the self-mobility of blobs in the semi-bounded oscillatory Stokes set-up, we need to introduce several useful quantities: the mobility in the unbounded domain, derived by Mazur et al. (Reference Mazur, Bedeaux and Mazur1974
b), and the wall reaction field tensor
$\boldsymbol {{\mathcal{R}}}$
, derived by Felderhof (Reference Felderhof2005).
6.4.1. Self-mobility in the free space
The self-mobility of a particle oscillating in free space is
$\boldsymbol {{\mathcal{M}}}_{\infty }= \mathcal{M}_{\infty }\,\boldsymbol {{\mathcal{I}}}$
, where
$\boldsymbol {{\mathcal{I}}}$
is the identity matrix and
\begin{align}{\mathcal{M}}_{\infty } &= \int\int S({\textbf r}-{\textbf q})\, {\boldsymbol G}^{(\infty)} ({\textbf r}-{\textbf r}^{\prime }) \, S({\textbf r}^{\prime } -{\textbf q}) \, {\mathrm d}^3{\textbf r} \, {\mathrm d}^3{\textbf r}^{\prime } \nonumber\\ &= \frac {1}{(2\pi)^3} \int |\hat {S}({\textbf k})|^2 \frac {1}{\eta (\alpha ^2 +k^2)}\left ({\textbf I}-\hat {{\textbf k}} \hat {{\textbf k}}\right) \, {\mathrm d}^3{\textbf k},\end{align}
where
${\textbf G}_{\infty }({\textbf r}-{\textbf r}^{\prime })$
is the Green function for the oscillatory Stokes flow in free space (Appendix H). Mazur et al. (Reference Mazur, Bedeaux and Mazur1974
b) derived the self-mobility of a solid spherical particle or radius
$a$
, corresponding to a Heaviside kernel,
$S(r) = \mathbb{V}^{-1} \, \Theta (r-a)$
with
$\mathbb{V}=4\pi a^3/3$
, and obtained
\begin{equation} \xi _{\text{MB}}=\mathcal{M}_{\infty }^{-1}= 6\pi \eta a [1+\alpha a + (1/3) (\alpha a)^2]. \end{equation}
This self-friction
$\xi _{\text{MB}}$
corresponds to a net induced force transferred to the fluid
${\textbf F}=\xi _{\text{MB}} \times ({\textbf u}-{\textbf u}^{(0)})$
, where
${\textbf u}^{(0)} = \int S({\textbf r}-{\textbf q}) {\textbf v}_0({\textbf r}) {\mathrm d}^3{\textbf r}$
is the volume averaged unperturbed velocity (i.e. the base flow velocity in the absence of the particle disturbance). Starting from (6.12), in Appendix G, we derive the self-mobility of a Gaussian blob (G10). Figure 5 compares the continuum-based analytical relation (G10) with the discrete scheme results, resulting in a perfect match. We also include in figure 5 mobility of a hard sphere of radius
$a$
with no-slip boundaries, given by (6.13) (Mazur et al. Reference Mazur, Bedeaux and Mazur1974
b). Notably, the self-mobility of the Gaussian blob is remarkably close to the hard-sphere results, with just a small deviation in its imaginary part for
$a/\delta \sim 1$
. The same comment applies for the 3-pt Peskin kernel, whose mobility is also shown.
Figure 5. Comparison between the self-mobility in free space given by the Mazur–Bedeaux relation for hard spheres (inverse of (6.12) (line)) and for the Gaussian blob (blue line) derived in (G10). The latter coincides with the numerical evaluations for the discrete set-up using the present code and Gaussian blobs (black dots) and 3-points Peskin kernel (green stars) for the free-space set-up using
$L_{\parallel } = 192a$
and
$H = 3a$
(varying
$a/\delta$
). The support and grid cell size for the Gaussian kernel were computed as explained in the Appendix, while for the Peskin kernel, we used
$a = h_g$
.
6.4.2. Self-mobility in the semi-bounded domain
Felderhof derived the self-mobility of an isolated point particle (
$S(r)=\delta _{Dirac}(r)$
) at a distance
$\mathrm{h}$
from the wall when submitted to an oscillatory force in a semibounded domain. It is illustrative to clarify the meaning of this mobility, which includes the reaction flow reflecting from the wall back to the particle. The semi-bounded Green tensor can be decomposed as
${\textbf G}({\textbf r},{\textbf r}^{\prime })=\boldsymbol {{\mathcal{R}}}({\textbf r},{\textbf r}^{\prime })+ {\textbf G}^{(\infty)}({\textbf r}-{\textbf r}^{\prime })$
, where
$\boldsymbol {{\mathcal{R}}}({\textbf r},{\textbf r}^{\prime }) \equiv {\textbf G}({\textbf r},{\textbf r}^{\prime }) -{\textbf G}^{(\infty)}({\textbf r}-{\textbf r}^{\prime })$
was called by Felderhof the ‘reaction field’ to indicate that it is the flow arising from the wall disturbance in response to the flow induced by a force applied at the blob location. In a dispersion, with some base flow
${\textbf v}^{(0)}$
and blobs imposing forces
$\text{F}_j$
to the fluid, the velocity of the blob
$i$
located at
${\textbf q}= {\textbf s} + \mathrm{h} \,\hat {{\textbf z}}$
would be
\begin{equation} {\textbf u}_i={\textbf u}_i^{(0)} +\sum _{i\ne j} \boldsymbol {{\mathcal{M}}}_{ij} {\textbf F}_j + \boldsymbol {{\mathcal{M}}}_{ii} {\textbf F}_i = {\textbf u}_i^{(0)} +\sum _{i\ne j} \boldsymbol {{\mathcal{M}}}_{ij} {\textbf F}_j + \boldsymbol {{\mathcal{R}}}_{S} {\textbf F}_i + \mathcal{M}_{\infty } {\textbf F}_i. \end{equation}
We have introduced
$\mathcal{M}_{\infty }$
in (6.12) and the smeared self-reaction field,
\begin{align} \boldsymbol {{\mathcal{R}}}_S(\mathrm{h}) &=\int\! \int S({\textbf r}-{\textbf q}) \boldsymbol {{\mathcal{R}}}({\textbf r},{\textbf r}^{\prime }) S({\textbf r}^{\prime }-{\textbf q}) {\mathrm d}^3{\textbf r} {\mathrm d}^3{\textbf r}^{\prime }= \nonumber \\ & =\int\! \int S_z(z-\mathrm{h})\, \left [\int \boldsymbol {{\mathcal{\hat {R}}}}({\textbf k},z,z^{\prime }) |\hat {S}_{\parallel }({\textbf k})|^2 {\mathrm d}^2 {\textbf k} \right] \, S_z(z^{\prime }-\mathrm{h}) {\mathrm d}z {\mathrm d}z^{\prime } \end{align}
Note that due to the system (and kernel) symmetries,
$\boldsymbol {{\mathcal{R}}}_S(\mathrm{h})=R_{\parallel }(\mathrm{h}) \, \hat {{\textbf s}}\, +R_{\perp }(\mathrm{h}) \, \hat {{\textbf z}}$
with
$R_{xx}=R_{yy} =R_{\parallel }$
and
$R_{\perp }=R_{zz}$
. Analytical expressions for the components
$\boldsymbol {\mathcal{R}}$
in the point-particle limit (Felderhof Reference Felderhof2005) (
$S({\textbf r})=\delta _{Dirac}({\textbf r})$
in (6.15)) are provided in (H6) and (H7) of Appendix H for completeness. In the case of a Gaussian blob, (6.15) becomes analytically intractable and, in this case, we evaluate
$\boldsymbol {\mathcal{R}}$
from (6.15) numerically (see Appendix H).
For an isolated particle (
$\boldsymbol {{\mathcal{M}}}_{ij}=0$
for
$i\ne j$
), one can simplify (6.14) (dropping the
$i$
index) to get
\begin{align} \left ({\textbf u}-{\textbf u}^{(0)}\right)= {\mathcal{\boldsymbol M}}_{\text{MBF}}(\mathrm{h})\,{\textbf F}, \end{align}
where we have defined
$\xi _{MB}=\mathcal{M}_{\infty }^{-1}$
to introduce what we called the ‘Mazur–Bedeaux–Felderhof’ (MBF) self-mobility,
\begin{equation} {\mathcal{\boldsymbol M}}_{\text{MBF}}= \xi ^{-1}_{\text{MB}}\,\left (\xi _{\text{MB}} \boldsymbol {{\mathcal{R}}}_S(\mathrm{h})+\boldsymbol {{\mathcal{I}}}\right). \end{equation}
Figure 6. Variation of the self-mobility obtained using a Gaussian kernel with varying hydrodynamic radius at different heights. Numerical results (dots) are compared with the analytical equation obtained using the predictions for the mobility for a Gaussian kernel near a wall (continous lines) and also with the commonly used MBF friction (6.17) (dashed lines). (a) Real part
$zz$
; (b) imaginary part
$zz$
; (c) real part
$xx$
; (d) imaginary part
$xx$
. The numerical calculations have been done using a box width
$L_x=L_y=243a$
so the finite-size effect is completely negligible and a box height
$2H = 2h + 4a$
is large enough to have the full kernel inside the box.
Again,
${\mathcal{\boldsymbol M}}_{\text{MBF}}$
is a diagonal tensor, with parallel
$\mathcal{M}_{\text{MBF}}^{\parallel }(\mathrm{h})$
and perpendicular
$\mathcal{M}_{\text{MBF}}^{\perp }(\mathrm{h})=\mathcal{M}_{\text{MBF}}^{(z)}(\mathrm{h})$
components. Figure 6 shows the excellent agreement between the MBF self-mobility obtained from the continuum setting ((6.15) and (6.17) and solid lines in figure 6) and from the discrete setting (dots, calculated with the present numerical scheme). The dashed line in figure 6 is the result of (6.17) combining the hard-sphere MB self-friction
$\xi _{MB}$
with the Felderhof point-particle reaction field tensor (H6 and H7). We call ‘standard-MBF’ this ‘hybrid’ combination of kernels (i.e. hard sphere for
$\xi _{MB}$
via (6.13) and point-particle reaction field
$\mathcal{R}$
via (H6), (H7)), which proves to be valid for
$a/\delta \lt 0.06$
. Interestingly, standard-MBF worsens as the height
$h$
is decreased. To understand the origin of this discrepancy, we compare in figure 7 the point-particle and Gaussian-blob result for the reaction-field tensor components
$R_{\parallel }(h)$
and
$R_{\perp }(h)$
. As it was also argued by Simha et al. (Reference Simha, Mo and Morrison2018), as the height is reduced, we find that the point-particle reaction field (solid line) deviates from the finite-radius result (dots). From figure 7, the point-particle limit ceases to be appropriate approximately for
$h \lt 2\,a$
, i.e. when the size of the particle becomes comparable to its distance to the wall, as one should expect.
Figure 8. Dependence of the self-mobility on the free-space set-up against the inverse of the scaled periodic box side
$a/L_{\parallel }$
and different values of the
$\delta /a = [2\nu /(a^2 \omega)]^{1/2}$
ratio. Solid lines indicate the mobility for
$L\rightarrow \infty$
using (G10) and the dotted lines correspond to (6.19). Black dashed line in
$\mathrm{Re}[\mathcal{M}_{zz}]$
corresponds to the Hasimoto result for 3-D lattices (6.18) and coloured dashed lines to (6.19).
6.4.3. Finite size effects in the self-mobility
We now study the effect of periodic images in the self-mobility, which is illustrated in figures 8 and 9. In the case of steady Stokes flow in a 3-D periodic box of side
$L$
, Hasimoto (Reference Hasimoto1959) analytically derived the drag force acting on any one of the small fixed spheres (obstacles) forming the 3-D periodic array. Due to Galilean invariance, this result can be applied to the mobility of an array of particles in a system with zero total momentum (more precisely, the sum of fluid plus particles’ momentum). Thus, it applies to systems with a vanishing net external force, e.g. where forces (particle–particle and particle–fluid) satisfy Newton’s third law. In such a case, the self-mobility of a particle in the 3-D array decreases with the lattice cell side
$L$
as
\begin{equation} \mathcal{M}^{3\text{-}D}(L)= \mathcal{M}(\infty) \left ( 1 - 2.8373 a/L+\ldots \right),\end{equation}
where
$\mathcal{M}(\infty)=6\pi \eta a$
is the Stokes mobility of an isolated particle and
$\mathcal{M}^{3\text{-}D}(L)$
the self-mobility in the cubic 3-D lattice. The system we are studying here differs from the Hasimoto situation in two ways. First, in the doubly periodic case, the images of one particle are placed in a 2-D lattice extending over the
$xy$
plane. Second, we impose zero net flow in the
$z$
direction (
${\textbf v}_{{wall}}\cdot \hat {{\textbf z}}=0$
in 4.7), while the total momentum in the plane does not vanish (see 4.17). Physically, this corresponds to a (infinite or semi-infinite) box which is fixed in the
$z$
direction, but it is allowed to shake in the plane. In terms of self-mobilities, one should expect that the
$zz$
component should decrease with the cell side
$L_{\parallel }$
, similarly to the Hasimoto picture. However, in the plane, the self-mobility is more resemblant to the settle velocity under jittering gravity and, as the box is enlarged, it should eventually increase with the number of particles per area in the periodic 2-D lattice (thus, as
$L_{\parallel }^{-2}$
). It is actually crucial for the present scheme to recover these strong in-plane hydrodynamic couplings which are present in applications, such as QCM (Delgado-Buscalioni Reference Delgado-Buscalioni2024). In the free-space doubly periodic set-up, the self-mobility corresponds to a very large box
$L_z \rightarrow \infty$
whose position is fixed in the vertical direction, but not in the plane. The effect of reducing
$L_{\parallel }$
in the free-space self-mobility is illustrated in figure 8. As predicted,
$\mathcal{M}_{zz}(L_{\parallel })$
is consistent with the results reported by Scalfi et al. (Reference Scalfi, Vitali, Kiefer and Netz2023) for a triply periodic, momentum-conserving molecular system. In their study, they examine the effect of periodic boundary conditions by summing (using the Ewald summation technique) the contribution of the interaction with each image to the velocity of the central particle. As shown in figure 8, at very low frequencies (
$\delta /a =48$
or
$\omega a^2/\nu \sim 10^{-3}$
),
$\mathrm{Re}[\mathcal{M}_{zz}(L_{\parallel })]$
converges to the Hasimoto result (6.18) within all the studied window
$5\lt L_{\parallel }/a\lt 200$
. As the frequency is increased, convergence takes place for more compact lattices (e.g. for
$L_{\parallel } \lt 12 a$
if
$\delta /a =5$
or
$\omega a^2/\nu \sim 10^{-1}$
). As
$L_{\parallel }$
increases, the imaginary
$zz$
-component of the mobility slowly converges to the mobility of an isolated particle (evaluated from (G10) for the Gaussian blob, see figure 6). Notably, the convergence is significantly slower at lower frequencies (large
$\delta /a$
values). Similarly, for the
$xx$
component, convergence to the
$L_{\parallel }\rightarrow \infty$
limit slows down with the frequency. However, opposite to the
$zz$
component,
$\mathcal{M}_{xx}(L_{\parallel })$
increases as one reduces the separation between images: as predicted, we find
$\mathcal{M}_{xx}\sim L_{\parallel }^{-2}$
due to the increase in the net force applied per unit area and the resulting hydrodynamic couplings between images.
Figure 9. Semibounded set-up: variation of the
$xx$
and
$zz$
terms of the mobility with the box size (points) compared with the predictions for an infinite box (solid lines) calculated using (6.17). For the real components, a comparison is also made with Hashimoto’s result for particles in a 3-D lattice, (6.18). The penetration length employed to make the calculations was
$\delta = 48 a$
. Dashed lines for
$\mathcal{M}_{xx}$
correspond to (6.19).
Finite-size effects in the semi-bounded case are particularly relevant to understand QCM signals, whose dissipation shift shows significant coverage effects due to hydrodynamic couplings between suspended analytes even at extremely low concentrations (Delgado-Buscalioni Reference Delgado-Buscalioni2024). Figure 9 presents
$\mathcal{M}_{xx}$
and
$\mathcal{M}_{zz}$
in the semi-bounded set-up against
$a/L_{\parallel }$
(we recall
$L_{\parallel } = L_x = L_y$
) and various
$h/\delta$
values. We present results for
$\delta /a =48$
which correspond to small particles within the point-particle limit, which (figure 7) is valid at least for
$h\gt 2a =\delta /24$
. In all cases, we consistently observe that the mobility converges to the
$L_{\parallel }\rightarrow \infty$
limit given by the MBF result in (6.17). As in the slit set-up, the
$zz$
component converges to the Hasimoto result (6.18) as the cell side
$L_{\parallel }$
is decreased. The imaginary part
$\mathcal{M}{zz}$
follows a different decay pattern (so far, not analytically predicted). Finite-size effects are noticeable up to significantly small particle coverage
$(a/L_{\parallel })^2\sim 10^{-4}$
, being more relevant as the particle height
$h$
is increased (which is consistent with QCM observations (Delgado-Buscalioni Reference Delgado-Buscalioni2024)). In the plane, the mobility increase with
$L_{\parallel }$
has the same origin as that explained for the slit set-up, being related to the hydrodynamic coupling between images in those directions where the net fluid momentum is not conserved. We have fitted the
$xx$
mobility as
\begin{equation} \mathcal{M}_{xx}(L,h) = \mathcal{M}_{xx}(\infty,h)\cdot \left [1 + C_1(h/\delta)\frac {a}{L_{\parallel }} + C_2(h/\delta)\left (\frac {a}{L_\parallel }\right)^2\right], \end{equation}
where
$C_1 \approx -4.94 s^{0.27} - 0.724 s$
and
$C_2 \approx 454 s^{0.71} + 352 s^{1.57}$
, where
$s \equiv h/\delta \in [0.1, 0.4]$
. Again, finite-size effects are more noticeable in the imaginary part of the mobility, reaching up to very small particle coverages
$(a/L_{\parallel })^2\sim 10^{-4}$
. As an aside, (6.19) also captures the large
$L_{\parallel }$
trend of
$\mathcal{M}_{xx}$
in the open slit set-up (figure 8).
7. Discussion
This work presents a novel spectral solver for oscillatory Stokes hydrodynamics in doubly periodic (open) boundaries. Similarly to the steady Stokes version (Hashemi et al. Reference Hashemi, Peláez, Natesh, Sprinkle, Maxian, Gan and Donev2023), the present method uses FFTs in the two periodic directions and Chebyshev transforms in the aperiodic direction. The method has been implemented in efficient implementation on GPUs by using a single 3-D FFT for each Fourier–Chebyshev transformation (Pérez Peláez Reference Peláez-Pérez2022; Hashemi et al. Reference Hashemi, Peláez, Natesh, Sprinkle, Maxian, Gan and Donev2023). Notably, the computational mesh is restricted to a finite relevant domain
$z\in [-H,H]$
, which is connected to the analytical solution of the exterior flow based on plane wave expansion. The scheme first derives the solution for the free space (i.e. without external boundaries at
$z=\pm H$
) and then adds a correction field to enforce the boundary conditions imposed by the boundaries (walls). These boundary conditions for the aperiodic direction
$z$
connect the exterior and interior flows, being derived from the Dirichlet-to-Neumann map. Here, we implement a no-slip bottom boundary (at
$z=-H$
) which oscillates in the plane with velocity
$v_{wall}$
. However the solver can be generalised in different ways, as we comment below. A relevant feature of the present scheme is that it is written in the frequency domain, i.e. solving the time-Fourier-transformed equations. Thus, fields are subjected to four different spectral transformations: ‘time’ evolution in the frequency domain (resulting in complex phasor fields), Fourier-based spectral description for the
$x, y$
directions (resulting, due to linearity, in independent equations for each in-plane wave vector) and finally, a spectral Chebyshev solver in the
$z$
direction for each wavevector. While this certainly complicates the solver structure, it translates into a computationally efficient algorithm. First, time-stepping is avoided, which brings a strong benefit in computational time, as the equations directly solve the (long-time) oscillatory regime, avoiding time-stepping and transient propagation of vorticity over the system. The gain in computational time can be estimated by considering the Courant number limitation
$\text{CFD}={\mathrm d}t \,\nu /h^2\sim 0.1$
on the time step
${\mathrm d}t$
arising from viscous momentum propagation across a given mesh size
$h$
. A time-stepping scheme will require a time
$T_{\nu }\sim L^2/\nu$
to propagate momentum across the box side
$L$
and a time
$T_{\omega }=2\pi /\omega$
to cover one oscillation. Consider a typical QCM scenario with
$\omega = 2\pi \times 10^6 \,\text{Hz}$
and penetration length
$\delta \sim 100\, \text{nm}$
. Taken as example the calculations carried out by Vázquez-Quesada et al. (Reference Vázquez-Quesada, Schofield, Tsortos, Mateos-Gil, Milioni, Gizeli and Delgado-Buscalioni2020), one would set a mesh of
$h\sim 4 \,\text{nm}$
to resolve a liposome of
$\sim 50\, \text{nm}$
radius placed in a tall enough box (to avoid finite size effects, i.e.
$L\sim 4 \delta \sim 400 \,\text{nm}$
). The number of time steps needed to cover the transient momentum propagation is
$T_{\nu }/{\mathrm d}t \sim 10^5$
, and this number also roughly correspond to
$T_{\omega }/{\mathrm d}t$
. By constrast, the present code will only require a single step to obtain the mobility of the immersed body. The time reduction is therefore very large. As an example, our previous simulations of
${\sim} 50 \, \mathrm{nm}$
spherical particles under wall oscillations using the standard time-stepping integration (Vázquez-Quesada et al. Reference Vázquez-Quesada, Schofield, Tsortos, Mateos-Gil, Milioni, Gizeli and Delgado-Buscalioni2020) required between 1 and 10 hours (depending on size and frequency), while calculations based on the present code (to be submitted), which required some iterations to take into account the viscoelastic stress of the structure in the fluid phase, required less than one minute. Even more importantly, simulations of proteins
${\sim} 3\, \text{nm}$
under unsteady Stokes flows with
$\delta \sim 100\, \text{nm}$
are just impossible in a standard code, due to the length scale disparities. Here, we can set
$H\sim 10\,\text{nm}$
and focus on the relevant domain surrounding the protein. Second, it is also important to state that the time for a single calculation of the present scheme scales linearly with the number of particles, and roughly linearly
$O(N\log N)$
with the number of fluid cells, thanks to the use of fast Fourier and Chebyshev transforms. Through careful bookkeeping and manipulation, our implementation only requires the application of a single 3-D FFT call to carry out the Fourier transform in the
$x$
and
$y$
directions plus the Chebyshev transform in the
$z$
direction. Details on the implementation of the fast Fourier and Chebyshev transform in similar boundary value problems can be found from Maxian et al. (Reference Maxian, Peláez, Greengard and Donev2021) and Pérez Peláez (Reference Peláez-Pérez2022, Appendix C) . The nature of the oscillatory complex equations impose several technical problems, complicating the Dirichlet–Neumann map, and also in the solution for the zero mode
$k=0$
, which required performing an integral involving the force density over the whole domain. This steps were efficiently solved by combining monomial expansion and Chebyshev integration (Appendix B).
We performed an extensive series of tests, with the scheme acting either as a free-space solver or for the semi-bounded domain (including a rigid no-slip wall at the bottom
$z=-H$
). The Green tensor for the semi-bounded domain in the reciprocal space was measured by exciting the flow with local force distributed over a smeared delta function (Gaussian kernel). The result is in excellent agreement with the smeared version of the analytical formulae derived by Felderhof (Reference Felderhof2005) (collected and corrected for the typo published in the scattered literature (Felderhof Reference Felderhof2006)). The spectral nature of the scheme translates in an excellent accuracy, presenting relative errors smaller than
$5\,\%$
with relatively large grid sizes
$h_{{grid}} \lt \delta /4$
, compared with
$h_{{grid}} \lt \delta /15$
for standard collocation discretisation (Vázquez-Quesada et al. Reference Vázquez-Quesada, Schofield, Tsortos, Mateos-Gil, Milioni, Gizeli and Delgado-Buscalioni2020; Delgado-Buscalioni Reference Delgado-Buscalioni2024) (we recall that
$\delta$
is the fluid penetration length). Self-mobilities and mutual mobilities were analysed using blobs defined by Gaussian kernels (we also test 3-point Peskin kernels) with a well-defined hydrodynamic radius
$a$
(Pérez Peláez Reference Peláez-Pérez2022). Importantly, using two ‘blobs’ (a force-emitter and a velocity-receiver), we verified that the reciprocal relations for the mobility tensor are satisfied to machine precision. Concerning the self-mobility, we found that the Gaussian blob mobility provides excellent agreement with the Mazur–Bedeaux mobility for a hard sphere oscillating in free space. In the semi-bounded domain, Felderhof’s reaction field for point-particles was found to be consistent with our results, provided the particle height over the surface is larger than its diameter. As previously expected (Simha et al. Reference Simha, Mo and Morrison2018), near to the surface, the point-particle limit breaks down. A more precise evaluation of the self-mobility of a finite particle near the wall would then require adding higher moments in the particle response (rotlet, stresslet couplings to translation) which is left for a future contribution. Finite-size effects due to hydrodynamic coupling between image particles in the doubly periodic system have been also investigated. In this set-up, the periodic images spread over the 2-D plane which introduce quantitative differences with respect to a 3-D periodic space. Finite-size effects propagate as
$L_{\parallel }^{-2}$
leading to strong hydrodynamic couplings particularly for the imaginary part of the in-plane mobility, which are present up to very small coverage
$\Theta = \pi a^2/L_{\parallel }^2\sim 10^{-4}$
. This is related to the reported sub-linear power-law scaling of the dissipation signal
$\Delta D$
in QCM of suspended particles (Delgado-Buscalioni Reference Delgado-Buscalioni2024), presenting
$\Delta D \sim \Theta ^{0.8}$
up to extremely small coverages.
The present fluid solver is the core of a new class of schemes to measure viscoelastic response of soft structures under oscillatory flow with applications, inter alia, to QCM of biomolecules, vibrational spectra of AFM or oscillatory excitation of vesicles, which will be presented in future contributions. The boundary condition at the bottom surface is here specified in the reciprocal space
\begin{align} \boldsymbol {\mathcal{B}}_{\pm H}[{\textbf v}] = {\textbf B}_{wall}({\textbf k},\pm H), \end{align}
so it can be generalised in several ways. Here,
$\boldsymbol {\mathcal{B}}_{\pm H}$
is a linear operator and
${\textbf B}_{wall}({\textbf k})$
is a phasor, in general, wavevector dependent. Here, we consider homogeneous surfaces (
${\textbf B}_{wall}=0$
for
$k\ne 0$
and
${\textbf B}_{wall}={\textbf v}_{wall}$
for
$k=0$
) with no-slip condition (
$\mathcal{B}_{-H}=1$
), but one can also impose partial slip using
$\mathcal{B}[v] = \partial _z v - \beta v$
, where
$\beta$
is the inverse of the slip-length. Additionally, the formalism allows further generalisations as the operator
$\mathcal{B}({\textbf k}, z=-H)$
may introduce partial derivatives with respect to
$z$
and in the plane, via the Fourier map
$\nabla _{{\textbf s}} \rightarrow \mathrm{i} {\textbf k}$
. This allows for general coarse-grained implementations of patterned surfaces, even active surfaces with a non-trivial wavevector dependent
${\textbf B}_{wall}({\textbf k})$
, or coupling with viscoelastic layers, via continuity of tangential stress and velocity. Slip velocity field due to ion transport in thin Debye layers (Schnitzer & Yariv Reference Schnitzer and Yariv2012) can be also easily implemented. In contrast, being monochromatic, the scheme is restricted to the linear response. Yet, it could be used as a fast tool to measure the drift forces arising from the time average of nonlinear terms such as those driving ultrasound forces (Bruus Reference Bruus2012; Balboa Usabiaga & Delgado-Buscalioni Reference Balboa Usabiaga and Delgado-Buscalioni2013) or in electrokinetics (Ramos et al. Reference Ramos, Morgan, Green and Castellanos1998). To this end, the scheme could be coupled to other transport equations (such as fluid compressibility and Poisson–Planck–Nerst, which has been already solved in real-time using open doubly periodic domains (Maxian et al. Reference Maxian, Peláez, Greengard and Donev2021)). Solvers for doubly periodic domains including transient effects in the time domain are also feasible and could be practical to regimes with strong aperiodic unsteadiness.
Acknowledgements
Funding
We acknowledge funding from the Spanish Agencia Estatal de Investigacion (AEI) with grant numbers PID 2020–117080RB-C51 and PDC 2021–121441-C21 (Prueba de Concepto program). Also AEI ‘María de Maeztu’ Program for Units of Excellence (grant R&D CEX 2023–001316-M) and funding from HORIZON-EIC-2023 PATHFINDER program (grant ID 101130615 ‘FASTCOMET’).
Declaration of interests
The authors report no conflict of interest.
Author contributions
R.P.P implemented the code, P.P-A derived analytical theory, performed the tests and error estimations, and R.D-B derived the numerical scheme and wrote the manuscript. All authors contributed equally to analysing data, reaching conclusions and revising the paper.
Data availability
The CUDA implementation of the present solver makes extensive use of the UAMMD library (Peláez et al. Reference Peláez, Ibáñez-Freire, Palacios-Alonso, Donev and Delgado-Buscalioni2025) and will be available at https://uammd.readthedocs.io/en/latest/.
Appendix A. Solution of (3.22), (3.29)
Equations (3.23) and (3.24) are ordinary differential equations for
$w(z)\equiv {\textbf v}^{\parallel }({\textbf k},z)$
, with generic form,
\begin{eqnarray} w^{\prime \prime }_1 -\tilde {k}^2 w_1 = \tilde {C}_1 e^{-kz}, \end{eqnarray}
\begin{eqnarray} w^{\prime \prime }_2 -\tilde {k}^2 w_2 = \tilde {C}_2 e^{kz},\end{eqnarray}
where we have defined
\begin{align} \tilde {C}_i =\mathrm{i} \frac {{\textbf k} C_i}{\eta }. \end{align}
Let us deal with
$w_1$
. The homogeneous solution is
$w_{1,h} = A_{1} e^{-\tilde {k} z} +A_{2} e^{\tilde {k} z}$
. A particular solution is
$w_{1,p} = A_{1,p} e^{-k z}$
because
\begin{align} w_{1,p}^{\prime \prime } - \tilde {k}^2 w_{1,p} = A_{1,p} \left(k^2-\tilde {k}^2 \right) e^{-kz} = \tilde {C}_1 e^{-kz}. \end{align}
Therefore,
\begin{equation} A_{p,1} = \frac {\tilde {C}_1}{k^2-\tilde {k}^2} = -\mathrm{i} \frac {{\textbf k} C_1}{\eta \alpha ^2}. \end{equation}
A similar result applies for
$w_2$
, for which the particular solution is
$w_{2,p}=A_{p,2} e^{kz}$
and
\begin{equation} A_{p,2} = -\mathrm{i} \frac {{\textbf k} C_2}{\eta \alpha ^2}. \end{equation}
In the case of (3.29), the homogeneous solution is
$v^{(h)}=A_{\pm } e^{\pm \tilde {k} z}$
and the particular solution of the non-homogeneous equation is now
$w(z) = A^{p}_{\pm } e^{\pm k z}$
, leading to
$(k^2 -\tilde {k}^2) A^{p}_{\pm } e^{\pm k z} = -\alpha ^2 A^{p}_{\pm } e^{\pm k z}$
and therefore
\begin{eqnarray} -\alpha ^2 A^{p}_{+} = \frac {k C_2}{\eta }, \end{eqnarray}
\begin{eqnarray} -\alpha ^2 A^{p}_{-} = -\frac {k C_1}{\eta }. \end{eqnarray}
The solution in (3.32) and (3.33) is indicated with
$A_{\perp,1}=A_{-}$
and
$A_{\perp,2}=A_{+}$
.
Appendix B. Green function formalism for the
$k=0$
parallel velocity
We particularise the problem for the domain
$z\in [-H,\infty]$
. The Green function of the problem satisfies
\begin{eqnarray} \partial _x^2 G(x,y)-\alpha ^2 G(x,y) = \delta (x-y), \end{eqnarray}
\begin{eqnarray} G(-H,y) =0, \end{eqnarray}
\begin{eqnarray} G(\infty,y) =0, \end{eqnarray}
\begin{eqnarray} \partial _x G(y_{+},y)-\partial _x G(y_{-},y) = 1, \end{eqnarray}
\begin{eqnarray} G(y_+,y) = G(y_{-},y) .\end{eqnarray}
The relation (B4) is the jump condition, obtained by integrating (B1) from
$x=y-\epsilon$
to
$x=y+\epsilon$
and then taking the limit
$\epsilon \rightarrow 0$
. The relation (B5) imposes continuity of
$G(x,y)$
at
$x=y$
.
For
$x\ne y$
, (4.8) is homogeneous and
\begin{align} G(x,y) = a_i e^{-\alpha x} + b_i e^{\alpha x} \quad \mathrm{for}\,x\ne y,\end{align}
where
$a_1$
and
$b_1$
are used for the domain
$x\lt y$
and
$a_2$
,
$b_2$
for the domain
$x\gt y$
. We now apply the boundary conditions:
-
(i) for
$x\lt y$
,
$G(-H,y)=0$
, so
$a_1 e^{\alpha H} + b_1 e^{-\alpha H}=0$
; -
(ii) for
$x\gt y$
,
$G(\infty,y)=0$
, so
$b_2=0$
.
This yields
\begin{align} b_1 = -a_1 e^{2\alpha H} \end{align}
and
\begin{eqnarray} a_1 = \frac {1}{2\alpha } \left (e^{-\alpha (y+2 H)}\right), \end{eqnarray}
\begin{eqnarray} a_2 = \frac {1}{2\alpha } \left (e^{-\alpha (y+2H)} -e^{\alpha y}\right). \end{eqnarray}
Therefore,
\begin{equation} G(x,y)= \left \{\begin{array}{cc} \dfrac {1}{2\alpha } \left (e^{-\alpha (x+y+2H)} -e^{\alpha (x-y)}\right) & x\lt y, \\ &\\ \dfrac {1}{2\alpha } \left (e^{-\alpha (x+y+2H)} -e^{-\alpha (x-y)}\right) & x\gt y. \end{array} \right . \end{equation}
As is customary, the solution is built as
\begin{align} v(x)=\int G(x,y) [-f(y)] {\mathrm d}y,\end{align}
and reassigning variables
$z\rightarrow x$
and
$z'\rightarrow y$
, we get
\begin{equation} {\textbf v}^{\parallel }(z)= -\int _{-H}^z [G(z,z')]_{z\gt z'} {\textbf f}^{\parallel }(z') {\mathrm d}z' - \int _z^{\infty } [G(z,z')]_{z\lt z'} {\textbf f}^{\parallel }(z') {\mathrm d}z'. \end{equation}
Or, specifically,
\begin{align} {\textbf v}^{\parallel }(z) & = - \frac {1}{2\eta \alpha } \left [\int _{-H}^z \left (e^{-\alpha (z+z'+2H)} -e^{-\alpha (z-z')}\right) {\textbf f}^{\parallel }(z') {\mathrm d}z' \right . \nonumber\\[4pt] & \quad \left . +\, \int _z^{\infty } \left (e^{-\alpha (z+z'+2H)} -e^{\alpha (z-z')}\right) {\textbf f}^{\parallel }(z') {\mathrm d}z' \right]. \end{align}
As
${\textbf f}(z)=0$
for
$z\gt H$
, this can be written as indicated in the main text (4.11).
Appendix C. Integrals via Chebyshev expansion.
The
$k=0$
mode requires the evaluation of the integral
\begin{align} \int _{-H}^H \sinh [\alpha (z'+H)] {\textbf f}^{\parallel }(z') {\mathrm d}z', \end{align}
where
${\textbf f}^{\parallel }(z)=\int\! \int \boldsymbol {I}_{\parallel } {\textbf f} {\mathrm d}x\,{\mathrm d}y$
is the projection of the net force in the plane and
$\boldsymbol {I}_{\parallel } \equiv \hat {{\textbf x}} \hat {{\textbf x}} +\hat {{\textbf y}} \hat {{\textbf y}}$
. To simplify the notation, we consider the integral for a generic function
$f(x)$
for which we have already obtained an expansion in Chebyshev polynomials,
\begin{align} f(x)=\sum _m f_m T_m(x). \end{align}
Using the transformation of variable
$x\rightarrow z/H$
, the integral becomes
\begin{align} H \int _{-1}^1 h(x)f(x) {\mathrm d}x,\end{align}
where
$h(x)\equiv \sinh [\alpha (x+1)]$
and now
$\alpha \rightarrow \alpha H$
. The weight function for the Chebyshev polynomials
\begin{align} w(x)=\frac {1}{\left (1-x^2\right)^{1/2}} \end{align}
is involved in their ortogonality conditions, which, in continuous form,
\begin{align} c_n \delta _{n m}\equiv \int _{-1}^1 T_n(x) T_m(x) w(x) {\mathrm d}x = \left \{ \begin{array}{cr} 0 & n\neq m, \\[2pt] \pi & n=m=0, \\[2pt] \dfrac {\pi }{2} & n=m\neq 0. \\ \end{array} \right . \end{align}
For any function
$h(x)$
, we define
$h^w(x)=h(x)/w(x)$
, so the integral can be then rewritten as
\begin{align} \int _{-1}^{1} h(x) f(x) {\mathrm d}x = \int _{-1}^{1} \frac {h(x)}{w(x)} f(x) {\mathrm d}z = \int _{-1}^{1} h^w(x) f(z) w(x) {\mathrm d}z .\end{align}
Expanding in Chebyshev and using the orthogonality conditions,
\begin{align} \int _{-1}^{1} h(x) f(x) {\mathrm d}x = \sum _{n,m} h^w_n f_m \int _{-1}^{1} T_n(x) T_m(x) w(x) {\mathrm d}x = \sum _{n} c_n h^w_n f_n,\end{align}
with
$c_0= \pi$
and
$c_{n\gt 0}=\pi /2$
. The coefficients of the Chebyshev expansion for
$h^w(x)$
are
\begin{align} h^w_n = \int _{-1}^{1} h^w(x) w(x) T_n(x) {\mathrm d}x= \int _{-1}^{1} h(x) T_n(x) {\mathrm d}x .\end{align}
The integrals can be evaluated by expanding
$h(x)$
and
$T_n(x)$
in monomials, e.g.
$h(x)=\sum _l H_l x^l$
and
$T_n(x)=\sum _i T_{n,i} x^i$
. As we use an expansion in Chebyshev of order
$N$
, without loss of accuracy, we just need to evaluate
\begin{align} h^w_n = \sum _{l=1}^{N-n} H_l \int _{-1}^{1} x^l T_n(x) {\mathrm d}x. \end{align}
The monomial expansion for the Chebyshev polynomials is well known,
\begin{align} T_n(x) = \sum _{k=0}^{\lfloor \frac {n}{2}\rfloor } a_{k,n}x^{n-2k}, \end{align}
where
\begin{align} a_{k,n} = \frac {n}{2}(-1)^{k}\frac {2^{n-2k}(n-k-1)!}{k!(n-2k)!} \end{align}
and the prefactor accompanying
$x^l$
is thus,
\begin{align} T_{n,l} = \begin{cases} a_{\frac {n-l}{2},n} & (n-l)\; \mathrm{even},\\ 0 &\; \mathrm{otherwise}.\end{cases} \end{align}
It should be noted that the factorials become too large to be computed when
$l\gt 20$
, but fortunately, the expression for
$a_{k,n}$
can be simplified using Stirling approximation,
\begin{align} a_{k,n} \approx \frac {n}{2k!}(-1)^{k}\exp \left [(n-2k)\log (2) + s_f(n-k-1) - s_f(n-2k)\right], \end{align}
where
\begin{align} s_f(n) = n\log (n) - n - 1 + \frac {\log (2\pi n)}{2} + \log \left (1+ \frac {1}{12n} + \frac {1}{288n^{2}} - \frac {139}{51\,840n^{3}}\right). \end{align}
Note that the
$k!$
factor can also be included in the exponential if necessary. Using Stirling approximation for
$l\gt 18$
yields a relative error smaller than
$10^{-10}$
in the computation of
$a_{k,n}$
. We also use Stirling approximation to evaluate the coefficients of the monomial expansion of
$h(x)$
. Since there is no approximation for the lower coefficients (up to
$l\leqslant 15$
), the overall computation of the integral is unaffected precision-wise. This procedure is much faster than another alternative based on the recurrent expressions
$T_{n}=2T_{n-1}-T_{n-2}$
Appendix D. Kernels
Kernels are smooth functions that we use to spread the forces acting on the blobs to the fluid and, once the fluid behaviour has been resolved, to interpolate the fluid’s velocity to the particles. There is a wide variety of kernels, each with its own specific characteristics. In this work, we have mainly used a Gaussian kernel, as it allows for achieving an arbitrary level of accuracy easily, as we will see later. Additionally, we have compared the results obtained using this kernel with those that would be achieved using the 3-point kernel by Peskin (Reference Peskin2002) which is less accurate but faster. The kernels depend on the distance
$\mathbf r=(x,y,z)$
of the particles to any point in the box, and usually they can be separated by dimensions,
\begin{equation} \delta (x,y,z)=\phi (x)\phi (y)\phi (z) .\end{equation}
D.1 Gaussian kernel
The Gaussian kernel, as its name suggests, employs a Gaussian function to transfer the forces from the particles to the fluid,
\begin{equation} \phi _G(r_\alpha) = \frac {1}{\left (2 \pi \sigma ^2\right)^{1/2}} \exp \left [-\frac {r_\alpha ^2}{2 \sigma ^2}\right] \quad \mathrm{for}\, r\lt r_c, \end{equation}
where
$ r_\alpha = x$
,
$ y$
or
$ z$
, and
$\sigma$
represents the width of the Gaussian, which is related to the hydrodynamic radius of the blob as
$ a = \sqrt {\pi } \sigma$
. In principle, the Gaussian function is greater than zero for any distance; however, to increase the efficiency of the spreading and interpolation functions, the kernel is truncated at a maximum distance
$ r_c$
. This distance is determined by a tolerance parameter
$\epsilon$
, which indicates that if
$\delta (r) \leqslant \epsilon$
, it can be approximated to 0. The cutoff radius of the kernel is related to the tolerance according to
\begin{equation} r_c(\epsilon) = \left (-2^{1/2} \log \left [-(2\pi \sigma ^2)^{d/2} \epsilon \right]\right)^{1/2} \, \sigma,\end{equation}
where
$ d$
indicates the number of dimensions of the system. The other important variable that determines the accuracy of the kernel is the fluid mesh
$ h = L_{box}/N_{cells}$
. In principle, this value could be set arbitrarily, but to ensure consistency in the accuracy with the tolerance that will be used, we define it using the expression
\begin{equation} \sigma (\epsilon) = \pi ^{1/2} \min [0.55-0.11 \log _{10}(3\epsilon), 1.65] \, h. \end{equation}
D.2 The 3-point Peskin kernel
The 3-point Peskin kernel uses a more compact function than the Gaussian kernel,
\begin{equation} \phi _{3pt}(\tilde {r}_\alpha) = \begin{cases} \dfrac {1}{3}\left ( 1 + \sqrt {1-3\tilde {r}_\alpha ^2}\right) & |\tilde {r}_\alpha | \lt 0.5, \\[4pt] \dfrac {1}{6}\left (5 - 3|\tilde {r}_\alpha | - \sqrt {1-3(1-|\tilde {r}_\alpha |)^2}\right) & 0.5 \leqslant |\tilde {r}_\alpha | \lt 1.5,\\[4pt] 0 & |\tilde {r}_\alpha | \gt 1.5, \end{cases}. \end{equation}
where
$ \tilde {r}_\alpha = r_\alpha / h$
. In this case, the support of the kernel always spans three cells and the hydrodynamic radius is
$a=0.91\,h_{\text{grid}}$
(Pérez Peláez Reference Peláez-Pérez2022). This allows the spreading and interpolation functions to be faster, at the cost of losing some accuracy.
Figure 10. Relative error obtained when calculating the xx (left column) and zz (right column) components of the mobility using: (a,b) a 3-point Peskin kernel; (c,d) a Gaussian kernel with
$\epsilon = 10^{-1}$
, resulting in a support of three cells; (e,f) a Gaussian kernel with
$\epsilon = 10^{-3}$
, resulting in a support of seven cells; and (g,h) a Gaussian kernel with
$\epsilon = 10^{-5}$
, resulting in a support of 11 cells.
To evaluate the accuracy of the kernels and analyse how the precision of the Gaussian kernel varies with tolerance, we studied the relative error of the mobility obtained depending on the position of the blob within the cell. The relative error was calculated as
\begin{equation} RE = \frac {||\mathcal{M}_{\alpha,\alpha }(x,y) - \langle \mathcal{M}_{\alpha,\alpha }\rangle ||}{||\langle \mathcal{M}_{\alpha,\alpha }\rangle ||}, \end{equation}
where
$\langle \ldots \rangle$
denotes the average over the entire cell and
$||\ldots ||$
represents the norm of the complex number. The results are shown in figure 10. The Peskin kernel has a relative error that varies between 2 % and 4 % in most of the cell, which is considerably better than the error obtained when using a Gaussian kernel with the same support (three cells), which is approximately 8 %. However, when we make
$\epsilon$
smaller, we achieve much lower errors, approximately 0.1 % when
$\epsilon = 10^{-3}$
and close to 0.001% when
$\epsilon = 10^{-5}$
. Throughout this work, a Gaussian kernel with
$\epsilon = 10^{-5}$
has been used.
Appendix E. Discretisation errors
To study the discretisation errors, we analyse how the relative error varies between the Felderhof’s smeared Green tensor (6.2) and the numerical calculations at each height z of the box when we apply a force on a Gaussian blob of radius
$a=10 h_{{grid}}$
at a height
$h=2.5 a$
in a box on length
$60 h_{{grid}}$
and height
$70 h_{{grid}}$
. We have studied the error by varying
$\delta$
between
$\delta = 0.004 h_{{grid}}$
and
$\delta = 2 h_{{grid}}$
.
Figure 11. Relative error between the analytical result of the Gaussian kernel smeared Green’s function,
$\boldsymbol {{\mathcal{G}}}_S({\textbf k}, z, \mathrm{h})$
(6.2), and the numerical results for the wavenumber
${\textbf n}=(1,1)$
. The smeared Green’s tensor is computed by applying a force to a Gaussian blob of radius
$a = 10 h_g$
, located at a height
$h = 2.5a$
above the wall. The relative error is shown as a function of the fluid cell position
$z$
and the penetration length. The left panels illustrate the relative error for the real part of Green’s tensor, while the right panels show the imaginary part. The top panels correspond to the xx component, the middle panels to the xz component and the bottom panels to the zz component.
In figure 11, we show the relative errors for the real (left panels) and imaginary (right) terms of the xx (top), xz (central) and zz (bottom) components of
$\boldsymbol {{\mathcal{G}}}_S({\textbf k}, z, \mathrm{h})$
. We can see that for
$\delta \gt 5 h_{{grid}}$
, the relative error is always smaller than a
$1\,\%$
and using a value of
$\delta \approx h_{{grid}}$
still leads to relative errors smaller than a
$10\,\%$
. It should be mentioned that the largest relative error for
$\delta /h_{{grid}} \simeq 1$
corresponds to
$z/a \sim 5$
in the real part of the xx mobility, due to the vanishing of the Green function component.
Appendix F. Zero finite-effects in vertical direction for arbitrary
${{H}}$
One of the main advantages of the present scheme is that the height of the inner domain (or computational domain)
$2H$
is an arbitrary parameter that does not influence the result, provided that all forces are exerted within it. This allows the use of very small computational boxes, with
$H$
much smaller than the fluid’s penetration length
$\delta$
. Moreover, for a given limitation in computer memory, this feature permits to significantly extend the area covered by the
$xy$
-plane (
$L_xLy)$
, to reach very small particle coverages
$\Theta$
. To prove the
$H$
-independence of the numerical scheme, we show in figure 12 the
$xx$
and
$ zz$
components of the self-mobility of a particle with radius
$ a = 0.1\delta$
located at height
$ h = 0.5\delta$
. Results for various values of
$H$
, in figure 12, clearly prove our claim showing a constant mobility for any
$H\gt 0.3\delta$
. The limit value of the computational height is determined by
$H\leqslant (h + a) / 2 = 0.3 \delta$
and corresponds to a portion of the particle kernel being outside the inner-domain (
$0\lt z\lt 2H$
). In such a case, not all forces spread to the fluid, resulting in the incorrect mobility. In conclusion, the only condition for
$H$
is that it should be sufficiently tall for all particles to be fully contained within the inner-domain
$0\lt z\lt 2H$
.
Figure 12. Dependence of the (a) xx and (b) zz components of the mobility on
$ H$
. The mobility has been computed for a particle with a radius
$ a = 0.1\delta$
located at a height
$ h = 0.5\delta$
. Continuous red and blue lines indicate the predictions of the MBF theory. The black vertical line indicates the value of
$ H$
at which
$ h + a = 2H$
. When
$ H$
is smaller than this value, a fraction of the particle is outside the box.
Appendix G. Calculation of the self-mobility of a Gaussian kernel
Equation (6.12) can be solved analytically by performing the integration in spherical coordinates:
\begin{align} k_x &= k \cos (\varphi) \sin (\theta), \nonumber \\ k_y &= k \sin (\varphi) \sin (\theta), \nonumber \\ k_z &= k \cos (\theta), \end{align}
where the azimuthal angle
$\varphi \in [0, 2\pi]$
and the polar angle
$\theta \in [0, \pi]$
. We also know that the Fourier transform of the Gaussian kernel is another Gaussian function:
\begin{equation} \hat {S}({\textbf k}) = e^{\frac {k^2\sigma ^2}{2}} .\end{equation}
In this way, the integral in (6.12) can be rewritten as
\begin{equation} {\mathcal{\boldsymbol M}}_{\infty } = \frac {1}{(2\pi)^3} \int ^\infty _0 \int ^\pi _0 \int ^{2\pi }_0 \frac {e^{k^2\sigma ^2}}{\eta (\alpha ^2+k^2)} \left ({\textbf I}-\hat {{\textbf k}} \hat {{\textbf k}}\right) k^2 \sin (\theta) {\mathrm d}\varphi {\mathrm d}\theta {\mathrm d}k. \end{equation}
First, we note that all the off-diagonal terms are zero because the integrals in
$\varphi$
that arise from the products
$k_\alpha k_\beta$
with
$\alpha \neq \beta$
vanish:
$\int _0^{2\pi } \sin (\varphi) \cos (\varphi) {\mathrm d}\varphi = \int _0^{2\pi } \sin (\varphi) {\mathrm d}\varphi = \int _0^{2\pi } \cos (\varphi) {\mathrm d}\varphi = 0$
. To compute the diagonal integrals, we start by calculating the angular integrals:
\begin{eqnarray} I_{xx} = \int ^\pi _0 \int ^{2\pi }_0 \left (1-\sin ^2(\theta)\cos ^2(\varphi)\right) \sin (\theta) {\mathrm d}\varphi {\mathrm d}\theta = \frac {8\pi }{3}, \end{eqnarray}
\begin{eqnarray} I_{yy} = \int ^\pi _0 \int ^{2\pi }_0 \left (1-\sin ^2(\theta)\sin ^2(\varphi)\right) \sin (\theta) {\mathrm d}\varphi {\mathrm d}\theta = \frac {8\pi }{3}, \end{eqnarray}
\begin{eqnarray} I_{zz} = \int ^\pi _0 \int ^{2\pi }_0 \left (1-\cos ^2(\theta)\right) \sin (\theta) {\mathrm d}\varphi {\mathrm d}\theta = \frac {8\pi }{3}. \end{eqnarray}
To compute these integrals, we have used
$\int _0^\pi \sin (\theta) {\mathrm d}\theta = 2$
,
$\int _0^\pi \sin ^3(\theta) {\mathrm d}\theta = {4}/{3}$
and
$\int _0^{2\pi } \sin ^2(\varphi) {\mathrm d}\varphi = \int _0^{2\pi } \cos ^2(\varphi) {\mathrm d}\varphi = \pi$
. In this way, we obtain that the self-mobility can be written as
${\mathcal{\boldsymbol M}}_{\infty } = {\mathcal{M}}_{\infty } {\textbf I}$
, where
\begin{equation} {\mathcal{M}}_{\infty } = \frac {1}{(2\pi)^3 \eta } \cdot \frac {8\pi }{3} \int _0^{\infty } e^{k^2\sigma ^2} \frac {k^2}{\alpha ^2 + k^2} {\mathrm d}k. \end{equation}
The fraction inside the integral can be expressed as
$({k^2}/{\alpha ^2 + k^2}) = 1 - ({\alpha ^2}/{\alpha ^2 + k^2})$
, so we can rewrite the integral as
\begin{equation} {\mathcal{M}}_{\infty } = \frac {1}{3\eta \pi ^2} \left [\int _0^{\infty } e^{-k^2\sigma ^2} {\mathrm d}k - \alpha ^2 \int _0^{\infty } \frac {k^2}{\alpha ^2 + k^2} {\mathrm d}k\right]. \end{equation}
The first integral is the Gaussian integral, whose value is known to be
${\sqrt {\pi }}/{2\sigma }$
. The second integral also has an analytical form (Gradshteuin et al. Reference Gradshteuin, Ryzhik, Jeffrey and Zwillinger2007, integral 3.466)
\begin{equation} \int _0^{\infty } \frac {k^2}{\alpha ^2 + k^2} {\mathrm d}k = \frac {\pi }{2\alpha } e^{\sigma ^2\alpha ^2} \text{erfc}(\sigma \alpha), \end{equation}
where
$\text{erfc}(x) = 1 - {2}/{\sqrt {\pi }} \int ^x_0 e^{-t^2} {\mathrm d}t$
is the complementary error function. Using this result and rearranging the terms, we arrive at
\begin{equation} {\mathcal{\boldsymbol M}}_{\infty } = \frac {1}{6\pi ^{3/2} \eta \sigma } \left [1 - \sqrt {\pi } \alpha \sigma \cdot \text{erfc}(\alpha \sigma) e^{\alpha ^2\sigma ^2}\right] {\textbf I} \end{equation}
and we recall that the effective radius is
$a=\sqrt {\pi } \sigma$
.
Appendix H. Calculation of the smeared self-reaction field
The calculation of the smeared self-reaction field involves solving a four-dimensional integral, as given in (6.15). Here, the tensor
$\boldsymbol {{\mathcal{\hat {R}}}}({\textbf k}, z, z^{\prime })$
is defined as
$\boldsymbol {{\mathcal{\hat {R}}}}({\textbf k}, z, z^{\prime }) = \boldsymbol {{\mathcal{G}}}({\textbf k}, z, z^{\prime }) - \boldsymbol {{\mathcal{G}}}^{(\infty)}({\textbf k}, z, z^{\prime })$
. Analytical expressions for all the components are available in the literature (Felderhof Reference Felderhof2005):
\begin{align} \mathcal{G}_{xx}({\textbf k},z,z^{\prime }) & = \frac {1}{8\pi ^2 \eta k \alpha ^2} \left[\frac {k}{k-s} (k (k+s) e^{-k z^{\prime }-k z}-k(k-s) e^{-k|z-z^{\prime }|} \right . \nonumber\\[3pt] & \quad \left . + s (k-s) e^{-s|z-z^{\prime }|}+ s(k+s) e^{-sz^{\prime }-sz} +2k s e^{-kz^{\prime }-sz} -2 k s e^{-sz^{\prime }-kz}) \right . \nonumber\\[3pt] & \quad \left . + \frac {k_y^2 (k+s)}{\alpha ^2 s} ( s (k+s) e^{-kz^{\prime }-kz} - s (k-s) e^{-k|z-z^{\prime }|}+ k (k-s) e^{-s|z-z^{\prime }|} \right . \nonumber\\[3pt] & \quad \left . - (k^2-ks-2s^2) e^{-sz^{\prime }-sz}-2 s^2 e^{-kz^{\prime }-sz}-2s^2 e^{-sz^{\prime }-kz})\right], \end{align}
\begin{equation} \mathcal{G}^{(\infty)}_{xx}({\textbf k}, z, z^{\prime }) = \frac {1}{8\pi ^2 \eta k s \alpha ^2} \left [ -k_x^2 s e^{-k |z - z^{\prime }|} + k(s^2 - k_y^2) e^{-s |z - z^{\prime }|} \right], \end{equation}
\begin{align} \mathcal{G}_{zz}({\textbf k},z,z^{\prime }) & = \frac {k}{8\pi ^2 \eta s (k-s) \alpha ^2} (s(k + s) e^{-k z^{\prime } - k z} + s(k - s) e^{- k|z-z^{\prime }|} \nonumber\\[4pt] & \quad - k(k - s) e^{-s|z-z^{\prime }|} + k ( k + s) e^{-sz^{\prime }-sz} \nonumber\\[4pt] & \quad - 2ks e^{-kz^{\prime }-sz} - 2ks e^{-sz^{\prime }-kz}), \end{align}
\begin{equation} \mathcal{G}^{(\infty)}_{zz}({\textbf k}, z, z^{\prime }) = \frac {k}{8\pi ^2 \eta s \alpha ^2} \left [ s e^{-k |z - z^{\prime }|} - k e^{-s |z - z^{\prime }|} \right],\end{equation}
with
$s = \sqrt {\alpha ^2+k^2}$
. Despite the availability of these analytical expressions, obtaining a closed-form expression for the integral in (6.15) is challenging. This complexity arises from the repeated integration of products involving fractions, polynomials and exponentials.
Given this difficulty, we opted to compute the integral numerically. However, direct numerical integration in four dimensions is computationally expensive. To address this, we first performed one of the integrals analytically, then evaluated the remaining integrals numerically. Specifically, we expressed the
${\textbf k}$
-integral in polar coordinates:
\begin{align} \mathcal{R}^{\alpha \alpha }_S(\mathrm{h}) & = \int ^\infty _0 S_z(z - \mathrm{h}) \, {\mathrm d}z \int ^\infty _0 S_z(z^{\prime } - \mathrm{h}) \, {\mathrm d}z^{\prime } \int ^\infty _0 |\hat {S}_{\parallel }(k)|^2 k \, {\mathrm d}k \nonumber\\[4pt] & \quad \times \int ^{2\pi }_0 \left ( \mathcal{G}_{\alpha \alpha }(k, \varphi, z, z^{\prime }) - \mathcal{G}^{(\infty)}_{\alpha \alpha }(k, \varphi, z, z^{\prime }) \right){\mathrm d}\varphi . \end{align}
The integral over the azimuthal angle
$\varphi$
is straightforward to solve analytically. For the
$zz$
component, there is no dependence on
$\varphi$
, so the integral simply evaluates to
$2\pi$
. In the case of the
$xx$
component, some terms involve
$k_x^2$
and
$k_y^2$
, which depend on
$\sin ^2(\varphi)$
or
$\cos ^2(\varphi)$
. The integrals of these terms are
$\pi$
, while the integrals of all other terms are
$2\pi$
.
After solving the angular integral analytically, the remaining three integrals were computed numerically using the ‘tplquad’ function from Python’s SciPy library.
We have validated the calculation of this integrals by comparing the numerical results computed for different particle radii with the analytical predictions for point particles (Felderhof Reference Felderhof2005),
\begin{align} \mathcal{R}^{(p)}_{\parallel }(h) & = -\frac {1}{192 \pi \eta h \nu ^4} \left[36 + 27 \nu + 6 \nu ^2 + 6 (6 + 12 \nu + 11 \nu ^2 + 6 \nu ^3 + 4 \nu ^4 ) e^{-2\nu } \right. \nonumber\\[3pt] &\quad - \left(144 + 144 \nu + 72 \nu ^2 + 24 \nu ^3 - 6 \nu ^4 + 2 \nu ^5 - \nu ^6 + \nu ^7 \right) e^{-\nu } \nonumber\\[3pt] & \quad + \nu ^8 E_1(\nu) + 12 \nu ^3 (2 \nu K_2(2\nu) + 3 K_3(2\nu)) \nonumber\\[3pt] &\quad \left. + 6 \pi \nu ^3 ( 2 \nu Y_2(2\nu) - 3 Y_3(2\nu) - 3 {\mathbf H}_{-3}(2\nu) - 2\nu {\mathbf H}_{-2}(2\nu) ) \right], \end{align}
\begin{align} \mathcal{R}^{(p)}_{\perp }(h) & = -\frac {1}{96 \pi \eta h \nu ^4} \left[6(6 + \nu ^2 + 8\nu ^3) + 6(6 + 12\nu + 9\nu ^2 + 2\nu ^3) e^{-2\nu } \right.\nonumber\\[3pt] &\quad - (144 + 144\nu + 48\nu ^2 - 6\nu ^4 + 10\nu ^5 + \nu ^6 - \nu ^7) e^{-\nu } \nonumber\\[3pt] &\quad + \nu ^6 (12 - \nu ^2) E_1(\nu) + 72 \nu ^2 K_0(2\nu) \nonumber\\[3pt] &\quad + 36\nu (2 + \nu ^2) K_1(2\nu) + 6\pi \nu ( 2\nu (3 - \nu ^2)(Y_0(2\nu) - {\mathbf H}_0(2\nu)) \nonumber\\ &\quad \left . - (6 - 5\nu ^2)(Y_1(2\nu) - {\mathbf H}_1(2\nu)) )\right] .\end{align}
In figure 7, we compare numerical calculations for a blob particle of radius
$a$
with the analytical predictions for point particles based on Felderhof’s reaction-field tensor against the particle height
$h$
. As the particles move further away from the bottom of the box, numerical results converge to the point-particle analytical prediction. However, significant discrepancies arise once the particle becomes close to the wall,
$h\lt 2 a$
.
Open access
