2.1 The fundamental periodic solution

We study the case of transport to a planar reactive surface with periodic planar arrays of inert disks of unit radius. The problem formulation and solution method are exactly the same as those in Shah & Shaqfeh (Reference Shah and Shaqfeh2015). Here, we reproduce brief relevant parts of Shah & Shaqfeh (Reference Shah and Shaqfeh2015) to introduce the notation used throughout this work. The heterogeneous surface ${\mathcal{B}}$ is subdivided into identical boxes ${\mathcal{B}}_{\boldsymbol{n}}$ (so that ${\mathcal{B}}=\cup {\mathcal{B}}_{\boldsymbol{n}}$ ), where $\boldsymbol{n}\in \mathbb{Z}^{2}$ . Each such unit has the same configuration of reactive and inert domains, denoted by ${\mathcal{P}}_{\boldsymbol{n}}$ and ${\mathcal{W}}_{\boldsymbol{n}}$ respectively (see figure 2). The entire reactive region is ${\mathcal{P}}=\cup {\mathcal{P}}_{\boldsymbol{n}}$ , while the inert region is ${\mathcal{W}}=\cup {\mathcal{W}}_{\boldsymbol{n}}$ . The area of the reactive region in each box is denoted by $P_{0}=|{\mathcal{P}}_{\boldsymbol{n}}|$ , while that for the inert region is $W_{0}=|{\mathcal{W}}_{\boldsymbol{n}}|$ . The centres of each box ${\mathcal{B}}_{\boldsymbol{n}}$ form a lattice structure,

(2.1) $$\begin{eqnarray}\boldsymbol{r}_{\boldsymbol{n}}=n_{1}\boldsymbol{L}_{1}+n_{2}\boldsymbol{L}_{2}\quad (n_{1},n_{2}=0,\pm 1,\pm 2,\ldots ),\end{eqnarray}$$

where $\boldsymbol{L}_{1}$ and $\boldsymbol{L}_{2}$ are steps along periodic directions on the heterogeneous surface. All periodic quantities have subscript ‘ $p$ ’. Following (1.2c ) and the asymptotic condition (1.6), the periodic concentration field $C_{p}$ satisfies

(2.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}_{x}^{2}C_{p}-Px_{3}\frac{\unicode[STIX]{x2202}C_{p}}{\unicode[STIX]{x2202}x_{1}}=0, & \displaystyle\end{eqnarray}$$

(2.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}C_{p}}{\unicode[STIX]{x2202}x_{3}}=\unicode[STIX]{x1D705}C_{p}\quad \text{at }\boldsymbol{x}\in {\mathcal{P}}_{\boldsymbol{n}},\boldsymbol{n}\in \mathbb{Z}^{2}, & \displaystyle\end{eqnarray}$$

(2.2c ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}C_{p}}{\unicode[STIX]{x2202}x_{3}}=0\quad \text{at }\boldsymbol{x}\in {\mathcal{W}}_{\boldsymbol{n}},\boldsymbol{n}\in \mathbb{Z}^{2}, & \displaystyle\end{eqnarray}$$

(2.2d ) $$\begin{eqnarray}\displaystyle & \displaystyle C_{p}\rightarrow 1+\frac{P_{0}}{B_{0}}S_{p}(\unicode[STIX]{x1D705},P,\unicode[STIX]{x1D719})x_{3}\quad \text{as }x_{3}\rightarrow \infty . & \displaystyle\end{eqnarray}$$

Figure 2.

A periodic system of reactive region ${\mathcal{P}}_{\boldsymbol{n}}$ with inert disks ${\mathcal{W}}_{\boldsymbol{n}}$ embedded in corresponding bounding boxes ${\mathcal{B}}_{\boldsymbol{n}}$ . Shear flow of strength $\dot{\unicode[STIX]{x1D6FE}}$ is present in the wall-normal direction. Also illustrated are the periodicity directions $\boldsymbol{L}_{1},\boldsymbol{L}_{2}$ in the lattice.

The Sherwood number, $S_{p}$ , is identical for each porous region ${\mathcal{P}}_{\boldsymbol{n}}$ , and is given by the expression

(2.3) $$\begin{eqnarray}S_{p}(\unicode[STIX]{x1D705},P,\unicode[STIX]{x1D719})=\frac{1}{|{\mathcal{P}}_{\boldsymbol{n}}|}\int _{{\mathcal{P}}_{\boldsymbol{n}}}\unicode[STIX]{x1D705}C_{p}(\boldsymbol{x}_{b})\,\text{d}A(\boldsymbol{x}_{b}),\end{eqnarray}$$

where the subscript ‘ $b$ ’ on the dummy variable indicates that point $\boldsymbol{x}_{b}$ is bound to the surface ${\mathcal{B}}$ . Shah & Shaqfeh (Reference Shah and Shaqfeh2015) devised a BEM that uses the Green’s function, $G_{p}$ , that satisfies the adjoint problem described in (2.4a–c ) below. It should be noted that the subscript ‘ $s$ ’ denotes the $_{1},_{2}$ coordinates.

(2.4a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}G_{p}+Px_{3}\frac{\unicode[STIX]{x2202}G_{p}}{\unicode[STIX]{x2202}x_{1}}=-\unicode[STIX]{x1D6FF}(x_{3}-y_{3})\left(\unicode[STIX]{x1D6FF}^{2}(\boldsymbol{x}_{s}-\boldsymbol{y}_{s,\boldsymbol{n}})-\frac{1}{B_{0}}\right), & \displaystyle\end{eqnarray}$$

(2.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}G_{p}}{\unicode[STIX]{x2202}x_{3}}=0\quad \text{at }\boldsymbol{x}\in {\mathcal{B}}_{\boldsymbol{n}}, & \displaystyle\end{eqnarray}$$

(2.4c ) $$\begin{eqnarray}\displaystyle & \displaystyle G_{p}\rightarrow 0\quad \text{as }x_{3}\rightarrow \infty . & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D6FF}^{2}(\boldsymbol{x}_{s}-\boldsymbol{y}_{s,\boldsymbol{n}})$

$\boldsymbol{y}_{\boldsymbol{n}}=(\boldsymbol{y}_{s,\boldsymbol{n}},y_{3})$

(2.5) $$\begin{eqnarray}\boldsymbol{y}_{s,\boldsymbol{n}}=\boldsymbol{y}_{s,0}+n_{1}\boldsymbol{L}_{1}+n_{2}\boldsymbol{L}_{2}\quad (n_{1},n_{2}=0,\pm 1,\pm 2,\ldots ).\end{eqnarray}$$

The renormalization term (Hasimoto Reference Hasimoto1959; Sangani & Behl Reference Sangani and Behl1989) $1/B_{0}$ balances the flux of singularity per periodic box. When the periodic singularities $\boldsymbol{y}_{b}$ (subscript $\boldsymbol{n}$ dropped for clarity) are bound to the surface ${\mathcal{B}}$ , we have (Shah & Shaqfeh Reference Shah and Shaqfeh2015)

(2.6) $$\begin{eqnarray}G_{p}(\boldsymbol{x}|\boldsymbol{y}_{b};P)=\frac{1}{B_{0}}\mathop{\sum }_{k\neq 0}\frac{-\text{Ai}(s)}{\text{Ai}^{\prime }(s_{0})(\text{i}Pk_{1})^{1/3}}\text{e}^{-\text{i}\boldsymbol{k}\boldsymbol{\cdot }(\boldsymbol{x}_{s}-\boldsymbol{y}_{s})},\end{eqnarray}$$

where

(2.7a,b ) $$\begin{eqnarray}s=\frac{k^{2}+\text{i}Pk_{1}x_{3}}{(\text{i}Pk_{1})^{2/3}}\quad \text{and }\quad s_{0}=\frac{k^{2}}{(\text{i}Pk_{1})^{2/3}},\end{eqnarray}$$

and $\boldsymbol{k}=(k_{1},k_{2})$ are vectors in the reciprocal lattice,

(2.8a ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{k}=2\unicode[STIX]{x03C0}n_{1}\boldsymbol{l}_{1}+2\unicode[STIX]{x03C0}n_{2}\boldsymbol{l}_{2}\quad (n_{1},n_{2}=0,\pm 1,\pm 2,\ldots ), & \displaystyle\end{eqnarray}$$

(2.8b ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{where}\quad \boldsymbol{l}_{1}=\frac{\boldsymbol{L}_{2}\times \hat{x}_{3}}{B_{0}},\quad \boldsymbol{l}_{2}=\frac{\hat{x}_{3}\times \boldsymbol{L}_{1}}{B_{0}}, & \displaystyle\end{eqnarray}$$

(2.8c ) $$\begin{eqnarray}\displaystyle & \displaystyle \text{so}\quad \boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{L}_{t}=2\unicode[STIX]{x03C0}n_{t}\quad (t=1,2). & \displaystyle\end{eqnarray}$$

Shah & Shaqfeh (Reference Shah and Shaqfeh2015) then derive a boundary integral equation for the periodic concentration field using the periodic Green’s function as the kernel,

(2.9) $$\begin{eqnarray}C_{p}(\boldsymbol{x})=1+\unicode[STIX]{x1D719}S_{p}x_{3}-\unicode[STIX]{x1D705}\int _{{\mathcal{P}}_{0}}G_{p}(\boldsymbol{y}_{b}|\boldsymbol{x};P)C_{p}(\boldsymbol{y}_{b})\,\text{d}A(\boldsymbol{y}_{b}).\end{eqnarray}$$

The Sherwood number $S_{p}$ is determined self-consistently by the following procedure.

1. (i) Solve (2.9) for $C_{p}(\boldsymbol{x}_{b})$ on ${\mathcal{P}}_{0}$ .

2. (ii) Obtain $S_{p}(\unicode[STIX]{x1D705},P,\unicode[STIX]{x1D719})$ using (2.3).

3. (iii) Generate $C_{p}(\boldsymbol{x})$ in the bulk as the left-hand side of (2.9).

The boundary integral equations are solved by meshing the reactive region over which the integral is defined and using the Galerkin approach to solve for nodal concentrations in the mesh. As the integral equations are periodic, we effectively solve by meshing the reactive region in only one periodic subunit.

2.2 Numerical implementation of the BEM

The Fourier series in (2.6) is evaluated at a grid point using Evald’s summation technique (Hasimoto Reference Hasimoto1959), as described in Shah & Shaqfeh (Reference Shah and Shaqfeh2015). The Green’s function at any point is then calculated using a bilinear interpolation between the values at the four grid points surrounding it. The integral equation (2.9) is solved numerically by converting to matrix equations using a Galerkin method on triangular elements for the reactive region, and this method differs from that presented in Shah & Shaqfeh (Reference Shah and Shaqfeh2015). Let the set of nodes be indexed by labels from the set ${\mathcal{N}}=\{i\}$ . The value of the concentration field at the node $i$ is denoted by $C_{i}$ , which is to be determined. We use piecewise continuous linear shape functions $\unicode[STIX]{x1D713}_{i}(\boldsymbol{x})$ centred at node $i$ . We solve the following set of linear equations for $C_{i}$ :

(2.10) $$\begin{eqnarray}\langle \unicode[STIX]{x1D713}_{j}(\boldsymbol{x}),1\rangle C_{j}+\unicode[STIX]{x1D705}\mathop{\sum }_{i\in {\mathcal{N}}}\langle \unicode[STIX]{x1D713}_{j}(\boldsymbol{x}),\langle G(\boldsymbol{x}_{0}|\boldsymbol{x};P),\unicode[STIX]{x1D713}_{i}(\boldsymbol{x}_{0})\rangle \rangle C_{i}=\langle \unicode[STIX]{x1D713}_{j}(\boldsymbol{x}),1\rangle ,\quad j\in {\mathcal{N}},\end{eqnarray}$$

where the standard $L_{2}$ inner product $\langle \cdot \,,\cdot \rangle$ is defined for periodic functions $f,g$ by an integral over the reactive region ${\mathcal{P}}_{0}$ inside one periodic box,

(2.11) $$\begin{eqnarray}\langle f(\boldsymbol{x}),g(\boldsymbol{x})\rangle =\int _{{\mathcal{P}}_{0}}f(\boldsymbol{x})g(\boldsymbol{x})\,\text{d}A(\boldsymbol{x}).\end{eqnarray}$$

The Green’s function has an integrable singularity at $\boldsymbol{x}=\boldsymbol{x}_{0}$ , and special care must be taken while evaluating the double inner product in (2.10) when $i=j$ . The surface concentration field and the Sherwood number are then approximately given by

(2.12) $$\begin{eqnarray}\displaystyle & \displaystyle C_{p}(\boldsymbol{x})\approx \mathop{\sum }_{j\in {\mathcal{N}}}C_{j}\unicode[STIX]{x1D713}_{j}(\boldsymbol{x}), & \displaystyle\end{eqnarray}$$

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle S_{p}(\unicode[STIX]{x1D705},P,\unicode[STIX]{x1D719})\approx \frac{\unicode[STIX]{x1D705}}{P_{0}}\mathop{\sum }_{j\in {\mathcal{N}}}\langle \unicode[STIX]{x1D713}_{j}(\boldsymbol{x}),1\rangle C_{j}. & \displaystyle\end{eqnarray}$$

Figure 3 shows the contours of the surface concentration field for a sample simulation with parameters $\unicode[STIX]{x1D705}=1,P=20$ and $\unicode[STIX]{x1D719}=0.95$ . One can observe that a concentration wake exists downstream of the inert patch. This wake is a concentration accretion region, similar in shape to the depletion region observed around the reactive disk in type ( ${\mathcal{R}}_{u}$ ) surfaces. The far-field surface concentration field is approximately $C_{p}\approx 1$ , while the concentration field inside the inert patch is significantly higher than everywhere else.

Figure 3.

Contour plot of surface concentration for the case of periodic inert disks embedded in a reactive wall; $\unicode[STIX]{x1D705}=1,P=20,\unicode[STIX]{x1D719}=0.95$ .

Before we investigate the results of the BEM simulations, we shall derive a scaling law for the Sherwood number at $\unicode[STIX]{x1D719}\rightarrow 1$ . This will provide physical intuition for the findings of the BEM simulations. Importantly, we note that since the radius of the inert disks is fixed to unity, the $\unicode[STIX]{x1D719}\rightarrow 1$ limit is approached by taking $B_{0}\rightarrow \infty$ .

4.1 Concentration field $C_{u}$ for a single inert patch in a reactive surface

At $P=0$ , the concentration field $C_{u}$ satisfies the axisymmetric Laplace equation and the solution may be expressed naturally in cylindrical coordinates ( $\unicode[STIX]{x1D70C}=(x_{1}^{2}+x_{2}^{2})^{1/2},x_{3}$ ). The solution is calculated as a superposition of three fields, each satisfying the Laplace equation but with different boundary and far-field asymptotic conditions that are summarized in table 2,

(4.5) $$\begin{eqnarray}\displaystyle C_{u} & = & \displaystyle \tilde{C_{u}}+(1+\unicode[STIX]{x1D705}x_{3}),\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle \tilde{C_{u}} & = & \displaystyle -\unicode[STIX]{x1D705}C_{1}+C_{2}.\end{eqnarray}$$

The term $C_{1}$ has a well-known closed-form solution (Gladwell, Barber & Olesiak Reference Gladwell, Barber and Olesiak1983) which can be expressed in Legendre polynomials and functions of spheroidal coordinates ( $\unicode[STIX]{x1D70C}=\cosh \unicode[STIX]{x1D709}\sin \unicode[STIX]{x1D702},x_{3}=\sinh \unicode[STIX]{x1D709}\cos \unicode[STIX]{x1D6FD}$ ),

(4.7) $$\begin{eqnarray}C_{1}=\frac{2}{\unicode[STIX]{x03C0}}Q_{1}(\text{i}\sinh \unicode[STIX]{x1D709})P_{1}(\cos \unicode[STIX]{x1D702}).\end{eqnarray}$$

On the heterogeneous surface, this simplifies to

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle C_{1}=\left\{\begin{array}{@{}ll@{}}-\displaystyle \frac{2}{\unicode[STIX]{x03C0}}(1-\unicode[STIX]{x1D70C}^{2})^{1/2}\quad & \text{at }x_{3}=0,\unicode[STIX]{x1D70C}<1,\\ 0\quad & \text{at }x_{3}=0,\unicode[STIX]{x1D70C}\geqslant 1,\end{array}\right. & \displaystyle\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}C_{1}}{\unicode[STIX]{x2202}x_{3}}}=\left\{\begin{array}{@{}ll@{}}0\quad & \quad \text{at }x_{3}=0,\unicode[STIX]{x1D70C}<1,\\ \displaystyle R(\unicode[STIX]{x1D70C})=\frac{2}{\unicode[STIX]{x03C0}}\left(\arcsin \left(\frac{1}{\unicode[STIX]{x1D70C}}\right)-\frac{1}{\unicode[STIX]{x1D70C}}\left(1-\frac{1}{\unicode[STIX]{x1D70C}^{2}}\right)^{-1/2}\right)\quad & \quad \text{at }x_{3}=0,\unicode[STIX]{x1D70C}\geqslant 1.\end{array}\right. & \displaystyle\end{eqnarray}$$

The above forms of $C_{1}$ and its normal gradient $\unicode[STIX]{x2202}C_{1}/\unicode[STIX]{x2202}x_{3}$ will be useful in developing correlations and understanding the structure of the concentration field (when $\unicode[STIX]{x1D705}\gg 1$ ) in the ensuing sections. The solution form for $C_{2}$ involves a tedious spectral expansion and is not presented here. The solution technique to obtain both $C_{1}$ and $C_{2}$ is outlined in detail in Gladwell et al. (Reference Gladwell, Barber and Olesiak1983). The total concentration profile as a function of the radial position from the centre of the inert disk is shown in figure 4 for the case of $\unicode[STIX]{x1D6FF}=0.05$ . One can see that the far-field surface concentration quickly approaches a constant value $C_{p}\approx \unicode[STIX]{x1D719}S_{p}/\unicode[STIX]{x1D705}$ for $\unicode[STIX]{x1D705}>1$ . The field inside the inert disk is seen to grow large at large $\unicode[STIX]{x1D705}$ while the far-field solution decays to zero. This is consistent with the fact that the average surface concentration is unity. The concentration profile as predicted by (4.2) matches the profile calculated using BEM simulations, especially at $\unicode[STIX]{x1D705}\gg 1$ . At $\unicode[STIX]{x1D705}=0.1$ , the relative error in the concentration profile is approximately $0.05$ , which is in fact $O(\unicode[STIX]{x1D6FF})$ . This validates the non-interacting periodic approximation (4.2) over a wide range of values of the Damköhler number, and confirms our error analysis of table 1. It is important to note that the Sherwood number $S_{p}$ used in figure 4 is the one obtained from the converged BEM simulations.

Table 2.

Summary of boundary conditions for various field variables.

Figure 4.

The surface concentration profiles as calculated using the ‘non-interacting periodic’ approximation (4.2) (dashed) and solving the BEM equations (2.9) (thin, black) are plotted for various values of $\unicode[STIX]{x1D705}$ and $\unicode[STIX]{x1D719}=0.95$ , i.e. $\unicode[STIX]{x1D6FF}=0.05$ , and $P=0.001$ . The Sherwood number in the multiplying constant $\unicode[STIX]{x1D719}S_{p}/\unicode[STIX]{x1D705}$ of equation (4.2) is calculated using the BEM simulations. Indeed, the non-interacting periodic approximation works very well as the Damköhler number becomes large. $\unicode[STIX]{x1D705}=0.1$ (a), 1 (b), 10 (c), 100 (d), 300 (e), 500 (f).

4.2 Boundary layer in the surface concentration field around the inert patch

As claimed before, (3.8a ) is not universally correct as there exists (see figure 4) a surface concentration boundary layer surrounding the inert disk, due to which the last integral in (3.4) can no longer be assumed to be small. We plot the surface concentration field, $C_{p}=\unicode[STIX]{x1D719}S_{p}C_{u}/\unicode[STIX]{x1D705}$ , as a function of the distance from the inert disk boundary, $\unicode[STIX]{x1D70C}-1$ , in figure 5 for large values of the Damköhler number. Again, it should be noted that here $S_{p}$ is obtained from BEM simulations. When $C_{p}$ is scaled by $\unicode[STIX]{x1D705}^{1/2}$ and the distance from the disk boundary is scaled by $\unicode[STIX]{x1D705}$ , all of the curves collapse to a single curve close to the disk boundary, as seen in figure 5. This proves that the boundary layer thickness scales with $1/\unicode[STIX]{x1D705}$ while the surface concentration field scales with $\unicode[STIX]{x1D705}^{-1/2}$ . Since the normal gradient is related to the surface concentration via (2.2b ), the local flux just outside the inert disk scales with $\unicode[STIX]{x1D705}^{1/2}$ . This deserves an in-depth discussion, for which we re-examine (4.5) and (4.6) along with the individual boundary conditions in table 2.

It should be noted that since $C_{1}=0$ on the reactive region (from (4.8)), the surface concentration field is given by $C_{u}=1+C_{2}$ . Thus, when $\unicode[STIX]{x1D705}\ll 1$ , we see that the surface concentration field is, as expected, $C_{u}\approx 1$ . When $\unicode[STIX]{x1D705}\gg 1$ , it should be noted that $C_{2}$ approximately satisfies the boundary condition $-C_{2}=\unicode[STIX]{x2202}C_{1}/\unicode[STIX]{x2202}x_{3}$ with an error of $O(1/\unicode[STIX]{x1D705})$ . Using the known solution (4.9), we obtain the following asymptotic expansion near the edge of the inert disk: $C_{u}\approx (2/\unicode[STIX]{x03C0})(\unicode[STIX]{x1D70C}-1)^{-1/2}$ at $\unicode[STIX]{x1D70C}\rightarrow 1^{+}$ . Consequently, from the non-interacting approximation, it follows that $C_{p}\sim (1/\unicode[STIX]{x1D705})(\unicode[STIX]{x1D70C}-1)^{-1/2}$ , where we used (1.10c ) to substitute $\unicode[STIX]{x1D719}S_{p}/\unicode[STIX]{x1D705}\sim 1/\unicode[STIX]{x1D705}$ . Equivalently, from the reaction condition (2.2b ), the local surface flux at $\unicode[STIX]{x1D705}\gg 1$ is given by $\unicode[STIX]{x2202}C_{p}/\unicode[STIX]{x2202}x_{3}\sim (\unicode[STIX]{x1D70C}-1)^{-1/2}$ . This form is reminiscent of the well-known singularity (Stone Reference Stone1989) in the surface flux at the disk boundary for an isothermal disk in an insulating surface, which is the type ( ${\mathcal{R}}$ ) counterpart to this problem. It follows immediately from the asymptotic forms of $C_{p}$ and $\unicode[STIX]{x2202}C_{p}/\unicode[STIX]{x2202}x_{3}$ that if the boundary layer thickness $\unicode[STIX]{x1D70C}-1$ scales with $1/\unicode[STIX]{x1D705}$ , the surface concentration field $C_{p}$ near the disk edge must scale with $\unicode[STIX]{x1D705}^{-1/2}$ . Indeed, the collapsed curves in figure 5 display a slope of $-1/2$ , corresponding to the $1/(\unicode[STIX]{x1D70C}-1)^{1/2}$ decay.

The boundary layer structure provides insight as to when the non-interacting approximation (4.2) breaks down. As the boundary layer thickness scales with $1/\unicode[STIX]{x1D705}$ , when $\unicode[STIX]{x1D705}\ll \unicode[STIX]{x1D6FF}^{1/2}$ , the boundary layer must scale with a length larger than the size of the periodic box. In such a case, the ‘inner’ surface concentration scales on the size of the periodic box and is not screened from the periodic influence of the domain. Thus, (4.2) then leads to errors of size $\unicode[STIX]{x1D6FF}$ , as noted before in figure 4(a).

Figure 5.

A rescaled surface concentration field is plotted versus the distance from the disk boundary for $P=0.001$ . Boundary layer existence and scaling is apparent from the fact that all of the curves approach a constant curve in the limit as $\unicode[STIX]{x1D705}\rightarrow \infty$ .

5.1 Case I: no flow, $P=0$

The non-interacting periodic inert disk approximation leads to an alternative more accurate approximation for the Sherwood number in type ( ${\mathcal{I}}_{p}$ ) surfaces (periodic inert disks). The key idea is to recognize that $\int _{{\mathcal{W}}_{0}}C_{p}(\boldsymbol{x}_{b})\,\text{d}A(\boldsymbol{x}_{b})+\int _{{\mathcal{P}}_{0}}C_{p}(\boldsymbol{x}_{b})\,\text{d}A(\boldsymbol{x}_{b})=B_{0}$ , which is a consequence of the average surface concentration being unity (see (1.5)). The integral over the reactive surface ${\mathcal{P}}$ is simply $P_{0}S_{p}/\unicode[STIX]{x1D705}$ (see (1.4)). Using the non-interacting periodic approximation and the relations $P_{0}=W_{0}\unicode[STIX]{x1D719}/\unicode[STIX]{x1D6FF},B_{0}=W_{0}/\unicode[STIX]{x1D6FF}$ , we can write $(\unicode[STIX]{x1D719}S_{p}/\unicode[STIX]{x1D705}W_{0})\int _{{\mathcal{W}}_{0}}C_{u}(\boldsymbol{x}_{b})\,\text{d}A(\boldsymbol{x}_{b})=1/\unicode[STIX]{x1D6FF}-\unicode[STIX]{x1D719}/\unicode[STIX]{x1D6FF}(S_{p}/\unicode[STIX]{x1D705})$ . This simplifies to

(5.1) $$\begin{eqnarray}S_{p}(\unicode[STIX]{x1D705},P=0,\unicode[STIX]{x1D719})=\frac{\unicode[STIX]{x1D705}}{\unicode[STIX]{x1D719}+\unicode[STIX]{x1D719}\unicode[STIX]{x1D6FF}\displaystyle \left[\frac{1}{W_{0}}\int _{{\mathcal{W}}_{0}}C_{u}(\boldsymbol{x}_{b})\,\text{d}A(\boldsymbol{x}_{b})\right]}.\end{eqnarray}$$

The surface integral above can be written as a sum of surface integrals of $C_{1}$ and $C_{2}$ (4.5)–(4.6). The surface integral of $C_{2}$ can be represented as a summation involving the spectral coefficients (Gladwell et al. Reference Gladwell, Barber and Olesiak1983), while $C_{1}$ can be explicitly integrated using the closed-form expression (4.8). The final form of the surface integral of $C_{u}$ in (5.1) is given by

(5.2) $$\begin{eqnarray}\frac{1}{W_{0}}\int _{{\mathcal{W}}_{0}}C_{u}(\boldsymbol{x}_{b})\,\text{d}A(\boldsymbol{x}_{b})=1+\frac{4\unicode[STIX]{x1D705}}{3\unicode[STIX]{x03C0}}+\text{corr}(\unicode[STIX]{x1D705}),\end{eqnarray}$$

where the integral involving $C_{2}$ contributes a ‘correction’ term $\text{corr}(\unicode[STIX]{x1D705})$ . We fit the following Padé approximant (figure 14, appendix B) using least squares with less than 1 % error:

(5.3) $$\begin{eqnarray}\frac{\text{corr}(\unicode[STIX]{x1D705})}{\unicode[STIX]{x1D705}}=\frac{4}{3\unicode[STIX]{x03C0}}\left(\frac{0.8837}{\unicode[STIX]{x1D705}^{0.6809}+0.8837}\right).\end{eqnarray}$$

The prefactor $4/3\unicode[STIX]{x03C0}$ has been taken out to emphasize that $\text{corr}(\unicode[STIX]{x1D705})/\unicode[STIX]{x1D705}\approx 4/3\unicode[STIX]{x03C0}$ at $\unicode[STIX]{x1D705}\ll 1$ . Using equations (5.1), (5.2) and (5.3), one can write the following improved approximation for the Sherwood number:

(5.4a ) $$\begin{eqnarray}\displaystyle S_{p}(\unicode[STIX]{x1D705},P=0,\unicode[STIX]{x1D719}) & = & \displaystyle \frac{\unicode[STIX]{x1D705}}{\unicode[STIX]{x1D719}+\unicode[STIX]{x1D719}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D705}\displaystyle \left[\frac{1}{\unicode[STIX]{x1D705}}+\frac{4}{3\unicode[STIX]{x03C0}}+\frac{\text{corr}(\unicode[STIX]{x1D705})}{\unicode[STIX]{x1D705}}\right]},\end{eqnarray}$$

(5.4b ) $$\begin{eqnarray}\displaystyle & = & \displaystyle \frac{\unicode[STIX]{x1D705}}{\unicode[STIX]{x1D719}+\unicode[STIX]{x1D719}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D705}\displaystyle \left[\frac{1}{\unicode[STIX]{x1D705}}+\frac{4}{3\unicode[STIX]{x03C0}}\left(1+\frac{0.8837}{\unicode[STIX]{x1D705}^{0.6809}+0.8837}\right)\right]}.\end{eqnarray}$$

Figure 6.

The Sherwood number $S_{p}(\unicode[STIX]{x1D705},0,\unicode[STIX]{x1D719})$ is plotted as a function of the Damköhler number at various values of the area fraction $\unicode[STIX]{x1D719}$ . The approximation (5.4b ) (thin lines) matches very well the ‘exact’ values of the Sherwood number obtained from BEM simulations (filled circles). The perturbation approximation (3.8a ) (dotted) is valid only when $\unicode[STIX]{x1D705}\ll 1$ .

The Sherwood number as predicted by (5.4b ) is plotted as a function of the Damköhler number in figure 6 for various values of $\unicode[STIX]{x1D719}$ . Certainly, the approximation (5.4b ) is very good in predicting the ‘exact’ Sherwood number obtained from BEM simulations. As expected, the approximation is better as we approach $\unicode[STIX]{x1D719}\rightarrow 1$ .

Figure 7.

The inverse of the Sherwood number is plotted at large $\unicode[STIX]{x1D705}$ ( $\unicode[STIX]{x1D705}=300$ ) to demonstrate the scaling with $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}$ as observed in (5.5). For both inert disks (filled circles) and squares (filled boxes), the same scaling is observed.

The form of (5.4b ) warrants a discussion. In particular, at $\unicode[STIX]{x1D705}=\infty$ , the term corresponding to the correction in the denominator is zero and we are left with an expression that is purely a function of the area fraction of disks,

(5.5) $$\begin{eqnarray}S_{p}(\unicode[STIX]{x1D705}=\infty ,P=0,\unicode[STIX]{x1D719})=\frac{3\unicode[STIX]{x03C0}}{4\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}}.\end{eqnarray}$$

Thus, equation (5.5) confirms the purely geometric nature of the problem in the diffusion-limited regime of type ( ${\mathcal{I}}_{p}$ ) surfaces. The prefactor $3\unicode[STIX]{x03C0}/4$ is a ‘shape factor’ corresponding to the disk-shaped inert patches. Indeed, for the case of square inert patches, we find evidence for a different prefactor in figure 7. Here, results for both disk- and square-shaped patches are shown at $\unicode[STIX]{x1D705}\gg 1$ . Intriguingly, this prefactor appears to be weakly shape-dependent. Sadhal & Tio (Reference Sadhal and Tio1991) showed that the contact resistance for the case of heat transfer between conducting solids with periodic circular adiabatic gaps (corresponding directly to the case of $P=0,\unicode[STIX]{x1D705}=\infty$ ) grows as $1/\unicode[STIX]{x1D6FF}$ to the leading order at $\unicode[STIX]{x1D6FF}\ll 1$ , which is a confirmation of (5.5).

Figure 8.

The term $1/(\unicode[STIX]{x1D719}S_{p})-1/\unicode[STIX]{x1D705}$ scales linearly with $1/(1-\unicode[STIX]{x1D719})$ at all values of the Damköhler number, as suggested by (5.7) and proved by the Sherwood number obtained from BEM simulations. The perturbation approximation (3.8a ) (dotted) is valid only when $\unicode[STIX]{x1D705}\ll 1$ . The Péclet number is set to zero here.

To demonstrate the importance of the boundary layer at $\unicode[STIX]{x1D705}\gg 1$ , observe that for the corresponding perturbation series approximation (3.8b ), the shape factor for the disk is $3\unicode[STIX]{x03C0}/8$ , which is half of the correct shape factor. Knowing the boundary layer structure from § 4.2, when $\unicode[STIX]{x1D705}\gg 1$ , we can write $C_{p}\approx (\unicode[STIX]{x1D719}S_{p}/\unicode[STIX]{x1D705})C_{u}\approx (\unicode[STIX]{x1D719}S_{p}/\unicode[STIX]{x1D705})(1-\unicode[STIX]{x2202}C_{1}/\unicode[STIX]{x2202}x_{3})$ . Matching this with the perturbation expansion (3.2), we get $C_{p}^{(1)}=-(\unicode[STIX]{x1D719}S_{p}/\unicode[STIX]{x1D705})(\unicode[STIX]{x2202}C_{1}/\unicode[STIX]{x2202}x_{3})$ . Using this relation for $C_{p}^{(1)}$ along with (4.9) and the unbounded Green’s function for $P=0$ , we can now analytically calculate the term neglected in the perturbation approximation in (3.4) (details are given in appendix D),

(5.6) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D705}\int _{{\mathcal{W}}_{0}}\int _{{\mathcal{P}}_{0}}G_{p}(\boldsymbol{y}_{b}|\boldsymbol{x}_{b};0)C_{p}^{(1)}(\boldsymbol{y}_{b})\,\text{d}A(\boldsymbol{y}_{b})\,\text{d}A(\boldsymbol{x}_{b})\nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D719}S_{p}\int _{{\mathcal{W}}_{0}}\int _{{\mathcal{P}}_{0}}\frac{-1}{2\unicode[STIX]{x03C0}\Vert \boldsymbol{x}_{b}-\boldsymbol{y}_{b}\Vert }{\displaystyle \frac{\unicode[STIX]{x2202}C_{1}(\boldsymbol{y}_{b})}{\unicode[STIX]{x2202}y_{3}}}\,\text{d}A(\boldsymbol{y}_{b})\,\text{d}A(\boldsymbol{x}_{b})\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{4}{3}\unicode[STIX]{x1D719}S_{p}.\end{eqnarray}$$

An immediate consequence of (5.6) is that (3.4) now reads as $P_{0}S_{p}/\unicode[STIX]{x1D705}=P_{0}-\unicode[STIX]{x1D719}S_{p}f(0)+4\unicode[STIX]{x1D719}S_{p}/3$ at $P=0$ . As $f(0)=8/3$ (see (3.8c )), the corrected version of (3.8b ) for circular disks is given by $1/S_{p}(\infty ,0,\unicode[STIX]{x1D719})=(4/3)\unicode[STIX]{x1D6FF}/\unicode[STIX]{x03C0}$ , which now has the correct shape factor. Thus, neglect of the last term in (3.4) leads to drastically different predicted values of the Sherwood number in figure 6 at $\unicode[STIX]{x1D705}\gg 1$ . For a different geometry of the inert patches, one can expect the Sherwood number $S_{p}(\unicode[STIX]{x1D705}=\infty ,P=0,\unicode[STIX]{x1D719})$ to scale with $1/(\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719})$ with a different shape factor due to the geometric nature of the limit $\unicode[STIX]{x1D705}\rightarrow \infty$ .

Figure 9.

The Sherwood number $S_{p}(\unicode[STIX]{x1D705},P,\unicode[STIX]{x1D719})$ is plotted as a function of the Damköhler number at various values of the area fraction $\unicode[STIX]{x1D719}$ for the cases (a) $P=0.001$ , (b) $P=2$ and (c) $P=20$ . The approximation (5.8) (thin lines) matches very well the ‘exact’ values of the Sherwood number obtained from BEM simulations (filled circles). The $P=0$ approximation (dashed lines, (5.4b )) is inaccurate at $P=20$ .

Casting (5.4b ) as

(5.7) $$\begin{eqnarray}\frac{1}{\unicode[STIX]{x1D719}S_{p}(\unicode[STIX]{x1D705},0,\unicode[STIX]{x1D719})}=\frac{1}{\unicode[STIX]{x1D705}}+(1-\unicode[STIX]{x1D719})\left[\frac{1}{\unicode[STIX]{x1D705}}+\frac{4}{3\unicode[STIX]{x03C0}}\left(1+\frac{0.8837}{\unicode[STIX]{x1D705}^{0.6809}+0.8837}\right)\right],\end{eqnarray}$$

we see that any approximation based on the method of resistances (for example, (3.8a )) fails as it is very different from (5.7). First, there is an additional multiplying factor $\unicode[STIX]{x1D719}$ for the Sherwood number on the left-hand side of (5.7). Second, the last term on the right-hand side now depends on the Damköhler number. The scaling of this last term is evident from the plot of $1/(\unicode[STIX]{x1D719}S_{p})-1/\unicode[STIX]{x1D705}$ versus $1/(1-\unicode[STIX]{x1D719})$ at various values of the Damköhler number in figure 8. Specifically, we see that the Sherwood number in (5.7), which is based on the non-interacting approximation (4.2), closely matches the values predicted by the BEM simulations at $\unicode[STIX]{x1D705}\gg 1$ and $\unicode[STIX]{x1D6FF}\ll 1$ . At small Damköhler number and $\unicode[STIX]{x1D6FF}>0.1$ , (5.7) deviates from the BEM data because of the errors of size $O(\unicode[STIX]{x1D6FF})$ or $O(\unicode[STIX]{x1D6FF}^{1/2}\unicode[STIX]{x1D705})$ in the non-interacting approximation (4.2) (see table 1). For the range of $\unicode[STIX]{x1D6FF}$ shown in figure 8, the maximum error in (5.7) occurs at $\unicode[STIX]{x1D705}=O(1)$ . Nonetheless, the relative error in the Sherwood number in (5.4b ) seen from the data in figure 6 is less than 2 %.

5.2 Case II: shear flow, $P>0$

To summarize, we first obtained a perturbation series based approximation for the Sherwood number, which resembles the method of resistances form in § 3. We then showed the existence of a surface concentration boundary layer (at $\unicode[STIX]{x1D705}\gg 1$ ) using the non-interacting periodic disk approach in § 4, which led to an improved correlation (5.4b ) that is valid only when $P=0$ . We observe that when the Damköhler number is small ( $\unicode[STIX]{x1D705}<1$ ), the correction term is $\text{corr}(\unicode[STIX]{x1D705})/\unicode[STIX]{x1D705}\approx 4/(3\unicode[STIX]{x03C0})$ . Thus, in this regime of the Damköhler number and at vanishingly small area fraction of disks, (5.7) is similar in form to (3.8a ) with $f(0)/W_{0}=8/(3\unicode[STIX]{x03C0})$ (see (3.8c )). Taking inspiration from this, to capture the effect of flow in our correlations, we bridge the approximation (3.8c ) for $f(P)$ in (3.8a ) with the last term in (5.7) and hypothesize the following universal approximation:

(5.8) $$\begin{eqnarray}S_{p}(\unicode[STIX]{x1D705},P,\unicode[STIX]{x1D719})=\frac{\unicode[STIX]{x1D705}}{\unicode[STIX]{x1D719}+\displaystyle \unicode[STIX]{x1D719}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D705}\left[\frac{1}{\unicode[STIX]{x1D705}}+\left(\left(\frac{4}{3\unicode[STIX]{x03C0}}+\frac{\text{corr}(\unicode[STIX]{x1D705})}{\unicode[STIX]{x1D705}}\right)^{-15/8}+\left(\frac{4.2075}{\unicode[STIX]{x03C0}}P^{-1/3}\right)^{-15/8}\right)^{-8/15}\right]},\end{eqnarray}$$

where $\text{corr}(\unicode[STIX]{x1D705})/\unicode[STIX]{x1D705}$ is given by (5.3). At $P=0$ , we recover (5.4b ), while the above correlation also matches the perturbation approximation when $\unicode[STIX]{x1D705}\ll 1$ by virtue of its construction. To understand the validity of this approximation, we plot the equivalent of figures 6 and 8 for the Péclet numbers $P=0.001,2$ and $20$ in figures 9 and 10 respectively. We again observe that the predictions of the new correlation are close to the values of the Sherwood number obtained from BEM simulations. The various correlations valid for the case of $\unicode[STIX]{x1D6FF}\ll 1$ are summarized in table 3. Other than having the correct behaviour in limiting cases, we also have a mathematical explanation for why (5.8) is a good match to BEM results for such a wide range of Péclet numbers. Noting that our error analysis of the non-interacting approximation (4.2) is valid for any Péclet number, we hypothesize that the area-averaging procedure leading to equation (5.1) must hold true and any errors observed must be introduced in the evaluation of the integral term $\int _{{\mathcal{P}}_{0}}C_{p}(\boldsymbol{x}_{b})\,\text{d}A(\boldsymbol{x}_{b})$ . These errors stem from leaving out the accretion wake that extends from the inert disk to the outside due to shear flow. These errors will be small as long as the accretion wake does not spread over a large area compared with the reactive region ${\mathcal{P}}_{0}$ .

Figure 10.

The term $1/(\unicode[STIX]{x1D719}S_{p})-1/\unicode[STIX]{x1D705}$ scales linearly with $1/(1-\unicode[STIX]{x1D719})$ at all values of the Damköhler number over a wide range of the Péclet number ( $P\in \{0.001,2,20\}$ ) (a–c), as suggested by (5.7) and (5.8) and proved by the Sherwood number obtained from BEM simulations.

This is indeed confirmed by the contour plots in figure 11 for Péclet numbers ranging from $0.001$ to $20$ and the corresponding data points in figure 9. Referring to figure 11(a,c,e), as the Péclet number increases, the accretion wake shifts downstream and becomes flatter, and its thickness scales with the boundary layer thickness. This is further shown by figure 11(b,d,f), where the wake stretches downstream of the inert patch with increasing Péclet number. This stretching and flattening of the wake region leads to increased Sherwood number. Even at $P=20$ , the wake size is comparable to the size of the inert region. Thus, the correlation (5.8) can be expected to be valid at large Péclet numbers as long as the inert disks are sufficiently dilute such that the wakes are not strongly interacting. Equation (5.8) also captures the fact that as the Péclet number is increased, the boundary layer thickness becomes the length scale for transport at the surface and the effect of shear must dominate the effect of reaction. However, at small Damköhler numbers, i.e. in the reaction-limited regime, the $1/\unicode[STIX]{x1D705}$ term in the denominator of (5.8) dominates and, as seen in figure 9, the Sherwood number is weakly dependent on the Péclet number. Thus, the rich interaction between the surface inert irregularities, reaction and flow is surprisingly well captured by (5.8).

Figure 11.

Contour plots of the concentration field $C_{p}-(1+\unicode[STIX]{x1D719}S_{p}x_{3})$ (a,c,e) in the vertical plane ( $x_{1}$ – $x_{3}$ ) and $C_{p}$ (b,d,f) on the surface ( $x_{1}$ – $x_{2}$ ) for a Damköhler number of $\unicode[STIX]{x1D705}=10$ and $\unicode[STIX]{x1D719}=0.95$ . Here, the boundary layer and wake structure are exposed. Panels (a,b) correspond to $P=0.001$ , (c,d) to $P=2$ , and (e,f) to $P=20$ .

It should be noted that for inert regions of different shapes, one would develop similar correlations in the spirit of (5.8). For instance, we expect the same asymptotic scaling of $f(P)/W_{0}$ with $P$ at $P\gg 1$ , and that of $S_{p}$ with $\unicode[STIX]{x1D6FF}$ at $\unicode[STIX]{x1D705}\gg 1$ . However, we expect that the constants in (3.8c ) and (5.3) will be shape-dependent. In our analysis, the inert patches are always aligned in the shear direction. While the Péclet number is small, the wakes between inert patches do not affect each other and the non-interacting approximation holds. This forms the central argument for the validity of (5.8). Thus, wake sharing (at $P\gg 1$ ) will be direction-dependent and thus influence the Sherwood number except in the limit $\unicode[STIX]{x1D6FF}\ll 1$ . This is supported by the previous study of Shah & Shaqfeh (Reference Shah and Shaqfeh2015), where the Sherwood number for disordered media was found to have the same scaling when $\unicode[STIX]{x1D719},\unicode[STIX]{x1D6FF}\sim O(0.1)$ .

Table 3.

Summary of correlations and dependence on various parameters. Equation (5.8) is the most general correlation, which can be used to obtain all others as special cases. It should be noted that in equation (5.8), the term involving $\text{corr}(\unicode[STIX]{x1D705})$ is given by (5.3).