Fluid permeation through a membrane with infinitesimal permeability under Reynolds lubrication

To understand the lubrication-dominated permeation through a membrane, numerical simulations of permeation through a moving corrugated permeable membrane is carried out with a fully validated numerical method. Through comparisons between the numerical results and the results of an asymptotic analysis of permeate ﬂux (under an inﬁnitesimal permeability condition) using Reynolds lubrication equation, the effect of permeation on lubrication and its inverse effect (i


Introduction
Mass transfer through a selective permeable membrane is typically observed in biological environments, for example, oxygen transport by red blood cells (RBCs), gas exchange in the alveoli, and nutrient absorption.
The mechanical aspects of oxygen transport through blood vessels may be characterised, for example, by the vessel geometry and RBC-RBC interactions.For the vessel geometry, the flow and mass transports are strongly dependent on the shape of the blood vessel, such as branching, merging, bending [1,2], and stenosis [3,4].For RBC-RBC interactions, pressure developed by lubrication could become significant in a narrow gap region between RBCs (or between an RBC and the vessel wall), considering that the lubrication force on a rigid spherical particle.For example, the lubrication pressure develops in inverse proportion and logarithmic proportion to the surface-to-surface distance [5,6,7,8] under normal-and tangential-mode motions, respectively.Despite its importance, lubrication-associated mass transfer has not received any considerable attention [9], and its modelling has largely lagged behind.
One of the reasons is that the lubrication pressure is sometimes difficult to reproduce under a limited spatial resolution.In view of handling the relative motion of spherical membranes (equivalent diameter D p ), a fixed fluid mesh system that is non-conforming to the membrane surface is often employed, and the lubrication pressure in a narrow gap of arbitrary non-spherical geometry is described by the Reynolds lubrication equation [10].The Reynolds-type lubrication is applicable when the ratio of the gap width to D p , , is much smaller than unity and the Reynolds number (based on D p and a reference velocity) multiplied by 2 is also much smaller than unity.The advantage of the Reynolds lubrication equation is that the dimension of the problem is locally reduced by one, as the pressure in the Reynolds lubrication regime is independent of the wall-normal distance.However, considering feasible computational mesh size, the above geometric and dynamic assumptions for the Reynolds lubrication equation strongly restricts the applicability of the equation to model particle-induced flows [11] in the inter-particle region, where non-Reynolds lubrication effect becomes predominant.
Another reason may be the difficulty in performing numerical simulation of the trans-membrane flux.
The permeate flux is dependent on the jump values in the hydrostatic and osmotic pressures (and therefore, it is proportional to the concentration jump under a dilute condition of a liquid solution), and a sharp representation of the membrane is important to reproduce the permeate flux [12,13]; in particular, it is essential to capture the discontinuous jumps across the membrane.The sharpness is not always guaranteed in a fixed mesh system non-conforming to the membrane geometry, and the problem would be even more pronounced in a lubrication-dominant situation where the flow is intrinsically under-resolved.Early efforts on simulating membrane permeation problems [14,15,16,17] employed a diffused interface, where a substantial thickness of the interface may be inevitable [18].The present authors have focused on the sharpness problem and developed an effective approach (based on the discrete-forcing immersed boundary method [19,20]) that can represent the sharp discontinuity of pressure across a permeable membrane [21] by strictly satisfying mass conservation in the vicinity of the permeable membrane.In addition, the mass transfer problems form a complex numerical system as they constitute a coupled system between the permeate solute/solvent, the membrane motion, and the motion of the ambient fluid [13].
The present paper focuses on the lubrication-dominated permeation of a solvent (or a pure fluid) through a membrane, and aims to establish the first step to understand membrane permeation under lubrication.To set up a simplified problem, permeation is solved through a membrane with a corrugated geometry moving at a constant velocity along the axis of corrugation (i.e., the tangential motion) in a two-dimensional parallel channel.A model of local permeate flux through a membrane is derived in the limit of infinitesimal permeability by analytically solving the Reynolds lubrication equation.The model prediction is compared to the numerical result; to isolate the resolution problem in the lubrication-dominant region in the present study, a fully resolved numerical simulation is carried out.The dependence of permeate flux on the lubrication pressure as well as the effect of permeation on lubrication are discussed for several corrugation amplitudes and permeation coefficients.

Governing equations
The fluid is an incompressible Newtonian fluid with a constant density (ρ) and constant viscosity (μ).The governing equations of the fluid are the equation of continuity and Navier-Stokes (N-S) equation: where u is the velocity, t is the time, p is the pressure and ν = μ/ρ.
The volumetric flux of the pure fluid (i.e., fluid without any solute) across the membrane is modelled by the following equation [22]: where n is the unit normal vector pointing from the rear-side of the membrane Ω 1 to the front-side Ω 2 , L p is the permeable coefficient and [ [p] ] is the hydrostatic pressure jump (i.e., discontinuity) across the membrane calculated with the limiting pressure values at the interface on the respective sides of the membrane, p 1 and

Problem statement and assumptions
We study the permeation across a corrugated rigid membrane in a parallel channel moving at a constant velocity (U 0 ), as schematically shown in Fig. 1.The corrugated membrane has an infinitesimal thickness.
The geometry of the corrugation and its prescribed motion are given by the following function: where k (= 2π/ ) is the wavenumber and δ is a parameter between 0 and 1.The corrugation has an infinite extension in the x direction.We assume that the amplitude of corrugation δh 0 is sufficiently small with respect to the half-channel width h 0 (i.e., δ 1) and that h 0 is much smaller than the wavelength of the corrugation (i.e., h 0 2π/k).The latter is equivalent to the aspect ratio ε (= h 0 / ) being sufficiently small, i.e., ε 1.
In the present study, the permeation is solved by transforming the frame to that attached on the membrane for simplicity: the top-and bottom-flat walls travel at −U 0 along the centre axis of the stationary corrugated membrane.When ε and the Reynolds number Re = U 0 h 0 /ν respectively satisfy ε 1 and ε 2 Re 1, the lubrication can be described by the Reynolds lubrication equation in both regions between the corrugated membrane and the flat walls.Assuming a small amplitude of corrugation (i.e., δ 1), the normal vector on the membrane surface is nearly parallel with the y axis.Therefore, in the limit of zero amplitude, only the y component of the permeation flux J n is effective, and it is denoted as where e y is the unit vector in the y direction (Fig. 1).By solving the lubrication equation, J y is modelled as a function of the membrane geometry.4 Theoretical: asymptotic permeation flux J y (x) With a membrane of no permeability (L p = 0, and therefore, J y = 0), the pressures on the lower and upper sides (Ω 1 and Ω 2 , respectively) of the corrugation, p (0) 1 (x) and p (0) 2 (x), follow separate Reynolds lubrication equations, where the superscript (0) stands for the no-permeable case.On the moving frame of reference x * (= x − U 0 t), the membrane geometry and pressure are transformed as follows: and p (0) * 1 follows the Reynolds lubrication equation: The derivation of the equation is detailed in Refs.[23,24].Following the lubrication analysis in Ref. [20], the pressure on the original (i.e., fixed) frame of reference is given as follows: From the geometric symmetry of Ω 1 and Ω 2 (see Fig. 1) and the one-dimensionality of the gap pressure, the pressure in Ω 2 is given as . Under an infinitesimal L p , we can reasonably assume that the permeate flux is Hereafter, we refer to J y (under L p → 0) as the asymptotic permeate flux.Using Eq.( 5) and h(x+π/k, t) = 2h 0 − h(x, t), the asymptotic permeate flux is given as where H (= 2h 0 ) is the full channel width (see Fig. 1) and L is the non-dimensional permeability defined as L p μ/H.

Outline of the numerical method
The numerical method for permeation employed herein is the same as that used in Takeuchi et al. [21], and it is briefly explained in this section.
The fluid motion is solved in an Eulerian frame, and a Cartesian fixed mesh of a uniform mesh size (Δ) is employed.The governing equations are discretised using a finite difference method.The fluid variables are defined on the collocated arrangement, and the spatial discretisation is performed by a second order central finite difference.The membrane is represented by connected Lagrangian marker points.
The present authors have developed a method that directly discretises the governing equation even at the grid points adjacent to the interface, while at the same time, ensuring consistency between the incompressible velocity and pressure fields [19,21].By their "consistent direct discretisation" for the discrete-forcing immersed boundary (DF-IB) approach, the non-slip condition on the interface was strictly imposed in a discrete sense while satisfying the mass and momentum conservations.This method enables capturing the sharp distribution of the velocity and pressure at the interface.Recently, this method was extended to allow permeate flux through the membrane as explained as follows.
For discretising on the cells not adjacent to the interface, ordinary discretisation of the second-order accuracy is adopted.On the other hand, in the interface vicinity (i.e., on cells adjacent to the interface), all the terms are discretised by considering the distance from the cell centre to the membrane (Eqs.(A4)∼ (A7) in Appendix A).The convective and viscous terms are time-updated by 4th-order Runge-Kutta method.A fractional step method is employed for coupling the velocity and pressure fields, and the velocity is corrected The above consistent direct discretisation for the DF-IB method is a suitable approach for membrane permeation, as the pressure fields on both sides reflect the local incompressible conservation at a cell level, which results in a sharp representation of the interface.This method has been fully validated for cases with a permeable membrane against a solvent [21] as well as for the non-permeable case (with stationary solid case [19] and moving/deforming membrane case [20]).The details of the numerical method can be found in Appendix A.

Results and discussion
In the present study, the domain length is fixed at 5H = 10h 0 , so that the aspect ratio of each flow region is sufficiently small: ε = h 0 / = 0.1.The Reynolds number is set at U 0 h 0 /ν = 0.5.Numerical simulation is carried out at the time increment Δt/(H/U 0 ) = 5 × 10 −5 and the spatial resolution h 0 /Δ = 20, unless specified otherwise.For this graph, the simulation is conducted at a higher spatial resolution of H/Δ = 50.
that the numerical result converges to Eq. ( 7) for all the L cases, even though the prediction assumes an infinitesimal L.
In the figure, the variation levels for L = 10 −5 ∼ 10 −3 are arranged at intervals of one order of magnitude, reflecting J y ∝ L (Eq.( 7)) in the infinitesimal condition.However, an irregular difference is observed for the case of L = 10 −2 , suggesting a non-trivial effect of membrane permeability on the lubrication pressure.This effect can be estimated with a higher-order term of pressure caused by the permeate flux.The following equation is solved for the induced pressure p (Jy) 1 caused by J y : The above equation is integrated twice with respect to x * by taking into account that both the average pressure and its gradient over 0 ≤ x ≤ are zero.Then, the effect of permeation on the lubrication pressure is found to be as follows: where The Taylor expansion of p (Jy) * 1 around δ = 0 is given by p (Jy) * 1 The above equation shows that, under the conditions of small δ and infinitesimal L, the permeate-induced pressure p (Jy) 1 is proportional to δ and L in the lowest order, and that the second-higher term (∝ δ 2 ) takes a mode of sin[2kx * ].Meanwhile, the involvement of the aspect ratio (ε) remains the same for the two terms, suggesting that the aspect ratio is not sensitive to the permeate-induced pressure.
Considering that ε = δ = O(10 −1 ) in the present study, Eq. (11) shows that p As explained above, the lubrication may not be explained by the ordinary Reynolds lubrication equation.
In the following, the effect of non-Reynolds lubrication is briefly studied.
Eq.( 7) suggests that ) is independent of L, which facilitates the evaluation of the convergence trends of the numerically-obtained pressure jump in the limit of L → 0. It is rather surprising that Eq.( 7) approximates the asymptotic permeate flux distributions well even at δ = 0.50, which is beyond the ideal Reynolds lubrication region.In Figure 4, the numerical solutions exhibit convergence to the respective asymptotic flux distributions for all the δ cases.
However, for all the δ cases, the maxima of the numerical asymptotic fluxes disagree with those of the analytical prediction, Eq.( 7), which is represented by the solid line in the figure.It might be expected that further higher-order flux components induced by p (Jy) * 1 alleviate the above disagreement in J y .However, considering that the order of magnitude of the flux L p p (Jy) 1 is L 2 δε −3 from Eq.( 11), the effect of p (Jy) 1 on J y is too small to fill up the disagreement.
- In the figure, "Analytical" means the prediction by Eq.( 7) under the Reynolds lubrication and infinitesimal permeability.
The effect of the x component of the flux is also insignificant as estimated in the following.The flux developing in the normal direction (n) is taken into account in the following lubrication equation for Ω 1 region: where (n x , n y ) are the x and y components of n, respectively, J x = J n • e x , and Here, p represents the pressure adjustment due to J x .The lowest-order mode of p : The normal vector n is also linearised as (n x , n y ) (−dh/dx, 1).Note that the magnitude of this n x is in the order of εδ; considering that the above lubrication equation with (n x , n y ) = (0, 1) reduces to Eq.( 8), Eq.( 12) describes the effect of n x on the pressure field.Then, p = O (δ 2 ).
Our preliminary study suggests that the disagreements in Fig. 4 which is only slightly smaller than J * y .

Conclusion
Under an infinitesimal permeability, the effective range of Reynolds lubrication on the membrane permeation was studied.The permeation through a moving corrugated permeable membrane was simulated with a fully-validated numerical method, and the results were compared with an analytical prediction based on the Reynolds lubrication equation.
The effects of lubrication on the permeation, permeation on the lubrication pressure, and the direction of flux on the permeation were independently studied, and the dependence of permeate flux and lubrication pressure on the geometry were identified.Even for a membrane with a simple geometry and restricted permeability condition, the mass transfer through the membrane cannot be fully explained by the Reynolds lubrication theory.Through the disagreement between the spatial distributions of numerical asymptotic flux and theoretical prediction, the applicable range of the Reynolds lubrication theory on permeation was evaluated.
The results suggest that the interaction between the permeation and lubrication is essential to reproduce the lubrication-associated mass transfer.The fully resolved simulation result of the lubrication pressure was influenced by the relaxation of the pressure jump due to the permeation flux, which in particular indicates an interesting implication of the effect of permeation on lubrication.To clarify this non-trivial effect of permeation on lubrication and modelling the mass transfer with lubrication, further studies are necessary in the non-Reynolds lubrication regime.

A.2 Governing equations
The governing equations of a fluid are Eqs.( 1) and ( 2), and these are non-dimensionalised by the reference velocity U and reference length H as follows: Here, the pressure is non-denationalised by ρU 2 , and Re is the Reynolds number.Note that the same symbols are used as the original equations (Eqs.( 1) and ( 2)) for representing the non-dimensional variables, as there would be no possibility of misunderstandings.
The volumetric flux of the pure fluid (i.e., fluid without any solute) across the membrane is shown in a non-dimensional form as follows: where L is the non-dimensionalised permeable coefficient defined as L p μ/H.

A.3 Spatial discretisation
The governing equations are discretised by a finite difference method.The flow variables are defined on the collocated arrangement, and the spatial discretisation is defined by the second-order central finite difference.
The convective and viscous terms are time-updated by the 4th-order Runge-Kutta method.A fractional step method is employed for coupling the velocity and pressure fields.In the following, the discretisations in the boundary cells (based on the configuration in Fig. A1) are presented to show the special treatments considering the distance between the surface and the cell centre.
The discretisations for the convective and viscous terms are as follows: • Convective term at (i − 1, j) • Convective term at (i, j) • Viscous term at (i − 1, j)

2 Fig. 1 :
Fig.1: A corrugated permeable membrane travelling at constant speed U 0 in +x direction in a parallel channel.
for a given pressure field at the next time step (Eqs.(A12)∼ (A15)).By substituting the corrected velocity field into the discretised mass conservation equation and replacing the velocity components at the interface vicinity with the pressure jump [[p] ] across the membrane (using Eq.(3)), pseudo-Poisson equation for the permeated pressure field is obtained (Eqs.(A10) and (A11)).The pressure equation takes a consistent form with the discretised velocity field by considering the distance to the interface, thereby ensuring consistency between the pressure and incompressible velocity fields.

Figure 2
Figure2shows the instantaneous velocity fields observed from the frame on the membrane for two different non-dimensional permeabilities (L) and three different corrugation amplitude parameters (δ).For both L cases, the locally developed Couette flows for the smallest amplitude case (δ = 0.1) suggest that the Reynolds lubrication (with a line source of permeate flux on the membrane surface) may hold reasonably in the regions away from the membrane.This observation is discussed later.On the other hand, for larger amplitudes of the corrugation δ ≥ 0.25, two-dimensional characteristics of the flows are remarkable: the non-negligible y-dependence of the pressure as well as the comparable magnitude of permeate flux to U 0 are due to the enhancement of the pressure discontinuity[  [p]  ] in comparison to the case of δ = 0.1.The result suggests that, for larger δ cases, the assumptions for the Reynolds lubrication equation in the entire domain is already invalidated; that is, non-Reynolds lubrication takes place.

Figure 3 Fig. 2 :
Figure 3 plots the maximum variation in J y (with respect to the analytical prediction) against δ for four different permeabilities at one-order magnitude intervals.Straight lines are also plotted in the figure to show the first-order converging trend.The constant decreasing trends of the variation with decreasing δ suggest

Fig. 3 :
Fig.3: The maximum variation in J y (with respect to the analytical prediction) plotted against δ for four different permeabilities L = 10 −2 , 10 −3 , 10 −4 and 10 −5 .The lines indicate the first-order converging trend.For this graph, the simulation is conducted at a higher spatial resolution of H/Δ = 50.
Figure 4 compares the longitudinal distributions of J y (x)/(U 0 L) for different L for the following three different δ values: (a) 0.10, (b) 0.25 and (c) 0.50.
solving the above equation with the lowest-order components of p (0) 1 , J y , and p (Jy) 1 are mainly caused by the effect of non-Reynolds lubrication[25] due to non-negligible gap width, even for the δ = 0.1 case in which the velocity fields in Figs.2(a) and 2(b) (as well as the pressure fields in Figs.B2(a) and B2(b) in Appendix B) look reasonably close to those of the Couette profile.While the Reynolds lubrication analysis, Eq.(14b), indicates J * y ∼ ε −1 δ, the non-Reynolds lubrication induces a flux with a magnitude of

FigureFig. A1 :
FigureA1shows the schematic of a membrane in an incompressible fluid in two dimensions.Hereafter, the computational cells partitioned by the object boundary are referred to as "boundary cells", as represented by the triangular symbol in the figure.For the boundary cell (i, j) in Fig.A1, the discretisations incorporating the boundary conditions are explained in two dimensions.In the following, (•) represents an interpolation operator of second-order accuracy, and δ x k is an operator of second-order central finite difference in the x k direction.The velocities at the cell face are denoted by U k (k = 1, 2) or (U, V ), and the fractional-step velocities (by excluding the pressure gradient) at the cell centre and cell face are represented by u * * and U * * k , respectively.The subscripts "b ∓ 0" stand for the (ε − + 1) Δx u b−0 − u i−1,j ε − Δx − (δ x u) i− 3