# The slow viscous flow around doubly-periodic arrays of infinite slender
cylinders

Lyndon<EMAIL_ADDRESS>Department of Mathematics, University of
Hull, Hull HU6 7RX, United Kingdom  Rohan Vernekar  Univ. Grenoble Alpes,
CNRS, LRP, 38000 Grenoble, France  Timm Krüger  School of Engineering,
Institute for Multiscale Thermofluids, University of Edinburgh, Edinburgh EH9
3FB, United Kingdom  Maciej Lisicki  Institute of Theoretical Physics, Faculty
of Physics, University of Warsaw, Warsaw, Pasteura 5, 02-093 Warsaw, Poland
David W. Inglis  School of Engineering, Faculty of Science and Engineering,
Macquarie University, Sydney, NSW 2109, Australia

###### Abstract

The slow viscous flow through a doubly-periodic array of cylinders does not
have an analytical solution. However, as a reduced model for the flow within
fibrous porous media, this solution is important for many real-world systems.
We asymptotically determine the flow around a doubly-periodic array of
infinite slender cylinders, by placing doubly-periodic two-dimensional
singularity solutions within the cylinder and expanding the no-slip condition
on the cylinder’s surface in powers of the cylinder radius. The asymptotic
solution provides a closed-form estimate for the flow and forces as a function
of the radius and the dimensions of the cell. The force is compared to results
from lattice-Boltzmann simulations of low-Reynolds-number flows in the same
geometry, and the accuracy of the no-slip condition on the surface of the
cylinder, predicted by the asymptotic theory, is checked. Finally, the
behaviour of the flow, flux, force and effective permeability of the cell is
investigated as a function of the geometric parameters. The structure of the
asymptotic permeability is consistent with other models for the flow parallel
to an array of rods. These models could be used to help understand the flows
within porous systems composed of fibres and systems involving periodic arrays
such as deterministic lateral displacement.

##  I Introduction

The slow viscous flow over multiple bodies is a notoriously complicated
problem with various applications [1]. For example, flowing colloids display
discontinuous shear thickening [2, 3], and programmable self-assembling
micromachines interact with each other through the flow to develop distinct
phases and shapes [4, 5]. Similarly, filter-feeding organisms use the flow
over a collection of microscopic fibres to capture food [6, 7, 8], and
periodic arrays of posts in microfluidics can sort particles in a process
called deterministic lateral displacement [9, 10, 11, 12].

These systems can be tricky to probe experimentally and are hard to model
theoretically because slow viscous flows have long-ranged interactions [1].
Even the relatively simplified geometries of a singularly or doubly-periodic
array of infinite cylinders have no exact solutions. Yet, such model arrays
have been studied since the late 1950s [13] and have played an important role
in understanding ordered or fibrous porous media [14]. Fibrous porous media,
like wool, hair, collagen, and fibreglass, can exist at much lower packing
fractions (below 1%) than granular porous media (60-70%) due to the large
aspect ratios of the bodies.

The lack of exact closed solutions means the dynamics in singularly or doubly-
periodic arrays of infinite cylinders are typically solved numerically [15,
16, 17, 18, 19, 20, 21, 22] or argued from geometries with known solutions
[23]. Through these methods, studies have explored the behaviour around
cylinders in different periodic domains [16, 18, 23], irregular domains [20],
the influence of slip [17], cylinder porosity [18], interacting fluid domains
[21, 22], and inertial effects [15, 13].

Singularly periodic arrays have also been studied asymptotically. Such work
provides a closed approximation to the solution that can be used when
numerical approaches struggle, provide insight into how the geometry
influences the behaviour and can be directly applied to new problems. Tamada
and Fujikawa [13] studied the drag on the cylinder when the periodic domain
was much greater than the cylinder radius, while Keller [24] considered the
behaviour in the lubrication limit. Barta and Weihs [19] later used slender-
body theory [25, 26, 27] to improve the accuracy and investigate the effects
of finite lengths and array size.

In contrast, relatively little asymptotic work has been done for doubly-
periodic arrays. In the mid-1980s, Drummond and Tahir estimated the
permeability of different periodic arrays through the matching of the flow
outside an infinite cylinder and a collection of singularities [28]. Their
tests suggest that the permeability, $k^{\prime}$, for an array of cylinders
typically has the form

$\frac{\mu k^{\prime}}{R^{\prime
2}}=\frac{1}{4\phi}\left(-\log(\phi)+\alpha-\beta\phi^{2}+2\phi\right)$ (1)

where $R^{\prime}$ is the radius of the cylinder, $\mu$ is the dynamic
viscosity, $\phi$ is the packing fraction (the volume occupied by the cylinder
divided by the total volume of the cell), and $\alpha$ and $\beta$ are
constants that depend on the geometry of the cell [28, 14]. The values of
$\alpha$ and $\beta$ need to be determined numerically for each geometry
considered. Drummond and Tahir’s models have been found to match several
experimental results [14]. Wang [16] later used the general solution to the
flow in a periodic box and the solution for flow outside a cylinder to
investigate the flow in the doubly-periodic arrays. They numerically enforced
the boundary conditions to the flow at discrete locations on the edge of the
domain to determine the unknown coefficients for each geometry they
investigated. Consistent with Eq. (1), they showed the leading logarithmic
behaviour of the drag and the permeability for doubly-periodic arrays in the
limit of small radius, but they also observed non-inertial vortices in front
of and behind the cylinder when the radius increases. However, to the authors’
knowledge, no closed-form asymptotic representation for the force and the flow
in a doubly-periodic domain has been proposed to date.

This paper determines the asymptotic flow around and the force from a doubly-
periodic array of cylinders, with dimensions $\ell^{\prime}\times h^{\prime}$,
in the limit of small cylinder radius, $R^{\prime}$. The solution is found by
constructing a complex-variable two-dimensional singularity representation for
the flow and satisfying the no-slip boundary condition on the cylinder surface
up to $\mathcal{O}(R^{\prime 4}/\ell^{\prime 4},R^{\prime 4}/H^{\prime 4})$.
The expansion provides a closed-form estimate for the flow and force in the
cell in terms of the scaled radius of the cylinder and the aspect ratio of the
two sides of the domain. The asymptotic force is compared to the force found
from lattice-Boltzmann simulations for the same geometry at small Reynolds
number, and the accuracy of the no-slip condition on the surface of the
cylinder, as predicted by the asymptotic model, is investigated. The
asymptotic solution allows us to create an analytical approximation for the
mean velocity through the cell, the pressure drop across the cell, and the
permeability of the system.

Section II introduces the doubly-periodic cylinder geometry considered. Sec.
III provides some background into the complex variable solutions of Stokes
flow and introduces the doubly-periodic singularities used. These
singularities are then employed in Sec. IV to determine the asymptotic flow
around the periodic cylinder system. The accuracy of these results is
investigated in Sec. V. Finally, the flows predicted are discussed and the
asymptotic permeability is determined in Sec. VI, before concluding the paper
in Sec. VII.

##  II The infinite cylinder in a doubly-periodic cell

Figure 1: A diagram of the doubly-periodic domain considered. The period along
the $x^{\prime}$-axis is $\ell^{\prime}$, and the period along the
$y^{\prime}$-axis is $h^{\prime}$. Any point in the domain on the complex
plane is represented by $z^{\prime}=x^{\prime}+\text{i}y^{\prime}$. The
cylinder has a radius of $R^{\prime}$ and is centred at $z_{0}^{\prime}$.
$\Delta h^{\prime}=h^{\prime}-2R^{\prime}$ is the gap between two adjacent
cylinders along the $y^{\prime}$-axis. The flow is taken along the
$x^{\prime}$-axis without any loss of generality, because of the linearity of
Stokes flow.

This paper considers the slow viscous flow over an infinite slender cylinder
of radius $R^{\prime}$ in a doubly-periodic domain with periods
$\ell^{\prime}$ and $h^{\prime}$, respectively (Fig. 1). The cylinder is
assumed to be stationary, and the background flow will be taken along the
$x$-axis. The dynamic viscosity of the fluid is $\mu$. In this geometry, the
gap between two adjacent cylinders along the $y$-axis is given by $\Delta
h^{\prime}=h^{\prime}-2R^{\prime}$. We scale all lengths by $\ell^{\prime}$,
velocities by the maximum velocity along the $x$-axis,
$u_{\text{m}}^{\prime}$, the force per unit length by $\mu
u_{\text{m}}^{\prime}$, the pressure by $\mu
u_{\text{m}}^{\prime}/\ell^{\prime}$, and the permeability by $\ell^{\prime
2}/\mu$. We note that, for the two-dimensional system, only the force per unit
length is defined.

The flow within the unit cell is assumed to satisfy the incompressible Stokes
equations

$\displaystyle-\nabla p+\mu\nabla^{2}\mathbf{u}$ $\displaystyle=$
$\displaystyle\mathbf{0},$ (2) $\displaystyle\nabla\cdot\mathbf{u}$
$\displaystyle=$ $\displaystyle 0$ (3)

where $p$ is the pressure, and $\mathbf{u}$ is the velocity. The flow must
also be periodic along the $x$\- and $y$-axis, and the no-slip condition
demands that the velocity should be zero on the surface of the cylinder.

As the geometry is two-dimensional, it is useful to solve for the flow on the
complex plane (Fig. 1). Any point in the the scaled domain on the complex
plane can be described by $z=x+\text{i}y$, where $x\in[0,1)$ and $y\in[0,h)$.
The surface of the cylinder is given by

$S=z_{0}+R\text{e}^{\text{i}\theta}$ (4)

where $h=h^{\prime}/\ell^{\prime}$ is the aspect ratio of the cell,
$R=R^{\prime}/\ell^{\prime}$ is the scaled cylinder radius, $z_{0}$ is the
location of the centre of the cylinder, and $\theta$ is the polar angle from
the real axis defined at the centre of the cylinder. The scaled gap between
adjacent cylinders along the $y$-axis is $\Delta h=\Delta
h^{\prime}/\ell^{\prime}=h-2R$. In the scaled coordinates,
$R\in[0,\min(1/2,h/2))$ and $R/h\in[0,\min(1/(2h),1/2))$.

##  III Background in complex solutions to Stokes flow

The incompressible Stokes equations are linear and time-independent. Solutions
to these equations only depend on the instantaneous configuration of the
system and can be constructed by the superposition of several flows. Solutions
to the Stokes equations, with an appropriate set of boundary conditions, are
also known to be unique. Even so, exact solutions to the Stokes equations are
only known in relatively simple geometries. As such, many flows are
approximated using numerical or asymptotic methods. These approaches often
exploit the Green’s function solution to the flow, called the Stokeslet. The
Stokeslet represents the flow from a point force and can be used to construct
solutions in two ways: the boundary integral method and the representation by
fundamental singularities. On the one hand, the boundary integral method uses
the properties of the Green’s function to convert the partial differential
equations into an integral equation over the boundaries of the domain. These
equations can then be inverted numerically to determine the solution. On the
other hand, the representation by fundamental singularities places the
Stokeslet and its derivatives outside the flow domain such that the boundary
conditions are satisfied [29]. In principle, such a distribution must exist
[1] and any suitable distribution must form the solution due to the uniqueness
of the flow. The latter approach is often employed to find asymptotic
solutions to common flow problems, such as the hydrodynamics of fibres [25,
26, 27].

In two dimensions, there is no solution to the Stokes equations for a point
force in an unbounded domain [1]. This so-called Stokes paradox is caused by
the flow in two dimensions growing as $\log r$ with the distance $r$ from the
point force. Solutions to Stokes flow in two dimensions, therefore, only exist
in bounded domains or force-free unbounded domains. In such a domain, the flow
can always be expressed in terms of a stream function, $\psi$, which is
related to the flow velocity through

$u=\frac{\partial\psi}{\partial y},\quad v=-\frac{\partial\psi}{\partial x}$
(5)

where $\mathbf{u}=(u,v)$. The stream function satisfies the equation

$\nabla^{4}\psi=0$ (6)

and has a general solution of the form

$\psi=\Im[\bar{z}f(z)+g(z)]$ (7)

where $\Im[f(z)]$ returns the imaginary part of $f(z)$ and the overbar denotes
the complex conjugate. $f(z)$ and $g(z)$ are analytic functions in the fluid
region and are referred to as Goursat functions [30]. The flow generated from
the above stream function is

$u(z)-\text{i}v(z)=-\bar{f}(z)+\bar{z}\frac{df}{dz}+\frac{dg}{dz}.$ (8)

Though conformal maps do not preserve the boundary conditions for Stokes flow,
this complex representation can be useful in determining the solutions to
various problems [31, 32, 33, 34]. However, no closed map exists for the flow
around a cylinder in a doubly-periodic domain and, thus, the flow around an
array of cylinders cannot currently be solved using conformal maps.

The two-dimensional flow from a point force per unit length in a doubly-
periodic domain was first determined by Hasimoto [35]. Luca and Crowdy [32]
later revisited this problem to determine higher-order singularities and
express it as a rapidly converging series. They showed that the flow from a
two-dimensional point force per unit length of strength $-8\pi
F=-8\pi(F_{x}+\text{i}F_{y})$ located at $z_{0}=x_{0}+\text{i}y_{0}$ in a
doubly-periodic cell with dimensions $x\in[0,1)$ and $y\in[0,h)$ [32] can be
written as

$\displaystyle u(z)-\text{i}v(z)=G_{S}(z-z_{0},F)$ $\displaystyle=$
$\displaystyle-\bar{F}\ln|P(\zeta,\rho)|^{2}+\Re[F]\ln|\zeta|^{2}-F\ln|\zeta|^{2}K(\zeta,\rho)$
$\displaystyle-F\ln(\rho^{2})\rho\frac{\partial\ln
P}{\partial\rho}-\frac{\Re[F]}{2\ln\rho}(\ln|\zeta|^{2})^{2}$

where $\zeta=\exp[2\pi\text{i}(z-z_{0})]$, $\rho=\exp(-2\pi h)$,
$z=x+\text{i}y$ is the location in space and $F_{x}$ ($F_{y}$) represent the
component $F$ in $x$ ($y$) direction. In Eq. (III), $P(\zeta,\rho)$ is the
Schottky–Klein prime function associated with the annulus $\rho<|\zeta|<1$
[32] and is given by

$\displaystyle P(\zeta,\rho)$ $\displaystyle=$
$\displaystyle(1-\zeta)\prod_{k=1}^{\infty}(1-\rho^{k}\zeta)(1-\rho^{k}\zeta^{-1})$
(10) $\displaystyle=$ $\displaystyle A(\rho)s(\zeta,\rho),$ (11)

where

$\displaystyle A(\rho)$ $\displaystyle=$
$\displaystyle\frac{\displaystyle\prod_{n=1}^{\infty}(1+\rho^{n})^{2}}{\displaystyle\sum_{n=1}^{\infty}\rho^{n(n-1)/2}},$
(12) $\displaystyle s(\zeta,\rho)$ $\displaystyle=$
$\displaystyle\sum_{n=-\infty}^{\infty}(-1)^{n}\rho^{n(n-1)/2}\zeta^{n}.$ (13)

Eq. (11) is a rapidly convergent series representation that is useful to
compute the flow [32]. From the Schottky–Klein prime function $P(\zeta,\rho)$
we can define

$\displaystyle K(\zeta,\rho)=\zeta\frac{\partial\ln P}{\partial\zeta}$
$\displaystyle=$ $\displaystyle\frac{\zeta s_{\zeta}}{s},$ (14)
$\displaystyle\rho\frac{\partial\ln P}{\partial\rho}$ $\displaystyle=$
$\displaystyle\frac{\rho A_{\rho}}{A}+\frac{\rho s_{\rho}}{s}$ (15)

where the subscript denotes the derivative with respect to the given variable.
Higher-order singularity solutions can be constructed from the Stokeslet by
taking the appropriate derivatives. For example, the first derivative creates
a force dipole, the second derivative creates a force quadrupole, etc.
Similarly, the flow from a source dipole can also be constructed by taking the
Laplacian of the Stokeslet flow [29].

The symmetries of the doubly-periodic cell mean that, for our model, only
singularities formed by an even number of derivatives of the Stokeslet will
contribute to the flow. Section IV will show that only the Stokeslet, force
quadrupole, source dipole, and source octupole are needed to solve the flow up
to $\mathcal{O}(R^{4},(R/h)^{4})$. These additional singularities are given by

$\displaystyle G_{Q}(z-z_{0},Q)$ $\displaystyle=$
$\displaystyle-\overline{QL(\zeta,\rho)}-2QL-Q\ln|\zeta|^{2}M(\zeta,\rho)-Q\ln(\rho^{2})\rho\frac{\partial
L}{\partial\rho},$ (16) $\displaystyle G_{D}(z-z_{0},D)$ $\displaystyle=$
$\displaystyle-DL(\zeta,\rho),$ (17) $\displaystyle G_{O}(z-z_{0},O)$
$\displaystyle=$ $\displaystyle ON(\zeta,\rho),$ (18)

where $G_{Q}(z-z_{0},Q)$ is the flow from a force quadrupole of the complex
strength $Q=Q_{x}+\text{i}Q_{y}$ located at $z_{0}$, $G_{D}(z-z_{0},D)$ is the
flow from a source dipole of the complex strength $D=D_{x}+\text{i}D_{y}$
located at $z_{0}$, and $G_{O}(z-z_{0},O)$ is the flow from a source octupole
of the complex strength $O=O_{x}+\text{i}O_{y}$ located at $z_{0}$:

$\displaystyle L(\zeta,\rho)=\zeta\frac{\partial K}{\partial\zeta}$
$\displaystyle=$ $\displaystyle\frac{\zeta[s(s_{\zeta}+\zeta
s_{\zeta\zeta})-\zeta s_{\zeta}^{2}]}{s^{2}},$ (19) $\displaystyle
M(\zeta,\rho)=\zeta\frac{\partial L}{\partial\zeta}$ $\displaystyle=$
$\displaystyle\frac{\zeta(s-\zeta s_{\zeta})[s(s_{\zeta}+3\zeta
s_{\zeta\zeta})-2\zeta
s_{\zeta}^{2}]+\zeta^{3}s^{2}s_{\zeta\zeta\zeta}}{s^{3}},$ (20) $\displaystyle
N(\zeta,\rho)=\zeta\frac{\partial M}{\partial\zeta}$ $\displaystyle=$
$\displaystyle\frac{\zeta\left(12\zeta^{2}s\left(s_{\zeta}+\zeta
s_{\zeta\zeta}\right)s_{\zeta}^{2}-\zeta
s^{2}\left(3\zeta^{2}s_{\zeta\zeta}^{2}+7s_{\zeta}^{2}+2\zeta\left(9s_{\zeta\zeta}+2\zeta
s_{\zeta\zeta\zeta}\right)s_{\zeta}\right)\right)}{s^{4}}$ (21)
$\displaystyle+\frac{\zeta\left(s^{3}\left(\zeta\left(\zeta^{2}s_{\zeta\zeta\zeta\zeta}+6\zeta
s_{\zeta\zeta\zeta}+7s_{\zeta\zeta}\right)+s_{\zeta}\right)-6\zeta^{3}s_{\zeta}^{4}\right)}{s^{4}},$
$\displaystyle\rho\frac{\partial L}{\partial\rho}$ $\displaystyle=$
$\displaystyle\frac{\zeta\rho\left(\left(s_{\zeta\rho}+\zeta
s_{\zeta\zeta\rho}\right)s^{2}-\left(2\zeta
s_{\zeta}s_{\zeta\rho}+s_{\rho}\left(s_{\zeta}+\zeta
s_{\zeta\zeta}\right)\right)s+2\zeta s_{\rho}s_{\zeta}^{2}\right)}{s^{3}}.$
(22)

We have left these functions in terms of the derivatives of $s$ because they
will be useful for the series expansion. A more compact representation for
each function can be found in Ref. [32].

##  IV The asymptotic flow around slender cylinders

The flow around a cylinder in a doubly-periodic array can be determined in the
slender limit ($R\ll 1$ and $R\ll h$) using the representation by fundamental
singularities. In this method, singularity solutions are placed within the
cylinder and their strength is chosen such that the boundary conditions are
satisfied [29]. If such a solution can be found, it must be the solution,
because of the uniqueness of Stokes flow.

Inspired by the seminal work of Chwang and Wu [29], we seek a singularity
representation in terms of a point force (Eq. (III)), source dipole (Eq.
(17)), force quadrupole (Eq. (16)) and source octupole (Eq. (18)) placed at
the centre of the cylinder, $z_{0}$. The flow from these singularities can be
written as

$u(z)-\text{i}v(z)=U+G_{S}(z-z_{0},F)+G_{D}(z-z_{0},R^{2}D)+G_{Q}(z-z_{0},R^{4}Q)+G_{O}(z-z_{0},R^{6}O)+\dots$
(23)

where $U$ is the unknown background flow speed in the absence of the cylinder,
$-8\pi F$ is the force per unit length on the fluid from the cylinder
(strength of the Stokeslet), $R^{2}D$ is the strength of the source dipole,
$R^{4}Q$ is the strength of the force quadrupole, and $R^{6}O$ is the strength
of the source octupole. The $R$ scaling in the strengths of the singularity is
chosen to simplify the analysis, and the dots represent the higher
singularities we are not considering. The proposed flow satisfies the doubly-
periodic nature of the cell, so only the no-slip condition at the surface of
the cylinder remains to be satisfied.

At the surface of the cylinder, the fluid velocity is zero and Eq. (23)
becomes

$0=U+G_{S}(R\text{e}^{\text{i}\theta},F)+G_{D}(R\text{e}^{\text{i}\theta},R^{2}D)+G_{Q}(R\text{e}^{\text{i}\theta},R^{4}Q)+G_{O}(R\text{e}^{\text{i}\theta},R^{6}O)+\dots$
(24)

Eq. (24) needs to be solved to determine the unknown strengths of each of the
singularities. Although a general solution cannot be obtained exactly, it is
possible to find a solution in the limit of small scaled radius ($R\ll 1$ and
$R\ll h$). In the small-$R$ limit, the singularities can be expressed as

$\displaystyle G_{S}(R\text{e}^{\text{i}\theta},F)$ $\displaystyle=$
$\displaystyle-F\left(G_{S}^{(0,0)}+R^{2}G_{S}^{(2,0)}\right)+F\text{e}^{-2\text{i}\theta}\left(1+R^{2}G_{S}^{(2,-2)}\right)$
(25)
$\displaystyle+FR^{2}\text{e}^{2\text{i}\theta}G_{S}^{(2,2)}+\mathcal{O}\left(R^{4},\frac{R^{4}}{h^{4}}\right),$
$\displaystyle G_{D}(R\text{e}^{\text{i}\theta},R^{2}D)$ $\displaystyle=$
$\displaystyle-D\frac{\text{e}^{-2\text{i}\theta}}{4\pi^{2}}-D\frac{R^{2}}{12}B(\rho)+\mathcal{O}\left(R^{4},\frac{R^{4}}{h^{4}}\right),$
(26) $\displaystyle G_{Q}(R\text{e}^{\text{i}\theta},R^{4}Q)$ $\displaystyle=$
$\displaystyle-Q\frac{R^{2}\text{e}^{2\text{i}\theta}}{4\pi^{2}}-Q\frac{R^{2}\text{e}^{-4\text{i}\theta}}{2\pi^{2}}+\mathcal{O}\left(R^{4},\frac{R^{4}}{h^{4}}\right),$
(27) $\displaystyle G_{O}(R\text{e}^{\text{i}\theta},R^{6}O)$ $\displaystyle=$
$\displaystyle-O\frac{3R^{2}\text{e}^{-4\text{i}\theta}}{8\pi^{4}}+\mathcal{O}\left(R^{4},\frac{R^{4}}{h^{4}}\right),$
(28)

where

$\displaystyle G_{S}^{(0,0)}$ $\displaystyle=$ $\displaystyle
1+\ln\left[4\pi^{2}R^{2}A^{2}(\rho)s_{\zeta}^{2}(1,\rho)\right]+2\rho\ln\rho\left[\frac{A^{\prime}(\rho)}{A(\rho)}+\frac{s_{\zeta,\rho}(1,\rho)}{s_{\zeta}(1,\rho)}\right],$
(29) $\displaystyle G_{S}^{(2,0)}$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}}{3}\left[B(\rho)+\frac{12}{\ln\rho}\right],$ (30)
$\displaystyle G_{S}^{(2,-2)}$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}}{6}\left[B(\rho)+\frac{12}{\ln\rho}\right],$ (31)
$\displaystyle G_{S}^{(2,2)}$ $\displaystyle=$
$\displaystyle\frac{\pi^{2}}{6}\left[3B(\rho)+\frac{12}{\ln\rho}+4\rho\ln(\rho)\frac{d}{d\rho}\left(\frac{s_{\zeta,\zeta,\zeta}(1,\rho)}{s_{\zeta}(1,\rho)}\right)\right],$
(32) $\displaystyle B(\rho)$ $\displaystyle=$ $\displaystyle
1-\frac{s_{\zeta,\zeta,\zeta,\zeta}(1,\rho)}{s_{\zeta}(1,\rho)}$ (33)

and we have used the properties $s(1,\rho)=0$, $s_{\zeta,\zeta}(1,\rho)=0$,
$s_{\zeta,\zeta,\zeta,\zeta}(1,\rho)+4s_{\zeta,\zeta,\zeta}(1,\rho)=0$. These
properties are proven by separating the summations over $n$ into even and odd
terms and noticing that the summation of the even terms is the negative of the
summation of the odd terms when $\zeta=1$. The above series representation
expresses the singularities as powers of $\text{e}^{\text{i}\theta}$. Since
different powers of $\text{e}^{\theta{i}\theta}$ are orthogonal, Eq. (24)
decomposes into four linear simultaneous equations for $F$, $D$, $Q$ and $O$.
The solutions of these equations give us

$\displaystyle F$ $\displaystyle=C_{F,U}U$
$\displaystyle=\frac{3}{3G_{S}^{(0,0)}+R^{2}\left[\pi^{2}B(\rho)+3G_{S}^{(2,0)}\right]}U+\mathcal{O}\left(R^{4},\frac{R^{4}}{h^{4}}\right),$
(34) $\displaystyle D$ $\displaystyle=C_{D,F}F$
$\displaystyle=4\pi^{2}\left(1+R^{2}G_{S}^{(2,-2)}\right)F+\mathcal{O}\left(R^{4},\frac{R^{4}}{h^{4}}\right),$
(35) $\displaystyle Q$ $\displaystyle=C_{Q,F}F$
$\displaystyle=4\pi^{2}G_{S}^{(2,2)}F+\mathcal{O}(R^{2}),$ (36) $\displaystyle
O$ $\displaystyle=C_{O,Q}Q$
$\displaystyle=-\frac{4\pi^{2}}{3}Q+\mathcal{O}(R^{2})$ (37)

where $C_{F,U}$, $C_{D,F}$, $C_{Q,F}$, and $C_{O,Q}$ are the linearity
coefficients that give the first subscript in terms of the second. The
coefficient relating the Stokelet strength and the background velocity,
$C_{F,U}$, is a scaled drag coefficient for the cylinder and behaves as
$1/(\ln(R)+c)$, where $c$ is a constant, to leading order in $R$. The drag on
slender rods often displays a similar $1/(\ln(R)+c)$ structure [36, 37, 29]
and is consistent with the limiting behaviour in previous studies of cylinders
in periodic arrays [28, 16].

Finally, the strength of $U$ is found by setting the maximum velocity in the
cell to 1. Conservation of mass dictates that the maximum velocity must lie on
the edge of the domain directly above the post, at $z-z_{0}=\text{i}h/2$.
Hence $U$ satisfies

$\displaystyle 1$ $\displaystyle=$ $\displaystyle
U+G_{S}\left(\frac{\text{i}h}{2},C_{F,U}\right)U+G_{D}\left(\frac{\text{i}h}{2},R^{2}C_{D,F}C_{F,U}\right)U+G_{Q}\left(\frac{\text{i}h}{2},R^{4}C_{Q,F}C_{F,U}\right)U$
(38)
$\displaystyle+G_{O}\left(\frac{\text{i}h}{2},R^{6}C_{O,Q}C_{Q,F}C_{F,U}\right)U+\mathcal{O}\left(R^{4},\frac{R^{4}}{h^{4}}\right)$

or

$\displaystyle U$ $\displaystyle=$
$\displaystyle\left\\{1+\left[G_{S}\left(\frac{\text{i}h}{2},1\right)+G_{D}\left(\frac{\text{i}h}{2},C_{D,F}\right)+G_{Q}\left(\frac{\text{i}h}{2},C_{Q,F}\right)\right.\right.$
(39)
$\displaystyle\left.\left.+G_{O}\left(\frac{\text{i}h}{2},C_{O,Q}C_{Q,F}\right)\right]C_{F,U}\right\\}^{-1}+\mathcal{O}\left(R^{4},\frac{R^{4}}{h^{4}}\right)$

where we have used the fact that the singularity strengths are real. Eqs.
(34), (35), (36),(37) and (39) uniquely determine all the unknown coefficients
in the proposed flow representation, Eq. (23), to
$\mathcal{O}(R^{4},R^{4}/h^{4})$. Hence, they complete the asymptotic model.
The asymptotic model can be used to predict the force per unit length from the
cylinder onto the fluid, $-8\pi F$ (Fig. 2a), the background flow in the
absence of the cylinders, $U$ (Fig. 2b), and the flow around the cylinder, Eq.
(23). We note that higher-order solutions can be constructed similarly but
require higher-order singularities to satisfy the boundary conditions.

Figure 2: Contour plots of (a) force per unit length, $-8\pi F$, and (b) the
background flow in the absence of posts, $U$. Dashed grey lines are lines of
constant cell aspect ratio, $h$.

##  V Validation tests of the asymptotic model

The asymptotic force per unit length on the fluid (Fig. 2a) decreases as the
scaled radius, $R$, decreases and the cell aspect ratio, $h$, increases, as
expected. As the cell aspect ratio decreases ($R/h$ increases), the cylinders
become closer together, increasing the hydrodynamic interactions between them
and the drag. Similarly, as the radius increases, the space between the
cylinders decreases and the drag increases.

The accuracy of the asymptotic force per unit length, $-8\pi F$ (Eq. (34)),
was quantified by comparing to the force per unit length determined from
lattice-Boltzmann (LB) simulations for the flow past a single stationary post
in a rectangular domain with periodic boundaries along both the $x$\- and
$y$-axis (Fig. 3) [38]. 200 LB nodes along the $x$-axis were used for all
simulated cases, while the number of nodes along the $y$-axis varied from 26
to 400. The flow was driven by a constant body force along the $x$-axis,
mimicking a constant pressure gradient. The resolution varied from 20 to 80 LB
nodes across the post diameter, with the stationary cylindrical boundary
resolved using a second-order accurate boundary condition [39]. We used the
9-velocity ‘compressible’ LB equilibrium model in the limit of low Mach
numbers ($Ma$) and Reynolds numbers ($Re$), with the two-relaxation time
collision operator [40]. The force on the post is computed with the momentum
exchange algorithm [41, 42]. The resulting flow is practically incompressible
and in the Stokes flow regime (since in all simulated cases, the LB Mach
number $Ma<1.9\times 10^{-5}$ and the Reynolds number $Re<1.7\times 10^{-2}$,
computed with the gap length, $\Delta h=h-2R$, and the maximum velocity
magnitude), and each simulation has been run until the maximum velocity
converges below $\leq 0.01\%$ relative change per time step. We did not see
that the simulation results change significantly when both $Ma$ and $Re$ are
varied by an order of magnitude.

The asymptotic force per unit length found in the asymptotic model, $-8\pi F$
(Eq. (34)), and the LB simulation, $F_{\text{LB}}$, agree well for small
scaled radii, $R$, and larger gaps between adjacent posts along the $y$-axis,
$\Delta h$ (Fig. 3a). The differences are harder to distinguish when plotted
against $R/h$ (Fig. 3b). The relative error, $E$, between the two results,
defined as

$\displaystyle E=\frac{F_{\text{LB}}-(-8\pi F)}{F_{LB}},$ (40)

is shown in Fig. 3c,d. Similarly to direct plots of force, the relative error,
$E$, increases with increasing scaled radius, $R$, and decreases with
increasing cell aspect ratio, $h$. The results suggest that when $R<0.1$ and
$R/h<0.1$ the asymptotic force displays less than 5% error (Eq. (40)). These
findings are consistent with the limits expected from the asymptotic
expansion. We note that for small gap sizes, $\Delta h<0.25$, the relative
error starts to decrease again. These small gap sizes tend to correspond to
$R/h>0.2$ and a rapid increase in $F$, indicating that it is well outside the
region of validity for the model. The apparent improvement is therefore likely
to be a coincidence.

Figure 3: (a,b) Comparison of the asymptotic force per unit length, $-8\pi F$
(lines), and the numerical results of lattice-Boltzmann simulations (dots).
(c,d) Relative error between the asymptotic force and the simulated force. The
relative error is defined by Eq. (40). (a) and (c) are plotted against the gap
between adjacent posts along the $y$-axis, $\Delta h$. b) and c) are plotted
against the scaled radius divided by the aspect ratio of the cell, $R/h$. In
all panels, $R=0.2$ is purple, $R=0.1$ is red, $R=0.085$ is green, $R=0.07$ is
yellow and $R=0.05$ is blue.

Unlike the asymptotic force, the background velocity within the cell without
the cylinder, $U$ (Fig. 2b), is less easy to understand. At small scaled
radii, $R$, the background velocity is positive. As $R$ increases, the
velocity decreases, ultimately changing its sign at some $R$. The rate of
decrease increases with smaller $h$. It is important to recognise that the
background velocity, $U$, is not the mean flow in the cell, but a mathematical
construct necessary to scale the flow by the maximum velocity in the cell. It
can never be observed numerically or experimentally. The mean flow in the cell
comes from the combination of the background velocity and the net flows
generated by the singularities (determined in Sec. VI). Therefore, the
background velocity changes its sign as the strengths of the singularities
involved increase. As the incompressible Stokes equations are linear, an
increase in the strengths of the singularities corresponds to an increase in
the flux they generate. The background flow, therefore, changes its sign to
reduce the mean flow and thereby keep the velocity everywhere within the cell
below 1.

While the background velocity within the cell without the cylinder, $U$, is
not useful for testing the accuracy of the predicted flow, the error on the
asymptotic flow, Eq. (23), can be estimated through the average error on the
no-slip condition on the surface of the cylinder. The average error on the no-
slip condition is defined as

$|\mathbf{u}|_{\theta}=\frac{1}{2\pi}\int_{0}^{2\pi}\left|u(z_{0}+R\text{e}^{\text{i}\theta})-\text{i}v(z_{0}+R\text{e}^{\text{i}\theta})\right|\,d\theta$
(41)

and varies with the scaled radius, $R$, and the cell aspect ratio, $h$. If the
asymptotic solution was exact, $|\mathbf{u}|_{\theta}=0$, since Stokes flows
are unique. Non-zero $|\mathbf{u}|_{\theta}$, therefore, indicates a
difference between the unknown exact solution and the asymptotic prediction.
The size of the difference must be proportional to $|\mathbf{u}|_{\theta}$
since Stokes flows are linear. Hence, small $|\mathbf{u}|_{\theta}$
corresponds to a small difference between the exact and the asymptotic flows.
Contour plots of the average error on the no-slip condition are shown in Fig.
4. The average error on the no-slip condition increases with increasing scaled
cylinder radius, $R$, and decreasing cell aspect ratio, $h$. On logarithmic
axes (Fig. 4b), the increase appears to occur at a constant rate, suggesting
power-law dependence on $R$ and $R/h$. The error of the asymptotic theory is
expected to increase with $R^{4}$ and $(R/h)^{4}$ and so is consistent with
the observed increase. A comparison of the relative error of the force per
unit length, $E$ (Fig. 3c,d), and the average error on the no-slip condition
suggests that when $|\mathbf{u}|_{\theta}<0.01$ the error of the force per
unit length is less than 6%. Similarly to the force, the 1% error appears to
occur around $R\sim 0.1$ and $R/h\sim 0.1$. For the remainder of the paper,
all phase diagrams will consider $R<0.2$ and $R/h<0.2$ to focus on the region
of validity of the asymptotic solution.

Figure 4: (a) Average error of the no-slip condition on the surface of the
cylinder predicted by the asymptotic flow, Eq. (23), as a function of the
scaled radius, $R$, and the scaled radius divided by the cell aspect ratio,
$R/h$. (b) Same as (a), but shown logarithmically. Dashed grey lines are lines
of constant cell aspect ratio, $h$.

##  VI The behaviour of the asymptotic flow

The asymptotic model developed in Sec. IV predicts the strengths of the
singularities (Eqs. (34), (35), (36), and (37)), the background velocity in
the absence of posts (Eq. (39)), and the flow throughout the domain. The
approximate flow at any point in the domain is given by Eq. (23) when the
strengths of the singularities and the background velocity are substituted
into the equation. The streamlines of the asymptotic flow are plotted in Fig.
5 for $z_{0}=(1+\text{i}h)/2$, and scaled radii of $R=0.05$, $0.1$, and $0.2$
and cell aspect ratios of $h=2$, $1$, and $0.5$. The scaled cylinder radius of
$R=0.2$ lies beyond the validity limit of the model, established in the
previous section, and is expected to show an error greater than 5%, especially
when the aspect ratio is $h=0.5$ ($R/h=0.4$). Near each cylinder, the flow
decreases. The size of the region with reduced velocity increases with scaled
radius, $R$, and cell aspect ratio, $h$, due to the periodic interactions. The
flows through cells with large aspect ratios also display a large portion of
the fluid travelling at almost the maximum-velocity cell, while cells with
smaller aspect ratios only reach the maximum velocity in a localised region
directly above the cylinder. The localisation of the maximum velocity is
caused by the conservation of mass which requires a faster velocity in the gap
to squeeze the same amount of fluid through the cell. We note that the flow
along the $x$-axis never reverses anywhere within the flow domain. Therefore,
the asymptotic solution does not predict any closed vortices. Non-inertial
vortices are typically found in Stokes flow at the leading and trailing edges
of cylinders with large radius [16]. The assumption that the cylinders are
slender ($R\ll 1$ and $R/h\ll 1$) in the expansion prevents the appearance of
these vortices in the asymptotic flow.

Figure 5: Streamlines for the flow predicted by the asymptotic solution (Sec.
V) for scaled cylinder radii $R=0.2$, $0.1$, and $0.05$ and cell aspect ratios
$h=2$, $1$, and $0.5$. The colour of the streamlines corresponds to the flow
speed.

The mean velocity in the periodic domain can be determined from the asymptotic
flow by integrating the flow over $y=[0,h)$ for any $x$. The incompressible
condition requires the flux through any plane normal to the $x$-axis to be the
same. Hence the mean fluid velocity is given by

$\displaystyle\langle u\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{h}\int_{0}^{h}u(x+\text{i}y)\,\text{d}y,$ (42)
$\displaystyle=$ $\displaystyle
U+\frac{1}{h}\int_{0}^{h}G_{S}(z-z_{0},F)\,\text{d}y+\frac{1}{h}\int_{0}^{h}G_{D}(z-z_{0},R^{2}D)\,\text{d}y$
$\displaystyle+\frac{1}{h}\int_{0}^{h}G_{Q}(z-z_{0},R^{4}Q)\,\text{d}y+\frac{1}{h}\int_{0}^{h}G_{O}(z-z_{0},R^{6}O)\,\text{d}y+\mathcal{O}\left(R^{4},\frac{R^{4}}{h^{4}}\right)$

where we chose $z_{0}=(1+\text{i}h)/2$. These integrals can be evaluated
exactly using the properties of the Schottky–Klein prime function and its
derivatives to find

$\langle u\rangle=U+\langle G_{S}\rangle F+\langle
G_{D}\rangle(R^{2}D+R^{4}Q)+\mathcal{O}\left(R^{4},\frac{R^{4}}{h^{4}}\right)$
(43)

where

$\displaystyle h\langle G_{S}\rangle$ $\displaystyle=$
$\displaystyle\frac{[\ln(\rho)]^{2}}{12\pi}+\frac{\ln\rho}{2\pi}\ln\left(\frac{P^{2}(-\sqrt{\rho},\rho)}{\sqrt{\rho}}\right)-\frac{\ln(\rho^{2})}{\pi}\sum_{k=1}^{\infty}k\ln\left(\frac{1+\rho^{k-1/2}}{1+\rho^{k+1/2}}\right),$
(44) $\displaystyle h\langle G_{D}\rangle$ $\displaystyle=$
$\displaystyle\frac{2K(-\sqrt{\rho},\rho)-1}{2\pi}$ (45)

are the fluid fluxes from a unit Stokeslet and source dipole, respectively.
The flux from the source dipole is the same as the flux from a force
quadrupole, while no flux is generated from the source octupole. The fluxes
from each singularity depend only on the cell aspect ratio, $h$, through
$\rho=\text{e}^{-2\pi h}$. Hence, the dependence in the mean velocity on the
scaled cylinder radius, $R$, comes from the background flow in the absence of
the post, $U$, the force on the flow from the post, $F$, the source dipole
strength, $R^{2}D$, and the force quadrupole strength, $R^{4}Q$. The
background flow in the absence of the post, $U$, and the force on the flow
from the post, $F$, depend on the scaled radius, $R$, logarithmically at
leading order. Hence the mean velocity is expected to display a weak
dependence on $R$ at leading order in the limit $R\ll 1$ and $R/h\ll 1$.

Figure 6: (a) Mean velocity in the cell, (b) pressure drop across the cell,
and (c) permeability of the cell as a function of the scaled radius, $R$, and
the scaled radius divided by the cell aspect ratio, $R/h$. (d) Same as (c),
but shown logarithmically. Dashed grey lines are lines of constant cell aspect
ratio, $h$.

Although the leading mean velocity in $R$ is dominated by $\ln(R)$, the full
asymptotic form, Eq. (43), shows a rich behaviour (Fig. 6a). For small $R$ and
$R/h$, the mean velocity is seen to increase as $R$ and $R/h$ decrease
(increasing $h$ for constant $R$) because smaller $R$ and $R/h$ increase the
separation between adjacent points, allowing more flow to travel near the
maximum velocity of the cell. However, when $R/h>0.1$ or $R>0.1$ changes to
the scaled radii tends to generate small changes in the mean velocity,
$\langle u\rangle$, as the contours become vertical lines. It is unclear what
the causes of this change are, but it may be related to issues in the
asymptotic models at larger $R$ and $R/h$.

The pressure drop across the cell can also be determined from the asymptotic
model (Fig. 6b). The pressure drop in a doubly-periodic cell is given by
$\Delta p=-4F/h$ where $F$ is the strength of the Stokeslet, Eq. (34).
Similarly to the force, the pressure drop increases as the cell aspect ratio,
$h$, decreases and the scaled radius, $R$, increases. The factor of $1/h$ in
the definition, however, changes the behaviour with $R$, for fixed $R/h$,
because $h$ must go to 0 to keep $R/h$ constant if $R\to 0$ and so the
pressure starts to diverge in this limit. The physical understanding of the
behaviour of the pressure drop is identical to that of the force.

The mean velocity in the cell and the pressure drop across the cell allow us
to calculate the scaled permeability for the doubly-periodic array, $k$:

$\displaystyle k=\frac{\langle u\rangle}{\Delta p}$ $\displaystyle=$
$\displaystyle-\frac{hU+h\langle G_{S}\rangle F+h\langle
G_{D}\rangle(R^{2}D+R^{4}Q)}{4F}+\mathcal{O}\left(R^{4},\frac{R^{4}}{h^{4}}\right)$
(46) $\displaystyle=$ $\displaystyle-\frac{h}{4C_{F,U}}-\frac{h\langle
G_{S}\rangle}{4}-\frac{h\langle
G_{D}\rangle}{4}(R^{2}C_{D,F}+R^{4}C_{Q,F})+\mathcal{O}\left(R^{4},\frac{R^{4}}{h^{4}}\right)$
$\displaystyle=$ $\displaystyle\frac{U}{\Delta p}-\frac{h\langle
G_{S}\rangle}{4}-\frac{h\langle
G_{D}\rangle}{4}(R^{2}C_{D,F}+R^{4}C_{Q,F})+\mathcal{O}\left(R^{4},\frac{R^{4}}{h^{4}}\right).$

The asymptotic permeability of the cell, Eq. (46), is dominated by the flux
from a unit Stokeslet, which depends only on the aspect ratio $h$, and
$U/\Delta p$, which depends on both the scaled radius $R$ and the aspect ratio
$h$. Hence, at very small $R$ ($R\ll 1$ and $R/h\ll 1$), the permeability
changes logarithmically with the radius, but it has a stronger dependence on
the aspect ratio of the cell. Interestingly, if the permeability is written in
terms of the packing fraction $\phi=\pi R^{2}/h$ instead of $R$, it can be
expressed as

$\frac{\mu k^{\prime}}{R^{\prime
2}}=\frac{k}{R^{2}}=\frac{1}{4\phi}\left(-\ln\phi+\alpha(h)+\beta(h)\phi^{2}+\gamma(h)\phi\right)+\mathcal{O}\left(R^{4},\frac{R^{4}}{h^{4}}\right)$
(47)

where

$\displaystyle\alpha(h)$ $\displaystyle=$ $\displaystyle-1-\langle
G_{S}\rangle-\ln\left[4\pi
hA^{2}(\rho)s_{\zeta}^{2}(1,\rho)\right]-2\rho\ln\rho\left[\frac{A^{\prime}(\rho)}{A(\rho)}+\frac{s_{\zeta,\rho}(1,\rho)}{s_{\zeta}(1,\rho)}\right],$
(48) $\displaystyle\gamma(h)$ $\displaystyle=$
$\displaystyle-\frac{hG_{S}^{(2,0)}}{\pi}-\frac{h\pi}{3}\left(B(\rho)+12\langle
G_{D}\rangle\right),$ (49) $\displaystyle\beta(h)$ $\displaystyle=$
$\displaystyle-4\langle
G_{D}\rangle\left(G_{S}^{(2,-2)}+G_{S}^{(2,2)}\right).$ (50)

The above form of the permeability is similar to the suggested form by
Drummond and Tahir for regular arrays [28], except all constants depend on the
aspect ratio of the cell. The apparent difference is likely due to the
different asymptotic approaches taken and the level of accuracy in the
expansions. The asymptotic permeability, found above, therefore gives a closed
form for these constants in doubly-periodic arrays.

The asymptotic permeability, Eq. (46), increases as the cell aspect ratio,
$h$, increases and the scaled radius, $R$, decreases (Fig. 6c,d). On a
logarithmic scale, the lines of constant $\log_{10}(k)$ are almost straight
lines (Fig. 6d), suggesting that $\log_{10}(k)\propto a\log(R)-b\log(R/h)$ or
$k\propto R^{a}(R/h)^{-b}$ where $a$ and $b$ are constants. The increase in
asymptotic permeability with increasing $h$ and decreasing $R$ is due to the
mean velocity in the cell, Eq. (43), having a much weaker dependence on $R$
and $h$ than the pressure drop across the cell, $\Delta p$ (Fig. 6b), over the
region of validity for the asymptotic model. Hence, when the pressure drop is
small, the permeability becomes large.

##  VII Conclusion

This paper asymptotically determines the slow viscous flow around a doubly-
periodic array of cylinders in the limit that the cylinder is slender. The
slender condition means that the scaled radius, $R$, is much smaller than 1
and the cell aspect ratio, $h$. The asymptotic solution is constructed using
the complex singularity solutions to two-dimensional Stokes flow and the
representation by fundamental singularities. The results provide an analytical
representation of the flow and the force per unit length as a function of the
scaled radius of the cylinder, $R$, and the aspect ratio of the domain, $h$,
and are accurate to $\mathcal{O}(R^{4},R^{4}/h^{4})$. The asymptotic force per
unit length on the flow from the cylinder is compared to lattice-Boltzmann
simulations for the same domain, and the accuracy of the asymptotically
applied no-slip condition on the surface of the cylinder is investigated. The
asymptotic results were then used to analyse the behaviour of the flow for
varying scaled cylinder radius, $R$, and cell aspect ratio, $h$, and a closed
representation of the permeability in the slender limit ($R\ll 1,h$) was
developed. In the future, our analysis could be extended beyond
$\mathcal{O}(R^{4},R^{4}/h^{4})$ through the addition of higher-order
singularities. Our results could be of use in modelling the flows within
porous systems composed of fibres and systems involving periodic arrays such
as deterministic lateral displacement.

##  VIII Acknowledgments

T.K. received funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program (803553). The
work of ML was supported by the National Science Centre of Poland grant Sonata
no. 2018/31/D/ST3/02408.

For the purpose of open access, the authors have applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission.

##  IX Data Management

The MATLAB Implementation of asymptotic formulae used within the text and the
data from the LB simulations plotted in fig 3 are available on GitHub [43].

## References

  * [1] S. Kim and S. J. Karrila, Microhydrodynamics: Principles and Selected Applications.  Boston: Courier Corporation, 2005.
  * [2] S. Jamali and J. F. Brady, “Alternative Frictional Model for Discontinuous Shear Thickening of Dense Suspensions: Hydrodynamics,” Phys. Rev. Lett., vol. 123, p. 138002, 2019.
  * [3] M. Wang, S. Jamali, and J. F. Brady, “A hydrodynamic model for discontinuous shear-thickening in dense suspensions,” J. Rheol. (N. Y. N. Y)., vol. 64, pp. 379–394, 2020.
  * [4] W. Wang, G. Gardi, P. Malgaretti, V. Kishore, L. Koens, D. Son, H. Gilbert, Z. Wu, P. Harwani, E. Lauga, C. Holm, and M. Sitti, “Order and information in the patterns of spinning magnetic micro-disks at the air-water interface,” Sci. Adv., vol. 8, p. 685, 2022.
  * [5] L. Koens, W. Wang, M. Sitti, and E. Lauga, “The near and far of a pair of magnetic capillary disks,” Soft Matter, vol. 15, pp. 1497–1507, 2019\.
  * [6] L. T. Nielsen, S. S. Asadzadeh, J. Dölger, J. H. Walther, T. Kiørboe, and A. Andersen, “Hydrodynamics of microbial filter feeding,” Proc. Natl. Acad. Sci., vol. 114, pp. 9373–9378, 2017.
  * [7] A. I. Lavrov, F. V. Bolshakov, D. B. Tokina, and A. V. Ereskovsky, “Fine details of the choanocyte filter apparatus in asconoid calcareous sponges (Porifera: Calcarea) revealed by ruthenium red fixation,” Zoology, vol. 150, p. 125984, 2022.
  * [8] J. R. Blake, S. R. Otto, and D. A. Blake, “Filter Feeding, Chaotic Filtration, and a Blinking Stokeslet,” Theor. Comput. Fluid Dyn., vol. 10, pp. 23–36, 1998.
  * [9] D. W. Inglis, J. A. Davis, R. H. Austin, and J. C. Sturm, “Critical particle size for fractionation by deterministic lateral displacement,” Lab Chip, vol. 6, p. 655, 2006.
  * [10] V. Biagioni, G. Balestrieri, A. Adrover, and S. Cerbelli, “Combining Electrostatic, Hindrance and Diffusive Effects for Predicting Particle Transport and Separation Efficiency in Deterministic Lateral Displacement Microfluidic Devices,” Biosensors, vol. 10, p. 126, 2020.
  * [11] S. Kim, B. H. Wunsch, H. Hu, J. T. Smith, R. H. Austin, and G. Stolovitzky, “Broken flow symmetry explains the dynamics of small particles in deterministic lateral displacement arrays,” Proc. Natl. Acad. Sci., vol. 114, p. 201706645, 2017.
  * [12] M. Jiang, A. D. Mazzeo, and G. Drazer, “Centrifuge-based deterministic lateral displacement separation,” Microfluid. Nanofluidics, vol. 20, p. 17, 2016.
  * [13] K. Tamada and H. Fujikawa, “The stready two-dimensional flow of viscous fluid at low Reynolds numbers passing through an infinite row of equal parallel circular cylinders.,” Q. J. Mech. Appl. Math., vol. 10, pp. 425–432, 1957\.
  * [14] G. W. Jackson and D. F. James, “The permeability of fibrous porous media,” Can. J. Chem. Eng., vol. 64, pp. 364–374, 1986.
  * [15] F. Ayaz and T. J. Pedley, “Flow through and particle interception by an infinite array of closely-spaced circular cylinders,” Eur. J. Mech. - B/Fluids, vol. 18, pp. 173–196, 1999.
  * [16] C. Y. Wang, “Stokes flow through a rectangular array of circular cylinders,” Fluid Dyn. Res., vol. 29, pp. 65–80, 2001.
  * [17] C. Y. Wang, “Stokes slip flow through a grid of circular cylinders,” Phys. Fluids, vol. 14, pp. 3358–3360, 2002.
  * [18] V. A. Kirsh, “Stokes flow past periodic rows of porous cylinders,” Theor. Found. Chem. Eng., vol. 40, pp. 465–471, 2006.
  * [19] E. Barta and D. Weihs, “Creeping flow around a finite row of slender bodies in close proximity,” J. Fluid Mech., vol. 551, p. 1, 2006.
  * [20] D. Shou, L. Ye, and J. Fan, “On the longitudinal permeability of aligned fiber arrays,” J. Compos. Mater., vol. 49, pp. 1753–1763, 2015.
  * [21] A.-T. Tran, H. Le-Quang, Q.-C. He, and D.-H. Nguyen, “Determination of the Effective Permeability of Doubly Porous Materials by a Two-Scale Homogenization Approach,” Transp. Porous Media, pp. 1–47, 2022.
  * [22] A.-T. Tran, H. L. Quang, Q.-C. He, and D.-H. Nguyen, “Solutions for two coupled Stokes and Darcy flows around and in spheroidal inclusions,” J. Eng. Math., vol. 132, p. 19, 2022.
  * [23] M. Maleki, R. J. Martinuzzi, W. Herzog, and S. Federico, “Orthotropic hydraulic permeability of arrays of parallel cylinders,” Phys. Rev. E, vol. 96, p. 033112, 2017.
  * [24] J. B. Keller, “Viscous flow through a grating or lattice of cylinders,” J. Fluid Mech., vol. 18, pp. 94–96, 1964.
  * [25] L. Koens and E. Lauga, “The boundary integral formulation of Stokes flows includes slender-body theory,” J. Fluid Mech., vol. 850, p. R1, 2018.
  * [26] J. B. Keller and S. I. Rubinow, “Slender-body theory for slow viscous flow,” J. Fluid Mech., vol. 75, pp. 705–714, 1976.
  * [27] R. E. Johnson, “An improved slender-body theory for Stokes flow,” J. Fluid Mech., vol. 99, pp. 411–431, 1979.
  * [28] J. Drummond and M. Tahir, “Laminar viscous flow through regular arrays of parallel solid cylinders,” Int. J. Multiph. Flow, vol. 10, pp. 515–540, 1984.
  * [29] A. T. Chwang and T. Y. Wu, “Hydromechanics of low-Reynolds-number flow. Part 2\. Singularity method for Stokes flows,” J. Fluid Mech., vol. 67, pp. 787–815, 1975.
  * [30] W. E. Langlois and M. O. Deville, Slow Viscous Flow.  Cham: Springer International Publishing, 2014.
  * [31] D. Crowdy and E. Luca, “Analytical solutions for two-dimensional singly periodic Stokes flow singularity arrays near walls,” J. Eng. Math., vol. 119, pp. 199–215, 2019.
  * [32] D. Crowdy and E. Luca, “Fast evaluation of the fundamental singularities of two-dimensional doubly periodic Stokes flow,” J. Eng. Math., vol. 111, pp. 95–110, 2018.
  * [33] E. Luca and D. G. Crowdy, “A transform method for the biharmonic equation in multiply connected circular domains,” IMA J. Appl. Math., vol. 83, pp. 942–976, 2018.
  * [34] D. G. Crowdy, “Treadmilling swimmers near a no-slip wall at low Reynolds number,” Int. J. Non. Linear. Mech., vol. 46, pp. 577–585, 2011.
  * [35] H. Hasimoto, “In:,” Lect. Notes Res. Sci., vol. 214, p. 17, 1974.
  * [36] J. Gray and G. J. Hancock, “The Propulsion of Sea-Urchin Spermatozoa,” J. Exp. Biol., vol. 32, pp. 802–814, 1955.
  * [37] J. Lighthill, Mathematical Biofluiddynamics.  SIAM, 1975.
  * [38] R. Vernekar, T. Krüger, K. Loutherback, K. Morton, and D. W. Inglis, “Anisotropic permeability in deterministic lateral displacement arrays,” Lab on a Chip, vol. 17, pp. 3318–3330, Sept. 2017.
  * [39] M. Bouzidi, M. Firdaouss, and P. Lallemand, “Momentum transfer of a Boltzmann-lattice fluid with boundaries,” Physics of Fluids, vol. 13, pp. 3452–3459, Oct. 2001.
  * [40] I. Ginzburg, F. Verhaeghe, and D. d’Humieres, “Two-relaxation-time Lattice Boltzmann scheme: About parametrization, velocity, pressure and mixed boundary conditions,” Communications in computational physics, vol. 3, pp. 427–478, Feb. 2008.
  * [41] A. J. C. Ladd, “Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation,” Journal of Fluid Mechanics, vol. 271, pp. 285–309, July 1994.
  * [42] L.-S. Luo, W. Liao, X. Chen, Y. Peng, and W. Zhang, “Numerics of the lattice Boltzmann method: Effects of collision models on the lattice Boltzmann simulations,” Physical Review E, vol. 83, p. 056710, May 2011.
  * [43] “Numerical implementation of asymptotic formulea in matlab.” https://github.com/LKoens/ADPSC.