Insect-Inspired Micropump : Flow in a Tube with Local Contractions

A biologically-inspired micropumping model in a three-dimensional tube subjected to localized wall constrictions is given in this article. The present study extends our previous pumping model where a 3D channel with a square cross-section is considered. The proposed pumping approach herein applies to tubular geometries and is given to mimic an insect respiration mode, where the tracheal tube rhythmic wall contractions are used/hypothesized to enhance the internal flow transport within the entire respiration network. The method of regularized Stokeslets-mesh-free computations is used to reconstruct the flow motions induced by the wall movements and to calculate the time-averaged net flow rate. The time-averaged net flow rates from both the tube and channel models are compared. Results have shown that an inelastic tube with at least two contractions forced to move with a specific time lag protocol can work as a micropump. The system is simple and expected to be useful in many biomedical applications.


Introduction
Bioinspired micromachines are considered to be one of the most novel directions toward miniaturizing equipment for several biomedical applications.This will be made possible by developing state-of-the-art microfabrication technologies that can help in building micro-and nano-scale structures.These microstructures can be then used as components in functional micromachines, such as pumps, valves and separators, where all of these components can be built and assembled on a single chip.Here, we focus our attention on a novel micropumping paradigm that is inspired by insect respiration physiology [1][2][3][4][5][6][7][8][9].
The flow transport within insect physiological systems is complex and non-stationary because of insects' use of body movements and rhythmic wall contractions along their tracheal tubes to enhance internal flow motions.According to the author's knowledge, to date, there is no detailed numerical simulation study given toward understanding or resolving insect internal flow motions.However, other multiple theoretical and numerical research efforts from the literature can actually be useful when modeling these bio-fluid problems.For example, a simplified solution that can describe the unsteady viscous flow motions in a semi-infinite pipe with wall contractions or expansion modes is given in [10].Another solution to the incompressible-viscous flow in an infinitely-long channel with prescribed pulsatile and sinusoidal wall motions is derived in [11].Furthermore, the unsteady squeezing of a viscous fluid from a shrinking or expanding tube is given analytically in [12][13][14].The flow in a tube due to wall contractions is also given by multiple studies [15][16][17].
The internal flow motions inside many physiological systems are normally driven by wave propagation of the muscular wall contractions, known as peristaltic pumping.Peristalsis-induced flow motions have been a main research area in classical biological fluid dynamics.There are several theoretical and numerical methods that have been used to study peristalsis-induced flow transport in two-dimensional and three-dimensional tubes [18].
In many of physiological systems that use wall contractions to transport blood or air, the kinematics of these wall motions have been shown to greatly impact the overall flow transport efficiency.Generally, there are two main types of localized wall contractions, namely stationary and progressive (local peristalsis) contractions can naturally exist in these systems, as reviewed in [19].For example, both types have been observed in the human duodenum.These contractions have been shown to produce and enhance fluid transport in the human intestine, as well.An important conclusion from this review [19] and supported by a basic theoretical analysis is that at least two stationary contractions operating with a slight time lag (phase lag) are needed in order to produce net flow transport [20].Additionally, insects' tracheal tubes are another example of physiological systems that use localized wall contractions as a natural pump to deliver oxygen to the entire proximal body cells [4][5][6][7].
In this study, we extend our bioinspired micropumping models [1][2][3][4][5][6][7][8] to simulate the flow transport in a tubular geometry.In particular, we present a three-dimensional numerical simulation based on the regularized Stokeslets method of fundamental solutions (MFS) to study flow motions within an insect tracheal-like tube.Results from this study are given to complement and to strengthen our previous pumping finding, which is inspired by the rhythmic contractions found in insect tracheal tubes.

Methods
In this section, we use the three-dimensional Stokeslets-mesh-free computational model based on the method of regularized Stokeslets [21] to study the contraction-induced flow pumping of viscous fluid in a tube with rhythmic contractions.Results from this analysis will be compared and validated against our 2D theoretical analysis derived previously [6].The flow structures and development in a tube with two wall contractions governed by a generic wall profile R(z, t) will be given in detail.The proposed wall profile permits contractions to move with various phases (time) lags with respect to each other.The details of the computational approach used herein can be found in [18,[21][22][23] and are summarized below.
The method of regularized Stokeslet [21] is based on finding exact solutions of the Stokes equations using body forces represented by smooth and localized elements, which satisfy the Stokes incompressible flow equations: where µ is the fluid viscosity, p * is the pressure, V * is the velocity and F * is the force per unit volume, which can be written in terms of Green's function and regularization (cut-off) expression as, where φ * (x * − x * 0 ) is a cut-off function with a property that φ * (x * ) dx = 1.In our computations, we use a specific cut-off function that is given as [18], where r * = x * − x * 0 .* is a small regularization parameter that is given to control the spread of the cut-off function.The velocity field and pressure solutions to the above governing equations can be given as: The indices i, j = 1, 2, 3 follow the Einstein summation convention.Furthermore, V * i = (w * , v * , u * ) is the flow velocity vector, and p * is the static pressure.x * i = (z * , y * , x * ) is the position vector.α * = (α * z , α * y , α * x ) is the force vector strength along z * , y * , x * directions.S ij is known as a regularized Stokeslet, and P j is the stress tensor, which can be given as: where For clarity, the velocity expression using the above equations can be easily rewritten in component form as follows, Similarly, the pressure is given as: Now, we consider the same non-dimensional parameters defined in [6]: The above velocity expressions in component form can be rewritten in a non-dimensional form as, Similarly, the non-dimensional pressure is given by: where

Results and Discussion
The MFS-Stokeslets mesh-free technique is used herein to obtain the velocity field and pressure within the tube pump under consideration.The solution is formed by summing all of the Stokeslets expressions with unknown intensities.The strengths of these Stokeslets singularities can be obtained by enforcing and satisfying the prescribed boundary conditions using direct collocation approaches [24] and [21].For the problem under consideration herein and shown in Figure 1, N Stokeslets (regularized source points) with unknown strengths α are distributed along the tube boundaries.As a result, the induced flow field by these singular points is then approximated by applying the principle of superposition and taking the effect from all Stokeslets points.The velocity field can be then evaluated as: Similarly, the pressure is given by: where 5/2 and α zj are the unknown force coefficients that represent the strengths of the Stokeslets in x, y and z-directions, respectively.x i = (x i , y i , z i ) is the position of any field point.x 0j = (x 0j , y 0j , z 0j ) are the locations of the Stokeslets source points.r ij = |x i − x 0j | is the distance between any field point and another source point.
In order to determine the unknown strengths α = (α xj , α yj , α zj ), the boundary conditions for velocity components and pressure are prescribed and collocated at the boundary field points, as shown in Figure 1a.Adding each fundamental solution using Equations ( 17)-( 20), the final solution can be obtained by solving a system of equations: where A is a matrix of size 3N × 3N with real entries formed by evaluating the right-hand side of the above expressions in Equations ( 17)- (20).b is a vector of size 1 × 3N of real entries formed by evaluating the left-hand sides of the same equations.In other words, all of the entries in both A and b are determined by enforcing the fundamental solutions to satisfy the boundary conditions by direct collocation.α represents the Stokeslets force intensities, which are to be determined by solving the above system of equations.
It should be mentioned here that, although the MFS is relatively easy to implement and does not require any meshing, it suffers from the ill-conditioning issue when finding the force strengths [21,[24][25][26].This method's drawback can lead to inaccurate Stokeslets strength calculations and, thus, inaccuracies in finding the induced flow field.A solution to the above linear system problem can be found by either collocating the Stokeslets sources at an optimal distance away from the domain boundaries as proposed by [24] or by using regularized Stokeslets expressions as derived by [21].Alternatively, similar results are accomplished by introducing a regularization to the matrix A itself [1]; this approach has been shown to provide an efficient and fast converging of the features when solving the linear system.Figure 1.Problem schematic given to mimic the insect's main tracheal tube segment with two contractions [1] and the Stokeslets-mesh-free numerical setup: (a) 3D tube with moving upper wall contraction profile R(z, t); (b) g 1 (t) and g 2 (t), the motion protocols assigned to the first and second contractions, respectively.
In this paper, we follow the algorithm given in [1] to solve for the Stokeslets strengths and to reconstruct the flow field in a 3D tube subjected to localized wall contractions that move according to a specific wall profile, as shown in Figure 1b.More details about implementing the Stokeslets-mesh-free method can be found in [2].This section is organized as follows: Firstly, the wall profile is prescribed in full.The algorithm that will be used to overcome the ill-conditioning linear algebra issue during the calculation of the Stokeslets strengths will be given next.Thirdly, the structure and the development of the flow field induced by only two walls' contractions that move with a non-zero phase lag will be given.Finally, the net flow rate and pumping effect will be given.The presented results are compared to our 2D theoretical solution for validation purposes [6].The details of the subsections are given below.

Wall Profile, R(z,t)
The tube wall profile R(z, t) and its contracting kinematics are prescribed using the following mathematical expressions, where f i (z) ∈ C r [0, 1] and g i (t) ∈ C r [0, T = 1/S t ] represent the spatial and the temporal distributions of the tube surface wall shape, respectively.N c defines the number of contractions, and A i is the amplitude assigned to each contraction.The spatial term of the above equation imitates the geometry of the wall contractions as: where α = 2π/δ, z i defines the beginning of each collapse region and d i ∈ (0, 1 − z i ] marks its end, as shown in Figure 1b.In this work, we only consider two contractions, i.e., N c = 2, where the first contraction moves in time according to the following: profile: while the second contraction moves according to: where the non-dimensional parameter β relates the phase lag between the first and second contractions. Although, during the compression regime, both contracting sites start to move at the same time, the first contraction reaches the maximum specified tube travel collapse (TC) distance faster than the second contraction.In other words, there will be a time lag, T g = (1−1/β)/(2S t ), between these two collapsing motions, which is equivalent to a phase lag θ 12 = π(1 − 1/β).Moreover, during the expansion phase, the first contraction returns back to the nominal (original) position and continues its period with zero amplitude until the second contraction completes its own cycle; then, both contractions start the second cycle together with the protocol shown in Figure 1b.It should be noted that both contractions have a similar time period T = 1/S t , i.e., S t = 1 and that is, if β = 1, there will be no phase lag and both contractions will be identical.Here, we used z 1 = 0.25, z 2 = 0.65, 3528, which render the maximum tube collapse ratio to be T C = 70%.This motion protocol has chosen, such that, after one cycle, the tube geometry will be returning to the initial position, and there will be no net flow due to any volume deformation/displacements.This proposed tube wall profile plays an important role in breaking the flow symmetry and producing a unidirectional net flow.

Algorithm for Finding Accurate Stokeslets Strengths
Finding an accurate solution to the ill-conditioned system given by Equation ( 21) is an important step and requires a careful implementation.The regularized parameter = 0.15 and the total number of source points of N = 128×64×64 have been optimized to resolve the flow details.In other words, a total number of Stokeslets N = 2 19 is obtained by finding the (N-independent) solution that indeed satisfies the wall boundary conditions.This numerical algorithm follows the Tikhonov-regularization technique given in [27] to solve an ill-conditioned linear system of equations.The algorithm is based on introducing a regularization parameter, h = 0.01, into the original matrix, A, which is then transformed into a modified symmetric matrix, A R .This modified matrix acts as a pre-conditioner that leads to a modified system of equations with an enhanced condition number.The algorithm is simple and only requires a priori calculations of coefficients A, vector b, regularization parameter h = 0.01 and a convergence criteria based on an assigned specific tolerance value tol = 10 −6 .The steps of this algorithm are listed in Algorithm 1.
Algorithm 1 Algorithm for finding the strengths of each Stokeslet-source point.
As a direct implementation of the above algorithm, results for both velocity field and pressure induced by tube wall motions are given and compared to our 2D analytical solution [6].More specifically, a tube with two walls collapses that move with a phase lag of θ 12 = 30°is considered.Flow field snapshots at compression time t = t/4 and expansion time t = 3T /4 are explained and used to carry out these comparisons.

Tube with Two Local Rhythmic Wall Contractions
The method of fundamental solution MFS along with the regularized Stokeslets approach are used to solve for the induced flow field in a three-dimensional tube with two moving wall contractions.The contractions are assigned to the tube wall profile R(x, t) in a specific protocol to grantee net flow production.A schematic to the problem setting, along with the distribution pattern of the source points, is given in Figure 1a.The motion protocol used to actuate wall contractions is shown in Figure 1b.The goal of this section is to perform 3D simulations and to compare the results to our previously-derived 2D analytical model [6].In addition, the net flow rate is also computed and compared to the 2D scenario.In other words, results from both the 3D Stokeslets mesh-free computations and from the 2D theoretical model are compared at different snapshots.The comparison process is given at two snapshots of times t = T /4 and 3T /4 to represent compression and expansion collapsing phases, respectively.A phase lag θ 12 = 30 • is used to set the contraction time between the collapsing sites.
The velocity field and axial pressure are computed and shown in Figures 2-4 using contour plots.Results have shown that, when θ 12 = 30 • , the stagnation zone between contraction sites (observed when θ 12 = 0 • , [6]) is removed, and there exists flow transport within the region between the two contractions, as shown by the axial and vertical velocity contour lines calculated analytically and numerically, as shown in Figures 2 and 3.Moreover, there are pressure gradients in the region between contraction spots, which indicate that there is a flow transport in this region, as shown in Figure 4.

Unidirectional Net Flow
The above results suggest that when the wall contractions move with different time (phase) lags, the system can be used to create a unidirectional flow and be used as a micropump.The velocity fields at different phase lags and at each time elapse during the contraction cycle are then used to calculate the time-averaged net flow rate Q T as a function of the contraction phase lag θ 12 parameter.To calculate the time averaged net flow rate, the instantaneous flow rate at a desired location along the axial direction is calculated firstly.The details of how the time averaged net flow rate is calculated can be also found in [18], which can be re-given as follows: where A(z, t) is the cross-sectional area along the tube length z and time t.The above integration is approximated numerically using the standard trapezoidal rule.Once we have computed the instantaneous flow rate q(z, t), we can average it over one period of time T = 1/S t .The time-averaged net flow rate is simply computed by averaging computed values at discrete times, since the flow rate is a periodic function of time.The time-averaged net flow rate per unit length of the channel is computed by: where T = 1/S t is the time period.In the discrete form, the time-averaged net flow rate can be computed using the quadrature trapezoidal rule, where ∆t = T /N T , N T is the number of quadrature points in a time period T = 1/S t .In Figure 5, we show the time-averaged net flow rate as a function of the phase lag parameter.In the same figure, results obtained by the 2D analytical solution [6] are compared to the 3D Stokeslets-mesh-free computations for the tube pumping model.Results have shown that, although there is not good agreement between 3D tube and the 2D analytical model, the net flow rate using the 3D computations has a similar trend as the 2D analytical solution.Furthermore, the maximum value of the net flow occurs at θ 12 = 110 • compared to θ 12 = 65 • , which is predicted previously by the 2D analytical model.This phase lag discrepancy between the 2D and 3D results might be related to the three-dimensional effect, where the net flow is expected to depend on all spatial dimensions and not only on the z-direction, as assumed in the 2D counterpart.Although the net flow rate distribution as a function of phase lag of both 3D Stokeslets-mesh-free computations and 2D analytical results are different, they still share similar physical interpretations about our pumping hypotheses, which is inspired by insects' respiration process and its rhythmic tracheal wall collapses.

Channel vs. Tube Simulations
The present study is a natural extension to our previously-published pumping model in a 3D channel with a square cross-section [7].Results have shown that the flow field obtained by both the tube and channel is analogous and exhibits similar structures during the compression and the relaxation phases.This is not surprising, since both conduits have similar undeformed volume and the wall contractions are driven by the same wall profile protocol.The flow rate obtained by the tube simulation has been found to be ∼ 4-fold greater than the channel flow rate.However, the maximum flow rate for both the tube and channel occurs at almost the same phase lag angle, i.e., at θ 12 105 • -110 • .
The observed increase in the tube net flow rate can simply be explained as follows: in the channel case, we only constrict from the upper walls, while in the tube case, contractions occur along the entire circumference.This means that the total deformation (volume change) over a collapsing cycle for the tube case is greater than the channel setup.We expect that when both the tube and the channel undergo similar deformations over the collapsing cycle, they will produce a similar net flow rate.
Finally, it should be noted that the tube simulation is more relevant when mimicking the physiological behavior of insect tracheal tubes.Additionally, the tube results can serve as a guide to design novel tubular micropumps for several biomedical applications.For example, a vessel-like (tubular conduit) one with multiple contractions can be designed, actuated properly and be used as a heart assist device.This might be made possible by the current revolution in nanotechnology, tissue engineering and artificial muscles.

Conclusions
Inspired by the insect respiration processes and their internal flow motions, a pumping model in a three-dimensional tube subjected to two rhythmic wall contractions is proposed.The flow field and time-averaged net flow rate induced by these wall contractions are reconstructed using the regularized Stokeslets-mesh-free method.Numerical results are compared to our previously-derived two-dimensional analytical pumping model.The flow field at different time snapshots has been found to be in good agreement with the 2D solution.The net flow rate over a complete cycle shows a distinct difference.In particular, the maximum net flow rate occurs at a phase lag that is different from what has been obtained using the 2D analytical model.In spite of the observed differences in the net flow rate calculations between the 3D simulations and the 2D model, the proposed system is capable of producing unidirectional flows and can be used as a micropump.This simple pump is hypothesized to be useful in building a bioinspired set of micromachines for the next generation of biomedical applications.

Figure 5 .
Figure 5. Time-averaged net flow rate comparisons between the 2D analytical solution and 3D Stokeslets-mesh-free computations.