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 firstorder 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 NavierStokes equations ^{li2001immersed, lee2003immersed}. The immersed interface method achieves secondorder accuracy in numerical approximation. For fluidstructure 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 secondorder 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 (IIMOT). By using the IIMOT, 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 tubeshaped 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 IIMOT for the NavierStokes equations and the Stokes equations achieved secondorder 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 infinitelythin, 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 
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 nonsmooth 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 nonsmooth 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 NavierStokes 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 steadystate 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 noslip 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.
Biperiodic 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 NavierStokes equations are derived in Ref. ^{li2013hybrid} and Ref. ^{li2014computing} respectively for closed interfaces. Both the Stokes equations and the NavierStokes 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 tubeshaped interfaces. This is shown in detail in Section 2.5 for the Stokes equations. The jump conditions for both the 3D Stokes and NavierStokes equations in cylindrical coordinates for both closed and open tubeshaped 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 NavierStokes 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 NavierStokes equations as shown by Ref. ^{leveque1997immersed} and Ref. ^{li2001immersed} respectively. This result occurs because the noslip condition implies that velocity and its material derivative are continuous at the interface ^{lai2001remark}.
It is sufficient to break the solution of the NavierStokes 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 NavierStokes equations with a singularlysupported 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 NavierStokes 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 NavierStokes 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 NavierStokes 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 IIMOT. 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 NavierStokes equations and the Stokes equations with the same interfacial forces are identical ^{beale2009velocity}, which is true because the noslip 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 IIMOT to the Stokes equations, it is only necessary to compute the corrections terms for the first and secondorder 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 semiLagrangian method
The regular part of the solution satisfies Eqs. 64–66. The advection terms and the body force are discretized using the backward difference formula, where the upstream values of velocity are found using the semiLagrangian method. A projection method is used to enforce the incompressibility condition.
3.3.1 SemiLagrangian method
The advection terms, and , can be thought of as material derivatives or the derivatives along the path of a fluid particle. In the semiLagrangian 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 ^{layton2003semiLagrangian, xiu2001semi}. Additionally, the nonlinear terms of the NavierStokes equations do not need to be discretized explicitly in the semiLagrangian method.
The semiLagrangian 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 predictorcorrector method
(72)  
(73) 
and
(74)  
(75) 
can be used to find the upstream positions at two previous time steps with secondorder 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 secondorder 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 NavierStokes 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 divergencefree vector fields so that This works because every vector field can be written as a Hodge decomposition where is divergencefree. This projected value satisfies the NavierStokes equations.
The regular part of the NavierStokes 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 divergencefree 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 Timestepping
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 secondorder method such as bilinear 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 onesided approximation for velocity could be used where all interpolation points are on the same side of the interface.
For the time stepping to be secondorder, 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 .

Update interface position .

Interpolate velocity at the interface from fluid grid using
(86) 
The interface positions is updated using the no slip condition
(87) 
Compute the interfacial force from the interface position .


Compute the Stokes solution, and .

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

Solve the following Poisson problem for pressure, :

Compute the discrete pressure gradient, using correction terms.

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


Compute regular solution, and .

Use semiLagrangian method to compute , , and

Solve momentum equation for intermediate regular velocity, .

Compute needed to project into divergentfree space.

Find and using .


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, secondorder 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. 88 – 90. 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  2e3 s  
Tether force constant  10 gm/s  
Elastic force constant  .1 gm/s  
Average Inlet Velocity  1 cm/s 
The solution was computed at various resolutions and then compared to a highresolution solution. The results of the convergence study shown in Table 2 indicate that the method is secondorder 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 
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 noslip 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.
The pressure in the computational domain at s, s, s and s are shown in Fig. 5. The IIMOT 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 NavierStokes 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 IIMOT for the NavierStokes 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.
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  1e3 s  
Number of time steps  16  
Time step size  6.25e5 s  
Tether force constant 
Comments
There are no comments yet.