Accepted Manuscript Title: Simulation of the motion of two elastic membranes in Poiseuille shear flow via a combined immersed boundary-lattice Boltzmann method Author: As’ad Alizadeh Abdolrahman Dadvand PII: DOI: Reference:
S1877-7503(15)30042-9 http://dx.doi.org/doi:10.1016/j.jocs.2015.11.008 JOCS 435
To appear in: Received date: Revised date: Accepted date:
29-6-2015 17-11-2015 22-11-2015
Please cite this article as: A. Alizadeh, Simulation of the motion of two elastic membranes in Poiseuille shear flow via a combined immersed boundary-lattice Boltzmann method, Journal of Computational Science (2015), http://dx.doi.org/10.1016/j.jocs.2015.11.008 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Simulation of the motion of two elastic membranes in Poiseuille shear flow via a combined immersed boundarylattice Boltzmann method As’ad Alizadeh1, Abdolrahman Dadvand2 Department of Mechanical Engineering, Urmia University of Technology, Urmia, Iran
ip t
1,2
Corresponding author, Email:
[email protected], Tel: (+98) 441-3554180
1
Email:
[email protected]
us
cr
2
Abstract
Ac ce pt e
d
M
an
In this study, the motion, deformation and rotation of two elastic membranes in a viscous shear flow in a microchannel with and without a groove are simulated utilizing a combined LBM-IBM. The membranes are considered as immersed elastic boundaries in the fluid flow. The membranes are represented in Lagrangian coordinates, while the fluid flow field is discretized by a uniform and fixed Eulerian mesh. The interaction of the fluid and membranes is modeled using an appropriate form of the Dirac delta function. Two geometrically different channels, namely, a simple channel and a channel with a groove are considered. In the simple channel case, when the membranes are placed on the symmetry axis of the channel, they continue to move and deform without any lift force and rotation induced. However, when the membranes are located off the symmetry axis, the pressure difference produced in the flow around the membranes would apply lift forces on them and expel them towards the center of the channel. In the case of channel with a groove, it was found that when a membrane tensile modulus was reduced its flexibility as well as rotational speed would increase. This would, in turn, result in a pressure difference in the flow around the membrane. The pressure difference would apply a lift force on the membrane and move it out of the groove. It is worth mentioning that the results of the present simulation have a good agreement with the available numerical and experimental results. Keywords: Lattice Boltzmann method (LBM); Immersed boundary method (IBM); Elastic membrane; Interaction; Grooved channel
1. Introduction
Most natural, engineering, biological and medical problems include the interaction of elastic solids and viscous fluids [1]. The movement of an immersed object in a fluid including the particle deposition process, the drug delivery through the blood, the red blood cells shape change in the capillaries and arteries, etc. have been studied by many researchers. The difficulty of dealing with such problems would be due to complexity of the grid generation and application of boundary conditions involved. The numerical methods developed for solving such problems are divided into two main categories: the dynamic- (moving-) mesh approach and the fixed-mesh approach. In the
Page 1 of 34
Ac ce pt e
d
M
an
us
cr
ip t
dynamic mesh approach, the computational domain is attached to the moving object and is updated along with the object movement in the fluid [2-4]. In the context of fixed mesh approach, various numerical methods have been proposed [5-8], among which one may refer to the immersed boundary method (IBM). At first this method was introduced for modeling and simulation of blood flow in the heart and heart valves dynamic [9-13]. In this approach, the Eulerian fixed grid and the moving Lagrangian grid are used for the fluid domain and the solid object, respectively. In this approach, the immersed object in the fluid is considered as an object with flexible boundary. The fluid motion over the object would result in small deformations to its boundary. On the other hand, the object’s tendency to return to the original state would create a force on the interface between the fluid and the object. Once this force is calculated on the Lagrangian grid points representing the object, the effect of this force on the Eulerian grid points is calculated and is added to the fluid domain equations to account for the effect of the presence of the object. Therefore, in whole computational domain the Navier-Stokes (N-S) equation must be solved with the forcing term resulting from the presence of the object in fluid [14]. There have been many studies carried out to simulate the elastic/rigid immersed body motion in a fluid using the IBM. In most of these studies, the N-S equation has been solved using conventional computational fluid dynamics (CFD) methods such as finite difference, finite element and finite volume discretization methods. Fadlun et al. [15] presented finite difference IBM for modeling fluid flow around a rigid object. The finite element-IBM has been developed to study the interaction of a submerged elastic solid and a fluid flow [16, 17]. During the past two decades, the lattice Boltzmann method (LBM) has matured to an alternative and efficient numerical scheme for the simulation of fluid flows and transport problems [18-20]. Unlike conventional numerical schemes based on discretization of macroscopic continuum equations, the LBM is based on microscopic models and mesoscopic kinetic equations. The LBM has been successfully used to simulate the single viscous fluid flow, the multicomponent and multi-phase flows [21, 22], flows under the complex geometric boundary conditions [23], suspension of fine particles and flow in micro-channels [24,25], a variety of complex flows of Newtonian and non-Newtonian fluids [26,27], diffusion-convection [28], diffusion-reaction [29], wave equation [30] and Poisson's equation [31]. The LBM is not limited to the hydrodynamic equations. This method can be used for simulation of processes that are difficult to be modeled using continuous methods. These include the blood flow in complex situations, the dynamics of blood components such as red blood cells [32], platelets and materials such as artificial blood. The LBM uses a fixed, regular grid rather than a time consuming adaptive lattice. It calculates fluid properties and velocity using the local values in a way that makes it suitable for massive parallel computing [20]. Unlike the conventional numerical methods such as finite difference method, finite element method and boundary element method that solve partial differential equations (PDEs) to calculate the physical quantities directly, in LBM these quantities are calculated using a distribution function. In addition, the corresponding PDE can be recovered directly by Chapman-Enskog expansion [33]. Lately, a combination of LBM and IBM is used for simulating the fluid-solid interaction involving rigid and/or flexible boundaries. Although many studies have been conducted on the motion of rigid and elastic objects within a fluid, most of these studies have focused on unbounded flows (infinite flow domain). In the limited work carried out to simulate the motion and deformation of elastic structures in shear flow, the numerical methods other than LBM have been used. Feng and Michaelides [34] proposed the first model of IBM-LBM. Their combined method involved most of the desirable properties of LBM and IBM. It uses an Eulerian grid (lattice) for flow field and a Lagrangian grid to trace the dynamics of solid
Page 2 of 34
Ac ce pt e
2. Governing equations
d
M
an
us
cr
ip t
particles (immersed object). This method not only represents the boundary smoother than the LBM, but also prevents fluctuations of velocities and forces acting on the object that can be seen in LBM. In general, it is so powerful in solving the problems with structure deformation [34]. This method was improved by Wu and Shu [14] and they implicitly calculated and applied the forces caused by the fluid acting on the immersed boundary, which in fact had a significant impact on the cost and accuracy. They simulated the motion of a particle in a typical shear flow and the motion of two particles in a channel using their proposed method. Dupuis et al. [35] studied the flow over an impulsively started cylinder at moderate Reynolds (Re) number. They investigated how the coupling method of the forcing term between the Eulerian and Lagrangian grids could affect the results. Zhang et al. [36, 37] studied the dynamic behavior of red blood cell in shear flow and channel flow and investigated several hemodynamic and rheological properties. Cheng et al. [38] have proposed a proper model to simulate the fast boundary movements and high pressure gradient occurred in the fluid-solid interaction. In their research, mitral valve jet flow considering the interaction of leaflets and fluid has been simulated. Navidbakhsh and Rezazadeh [39] carried out a numerical study on the behavior of malaria-infected red blood cell. Vahidkhah and Abdollahi [40] simulated the motion of a massless elastic object in a two-dimensional viscous channel flow numerically using IBM-LBM. Dadvand et al. [41] investigated numerically the motion and deformation of a red blood cell in a viscous shear flow utilizing a combined LBM-IBM. Le et al. [42] used immersed interface method (IIM) to simulate the membrane motion in a two-dimensional channel flow. In the present work, the motion, deformation and rotation of two elastic membranes in a viscous shear flow in a microchannel with rectangular cross section are simulated utilizing a combined LBM-IBM. Two geometrically different channels, namely, a simple channel and a channel with a groove are considered. In the simple channel case, the membranes are placed both on the symmetry axis of the channel and off the symmetry axis.
It was mentioned in introduction that in the IBM the fluid is represented on an Eulerian coordinate and the structure is represented on a Lagrangian coordinate. A typical twodimensional example of an elastic solid membrane with curved boundary has been shown in Fig. 1. Consider a flexible solid membrane with the curved boundary immersed in the twodimensional incompressible viscous fluid domain . The membrane boundary is characterized by the Lagrangian parameter s, and the fluid domain is represented by Eulerian coordinates . Hence any point on the membrane can be written as where is arc length, and is time. Hence, the equations governing the combination of fluid and solid motions are as following: (1) (2) (3) To satisfy the no-slip boundary condition on the fluid-solid interface, the velocity of any point on the solid surface must be equal to that of the adjacent fluid particle, i.e.,
Page 3 of 34
(4) In the above equations, respectively. In addition,
and and
are the mass density and dynamic viscosity of the fluid, indicate the velocity and pressure fields, respectively. The
us
cr
ip t
term on the right-side of Eq. (2) denotes the membrane forces (tensile and bending) due to the elastic boundary immersed in the fluid. Equation (3) indicates that the force density of the fluid is obtained from the force density of the membrane . Equation (4) represents the no-slip condition at the fluidsolid interface, as the solid boundary moves with the same velocity as that of the surrounding fluid. There are various methods developed to simplify and smooth the Dirac delta function, among which the model proposed by Peskin is the most popular one [43].
M
an
where
(5)
Ac ce pt e
d
In Eq. (5), is the distance between the Eulerian nodes and in Eq. (6), is the distance between the Eulerian and Lagrangian nodes. In the IBM, to calculate the momentum exchange, the following collision function is used: (6)
where is the density distribution function of particles with the velocity located at position at time . is time step, is the equilibrium distribution function, indicates the dimensionless relaxation time and denotes the body force associated with the immersed body. The membranes are assumed to be far enough away from the front and back walls of the channel (i.e., the channel depth is very large) and hence the problem is assumed to be twodimensional (2D). Therefore, in the present research, the 2D-LBM with D2Q9 model has been used. In D2Q9 model, the particle velocity in the corresponding nine directions can be written as follows:
where Here, is the distance between two successive nodes in the Euler grid. . The equilibrium density distribution In the present research it is assumed that function is written as follow:
Page 4 of 34
(7)
cr
ip t
The fluid pressure is calculated via an isothermal equation of state ( , where is the speed of sound and is density. In addition, are weight coefficients with the following values,
The term
in Eq. (6) is the so-called BGK approximation of the collision term
us
in the Boltzmann equation. The elastic force in the lattice Boltzmann equation as,
(8)
an
.
is defined
M
In addition, the density and microscopic fluid velocity are calculated from the following relations, (9)
Ac ce pt e
d
(10)
The kinematic viscosity as follows:
in the D2Q9 model is related to the dimensionless relaxation time
.
The Lagrangian force density forces, i.e.,
(11)
comprises two parts of tension-compression
This force is related to the elastic potential energy density work theorem),
and bending
(12) as follows (thanks to the virtual
(13) where (14)
Page 5 of 34
and .
(15)
ip t
Here, and are elastic modulus (tension/compression constant) and bending modulus, respectively. The discretized form of Lagrangian force density , and elastic potential energy density will become,
In Eqs. (16)-(19) and
(17)
(18)
(19)
is the total number of Lagrangian nodes on the
are elastic Lagrangian forces associated with the node
on the
d
membrane),
(
M
an
us
cr
(16)
Ac ce pt e
membrane and is the Kronecker delta function. When the Lagrangian forces on the membrane are calculated, all the translational and rotational speeds are updated explicitly. It should be noted that the solid membrane moves continuously based on Newtonian dynamics and finally the new position of the membrane is obtained.
3 Validation
In this section, the accuracy of the current numerical results is substantiated through the comparison with the existing numerical and experimental results. Two validation cases are considered.
3.1 Validation I A conventional phenomenon that an elastic membrane may encounter when it is placed in a shear flow is the so-called tank treading (TT) motion, i.e. the membrane rotates about its own perimeter. The rotation of membrane results from fluid shear force exerted on it. The TT motion is accompanied by a lift force applied on the deformed membrane. This would cause the membrane to move toward the center of the channel provided other external forces such as gravity force (as in the present work) are absent. It should be noted that in the TT motion, after an initial deformation, the shape of the membrane would remain unchanged.
Page 6 of 34
Ac ce pt e
d
M
an
us
cr
ip t
In the first test case (validation I), an elastic membrane is placed in a Poiseuille shear flow near the bottom wall of a channel with rectangular cross section. The current results are compared with the experimental results of Fischer and Schmid-shonbin [44] and the numerical results of Frank et al. [45] (see Fig. 2). Following Fischer and Schmid-shonbin [44] and Franke et al. [45], the fluid viscosity , the membrane elastic modulus , the membrane bending modulus and the number are considered to be , , and , respectively. It may be noted that the Re number is defined as , where is the maximum velocity at the channel center without membrane and is the width of the channel. As illustrated in Fig. 2, the membrane approaches the center of channel during its translational and rotational motion in the shear flow. The circular symbol showing a pint on the membrane clearly illustrates the rotation. There is a reasonable agreement between the present numerical results (Fig. 2C) and experimental results of [44] (Fig. 2A). In addition, the present numerical result (Fig. 2C) compared to the numerical result of [45] obtained using the finite element-immersed boundary method (Fig. 2B) seems to be more accurate (the membrane deformation in the present research is closer to the experimental results [44]). The reason may be related to the fact that in this study non-uniform Lagrangian points have been used for discretization of the membrane. That is, higher density of Lagrangian points has been used on the part of the membrane located near the channel wall. This part of the membrane would deform further due to greater shear force. Due to the viscosity of the flow, boundary layers develop in the channel. In fully developed flow both the upper and lower boundary layers coincide at a point on the symmetry axis of the channel and hence the shear flow fills the whole cross-sectional area of the channel. Under such condition, the upper part of the membrane will be affected by a higher stream wise velocity than the lower part located near the lower wall of the channel. This would cause the membrane to rotate clockwise (see Fig. 2C). After a definite time, the membrane will no longer rotate and the angle will remain fixed. On the other hand, according to the relation , the lower part of the membrane is affected by a greater shear force than the upper part. This will create a lift force, which is applied to the membrane and move it towards the center of channel (see Fig. 2C).
3.2 Validation II
In the second test case (validation II), the motion and deformation of an elastic membrane placed on the centerline of an 8µm channel (see Fig. 3A) are studied. Flow rate is adjusted so that membrane velocity is approximately . The current numerical results for a single membrane (Fig. 3B) are compared with experimental observations in glass tubes [46] (Fig. 3C) and numerical result of Secomb et al. [47] (Fig. 3D). According to the numerical work of Secomb et al. [47], the initial radius of the membrane is set equal to 2.66 . The suspending fluid viscosity was set to 0.001 Pa.s, which is close to plasma viscosity. The membrane elastic shear modulus and bending modulus are taken to be and , respectively. Good qualitative agreement with the experimental configuration is seen. Since the initial position of the membrane is on the center-line of the channel ( ), it remains on the centre-line and its shape evolves to a two-dimensional parachute-like shape (convex in front and concave at the rear) over a period of about 100 ms. This change in shape is accompanied by an increase in the membrane length from the initial circular shape and a slight decrease in width (Fig. 3E).
Page 7 of 34
4 Results and discussion
ip t
Consider the Poiseuille (pressure-driven) fully developed flow in a horizontal microchannel with rectangular cross section. The objective of the present work is to study the dynamics of two flexible membranes when moving in the channel. Two geometrically different channels, namely, a simple channel and a channel with a groove are considered (see Fig. 4). In the simple channel case, the initial position of the membranes on their motion and deformation are assessed. In the channel with a groove, the membranes tensile coefficients on their behavior are investigated.
cr
4.1 Dynamics of two elastic membranes moving in a simple microchannel
us
The grid dimensions and the membrane characteristics for the simple channel case are given in Table 1.
an
4.1.1 Membranes are placed at the center of the microchannel
Ac ce pt e
d
M
In this section, two elastic membranes are placed at the center of a simple channel (see Fig. 5A). The initial position of the membrane I is (0.9, 0.75) and that of the membrane II is (0.25, 0.75). Figure 5 depicts the initial shape of the membranes along with the flow velocity vectors (see Fig. 5A) and the membranes shapes along with the pressure contours at two different times (see Figs. 5B and 5C). According to the pressure contours it can be observed that the pressure is symmetric about the axis of symmetry (the centerline of the channel). Therefore, the membranes do not migrate away from the center of the channel. However, the very small pressure difference in the flow in the direction would cause the membranes to move from left to right. It can also be seen that although the membranes experience axisymmetric pressures they undergo different deformations. This may be attributed to the variations of the pressure along the flow direction. The overall motion and deformation of the membrane I and the membrane II in the channel have been shown in Figs. 6A and 6B, respectively. It can be observed that membrane II experiences greater deformation than membrane I due to the (relatively) higher pressure applied to the rear part of the membrane II. In order to quantitatively justify this statement, variations with time of the length/width ratio (a/b) representing deformation intensity of the membrane are plotted in Fig. 7. It can be seen from this figure that at all times membrane II has a lower length/width ratio than membrane I, indicating that membrane II experiences greater deformation than membrane I. This is again attributed to the (relatively) smaller pressure applied to the rear part of the membrane I and (relatively) higher pressure applied to the rear part of the membrane II.
4.1.2 Membranes are placed asymmetrically near the walls of a simple microchannel In this section, the two elastic membranes are placed asymmetrically near the walls of a simple microchannel (see Fig. 8A). The initial positions of the membrane I (lower membrane) and membrane II (upper membrane) are respectively (0.25, 0.25) and (0.25, 1.26). The results are represented in terms of flow velocity magnitude and pressure contours.
Page 8 of 34
us
cr
ip t
It can be seen that at time , the upper membrane (membrane II) has been elongated more than the membrane I. This may be due to the fact that a larger portion of the membrane II has been placed in the low velocity region near the upper wall of the channel (see velocity contours in Fig. 8B) and hence it demands greater force to move as compared to the membrane I (see pressure contours in Fig. 8B). The higher velocity magnitude of the flow around the lower membrane would cause it to move faster than the upper one (see Fig. 8C). The vertical movements of the centroid of the membranes are depicted in Fig. 9. It can be observed that the membrane I migrates toward the centerline of the channel faster than the membrane II. The overall motion and deformation of the membrane I and the membrane II in the channel have been shown in Fig. 10. It can be observed that the membrane I has moved a longer distance and has further migrated toward the centerline of the channel as compared with the membrane II. This may be attributed to the larger initial distance of the membrane I from the wall.
4.1.3 Membranes are placed symmetrically near the walls of a simple microchannel
Ac ce pt e
d
M
an
In this section, the two elastic membranes are located symmetrically near the walls of a simple microchannel (see Fig. 11A). The initial positions of the membrane I and membrane II are respectively (0.25, 0.25) and (0.25, 1.25). The results are represented in terms of and pressure contours. dimensionless flow velocity It is observed that the two membranes migrate toward the centerline of the channel simultaneously. Their translational and rotational speeds as well as their deformation are the same. This may be due to the symmetric fashion of the pressure field exerted on them by the flow field. The vertical movements of the centroid of the membranes associated with Fig. 11 are depicted in Fig. 12. It can be observed that the two membranes migrate toward the centerline of the channel simultaneously and sweep the same distance measured from the channel inlet. The overall motion and deformation of the membrane I and the membrane II in the channel associated with symmetric case have been shown in Fig. 13. It can be observed that the membrane I has moved the same distance and has migrated toward the centerline of the channel as equally as the membrane II. This may be attributed to the same initial distance of the membranes from the centerline of the channel. Comparison of the asymmetric and symmetric cases discussed respectively in sections 4.1.2 and 4.1.3 reveals that the membrane II in the asymmetric case has less influence on the membrane I as it moves a distance away from it.
4.2 Dynamics of two elastic membranes when moving in a grooved microchannel In this section, two elastic membranes are placed in the groove of a grooved channel (see Fig. 4C). The flow through the channel can under certain conditions move the membranes out of the groove. These conditions encompass the initial positions of the membranes, flow rate, depth of the groove and the elasticity of the membranes. Two cases are considered: in the first case (section 4.2.1), the depth and width of the groove is considered to be 0.9 and 1.4 , respectively. Besides, the initial positions of the membrane I and membrane II are considered to be (2.5, 0.22) and (1.8, 0.22), respectively. In the second case (section 4.2.2), all the conditions are considered to be the same as the first case, but the tensile modulus of
Page 9 of 34
the membrane II is changed. In both cases the tensile modulus of the membrane I is set to be 0.1.
4.2.1 Case I: the tensile modulus of the membrane II is equal to 0.1
cr
4.2.2 Case II: the tensile modulus of the membrane II is equal to 0.3
ip t
In this Case I, the tensile modulus of the membrane II is equal to 0.1, which is equal to that of the membrane I. Other physical features of the membranes are identical to those given in Table 1. In this case, both of the membranes are able to move out of the groove due to the proper pressure differences made in the flow around the membranes (see Fig. 14).
5 Conclusions
M
an
us
In this Case II, the tensile modulus of the membrane II is increased to 0.3, which is greater than that of the membrane I. Other conditions are identical to the Case I. This would result in a decrease in the elasticity of the membrane II. Therefore, the membrane II only rotates in the groove and is unable to move out of the groove due to the inappropriate pressure differences made in the flow around it (see Fig. 15). This may be due to the fact that the Lagrangian forces exerted by the membrane II on the flow field are too small to create proper lift force and move the membrane II out of the groove.
Ac ce pt e
d
The motion, deformation and rotation of two elastic membranes in a viscous shear flow in a microchannel with and without a groove are simulated utilizing a combined LBM-IBM. The membranes are considered as immersed elastic boundaries in the fluid flow. The membranes are represented in Lagrangian coordinates, while the fluid flow field is discretized by a uniform and fixed Eulerian mesh. The interaction of the fluid and membranes is modeled using an appropriate form of the Dirac delta function. Two geometrically different channels, namely, a simple channel and a channel with a groove are considered. In the simple channel case, when the membranes are placed on the symmetry axis of the channel, they continue to move and deform without any lift force and rotation induced. However, when the membranes are located off the symmetry axis, the pressure difference produced in the flow around the membranes would apply lift forces on them and expel them towards the center of the channel. In the case of channel with a groove, it was found that when a membrane tensile modulus was reduced its flexibility as well as rotational speed would increase. This would, in turn, result in a pressure difference in the flow around the membrane. The pressure difference would apply a lift force on the membrane and move it out of the groove. The results of the present simulation were compared with existing numerical and experimental results, which had a good accordance with them. Moreover comparing the current numerical results with LBM-FE results of [45] indicates that the present results are more accurate. The reason may be related to the fact that in present study the non-uniform Lagrangian points have been used for discretization of the immersed boundaries.
References
Page 10 of 34
Ac ce pt e
d
M
an
us
cr
ip t
[1] S. Nadeem, F. Naeem, Thin film flow of a second grade fluid over a stretching/shrinking sheet with variable temperature-dependent viscosity, Chinese Physics. 27(3) (2010) 34704-34707 [2] H.H. Hu, Direct simulation of flows of solid-liquid mixtures, International Journal of Multiphase Flow. 22(2) (1996) 335-352 [3] B. Maury, Direct simulations of 2D fluid-particle flows in biperiodic domains, Journal of computational physics. 156 (2) (1999) 325-351 [4] H.H. Hu, N.A. Patankar, M. Zhu, Direct numerical simulations of fluid–Solid systems using the arbitrary Lagrangian–Eulerian technique, Journal of Computational Physics. 169(2) (2001) 427- 462 [5] R. Glowinski, T.W. Pan, T.I. Hesla, D.D. Joseph, J. Periaux, Distributed Lagrange multiplier/fictitious domain method for flows around moving rigid bodies: application to particulate flow, International Journal for Numerical Methods in Fluids. 30(8) (1999)1043-1066 [6] R. Glowinski, T.W. Pan, T.I. Hesla, D.D. Joseph, Distributed Lagrange multiplier/fictitious domain method for particulate flows, International Journal of Multiphase Flow. 25(5) (1999)755-794 [7] Z. Yu, N. Phan-Thien, Y. Fan, R. I. Tanner, Viscoelastic mobility problem of system of particles, Journal of non-Newtonian fluid mechanics. 104(2) (2002)87-124 [8] J. Wu, C. Shu, Particulate flow simulation via boundary condition-enforced immersed boundary-lattice Boltzmann scheme, Communications in Computational Physics. 7(4) (2010)793 [9] C.S. Peskin, Flow patterns around heart valves: a digital computer method for solving the equations of motion, PhD thesis, Albert Einstein College of Medicine, Chicago (1972) [10] C.S. Peskin, Numerical analysis of blood flow in the heart, Journal of Computational Physics. 25(3) (1977) 220-252 [11] C.S. Peskin, The fluid dynamics of heart valves, experimental, theoretical and computational methods, Annual Review of Fluid Mechanics. 14(1) (1982) 235–259 [12] D.M. McQueen, C.S. Peskin, E.L. Yellin, Fluid dynamics of the mitral valve, physiological aspects of a mathematical model, American Journal of Physiology-Heart and Circulatory Physiology. 242(6) (1982)1095–1110 [13] D.M. Mc Queen, C.S. Peskin, Computer-assisted design of butterfly bileaflet valves for the mitral position, Scandinavian Cardiovascular Journal. 19(2) (1985)139–148 [14] J. Wu, C. Shu, Particulate flow simulation via boundary condition-enforced immersed boundary-lattice Boltzmann scheme, Communications in Computational Physics. 7(4) (2010) 793 [15] E. Fadlun, R. Verzicco, P. Orlandi, J. Mohd-Yusof, Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations, Journal of Computational Physics. 161(1) (2000) 35-60 [16] X. Wang, W.K. Liu, Extended immersed boundary method using FEM and RKPM, Computer Methods in Applied Mechanics and Engineering. 193(12) (2004) 1305-1321 [17] L. Zhang, A. Gerstenberger, X. Wang, W.K. Liu, Immersed finite element method, Computer Methods in Applied Mechanics and Engineering. 193(21) (2004) 2051-2067 [18] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annual Review of Fluid Mechanics. 30(1) (1998) 329-364 [19] S. Succi, The lattice Boltzmann equation for fluid dynamics and beyond: Oxford University press. 2001 [20] M.C. Sukop, D.T. Thorne Jr., Lattice Boltzmann modeling: an introduction for geoscientists and engineers, New York, 2007
Page 11 of 34
Ac ce pt e
d
M
an
us
cr
ip t
[21] H. Liu, Y. Ju, N. Wang, G. Xi, Y. Zhang, Lattice Boltzmann modeling of contact angle and its hysteresis in two-phase flow with large viscosity difference, Physical Review E. 92 (2015) 033306 [22] H. Huang, L. Wang, X.Y. Lu, Evaluation of three lattice Boltzmann models for multiphase flows in porous media, Computers and Mathematics with Applications. 61(12) (2011) 3606-3617 [23] R.W. Nash, H.B. Carver, M.O. Bernabeu, J. Hetherington, D. Groen, T. Krüger, P.V. Coveney, Choice of boundary condition for lattice-Boltzmann simulation of moderateReynolds-number flow in complex domains, Physical Review E. 89(2) (2014) 023303. [24] H. Hassanzadeh Afrouzi, K. Sedighi, M. Farhadi, A. Moshfegh, Lattice Boltzmann analysis of micro-particles transport in pulsating obstructed channel flow, Computers and Mathematics with Applications. 70(5) (2015) 1136-1151 [25] Y. Bakhshan, A. Omidvar, Calculation of friction coefficient and analysis of fluid flow in a stepped micro-channel for wide range of Knudsen number using Lattice Boltzmann (MRT) method, Physica A: Statistical Mechanics and its Applications. 440(15) (2015) 161-175 [26] C. Cao, S. Chen, J. Li, Z. Liu, L. Zha, S. Bao, C. Zheng, Simulating the interactions of two freely settling spherical particles in Newtonian fluid using lattice-Boltzmann method, Applied Mathematics and Computation. 250(1) (2015) 533-551 [27] M. Yoshino, Y. Hotta, T. Hirozane, M. Endo, A numerical method for incompressible non-Newtonian fluid flows based on the lattice Boltzmann method, J. Non-Newtonian Fluid Mech. 147(1) (2007) 69-78 [28] J. Huang, W.A. Yong, Boundary conditions of the lattice Boltzmann method for convection–diffusion equations, Journal of Computational Physics. 300(1) (2015) 70-91 [29] L. Chen, Q. Kang, B. Carey, W.Q. Tao, Pore-scale study of diffusion–reaction processes involving dissolution and precipitation using the lattice Boltzmann method, International Journal of Heat and Mass Transfer. 75 (2014) 483-496 [30] G.W. Yan, A lattice Boltzmann equation for waves, J. Computational Physics. 161 (2000) 61–69 [31] H. Wang, G. Yan, B. Yan, Lattice Boltzmann model based on the rebuilding-divergency method for the Laplace equation and the Poisson equation, Journal of Scientific Computing. 46(3) (2011) 470-484 [32] P. Bagchi, Mesoscale simulation of blood flow in small vessels, Biophysical Journal. 92(6) (2007) 1858–1877 [33] A.A. Mohamad, Lattice Boltzmann Method: Fundamentals and Engineering Applications with Computer Codes, New York, 2011 [34] Z.G. Feng, E.E. Michaelides, The immersed boundary-lattice Boltzmann method for solving fluid–particles interaction problems, Journal of Computational Physics. 195(2) (2004)602-628 [35] A. Dupuis, P. Chatelain, P. Koumoutsakos, An immersed boundary-lattice-Boltzmann method for the simulation of the flow past an impulsively started cylinder, Journal of Computational Physics. 227(9) (2008) 4486-4498 [36] J.F. Zhang, P.C. Johnson, A.S. Popel, An immersed boundary lattice Boltzmann approach to simulate deformable liquid capsules and its application to microscopic blood flows, Physical Biology. 4(4) (2007) 285-295 [37] J.F. Zhang, P.C. Johnson, A.S. Popel, Red blood cell aggregation and dissociation in shear flows simulated by lattice Boltzmann method, Journal of Biomechanics. 41(1) (2008) 47-55 [38] Y. Cheng, H. Zhang, Immersed boundary method and lattice Boltzmann method coupled FSI simulation of mitral leaflet flow, Computers and Fluids. 39(5) (2010) 871-881
Page 12 of 34
Ac ce pt e
d
M
an
us
cr
ip t
[39] M. Navidbakhsh, M. Rezazadeh, An immersed boundary-lattice Boltzmann model for simulation of malaria-infected red blood cell in micro-channel, Scientia Iranica. 19(5) (2012) 1329-1336 [40] K. Vahidkhah, V. Abdollahi, Numerical simulation of a flexible fiber deformation in a viscous flow by the immersed boundary-lattice Boltzmann method, Communications in Nonlinear Science and Numerical Simulation. 17(3) (2012) 1475-1484 [41] A. Dadvand, M. Baghalnezhad, I. Mirzaee, B.C. Khoo, S. Ghoreishi, An immersed boundary–lattice Boltzmann approach to study the dynamics of elastic membranes in viscous shear flows, Journal of Computational Science, 5(5) (2014), 709-718 [42] D.V. Le, B.C. Khoo, J. Peraire, An immersed interface method for viscous in compressible flows involving rigid and flexible boundaries, Journal of Computational Physics. 220(1) (2006) 109-138 [43] C.S. Peskin, The immersed boundary method, Acta Numerica. 11 (2002) 479-517 [44] T. Fischer, H. Schmid-Schönbein, Tank-trading motion of red blood cell membranes in viscometric flow: behavior of intracellular and extracellular markers, Blood Cells. 3 (1977) 351-365 [45] T. Franke, R.H. Hoppe, W. Linsenmann, C. Schmid, L. Willbold, A. Wixforth, Numerical simulation of the motion and deformation of red blood cells and vesicles in microfluidic flows, Computing and Visualization in Science. 14(4) (2011) 167-180 [46] T.W. Secomb, Mechanics and computational simulation of blood flow in microvessels, Medical engineering & physics. 33(7) (2011) 800-804 [47] T.W. Secomb, B. Styp-Rekowska, A.R. Pries, Two-dimensional simulation of red blood cell deformation and lateral migration in microvessels, Annals of Biomedical Engineering. (35) (2007) 755-765
Page 13 of 34
ip t cr us
Table 1 Grid dimensions and membranes characteristics
300 100 60 0.0000004 0.1 N/m 0.001 N/m 0.22 µm 0.00001
Ac ce pt e
d
M
an
Number of Eulerian points in the x-direction Number of Eulerian points in the y-direction Number of Lagrangian points Time step Membrane tensile coefficient Membrane bending coefficient Initial radius of the membrane Gravity force
Velocity interpolation
Spread of Lagrangian force
Page 14 of 34
an
us
cr
ip t
Fig. 1 Schematic representation of immersed boundary (Lagrangian coordinates) and Eulerian mesh for fluid and spreading the Lagrangian force from the membrane boundary points to the Eulerian nodes.
Ac ce pt e
d
M
50 µm
(A)
y(µm)
150
100 80 60 40 20 0
0
100
300
500
x(µm)
(B)
Page 15 of 34
ip t
(C)
Ac ce pt e
d
M
an
us
cr
Fig. 2 Tank-treading motion of an elastic membrane in a channel shear flow; (A) Experimental results of [44], (B) Numerical results of [45] and (C) Present numerical results obtained using the same code as the one employed by Dadvand et al. [41] but with different Dirac delta function
(A) (B)
(C)
(D)
Page 16 of 34
8
7 Length
4
Width
0
20
40
t (ms)
60
80
cr
5
ip t
6
100
us
Length,Width (µ m)
Numerical [47] Present study ZONE 001 ZONE 001
d
M
an
(E) Fig. 3 Motion and deformation of a membrane (initial radius, when placed on the centerline of an 8-µm channel (A); Comparison between the current results (B), the experimental data of Secomb [46] (C) and numerical results of Secomb et al. [47] (D); E: Variations of length and width of the membrane with time.
I
Ac ce pt e
II
H
(A)
II
I
(B)
Page 17 of 34
II
I
D
ip t
W (C)
Ac ce pt e
d
M
an
us
cr
Fig. 4 Simple channels (A) and (B) with different positioning of the membranes, and channel with a groove (C)
(5A)
Page 18 of 34
us
cr
ip t
(5B)
Ac ce pt e
d
M
an
(5C) Fig. 5 Pressure contours associated with the motion and deformation of two elastic membranes at the center of a fully developed Poiseuille flow in a channel
(6A)
(6B)
Page 19 of 34
M
an
us
cr
ip t
Fig. 6 Motion and deformation of the membrane I (A) and the membrane II (B) in the channel when they are placed on the symmetry axis of the channel
1
Ac ce pt e
d
0.96
a/b
0.92
0.88
Membrane II Membrane I
0.84
0.8
0
0.001
0.002
0.003
t
0.004
0.005
0.006
Fig. 7 Variation of length/width ratio with time representing membrane deformation intensity associated with the membranes considered in Fig. 6
Page 20 of 34
d
Ac ce pt e (8A)
Page 21 of 34
us
an
M
cr
ip t
d
Ac ce pt e (8B)
Page 22 of 34
us
an
M
cr
ip t
ip t cr us an M d Ac ce pt e
(8C)
Fig. 8 Motion and deformation of two membranes placed asymmetrically near the walls of a simple microchannel
Page 23 of 34
1.25
Membrane I Membrane II
ip t
1
cr
ycenter
0.75
us
0.5
0 0.002
M
0
an
0.25
t
0.004
0.006
Ac ce pt e
d
Fig. 9 Vertical movement of the centroid of the membranes associated with Fig. 8
Fig. 10 Sequential frames of the membranes motion when they are placed asymmetrically in the channel
Page 24 of 34
ip t
Ac ce pt e
d
M
an
us
cr
(11A)
(11B)
Page 25 of 34
ip t cr us an M d Ac ce pt e
(11C)
Fig. 11 Motion and deformation of two membranes placed symmetrically near the walls of a simple microchannel
Page 26 of 34
1.25
Membrane I Membrane II
1
ip t
ycenter
0.75
us
cr
0.5
0 0
0.002
an
0.25
t
0.004
0.006
Ac ce pt e
d
M
Fig. 12 Vertical movement of the centroid of the membranes associated with Fig. 11
Fig. 13 Sequential frames of the membranes motion when they are placed symmetrically in the simple channel.
Page 27 of 34
ip t cr us
Ac ce pt e
d
M
an
(14A)
(14B)
Page 28 of 34
ip t cr us
Ac ce pt e
d
M
an
(14C)
(14D)
Fig. 14 Sequential snapshots of the motion of two identical elastic membranes in a channel with a groove
Page 29 of 34
ip t cr us
Ac ce pt e
d
M
an
(15A)
(15B)
Page 30 of 34
ip t cr us
(15C)
Ac ce pt e
d
M
an
Fig. 15 Sequential snapshots of the motion of two different elastic membranes in a channel with a groove
Page 31 of 34
*Biographies (Text)
Asad Alizadeh is from Bookan city, Iran. He is currently M.Sc. student at the Urmia University
ip t
of Technology, Urmia, Iran in Mechanical Engineering (energy conversion). His thesis is on the combined immersed boundary-lattice Boltzmann method to simulate elastic membranes in shear flows.
Ac ce p
te
d
M
an
us
cr
Abdolrahman Dadvand, assistant professor, received his PhD from the Department of Mechanical Engineering at the University of Tabriz, Tabriz, Iran with his dissertation on spark bubble micro-droplet generation. He continued his PhD research at the National University of Singapore (NUS) as an attachment student. He joined Urmia University of Technology (UUT), Urmia, Iran in 2010 and since then, he has worked on CFD methods (such as Boundary Element, Finite Volume, Lattice Boltzmann and Immersed Boundary methods). He has published more than 50 peer reviewed Journal and conference papers. His current projects focus on the hybrid computational methods.
Page 32 of 34
*Biographies (Photograph)
Ac ce p
te
d
M
an
us
cr
ip t
Abdolrahman Dadvand Asad Alizadeh
Page 33 of 34
*Highlights (for review)
Highlights o The motion, deformation and rotation of two elastic membranes in a Poiseuille shear flow are simulated utilizing a combined LBM-IBM.
ip t
o Two geometrically different channels, namely, a simple channel and a channel with a groove are considered.
cr
o In the simple channel case, when the membranes are placed on the symmetry axis of the channel, they continue to move along the axis. However, when the membranes are located off the symmetry axis, they migrate towards the center of the channel. In the case of channel with a groove, it was found that when the membrane tensile modulus was reduced it moves out of the groove.
o
The results of the present simulation were compared with existing numerical and experimental results, which had a good accordance with them.
Ac ce p
te
d
M
an
us
o
Page 34 of 34