On the boundary layer structure near a highly permeable porous interface

The method of matched asymptotic expansions is used to study the canonical problem of steady laminar flow through a narrow two-dimensional channel blocked by a tight-fitting finite-length highly permeable porous obstacle. We investigate the behaviour of the local flow close to the interface between the single-phase and porous regions (governed by the incompressible Navier--Stokes and Darcy flow equations, respectively). We solve for the flow in these inner regions in the limits of low and high Reynolds number, facilitating an understanding of the nature of the transition from Poiseuille to plug to Poiseuille flow in each of these limits. Significant analytic progress is made in the high-Reynolds-number limit, and we explore in detail the rich boundary layer structure that occurs. We consider the three-dimensional generalization to unsteady laminar flow through and around a tight-fitting highly permeable cylindrical porous obstacle within a Hele-Shaw cell. For the high-Reynolds-number limit, we give the coupling conditions and interfacial stress in terms of the outer flow variables, allowing information from a nonlinear three-dimensional problem to be obtained by solving a linear two-dimensional problem. Finally, we illustrate the utility of our analysis by considering the specific example of time-dependent forced far-field flow in a Hele-Shaw cell containing a porous cylinder with a circular cross-section. We determine the internal stress within the porous obstacle, which is key for tissue engineering applications, and the interfacial stress on the boundary of the porous obstacle, which has applications to biofilm erosion. In the high-Reynolds-number limit, we demonstrate that the fluid inertia can result in the cylinder experiencing a time-independent net force, even when the far-field forcing is periodic with zero mean.

The method of matched asymptotic expansions is used to study the canonical problem of steady laminar flow through a narrow two-dimensional channel blocked by a tight-fitting finite-length highly permeable porous obstacle. We investigate the behaviour of the local flow close to the interface between the single-phase and porous regions (governed by the incompressible Navier-Stokes and Darcy flow equations, respectively). We solve for the flow in these inner regions in the limits of low and high Reynolds number, facilitating an understanding of the nature of the transition from Poiseuille to plug to Poiseuille flow in each of these limits. Significant analytical progress is made in the high Reynolds number limit, and we explore in detail the rich boundary layer structure that occurs. We derive general results for the interfacial stress and for the conditions that couple the flow in the outer regions away from the interface. We consider the three-dimensional generalization to unsteady laminar flow through and around a tight-fitting highly permeable cylindrical porous obstacle within a Hele-Shaw cell. For the high Reynolds number limit, we give the coupling conditions and interfacial stress in terms of the outer flow variables, allowing information from a nonlinear three-dimensional problem to be obtained by solving a linear two-dimensional problem. Finally, we illustrate the utility of our analysis by considering the specific example of time-dependent forced far-field flow in a Hele-Shaw cell containing a porous cylinder with a circular cross-section. We determine the internal stress within the porous obstacle, which is key for tissue engineering applications, and the interfacial stress on the boundary of the porous obstacle, which has applications to biofilm erosion. In the high Reynolds number limit, we demonstrate that the fluid inertia can result in the cylinder experiencing a time-independent net force, even when the far-field forcing is periodic with zero mean.

Introduction
Flow through a narrow channel containing a porous blockage is a canonical problem with numerous practical applications. For example, this flow configuration is used: to purify colloids in separation science, where it is also known as dead-end filtration (McCarthy et al. 1998; Van der Bruggen et al. 2003;Bessiere et al. 2005); to deliver nutrient to cells growing in a porous tissue construct in tissue engineering (El Haj et al. 1990;Jaasma, Plunkett & O'Brien 2008;O'Dea, Waters & Byrne 2009); and to erode porous biofilms that have grown within a pipe, where the erosion is dependent on the interfacial shear stress (Picioreanu, Van Loosdrecht & Heijnen 2001;Telgmann, Horn & Morgenroth 2004;Duddu, Chopp & Moran 2009). In many of these problems, the flow near the interface between the fluid and the porous blockage is of experimental or industrial interest. For example, the stress acting on the interface may be important because it is often coupled to some interesting physical phenomenon, such as erosion, mechanotransduction or movement of the entire porous blockage. However, in many models the behaviour of the flow near the interface is ignored, leading to incomplete information about the flow. As the lubrication equations used to approximate the flow away from the interface are reduced in order, approximated interfacial boundary conditions must be applied to couple the flows away from the interface (Waters et al. 2006;Cummings & Waters 2007;Cummings et al. 2009;O'Dea et al. 2009O'Dea et al. , 2013. A detailed understanding of the flow behaviour close to the interface will enable the correct coupling conditions to be applied and (dynamic) interfacial effects to be accurately included. It is, therefore, of fundamental interest to understand the flow near the interface of a porous blockage.
Our main motivation arises from a tissue engineering problem. One approach to in vitro tissue engineering is to seed cells onto a porous biomaterial scaffold, which is then cultured within a bioreactor. The combination of cells and scaffold is referred to as a tissue construct and one particular aim of tissue engineering is to make porous tissue constructs as permeable as possible whilst maintaining their structural integrity. This is to enhance nutrient delivery to the cells residing in the interior of the construct via advection. Moreover, many cells are mechanosensitive, and their proliferation rate depends on the stress that they experience. In a high aspect ratio vessel bioreactor, the porous tissue construct is placed within a bioreactor shaped like a Petri dish turned on its side, and the entire set-up is saturated with a nutrient-rich fluid, and rotated around the bioreactor axis (figure 1). The construct moves according to the force it experiences and, in the short term, the construct undergoes periodic motion. However, over a long time, the tissue construct drifts from its periodic orbit, an effect that can be attributed to weak inertia. To predict the construct trajectory and, ultimately, the nutrient transport in such a bioreactor, we must determine the forces acting on the construct. Thus, determining the flow through the porous construct is of particular importance. To gain insights into this problem, we consider the flow past a tissue construct held at a fixed position, with a particular view to investigating the effect of weak inertia. Such an analysis can then be used to determine the construct trajectory when it is free to move. Another application of interest is the interfacial erosion of a porous biofilm, which is proportional to the square root of the fluid shear stress evaluated at the interface (Duddu et al. 2009). Thus, it is of interest to determine both the internal stress within the porous obstacle for tissue engineering applications and the interfacial stress acting on the obstacle boundary for biofilm erosion problems.
The main aim of this paper is to investigate steady laminar flow through a narrow two-dimensional channel blocked by a tightly fitting finite-length porous obstacle, and then to generalise this analysis to investigate a time-dependent three-dimensional flow past a porous cylinder with an arbitrary smooth cross-section within a Hele-Shaw cell, as illustrated in figure 2. In particular, we are interested in the behaviour of the flow near the interface between the single-phase and porous regions, governed by the Navier-Stokes and Darcy equations, respectively. Transition from plug to Poiseuille flow in a two-dimensional channel is a classic problem, with the asymptotic structure examined by Van Dyke (1970) and Wilson (1971), and numerical solutions given by, for example, Bodoia & Osterle (1961) and Brandt & Gillis (1966). Although Van Dyke (1970) mentions that a porous mesh would have to be used experimentally to induce uniform flow at the entry to a channel, the nature of the coupled flow between single-phase and porous regions remains an open question -one which we will answer in this paper. Thompson (1968) used the method of matched asymptotic expansions to analyse the three-dimensional flow near a solid circular cylinder within a Hele-Shaw cell, exploiting the small aspect ratio of the cell. Due to the no-flux condition on the solid cylinder, the normal flow near the interface is small. We are motivated by the tissue engineering application in which there is flow within a highly permeable porous obstruction of arbitrary smooth cross-section. The normal flow, in general, is then significant near the interface and the flows inside and outside the obstacle are coupled through their boundary conditions at the interface. Although we are interested in a highly permeable obstacle, we do not consider the Brinkman extension to the Darcy equations in this paper because the porous construct must maintain its structural integrity. That is, we consider a porous medium whose underlying solid matrix is connected, so that Brinkman's equations do not apply (Auriault 2009).
We shall make extensive use of the method of matched asymptotic expansions by exploiting the small aspect ratio of the channel/Hele-Shaw cell -one of our goals is to understand the flow near the interface (in an 'inner' region) given a general flow away from the interface (in an 'outer' region). Although the outer problems satisfy reduced equations, the inner problems are quasi-two-dimensional in each plane perpendicular to the interface. In particular, we determine the behaviour of inherently local properties, such as the stress acting on the interface. Additionally, we derive systematically the conditions that couple the outer equations for a general far-field forcing and for a general smooth cross-section of the porous obstacle. Due to their generality, the results obtained in this paper can significantly reduce the computational expense required to solve flow problems with coupled single-phase and porous regions, whilst retaining important interfacial information.
There is a large literature addressing the question of interfacial conditions on the boundary of a porous obstacle (and it remains an area of active research); see, for example, Nield & Bejan (2006) for a comprehensive review. We do not address this question here, and use as our interfacial boundary conditions: continuity of flux and continuity of pressure, together with a no tangential slip condition. The first two are derived in Levy & Sanchez-Palencia (1975), the last is a special case of the general tangential slip condition derived in Carraro et al. (2015).
The structure of this paper is as follows. In § 2, we consider the steady twodimensional case of flow through a channel containing a tight-fitting highly permeable porous obstacle. In § 2.1, we formulate the mathematical problem. In § 2.3, we describe the asymptotic structure in the small aspect ratio limit and present the problems to be solved to couple the outer single-phase and porous regions. In § § 2.4 and 2.5, we investigate the asymptotic behaviour of these coupling conditions in the small and large Reynolds number limits, respectively. In particular, we find that the large Reynolds number limit induces a further boundary layer structure, and we investigate this in full. In § 3 we extend our analysis to an unsteady three-dimensional flow in a Hele-Shaw cell past a tight-fitting highly permeable cylindrical porous obstacle. If the curvature of the cross-sectional boundary of the porous obstacle is not large (in a sense to be made precise later), the boundary layer structure in each plane perpendicular to the interface is equivalent to the two-dimensional case considered in § 2, and we are able to generalize the two-dimensional results to the three-dimensional case. In § 4 we apply the general results of § 3 to a time-dependent forced far-field flow in a Hele-Shaw cell containing a porous cylinder with circular cross-section. We determine the internal and interfacial stress within and on the porous construct, and show that fluid inertia can result in the generation of a net force on the cylinder, even when the far-field flow is periodic with zero mean. If the cylinder were allowed to move, this would cause a long-term drift, thus highlighting the singular perturbative nature of fluid inertia, and how it can be important in dynamical problems arising in lubrication flow.

Model formulation
We consider the steady two-dimensional flow of an incompressible Newtonian fluid with density ρ and viscosity µ in an infinitely long channel of height h. The channel contains a fully saturated porous obstacle with rectangular cross-section, height h and length 2l, as illustrated in figure 2(a). The aspect ratio = h/l is assumed to be small. We use Cartesian coordinates (x, z) along and across the channel, as illustrated in figure 2(a), such that in the channel z ∈ [0, h] and the plug lies in x ∈ [−l, l], z ∈ [0, h]. The flow is driven by a prescribed upstream pressure gradient, such that in the far field we have Poiseuille flow with cross-sectionally averaged velocity of magnitude U av .
The flow exterior to the porous plug (the light regions in figure 2a), which we refer to as the 'exterior flow', is governed by the incompressible steady Navier-Stokes equations with no body force. The exterior fluid velocity, pressure, and stress tensor are denoted by u = ue x + we z , p, and σ , respectively, where e x and e z are the unit vectors in the x-and z-directions, respectively. The fluid inside the porous plug (the dark region in figure 2a), which we refer to as the 'interior flow', is governed by the incompressible Darcy equations, again with no body force. The Darcy (volumeaveraged) velocity and (interstitial) pressure within the porous plug are denoted by Q = Ue x + We z and P, respectively. The porous plug has constant permeability k. We do not use the porosity of the plug as a parameter as it is absorbed into the definition of the volume-averaged velocity, but note that it is assumed to be constant.
The dimensionless governing equations in the exterior region |x| > 1, 0 < z < 1 are the scaled Navier-Stokes equations given by where Re = ρU av h/µ is the modified Reynolds number based on the channel height. The governing equations in the porous plug |x| < 1, 0 < z < 1 are the scaled Darcy equations given by where the dimensionless permeability K = k/h 2 . For the Darcy equations to be valid, we must have K 1. We discuss the typical sizes of the three dimensionless parameters appearing in our model, , Re and K, in § 2.2.

Boundary conditions
On the channel walls we impose no-flux and no-slip conditions on the exterior flow and a no-flux condition on the Darcy flow, so that u = 0 on z = 0, 1 for |x| > 1; W = 0 on z = 0, 1 for |x| < 1. (2.3a,b) On the inflow (x = −1) and outflow (x = 1) interfaces, continuity of flux, continuity of pressure and a no tangential slip condition (Levy & Sanchez-Palencia 1975;Carraro et al. 2015) are given by The flow is driven by imposing a constant, upstream pressure gradient in the far field, with dp dx → −12 as x→ −∞,   (table 1). We will investigate the dimensionless problem defined by (2.1)-(2.5) in specific regions of this parameter space.
The dimensionless problem is characterised by two length scales: the obstacle halflength, 1, and the channel height, . We proceed by considering the case in which 1, the channel height then being much smaller than the half-length of the obstacle. The value of K must be small enough for the underlying solid matrix to be connected, resulting in Darcy's equations governing the flow within the porous medium (Auriault 2009), but large enough to enhance nutrient transfer via fluid advection to the cells within the plug. There is an absolute upper bound of K < 1/12, and this can be seen by considering unobstructed Poiseuille flow through a channel. In reality, K 1, and we work in this limit. We see in the three-dimensional case, considered in § 3, that K = O( ) corresponds to an obstacle which is impermeable to the flow at leading order. Hence, this limit reduces to the leading-order problem given in Thompson (1968), who considered the three-dimensional flow near a solid circular cylinder within a Hele-Shaw cell. We consider instead the limit in which K 1, corresponding to the regime in which the plug is long, thin and permeable at leading order in . Examples of long, thin tissues include skin (Lei et al. 2011), urothelial tissue (Gabouev et al. 2003 and the cornea (Su et al. 2003). Although we do assume that K 1 throughout this paper, we effectively treat K = O(1) as → 0 in our following asymptotic analysis for mathematical expediency.
We must also define the relative scaling of Re. There is a distinguished limit when Re = O(1), and the investigation of this limit is the subject of the next section. Further analytic progress can be made when Re is either small or large: we consider the asymptotic sublimit Re 1 in § 2.4, and the asymptotic sublimit 1 1/K Re 1/ in § 2.5. We note that Darcy's law (2.2) and the continuity of pressure condition (2.4b) would require modification if we were to generalize this work to significantly larger Reynolds numbers (Kaviany 2012). 2.3. Asymptotic structure for Re = O(1) We seek a solution using the method of matched asymptotic expansions with small parameter . The asymptotic structure is shown in figure 3(a), where each of seven asymptotic regions are labelled I-VII. We describe the regions as follows, noting that in each case the relevant transverse length scale is the channel height, . There are three 'outer' regions characterised by an O(1) axial length scale: two for the exterior fluid (regions I and VII) in which we have Poiseuille flow at leading order, and one for the interior fluid (region IV) in which we have plug flow at leading order. Thus, there is a transition from Poiseuille to plug to Poiseuille flow, between regions I, IV and VII, respectively, as illustrated schematically in figure 3(b). The inflow transition from Poiseuille to plug flow occurs near the interface x = −1, where there are two 'inner' regions, one on each side of the interface, each characterised by an axial length scale of O( ). The exterior inflow inner region is region II, and the porous inflow inner region is region III. In a similar manner, the outflow transition from plug to Poiseuille flow occurs near the interface x = 1, where again there are two inner regions, one on each side of the interface, and each characterised by an axial length scale of O( ). The porous outflow inner region is region V and the exterior outflow inner region is region VI.
(2.6) Substituting (2.6) into (2.1), we find that the leading-order governing equations for the exterior flow are simply the lubrication equations 0 = −p 0x + u 0zz , 0 = −p 0z , 0 = u 0x + w 0z for |x| > 1, 0 < z < 1. (2.7a−c) Substituting (2.6) into (2.2), the leading-order governing equations for the interior flow are At leading-order, the channel wall boundary conditions (2.3) become u 0 = 0 on z = 0, 1 for |x| > 1, (2.9a) W 0 = 0 on z = 0, 1 for |x| < 1, (2.9b) while the far-field condition (2.5) becomes dp 0 dx → −12 as x→ −∞. (2.10) We are unable to impose at leading-order the interfacial boundary conditions (2.4) in the outer regions; these conditions are satisfied within the inner regions described in § 2.3.2. For the outer problems, we must instead impose matching conditions (using the matching law given in Van Dyke (1975)). These are determined in § 2.3.2, and we state them here to close the leading-order outer problem. As the outer pressure is independent of z and the fluid velocity is proportional to the pressure gradient, we find that we may impose continuity of total flux and pressure at each interface (corresponding to Neumann and Dirichlet conditions on the pressure, respectively). That is, We note that whilst the conditions (2.11) close the outer problem, they give no information about, for example, the leading-order stress acting on the interfaces. The leading-order solutions are readily obtained as follows: in region I, in region IV, and in region VII, The pressure is determined up to an arbitrary constant which (at this order) we fix by supposing p 0 + 12x→ 0 as x→ −∞, without loss of generality. We remark that the presence of the porous plug introduces a jump in pressure of 24 − 2/K ∼ −2/K (as K is small). Equivalently, an O(1) pressure drop creates an O(K) flow through the porous plug. We note that the horizontal components of velocity (2.12a) and (2.13a) are incompatible with the continuity of flux condition (2.4a) on the interfacial inflow boundary for the full problem. Similarly for outflow, we find that (2.13a) and (2.14a) are incompatible with (2.4a). In the next section we introduce inner regions to resolve this issue, but for the remainder of this section we consider higher-order terms in the outer regions as the higher-order pressure is coupled to the leading-order flow in the inner regions.
Using the leading-order solutions (2.12) and (2.14), the O( ) terms in the governing equations for the exterior flow (2.1) become Similarly, using the leading-order solutions (2.13), the O( ) terms in the governing equations for the interior flow (2.2) become The channel wall boundary conditions (2.3) at this order are as follows while the far-field condition (2.5) implies dp 1 dx → 0 as x→ −∞. (2.18) The matching conditions across each interface are determined in § 2.3.2, and we give them here to close the problem. For inflow, we have 1 0 (2.19a,b) and for outflow we have where the pressure jumps Π − and Π + are functions of the parameters Re and K. In § 2.3.2, we specify the problems that need to be solved to determine these functions, while in § § 2.4 and 2.5 we determine their asymptotic behaviour for small and large Re.
Given Π ± , the solution to (2.15)-(2.20) is readily determined: in region I, in region IV, and in region VII, As with the leading-order solutions (2.12)-(2.14), the pressure field is unique up to an arbitrary constant. For convenience, we choose the constant such that P 1 (0, 0.5) = 0, without loss of generality. We proceed by considering the inner regions, highlighting the problems that would need to be solved numerically to determine the functions Π − and Π + .

Inner regions
The appropriate scalings in the inflow and outflow inner regions are x = ∓1 + X ∓ , respectively, with (w, W) = (w, W)/ . To make it clear that we are working in the inner regions we introduce overbars on the other dimensionless variables, writing (u, U) = (u, U), (p, P) = (p, P) and (u, Q) = (u, U)e x + (w, W)e z . Using vector operators in terms of X ∓ and z as appropriate, the exterior fluid governing equations (2.1) become Re(u · ∇)u = − −1 ∇p + ∇ 2 u, 0 = ∇ · u in regions II and VI, (2.24a,b) which are valid when X − < 0 for inflow (region II) and X + > 0 for outflow (region VI), with 0 < z < 1 in both cases. The interior fluid governing equations (2.2) become which are valid for X − > 0 for inflow (region III) and X + < 0 for outflow (region V), with 0 < z < 1 in both cases. Similarly to (2.6), we form inner expansions in powers of as follows: At leading order, it follows from (2.24a) that 0 = −∇p 0 in regions II and VI, (2.27) while (2.25a) gives, similarly, 0 = −K∇P 0 in regions III and V. (2.28) The appropriate interfacial conditions on X ∓ = 0 for the leading-order pressure are given by the leading-order terms in (2.4b), which become p 0 = P 0 on X ∓ = 0. (2.29) Thus p 0 and P 0 are constant and equal in regions II and III and in regions V and VI, leading to the matching condition (2.11b). Using the leading-order outer solutions (2.12)-(2.14) for the pressure, we find p 0 = 12 in region II, P 0 = 12 in region III, (2.30a,b) for inflow, and for outflow. We note that these leading-order pressures are constant on the interface, as predicted by Levy & Sanchez-Palencia (1975).
As the leading-order fluid pressures are constant within each inner region, the leading-order inner flow is driven by the first correction to pressure, and these equations are given by the O(1) terms in (2.24) and (2.25). We illustrate the leading-order flow problems in figure 4 for inflow, and in figure 5 for outflow. In both cases, the inner flow transition problem is governed by the two-dimensional steady Navier-Stokes equations coupled to the Darcy equations. To fully determine the inflow and outflow pressure jumps Π − and Π + , we would have to solve the problems indicated in figures 4 and 5 numerically for two parameters: Re and K. However, to couple the outer regions, we require another condition on the pressure. We note that we can analytically determine suitable coupling conditions for the horizontal components of velocity averaged over the channel height by appealing to global conservation of mass, rather than solving the full problem. By integrating the continuity equations across the channel height, using the conditions for the flow in the normal direction at each boundary and matching with the outer solutions, we deduce the following coupling conditions for the outer horizontal components of velocity (2.32a,b) We proceed by determining the inner flow in the asymptotic sublimits Re 1 and 1 1/K Re 1/ , starting with the former. This will allow us to calculate quantities that are unobtainable from sole consideration of the outer problem, such as the stress acting on the interface and the limiting behaviour of the functions Π − , Π + .
2.4. Small Reynolds number: Re 1 In the sublimit in which Re 1, the leading-order problems outlined in figures 4 and 5 become Stokes flow coupled with Darcy flow. Therefore, the flow equations are now reversible, and it suffices to solely consider inflow. We therefore use X − = X for ease of notation. Further, we note that the inflow and outflow pressure jumps are now linked via the expression Π + = −Π − .
We solve the resulting leading-order system numerically (for the streamfunction within the exterior flow and for the pressure within the porous flow) using a second-order accurate central finite-difference scheme. Due to the elliptic nature of the governing equations, we use an iterative method to bypass some of the cost of solving the fully coupled elliptic problems in one step. That is, we iterate between solving Laplace's equation for the pressure in region III, and the biharmonic equation for the streamfunction in region II, and use the solution found in the previous iteration to determine the boundary conditions for the current iteration. We truncate the infinite domain to [−1, 0] for X in region II and [0, 1] for X in region III (we note that an extended domain adds little accuracy to the scheme), using 300 grid points for X in each region and 300 grid points for z. The simulation is halted once the difference between successive iterations at each grid point is less than 10 −4 .
The normal and shear components of stress acting on the interface are respectively, as , Re→ 0. We show the O( ) stress results for K = 0.01 in figure 6, and calculate Π − = 1.92 to three significant figures. The leading-order normal stress term in (2.33a) takes the specific value of −12 because we have imposed P(0, z) = 0 to uniquely define the pressure. In fact, it is the O( ) terms in (2.33) which are more interesting, as these are the dominant terms which vary in z and, unlike the leading-order pressure, cannot simply be determined from the outer solution. Hence, the boundary layer analysis is essential to obtain these results. The reason the outer solutions cannot make good predictions of the z-variation of interfacial stress is because of the tendency of the exterior flow to move towards the channel edges in transition between Poiseuille and plug flow in the boundary region. The interfacial normal stress is greater/smaller than predicted by the outer flow towards the centre/edges of the channel, whilst the magnitude of interfacial shear stress is greater than predicted by the outer solution.
2.5. Large Reynolds number: 1 1/K Re 1/ We now consider the sublimit in which 1 1/K Re 1/ . We do not consider further the sublimits in which ReK = O(1) or 1 Re 1/K 1/ . This is because we are interested in the tissue engineering problem with porous plugs whose values of K are not too small. Moreover, the above sublimits yield coupled nonlinear problems which are not amenable to analytic study. We emphasize that the limit we consider is when the external fluid has a large Reynolds number, but the pore Reynolds number is small, so that the Darcy equations and continuity of pressure conditions are valid.
In this section we perform an asymptotic analysis using the small parameter δ = 1/Re 1 for ease of notation. Therefore, when we refer to, for example, 'leading order' in this section, it is meant in reference to δ. In the limit δ→ 0, there is a significant change in boundary layer structure within the inner regions II and VI compared to the O(1) Reynolds number case. The limit as δ→ 0 is singular, and we expect viscous effects to become important in regions near each boundary. Furthermore, as fluid momentum dominates within regions II and VI, the flow direction is important. Indeed, we obtain a different boundary layer structure for inflow and outflow, and therefore split the analysis into these two cases, starting with inflow.
In the limit of small δ, the first-correction terms u 1 , Q 1 , p 1 and P 1 from the asymptotic series in the outer regions (2.6) are now scaled with δ −1 to be consistent with the new inner scalings, as follows: (u 1 , Q 1 , p 1 , P 1 ) = δ −1 (u 10 , Q 10 , p 10 , P 10 ) + (u 11 , Q 11 , p 11 , P 11 ) + O(δ). (2.34) The scalings for the pressure jumps Π − and Π + will turn out to be our objective being to determine Π − 1 and Π + 0 . Though we have not yet fully described the inner asymptotic structure, we give the form of the asymptotic expansions within each region in table 2 for ease of reference.

Inflow
The inflow inner problem (within regions II and III) is shown in figure 4. The limit as δ→ 0 yields leading-order inviscid flow in region II away from the boundaries, with five additional boundary layers being required to satisfy all of the boundary conditions. Two of these viscous boundary layers are on the channel walls, one is on the porous interface, whilst the remaining boundary layers occur in the corners between the channel wall and porous interface. The boundary layer structure for the inflow problem is shown in figure 7.
We start by considering the coupled regions II and III, and show that the transition to Poiseuille flow occurs solely within region III at leading order in δ.

Region
Asymptotic expansions Relative scalings Asymptotic expansions within boundary layers in the large Reynolds number limit.

Regions II and III
Within regions II and III, we make the following asymptotic expansions as δ→ 0, where the non-integer powers of δ come from corrections to the flow that arise due to the presence of a boundary layer near the channel walls (regions IIb in figure 7). These terms are used in appendix A, but are not required further in this section. Substituting the asymptotic expansions (2.36) into the governing equations presented in figure 4, the O(1/δ) terms in region II are trivially satisfied, but so far Q 00

II
IIa III

IIb IIc
IIb IIc FIGURE 7. Schematic diagram of the inner regions IIa-c within the inflow boundary region II. The light region is the exterior fluid and the dark region is the porous obstacle. The size of the boundary layers has been exaggerated for illustrative purposes.
is unknown. Proceeding with the expansions, we find at O(1) that u 01 = (u 01 , w 01 ) is governed by The boundary conditions we impose at this order in region II are again obtained via matching. On the channel walls (via matching with region IIb in appendix A), we have no flux, so that while in the far field we match with region I, to obtain Similarly, the boundary conditions in region III are while matching with region IV gives We note that the far-field matching condition Q 00 ∼ e x as X − → ∞ follows consistently from (2.42) using (2.38a).
The interfacial conditions at X − = 0 for (2.37)-(2.38) are formally derived from consideration of the inner-inner boundary layer (region IIa) in § 2.5.3, within which we are able to satisfy the no tangential slip condition, but additionally find that the leading-order flow and the first correction to the pressure p 11 are unchanged. From § 2.5.3, the correct conditions to couple regions II and III are continuity of flux and pressure, with The system to solve in regions II and III (up to the O(δ) terms) is given by (2.37)-(2.43). We see that regions II and III have decoupled: the boundary condition (2.43a) allows us to solve for Q 00 and P 11 using (2.38) with (2.41)-(2.42), and then we can use the boundary condition (2.43b) to determine u 01 and p 11 using (2.37) with (2.39)-(2.40). The governing equations (2.38), with boundary conditions (2.41)-(2.42) and (2.43a), are solved by Since u 00 is independent of X − within region II, the leading-order transition from Poiseuille to plug flow occurs entirely within the porous medium in region III, and is governed by (2.44).
We now complete the solution in region II by determining u 01 and p 11 . As well as allowing us to determine Π − 1 , knowledge of p 11 gives the normal stress acting on the interface up to O( ), the first order at which the stress varies in the z-direction. A rearrangement of the exterior flow (2.37) allows us to decouple the region II system (where −∞ < X − < 0 and 0 < z < 1) into (2.45b) We solve (2.45a) for w 01 first, using boundary conditions (2.39) on z = 0, 1 and (2.40a) as X − → −∞. We derive a consistent boundary condition for w 01 on X − = 0 using the boundary condition (2.43b) for the pressure on X − = 0, in combination with the governing equation (2.37b), giving (2.46) The boundary condition (2.46) ensures that w 01 does not vanish on X − = 0, inducing a slip velocity on this interface within region II. This issue is taken care of in an inner-inner boundary layer closer to the interface (region IIa in figure 7) in § 2.5.3. The linearity of the governing equation (2.45a), with boundary conditions (2.39), (2.40a) and (2.46), also allows us to deduce that w 01 scales with K −1 . Further, from the continuity equation (2.37c) it follows that u 01 also scales with K −1 . We solve for w 01 using a standard central finite-difference scheme. By numerically integrating w 01 with respect to X − , we are able to plot the streamlines in figure 8. In figure 8(a), we show the streamlines for the first correction to the velocity, and in figure 8(b), we show the streamlines for the full velocity up to this first correction. We see that the fluid moves from the centre of the channel towards the channel walls, as expected. Therefore, adjustment from Poiseuille to plug flow occurs in region II, but at higher order. This movement induces a horizontal slip velocity on the channel walls, which is resolved in appendix A by introducing a boundary layer near the channel walls in which z = O(δ 1/3 ) on the bottom wall, with a similar layer near the top wall. These regions are labelled as IIb in figure 7. In appendix A, we show that the flow near z = 0 has the form where C k , µ k , and g k are described in (A1)-(A 2), Ai is the Airy function of the first kind and a k is defined below (A 10). From (2.47), we see that the leading-order solution is preserved near the channel wall. In fact, the slip is resolved within higherorder terms in the boundary layer, and does not affect the leading-order solution in region II. The problem for p 11 is given by (2.45b) with boundary conditions (2.40b) and (2.43b). We obtain consistent boundary conditions for p 11 on the channel walls by substituting the leading-order flow solution into the governing equation (2.37b) to obtain p 11z = 0 on z = 0, 1 for X − < 0. (2.48) We solve the full system for p 11 , given by (2.45b) with boundary conditions (2.40b), (2.43b), and (2.48) numerically using a central finite-difference scheme. Using a similar scaling argument as for w 01 , we find that Π − 1 scales with 1/K. We are able to calculate Π − 1 by plotting (p 11 + 12X − )K, and using the matching condition (2.40b). As Π − 1 is independent of z, and to more easily show the pressure drop on a graph, we plot ( 1 0 p 11 dz + 12X − )K for varying X − in figure 9. We find thereby that . The function ( 1 0 p 11 dz + 12X − )K for varying X − in region II, which allows us to determine the outer pressure drop for inflow.

Inner-inner region IIa
We now investigate the inner-inner boundary layer closer to the interface, region IIa in figure 7. This boundary layer is required to satisfy the no tangential slip condition on the interface, and captures a balance between the relevant inertial terms and the viscous terms induced by the presence of the porous interface.
The relevant scalings are X − = δ X and w 0 = δ w 0 , and we introduce the new variables u 0 = u 0 and p 1 = p 1 to signify that we are working in region IIa. Under these scalings, in region IIa (where −∞ < X < 0 and 0 < z < 1) the exterior flow equations in figure 4 become (2.49c) and the interfacial boundary conditions become The purpose of this boundary layer is to determine the tangential velocity w 0 , and to ensure that it is possible to satisfy the no tangential slip boundary condition on X − = 0. Therefore, we form asymptotic expansions as follows: The relevant far-field condition arises via matching with region II, giving where w 01 is determined in § 2.5.2. The interfacial conditions (2.53a,b) are used to couple regions II and III in § 2.5.2. The remaining terms in the system (2.52)-(2.54) yield w 01 , and are solved by ( 2.55) We see that (2.55) can only be a solution for inflow, as we require 6z(1 − z) X < 0 as we move back into region II. For outflow, this inner-inner region is not a valid boundary layer, and we must tackle the problem in a different manner. We note that this boundary layer persists for alternate tangential slip conditions. We can showcase the robust nature of the boundary layer via two examples. The first example uses the general tangential slip condition derived in Carraro et al. (2015). That is, instead of (2.53c), we consider a condition of the form w = λu on X − = 0, for 0 < z < 1, (2.56) which reduces to (2.53c) when λ = 0. Here, λ is a slip coefficient, obtained by solving the cell problem given in Carraro et al. (2015). Using this condition, the transverse velocity is in region IIa. In the second example, we use a Beavers and Joseph slip condition (Beavers & Joseph 1967) of the form whereλ is the experimentally determined slip coefficient. In this case, the transverse velocity would be Hence the boundary layer structure persists if the no-slip condition is modified in either of these two ways. Physically, we can conclude that the thickness of each boundary layer involved in inflow is unchanged for any of the three boundary conditions applied at the interface. Moreover, the flow induced by each of the three (a) The normal stress (we plot Kp 11 (0, z)) and (b) the shear stress (we plot 6z(1 − z)Kw 01 (0, z)). As both the normal and shear stress scale with 1/K at this order, both plotted terms are independent of K.
boundary conditions is only different in region IIa, the boundary layer closest to the interface.
Returning to the original no-slip boundary condition (2.50c), it is straightforward to calculate that the components of stress on the interface at x = −1 are given by as , δ→ 0, with δ 1. Whilst the normal stress (2.60a) is nominally O(1), the leading-order term comes from the decision to impose P(0, 0.5) = 0 to uniquely define the pressure. Thus, the more interesting result occurs at O( ), which is the lowest order at which the normal stress varies in the z-direction. For the shear stress (2.60b), the leading-order terms are O( ), and this varies in z. Those terms which vary in z cannot be determined solely from the outer solution: they are obtained from our boundary layer analysis and are each plotted in figure 10. We show in the previous section that p 11 and w 01 both scale with 1/K, thus the stress also scales with 1/K at O( ).
Additionally, since we have imposed P(0, z) = 0, we have no O( /δ) constant term from the pressure in the normal stress (2.60a). In three dimensions, the tangential (that is, in the (e x × e z )-direction) variation in pressure ensures that there will be an O( /δ) term in the corresponding normal stress component no matter how we define the pressure (as we show in § 3.6.1).
We have now solved for the flow variables within the inflow boundary layers for the regions that affect the stress and the pressure drop. As shown in figure 7, there are further boundary layers which are important for inflow, and we discuss these in appendix A.

Summary
The main part of the leading-order inflow may be summarized as follows. The leading-order axial flow and pressure are unchanged until they reach the interface, after which they transition to plug flow inside the obstacle in region III, whose length is comparable to the channel height. This has an effect on the correction to normal flow and pressure in region II. Meanwhile, the tangential velocity satisfies the no tangential slip condition on the obstacle surface in region IIa -a boundary layer of even smaller width, which is comparable to the length of the channel height multiplied by the reciprocal of the Reynolds number.
Combining the original asymptotic expansion (2.26) with the asymptotic expansion in region II, namely (2.36c), we deduce that the pressure jump between outer inflow regions I and IV is O( ). Specifically, we have determined that where α ≈ −0.117. The argument presented in this section will not hold if the fluid is leaving the obstacle rather than entering it. This is because the boundary layer region IIa is only valid for inflow (as can be seen from (2.55)). We reconsider our approach for outflow in the next section.

Outflow
The outflow inner problem (within regions V and VI) is shown in figure 5. We again seek a solution as δ→ 0 with δ 1 using the method of matched asymptotic expansions.
The asymptotic structure for outflow is shown schematically in figure 11. The flow in region V is unchanged at leading order and reaches the interface as uniform plug flow with magnitude 1. There is an inviscid core in the centre of region VI, flanked by Prandtl boundary layers near z = 0, 1 whose thicknesses are of O((δX + ) 1/2 ). These boundary layers (region VIb) grow downstream, and the transition to Poiseuille flow occurs when their thicknesses become of O(1). This occurs on a horizontal length scale of O(1/δ) in a transition region we denote region VIa. There is a third boundary layer, denoted region VIc, which occurs in the corner between the channel wall and the porous interface. This region deals with a singularity at the corner that arises in region VIb, as in the classic Prandtl problem.
In this section, we investigate the regions VI and VIa to understand the transition from plug to Poiseuille flow and to calculate the pressure jump between the outer regions IV and VII. Regions VIb and VIc are discussed in appendix B.1 since they are not important to the transition. Recall that, from (2.61), we calculated an O( ) pressure jump between outer regions I and IV for inflow. We show that the pressure jump between outer regions IV and VII is of O( /δ) for outflow.
We start by investigating how the fluid moves through region VI. The leading-order behaviour is plug flow, and we determine the effect of the interfacial no tangential slip condition on the correction to the plug flow. In particular, we determine the far-field behaviour of the flow, which allows us to impose the correct matching conditions for the problem of transition to Poiseuille flow in region VIa.

Inner regions V and VI
To determine the outflow pressure jump (2.35b), we must consider first-order flow terms. Therefore, we pose the asymptotic expansions 62a−d) as δ→ 0, where the non-integer powers of δ arise due to the correction terms from the boundary layers in region VIb near the channel walls. We determine the leading-order behaviour within the coupled regions V and VI, but do not fully solve for the higher-order terms due to their complexity. Instead, we analyse the higher-order problems to obtain the information required for the matching condition with the transition region VIa.
Substituting the asymptotic expansions (2.62) into the system outlined in figure 5, the leading-order equations are trivially satisfied. The O(δ 1/2 ) exterior flow equations are u 0aX + = −p 1aX + , w 0aX + = −p 1az , 0 = u 0aX + + w 0az in region VI, (2.63a−c) and the O(δ 1/2 ) interior flow equations are In region V, the first correction to the boundary condition on the channel walls is given by while the far-field conditions are obtained by matching with region IV, to yield In region VI, the boundary conditions on the channel walls are obtained by matching with region VIb. We show in appendix B.1 that, at this order, the flow problem within region VIb is equivalent to the standard Prandtl problem of uniform flow past a flat plate (Prandtl 1905). The matching conditions are therefore given by w 0a =β(4X + ) −1/2 on z = 0 for X + > 0, (2.67a) w 0a = −β(4X + ) −1/2 on z = 1 for X + > 0, (2.67b) whereβ (≈ 1.721) is a constant that is determined numerically from a nonlinear ordinary differential equation (discussed further in appendix B.1); the parameter β  (1970) and Wilson (1971) is related toβ viaβ = 2 1/2 β. Finally, the interfacial conditions on X + = 0 are u 0a = U 0a , p 1a = 0, w 0a = 0 on X + = 0 for 0 < z < 1. (2.68a−c) We do not solve the coupled system of equations (2.64)-(2.68) here, though we note that it is possible to obtain an implicit representation for w 0a via a Fourier transformation. Instead, we acquire the necessary information for the matching with region VIa, which allows us to determine Π + 0 in the next section. We note that (2.63) can be rearranged to deduce the following three results: firstly, that ∇ 2 w 0a = 0 in region VI; (2.69) secondly, that p 1a is the harmonic conjugate of w 0a ; and, thirdly, that p 1a + u 1a is independent of X + . We deduce the leading-order far-field behaviour of w 0a by considering the governing equation (2.69) together with the channel wall boundary conditions (2.67) as X + → ∞, to obtain the expansion Using the far-field condition (2.70) in the governing equations (2.63a,c), we further deduce that u 0a ∼β(4X + ) 1/2 , p 1a ∼ −β(4X + ) 1/2 as X + → ∞ for 0 < z < 1. (2.71a,b) We note that the far-field expansions (2.70)-(2.71) have no dependence on the permeability K. As these far-field expansions are used to form matching conditions, which transfer information into the transition region VIa (which is where we calculate the pressure drop Π + 0 ), we note that Π + 0 will be independent of K. We carry out this matching in appendix B.2.
The components of stress on the interface at x = 1 are given by as , δ→ 0 with δ 1, where we have not explicitly calculated p 1a , u 0a , or w 0a . In terms of the first order at which the stress varies in z, we find that the normal stress (2.72a) for outflow is a factor of O(δ −1/2 ) larger than for inflow (2.60a), while the shear stress (2.72b) for outflow is a factor of O(δ 1/2 ) smaller than for inflow (2.60b).

Transition region VIa
We now consider the region in which the flow transitions from plug to Poiseuille flow. This occurs when the growing Prandtl boundary layers with thickness of O((δX + ) 1/2 ) become of O(1), that is, when X + = O(1/δ), bringing about a balance in the exterior flow equations between fluid momentum and viscosity across the full channel width. Recall that we have called this region VIa, as illustrated in figure 11.
The relevant scalings are X + = δ −1 X and w 0 = δ w 0 , and we introduce the new variables u 0 = u 0 and p 1 = p 1 to signify that we are working in region VIa. Then, in region VIa (where 0 < X < ∞ and 0 < z < 1) the exterior flow equations given in figure 5 become (2.73c) We pose the asymptotic expansions as δ→ 0, and substitute them into the governing equations (2.73) to yield the following leading-order governing equations in region VIa: The boundary conditions on the channel walls become u 00 = 0 on z = 0, 1 for X > 0. (2.76) The leading-order matching conditions are derived in appendix B.2. They are obtained by forming composite expansions between regions VI and VIb, then using Van Dyke's matching rule (Van Dyke 1975). We use the multiplicative composite expansions (Van Dyke 1975) to ensure that the matching condition for the velocity satisfies the no-slip condition on the channel walls. This allows a greater accuracy when we solve numerically the resulting system. The matching conditions are given by u 00 ∼ (1 + 2β X 1/2 )f ( η (0) )f ( η (1) ), (2.77a) as X ↓ 0, where f is the standard Blasius similarity solution (defined in (B 7)), η (0) = z/ X 1/2 and η (1) = (1 − z)/ X 1/2 . The matching conditions (2.77) are consistent with a small X expansion (for fixed z) of the governing equations (2.75). The matching conditions between regions VIa and VII are given by where Π + 0 is determined from the solution in region VIa. The full system for the flow in region VIa is given by (2.75)-(2.78).
In Bodoia & Osterle (1961), this same transition region was considered, though was not derived formally through asymptotic methods. The authors used the boundary conditions u 00 = 1, w 01 = 0, p 10 = 0 on X = 0 for 0 < z < 1, (2.79a−c) instead of the matching conditions (2.77). We note that, unlike (2.77b), the condition (2.79b) is not consistent with the small X approximation to the governing equations (2.75), from which it can be deduced that w 01 is singular at the corners at X = 0, z = 0, 1. With the boundary conditions (2.79), it was found in Bodoia & Osterle (1961) that the pressure jump Π + 0 ≈ −0.338. We now calculate Π + 0 numerically using the formal matching conditions (2.77). We start the numerical simulation at X = 10 −4 to avoid the inverse square root singularities in w 01 and in p 10 X , and use the finite-difference scheme outlined in Bodoia & Osterle (1961). In figure 12, we show the numerical solutions for p 10 + 12 X as a function of X, and we emphasise that p 10 is independent of z. From the matching condition (2.78c), we determine that the value of the leading-order outflow pressure jump between outer regions IV and VII is given by Π + 0 = −0.331 to 3 significant figures. We note that the value of Π + found by Bodoia & Osterle (1961) using their inconsistent boundary conditions (2.79) only differs by approximately 2 % from the value of Π + obtained using formal matching conditions.
Combining the original asymptotic expansion (2.26) with the asymptotic expansion in region VIa (2.74c), we deduce that the pressure jump between outer outflow regions IV and VII is O( /δ). Specifically, we have determined that where Π + 0 ≈ −0.331.

Summary
We have thus far determined the boundary layer structure for two-dimensional flow through a channel with stationary walls containing a porous obstacle. We formulated the full problem for an O(1) Reynolds number, then considered the asymptotic sublimits in which Re 1 and 1 Re 1/ . In the first sublimit, the inflow and outflow boundary layer structure was found to be the same, whereas in the latter sublimit, a rich boundary layer structure was unveiled, which differed for inflow and outflow. We determined how flow transitions from Poiseuille to plug for inflow and from plug back to Poiseuille for outflow, which allowed us to determine the leading-order asymptotic values of the pressure jumps Π − and Π + across regions II/III and regions V/VI, respectively. For inflow, where the flow problems in regions II and III decouple, we calculated the stress acting on the interface up to the order at which there is a variation in z. For outflow, where the relevant flow problems in regions V and VI are coupled, we determined the scalings of the stress acting on the interface up to the order at which there is a variation in z.
In the next section, we generalize the problem and consider unsteady threedimensional flow in a Hele-Shaw cell past a tight-fitting, highly permeable cylindrical porous obstacle, with a smooth boundary. As the fluid can now travel around the porous obstacle, the high Reynolds number limit has an entrainment effect for outflow. We show that this entrainment can lead to a net force acting on the obstacle, even when the unidirectional far-field forcing is periodic with a zero mean.

Unsteady three-dimensional flow
A schematic of the three-dimensional problem is illustrated in figure 2(b), and a two-dimensional plan view is illustrated in figure 13. We initially work in a Cartesian coordinate system (x, y, z) with origin at the obstacle centre of mass projected onto the planar plate of the Hele-Shaw cell in the plane z = 0, the other plate lying in the plane z = h. The obstacle axis is normal to the two parallel flat sides of the Hele-Shaw cell, and is parallel to the z-axis. The typical cross-sectional extent of the obstacle is l, the Hele-Shaw cell height is h, and we shall also assume that = h/l 1.
The exterior and interior fluid velocities are denoted by u s = u s e x + v s e y + w s e z and Q s = U s e x + V s e y + W s e z , respectively. The exterior and interior fluid pressures are denoted by p s and P s , respectively. We use the superscript s to denote variables in the three-dimensional problem, and we later reduce the three-dimensional problem by relating these variables to their two-dimensional counterparts in the three-dimensional analogue of the inner regions considered in § § 2.3-2.5. The typical magnitude of the unidirectional far-field flow is U ∞ . We non-dimensionalize by scaling the variables as follows: (x, y) ∼ l, z ∼ l, t ∼ l/U ∞ (corresponding to an O(1) Strouhal number), (p s , P s ) ∼ µU ∞ /( h), (u s , v s , U s , V s ) ∼ U ∞ and (w s , W s ) ∼ U ∞ . With these scalings, the dimensionless Navier-Stokes equations, which hold inside the Hele-Shaw cell but outside the porous obstacle, are given by where Re = ρhU ∞ /µ is the global Reynolds number, ∇ is the three-dimensional gradient operator and ∇ 2 ⊥ = ∂ xx + ∂ yy is the two-dimensional Laplacian operator. The dimensionless Darcy equations, which hold inside the porous obstacle, are given by where K = k/h 2 1, as before.

Boundary conditions
The boundary conditions on the plates of the Hele-Shaw cell at z = 0, 1 are those of no slip and no flux for the exterior flow and of no flux for the interior flow, so that u s = 0 on z = 0, 1 outside the obstacle, (3.3a) W s = 0 on z = 0, 1 inside the obstacle. (3.3b) On the interface, we impose three-dimensional equivalents of the two-dimensional interfacial conditions given by (2.4). These are continuity of normal flux, continuity of pressure and a no tangential slip condition, as follows evaluated at the interface, where n is the unit normal pointing out of the porous obstacle.
The far-field condition is where G(t) is a dimensionless function of time. We take |G(t)| and |G (t)| to be of O(1), so that the far-field forcing does not have a large velocity or acceleration.
3.2. Asymptotic structure The dimensionless problem (3.1)-(3.5) is characterised by two length scales: the typical cross-sectional extent of the obstacle, 1, and the Hele-Shaw cell height, 1. We again seek a solution using the method of matched asymptotic expansions in terms of the small parameter . We note that, in contrast to § 2, not all of the fluid has to pass through the porous obstacle as it is now able to flow around the porous obstacle. In this section, we show that this difference results in a Darcy velocity of O(K) with an O(1) pressure drop, in comparison to the two-dimensional case where the Darcy velocity was of O(1) with an O(1/K) pressure drop. As the continuity of flux condition (3.4a) ensures that u s · n = O(K) near the interface, we are able to deduce that K = O( ) corresponds to an obstacle which is impermeable at leading order in . As discussed previously, we wish to investigate the case where the obstacle is as permeable as possible, and we therefore consider the limit where K 1. The local Reynolds number (defined later) characterises the flow close to the interface. We initially consider the distinguished limit whereby the local Reynolds number is of O(1), then consider the asymptotic sublimits of a small and large local Reynolds number in more detail, as before.
The outer problems are characterised by an O(1) length scale in the (x, y)-plane (with the Hele-Shaw cell height, , being the length scale in the z-direction), where we have at leading order Poiseuille and plug flow in the exterior and interior regions, respectively. We assume that the obstacle is a cylinder whose cross-section has a curvature of O(1), so that the surface near a point on the obstacle boundary is well approximated by the tangent plane at that point. The inner regions are then characterised by an O( ) extension of each local normal plane to the interface and occur on either side of the interface, so that the inner problems are similar to those considered in § 2, and the boundary layer structure is its three-dimensional equivalent.

Outer regions
In the outer regions, we take asymptotic expansions in the limit as → 0, of the form for f ∈ {u s , v s , w s , U s , V s , W s , p s , P s }. At leading-order we find that, exterior to the obstacle, where g 1 (z) = z(z − 1)/2, while interior to the obstacle, so that both the interior and exterior pressure satisfy Laplace's equation in two dimensions. We define the two-dimensional cross-section of the obstacle as D and the one-dimensional boundary of this as ∂D. Thus, the governing equations are The far-field condition is In § 3.4, we find that the leading-order coupling conditions on the interface are where (3.11b) is equivalent to the continuity of total flux condition given in (2.11a). The leading-order problem is defined by (3.9)-(3.11), and can be solved for a given obstacle boundary ∂D. For future reference, we deduce from (3.8) that the segments of the boundary where we have inflow/outflow at leading order occur when ∂P s 0 /∂n evaluated on the boundary is positive/negative.
For the O( ) correction terms we find that, outside D, while in the obstacle, so that the pressures satisfy The time derivative of ∇ 2 p s 0 does not appear in (3.14a) due to the governing equation (3.9a). The far-field condition is The first-order coupling conditions on the interface are determined in § 3.4 and appendix C, and are given by (3.16b) for u 0 · n = 0, where a superscript ending in −/+ refers to an inflow/outflow quantity, and H(x) is the Heaviside step function. These inflow/outflow segments of ∂D are determined from (3.7)-(3.11) as previously discussed. The presence of the Heaviside functions in (3.16) introduces the possibility of pressure singularities appearing at points on the boundary at which u 0 · n = 0. We discuss this in more detail once we determine the asymptotic approximations of the coefficients premultiplying these functions. The coupling condition (3.16b) represents conservation of mass across the interface at this order and, in contrast to the two-dimensional versions (2.19a), (2.20a), has additional terms which correspond to an entrainment effect at this order. That is, the jump in the O( ) average normal velocity across the interface (between outer regions) has a contribution from the variation of the leading-order tangential velocity in the inner regions.

Inner regions
As the inner regions are close to the obstacle boundary, it will be convenient to work in general curvilinear coordinates (n, s, z) such that n = 0 on the boundary, with n > 0 corresponding to the exterior of D, and s being the arc length along the boundary measured anticlockwise, as illustrated in figure 13.
The relevant inner scalings are n = N and (w s , W s ) = (w s , W s )/ , and we define the components of fluid velocity in curvilinear coordinates by u s = u s n + v s t + w s e z and Q s = U s n + V s t + W s e z in the exterior and interior regions, respectively. Under these scalings, the exterior flow equations (3.1) become (see, for example, Schlichting & Gersten (2000)) where κ(s) denotes the curvature of the obstacle boundary (positive if the centre of the osculating circle lies in the region in which n < 0 for a given s), and we assume that |κ| = O(1). Similarly, the interior flow equations (3.2) become We form inner expansions in powers of as follows: substituting these expansions into the inner exterior flow equations (3.17a,c) and the interior flow equations (3.18a,c), then equating powers of , yields the leading-order pressure equations This leading-order system is the same as in the two-dimensional case, and we proceed in exactly the same way. Matching with the outer pressures, and using the leading-order version of the continuity of pressure condition (3.20b), we deduce that the leading-order pressures are unchanged through the inner regions and we are justified in writing the coupling condition (3.11a). In a similar manner, integrating the leading-order version of the continuity equations (3.17d) and (3.18d) over the cell height, using the cell wall boundary conditions (3.19), and then matching with the outer regions, justifies the coupling condition (3.11b). The entrainment effect does not occur at this order. The first-order flow equations in the (N, z)-plane decouple from the tangential velocities v s and V s , and the tangential coordinate s can be treated as a parameter in the problem. We can further reduce these problems to the two-dimensional cases considered previously using the scalings where the variables on the right-hand side without a superscript s are variables from the two-dimensional problem, and ν(s, t) = − 1 0 u 0 · n dz = p s 0n (0, s, t)/12 = KP s 0n (0, s, t) is positive for inflow and negative for outflow. With N = −X − for inflow and N = X + for outflow we obtain the systems shown in figures 4 and 5, but with Re replaced by the appropriate local Reynolds number, namely Re(s, t) = |ν|Re. The three-dimensional pressure jump functions in (3.16a) are related to their two-dimensional equivalents via the expression (3.24) We have, therefore, already calculated the asymptotic behaviour of Π s− and Π s+ for the pressure coupling condition (3.16a) in § § 2.5.1 and 2.5.5. We see that Π s− , Π s+ → 0 as ν→ 0 and thus (3.16a) does not impose a pressure singularity at points where ν = 0. However, the possibility remains that (3.16b) could still impose a pressure singularity when ν = 0, in which case an additional inner region would be present. We proceed by assuming that |ν| = O(1), and discuss any pressure singularities when they occur. We now calculate the behaviour of the tangential velocities v s 0 and V s 0 within the inner regions (which are new to the three-dimensional problem). This will allow us to determine local properties near the interface, for example, the stress, and to couple the outer regions by determining the functions Λ s− and Λ s+ in the final coupling condition (3.16b).
Within the obstacle, the tangential velocity is given at leading order by V s 0 ≡ −KP s 0s (0, s, t). (3.25) The scaling v s 0 (N, s, z, t) = −(p s 0s (0, s, t)/12)v 0 (N, z; s, t), (3.26) where we reiterate that s and t are now parameters in the inner three-dimensional problem, results in the exterior flow equation Here, we introduce X such that X = X − < 0 for inflow, and X = X + > 0 for outflow, thus remaining consistent with the two-dimensional analysis. We emphasise that, as with the two-dimensional case, we expect the behaviour of the tangential velocity to vary depending on whether we have inflow or outflow. The cell wall boundary conditions for v 0 are given by v 0 = 0 on z = 0, 1 for X − < 0 and X + > 0, (3.28) the interfacial conditions are v 0 = 0 on X = 0 for 0 < z < 1, (3.29) and the matching conditions are v 0 → 6z(1 − z) as X − → −∞ and X + → ∞ for 0 < z < 1. (3.30) In appendix C, we integrate the O( ) continuity conditions (3.17d), (3.18d) over the channel height to obtain the expression encompassing the entrainment effect which occurs at this order. Using the outer solutions (3.7) and the tangential velocity scaling (3.26) in (3.31), we deduce that the inflow and outflow flux jumps Λ s− and Λ s+ in (3.16b) are given by where v 0 satisfies the system (3.27)-(3.30). We proceed by determining the leadingorder asymptotic behaviour of Λ s− and Λ s+ in the sublimits in which Re 1 and 1 1/K Re 1/ . 120 M. P. Dalwadi, S. J. Chapman, S. L. Waters and J. M. Oliver 3.5. Small local Reynolds number: Re 1 For the small local Reynolds number case, Re 1, the inertial terms in the governing equation (3.27) are not present at leading order in Re. Therefore, v 0 has no dependence on s or t (which come from Re(s, t)), nor on whether we have inflow or outflow. The local Reynolds number uses the local velocity component normal to the interface averaged over the channel height. Thus, the problem we present in this section governs the case when the normal local velocity is small, as well as the case when the global Reynolds number is small. The former case can occur, for example, in the transition between inflow and outflow for a large global Reynolds number, or when the obstacle is close to impermeable. For the latter reason, this problem also appears in the impermeable obstacle case considered in Thompson (1968).
We consider outflow without loss of generality, the O(1) terms in (3.27) yielding the governing equation (3.33) It follows from (3.33) and the outflow versions of the boundary conditions given in (3.28)-(3.30) that where α k = (2k + 1)π for non-negative integers k. Using this flow result, we can deduce that the shear stress in the tangential direction is where ζ is the Riemann zeta function. As inflow and outflow are reversible in the Stokes equations, the Heaviside functions in (3.16b) disappear and we can deduce that, in the limit in which Re 1, no pressure singularities exist up to O( ).
3.6. Large local Reynolds number: 1 1/K Re 1/ Now we consider the local asymptotic sublimit 1 1/K Re 1/ and, as in § 2.5, we find it useful to work with the inverse local Reynolds number, δ s = Re −1 1, for ease of notation. The asymptotic structure here is the same as in § 2.5; our goal is to determine the behaviour of the tangential velocity v 0 from the system (3.27)-(3.30), and use this to determine the asymptotic behaviour of the flux jump functions from (3.32).
We scale the first correction terms in the outer regions with the global Reynolds number as follows (u s 1 , Q s 1 , p s 1 , P s 1 ) = Re(u s 10 , Q s 10 , p s 10 , P s 10 ) + (u s 11 , Q s 11 , p s 11 , P s 11 ) + O(Re −1 ), (3.37) but the inner regions will be asymptotic expansions in the local Reynolds number. We split the analysis into inflow and outflow, starting with the former.

Inflow
The inflow behaviour is similar to the flow considered in § 2.5.1, for which the asymptotic structure is given in figure 7. We show that the leading-order tangential velocity is independent of X − through region II and region III, but the next-order flow depends on X − in region II. In this section we also show that the no tangential slip condition (3.29) is satisfied in region IIa, and we calculate the asymptotic behaviour of Λ s− as δ s → 0.
In region II, the tangential flow is unchanged at leading order and satisfies the asymptotic expansion Substituting (3.38) into the governing equation (3.27) and equating terms of O(1), gives with the matching condition v 01 → 0 as X − → −∞.
(3.40) Equation (3.39), with far-field condition (3.40), can be integrated numerically with respect to X − to determine v 01 . The resulting correction to the tangential velocity induces a slip on the cell walls, which is remedied by a boundary layer in region IIb (as illustrated in figure 7), the details of which are omitted from this paper for the sake of brevity. The coordinate scaling in region IIa is X − = δ s X, and we use hats to denote all variables in this region. We use the flow variables calculated in § 2.5.3, and pose an asymptotic expansion for the tangential velocity of the form v 0 ( X, z; s, t) = v 00 ( X, z; s, t) + O(δ s ), (3.41) where the leading-order term satisfies the governing equation (3.42) with the no tangential slip condition (3.29) yielding v 00 = 0 on X = 0, 0 < z < 1, (3.43) and a matching condition with region II given by v 00 ( X, z; s, t) ∼ 6z(1 − z) as X→ −∞, for 0 < z < 1. (3.44) The system (3.42)-(3.44) is solved by v 00 = 6z(1 − z)(1 − exp(6z(1 − z) X)). (3.45) The leading-order contribution to Λ s− in (3.32a) is of O(δ s ), and comes from both the O(δ s ) term in region II and the O(1) term in region IIa, as region IIa is O(δ s ) smaller than region II. Therefore, substituting the asymptotic expansion (3.38) into (3.32a), and using the governing equation (3.39) to relate the velocity v 01 to w 01 (which has already been calculated numerically in § 2.5.2), along with the leading-order solution for the tangential velocity in region IIa (3.45), the leading-order contribution to (3.32a) becomes as , 1/Re→ 0 with 1/Re K 1. We note that the integral is finite because w 01 → 0 as z→ 0, 1.
In contrast to the two-dimensional case, the pressure now varies in the s-direction. Therefore, the component of the interfacial stress in the normal direction is different from the two-dimensional case (2.60a), and in three dimensions is given by as , 1/Re→ 0 with 1/Re K 1. We note that P s 0 and P s 10 have no z-dependence. Therefore, as in two dimensions, the first order at which the normal stress varies in z is at O( ). Finally, we determine the leading-order term in the component of the shear stress on the obstacle boundary in the s-direction for inflow, σ ns . From the tangential flow solution in region IIa given by (3.45), and recalling the scaling (3.26), we deduce that as , 1/Re→ 0 with 1/Re K 1.

Outflow
Outflow occurs in a similar manner to § 2.5.5, for which the asymptotic structure is illustrated in figure 11. The tangential velocity within region V is plug flow (as given by (3.25)) and the tangential flow in region VI is of O(δ s ). The tangential flow then transitions to Poiseuille flow in region VIa, where it becomes of O(1). In this section we solve for the leading-order tangential flow for outflow and calculate the asymptotic behaviour of Λ s+ as δ s → 0.
In region VI, we use the flow variables calculated in § 2.5.6 as well as the asymptotic expansion which satisfies the no tangential slip condition (3.29). The no-slip condition on the cell walls is satisfied in region VIb (as illustrated in figure 11) via a similarity solution, as described in appendix B.1. As with the three-dimensional inflow case considered in the previous section, the component of the interfacial stress in the normal direction is now σ nn = −(P s 0 (0, s, t) + ReP s 10 (0, s, t)) + O( /δ 1/2 s ), (3.50) as , 1/Re→ 0 with 1/Re 1, and where P s 0 and P s 10 have no z-dependence. Therefore, as in two dimensions, the first order at which the normal stress varies in z is at O( /δ 1/2 s ). From (3.49) and the scaling (3.26), we deduce that the stress on the interface in the tangential direction is σ ns ∼ − p s 0s (0, s, t) 12 v 0X + (0, z; s, t) = ν(s, t)Re p s 0s (0, s, t), (3.51) as , 1/Re→ 0 with 1/Re K 1. The transition to Poiseuille flow occurs in region VIa, within which the relevant coordinate scaling is X + = δ −1 s X. We further use v 0 = v 0 to indicate that we are in this intermediate region, and use the variables introduced in § 2.5.7. We substitute the asymptotic expansion (3.52) into the governing equation (3.27) to obtain the leading-order governing equation (3.53) The cell wall boundary conditions (3.28) become v 00 = 0 on z = 0, 1 for X > 0, (3.54) the matching condition with region VI is v 00 → 0 as X ↓ 0 for 0 < z < 1, (3.55) and the far-field behaviour is determined via the matching condition (3.30), and is given by v 00 ∼ 6z(1 − z) as X→ ∞ for 0 < z < 1. (3.56) The system (3.53)-(3.56) details the transition from plug to Poiseuille tangential flow. To determine the flux jump Λ s+ at leading order, the system must be solved numerically. We note that the system (3.53)-(3.56) is independent of K, and hence only needs to be solved once to determine Λ s+ at leading order.
Accounting for the change of variable X + = δ −1 s X being dependent on s, we find and it follows that the leading-order behaviour of Λ s+ is given by (3.58b) We calculate v 00 numerically using a second-order central finite-difference scheme, then determine Λ a using the trapezium rule. We find that Λ a ≈ 0.125 to three significant figures. Combining the outflow flux jump (3.58) with its inflow equivalent (3.46), and recalling that ν = KP s 0n (0, s, t), we find that the leading-order (in δ) coupling condition (3.16b) in the large Reynolds number asymptotic sublimit is given by (3.59a) The leading-order pressure jump condition (3.16a) in the large Reynolds number sublimit is given by We see that the coefficient premultiplying the Heaviside function in (3.59a) does not necessarily vanish when u 0 · n = 0, whereas the coefficient premultiplying the Heaviside function in (3.59b) does vanish. Thus, at these points the pressure is continuous but the pressure derivative along the boundary will have a log singularity.

Summary
We now summarise the main results for the three-dimensional case. The general outer problem is governed by the two-dimensional linear equations: outside D The coupling conditions (in outer variables) on the boundary of the obstacle crosssection ∂D are given by Here, ∂/∂n denotes the outward normal derivative to the obstacle, and we define the functions Π (s, t; Re, K) = Π − (s, t; Re, K)H(P s 0n ) + Π + (s, t; Re, K)H(−P s 0n ), (3.63a) Λ(s, t; Re, K) = Λ s− (s, t; Re, K)H(P s 0n ) + Λ s+ (s, t; Re, K)H(−P s 0n ), (3.63b) where Π − and Π + are outputs to the nonlinear systems presented in figures 4 and 5, Λ s− and Λ s+ are solutions to (3.32a) and (3.32b) respectively and H(x) denotes the Heaviside step function.
On the boundary layer structure near a highly permeable porous interface

125
For , Re→ 0, we found that Π − = −Π + and these functions can be numerically determined by solving a system of linear equations, as described in § 2.4. Additionally, from (3.36), we found that Λ s− , Λ s+ ∼ 93ζ (5) 2π 5 p s 0ss (0, s, t) as , Re→ 0, (3.64) where ζ is the Riemann zeta function. For , 1/Re→ 0 with 1/Re K 1, we found in § § 2.5.1 and 2.5.5, respectively, that recalling that Re(s, t) = K|P s 0n |Re, where α ≈ −0.117 and Π + 0 ≈ −0.331. In the same limit, we also found that Λ s− = O(1/(ReK)) (the prefactor at this order of magnitude is given in (3.46)), and that Λ s+ ∼ −ReKΛ a ∂ ∂s (P s 0n (0, s, t)p s 0s (0, s, t)) (3.66) as , 1/Re→ 0 with 1/Re K 1, where Λ a ≈ 0.125 (defined in (3.58b)). In the limit in which , Re→ 0, we found that the interfacial stress σ (n)| ∂D = −p s 0 n| ∂D + O( ), where n is the unit normal vector on the interface directed towards the exterior fluid. The prefactors in the error term can be determined by solving a linear system of equations. For these stress components in the n-and z-directions, we must solve the problem of Stokes flow coupled to Darcy flow, and present the system to be solved numerically in § 2.4. The error term due to the stress component in the s-direction is given in (3.35). For the stress components in the n-and z-directions, we showed in figure 6 the O( ) stress acting on the interface, the first order at which the stress varies in z.
In the limit in which , 1/Re→ 0 with 1/Re K 1, the interfacial stress is given by σ (n)| ∂D ∼ −(P s 0 (0, s, t) + ReP s 10 (0, s, t))n − ReK(3P s 0n (0, s, t)p s 0s (0, s, t)(z(1 − z)) 2 )H(P s 0n )t, ( 3.67) where t = e z × n is the unit tangent vector on the interface in the anticlockwise direction. The given stress components in the normal direction in (3.67) do not vary in z, and the order of magnitude of the error is different for inflow and for outflow. with dimensionless radius 1. Such a set-up can arise in tissue engineering, where the interior of the porous obstacle is seeded with cells. Another application is in modelling the erosion of a porous biofilm, for example, the attempted removal of dental plaque from between teeth. In these applications, quantities of physical interest can be determined from our analysis. Firstly, the net force acting on the porous obstacle determines the movement of the porous obstacle were it free to move. Although a periodic forcing in the low Reynolds number regime would lead to no net movement over one oscillation (Purcell 1977), the same is not true for a large Reynolds number regime, and we investigate the latter case here. Secondly, in tissue engineering applications, cell growth is often coupled to the shear stress that cells experience, which is one part of a process known as mechanotransduction. Thus, to help with cell placement within a porous scaffold, it is important to know how the internal shear stress varies within an obstacle. We use results from Whittaker et al. (2009) to estimate the spatial variation of shear stress within the porous obstacle, and thus the spatial variation of shear stress experienced by cells placed within such an obstacle. Finally, biofilm erosion is often taken to be proportional to the square root of the shear stress acting on the interface (Duddu et al. 2009) and, to this end, we determine the interfacial shear stress.
We use cylindrical polar coordinates (r, θ, z), where r = 0, z = 1/2 corresponds to the centre of the porous obstacle, and define the components of the exterior and interior fluid velocities as u s = u r e r + u θ e θ + u z e z and Q s = U r e r + U θ e θ + U z e z , respectively.
We consider the (more interesting) case where , 1/Re→ 0 with 1/Re K 1, and use the corresponding outer system as outlined in § 3.7. That is, we consider asymptotic expansions of the form u s ∼ u s 0 + Reu s 10 as , 1/Re→ 0 with 1/Re K 1, (4.1) with similar expressions for Q s , p s and P s . The governing equations for the exterior fluid are |∇p s 0 | 2 = 0 for r > 1, 0 < θ 2π, 0 < z < 1; (4.2a,b) and, for the interior flow, The coupling conditions on r = 1 are given by where K is the dimensionless permeability, Π + 0 , Λ a and Λ b are constants defined in § 3.7, and H(x) is the Heaviside step function. The far-field conditions are p s 0 ∼ −12xG(t), p s 10 → 0 as r→ ∞. (4.5a,b) The leading-order problem is solved by 12K). The points on the obstacle boundary at which the flow transitions between inflow and outflow are when P s 0r = 0 which, for this specific problem, is when θ = π/2, 3π/2, and the region θ ∈ (π/2, 3π/2) admits inflow when G > 0, and outflow when G < 0. The first-order problem is solved by a n (t)r −n cos nθ , (4.7a) where the coefficients a n , b n can be determined via a standard application of Fourier series, and we give b n (which are used in the following analysis) in appendix D. We note that the even modes (apart from n = 0, 2) vanish, and (a 2n+1 , b 2n+1 ) ∼ ((−1) n /(2n + 1) 2 )(α, β) as n→ ∞ (where α and β are constant). Thus, at the points θ = π/2, 3π/2 (where there is a transition between inflow and outflow) the O( Re) pressures are continuous but their derivatives with respect to θ are singular, as expected. Using (3.67), the total force acting on the obstacle, F, is given by as , 1/Re→ 0 with 1/Re K 1. Thus the only contribution to the total force acting on the obstacle from P s 10 is from the coefficient b 1 (t), given by Thus, the force (4.8) is given by as , 1/Re→ 0 with 1/Re K 1.
If the far-field forcing is periodic, with zero mean, the force acting on the obstacle does not necessarily vanish. That is, if there exists T > 0 such that G(t + T) ≡ G(t) for all t, and where T 0 G(t) dt = 0, the net force experienced by the obstacle over one oscillation is given by .11) as , 1/Re→ 0 with 1/Re K 1. The integral on the right-hand side of (4.11) vanishes for single-mode oscillations (for example, G(t) = cos t) but, in general, not for more complicated periodic functions (for example, G(t) = cos t − cos 2t). Even though the forcing is periodic, the effect of fluid inertia is to impart a net force over one oscillation. This would cause the obstacle to drift were it free to move. We note that the right-hand side of (4.11) vanishes as K→ 0, and hence the net force would not appear at this order for impermeable obstacles. We can attribute the appearance of the net force to the entrainment effect that occurs for outflow, which we were able to capture via a formal boundary layer analysis.
In tissue engineering applications, the shear stress S experienced within the porous obstacle is of interest, as cells are placed within the porous obstacle and their growth may be dependent on the shear stress. As different cells will require different levels of shear stress to grow optimally, it is important to understand how the shear stress experienced by cells will vary across the porous obstacle. In Whittaker et al. (2009), it was shown that the shear stress experienced by cells within the individual porous scaffold pores can be estimated from the Darcy velocity. In particular, the shear stress is proportional to the magnitude of the Darcy velocity, which can be obtained by the pressure solutions (4.6b) and (4.7b), to obtain (n + 1)b n+1 (t)r n cos nθ . (4.12) We can see from (4.12) that the spatial variation of shear stress through the obstacle occurs at O( ReK). We show that there can be significant spatial variation of shear stress within the obstacle (figure 14). From the coefficients b n (t), given in appendix D, we can deduce that the spatial variation is symmetric across x = 0 when T 0 G|G| dt = 0. Thus, all single-mode oscillations will cause an internal shear stress which is symmetric across x = 0 at O( Re), for example G(t) = (cos t)/ √ π (figure 14(a)), but this will not be the case for more complicated forcing functions, such as G(t) = (cos t − cos 2t)/ √ 2π ( figure 14(b)). The prefactors of these functions are chosen such that 2π 0 G 2 (t) dt = 1. We see that there are apparent singularities in the internal shear stress at (x, y) = (0, ±1), these are due to the logarithmic singularities that we have previously discussed. Although the shear stress will increase near these points, their singular nature is unphysical and the shear stress will be bounded; we expect that a formal realization of this could be achieved by considering boundary regions near these points if required.
We can also apply the results of our boundary layer analysis to predict the erosion of a porous biofilm, by calculating the shear stress acting on the interface. The square root of this shear stress is proportional to the interfacial erosion due to the flow (Duddu et al. 2009). Using the flow solution (4.6) in (3.67), we see that the leading-order shear stress, τ , is On the boundary layer structure near a highly permeable porous interface . This term is proportional to the shear stress experienced within the porous obstacle. We use far-field forcing functions of (a) G(t) = cos(t)/ √ π, (b) G(t) = (cos t − cos 2t)/ √ 2π. The functions are chosen such that 2π 0 G 2 (t) dt = 1. In both figures, we take K = 0.01, we use contour spacings of 0.05 and, as there are logarithmic singularities at (x, y) = (0, ±1), we do not plot contours of values larger than 0.3. We see that the symmetry exhibited in (a) is not present in (b). This is due to the effect of fluid inertia.
where t = e z × n is the tangential vector on the interface in the direction of increasing θ. The Heaviside step function in (4.13) means that τ vanishes for θ ∈ (π/2, 3π/2) when G < 0, and for θ ∈ (−π/2, π/2) when G > 0. For a periodic far-field forcing with zero mean, the temporal average over one oscillation of the square root of the magnitude of shear stress is 1 T T 0 |τ | 1/2 dt = Γ z(1 − z)| sin 2θ| 1/2 , (4.14a) |G|H(G(t)) dt for θ ∈ (π/2, 3π/2), T 0 |G|H(−G(t)) dt for θ ∈ (−π/2, π/2). (4.14b) Here, (4.14a) is proportional to the initial biofilm erosion. Thus, we have deduced that this erosion has a shape of z(1 − z)| sin 2θ| 1/2 , repeated in each quadrant of the interface. The magnitude of this erosion will be the product of Γ and the constant of proportionality linking erosion to the square root of shear stress. From the form of Γ , given in (4.14b), we see that the magnitude of erosion is not necessarily symmetric across x = 0 and the asymmetry will depend on G(t). The interfacial erosion is locally maximal at (θ , z) = (π(1 + 2n)/4, 0.5) for n ∈ Z and, in each quadrant, decreases from these points of maxima with two lines of symmetry, one with constant z and the other with constant θ ( figure 15). . This shape is repeated in each quadrant of the circular cylinder interface, but the magnitude will vary, as discussed in the text. We use contour spacings of 0.05.

Discussion
We make extensive use of the method of matched asymptotic expansions to investigate the laminar flow around and through a porous obstacle in both a channel and a Hele-Shaw cell. In particular, we determine how the flow behaves close to an interface (in a region whose width is of the order of the small aspect ratio) between single-phase and porous flows (governed by the Navier-Stokes and Darcy equations, respectively) for both small and large Reynolds numbers.
Our analysis allows us to resolve everywhere the leading-order fluid velocity, a necessary condition to start investigating the nutrient transport in the system, with the eventual goal to optimize nutrient delivery to cells within the porous obstacle. Furthermore, we obtain important characteristics of the inner flow (in a sense made formal within the main text), such as the interfacial stress, and determine suitable conditions to couple the (simpler) outer flows. The analytical approach reveals the dependence of the flow on the Reynolds number and the dimensionless permeability, without resorting to numerically expensive parameter sweeps for an inherently nonlinear problem. Importantly, we note that any variation of the interfacial stress in the direction transverse to the channel/cell cannot be determined solely from the outer solutions, and requires an in-depth boundary layer analysis. In particular, we must continue past the leading-order terms in the asymptotic expansions to obtain the relevant contributions to the stress.
In the first part of this paper we consider two-dimensional steady flow through a channel fully blocked by a finite-length porous obstacle and, in particular, examine the regions near the interface between single-phase and porous flow. These inner regions satisfy the full two-dimensional Navier-Stokes equations, and we calculate the pressure and flux jumps between outer regions in the asymptotic limits of small and large Reynolds number.
Whilst we include results for the small Reynolds number limit for comparison and completeness, the case that is more physically relevant to tissue engineering applications is the large Reynolds number limit. In this paper, we show that the latter case also exhibits more interesting flow behaviour, and we unveil a rich boundary layer structure. The leading-order (in terms of the large Reynolds number) transition between Poiseuille and plug flow occurs within the porous medium for inflow, and outside the porous medium for outflow. Therefore, the boundary layer structure is different for inflow and outflow, reflecting the irreversible nature of the large Reynolds number flow. We determine that the pressure jump is an order of the Reynolds number larger in magnitude for outflow than for inflow. There is no flux jump between outer regions, due to the restriction of fluid movement to a two-dimensional plane.
In the second part of this paper we extend the two-dimensional work to consider unsteady three-dimensional flow in a Hele-Shaw cell containing a cylindrical porous obstacle (which fits tightly between the parallel planar plates of the cell) whose cross-sectional boundary is smooth. Due to the smooth boundary, the flow behaviour and boundary layer structure close to the interface are similar to the two-dimensional channel flow considered previously and we can reduce much of the three-dimensional problem to the two-dimensional problem considered in the first part of the paper. However, in contrast to the two-dimensional channel flow case, a small flux jump between the outer flows is now present, as the tangential flow allows a small amount of fluid to move around, instead of through, the porous obstacle. We summarise the main results of the three-dimensional problem in § 3.7, including the results for the interfacial stress and the derived conditions required to couple the outer flows.
In § 4 we apply the general results to a three-dimensional cylinder with a circular cross-section in a flow with a far-field forcing. From this, we are able to make several physical predictions. We deduce that, due to the fluid inertia, a periodic forcing with zero mean could cause a drift on the cylinder were it free to move. This is relevant to the tissue engineering application we discussed in § 1, where it is important to control the movement of the porous obstacle, which has implications for the delivery of nutrients to cells within the construct. Moreover, knowledge of this long-term drift may be helpful in flushing blockages from pipes. We also calculate the spatial variation in internal shear stress through the obstacle. As cell growth can be coupled to shear stress, this has important implications for cell placement in tissue engineering applications. We show that the spatial variation in shear stress is symmetric through the obstacle for a single-mode far-field oscillatory forcing, but can break symmetry for multiple-mode far-field oscillatory forcings. We are also able to determine the interfacial stress and use this result to obtain the initial interfacial erosion of a porous biofilm, assuming that interfacial erosion is proportional to the square root of the modulus of shear stress with a given constant of proportionality.
In the large Reynolds number limit, we apply three different tangential slip conditions on the interface for inflow to see the effect of varying boundary condition. We find that any of the three tangential slip conditions we apply produce the same leading-order coupling conditions for inflow. Thus, knowledge of the tangential velocity condition is important for inflow if local information about, say, the interfacial stress were required, and it is important for outflow if global information about the outer problem is required. Additionally, we note that although we imposed continuity of pressure on the full problem, the coupling condition (3.62a) would remain the same if we had imposed continuity of stress instead. This is because the viscous term in the normal stress does not contribute to the stress at leading order. The higher-order pressure coupling condition (3.62b) would remain unchanged in the large Reynolds number limit for the same reason, but would be different in the small Reynolds number limit.
In this work, we have considered a porous medium within which the flow is governed by Darcy's equations. We restrict ourselves to these governing equations because the tissue engineering application we are interested in requires a porous scaffold that maintains its structural integrity whilst moving within the flow. Thus, at the very least, the underlying solid matrix of which the porous medium consists must be connected, and Brinkman's equations will not apply (Auriault 2009). However, there may be applications that are not as restrictive, where Brinkman's equations may apply. In such cases, the coupling conditions for the outer problem would be different to those derived in this paper, and a new analysis would be necessary to determine these. Whilst the issue of coupling the plain and porous flow regions is numerically simpler (as Brinkman's equations are able to be applied to both plain and porous fluid domains), the presence of the Laplacian viscous term means that it is more difficult to make analytic progress with the resulting partial differential equations. For example, the boundary layer structure may be significantly more complex within the porous medium.
Although we have considered a porous obstacle whose interfacial boundary is normal to the channel or bioreactor walls, this may not be the case. We show in § 4 that erosion of an initially straight wall is not uniform in any tangent direction to the interface. More generally, manufacturing constraints could cause the porous insert to have non-straight walls. Such imperfections will only change the flow problem if the gradient of the walls is significant in one of the boundary layers. If this occurs, the change in problem geometry yields a complicated extension to the problem we have presented. However, this extension is vital if the full dynamic problem of erosion is to be considered in, for example, the shape evolution of a biofilm.
The tissue engineering experiment discussed in § 1 involves a moving porous obstacle, which is possible as there are small gaps between the flat sides of the porous obstacle and the flat sides of the bioreactor (figure 1). Whilst our work has considered unsteady flow, we have restricted ourselves to a pinned obstacle in this paper, and neglected the small gaps. To consider the dynamic problem fully, we must investigate the role of these gaps. The presence of gaps would have no effect on the flow up to the asymptotic orders considered if these gaps were smaller than the size of the boundary layers presented (Dalwadi 2014). However, a gap height of the same order as the width of one of these boundary layers would significantly complicate the boundary layer problem, as a result of the significant change in the geometry of one or more domains in which we solve various submodels arising in the boundary layer analysis. Knowledge of the work presented here allows the porous obstacle trajectory to be calculated once the effect of the small gaps has been calculated. As the obstacle movement would affect the flow, this would allow a deeper understanding of the dynamics of nutrient transport and, ultimately, tissue growth.
The analytic nature of our results can significantly reduce the numerical expense of solving such dynamical problems. This is because the coupling conditions and stress components that we have determined are all in terms of the outer variables. Hence, we have reduced a nonlinear three-dimensional problem to a linear two-dimensional problem. More generally, this work highlights how the exploitation of an underlying separation of scales can significantly reduce the computational complexity of a physical problem.