3.1. Limit of small $h$

The following analysis applies in the limit where the force is close to the boundary such that $\alpha h \ll 1$ and $h \ll r$ , where $r=\sqrt {\rho ^2+z^2}$ . As in the case of the Stokes fluid and a force perpendicular to the wall (Blake Reference Blake1971; Blake & Chwang Reference Blake and Chwang1974), the zeroth and first-order terms in $h$ vanish. The Green’s function can therefore be expanded to the leading order as

(3.11) \begin{align} \boldsymbol{G} = h^2 \boldsymbol{G}^{(2)} + \mathcal{O} ( h^3) . \end{align}

The same expansion also holds for the pressure $P$ and the function $\varPsi$ . Specifically,

(3.12) \begin{align} P = h^2 P^{(2)} + \mathcal{O} ( h^3) , \quad \varPsi = h^2 \varPsi ^{(2)} + \mathcal{O} ( h^3) , \end{align}

with

(3.13) \begin{align} P^{(2)}=2\int _0^\infty k^2(k+K ) \, \textrm{e}^{-kz}\, J_0(k\rho ) \, \mathrm{d} k \end{align}

and

(3.14) \begin{align} \varPsi ^{(2)}=-2\int _0^\infty k(K+k ) \, \textrm{e}^{-Kz}\, J_0(k\rho ) \, \mathrm{d} k . \end{align}

The first term in each of the integrals can easily be recognised as the second $z$ -derivative of the free-space solution (3.6). The solution therefore consists of a force quadrupole and further corrections resulting from the mixed terms mentioned above. We decompose the solutions as

(3.15) \begin{align} P^{(2)}=P^{(2)}_1 +P^{(2)}_2 , \quad \varPsi ^{(2)}=\varPsi ^{(2)}_1+ \varPsi ^{(2)}_2, \end{align}

where

(3.16) \begin{align} P^{(2)}_1 = \partial _z^2 P^\infty |_{h=0} , \quad \varPsi ^{(2)}_1= \partial _z^2 \varPsi ^\infty |_{h=0} . \end{align}

In the following, all expressions are evaluated at $h=0$ unless $h$ appears explicitly.

3.2. Pressure field

The first contribution to the pressure in (3.15) is obtained as the second derivative of the pressure from (3.3) and reads

(3.17) \begin{align} P^{(2)}_1 = 6 \, \frac {5z^3-3r^2z}{r^7}. \end{align}

It corresponds to an axisymmetric source octupole, which can be expressed in terms of a Legendre polynomial as

(3.18) \begin{align} P_1^{(2)} = \frac {12}{ r^4} \, P_3(z/r) . \end{align}

For the second contribution, we are not aware of a closed-form analytical solution. By substituting $k = \alpha t$ , the second term from the integral in (3.13) can be written as

(3.19) \begin{align} P^{(2)}_2 = 2\alpha ^4 \int _0^\infty t^2 \sqrt {t^2+1} \, \textrm{e}^{-\alpha z t} \, J_0(\alpha \rho t) \, \mathrm{d} t . \end{align}

In the far-field limit, we have $\alpha z \gg 1$ . Therefore, only small values of the integration variable $t$ contribute to the integral. We can expand the expression $\sqrt {t^2+1}$ in a Taylor series

(3.20) \begin{align} \sqrt {t^2+1}=\sum _{n=0}^{\infty } \varphi _n t^{2n}, \end{align}

where

(3.21) \begin{align} \varphi _n = \frac {(-1)^{n+1} (2n-3)!!}{ 2^{n} n! } , \end{align}

with $n!!$ denoting a double factorial. The leading coefficients have the values $\varphi _0 = 1$ , $\varphi _1 = 1/2$ , $\varphi _2 = -1/8$ , $\varphi _3 = 1/16$ , $\varphi _4 = -5/128$ , etc.

Inserting (3.20) into (3.19), each term corresponds to a derivative of the free-space pressure Green’s function $P^\infty$ with respect to $z$ . This yields the series representation

(3.22) \begin{align} P^{(2)}_2=-\sum _{n=0}^{{N}} \frac {\varphi _n} {\alpha ^{2n-1}} \frac {\partial ^{2n+1}}{\partial z^{2n+1}} P^\infty . \end{align}

Here, $N$ is the maximum number of terms in the series required for the best approximation. This upper limit is necessary because after the initial convergence, the general term diverges for large $n$ . Equation (3.22) can be expressed in terms of Legendre polynomials as

(3.23) \begin{align} P_2^{(2)} = 2\sum _{n=0}^N (2n+2)! \, \frac { \varphi _n}{\alpha ^{2n-1}} \, \frac {P_{2n+2} (z/r)}{r^{2n+3}} . \end{align}

For $N=0$ , the leading-order term is given by

(3.24) \begin{align} P_2^{(2)} \simeq \frac {2 \alpha }{r^3} \left ( 3 \, \frac {z^2}{r^2}-1\right ) \end{align}

and has the properties of an axisymmetric source quadrupole. The final asymptotic expression for the far-field pressure is

(3.25) \begin{align} P=h^2\big(P^{(2)}_1+P^{(2)}_2\big) \simeq h^2 \left ( 2\alpha \, \frac {3z^2-r^2}{r^5}+ 6 \, \frac {5z^3-3r^2z}{r^7} \right ) \! . \end{align}

In the near field, $\alpha r \ll 1$ , we evaluate (3.19) by expanding $\sqrt {t^2+1}$ for large $t$ as

(3.26) \begin{align} \sqrt {t^2+1} = t + \frac {1}{2t} - \frac {1}{8t^3} + \mathcal{O} \left ( \frac {1}{t^5} \right )\!. \end{align}

We obtain

(3.27) \begin{align} P_2^{(2)} \simeq 6 \, \frac {5z^3-3r^2z}{r^7} +\alpha ^2\, \frac {z}{r^3} + \frac {\alpha ^4}{4} \log \left ( \alpha (r+z)\right ) \quad (\alpha r \ll 1). \end{align}

The first term is identical to $P_1^{(2)}$ from (3.17). Together, they give the pressure $P=12 h^2 (5z^3-3r^2z)r^{-4}$ , in agreement with the pressure from Blake’s solution (Blake Reference Blake1971) in the limit $h\to 0$ . The second term represents the leading-order correction in a Brinkman fluid. Note that the improper integral in (3.19) diverges in the third term, but it can still be evaluated up to a constant by first computing $\boldsymbol {\nabla }\!P_2^{(2)}$ and then integrating it in space.

A comparison between the numerically evaluated $P^{(2)}_2$ and both asymptotic approximations is shown in figure 1. The far-field prediction from (3.24), shown in figure 1(a), agrees well with the numerical integration for $\alpha r \gg 1$ . Overall, the best agreement is achieved with $N=1$ . Adding additional terms improves the accuracy at $\alpha r \gg 1$ , but worsens the divergence in the intermediate range ( $\alpha r \lesssim 2$ ). The near-field behaviour predicted by (3.27), illustrated in figure 1(b), gives good agreement for $\alpha r \lesssim 1$ . Because $P^{(2)}_2$ is just one of the contributions, and the others will be given by closed-form analytical expressions, the overall error of any of the approximations in the final result will be less pronounced.

Figure 1. Comparison of the numerically evaluated $P^{(2)}_2$ (symbols) with the asymptotic approximations (solid lines) from (3.24) in the far-field limit $\alpha r \gg 1$ and (3.27) (up to order $\mathcal O(\alpha ^2)$ ) in the near-field limit $\alpha r \ll 1$ .

3.3. Shear-driven flow

The first contribution to $\varPsi ^{(2)}$ in (3.15) is obtained as the second derivative of $\varPsi ^\infty$ from (3.3), and reads

(3.28) \begin{align} \varPsi ^{(2)}_1= 2\left ( -(1+\alpha r) + \frac {z^2}{r^2}\big(3 + 3\alpha r +(\alpha r)^2 \big) \right ) \frac {\textrm{e}^{-\alpha r}}{r^3}. \end{align}

The second contribution, expressed with a dimensionless integrand, is

(3.29) \begin{align} \varPsi ^{(2)}_2 = -2\alpha ^3 \int _0^\infty t^2 \, \textrm{e}^{-\alpha z \sqrt {t^2+1}} \, J_0(\alpha \rho t) \, \mathrm{d} t , \end{align}

which can be evaluated exactly. To do so, we introduce the following convergent improper integral

(3.30) \begin{align} S_0 = \int _0^\infty \frac {1}{\sqrt {t^2+1}} \, \textrm{e}^{-a \sqrt {t^2+1}}\, J_0(b t) \, \mathrm{d} t , \end{align}

which can be found in Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2014), equation (6.645.1), by making the substitution $x = \sqrt {t^2 + 1}$ . Defining $\lambda _\pm = ( \sqrt {a^2+b^2} \pm a ) / 2$ , we obtain

(3.31) \begin{align} S_0 = I_0\kern-1pt\left ( \lambda _- \right ) K_0\kern-1pt\left ( \lambda _+ \right ) \!, \end{align}

with $K_0$ and $I_0$ denoting the zeroth-order modified Bessel functions of the second and first kinds, respectively. Using higher-order derivatives of $S_0$ with respect to $a$ , denoted by $S_0^{(m)} = \partial ^m S_0 / \partial a^m$ , the solution of the integral in (3.29) is given by

(3.32) \begin{align} \varPsi ^{(2)}_2 = 2\alpha ^3 \big(S_0^{(3)} - S_0^{(1)}\big) , \end{align}

evaluated at $a=\alpha z$ , $b=\alpha \rho$ and $\lambda _\pm =\alpha (r\pm z)/2$ . By simplifying the resulting expressions, we obtain

(3.33) \begin{align} \varPsi _2^{(2)} = I_0(\lambda _-) \left [ A_0 K_0(\lambda _+) + \varGamma _+ K_1(\lambda _+) \right ] +I_1(\lambda _-) \left [ A_1 K_1(\lambda _+) + \varGamma _- K_0(\lambda _+) \right ] \! , \end{align}

where

(3.34) \begin{align} A_0 = \frac {\alpha ^2 z}{r^2} \left ( 1-\frac {3z^2}{r^2} \right ) \! , \quad A_1 = - \frac {3\alpha ^2 z\rho ^2}{r^4} , \end{align}

and

(3.35) \begin{align} \varGamma _\pm = \frac {\alpha }{ r^2} \left ( 1 \pm \big (1+\alpha ^2\rho ^2\big ) \frac {z}{r} -\frac {3z^2}{r^2} \mp \frac {3z^3}{r^3} \right ) \! . \end{align}

From the asymptotic form of the modified Bessel functions ( $K_0(x)\sim \exp (-x)$ and $I_0(x) \sim \exp (x)$ ), it is clear that $\varPsi ^{(2)}_2$ decays $\sim\!\exp (-\alpha z)$ , and the same then holds for the velocity field. Therefore, $\varPsi ^{(2)}_2$ describes the short-range flows induced by the tractions on the no-slip boundary.

3.4. Velocity field

With the knowledge of the two scalar fields, the velocity Green’s function is determined as

(3.36) \begin{align} \boldsymbol{G}^{(2)}= \partial _z^2\boldsymbol{G}^\infty +\alpha ^{-2} \boldsymbol{\nabla }\big({-}P^{(2)}_2 + \partial _z \varPsi ^{(2)}_2 \big) -\hat {\boldsymbol{e}}_z \varPsi ^{(2)}_2 . \end{align}

Due to the length of the resulting expression, we do not present the full result here. It can, however, be readily obtained using computer algebra systems or symbolic computation tools.

3.4.1. Far field

We now assume $\alpha r \gg 1$ , retaining only the terms with $\textrm{e}^{-\alpha z}$ , but not those with $\textrm{e}^{-\alpha r}$ . The field $\varPsi _1^{(2)}$ does not contribute to the far field, and $\varPsi _2^{(2)}$ has the asymptotic expansion

(3.37) \begin{align} \varPsi _2^{(2)} \simeq \frac {\textrm{e}^{-\alpha z}}{r^3} \left ( 2 + \frac {3\alpha z (3 + \alpha z)}{(\alpha r)^2} + \frac {15 \alpha z \big(15 + 15\alpha z + 6 (\alpha z)^2 + (\alpha z)^3 \big)}{4(\alpha r)^4} \right ) \! . \end{align}

The Green’s function for the velocity can now be split up as

(3.38) \begin{align} \boldsymbol{G}^{(2)} = \boldsymbol{W} + \boldsymbol{H} \, , \end{align}

where the first term $\boldsymbol{W}$ represents contributions from $\boldsymbol{G}^\infty$ and $P_2^{(2)}$ that correspond to a series of source multipoles. Defining

(3.39) \begin{align} \boldsymbol{g} = - \boldsymbol{\nabla } \frac 1 r = \frac {\rho \hat {\boldsymbol{e}}_\rho + z \hat {\boldsymbol{e}}_z}{(\rho ^2+z^2)^{3/2}} \end{align}

as the flow of a source singularity, this can be written compactly as

(3.40) \begin{align} \boldsymbol{W} = \frac {2}{\alpha ^2} \frac {\partial ^3 \boldsymbol{g}}{\partial z^3} -2 \sum _{n=0}^{N} \frac { \varphi _n}{\alpha ^{2n+1}} \, \frac {\partial ^{2n+2} \boldsymbol{g}}{\partial z^{2n+2}}, \end{align}

with $\varphi _n$ defined in (3.21). In terms of Legendre polynomials, the radial component of $\boldsymbol{W}$ can be expressed as

(3.41a) \begin{align} W_\rho = -\frac {12}{\alpha ^2 r^5} P_4^1(z/r) -2 \sum _{n=0}^N (2n+2)! \, \frac {\varphi _n}{\alpha ^{2n+1}} \frac {P_{2n+3}^1 (z/r)}{r^{2n+4}} , \end{align}

with $P_n^1(z/r) = -(\rho /r)\, P_n'(z/r)$ denoting the associated Legendre function of degree $n$ and order 1, where $P_n'(x)={\mathrm{d}} P_n(x)/{\mathrm{d}} x$ . The axial component is given by

(3.41b) \begin{align} W_z = \frac {48}{\alpha ^2 r^5 } \, P_4(z/r) + 2 \sum _{n=0}^N (2n+3)! \, \frac {\varphi _n}{\alpha ^{2n+1}} \frac {P_{2n+3}(z/r)}{r^{2n+4}} . \end{align}

By retaining only the first three terms of the series ( $N=2$ ), we obtain

(3.42a) \begin{align} W_\rho &\simeq \frac {3\rho }{\alpha r^5} \left ( \frac {10z^2}{r^2}-2 +\frac {10 z}{\alpha r^2} \left ( \frac {7z^2}{r^2}-3 \right ) +\frac {15}{(\alpha r)^2} \left ( \frac {21z^4}{r^4}-\frac {14z^2}{r^2}+1 \right ) \right )\!, \\[-12pt] \nonumber \end{align}

(3.42b) \begin{align} W_z &\simeq \frac {3}{\alpha r^5} \left ( \frac {10z^3}{r^2}-6z + \frac {2}{\alpha } \left ( \frac {35z^4}{r^4}-\frac {30z^2}{r^2} + 3 \right ) + \frac {5z}{(\alpha r)^2} \left ( \frac {63z^4}{r^4}-\frac {70z^2}{r^2}+15 \right ) \right )\!. \\[4pt] \nonumber \end{align}

The second term $\boldsymbol{H}$ in (3.38) contains the far-field contributions stemming from $\varPsi _2^{(1)}$ , which can be evaluated exactly, with radial and axial components given in the far field by

(3.43a) \begin{align} H_\rho &\simeq \frac {\textrm{e}^{-\alpha z}\rho }{\alpha r^5} \left ( 6+\frac {15}{(\alpha r)^2} \big ( (\alpha z)^2 + 3 \alpha z -3\big ) \right )\! , \\[-12pt] \nonumber \end{align}

(3.43b) \begin{align} H_z &\simeq -\frac {\textrm{e}^{-\alpha z}}{\alpha ^2r^5} \left ( 18 + \frac {45 \alpha z}{(\alpha r)^2} (\alpha z+5) \right )\! . \\[9pt] \nonumber \end{align}

Finally, the far-field expression for the velocity is $\boldsymbol{G}=h^2(\boldsymbol{W} +\boldsymbol{H})$ , with $\boldsymbol{W}$ given by (3.42), and $\boldsymbol{H}$ given by (3.43).

A comparison between the far-field flow and the result of numerical integration of (3.9) and their derivatives is shown in figure 2, showing near-perfect agreement in the far-field regime.

Figure 2. Streamline plot of the axisymmetric flow induced by a Brinkmanlet near a no-slip wall, obtained using (a) the far-field solution from (3.42) and (3.43), and (b) numerical integration for $\alpha h = 0.1$ .

3.4.2. Near field

The velocity can also be analytically evaluated in the near-field limit, where $\alpha r \ll 1$ , while still maintaining the condition $r \gg h$ . We obtain

(3.44) \begin{align} \varPsi _2^{(2)} \simeq \frac {2r^2-6z^2}{r^5} + \alpha ^2 \frac {z^2}{r^3} +\frac {\alpha ^4}{4} \left ( z\left ( \frac {1}{4} + \gamma + \log \left ( \frac {\alpha (r+z)}{4} \right ) \right ) - \frac {z^2}{r} \right ) \! . \end{align}

Substituting the expressions for $\varPsi _2^{(2)}$ from (3.44) and $P_2^{(2)}$ from (3.27), together with $\boldsymbol{G}^\infty$ from (3.4), into (3.36), we obtain

(3.45a) \begin{align} G_\rho ^{(2)} &\simeq \frac {6\rho z}{r^7}( 3z^2-2\rho ^2) + \frac {\alpha ^2 \rho z }{2r^5} ( 2\rho ^2-z^2), \\[-12pt] \nonumber \end{align}

(3.45b) \begin{align} G_z^{(2)} &\simeq \frac {6z^2}{r^7} ( 2z^2-3\rho ^2) + \frac {\alpha ^2 z^2}{2r^5} ( \rho ^2-2z^2). \\[9pt] \nonumber \end{align}

The terms of order $\mathcal{O}(\alpha ^0)$ coincide with Blake’s tensor for a normal force in the limit $r\gg h$ (Blake & Chwang Reference Blake and Chwang1974). Terms of order $\alpha ^2$ then represent the leading correction in the Brinkman flow.

3.5. Discussion

A complete solution of the flow problem in the presence of a no-slip wall will also require the Green’s function for a force acting parallel to the wall, in which case the height dependence will be $\sim\!h$ . Although the absence of axial symmetry makes the problem more complex, we hope that the solution method developed here can be generalised to parallel forces as well.

Another open question is whether the same problem can be solved in a two-dimensional Brinkman fluid, which can also be used as an approximation for the Hele-Shaw flow in a thin channel.

Out results can be used to model hydrodynamic interactions in porous environments where particles or microorganisms move near boundaries. A force-free swimmer can, for example, be modelled by placing two oppositely directed point forces asymmetrically relative to the swimmer body (Menzel et al. Reference Menzel, Saha, Hoell and Löwen2016; Daddi-Moussa-Ider & Menzel Reference Daddi-Moussa-Ider and Menzel2018; Hosaka et al. Reference Hosaka, Golestanian and Daddi-Moussa-Ider2023). In biological contexts near boundaries, such as microorganism motion in gels or tissues, these solutions describe long-range hydrodynamic disturbances. Moreover, they serve as fundamental building blocks for modelling suspensions near walls or active matter in confined or partially permeable environments.

Although we introduced the Brinkman equation for porous media, our solution remains valid for unsteady flows in the presence of inertia, where we use $\alpha ^2= - {\rm i} \omega /\nu$ , with the angular frequency $\omega$ , and the kinematic viscosity $\nu =\eta /\rho$ (Kim & Karrila Reference Kim and Karrila2013). For instance, the solution can then be applied for any oscillatory assembly of forces in the proximity of a no-slip boundary. For example, it has been shown that for sufficiently far distances, inertial effects become relevant in the flows driven by a beating cilium (Wei et al. Reference Wei, Dehnavi, Aubin-Tam and Tam2019, Reference Wei, Dehnavi, Aubin-Tam and Tam2021). A straightforward extension of our solution could also apply to source singularities near a wall, such as an oscillating gas bubble in the linear regime, or more complex multipole singularities.