Computing Viscous Flow Along a 3D Open Tube Using the Immerse Interface Method

12/24/2021
by   Sarah E Patterson, et al.
0

In a companion study <cit.>, we present a numerical method for simulating 2D viscous flow through an open compliant closed channel, drive by pressure gradient. We consider the highly viscous regime, where fluid dynamics is described by the Stokes equations, and the less viscous regime described by the Navier-Stokes equations. In this study, we extend the method to 3D tubular flow. The problem is formulated in axisymmetric cylindrical coordinates, an approach that is natural for tubular flow simulations and that substantially reduces computational cost. When the elastic tubular walls are stretched or compressed, they exert forces on the fluid. These singular forces introduce unsmoothness into the fluid solution. As in the companion 2D study <cit.>, we extend the immersed interface method to an open tube, and we compute solution to the model equations using the resulting method. Numerical results indicate that this new method preserves sharp jumps in the solution and its derivatives, and converges with second-order accuracy in both space and time.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

page 16

page 17

page 22

page 23

12/23/2021

Analysis of the diffuse interface method for the Stokes-Darcy coupled problem

We consider the interaction between a free flowing fluid and a porous me...
06/17/2021

The Biot-Stokes coupling using total pressure: formulation, analysis and application to interfacial flow in the eye

We consider a multiphysics model for the flow of Newtonian fluid coupled...
01/22/2021

The art of coarse Stokes: Richardson extrapolation improves the accuracy and efficiency of the method of regularized stokeslets

The method of regularised stokeslets is widely used in microscale biolog...
07/02/2019

Parallel time-stepping for fluid-structure interactions

We present a parallel time-stepping method for fluid-structure interacti...
07/15/2019

Hydrodynamic Simulations using GPGPU Architectures

Simulating the flow of different fluids can be a highly computational in...
09/17/2020

Self-attenuation of extreme events in Navier-Stokes turbulence

Turbulent fluid flows are ubiquitous in nature and technology, and are m...
07/23/2020

An integral-based spectral method for inextensible slender fibers in Stokes flow

Every animal cell is filled with a cytoskeleton, a dynamic gel made of i...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

A detailed description of 3D viscous fluid flow through a compliant or actively moving tube is of interest in many biological applications including pumping via peristalsis in a valveless heart santhanakrishnan2011fluid, food mixing in the intestine tharakan2010mass, and blood flow through a vessel tang2002simulating. Computer simulations that assume rigid walls often fail to predict some essential characteristics of tubular flow, such as pressure wave propagation. Therefore, such simulations cannot be considered reliable in every situation, such as the case when the vessels undergo relatively large displacements.

A natural way to model flow in compliant vessels is to frame it as an immersed boundary problem arthurs1998modeling, peskin2002immersed. Immersed boundary problems are a subset of FSI problems in which a thin structure or physical boundary is present in the fluid hou2012review. The immersed boundary formulation was developed by Charles Peskin to study blood flow in the heart but has been extended to a variety of biological applications including blood clotting, aquatic animal locomotion, fluid dynamics in the inner ear, arteriolar flow, and flow in collapsible tubes arthurs1998modeling, peskin1995general. The popularity of this method is due to its ability to model fluid interactions with complex, passive, or active elastic material with relative ease compared to traditional body fitted approaches. The fluid solution is computed on a fixed Cartesian grid which does not need to conform to the geometry of the structure. The structure can then be represented using Lagrangian variables that can move over the fluid grid unimpeded. The Eulerian and Lagrangian variables are related by interaction equations that contain the delta Dirac function mittal2005immersed, peskin2002immersed.

The immersed boundary method generally only converges with first-order accuracy and does not capture the discontinuities in the velocity gradient and pressure near the interface. Since the vessel movement is determined from the local fluid velocity, these inaccuracies cause the movement of the interface to deteriorate over time and lose volume. Although progress has been made on improving the volume conservation of the immersed boundary method, this problem has not been eliminated peskin1993improved. Additionally, the pressure and velocity near the interface are needed for some applications. A sharp interface method, like the immersed interface method, is required in order to determine these variables accurately.

The immersed interface method overcomes the drawbacks in immersed boundary methods by sharply capturing the discontinuities of fluid solutions and the movement of the immersed structures. It is based on a method used to compute solutions to Poisson’s equation in irregular domains developed by Anita Mayo mayo1984fast. The immersed interface method has been applied to Poisson’s equation leveque1994immersed, the Stokes equations leveque1997immersed, and the Navier-Stokes equations li2001immersed, lee2003immersed. The immersed interface method achieves second-order accuracy in numerical approximation. For fluid-structure interactions problems, the immersed interface method has been used to compute the coupled motion of a viscous fluid and a simple closed elastic interface with second-order accuracy le2006immersed, lee2003immersed, li2014computing, li2001immersed, tan2009immersed, xu2006immersed.

In the original derivation of the jump conditions, the interface was a closed surface. Therefore, blood vessels have been modeled as a closed interface in the shape of a tube with capped ends arthurs1998modeling, li2014computing, li2013hybrid, rosar2001fluid, smith2003advective, rosar1994three. The closed tube was immersed inside a rectangular fluid domain. The flow in the closed tube is driven by adding a fluid source and sink to opposing ends of the tube. These modifications create unrealistic flow in biological models of blood vessels, especially near the source and sink. Typically, a large computational domain is used to compute solutions, but only the results in the center of the tube away from the source and sink are considered. We have created a novel extension of the immersed interface method to interfaces that are not closed, but instead, are shaped like an open tube that spans from one end of the fluid domain to the other. This method will be referred to as the immersed interface method for open tubes (IIM-OT). By using the IIM-OT, we can create a more natural fluid profile using a smaller computational domain.

Section 2 details the problem formation. In particular, the immersed interface method requires the magnitudes of the discontinuities or jump conditions of the primitive variables, and their derivatives are determined a priori. We show that the jump conditions for the immersed interface method for both closed and open tube-shaped interface can be computed in the same manner from the force strength function. Section 3 describes the numerical method. Section 4 includes numerical simulations to show that the novel IIM-OT for the Navier-Stokes equations and the Stokes equations achieved second-order accuracy in space for both 2D simulations in rectangular coordinates and 3D simulations in axisymmetric cylindrical coordinates.

2 Problem formulation

2.1 Computational domain and immersed interface

We formulate a model that simulates fluid flow through an open tube with compliant walls, extending from one boundary of the computational domain to the opposite boundary. We formulate the 3D model in cylindrical coordinates in the fluid domain which contains an immersed interface in the shape of an infinitely-thin, compliant tube with a nonuniform diameter that extends from one end of the domain to the opposing end, shown in Fig. LABEL:fig:cyl3Ddomain. Additionally, we assume that the flow is axisymmetric, thereby allowing us to reduce the computational costs from 3D to 2D. The 2D computational domain is a slice of the 3D domain where or , shown in Fig. LABEL:fig:cyl2Ddomain. The immersed interface surface can be represented as a 2D curve, which is also be called since the structure of interest is clear from the context. Assume intersects the computational domain boundary at and . Let

be the unit normal vector oriented towards the outside of the tube.


a b
Figure 1: (a) The 3D fluid domain is radially symmetric and represented in cylindrical coordinates . The open interface is a surface in the shape of an irregular tube that spans . (b) The 2D computational domain is a slice of the 3D domain where or . The computational slice contains an open immersed interface curve . Note intersects the computational domain boundary at and . is the unit normal oriented towards the outside of the tube.

2.2 Boundary condition

Appropriate boundary conditions must be chosen to drive the fluid through the tube, which can be done by creating an axial pressure gradient inside the vessel. One way to accomplish this is to specify the value of pressure on the inlet and outlet of the tube using Dirichlet boundary conditions.

At the inlet, it is assumed that has a parabolic profile consistent with 3D Poiseuille flow. In 3D axisymmetric cylindrical coordinates, the pressure gradient and velocity are related by

(1)
(2)

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. The boundary conditions for the full solution are non-smooth for and discontinuous for at the inlet of the tube. The boundary conditions for the full solution are imposed on the Stokes part since is non-smooth and is discontinuous near the boundary. Homogeneous boundary conditions of the same type as the full solution are imposed for the regular solution. Below are the boundary conditions used for the two simulations in Section LABEL:sec:NSCylResults where the average inlet velocity is denoted .

  • Inlet Boundary ()

    (3)
    (4)
    (5)
    (6)
  • Outlet Boundary ()

    (7)
    (8)
    (9)
  • Top Boundary ()

    (10)
    (11)
    (12)
  • Bottom Boundary ()

    (13)
    (14)
    (15)

2.3 Fluid structure interactions

The Navier-Stokes equations take the form

(16)
(17)

where , where , , and denote velocity components in the -, -, and -directions, respectively, is the pressure, and is viscosity; and is the interfacial force, which is singularly supported along (see below). We connsider also the zero Reynolds number regimes given by the Stokes equations:

(18)

and the continuity equation (Eq. 17)

The force exerted by the interface can be written as

(19)

where denotes the position of the interface, in 3D is the material coordinate(s) that parameterize the interface curve or surface at time , is the force strength at the point , and is the Dirac delta function. Since all 3D cases are axisymmetric, the interface is effectively a 2D curve, and it is sufficient to only consider .

The force strength from Eq. 19

(20)

is comprised of a two major components: an elastic force and a tether force . Since the interface has elastic properties, deviation from its resting configuration generates a restorative force. In 3D cylindrical coordinates, the elastic force derived in Ref.  lai2011simulating is given by

(21)

where is the mean curvature and the tension is given by

(22)

In Eq. 22, controls the stiffness of the interface. In Eq. 21, the unit tangent vector to is given by

(23)

We assume that at steady state with the tethers at their equilibrium positions, tubular flow is that of Poiseuille flow, with fluid outside of the channel at rest. Thus, at steady state, one expects a jump discontinuity in p and in

across . To generate this steady-state flow profile, we assignn to the tether forces two components:

(24)

where is the force component needed to support Poiseuille flow along the channel (by generating the necessary jumps in the solution and its derivative).

The second component arises from the displacement of the tethers from their equilibrium positions or anchor points. The interface control knots are tethered to anchor points in the fluid domain by a spring with resting length 0. Suppose the interface is anchored to . Let be the spring force constant. Then tether force is

(25)

If the interface knots move away from their anchor points, restorative forces are generated. Interface movement can either be restricted by using stationary anchor points or induced by moving the position of the anchor points in time.

The interface is deformable and is assumed to move at the same speed as the local fluid. The no-slip condition

(26)

describes this motion.

2.4 Boundary conditions

Appropriate boundary conditions must be chosen to drive the fluid through the tube, which can be done by creating an axial pressure gradient in the vessel. One way to accomplish this is to specify the value of pressure on the inlet and outlet of the tube using Dirichlet boundary conditions. Another way is to impose inhomogeneous periodic boundary conditions. Specific boundary conditions used for the 2D and 3D cases simulations are discussed before the results.

Bi-periodic boundary conditions are imposed for velocity, which implies that the volume of the fluid in the tube remains constant in time. To drive flow, we prescribe a pressure gradient inside the tube, by requiring there to be a constant difference in pressure at the inlet and outlet of the tube

(27)

where is the radius of the tube. The derivatives of pressure on the and boundaries are required to be equal

(28)

This is referred to as an inhomogeneous periodic boundary condition. Not only does this create a pressure gradient in the tube, but it also forces the pressure gradient to be periodic, which is needed since periodic boundary conditions are imposed for velocity. Periodic boundary conditions are imposed for pressure at the boundaries.

2.5 Derivation of jump conditions

The jump conditions for the Stokes equations and the Navier-Stokes equations are derived in Ref.  li2013hybrid and Ref. li2014computing respectively for closed interfaces. Both the Stokes equations and the Navier-Stokes equations have the same jump conditions given the same interfacial forces li2013hybrid, li2014computing. By creating a fictitious closed interface by adding segments to close the open tube interface, the same derivation for closed interfaces from Ref. li2014computing can be applied to open tube-shaped interfaces. This is shown in detail in Section 2.5 for the Stokes equations. The jump conditions for both the 3D Stokes and Navier-Stokes equations in cylindrical coordinates for both closed and open tube-shaped interfaces are

(29)
(30)
(31)
(32)

where is the angle between the normal and direction li2013hybrid, li2014computing.

Next, we derive the jump conditions across an open tube for the 3D Stokes equations in axisymmetric cylindrical coordinates. We show that the jump conditions in cylindrical coordinates for an interface that is in the shape of an open tube have the same dependence on force and interface position as in the closed tube case.

Again, the interface must be shaped like a closed surface for the derivation. Therefore, consider the fictitious closed surface

(33)

which is formed by the interface and the inlet and outlet of the tube which intersect the fluid domain at and respectively.

Let be an extended domain that contains the surface interface where the distance between and shrinks to zero as approaches zero. Note that extends outside of the original fluid domain .

Although symmetry is imposed later, first consider the fully 3D Stokes equations in cylindrical coordinates, given by

(34)
(35)
(36)
(37)

where , , , and are the velocities in the , , and directions respectively, and with

(38)

for , where is the extension of which is supported on and agrees with on . Even if the correction terms are computed for the fictitious closed domain there is no contribution from the extensions of the interface. Therefore, the hat symbol is dropped for the rest of the derivation. Now the derivation mirrors the procedure found in Ref. li2013hybrid.

Taking the divergence of Eqs. 34, 35, and 36 gives Poisson’s equation for pressure

(39)

Let be an arbitrary twice continuously differentiable test function. Multiply the left hand side of Eq. 39 by and integrate over

(40)
(41)
(42)

The first and second lines are by Green’s first Identity and the divergence theorem, respectively. The last line is found by taking the limit as goes to zero. The last term goes to zero since pressure is bounded. The surface element is . By imposing the axisymmetric assumptions on and we get that

(43)

where denotes the curve formed by restricting the surface to a fixed . Multiply the right hand side of Eq. 39 by and integrate over .

(44)
(45)
(46)

The last line is found using the divergence theorem and noting that the has support on and is zero on . Applying the axisymmetric assumption that and are independent of to Eq. 46 yields

(47)
(48)
(49)

where denotes the curve formed by restricting to a fixed . The second line is found by breaking the force into components normal and tangential to the interface. The last line is found using integration by parts. Since is arbitrary, we have

(50)
(51)

To find the jump conditions for , multiply Eq. 34 by and integrate

(52)

The first term on the right hand side of the Eq. 52 approaches the following as

(53)
(54)

The second term on the right hand side of the Eq. 52

(55)

approaches 0 as . This is true because by the axisymmetric condition when , approaches zero. Additionally , which is finite.

The first term on the left hand side of the Eq. 52 approaches the following as .

(56)
(57)

where is the angle between the normal and direction. Therefore

(58)

Similarly,

(59)

3 Numerical method

To solve the immersed boundary problem, we compute the fluid velocity and pressure on a fixed Eulerian grid. A moving Lagrangian frame of reference is used to track the location of the interface over time. In 2D, the Stokes equations are solved on the fluid computational grid,

where is the grid spacing and and are the number of grid points in the and directions respectively. In 3D, the Stokes equations are solved on the fluid computational grid,

where is the grid spacing and and are the number of grid points in the and directions respectively. In both the 2D and 3D case, the interface position at time is tracked by boundary markers that are connected by a cubic spline.

3.1 Velocity decomposition method

The velocity decomposition method, developed by Ref.  beale2009velocity, can be used to avoid computing most of the correction terms that are needed for the immersed interface method for the Navier-Stokes equations. This method leverage the fact that, given the same singular interfacial force, the jumps in the solutions are identical for both the Stokes equations and the Navier-Stokes equations as shown by Ref.  leveque1997immersed and Ref.  li2001immersed respectively. This result occurs because the no-slip condition implies that velocity and its material derivative are continuous at the interface lai2001remark.

It is sufficient to break the solution of the Navier-Stokes equations into a part that satisfies the Stokes equations with the singular interfacial forces and a regular part that solves the remaining equations as shown below.

Suppose and solve the Navier-Stokes equations with a singularly-supported interfacial force . Let

(60)
(61)

where the subscripts and are used for the Stokes and regular part, respectively. The Stokes part, and , satisfies the Stokes equations with the interfacial force

(62)
(63)

The regular part, and , satisfies the remaining terms in the full Navier-Stokes equations

(64)
(65)

where

(66)

The Stokes part is found first, so is known and treated as a body force when solving for the regular part. The force is a continuous function in the fluid domain since and its material derivative are continuous. Notice that the full velocity is used in the transport of and in .

This method has several advantages over the immersed interface method for the Navier-Stokes equations. First, correction terms only need to be computed for the Stokes part of the solution, saving considerable time in implementation since numerous correction terms are needed for the additional terms in the Navier-Stokes equations. Secondly, the velocity decomposition method can also be implemented with fractional time steps that help overcome some of the challenges of simulating stiff interfaces. Small fractional time steps are used to move the interface, which is done by only computing the velocity on the interface. This computation is done using Stokeslets to compute the Stokes part of the velocity at the interface and interpolating the regular part of velocity on the interface in time

beale2009velocity. Therefore, the interface position can be updated without computing the full solution at each fractional time step.

3.2 Stokes part: immersed interface method for open tube interfaces

The Stokes part of the equation is solved using the IIM-OT. This method requires the jump conditions for the solution to be known a priori. The velocity decomposition method is based on the fact that the jump conditions for the Navier-Stokes equations and the Stokes equations with the same interfacial forces are identical beale2009velocity, which is true because the no-slip condition implies that

(67)

since the interface moves at the velocity of the local fluid lai2001remark. Taking the material derivative of velocity gives

(68)

Therefore, the advection term does not contribute to the jump conditions lai2001remark.

Since the jump conditions are known, the immersed interface method for open tube interfaces needs to only correct the incompressible Stokes equations. First, the incompressible Stokes equations are converted into a series of uncoupled Poisson problems. The Poisson problem for pressure, is found by taking the divergence of Eq. 62 and use the incompressibility condition

(69)
(70)

Once is known, Eq. 62 reduces to a Poisson problem for each component of velocity, . Therefore, to solve the Stokes part it is sufficient to use the immersed interface method for Poisson problems.

The immersed interface method for Poisson problems was developed by Ref.  leveque1994immersed. The Poisson problems are discretized using the finite difference method. Therefore, to apply the IIM-OT to the Stokes equations, it is only necessary to compute the corrections terms for the first- and second-order spatial derivatives. Once the jump conditions for the open tube interface are known, the correction terms can be computed using the immersed interface method for Poisson problems.

3.3 Regular part: projection method and semi-Lagrangian method

The regular part of the solution satisfies Eqs. 6466. The advection terms and the body force are discretized using the backward difference formula, where the upstream values of velocity are found using the semi-Lagrangian method. A projection method is used to enforce the incompressibility condition.

3.3.1 Semi-Lagrangian method

The advection terms, and , can be thought of as material derivatives or the derivatives along the path of a fluid particle. In the semi-Lagrangian method, the material derivative is computed at Eulerian grid points, which combines the Lagrangian and Eulerian variable frameworks and has several benefits over a purely Lagrangian or Eulerian framework.

The material derivative is computed at Eulerian grid points. Therefore, the resolution of the grid does not deteriorate or get distorted with the fluid flow, which can occur in Lagrangian methods.

Particles cannot physically cross the impermeable immersed interface so that the material derivatives will be smooth, which reduces the number of correction terms needed in the immersed interface method. For example, correction terms would be needed for the time derivatives computed on an Eulerian grid when the interface crosses over a grid point in a time step. Additionally, is discontinuous at , so the terms, and , would need to be corrected near the interface if discretized on the Eulerian grid layton2003semi-Lagrangian, xiu2001semi. Additionally, the nonlinear terms of the Navier-Stokes equations do not need to be discretized explicitly in the semi-Lagrangian method.

The semi-Lagrangian method is used to compute the velocity along the path of a fluid particle that passes through the grid point at time , which is used to compute the material derivative of the velocity of the particles. This is done by first computing the upstream positions of the particle, and , at two previous times, and , respectively. An example of a 1D particle path is shown in Fig. LABEL:fig:particlePath.

The upstream position of the particle located on the Eulerian fluid grid point, , at time can be found by integrating

(71)

backwards in time. The predictor-corrector method

(72)
(73)

and

(74)
(75)

can be used to find the upstream positions at two previous time steps with second-order temporal accuracy. Notice that is approximated using the time extrapolation .

It is unlikely that the upstream positions coincide with grid points. Therefore the velocity at upstream values

(76)
(77)

are interpolated in space from the velocity at grid points using cubic Lagrangian interpolation. This interpolation has 4th order accuracy in space. Once the upstream values of velocity are known, the material derivative is found using the second-order backwards difference formula. Therefore, Eq. 64 can be discretized using

(78)

3.3.2 Projection method

The coupling of the velocity components and pressure in the Navier-Stokes equations would create a large finite difference system if all components were solved simultaneously. Instead, the projection method is used to solve the regular part. In this procedure, an intermediate value of velocity and pressure are found that satisfy

(79)

The intermediate values are then projected, , into the subspace of divergence-free vector fields so that This works because every vector field can be written as a Hodge decomposition where is divergence-free. This projected value satisfies the Navier-Stokes equations.

The regular part of the Navier-Stokes equations can be numerically solved using the following projection method. An intermediate velocity value, is found by solving

(80)

where is the value of pressure from the previous time step. The intermediate value is approximately projected into divergence-free space by defining as

(81)

which can be found by solving

(82)

Then, velocity and the gradient of pressure can be found by

(83)
(84)

3.4 Time-stepping

The interface movement is govern by Eq. LABEL:eq:noSlip. Suppose the interface position at time is represented by boundary markers . The velocity at each marker can be found by interpolating the velocity from the fluid grid to the markers using a second-order method such as bi-linear interpolation. Care must be taken during this interpolation since the velocity near the interface is not smooth. The velocity could be corrected near the interface or a one-sided approximation for velocity could be used where all interpolation points are on the same side of the interface.

For the time stepping to be second-order, the velocity on the interface at two previous time levels, and , are needed. The interface position is then updated as follows:

(85)

To move the interface initially, forward Euler’s method was used.

3.5 Velocity decomposition method algorithm outline

At time , the following are known: , , , , and for and . The following procedure can be used to compute the solution at .

  1. Update interface position .

    1. Interpolate velocity at the interface from fluid grid using

      (86)
    2. The interface positions is updated using the no slip condition

      (87)
    3. Compute the interfacial force from the interface position .

  2. Compute the Stokes solution, and .

    1. Use force to find the correction terms, and , for pressure and velocity respectively using the immersed interface method.

    2. Solve the following Poisson problem for pressure, :

    3. Compute the discrete pressure gradient, using correction terms.

    4. Solve the following Poisson problems for each component of velocity, .

  3. Compute regular solution, and .

    1. Use semi-Lagrangian method to compute , , and

    2. Solve momentum equation for intermediate regular velocity, .

    3. Compute needed to project into divergent-free space.

    4. Find and using .

  4. Compute the full solution at time by adding the Stokes and regular part.

Forward Euler’s method is used for the first step of the time stepping methods. The local truncation error is and this is only used for one time step, second-order temporal accuracy is still maintained.

4 Numerical results

4.1 Stokes Results

A grid refinement study was conducted to test the spatial convergence rate for the immersed interface method in an open tube for the Stokes equations in 3D axisymmetric cylindrical coordinates shown in Eqs. 8890. In this test, the tube was displaced and allowed to return to its resting position as a tube with a constant radius. In addition to confirming the order of the method, this test also aims to validate that approximate 3D Poiseuille flow is achieved when the tube has near constant radius.

(88)
(89)
(90)

For the cylindrical coordinate test cases, the domain is assumed to be radially symmetric. This allows three dimensional kinetics with two dimensional computational costs. Although, the fluid domain is , the solution is only computed on the computational domain which is the slice of the fluid domain. cm is the resting radius of the tube.

For this case, the interface was initially displaced in the direction according to the function

(91)

where is the resting radius of the tube. The interface was then allowed to return to its resting position as a straight tube.

Flow is driven through the tube by prescribing a pressure gradient over the computational domain using Dirichlet boundary conditions. At the inlet, it is assumed that has a parabolic profile consistent with 3D Poiseuille flow. In 3D axisymmetric cylindrical coordinates, the pressure gradient and velocity are related by

(92)
(93)
  • Inlet Boundary ()

  • Outlet Boundary ()

  • Top Boundary ()

  • Bottom Boundary ()

The remainder of the parameters used for this experiment are shown in Table 1.

Parameter Symbol Value
Viscosity 0.1 gm/(cms)
Density 1 gm/cm
Domain length 4.25 cm
Domain height 1.625 cm
Initial tube radius 0.7 cm
Interface control points 100
Simulation length 4 s
Time steps 2000
Time step size 2e-3 s
Tether force constant 10 gm/s
Elastic force constant .1 gm/s
Average Inlet Velocity 1 cm/s
Table 1: Parameters used for the spatial convergence study of the IIM-OT for 3D Stokes flow in axisymmetric cylindrical coordinates.

The solution was computed at various resolutions and then compared to a high-resolution solution. The results of the convergence study shown in Table 2 indicate that the method is second-order accurate in space.

Grid Size p w u
Order Order Order
35 14 0.03741 0.25162 0.10789
69 27 0.00796 2.23 0.10349 1.28 0.02643 2.03
137 53 0.00189 2.07 0.01862 2.47 0.00628 2.07
273 105 0.00044 2.12 0.00456 2.03 0.00145 2.12
545 209 0.00008 2.37 0.00027 4.07 0.00029 2.34
Table 2: Spatial convergence results of the IIM-OT for Stokes flow in cylindrical coordinates. Solutions are compared to a high resolution, solution at time .

Fig. 3 shows the interface position overlaid on a quiver plot that show the velocity as an arrow with components various times. Since the most dramatic changes occur in the beginning of the simulation, the Fig. 1(a), 1(b), and 2(a) show the results at , , and respectively. Recall that because of the no-slip condition on the interface, the interface moves at the local fluid velocity. Therefore, the velocity slope field in Fig. 1(a) indicates that the interface is moving towards its resting position as a straight tube (shown in Fig. 2(b)). Initially in Fig. 1(a), some of the fluid is pushed to the left towards the proximal end of the tube as the interface relaxes. As the interface gets closer to its resting configuration, it exerts less force on the fluid. Therefore, in Figs. 1(b)2(b), all fluid in the tube moves to the right towards the distal end of the tube. Fig. 2(b) shows the resting position of the interface, in which the velocity exhibits a parabolic profile inside the interface and dissipates to zero outside of the interface.

((a)) The interface position and velocity quiver plot at s shows the interface is moving towards cm.
((b)) The interface position and velocity quiver plot at s shows the interface is moving towards cm.
Figure 2: The interface (black line) is shown at various times as it relaxes from its displaced in the axisymmetric cylindrical computational domain. The interface position is superimposed on a quiver plot that show the velocity as an arrow with components .
((a)) The interface position and velocity quiver plot at s shows the interface is moving towards cm.s
((b)) The interface position and velocity quiver plot at s shows the interface at its resting position where cm.
Figure 3: The interface (black line) is shown at various times as it relaxes from its displaced in the axisymmetric cylindrical computational domain. The interface position is superimposed on a quiver plot that show the velocity as an arrow with components .
((a)) The initial pressure in the computational domain at s.
((b)) The pressure at s as the interface is returning to its resting position.
Figure 4: Pressure in the axisymmetric cylindrical computational domain at various times as the interface relaxes from its displaced position to rest.
((a)) The pressure at s as the interface is returning to its resting position.
((b)) The pressure is linearly decreasing inside the tube and zero outside at s when the interface is in its resting position.
Figure 5: Pressure in the axisymmetric cylindrical computational domain at various times as the interface relaxes from its displaced position to rest.

The pressure in the computational domain at s, s, s and s are shown in Fig. 5. The IIM-OT captures the discontinuity in pressure across the interface as can be seen in Fig. 5. Additionally, pressure decreases along the length of the domain inside of the interface, which creates a pressure gradient that drives fluid through the open tube. Another pressure gradient outside of the tube can be seen in Fig. 3(a) which drives the fluid to restore the interface to its resting position. Finally, Fig. 4(b) shows a linearly decreasing pressure that is typical of Poiseuille flow in a cylinder.

4.2 Navier-Stokes Equations

4.2.1 Initial solutions

For the velocity decomposition method, an initial solution is needed for both the full solution and the regular part. In the two simulations given in Section LABEL:sec:NSCylResults, the interface is initially a straight tube with a constant radius. Therefore, the initial solution of the full solution is taken to be approximately Poiseuille flow. Since the interfacial force is in the Stokes part, the Stokes part will approximate Poiseuille flow. Therefore, the Stokes part and the full solution will have identical initial solutions. The initial solution of the regular part will be identically zero. The initial solutions are

(94)
(95)
(96)
(97)

where is the average inlet velocity.

4.2.2 Spatial convergence study: tethers with sinusoidal movement

A grid refinement study was conducted to test the spatial convergence rate for the IIM-OT for the Navier-Stokes equations in 3D axisymmetric cylindrical coordinates. In this simulation, the initially straight interface is actively moved by displacing the tether anchor positions in time. The -position of the tether anchor at position at time is given by

(98)

where is the initial interface radius, the period is and the amplitude is m. The positions of the interface tethers at ms are shown in Fig. 6.

Figure 6: Interface tether position at ms.

Even though the fluid domain is , the solution is only computed on the computational domain which is the slice of the fluid domain. The initial radius of the tube is m.

The boundary conditions and initial solutions for the full solution and the Stokes and regular parts are given in Section LABEL:subsec:NScylBC. The parameters for this simulation are shown in Table 3.

Parameter Symbol Value
Viscosity 0.0175 gm/(cms)
Density 1.055 gm/cm
Domain length 64 m
Domain height 32 m
Initial tube radius 14.14 m
Simulation length 1e-3 s
Number of time steps 16
Time step size 6.25e-5 s
Tether force constant