A Model for Stokes Flow in Domains with Permeable Boundaries

Abstract

We derive a new computational model for the simulation of viscous incompressible flows bounded by a thin, flexible, porous membrane. Our approach is grid-free and models the boundary forces with regularized Stokeslets. The flow across the porous membranes is modeled with regularized source doublets based on the notion that the flux velocity across the boundary can be viewed as the flow induced by a fluid source/sink pair with the sink on the high-pressure side of the boundary and magnitude proportional to the pressure difference across the membrane. Several validation examples are presented that illustrate how to calibrate the parameters in the model. We present an example consisting of flow in a closed domain that loses volume due to the fluid flux across the permeable boundary. We also present applications of the method to flow inside a channel of fixed geometry where sections of the boundary are permeable. The final example is a biological application of flow in a capillary with porous walls and a protein concentration advected and diffused in the fluid. In this case, the protein concentration modifies the pressure in the flow, producing dynamic changes to the flux across the walls. For this example, the proposed method is combined with finite differences for the concentration field.

1. Introduction

We present a new model for a type of problem in which a thin flexible porous boundary, such as a membrane, interacts with an incompressible viscous fluid. Problems of this nature arise in numerous biological applications such as in modeling cell membranes [1,2], microvasculature [3,4], and specialized anatomical structures such as the renal tubule [5,6,7]. Membranes are used in numerous medical applications such as drug delivery [2,8], biosensors [9,10] and tissue engineering [11]. Other more general engineering applications include membrane distillation processes [12], water filtration [8], and wastewater treatment [13]. There are widespread applications of microfiltration processes that take advantage of selective tuning of membrane permeability [14,15], and many theoretical and experimental investigations have focused on the swelling and permeability of gel dispersions in the microfluidic domain [16]. Theoretical studies have investigated the dynamics of particle interactions with deformable permeable membranes [17] as well as the behaviors of these membranes in complex flows [18].

There are numerical methods available for fluid flows bounded by a porous domain where the porous boundary does not deform [19,20]. Our goal is to address the case when an infinitely thin, flexible porous membrane divides the fluid domain. Some of the most closely related methods to our approach are described in [21,22]. Stockie [21] uses the immersed boundary method to model the Navier–Stokes equations with immersed elastic membranes and incorporates the porosity of the membrane by introducing a flux velocity across the boundary proportional to the hydrostatic pressure jump. A similar strategy is employed in [23] where the immersed boundary method and Darcy’s law are used to simulate the aerodynamics of a 2-D parachute, and in [24] to investigate effects of porosity on the motion of a flapping filament. Layton [22] uses a Cartesian grid method combined with Mayo’s technique [25], similar to the immersed interface method [26], to incorporate jump conditions into the finite-difference stencils. This formulation also incorporates a solute concentration into the model so that the flux velocity across the membrane can be a function of the hydrostatic pressure jump or the solute concentration differences across the boundary, a condition often encountered in biological settings. Using a coupled immersed boundary–lattice Boltzmann method, Pepona calculates fluid flow through and around a volumetric porous body that deforms elastically [27].

Exact solutions of Stokes flow through a permeable tube have been considered in several works. In particular, the seminal work of Berman [28] considered laminar Stokes flow near a porous boundary and has since been expanded to model flows in permeable cylinders [29,30,31]. These models assume low permeability and a small ratio between transverse and axial velocities. More recent studies have relaxed this assumption by deriving an exact solution for Stokes flow through a pipe of arbitrary permeability [19]. This work was later expanded to consider Stokes flow through a pipe of temporally dynamic radius to simulate “pumping” mechanisms encountered in numerous biological systems [20]. These analytical solutions assume that the position of the porous boundary is independent of the flow, as such they are useful for validation of numerical methods with a similar assumption.

Here, we present a grid-free method to model the motion of Stokes flows bounded by a flexible, elastic, porous membrane. The formulation of the model, as presented in the next section, applies to both two and three dimensions. The rest of the article focuses on the implementation and examples in two dimensions. The approach is based on the method of regularized Stokeslets, augmented with regularized source doublets on the porous boundaries. The latter is based on the notion that the flux velocity across the boundary can be considered as the flow induced by a fluid source/sink pair with the source on the high-pressure side of the boundary and magnitude proportional to the pressure difference across the membrane.

2. The Model Formulation

We consider a closed elastic membrane $\Gamma$ that encloses a two-dimensional domain $\mathsf{\Omega }$. The membrane is given parametrically by $\mathbf{X}\left(s\right)$ where s is arc length, and is permeable to the fluid, whose motion is assumed to be described by the Stokes equations in ${\mathbb{R}}^2/\Gamma$. Since the membrane is elastic, it supports a force density given by $\mathbf{f}\left(s\right)$. The equations of motion are

$\begin{array}{ccc}\hfill \mathsf{\mu }\Delta \mathbf{u}=\nabla p-\mathbf{F},\nabla ·\mathbf{u}=0& & \mathrm{i}\mathrm{n}{\mathbb{R}}^2/\Gamma \hfill \\ \hfill \mathbf{F}\left(\mathbf{x}\right)={\int }_{\Gamma }\mathbf{f}\left(s\right)\delta (\mathbf{x}-\mathbf{X}\left(s\right))ds& & \mathrm{f}\mathrm{o}\mathrm{r}\mathbf{x}\in {\mathbb{R}}^2\hfill \end{array}$

A fluid particle moves with velocity $\mathbf{u}$ so that $d\mathbf{x}/dt=\mathbf{u}\left(\mathbf{x}\right)$. Due to the permeability, the membrane moves with a different velocity field, namely $d\mathbf{X}/dt=\mathbf{u}\left(\mathbf{X}\right)+\mathbf{U}$, where $\mathbf{U}$ is normal to the membrane and depends on the jump across the membrane of pressure or some other quantity.

2.1. A Different Formulation of the Problem

Given that the velocity of the fluid and the membrane differ by $\mathbf{U}$, we can write a general expression for the velocity field in the fluid domain including material points on the membrane as

$\frac{d\mathbf{x}}{dt}=\mathbf{u}\left(\mathbf{x}\right)+{\int }_{\Gamma }\mathbf{U}\left(s\right)\delta (\mathbf{x}-\mathbf{X}\left(s\right))ds$

where the last term on the right represents the membrane velocity. This term can also be interpreted as a velocity field whose divergence is ${\int }_{\Gamma }\mathbf{U}\left(s\right)·\nabla \delta (\mathbf{x}-\mathbf{X}\left(s\right))ds$. This leads to the formulation

(1)$\begin{array}{ccc}\hfill \mathsf{\mu }\Delta \mathbf{u}=\nabla p-\mathbf{F},\nabla ·\mathbf{u}=A& & \mathrm{i}\mathrm{n}{\mathbb{R}}^2\hfill \end{array}$

(2)$\begin{array}{ccc}\hfill \mathbf{F}\left(\mathbf{x}\right)={\int }_{\Gamma }\mathbf{f}\left(s\right)\delta (\mathbf{x}-\mathbf{X}\left(s\right))ds& & \mathrm{f}\mathrm{o}\mathrm{r}\mathbf{x}\in {\mathbb{R}}^2\hfill \end{array}$

(3)$\begin{array}{ccc}\hfill A\left(\mathbf{x}\right)={\int }_{\Gamma }\mathbf{U}\left(s\right)·\nabla \delta (\mathbf{x}-\mathbf{X}\left(s\right))ds& & \mathrm{f}\mathrm{o}\mathrm{r}\mathbf{x}\in {\mathbb{R}}^2\hfill \end{array}$

2.2. The Proposed Regularized Formulation of the Problem

We use the linearity of the Stokes equations to separate the formulation in Equations (1)–(3) into two different problems whose solutions are added together. Let $(\mathbf{u},p)$ satisfy the incompressible Stokes equations with regularized external forcing

(4)$\begin{array}{cc}\hfill \mathsf{\mu }\Delta \mathbf{u}=\nabla p-\mathbf{F},\nabla ·\mathbf{u}=0& \mathrm{i}\mathrm{n}{\mathbb{R}}^2/\Gamma \hfill \\ \hfill \mathbf{F}\left(\mathbf{x}\right)={\int }_{\Gamma }\mathbf{f}\left(s\right){\varphi }_{\delta }(\mathbf{x}-\mathbf{X}\left(s\right))ds& \mathrm{f}\mathrm{o}\mathrm{r}\mathbf{x}\in {\mathbb{R}}^2\hfill \end{array}$

and let $(\mathbf{v},q)$ satisfy the Stokes equations without external forcing but with divergence given by a regularized source doublet distribution on the membrane

(5)$\begin{array}{ccc}\hfill \mathsf{\mu }\Delta \mathbf{v}=\nabla q,\nabla ·\mathbf{v}=A& \hfill & \hfill \mathrm{i}\mathrm{n}{\mathbb{R}}^2\\ \hfill A\left(\mathbf{x}\right)=-{\int }_{\Gamma }\frac{b\left(s\right)}{\mathsf{\mu }}\widehat{\mathbf{n}}\left(s\right)·\nabla {\varphi }_{\delta }(\mathbf{x}-\mathbf{X}\left(s\right))ds& \hfill & \hfill \mathrm{f}\mathrm{o}\mathrm{r}\mathbf{x}\in {\mathbb{R}}^2\end{array}$

Then the fluid particles move with $d\mathbf{x}/dt=\mathbf{u}\left(\mathbf{x}\right)$ and the membrane moves with $d\mathbf{X}/dt=\mathbf{u}\left(\mathbf{X}\right)+\mathbf{v}\left(\mathbf{X}\right)$. The divergence of $\mathbf{v}$ may be defined by the density function $b\left(s\right)=\beta (\mathbf{f}·\widehat{\mathbf{n}})$, where $\widehat{\mathbf{n}}$ is the outward unit normal vector and $\beta$ is a coefficient that determines the membrane’s permeability, and is subject to the membrane properties.

Figure 1 shows a generic situation of the type considered in the next section. The figure shows the velocity field in a neighborhood of the membrane (top panel) as well as the divergence of the velocity (bottom panel) and the source doublet distribution that generates the flow. The bottom panel shows the vectors $-b\left(s\right)\widehat{\mathbf{n}}\left(s\right)$ (black arrows) and the divergence of the velocity (colormap). The latter is negative just inside the membrane, indicating a distribution of sinks, and positive just outside, indicating a distribution of sources.

The function ${\varphi }_{\delta }\left(\mathbf{x}\right)$ is used in the method of regularized Stokeslets [32,33] and is characterized by being smooth, radially symmetric functions concentrated at its center, like narrow Gaussians. The total integral is 1 and the parameter $\delta$ controls the spread of the function. In the formulation in (4) and (5), one could choose two different functions ${\varphi }_{\delta }\left(r\right)$ or use different values of $\delta$ even for the same function. In practice, the size of $\delta$ is influenced by the singularity of the corresponding singular kernel.

The solution of the problem in (4) is a regularized Stokeslet, given by [32,33]

(6)$\begin{array}{ccc}\hfill \mathbf{u}\left(\mathbf{x}\right)& =& \frac{1}{\mathsf{\mu }}{\int }_{\Gamma }(\mathbf{f}\left(s\right)·\nabla )\nabla B_{\delta }(\mathbf{x}-\mathbf{X}\left(s\right))-\mathbf{f}\left(s\right)G_{\delta }(\mathbf{x}-\mathbf{X}\left(s\right))ds\hfill \end{array}$

(7)$\begin{array}{ccc}\hfill p\left(\mathbf{x}\right)& =& {\int }_{\Gamma }\mathbf{f}\left(s\right)·\nabla G_{\delta }(\mathbf{x}-\mathbf{X}\left(s\right))ds\hfill \end{array}$

where the regularized Green’s and biharmonic functions, $G_{\delta }\left(r\right)$ and $B_{\delta }\left(r\right)$, corresponding to a radially symmetric blob ${\varphi }_{\delta }\left(r\right)$ are defined as the solutions of $\Delta G_{\delta }\left(r\right)={\varphi }_{\delta }\left(r\right)$ and $\Delta B_{\delta }\left(r\right)=G_{\delta }\left(r\right)$ in ${\mathbb{R}}^2$.

Derivation of the Source Doublet Solution

Take the divergence of the Stokes Equation (5) to get $\mathsf{\mu }\Delta (\nabla ·\mathbf{v})=\Delta q$, and substitute the divergence expression from Equation (5), from which we find that the pressure is given by

(8)$q\left(\mathbf{x}\right)=-\mathsf{\mu }\frac{1}{\mathsf{\mu }}{\int }_{\Gamma }b\left(s\right)\widehat{\mathbf{n}}\left(s\right)·\nabla {\varphi }_{\delta }(\mathbf{x}-\mathbf{X}\left(s\right))ds.$

Now the Stokes equation becomes

$\mathsf{\mu }\Delta \mathbf{v}=-\mathsf{\mu }{\int }_{\Gamma }\frac{b\left(s\right)}{\mathsf{\mu }}(\widehat{\mathbf{n}}\left(s\right)·\nabla )\nabla {\varphi }_{\delta }(\mathbf{x}-\mathbf{X}\left(s\right))ds$

which leads to the velocity

(9)$\mathbf{v}\left(\mathbf{x}\right)=-{\int }_{\Gamma }\frac{b\left(s\right)}{\mathsf{\mu }}(\widehat{\mathbf{n}}\left(s\right)·\nabla )\nabla G_{\delta }(\mathbf{x}-\mathbf{X}\left(s\right))ds$

For convenience we set $\widehat{\mathbf{x}}=\mathbf{x}-\mathbf{X}\left(s\right)$, $r=|\widehat{\mathbf{x}}|$ and summarize the solutions above as

$\begin{array}{ccc}\hfill \mathrm{Stokeslet}\text{:}\mathbf{u}\left(\mathbf{x}\right)& =& \frac{1}{\mathsf{\mu }}{\int }_{\Gamma }\mathbf{f}\left(s\right)H_1\left(r\right)+(\mathbf{f}\left(s\right)·\widehat{\mathbf{x}})\widehat{\mathbf{x}}{H}_2\left(r\right)ds\hfill \\ \hfill p\left(\mathbf{x}\right)& =& {\int }_{\Gamma }(\mathbf{f}\left(s\right)·\widehat{\mathbf{x}})P\left(r\right)ds\hfill \\ \hfill \mathrm{Source}\mathrm{Doublet}\text{:}\mathbf{v}\left(\mathbf{x}\right)& =& -\frac{1}{\mathsf{\mu }}{\int }_{\Gamma }b\left(s\right)[\widehat{\mathbf{n}}\left(s\right)S_1\left(r\right)+(\widehat{\mathbf{n}}\left(s\right)·\widehat{\mathbf{x}})\widehat{\mathbf{x}}{S}_2\left(r\right)]ds\hfill \\ \hfill q\left(\mathbf{x}\right)& =& -{\int }_{\Gamma }b\left(s\right)(\widehat{\mathbf{n}}\left(s\right)·\widehat{\mathbf{x}})Q\left(r\right)ds\hfill \end{array}$

where the regularizing functions are given by

(10)$\begin{array}{ccc}\hfill H_1\left(r\right)& =& \frac{{B_{\delta }}^{\prime }\left(r\right)}{r}-G_{\delta }\left(r\right)+\frac{1}{{8\pi }^{\prime }}{H}_2\left(r\right)=\frac{r{B_{\delta }}^{″}\left(r\right)-{B_{\delta }}^{\prime }\left(r\right)}{r^3}\hfill \\ \hfill S_1\left(r\right)& =& \frac{{G_{\delta }}^{\prime }\left(r\right)}{r},S_2\left(r\right)=\frac{r{G_{\delta }}^{″}\left(r\right)-{G_{\delta }}^{\prime }\left(r\right)}{r^3}\hfill \\ \hfill P\left(r\right)& =& \frac{{G_{\delta }}^{\prime }\left(r\right)}{r},Q\left(r\right)=\frac{{{\varphi }_{\delta }}^{\prime }\left(r\right)}{r}\hfill \end{array}$

2.3. Choice of Blobs

For a radially symmetric blob ${\varphi }_{\delta }\left(r\right)$ the minimum requirement is to have total integral equal to one

$1=2\pi {\int }_0^{\infty }r{\varphi }_{\delta }\left(r\right)dr.$

We use conditions derived in the first numerical example as part of the process of calibrating the source doublet coefficient $\beta$. In all cases, we choose a blob ${\varphi }_{\delta }\left(r\right)$ and we use Equation (10) to find all necessary regularizing functions. Table 1 summarizes the blob functions used in this article.

3. Numerical Examples

In order to assess the method, we present first two numerical examples whose exact solutions are known.

3.1. Example 1: A Circular Permeable Membrane under Tension

Consider a circular permeable membrane $\Gamma$ given initially by $\mathbf{X}\left(s\right)=(R_{0}cos\left(\frac{s}{R_0}\right),$ $R_{0}sin\left(\frac{s}{R_0}\right))$, where s is arclength. The force density is given by $\mathbf{f}\left(s\right)=\kappa \left(s\right)\widehat{\mathbf{n}}\left(s\right)$, where $\kappa$ is curvature. The equations of motion are

$\begin{array}{ccc}\hfill \mathsf{\mu }\Delta \mathbf{u}=\nabla p-\mathbf{F},\nabla ·\mathbf{u}=0& & \mathrm{i}\mathrm{n}{\mathbb{R}}^2/\Gamma \hfill \\ \hfill \mathbf{F}\left(\mathbf{x}\right)={\int }_{\Gamma }\mathbf{f}\left(s\right)\delta (\mathbf{x}-\mathbf{X}\left(s\right))ds& & \mathrm{f}\mathrm{o}\mathrm{r}\mathbf{x}\in {\mathbb{R}}^2\hfill \\ \hfill \mathbf{U}=-\frac{\alpha }{\mathsf{\mu }}〚p〛\widehat{\mathbf{n}}& & \mathrm{f}\mathrm{o}\mathrm{r}\mathbf{x}\in \Gamma \hfill \end{array}$

Here, $\alpha$ is a permeability parameter and $〚p〛$ represents the jump in pressure across the membrane, $〚p〛=p\left({\mathbf{X}}^{+}\right)-p\left({\mathbf{X}}^{-}\right)$. The parameters in this example were assumed to be dimensionless and given by $R_0=1$, $\alpha =0.1$ and $\mathsf{\mu }=1$.

Exact solution. In this simple example, the membrane will remain circular for all times as it shrinks to a point in finite time. Due to incompressibility, the tension in the membrane serves to increase the pressure inside the circle but produces no fluid motion, so $\mathbf{u}=0$. The added membrane velocity is

$\mathbf{U}=-\frac{\alpha }{\mathsf{\mu }}〚p〛\widehat{\mathbf{n}}$

where the pressure jump is computed from $〚p〛=\mathbf{f}\left(s\right)·\widehat{\mathbf{n}}\left(s\right)=\kappa$. Therefore $\mathbf{U}=-\frac{\alpha }{\mathsf{\mu }}\kappa \widehat{\mathbf{n}}$. For a circular membrane, this reduces to an equation for its radius

(11)$\frac{dR}{dt}=-\frac{\alpha }{\mathsf{\mu }R},R\left(0\right)=R_0$

with solution

$R_{exact}\left(t\right)=\sqrt{R_0^2-\frac{2\alpha }{\mathsf{\mu }}t}\mathrm{f}\mathrm{o}\mathrm{r}0\le t\le \frac{R_0^2\mathsf{\mu }}{2\alpha }.$

Calibrating $\mathbf{\beta }$. The membrane velocity $\mathbf{v}\left(s\right)$ is normal to the membrane and $\mathbf{v}\left(s\right)·\widehat{\mathbf{n}}\left(s\right)$ is independent of position s along the membrane due to the circular symmetry. Thus, we compute it at $s=0$, corresponding to $\mathbf{X}\left(0\right)=\left(r\right(t),0)$. At this location, Equation (9) and $b\left(s\right)=\beta \kappa$ give

$\begin{array}{ccc}\hfill \frac{dR}{dt}=\mathbf{v}·\widehat{\mathbf{n}}& =& -\frac{1}{\mathsf{\mu }}{\int }_{\Gamma }b\left(s\right)(\widehat{\mathbf{n}}\left(s\right)·\nabla )\frac{\partial }{\partial x}{G}_{\delta }(\mathbf{x}-\mathbf{X}\left(s\right))ds\hfill \\ & =& -\frac{\beta }{\mathsf{\mu }}\kappa {\int }_{\Gamma }(\widehat{\mathbf{n}}\left(s\right)·\nabla )\frac{\partial }{\partial x}{G}_{\delta }(\mathbf{x}-\mathbf{X}\left(s\right))ds\hfill \\ & =& -\frac{\beta }{\mathsf{\mu }}\kappa R{\int }_0^{2\pi }(\widehat{\mathbf{n}}\left(\theta \right)·\nabla )\frac{\partial }{\partial x}{G}_{\delta }(\mathbf{x}-\mathbf{X}\left(\theta \right))d\theta \hfill \end{array}$

Using the particular blob ${\varphi }_{\delta }\left(r\right)=\frac{2{\delta }^4}{\pi {(r^2+{\delta }^2)}^3}$ we find that

(12)$\begin{array}{ccc}\hfill \mathbf{v}·\widehat{\mathbf{n}}& =& -\frac{\beta }{\mathsf{\mu }}\kappa R{\int }_0^{2\pi }\frac{3R^2(R^2+{\delta }^2)+({\delta }^4-3{\delta }^2{R}^2-4R^4)cos\left(\theta \right)+R^{4}cos\left(2\theta \right)}{\pi {({\delta }^2+2R^2-2R^{2}cos\left(\theta \right))}^3}d\theta \hfill \\ & =& -\frac{\beta }{\mathsf{\mu }}\kappa \frac{12R^3(2R^2+{\delta }^2)}{\delta {(4R^2+{\delta }^2)}^{5/2}}=-\frac{\beta \kappa }{\delta \mathsf{\mu }}\left(\frac{3}{4}-\frac{3{\delta }^2}{32R^2}-\frac{15{\delta }^4}{512R^4}+O({\delta }^6/R^6)\right).\hfill \end{array}$

In order for this velocity to be consistent with the exact velocity $\mathbf{v}·\widehat{\mathbf{n}}=-\alpha \kappa /\mathsf{\mu }$ in Equation (11), we choose

$\beta =\frac{4}{3}\delta \alpha .$

The coefficient $4/3$ depends on the blob used.

Based on the calculation of $\mathbf{v}·\widehat{\mathbf{n}}$, we set $\beta =(4/3)\delta \alpha$ and computed numerically the solution of the single differential equation $R^{\prime }\left(t\right)=-\beta \frac{12R^2(2R^2+{\delta }^2)}{\mathsf{\mu }\delta {(4R^2+{\delta }^2)}^{5/2}}$ with $R\left(0\right)=1$, corresponding to the blob ${\varphi }_{\delta }\left(r\right)$. The solution using $\delta =0.1$ and $\mathsf{\mu }=1$ is compared with the exact solution in Figure 2.

Notice that the two curves in the figure are nearly indistinguishable until the circle has shrunk to a radius of about 0.1. Beyond this point, the blob size $\delta$ is equal to or larger than the radius of the circle and we do not expect the regularized velocity to approximate well the exact velocity. To ameliorate this issue we also computed the solution by scaling $\delta$ dynamically with the circle radius and reduced the error from a maximum of 0.078 to a maximum of 0.016.

We also computed the solution using a full simulation based on the velocity field

$\mathbf{v}\left(\mathbf{x}\right)=-\frac{1}{\mathsf{\mu }}{\int }_{\Gamma }b\left(s\right)[\widehat{\mathbf{n}}\left(s\right)S_1\left(r\right)+(\widehat{\mathbf{n}}\left(s\right)·\widehat{\mathbf{x}})\widehat{\mathbf{x}}{S}_2\left(r\right)]ds$

with $b\left(s\right)=\beta \kappa \left(s\right)$. We discretized the initial circle using $N=400$ points equally spaced. The first and second derivatives of the membrane discretization with respect to the arc length parameter $s\in [0,2\pi R_0)$, are approximated by

$D{\mathbf{X}}_n=\frac{({\mathbf{X}}_{n+1}-{\mathbf{X}}_{n-1})}{2\Delta s_0},D^2{\mathbf{X}}_n=\frac{{\mathbf{X}}_{n+1}-2{\mathbf{X}}_n+{\mathbf{X}}_{n-1}}{\Delta s_0^2}$

where $\Delta s_0=2\pi R_0/N$. They are used to approximate the outer normal vectors and curvature

${\mathbf{N}}_n=\frac{(D{Y}_n,-D{X}_n)}{\left|D\mathbf{X}\right|},K_n=\frac{D{X}_n{D}^2{Y}_n-D{Y}_n{D}^2{X}_n}{{\left|D\mathbf{X}\right|}^3}.$

For $t>0$, we set $\Delta s\left(t\right)=2\pi R\left(t\right)/N$, where $R\left(t\right)=\frac{1}{N}{\sum }_n\left|{\mathbf{X}}_n\right|$. Using the notation $\widehat{\mathbf{x}}={\mathbf{X}}_j-{\mathbf{X}}_n$ and $r_{jn}=|{\mathbf{X}}_j-{\mathbf{X}}_n|$, the membrane points are evolved with the equation

(13)$\frac{d{\mathbf{X}}_j}{dt}=-\frac{1}{\mathsf{\mu }}\sum _{n=1}^N\beta K_n\left({\mathbf{N}}_n{S}_1\left(r_{jn}\right)+({\mathbf{N}}_n·\widehat{\mathbf{x}})\widehat{\mathbf{x}}{S}_2\left(r_{jn}\right)\right)\Delta s\left(t\right),j=1,2,\dots N$

The calibration method for $\beta$ allows for the introduction of two parameters $a_1$ and $a_2$ as unknown coefficients in the definition of the blob. We define a blob with integral equal to 1

$\psi =\frac{2{\delta }^4((6-a_1-3a_0)r^4+a_1{r}^2{\delta }^2+a_0{\delta }^4)}{\pi {(r^2+{\delta }^2)}^5}$

Following the same calibration procedure as before leads to a version of Equation (12) of the form

$\mathbf{v}·\widehat{\mathbf{n}}=-\frac{\beta \kappa }{\mathsf{\mu }\delta }\left(\frac{9+13a_0+a_1}{32}+\frac{3(-15+5a_0+a_1)}{256}\frac{{\delta }^2}{R^2}+\frac{15(-35+17a_0+5a_1)}{4096}\frac{{\delta }^4}{R^4}+O\left(\frac{{\delta }^6}{R^6}\right)\right)$

and we select the parameters $a_1=-10$ and $a_0=5$ to make $-15+5a_0+a_1=-35+17a_0+5a_1=0$ and reduce the expression to

(14)$\mathbf{v}·\widehat{\mathbf{n}}=-\frac{\beta \kappa }{\mathsf{\mu }\delta }\left(2+O\left(\frac{{\delta }^6}{R^6}\right)\right)$

which leads to $\beta =\frac{1}{2}\delta \alpha$ and the blob

${\psi }_{\delta }\left(r\right)=\frac{2{\delta }^4(r^4-10{\delta }^2{r}^2+5{\delta }^4)}{\pi {(r^2+{\delta }^2)}^5}.$

We mention that Beale and collaborators [34,35] have derived conditions for regularized single and double layer potentials that increase the order of the regularization error. The conditions are that even moments of a certain function must be zero. Specifically, for each $k=0,1,2,\dots$, the condition

(15)${\int }_0^{\infty }(F\left(r\right)-1)r^{2k}dr=0,\text{where }F\left(r\right)=2\pi {\int }_0^{r}s{\varphi }_{\delta }\left(s\right)ds$

increases the order of accuracy by two. Using these conditions (for $k=0$ and $k=1$) leads to the same two equations for $a_0$ and $a_1$, so that the two procedures are equivalent.

The results are shown in Figure 3 using both blobs ${\varphi }_{\delta }\left(r\right)=\frac{2{\delta }^4}{\pi {(r^2+{\delta }^2)}^3}$ and the more accurate ${\psi }_{\delta }$ reducing the maximum error from 0.0827 to 0.0216 with fixed $\delta =0.1$.

3.2. Example 2: A Circular Permeable Membrane with Circular Equilibrium

We consider the same circular membrane as in the previous example except that the force density along the membrane is given by $\mathbf{f}\left(s\right)=(\frac{R_{eq}}{R\left(t\right)}-1)\widehat{\mathbf{n}}\left(s\right)$, where $R_{eq}$ is the equilibrium radius assumed to satisfy $0\le R_{eq}<R_0$. In this case, the circle will shrink from initial radius $R_0$ to an asymptotic value of $R_{eq}$. The exact solution of the problem is given by [21]

$R_{exact}\left(t\right)=R_{eq}\left(1+W\left(c{e}^{-\alpha t/\mathsf{\mu }{R}_{eq}}\right)\right),\text{where }c=\left(\frac{R_0}{R_{eq}}-1\right)exp\left(\frac{R_0}{R_{eq}}-1\right)$

and $W\left(x\right)$ is the Lambert W-function. Using the parameter values $\alpha =0.1$, $\mathsf{\mu }=1$, $R_0=1$ and $R_{eq}=0.75$, we discretized the initial circle with N points ${\mathbf{X}}_n$, as before. The normal vectors, particle spacing $\Delta s\left(t\right)$ and circle radius $R\left(t\right)$ were computed as in the previous example so that the membrane motion is given by

(16)$\frac{d{\mathbf{X}}_j}{dt}=-\sum _{n=1}^N\beta \left(\frac{R_{eq}}{R\left(t\right)}-1\right)\left({\mathbf{N}}_n{S}_1\left(r_{jn}\right)+({\mathbf{N}}_n·\widehat{\mathbf{x}})\widehat{\mathbf{x}}{S}_2\left(r_{jn}\right)\right)\Delta s\left(t\right),j=1,2,\dots N$

where $\widehat{\mathbf{x}}={\mathbf{X}}_j-{\mathbf{X}}_n$ and $r_{jn}=|{\mathbf{X}}_j-{\mathbf{X}}_n|$. This system of equations is solved using a 4th order Runge–Kutta method and the blob ${\varphi }_{\delta }\left(r\right)$ in Table 1. The computed radius was compared to the exact solution resulting in the maximum errors shown in Table 2. The initial arc length parameter is computed as $\Delta s_0=2\pi R_0/N$ where $R_0=1$ and N is the number of membrane points. There are several observations we can make about the errors. Notice that the value of the regularization parameter $\delta$ is not constant in each column of the table since it depends on N. The errors for constant $\delta$ are along diagonals of the table; for example, the underlined values correspond to $\delta =\pi /50$. With $\delta$ held constant, the errors decrease as $\Delta s_0$ decreases until the discretization error is small enough that the dominant error is due to the regularization $\delta$ (lower right corner of the table). If $\delta$ is chosen proportional to $\Delta s_0$ with sufficiently large proportionality constant (rightmost column of the table), the errors decrease by a factor of 4 as N doubles, indicating second-order convergence as expected from the choice in Equation (12). If $\delta$ is too small compared to $\Delta s_0$, the errors reach a plateau (leftmost columns). The corresponding results using the blob $\psi \left(r\right)$ in Table 1 are displayed in Table 3. Notice that the errors are smaller by several orders of magnitude.

Each row corresponds to a fixed discretization of the membrane and we see that for fixed N, as $\delta$ increases from $0.5\Delta s_0$ to $8\Delta s_0$, the errors decrease to a minimum and eventually increase. The minimum values are smaller for finer discretizations (i.e., for larger N) and they occur at values of $\delta$ that are not the same multiple of $\Delta s_0$. In other words, the errors decrease as both $\delta$ and $\Delta s_0$ are decreased but the optimal value of $\delta$ is not a fixed multiple of $\Delta s_0$. This can be observed in Figure 4 left, which shows a log-log plot of the error versus $\delta$. When $\delta$ is too small compared to $\Delta s_0$, the errors are large since the blobs do not overlap sufficiently and fluid may leak between discretization nodes regardless of the value of $\beta$. On the other hand, when $\delta$ is too large compared to $\Delta s_0$, the regularization error is the largest contributor and the errors collapse onto a single line (in log-log scale) as long as the discretization is fine enough to approximate the smooth integral. Notice that the slope of the line corresponds to the sixth order leading error term in the value of $\beta$ found from the asymptotic expansion in Equation (14).

Figure 4 right shows the same errors except that they are plotted as a function of $(\delta /\Delta s)$. In this scaling, the errors on the left side of the figure collapse onto a curve $Error\propto exp(-6\delta /\Delta s)$, indicating that in this example the discretization error decreases exponentially for the values tested when $\delta <3\Delta s$.

3.3. Example 3: Computing the Pressure

We consider an initially elliptical permeable membrane $\Gamma$ given by $\mathbf{X}=(cos\theta ,1/2sin\theta )$. As before, the force density is given by $\mathbf{f}\left(s\right)=\kappa \left(s\right)\widehat{\mathbf{n}}\left(s\right)$ and the equations of motion are

$\begin{array}{ccc}\hfill \mathsf{\mu }\Delta \mathbf{u}=\nabla p-\mathbf{F},\nabla ·\mathbf{u}=0& & \mathrm{i}\mathrm{n}{\mathbb{R}}^2/\Gamma \hfill \\ \hfill \mathbf{F}\left(\mathbf{x}\right)={\int }_{\Gamma }\mathbf{f}\left(s\right)\delta (\mathbf{x}-\mathbf{X}\left(s\right))ds& & \mathrm{f}\mathrm{o}\mathrm{r}\mathbf{x}\in {\mathbb{R}}^2\hfill \\ \hfill \mathbf{U}=-\frac{\alpha }{\mathsf{\mu }}〚p〛\widehat{\mathbf{n}}& & \mathrm{f}\mathrm{o}\mathrm{r}\mathbf{x}\in \Gamma \hfill \end{array}$

The parameters used were $\alpha =0.01$, $N=400$, and $\mathsf{\mu }=1$. The exact solution is unknown. The numerical solution at selected times is shown in Figure 5. Since the force density is proportional to the local curvature, the points along the major axis, where the curvature is largest, move in faster than other parts of the membrane. The two points along the minor axis move outward, changing the membrane’s shape into a more circular one. If the permeability is zero, no fluid can escape through the membrane and the membrane approaches a circle with area equal to the area bounded by the initial condition. When the permeability is nonzero, $\alpha >0$, fluid crosses the membrane as a function of the pressure drop across the boundary. This leads to an ever shrinking shape until it approaches a single point in finite time.

The computation proceeds as follows: for $\widehat{\mathbf{x}}=\mathbf{x}-\mathbf{X}\left(s\right)$ and $r=|\widehat{\mathbf{x}}|$, we compute the fluid velocity

$\mathbf{u}\left(\mathbf{x}\right)={\int }_{\Gamma }\mathbf{f}\left(s\right)H_1\left(r\right)+(\mathbf{f}\left(s\right)·\widehat{\mathbf{x}})\widehat{\mathbf{x}}{H}_2\left(r\right)ds,\mathbf{f}\left(s\right)=\kappa \left(s\right)\widehat{\mathbf{n}}\left(s\right)$

using the discretization

(17)$\mathbf{u}\left(\mathbf{x}\right)=\sum _{n=1}^N{K}_n\left[{\mathbf{N}}_n{H}_1\left(r\right)+({\mathbf{N}}_n·(\mathbf{x}-{\mathbf{X}}_n))(\mathbf{x}-{\mathbf{X}}_n)H_2\left(r\right)\right]\Delta s\left(t\right),$

with $r=|\mathbf{x}-{\mathbf{X}}_n|$. The membrane velocity at ${\mathbf{X}}_j$ is given by

(18)$\frac{d{\mathbf{X}}_j}{dt}=\mathbf{u}\left({\mathbf{X}}_j\right)-\sum _{n=1}^N\beta K_n\left({\mathbf{N}}_n{S}_1\left(r_{jn}\right)+({\mathbf{N}}_n·\widehat{\mathbf{x}})\widehat{\mathbf{x}}{S}_2\left(r_{jn}\right)\right)\Delta s\left(t\right),$

where $\widehat{\mathbf{x}}={\mathbf{X}}_j-{\mathbf{X}}_n$ and $r_{jn}=|{\mathbf{X}}_j-{\mathbf{X}}_n|$. In this example, the formulas are derived for the blob ${\psi }_{\delta }\left(r\right)$ from Table 1, where we used the value $\beta =\frac{1}{2}\delta \alpha$ based on the calibration with this blob.

The pressure is given by

$\begin{array}{ccc}\hfill p\left(\mathbf{x}\right)& =& {\int }_{\Gamma }\kappa \left(s\right)(\widehat{\mathbf{n}}\left(s\right)·\nabla )G_{\delta }\left(\right|\mathbf{x}-\mathbf{X}\left(s\right))ds\hfill \\ & =& {\int }_{\Gamma }(\kappa \left(s\right)-\kappa \left(s_0\right))(\widehat{\mathbf{n}}\left(s\right)·\nabla )G_{\delta }\left(\right|\mathbf{x}-\mathbf{X}\left(s\right))ds+\kappa \left(s_0\right)\chi \left(\mathbf{x}\right)\hfill \end{array}$

where $\chi \left(\mathbf{x}\right)=-1$ when $\mathbf{x}$ is inside the ellipse and 0 otherwise. A surface plot of the pressure, along with its cross-section in the $xz$-plane and contours in the $xy$-plane are shown in Figure 6.

3.4. Example 4: Flow in a Channel with Permeable Walls

We consider flow in a channel with walls that are permeable or have permeable sections, as depicted schematically in Figure 7. We assume that the flow at the inlet and at impermeable sections of the channel walls is given. The goal is to determine the flow within the channel for given boundary permeability. If the fluid velocity on the boundary of the channel were known, the flow in the channel could be computed using Stokeslets by solving for the boundary force. In the current case however, the fluid flow along the permeable boundary is unknown and satisfies $\mathbf{u}·\widehat{\mathbf{n}}=-\frac{\alpha }{\mathsf{\mu }}〚p〛$, as in previous examples.

As a first step we use Stokeslets and source doublets to compute the velocity boundary condition in the permeable sections of the channel walls. Let $\Gamma$ represent the inlet, bottom and top sides of the channel, which constitute the boundary of the computational domain where forces are applied. Note that the channel outlet is not part of $\Gamma$. Let ${\Gamma }_{\beta }$ be the permeable portion of $\Gamma$, where $\beta \ne 0$, and set $\mathbf{X}\left(s\right)$ to be a parametrization of $\Gamma$. We look for a force distribution $\mathbf{g}\left(s\right)$ on the channel walls along with a distribution of source doublets $\beta \left(s\right)(\mathbf{g}\left(s\right)·\widehat{\mathbf{n}}\left(s\right))$ in the permeable region, with $\beta \left(s\right)$ given, in order to determine the velocity on ${\Gamma }_{\beta }$.

The procedure starts by setting ${\mathbf{u}}_b$ on $\Gamma$ as

${\mathbf{u}}_b\left(s\right)=\left\{\begin{array}{cc}{\mathbf{U}}_0\left(s\right),\hfill & \mathbf{X}\left(s\right)\in \Gamma \setminus {\Gamma }_{\beta }\hfill \\ \text{unknown},\hfill & \mathbf{X}\left(s\right)\in {\Gamma }_{\beta }\hfill \end{array}\right.$

where ${\mathbf{U}}_0\left(s\right)$ is known at the inlet (parabolic) and along impermeable boundaries (zero). Setting $\widehat{\mathbf{x}}=\mathbf{X}\left(s_0\right)-\mathbf{X}\left(s\right)$ and $r=|\widehat{\mathbf{x}}|$ and using the following notation for the regularized Stokeslet and source doublet kernels

$\begin{array}{ccc}\hfill St(\widehat{\mathbf{x}},s,s_0)& =& \mathbf{g}\left(s\right)H_1\left(r\right)+(\mathbf{g}\left(s\right)·\widehat{\mathbf{x}})\widehat{\mathbf{x}}{H}_2\left(r\right)\hfill \\ \hfill SD(\widehat{\mathbf{x}},s,s_0)& =& -\beta \left(s\right)(\mathbf{g}\left(s\right)·\widehat{\mathbf{n}}\left(s\right))\left[\widehat{\mathbf{n}}\left(s\right)S_1\left(r\right)+(\widehat{\mathbf{n}}\left(s\right)·\widehat{\mathbf{x}})\widehat{\mathbf{x}}{S}_2\left(r\right)\right]\hfill \end{array}$

we enforce the boundary conditions

(19)$\left\{\begin{array}{cc}{\int }_{\Gamma }St(\widehat{\mathbf{x}},s,s_0)+SD(\widehat{\mathbf{x}},s,s_0)ds=\mathsf{\mu }{\mathbf{U}}_0\left(s_0\right),\hfill & \mathbf{X}\left(s_0\right)\in \Gamma \setminus {\Gamma }_{\beta }\hfill \\ {\int }_{\Gamma }St(\widehat{\mathbf{x}},s,s_0)ds=0,\hfill & \mathbf{X}\left(s_0\right)\in {\Gamma }_{\beta }\hfill \end{array}\right.$

These are based on the assertion that the known boundary velocity is the sum of contributions from the Stokeslets and source doublets, while the unknown velocity across permeable channel walls is entirely due to the source doublets. The Stokeslets velocity contribution in those sections is zero.

Equation (19) must be inverted to find $\mathbf{g}\left(s\right)$. Once the boundary force is known, the fluid velocity at the permeable boundary ${\Gamma }_{\beta }$ is set to ${\int }_{\Gamma }SD(\widehat{\mathbf{x}},s,s_0)ds$. Now that all boundary velocities ${\mathbf{u}}_b$ are known, the second step is to invert the Stokeslet kernel

${\int }_{\Gamma }{H}_1\left(r\right)+(\mathbf{f}\left(s\right)·\widehat{\mathbf{x}})\widehat{\mathbf{f}}{H}_2\left(r\right)ds=\mathsf{\mu }{\mathbf{u}}_b$

to find $\mathbf{f}\left(s\right)$ and use it to compute the flow in the channel.

Figure 8 shows the result of the velocity boundary condition computation as well as the computed velocity near the outlet of the channel (right side) using the blob ${\psi }_{\delta }\left(r\right)$ in Table 1. The figures are for three different values of $\beta$ while the inlet flow is the same in all cases. The viscosity was set to $\mathsf{\mu }=1$ and the numerical parameters used were $\delta =0.0224$ and wall discretization size $h=1/160$. The figure shows that as $\left|\beta \right|$ increases, the permeability of the wall segment also increases, leading to more fluid escaping the channel at the top wall. When $\beta =-0.00024$, more fluid escapes through the permeable wall than comes in through the inlet, leading to flow reversal at the right end of the channel. Figure 8 shows the velocity field in the channel computed from the Stokeslets $\mathbf{f}\left(s\right)$ that result from imposing the velocity boundary conditions. For $\beta =-0.0001$, particles entering the channel at the inlet will flow to the outlet except those that are near the top wall. These get ejected out of the channel through the permeable top wall segment. The same is true for $\beta =-0.00018$ although particles in a larger region of the channel will flow out of the top. For $\beta =-0.00024$, particles entering the channel at the inlet will necessarily flow out of the top permeable wall segment due to the flow reversal on the right side of the channel.

The flow rates shown in Table 4 were computed at the inlet, permeable wall, and outlet by a midpoint rule approximation of the integral

$\mathrm{F}\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}={\int }_{\Gamma }\mathbf{u}\left(s\right)·\widehat{\mathbf{n}}\left(s\right)ds$

Note that the exact flow rate at the inlet is ${\int }_0^{1}4y(1-y)dy=-2/3$.

3.5. Example 5: Flow in a Capillary with Protein Concentration

We consider a biological application of our formulation by modeling a single rat glomerular capillary using source doublets to represent pores in the membrane walls. The renal glomerulus is the structure responsible for filtering the blood in the kidney, and is composed of a tortuous network of over 300 capillary segments, wrapped into a sphere and surrounded by an outer capsule filled with the fluid filtrate from the capillaries [36]. Since this space is made up of fluid and not tissue, it is reasonable to model a glomerular capillary as a channel filled with and surrounded by fluid. Filtration in the glomerular capillaries is driven by a hydrostatic pressure $P_{hyd}$ [37]. The glomerular capillary walls are selectively permeable to fluid and solutes such as sodium and potassium, but impermeable to larger macromolecules including plasma proteins such as albumin. Due to the plasma protein concentration differential across the capillary wall between the blood plasma and filtrate, a colloid osmotic pressure, $\mathsf{\Pi }$, is exerted in the direction of highest concentration, in this case back into the capillary and thus in opposition to $P_{hyd}$. An expression for $\mathsf{\Pi }$ has been derived experimentally as a function of plasma protein concentration $C_p$:

$\mathsf{\Pi }=a_1{C}_p+a_2{C}_p^2$

where $C_p$ is measured in g/dL and $a_1=1.629$ mmHg g${}^{-1}$ dL, $a_2=0.2935$ mmHg${}^2$ g${}^{-2}$ dL${}^2$ [37]. Many 1-D models have been developed to estimate fluid filtration along the length of a glomerular capillary or the glomerulus as a whole [38,39,40,41,42]. Some computational fluid dynamics models have been developed to investigate glomerular filtration, including a model which used the immersed boundary method and finite element analysis to consider solute filtration in an idealized glomerulus in 3-D [43]. This study, while robust and detailed, did not consider the anatomy of the glomerulus and instead condensed the glomerulus into five cylinders that filtered in parallel. A 2-D model of fluid and solute filtration across a single portion of the glomerular membrane [44] studied intensely the sieving of solutes through the membrane, but did not consider the fluid dynamics within the capillary. To our knowledge the application of computational fluid dynamics to consider the dynamics of fluid filtration in an anatomically-accurate glomerular capillary in 2-D has not been performed.

To accomplish this task, we modify the previously-described channel algorithm (example 4) by introducing a concentration of plasma proteins $C_p$ that are transported by the fluid within the channel domain and cannot pass through the membrane:

$C_p(x,y)=0\mathrm{for}(x,y)\notin \{0<x<L,0<y<H\}$

$\frac{d{C}_p}{dt}+\mathbf{u}·\nabla C_p=\sigma \Delta C_p,$

with homogeneous Neumann conditions on the boundaries of the domain. At the inlet and outlet, the concentration is advected by the flow, so there is a need for an assumption on the concentration upstream of the inlet. We assume that $C_p(0^{-},y)$ is fixed at a high value so that plasma proteins are continuously introduced into the channel at the inlet.

To obtain the fluid velocity $\mathbf{u}$ used to advect the concentration, we follow the algorithm outlined in example 4 with some modifications. Namely, the initial boundary flow ${\mathbf{u}}_b$ must be modified to account for the influence of ${\mathsf{\Pi }}_g$ driving fluid back into the capillary and $P_{hyd}$ driving fluid out. ${\mathsf{\Pi }}_g$ is computed based on the concentration nearest to the permeable wall:

$\begin{array}{ccc}\hfill \mathsf{\Pi }(x,0)& =& a_1{C}_p(x,0)+a_2{C}_p{(x,0)}^2\hfill \\ \hfill \mathsf{\Pi }(x,H)& =& a_1{C}_p(x,H)+a_2{C}_p{(x,H)}^2\hfill \end{array}$

Assuming a constant $P_{hyd}$ along the length of the capillary [44] the effective filtration pressure difference $P_{filt}$ across the channel walls is:

$\begin{array}{ccc}\hfill P_{filt}(x,0)& =& \mathsf{\Pi }(x,0)-P_{hyd}\hfill \\ \hfill P_{filt}(x,H)& =& -(\mathsf{\Pi }(x,H)-P_{hyd})\hfill \end{array}$

The resultant forces on the channel walls $\mathbf{h}$ are solved for by inverting the source doublet kernel for the pressure:

$P_{filt}\left(\mathbf{x}\right)=-{\int }_{\Gamma }(\mathbf{h}·\widehat{\mathbf{n}})(\widehat{\mathbf{n}}·\widehat{\mathbf{x}})Q\left(r\right)ds$

Unlike previously described, this step assumes the constant $\beta =1$ as part of the formulation for pressure. This is because, in this case, the pressure-to-force relationship is independent of the membrane permeability. The constant $\beta$ is incorporated in the next step using the source doublet notation:

$\begin{array}{ccc}\hfill S{D}_h(\widehat{\mathbf{x}},s,s_0)& =& -\beta \left(s\right)(\mathbf{h}\left(s\right)·\widehat{\mathbf{n}}\left(s\right))\left[\widehat{\mathbf{n}}\left(s\right)S_1\left(r\right)+(\widehat{\mathbf{n}}\left(s\right)·\widehat{\mathbf{x}})\widehat{\mathbf{x}}{S}_2\left(r\right)\right]\hfill \end{array}$

and we enforce the boundary conditions

(20)$\left\{\begin{array}{cc}{\int }_{\Gamma }St(\widehat{\mathbf{x}},s,s_0)+SD(\widehat{\mathbf{x}},s,s_0)ds=\mathsf{\mu }{\mathbf{U}}_0\left(s_0\right)-{\int }_{\Gamma }S{D}_h(\widehat{\mathbf{x}},s,s_0)ds,\hfill & \mathbf{X}\left(s_0\right)\in \Gamma \setminus {\Gamma }_{\beta }\hfill \\ {\int }_{\Gamma }St(\widehat{\mathbf{x}},s,s_0)ds=0.\hfill & \mathbf{X}\left(s_0\right)\in {\Gamma }_{\beta }\hfill \end{array}\right.$

The additional term on the right hand side of the equation is known, as such Equation (20) can be inverted to obtain $\mathbf{g}\left(s\right)$ and the boundary velocities in the permeable region ${\Gamma }_{\beta }$ are set to ${\int }_{\Gamma }SD(\widehat{\mathbf{x}},s,s_0)ds$. The boundary velocities ${\mathbf{u}}_b$ are known and the Stokeslet kernel is subsequently inverted for the forces $\mathbf{f}$ which are used to find the channel fluid velocity. This velocity is then used to advect the plasma protein concentration. In Figure 9, we perform a temporal simulation in which an initial Gaussian concentration profile of protein (yellow circles in top left panel) is advected by the fluid while the fluid is filtered at the wall. As additional concentration is introduced at the inlet (yellow at the left wall in the bottom panels), the concentration develops local maxima at different locations in the channel. This simulation demonstrates the versatility of our method for predicting distribution of proteins throughout the vessel, as the protein concentration can vary with both time and location within the vessel. This is a significant improvement on previous models, which only consider axial changes in plasma protein concentration [38,39,40,41,42].

The velocity at the inlet is defined as

(21)${\mathbf{U}}_0=u_{0}y(H-y)$

where $u_0$ is chosen based on experimental data (Table 5). Parameters are chosen for the model based on measured dimensions of rat glomerular capillaries, red blood cell velocities and plasma viscosity [45].

As explained in example 3, we use the blob ${\psi }_{\delta }\left(r\right)$, $\beta =\frac{1}{2}\frac{\delta \alpha }{\mathsf{\mu }}$. In this example, by defining $\alpha =k\mathsf{\mu }$ we use the same equation for $\beta$ and the dimensions of $\beta$ remain consistent. We calculate the fraction of filtered fluid to fluid entering the left end of the channel ($FF$) by integrating the normal velocity through the permeable section of the membrane and dividing by the left velocity boundary condition:

(22)$FF=\frac{{\int }_{{\Gamma }_{\beta }}\mathbf{u}\left(s\right)·\widehat{\mathbf{n}}\left(s\right)ds}{{\mathbf{U}}_0\left(s_0\right)}$

Additional analysis elucidates a nonlinear relationship between the hydraulic conductivity k and the filtered fluid volume (Figure 10). This nonlinearity is transferred to the forces exerted on the wall in the normal and tangential directions (Figure 11). For forces in the tangential direction, a loss of fluid due to permeability reduces the force magnitude on the length of the vessel. The difference between these tangential stress profiles is nonlinear, corresponding to the nonlinear relationship between k and $FF$. For forces in the normal direction, a linear drop in force across the vessel is seen for all values of k. This is expected due to the necessary pressure drop across the vessel to facilitate the flow imposed at the left boundary. However with increasing k the force in the normal direction is reduced nonlinearly, corresponding to the nonlinear relationship between k and $FF$.

Inherent to our current application of Stokeslets and source doublets to compute the fluid dynamics within a permeable glomerular capillary is the assumption that the walls are not deflected by flow through the capillary such that the capillary diameter remains constant in time. In reality, the glomerular capillary walls change in time due to the pulsatility of blood pressure. We have previously estimated the magnitude of strain of the glomerular capillary walls using a mathematical model of blood flow through an entire anatomically accurate glomerular capillary network [48], and found that the magnitude of glomerular capillary wall strain is, under physiological conditions, less than 1%. Future models will investigate dynamics of adding flexibility to the vessel wall as well as changing the diameter of the vessel as a function of length.

4. Discussion and Conclusions

We have presented a model for computing flows bounded by thin permeable membranes in which the amount of fluid that is filtrated across the membrane depends on local flow properties such as the pressure drop across the membrane, the concentration of a solute, and the membrane permeability.

The derivation was based on the concept of placing fluid sources on one side of the membrane and corresponding sinks on the other side, creating a source doublet field along the permeable parts of the membrane. We note that the resulting source doublet kernel is the same as that of a flow proportional to the pressure gradient due to a force normal to the boundary. Indeed, decomposing a force $\mathbf{f}{\varphi }_{\delta }$ into a gradient $\nabla q$ and a zero-divergence vector $\mathbf{u}$ field gives

$\mathbf{f}{\varphi }_{\delta }=\nabla q+\mathbf{u}\Rightarrow (\mathbf{f}·\nabla )\nabla G_{\delta }=\nabla q.$

The kernel in Equation (9) is the same as this expression for $\nabla q$.

Our model builds on the method of regularized Stokeslets through the use of regularized sources. The resulting boundary integral formulation is free of singularities and the regularization parameter $\delta$ provides a length scale for the thickness of the membrane. In addition, the strength of the source doublet is proportional to the regularization parameter with proportionality constant that depends on the ${\varphi }_{\delta }$ used. Our numerical results demonstrate that designing the regularization function ${\varphi }_{\delta }$ to satisfy the moment condition in Equation (15) [34,35] leads to higher accuracy, in practice.

The last two examples present applications to channel flow. Example 4 involves a channel with only a section of the boundary that is permeable to the fluid, while the velocity field at the Inlet is assumed to be given and the no slip condition is enforced in the impermeable portions of the boundary, the velocity boundary condition is unknown in the permeable section, and therefore, unknown also at the outlet. We use the model proposed here specifically to determine the velocity boundary condition at the permeable section of the channel in a way that is consistent with the normal stress along the boundary. Once this is known, the flow in the channel is computed using the method of regularized Stokeslets. The example shows that if the permeability is high enough, more fluid escapes through the permeable membrane than fluid comes in through the inlet, forcing fluid into the channel at the outlet. The final example provides an application in which a concentration $C_p$ flows and diffuses in the channel, whose walls are permeable to the fluid but not to $C_p$, and the concentration level affects the permeability locally.

There are extensions to this model that we expect to pursue in the future. Each wall of the channel can be treated as a set of connected linear segments with a force density that varies linearly on each segment. Then the integral in Equations (6) and (9) can be computed analytically, as explained in [49]. This would allow the use of substantially smaller values of the regularization $\delta$ without increasing (in fact, decreasing) the number of nodes discretizing the boundary. Finally, we note that our boundary integral model describing the flow across a permeable membrane is valid in three dimensions as well, so that extending the method to 3D is straight forward.

Author Contributions

Conceptualization, R.C. and M.H.-V.; methodology, R.C.; software, R.C., M.H.-V. and O.R.; formal analysis, R.C., M.H.-V. and O.R.; writing—original draft preparation, R.C. and M.H.-V.; writing—review and editing, R.C. and O.R.; visualization, R.C., M.H.-V. and O.R.; supervision, R.C.; project administration, R.C. All authors have read and agreed to the published version of the manuscript.

Funding

At the time this study was conducted, O.R. was a graduate student in the Tulane University Bioinnovation Program, supported by NSF IGERT Grant NSF DGE-1144646. O.R. received additional funding from the NIH fellowship grant NIH F31 DK121445.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. General setting of the problems considered in this article. The figure shows a section of a circular permeable membrane. The top panel shows a neighborhood of the membrane and the velocity field across the membrane from the inside to the outside. This flow is generated by a distribution of source doublets on the membrane, [Forumla omitted. See PDF.] with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] is the outward unit normal. The bottom panel shows the vectors [Forumla omitted. See PDF.] (black arrows) and the divergence of the velocity (colormap).

Figure 2. Numerical solution of the ODEs for [Forumla omitted. See PDF.] for the blobs [Forumla omitted. See PDF.] (black) and [Forumla omitted. See PDF.] (blue). The dashed curve is the exact solution of the problem. The left panel shows results for a constant regularization [Forumla omitted. See PDF.]. The right panel shows results when the regularization size is varied dynamically according to [Forumla omitted. See PDF.].

Figure 3. Computation of the radius of the circle using the method in Equation (13) and the two blobs [Forumla omitted. See PDF.] (dash-dot) and [Forumla omitted. See PDF.] (solid) from Table 1. The dashed curve is the exact solution of the problem. The radius is shown in logarithmic scale to appreciate the difference between the solutions. The inset shows the results in linear scale near [Forumla omitted. See PDF.], when the exact solution shrinks to a point.

Figure 4. Errors in the circle radius in example 2 using the blob [Forumla omitted. See PDF.]. The (left) panel shows the errors plotted versus the discretization size [Forumla omitted. See PDF.]. The (right) panel shows the errors plotted versus [Forumla omitted. See PDF.].

Figure 5. Numerical solution of example 3 for the blob [Forumla omitted. See PDF.] in Table 1 using permeability parameter [Forumla omitted. See PDF.], regularization parameter [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] = 1. The initial membrane shape was the ellipse of major axis 1 and minor axis 1/2. The figure shows the computed solution at times [Forumla omitted. See PDF.].

Figure 6. Numerical solution of the pressure in example 3 at time [Forumla omitted. See PDF.] for the blob [Forumla omitted. See PDF.] in Table 1 using permeability parameter [Forumla omitted. See PDF.], regularization parameter [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.]. The left panel shows the pressure in a region containing the membrane. The middle panel shows the pressure along the line [Forumla omitted. See PDF.]. The right panel shows contours of the pressure.

Figure 7. Schematic of the channel. Parabolic flow given by [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.] enters the channel at the left inlet. The height of the channel is [Forumla omitted. See PDF.] and the length is [Forumla omitted. See PDF.]. [Forumla omitted. See PDF.], the portion of the top wall within [Forumla omitted. See PDF.], is permeable and is represented by a dashed boundary; the rest of the top wall and the entire bottom wall ([Forumla omitted. See PDF.]) are solid. The right outlet is open. Dotted lines show potential fluid particle paths as the fluid interacts with the channel walls.

Figure 8. Computed velocity boundary conditions (arrows) and velocity field streamlines for different values of [Forumla omitted. See PDF.]. The colors represent the magnitude of the velocity in the channel. The numerical parameters used are [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 9. Advection of concentration of plasma protein, [Forumla omitted. See PDF.], in a permeable channel with a constant velocity profile at the left boundary. Sources of concentration include an initial gaussian of concentration in the channel which is then added to by a constant protein concentration introduced at the left boundary. Red arrows indicate the velocity due to the hydrostatic pressure, [Forumla omitted. See PDF.], while green arrows indicate that due to the colloid osmotic pressure [Forumla omitted. See PDF.]. Thus, these velocities are greatest where the concentration is largest, and change with the transport of the concentration. Cyan indicates the net flow due to both pressures in addition to the force exerted by the fluid moving with the velocity at the left boundary. The net velocity is scaled by 5 so as to be visible, when it is dwarfed in comparison to the velocities due to [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. All three velocities are scaled by 2 for visibility.

Figure 10. Filtration fraction, [Forumla omitted. See PDF.] as a function of the hydraulic conductivity k of the vessel wall which ultimately dictates the value of [Forumla omitted. See PDF.].

Figure 11. Tangential forces (left) and normal forces (right) on the vessel length for varied values of hydraulic conductivity k.

Table 1

Summary of regularizing functions. The blob $${\varphi }_{\delta }\left(r\right)$$ satisfies only the condition that its integral be one. The function $$\psi \left(r\right)$$ satisfies two more conditions for higher accuracy as described in example 1.

$${\displaystyle {\mathit{\varphi }}_{\mathit{\delta }}\left(\mathit{r}\right)=\frac{2{\mathit{\delta }}^4}{\mathit{\pi }{({\mathit{r}}^2+{\mathit{\delta }}^2)}^3}}$$	$${\displaystyle {\mathit{\psi }}_{\mathit{\delta }}\left(\mathit{r}\right)=\frac{2{\mathit{\delta }}^4({\mathit{r}}^4-10{\mathit{\delta }}^2{\mathit{r}}^2+5{\mathit{\delta }}^4)}{\mathit{\pi }{({\mathit{r}}^2+{\mathit{\delta }}^2)}^5}}$$

Table 2

Errors in the circle radius in example 2 using the blob $${\varphi }_{\delta }$$. The errors were computed as $$ma{x}_t|R\left(t\right)-R_{exact}\left(t\right)|$$. The initial arc length discretization is defined as $$\Delta s_0=2\pi R_0/N$$. The underlined values correspond to errors for a fixed value of $$\delta$$.

N	$$\mathit{\delta }=0.5\Delta {\mathit{s}}_0$$	$$\mathit{\delta }=\Delta {\mathit{s}}_0$$	$$\mathit{\delta }=2\Delta {\mathit{s}}_0$$	$$\mathit{\delta }=4\Delta {\mathit{s}}_0$$	$$\mathit{\delta }=8\Delta {\mathit{s}}_0$$
50	$$1.3241\times {10}^{-2}$$	$$5.2417\times {10}^{-4}$$	$$9.9343\times {10}^{-4}$$	$$4.2410\times {10}^{-3}$$	$$1.8927\times {10}^{-2}$$
100	$$1.3265\times {10}^{-2}$$	$$6.6654\times {10}^{-4}$$	$$2.4256\times {10}^{-4}$$	$$9.9499\times {10}^{-4}$$	$$4.2410\times {10}^{-3}$$
200	$$1.3271\times {10}^{-2}$$	$$7.0338\times {10}^{-4}$$	$$5.9440\times {10}^{-5}$$	$$2.4380\times {10}^{-4}$$	$$9.9499\times {10}^{-4}$$
400	$$1.3272\times {10}^{-2}$$	$$7.1265\times {10}^{-4}$$	$$1.3986\times {10}^{-5}$$	$$6.0624\times {10}^{-5}$$	$$2.4380\times {10}^{-4}$$
800	$$1.3273\times {10}^{-2}$$	$$7.1497\times {10}^{-4}$$	$$2.7275\times {10}^{-6}$$	$$1.5135\times {10}^{-5}$$	$$6.0624\times {10}^{-5}$$

Table 3

Errors in the circle radius in example 2 using the blob $${\psi }_{\delta }$$. The errors were computed as $$ma{x}_t|R\left(t\right)-R_{exact}\left(t\right)|$$. The initial arc length discretization is defined as $$\Delta s_0=2\pi R_0/N$$. The underlined values correspond to errors for a fixed value of $$\delta$$.

N	$$\mathit{\delta }=0.5\Delta {\mathit{s}}_0$$	$$\mathit{\delta }=\Delta {\mathit{s}}_0$$	$$\mathit{\delta }=2\Delta {\mathit{s}}_0$$	$$\mathit{\delta }=4\Delta {\mathit{s}}_0$$	$$\mathit{\delta }=8\Delta {\mathit{s}}_0$$
50	$$4.4493\times {10}^{-2}$$	$$6.4214\times {10}^{-3}$$	$$5.0403\times {10}^{-5}$$	$$7.1866\times {10}^{-5}$$	$$1.2301\times {10}^{-3}$$
100	$$4.4493\times {10}^{-2}$$	$$6.3945\times {10}^{-3}$$	$$4.5526\times {10}^{-5}$$	$$1.5954\times {10}^{-6}$$	$$7.1866\times {10}^{-5}$$
200	$$4.4493\times {10}^{-2}$$	$$6.3877\times {10}^{-3}$$	$$4.4542\times {10}^{-5}$$	$$2.7675\times {10}^{-8}$$	$$1.5951\times {10}^{-6}$$
400	$$4.4493\times {10}^{-2}$$	$$6.3861\times {10}^{-3}$$	$$4.4301\times {10}^{-5}$$	$$1.0104\times {10}^{-9}$$	$$2.7490\times {10}^{-8}$$
800	$$4.4493\times {10}^{-2}$$	$$6.3856\times {10}^{-3}$$	$$4.4240\times {10}^{-5}$$	$$7.8814\times {10}^{-10}$$	$$5.6394\times {10}^{-10}$$

Table 4

Computed flow rates at inlet, top wall, and outlet for different values of $$\beta$$. The numerical parameters used are $$h=1/160$$ and $$\delta =0.0224$$.

Flow Rate	$$\mathit{\beta }=-0.0001$$	$$\mathit{\beta }=-0.00018$$	$$\mathit{\beta }=-0.00024$$
inlet	$$-0.66667$$	$$-0.66667$$	$$-0.66667$$
top	$$0.36402$$	$$0.65678$$	$$0.87733$$
outlet	$$0.30265$$	$$0.009889$$	$$-0.21067$$

Table 5

Parameters for glomerular capillary simulations.

Parameter	Value	Units	Ref
Plasma velocity, $$u_0$$	0.09	cm s$${}^{-1}$$	[45]
Capillary length, L	30	$$\mathsf{\mu }$$m	[45]
Capillary diameter, H	10	$$\mathsf{\mu }$$m	[45]
Plasma viscosity, $$\mathsf{\mu }$$	1.24	cP	[45]
Plasma protein diffusivity, $$\sigma$$	68	$$\mathsf{\mu }$$m$${}^2$$ s$${}^{-1}$$	[46]
Hydraulic conductivity, k	$$5\times {10}^{-3}$$	nL min$${}^{-1}$$ mmHg$${}^{-1}$$ $$\mathsf{\mu }$$m$${}^{-2}$$
Hydrostatic pressure, $$P_{hyd}$$	39	mmHg	[47]

References

1. Stillwell, W.; Ehringer, W.; Jenski, L.J. Docosahexaenoic acid increases permeability of lipid vesicles and tumor cells. Lipids; 1993; 28, pp. 103-108. [DOI: https://dx.doi.org/10.1007/BF02535772] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/8441334]

2. Yang, N.J.; Hinner, M.J. Getting across the cell membrane: An overview for small molecules, peptides, and proteins. Methods Mol. Biol.; 2015; 1266, pp. 29-53. [DOI: https://dx.doi.org/10.1007/978-1-4939-2272-7-3] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25560066]

3. Michel, C.; Curry, F. Microvascular permeability. Physiol. Rev.; 1999; 79, pp. 703-761. [DOI: https://dx.doi.org/10.1152/physrev.1999.79.3.703] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10390517]

4. Nagy, J.A.; Benjamin, L.; Zeng, H.; Dvorak, A.M.; Dvorak, H.F. Vascular permeability, vascular hyperpermeability and angiogenesis. Angiogenesis; 2008; 11, pp. 109-119. [DOI: https://dx.doi.org/10.1007/s10456-008-9099-z] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18293091]

5. Pennell, J.P.; Lacy, F.B.; Jamison, R.L. An in vivo study of the concentrating process in the descending limb of Henle’s loop. Kidney Int.; 1974; 5, pp. 337-347. [DOI: https://dx.doi.org/10.1038/ki.1974.49] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/4427415]

6. Knepper, M.A.; Danielson, R.A.; Saidel, G.M.; Post, R.S. Quantitative analysis of renal medullary anatomy in rats and rabbits. Kidney Int.; 1977; 12, pp. 313-323. [DOI: https://dx.doi.org/10.1038/ki.1977.118] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/604620]

7. Layton, A.T. Mathematical modeling of kidney transport. Wiley Interdiscip. Rev. Syst. Biol. Med.; 2013; 5, pp. 557-573. [DOI: https://dx.doi.org/10.1002/wsbm.1232]

8. Jackson, E.A.; Hillmyer, M.A. Nanoporous membranes derived from block copolymers: From drug delivery to water filtration. ACS Nano; 2010; 4, pp. 3548-3553. [DOI: https://dx.doi.org/10.1021/nn1014006]

9. Baronas, R.; Ivanauskas, F.; Kaunietis, I.; Laurinavicius, V. Mathematical Modeling of Plate- gap Biosensors with an Outer Porous Membrane. Sensors; 2006; 6, pp. 727-745. [DOI: https://dx.doi.org/10.3390/s6070727]

10. Hurk, R.v.d.; Evoy, S. A review of membrane-based biosensors for pathogen detection. Sensors; 2015; 15, pp. 14045-14078. [DOI: https://dx.doi.org/10.3390/s150614045]

11. Stamatialis, D.F.; Papenburg, B.J.; Girones, M.; Saiful, S.; Bettahalli, S.N.; Schmitmeier, S.; Wessling, M. Medical applications of membranes: Drug delivery, artificial organs and tissue engineering. J. Membr. Sci.; 2008; 308, pp. 1-34. [DOI: https://dx.doi.org/10.1016/j.memsci.2007.09.059]

12. Ugrozov, V.V.; Elkina, I.B. Mathematical modeling of influence of porous structure a membrane on its vapour-conductivity in the process of membrane distillation. Desalination; 2002; 147, pp. 167-171. [DOI: https://dx.doi.org/10.1016/S0011-9164(02)00528-3]

13. Yun, M.A.; Yeon, K.M.; Park, J.S.; Lee, C.H.; Chun, J.; Lim, D.J. Characterization of biofilm structure and its effect on membrane permeability in MBR for dye wastewater treatment. Water Res.; 2006; 40, pp. 45-52. [DOI: https://dx.doi.org/10.1016/j.watres.2005.10.035] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16360189]

14. Kim, W.K.; Kanduč, M.; Roa, R.; Dzubiella, J. Tuning the Permeability of Dense Membranes by Shaping Nanoscale Potentials. Phys. Rev. Lett.; 2019; 122, 108001. [DOI: https://dx.doi.org/10.1103/PhysRevLett.122.108001]

15. Kim, W.K.; Chudoba, R.; Milster, S.; Roa, R.; Kanduč, M.; Dzubiella, J. Tuning the selective permeability of polydisperse polymer networks. Soft Matter; 2020; 16, pp. 8144-8154. [DOI: https://dx.doi.org/10.1039/D0SM01083A]

16. Holmqvist, P.; Mohanty, P.S.; Nägele, G.; Schurtenberger, P.; Heinen, M. Structure and Dynamics of Loosely Cross-Linked Ionic Microgel Dispersions in the Fluid Regime. Phys. Rev. Lett.; 2012; 109, 048302. [DOI: https://dx.doi.org/10.1103/PhysRevLett.109.048302] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23006114]

17. Daddi-Moussa-Ider, A.; Gekle, S. Brownian motion near an elastic cell membrane: A theoretical study. Eur. Phys. J. E; 2018; 41, 19. [DOI: https://dx.doi.org/10.1140/epje/i2018-11627-6]

18. Bächer, C.; Gekle, S. Computational modeling of active deformable membranes embedded in three-dimensional flows. Phys. Rev. E; 2019; 99, 062418. [DOI: https://dx.doi.org/10.1103/PhysRevE.99.062418]

19. Herschlag, G.; Liu, J.G.; Layton, A.T. An Exact Solution for Stokes Flow in a Channel with Arbitrarily Large Wall Permeability. SIAM J. Appl. Math.; 2014; 75, pp. 2246-2267. [DOI: https://dx.doi.org/10.1137/140995854]

20. Herschlag, G.; Liu, J.G.; Layton, A. Fluid extraction across pumping and permeable walls in the viscous limit. Phys. Fluids; 2016; 28, 041902. [DOI: https://dx.doi.org/10.1063/1.4946005]

21. Stockie, J.M. Modelling and simulation of porous immersed boundaries. Comput. Struct.; 2009; 87, pp. 701-709. [DOI: https://dx.doi.org/10.1016/j.compstruc.2008.11.001]

22. Layton, A.T. Modeling Water Transport across Elastic Boundaries Using an Explicit Jump Method. SIAM J. Sci. Comput.; 2006; 28, pp. 2189-2207. [DOI: https://dx.doi.org/10.1137/050642198]

23. Kim, Y.; Peskin, C. 2–D Parachute Simulation by the Immersed Boundary Method. SIAM J. Sci. Comput.; 2006; 28, pp. 2294-2312. [DOI: https://dx.doi.org/10.1137/S1064827501389060]

24. Natali, D.; Pralits, J.O.; Mazzino, A.; Bagheri, S. Stabilizing effect of porosity on a flapping filament. J. Fluids Struct.; 2016; 61, pp. 362-375. [DOI: https://dx.doi.org/10.1016/j.jfluidstructs.2015.11.016]

25. Mayo, A. The fast solution of Poisson’s and the biharmonic equations on irregular regions. SIAM J. Numer. Anal.; 1984; 21, pp. 285-299. [DOI: https://dx.doi.org/10.1137/0721021]

26. LeVeque, R.J.; Li, Z. The Immersed Interface Method for Elliptic Equations With Discontinuous Coefficients and Singular Sources. SIAM J. Numer. Anal.; 1992; 31, pp. 332-364.

27. Pepona, M.; Favier, J. A coupled Immersed Boundary–Lattice Boltzmann method for incompressible flows through moving porous media. J. Comput. Phys.; 2016; 321, pp. 1170-1184. [DOI: https://dx.doi.org/10.1016/j.jcp.2016.06.026]

28. Berman, A.S. Laminar flow in channels with porous walls. J. Appl. Phys.; 1953; 24, pp. 1232-1235. [DOI: https://dx.doi.org/10.1063/1.1721476]

29. Yuan, S. Laminar Pipe Flow with Injection and Suction through a Porous Wall; Technical Report James Forrestal Research Center, Princeton University: Princeton, NJ, USA, 1955.

30. Wah, T. Laminar flow in a uniformly porous channel. Aeronaut. Q.; 1964; 15, pp. 299-310. [DOI: https://dx.doi.org/10.1017/S0001925900010908]

31. Terrill, R.M.; Shrestha, G.M. Laminar flow through parallel and uniformly porous walls of different permeability. Z. Angew. Math. Phys. ZAMP; 1965; 16, pp. 470-482. [DOI: https://dx.doi.org/10.1007/BF01593923]

32. Cortez, R. The method of regularized Stokeslets. SIAM J. Sci. Comput.; 2001; 23, pp. 1204-1225. [DOI: https://dx.doi.org/10.1137/S106482750038146X]

33. Cortez, R.; Fauci, L.; Medovikov, A. The method of regularized Stokeslets in Three Dimensions: Analysis, Validation and Application to Helical Swimming. Phys. Fluids; 2005; 17, 031504. [DOI: https://dx.doi.org/10.1063/1.1830486]

34. Beale, J.T. A Convergent Boundary Integral Method for Three-Dimensional Water Waves. Math. Comput.; 2001; 70, pp. 977-1029. [DOI: https://dx.doi.org/10.1090/S0025-5718-00-01218-7]

35. Tlupova, S.; Beale, J.T. Regularized single and double layer integrals in 3D Stokes flow. J. Comput. Phys.; 2019; 386, pp. 568-584. [DOI: https://dx.doi.org/10.1016/j.jcp.2019.02.031]

36. Shea, S.M. Glomerular hemodynamics and vascular structure: The pattern and dimensions of a single rat glomerular capillary network reconstructed from ultrathin sections. Microvasc. Res.; 1979; 18, pp. 129-143. [DOI: https://dx.doi.org/10.1016/0026-2862(79)90023-2]

37. Arendshorst, W.; Navar, L. Renal circulation and glomerular hemodynamics. Diseases of the Kidney; 8th ed. Schrier, R. Walters Kluwer/Lippincott Williams and Wilkins: Philadelphia, PA, USA, 2007; Chapter 2 1, pp. 54-95.

38. Deen, W.; Robertson, C.; Brenner, B. A model of glomerular ultrafiltration in the rat. Am. J. Physiol.-Leg. Content; 1972; 223, pp. 1178-1183. [DOI: https://dx.doi.org/10.1152/ajplegacy.1972.223.5.1178] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/4654350]

39. Deen, W.M.; Bohrer, M.P.; Brenner, B.M. Macromolecule transport across glomerular capillaries: Application of pore theory. Kidney Int.; 1979; 16, pp. 353-365. [DOI: https://dx.doi.org/10.1038/ki.1979.138]

40. Deen, W.M.; Lazzara, M.J.; Myers, B.D. Structural determinants of glomerular permeability. Am. J. Physiol.-Ren. Physiol.; 2001; 281, pp. F579-F596. [DOI: https://dx.doi.org/10.1152/ajprenal.2001.281.4.F579]

41. Lambert, P.; Aeikens, B.; Bohle, A.; Hanus, F.; Pegoff, S.; Van Damme, M. A network model of glomerular function. Microvasc. Res.; 1982; 23, pp. 99-128. [DOI: https://dx.doi.org/10.1016/0026-2862(82)90035-8]

42. Remuzzi, A.; Brenner, B.M.; Pata, V.; Tebaldi, G.; Mariano, R.; Belloro, A.; Remuzzi, G. Three-dimensional reconstructed glomerular capillary network: Blood flow distribution and local filtration. Am. J. Physiol.-Ren. Physiol.; 1992; 263, pp. F562-F572. [DOI: https://dx.doi.org/10.1152/ajprenal.1992.263.3.F562]

43. de Sousa, D.; Ferreira, M.C. Ultrafiltration in Renal Glomerular Capillaries: Theoretical Effects of Ultrastructure. Ph.D. Thesis; Massachusetts Institute of Technology: Cambridge, MA, USA, 1994.

44. Layton, A.T.; Edwards, A. Glomerular Filtration. Mathematical Modeling in Renal Physiology; Springer: Berlin/Heidelberg, Germany, 2014; pp. 7-41.

45. Ferrell, N.; Sandoval, R.M.; Bian, A.; Campos-Bilderback, S.B.; Molitoris, B.A.; Fissell, W.H. Shear stress is normalized in glomerular capillaries following 5/6 nephrectomy. Am. J. Physiol.-Ren. Physiol.; 2015; 308, pp. F588-F593. [DOI: https://dx.doi.org/10.1152/ajprenal.00290.2014]

46. Wakeham, W.; Salpadoru, N.; Caro, C. Diffusion coefficients for protein molecules in blood serum. Atherosclerosis; 1976; 25, pp. 225-235. [DOI: https://dx.doi.org/10.1016/0021-9150(76)90029-0]

47. Zatz, R.; Dunn, B.R.; Meyer, T.W.; Anderson, S.; Rennke, H.G.; Brenner, B.M. Prevention of diabetic glomerulopathy by pharmacological amelioration of glomerular capillary hypertension. J. Clin. Investig.; 1986; 77, pp. 1925-1930. [DOI: https://dx.doi.org/10.1172/JCI112521] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/3011862]

48. Richfield, O.; Cortez, R.; Navar, L.G. Simulations of increased glomerular capillary wall strain in the 5/6-nephrectomized rat. Microcirculation; 2021; 28, e12721. [DOI: https://dx.doi.org/10.1111/micc.12721] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34192389]

49. Cortez, R. Regularized Stokeslet segments. J. Comput. Phys.; 2018; 375, pp. 783-796. [DOI: https://dx.doi.org/10.1016/j.jcp.2018.08.055]