A homogenised model for flow, transport and sorption in a heterogeneous porous medium
Lucy C. Auton, Satyajit Pramanik, Mohit P. Dalwadi, Christopher W. MacMinn, Ian M. Griffiths
AA homogenised model for flow, transport and sorption in a heterogeneousporous medium
Lucy C. Auton,
1, 2
Satyajit Pramanik,
2, 3
Mohit P. Dalwadi, Christopher W. MacMinn, and Ian M. Griffiths ∗ Mathematical Institute, University of Oxford, Oxford, OX2 6GG, UK Department of Engineering Science, University of Oxford, Oxford, OX1 3PJ, UK Department of Mathematics, Indian Institute of Technology Gandhinagar,Palaj, Gandhinagar – 382355, Gujarat, India (Dated: January 20, 2021)
Abstract
A major challenge in flow through porous media is to better understand the link between pore-scale mi-crostructure and macroscale flow and transport. For idealised microstructures, the mathematical frameworkof homogenisation theory can be used for this purpose. Here, we consider a two-dimensional microstruc-ture comprising an array of circular obstacles, the size and spacing of which can vary along the length ofthe porous medium. We use homogenisation via the method of multiple scale to systematically upscale anovel problem that involves cells of varying area to obtain effective continuum equations for macroscaleflow and transport. The equations are characterized by the local porosity, an effective local anisotropicflow permeability, and an effective local anisotropic solute diffusivity. These macroscale properties dependnon-trivially on both degrees of microstructural geometric freedom (obstacle size and spacing). We takeadvantage of this dependence to compare scenarios where the same porosity field is constructed with dif-ferent combinations of obstacle size and spacing. For example, we consider scenarios where the porosityis spatially uniform but the permeability and diffusivity are not. Our results may be useful in the design offilters, or for studying the impact of deformation on transport in soft porous media. ∗ ian.griffi[email protected] a r X i v : . [ phy s i c s . f l u - dyn ] J a n . INTRODUCTION Fluid flow and solute transport in porous media occurs in a wide variety of situations, includingcontamination transport [5, 18], lithium-ion batteries [14], hydrogeology [12], biofilms [11], bones[13] and soil transport [9]. Macroscopic flow, transport and sorption in porous media dependcritically on pore-scale quantities such as pore structure, connectivity of the fluid domain and thelocalised fluid–solid interactions. Due to the geometric complexity of the domains of the fluid andsolid phases, macroscopic transport and sorption equations that are uniformly defined throughoutthe entire porous medium are desirable. Here, we aim to improve the understanding of the rolethat microstructural heterogeneity and anisotropy has on macroscale models for fluid flow, solutetransport and sorption within porous materials with a simple two-dimensional model.Many porous materials are intrinsically heterogeneous at the pore-scale, including soils, rocksand biological tissues, however, the effect of these pore-scale heterogeneities on the macroscalebehaviour is often overlooked. Peat soil is intrinsically heterogeneous and anisotropic, whichleads to complex flows [1]. Beckwith et al. [1] state in their study that it is heterogeneity and notanisotropy that is the main cause of complex patterns of flow. Wang et al. [24], however, statethat soil anisotropy is crucial in peatland hydrology and highlight that additional research is re-quired to fully understand effects of soil anisotropy on solute transport. Clavaud et al. [6] studythe effect of the pore space geometry on the anisotropy of permeability on sandstone, limestone,and volcanic rocks through imaging. They find that, for each type of rock, there are differentrelationships between the pore-structure and permeability anisotropy and state that further workis needed to fully understand this relationship in limestone and volcanic rocks. O’Dea et al. [16]present a model broadly applicable to tissue-engineering applications, finding that the microstruc-ture induces anisotropy in the flow properties of the porous material and highlight the importanceof including microstructural information to enable the flow dynamics and nutrient delivery to becorrectly determined. Changes on the pore-scale can lead to large deviations from macroscopicmodels derived with a homogeneous microstructure. Rosti et al. [19] find that microstructuralchanges due to deformation of the solid phase constituting the porous matrix caused a breakdownof Darcy’s law.Three-dimensional structures are commonly approximated by a simplified two-dimensionalanalogue [3, 17]. Blum et al. [3] consider the effect of cytoskeletal structure on diffusion of parti-cles that do not bind to the cytoskeletal structure. A model they use for the cytoskeletal structure is2 two-dimensional periodic lattice composed of circles; they find that this two-dimensional modelgives generally the same result as the three-dimensional model. Non-woven filters composed offibres form a major part of the filtration industry and are found in numerous filters including airpurifiers and vacuum cleaners [22]. Printsypar et al. [17] consider a non-woven filter-mediumwith a microstructure comprising unidirectional fibres and model this with a two-dimensional ar-ray of circles. Arrays of micropillars are common in microfluidics and, provided each pillar ismuch taller than both its width and its separation from other pillars (such as in 2 and 25), theflow through a cross-section can also be modelled via a two-dimensional array. Similarly, flowthrough magnetic-separation filters, composed of wire wool, can be simply modelled via a two-dimensional cross-section [15].Porous media, and two-dimensional cross-sections thereof, are characterised by at least twodistinct length scales: the characteristic length of each pore/solid grain (pore-scale) and the char-acteristic length of the porous medium itself (macroscale) [23]. Obtaining a full understandingof the role of the pore structure on flow, transport and sorption via direct numerical solution ina complex geometry is computationally expensive, and prohibitively so when the pore-scale andmacroscale have very different lengths. In the latter case, one way to deal with such disparities isto systematically derive an upscaled macroscale model that is uniformly valid on the entire porousmedium, and contains pertinent pore-scale information in terms of an effective permeability, aneffective diffusivity and a source/sink term.There are many common methods for upscaling equations including the method of moments,renormalisation group theory, and homogenisation via volume averaging or the method of mul-tiple scales (MMS) [20, 28]. Both these methods of homogenisation yield the same macroscaleequations, but the routes to obtaining these are different. In essence, both methods identify thegoverning equations on the pore-scale, which are subject to closure conditions, and use this pore-scale problem to derive a system of equations over the macroscale. The MMS classically requiresthe medium to be strictly periodic whereas the volume-averaging method assumes that approxi-mate effective properties can be obtained from consideration of a periodic case [10]. The formalnature of the MMS enables higher-order corrections to the leading-order macroscale equations tobe determined. Conversely, the volume-averaging method can be more physically intuitive [see,for example, 10, 26–28] but no higher-order corrections nor precise quantification of errors can bedetermined. While classic homogenisation techniques that use the MMS assume strictly periodicmicrostructure [20, 21], Bruna & Chapman [4] and Dalwadi et al. [7, 8] exploit local periodicity3n the microstructure. Dalwadi et al. [7, 8] introduce one slowly varying change between adjacentcells, but the total area/volume of each cell remains constant throughout the material.In this paper, we derive a mathematical model to allow us to examine the role that the pore-structure plays on the macroscopic flow, transport, and sorption within a porous medium. Weconsider steady flow through a two-dimensional heterogeneous filter comprising an array of circu-lar obstacles arranged on a rectangular lattice. We allow for heterogeneities in the microstructurenot only by variations in the obstacle size [as in 7, 8] but also in the longitudinal spacing betweenobstacles, over the length of the porous medium. A novelty that this model offers over pre-existinghomogenised models is that the two routes through which the microstructure may vary affords amuch richer parameter space for exploration, including the ability to have a heterogeneous mi-crostructure while maintaining a constant porosity, and allowing for a specific study of anisotropy.On the complex fluid domain, the flow is modelled via the Stokes equations, with no-slip andno-penetration conditions on the surface of the obstacles, and the solute transport is modelled viaa time dependent advection-diffusion equation, with removal on the surface of the obstacle (§II).Following Bruna & Chapman [4] and Dalwadi et al. [7, 8], we exploit the local periodicity of thepore-scale geometry to homogenise the pore-scale problem via the use of the MMS (§III). Sincewe consider a microstructure that slowly varies in both the length of period and in size of solid ob-stacles, the total areas of each cell we consider varies along the length of the filter. The homogeni-sation method provides equations for fluid flow, solute transport and sorption that are uniformlyvalid throughout the heterogeneous porous medium with the effective permeability and diffusivitytensors determined numerically. These tensors are strongly anisotropic and show that porosityalone is an insufficient measure of the pore-structure (§IV). We use the homogenised model toinvestigate the effects of heterogeneous pore-structure in a simple one-dimensional steady-statefiltration problem (§V). Finally, we discuss the merits and limitations of the model (§VI).
II. MODEL PROBLEM
We consider the steady flow of fluid carrying a passive solute through a rigid porous medium intwo dimensions. The solute advects, diffuses and is removed via adsorption to the solid structure.The spatial coordinate is ˜ x := ˜ x e + ˜ x e , with ˜ x and ˜ x the dimensional longitudinal andtransverse coordinates, respectively, and e and e the longitudinal and transverse unit vectors,respectively. The fluid enters the porous medium uniformly through the inlet at the left ( ˜ x = 0 )4 IG. 1. We consider the flow of fluid carrying solute through a heterogeneous porous filter in two dimen-sions. The porous medium has length ˜ L and is formed of an array of circular obstacles of radius ˜ R (˜ x ) , eachlocated in the centre of a rectangular cell of transverse height ˜ h and longitudinal width A (˜ x )˜ h . The porousmedium is thus uniform in the transverse ( ˜ x ) direction but heterogeneous in the longitudinal ( ˜ x ) direc-tion. We assume that the spacing between obstacles is small relative to the length of the porous medium: (cid:15) := ˜ h/ ˜ L (cid:28) . and exits the porous medium through the outlet at the right ( ˜ x = ˜ L ) (Figure 1). We denotedimensional quantities with a tilde.The entire domain of the porous medium, denoted ˜Ω , comprises both the fluid and the solidstructure of the domain. The latter constitutes an array of circular solid obstacles, as discussed inmore detail below. We assume that the solute particles are small relative to the solid obstacles, andwe measure the local density of solute (amount of solute per volume of fluid) via the concentrationfield ˜ c ( ˜ x , ˜ t ) , where ˜ t is dimensional time. This concentration field is defined within the fluid phaseof the porous medium, denoted ˜Ω f . Note that we do not track any solute that is adsorbed to thesolid surface. The porous medium can be partitioned into an array of rectangular cells of fixedheight ˜ h and varying width A (˜ x )˜ h , where A is the dimensionless aspect ratio. Each cell containsa fixed and rigid circular obstacle of radius ˜ R (˜ x ) at its centre. The solid domain is the union5f these obstacles, and is denoted ˜Ω s := ˜Ω \ ˜Ω f . To prevent the obstacles from overlapping, werequire that R ≤ min( A ˜ h, ˜ h ) . This construction leads to a porous medium whose properties varyin the longitudinal direction but not in the transverse direction (see Figure 1). We further assumethat the porous medium is composed of a large number of obstacles in the longitudinal direction,which requires (cid:15) := ˜ h/ ˜ L (cid:28) with A = O (1) .We assume that the fluid is incompressible and Newtonian, and that the flow is steady anddominated by viscosity. As such, the fluid velocity ˜ v ( ˜ x ) and pressure ˜ p ( ˜ x ) satisfy the Stokesequations, subject to no-slip and no-penetration boundary conditions on the solid obstacles, − ˜ ∇ ˜ p + ˜ µ ˜ ∇ ˜ v = , ˜ x ∈ ˜Ω f , (1a) ˜ ∇ · ˜ v = 0 , ˜ x ∈ ˜Ω f , (1b) ˜ v = , ˜ x ∈ ∂ ˜Ω s , (1c)where ˜ µ is the dynamic viscosity of the fluid, ∂ ˜Ω s denotes the fluid–solid interface and ˜ ∇ is thegradient operator with respect to ˜ x .We model solute transport and sorption via the standard advection–diffusion equation with alinear, partially adsorbing condition at the fluid–solid interface: ∂ ˜ c∂ ˜ t = ˜ ∇ · (cid:16) ˜ D ˜ ∇ ˜ c − ˜ v ˜ c (cid:17) , ˜ x ∈ ˜Ω f , (2a) − ˜ γ ˜ c = ˜ n s · (cid:16) ˜ D ˜ ∇ ˜ c − ˜ v ˜ c (cid:17) , ˜ x ∈ ∂ ˜Ω s , (2b)where ˜ D is the coefficient of molecular diffusion, ˜ n s is the outward-facing unit normal to ∂ ˜Ω s ,and ˜ γ ≥ is the constant adsorption coefficient. Note that the second term on the right-handside of Equation (2b) vanishes due to Equation (1c). Further, note that ˜ γ = 0 corresponds to noadsorption and ˜ γ → ∞ corresponds to instantaneous adsorption, where the latter is equivalent toimposing ˜ c = 0 on ∂ ˜Ω s . We parameterise the fluid–solid interface ∂ ˜Ω s by the line ˜ f s ( ˜ x ) = 0 ,where the sign of ˜ f s is chosen so that ˜ n s ( ˜ x ) := ˜ ∇ ˜ f s (cid:12)(cid:12)(cid:12) ˜ ∇ ˜ f s (cid:12)(cid:12)(cid:12) , (3)is the outward-facing normal to the fluid domain.We make Equations (1)–(2) dimensionless via the scalings ˜ x = ˜ L ˆ x , ˜ v = ˜ V ˆ v , ˜ p = (cid:32) ˜ µ ˜ V (cid:15) ˜ L (cid:33) ˆ p, ˜ c = ˜ C ˆ c, and ˜ t = (cid:32) ˜ L ˜ D (cid:33) t, (4)6here ˜ V and ˜ C are the average inlet velocity and the average inlet concentration, respectively; ˆ x and t denote the dimensionless spatial and temporal coordinates, respectively; and ˆ v = ˆ v ( ˆ x ) , ˆ p = ˆ p ( ˆ x ) and ˆ c = ˆ c ( ˆ x , t ) denote the dimensionless velocity, pressure and concentrations fields,respectively. This pressure scale balances the macroscopic pressure gradient against viscous dis-sipation at the pore-scale, as is standard in lubrication problems. Employing the scalings in Equa-tion (4), the flow problem (Eqs. 1) becomes − ˆ ∇ ˆ p + (cid:15) ˆ ∇ ˆ v = , ˆ x ∈ ˆΩ f , (5a) ˆ ∇ · ˆ v = 0 , ˆ x ∈ ˆΩ f , (5b) ˆ v = , ˆ x ∈ ∂ ˆΩ s , (5c)where ˆ ∇ is the gradient operator with respect to ˆ x . Similarly, the transport problem (Eqs. 2)becomes ∂ ˆ c∂t = ˆ ∇ · (cid:16) ˆ ∇ ˆ c − Pe ˆ v ˆ c (cid:17) , ˆ x ∈ ˆΩ f , (6a) − (cid:15)γ ˆ c = ˆ n s · (cid:16) ˆ ∇ ˆ c − Pe ˆ v ˆ c (cid:17) , ˆ x ∈ ∂ ˆΩ s , (6b)where ˆ n s ( ˆ x ) , a function of ˆ x , is the outward-facing normal to ˆΩ f and the P´eclet number Pe :=˜ L ˜ V / ˜ D measures the rate of advective transport relative to that of diffusive transport and the dimen-sionless adsorption rate γ := ˜ γ ˜ L/ ( (cid:15) ˜ D ) measures the rate of adsorption relative to that of diffusivetransport. Note that γ ≡ Da /(cid:15) , where Da = ˜ γ/ ˜ V is the Damk¨ohler number of the first kind.As discussed in more detail below, the subsequent analysis requires that Pe , γ ∼ are constantsindependent of (cid:15) ; this represents a distinguished limit as highlighted below.Finally, the dimensionless fluid–solid interface becomes ˆ f s ( ˆ x ) = 0 and Equation (3) becomes ˆ n s ( ˆ x ) := ˆ ∇ ˆ f s (cid:12)(cid:12)(cid:12) ˆ ∇ ˆ f s (cid:12)(cid:12)(cid:12) . (7) III. HOMOGENISATION
Here, we approach the problem with homogenisation by the method of multiple scales (MMS).Homogenisation via the MMS is an asymptotic technique for domains that can be representedas the union of a large number of locally periodic cells. Following the MMS, we isolate andsolve the problem of flow and solute transport in an individual cell. We then construct a modelfor macroscopic flow and transport through the entire porous medium from the solution to these7ndividual cell problems via local averaging. The result is a system of equations that are uniformlyvalid for all ˆ x ∈ ˆΩ . A. Two spatial scales
Applying the MMS, we consider the spatial domain on two distinct length scales: the macroscale x := ˆ x , relative to which the length of the porous medium is precisely one, and the microscale y := ˆ x /(cid:15) , relative to which the transverse distance between the centres of the obstacles is pre-cisely one. The scalings defined in Equation (4) provide the macroscale and microscale geometriesshown in Figure 2, where A (˜ x ) = a ( x ) and ˜ R (˜ x ) = ˜ hR ( x ) . (8)Note that any arbitrary distribution of obstacles in the longitudinal directions ( i.e., arbitrary longi-tudinal heterogeneity) can be imposed by fixing the functions a ( x ) and R ( x ) .On the macroscale, the domain of the porous medium, x ∈ Ω , comprises a fluid domain Ω f and a complementary solid domain Ω s . On the microscale, we discretise the porous medium intorectangular cells of width a ( x ) and unit height centred around each solid obstacle; the origin ofthe microscale coordinate y is the centre of an arbitrary cell. Each cell ω ( x ) comprises a fluidphase ω f ( x ) and a solid obstacle ω s ( x ) := ω ( x ) \ ω f ( x ) . The fluid–solid interface ∂ω s ( x ) is the boundary of the obstacle. Each cell has four additional boundaries that separate it fromneighbouring cells. We denote the bottom and top boundaries ∂ω i = for i = 1 , , respectively,and the left and right boundaries ∂ω i || for i = 1 , , respectively (Figure 2). To proceed withhomogenisation, we need to construct and solve an isolated cell problem for arbitrary values of a and R .Following the MMS, we take x and y to be independent spatial parameters. We thereforerewrite all functions of ˆ x as functions of x and y : ˆ v ( ˆ x ) := v ( x , y ) , ˆ p ( ˆ x ) := p ( x , y ) , and ˆ c ( ˆ x , t ) := c ( x , y , t ) . Additionally, spatial derivatives become ∂∂ ˆ x i = ∂∂x i + 1 (cid:15) ∂∂y i , for i = 1 , . (9)For a given quantity Z ( x , y , t ) , there are two different averages of interest: the intrinsic (fluid)average (cid:104) Z (cid:105) ( x , t ) := 1 | ω f ( x ) | (cid:90) ω f ( x ) Z ( x , y , t )d S y , (10)8 IG. 2. The dimensionless problem on two spatial scales: the macroscale (left) and the cell problem onthe microscale (right). Here we have introduced the notation O x (1) to signify a quantity that is O (1) whenmeasured with respect to the x variable. where | ω f | = a − πR , and the volumetric average | ω ( x ) | (cid:90) ω ( x ) Z ( x , y , t )d S y , (11) | ω | = a , the latter of which is the average of choice in Dalwadi et al. [8]. Here, d S y := d y d y isan area element of the fluid region and the porosity φ is φ ( x ) = | ω f ( x ) || ω ( x ) | ≡ − πR ( x ) a ( x ) . (12)Thus, (cid:104) c (cid:105) is the amount of solute per unit fluid area within the porous medium, while φ (cid:104) c (cid:105) , thevolumetric average of the concentration, is the amount of solute per unit total area.We define the average velocity, pressure, and concentration as V ( ˆ x ) ≡ V ( x ) := (cid:104) v (cid:105) , P ( ˆ x ) ≡ P ( x ) := (cid:104) p (cid:105) , and C ( ˆ x , t ) ≡ C ( x , t ) := (cid:104) c (cid:105) , (13)respectively. Note that φ V is the standard Darcy flux.9 . Flow problem For a passive tracer as described in this paper, the flow problem (Eqs. 5a) does not depend on c .Using Equation (9), Equations (5) in an arbitrary cell become − (cid:18) (cid:15) ∇ x + 1 (cid:15) ∇ y (cid:19) p + (cid:18) (cid:15) ∇ x + ∇ x · ∇ y + ∇ y · ∇ x + 1 (cid:15) ∇ y (cid:19) v = , y ∈ ω f ( x ) , (14a) (cid:18) ∇ x + 1 (cid:15) ∇ y (cid:19) · v = 0 , y ∈ ω f ( x ) , (14b) v = , y ∈ ∂ω s ( x ) , (14c)where ∇ x and ∇ y are the gradient operators with respect to the coordinates x and y , respectively.Note that, for clarity of presentation in what follows, we have divided by (cid:15) during the derivationof Equation (14a). To close this cell problem, we must also impose boundary conditions at theexternal boundaries of the cell. Recall that the MMS requires periodicity of v , p and c at leadingorder ( i.e., local periodicity). To achieve this, the microscale geometry must also be locally peri-odic. In the transverse direction, each cell is geometrically identical to adjacent cells by definition;in the longitudinal direction, we must constrain the variation in R and a so that, at leading order,each cell is the same as adjacent cells. Thus, we enforce that the difference in R ( x ) and a ( x ) between two consecutive cells is O x ( (cid:15) ) ≡ O y ( (cid:15) ) , where O (cid:63) (1) means O (1) when (cid:63) is strictlyorder one. In other words, the spacing scales with (cid:15) when measured with respect to the x variableand scales with (cid:15) when measured with respect to the y variable. This is equivalent to saying that | ∇ x R | = O x (1) and | ∇ x a | = O x (1) ( i.e., that these quantities only vary on the macroscale). En-forcing periodicity of all quantities at both the top and bottom ( ∂ω = ( x ) := ∂ω ( x ) ∪ ∂ω ( x ) )and left and right ( ∂ω || ( x ) := ∂ω || ( x ) ∪ ∂ω || ( x ) ) cell boundaries, respectively, leads to v and p periodic on y ∈ ∂ω = ( x ) and ∂ω || ( x ) . (14d)We now seek an asymptotic solution to Equations (14) by expanding v and p in powers of (cid:15) : v ( x , y ) = v (0) ( x , y ) + (cid:15) v (1) ( x , y ) + · · · as (cid:15) → , (15a) p ( x , y ) = p (0) ( x , y ) + (cid:15)p (1) ( x , y ) + · · · as (cid:15) → . (15b)Considering terms of O y (1 /(cid:15) ) in Equations (14a) gives ∇ y p (0) = , (16)10rom which we conclude the standard result that, at leading order, the pressure is uniform on themicroscale: p (0) = p (0) ( x ) [8, 21].Considering terms of O y (1 /(cid:15) ) in Equation (14) gives − ∇ x p (0) − ∇ y p (1) + ∇ y v (0) = , y ∈ ω f ( x ) , (17a) ∇ y · v (0) = 0 , y ∈ ω f ( x ) , (17b) v (0) = , y ∈ ∂ω s ( x ) , (17c)with v (0) and p (1) periodic on y ∈ ∂ω = ( x ) and ∂ω || ( x ) . (17d)The form of Equations (17) suggest that we can scale ∇ x p (0) out of the problem via the substitu-tions v (0) = − K ( x , y ) · ∇ x p , (18a) p (1) = − Π ( x , y ) · ∇ x p + ˘ p ( x ) , (18b)where ˘ p ( x ) is a scalar function, K ( x , y ) is a tensor function, and Π ( x , y ) is a vector function.Using Equations (18) and the fact that p (0) is independent of y , Equations (17) become (cid:0) I − ∇ y ⊗ Π + ∇ y K (cid:1) · ∇ x p (0) = , y ∈ ω f ( x ) , (19a) ( ∇ y · K ) · ∇ x p (0) = , y ∈ ω f ( x ) , (19b) K · ∇ x p (0) = , y ∈ ∂ω s ( x ) , (19c)where I is the identity tensor. Equations (19) must hold for arbitrary ∇ x p (0) , hence K ( x , y ) and Π ( x , y ) must satisfy the system I − ∇ y ⊗ Π + ∇ y K = , y ∈ ω f ( x ) , (20a) ∇ y · K = , y ∈ ω f ( x ) , (20b) K = , y ∈ ∂ω s ( x ) , (20c)with K ij and Π i periodic on y ∈ ∂ω = ( x ) and ∂ω || ( x ) , (20d)and note that ( ∇ y ⊗ Π ) ij = ∂ Π j ∂y i and ( ∇ y · K ) i = (cid:88) j =1 ∂ K ji ∂y j . (20e)11quations (20) must then be solved numerically for each desired cell geometry ( i.e., pairs of a and R ). Note that Equations (20) do not depend on ∇ x p (0) .To derive a macroscale relationship between velocity and pressure from Equation (18a), weexpand the averaged quantities defined in Equation (13) in powers of (cid:15) : V ( ˆ x ) = V (0) ( ˆ x ) + (cid:15) V (1) ( ˆ x ) + · · · as (cid:15) → , (21a) P ( ˆ x ) = P (0) ( ˆ x ) + (cid:15)P (1) ( ˆ x ) + · · · as (cid:15) → . (21b)Note that P (0) ( ˆ x ) = (cid:104) p (0) ( x ) (cid:105) ≡ p (0) ( x ) , (22)since p (0) is independent of y . We then take the intrinsic average of Equation (18a) to determinethat the leading-order macroscale velocity depends on gradients in the leading-order macroscalepressure according to Darcy’s law: φ V (0) = − K ( φ, R ) · ˆ ∇ P (0) , (23a)where we have introduced the effective macroscale permeability tensor K ( φ, R ) := φ (cid:104) K (cid:105) , (23b)where φ and R are known functions of ˆ x . Prescribing both φ and R determines a via Equa-tion (12). Being averaged in y , Equation (23a) depends on ˆ x = x only and we have thereforereplaced ∇ x with ˆ ∇ . The symmetry of the geometry and the boundary conditions ( i.e., reflec-tional symmetry along both the y and y axes) suggests that K is diagonal in general, and that itreduces to a scalar multiple of I for a = 1 .Equation (23a) provides two equations for three unknowns. To develop another constraintin terms of V (0) and P (0) , we take the intrinsic average of Equation (17b). We manipulate theintrinsic average of Equation (17b) using a suitable transport theorem (see Appendix A) to arriveat ˆ ∇ · ( φ V (0) ) = 0 , (23c)which closes the system defined in Equations (23a). The details of this derivation are shown inAppendix B.To evaluate K , we solve Equations (20) in COMSOL Multiphysics ® using the ‘Laminar Flow(spf)’ interface (‘Fluid Flow’ → ‘Single Phase Flow’ → ’Laminar flow (spf)’). The domain isdiscretised using the ‘Physics-controlled mesh’ with the element size set to ‘Extremely fine’. In§IV we graphically present K ( φ, R ) and discuss its implications.12 . Transport problem We now perform a similar procedure for the solute-transport problem (Eqs. 6) which, under thespatial transformations (Eq. 9), become (cid:15) ∂c∂t = ( (cid:15) ∇ x + ∇ y ) · (cid:20)(cid:18) ∇ x + 1 (cid:15) ∇ y (cid:19) c − v c (cid:21) , y ∈ ω f ( x ) , (24a) − (cid:15)γc + O y ( (cid:15) ) = ( n ys + (cid:15) ∇ x R ) · (cid:20)(cid:18) ∇ x + 1 (cid:15) ∇ y (cid:19) c − v c (cid:21) , y ∈ ∂ω s ( x ) , (24b)with v , c, periodic on y ∈ ∂ω = and ∂ω || ( x ) . (24c)Note that we have extended the definition of the fluid–solid interface as introduced in Equation (7)to ˆ f s ( ˆ x ) = f s ( x , y ) , which vanishes on the interface: f s ( x , y ) ≡ R ( x ) − | y | = 0 , (24d)where | y | := (cid:112) y + y . Using Equation (7) and Equation (9), the outward-facing unit normal n s ( x , y ) = ˆ n s ( ˆ x ) becomes n s ( x , y ) = (cid:0) ∇ x + (cid:15) ∇ y (cid:1) f s (cid:12)(cid:12)(cid:0) ∇ x + (cid:15) ∇ y (cid:1) f s (cid:12)(cid:12) = n ys + (cid:15) ∇ x R | n ys + (cid:15) ∇ x R | , (24e)where n ys := − y / | y | is the outward-facing unit normal to the circular obstacle in a given cell [8].Note also that, for clarity of presentation in what follows, we have multiplied by (cid:15) when derivingEquation (24a) from Equation (6a).We now consider an expansion of the concentration field of the form c ( x , y , t ) = c (0) ( x , y , t ) + (cid:15)c (1) ( x , y , t ) + (cid:15) c (2) ( x , y , t ) + · · · as (cid:15) → . (25)Note that we take Pe , γ ∼ to be constants independent of (cid:15) ; this corresponds to a distinguishedlimit where all the transport mechanisms balance over the macroscale ( cf., Equation (40a)). Con-sidering Equation (24) at leading order — that is, O (1 /(cid:15) ) — we obtain ∇ y c (0) = 0 , y ∈ ω f ( x ) , (26a) n ys · ∇ y c (0) = 0 , y ∈ ∂ω s ( x ) , (26b) c (0) periodic on y ∈ ∂ω = ( x ) and ∂ω || ( x ) . (26c)13y inspection, c (0) ( x , t ) is a non-trivial solution to this system. By linearity, this solution is uniqueand independent of y ( i.e., ∇ y c (0) ≡ ).Considering Equations (24) at O y (1) , we obtain ∇ y c (1) = 0 , y ∈ ω f ( x ) , (27a) n ys · ∇ y c (1) = − n ys · ∇ x c (0) , y ∈ ∂ω s ( x ) , (27b) c (1) periodic on y ∈ ∂ω = and ∂ω || ( x ) , (27c)where we have used ∇ y c (0) ≡ , microscale incompressibility (Eq. 17b), and the no-slip and no-penetration conditions (Eq. 17c) on the solid surface. The form of Equations (27) suggest that wecan scale ∇ x c (0) out of the problem via the substitution c (1) ( x , y , t ) = − Γ ( x , y ) · ∇ x c (0) ( x , t ) + ˘ c ( x , t ) , (28)where ˘ c is a scalar function and Γ is a vector function. Substituting Equation (28) into Equation(27), we obtain a system of equations for Γ : ∇ y Γ i = 0 , y ∈ ω f ( x ) , (29a) n ys · ∇ y Γ i = n i , y ∈ ∂ω s ( x ) , (29b) Γ i periodic on y ∈ ∂ω = and ∂ω || ( x ) , (29c)where i ∈ { , } and n i := ( n ys ) i . Equations (29) must then be solved numerically for each i andeach desired cell geometry ( i.e., pairs of a and R ). Note that Equations (29) do not depend on ∇ x c (0) and are only functions of the cell geometry.The goal of this analysis remains to determine a macroscale equation for the concentration.Since there are no macroscopic transport mechanisms present at this order, there is not enoughinformation to determine a macroscale governing equation for the concentration. Hence, we musttherefore proceed to the next order in Equations (24), which yields ∂c (0) ∂t = ∇ y · A + ∇ x · B , y ∈ ω f ( x ) , (30a) − γc (0) = n ys · A + B · ∇ x R y ∈ ∂ω s ( x ) , (30b) c (2) periodic on y ∈ ∂ω = ( x ) and ∂ω || ( x ) , (30c)where A := (cid:0) ∇ y c (2) + ∇ x c (1) (cid:1) − Pe (cid:0) v (0) c (1) − v (1) c (0) (cid:1) , (30d)14nd B := ∇ y c (1) + ∇ x c (0) − Pe v (0) c (0) . (30e)Note that B · ∇ x R ≡ B ∂R/∂x since ∂R/∂x ≡ . Taking the intrinsic average of Equation(30a) and multiplying by aφ gives aφ ∂C (0) ∂t = (cid:90) ω f ∇ y · A d S y + (cid:90) ω f ∇ x · B d S y , (31)where C (0) is the leading-order term in the expansion of C , defined in Equation (13): C ( ˆ x , t ) = C (0) ( x , t ) + (cid:15)C (1) ( x , t ) + (cid:15) C (2) ( x , t ) + · · · as (cid:15) → . (32)Note that C (0) = c (0) by the definition of C (Eq. 13), since c (0) is independent of y . Applying thedivergence theorem to the first integral on the right-hand side of Equation (31) yields (cid:90) ω f ∇ y · A d S y = (cid:90) ∂ω n ys · A d s y = (cid:90) ∂ω s n ys · A d s y + (cid:90) ∂ω (cid:3) n y (cid:3) · A d s y , (33)where d s y signifies a scalar line integral and n y (cid:3) is the outward-facing unit normal to ∂ω (cid:3) := ∂ω = ( x ) ∪ ∂ω || ( x ) . The last term of the right-hand side of Equation (33) vanishes due to period-icity. Using Equation (30b), (cid:90) ω f ∇ y · A d S y = − (cid:90) ∂ω s B ∂R∂x d s y − (cid:90) ∂ω s γc (0) d s y , (34a) = − (cid:90) ∂ω s B ∂R∂x d s y − aγ (1 − φ ) R C (0) , (34b)where the last equality is a consequence of Equation (12).To manipulate the second integral on the right-hand side of Equation (31), we apply the trans-port theorem (Eq. A11): (cid:90) ω f ∇ x · B d S y = ∇ x · (cid:90) ω f B d S y − ¯ B ∂a∂x + (cid:90) ∂ω s B ∂R∂x d s y = a (cid:34) ∇ x · (cid:32) a (cid:90) ω f B d S y (cid:33)(cid:35) + 1 a ∂a∂x (cid:90) ω f B d S y − ¯ B ∂a∂x + (cid:90) ∂ω s B ∂R∂x d s y , (35)where ¯ B := (cid:90) − B d y . (36)The second and third terms on the right-hand side of Equation (35) cancel, since (cid:90) ω f B d S y = (cid:90) a − a (cid:90) − B d y d y = (cid:90) a − a ¯ B d y = a ¯ B , (37)15ecause ¯ B is independent of y (see Appendix C). Thus, combining Equation (37) with Equation(35), we obtain (cid:90) ω f ∇ x · B d S y = a (cid:34) ∇ x · (cid:32) a (cid:90) ω f B d S y (cid:33)(cid:35) + (cid:90) ∂ω s B ∂R∂x d s y . (38)Finally, combining Equations (31), (34b) and (38) we obtain φ ∂C (0) ∂t = (cid:34) ∇ x · (cid:32) a (cid:90) ω f ∇ y c (1) + ∇ x c (0) − Pe v (0) c (0) d S y (cid:33)(cid:35) − γ (1 − φ ) R C (0) , (39)where we note that the final term on the right-hand side of Equation (38) cancels with the firstterm on the right-hand side of Equation (34b). Hence, using the definitions of c (1) (Eq. 28), C (0) (Eq. 32) and V (0) (Eq. 21a), Equation (39) can be written as φ ∂C (0) ∂t = ˆ ∇ · (cid:104) φ D ( φ, R ) · ˆ ∇ C (0) − Pe V (0) φC (0) (cid:105) − F ( φ, R ) C (0) , (40a)where the effective diffusivity tensor D ( φ, R ) is D ( φ, R ) := (cid:28) − aφ (cid:90) ω f ( x ) ∇ y ⊗ Γ ; d S y , (40b)and where the effective adsorption rate F ( φ, R ) is F ( φ, R ) = 2 γ (1 − φ ) R = 2 πRγa . (40c)Note that φ and R are functions of x and prescribing both φ and R determines a via Equation (12).Being averaged in y , Equation (40a) now depends on x and t only and we have therefore replaced ∇ x with ˆ ∇ . In the limit a = 1 , Equations (40) become the same system as Equation (3.22) inDalwadi et al. [8] in two dimensions ( i.e., d = 2 ), but written in terms of the intrinsic average C (0) ( ˆ x , t ) ≡ (cid:104) c (0) (cid:105) rather than the volumetric average φ (cid:104) c (0) (cid:105) ≡ φC (0) .To evaluate D , we solve Equations (29) numerically using COMSOL Multiphysics ® usingthe ‘Laplace Equation (lpeq)’ interface (‘Classical PDEs’ → ‘Mathematics branch’ → ‘LaplaceEquation (lpeq)’). As for the flow problem, the domain is discretised using the ‘Physics-controlledmesh’ with the element size set to‘Extremely fine’.Equation (40a) describes macroscopic transport by advection and diffusion in a porous mediumwith chemical sorption such that φ D · ˆ ∇ C (0) is the diffusive flux per unit area of porous mediumand D · ˆ ∇ C (0) is the diffusive flux per unit area of fluid. Thus in §IV we graphically present φ D ( φ, R ) and discuss its implications. 16 IG. 3. The porosity φ of a rectangular cell increases with aspect ratio a and decreases with obstacle radius R according to Equation (12). (a) a versus φ for R ∈ { . , . , . , . , . , . , . , . , . , . } (dark to light). The attainable region of the φ - a plane (shaded grey) is bounded below by φ min given byEquation (12) with a = 2 R (dot-dashed line). The cell geometry for two distinct points with R = 0 . is shown pictorially. (b) a versus R for φ ∈ { . , . , . , . , . , . , . , . , . , . } (dark to light). The attainable region of the R – a plane (shaded grey) is also bounded below by a min := 2 R (dot-dashed line). Note that the smallest attainable φ for any R , a combination is φ min ( R = 1 /
2) = 1 − π/ . IV. EFFECTIVE PERMEABILITY AND DIFFUSIVITY TENSORS
In this section, we explore the impact of microstructure on macroscopic flow and transport byanalysing the effective permeability and effective net diffusivity tensors, K and φ D , respectively.Recall that these tensors depend on microstructure via a , R , and φ , any two of which are indepen-dent (Figure 3). We therefore have one additional degree of microstructural freedom relative toDalwadi et al. [8] and this allows us to explore the anisotropy in the system. We explore the effectof the a , R and φ parameter space on K and φ D in Figures 4 and 5, respectively. The effectivediffusivity D is shown for reference in Figure 8 (top and middle row; Appendix D).Increasing φ at fixed R is achieved by increasing a (Figure 3 (a)), such that the obstacles movefurther apart in the longitudinal direction only; as a result, K , K , φD and φD all increase(Figure 4 (a) and (c) and Figure 5 (a) and (c)). As φ → ( a → ∞ ), both K and K diverge17 IG. 4. The effective longitudinal permeability K (top row), transverse permeability K (middlerow) and the permeability–anisotropy ratio K /K (bottom row) depend strongly on microstructure.Left column: K , K and K /K against φ for fixed values of R ∈ { . , . , . , . , . } , with a varying according to Equation (12). Right column: the same quantities against R for fixed values of φ ∈ { . , . , . , . , . } , with a varying according to Equation (12). For a given value of R , the min-imum porosity φ min ( R ) is given by Equation (12) with a = 2 R . Note that K is non-zero at φ min ( R ) forall values of φ (dot-dashed curve, top row), whereas K vanishes at φ min ((c) and (e) red vertical asymp-totes). There exists a smallest possible R for any given φ ((d) and (f) red vertical asymptotes). In all cases, K and K are as defined in Equation (23b) and calculated using COMSOL Multiphysics ® . The effectivepermeability is isotropic when a ≡ (solid black curves; Dalwadi et al. [8]). The limit R → and a → corresponds to a set of parallel but disconnected channels with unit transverse width, for which φ → , K → , and K → / (black diamonds in top row).
18s the resistance to flow vanishes (Figure 4 (a) and (c)), and both φD and φD tend to 1 asmolecular diffusion becomes unobstructed (Figure 5 (a) and (c)). As φ → φ min ( R ) ( a → R )at fixed R , the obstacles move closer together in the longitudinal direction and the pore spacebecomes disconnected in the transverse direction, so that K and φD vanish; K and φD areminimised but do not vanish. Further taking R → , the longitudinal problem reduces to a set ofdisconnected parallel channels of unit transverse width, for which K = 1 / (Figure 4 top row,black diamond).Increasing R at fixed φ is similarly achieved by increasing a (Figure 3 (b)), in which casethe transverse channels between obstacles grow wider while the longitudinal channels betweenobstacles grow narrower. As a result, K and φD decrease while K and φD increase. As R → / at fixed φ , the longitudinal channels close and K and φD vanish, but the transversechannels become wider and K and φD are maximised. The longitudinal permeability, K , isweakly non-monotonic in R for larger values of φ (Figure 4 (b)), which means that the longitudinalpermeability of a high porosity filter can be maximised for a given φ by appropriately varying R and a .When a ≡ , equivalent to the case considered in Dalwadi et al. [8], K and φ D becomeisotropic. Increasing φ corresponds to decreasing R , in which case both the longitudinal and trans-verse spacing between obstacles decreases (Figure 3) which decreases K = K and φD = φD (Figures 4 and 5, solid black lines). For a (cid:54) = 1 , our microstructure is inherently anistropic( K (cid:54) = K , φD (cid:54) = φD ). For a < , the longitudinal channels are wider than the trans-verse channels, such that K /K < and D /D < and both ratios vanish as a → R ( φ → φ min where φ min is given by Equation (12); Figure 4 (e) and (f) and Figure 5 (e) and (f),respectively). For a > , the longitudinal channels are narrower than the transverse channels,such that K /K > and D /D > . The permeability–anisotropy ratio, K /K , increasesmonotonically with both φ and R , and diverges as φ → at fixed R ( K diverges faster than K because the obstacles never get further apart in the transverse direction) and as R → / atfixed φ ( K vanishes; Figure 4 (e) and (f)). The diffusivity–anisotropy ratio, D /D , increasesmonotonically with R for all φ ∈ ( φ min , , diverging as R → / (Figure 5 (f)). For fixed R , thisratio increases monotonically with φ for a ≤ , is equal to unity for a = 1 (isotropic geometry),and must approach unity as φ → ( a → ∞ ; unobstructed molecular diffusion; Figure 5(e)).These bounds require that D /D has an intermediate maximum in φ (or in a ) at fixed R , theamplitude of which diverges as R → / . Specifically, the non-monotonicity in the ratio D /D IG. 5. The net longitudinal diffusivity φD (top row), net transverse diffusivity φD (middle row)and the diffusivity–anisotropy ratio D /D (bottom row) depend strongly on microstructure. Left col-umn: φD , φD and D /D against φ for fixed values of R ∈ { . , . , . , . , . } , with a varying according to Equation (12). Right column: the same quantities against R for fixed values of φ ∈ { . , . , . , . , . } , with a varying according to Equation (12). For a given value of R , theminimum porosity φ min ( R ) is given by Equation (12) with a = 2 R . Note that φD is non-zero at φ min ( R ) for all values of φ (dot-dashed curve, top row), whereas φD vanishes at φ min . In all cases, D and D are as defined in Equation (40b) and calculated using COMSOL Multiphysics ® . Both D and D tendto 1 as φ → , which corresponds to the limit of free-space diffusion. The net diffusivity, φ D , is isotropicwhen a = 1 (solid black curves), in agreement with the results presented in Dalwadi et al. [8]. φD and φD . For φ = φ min ( a = 2 R ) there isno transverse connectivity thus separating the obstacles slightly (a small increase in a ) leads toa sharp increase in φD but only a slight increase in φD since the longitudinal connectivityis unchanged and most longitudinal mixing occurs in the longitudinal channels. Conversely, as φ → ( a → ∞ ) the longitudinal spacing between obstacles diverges which means that φD isvery sensitive to changes in a as longitudinal mixing occurs predominantly between longitudinallyadjacent obstacles (in the transverse channels), thus φD approaches unity rapidly. However, sig-nificant transverse connectivity is preserved for large a so increasing a further has minimal effecton φD since most transverse mixing occurs in the transverse channels in this limit. Note that,the longitudinal diffusivity D is non-monotonic in φ for each R as a varies (Appendix D, Figure8). V. SIMPLE ONE-DIMENSIONAL FILTER
We now use the homogenised model to understand the effect of microstructure and P´eclet num-ber on filter efficiency in the context of a simple one-dimensional steady-state filtration problem.Specifically, we consider Equations (23) and (40) at steady state, with imposed flux and concen-tration at the inlet, φ V (0) = e at ˆ x = 0 , (41a) C (0) = 1 at ˆ x = 0 , (41b)and passive outflow at the outlet, ∂C (0) ∂ ˆ x = 0 at ˆ x = 1 . (41c)Since these boundary conditions (Eq. 41a) are compatible with unidirectional flow, we take V (0) ( ˆ x ) = V (0)1 (ˆ x ) e and C (0) ( ˆ x , t ) = C (0) (ˆ x ) . Thus, Equation (23c) leads to ddˆ x (cid:16) φV (0)1 (cid:17) = 0 , (42)which, on application of the inlet condition (Eqs. 41a), gives the macroscale flux φV (0)1 ≡ for all ˆ x . Hence, the homogenised governing equation for the steady concentration distribution C (0) (ˆ x ) (Eq. 40a) becomes ddˆ x (cid:20) φD ( φ, R ) d C (0) dˆ x − Pe C (0) (cid:21) = F ( φ, R ) C (0) for ˆ x ∈ (0 , , (43)21here F ( φ, R ) is defined in Equation (40c). We solve Equation (43) subject to Equations (41b)and (41c) numerically using a finite-difference scheme. We discretise the interval [0, 1] using auniform mesh of size ∆ x = 1 /N and we approximate derivatives using a second-order accuratecentral-difference formula. The results presented below were obtained using N = 500 . Below,we consider filters with varying porosity and filters with uniform porosity, but with heterogeneousmicrostructure in both cases. A. Porosity gradients
We first consider filters with varying porosity. Recall that a given porosity φ may be achievedin two different ways: by fixing a (ˆ x ) and varying R (ˆ x ) , as considered in Dalwadi et al. [8] for a ≡ ; or by fixing R (ˆ x ) and varying a (ˆ x ) . We consider these options in Figure 6, for thesame three porosity fields in both cases: linearly increasing with ˆ x , uniform in ˆ x , or linearlydecreasing with ˆ x . Specifically, we take φ (ˆ x ) = φ + m φ (ˆ x − . , where φ is the averageporosity (also the mid-point porosity) and m φ is the porosity gradient. We take φ = 0 . and m φ = 0 . (increasing in ˆ x ), m φ = 0 (uniform in ˆ x ), or m φ = − . (decreasing in ˆ x ). Note that, | m φ | defines the filter microstructure and sgn( m φ ) is simply the orientation of the filter.When varying φ by varying R at fixed a ≡ , as considered by Dalwadi et al. [8], the sign of theporosity gradient has a modest impact on the concentration distribution within the filter: m φ > leads to a steeper gradient in C (0) near the inlet and a shallower gradient in C (0) near the outlet,whereas m φ < leads to a more uniform gradient in C (0) throughout the filter (Figure 6(a)).However, the outlet concentration C (0) out := C (0) (1) is remarkably insensitive to m φ . The outletconcentration is slightly lower for m φ < , and this slight difference decreases as Pe increases(Figure 6(b)). As Pe increases, advection becomes stronger causing more solute to be sweptthrough the filter; as a result, C (0) (ˆ x ) increases with Pe for all ˆ x , and C (0)out more than doubles as Pe increases from 0 to 10. The case when a ≡ is considered in more detail in Dalwadi et al. [8].Varying φ by varying a at fixed R ≡ . leads to qualitatively similar results, but C (0) ( x ) and C (0) out are more sensitive to m φ (Figure 6, bottom row). For all Pe , attaining a desired porosity gradientvia varying R leads to a more efficient filter than varying a , in the sense that C (0) out is lower for thesame φ ( x ) . 22 IG. 6. Steady-state concentration field C (0) (ˆ x ) for Pe ∈ { , , } (left column) and outlet concentration C (0)out := C (0) (1) as a function of Pe (right column). Top row: a ≡ and φ (ˆ x ) = 0 . m φ (ˆ x − . .Bottom row: R ≡ . and φ (ˆ x ) = 0 . m φ (ˆ x − . . In both cases, m φ = 0 . (solid), m φ = 0 (dotted),and m φ = − . (dashed). B. Microstructural gradients with uniform porosity
We now consider filters with uniform porosity but gradients in microstructure. We thereforefix φ and simultaneously vary R and a with ˆ x , recalling that a is related to R via Equation (12).We consider two types of variation: an imposed gradient in R with a varying to maintain constant φ (via Equation (12)) or an imposed gradient in a with R varying to maintain constant φ (viaEquation (12)). We consider these options in Figure 7.23 IG. 7. Steady-state concentration field C (0) (ˆ x ) for φ ∈ { . , . , . } (left column), and outletconcentration C (0)out := C (0) (1) as a function of φ (right column). Top row: Pe = 1 and R (ˆ x ) =0 .
36 + m R (ˆ x − . for m R = 0 . (solid), m R = 0 (dotted), m R = − . (dashed). The aspectratio a varies following Equation (12). For R = 0 . and m R = ± . the minimum attainable porosityis − . π ≈ . ((b); black dots). Bottom row: Pe = 1 and a (ˆ x ) = 1 .
34 + m a (ˆ x − . for m a = 1 . (solid), m a = 0 (dotted), and m a = − . (dashed). The radius R varies following Equation (12).For a = 1 . and m a = ± . the minimum attainable porosity is − . π ≈ . ((d); black dots). We first consider R (ˆ x ) = R + m R (ˆ x − . , where R is the average obstacle radius (also themid-point) radius and m R is the gradient (Figure 7, top row). We take R = 0 . and m R = 0 . (increasing in ˆ x ), m R = 0 (uniform in ˆ x ) or m R = − . (decreasing in ˆ x ). When varying R linearly, the sign of m R has a modest impact on the concentration distribution: for any uniform24orosity, m R > leads to a shallower gradient in C (0) and a higher concentration at every pointwithin the filter, including the outlet. Thus, for any uniform porosity, m R < is always a moreefficient filter than m R > .We next consider a (ˆ x ) = a + m a (ˆ x − . , where a is the average cell width (also thecell width at the mid-point) and m a is the gradient of a over ˆ x (Figure 7 bottom row). We take a = 1 . and m a = 1 . (increasing in ˆ x ), m a = 0 (constant in ˆ x ), or m a = − . (decreasingin ˆ x ). From Equation (12) for fixed φ , it can be seen that a ∝ R , thus, a decrease in a mustbe mirrored by a decrease in R to maintain a uniform φ . Thus, we expect the same qualitativebehaviour for a linear gradient in R (Figure 7, top row) as for a linear gradient in a (Figure 7,bottom row). We find that m a < leads to a more efficient filter for all φ .Note that, for any φ ∈ (1 − . π, — that is, the range of porosities attainable for imposedlinear gradients in both R and a (see Figure 7 caption) — prescribing a and taking m a < predictsmost efficient filter considered here (comparing Figure 7 (b) and Figure 7 (d)). Similarly, for large φ (cid:38) . prescribing a and taking m a > predicts the least efficient filter, whereas, for small φ (cid:46) . prescribing R and taking m R > predicts the least efficient filter. VI. CONCLUSIONS
We have systematically derived a macroscopic model for flow, transport and sorption duringsteady flow in a two-dimensional heterogeneous and anisotropic porous medium by generalis-ing standard homogenisation theory for locally periodic systems [8]. The heterogeneity origi-nates from slowly varying obstacle size and/or obstacle spacing along the length of the porousmedium, the latter also induces strong anisotropy within the problem. For the flow problem, weobtain Darcy’s law with an anisotropic effective permeability tensor, determined numerically, andfor the solute concentration problem we obtain an advection–diffusion–reaction equation with ananisotropic effective diffusivity tensor, also determined numerically. The effective permeability,effective diffusivity and the removal term are functions of the porosity, obstacle radius and ob-stacle spacing, two of these being free choices and the third prescribed by simple cell geometry.This work illustrates and quantifies how the permeability and diffusivity of a porous medium notonly depend on the porosity of the medium, but also depend strongly on the microstructure of themedium.The model derived in §III is two-dimensional and a direct physical analogue would be a quasi-25wo-dimensional filter composed of solid circular pillars, centered on a rectangular grid, that aresufficiently tall that the boundary effects of the fluid at the top and bottom of this medium may beneglected. This is a simple but appropriate model for non-woven fibrous filters which form a majorpart of the filtration industry ( e.g., those in air purifiers and vacuum cleaners), magnetic separationfilters that are composed of wire wool, and microfluidic devices containing tall micropillars. Thestrong anisotropy in the problem could be useful for filter design; it is achieved while maintainingthe circular shape of the obstacles and the principal directions of the permeability and diffusivitytensors are fixed as the longitudinal and transverse directions. Furthering our understanding of theimpacts of microstructural heterogeneity and anisotropy in general, is of use to many other areasof research including hydrology and biology ( e.g.,
16, 24).The homogenisation procedure we used allows for slowly varying changes to the cell surround-ing each circle that comprises the filter. This means that the total area of each individual cell maydiffer between cells. This is a new aspect to homogenisation and we have carefully derived a trans-port theorem to account for how these microstructural changes affect the macroscale transport.Using this transport theorem we have shown that macroscale incompressibility is preserved (thedivergence of the Darcy flux vanishes) and that this is independent of the individual cell size. Thetwo degrees of microstructural freedom (varying obstacle size and spacing) enable us to considera wide range of heterogeneities on the microscale, for example, to maintain a uniform porositywhile systematically varying the microstructure. These macroscale equations are computationallyinexpensive to solve, allowing for optimisation of parameters through large sweeps, which wouldnot be possible with direct numerical simulations.In the analysis we conducted, we considered a regime in which advection, diffusion and re-moval balance in the final macroscale equations, allowing sub-limits to be taken without repetitionof the interim analysis. Note that if the P´eclet number is large (comparable in magnitude to theratio of filter size to obstacle spacing) then dispersive terms — that is, terms that are proportionalto the product of velocity and concentration gradient — arise in the final homogenised equations.This is the subject of future work.We considered a simple model problem for a one-dimensional filter with chemical adsorptionat steady state. Measuring efficiency as the amount of solute removed by the filter per unit time,we found that negative porosity gradients lead to a more efficient filter than positive porositygradients or filters of uniform porosity. Further, for a fixed porosity, decreasing obstacle size ordecreasing obstacle spacing lead to more efficient filters than their respective constant or increasing26ounterparts. For a given porosity decreasing the obstacle spacing linearly leads to a more efficientfilter than linearly decreasing the obstacle radius.While we have defined efficiency to mean instantaneous performance, there are further consid-erations to a filter’s efficiency. Factors such as manufacturing costs, filter lifetime and fluid fluxoutput may also need to be considered. For example, if the coating on the solid obstacles wasvery expensive then we may wish to minimise the amount of surface area of the solid obstacleswhile maximising performance. Further, we assumed that the solid surface never saturates withsolute. However, in practice, the number of active sites where the solute can attach to the solidwill decrease as solute adsorbs, which may reduce the efficiency. In this case, this effect maybe mitigated by ensuring that there are active sites throughout the full length of the filter so thatthe chance a solute particle comes into contact with an active site is maximised. The simple one-dimensional filter model we considered only predicts initial or instantaneous filter efficiency andwill therefore not predict the total amount of contaminant filtered out over the life span of a filterif properties were to change with time. However, the equations derived in this paper can readilybe generalised to describe such a case. All of these additional considerations to filter design leadto multiple optimisation problems, requiring large parameter sweeps, for which a computationallyinexpensive model, such as this, is vital.In our analysis we have assumed that the solute particles are negligibly small; for particlesthat are not negligible in size relative to the smallest distances between adjacent obstacles (chokepoints), we would also need to consider the effects of choking of the filter due to particle build-up.Avoidance of such filter blockages requires sufficiently wide longitudinal connectivity. Hence,a filter comprising obstacles whose radii increase with depth is desirable, since such a gradientallows for more build-up of solute on the solid obstacles near the inlet without choking the filter.This scenario was considered in Dalwadi et al. [7]. However, in our case, we also have the pos-sibility of varying the spacing between obstacles. This additional degree of freedom allows us torespect a positive gradient in the obstacle radii to mitigate the risks of blockages while also havingeither a negative gradient in the porosity or obstacle spacing to enable more efficient filters.While it has been shown that the effective diffusivity for a porous medium with obstacles ona uniform square grid was qualitatively similar to a porous medium with obstacles on a uniformhexagonal grid for all porosities [4], we expect that the addition of anisotropy to the hexagonalproblem, obtained by varying the longitudinal obstacle spacing, will cause the effective perme-abilities and diffusivites to diverge from those determined here. For example, in certain limits, the27exagonal problem reduces to a series of longitudinal channels while in other limits, the hexag-onal problem reduces to a series of transverse channels. Consequently in the latter limit, for thehexagonal structure, the longitudinal permeability and diffusivity must vanish, while for the rect-angular structure the longitudinal permeability and diffusivity remain non-zero for the all param-eter combinations. In general, the hexagonal structure of obstacles will mean that the longitudinalpermeability will be more sensitive to longitudinal obstacle spacing than it is for obstacles in arectangular structure. This is because with a hexagonal grid, altering the longitudinal spacing al-ters both the longitudinal and transverse distances between neighbouring obstacles, while for arectangular grid, altering the longitudinal spacing does not alter the transverse distance betweenobstacles. This illustrates that, when anisotropy is introduced into a problem, the microstructurebecomes more significant than for isotropic problems.It would be straightforward to generalise our approach to a three-dimensional porous mediumcomprising spherical obstacles centred on a cuboid grid that is homogeneous in two directions,but again allowing for arbitrary variation of both obstacle radius and obstacle spacing in thelongitudinal direction. We would expect that the results would be qualitatively similar to thetwo-dimensional problem considered here, however connectivity does not vanish when obstaclestouch. This would then mean that we have non-zero permeability and diffusivity in all directionsthroughout the entire parameter space.A final point to note is that the spacing between obstacles may change when a filter is subjectto an effective stress. By coupling the model presented here to a law that relates the spacingof the obstacles to the strain of the porous medium, we can derive homogenised equations for afilter undergoing longitudinal deformation. Modelling the filter as a series of circles on a varyinghexagonal grid will better describe granular materials and this is the focus of future work.In summary, the results presented in this manuscript form a comprehensive framework for de-scribing the transport and adsorption properties through heterogeneous porous media. The modelcan be used to answer questions on the filtration performance of such porous media as well asbeing well-equipped for the generalisation to more complicated scenarios.28 ppendix A: The Transport Theorem1. The Generalised Transport Theorem
Firstly, we present a generalised form of the transport theorem which allows us to interchange ∇ x with integration over a cell of arbitrary geometry. Consider the region α bounded by thesurface ∂α , and suppose that this region moves and/or deforms with time t . Denote the position ofpoints on ∂α ( t ) by y b ( t ) . The Reynolds Transport Theorem states that dd t (cid:90) α ( t ) ζ d V = (cid:90) α ( t ) ∂ ζ ∂t d V + (cid:90) ∂α ( t ) (cid:18) ∂ y b ∂t · n (cid:19) ζ d A, (A1)for an arbitrary vector field ζ ( ˆ x , t ) , where d V signifies a volume integral, d S signifies a surfaceintegral, n is the outward normal to ∂α and the time derivative ∂ y b /∂t can be identified as thelocal velocity of ∂α ( t ) .In Equation (A1), t plays the role of an arbitrary scalar parameter. In other words, Equa-tion (A1) remains valid if we suppose that the region moves and/or deforms according to someother scalar parameter ξ , in which case we have that dd ξ (cid:90) β ( ξ ) z d V = (cid:90) β ( ξ ) ∂ z ∂ξ d V + (cid:90) ∂β ( ξ ) (cid:18) ∂ y b ∂ξ · n (cid:19) z d A, (A2)for an arbitrary vector field z ( ˆ x , ξ ) , and where the domain β is a function of ξ Note that thederivative ∂ y b /∂ξ can no longer be identified as a velocity in the traditional sense.Now, consider several independent parameters as a vector ξ = ξ i e i , where here we are adoptingthe Einstein summation convention and e i is the unit normal in the i th direction. The correspondingdivergence with respect to this vector is then ∇ ξ · (cid:90) β ( ξ ) z d V = (cid:18) e i ∂∂ξ i (cid:19) · (cid:90) β ( ξ ) ( z j e j ) d V = ∂∂ξ i (cid:90) β ( ξ ) z i d V. (A3)Equation (A2) then indicated the following expression for the right-hand side of Equation (A3) ∂∂ξ i (cid:90) β ( ξ ) z i d V = (cid:90) β ( ξ ) ∂z i ∂ξ i d V + (cid:90) ∂β ( ξ ) (cid:18) ∂ y b ∂ξ i · n (cid:19) z i d A. (A4)Returning to vector notation, we can rewrite this result as ∇ ξ · (cid:90) β ( ξ ) z d V = (cid:90) β ( ξ ) ∇ ξ · z d V + (cid:90) ∂β ( ξ ) n · G · z d A, (A5a)29here G = ∂y bi ∂ξ j e i e j = (cid:18) ∂ y b ∂ ξ (cid:19) (cid:124) = ( ∇ ξ y b ) (cid:124) , (A5b)is the Jacobian of the dependence of ∂β on ξ . Thus, Equation (A5) defines the generalised trans-port theorem.
2. The transport theorem for cell geometry of varying a and R We use the generalised the Reynolds transport theorem (Eq. A5) leading to the generalisedtransport theorem for a vector field z ( x , y ) , over the region ω f , which states that ∇ x · (cid:90) ω f ( x ) z d S y = (cid:90) ω f ( x ) ∇ x · z d S y + (cid:90) ∂ω ( x ) n y · G · z d s y , (A6)where ∂ω ( x ) := ∂ω = ( x ) ∪ ∂ω || ( x ) ∪ ∂ω s ( x ) , n y is the outward-facing unit normal to ∂ω ( x ) ,and G := (cid:18) ∂ y b ∂ x (cid:19) (cid:124) ≡ ( ∇ x ⊗ y b ) (cid:124) , (A7)is the Jacobian of the dependence of ∂ω ( x ) on x where y b denotes the points on ∂ω ( x ) .The outward-facing normal to the fluid domain n y comprises n y (cid:3) and n ys ; n y (cid:3) can be furtherdecomposed into n y = and n y || which are the outward-facing normals to the boundaries ∂ω = ( x ) and ∂ω || ( x ) respectively. Note that n ys is the inward-facing normal of a circle centred at y = with fixed radius R — that is, n ys := − y / | y | . We parameterise y b with the variable η and simi-larly decompose y b ( η ) into y b = ( η ) , y b || ( η ) and y bs ( η ) , which represent the points on the boundaries ∂ω = ( x ) , ∂ω || ( x ) and ∂ω s ( x ) , respectively. Thus, y b = ( η ) := (cid:26) y = η, y = ±
12 : η ∈ (cid:18) − a ( x )2 , a ( x )2 (cid:19)(cid:27) , (A8a) y b || ( η ) := (cid:26) y = ± a ( x )2 , y = η : η ∈ (cid:18) − , (cid:19)(cid:27) , (A8b)and y bs ( η ) := { y = R ( x ) cos η, y = R ( x ) sin η : η ∈ (0 , π ) } . (A8c)Since R is independent of x , the tensor G has only two non-zero components, G and G . On ∂ω = and ∂ω || , G has at most one non-zero component: ∂y b = , /∂x or ∂y b || , /∂x , respectively.On ∂ω = , n y = = (0 , ± so that n y = · G = and hence does not contribute to the transporttheorem (Eq. A6). For ∂ω || , n y || = ( − , and on ∂ω || , n y || = (1 , so that on both ∂ω || and ∂ω || , n y || · G = (cid:18) ∂a∂x , (cid:19) . (A9)30e parameterise n ys with θ ∈ (0 , π ) so that n ys = − (cos θ, sin θ ) , (A10a)which gives n ys · G = (cid:18) − ∂R∂x , (cid:19) . (A10b)Thus, the transport theorem (Eq. A6) becomes ∇ x · (cid:90) ω f z ( x , y ) d S y = (cid:90) ω f ∇ x · z ( x , y ) d S y + (cid:90) − ∂a∂x z d y − (cid:90) ∂ω s ∂R∂x z d s y , (A11)where the last integral is straightforward to parameterise in terms of η (Eq. A8c). Appendix B: Incompressibility on the macroscale
In this Appendix, we show that ˆ ∇ · ( φ V ) vanishes. On application of the transport theorem(Eq. A11) and the product rule, the divergence of the flux is ˆ ∇ · ( φ V (0) ) = 1 a (cid:90) ω f ∇ x · v (0) d S y + 1 a (cid:90) − ∂a∂x v (0)1 d y − a ( ∇ x a ) · (cid:90) ω f v (0) d S y − (cid:90) ∂ω s ∂R∂x v (0)1 d s y , (B1)where the last term will vanish due to the no-slip and no-penetration conditions on ∂ω s (Eq. 17c).Considering O y (1) terms in Equation (14b) and the divergence theorem, we rewrite the first termon the right-hand side of Equation (B1) as a (cid:90) ω f ∇ x · v (0) d S y = − a (cid:90) ∂ω v (1) · n y d s y = − a (cid:32)(cid:90) ∂ω = v (1) · n y = d s y + (cid:90) ∂ω || v (1) · n y || d s y + (cid:90) ∂ω s v (1) · n ys d s y (cid:33) . (B2)The first two terms on the right-hand side of Equation (B2) vanish due to periodicity of v (1) on ∂ω = and ∂ω || , respectively and the last term on the right-hand side vanishes due to the no-slip andno-penetration boundary condition (Eq. 17c). Thus a (cid:90) ω f ∇ x · v (0) d S y = 0 . (B3)31he second term on the right-hand side of Equation (B1) becomes a (cid:90) − ∂a∂x v (0)1 d y = 1 a ∂a∂x (cid:90) − v (0)1 d y = ¯ v (0)1 a ∂a∂x , (B4)since a is independent of y and where ¯ v (0)1 ( x ) := (cid:90) / − / v (0)1 ( x , y ) d y . (B5)We interpret ¯ v to be the total longitudinal flux through any transverse cross-section of the mi-croscale problem. Note that this quantity must be independent of y due to conservation of mass:recall that Equation (17b) gives ∇ y · v (0) = 0 , so that (cid:90) − ∇ y · v (0) d y = ∂ ¯ v (0)1 ∂y + (cid:90) − ∂v (0)2 ∂y d y = ∂ ¯ v (0)1 ∂y + (cid:104) v (0)2 (cid:105) − = ∂ ¯ v (0)1 ∂y . (B6)The third term on the right-hand side of Equation (B1) is − a ∇ x a · (cid:90) ω f ( x ) v (0) d S y = − a (cid:90) ω f ( x ) ∇ x a · v (0) d S y = − a ∂a∂x (cid:90) ω f ( x ) v (0)1 d S y = − a ∂a∂x (cid:90) a − a (cid:90) − v (0)1 d y d y = − a ∂a∂x (cid:90) a − a ¯ v (0)1 d y . (B7)Hence Equation (B7) becomes − a ∇ x a · (cid:90) ω f ( x ) v (0) d S y = − a ∂a∂x ¯ v (0)1 (cid:90) a − a d y = − ¯ v (0)1 a ∂a∂x . (B8)Thus, combining the results from Equations (B4) and (B8) we find that, at leading order, ˆ ∇ · (cid:0) φ V (0) (cid:1) = 0 , (B9)so that total mass of fluid entering the porous medium at the inlet is equal to the total mass of fluidat the outlet. Note that it is the divergence of the leading-order Darcy flux Q (0) := φ V (0) thatvanishes. Appendix C: Independence of ¯ B with respect to y In this Appendix we show that ¯ B is also independent of y . Equation (27a) gives ∇ y c (1) = ∇ y · (cid:0) ∇ y c (1) − v (0) c (0) (cid:1) = ∇ y · (cid:0) B − ∇ x c (0) (cid:1) = ∇ y · B = 0 , (C1)32ecalling that c (0) is independent of y . Hence (cid:90) − ∇ y · B d y = (cid:90) − ∂ B ∂y + ∂ B ∂y d y = ∂∂y (cid:32)(cid:90) − B d y (cid:33) + [ B ] − = ∂ ¯ B ∂y , (C2)where the last equality is due to the periodicity of c (1) , c (0) and v (0) and their gradients. Thus, ¯ B is independent of y . Appendix D: Non-monotonicity of D In this Appendix, we show the diffusivity tensor D (Figure 8 (a)–(d)) and examine the non-monotonicity of the longitudinal diffusivity D for fixed R as φ varies. We consider the minimumlongitudinal diffusivity D (cid:63) and the unique value of φ (cid:63) to which it corresponds (Figure 8 (e)). Eachvalue φ (cid:63) corresponds to a particular pair a (cid:63) and R (cid:63) (Figure 8 (f)). Note that a (cid:63) is approximatelyrelated to R (cid:63) via a (cid:63) = 2 R (cid:63) √ π (Figure 8 (f), blue dot-dashed line). This linear relationship is agood fit for small R (cid:63) , but slightly overestimates the true value of a (cid:63) for larger values of R (cid:63) .33 IG. 8. (a) and (c) D and D , respectively, against φ for fixed values of R ∈ { . , . , . , . , . } ,with a varying according to Equation (12). (b) and (d) D and D , respectively, against R for fixed valuesof φ ∈ { . , . , . , . , . } , with a varying according to Equation (12). For a given value of R , theminimum porosity φ min ( R ) is given by Equation (12) with a = 2 R . Note that D is non-zero at φ min ( R ) for all values of φ (dot-dashed curve, top row), whereas D vanishes at the corresponding φ min . Both D and D are as defined in Equation (40b) and calculated using COMSOL Multiphysics ® . Both D and D tend to 1 as φ → , which corresponds to the limit of free-space diffusion. The effective diffusivity isisotropic when a = 1 (solid black curves, top 2 rows), in agreement with the results presented in Dalwadi et al. [8]. (e) The minimum longitudinal diffusivity D = D (cid:63) for each value of R in Figure 5(a) againstthe corresponding porosity φ = φ (cid:63) and (f) the corresponding value of a = a (cid:63) against R (cid:63) . The latter iswell predicted by the line a (cid:63) = 2 R (cid:63) √ π (blue dot-dashed line) for small R (cid:63) , becoming non-monotonic near R (cid:63) = 0 . .
1] B
ECKWITH , C. W, B
AIRD , A. J. & H
EATHWAITE , A. L. 2003 Anisotropy and depth-related het-erogeneity of hydraulic conductivity in a bog peat. II: modelling the effects on groundwater flow.
Hydrological processes (1), 103–113.[2] B EN ´ ITEZ , J. J., T
OPOLANCIK , J., T
IAN , H. C., W
ALLIN , C. B, L
ATULIPPE , D. R., S
ZETO , K.,M
URPHY , P. J., C
IPRIANY , B. R., L
EVY , S. L., S
OLOWAY , P. D. & C
RAIGHEAD , H. G. 2012Microfluidic extraction, stretching and analysis of human chromosomal DNA from single cells.
Labon a Chip , 4848–4854.[3] B LUM , J. J., L
AWLER , G., R
EED , M. & S
HIN , I. 1989 Effect of cytoskeletal geometry on intracel-lular diffusion.
Biophysical Journal , 995–1005.[4] B RUNA , M. & C
HAPMAN , S. J. 2015 Diffusion in spatially varying porous media.
SIAM Journal onApplied Mathematics (4), 1648–1674.[5] B RUSSEAU , M. L. 1994 Transport of reactive contaminants in heterogeneous porous media.
Reviewsof Geophysics (3), 285–313.[6] C LAVAUD , J.-B., M
AINEULT , A., Z
AMORA , M., R
ASOLOFOSAON , P. & S
CHLITTER , C. 2008 Per-meability anisotropy and its relations with porous medium structure.
Journal of Geophysical Research:Solid Earth (B1).[7] D
ALWADI , M. P., B
RUNA , M. & G
RIFFITHS , I. M. 2016 A multiscale method to calculate filterblockage.
Journal of Fluid Mechanics , 264–289.[8] D
ALWADI , M. P., G
RIFFITHS , I. M. & B
RUNA , M. 2015 Understanding how porosity gradients canmake a better filter using homogenization theory.
Proceedings of the Royal Society A , 20150464.[9] D
ALY , K. R. & R
OOSE , T. 2015 Homogenization of two fluid flow in porous media.
Proceedings ofthe Royal Society A , 20140564.[10] D
AVIT , Y., B
ELL , C. G., B
YRNE , H. M., C
HAPMAN , L. A. C., K
IMPTON , L. S., L
ANG , G. E.,L
EONARD , K. H. L., O
LIVER , J. M., P
EARSON , N. C., S
HIPLEY , R. J., L., W
ATERS
S., P.,W
HITELEY
J., D., W
OOD
B. & M., Q
UINTARD
Advances in Water Resources ,178–206.[11] D AVIT , Y., B
YRNE , H., O
SBORNE , J., P
ITT -F RANCIS , J., G
AVAGHAN , D. & Q
UINTARD , M. 2013Hydrodynamic dispersion within porous biofilms.
Physical Review E (1), 012718.
12] D
OMENICO , P. A. & S
CHWARTZ , F. W. 1990
Physical and Chemical Hydrogeology . Wiley, NewYork.[13] F
RITTON , S. P. & W
EINBAUM , S. 2009 Fluid and solute transport in bone: flow-induced mechan-otransduction.
Annual Review of Fluid Mechanics , 347–374.[14] L I , N., W EI , W., X IE , K., T AN , J., Z HANG , L., L UO , X., Y UAN , K., S
ONG , Q., L I , H., S HEN ,C., R
YAN , E. M., L
ING , L. & B
INGQING , W. 2018 Suppressing dendritic lithium formation usingporous media in lithium metal-based batteries.
Nano Letters , 2067–2073.[15] M ARIANI , G., F
ABBRI , M., N
EGRINI , F. & R
IBANI , P. L. 2010 High-Gradient Magnetic Separationof pollutant from wastewaters using permanent magnets.
Separation and Purification Technology ,147–155.[16] O’D EA , R. D., N ELSON , M. R., E L H AJ , A. J., W ATERS , S. L. & B
YRNE , H. M. 2015 A multi-scale analysis of nutrient transport and biological tissue growth in vitro.
Mathematical Medicine andBiology: A Journal of the IMA (3), 345–366.[17] P RINTSYPAR , G., B
RUNA , M. & G
RIFFITHS , I. M. 2019 The influence of porous-medium mi-crostructure on filtration.
Journal of Fluid Mechanics , 484–516.[18] Q
UINTARD , M. & W
HITAKER , S. 1994 Convection, dispersion, and interfacial transport of contami-nants: Homogeneous porous media.
Advances in Water Resources , 221–239.[19] R OSTI , M. E., P
RAMANIK , S., B
RANDT , L. & M
ITRA , D. 2020 The breakdown of Darcy’s law in asoft porous material.
Soft Matter , 939–944.[20] S ALLES , J., T
HOVERT , J.-F., D
ELANNAY , R., P
REVORS , L., A
URIAULT , J.-L. & A
DLER , P. M.1993 Taylor dispersion in porous media. Determination of the dispersion tensor.
Physics of Fluids A:Fluid Dynamics (10), 2348–2376.[21] S HIPLEY , R. J. & C
HAPMAN , S. J. 2010 Multiscale modelling of fluid and drug transport in vasculartumours.
Bulletin of Mathematical Biology , 1464–1491.[22] S PYCHAŁA , M. & S
TARZYK , J. 2015 Bacteria in non-woven textile filters for domestic wastewatertreatment.
Environmental Technology (8), 937–945.[23] T OMIN , P. & L
UNATI , I. 2016 Investigating Darcy-scale assumptions by means of a multiphysicsalgorithm.
Advances in Water Resources , 80–91.[24] W ANG , M., L IU , H., Z AK , D. & L ENNARTZ , B. 2020 Effect of anisotropy on solute transport indegraded fen peat soils.
Hydrological Processes (9), 2128–2138.[25] W ANG , Z., W U , H.-J, F INE , D., S
CHMULEN , J., H U , Y., G ODIN , B., Z
HANG , J. X. J. & L IU , X.
013 Ciliated micropillars for the microfluidic-based isolation of nanoscale lipid vesicles.
Lab on aChip , 2879–2882.[26] W HITAKER , S. 1986 Flow in porous media I: A theoretical derivation of Darcy’s law.
Transport inPorous Media , 3–25.[27] W HITAKER , S. 2013
The Method of Volume Averaging . Springer Science & Business Media, B.V.[28] W
OOD , B. D., C
HERBLANC , F., Q
UINTARD , M. & W
HITAKER , S. 2003 Volume averaging fordetermining the effective dispersion tensor: closure using periodic unit cells and comparison withensemble averaging.
Water Resources Research (8), 1210.(8), 1210.