Computing viscous flow along a 2D open channel using the immersed interface method

We present a numerical method for simulating 2D flow through a channel with deformable walls. The fluid is assumed to be incompressible and viscous. We consider the highly viscous regime, where fluid dynamics are described by the Stokes equations, and the less viscous regime described by the Navier–Stokes equations. The model is formulated as an immersed boundary problem, with the channel defined by compliant walls that are immersed in a larger computational fluid domain. The channel traverses through the computational domain, and the walls do not form a closed region. When the walls deviate from their equilibrium position, they exert singular forces on the underlying fluid. We compute the numerical solution to the model equations using the immersed interface method, which preserves sharp jumps in the solution and its derivatives. The immersed interface method typically requires a closed immersed interface, a condition that is not met by the present configuration. Thus, a contribution of the present work is the extension of the immersed interface method to immersed boundary problems with open interfaces. Numerical results indicate that this new method converges with second‐order accuracy in both space and time, and can sharply capture discontinuities in the fluid solution.


INTRODUCTION
A detailed and accurate description of fluid flow through a channel with compliant or actively moving walls is of interest in many biological applications. Small changes in the diameter of the channel can cause drastically different fluid behavior. For example, the main mechanism for blood pumping in a valveless heart is thought to be peristalsis, a wave-like contraction. [1][2][3] Deformable channels or tubes are also used to model obstructed ureters, 4 food mixing in the intestine, 5 and blood flow through a vessel 6,7 or renal tubules. 8,9 A natural way to model flow along a compliant channel is to frame it as an immersed boundary problem. [10][11][12][13][14] Immersed boundary problems are typically a subset of fluid-structure interaction problems in which a thin structure or physical boundary is present in the fluid. 15 The immersed boundary applies a singular force to the underlying fluid, which is often described by Navier-Stokes flow, but can be described by Stokes flow for highly viscous fluids. The immersed boundary method transfers the singular boundary forces onto the underlying fluid using approximate (smooth) Dirac delta functions typically with (h) support. Because the delta function is smoothed, this approach does not capture the jump discontinuity in the solution (e.g., pressure) at the immersed boundary but rather approximates the solution as a continuous function with a large gradient. In general, the immersed boundary method computes approximations with first-order spatial accuracy.
If better accuracy is desired, especially near the immersed boundary, one may use the immersed interface method developed by LeVeque and Li. 16,17 The immersed interface method captures the sharp jumps in the solution and its derivatives and generates approximations with second-order accuracy. The key idea in the immersed interface method is the incorporation of known jumps in the solution or its derivatives into the finite difference schemes. The immersed interface method has been applied to a variety of applications to compute the coupled motion of a viscous fluid and a thin closed interface with second-order accuracy. [18][19][20][21][22] The immersed interface method requires that the immersed interfaces be closed. Thus, to simulate flow through a channel, one typically approximates the channel as a closed, elongated interface with capped ends, immersed in a larger computational domain. 10,21,[23][24][25][26] The flow is then driven by a pair of fluid source and sink located at the opposing ends of the channel. A major downside of this setup is the unrealistic flow near the source and sink. As a remedy, a longer channel is modeled, and only the flow sufficiently far from the source and sink is considered.
This study develops a modified version of the immersed interface method that can be applied directly to an open channel without the representation of (unrealistic) fluid source and sink. To accomplish this, we derive the jump conditions for the open interface and apply the method to the Stokes equations as well as the Navier-Stokes equations. Numerical examples indicate that the method captures the sharp jumps in the solutions and achieves second-order accuracy in both time and space.

Computational domain and immersed interface
We formulate a model that simulates fluid flow through an open channel with compliant walls extending from one side of a rectangular computational domain to the opposite side. The model is formulated in 2D rectangular coordinates. To define the model, consider a rectangular domain given by Ψ = [0, L] × [−H, H] with length L and width 2H. Let Γ represent the compliant walls of the channel immersed in Ψ. Specifically, Γ is an interface immersed in the fluid formed by two distinct smooth curves X = G 1 (s) for s ∈ [0, L 1 ] such that L 1 is the resting length of the curve G 1 which intersects Ψ at (0, a 1 ) and (L, b 1 ) and X = G 2 (s) for s ∈ [0, L 2 ] such that L 2 is the resting length of G 2 which intersects Ψ at (0, a 2 ) and (L, b 2 ), as shown in Figure 1. We are interested in tracking flow through the channel (i.e., between G 1 and G 2 ), driven by a pressure gradient between x = 0 and x = L. Let n be the unit vector normal to Γ oriented toward Ψ away from the channel. The interface Γ consists of two separate curves given by G 1 and G 2 that intersect the x = 0 and x = L on Ψ at (x, y) = (0, a 1 ), (x, y) = (0, b 1 ), and (x, y) = (0, a 2 ), (x, y) = (0, b 2 ), respectively. The unit normal, defined to be n, is oriented toward Ψ

Boundary conditions
To properly state the problem, appropriate boundary conditions must be imposed on fluid velocity and pressure. Biperiodic boundary conditions are assumed for fluid velocity, which implies that the volume of the fluid in the channel remains constant in time. To drive flow, we prescribe a pressure gradient inside the channel. For simplicity, we assume that a 1 = b 1 and a 2 = b 2 . Thus, the width of the inlet and outlet is the same. We require there to be a constant difference in pressure, denoted P diff , at the inlet and outlet of the channel The derivatives of pressure on the x = 0 and x = L boundaries are required to be equal We will refer to this setup as "inhomogeneous periodic boundary conditions." Not only does this setup generate a pressure gradient along the channel, but it also forces the pressure gradient to be periodic in x, which is consistent with the periodic boundary conditions imposed for velocity. Periodic boundary conditions are imposed for pressure at the y = ±H boundaries.
Note that boundary conditions at the inflow and outflow of a channel are a rich research area. [27][28][29] These boundary conditions were chosen for simplicity, but this method can be used with a wide range of boundary conditions that would allow for the accumulation and release of fluid in the channel.

Fluid-structure interactions
The Navier-Stokes equations are given by where u = (u, v), where u and v denote fluid velocity components in the x and y directions, respectively; p is the pressure; is viscosity; and F is the interfacial force, which is singularly supported along Γ (see below). We consider also the zero Reynolds number regimes given by the Stokes momentum equation: and the continuity equation (Equation (4)). The force exerted by the interface Γ can be written as where X denotes the position of the interface, = (s, t) is the material coordinate(s) that parameterize the interface curve at time t, f( ) is the force strength at the point X( ), and is the Dirac delta function.
The force strength f( ) from Equation (6) has two major components, an elastic force f E and a tether force f T : Since the interface is elastic, any deviation from its resting configuration generates a restorative force. The elastic force can be modeled using Hooke's law where the unit tangent vector to Γ is given by The tension T(s, t) in Equation (8) is given by where a E controls the stiffness of the interface. The interface is tethered to a set of anchor points evenly distributed along the length of the interface. Interface movement can either be restricted by using stationary anchor points or be induced by moving the position of the anchor points.
The second force component f T arises from the displacement of the interface control knots from their anchor points.
Suppose the interface X is discretized by a set of interface control knots by a spring with resting length 0. Figure 2 shows three interface control knots X i − 1 , X i , and X i + 1 connected to the tether points X i−1 , X i , and X i+1 , respectively, by a spring. Let a T be the spring force constant. Then the tether force is given by Thus, if the interface knots move away from their anchor points, restorative forces are generated. The interface is deformable and moves at the same speed as the local fluid. The no-slip condition describes this motion.

Jump conditions
The derivation of the jump conditions normally requires that the interface be a closed curve. Here, we extend the derivation to an open channel. Details are explained in the Appendix. Briefly, we extend Γ to a fictitious closed and piecewise smooth curve Γ c by taking the union of Γ, the line joining the points {(0, a 1 ), (0, a 2 )}, and the line joining the points {(L, b 1 ), (L, b 2 )}. The rest of the derivation is similar to the closed curve case. The jump conditions given in Equations (13)- (16) are the same for the Stokes equations and for the Navier-Stokes equations. This follows from the no-slip condition assumed for the interfacial motion, which generates a continuous velocity field and material derivative across the interface. 30 [p] = f ⋅ n (13) F I G U R E 2 A section of Γ is discretized using the interface control knots X i − 1 , X i , and X i + 1 . The points X i − 1 , X i , and X i + 1 are connected to the tether points X i−1 , X i , and X i+1 , respectively, by springs (waving lines). These springs generate the tether force f T where the jump [⋅] of a quantity, say q(X, t), at time t and point X ∈ Γ is Below we refer to this extended immersed interface method as simply the immersed interface method.

Stokes problem
To solve the Stokes equations with a singular force (Equations (4) and (5)), we compute the fluid velocity and pressure on a fixed Eulerian grid: where h is the grid spacing and N x = L h + 1 and N y = H h + 1 are the number of grid points in the x and y directions, respectively. A moving Lagrangian frame of reference is used to track the location of the interface Γ. The interface position that are connected by a cubic spline. 31 The Stokes solution is computed by applying the immersed interface method to solve a series of three Poisson equations: one obtained by taking the divergence of Equation (5) to yield the other two being the two velocity components of Equation (4). Specifically, Equations (4) and (18) are discretized using second-order finite difference, with jump conditions Equations (13)-(16) incorporated into the finite difference stencil as correction terms. For details on computing the correction terms, see Reference 16.

Navier-Stokes problem
A large number of correction terms are involved in the application of the immersed interface method directly to the Navier-Stokes equations. As an alternative, we adopt the velocity decomposition approach. 32 This approach leverages the fact that, for a given singular interfacial force, the jumps in the solutions (see Equations (13)-(16)) are identical for both the Stokes equations and the Navier-Stokes equations. 16,17 To apply the velocity decomposition approach, we split the Navier-Stokes equations into two parts: a singular part that satisfies the Stokes equations including the singular force (denoted by u s and p s ), and a regular part ( u r and p r ): Taking the difference between Equations (3) and (5), we determine that the regular part, u r and p r , satisfies where F b is a body force Note that the full velocity u is used in the transport of u r and in F b . Then substituting Equation (19) into Equation (4) yields To advance the full solution from t n to t n + 1 , the Stokes part (u n+1 s , p n+1 s ) is first computed, following the procedures described in Sect. 3.1. This allows us to update F n+1 b by integrating along the fluid trajectory backward in time (more below). Then the regular solution (u n+1 r , p n+1 r ) is computed. Equations 20 and 22 are essentially the Navier-Stokes equations with a body force and are solved using the projection method. That is, we first compute the intermediate solution u * ,n+1 r using Equation (20) alone and then project u * ,n+1 r onto the divergence free space to yield u n+1 r and p n+1 r . To compute u * ,n+1 r , we first rewrite Equation (20) in terms of the material derivative as follows Note that u r and its material derivatives are smooth along the fluid trajectories. This motivates us to use the semi-Lagrangian time-discretization method. The semi-Lagrangian method computes the solution at fixed Eulerian grid points x by integrating the solution backward in time along the trajectories of fluid particles that pass through those grid points at t n + 1 . In our case, the semi-Lagrangian discretization is applied to the material derivative Du r /Dt. To proceed, we first compute the upstream positions of the particle, x n at t n and x n − 1 at t n − 1 , which can be performed by integrating backward in time. Once x n and x n − 1 are found, we compute u r at these locations and times. It is unlikely that the upstream positions coincide with grid points. Therefore, these upstream velocity valuesũ n r = u r (x n , t n ) andũ n−1 r = u r (x n−1 , t n−1 ) are approximated using spatial interpolation. The material derivative is then approximated using the second-order backward difference formula, The intermediate solution u * ,n+1 r is then projected, P(u * ,n+1 r ) = u n+1 r , into the subspace of divergence-free vector fields so that ∇ ⋅ u n+1 r = 0. Specifically, we solve and update u n+1 r and p n+1 The full solution is given by Equation (19). Details for the velocity decomposition method can be found in Reference 32.

Boundary condition and solution decomposition
In the velocity decomposition method, boundary conditions must be imposed for both the Stokes and regular parts such that the boundary conditions for the full solution are still satisfied. Recall that biperiodic boundary conditions are prescribed for the full velocity solution. This is achieved by imposing biperiodic boundary conditions for both u s and u r as well as the intermediate solution u * r .
The full pressure solution must satisfy the inhomogeneous periodic boundary conditions given by Equations (1) and (2). This is achieved by imposing the inhomogeneous periodic conditions on p s along x and periodic boundary conditions along y, and by imposing the regular biperiodic boundary conditions on p r . Biperiodic boundary conditions are assumed for the auxiliary variable .

Stokes results
We first conduct a grid refinement study on the Stokes problem to test the spatial convergence of the method in the zero Reynolds number regime. In this test, the interface is initially displaced and allowed to return to its resting position as respectively. The pressure drop along the channel (in x) is taken to be P diff = 2.544 gm/(s 2 μm). Model parameters for this simulation are presented in Table 1.
For the Stokes problem, the fluid solutions are steady and only indirectly depend on time because the position of the interface may be time dependent. Therefore, we conducted a spatial grid refinement study at t = 0. This goal of this study is to test the spatial convergence of method without added temporal errors (the spatial convergence with temporal errors are presented in Table 4). Computed solutions were compared with an N x = 1089 by N y = 1666 high-resolution numerical solution. The order of convergence for a solution component w is found using the formula where w* is the N x = 1089 by N y = 1666 high-resolution solution, w h is the solution computed with spatial grid size h, w h 2 is the solution computed with spatial grid size h 2 , and || ⋅ || ∞ is the infinity norm. Results in Table 2 indicate that second-order spatial convergence was achieved for pressure and both components of velocity. Furthermore, the jump discontinuities in p and in the normal derivative of u are captured robustly. Solution profiles are not shown for these results because the sharp-interface feature of the method is also captured and illustrated in the Navier-Stokes example below.

Navier-Stokes results
For the second example, we solve the Navier-Stokes equations to assess the spatial and temporal accuracy in the nonzero Reynolds number regime.

Spatial convergence and steady-state profile
We start with a fluid that is at rest throughout the domain and a channel with uniform width W = 10.1. The axial pressure gradient is then gradually increased during the first 0.125 s and then held constant according to the function: where R = W/2. The channel boundaries are tethered to its initial and resting position, that is, a channel with a width of W. Table 3 shows the remainder of the parameters.

TA B L E 3 Parameters used for the spatial convergence study of the Navier-Stokes problem
At steady state, pressure gradient and velocity are described by 2D Poiseuille flow, that is, It follows that pressure decreases linearly along the channel and u has a parabolic profile given by In Figure 3 a solid line shows the expected Poiseuille velocity profile given the parameters of the simulation. The computed u-component of the velocity profile is shown on three slices along the length of the channel at time t = 0.05 s. The velocity profile on the slice where x = 1, x=16, and x = 29 μm are denoted by the ○, △, and □ symbols, respectively. While the computed u exhibits an approximate parabolic profile along the entire length of the channel, it is particularly noteworthy that the results near the inlet and outlet of the channel exhibit this expected behavior.
The computed p is shown in Figure 4. The interface segments are located at approximately y = ±5.05 μm. It can be seen that p decreases linearly along the channel. Outside of the channel, the pressure is approximately zero. There is a sharp discontinuity in pressure across the interface.
To test the spatial convergence of the method, the solution was computed at various spatial resolutions. The number of interface control points for each segment of the interface is N b = where w* is the N x = 1025 by N y = 513 high-resolution solution, w h is the solution computed with spatial grid size h, w h 2 is the solution computed with spatial grid size h 2 , and || ⋅ || ∞ is the infinity norm. Convergence results indicate that the method is approximately second-order accurate in space; see Table 4.

Temporal convergence
The above configuration is not ideal for assessing temporal convergence in part because the interface position undergoes little deformation. Thus, we consider a different configuration in which the channel walls move. Specifically, the wall tether anchors are displaced in time, causing the walls to actively deform. The y-position of the tethers at position x for the upper and lower interface, denoted Y top, tether (x) and Y bottom, tether (x), respectively, are given by The fluid solutions are computed in the domain with a length of 61 μm and width 21 μm. The remainder of the parameter values for this simulation are given in Table 5. The temporal accuracy of the method is assessed by refining the time

TA B L E 5 Parameters used for the temporal convergence study of
the Navier-Stokes problem step and comparing the solutions at time t = 12.5 μs to a high-resolution solution computed after 64 time steps of size Δt = 1.8e − 7 when the solution is changing at a sufficiently high rate. The order of convergence for a solution component w is found using the formula where w* is the high-resolution solution computed with time steps of size Δt = 1.8e − 7, w Δt is the solution computed with time steps of size Δt, w Δt 2 is the solution computed with time steps of size Δt 2 , and || ⋅ || ∞ is the infinity norm. Convergence results, given in Table 6, indicate that this method can achieve second-order accuracy in time.

Stenosed channel
Finally, flow through a stenosed channel is presented to showcase a potential application of this method. The fluid solutions are computed in the domain with a length of 32 μm and width 16 μm. In this example, the interface walls are tethered in their initial configuration as a straight channel with an indentation centered at x = 21 μm which represents a stenosed channel. The y-position of the tethers at position x at time t for the upper and lower interface, denoted Y top, tether (x, t) and Y bottom, tether (x, t), respectively, are given by The parameters for this application are identical to those presented in Table 3.
To test the spatial convergence of the method, the solution was computed at various spatial resolutions. The number of interface control points for each segment of the interface is N b = N x +1 2 . The solutions are compared with high-resolution solutions computed on a 513 × 257 grid at time t = 0.1 s. The order of convergence for a solution component w is found using the formula where w* is the N x = 513 by N y = 257 high-resolution solution, w h is the solution computed with spatial grid size h, w h 2 is the solution computed with spatial grid size h 2 , and || ⋅ || ∞ is the infinity norm. Convergence results indicate that the method is approximately second-order accurate in space; see Table 7.  The results for this study indicate that even when the interface deviates substantially from a straight channel, the spatial convergence remains above first order, although second-order convergence was not quite achieved.

TA B L E 6
The following images were generated from 129 × 65 resolution solution at time t = 0.1 s. The position of the interface at time t = 0.1 s is shown as a black line in Figure 5. The blue arrows indicate the direction and magnitude of the fluid velocity through the channel at time t = 0.1 s. Notice that flow is fastest near the center of the channel. Above the stenosed channel, the fluid pushes slightly toward the wall then slows down through the stenosed section of the channel. There is very little flow outside of the interface.
The pressure in the computational domain is shown in Figure 6. It is worth noting that the pressure decreases almost linearly throughout the channel. This constant pressure gradient drives the fluid flow through the channel. The pressure is approximately zero outside of the channel and exhibits a sharp jump across the interface.

DISCUSSION
We have presented a numerical method for simulating viscous fluid flow through an open channel with deformable walls. The model is formulated as an immersed boundary problem, with the channel spanning from one end of the computational domain to the other. We apply the method to the Stokes equations and the Navier-Stokes equations. The Stokes equations are solved using a method that is an extension of the immersed interface method, which originally requires the immersed interface to be closed. This method gives second-order accurate values by incorporating known jumps for the solution and its derivatives into a finite difference method. The Navier-Stokes equations are solved using the velocity decomposition approach, which decomposes the velocity into a "Stokes" part and a "regular" part. The first part is determined by the Stokes equations and the singular interfacial force. The regular part of the velocity is given by the Navier-Stokes equations with a body force resulting from the Stokes part. The regular velocity is obtained using a time-stepping method that combines the semi-Lagrangian method with the backward difference formula. Numerical examples are presented to demonstrate the method's performance for both the Stokes and Navier-Stokes problems. In the Stokes model, the interface is displaced from and allowed to return to its resting position as a straight channel. Three numerical cases were run for the Navier-Stokes model: a straight channel, a straight channel that is actively deformed by moving the tether anchor points, and a stenosed channel. These numerical cases indicate that the method converges with second-order spatial and temporal accuracy when the interface undergoes small displacements and near second-order spatial accuracy then the interface substantially deviates from a straight channel (as was the case with the stenosed channel). In that study, there were points of high pressure near the intersection of the interface and the inlet and outlet walls. This may be the cause of the loss of accuracy in the case of large interface displacement. Regardless, this method has many potential applications for flow in biological tubes where the tube is pliable or actively moving, but does not undergo drastic deformation. One such application is the myogenic response of vessels in the microcirculation to sudden changes in blood pressure. 33 All experiments in this article were in the low Reynolds number regime. Future work would also include optimizing the implementation of this method and running experiments with higher Reynolds numbers.
The development of the present method is motivated by our interest in simulating biological problems with flows through biological tubes 34,35 or microfluidic devices. 36 The present model is formulated for 2D. Hence, the flow that it describes may differ significantly from flow through a tube. In future studies, the immersed interface method may be extended to compute fluid flows through a 3D open tube. Additional work is also being conducted on using this method to model congestion in the microcirculation and the myogenic response of vessels in the microcirculation to sudden changes in blood pressure.

PEER REVIEW INFORMATION
Engineering Reports thanks Mehdi Jabbarzadeh and other anonymous reviewers for their contribution to the peer review of this work.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.

CONFLICT OF INTEREST
Authors have no conflict of interest relevant to this article.

R E F E R E N C E S APPENDIX A. DERIVATION OF JUMP CONDITIONS
The derivation of the jump conditions normally requires that the interface be a closed curve. Recall that the interface Γ was formed by two distinct smooth curves that intersect the computational domain boundary Ψ at (0, a 1 ) and (L, b 1 ) and at (0, a 2 ) and (L, b 2 ), respectively. Since Γ is not closed, we extend Γ to a fictitious closed and piecewise smooth curve Γ c by taking the union of Γ, the line joining the points {(0, a 1 ), (0, a 2 )}, and the line joining the points {(L, b 1 ), (L, b 2 )}. See Figure A1. Note that the jump conditions will be the same for the Stokes equations and for the Navier-Stokes equations.
In this section, we use the notation [⋅] to represent the jump discontinuity of a quantity, say q(X, t), at time t and point X ∈ Γ c , where [q(X, t)] = lim →0 + q(X + n, t) − lim →0 − q(X + n, t).
We first derive jump conditions for p. Let (x) be an arbitrary twice continuously differentiable test function defined on the extended domain Ψ = [− , L + ] × [− , H + ] (see Figure A1). Extend F to Γ c by defining where f is a piecewise smooth extension to Γ c . Integrating the product of the divergence of the boundary force F and over the extended fluid domain Ψ , The last line can be obtained via integration by parts and noting that is zero away from Γ c . This calculation holds on any subset of Ψ which contains Γ c . In particular, it holds in the belt domain Ω = {x ∈ Ψ |min y∈Γ c ||x − y|| < }, which encloses the interface Γ c . Let Γ + and Γ − denote the inner and outer boundary of Ω , respectively, as shown in Figure A2. By where is the angle between the outward normal and the x-axis, n = (cos , sin ), and = (− sin , cos ). Solving this linear equation for x and y yields x = n cos − s sin (A13) y = n sin + s cos . Then For any twice continuously differentiable test function , We can conclude that Next, we derive the jump conditions for the u component of velocity. The derivation for v is similar. The velocity is continuous, so we only need to derive the jump conditions for the normal derivative, u n . We start by multiplying the u component of Equation (5) by test function and integrating.