**Surfactant Analysis of Thin Liquid Film in the Human Trachea via Application of Volume of Fluid (VOF)**

Sujudran Balachandran *Bumi Armada Berhad Malaysia/Singapore* 

#### **1. Introduction**

448 Fluid Dynamics, Computational Modeling and Applications

H.W. Oh (Ed.), 391-404. INTECH Education and Publishing, Vienna. Marinho, D.A.; Reis, V.M.; Vilas-Boas, J.P.; Alves, F.B.; Machado, L.; Rouboa, A.I. & Silva, A.J.

Fluid Dynamics. *Brazilian Archives of Biology and Technology*, 5(2), 437-442. Marinho, D.A.; Barbosa, T.M.; Reis, V.M.; Kjendlie, P.L.; Alves, F.B.; Vilas-Boas, J.P.;

swimming. *Journal of Biomechanics*, 42, 13, 2188-2190

*Journal of Human Sport and Exercise*, 6, 1, 87-93.

fluid dynamics. *Journal of Biomechanics*, 39, 1239-1248.

*Journal of Sports Science and Medicine*, 7, 1, 60-66.

swimming. *Animal Biology*, 55, 1, 17-40.

but first think about it. *American Swimming Magazine*, 5, 23-32.

*Sport Sciences*, 6, Suppl. 2, 185-189.

Oxford.

35, 519-524.

*Biomechanics*, 15*,* 3-26.

enhanced by a small finger spread. *Journal of Applied Biomechanics*, 26, 87-92. Minetti, A.E.; Machtsiras, G. & Masters, J.C. (2009). The optimum finger spacing in human

Moreira, A.; Rouboa, A.; Silva, A.; Sousa, L.; Marinho, D.; Alves, F.; Reis, V.; Vilas-Boas, J.P.;

Pallis, J.M.; Banks, D.W. & Okamoto, K.K. (2000). 3D computational fluid dynamics in

Pendergast, D.R.; Capelli, C.; Craig Jr, A.B.; di Prampero, P.E.; Minetti, A.E.; Mollendorf1, J.;

Roberts, B.S.; Kamel, K.S.; Hedrick, C.E.; MLean, S.P. & Sharpe, R.L. (2003). Effect of Fastskin

Rouboa, A.; Silva, A.; Leal, L.; Rocha, J. & Alves, F. (2006). The effect of swimmer's

Sanders, R.H. (1999). Hydrodynamic characteristics of a swimmer's hand. *Journal of Applied* 

Sanders, R.B.; Rushall, H., Toussaint, H., Steager, J. & Takagi, H. (2001). Bodysuit yourself

Sato, Y & Hino, T. (2002). Estimation of thrust of swimmer's hand using CFD. In: *Proceedings of 8th symposium on nonlinear and free-surface flows*, 71-75. Hiroshima. Schleihauf, R.E. (1979). A hydrodynamic analysis of swimming propulsion. In: *Swimming III*, J. Terauds & E.W. Bedingfield (Eds.), 70-109. University Park Press, Baltimore. Silva, A.J.; Rouboa, A.; Moreira, A.; Reis, V.; Alves, F.; Vilas-Boas, J.P. & Marinho, D. (2008).

Takagi, H.; Shimizu, Y.; Kurashima, A. & Sanders, R. (2001). Effect of thumb abduction and

around a cylinder. *Portuguese Journal of Sport Sciences*, 6(Suppl. 1), 105. Neiva, H.P.; Vilas-Boas, J.P.; Barbosa, T.M.; Silva, A.J. & Marinho, D.A. (2011). 13th FINA

swimming using computacional fluid dynamics, In: *Computational Fluid Dynamics*,

(2010b). Design of a three-dimensional hand/forearm model to apply Computational

Machado, L.; Silva, A.J. & Rouboa, A.I. (2010c). Swimming propulsion forces are

Carneiro, A. & Machado, L. (2006). Computational analysis of the turbulent flow

World Championships: Analysis of Swimsuits Used By Elite Male Swimmers.

competitive sail, yatch and windsurfer design, In: *The Engineering of Sport: Research, Development and Innovation*, F. Subic & M. Haake (Eds.), 75-79. Blackwell Science,

Termin, A. & Zamparo, P. (2006). Biophysics in swimming. *Portuguese Journal of* 

suit on submaximal freestyle swimming. *Medicine and Science in Sports and Exercise*,

hand/forearm acceleration on propulsive forces generation using computational

Analysis of drafting effects in swimming using computational fluid dynamics.

adduction on hydrodynamic characteristics of a model of the human hand. In: *Proceedings of Swim Sessions of the XIX International Symposium on Biomechanics in Sport*, J. Blackwell & R. Sanders (Eds.), 122-126. University of San Francisco, San Francisco. Tamura, H.; Nakazawa, Y.; Sugiyama, Y.; Nomura, T. & Torii, N. (2002). Motion analysis and shape evaluation of a swimming monofin. *The Engineering Sports*, 4, 716-724. Toussaint, H. & Truijens, M. (2005). Biomechanical aspects of peak performance in human The human tracheobronchial tree is a complex branched distribution system in charge of renewing the air inside the acini which are the gas exchange units. The surfactant factor existing in the acini of the human tracheobronchial tree is exposed to thin liquid film flow where the application of volume of fluid (VOF) can significantly determine the membrane effects on the branching asymmetry. The thin film application in this chapter will focus on the breathing airway which is commonly known as the windpipe (trachea) of an adult human. The upper human airway is the primary conduit for inspiration in the breathing process. Air entering the mouth passes through pharynx and flows into the trachea via the glottal region. The air which enters the windpipe applies a surfactant pressure at the lipoprotein complex, where this complex is developed as a thin layer on the trachea (Caro et al. 2002). The surfactant is a lipo-protein complex, which is a highly surface-active material found in the fluid lining of the air-liquid interface in the trachea surface. This surfactant plays a dual function of preventing sudden collapse during the breathing cycle and protection from injuries and infections caused by foreign bodies and pathogens. The varying degrees of structure-function abnormalities of surfactant have been associated with obstructive trachea diseases, respiratory infections, respiratory distress syndromes, interstitial lung diseases, pulmonary alveolar proteinosis, cardiopulmonary bypass surgery and smoking habits. For some of the pulmonary conditions, especially respiratory distress syndrome, surfactant therapy is on the horizon. In order to understand the behaviour and relevant condition of the surfactant in the human trachea, it is important to apply the volume of fluid method on these surfaces. The phenomena that occur on the trachea will ensure that the surfactant responsibility in resolving the potential obstruction of breathing. The surfactant factor may occur with non-lateral conditions in space as well as during the inspiration of breathing in the human body. The fluid interaction in the thin surface results in serious impairment by obstructive trachea diseases as mentioned earlier. The pulmonary surfactant is essential for normal breathing, alveolar stability and as a host defence system in the lungs. The interface of surfactant films reduces the surface tension to extremely low values when it is compressed during expiration. This protects our lungs from collapse during breathing out. Thus application of the Volume of Fluid (VOF) method is introduced in this paper to study the behaviour of the pulmonary surfactant in the human trachea.

Surfactant Analysis of Thin Liquid Film in the

the boundary layer becomes significant in this region.

Fig. 3a. Mesh concentration on lung model

Fig. 3b. Mesh concentration on lung model (isometric view)

The categorization of the anatomy of tracheobronchial tree airways constitutes the first step in the examination of the respiratory flow field and surfactant impact with volume of fluid

**3. CFD simulation** 

Human Trachea via Application of Volume of Fluid (VOF) 451

been included. First, the airway dimensions at the first generations (G0-G0 and G1-G1) are specific to the human anatomy and they are essentially independent of physiological variability (Martonen et al. 2003). Second, for higher generations (G2-G2 and G3-G3), a systematic branching asymmetry has been modelled in the different tree bifurcations. Figure 3a and 3b illustrate the overall mesh of the model selected. The meshes are generated using a tetrahedral unstructured mesh format. The total number of grids generated is approximately 0.62 million. Concentrated mesh is applied near the branches of the bifurcation to increase the result interpretation. The flow details are critical as the effect of

## **2. Model selection**

Figure 1 shows the three-dimensional model of the lung. Detailed geometries of this model were extracted from the anatomical model by Schmidt et al. (2004). This explicit human lung was derived from High-Resolution Computed Tomography (HRCT) imaging of an in vitro preparation. This model was a selection of data extracted from the anatomy of a healthy human lung of an adult male. This lung model is free from pathological alteration (Gemci et al. 2008). In Figure 2, an isometric view of the human lung is presented. The notation of G0, G1, G2 and G3 are the cross-sectional gate to facilitate the model drawn for indication segments. This model was selected based on the actual structure of the human lung. No additional tissue configuration was created due to investigation which only focused on the airway blockage area. At each generation, the branching is essentially dichotomous; each airway being divided into two smaller daughter airways (Farkas et al. 2007). The tree starts at the trachea (G0-G0) whose average diameter and length are respectively D0 = 1.8 cm and L0 = 12 cm in the healthy human adult (Allen et al. 2004), and ends in the terminal bronchioles. From the trachea to terminal bronchioles, this element is located on average in generation G1-G1, these airways are purely conducting pipes. Two important features have

Fig. 1. General three-dimensional lung airway model

Fig. 2. Isometric view of lung airway

Figure 1 shows the three-dimensional model of the lung. Detailed geometries of this model were extracted from the anatomical model by Schmidt et al. (2004). This explicit human lung was derived from High-Resolution Computed Tomography (HRCT) imaging of an in vitro preparation. This model was a selection of data extracted from the anatomy of a healthy human lung of an adult male. This lung model is free from pathological alteration (Gemci et al. 2008). In Figure 2, an isometric view of the human lung is presented. The notation of G0, G1, G2 and G3 are the cross-sectional gate to facilitate the model drawn for indication segments. This model was selected based on the actual structure of the human lung. No additional tissue configuration was created due to investigation which only focused on the airway blockage area. At each generation, the branching is essentially dichotomous; each airway being divided into two smaller daughter airways (Farkas et al. 2007). The tree starts at the trachea (G0-G0) whose average diameter and length are respectively D0 = 1.8 cm and L0 = 12 cm in the healthy human adult (Allen et al. 2004), and ends in the terminal bronchioles. From the trachea to terminal bronchioles, this element is located on average in generation G1-G1, these airways are purely conducting pipes. Two important features have

G3-G3

G0-G0

**2. Model selection** 

Fig. 1. General three-dimensional lung airway model

G1-G1

G2- G2

Fig. 2. Isometric view of lung airway

been included. First, the airway dimensions at the first generations (G0-G0 and G1-G1) are specific to the human anatomy and they are essentially independent of physiological variability (Martonen et al. 2003). Second, for higher generations (G2-G2 and G3-G3), a systematic branching asymmetry has been modelled in the different tree bifurcations.

Figure 3a and 3b illustrate the overall mesh of the model selected. The meshes are generated using a tetrahedral unstructured mesh format. The total number of grids generated is approximately 0.62 million. Concentrated mesh is applied near the branches of the bifurcation to increase the result interpretation. The flow details are critical as the effect of the boundary layer becomes significant in this region.

Fig. 3a. Mesh concentration on lung model

Fig. 3b. Mesh concentration on lung model (isometric view)

## **3. CFD simulation**

The categorization of the anatomy of tracheobronchial tree airways constitutes the first step in the examination of the respiratory flow field and surfactant impact with volume of fluid

Surfactant Analysis of Thin Liquid Film in the

be computed through the following constrains:

each cell is given by:

following form:

through the density

The value of

values:

the following equation:

on the volume condition 2 1

**3.1.2 The momentum equation** 

domain. The momentum equation is:

equation dependant on the volume fraction

 

, the following equations are obtained since

and viscosity

Human Trachea via Application of Volume of Fluid (VOF) 453

. *<sup>q</sup> <sup>q</sup> q*

The primary phase of the fluid will not be solved by the volume fraction equation, but will

*v t* 

> 1 <sup>1</sup> *<sup>n</sup> q q*

For each of the phases, the properties will be presented in control volume of the respective component phase. For the two phase system, for example, if the phases are represented by the subscripts 1 and 2, and if the volume fraction of these is being tracked, the density in

22 21

and for every n – phase system , the volume fraction averaged density takes on the

As mentioned earlier, both fluids will share the same momentum equation throughout the

*j i j j i j ij i u u P u u u g B*

Velocity component and direction denoted by *uj* and *<sup>j</sup> x* in j direction and t is the time. The summation of convention is applied when an index is repeated and i, j = 1, 2 for this problem. Pressure and gravitational acceleration are denoted by P and g respectively and B is a source term due to any additional force applied (per unit volume). This momentum

22 2 0 *<sup>i</sup>*

*u Dt t x*

*i*

1. By applying the properties of density

 and 

<sup>1</sup> for the first phase is not solved directly at any stage, but it is solved based

 

 

, *j j <sup>i</sup>*

 

*t u x xx x*

 

*D*

 

 

*q*

(1)

(2)

( ) 1 , (3)

*q q* (4)

*<sup>q</sup>* where in this case *q* = 1, 2 of the phases

. The computational domain is calculated through

(6)

2 2 11 , (7)

depend on the volume fraction

and viscosity

(5)

*S*

(VOF). The intention of this study is to understand the behaviour of the membrane which acts as a thin liquid film and the role of the surfactant on the model. In this study, the homogenous airflow is a comprehensive digital reference model with a maximum of four sections of pulmonary airway tree computed using the commercial CFD code FLUENT® (version 6.2). The computational mesh was generated using the FLUENT mesh generation code GAMBIT®. FLUENT employs a finite-volume method to solve the Navier-Stokes and continuity equations on an arbitrarily shaped flow domain with appropriate boundary conditions. The steady-state solution of the flow fields assumed converged when the residuals reduced to less than 10-4. A typical run time was approximately 38 hours on a normal processor. In this CFD study of surfactant in the human pulmonary tree model, computations were performed at 28.3 L/min for a quasi steady-state volumetric adult inhalation flow rate for pulmonary surfactant.

#### **3.1 Volume of Fluid (VOF) method in Fluent 6.2**

The Volume of Fluid (VOF) model has been widely used in many research fields and commercial software. Known for its wide range of usage, implementation of this model in commercial code has been popular and FLUENT 6.2 is one of the software where their codes have been implemented. Isothermal VOF formulation was considered in this project with heat transfer and is assumed negligible due to analysis being based on surface interface. Formulation of laminar flow is presented here. For the turbulence effect, it can be done by averaging the Navier Stokes equation and detailed formulation can be found in Versteeg and Malalsekera (1995).

#### **3.1.1 The basic of VOF for interface development**

The VOF model is based on a volume fraction denoted by where this volume of fraction is using a two fluids mixture in a fixed computational grid. The transport of the two fluids results in solving the multiphase fluid equation using the values and interface location is reconstructed from the volume fraction field. The following procedures are used in solving the VOF model:


All these steps are repeated until a convergence criterion is achieved. In VOF, especially in FLUENT 6.2, one can conduct a simulation for two different fluids without mixing the respective fluids. There will be an interface region between the both liquids and for each region *<sup>q</sup>* of 1 and 0 values will be given.

*<sup>q</sup>* = 0: the cell is empty

*<sup>q</sup>* = 1: the cell is full

0 < *<sup>q</sup>* < 1: the cell contains the interface between the qth fluid and one or more other fluids.

Using the same momentum equation, both fluids shared the same velocity field and it is solved throughout their domain. The tracking interface between both phases is accomplished by the solution of a continuity equation for the volume fraction of one of the phases. The equation is:

(VOF). The intention of this study is to understand the behaviour of the membrane which acts as a thin liquid film and the role of the surfactant on the model. In this study, the homogenous airflow is a comprehensive digital reference model with a maximum of four sections of pulmonary airway tree computed using the commercial CFD code FLUENT® (version 6.2). The computational mesh was generated using the FLUENT mesh generation code GAMBIT®. FLUENT employs a finite-volume method to solve the Navier-Stokes and continuity equations on an arbitrarily shaped flow domain with appropriate boundary conditions. The steady-state solution of the flow fields assumed converged when the residuals reduced to less than 10-4. A typical run time was approximately 38 hours on a normal processor. In this CFD study of surfactant in the human pulmonary tree model, computations were performed at 28.3 L/min for a quasi steady-state volumetric adult

The Volume of Fluid (VOF) model has been widely used in many research fields and commercial software. Known for its wide range of usage, implementation of this model in commercial code has been popular and FLUENT 6.2 is one of the software where their codes have been implemented. Isothermal VOF formulation was considered in this project with heat transfer and is assumed negligible due to analysis being based on surface interface. Formulation of laminar flow is presented here. For the turbulence effect, it can be done by averaging the Navier Stokes equation and detailed formulation can be found in Versteeg

is using a two fluids mixture in a fixed computational grid. The transport of the two fluids

reconstructed from the volume fraction field. The following procedures are used in solving

ii. From the existing interface in the model, volume of fraction is obtained. (Truncated volumes calculation was done at each interfacial cell, this discrete volume fraction is

iv. Finally, the interface geometry and location are obtained using the new volume

All these steps are repeated until a convergence criterion is achieved. In VOF, especially in FLUENT 6.2, one can conduct a simulation for two different fluids without mixing the respective fluids. There will be an interface region between the both liquids and for each

*<sup>q</sup>* < 1: the cell contains the interface between the qth fluid and one or more other fluids. Using the same momentum equation, both fluids shared the same velocity field and it is solved throughout their domain. The tracking interface between both phases is accomplished by the solution of a continuity equation for the volume fraction of one of the

where this volume of fraction

values and interface location is

inhalation flow rate for pulmonary surfactant.

and Malalsekera (1995).

the VOF model:

fraction.

*<sup>q</sup>* = 0: the cell is empty

*<sup>q</sup>* = 1: the cell is full

phases. The equation is:

region 

0 < 

**3.1 Volume of Fluid (VOF) method in Fluent 6.2** 

**3.1.1 The basic of VOF for interface development**  The VOF model is based on a volume fraction denoted by

results in solving the multiphase fluid equation using the

used instead of the exact interface information.) iii. The volume fraction field is adverted. (Transport of fluid.)

*<sup>q</sup>* of 1 and 0 values will be given.

i. Initial interface shape is modelled in computation domain.

$$\frac{\partial \alpha\_q}{\partial t} + \vec{v} . \nabla \alpha\_q = \frac{S\_{\alpha\_q}}{\rho\_q} \tag{1}$$

The primary phase of the fluid will not be solved by the volume fraction equation, but will be computed through the following constrains:

$$\sum\_{q=1}^{n} \alpha\_q = 1 \tag{2}$$

For each of the phases, the properties will be presented in control volume of the respective component phase. For the two phase system, for example, if the phases are represented by the subscripts 1 and 2, and if the volume fraction of these is being tracked, the density in each cell is given by:

$$
\rho = a\_2 \,\, \rho\_2 + (1 - a\_2)\rho\_1. \tag{3}
$$

and for every n – phase system , the volume fraction averaged density takes on the following form:

$$
\rho = \Sigma \alpha\_q \rho\_q \tag{4}
$$

#### **3.1.2 The momentum equation**

As mentioned earlier, both fluids will share the same momentum equation throughout the domain. The momentum equation is:

$$\frac{\partial}{\partial t}\rho u\_j + \rho u\_i \frac{\partial u\_j}{\partial u\_i} = -\frac{\partial P}{\partial \mathbf{x}\_j} + \mu \frac{\partial}{\partial \mathbf{x}\_i} \left(\frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i}\right) + \rho \mathbf{g}\_j + B\_{j,i} \tag{5}$$

Velocity component and direction denoted by *uj* and *<sup>j</sup> x* in j direction and t is the time. The summation of convention is applied when an index is repeated and i, j = 1, 2 for this problem. Pressure and gravitational acceleration are denoted by P and g respectively and B is a source term due to any additional force applied (per unit volume). This momentum equation dependant on the volume fraction *<sup>q</sup>* where in this case *q* = 1, 2 of the phases through the density and viscosity . The computational domain is calculated through the following equation:

$$\frac{D\alpha\_2}{Dt} = \frac{\partial \alpha\_2}{\partial t} + \mu\_i \frac{\partial \alpha\_2}{\partial x\_i} = 0 \tag{6}$$

The value of <sup>1</sup> for the first phase is not solved directly at any stage, but it is solved based on the volume condition 2 1 1. By applying the properties of density and viscosity , the following equations are obtained since and depend on the volume fraction values:

$$
\rho = a\_2 \rho\_2 + a\_1 \rho\_1 \tag{7}
$$

Surfactant Analysis of Thin Liquid Film in the

surface normal and the gradient is

the mean curvature.

by:

the cell.

can simplified to:

**3.2 Boundary conditions** 

hence, *p*

where ˆ *<sup>n</sup>*

In equation *p* is the pressure drop across the surface,

And the curvature represented by divergence of the unit normal, *n*ˆ

*vol ij*

*F*

And if only two phases are present in a cell, then

The liquid was considered Newtonian with viscosity

develop consistent laminar flows throughout the geometry.

*n*

Human Trachea via Application of Volume of Fluid (VOF) 455

The formulation of CSF model is used in surface curvature computation and it is located from the local gradient in the surface normal at the interface. For example, let *n* be the

> *n*

With this equation, the pressure variation across the interface is assumed linear and given

2 1 ( ) *<sup>q</sup> p p* 

*pairsij i j i j*

1 <sup>2</sup> , ( )

<sup>1</sup>

 

> 

 is the volume averaged density computed using Equation 4 . The simplified volume force shows that the surface tension source term for a cell is proportional to the average density in

flows in the region. The film thicknesses were obtained from the lipo lipid thickness in the human airway and it was set to be constant from the entrance of the tracheobronchial tree and uniform to its outlet. The membranes were patched with a similar thickness of the film, in order to reduce the computation period. A very fine mesh was created at the edge of the lung walls. The number of cells used was 198200 and a finer grid was placed near the thin film where the location of the interface was used. Two velocities inlets and a pressure outlet were used as boundary condition inputs. Inlet velocities were separated into two parts; small inlet and big inlet. The small inlet was selected for the thin film boundary condition and the flows maintained at 0.0001 m/s. Meanwhile the big inlet was set to 0.001 m/s to

*i i*

*i j*

, density

2

 

 *<sup>i</sup> <sup>j</sup>* and

*i i j j jj i i*

The force at the surface can be expressed as a volume force using divergence theorem:

 

*vol ij*

*F*

(11)

*<sup>q</sup>* (12)

(14)

(16)

and surface tension

*<sup>j</sup>* , this equation

.*n*ˆ

*<sup>n</sup>* (13)

(15)

*<sup>i</sup>* 

is the surface tension and

is

*<sup>q</sup>* . Then, the volume fraction will be:

$$
\mu = \alpha\_2 \mu\_2 + \alpha\_1 \mu\_1 \tag{8}
$$

The volume transport equation then, yields an updated volume fraction field with discrete values in computational cell. The values will later be used in reconstructing the new interface location.

#### **3.1.3 Interface reconstruction**

Reconstruction of new interface needs formulation from the discrete values of volume fraction inside the domain. In FLUENT 6.2, there are four interface reconstruction schemes, which are geometric reconstruction, donor-acceptor, Euler explicit and Euler implicit. Both Euler explicit and Euler implicit schemes treat cells that lie near the interface using the same interpolation as other cells. Meanwhile, the geometric reconstruction and the donor acceptor scheme treat the interfaces in an advanced interpolation. Full description of these schemes can be found in the FLUENT 6.2 User Manual (2002). A brief explanation on geometric reconstruction is explained here as it is recommended in FLUENT 6.2.

In the FLUENT 6.2 User Manual, the geometric reconstruction scheme approach standard interpolation to obtain face fluxes as the cell fills with one phase or another. Only at the cell near the interface between two phases is this method widely applied. This method is generalized for unstructured meshes and obtained from the work of Young (Piecewise linear 1982). It is assumed that for the calculation of the advection of fluid through the cell faces, the interface between two fluids has a linear slope within each cell. The following are the procedures of reconstruction on interface using geometric reconstruction scheme:


Aside from the steady state problem, the geometric reconstruction scheme also available for transient cases. User defined functions can be used in order to obtain the desired solution of interest.

#### **3.1.4 Effect of surface tension in volume of fraction**

In VOF, attraction between the molecules in a fluid causes the rises of surface tension effect. This surface tension is a force only acting on the surface. The surface tension model in FLUENT is the continuum surface force (CSF) model proposed by Brackbill et al. 1992. In this model, the addition of surface tension to the VOF calculation results in a source term in the momentum equation. The pressure drop across the surface depends upon the surface tension coefficient and the surface curvature as measured by two radii in an orthogonal direction, R1 and R2:

$$
\Delta p = p\_2 - p\_1 = \sigma \left( \frac{1}{R\_1} + \frac{1}{R\_2} \right) \tag{9}
$$

$$\mathbf{x} = \left(\frac{1}{R\_1} + \frac{1}{R\_2}\right),\tag{10}$$

where,

hence, *p*

454 Fluid Dynamics, Computational Modeling and Applications

The volume transport equation then, yields an updated volume fraction field with discrete values in computational cell. The values will later be used in reconstructing the new

Reconstruction of new interface needs formulation from the discrete values of volume fraction inside the domain. In FLUENT 6.2, there are four interface reconstruction schemes, which are geometric reconstruction, donor-acceptor, Euler explicit and Euler implicit. Both Euler explicit and Euler implicit schemes treat cells that lie near the interface using the same interpolation as other cells. Meanwhile, the geometric reconstruction and the donor acceptor scheme treat the interfaces in an advanced interpolation. Full description of these schemes can be found in the FLUENT 6.2 User Manual (2002). A brief explanation on geometric

In the FLUENT 6.2 User Manual, the geometric reconstruction scheme approach standard interpolation to obtain face fluxes as the cell fills with one phase or another. Only at the cell near the interface between two phases is this method widely applied. This method is generalized for unstructured meshes and obtained from the work of Young (Piecewise linear 1982). It is assumed that for the calculation of the advection of fluid through the cell faces, the interface between two fluids has a linear slope within each cell. The following are the procedures of reconstruction on interface using geometric reconstruction scheme: i. First, calculation of position of the linear interface relative to the centre or each partial

ii. Second, the advection amount of fluid through each face is calculated using computed interface representation. Normal and tangential velocity distribution on the faces is also

iii. Finally, each cell of volume fraction is calculated by computing the balance fluxes from

Aside from the steady state problem, the geometric reconstruction scheme also available for transient cases. User defined functions can be used in order to obtain the desired solution of

In VOF, attraction between the molecules in a fluid causes the rises of surface tension effect. This surface tension is a force only acting on the surface. The surface tension model in FLUENT is the continuum surface force (CSF) model proposed by Brackbill et al. 1992. In this model, the addition of surface tension to the VOF calculation results in a source term in the momentum equation. The pressure drop across the surface depends upon the surface tension coefficient and the surface curvature as measured by two radii in an orthogonal

2 1

*pp p R R* 

> 1 2 1 1 *R R*

 

1 2 1 1

(9)

, (10)

2 2 11 (8)

 

reconstruction is explained here as it is recommended in FLUENT 6.2.

cell is done based on values of volume fraction.

**3.1.4 Effect of surface tension in volume of fraction** 

used as information in calculating the advection of fluid.

interface location.

**3.1.3 Interface reconstruction** 

the previous step.

direction, R1 and R2:

interest.

where,

$$
\Delta p = \sigma \kappa \tag{11}
$$

In equation *p* is the pressure drop across the surface, is the surface tension and is the mean curvature.

The formulation of CSF model is used in surface curvature computation and it is located from the local gradient in the surface normal at the interface. For example, let *n* be the surface normal and the gradient is*<sup>q</sup>* . Then, the volume fraction will be:

$$m = \nabla \alpha\_q \tag{12}$$

And the curvature represented by divergence of the unit normal, *n*ˆ .*n*ˆ

*n*

$$\text{where}\tag{13}$$

$$\hat{n} = \frac{n}{|n|}\tag{13}$$

With this equation, the pressure variation across the interface is assumed linear and given by:

$$p\_2 = p\_1 + \sigma \kappa (\Delta \alpha\_q) \tag{14}$$

The force at the surface can be expressed as a volume force using divergence theorem:

$$F\_{\rm vol} = \sum\_{\rm pairs; i, i < j} \sigma\_{ij} \frac{\alpha\_i \rho\_i \kappa\_j \nabla \alpha\_j + \alpha\_j \rho\_j \kappa\_i \nabla \alpha\_i}{\frac{1}{2} (\rho\_i + \rho\_j)} \tag{15}$$

And if only two phases are present in a cell, then *<sup>i</sup> <sup>j</sup>* and *<sup>i</sup> <sup>j</sup>* , this equation can simplified to:

$$F\_{\rm vol} = \sigma\_{ij} \frac{\rho \kappa\_i \nabla \alpha\_i}{\frac{1}{2} (\rho\_i + \rho\_j)} \tag{16}$$

 is the volume averaged density computed using Equation 4 . The simplified volume force shows that the surface tension source term for a cell is proportional to the average density in the cell.

#### **3.2 Boundary conditions**

The liquid was considered Newtonian with viscosity , density and surface tension flows in the region. The film thicknesses were obtained from the lipo lipid thickness in the human airway and it was set to be constant from the entrance of the tracheobronchial tree and uniform to its outlet. The membranes were patched with a similar thickness of the film, in order to reduce the computation period. A very fine mesh was created at the edge of the lung walls. The number of cells used was 198200 and a finer grid was placed near the thin film where the location of the interface was used. Two velocities inlets and a pressure outlet were used as boundary condition inputs. Inlet velocities were separated into two parts; small inlet and big inlet. The small inlet was selected for the thin film boundary condition and the flows maintained at 0.0001 m/s. Meanwhile the big inlet was set to 0.001 m/s to develop consistent laminar flows throughout the geometry.

Surfactant Analysis of Thin Liquid Film in the

Fig. 6. Velocity profile of contour Figure 6

region taken approximately F=0.5 from the VOF calculation.

tension was neglected.

Human Trachea via Application of Volume of Fluid (VOF) 457

shown in Figure 5. The contours are from the calculation in which the effect of surface

Fig. 5. Contour of volume of fraction with surface tension neglected of airway mode

There are only small effects of interface changes shown for the simulation above. Near the edge of both membranes bordering adhesive to the human lung wall, the region bordering both fluids shows a minimum deformation over each fluid. Conversely there is some numerical diffusion involved where smeared effects can be seen at the interface resolution. In order to improve this factor, finer grids are able to solve this problem, but it occupies more computation time and it was not attempted here. This simulation was repeated by adding surface tension factor. The surface tension coefficient was taken to be 2 Nm as per lipo lipid configuration in the human body to obtain the effect of surface tension at the interface. The contours of volume of fraction are shown in Figure 7 based on the interface

Fig. 4. Patch region in the semi symmetrical airway model

These velocities were applied as the Reynolds number for thin films considered to be 1 . The density of water and viscosity of the lipo liquid was fixed to be 1023 kg/m3 and 0.00093 kg/ms and for air the density was 1.225 kg/m3 and viscosity 1.7894 x 10-5 kg/ms. The calculations were performed transiently using FLUENT 6.2 until the steady state system was reached. Initial the study was done to anticipate that the solution could be run in a steady state so as to avoid time consuming transient calculations, but for better results for the interface it was run under transient mode. A grid dependency study was done until it was found that the location of the interface was not changing with the grid refinement. The time step was set to be 0.001s and simulation was done until 6750 iteration for the solution to converge. Later the iteration was reduced to 4600 iteration as the convergence residual was achieved sooner than expected. After total time of 38 hours, all the simulations reached convergence state. In order to achieve the state of convergence, severe computational time step needed to be selected. Initially time step was set to be 0.00001s and it not only prolonged the state of convergence, but also required more computational time. An assumption was made to set the time step to be 0.001s after 2500 iterations with results not showing prominent converging state of residuals. It was important to select the correct values of velocity based on thin film liquid flow. As per the initial guesses done with different velocity magnitude, the results show the region of multiphase fluids washed by the greater velocity of air. Based on several recalculations, attempts were made to maintain the thin film velocity to be as low as 1 (Reynolds number) in order to achieve the desired profile.

#### **4. Results and discussion**

The simulations on lipo lipid with interface were done with two different selections of velocity intake based on nominal inhalation in the human body. To obtain more detailed studies on each of these cases, the effect of surface tension were also included later to the simulation to provide better understandings of the membrane function of the pulmonary surfactant (lipo lipid in the human airway). A contour plot of lipid air volume fraction is

These velocities were applied as the Reynolds number for thin films considered to be 1 . The density of water and viscosity of the lipo liquid was fixed to be 1023 kg/m3 and 0.00093 kg/ms and for air the density was 1.225 kg/m3 and viscosity 1.7894 x 10-5 kg/ms. The calculations were performed transiently using FLUENT 6.2 until the steady state system was reached. Initial the study was done to anticipate that the solution could be run in a steady state so as to avoid time consuming transient calculations, but for better results for the interface it was run under transient mode. A grid dependency study was done until it was found that the location of the interface was not changing with the grid refinement. The time step was set to be 0.001s and simulation was done until 6750 iteration for the solution to converge. Later the iteration was reduced to 4600 iteration as the convergence residual was achieved sooner than expected. After total time of 38 hours, all the simulations reached convergence state. In order to achieve the state of convergence, severe computational time step needed to be selected. Initially time step was set to be 0.00001s and it not only prolonged the state of convergence, but also required more computational time. An assumption was made to set the time step to be 0.001s after 2500 iterations with results not showing prominent converging state of residuals. It was important to select the correct values of velocity based on thin film liquid flow. As per the initial guesses done with different velocity magnitude, the results show the region of multiphase fluids washed by the greater velocity of air. Based on several recalculations, attempts were made to maintain the thin film velocity to be as low as 1 (Reynolds

The simulations on lipo lipid with interface were done with two different selections of velocity intake based on nominal inhalation in the human body. To obtain more detailed studies on each of these cases, the effect of surface tension were also included later to the simulation to provide better understandings of the membrane function of the pulmonary surfactant (lipo lipid in the human airway). A contour plot of lipid air volume fraction is

Fig. 4. Patch region in the semi symmetrical airway model

number) in order to achieve the desired profile.

**4. Results and discussion** 

shown in Figure 5. The contours are from the calculation in which the effect of surface tension was neglected.

Fig. 5. Contour of volume of fraction with surface tension neglected of airway mode

#### Fig. 6. Velocity profile of contour Figure 6

There are only small effects of interface changes shown for the simulation above. Near the edge of both membranes bordering adhesive to the human lung wall, the region bordering both fluids shows a minimum deformation over each fluid. Conversely there is some numerical diffusion involved where smeared effects can be seen at the interface resolution. In order to improve this factor, finer grids are able to solve this problem, but it occupies more computation time and it was not attempted here. This simulation was repeated by adding surface tension factor. The surface tension coefficient was taken to be 2 Nm as per lipo lipid configuration in the human body to obtain the effect of surface tension at the interface. The contours of volume of fraction are shown in Figure 7 based on the interface region taken approximately F=0.5 from the VOF calculation.

Surfactant Analysis of Thin Liquid Film in the

Fig. 9. Model with surface tension 2 Nm

Fig. 10. Schematic of the interface

corresponding value of y along coordinate distance x as:

differentiated to obtain the interface slope:

can be obtained.

Human Trachea via Application of Volume of Fluid (VOF) 459

Smooth surface

The iso-surface obtained above was computed in a polynomial condition to obtain a better understanding between the two simulations. It is understood that the two phases are separated by an interface and the normal velocity of the interface must denote the value of zero. From the results obtained, the interface region of the volume fraction varies from 1 (lipo lipid) to zero (air). It can be concluded that the interface can be located along this region. Starting from the volume fraction contour plot, such as the one shown in Figure 10, the coordinates and velocities for each point (defined by the grid) on the assumed interface

Considering a point on the interface, see Figure 10, the CFD calculated velocity components (Vx, Vy) correspond to the x, y values of the coordinate system. Normal and tangential velocities need to be calculated to evaluate the interface boundary conditions. An interface curve was fitted through sets of computed points corresponding to each of the volume fraction values investigated. The height of the interface can be expressed through the

This interface function of y(x) defined, which is a continuous function, can then be

'( ) *dy <sup>f</sup> <sup>x</sup>*

From the geometry, values of the unit normal vector can be obtained as follows:

*y f x*( ) (17)

*dx* (18)

Fig. 7. Interface section F=0.5 without surface tension, 0 Nm

In order identify the effect of the pulmonary surfactant, the edge perimeter of each simulation was selected to obtain the interface section for both, with and without surface tension. Figure 8 and Figure 9 will provide details of the edge perimeter membrane interface of pulmonary surfactant. The figure was studied from the entrance of the airway up until the first bifurcation section and the pulmonary effects were uniform throughout the model. As can be seen from the figure, the effect of surface tension plays a vital role in surface reconstruction. It is very prominent to see that as the value of surface tension increases from 0 to 2 (0 = no surface tension), the effect on the interface is to deliver a smoother curve on the edge of bifurcation and wall. It can be concluded that with surface tension effect, the volume fraction creates an even structure when surface tension increases. This proves the understanding or surface tension effect where the results show smoother effects as surface tension is applied. This also agreed with human lipo lipid membrane function where the existence of the membrane as thin liquid film provides a better breathing cycle and protect from injuries by preventing the collapse or absences of surfactant.

Fig. 8. Model without surface tension 0 Nm

Fig. 9. Model with surface tension 2 Nm

In order identify the effect of the pulmonary surfactant, the edge perimeter of each simulation was selected to obtain the interface section for both, with and without surface tension. Figure 8 and Figure 9 will provide details of the edge perimeter membrane interface of pulmonary surfactant. The figure was studied from the entrance of the airway up until the first bifurcation section and the pulmonary effects were uniform throughout the model. As can be seen from the figure, the effect of surface tension plays a vital role in surface reconstruction. It is very prominent to see that as the value of surface tension increases from 0 to 2 (0 = no surface tension), the effect on the interface is to deliver a smoother curve on the edge of bifurcation and wall. It can be concluded that with surface tension effect, the volume fraction creates an even structure when surface tension increases. This proves the understanding or surface tension effect where the results show smoother effects as surface tension is applied. This also agreed with human lipo lipid membrane function where the existence of the membrane as thin liquid film provides a better breathing cycle and protect from injuries by preventing the collapse or

Distored surface

Fig. 7. Interface section F=0.5 without surface tension, 0 Nm

absences of surfactant.

Fig. 8. Model without surface tension 0 Nm

The iso-surface obtained above was computed in a polynomial condition to obtain a better understanding between the two simulations. It is understood that the two phases are separated by an interface and the normal velocity of the interface must denote the value of zero. From the results obtained, the interface region of the volume fraction varies from 1 (lipo lipid) to zero (air). It can be concluded that the interface can be located along this region. Starting from the volume fraction contour plot, such as the one shown in Figure 10, the coordinates and velocities for each point (defined by the grid) on the assumed interface can be obtained.

Fig. 10. Schematic of the interface

Considering a point on the interface, see Figure 10, the CFD calculated velocity components (Vx, Vy) correspond to the x, y values of the coordinate system. Normal and tangential velocities need to be calculated to evaluate the interface boundary conditions. An interface curve was fitted through sets of computed points corresponding to each of the volume fraction values investigated. The height of the interface can be expressed through the corresponding value of y along coordinate distance x as:

$$y = f(\mathbf{x})\tag{17}$$

This interface function of y(x) defined, which is a continuous function, can then be differentiated to obtain the interface slope:

$$\frac{dy}{dx} = f'(x) \tag{18}$$

From the geometry, values of the unit normal vector can be obtained as follows:

Surfactant Analysis of Thin Liquid Film in the

**5. Conclusion** 

**6. Acknowledgment** 

attention and time.

**7. References** 

Human Trachea via Application of Volume of Fluid (VOF) 461

The comparison was done between the plot of normal velocity without surface tension and with surface tension. From Figure 11 the surface tension creates a smoother surface. The values of normal velocities of each case provide sufficient value of zero. It is understood that the interface region is zero normal velocity and with the application of the surface tension in the simulation it denotes as the numerical value. The results obtained clearly indicate that using VOF in FLUENT facilitate a good approach in determining the interface of new surface reconstruction. In the cases of using VOF method in the thin liquid film flow, this

In this chapter, a 3D model of the human tracheobronchial tree was investigated by using commercial software FLUENT 6.2. The major study was done to understand the behaviour of pulmonary surfactant in the human airway whereby the surface tension presences in airways are significant and without the existence of surfactant, it will severely affect the condition of inhalation. Apart from this, presenting the effects of surface tension in VOF provide significant results with application in thin liquid film studies. Finally, the use of the VOF model to calculate thin film was considered. It is a known fact that in a thin film flow, surface tension plays a vital role and is prominent in controlling the film surface. The results of VOF show that application of this method in the human airway with thin film flow accommodates a high level of accuracy. The conclusion was made that use of VOF in thin film flow is appropriate in identifying the real interface using volume fraction of VOF and

provides significant results on pulmonary surfactant condition in thin film flow.

the pulmonary surfactant in the human airway plays a major role for inhalation.

I have taken great effort in writing this book chapter, however, it would not have been possible without the kind support and help of many individuals and my organization. I would like to extend my sincere thanks to all of them. I am highly indebted to Bumi Armada Berhad and my family members for their guidance and constant motivation as well as for providing necessary support in completing this chapter. I would like to express my special gratitude and thanks to InTech and the University of Rijeka, Croatia, for giving me such

Allen G.M, Shortall B.P, Gemci T, Corcoran T.E, and Chigier N.A, (2004). Computational

Brackbill, J.U., Kothe, D.B. and Zemach, C. (1992). A continuum method for modelling

Caro C.G, Schroter R.C, Watkins N, Sherwin S.J, and Sauret V, (2002). Steady inspiratory

Cerne, G., Petelin, S. and Tiselj, I. (2000). Upgrade of the VOF method for the simulation of

surface tension, *Journal of Computational Physics* 100, pp 335–354

*Journal of Biomechanical Engineering* 126, pp 604–613

*the Royal Society*, pp 458:791 - 809.

*Meeting* Boston, Massachusetts

simulations of airflow in an in vitro model of the pediatric upper airways, *ASME* 

flow in planar and non planar models of human bronchial airways, *Proceedings of* 

the dispersed flow, *Proceedings of the ASME 2000, Fluid Engineering Division Summer* 

$$\hat{n}\_x = -\frac{dy/dx}{\sqrt{1 + \left(\frac{dy}{d\infty}\right)^2}}, \quad \hat{n}\_y = \frac{1}{\sqrt{1 + \left(\frac{dy}{d\infty}\right)^2}}\tag{19}$$

Similarly the unit tangential vectors are given by:

$$\hat{m\_y} = -\frac{dy/dx}{\sqrt{1 + \left(\frac{dy}{dx}\right)^2}}, \quad \hat{n\_x} = \frac{1}{\sqrt{1 + \left(\frac{dy}{dx}\right)^2}}\tag{20}$$

Having obtained these unit vectors, the normal and tangential velocity components of the interface are respectively given by:

$$V\_{\stackrel{\frown}{n}} = V\_{\text{x}} \stackrel{\frown}{n}\_{\text{x}} + V\_{\text{y}} \stackrel{\frown}{n}\_{\text{y}} \tag{21}$$

$$V\_{\stackrel{\frown}{t}} = V\_{\text{x}} \stackrel{\frown}{t}\_{\text{x}} + V\_{\text{y}} \stackrel{\frown}{\mathbf{n}}\_{\text{y}} \tag{22}$$

The calculations of first derivative *dy dx* ensure the accuracy of the result. Here *dy dx* represent the slope at any point x on the curve. Consequently, it is important to determine the precise interface curve representation *f x*( ) that passes through all data points.

Computation of normal velocity at the interface was done based on equation 21 and was repeated for both cases of volume of fraction on the conducted simulation. The objective of the investigation was to verify the accurateness of the obtained interface and its sensitivity to selection of significant volume fraction. The value of 0.5 of volume of fraction used here was based on previous work that shows at this particular value the interface can be presented more significantly. Using equation 21, the normal velocity along the x direction of the interface was calculated. The plot of normal velocity with respect to x direction is shown in Figure 11.

Fig. 11. Graph of normal velocity of iso -0.5 of both simulations

The comparison was done between the plot of normal velocity without surface tension and with surface tension. From Figure 11 the surface tension creates a smoother surface. The values of normal velocities of each case provide sufficient value of zero. It is understood that the interface region is zero normal velocity and with the application of the surface tension in the simulation it denotes as the numerical value. The results obtained clearly indicate that using VOF in FLUENT facilitate a good approach in determining the interface of new surface reconstruction. In the cases of using VOF method in the thin liquid film flow, this provides significant results on pulmonary surfactant condition in thin film flow.

## **5. Conclusion**

460 Fluid Dynamics, Computational Modeling and Applications

*ny*

*nx*

 

 

, <sup>2</sup> 1

*dy*

(19)

(20)

 

*dx*

*dy dx*

(21)

(22)

Surface Tension 2 Nm

 

1

, <sup>2</sup> 1

1

2

2

Having obtained these unit vectors, the normal and tangential velocity components of the

*x y x y <sup>n</sup> V Vn Vn* 

*x y x y <sup>t</sup> V Vt Vn* 

The calculations of first derivative *dy dx* ensure the accuracy of the result. Here *dy dx* represent the slope at any point x on the curve. Consequently, it is important to determine

Computation of normal velocity at the interface was done based on equation 21 and was repeated for both cases of volume of fraction on the conducted simulation. The objective of the investigation was to verify the accurateness of the obtained interface and its sensitivity to selection of significant volume fraction. The value of 0.5 of volume of fraction used here was based on previous work that shows at this particular value the interface can be presented more significantly. Using equation 21, the normal velocity along the x direction of the interface was calculated. The plot of normal velocity with respect to x direction is shown

0.00 0.20 0.40 0.60 0.80 1.00 1.20

Surface Tension 0 Nm

**x**

the precise interface curve representation *f x*( ) that passes through all data points.

Fig. 11. Graph of normal velocity of iso -0.5 of both simulations

*dy dx*

 

*dy*

 

*dx*

1

1

*dy dx <sup>n</sup>*

*dy dx <sup>n</sup>*

*x*

*y*

 

Similarly the unit tangential vectors are given by:

interface are respectively given by:

in Figure 11.

**Normal velocity (m/s)**


> In this chapter, a 3D model of the human tracheobronchial tree was investigated by using commercial software FLUENT 6.2. The major study was done to understand the behaviour of pulmonary surfactant in the human airway whereby the surface tension presences in airways are significant and without the existence of surfactant, it will severely affect the condition of inhalation. Apart from this, presenting the effects of surface tension in VOF provide significant results with application in thin liquid film studies. Finally, the use of the VOF model to calculate thin film was considered. It is a known fact that in a thin film flow, surface tension plays a vital role and is prominent in controlling the film surface. The results of VOF show that application of this method in the human airway with thin film flow accommodates a high level of accuracy. The conclusion was made that use of VOF in thin film flow is appropriate in identifying the real interface using volume fraction of VOF and the pulmonary surfactant in the human airway plays a major role for inhalation.

#### **6. Acknowledgment**

I have taken great effort in writing this book chapter, however, it would not have been possible without the kind support and help of many individuals and my organization. I would like to extend my sincere thanks to all of them. I am highly indebted to Bumi Armada Berhad and my family members for their guidance and constant motivation as well as for providing necessary support in completing this chapter. I would like to express my special gratitude and thanks to InTech and the University of Rijeka, Croatia, for giving me such attention and time.

## **7. References**


**20** 

*Japan* 

**3D Particle Simulations of** 

**in Micro-Capillary Vessel** 

*Kyushu Institute of Technology* 

*Hitachi Cooperation* 

Katsuya Nagayama1 and Keisuke Honda2

**Deformation of Red Blood Cells** 

With the increase in arteriosclerosis, thrombosis, etc., in order to find out the cause, research of the flow characteristic of blood attracts attention. As for the analysis of the flow phenomenon of the RBC (Red Blood Cell or Erythrocyte), the numerical simulation (Wada et al., 2000, Tanaka et al., 2004) as well as experiment observation (Gaehtgens et al., 1980) is becoming a strong tool. Particle methods, such as SPH method (Monaghan J., 1992) and the MPS method (Koshizuka, 1997), treats both solid and liquid as particles, and can be applied to complicated flow analysis. When applying a particle method to the flow analysis of RBC, RBC is divided into the elastic film and internal liquid, and its deformation was analyzed in

The RBC which is actually flowing in our body occupies 40-60% by volume ratio of blood (hematocrit), and is numerous. The objective of our research is clarifying the flow characteristic of the blood flow containing many RBCs. We reported preliminarily simulation of 2D blood flow (Nagayama et al., 2004), where many RBCs were simply treated as a lump of an elastic particle, the flow was analyzed qualitatively. Moreover, three dimensional RBC was modelled with double structure, and the RBC shape in flow was more realistic (Nagayama et al., 2005). The relation of the blood vessel diameter and the bloodflow with many RBCs was studied by 2D model (Nagayama, 2006) and by 3D model (Nagayama et al., 2008a). The model was also applied to 3D blood flow in capillary bend

The objective is to understand the fundamental flow phenomenon in a blood vessel. In this

In Section 2, simulation model was described. And the shape of single red blood cell in static

In Section 3, blood flows with RBCs in capillary straight tube were simulated. And the relations of the blood vessel diameter and the hematocrit were investigated. Furthermore, transient phenomena of interacting red blood cells and their shape were investigated. In Section 4, the model is applied to the capillary vessel flow at finger tip edge. The capillary vessel is modelled as two cases. One case is bent tube and another is bent and twisted tube,

paper, 3 dimensional blood flows with RBCs in capillary tube were simulated.

**1. Introduction** 

detail (Tanaka et al., 2004, Tsubota et al., 2006).

tube (Nagayama et al., 2008b).

and RBC deformation were investigated.

fluid was shown.

Farkas A, and Balashazy I, (2007). Simulation of the effect of local obstructions and blockage on airflow and aerosol deposition in central human airways, *Journal of Aerosol Science* 38: pp 865–884

FLUENT 6.2 User Manual, (2002)

