+,*''".*+",0+,*''".*+",0
!(%*+!")+,*'!(%*+!")+,*'
!'"%',*"%+' "'*"'
-%","('+
!'"%',*"%+' "'*"'
)*,&',
('#- ,1-")(*(-+))*(!(*+"&-%,"' "*1(/"'('#- ,1-")(*(-+))*(!(*+"&-%,"' "*1(/"'
*%"+," (&,*"*)*+',,"('+(,!!-&'*+)"*,(*0*%"+," (&,*"*)*+',,"('+(,!!-&'*+)"*,(*0
+0+,&+0+,&
!*"+,()!**((,*
+,*''".*+",0
',!('0,*,&'
-/(
+,*,-/(
(%%(/,!"+'","('%/(*$+,!,,)+"*%"-/(&!'"%)-
*,(,!,*"%+"''' "'*"' (&&('+',!!'"%' "'*"' (&&('+
",,"('(,!"+))*",,"('(,!"+))*
*((,!*"+,()!**',*,&'',!('0('#- ,1-")(*(-+))*(!(*+"&-%,"'
"*1(/"'*%"+," (&,*"*)*+',,"('+(,!!-&'*+)"*,(*0+0+,&
!'"%'
,*"%+' "'*"' -%","('+
!,,)+"*%"-/(&!'"%)-
A Conjugate Fluid–Porous Approach
for Simulating Airflow in Realistic
Geometric Representations of the
Human Respiratory System
Christopher T. DeGroot
1
Mem. ASME
Department of Mechanical and Materials Engineering,
Western University,
London, ON N6A 5B9, Canada
e-mail: cdegroo5@uwo.ca
Anthony G. Straatman
Mem. ASME
Department of Mechanical and Materials Engineering,
Western University,
London, ON N6A 5B9, Canada
e-mail: astraatman@eng.uwo.ca
Simulation of flow in the human lung is of great practical interest
as a means to study the detailed flow patterns within the airways
for many physiological applications. While computational simula-
tion techniques are quite mature, lung simulations are particu-
larly complicated due to the vast separation of length scales
between upper airways and alveoli. Many past studies have pre-
sented numerical results for truncated airway trees, however,
there are significant difficulties in connecting such results with re-
spiratory airway models. This article presents a new modeling
paradigm for flow in the full lung, based on a conjugate
fluid–porous formulation where the upper airway is considered as
a fluid region with the remainder of the lung being considered as
a coupled porous region. Results are presented for a realistic lung
geometry obtained from computed tomography (CT) images,
which show the method’s potential as being more efficient and
practical than attempting to directly simulate flow in the full lung.
[DOI: 10.1115/1.4032113]
1 Introduction
The ability to numerically simulate flow in the human lung is of
great interest to the medical community because of the potential
advancements in respiratory drug delivery, particle deposition,
etc., that can be attained with detailed knowledge of the flow pat-
terns within the lung [1]. It is not, however, a simple task to simu-
late flow in the whole lung because of the many orders of
magnitude separating the length scales of the upper airways and
the alveoli. The air-filled portion of the lung can be divided into
two main regions: (i) the conducting region, which consists of a
network of bifurcating airways that transport air to the respiratory
region, and (ii) the respiratory region, where gas exchange occurs
by passive diffusion through the thin walls of the alveolar sacs
which surround the respiratory airways [2]. At the trachea, the first
airway generation in the conducting region, the diameter is of the
order of centimeters; at the alveolar sacs, the terminal generation
in the respiratory region, the diameter is of the order of a fraction
of a millimeter [2]. In addition to the large separation of length
scales, there are approximately 300 10
6
alveolar sacs in the
human lung [2], which cannot possibly be considered individually
when simulating processes in the lung using the computational
technology of today or even the foreseeable future.
Computational fluid dynamics (CFD) simulations have been
conducted in both the upper airways [1,3–19] and the alveolated
ducts within the respiratory region [20–31], however, there is
great difficulty in connecting the results at the two levels in order
to simulate processes in the whole lung. As computational power
has increased over the years, the total number of airway segments
considered in upper airway simulations has increased from seven
[3] to over 1400 [11] (although the computational mesh was not
sufficiently refined in this case to have grid-independent results).
That being said, the total number of airway segments in a full 23-
generation airway tree would be more than 1.6 10
7
and would
require hundreds of billions of mesh elements to adequately dis-
cretize the domain for simulation [15,18]. Since simulations of
this size appear to be out of reach for quite some time, other inno-
vative approaches need to be developed to capture the fundamen-
tal information of a full lung simulation model with less
computational effort.
Several innovative methods for simulating more branches of
the airway tree with less effort have been proposed including
those that simulate small subsections of the airway tree and use
the outlet condition of one subunit as the input to the next subunit
[7,10], the use of partially resolved airway trees [11,15], and the
coupling of three-dimensional CFD models for the upper airways
with one-dimensional resistance models for the lower airways
[1,14]. In addition to saving computational time, the use of ap-
proximate models for smaller airways is quite pragmatic given
that the geometric models at such small scales are approximate as
well. As a result, it may be argued that from a practical standpoint
it is a better use of resources to solve a simplified model in the
small airways rather than solving the flow in an approximate
geometry with high accuracy.
It has been proposed by Owen and Lewis [32], in their theoreti-
cal work on high-frequency ventilation, that the lung parenchyma
can be modeled as a porous continuum. Using homogenization
and volume-averaging, they developed a model to describe the
flow and tissue deformation for small uniform samples of lung tis-
sue. While the theoretical development of their model was well-
founded, it depended on several effective properties of the porous
lung parenchyma that were only roughly estimated and only select
one-dimensional results were obtained for the flow and tissue
deformations. Lande and Mitzner [33] used the equations of
Owen and Lewis [32] to develop a model for lung impedance
based on parameters describing the properties of the lung.
In light of the preceding discussion, it is the goal of this work to
propose a new method for simulating flow in the full lung by treat-
ing it as a conjugate fluid–porous domain. We hypothesize that
the large upper airways can be treated as a fluid region and the
remaining airways and alveoli can be treated as a connected,
coupled porous region where the flow is driven by the moving
boundary of the lung. The focus of this work is on the develop-
ment and implementation of a robust mathematical model of a
conjugate lung model based on a realistic representative lung ge-
ometry, as well as a demonstration of the potential applications of
the model. In future studies, it is envisioned that the proposed
model could be used for more comprehensive parameter studies,
alongside experimental measurements.
In this article, the governing equations and the associated nu-
merical methods will first be outlined, followed by the setup of
the problem, including the method for obtaining the lung geome-
try from CT images, details of creating a meshed lung model, pa-
rameter estimation, and boundary conditions. Finally, results are
presented using the coupled fluid–porous model of the lung solved
using the proposed numerical method. Ultimately, it is shown that
the proposed modeling paradigm has potential as an efficient tool
for simulating processes in the full human lung.
2 Methods
2.1 Governing Equations and Numerical Methods. It is
proposed that the lung may be considered as a conjugate
1
Corresponding author.
Manuscript received January 28, 2015; final manuscript received November 27,
2015; published online January 29, 2016. Assoc. Edito r: Naomi Chesler.
Journal of Biomechanical Engineering MARCH 2016, Vol. 138 / 034501-1Copyright
V
C
2016 by ASME
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/29/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use
fluid–porous domain, where the fluid region is a truncated airway
tree and the porous region is the remainder of the lung. In the fluid
region, the flow is governed by the continuity and Navier–Stokes
equations for an incompressible fluid, given as
ru ¼ 0 (1)
and
q
f
@u
@t
þr uu
ðÞ
¼rp þ l
f
r
2
u (2)
respectively. In the porous region, flow is governed by the
volume-averaged counterparts of Eqs. (1) and (2), given as
rhui¼
@e
@t
(3)
and
q
f
@hui
@t
þr
1
e
huihui
¼erhpi
f
þ l
f
r
2
hui
el
f
K
hui
(4)
respectively, where a detailed derivation of Eqs. (3) and (4) can
be found in Ref. [34]. For the purposes of this work, it is assumed
that the orientation of the ducts is random such that the flow has
no preferred direction within the parenchyma, such that the per-
meability can be taken as a scalar. The Reynolds number is taken
to be sufficiently low that non-Darcy effects are negligible, which
is certainly reasonable in the respiratory airways. To further sim-
plify the governing equations, spatial variation of porosity is not
included, although its addition is relatively straightforward to
implement numerically after expansion of the convection term by
the product rule. The temporal variation in porosity is retained in
the governing equations, although porosity is taken as being con-
stant in time for the case considered herein.
At the interface between the fluid and porous regions, the equa-
tions are coupled by appropriate conditions which ensure a bal-
ance of fluxes, as well as a balance of viscous and pressure forces.
The basic numerical approach used in this work follows that pre-
sented previously by the authors in Ref. [35], which has been
shown to be an effective method of coupling fluid and porous
regions where general unstructured grids are required to discretize
the domain of interest.
Since breathing is driven by the motion of the boundary of the
lung, the computational mesh must be considered to be in motion.
As a result, some extensions to the original numerical method pre-
sented in Ref. [35] must be made in order to reformulate the
method in an arbitrary Lagrangian–Eulerian (ALE) framework,
which is briefly described below. Since the fluid is incompressi-
ble, the mesh motion does not impact the continuity equations
[36], so only the momentum equations will be discussed in light
of ALE framework. Rewriting the momentum equations in inte-
gral form for an arbitrary control volume in space, X, which is
bounded by the moving control surface @X, results in [36,37]
@
@t
ð
X
q
f
udV þ
ð
@X
q
f
uu u
s
ðÞ
ndS ¼
ð
@X
pndS þ
ð
@X
l
f
ru n dS
(5)
and
@
@t
ð
X
q
f
huidV þ
ð
@X
q
f
e
huihuiu
s
ðÞ
ndS ¼
ð
@X
ehpi
f
ndS
þ
ð
@X
l
f
rhuindS
ð
X
el
f
K
huidV (6)
where the unit-normal vector to the surface @X is denoted n, and
u
s
is the velocity of the control surface, which is used to account
for mesh motion.
It is important when calculating flows on moving meshes that
the surface velocity, u
s
, is selected in such a way that volume is
conserved in order to avoid artificial mass sources in the domain.
This concept is expressed through the “geometric conservation
law” [36–40], which is taken into account in the discretization of
the governing equations for this work. Discretization of the gov-
erning equations is carried out using a spatially second-order
finite-volume method, as described in Ref. [35], with the appropri-
ate modifications to the transient and convection terms to account
for the changing cell volume and control surface velocity [36,37].
In addition to discretization of the governing equations, a pro-
cedure for updating the locations of the mesh nodes, based on pre-
scribed boundary motions, is required. In this work, the motion of
the mesh nodes is determined by numerically solving Laplace’s
equation with a variable diffusion coefficient and Dirichlet condi-
tions on all domain boundaries [41], which is reasonably robust
for large deformations, is relatively straightforward to implement,
and is not overly expensive to compute. The Laplace equation is
solved using a cell-centered finite-volume method, similar to that
used to solve all other transport equations in this work, where the
mesh stiffness coefficient is taken to be inversely proportional to
the cell volume such that larger cells absorb more of the motion
and smaller cells move more like rigid bodies.
2.2 CT Image Segmentation. Low-dose whole lung CT
images obtained from a single breath hold with a 1.25 mm slice
thickness were obtained from the ELCAP Public Lung Image
Database.
2
From these images, a three-dimensional representation
of the airway tree, up to a certain bifurcation, and a representation
of the remaining lung volume were obtained. The truncated air-
way tree was extracted using the segmentation software
ITK-SNAP
[42], however, only a limited number of bifurcations were cap-
tured to prevent the segmented volume from leaking into regions
of the lung outside of the airways. The resulting segmentation was
exported as a stereolithography (STL) file and was smoothed and
decimated, to reduce the total number of facets, using
MESHLAB.
3
The remaining lung volume was segmented by first using the
automatic threshold segmentation feature in
OSIRIX [43], however,
due to limitations in the segmentation algorithm, it was not possi-
ble to exclude the airway tree from the segmented volume. Thus,
the airways were manually removed from the segmented volume
and any geometric problems were manually repaired. The seg-
mentation obtained in
OSIRIX was then exported as a new image se-
ries and segmented using
ITK-SNAP [42], which produced
significantly smoother surfaces than
OSIRIX. Again, MESHLAB was
used to smooth and decimate the exported faceted surface.
2.3 Geometric Model and Meshing. After segmentation, the
STL surfaces generated for the airway tree and the remaining lung
volume were imported into
ANSYS
V
R
ICEM CFD, Release 13.0, to
be combined, further cleaned up, and meshed for CFD simula-
tions. The airway and lung volume surfaces were first intersected
and the portion of the lung surface penetrated by the airway tree
was removed. Additionally, the ends of the airway tree branches
and trachea were modified such that they formed flat surfaces.
The resulting surface model is shown in Fig. 1. From this figure, it
is obvious that the CT images did not include details of the mouth
and larynx. If one were to conduct simulations of purely inspira-
tory flow in a truncated airway tree, this configuration would not
necessarily pose a problem since a standard inlet boundary condi-
tion could be applied where the trachea is truncated. However, in
this case, a full breathing cycle with both inspiratory and expira-
tory flow is to be considered and the flow is to be driven by the
motion of the boundary. Therefore, the opening where the trachea
is truncated must allow both inflow and outflow, possibly simulta-
neously when the bulk flow is changing direction. As a result, the
2
http://www.via.cornell.edu/databases/lungdb.html
3
http://meshlab.sourceforge.net/
034501-2 / Vol. 138, MARCH 2016 Transactions of the ASME
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/29/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use
trachea was extended somewhat to allow the flow to develop more
fully before exiting the domain, and a hemispherical cap was
added to mimic expansion to an open environment. Treatment of
this boundary will be discussed further in Sec. 2.6.
The resulting geometric model was meshed in
ANSYS
V
R
ICEM
CFD, Release 13.0, using tetrahedral volumes that were refined
near fluid–porous interfaces as well as at the airway walls and the
lung boundary. The meshed model is shown in Figs. 2(a) and 2(b)
which show the main lung volume and the trachea extension,
respectively.
2.4 Grid Motion. To provide boundary conditions for the
Laplace equation governing the mesh motion, an estimate of the
breathing rate and the nature of the lung motion are required. For
the test cases considered, a breathing rate of 12 breaths per minute
was chosen, which is within the normal range of breathing rates
[44]. The motion of the lung, which drives breathing, is a result of
the motion of the diaphragm, a thin sheet of muscle located along
the bottom surface of the lung [44]. Although there is also motion
of the ribcage, it has been shown that the magnitude of the motion
of the diaphragm is somewhat larger [45]. Thus, for simplicity, it
is assumed for the purposes of this work that all motions are
driven by the diaphragm and that the magnitude of its motion is
1.5 cm, consistent with normal breathing [44]. It is found that this
specification of the lung motion results in a tidal volume of
approximately 400 mL, which is slightly smaller than the typical
tidal volume of 500 mL given by West [44], thus represents light
breathing. It should be noted at this point that the lung motion is
simply a boundary condition to the mesh motion algorithm and is
thus easily modified if one wishes to explore more elaborate spec-
ifications of its motion to reflect a specific breathing pattern.
Although it has been stated that the mesh motion at the dia-
phragm should be 1.5 cm, the mesh motion must be specified on
all boundaries in order to solve for the motion of the interior
nodes. In this case, the motion of all boundaries nodes is taken to
be in the z-direction, as shown in Fig. 1. The magnitude of the dis-
placement should be 1.5 cm at the base of the lung, whereas at the
top of the lung it should be zero. Thus, an appropriate variation
between these two values must be specified. To avoid motion of
the truncated airways, the motion is blended between a value of
zero and 1.5 cm in the region below the truncated airways.
2.5 Parameter Estimation. According to Weibel’s idealized
“model A” lung geometry, more than three quarters of the volume
of air in the lung resides in generations 20–23, which are fully
alveolated ducts or alveolar sacs [2]. Further, almost 90% of the
air resides in generations 17–23 which have at least some degree
of alveolation [2]. Thus, a logical approximation for the perme-
ability of the lung parenchyma is the permeability of an alveolated
duct, which has been determined numerically in Ref. [46]. Taking
a volume-weighted average of the airway diameters in generations
17–23 given by Weibel [2] yields an average diameter of
0.43 mm. Using this average diameter and the result for the
dimensionless permeability given in Ref. [46] yields an estimate
for the average permeability of the lung parenchyma of K/
e ¼ 1.75 10
9
m
2
, where the porosity remains to be determined.
According to the measurements of Gehr et al. [47], the lung is
composed of 86.5% air, 5.7% capillary blood, and 7.8% tissue.
The same study also found that this composition does not change
significantly in space, justifying the previous assumption of zero
spatial porosity gradient. Results of Kampschulte et al. [48] indi-
cate significant variations in air volume fractions between sub-
jects, with measurements ranging from approximately 81–93%.
Therefore, a rough estimate of e ¼ 0.9 is used for the purposes of
this work and is taken as constant in time for simplicity, noting
that extension to a time-varying porosity is included in the
formulation.
Fig. 2 Plots of the computational meshes for the lung geome-
try showing (a) the main lung mesh and (b) the trachea
extension
Fig. 1 An illustration of the combined airway tree (inner struc-
ture) and the lung surfaces (remaining structure)
Journal of Biomechanical Engineering MARCH 2016, Vol. 138 / 034501-3
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/29/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use
Basedonthevaluee ¼ 0.9, the permeability is K ¼ 1.58 10
9
m
2
.
Note that this permeability estimate is for a representative airway
diameter, which will change with expansion and contraction of
the lung. For the purposes of this work, it will be assumed that the
hydraulic diameter of the duct varies sinusoidally with time,
according to
D
H
¼ D
H;0
½1 þ K
D
H
cosðxtÞ (7)
where D
H,0
is the hydraulic diameter at mean inflation, and K
DH
is
the dimensionless amplitude of the change in hydraulic diameter.
This leads to the permeability as a function of time, given as
K ¼ K
0
½1 þ K
D
H
cosðxtÞ
2
(8)
where K
0
is the permeability at mean inflation, i.e.,
K
0
¼ 1.58 10
9
m
2
.
According to the studies by Sznitman et al. [24] and Harding
and Robinson [28], the volume of an alveolated duct changes by
approximately 15% from minimum to maximum inflation, indi-
cating that the length scale of the duct changes by about 5%.
Accordingly, the amplitude of the change in duct diameter is
taken to be K
D
H
¼ 0:025, such that D
H
varies between 0.975D
H
and 1.025D
H
.
It should be noted that the parameter estimates described in this
section are meant simply to provide the model with data that are
of approximately the right order of magnitude in order to assess
the reasonableness of the final results. For more detailed studies,
or those specific to a particular patient, more accurate values for
these parameters should be obtained.
2.6 Boundary Conditions. The walls of the airway tree are
considered as fixed no-slip, impenetrable walls, such that all ve-
locity components may be set to zero and the pressure is extrapo-
lated from the interior of the domain to the boundary. Note that
the walls of the extended trachea and the flat lower surface of the
hemispherical cap are treated in the same way as fixed walls. The
lung boundary is considered to be a moving no-slip, impenetrable
wall and the velocity components are selected such that there is
no mass flux through any of the control surfaces on the moving
walls. On the moving walls, the pressure is also extrapolated from
the interior of the domain.
At the open boundary, i.e., the curved portion of the hemispher-
ical cap, the derivatives of all velocity components normal to the
boundary are set to zero. When the fluid is exiting the domain, a
constant static pressure is specified, whereas when the fluid enters
the domain, a constant dynamic pressure is specified, as recom-
mended by Mathur and Murthy [49] for open boundaries.
3 Results and Discussion
3.1 Scope. The model for air flow in the human lung
described in Sec. 2 was run for five full breath cycles to ensure
that any initial transient behavior was dissipated and results were
extracted from the final full breath cycle that was computed. For
the purposes of the results to be presented, the time t ¼ 0 corre-
sponds to the beginning of the final breath cycle. The computa-
tional mesh shown previously in Fig. 2 contained 911,156
tetrahedral control volumes. In this section, results will be pre-
sented for the maximum pressure difference across the domain,
representing the pressure drop from the trachea to the most distal
alveolus, as well as pressure contours on the surface of the lung
and velocity vectors in a plane intersecting the bifurcation from
the trachea to the main bronchi. Although detailed experimental
Fig. 3 A plot of the maximum pressure difference within the
domain as a function of time, where t 5 0 represents the begin-
ning of the breath cycle at maximum inflation
Fig. 4 Contour plots of the pressure in Pascals for the times:
(a) 0.75, (b) 1.25, (c) 1.75, (d) 2.5, (e) 3.25, (f) 3.75, (g) 4.25, and
(h) 5.0 s from the beginning of the breath cycle
034501-4 / Vol. 138, MARCH 2016 Transactions of the ASME
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 01/29/2016 Terms of Use: http://www.asme.org/about-asme/terms-of-use