Transport of solute and solvent driven by lubrication pressure through non-deformable permeable membranes

Abstract

A discrete-forcing immersed boundary method with permeable membranes is developed to investigate the effect of lubrication on the permeations of solute and solvent through membrane. The permeation models are incorporated into the discretisation at the fluid cells including the membrane, and discretised equations for the pressure Poisson equation and convection–diffusion equation for the solute are represented with the discontinuities at the membrane. The validity of the proposed method is established by the convergence of the numerical results of the permeate fluxes (solute and solvent) to higher-order analytical models in a lubrication-dominated flow field. As a model of the mass exchange between inside and outside of a biological cell flowing in a capillary, a circular membrane is placed between parallel flat plates, and the effect of lubrication is investigated by varying the distance between the membrane and the walls. The pressure discontinuity near the wall is larger than that at the stagnation point, which is a highlighted effect of lubrication. In the case of a small gap, the solute transport is dominated by convection inside the circular membrane and by diffusion outside. Through the time variation of the concentration in the circular membrane, lubrication is shown to enhance mass transport from/to inside and outside the membrane.

1 Introduction

In biological environments and industrial applications, mass transport through permeable membranes takes place in various ways. Exchange of solute and water through microvascular wall is largely passive (Michel and Curry 1999), and the relationship between the structural elements of the capillary wall and the permeability coefficients for solutes of various sizes has been determined by systematic studies (Sugihara-Seki and Fu 2005). As an example, the filtration performance of the kidney is significantly affected by the narrowing and occlusion of the vascular lumen (Cannon et al. 1974), which is normally held in tension by intravascular pressure. On the contrary, it has been indicated that an increased permeability of glomerular capillary may lead to capillary occlusion owing to protein deposition (Purkerson et al. 1976), and the lubrication at narrow capillaries may play an important role in those processes. Using a lubrication-based model, Secomb et al. (1998, 2001) demonstrated that the hydrostatic pressure generated within the endothelial surface layer alters the shape of the red blood cells and the wall-cell distance depending on the flow velocity in the capillary [as well as the geometry of the vessel wall (Secomb and Hsu 1996, 1997)]. By solving a coupled problem of hydrostatic and osmotic pressures across the endothelial surface layer, Hu and Weinbaum (1999) reported a non-uniform distribution of mass concentrations and the corresponding non-uniform distribution of effective osmotic pressure.

The above processes are commonly characterised by transport phenomena under lubrication. In the ideal lubrication state in a negligibly small gap region (between the interfaces), the pressure increases locally. However, there may be a number of cases in biological environments where the ideal condition for the lubrication theory is violated (Takeuchi et al. 2021). For example, in relatively large gaps, the theory could deviate from the conventional Reynolds lubrication theory (Takeuchi and Gu 2019), resulting in an underestimation of solute and solvent permeations driven by the pressure difference. One of the difficulties associated with the numerical analysis of lubrication is that by introducing $\epsilon$ as the ratio of the gap width and reference length, the minimum number of grid points that are required to capture the pressure increase generated by lubrication is $\epsilon ^{-1/2}$ (Takeuchi et al. 2021), and, the resolution, therefore, becomes insufficient on a coarse grid system; meanwhile, the fine grid system becomes computationally demanding. In addition, the lubrication pressure decreases owing to permeation, which makes the analysis or prediction difficult (Takeuchi et al. 2021).

Another difficulty in the numerical simulation of permeation for a two-component fluid (i.e. solute and solvent) is the accuracy of the flow around the membrane. For this problem, a number of numerical methods have been proposed. The immersed boundary (IB) method proposed by Peskin was first applied to the analysis of blood flow in a heart (Peskin 1972). The interaction force from the object is incorporated into the external force term in the equation of motion of the fluid, and an approximate delta function is used to distribute the interaction force from the Lagrange marker to the grid points of the surrounding fluid. Numerical methods in the framework of taking Eulerian and Lagrangian approaches for fluid and solid (such as the IB method), respectively, have been developed for heat and mass transport analysis. Gong et al. (2014) and Wang et al. (2020) proposed a numerical method for the analysis of oxygen transport with red blood cells by treating the diffusion flux as an independent variable and using an approximate delta function to incorporate the effect of the diffusion flux in the membrane into the source term of the advection–diffusion equation. On the other hand, the ghost-cell method, the immersed interface method and the direct forcing immersed boundary method are typical methods for imposing interface conditions on discretised equations. The ghost-cell method implicitly assigns interface conditions to a virtual cell in the solid phase adjacent to the fluid cell (called the ghost point) by an interpolation function using the values of the mirror point of the interface and the nearby cells in the fluid phase. The ghost-cell method has been applied, for example, to particle multiphase flows with reactions at the particle surface (Lu et al. 2018). The immersed interface method (LeVeque and Li 1994; Layton 2006; Jayathilake et al. 2010) reproduces a sharp interface by applying a finite difference discretisation that takes into account jumps at the interface without using interpolation functions. The pressure jump is calculated by the singular force acting on the interface and the concentration jump by interpolation from the fluid cell near the membrane (Jayathilake et al. 2010). Miyauchi et al. (2015, 2017) proposed a finite element formulation for fluid permeation through a deformable membrane in a two-component fluid by incorporating the discontinuities of pressure and concentration into the discretised equations, and they showed that the reproduction of the sharpness of the discontinuities in pressure and concentration fields at the membrane is important for the accurate prediction of permeate fluxes.

In the present study, to investigate mass transport induced by lubrication pressure (hereafter, lubrication-induced mass transport) in a two-component fluid separated by a membrane, we develop a discrete-forcing (DF) IB method to capture the concentration and pressure distributions sharply along the membrane and to analyse the membrane permeation accurately. In the proposed method, the membranes are represented by Lagrange markers to maintain a high resolution in the vicinity of the membranes, and a Cartesian mesh fixed in space is used to solve fluid flow along the arbitrary geometry of the membrane. The boundary condition on the fluid-membrane interface is enforced by directly specifying the boundary values into the discretised equations, thereby enabling the sharp representation of the object interface as well as tight conservations of mass and momentum (and, therefore, the sharp distributions of the pressure and mass concentration), which distinguishes the method from previous DF-IB methods.

In the present study, based on the concept of the above numerical method with membrane permeations of solute and solvent, we analyse the effect of lubrication on the mass transport through a membrane in two-dimensional space. For this purpose, a new pressure Poisson equation is derived by incorporating the concentration jump on the membrane. The numerical method is validated through a comparison with analytical predictions of permeate fluxes for the case of a moving corrugated membrane. To study the mass transport in a lubrication-dominant environment, a system with a circular membrane placed near a solid wall is set up for different wall-membrane distances, and the time development of concentration field is discussed by decomposing the fluxes into the components of convection, diffusion, and permeation.

2 Governing equations

The fluid is an incompressible Newtonian fluid governed by the continuity equation and the Navier–Stokes (N–S) equation:

$$\begin{aligned}&\nabla \cdot \varvec{u} = 0 \ , \end{aligned}$$

(1)

$$\begin{aligned}&\rho _{\mathrm{f}}\left( \frac{\partial ^{} \varvec{u}}{\partial t^{}} + \varvec{u} \cdot \nabla \varvec{u} \right) = -\nabla p + \mu _{\mathrm{f}}\nabla ^2 \varvec{u} \ , \end{aligned}$$

(2)

where $\varvec{u}$ is the fluid velocity, $\rho _{\mathrm{f}}$ is the fluid density, t is the time, p is the pressure, and $\mu _{\mathrm{f}}$ is the viscous coefficient.

The transport of solute is governed by the following unsteady convection-diffusion equation:

$$\begin{aligned} \frac{\partial ^{} c}{\partial t^{}} + \nabla \cdot \left( c \varvec{u} \right) = D \nabla ^2 c \ , \end{aligned}$$

(3)

where c and D are the solute concentration and diffusion coefficient, respectively.

The permeate fluxes of solvent and solute through membrane are modelled as follows (Katchalsky and Curran 1961):

$$\begin{aligned} \varvec{J}_v= L_{\mathrm{p}}\big (\llbracket \,p\,\rrbracket - \sigma \llbracket \,\Pi \,\rrbracket \big ) \varvec{n} \ , \end{aligned}$$

(4a)

$$\begin{aligned} \varvec{J}_s= L_{\mathrm{p}}{\overline{c}} \big \{ \left( 1 - \sigma \right) \llbracket \,p\,\rrbracket + \left( \mu _{\mathrm{L}}- \sigma \right) \llbracket \,\Pi \,\rrbracket \big \} \varvec{n} \ , \end{aligned}$$

(4b)

where $L_{\mathrm{p}}$ is the permeability, $\sigma$ is the repulsion coefficient, and $\mu _{\mathrm{L}}= L_{\mathrm{d}}/L_{\mathrm{p}}$, with $L_{\mathrm{d}}$ being the phenomenological coefficients relating the diffusion flow to the osmotic pressure $\Pi$. Denoting the limiting values of $\varphi$ on the membrane as $\varphi ^-$ and $\varphi ^+$ in the rear and front sides, respectively, the jump and average values of $\varphi$ are defined as $\llbracket \,\varphi \,\rrbracket = \varphi ^- -\varphi ^+$ and ${\overline{\varphi }}=\left( \varphi ^- +\varphi ^+ \right) /2$, respectively. The unit normal vector $\varvec{n}$ on the membrane is defined in the direction from the rear to front sides.

In this study, we assume that the two-component fluid (i.e. solute and solvent) is a dilute solution, and van’t Hoff’s equation $\Pi = R T c$ (R: the gas constant, T: the temperature) is used to convert $\llbracket \,\Pi \,\rrbracket$ in Eq. (4) into the concentration jump $\llbracket \,c\,\rrbracket$. Throughout the study, no deformation is considered for membranes.

3 Numerical method

In this section, discretisations at the boundary cells are explained using a DF-IB method. Boundary cells are fluid cells that are separated into two regions by the membrane, and are represented by the triangular symbol for the configuration in Fig. 1.

For a single-component fluid (i.e. no solute), the DF-IB method proposed by Sato et al. (2016) and Takeuchi et al. (2018) directly discretises the N-S equation even at the boundary cells, while at the same time, their method guarantees the consistency between the discretised equations for the incompressible velocity and pressure fields. Using their “consistent direct discretisation” for the DF-IB approach, the no-slip and impermeable conditions on the interface were strictly imposed in a discrete sense, while satisfying the mass and momentum conservations, which enables capturing the sharp distribution of the velocity and pressure at the interface. The DF-IB method was extended to enable permeation of the solvent through the membrane by Takeuchi et al. (2019) and Tazaki et al. (2020); they showed discretisations that consider the permeable condition for the solvent at the interface. In the present study, we further deal with membranes which also allow solute permeation (i.e. bi-permeation for a two-component fluid), and the discretisation of the governing equations considering solvent and solute permeation is explained in two-dimensional space based on the DF-IB method for the configuration shown in Fig. 1.

In the following, the Cartesian coordinate system is adopted, and uniform fluid cells are arranged in the domain. The variables are defined on the collocated grid points; the primary variables (i.e. velocities u and v, pressure p and concentration c) are defined at the centres of the fluid cells, and the gradients and contra-variant velocity components are on the cell faces. For each time step, the concentration field is updated, followed by the time-marching procedure of the flow field.

Fig. 1

Schematic of the membrane and the Cartesian grid. The membrane separates the fluid cells $(i{-}1,j)$ and (i, j), and these cells are referred to as “boundary cells”. The distance between the membrane and the centre of a boundary cell is denoted by $\omega ^{\mp }\,\Delta x$

3.1 Discretisation of fluid equations with solute permeation

Based on the approach in Takeuchi et al. (2018, 2019), the discretisation for the N–S equation is briefly explained, with a focus on the improvements to consider the concentration jump at the boundary cells $(i{-}1,j)$ and (i, j) in Fig. 1.

The incompressible velocity and pressure fields are coupled by a fractional-step method. The discretisations of the convective and viscous terms are the same as those provided in Takeuchi et al. (2018). The major differences from the impermeable case (Takeuchi et al. 2018) and permeable case for the solvent (Takeuchi et al. 2019) appear to incorporate the concentration jump at the membrane into discretisation. In the following, the explanation of the numerical method focuses on the pressure Poisson equation and velocity correction procedure.

In the numerical study, the governing equations are non-dimensionalised using the reference velocity U, the reference length H, the reference pressure $\rho _{\mathrm{f}}U^2$, the reference time H/U, and the reference concentration C, and the non-dimensional variables are denoted with a tilde $\displaystyle { \widetilde{(\cdot )} }$. For example, the non-dimensional forms of the permeate fluxes are as follows:

$$\begin{aligned} \frac{\varvec{J}_v}{U}= \mathrm {Re}L\big \{ \llbracket \,{\widetilde{p}}\,\rrbracket - \sigma \mu _{\mathrm{p}}\llbracket \,{\widetilde{c}}\,\rrbracket \big \} \varvec{n} \ , \end{aligned}$$

(5a)

$$\begin{aligned} \frac{\varvec{J}_s}{CU}= \mathrm {Re}L{\overline{{\widetilde{c}}}} \big \{ (1-\sigma ) \llbracket \,{\widetilde{p}}\,\rrbracket + \mu _{\mathrm{p}}(\mu _{\mathrm{L}}- \sigma ) \llbracket \,{\widetilde{c}}\,\rrbracket \big \}\varvec{n} \ , \end{aligned}$$

(5b)

where $\mathrm {Re}$ is the Reynolds number defined as $\rho _{\mathrm{f}}{U H}/\mu _{\mathrm{f}}$, $L={\mu _{\mathrm{f}}L_{\mathrm{p}}}/{H}$ is the non-dimensional permeability of the solvent, $\mu _{\mathrm{p}}={RTC}/{\rho _{\mathrm{f}}U^2}$ is the pressure ratio, and ${\overline{{\widetilde{c}}}}=\left( {\widetilde{c}}^{\,-} + {\widetilde{c}}^{\,+} \right) /2$.

At the boundary cells $(i-1,j)$ and (i, j), the pressure equations considering the discontinuities are expressed as follows:

• at the boundary cell $(i-1,j)$

  $$\begin{aligned}- \frac{\Delta {\widetilde{t}}}{\left( \omega ^- + 0.5 \right) \Delta {\widetilde{x}}} \left\{ \left( \delta _{{\widetilde{x}}} {\widetilde{p}}^{\,n+1} \right) _{i-\frac{3}{2},j} +\frac{\mathrm {Re}Ln_x\llbracket \,{\widetilde{p}}\,\rrbracket ^{n+1}}{ \Delta {\widetilde{t}}} \right\} + \Delta {\widetilde{t}}\left( \delta _{{\widetilde{y}}}\delta _{{\widetilde{y}}} {\widetilde{p}}^{\,n+1} \right) _{i-1,j} \nonumber \\\quad = \frac{\left( {\widetilde{u}}_m- \mathrm {Re}Ln_x \mu _{\mathrm{p}}\llbracket \,{\widetilde{c}}\,\rrbracket ^{n+1} \right) - {\widetilde{U}}^{\star \star }_{i-\frac{3}{2},j}}{\left( \omega ^- + 0.5 \right) \Delta {\widetilde{x}}} + \left( \delta _{{\widetilde{y}}} {\widetilde{V}}^{\star \star } \right) _{i-1,j} \ , \end{aligned}$$

  (6)

• at the boundary cell (i, j)

  $$\begin{aligned}\frac{\Delta {\widetilde{t}}}{\left( \omega ^+ + 0.5 \right) \Delta {\widetilde{x}}} \left\{ \left( \delta _{{\widetilde{x}}} {\widetilde{p}}^{\,n+1} \right) _{i+\frac{1}{2},j} +\frac{\mathrm {Re}Ln_x \llbracket \,{\widetilde{p}}\,\rrbracket ^{n+1}}{ \Delta {\widetilde{t}}} \right\} + \Delta {\widetilde{t}}\left( \delta _{{\widetilde{y}}}\delta _{{\widetilde{y}}} {\widetilde{p}}^{\,n+1} \right) _{i,j} \nonumber \\\quad = \frac{ {\widetilde{U}}^{\star \star }_{i+\frac{1}{2},j} - \left( {\widetilde{u}}_m- \mathrm {Re}Ln_x \mu _{\mathrm{p}}\llbracket \,{\widetilde{c}}\,\rrbracket ^{n+1} \right) }{\left( \omega ^+ + 0.5 \right) \Delta {\widetilde{x}}} + \left( \delta _{{\widetilde{y}}} {\widetilde{V}}^{\star \star } \right) _{i,j} \ , \end{aligned}$$

  (7)

where $\Delta {\widetilde{x}}$ is the grid spacing and $\Delta {\widetilde{t}}$ is the time increment, $\delta _{{\widetilde{x}}}$ and $\delta _{{\widetilde{y}}}$ are the second-order central differences, $n_x$ is the x component of the normal vector, $\omega ^{\mp }\Delta {{\widetilde{x}}}$ are the distances between the membrane and the cell centre (see Fig. 1), ${\widetilde{U}}$ and ${\widetilde{V}}$ are the velocities in the x and y directions at the cell face, respectively, and ${\widetilde{u}}_m$ is the membrane velocity in the x direction. The superscripts “${\star \star }$” and $(n+1)$ represent the fractional-step velocity and the time level, respectively. By linearly extrapolating the pressure values from both sides of the membrane to determine $\llbracket \,{\widetilde{p}}\,\rrbracket$, Eqs. (6) and (7) (together with the discretised equations at non-boundary cells) constitute a closed set of simultaneous equations.

Note that the discretisations of the convective and viscous terms may include $1/\omega ^{\mp }$, and the present study employs a wider stencil when $\omega ^{\mp }\rightarrow 0$ to cope with the singular behaviours. For more details, refer to Takeuchi et al. (2018).

The velocities at the fluid cell centre and cell face ${\widetilde{u}}_{i-1,j},{\widetilde{u}}_{i,j},{\widetilde{U}}_{i-\frac{1}{2},j}$ are corrected as follows:

• at the boundary cell $(i-1,j)$

  $$\begin{aligned}&{\widetilde{u}}^{\,n+1}_{i-1,j} = {\widetilde{u}}^{\star \star }_{i-1,j} - \Delta {\widetilde{t}}\left[ \phi _3^- \left( \frac{\delta {\widetilde{p}}^{\,n+1}}{\delta {\widetilde{x}}} \right) \right] _{i-1,j} , \end{aligned}$$

  (8a)
  $$\begin{aligned}&{\widetilde{U}}^{n+1}_{i-\frac{1}{2},j} = \left[ \phi _4^- \left( {\widetilde{U}}^{\star \star } \right) \right] _{i-\frac{1}{2},j} - \Delta {\widetilde{t}}\left[ \phi _5^- \left( \frac{\delta {\widetilde{p}}^{\,n+1}}{\delta {\widetilde{x}}} \right) \right] _{i-\frac{1}{2},j} \ , \end{aligned}$$

  (8b)

• at the boundary cell (i, j)

  $$\begin{aligned}&{\widetilde{u}}^{\,n+1}_{i,j} = {\widetilde{u}}^{\star \star }_{i,j} - \Delta {\widetilde{t}}\left[ \phi _3^+ \left( \frac{\delta {\widetilde{p}}^{\,n+1}}{\delta {\widetilde{x}}} \right) \right] _{i,j} , \end{aligned}$$

  (9a)
  $$\begin{aligned}&{\widetilde{U}}^{n+1}_{i-\frac{1}{2},j} = \left[ \phi _4^+ \left( {\widetilde{U}}^{\star \star } \right) \right] _{i-\frac{1}{2},j} - \Delta {\widetilde{t}}\left[ \phi _5^+ \left( \frac{\delta {\widetilde{p}}^{\,n+1}}{\delta {\widetilde{x}}} \right) \right] _{i-\frac{1}{2},j} \ , \end{aligned}$$

  (9b)

where the interpolation functions $\phi ^{\mp }_{3,4,5}$ are expressed with the discontinuities in the following form:

$$\begin{aligned}\left[ \phi _3^- \left( \frac{\delta {\widetilde{p}}^{\,n+1}}{\delta x} \right) \right] _{i-1,j} = \frac{1}{\omega ^- + 0.5} \left\{ \omega ^- \left( \delta _{{\widetilde{x}}} {\widetilde{p}}^{\,n+1} \right) _{i-\frac{3}{2},j} + 0.5 \left( \frac{\delta {\widetilde{p}}^{\,n+1}}{\delta x} \right) ^- \right\} ,\\\left[ \phi _3^+ \left( \frac{\delta {\widetilde{p}}^{\,n+1}}{\delta x} \right) \right] _{i,j} = \frac{1}{\omega ^+ + 0.5} \left\{ \omega ^+ \left( \delta _{{\widetilde{x}}} {\widetilde{p}}^{\,n+1} \right) _{i+\frac{1}{2},j} + 0.5 \left( \frac{\delta {\widetilde{p}}^{\,n+1}}{\delta x} \right) ^+ \right\} ,\\\left[ \phi _4^- \left( {\widetilde{U}}^{\star \star } \right) \right] _{i-\frac{1}{2},j} = \frac{1}{\omega ^- + 0.5} \left\{ \left( \omega ^+ - 0.5 \right) {\widetilde{U}}^{\star \star }_{i-\frac{3}{2},j} + u_{m} - \mathrm {Re}Ln_x \mu _{\mathrm{p}}\llbracket \,{\widetilde{c}}\,\rrbracket ^{n+1} \right\} ,\\\left[ \phi _4^+ \left( {\widetilde{U}}^{\star \star } \right) \right] _{i-\frac{1}{2},j} = \frac{1}{\omega ^+ + 0.5} \left\{ \left( \omega ^+ - 0.5 \right) {\widetilde{U}}^{\star \star }_{i+\frac{1}{2},j} + u_{m} - \mathrm {Re}Ln_x \mu _{\mathrm{p}}\llbracket \,{\widetilde{c}}\,\rrbracket ^{n+1} \right\} ,\\\left[ \phi _5^- \left( \frac{\delta {\widetilde{p}}^{\,n+1}}{\delta x} \right) \right] _{i-\frac{1}{2},j} = \frac{1}{\omega ^- + 0.5} \left\{ \left( \omega ^+ - 0.5 \right) \left( \delta _{{\widetilde{x}}} {\widetilde{p}}\right) _{i-\frac{3}{2},j} - \frac{\mathrm {Re}Ln_x}{\Delta {\widetilde{t}}} \llbracket \,{\widetilde{p}}\,\rrbracket ^{\,n+1} \right\} ,\\\left[ \phi _5^+ \left( \frac{\delta {\widetilde{p}}^{\,n+1}}{\delta x} \right) \right] _{i-\frac{1}{2},j} = \frac{1}{\omega ^+ + 0.5} \left\{ \left( \omega ^+ - 0.5 \right) \left( \delta _{{\widetilde{x}}} {\widetilde{p}}\right) _{i+\frac{1}{2},j} - \frac{\mathrm {Re}Ln_x}{\Delta {\widetilde{t}}} \llbracket \,{\widetilde{p}}\,\rrbracket ^{\,n+1} \right\} \ . \end{aligned}$$

3.2 Discretisation of solute equation with permeation

The x and y components of the solute flux $(c\varvec{u} - D \nabla c)$ are denoted as ${\widetilde{j}}_{sx}$ and ${\widetilde{j}}_{sy}$, respectively. Using ${\widetilde{u}}_m$ and Eq. (4), ${\widetilde{j}}_{sx}$ on the membrane is given as follows:

$$\begin{aligned} {\widetilde{j}}_{sx}^{\mp }= {\widetilde{c}}^{\,\mp } {\widetilde{u}}^\mp - \frac{1}{\mathrm {Pe}} \left( \delta _{{\widetilde{x}}}{{\widetilde{c}}} \right) ^\mp \ , \nonumber \\= {\widetilde{c}}^{\,\mp } {\widetilde{u}}_m+ \mathrm {Re}L{\overline{{\widetilde{c}}}} \left\{ (1-\sigma ) \llbracket \,{\widetilde{p}}\,\rrbracket + \mu _{\mathrm{p}}(\mu _{\mathrm{L}}- \sigma )\llbracket \,{\widetilde{c}}\,\rrbracket \right\} n_x \ , \end{aligned}$$

(10)

where $\mathrm {Pe}$ is the Peclet number. Then, using the Crank–Nicolson method for the time update of the diffusion term, the discretisation of the unsteady convection–diffusion equation (3) is given as follows:

• at the boundary cell $(i-1,j)$

  $$\begin{aligned}{\widetilde{c}}^{\,n+1}_{i{-}1,j} + \frac{\Delta {\widetilde{t}}}{2} \frac{ \left( ({\widetilde{c}}^{\,-})_{i,j}^{n+1} {\widetilde{u}}_m^{\,n} + \mathrm {Re}L{\overline{{\widetilde{c}}}}^{\,n} n_x \mu _{\mathrm{p}}(\mu _{\mathrm{L}}- \sigma ) \llbracket \,{\widetilde{c}}\,\rrbracket ^{n+1} \right) - \left( {\widetilde{j}}_{sx} \right) ^{n+1}_{i{-}\frac{3}{2},j} }{\left( \omega ^{-} + 0.5 \right) \Delta {\widetilde{x}}} + \frac{\Delta {\widetilde{t}}}{2} \left( \delta _{{\widetilde{y}}}\,{\widetilde{j}}_{sy} \right) ^{n+1}_{i{-}1,j} \nonumber \\\quad = {\widetilde{c}}^{\,n}_{i{-}1,j} - \frac{\Delta {\widetilde{t}}}{2} \frac{ \left( {\widetilde{j}}_{sx}^- \right) ^{n} - \left( {\widetilde{j}}_{sx} \right) ^{n}_{i{-}\frac{3}{2},j} + \mathrm {Re}L{\overline{{\widetilde{c}}}}^{\,n} n_x(1-\sigma ) \llbracket \,{\widetilde{p}}\,\rrbracket ^n }{\left( \omega ^{-} + 0.5 \right) \Delta {\widetilde{x}}} - \frac{\Delta {\widetilde{t}}}{2} \left( \delta _{{\widetilde{y}}}\, {\widetilde{j}}_{sy} \right) ^{n}_{i{-}1,j} \ , \end{aligned}$$

  (11)

• at the boundary cell (i, j)

  $$\begin{aligned}{\widetilde{c}}^{\,n+1}_{i,j} + \frac{\Delta {\widetilde{t}}}{2} \frac{ \left( {\widetilde{j}}_{sx} \right) ^{n+1}_{i+\frac{1}{2},j} - \left( ({\widetilde{c}}^{\,+})_{i,j}^{n+1} {\widetilde{u}}_m^{\,n} + \mathrm {Re}L{\overline{{\widetilde{c}}}}^{\,n} n_x \mu _{\mathrm{p}}(\mu _{\mathrm{L}}- \sigma ) \llbracket \,{\widetilde{c}}\,\rrbracket ^{n+1} \right) }{\left( \omega ^{+} + 0.5 \right) \Delta {\widetilde{x}}} + \frac{\Delta {\widetilde{t}}}{2} \left( \delta _{{\widetilde{y}}} {\widetilde{j}}_{sy} \right) ^{n+1}_{i,j} \nonumber \\\quad = {\widetilde{c}}^{\,n}_{i,j} - \frac{\Delta {\widetilde{t}}}{2} \frac{ \left( {\widetilde{j}}_{sx} \right) ^{n}_{i+\frac{1}{2},j} - \left( {\widetilde{j}}_{sx}^+ \right) ^{n} - \mathrm {Re}L{\overline{{\widetilde{c}}}}^{\,n} n_x(1-\sigma ) \llbracket \,{\widetilde{p}}\,\rrbracket ^n }{\left( \omega ^{+} + 0.5 \right) \Delta {\widetilde{x}}} - \frac{\Delta {\widetilde{t}}}{2} \left( \delta _{{\widetilde{y}}} {\widetilde{j}}_{sy} \right) ^{n}_{i,j} \ . \end{aligned}$$

  (12)

Note that to linearise the term ${\overline{{\widetilde{c}}}}\llbracket \,{\widetilde{c}}\,\rrbracket$ in Eq. (10), the average concentration on the membrane ${\overline{{\widetilde{c}}}}$ is evaluated at the time level n (i.e. it is treated as a known value). By linearly extrapolating for ${\widetilde{c}}^{\,\mp }$ from the respective sides, Eqs. (11) and (12) constitute a closed set of simultaneous equations for c, which is solved using the pre-conditioned BiCGSTAB method.

Fig. 2

Schematic of a corrugated permeable membrane travelling in a parallel channel at a constant speed $U_0$ in the $+x$ direction

4 Validation

4.1 Problem statement

We consider a problem with a corrugated permeable membrane translating in a region between no-slip parallel plates, as illustrated in Fig. 2. The geometry and motion of the corrugation are given by $h(x,t) = h_0 + h_0 \delta \cos [k(x-U_0 t)]$, where $h_0 = H_0/2$ is the half channel height, $k=2\pi /l$ is the wavenumber, $\delta$ is a dimensionless parameter between 0 and 1, and $U_0$ is a constant velocity. Periodic boundary conditions are imposed on the left and right boundaries. No solvent permeation is considered on the top and bottom boundaries, whereas the following boundary concentrations are prescribed on the top and bottom boundaries, respectively: $c_{\mathrm{top}}(x) = c_0 [1 + \sin (kx)]$ and $c_{\mathrm{bot}}(x) = c_0 [1 - \sin (kx)]$. In this problem, $2c_0$ is taken as the reference concentration (C).

The lower and upper regions of the membrane are denoted as $\Omega _1$ and $\Omega _2$, respectively. Hereafter, the variables in those regions are distinguished by the subscripts 1 and 2, and jump values are defined on the membrane as $\llbracket \,\varphi \,\rrbracket =\varphi _1-\varphi _2$.

The narrow gaps between the corrugation and the flat walls induce the lubrication. Therefore, in this problem, the permeate fluxes are driven by the hydrostatic pressure difference developed by lubrication, and at the same time, osmotic pressure difference owing to the difference in concentration either promotes or impede the permeations by Eq. (4).

4.2 Analytical models of lubrication-induced permeate fluxes

By introducing $\epsilon$ as the ratio of the channel width to the channel length $h_0/l$, we assume $\epsilon \ll 1$ and $\epsilon \mathrm {Re}\ll 1 \left( \mathrm {Re}=\rho _{\mathrm{f}}U_0h_0/\mu _{\mathrm{f}}\right)$. The hydrostatic pressure is described by the Reynolds lubrication equation for a narrow gap between the corrugation and the flat plate. We further assume an infinitesimal limit for $L_{\mathrm{p}}$ to isolate the effect of permeability to develop an asymptotic analytical model. Then, the pressure in $\Omega _1$ is expressed as follows (Takeuchi and Gu 2019; Tazaki et al. 2020):

$$\begin{aligned} & p_1^{(0)}(x,t) = -\frac{6}{\pi \epsilon ^2}\frac{U_0\mu _{\mathrm{f}}}{l} \frac{\delta }{2+\delta ^2}\\ &\quad\frac{ \big \{2 + \delta \cos [k(x-U_0 t)]\big \} \sin [k(x-U_0 t)]}{\big \{(1 + \delta \cos [k(x-U_0 t)]\big \}^2} \ , \end{aligned}$$

(13)

where the superscript (0) represents the 0th-order pressure (i.e. the pressure satisfying the Reynolds lubrication equation). Owing to the symmetry of the computational domain, the pressure in $\Omega _2$ is given as $p_2^{(0)}(x,t) = p_1^{(0)}(x+\pi /k,t)$.

When $\epsilon \ll 1$ and $\epsilon \mathrm {Pe}\ll 1 \left( \mathrm {Pe}=U_0h_0/D\right)$, the mass transport of the solute is dominated by diffusion in the y direction. Then, the wall-normal concentration distribution is assumed to be uniform, and the wall-tangential distribution is approximated to be equivalent to the boundary concentrations: $c_1(x, {}^\forall y) \simeq c_{\mathrm{bot}}(x)$ and $c_2(x, {}^\forall y) \simeq c_{\mathrm{top}}(x)$.

The pressure jump, concentration jump and the average concentration on the membrane are obtained as follows:

$$\begin{aligned}\llbracket \,p^{(0)}\,\rrbracket = p_1^{(0)}(x,t) - p_2^{(0)}(x,t) = -\frac{24}{\pi \epsilon ^2}\frac{U_0\mu _{\mathrm{f}}}{l}\frac{\delta }{2+\delta ^2}\frac{\sin {[k(x-U_0 t)]}}{(1-\delta ^2\cos ^2{[k(x-U_0 t)]}^2}, \end{aligned}$$

(14a)

$$\begin{aligned}\llbracket \,c\,\rrbracket &= c_1(x) - c_2(x) = c_{\mathrm{bot}}(x) - c_{\mathrm{top}}(x)\\& = - 2 c_0 \sin {(kx)}= - C \sin {(kx)}, \end{aligned}$$

(14b)

$$\begin{aligned}{\overline{c}} = \frac{c_1(x) + c_2(x)}{2}=c_0 = \frac{C}{2} \ . \end{aligned}$$

(14c)

Then, the y components of the asymptotic permeate fluxes are approximated as follows:

$$\begin{aligned}&\varvec{J}_{v}(x,t) \cdot \varvec{e}_y = - \frac{24 LU_0}{\pi \epsilon }\frac{\delta }{2+\delta ^2}\frac{\sin {[k(x-U_0 t)]}}{\left( 1-\delta ^2\cos ^2{[k(x-U_0 t)]} \right) ^2} + \mathrm {Re}LU_0 \sigma \mu _{\mathrm{p}}\sin {(kx)} \ , \end{aligned}$$

(15a)

$$\begin{aligned}&\varvec{J}_{s}(x,t) \cdot \varvec{e}_y = -\frac{C}{2} \left\{ (1-\sigma )\frac{24 LU_0}{\pi \epsilon }\frac{\delta }{2+\delta ^2}\frac{\sin {[k(x-U_0 t)]}}{\left( 1-\delta ^2\cos ^2{[k(x-U_0 t)]}\right) ^2} + \mathrm {Re}LU_0 (\mu _{\mathrm{L}}-\sigma ) \mu _{\mathrm{p}}\sin {(kx)} \right\} \ . \end{aligned}$$

(15b)

Here, the permeability is normalised as $L\, (=L_{\mathrm{p}}\mu _{\mathrm{f}}/H_0)$.

To improve the prediction of lubrication-induced flow over a larger $\epsilon$ range (i.e. $\epsilon \lesssim 1$), a higher-order lubrication model (Takeuchi and Gu 2019) is applied to describe the wall-normal distribution of pressure. Takeuchi et al. (2021) showed the 2nd-order pressure for the corrugated membrane as:

$$\begin{aligned} p_1^{(2)*} = -4\pi \frac{U_0\mu _{\mathrm{f}}}{H_0} \frac{y}{h_0} \frac{\epsilon ^2 \delta }{2+\delta ^2}\\ \quad\times \frac{ \left( 4 \delta ^2+\left( \delta ^2+2\right) \delta \cos (k x^*)-1\right) (2 \delta \cos (k x^*)+2 -3 y/h_0)}{ (1+\delta \cos (k x^*))^4} \sin (k x^*) \ . \end{aligned}$$

(16)

where the superscript “$*$” represents the value observed on the frame fixed at the membrane; $x^* = x - U_0t$. Using $p_2^{(2)*}(x^*)= p_1^{(2)*}(x^*+\pi /k)$, the pressure jump (with the higher-order correction) is given as $\llbracket \,p^*\,\rrbracket =\llbracket \,p^{(0)*}+p^{(2)*}\,\rrbracket$, and the corresponding fluxes (denoted as $\varvec{J}_v^{(0+2)*}$ and $\varvec{J}_s^{(0+2)*}$) are calculated from Eq. (4). The explicit forms of the fluxes are not presented here because of the long mathematical expressions.

4.3 Simulation conditions

The simulation parameters are set as follows: the channel length $l=5H_0$, the grid resolution is $H_0/\Delta {x}=40$, the time increment $\Delta t/(H_0/U_0)=5\times 10^{-5}$, $\sigma =0.5$, and $\mu _{\mathrm{L}}= \mu _{\mathrm{p}}= 1$. The effect of the amplitude of corrugation is investigated at $\delta = 0.1$ unless specified otherwise. At the above spatial resolution, the y variation of the corrugation $2h_0\delta$ is covered by 4 grid points, which is sufficient from our previous study (Takeuchi et al. 2018). The Reynolds number and the Peclet number are fixed at $\mathrm {Re}=0.5$ and $\mathrm {Pe}=0.5$, respectively. The permeability $L$ is varied in the following range: $L=10^{-5},10^{-4},10^{-3},10^{-2}, 10^{-1},10^{0}$.

By substituting the above values of $\mu _{\mathrm{L}}$, $\mu _{\mathrm{p}}$ and $\sigma$ into Eq. (5), the permeate fluxes are simplified as follows:

$$\begin{aligned}\frac{\varvec{J}_v}{ U_0} = \mathrm {Re}L\frac{2\llbracket \,{\widetilde{p}}\,\rrbracket -\llbracket \,{\widetilde{c}}\,\rrbracket }{2} \varvec{n}, \end{aligned}$$

(17a)

$$\begin{aligned}\frac{\varvec{J}_s}{CU_0} = \mathrm {Re}L\,{{\overline{{\widetilde{c}}}}} \,\frac{\llbracket \,{\widetilde{p}}\,\rrbracket+\llbracket \,{\widetilde{c}}\,\rrbracket }{2} \varvec{n} \ . \end{aligned}$$

(17b)

Considering that the base functions of $\llbracket \,p\,\rrbracket$ and $\llbracket \,c\,\rrbracket$ are $\sin [k(x-U_0t)]$ and $\sin (kx)$ (see Eq. 14), there are the moments when $\llbracket \,p\,\rrbracket$ and $\llbracket \,c\,\rrbracket$ weaken and strengthen each other; at $t/(H_0/U_0)=(2m-1)\pi /kU_0$ $(m=1,2,\cdots )$, $\llbracket \,p\,\rrbracket$ and $\llbracket \,c\,\rrbracket$ weaken and strengthen each other for $\varvec{J}_v$ and $\varvec{J}_s$, respectively, while at $t/(H_0/U_0)=2m\pi /kU_0$, $\llbracket \,p\,\rrbracket$ and $\llbracket \,c\,\rrbracket$ strengthen and weaken each other for $\varvec{J}_v$ and $\varvec{J}_s$, respectively. In the following, $m=3$ is taken (i.e. $t/(H_0/U_0)=12.5$ and 15.0). At $t/(H_0/U_0)=12.5$, the solvent and solute permeations reach the respective minimum and maximum strengths, while the permeations of the solvent and solute are strongest and weakest, respectively, at $t/(H_0/U_0)=15.0$.

4.4 Simulation results

Figures 3, 4, 5 and 6 show the pressure and concentration fields at time $t/(H_0/U_0)=12.5$ and 15.0 for the following permeabilities: $L=10^{-5},10^{-2},10^{-1},10^{0}$. The results are visualised on the coordinate system fixed on the corrugated membrane (i.e. $x^*$ and y).

In Figs. 3 and 5, the pressure distributions at $L=10^{-5}$ and $10^{-2}$ tend to be insensitive to y in both the upper and lower regions (as predicted by lubrication theory). The y-insensitive distribution is more pronounced at $t/(H_0/U_0)=12.5$ (Fig. 3a, b) than the cases at $t/(H_0/U_0)=15.0$ (Fig. 5) because for this weak $\varvec{J}_v$ condition, the apparent permeability $L$ is supposed to be low. However, for $L=10^{-1}$ and $10^{0}$ at both times (Figs. 3c, d, 5c, d), stronger y-dependent trend is observed in the pressure field; in particular, high- and low-pressure regions near the flat plates are highlighted at $t/(H_0/U_0)=15.0$. The result indicates the presence of a regime that is not described by Reynolds lubrication equation, which is referred to as the non-Reynolds lubrication regime (Takeuchi et al. 2021; Takeuchi and Gu 2019), hereafter.

Fig. 3

Longitudinal distributions of pressure p at $\delta =0.1$ for the four $L$ values at the time of the maximum solvent flux $t/(H_0/U_0)=12.5$. The pressure is shown normalised by $\rho _{\mathrm{f}}U_0^2$

Fig. 4

Longitudinal distributions of concentration c at $\delta =0.1$ for the four $L$ values at the time of the minimum solute permeate flux $t/(H_0/U_0)=12.5$

Fig. 5

Longitudinal distributions of pressure p at $\delta =0.1$ for the four $L$ values at the time of the minimum solvent permeate flux $t/(H_0/U_0)=15.0$. The pressure is shown normalised by $\rho _{\mathrm{f}}U_0^2$

Fig. 6

Longitudinal distributions of concentration c at $\delta =0.1$ for the four $L$ values at the time of the maximum solute permeate flux $t/(H_0/U_0)=15.0$

The concentration field tends to show the y-insensitive distributions in all cases (Figs. 4,  6), which strongly reflects the effect of the small Peclet number and small aspect ratio of the channel. As a result, $\displaystyle { \max _{x^*}\llbracket \,c\,\rrbracket/C }$ is nearly equal to 1, which will be mentioned again.

Fig. 7

Comparison between the numerical result and the analytical model for the solvent permeate flux $(\varvec{J}_{v} \cdot \varvec{e}_y)$ in the y-direction at $\delta =0.1$ for the four $L$ cases at the two instants when $\varvec{J}_{v}$ is a weakest and b strongest

Fig. 8

Comparison between the numerical result and the analytical model for the solute permeate flux $(\varvec{J}_{s} \cdot \varvec{e}_y)$ in the y-direction with $\delta =0.1$ for the four $L$ cases at the instants when $\varvec{J}_{s}$ is a strongest and b weakest

To visualise the permeate fluxes in the $L$ range of the y-insensitive pressure distribution, Figs. 7 and 8 show the profiles of the y components of the permeate fluxes of the solvent and solute in the range of $L$ between $10^{-5}$ and $10^{-2}$. The dashed lines represent $\varvec{J}_v^{(0)}$ (Fig. 7) and $\varvec{J}_s^{(0)}$ (Fig. 8), and the solid lines represent $\varvec{J}_v^{(0+2)}$ (Fig. 7) and $\varvec{J}_s^{(0+2)}$ (Fig. 8). The vertical axis is the permeate flux divided by $L$. Note that $\varvec{J}_v/L$ and $\varvec{J}_s/CL$ are essentially the same as $2\llbracket \,p\,\rrbracket -\llbracket \,c\,\rrbracket$ and $\llbracket \,p\,\rrbracket +\llbracket \,c\,\rrbracket$, respectively, which facilitate the study of the convergence of the discontinuities of pressure and concentration in the limit of $L\rightarrow 0$. The figures show that the numerical results exhibit converging trends towards the asymptotic solutions at $L\rightarrow 0$, and $\varvec{J}_v^{(0+2)}$ and $\varvec{J}_s^{(0+2)}$ describe the permeation physics better than $\varvec{J}_v^{(0)}$ and $\varvec{J}_s^{(0)}$. The results indicate that the effect of the non-Reynolds lubrication regime needs to be corrected.

Fig. 9

Deviation of the numerical result of the a solvent and b solute permeate fluxes in the $\text {L}^2$ norm from the analytical models of $\varvec{J}_v^{(0+2)}\cdot \varvec{e}_y$ and $\varvec{J}_s^{(0+2)}\cdot \varvec{e}_y$, respectively, plotted against $L$ for two amplitude parameters $\delta$

Figure 9 summarises the convergence of the numerical permeate fluxes towards $\varvec{J}_v^{(0+2)}\cdot \varvec{e}_y$ and $\varvec{J}_s^{(0+2)}\cdot \varvec{e}_y$ in the $\text {L}^2$ norm for $\delta =0.1$ and 0.5. The convergence at about the first-order rate of $L$ is observed for both $t/(H_0/U_0)=12.5$ and 15.0. In addition, the error levels are improved with decreasing $\delta$. The results indicate the validity of the developed simulation method for lubrication-induced permeation.

Fig. 10

Numerical result of the solvent permeate flux $(\varvec{J}_{v} \cdot \varvec{e}_y)$ at $\delta =0.1$ for $L=10^{-1},10^{0}$ at the same instants as in Fig. 7

Fig. 11

Numerical result of the solute permeate flux $(\varvec{J}_{s} \cdot \varvec{e}_y)$ at $\delta =0.1$ for $L=10^{-1},10^{0}$ at the same instants as in Fig. 8

At the end of this section, a characteristic situation for the balance between $\llbracket \,p\,\rrbracket$ and $\llbracket \,c\,\rrbracket$ is investigated. Figures 10 and 11 show the numerical results of the y components of the permeate fluxes of the solute and solvent at $L=10^{-1}$ and $10^{0}$. From the graphs, the distribution of the solvent flux (Fig. 10) changes little with time, whereas the phase difference in the distribution of the solute permeate flux is significant (Fig. 11). This is because the maximum concentration jump on the membrane is approximately unity in all $L$ cases for the small Peclet number employed, while the maximum pressure jump varies as $\llbracket \,{\widetilde{p}}\,\rrbracket =5.0, 3.0, 0.8, 0.6$ for $L=10^{-5}, 10^{-2}, 10^{-1}, 10^{0}$, respectively. Therefore, at $L=10^{-1}$ and $10^{0}$, the concentration jump on the membrane is more influential than the pressure jump, and the diffusion of the solute (i.e. osmotic pressure difference) has the predominant effect, which renders an interesting implication of the lubrication-induced permeation of both the solute and solvent at large permeabilities.

In summary, the pressure-dominant permeation in a small $L$ range (i.e. $L=10^{-5}$ and $10^{-2}$; ideally $L\rightarrow 0$) shows good agreement of the numerical permeate fluxes with the higher-order models $\varvec{J}_v^{(0+2)}$ and $\varvec{J}_s^{(0+2)}$, whereas over a finite $L$ range, a different type of higher-order model may be necessary to describe the nonlinearity owing to the above comparable effects of $\llbracket \,p\,\rrbracket$ and $\llbracket \,c\,\rrbracket$.

Fig. 12

Schematic of a transport problem of solute from a circular permeable membrane between parallel plates moving at a constant speed $U_0$

5 Lubrication-induced mass transport through circular membrane

As a simplified model of mass transport inside and outside a biological cell in a capillary, we set up a problem involving a circular permeable membrane ($D_{\mathrm{p}}$ in diameter) between parallel plates ($H_0$ in distance), as shown in Fig. 12. The top and bottom walls are moved in the $-x$ direction at a velocity $U_0$, and the non-deformable circular membrane is fixed at the centre level of the channel. No-slip and impermeable conditions are imposed on the top and bottom walls, and periodic boundary conditions are imposed on the left and right sides. The initial concentrations $C_0$ and 0 are uniformly given inside and outside the circular membrane ($\Omega _1$ and $\Omega _2$), respectively. Then, to analyse the effect of lubrication on permeation, simulations are carried out by varying the distance $H_0$ between the parallel plates.

In the simulations, the domain length and the grid spacing are set at $l=3D_{\mathrm{p}}$ and $\Delta {x}=D_{\mathrm{p}}/40$, respectively, and the time increment is fixed at $\Delta t/(D_{\mathrm{p}}/U_0)=1\times 10^{-5}$. The dimensionless numbers for the permeate flux models are as follows: $\mathrm {Re}=\rho _{\mathrm{f}}U_0D_{\mathrm{p}}/\mu _{\mathrm{f}}=1$, $\mathrm {Pe}=U_0D_{\mathrm{p}}/D=1$, $\sigma =0.5$ and $\mu _{\mathrm{L}}= \mu _{\mathrm{p}}= 1$. The permeability is set to $L=10^{-2}$ unless specified otherwise.

The channel width is varied in the following range: $H_0/D_{\mathrm{p}}=1.5, 2.0, 3.0$, which correspond to the narrowest width (i.e. the distance between the lowest/highest point of the circular membrane and the bottom/top wall) having a value of $(H_0-D_{\mathrm{p}})/2=0.25D_{\mathrm{p}}, 0.50D_{\mathrm{p}}, 1.0D_{\mathrm{p}}$, respectively. Although the case of the smallest gap is already out of the range of the ideal lubrication condition $\epsilon \ll 1$, our previous studies (Takeuchi and Gu 2019; Takeuchi et al. 2021) showed that the lubrication phenomenon can still be described by including a higher-order correction for a non-negligible gap. Note that the effective longitudinal length scales of the lubrication region can be estimated to be $\sqrt{\varepsilon }D_p$ (Takeuchi et al. 2021), and the effective scales for the above three cases are $0.5D_p$, $0.71D_p$, and $1.0D_p$, respectively, which are sufficiently smaller than the domain length l, indicating that the permeation induced by lubrication pressure is influenced little by the mirror images of the circular membranes due to the periodic boundary condition.

The pressure field in the entire domain at the time $t/(D_{\mathrm{p}}/U_0)=3.0\ (=l/D_{\mathrm{p}},$ see Fig. 12) is shown in Fig. 13. The figure shows that lubrication in the wall-membrane gap causes an increase in the pressure in the right-half side of $\Omega _2$ with decreasing wall-membrane distance.

The concentration fields in $\Omega _1$ and $\Omega _2$ are visualised in Fig. 14. Recalling that the initial concentration is zero in $\Omega _2$, more solute goes out of the membrane as the wall-membrane gap becomes narrower. The x variations of the concentration field in $\Omega _1$ are evident for all cases, and the largest gradient for the case of $H_0/D_{\mathrm{p}}=1.5$ indicates that the unsteady transport is predominant. However, the pressure shows different trends. The pressures in $\Omega _1$ remain nearly uniform for all the cases, whereas the pressure discontinuity varies along the membrane and takes the maximum and minimum values in the near-wall regions. These positions of the largest $|\llbracket \,p\,\rrbracket |$ are predicted by the higher-order lubrication model (Takeuchi and Gu 2019), and it is worth noting that the stagnant point (i.e. the right-hand side intersection of the centreline and membrane) does not show the largest pressure for lubrication flows.

Therefore, the strong variations of c and p in $\Omega _1$ and $\Omega _2$, respectively, indicate characteristic distributions of the lubrication-induced permeate fluxes along the membrane. In the following, the effect of lubrication on permeate fluxes of the solvent and solute is discussed.

Fig. 13

Pressure field p(x, y) at $t/(D_{\mathrm{p}}/U_0)=3.0$. The pressure is shown normalised by $\rho _{\mathrm{f}}U_0^2$. The colour ranges are varied to adjust the maximum and minimum values of the respective cases

Fig. 14

Concentration field c(x, y) inside and outside the circular membrane at $t/(D_{\mathrm{p}}/U_0)=3.0$ for three $H_0/D_{\mathrm{p}}$ values. a, b $H_0/D_{\mathrm{p}}=1.5$, c, d $H_0/D_{\mathrm{p}}=2.0$, e, f $H_0/D_{\mathrm{p}}=3.0$

Figure 15 shows the velocity field and solvent fluxes through the membrane at time $t/(D_{\mathrm{p}}/U_0)=3.0$. Larger permeate fluxes are observed for smaller membrane-wall distance cases, indicating the active exchange of solvent between $\Omega _1$ and $\Omega _2$. The velocity in $\Omega _1$ is especially large for $H_0/D_{\mathrm{p}}=1.5$ (Fig. 15a). Considering that the concentration field could be strongly influenced by the local velocity, it may be interesting to decompose the mass flux into the convective and diffusive fluxes of the solute.

Figures 16 and 17 compare the fluxes of convection, diffusion, and permeation of the solute inside and outside the membrane, respectively, at the same instant as Fig. 15. In the case of the small gap ($H_0/D_{\mathrm{p}}=1.5$), the solute transport in $\Omega _1$ is dominated by convection (Fig. 16a), which indicates that the solvent permeation by lubrication enhances the transport of the solute in $\Omega _1$. Interestingly, diffusion is predominant in $\Omega _2$ (Fig. 17a), particularly in the narrowest gap regions between the membrane and flat walls. This is because the value of c varies largely in the x direction in this region, as shown in Fig. 14b.

For the other distance cases (i.e. $H_0/D_{\mathrm{p}}=2.0$ and 3.0), the effects of convective and diffusive fluxes are comparable, indicating that the lubrication is less active.

Fig. 15

Solvent velocity and the permeate flux at $t/(D_{\mathrm{p}}/U_0)=3.0$. Both are shown normalised by $U_0$

Fig. 16

Diffusive, convective and permeate fluxes of the solute in a circular membrane at $t/(D_{\mathrm{p}}/U_0)=3.0$. The fluxes are shown normalised by $C_0U_0$

Fig. 17

Diffusive, convective and permeate fluxes of the solute outside the circular membrane at $t/(D_{\mathrm{p}}/U_0)=3.0$. The fluxes are shown normalised by $C_0U_0$

Fig. 18

Temporal evolution of the total amount of mass inside the membrane of $L=10^{-2}$

Figure 17 shows that the outgoing flux (on the left-hand side of the membrane) is greater than the incoming flux for all of the $H_0/D_{\mathrm{p}}$ cases. To compare the effect of the wall-membrane distance on the net permeate flux, the time evolution of the amount of solute inside the circular membrane is given as follows:

$$\begin{aligned} \frac{1}{\pi C_0D_{\mathrm{p}}^2/4}\int _{\Omega _1} c(x,y,t) \ \mathrm{d}{x}\mathrm{d}{y} \ \end{aligned}$$

(18)

and is plotted in Fig. 18. The graph shows that for $H_0/D_{\mathrm{p}}=1.5$, the solute permeation from the inside to the outside is faster than the other two cases. As observed previously, the solute permeation is enhanced by the lubrication; the increase in lubrication pressure at the near-wall region of the membrane (Fig. 13a) causes a local increase in the solvent permeate flux into the membrane, resulting in the enhancement of convection and permeation of solute.

Table 1 Slope of Eq. (18) at $t=0$ for $H_0/D_{\mathrm{p}}=1.5$

Table 1 summarises the values of the slope (i.e. time derivative) of Eq. (18) evaluated at $t=0$ for the six $L$ cases at $H_0/D_{\mathrm{p}}=1.5$. The initial decrease rates of the mass exhibit approximately proportional to $-L$. Although Eq. (15) is for different geometry, the solute permeation being proportional to $L$ may support the general trend of lubrication-induced permeation.

6 Conclusion

To study the effect of lubrication on the permeations of solute and solvent through membrane, DF-IB method was proposed with a permeable membrane, and the contribution of a higher-order mode of flux to permeation was highlighted through mathematical modelling.

In the numerical study, the permeate flux models for the solute and solvent were incorporated into the DF-IB method by considering direct discretisation and the consistent coupling (between the incompressible velocity and pressure fields) in the immediate vicinity of the membrane, and discretised equations of the pressure Poisson equation and convective-diffusion equation for the solute were expressed with the discontinuities at the membrane.

The proposed method was validated by showing the convergence of the numerical result to the higher-order analytical models of permeate fluxes (of solute and solvent) in a lubrication-dominated flow field for a problem involving the movement of corrugated membrane between parallel flat plates. The lubrication was found to promote the permeation of the solvent and solute. However, this effect depends on the permeability of the membrane, and for highly permeable membranes, the nonlinearity in the solute permeation through the membrane becomes non-negligible.

As a model of mass release from a biological cell flowing in capillary, a flow problem with a circular membrane placed between parallel plates was simulated, and the effect of lubrication was investigated by varying the distance between the membrane and plates. The pressure discontinuity on the membrane in a near-wall region is larger than at the stagnant point, which highlights the effect of lubrication on the permeate fluxes of solute and solvent. In particular, for a small gap case, the solute transport was dominated by convection inside the circular membrane and by diffusion outside. The temporal evolution of the concentration in the circular membrane indicated the enhancement of solute release by lubrication.

In the present study, the circular membrane was assumed to be non-deformable, and the effect of lubrication on the membrane permeation was found to be predominant in a narrow channel. For general deformable membranes, considering that lubrication forces promote the deformation of the membrane (Secomb et al. 2001), mass transport phenomena may be influenced by time-dependent variations of the membrane shape. Therefore, studies to investigate the effect of deformation under lubrication on mass transport through membranes are the subject of future work by the present authors.