3.1. Net force on a fixed sphere

We consider the special case of a fixed particle, $\boldsymbol {U} =\boldsymbol {0}$, which simplifies (3.6) to

(3.9)\begin{equation} \boldsymbol{U}=\boldsymbol{0}:\quad \boldsymbol{F}=\boldsymbol{\mathcal{M}}_{f} \boldsymbol{\cdot} \boldsymbol{f}_{\!\!f}^{{\circ}} + \boldsymbol{\mathcal{M}}_{n} \boldsymbol{\cdot} \boldsymbol{f}_{\!\!n}^{{\circ}}. \end{equation}

It is useful to decompose the point forces and the force on the particle in parallel and perpendicular directions of the separation vector. Using $\boldsymbol {f}_{\{\,f,n\}}^{\circ } = {f}_{\{\,f,n\}}^{\parallel } \hat {\boldsymbol {R}} + {f}_{\{\,f,n\}}^{\perp } \hat {\boldsymbol {R}}^{\perp }$ and $\boldsymbol {F}=F^\parallel \hat {\boldsymbol {R}}+F^\perp \hat {\boldsymbol {R}}^{\perp }$, we get

(3.10a,b)\begin{equation} {F}^{{\parallel}} = A_f (R) {f}_f^{{\parallel}} + A_n (R) {f}_n^{{\parallel}} ,\quad {F}^{{\perp}} = B_f (R) {f}_f^{{\perp}} + B_n (R) {f}_n^{{\perp}} . \end{equation}

We begin by analysing the limiting behaviours of $A_{\{\,f,n\}}$ and $B_{\{\,f,n\}}$ at early ($t\to 0$) and long ($t\to \infty$) times. These limits can be computed analytically using the following relationships: $f(t\to 0^+)= \lim _{s\to \infty } s f(s)$ and $f(t\to \infty )= \lim _{s\to 0} s f(s)$ and the results are

(3.11a) \begin{align} {A}_{ \{\,f , n\} } (t\to 0)&= \frac{3}{2} \left( \frac{a}{R} \right) -\frac{1}{2} {\left( \frac{a}{R} \right)}^3 , \quad {B}_{ \{\,f , n\} } (t\to 0)= \frac{3}{4} \left( \frac{a}{R} \right) + \frac{1}{4} {\left( \frac{a}{R} \right)}^3 , \end{align}

(3.11b) \begin{align} A_f (t\to \infty) &= \frac{6(1-\nu)}{5-6\nu} \left( \frac{a}{R} \right) - \frac{1}{5-6\nu} {\left( \frac{a}{R} \right)}^3 \nonumber\\ & \quad + \frac{3a(1-2\nu)}{{\beta}_{{\circ}}^2 R^3 (5-6\nu)} [ 1+a {\beta}_{{\circ}} -(1+ {\beta}_{{\circ}} R )\times {\mathrm{e}}^{-{\beta}_{{\circ}} (R-a)} ] , \end{align}

(3.11c) \begin{align} B_f(t\to \infty) &= \frac{3(1-\nu)}{5-6\nu} \left( \frac{a}{R} \right) + \frac{1}{2(5-6\nu)} {\left( \frac{a}{R} \right)}^3 \nonumber\\ & \quad - \frac{3a(1-2\nu)}{2 {\beta}_{{\circ}}^2 R^3 (5-6\nu)} [ 1+a {\beta}_{{\circ}} -(1+ {\beta}_{{\circ}} R+ {\beta}_{{\circ}}^2 R^2 ) \times {\mathrm{e}}^{-{\beta}_{{\circ}} (R-a)} ] , \end{align}

(3.11d)\begin{align} A_n (t\to \infty) &= \frac{6(1-\nu)}{5-6\nu} \left( \frac{a}{R} \right) - \frac{1}{5-6\nu} {\left( \frac{a}{R} \right)}^3 , \end{align}

(3.11e)\begin{align} B_n (t\to \infty) &= \frac{3(3-4\nu)}{2(5-6\nu)} \left( \frac{a}{R} \right) + \frac{1}{2(5-6\nu)} {\left( \frac{a}{R} \right)}^3. \end{align}

As shown in (3.11a), at short times the response functions for both fluid and network asymptote to their well-known respective forms in Stokes flow (Kim & Karrila Reference Kim and Karrila2013). This can be explained by noting that at early times the network co-moves with the fluid phase, $\boldsymbol {v}_n(\boldsymbol {r},t\to 0)=\boldsymbol {v}_f(\boldsymbol {r},t\to 0)$ and $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {v}_{\{\,f,n\}}(\boldsymbol {r},t\to 0)=0$. Furthermore, at long times the response functions of the network phase ($A_n$ and $B_n$), approach their well-known form in linear elasticity (Phan-Thien & Kim Reference Phan-Thien and Kim1994); see (3.11d) and (3.11e). The long time forms of the fluid response functions ($A_f$ and $B_f$) include extra $\beta _\circ$ dependent terms, which, as expected, become identically zero for an incompressible network.

Variations of the net force with time and the distance from the point forces in the parallel direction are determined by the behaviour of $A_{\{\,f,n\}}$. Since we have the analytical form of these functions in short and long times in (3.11), we study the time-dependent behaviour of the following quantities:

(3.12a,b)\begin{align} \hat{A}_{\{\,f,n\}}= \frac{A_{\{\,f,n\}}(t,R)-A_{\{\,f,n\}}(t\to \infty,R)}{A_{\{\,f,n\}}(0,R)-A_{\{\,f,n\}}(t\to \infty,R)},\quad \hat{B}_{\{\,f,n\}}=\frac{B_{\{\,f,n\}}(t,R)-B_{\{\,f,n\}}(t\to \infty,R)}{B_{\{\,f,n\}}(0,R)-B_{\{\,f,n\}}(t\to \infty,R)}.\end{align}

Here $\hat {A}$ and $\hat {B}$, by construction relax over time from $\hat {A}=\hat {B}=1$ at $t=0$ to $\hat {A}=\hat {B}\to 0$ as $t\to \infty$. Next, we explore how the relaxation dynamics changes with the separation distance, $R$, and material properties of the fluid and network, $\eta, G, \beta _\circ$ and $\nu$.

The expressions for $A_{\{\,f,n\}}$ and $B_{\{\,f,n\}}$ in (3.2) contain two types of $s$ and $R$ dependencies. One obvious relaxation mechanism is shear relaxation, $\tau _\circ =\eta /G$, which is independent of $\beta _\circ$ and $\nu$ and the separation distance, $R$. This is the relaxation time measured in the standard particle tracking microrheology. The terms involving $S_1$, $S_2$ and $S_3$ can be written as a product of only $s$-dependent and only $r$-dependent functions, $\tilde {\mathcal {Q}}(s) \mathcal {P}(r)$, resulting in a $\mathcal {Q}(t) \mathcal {P}(r)$ form after Laplace inversion. As a result, the relaxation behaviour becomes independent of $R$. In comparison, the terms that include $\exp (-\beta (s) r)$ and $\exp (-\alpha (s) r)$ cannot be written as a product of $s$- and $r$-dependent functions, and so the relaxation time of the force on the sphere becomes a function of the distance from the point forces.

As shown in our earlier work (Moradi et al. Reference Moradi, Shi and Nazockdast2022) and others (Detournay & Cheng Reference Detournay and Cheng1993; Doi Reference Doi2009), the divergence of the displacement field in a linear elastic network, $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {v}_n$, is described by a diffusion equation, where the diffusion coefficient is $D=\tau _\circ \alpha _\circ ^{-2}$. Based on this, we can introduce two extra time scales: $\tau _D^a=(\alpha _\circ a)^2\tau _\circ$, which is the time scale for the network compressibility to diffuse over a particle radius, $a$; and $\tau _D^{R}=(\alpha _\circ R)^2\tau _\circ$ is the time for diffusing a distance $R$ from the point forces.

If the $S_1$, $S_2$ and $S_3$ terms dominate the response, we expect the force relaxation time to be independent of $R$ and determined by $\tau _D^ a$; alternatively, if $\exp (-\alpha (s) r)$ and $\exp (-\beta (s) r)$ terms determine the behaviour of $A_{\{\,f,n\}}$ and $B_{\{\,f,n\}}$, we expect the force relaxation to be distance dependent and in part determined by $\tau _D^{R}$.

We begin with presenting the results for the relaxation dynamics of fluid phase functions $\hat {A}_f$ and $\hat {B}_f$ versus time in figure 3 for different values of Poisson ratio, $\nu$ (figure 3a), permeability, $\beta _\circ$ (figure 3b) and distance from the point force $R$ (figure 3c). The solid lines that go through the markers are the computed values of $\widehat {(S_1/\varDelta )}$, which is simply the values of $S_1/\varDelta$ made dimensionless in the exact manner as $A_f$ and $B_f$ in (3.2). The relaxation time is made dimensionless with $\tau _D^a$ in all figures. When $R$ is fixed and $\nu$ and $\beta _\circ$ are varied (figure 3a,b), $\tau _D^a$ and $\tau _D^R$ are linearly related by a constant factor, $\tau _D^R/\tau _D^a = (R/a)^2$, and using $\tau _D^R$ (instead of $\tau _D^a$) would only shift the plots in the $x$ axis and would not change the overall superposition of the curves. In comparison, when $R$ is varied, the two time scales become distinct. Thus, to test if the relaxation is distance dependent or not, we present the results of figure 3(c) as a function of $t/\tau _D^R$ in the insets of those figures.

Figure 3.

Time variations of the fluid hydrodynamic functions, $\hat {A}_f$ (top row) and $\hat {B}_f$ (bottom row). (a) Different values of Poisson ratio, $\nu$, while setting $\beta _\circ =5$ and $R=2$. (b) Different values of permeability, $\beta _\circ$, while setting $\nu =0.3$ and $R=2$. (c) Different separation distances, $R$, while setting $\nu =0.3$ and $\beta _\circ =5$. The time axis is made dimensionless with $\tau _D^a=a^2\alpha _\circ ^2\tau _\circ$ in all the main plots. The values of $R$ and $\beta _\circ$ are made dimensionless with the sphere radius, which is equivalent to taking $a:=1$. The inset figures in figure 3(c) show the same results as the main figure, but the time axes are made dimensionless with $\tau _D^R=R^2\alpha _\circ ^2\tau _\circ$.

As can be seen, in all cases $\widehat {(S_1/\varDelta )}$ (solid lines) closely match the predictions that include all the terms (symbols). This suggests that the relaxation is distance independent and that $\tau _D^a$ (and not $\tau _D^R$) is a more appropriate choice for relaxation time of $\hat {A}_f$ and $\hat {B}_f$. This can be clearly seen by comparing the nearly perfect superposition of the plots in the main figure against the inset plots in figure 3(c). Furthermore, the relaxation curves for different values of $\nu$ superimpose, when time is made dimensionless with $\tau _D^a$. The superposition is not as clean at shorter times when permeability is varied (see figure 3b), since shear relaxation time ($\tau _\circ$) becomes increasingly more important at shorter time scales, where the fluid and network phase co-move. However, at longer times ($t/\tau _D^a\ge 10$) the plots superimpose, showing that the long time behaviour of these hydrodynamic functions is determined by $\tau _D^a$. Collectively, these results show that the relaxation dynamics of $\hat {A}_{f}$ and $\hat {B}_f$ is not distance dependent, and is dominated by the terms involving $(S_1/\varDelta )(a/r)$ and the relaxation time, $\tau _D^a$.

Figure 4 shows the relaxation dynamics of the network functions, $\hat {A}_n$ and $\hat {B}_n$. In this case, the solid lines are just showing the detailed predictions and are not showing $\widehat {(S_1/\varDelta )}$, since $\widehat {(S_1/\varDelta )}$ values did not match the detailed predictions. The time axis is made dimensionless by $\tau _D^{R}$ in all the main plots. Again, when the separation distance is varied, we replot the results of figure 4(c) in the inset figure, but with the dimensionless time being defined as $t/\tau _D^a$. As can be seen in figure 4(c), using $t/\tau _D^R$ gives a better collapse of the relaxation plots at long times and at different $R$ in the main plot, compared with the inset figures, with $t/\tau _D^a$ as the $x$ axis. However, the collapse is not as clean as the results for $\hat {A}_f$ and $\hat {B}_f$. For both cases of varying permeability and distance (figures 4b and 4c), the plots fail to collapse at shorter times, especially for $\hat {B}_n$, where the shape of the plots qualitatively changes with varying $\beta _\circ$ and $R$.

Figure 4.

Time variations of the network hydrodynamic functions, $\hat {A}_n$ (top row) and $\hat {B}_n$ (bottom row). (a) Different values of Poisson ratio, $\nu$, while setting $\beta _\circ =5$ and $R=2$. (b) Different values of permeability, $\beta _\circ$, while setting $\nu =0.3$ and $R=2$. (c) Different separation distances, $R$, while setting $\nu =0.3$ and $\beta _\circ =5$. The time axis is made dimensionless with $\tau _D^R=R^2\alpha _\circ ^2\tau _\circ$ in all the main plots. The inset figures in figure 4(c) show the same results as the main figure, but the time axes are made dimensionless with $\tau _D^a=a^2\alpha _\circ ^2\tau _\circ$.

All together these results suggest that the relaxation behaviour of the network functions $\hat {A}_n$ and $\hat {B}_n$ at long times is distance dependent and best described by $\tau _D^{R}$. However, at shorter times the fluid shear stress plays a more central role and the relaxation dynamics is a function of all three time scales ($\tau _\circ$, $\tau _D^a$ and $\tau _D^R$). Including the shear stresses in the conservation equations allows us to study the complex variations in the relaxation behaviour from short to long time scales.

4.1. Slip BC for the network phase

In deriving the expressions for the sphere's response to point forces, we assumed a no-slip BC for the network phase at the sphere's surface. This assumption generally holds when the size of the sphere is considerably larger than the network and there is no active interactions between the network and the sphere. As an example, the mesh size for actin networks in the cell is roughly in the range 20–200 nm (Charras et al. Reference Charras, Yarrow, Horton, Mahadevan and Mitchison2005; Keren et al. Reference Keren, Yam, Kinkhabwala, Mogilner and Theriot2009), and so the no-slip BC is appropriate for micron-sized endogenous bodies and injected probes used in microrheological measurements. But in many cases, such as transport of proteins and nanoscopic particles, the size of the inclusion can become comparable or smaller than the mesh size. In such conditions, we must account for slip velocity between the inclusion and the network phase (a no-slip BC generally holds for the fluid phase down to atomistic lengths).

One method of modelling the slip velocity is to use the equivalent of the Navier equation, where the tangential slip velocity is proportional to the traction in that direction, i.e.

(4.1)\begin{equation} \boldsymbol{v}_n =\boldsymbol{U}+\frac{b}{\mu}(\boldsymbol{I}-\boldsymbol{n} \boldsymbol{n}) \boldsymbol{\cdot} (\boldsymbol{\varSigma}_n \boldsymbol{\cdot} \boldsymbol{n}), \end{equation}

where $b$ is the slip length and $\mu$ is the term with dimensions of viscosity. Recall that the above equation is written in Laplace space. The time dependence of $\mu$ depends on the mechanical coupling of the network with the sphere's surface in microscopic scale. For example, if the network is not tethered to the boundary, we may assume that the slip dynamics is viscous, and take $\mu (t)=\mu _\circ \to \mu (s)=\mu _\circ /s$; in the opposite case of the network tethered to the boundary, the slip dynamics may be predominantly elastic, leading to slip displacements rather than slip velocities: $\mu (t)=\mu _\circ \delta (t) \to \mu (s)=\mu _\circ$. In any case, we recover a no-slip BC when $b=0$.

Our goal is to find the leading-order correction to the force on a sphere moving with velocity $\boldsymbol {U}$, subject to this BC. We would like to avoid solving this boundary value problem by using the reciprocal theorem. We use the sphere moving with velocity $\hat {\boldsymbol {U}}$, with a no-slip BC for both phases as our auxiliary problem, and use the slip BC for the network as our target solution. Upon substitution we find that

(4.2)\begin{equation} \boldsymbol{F} \boldsymbol{\cdot} {\hat{\boldsymbol{U}}}= {\hat{\boldsymbol{F}}} \boldsymbol{\cdot} \boldsymbol{U}+\left(\frac{b}{\mu}\right) \int_{S_p} ( {\hat{\boldsymbol{\varSigma}}}_n \boldsymbol{\cdot} \boldsymbol{n})\boldsymbol{\cdot} [(\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot} (\boldsymbol{\varSigma}_n\boldsymbol{\cdot} \boldsymbol{n})] \,\mathrm{d}S.\end{equation}

Next, we assume a small slip length (${b}/{a}\ll 1$), and so we can write down the following perturbation expansions:

(4.3a,b)\begin{equation} \boldsymbol{v}_n = {\hat{\boldsymbol{v}}}_n +\left(\frac{b}{a}\right)\boldsymbol{v}_{n}^{(1)} +O\left[\left(\frac{b}{a}\right)^2\right],\quad \boldsymbol{\varSigma}_n ={\hat{\boldsymbol{\varSigma}}}_n + \left(\frac{b}{a}\right) \boldsymbol{\varSigma}_n^{(1)} +O\left[\left(\frac{b}{a}\right)^2\right]. \end{equation}

Substituting these expansions into (4.2), replaces $\boldsymbol {\varSigma }_n$ with ${\hat {\boldsymbol {\varSigma }}}_n$, and so we get

(4.4)\begin{equation} \boldsymbol{F} \boldsymbol{\cdot} {\hat{\boldsymbol{U}}} = {\hat{\boldsymbol{F}}} \boldsymbol{\cdot} \boldsymbol{U}+\left(\frac{b}{\mu}\right) \int_{S_p} [ ({\hat{\boldsymbol{\tau}}}_n\boldsymbol{\cdot} {\hat{\boldsymbol{\tau}}}_n )-({\hat{\boldsymbol{\tau}}}_n \boldsymbol{\cdot} \boldsymbol{n})^2] \,\mathrm{d}S +O\left[\left(\frac{b}{a}\right)^2\right],\end{equation}

where ${\hat {\boldsymbol {\tau }}}_n={\hat {\boldsymbol {\varSigma }}}_n \boldsymbol {\cdot } \boldsymbol {n}$ is the network traction of the auxiliary problem. We see that the reciprocal theorem can be used to compute the first-order correction to the drag due to slip in the network phase. Almost the exact same derivation can be used to model the slip BC in Stokes flow, as outlined in Masoud & Stone (Reference Masoud and Stone2019).

When the inclusion is considerably smaller than the mesh ($b/a\gg 1$), it can move through the network pores and channels without experiencing any force from the network. In this limit, no network traction, $\boldsymbol {\varSigma }_n \boldsymbol {\cdot } \boldsymbol {n}=0$ at $|\boldsymbol {r}|=a$ would be a more appropriate BC. In this scenario and under a constant external force, the sphere will move with a steady-state velocity (viscous response) with the fluid velocity and the sphere's drag determined by the Brinkman equation. In comparison, when a no-slip BC is assumed for the network, as we showed earlier, the long time behaviour of the sphere under a constant force will be elastic.

Let us now consider the auxiliary problem of a sphere moving with velocity $\hat {\boldsymbol {U}}$, but this time with a no-slip BC only applied to the fluid domain and zero traction BCs applied for the network phase. It is straightforward to find the displacement fields for this problem, using the general solutions of PVE equations (Moradi et al. Reference Moradi, Shi and Nazockdast2022). Fu, Shenoy & Powers (Reference Fu, Shenoy and Powers2008) also provide a closed-form solution of this problem. Assuming that we have the solution for this auxiliary problem, we would like to find the net force on the sphere moving with velocity $\boldsymbol {U}$ and subjected to point forces in the fluid and network phases, using no traction from the network as a BC. After substituting the auxiliary solution to the reciprocal theorem, and applying the BC, we recover (3.6) and (3.7a,b), with the distinction that $A_{\{\,f,n\}}$ and $B_{\{\,f,n\}}$ now correspond to the auxiliary problem with zero network traction as the BC.

Finally, we ask if the combination of perturbation methods and the reciprocal theorem can be used to find the sphere's response functions to point forces, subject to the slip BC given by (4.1). Again, we use the sphere moving with velocity $\hat {\boldsymbol {U}}$ with no-slip BCs for both phases as the auxiliary problem, and seek to find the net force on the sphere in response to point forces $\boldsymbol {f}_{\{\,f,n\}} =\boldsymbol {f}_{\{\,f,n\}}^{\circ } \delta (\boldsymbol {r} - \boldsymbol {R})$, with the network BC given by (4.1). After substitution and using $b/a\ll 1$ we get

(4.5)\begin{align} \boldsymbol{F} \boldsymbol{\cdot} {\hat{\boldsymbol{U}}} &= {\hat{\boldsymbol{F}}} \boldsymbol{\cdot} \boldsymbol{U}+{\hat{\boldsymbol{v}}}_f (\boldsymbol{R}) \boldsymbol{\cdot} \boldsymbol{f}_{\!\!f}^{{\circ}}+{\hat{\boldsymbol{v}}}_n(\boldsymbol{R})\boldsymbol{\cdot} \boldsymbol{f}_{\!\!n}^{{\circ}} \nonumber\\ &\quad +\left(\frac{b}{\mu}\right) \int_{S_p} ( {\hat{\boldsymbol{\varSigma}}}_n \boldsymbol{\cdot} \boldsymbol{n})\boldsymbol{\cdot} [(\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot} (\boldsymbol{\varSigma}_n^{(0)} \boldsymbol{\cdot} \boldsymbol{n})] \,\mathrm{d}S +O\left[\left(\frac{b}{a}\right)^2\right]. \end{align}

The first line of (4.5) is exactly the expression we get using the no-slip BC in both phases. The second line shows the correction to the net force on the sphere due to network slip. Here, $\boldsymbol {\varSigma }_n^{(0)}$ is the network stress in the target problem with a no-slip BC in both phases. Thus, we see that using the reciprocal theorem comes with the caveat of having to solve the problem of a sphere subjected to point forces, with a no-slip BC in both phases. Of course, the solution can be obtained, using the general solutions and the fundamental solutions to point forces in the fluid and network phase, but the calculations are lengthy and will not be pursued here.

4.2. Nonlinearities in the fluid and network forces and stresses

Another utility of the reciprocal theorem in fluid mechanics is in calculating the leading-order effects of nonlinear viscoelastic stresses (Leal Reference Leal1980; Hu & Joseph Reference Hu and Joseph1999; Koch & Subramanian Reference Koch and Subramanian2006; Khair & Squires Reference Khair and Squires2010; Elfring Reference Elfring2015) and finite inertia (Ho & Leal Reference Ho and Leal1974) in the net force and velocity of a particle, without the need to directly solve the nonlinear problem. See Masoud & Stone (Reference Masoud and Stone2019) for a comprehensive review of the topic.

The same ideas can be directly applied to study the nonlinear effects in biphasic systems. To demonstrate this point, we reconsider the problem of a sphere moving a prescribed velocity $\boldsymbol {U}$ in an isotropic elastic network and a Newtonian fluid (with a no-slip BC in both phases), and derive an expression that allows computing the leading-order correction to the sphere's drag due to nonlinear terms that appear in the equations in finite deformations. One possible source of nonlinearity in the problem is the contribution of hyperelastic stresses in finite deformations of the elastic network. The simplest hyperelastic model is the Saint Venant–Kirchhoff model (Nezamabadi, Zahrouni & Yvonnet Reference Nezamabadi, Zahrouni and Yvonnet2011), where a linear stress–strain relationship is still assumed and only geometric nonlinearities are accounted for. We choose this model for its simplicity, but the methodology can be directly used for other choices of hyperelastic models that reduce to linear elasticity in small deformations. In the Saint Venant–Kirchhoff model the elastic stress in Eulerian coordinates is given by

(4.6a,b)\begin{equation} \boldsymbol{\sigma}_n = 2G \boldsymbol{e} +\lambda \,\text{trace}(\boldsymbol{e}),\quad \boldsymbol{e}=\tfrac{1}{2}[\boldsymbol{\nabla} \boldsymbol{u}_n+(\boldsymbol{\nabla} \boldsymbol{u}_n)^T] - \boldsymbol{\nabla} \boldsymbol{u}_n\boldsymbol{\cdot} (\boldsymbol{\nabla} \boldsymbol{u}_n)^T, \end{equation}

where $\boldsymbol {e}$ is the Eulerian finite strain tensor, which also includes the nonlinear term that is the product of the displacement gradient and its transpose. Another source on nonlinearity is in the term that models the relative frictional forces between the fluid and network, $\xi (\partial \boldsymbol {u}_n/\partial t -\boldsymbol {v}_f)$. Specifically, the partial time derivative of the displacement field ($\partial \boldsymbol {u}_n/\partial t$) must be replaced with the velocity of material elements (Fu et al. Reference Fu, Shenoy and Powers2010), which is defined as $\textrm {d}\boldsymbol {u}_n/\textrm {d}t =\partial \boldsymbol {u}_n/\partial t + (\textrm {d}/\textrm {d}t)(\boldsymbol {u}_n(t) \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u}_n(t))$.

We use perturbation expansion of the displacement and stress fields and the net force, i.e.

(4.7a–c)\begin{equation} \boldsymbol{v}_{\{\,f,n\}}=\boldsymbol{v}^{(0)}_{\{\,f,n\}} + \boldsymbol{v}^{(1)}_{\{\,f,n\}}, \quad \boldsymbol{\varSigma}_{\{\,f,n\}}=\boldsymbol{\varSigma}^{(0)}_{\{\,f,n\}} + \boldsymbol{\varSigma}^{(1)}_{\{\,f,n\}},\quad \boldsymbol{F}=\boldsymbol{F}^{(0)}+\boldsymbol{F}^{(1)}, \end{equation}

and assume a no-slip BC, which reduces to the following BCs for $\boldsymbol {v}^{(0)}$ and $\boldsymbol {v}^{(1)}$:

(4.8a,b)\begin{equation} \boldsymbol{v}_{\{\,f,n\}} ^{(0)} (|\boldsymbol{r}|=a)=\boldsymbol{U}, \quad \boldsymbol{v}_{\{\,f,n\}}^{(1)} (|\boldsymbol{r}|=a) =\boldsymbol{0}. \end{equation}

The solution to the linear problem is known: $\boldsymbol {F}^{(0)}=-\mathcal {R}\boldsymbol {U}$. Importantly, $\boldsymbol {v}^{(1)}_{\{\,f,n\}}$ are also described by linear biphasic equations with the difference that the equations also include body forces that are only dependent on the solution of the linear problem:

(4.9a)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\varSigma}_f^{(1)}+ \xi (\boldsymbol{v}^{(1)}_n - \boldsymbol{v}^{(1)}_f ) ={-}\boldsymbol{f}_{\!\!f},\quad \boldsymbol{f}_{\!\!f}=s\mathcal{L}[\boldsymbol{u}_n^{(0)} (t)\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}^{(0)}_n(t)], \end{gather}$$

(4.9b)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\varSigma}_n^{(1)} - \xi (\boldsymbol{v}^{(1)}_n - \boldsymbol{v}^{(1)}_f ) ={-}\boldsymbol{f}_{\!\!n},\quad \boldsymbol{f}_{\!\!n}=\mathcal{L}[\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma}^{(0)}_{n,{n.l}}(t)]- s\mathcal{L}[\boldsymbol{u}_n^{(0)}(t)\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_n^{(0)}(t)]. \end{gather}$$

Here $\boldsymbol {\sigma }^{(0)}_{n, {n.l}}=-2G[\boldsymbol {\nabla } \boldsymbol {u}_n^{(0)} \boldsymbol {\cdot } (\boldsymbol {\nabla } \boldsymbol {u}_n^{(0)})^T]-\lambda \, \textrm {trace}[\boldsymbol {\nabla } \boldsymbol {u}_n^{(0)} \boldsymbol {\cdot } (\boldsymbol {\nabla } \boldsymbol {u}_n^{(0)})^T]$ is the nonlinear piece of the elastic stress evaluated from the linear solution.

Next, we use the reciprocal theorem to find the drag contribution for the nonlinear terms, $\boldsymbol {F}^{(1)}$, without directly solving (4.9). We use the solution to the linear problem with no body force as our auxiliary problem. Upon substitution we get

(4.10)\begin{equation} {\boldsymbol{F}} \boldsymbol{\cdot} {\hat{\boldsymbol{U}}} =\int_V ({\hat{\boldsymbol{v}}}_f\boldsymbol{\cdot} \boldsymbol{f}_{\!\!f}+{\hat{\boldsymbol{v}}}_n \boldsymbol{\cdot} \boldsymbol{f}_{\!\!n} )\, \mathrm{d}V, \end{equation}

where ${\hat {\boldsymbol {U}}} := \boldsymbol {U}$, ${\boldsymbol {F}} :=\boldsymbol {F}^{(1)}$ and $\hat {\boldsymbol {v}}_{\{\,f,n\}} := \boldsymbol {v}_{\{\,f,n\}}^{(0)}$. Thus, $\boldsymbol {F}^{(1)}$ can be computed by evaluating the integral involving only known terms. Analytical evaluation of this integral is tedious, but the numerical integration is straightforward.

The same methodology can be used to account for weak nonlinear stresses in the fluid phase, due to viscoelasticity (small Deborah numbers), or shear rate-dependent viscosity (Leal Reference Leal1980; Elfring & Lauga Reference Elfring and Lauga2015; Boyko & Stone Reference Boyko and Stone2021).

4.3. Self-propulsion in biphasic systems

Another major application of the reciprocal theorem is in modelling the self-propulsion of active biological or synthetic entities at low Reynolds flows (Becker et al. Reference Becker, Koehler and Stone2003; Lauga & Powers Reference Lauga and Powers2009; Masoud & Stone Reference Masoud and Stone2014; Elfring Reference Elfring2015; Li et al. Reference Li, Lauga and Ardekani2021). As with the previous examples, the framework in Stokes flow can be directly extended to biphasic systems. We focus on the general case, where the active interactions between the motile object with the biphasic medium leads to slip velocities on the surface of the object with respect to the network and the fluid phases. The exact form of this slip velocity depends very much on the specifics of the involved active processes in smaller scales. For the purpose of this discussion, we assume the slip velocities in both phases, $\boldsymbol {v}_{\{\,f,n\}}^s$, are known. Our goal is to find the net propulsion speed of the force-free particle, $\boldsymbol {U}$, due to these surface activities. Hence, the velocity on the particle's surface is $\boldsymbol {v}_{\{\,f,n\}}=\boldsymbol {U}+\boldsymbol {v}^s_{\{\,f,n\}}$.

We use the solution to a particle moving with ${\hat {\boldsymbol {U}}}$ under a no-slip BC as the auxiliary solution, and seek to find $\boldsymbol {U}$. After substitution we get

(4.11)\begin{equation} {\hat{\boldsymbol{F}}} \boldsymbol{\cdot} \boldsymbol{U}={-}\int_{S_p} (\boldsymbol{v}^s_{f}\boldsymbol{\cdot} {\hat{\boldsymbol{\varSigma}}}_f + \boldsymbol{v}^s_{n}\boldsymbol{\cdot} {\hat{\boldsymbol{\varSigma}}}_n )\boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d}S, \end{equation}

where we used the force-free condition on the motile particle ($\int _{S_p} (\boldsymbol {\varSigma }_f +\boldsymbol {\varSigma }_n )\boldsymbol {\cdot } \boldsymbol {n} \, \mathrm {d} S=\boldsymbol {0}$) to derive the final expression. This expression holds for any particle shape, but is only useful if the solution to the auxiliary problem is either known or can be obtained more easily using analytical and numerical methods. This formulation is also extendable to a collection of particles (Elfring Reference Elfring2015; Papavassiliou & Alexander Reference Papavassiliou and Alexander2015).