Dynamics of suspended rigid aggregating particles in flowing medium: theory, analysis and scientific computing

Figure 1a illustrates the adhesion of two particles which are modelled as two rigid spheres, $s_{i}$ (with radius $r_{i}$) and $s_{j}$ (with radius $r_{j}$), respectively. To simplify the visualisation of the micro-level dynamics, consider a moving frame, $R_{m}$, with origin $O_{m}$ (sometimes referred as the Lagrangian frame of reference) fixed on the surface of the sphere $s_{i}$ at a point equidistant from the edge of the separation gap. The unit vectors for this frame of reference are $\text{\boldmath$e$}_{x}$, $\text{\boldmath$e$}_{y}$ and $\text{\boldmath$e$}_{z}$. For a given spatial point ${\text{\boldmath$x$}}=(x,y,z)$ in this moving frame, the total relative velocity (of sphere $s_{i}$ with respect to $s_{j}$) is ${\text{\boldmath$v$}}=u(\text{\boldmath$U$})\text{\boldmath$e$}_{x}+r_{i}\omega(\text{\boldmath$U$})\text{\boldmath$e$}_{y}$, where the relative translational velocity, $u$, and the rotational velocity $\omega$, are functions of the fluid velocity, $U$ (§2.2.1). Let $d(x)$ be the vertical separation gap between the two spheres. Define $A_{\text{Tot}}g(x)\,dA$ as the number of bonds that are attached between the surfaces $A$ and $A+dA$ at time $t$. $A_{\text{Tot}}$ is the total number of binding ligands. The function $g$ is synonymous with the term sticking probability [46]. The total number of bonds formed is $\int_{a_{s}}A_{\text{Tot}}g(x)\,dA$ , where $a_{s}=\pi r^{2}_{s}$ is the area of adhesion (refer §2.1.3 for details on the radius of adhesion, $r_{s}$). In further description of the model we denote $d(x)\equiv d$, without loss of generalisation.

The forward and reverse reaction rates for the binding ligands are then written as Boltzmann distributions, allowing highly stretched bonds to be readily broken by thermal energy fluctuations. The kinetics are also influenced by the surface potential of the two charged surfaces. Further, we cater for the ligands tilting by a finite angle $\alpha_{0}$ with respect to the vertical direction. This tilt is again expressed as a Boltzmann distribution, $\mathcal{D}(\alpha_{0})$, such that a bond may form between the two spheres for a given angle $\alpha_{0}\in(-\frac{\pi}{2},\frac{\pi}{2})$. With these degrees of freedom, the bond attachment/detachment rates are

	$\displaystyle K_{\text{on}}(x)$	$\displaystyle=K_{\text{on,eq}}\exp\left[\frac{-\lambda_{s}(l(x)-{l}_{0})^{2}+W(d)}{2{k}_{B}T}\right]\mathcal{D}(\alpha_{0}),$
	$\displaystyle K_{\text{off}}(x)$	$\displaystyle=K_{\text{off,eq}}\exp\left[\frac{(\lambda_{0}-\lambda_{s})(l(x)-{l}_{0})^{2}+W(d)}{2{k}_{B}T}\right],$		(1)

where ${k}_{B}$ is the Boltzmann constant, $T$ is the temperature, $l_{0}$ is the mean rest length of the binders, $\lambda_{0}$ is the binder stiffness coefficient, and $\lambda_{s}$ is the spring constant of the transition state used to distinguish catch ($\lambda_{0}<\lambda_{s}$) from slip ($\lambda_{0}>\lambda_{s}$) bonds [7]. $W(d)$ is the total surface potential described in §2.1.2. As depicted in Figure 1a, $l=\sqrt{d^{2}+x^{2}}$ is the length of a bond in a stretched configuration. The energy associated with tilting a bond from its vertical position is $(1/2)\lambda_{\theta}\alpha_{0}^{2}$, ($\lambda_{\theta}$ being the torsion constant) and the corresponding Boltzmann distribution is

$$\mathcal{D}(\alpha_{0})=\exp\left(-\frac{\lambda_{\theta}\alpha_{0}^{2}}{2{k}_{B}T}\right)\frac{1}{D_{0}}\,,\quad\alpha_{0}=\tan^{-1}\frac{x}{d}\,,$$

(2)

where $D_{0}=\int^{{\pi}/{2}}_{-{\pi}/{2}}\exp\big{[}{-\frac{\lambda_{\theta}\alpha^{2}_{0}}{2{k}_{B}T}}\big{]}\,d\alpha_{0}=\sqrt{\frac{2\pi k_{B}T}{\lambda_{\theta}}}\text{erf}\left(\frac{\pi}{2}\sqrt{\frac{\lambda_{\theta}}{2k_{B}T}}\right)$ is the normalization constant for all possible tilt orientations along the flow-direction. Ignoring non-equilibrium binding kinetics (i.e., $\frac{\partial g}{\partial t}=0$), in the limit of small binding affinity and abundant ligands on the binding surface (i.e., $K_{\text{eq}}=\nicefrac{{(A_{\text{Tot}}K_{\text{on,eq}})}}{{K_{\text{off,eq}}}}\ll 1$), the evolution equation for the sticking probability is [7, 35],

$${\text{\boldmath$v$}}\cdot\nabla g=A_{\text{Tot}}K_{\text{on}}-K_{\text{off}}g\,,\quad g=0\quad\text{for }x\geq r_{s}\,.$$

(3)

For two charged spheres, of radii $r_{i},r_{j}$, the potential due to the repulsive Coulombic forces in the gap of size $d$ is [13]

$$W_{\text{C}}(d)=2\pi\epsilon_{0}\epsilon\psi^{2}_{0}\Big{(}\frac{2r_{i}r_{j}}{r_{i}+r_{j}}\Big{)}e^{-\delta d},$$

(4)

where $\delta$ is the Debye length, $\epsilon$ and $\epsilon_{0}$ the dielectric constant of vacuum and the medium, respectively, and $\psi_{0}$, the average zeta potential of the diffuse cloud of charged counterions. The potential due to the Van der Waal forces for these spheres in the regime of close contact is

$$W_{\text{VW}}(d)=-\frac{A}{6d}\frac{r_{i}r_{j}}{r_{i}+r_{j}},$$

(5)

where A is the Hamaker constant, measuring the Van der Waal ‘two-body’ pair-interaction for macroscopic spherical objects. The total surface potential is $W(d)=W_{\text{C}}(d)+W_{\text{VW}}(d)$, which is pair-wise attractive over short and long range, and pair-wise repulsive over intermediate range.

First, the mesoscale hydrodynamic forces, $F$, and torques, $T$, on two rigid spheres in Stokes’ flow are listed. These were derived by O’Neill  [29] using lubrication theory. At leading order in $\epsilon=\frac{d}{r_{i}}(\ll 1)$, these expressions are

		$\displaystyle{\text{\boldmath$F$}}^{t}=$	$\displaystyle-\frac{8\pi\mu r_{i}u}{5}\frac{\zeta(2+\zeta+2\zeta^{2})}{(1+\zeta)^{3}}\ln(\epsilon){\bf e}_{x},$		(6)
		$\displaystyle{\text{\boldmath$T$}}^{t}=$	$\displaystyle-\frac{4\pi\mu r^{2}_{i}u}{5}\frac{\zeta(4+\zeta)}{(1+\zeta)^{2}}\ln(\epsilon){\bf e}_{y},$
		$\displaystyle{\text{\boldmath$F$}}^{r}=$	$\displaystyle-\frac{4\pi\mu r^{2}_{i}\omega}{5}\frac{\zeta(4+\zeta)}{(1+\zeta)^{2}}\ln(\epsilon){\bf e}_{x},$
		$\displaystyle{\text{\boldmath$T$}}^{r}=$	$\displaystyle-\frac{16\pi\mu r^{3}_{i}\omega}{5}\frac{\zeta}{(1+\zeta)}\ln(\epsilon){\bf e}_{y},$
		$\displaystyle{\text{\boldmath$F$}}^{f}=$	$\displaystyle\frac{3\pi\mu r_{i}U(25\zeta+16)}{4\zeta}\Big{(}\frac{\zeta}{1+\zeta}\Big{)}^{2}{\bf e}_{x},$
		$\displaystyle{\text{\boldmath$T$}}^{f}=$	$\displaystyle 8\pi\mu r^{2}_{i}U\left(\frac{\zeta}{1+\zeta}\right)^{3}\left(1-\frac{3}{16}\left(\frac{\zeta}{1+\zeta}\right)^{3}\right){\bf e}_{y},$

where $U$ is the cross sectionally averaged fluid flow (§2.2.1), $\zeta=\frac{r_{i}}{r_{j}}$ and $\mu$ is the dynamic viscosity of the fluid. The superscripts $(t,r,f)$ denote components arising from relative translation, rotation and the flow of the surrounding fluid, respectively. These hydrodynamic forces must balance the forces arising from all the bond interactions,

	$\displaystyle{\text{\boldmath$F$}}_{E}(u,\omega,d)=$	$\displaystyle-\int_{a_{s}}g(1-1/l)\left[x\text{\boldmath$e$}_{x}+d\text{\boldmath$e$}_{z}\right]dA,$
	$\displaystyle{\text{\boldmath$F$}}_{C}(u,\omega,d)=$	$\displaystyle\int_{a_{s}}g\frac{r_{i}r_{j}}{r_{i}+r_{j}}\left(-\frac{4\pi\epsilon_{0}\epsilon\psi_{0}}{\delta}\exp^{-\delta d}+\frac{A}{6{d}^{2}}\right){\text{\boldmath$e$}}_{z}dA,$
	$\displaystyle{\text{\boldmath$F$}}_{T}(u,\omega,d)=$	$\displaystyle(\nicefrac{{\lambda_{\theta}}}{{(2\lambda_{0}{\it l}_{0}^{2})}})\int_{{a_{s}}}({g}/{l}^{2})\alpha_{0}\left[-{d}{\text{\boldmath$e$}}_{x}+{x}{\text{\boldmath$e$}}_{z}\right]dA.$			(7)

The subscripts $i=E,C,T$ denote extension, surface charges and torsion, respectively. The adhesion area of the circular patch is given by $a_{s}=\pi r^{2}_{s}$, where the adhesion radius, $r_{s}$, is found by Mahadevan (details in supplementary material, [21]) using a scaling law argument of the ‘settling phase of the particles’,

$$r_{s}=2\left(\frac{{\it k}_{B}T{\it l}_{0}}{\lambda_{0}}\right)^{1/2}\left(\frac{1}{r_{i}}+\frac{1}{r_{j}}\right)(r_{i}+r_{j})^{1/2}.$$

(8)

The coupling across the scales is obtained via the global force balance on the spheres in the horizontal and vertical directions, and the torque balance about the center of mass of the spheres (i.e., forces and torques from eqn. (6) with forces from eqn. (7)), which solves for the unknown variables, $u,\omega,d$, i.e.,

	$\displaystyle({\bf F}_{E}+{\bf F}_{C}+{\bf F}_{T})\cdot{\bf e}_{z}=0,$
	$\displaystyle({\bf F}_{E}+{\bf F}_{C}+{\bf F}_{T})\cdot{\bf e}_{x}+(F^{r}+F^{t}+F^{s})=0,$
	$\displaystyle\frac{\beta r_{i}}{\zeta+1}({\bf F}_{E}+{\bf F}_{C}+{\bf F}_{T})\cdot{\bf e}_{x}+(T^{r}+T^{t}+T^{s})=0.$		(9)

Finally, in a Stokes regime, the collision frequency function, $\beta$ (eqn. (15) in §2.2), is proportional to the total microscale force, ${\text{\boldmath$F$}}_{\text{Tot}}={\text{\boldmath$F$}}_{E}+{\text{\boldmath$F$}}_{C}+{\text{\boldmath$F$}}_{T}$, given by [41]

$$\beta=\gamma_{A}\|{\text{\boldmath$F$}}_{\textrm{Tot}}\|,$$

(10)

where $\gamma_{A}=\frac{\lambda_{0}}{A_{\text{Tot}}K_{eq}A\pi\mu}$ is known as the aggregation contact efficiency parameter [46]. The operator $\|\cdot\|$ denotes the magnitude of a vector. We remark that in the case of static equilibrium ($u=\omega=0$) with no bond tilting effects on the binders (${\lambda_{\theta}}^{-1}\rightarrow 0$) the evolution of $g$ (eqn. (3)) trivially reduces to

$$g=\nicefrac{{A_{\text{Tot}}K_{\text{on}}}}{{K_{\text{off}}}}=K_{eq}\exp^{-\lambda_{0}(d-{\it l}_{0})^{2}},$$

(11)

and further, with perfect binding ($l=d\rightarrow{\it l}_{0}$ so that $g\rightarrow K_{eq}$), the collision frequency reduces to

	$\displaystyle\beta(\mathit{u},\mathit{v})_{d\rightarrow{\it l}_{0}}$	$\displaystyle=\gamma_{A}\|\sum_{i=E,C}F_{i}\|\pi r^{2}_{s}$
		$\displaystyle=\frac{4k_{B}T}{\mu}\frac{(\mathit{u}^{1/3}+\mathit{v}^{1/3})^{2}}{\mathit{u}^{1/3}\mathit{v}^{1/3}}\frac{1}{{\it l}_{0}A}\left(-\frac{4\pi\epsilon_{0}\epsilon\psi^{2}_{0}}{\delta{\it l}_{0}}\exp^{-\delta{\it l}_{0}}+\frac{A}{6{\it l}^{2}_{0}}\right),$		(12)

where $\mathit{u}=\frac{4\pi}{3}r^{3}_{i}$ and $\mathit{v}=\frac{4\pi}{3}r^{3}_{j}$ are the respective volumes of the spheres. In a neutral solution with no dissolved ions (i.e., the Debye length, $\delta^{-1}\rightarrow 0$), we have

$$\beta(\mathit{u},\mathit{v})_{d\rightarrow{\it l}_{0},\delta^{-1}\rightarrow 0}=\beta^{Br}=\frac{2{\it k}_{B}T}{3\mu{\it l}^{3}_{0}}\frac{(\mathit{u}^{1/3}+\mathit{v}^{1/3})^{2}}{\mathit{u}^{1/3}\mathit{v}^{1/3}},$$

(13)

which is the collision frequency function for Brownian diffusion for hard spheres commonly used in many theoretical and experimental studies of aggregation [10, 13]. To summarize, the micro-scale surface adhesion model of two particles has the following features:

• •

  Each particle constitutes a rigid spherical core onto which linear, hookean, spring-like binding ligands are attached. The surface of the coalescing particles are linked through these ligands. The ligand kinetics is mediated by the bond formation/breakage rates and modelled using eqn. (3).

• •

  The rigid core of the particles, is charged and suspended in an ionic medium. These effects are modelled via the repulsive Coulombic interactions (eqn. (4)) and the attractive Van der Waal interactions (eqn. (5)).

• •

  The coupling between dynamics at the micro(particle)-level and the meso(cluster)-level, is achieved through the global balance of the meso-scale drag forces and couples (eqn. (6)) with the micro-scale forces associated with the binder kinetics (eqn. (7)).

This section justifies and constructs a model of the long-term dispersion of the particle aggregates along a long pipe that is carried by Poiseuille flow down the pipe [23]. This approach, supported by center manifold theory, can be readily extended to cater for reactive [49] or sedimenting material [47] in long pipes of complex and varying geometry [23].

Consider the particle aggregates being carried by a fluid flow in a round pipe of radius $a$ (Figure 1c). We assume that the particles are neutrally buoyant, i.e., the particles have no effect on the density of the fluid so that concentration differences do not cause density differences that could generate secondary flows. The along pipe flow is the Poiseuille flow, $\text{\boldmath$U$}(R)=2U(1-R^{2}/a^{2})\text{\boldmath$e$}_{X}$, where $U$ is the cross-sectionally averaged downstream velocity (defined later). Using the pipe radius, $a$, and the time associated with cross-pipe diffusion, $\tau=a^{2}/D$, as the characteristic length and time scale, the particle conservation pde (15) reduces to

	$\displaystyle\frac{\partial n}{\partial t}+U(R)\frac{\partial n}{\partial X}$	$\displaystyle=\frac{1}{R}\frac{\partial}{\partial R}\left(R\frac{\partial n}{\partial R}\right)+\frac{1}{R^{2}}\frac{\partial^{2}n}{\partial\theta^{2}}+\frac{\partial^{2}n}{\partial X^{2}}+\widetilde{\phi(n)}$
		$\displaystyle=\mathcal{L}(n)+\frac{\partial^{2}n}{\partial X^{2}}+\widetilde{\phi(n)},$		(16)

for pouiseuille flow. $\widetilde{\phi(n)}=\frac{a^{2}}{D}\phi(n)$. The boundary condition is that of zero diffusive flux on the wall,

$$\frac{\partial n}{\partial R}=0\quad\text{on}\quad R=1.$$

(17)

Note that the equilibria, $n^{*}$, of the pde (16) satisfies $\widetilde{\phi(n^{*})}=0$. We find a slow manifold about such an equilibria and hence deduce a model global in size distribution and local to small $\partial/\partial x$, corresponding to the low wavenumber solution in Fourier space. Choosing to write the model in terms of the cross-pipe average, to an initial approximation, the slow manifold model is the trivial

$$n^{0}\approx\overline{n}(X,t)\quad\text{such that}\quad\frac{\partial\overline{n}}{\partial t}\approx\widetilde{\phi(\overline{n}}),$$

(18)

where the amplitude $\overline{n}$ is defined via the projection, $\overline{n}(X,t)\!=\!\!\int_{0}^{1}\!\!\!\int_{0}^{2\pi}\!\!\!n(R,\!\theta,\!X,\!t)\!Rd\theta dR$​, which is the cross-pipe average. Now seek to iteratively refine this model (18), using the corrections ${\hat{n}}(N,R,\theta)$ and ${\hat{g}}(N)$, where

$$n^{1}\approx n^{0}+{\hat{n}}\quad\text{such that}\quad\frac{\partial\overline{n}}{\partial t}\approx\widetilde{\phi}+{\hat{g}}(\overline{n})\,,$$

(19)

is a more refined model. Substituting into the advection diffusion pde (16) and neglecting products of the small corrections (i.e., $-U\partial\hat{n}/\partial X,\partial^{2}\hat{n}/\partial X^{2}$) gives

$$\frac{\partial n^{1}}{\partial t}=\mathcal{L}(n^{1})-U(r)\overline{n}_{X}+\overline{n}_{XX}+\widetilde{\phi},$$

(20)

where the subscript, $X$, denotes the longitudinal derivatives (see Figure 1c). Using the approximation (19), the chain rule, ${\partial n^{1}}/{\partial t}=\left({\partial n^{1}}/{\partial N}\right)\left({\partial N}/{\partial t}\right)\approx\widetilde{\phi}+\hat{g}$ and the identity $\mathcal{L}(\overline{n})=0$, we have

$$\mathcal{L}(\hat{n})-\hat{g}=U(r)\overline{n}_{X}-\overline{n}_{XX}.$$

(21)

Note that $\mathcal{L}$ is self-adjoint, upon defining the inner product $\langle v,w\rangle=\int_{0}^{1}\int_{0}^{2\pi}vw\,Rd\theta\,dR$. An adjoint eigenvector of this operator is $z(R,\theta)\!=$constant. Projecting eqn. (21) in the space spanned by this eigenvector, gives

$$-\pi{\hat{g}}=\overline{n}_{X}\int\!\!\!\int U(R)\,R\,d\theta\,dR-\pi\overline{n}_{XX}\,,$$

which simplifies to

$${\hat{g}}=-U\overline{n}_{X}+\overline{n}_{XX}\,,$$

(22)

as $U(R)=2U(1-R^{2})$ has cross-sectional average of $U$ by the non dimensionalisation. Using this update ${\hat{g}}$ and the fact that there are no $\theta$ variations in the right-hand side of eqn. (21), the first correction ${\hat{n}}$ can be found as ${\hat{n}}=w_{1}(R)U\overline{n}_{X}$ where $w_{1}(R)=\frac{1}{24}\left(-2+6R^{2}-3R^{4}\right)$. The arbitrary constant of integration is found from the constraint that the cross-sectional average of every correction ${\hat{n}}$ has to be zero. Correcting the initial model (18) with these ${\hat{n}}$ and ${\hat{g}}$ we deduce a refined model

$$n^{1}\approx\overline{n}+w_{1}(R)U\overline{n}_{X}\quad\text{such that}\quad\frac{\partial\overline{n}}{\partial t}\approx-U\overline{n}_{X}+\overline{n}_{XX}+{\tilde{\phi}}\,.$$

(23)

The Approximation Theorem suggests that the $N_{XX}$ term needs correction because the residuals in the equation resulting from this approximate model, (see the iterate, eqn. (24b)), are of order 2 [23]. Hence, following the above procedure, a second level of refinement leads to

	$\displaystyle n\approx\overline{n}+w_{1}(R)U\overline{n}_{X}+w_{2}(R)U^{2}\overline{n}_{XX}$		(24a)
	$\displaystyle\frac{\partial\overline{n}}{\partial t}\approx-U\overline{n}_{X}+\left(1+\frac{U^{2}}{48}\right)\overline{n}_{XX}+{\widetilde{\phi}}$		(24b)

where $w_{2}(R)=\frac{1}{11520}\left(31-180R^{2}+300R^{4}-200R^{6}+45R^{8}\right)$. Higher order corrections of the estimate in eqn. (24b) can be easily derived using computer algebra.¹¹1http://reduce-algebra.com/ The reaction diffusion eqn. (24b) has an effective diffusion coefficient of $1+U^{2}/{48}$, i.e., the effective rate of dispersion along the pipe is faster than what is predicted by simple cross-sectional averaging. The dimensional form of the model (24b), is

$$\frac{\partial\overline{n}}{\partial t}\approx-U\overline{n}_{X}+\left(D+\frac{U^{2}a^{2}}{48D}\right)\overline{n}_{XX}+\phi(\overline{n})\,.$$

(25)

The extra dispersion (with coefficient $U^{2}a^{2}/(48D)$), also called the shear dispersion, is proportional to the square of the velocity and inversely proportional to the cross-pipe diffusivity and implies that the smaller the molecular diffusion, the larger the shear dispersion, a paradox which is well established [38].

We construct $j$th patch of width $(2k+1)h$ and denote it with macro-scale lattice point $X_{j}$ (Figure 2a). $k$ is the patch half-width and $h$ is the meso-scale lattice spacing. First step, for each point within a patch and for every time-step, the non-linear aggregation term, $\phi(n)$ (eqn. (14)), is calculated. The integrals (eqn. (7)) in the collision frequency function (eqn. (10)) using the adaptive Lobatto quadrature (via Matlab function quadl). The function, $\phi(n)$, is numerically evaluated by employing the Galerkin approximation scheme developed by Banks [18] and adopted by others [42]. While evaluating the integrals in eqn. (14), we chose $\mathit{u}_{\text{min}}=1$ femtoliter (fL) as a lower bound in our simulations, since the radius of the smallest aggregate (consisting of a single particle) was of the order of a micron. In our discrete aggregation model we chose the upper bound, $\mathit{u}_{\text{max}}=2000$ fL and the volume resolution of $N=100$ bins per fL. Other parameters used in the simulation are listed in Table 1 (see §3). As reported in [18], the convergence of this first order approximation scheme was tested and a linear (first order) relationship was found between the $l^{\infty}$-error of the solution and the mesh-size, $\delta u$ ($\delta u$ being the discretisation in the aggregate volume). The non-dimensional advection-reaction-diffusion equation (24b) is discretised using the standard Crank-Nicolson scheme. The non-linear reaction term is treated explicitly. Although the size of the domain is arbitrarily fixed at $K=100$ patches, the reason for the choice the patch half-width, $k$, buffer half-width, $b$, and the lattice spacings, $h,H$, is explained in §3. Physical inlet and exit conditions are needed for the macro-scale solution in the first patch, $X_{1}$, and the last patch, $X_{K}$, respectively. The conditions are [37],

	$\displaystyle N_{1}(t)$	$\displaystyle\approx\frac{1}{\pi\kappa}e^{-\frac{t}{\kappa}}\!\ast\!\overline{U(r)n_{1}(t)}+\frac{1}{\pi\kappa}\left[\overline{U(r)r_{1}n_{1}(t)}-\frac{1}{\kappa}e^{-\frac{t}{\kappa}}\!\ast\!\overline{U(r)r_{1}n_{1}(t)}\right]$		(26a)
	$\displaystyle\frac{\partial N}{\partial X}$	$\displaystyle\approx e^{-\frac{t}{\kappa}}\ast\frac{\partial N}{\partial t},\qquad\text{at}~j=K$		(26b)

where $\kappa=2/105$, $n_{1}(t)$ is the prescribed particle size distribution across the domain entrance. The overbar ( $\overline{\phantom{w}}$ ), denotes the cross-stream average and the operator $\ast$ denotes the convolution, $f(t)\ast g(t)=\int_{0}^{t}f(\tau)g(t-\tau)d\tau$. The functions $r_{1}$ is given by

$$r_{1}(R)=\frac{1}{280}\left(35R^{4}-70R^{2}+11\right),$$

(27)

The inlet condition (26a) ensures that the second order shear dispersion model (24) is accurate after the entry transients have disappeared and the exit condition (26b) ensures that the exit does not have an upstream influence (via diffusion) [37]. In following discussions, we drop the overbar over the variables, so that $n=n(x,t)$ is the cross-sectionally averaged particle size distribution function.

After eqn. (25) is solved within the $j$th patch for $n_{ij}$, together with the coupling conditions (§2.3.3), we define a macro-scale solution field within each of these patches. We propose a macro-scale solution, referred as the patch amplitude from hereon, as an average over the meso-scale field within the patch ‘core’ (Figure 2b),

$$N_{j}(t)=N(X_{j},t)=\sum^{i=b}_{i=-b}\frac{n_{ij}}{2b+1},$$

(28)

where $b$ is the patch-buffer half width, such that $0\leq b<\,$k. Some patch dynamics methods extract the macro-scale solution by trivially using the meso-scale solution at the patch center, corresponding to $b=0$ [36]. However, unlike the ‘rough’ micro-scale structure in the present setup, these methods are for ‘smooth’ micro-scale dynamics with only slight variations within patches.

Two coupling conditions are required to solve eqn. (25), at each end within the $j$th patch. These conditions are derived from an interpolation of the macro-scale field, $N_{j}$, from its nearest macro-scale neighbours, $N_{j\pm 1},N_{j\pm 2},\ldots$ To write these conditions, we define a fractional step $r$, connecting the spacing between the meso-scale domain, h, and the macro-scale domain, H, such that $rH=(k-b)h$ (Figure 2b). We also define the meso-scale and the macro-scale step operators, $\varepsilon,\overline{\varepsilon}$, such that $\varepsilon^{\pm 1}n_{ij}=n_{i\pm 1j}$ and $\overline{\varepsilon}^{\pm 1}N_{j}=N_{j\pm 1}$, respectively. Thus, $\varepsilon^{\pm(k-b)}N_{j}=\overline{\varepsilon}^{\pm r}N_{j}$. Finally, we introduce the macro-scale mean and difference operators, $\overline{\mu}=(\overline{\varepsilon}^{1/2}+\overline{\varepsilon}^{-1/2})/2$ and $\overline{\delta}=(\overline{\varepsilon}^{1/2}-\overline{\varepsilon}^{-1/2})/2$, respectively. Using the condition (28), note that

	$\displaystyle\varepsilon^{\pm(k-b)}N_{j}$	$\displaystyle=\sum^{i=k}_{i=k-2b}\frac{n_{\pm ij}}{2b+1}=\overline{\varepsilon}^{\pm r}N_{j}$
		$\displaystyle=(\overline{\mu}\overline{\delta}+\frac{1}{2}\overline{\mu}^{2}+1)^{\pm r}N_{j}$
		$\displaystyle\approx N_{j}+\sum_{q=1}^{\Gamma}\prod_{l=0}^{q-1}(r^{2}-l^{2})\frac{\pm(2q/r)\overline{\mu}\overline{\delta}^{2q-1}+\overline{\delta}^{2q}}{(2q)!}N_{j}+\mathcal{O}({\Gamma+1}),$		(29)

which are the required coupling conditions. The $\approx$ can be replaced with an equality provided $\Gamma\rightarrow\infty$. Detailed discussion on the theoretical support for patch dynamics are provided elsewhere [36]. Finally, it has been shown that it is possible define arbitrary patch coupling conditions, provided suitably large buffers are chosen to shield the patch ‘core’ from unwanted consequences of these coupling conditions over short evaluation times [39]. To summarise, given the macro-scale solution, $N_{j}(t_{n-1})$ at time $t_{n-1}$, the patch dynamics algorithm for finding the macroscale solution , $N_{j}(t_{n})$ at time $t_{n}$, is

1. 1.

   Using the micro-scale field variables evaluated at time $t_{n-1}$, setup the non-linear reaction term, $\phi(n)$ (eqn. (14)), within each patch $j$.

2. 2.

   Numerically solve eqn. (25), together with coupling conditions (29) within all patches $j$, such that $N_{j}=N_{j}(t_{n-1})$.

3. 3.

   From the amplitude condition (28), determine the macro-scale solution at time $t_{n}$ for all macro-scale points $X_{j}$ (Figure 2a), such that $n_{ij}=n_{ij}(t_{n}),N_{j}=N_{j}(t_{n})$.

4. 4.

   Using the macro-scale solution $N_{j}(t_{n})$ evaluated in Step 3, return to Step 1 to determine the solution at time step $t_{n+1}$.

The computational technique is ‘equation-free’ since no attempt is made to derive an explicit macro-scale model. The procedure simply evaluates ‘on the fly’ numerical solutions of the model at each macro-scale point, $X_{j}$, and time-step, $t_{n}$.