CET Volume 86 DOI: 10.3303/CET2186193 Paper Received: 5 September 2020; Revised: 1 March 2021; Accepted: 10 May 2021 Please cite this article as: Procopio G., Adrover A., Giona M., 2021, Generalized Reflection Method for the Stokes Equation in Confined Geometries: Applications to Microfluidics and Particle Transport, Chemical Engineering Transactions, 86, 1153-1158 DOI:10.3303/CET2186193 CHEMICAL ENGINEERING TRANSACTIONS VOL. 86, 2021 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Sauro Pierucci, Jiří Jaromír Klemeš Copyright © 2021, AIDIC Servizi S.r.l. ISBN 978-88-95608-84-6; ISSN 2283-9216 Generalized Reflection Method for the Stokes Equation in Confined Geometries: Applications To Microfluidics and Particle Transport Giuseppe Procopio, Alessandra Adrover, Massimiliano Giona* DICMA, Università degli Studi di Roma ``La Sapienza'', via Eudossiana 18, 00184, Roma, Italy massimiliano.giona@uniroma1.it A generalized reflection method is proposed for determining the hydrodynamic resistance matrix of a rigid particle translating in a confined Stokesian fluid, by reducing its estimate to the solution of two simpler problems: (i) the singular expansion of the hydrodynamic variables due to the particle motion in the free space, represented by the Faxén operator associated with the particle geometry, and (ii) the determination of the Green function associated with the specific channel geometry. This approach extends the Stokeslet approximation to a Faxén approximation that is reliable to ratios of the particles characteristic length to the characteristic distance from the walls higher than the existing literature approximations, the error of which is order of O /d . Given the solution of the two above mentioned subproblems, the proposed method applies to any particle and channel geometry. The application to prolate spheroidal particles translating parallel to a plane is discussed and compared against numerical simulations. 1. Introduction The drag force acting on a spherical rigid particle moving with velocity in a Newtonian unconstrained fluid at low Reynolds numbers is provided by the celebrated Stokes law 6πμr , where is the fluid viscosity and r the particle radius. However, a great deal of applications involve particles possessing generic geometries and flow systems confined by rigid walls. Typical applications concern separation in microfluidics devices (Adrover et al., 2019), chromatography (Edam et al., 2011), hydrodynamic chromatography (Adrover et al., 2018), haemodynamics (Goldsmith and Skalak, 1975), the modelling complex fluids (Happel and Brenner 2012) and, in principle, any low-Reynolds number hydrodynamic problem in which the characteristic length of the device (L) is comparable with the characteristic linear particle size (Happel and Brenner, 2012). The main effect of a confined hydrodynamics as regards the interaction with a moving rigid particle is that the drag forces become position dependent and increase either as → L or as the distance of the particle to the walls decreases. This poses interesting transport issues as, owing to the fluctuation-dissipation theorem in isothermal conditions, the effective transport parameters (diffusivity tensor) do depend on particle position within the hydrodynamic device. Closed-form expressions for the drag forces have been determined only in a handful of axial-symmetric cases, such for a sphere near a plane (Brenner, 1961), or at the center of a cylindrical channel (Haberman and Sayre, 1958). However, also in these simple cases, the exact functional dependence on the particle position is not available. Explicit approximated relations can be derived only under certain assumptions, e.g., whenever the particle distance to the walls is much greater than the characteristic linear particle size. These approximations involve the determination of the lower-order terms in the expansion of the solution in series of /d, in which the error is order of O /d . A pioneering work on the drag forces in confined geometries has been carried out by Lorentz in the case of a sphere near a plane by using a method of images in order to determine the lower-order terms of the exact solution. Subsequently, Faxén, see (Happel and Brenner, 2012), analyzed the problem of a sphere between two parallel planes by considering two special geometrical configurations: where the sphere is placed at the midpoint and at the quarter distance from one of the plates. 1153 The problem of a sphere moving in a cylindrical channel has been addressed by several authors see e.g. (Hasimoto, 1976). As regards non-spherical particles, attention has been focused on slender bodies (Barta and Liron, 1988), whereas few contributions have considered ellipsoids (Wakiya, 1957). A reflection method of general applicability has been developed by Brenner (Brenner, 1962), with the aid of which a general expression for the drag forces, correct up to terms O /d , has been obtained (Cox and Brenner, 1967). Nevertheless, a general relation, explicit with respect to the geometric parameters of the particle and of channel and correct up to higher-order terms in l/d is still not available. In this work we propose a new approach, based on the reflection method, for determining the forces acting on a rigid particle within a confined geometry from the solution of two simpler problems: i) the general Faxén operator, associated with the motion of the body in the free space (unbounded case), and ii) the Green's function of the Stokes flow within the channel under consideration in the absence of the particle, in the presence of an impulsive perturbation. This approach applies under the hypothesis that d ≫ and proves to be correct up to terms order of O /d . The two above mentioned sub-problems, upon which the present approach is based, are significantly easier to tackle than the original problem of a rigid particle within a constrained geometry, and their solutions are, in many cases, available in the literature, at least in series expansion to the desired degree of accuracy. Due to space limitations of the present communication, the derivation of the method is only briefly outlined. A complete derivation can be found in (Procopio and Giona, 2021). As a case study we consider the estimate of the drag acting on a prolate spheroid translating parallel to a rigid plane in three different configurations, comparing the theoretical results to the numerical hydrodynamic simulations. 2. The reflection method To begin with, let us briefly review the basic steps in the reflection method (Brenner, 1962). This method is grounded on the linearity of the Stokes problem for an incompressible flow μ∇ − ∇p = 0 , ∇ ⋅ = 0 (1) where and p are the velocity pressure fields at within the flow domain. No-slip conditions apply both at the external particle surface S , moving at velocity and at fixed solid walls S of the domain = 0   for    ∈ S , = 0   for    ∈ S (2) Owing to the linearity of the Stokes problem, the solution for the velocity and pressure fields can be expressed as the superposition of a countable system of fields = + + ⋯ + ⋯ , p = p + p + ⋯ + p + ⋯ (3) such that: (i) the first field is due to the motion of the particle in an unbounded domain (free space) with = for ∈ S (ii) any even velocity field cancels the preceding one in the expansion at the walls S , and (iii) any odd field cancels the preceding one at the particle surface S = − for ∈ S , = − for ∈ S (4)i = 1,2... As a consequence, the total stress tensor admits the analogous series expansion = + + ⋯ + + ⋯ , = p − μ  ∇ ∇ (5) being the identity matrix, so that the drag force acting on the particle becomes = − ⋅ d = + + ⋯ + + ⋯ (6) Owing to that properties that the even fields are smooth on the surface of the particle and that the stress tensor is divergence-free (by definition of the Stokes flow), it is possible to show, using the Gauss divergence theorem, that their contribution to the total force is vanishing, = − ⋅ d = − ∇ ⋅ dV = 0 (7) 1154 V being the particle domain, and thus = + + ⋯ + + ⋯ (8) 3. Estimate of the forces This paragraph provides a relatively simple way to estimate the forces entering eq. (8) from the solution of two simpler problems. Below, the conceptual steps in the formulation of the method are outlined, leaving all the (lengthy) calculations to (Procopio and Giona, 2021). We make use of the typical notation used in Stokesian mechanics as developed e.g., in classical monographs in the field (Kim and Karrila, 2005). The drag force acting on the particle moving in an unbounded fluid depends linearly both on the viscosity and on its velocity . Therefore, the force related to the unbounded problem can be expressed as = −8πμ ⋅ (9) where represents the resistance matrix which depends solely on the particle geometry. On the other hand, the force , which is of order O /d , is that experienced by the particle at a point as it were immersed within the field . This is because at the surface of the particle + = . Therefore, we can use the generalized Faxén theorem (Kim, 1985) to expressed = −8πμ , ⋅ (10) where , is an integro-differential operator referred to as the generalized Faxén operator (Kim and Karrila, 2013), which, applied to a field returns the components of the force experienced by a body placed at . As shown by Kim (Kim, 1985), the Faxén operator corresponds to the operator that once applied to the poles of the Oseen's tensor , , returns the velocity field of the unbounded problem, i.e., = ⋅ , ⋅ , (11) as the Oseen's tensor , represents the Green's function of the unbounded Stokes problem ∇ , − ∇ , = −8πμ δ − , ∇ ⋅ , = (12) see (Happel and Brenner, 2012). As known (Pozrikidis, 1992), any solution of a Stokes problem in a domain can be expressed in terms of singularities placed outside the domain. Specifically, the Green's function with a given domain with boundaries S can be written as the sum of the free-space Oseen's tensor plus a contribution , due to these singularities , = , + , (13) where , is such that , = − , for ∈ S (14) and it is singular at poles outside flow domain and regular inside. It is possible to show (Procopio and Giona, 2021), that can be expressed as = ⋅ (15) where the operator is given by = , ⋅ , ⋅ , ⋅ (16) where ’ are the poles of , . Thus, the operator can be obtained from the knowledge of the unbounded Stokes problem for the rigid particle (which provide the expression for the Faxen's operator and for ) and of 1155 the Green's function , associated with the constrained geometry. Under the hypothesis that d ≫ l, the iterative application of the reflection method provides the following approximate relation (Procopio and Giona, 2021) ≈ ⋅ (17) In this way, the substitution into the series expansion for the drag force yields the following compact result = −8πμ ⋅ ⋅ (18) where = (19) 3.1. A case study: spheroidal particles near a plane The theory outlined in the previous paragraph, and specifically eqs. (18)-(19), can be tested in a non trivial case represented by the motion of a spheroidal particle close to a plane. Let h be the distance between the plane and the centroid of the spheroidal particle. We consider a prolate spheroid, geometrically described by equation + = 0 (20) where and are the half-lenghts of the major and minor axes respectively. Given the eccentricity of the spheroid, and are related by = √1 − while the focal length equals = . The Faxén operator for this class of spheroids can be found in the literature (Kim and Karrila, 2013) and takes the form , = 2 d 1 + 1 − 1 − ∇4 (21) where the integration variable is the dimensionless coordinate spanning the major axis centered at the position . In this case the resistance matrix can be expressed as = 2c , where is a symmetric matrix, depending exclusively on the eccentricity and possessing an eigenvalue A = −2 + 1 + log 1 + / 1 − (22) related to the motion along the major axis, and the double eigenvalue A = − 2−2 + 1 − 3 log 1 + / 1 − (23) associated with the motion along the minor axis. The regular part of the Green's function tensor , for the Stokes problem confined by a plane with normal directed towards the fluid is (Blake, 1971) , = − , + 2h ⋅ ⋅ ∇ , + ∇ , (24) where = − 2 ⋅ . Therefore, by applying equation (16) at = − ’ = h , the matrix entering eq. (18) can be evaluated for any orientation and translation of the spheroid. Three particular configurations are considered, where the spheroid can move parallel to the plane, thus with ⊥ : (a) the major axis of the spheroid is orientated along the direction of the motion ( ∥ ), (b) the major axis is orientated along the normal of the plane ( ∥ ), (c) the major axis is orientated perpendicular to both the direction of the motion and the normal ( ⊥ , ). 1156 Figure 1: Panel (a): Normalized velocity field coplanar to the normal to the plane and the velocity of the spheroid moving in the configuration (a). Panel (b): Dimensionless drag acting on the spheroid in the configuration (a) for several values of the eccentricity deriving from theoretical prediction (lines) and from FEM numerical simulation (symbols). For any of these configurations, it is possible to define the dimensionless drag λ = F / 8πμaU with j = a, b, c, where the drag forces F can be evaluated by applying the theory developed in the previous section. It is possible to provide closed-form expressions for λ as a function of the geometric parameters of each configuration, but the analytic expressions are not reported for the sake of length limitations of the present communication. In order to check the theoretical results, we have performed FEM numerical simulations by means of the software Comsol Multiphysics 5.5. Figure 1 panel (a) shows the velocity profile obtained numerically in the configuration (a). Figure 1 panel (b) shows the comparison between the values of λ obtained by the equation (18) and the numerical simulation for some values of the eccentricity. A satisfactory agreement between the theoretical approach proposed and hydrodynamic simulations can be observed. Analogous results have been obtained for λ and λ , but not reported for space limitations. With the aid of these results, it is possible evaluate which configuration (amongst cases a, b and c) experiences a minor drag, i.e., which of the three configuration is the most energetically advantageous. As shown in Figure 2 panel (a), where the dimensionless drags are plotted as function of / , the most advantageous configuration, at a given distance between the centroid and the plane, is always the configuration (a), which is also the most favourable configuration in unbounded conditions since the equivalent hydrodynamics radius is the half-length of the minor axis. Figure 2: Dimensionless drag acting on a spheroid as a function of the normalized distance of the centroid to the plane / (panel a) and of the gap / (panel b) in the three configurations (a) (solid lines), (b) (dashed lines) and (c) (dot-dashed lines) for two different values of the eccentricity. The other two configurations, that are equivalent in unbounded conditions, provide almost the same drag especially for small values of eccentricity . However, configuration (c) experiences slightly less resistance than configuration (b) and this becomes evident only for values of ≈ 1 and ≈ . On the other hand, if the dimensionless drag forces are compared respect to the same gap δ between the surface of the spheroid and the plane results provide a slightly different interpretation. In this case, the fluid 1157 volume confined by configurations (a) and (c) is always greater than in case (b). Therefore, as depicted in Figure 2 panel (b) the configuration experiencing greatest drag is case (b) for any value of / , whereas the most advantageous one is case (a) for high values of / and case (c) for spheroids very close to the plane. 4. Conclusions In this contribution, we have proposed a method to reduce the estimate of the resistance matrix for a rigid particle in a confined geometry to the solution of two simpler problems, related to the geometry of the particle (in unbounded conditions), and to the geometry of the hydrodynamic device in the presence of a pointwise forcing. As a case study the drag acting on a prolate spheroidal particle moving parallel to a plane in three different configurations has been considered, and the theoretical predictions agree perfectly with numerical simulations for any value of the eccentricity of the spheroid. Moreover, the theoretical predictions are accurate up to terms of O /d , while the existing literature provide extimate accurate up to O /d . A further extension of the method can be performed in order to include rigid rotations and torques acting on rigid particles. References Adrover A., Cerbelli S., & Giona M. (2018). Taming axial dispersion in hydrodynamic chromatography columns through wall patterning. Physics of Fluids, 30(4), 042002. Adrover A., Venditti C., & Giona M., 2019, Laminar dispersion at low and high Peclet numbers in a sinusoidal microtube: Point-size versus finite-size particles. Physics of Fluids, 31(6), 062003. Barta E. & Liron N., 1988, Slender body interactions for low Reynolds numbers—part I: body-wall interactions, SIAM Journal on Applied Mathematics, 48(5), 992-1008. Blake J. R., 1971, A note on the image system for a stokeslet in a no-slip boundary, Mathematical Proceedings of the Cambridge Philosophical Society 70, 303–310. Brenner H., 1961, The slow motion of a sphere through a viscous fluid towards a plane surface, Chemical Engineering Science 16, 242–251. Brenner H., 1962, Effect of finite boundaries on the Stokes resistance of an arbitrary particle, Journal of Fluid Mechanics 12, 35–48. Edam R., Eeltink S., Vanhoutte D. J., Kok W. T. & Schoenmakers P. J., 2011, Hydrodydinamic chromatography of macromolecules using polymer monolithic columns, Journal of Chromatography A 1218, 8638–8645. Goldsmith H. & Skalak R., 1975, Hemodynamics, Annual Review of Fluid Mechanics 7, 213–247. Haberman W. L. & Sayre R. M., 1958, Motion of rigid and fluid spheres in stationary and moving liquids inside cylindrical tubes tech. Rep., David Taylor Model Basin, Washington DC. Happel J. & Brenner H., 2012, Low Reynolds number hydrodynamics: with special applications to particulate media, Springer Science & Business Media, New York. Hasimoto H., 1976, Slow motion of a small sphere in a cylindrical domain. Journal of the Physical Society of Japan, 41(6), 2143-2144. Kim S., 1985, Note on faxen laws for nonspherical particles. International Journal of Multiphase Flow 11, 713– 719. Kim S. & Karrila S. J., 2013, Microhydrodynamics: principles and selected applications, Dover Publications, New York. Pozrikidis C., 1992, Boundary integral and singularity methods for linearized viscous flow, Cambridge University Press, Cambridge. Procopio G. and Giona M., 2021 Resistance on a generic particle moving in a confined Stokesian fluid (in preparation). Wakiya S., 1957, Viscous flows past a spheroid. Journal of the Physical Society of Japan 12, 1130–1141. 1158