UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Numerical simulation of the air flow and particulate deposition in emphysematous human acini Dutta, Amitvikram 2016

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2017_february_Dutta_Amitvikram.pdf [ 4.33MB ]
JSON: 24-1.0340526.json
JSON-LD: 24-1.0340526-ld.json
RDF/XML (Pretty): 24-1.0340526-rdf.xml
RDF/JSON: 24-1.0340526-rdf.json
Turtle: 24-1.0340526-turtle.txt
N-Triples: 24-1.0340526-rdf-ntriples.txt
Original Record: 24-1.0340526-source.json
Full Text

Full Text

Numerical Simulation of the Air Flowand Particulate Deposition inEmphysematous Human AcinibyAmitvikram DuttaB.E. Hons., BITS Pilani, 2014A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinTHE COLLEGE OF GRADUATE STUDIES(Mechanical Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan Campus)December 2016c© Amitvikram Dutta, 2016The undersigned certify that they have read, and recommend to theCollege of Graduate Studies for acceptance, a thesis entitled: NumericalSimulation of the Air Flow and Particulate Deposition in Em-physematous Human Acini submitted by Amitvikram Dutta in partialfulfilment of the requirements of the degree of Master of Applied ScienceDr. Joshua Brinkerhoff, Applied Science/School of EngineeringSupervisor, Professor (please print name and faculty/school above the line)Dr. Andre Phillion, Materials Science and Engineering, McMaster UniversityCo-Supervisor, Professor (please print name and faculty/school above the line)Dr. Mina Hoorfar, Applied Science/School of EngineeringSupervisory Committee Member, Professor (please print name and faculty/school abovethe line)Dr. Liwei Wang, Applied Science/School of EngineeringUniversity Examiner, Professor (please print name and faculty/school above the line)Dr. Neil Eves, Health and Exercise ScienceExternal Examiner, Professor (please print name and faculty/school above the line)12/16/2016(Date Submitted to Grad Studies)iiAbstractEmphysema, is a destructive process that leads to the permanent en-largement of air spaces within the parenchyma of the lung. Along withchronic bronchitis, emphysema forms one of the two components of ChronicObstructive Pulmonary Disease (COPD), a serious condition that is respon-sible for severe limitation of expiratory airflow in its victims. The earlystages of emphysema are charecterized by the destruction of tissue in thepulmonary acinus - the part of the lung airway tree responsible of gas ex-change with the bloodstream. Little is known how emphysema affects airflowwithin the acinus especially in the early stages of the disease. In this the-sis computational fluid dynamics simulations are performed of airflow in amathematically-derived model of a section of the pulmonary acinus. Thecomputational domain consists of two generations of the acinus with alve-olar geometries approximated as closely-packed, fourteen-sided polygons.Physiologically realistic flow rates and wall motions are used to capture theacinar flow during the inspiratory and expiratory phases of the breathingcycle. The effects of emphysema on the airway wall motion, flow rates, andseptal destruction are simulated at various stages of the disease’s progressionto identify the effect on the flow in the acinar region. Parametric studies arepresented to independently assess the relative influence of septal destructioniiiAbstractand the emphysematous degradation of airway motion and flow rates. Theresults illustrate that septal destruction lowers the flow resistance throughthe alveolar ducts but has little influence on the mass transport of oxygeninto the alveoli. Septal destruction has a net effect on the flow field byfavouring the development of recirculatory flow patterns in individual alve-oli. The effects of the gradually advancing emphysema on the deposition ofmicron-sized particles in the acinus are also studied. The simulations arecategorized according to particle size and the relative orientation of the grav-itational vector to the incoming flow. Emphysematous destruction increasesthe deposition of particles in affected ducts, with the greatest increase occur-ring for the larger particle size when the gravity vector is oriented tangentialto the incoming flow.ivPrefaceThe work outlined in this thesis was conducted at UBC’s Okanagan CFDLab by Amitvikram Dutta under the supervision of Dr. Joshua Brinkerhoffand Dr. Andre Phillion. Frequent discussions on the physiological aspects ofthe simulations were conducted with Dr. Dragos Vasilescu, Dr. Tillie Hack-ett and Dr. James C. Hogg of the Heart and Lung Institute (HLI) in UBCVancouver. I was responsible for designing and carrying out the simulations,while Dr. Brinkerhoff and Dr. Phillion helped with the interpretation of theresults.Portions of Chapter 2, 4 and 5 have been presented at a conference.Dutta A., Vasilescu D.M., Hogg J.C., Phillion A.B., Brinkerhoff J. (2016June). Simulation of Airflow in Idealized Emphysematous Human Aci-nus.24th Annual Conference of the CFD Society of Canada, Kelowna.vTable of ContentsThesis Committee . . . . . . . . . . . . . . . . . . . . . . . . . . iiAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . .xviiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .xviiiChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2.1 Human lung anatomy and physiology . . . . . . . . . 21.2.2 Chronic obstructive pulmonary disease (COPD) . . . 21.2.3 Diagnosis of COPD . . . . . . . . . . . . . . . . . . . 4viTABLE OF CONTENTS1.2.4 Emphysema . . . . . . . . . . . . . . . . . . . . . . . . 61.3 Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.4 Overview of the thesis . . . . . . . . . . . . . . . . . . . . . . 9Chapter 2: Literature Review . . . . . . . . . . . . . . . . . . . 102.1 Anatomy of the human lung . . . . . . . . . . . . . . . . . . . 102.1.1 Airway structure . . . . . . . . . . . . . . . . . . . . . 102.1.2 Mechanism of pulmonary ventilation . . . . . . . . . . 122.1.3 Effects of emphysema on mechanical behavior of lungairways . . . . . . . . . . . . . . . . . . . . . . . . . . 142.1.4 The pulmonary acinus . . . . . . . . . . . . . . . . . . 172.2 Numerical simulations of alveolar airflows . . . . . . . . . . . 202.2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . 202.2.2 Governing equations . . . . . . . . . . . . . . . . . . . 212.2.3 The nature of alveolar flow . . . . . . . . . . . . . . . 232.2.4 Numerical simulations of particulate transport in alveoli 262.3 CFD simulation of diseased acini . . . . . . . . . . . . . . . . 312.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32Chapter 3: Scope and Objectives . . . . . . . . . . . . . . . . . 34Chapter 4: Methodology . . . . . . . . . . . . . . . . . . . . . . 364.1 Computational domain . . . . . . . . . . . . . . . . . . . . . . 364.2 Boundary conditions and wall motion . . . . . . . . . . . . . 384.3 Spatial mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.4 Solution approach . . . . . . . . . . . . . . . . . . . . . . . . 46viiTABLE OF CONTENTS4.4.1 Continuity and momentum equations . . . . . . . . . 464.4.2 Finite volume method . . . . . . . . . . . . . . . . . . 464.4.3 Scalar transport . . . . . . . . . . . . . . . . . . . . . 504.4.4 Particle transport . . . . . . . . . . . . . . . . . . . . 534.4.5 Initial conditions . . . . . . . . . . . . . . . . . . . . . 544.4.6 Solution hardware . . . . . . . . . . . . . . . . . . . . 54Chapter 5: Results and Discussion . . . . . . . . . . . . . . . . 555.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 555.2 Flow field and oxygen transport . . . . . . . . . . . . . . . . . 555.2.1 Oxygen transport . . . . . . . . . . . . . . . . . . . . . 555.2.2 Hydraulic losses . . . . . . . . . . . . . . . . . . . . . 575.2.3 Flow field visualization . . . . . . . . . . . . . . . . . . 665.3 Particle transport . . . . . . . . . . . . . . . . . . . . . . . . . 715.3.1 Simulation conditions and boundary conditions . . . . 715.3.2 Temporal variation of deposition fraction . . . . . . . 725.3.3 Variation of steady state deposition fraction . . . . . . 79Chapter 6: Summary and Conclusions . . . . . . . . . . . . . . 856.1 Model development . . . . . . . . . . . . . . . . . . . . . . . . 856.2 Oxygen transport and flow-field . . . . . . . . . . . . . . . . . 866.3 Particle transport . . . . . . . . . . . . . . . . . . . . . . . . . 876.4 Summary of conclusions . . . . . . . . . . . . . . . . . . . . . 896.5 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91viiiTABLE OF CONTENTSAppendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102ixList of TablesTable 4.1 Variation of parameters governing flow-rate and air-flow motion due to emphysema . . . . . . . . . . . . . 40Table 4.2 Mesh independence study, with the observed variablebeing the static pressure drop across Duct I . . . . . . 45Table 5.1 Simulation matrix for particle deposition study . . . . 72xList of FiguresFigure 1.1 The pulmonary airways. Adapted with permissionfrom [1]. . . . . . . . . . . . . . . . . . . . . . . . . . 3Figure 1.2 Spirometer . . . . . . . . . . . . . . . . . . . . . . . . 5Figure 1.3 Destruction of acinar parenchyma in emphysema. (Re-produced with permission from [2], Copyright Mas-sachusetts Medical Society) . . . . . . . . . . . . . . . 7Figure 2.1 The pulmonary airways. Adapted with permissionfrom [1]. . . . . . . . . . . . . . . . . . . . . . . . . . 11Figure 2.2 The pulmonary acinus. Adapted with permission from[3]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12Figure 2.3 Mechanical analogy of airways in human lung. Adaptedwith permission from [4]. . . . . . . . . . . . . . . . . 13Figure 2.4 Effect of COPD on FEV1 and FV C (Image repro-duced with permission from HSE Digital Communi-cations) . . . . . . . . . . . . . . . . . . . . . . . . . . 15Figure 2.5 Variation of Flow Rate vs Volume for healthy andCOPD afflicted lungs (adapted with permission from[5]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16xiLIST OF FIGURESFigure 2.6 Detailed cutaway of pulmonary acinus (Shared underCreative Commons Attribution 2.5 Generic license.) . 18Figure 2.7 Alveolar shapes proposed by Weibel (Adapted withpermission from [6]) . . . . . . . . . . . . . . . . . . 19Figure 2.8 SEM image of terminal alveolar regions showing aclose-packed honeycomb-like structure (Adapted withpermission from [3]) . . . . . . . . . . . . . . . . . . 19Figure 2.9 Close packed acinar geometry developed by Fung (Adaptedwith permission from [7]) . . . . . . . . . . . . . . . 21Figure 2.10 ’Seperatrix’ in alveolar flow (Adapted with permis-sion from [8]) . . . . . . . . . . . . . . . . . . . . . . 24Figure 2.11 Recirculatory flow in the pulmonary alveolus (Adaptedwith permission from [9]) . . . . . . . . . . . . . . . 25Figure 2.12 Relative streamlines showing the flow- field differ-ences in proximal, medial and distal generations ofthe acinus (Adapted with permission from [10]) . . . 26Figure 2.13 Deposition patterns of 3µm particle trajectories un-der gravity, coloured by velocity magnitude (m/s)(Adapted with permission from [11]) . . . . . . . . . 29Figure 2.14 Particle deposition efficiencies under differing gravi-tational orientations (Adapted with permission from[11]) . . . . . . . . . . . . . . . . . . . . . . . . . . . 30xiiLIST OF FIGURESFigure 4.1 Computational domain of section of pulmonary aci-nus for present study. The blue and yellow arrowsindicate the flow direction during inspiration and ex-piration, respectively. . . . . . . . . . . . . . . . . . . 37Figure 4.2 Progressive destruction of alveolar septa in the com-putational domain. (a) Healthy case (b) Case I (c)Case II (d) Case III. . . . . . . . . . . . . . . . . . . . 38Figure 4.3 (a) Displacement of the airway wall from rest duringinspiriation and expiration for the healthy and dis-eased cases. (b) Temporal variation of flow rate intothe pulmonary acinus during inspiration and expira-tion for the healthy and diseased cases. . . . . . . . . 42Figure 4.4 Computational mesh composed of tetrahedral elements.Note refined areas at inlet and outlet . . . . . . . . . 45Figure 4.5 Solution Algorithm . . . . . . . . . . . . . . . . . . . 51Figure 5.1 Temporal variation of the percentage increase in oxy-gen concentration within a duct (C) for cases withboth septal destruction and emphysematous wall mo-tion/inlet flow rates. . . . . . . . . . . . . . . . . . . . 58Figure 5.2 Temporal variation of the percentage increase in oxy-gen concentration within a duct (C) for cases withonly emphysematous septal destruction. . . . . . . . . 59xiiiLIST OF FIGURESFigure 5.3 Temporal variation of the percentage increase in oxy-gen concentration within a duct (C) for cases withonly emphysematous wall motion/inlet flow rates. . . 60Figure 5.4 Temporal variation in the static pressure drop acrossthe alveolar ducts for cases with both septal destruc-tion and emphysematous wall motion/inlet flow rates. 62Figure 5.5 Temporal variation in the static pressure drop acrossthe alveolar ducts for cases with only emphysematousseptal destruction. . . . . . . . . . . . . . . . . . . . . 63Figure 5.6 Temporal variation in the static pressure drop acrossthe alveolar ducts for cases with only emphysematouswall motion/inlet flow rates. . . . . . . . . . . . . . . 64Figure 5.7 Velocity contours at peak inspiration (a) Healthy case(b) Case I (c) Case II (d) Case III. Arrows indicategeneral flow direction in the acinus. . . . . . . . . . . 67Figure 5.8 Velocity contours at peak expiration (a) Healthy case(b) Case I (c) Case II (d) Case III. Arrows indicategeneral flow direction in the acinus. . . . . . . . . . . 68Figure 5.9 Streamlines in an individual alveolus from Duct 1near the end of inspiration for the healthy and dis-eased cases. Streamlines are coloured according tovelocity magnitude. . . . . . . . . . . . . . . . . . . . 69Figure 5.10 Orientation of the gravity vector with respect to thecomputational geometry for the normal and tangen-tial cases respectively . . . . . . . . . . . . . . . . . . 73xivLIST OF FIGURESFigure 5.11 Deposition fraction of 1 µm particles, with the gravityvector oriented in the normal direction . . . . . . . . 74Figure 5.12 Deposition fraction of 3 µm particles, with the gravityvector oriented in the normal direction . . . . . . . . 75Figure 5.13 Deposition fraction of 1 µm particles, with the gravityvector oriented in the tangential direction . . . . . . . 76Figure 5.14 Deposition fraction of 3 µm particles, with the gravityvector oriented in the tangential direction . . . . . . . 77Figure 5.15 Particle deposition efficiencies under differing gravi-tational orientations (Adapted with permission from[11]) . . . . . . . . . . . . . . . . . . . . . . . . . . . 78Figure 5.16 Plot of final deposition fraction after five breath cy-cles for 1 µm particles, gravity vector in normal di-rection. The filled circles indicate healthy acinar ductwhile empty circles indicate emphysematous destruc-tion. The labels indicate the duct number. . . . . . . 80Figure 5.17 Plot of final deposition fraction after five breath cy-cles for 3 µm particles, gravity vector in normal di-rection. The filled circles indicate healthy acinar ductwhile empty circles indicate emphysematous destruc-tion. labels indicate the duct number. . . . . . . . . . 82xvLIST OF FIGURESFigure 5.18 Plot of final deposition fraction after five breath cy-cles for 1 µm particles, gravity vector in tangentialdirection. The filled circles indicate healthy acinarduct while empty circles indicate emphysematous de-struction. labels indicate the duct number. . . . . . . 83Figure 5.19 Plot of final deposition fraction after five breath cy-cles for 3 µm particles, gravity vector in tangentialdirection. The filled circles indicate healthy acinarduct while empty circles indicate emphysematous de-struction. labels indicate the duct number. . . . . . . 84xviAcknowledgementsFirst and foremost I would like to gratefully acknowledge the guidanceand help provided by my principal supervisor, Dr. Joshua Brinkerhoff. Noless important was the direction provided by Dr. Andre Phillion, who de-fined the extent of the thesis.I would also like to acknowledge the help provided to me by Dr.JamesC. Hogg, Dr. Dragos Vasilescu and Dr. Tillie Hackett of the the Heart andLung Institute in Vancouver, for providing necessary and valuable insightinto the nature of COPD.My colleagues and friends at the Okanagan CFD Lab have provided in-valuable assistance in the form of ideas and feedback and this thesis couldhave been borne to fruition without their help. I would also like to ac-knowledge other colleagues and friends at the School of Engineering at UBCOkanagan, particularly Mr. Oleg Shabarchin.xviiDedicationTo my parents, for the start they gave me in life.xviiiChapter 1Introduction1.1 MotivationThis thesis studies the effects of early onset emphysema in human lungsby means of numerical simulation of the airflow in diseased and healthyportions of the lungs. Emphysema is one of the component conditions ofChronic obstructive pulmonary disease (COPD), which is defined as “... apersistent airflow limitation that is usually progressive and associated withenhanced chronic inflammatory response within the lung to noxious particlesor gases” [12]. Thus COPD is an umbrella term for a host of conditions whichmay produce the symptoms that are outlined in its definition.COPD is projected to become the third leading cause of death worldwideby 2020 [12]. In addition to direct economic costs associated with the treat-ment of the disease (estimated to be $1.5 billion per year in Canada alone[13]), COPD also has significant impacts on quality of life and productivityin the workplace [12]. Clearly, a deeper understanding of the componentsof this disease, its genesis, progression and its attendant effects is of theessence.11.2. Background1.2 Background1.2.1 Human lung anatomy and physiologyThe human lungs are a pair of air-filled organs that are primarily respon-sible for gas-exchange between the atmosphere and the bloodstream. Air istransported into each lung through a network of dichotomously branchingairways which branch off from the trachea which is the primary airway andthese airways decrease in diameter with each distal generation. The airwaysare divided into two distinct generations, the conducting and respiratoryairways. The former transports the air into the respiratory zone which isresponsible for the gas-exchange process. While the conducting airways aresmooth walled, the respiratory zone has airways whose walls are populatedwith structures known as known as alveoli. These walls of these alveoliare thin enough to permit diffusive exchange of gases between the lung airand the blood-stream thereby supplying fresh oxygen while carrying awaycarbon-dioxide.1.2.2 Chronic obstructive pulmonary disease (COPD)The chronic airflow limitation in COPD is caused by two specific syn-dromes: emphysema and small airways disease (obstructive bronchiolitis)[12]. Emphysema is the deterioration of pulmonary tissue, called parenchyma,and it results in loss of respiratory airways for gas exchange and elastic re-coil. Small Airways Disease affects airways of diameters less than 2 mm andcauses structural changes including remodelling and narrowing. Chronicbronchitis is also frequently included as one of the component entities of21.2. BackgroundFigure 1.1: The pulmonary airways. Adapted with permission from [1].COPD [12, 14]. It is defined as “...the presence of cough and sputum pro-duction for at least three months of two successive years” [12]. However,chronic bronchitis may exist independently of the airflow limitation, consid-ered to be the defining feature of COPD [15], and is thus not considered inthe present study. The extent of contribution of emphysema, small airwaysdisease and chronic bronchitis and indeed their presence is highly variablewithin the population of COPD patients. Thus, in some cases, emphy-sema accounts for a significant part of the airflow limitation, while in otherinstances it might only play a secondary role. Additionally, COPD oftenco-exists with other diseases, which may complicate its diagnosis [12].The most well-known risk factor for COPD is inhaled cigarette smokewith the amount of and duration of smoking contributing directly to the rela-31.2. Backgroundtive severity of the disease [12]. However, there is evidence that non-smokersmay also develop chronic airflow limitation, and among people with the samesmoking history, not all will develop COPD [12]. Globally, COPD has alsobeen linked to sustained exposure to combustion products of biomass fuelssuch as coal, wood and straw [16], and inhalation of organic and inorganicdusts, chemical agents and fumes [17]. Additional factors, including genet-ics, gender, age, lung growth, socio-economic status and childhood infectionshave also been acknowledged as having a role in the progression of this dis-ease [12, 16]. The most well established genetic risk factor is the deficiencyof alpha—1 antitrypsin. This serum acts as an inhibitor and is releasedduring lung inflammation in order to protect the body from harmful effectsof other enzymes [18].1.2.3 Diagnosis of COPDThe three primary symptoms of COPD are chronic cough, breathless-ness and sputum production. Tightness of the chest and wheezing soundsmay also be present [12]. A diagnosis of COPD is usually considered for anypatient that displays these three characteristic symptoms. However, clin-ical diagnosis is only made after performing a Lung Function Test (LFT)[12, 19]. This operation is carried out in two identical parts, before andafter administration of an inhaled bronchodilator. Specifically, the patientis asked to take in a full breath and then exhale the air in a forced manneruntil the lungs are emptied and with the aid of a recording device known as aspirometer (Fig. 1.2), two parameters are measured. The Forced ExpiratoryVolume in one second (FEV1) is the volume of air exhaled in the first sec-41.2. Backgroundond, while the Forced Vital Capacity (FV C) is the total volume of air thatcan be exhaled with no time constant. In both parts, the ((FEV1)/FV C)ratio is calculated and if this value is less than 0.7 after the inhalation ofthe bronchodilator, then a diagnosis of COPD is usually considered [12].Figure 1.2: SpirometerStudies have shown that early detection and treatment of COPD yieldssubstantially better results for patients [12, 20]. However, the early detec-tion of COPD is problematic, since in a majority of the cases patients areeither asymptomatic or attribute the early symptoms of COPD, such asincreasing lack of breath, to advancing age [20]. Airway conditions in thelower airways are also difficult to diagnose [21] with the standard lung func-tion tests and spirometry. This results in many patients being diagnosedonly after considerable and irreversible damage has occurred to the lungs.Additionally, commonly published spirometry data has significant discor-51.2. Backgrounddance amongst various references [22] and the American Thoracic Societyhas recommended against using spirometric parameters to diagnose lowerairway disease [23]. Therefore, there exists a gap in understanding withinthe medical community of the effects of the progression of emphysematousdestruction in the earliest stages of COPD, which in turn hinders access topossible routes towards early diagnosis. While it has been confirmed thatthe earliest stages of emphysema are characterised by the destruction of air-ways in the pulmonary acinus, the effects of their destruction on the airflowand particle deposition in the alveolar spaces are unknown. In addition,while this destruction has been correlated to the increase in resistance toairflow, the effects of this decreased airflow and indeed the effects of theprogressive nature of the destruction have yet to be quantified.1.2.4 EmphysemaThe word emphysema originated from the Greek term “emphysan” ,meaning “blow into” [24]. When applied to the lungs, it implies excessiveair within the parenchymal tissue of the lung [24]. The definition of emphy-sema as it stands today originated in a 1984 workshop held by the Divisionof Lung Disease at the National Heart and Blood Institute, which statedthat emphysema is “... a condition of the lung, characterized by abnor-mal, permanent enlargement of air-spaces, distal to the terminal bronchiole,accompanied by destruction of their walls and without obvious fibrosis”[25, 26].The emphysematous condition of the lung is characterized by loss ofalveolar walls, as depicted in Fig. 1.3, which in the case of closely packed61.3. ApproachFigure 1.3: Destruction of acinar parenchyma in emphysema. (Reproducedwith permission from [2], Copyright Massachusetts Medical Society)lower generations of the acinus also forms the with subsequent destructionof the capillary bed. The air-spaces are dilated and the small airways arenarrowed and have atrophied walls [5]. This implies that pulmonary airwayslose their connection with the inner surface of the lung, which hitherto hadbeen provided by the elastic parenchymal mesh.1.3 ApproachTraditional experimental methods of direct observations on the effectsof this disease are not applicable for the study of early stage COPD. WithComputed Tomography (CT), airways are only visible for approximately 671.3. Approachor 7 generations, while the airways which form the pulmonary acinus areusually below the resolution of conventional CT [27]. Magnetic ResonanceImaging (MRI), which is also frequently used for pulmonary diagnostics ishampered by a low signal return due to lack of tissue in the alveolar spaces.This is further excarbated in emphysema due to loss of tissue [28]. As aresult, while excised lung tissue from deceased human subjects may later bestudied ex vivo under high resolution micro-CT studies, the same cannotbe said for living tissue under in vivo conditions. This shortcoming in theexperimental method is hereby overcome by the use of numerical simulation.A simulation-based approach is also justified due to the ability to artifi-cially control the conditions of the simulation by the modeller. This allowsfor identifying the relative importance of various effects that cannot be inde-pendently controlled in a laboratory setting. Not only does ComputationalFluid Dynamics (CFD) simulation allow observation of the effects of the dis-ease on the airflow, but it also allows the progressive nature of the diseaseto be modelled with a greater degree of control than can be obtained fromex vivo samples. Additionally CFD simulation techniques allow real-timevalues of variables such as particle deposition to be simulated under a vari-ety of controlled conditions that would be extremely hard to replicate underexperimental conditions.This thesis studies the effects of destruction of parenchymal tissue in theacinar airways due to emphysema by building a model which incorporatesthe effects of the disease on the airway wall motion and airflow in the pul-monary acinus. Specifically, numerical simulations of the airflow are thencarried out in order to determine the effects of the disease, the severity of81.4. Overview of the thesiswhich is progressively increased over successive iterations. The results ofthese simulations are then subsequently used to determine consequences ofthe progressive nature of the disease on oxygen transport, pressure drop inthe alveolar ducts and particle deposition in the pulmonary acinus.1.4 Overview of the thesisThe remainder of the thesis is divided into five parts. Chapter 2 includesa detailed literature review covering COPD in general, emphysema in par-ticular, and CFD approaches that have been employed to date in order tosimulate airflow in the human lungs. This is followed by a detailed explana-tion of the scope and objectives of the thesis in Chapter 3. Chapter 4 detailsthe numerical methods that have been employed in the study, followed bythe results and their discussions in Chapter 5. A summary of the thesis, itsconclusions, and recommended future work is given in Chapter 6.9Chapter 2Literature Review2.1 Anatomy of the human lung2.1.1 Airway structureThe root of the human airway tree is the trachea or windpipe. Thisstructure, which measures 1-2 cm in diameter in normal adult human lungs[29], branches in the chest cavity with one daughter branch per lung. Theleft and right bronchi subsequently undergo multiple dichotomous subdivi-sions leading to daughter bronchioles and finally terminal bronchioles. Thisdichotomous branching was detailed in Fig 1.1 and is repeated in Fig 2.1.The detailed anatomical structure of airways was first characterisedthrough the study of preserved cadaver lung tissue [6, 30]. These studiesestablished the dichotomous nature of the airway branching. These airwayswere later shown to span up to an average of 23 generations [31]. The air-way generations are numbered from the trachea which is generation zeroand each subsequent branching increases the generation number by one asseen in Fig 2.1. Also depicted in Fig. 2.1 is the functional classification ofthe pulmonary airway system. The airways may be broadly divided intotwo regions. The first region -the conducting airways- has the trachea as its102.1. Anatomy of the human lungFigure 2.1: The pulmonary airways. Adapted with permission from [1].root and consists of smooth-walled airways whose sole function is to con-duct the flow deeper into the lung [6]. This region spans approximately 14generations of dichotomous branching. The final bronchioles of this gener-ation are termed as the terminal bronchioles. The subsequent generations-the respiratory airways- comprise the respiratory zone. From herein, asthe flow descends, the walls of the airways begin to be populated by indi-vidual, isolated structures known as alveoli. These transitional bronchioleswith isolated alveoli form the beginning of the respiratory unit known as thepulmonary acinus. The acinus represents the basic unit responsible for ex-change of gases between the bloodstream and the lung. The alveoli continueto increase in size among the more distal airways [32] until finally they coverthe entire surface of the airway with each alveolus having an opening into112.1. Anatomy of the human lungthe common duct [33]. Fig. 2.2 depicts a single pulmonary acinus spanningeight generations; the respiratory bronchioles showonly isolated alveoli whilealveolar ducts are completely covered with alveolar outgrowths.Figure 2.2: The pulmonary acinus. Adapted with permission from [3].2.1.2 Mechanism of pulmonary ventilationThe physiological function of pulmonary ventilation — or what is morecommonly known as the act of breathing — is to ensure a gas exchange thatmeets the requirements of tissue metabolism in the body, thereby ensuringhomeostasis [34]. This requires airflow through the pulmonary airways, tothe alveoli, which are the sites of gas exchange. The presence of such aflow may be achieved by creating a pressure difference between the two endsof the airways. Therefore, during inspiration, the pressure in the airwaysmust be lower than the ambient pressure at the mouth and vice versa duringexpiration.Inspiration is an active process, whereby the contraction of the diaphragmand the external intercostal muscles expand the tissue matrix, inflating the122.1. Anatomy of the human lungFigure 2.3: Mechanical analogy of airways in human lung. Adapted withpermission from [4].airways and generating the suction pressure required [34]. This is illustratedin Fig. 2.3 where the springs attached to the tube representing the airwaymimic the action of the elastic tissue matrix during the breath cycle. As theair flows into the lungs, the pressure equalises with the atmosphere, untilthere is no more flow, which signals the end of inspiration. Normal or tidalexhalation by contrast is a passive process. As the intercostal muscles relax,the natural elasticity of the parenchymal tissue recoils the fibres attached tothe airways thereby compressing them. This phenomenon, known as elasticrecoil, is responsible for compressing the airways, increasing the pressure inthem, and driving the air out.132.1. Anatomy of the human lung2.1.3 Effects of emphysema on mechanical behavior of lungairwaysThe primary effect of emphysema on the mechanics of the lung is theloss of elastic recoil. This in turn reduces the expiratory driving pressure[35]. This translates to a reduced airflow during expiration in a conditionthat is termed as Expiratory Flow Limitation (EFL). Under EFL conditionsthe expiratory flow rate becomes independent of expiratory muscle effort.Instead the flow rate is determined by the static lung recoil pressure andresistance of the airways up-stream from the affected segments [36, 37]. TheEFL conditions also increase the time required for a complete expiration ofthe air inhaled during the inspiratory phase. This fact coupled with the lossof elastic recoil directly leads to an increase in the End-Expiratory-LungVolume (EELV)- the residual volume of the lung at the end of expiration[35]. This is in turn quantified by the variable forced vital capacity (FV C).FV C is defined as the amount of air which can be forcibly exhaled from thelungs after taking the deepest breath possible. This can be measured usingspirometric techniques discussed in Chapter 1 and as depicted in Fig 2.4COPD also results in a decrease in FV C. This is directly linked to the lossof elastic recoil [5] and results in air-trapping in the lung. This adds to thepermanent enlargement of the air-spaces, parts of which had been previouslyrendered void, by the loss of parenchymal tissue.As a result the expiratory flow is limited and significantly hampered. Ad-ditionally, during inspiration the maximum flow rate is observed to decreasein patients with emphysema. As a result, both during inspiration and expi-142.1. Anatomy of the human lungFigure 2.4: Effect of COPD on FEV1 and FV C (Image reproduced withpermission from HSE Digital Communications)ration the natural self-similar motion of the alveolar airways is truncated andas a direct consequence the magnitude of airflow through them, decreases.The classic description of how the flow rate of air in and out of the lungsvaries with volume is given in Fig. 2.5. During inspiration, the flow ratedependence on volume is appreciably sinusoidal (COA). During expiration,the flow rate rises with decreasing volume linearly (OB), until it reaches amaximum expiratory flow rate (B). Thereafter it falls to zero linearly (BC),signalling the end of complete expiration. The values of the correspondingquantities for a patient suffering from COPD are shown in red on the same152.1. Anatomy of the human lungFigure 2.5: Variation of Flow Rate vs Volume for healthy and COPD af-flicted lungs (adapted with permission from [5])162.1. Anatomy of the human lungfigure. As mentioned before the maximum flow rate during inspiration andexpiration for the diseased case (D and E respectively) is reduced, whencompared to the healthy case.Although these values in Fig. 2.5 reflect the variation of the wholelung volume and the flow rate measured at the mouth, it is assumed theywould reflect the volume flow rate relationship in an alveolar duct. It is alsoimportant to note that the volume values plotted on the horizontal axis donot represent the actual volume of the lung or alveolus itself, but rather,the variation in volume from the resting position prior to inspiration as isrepresented by points O and C in Fig The pulmonary acinusThe dichotomously branching airways of the pulmonary airway tree ter-minate in structures that are termed as acini. The pulmonary acinus formsthe respiratory zone of the airway tree and is termed as such since it is theprincipal location of gas exchange between the inhaled air and blood vessels.This complexity of structures necessary to achieve efficient gas-exchange isdetailed in Fig 2.6. The velocity of airflow in these airways is much lowerthan that of the higher conducting airways and this serves to promote thediffusive exchange of gases.The geometrical shape of the alveoli is a matter of some debate. Theseminal work by Weibel [6] mentions three distinct possibilities based onsurface-area and volumetric measurements. Two are trapezoids with thethird being a spherical cap as shown in Fig. 2.7. In recent years the modelof the spherical cap with parameters established by [32] has been popular172.1. Anatomy of the human lungFigure 2.6: Detailed cutaway of pulmonary acinus (Shared under CreativeCommons Attribution 2.5 Generic license.)182.1. Anatomy of the human lungFigure 2.7: Alveolar shapes proposed by Weibel (Adapted with permissionfrom [6])in studies relating to flow in the respiratory zone [9, 38, 39]. However,this geometry has a number of shortcomings. In particular, it does notrepresent the close packed nature of the alveoli that is usually observed inthe terminal regions of the respiratory zone (Fig. 2.8). In order to addressFigure 2.8: SEM image of terminal alveolar regions showing a close-packedhoneycomb-like structure (Adapted with permission from [3])the shortcomings of the alveolar cap models, an alternative arrangement of192.2. Numerical simulations of alveolar airflowsalveoli in the acinar regions was proposed by Fung [7]. The proposed modelwas envisaged on the basis of three assumptions; 1) all alveoli are equaland space-filling before they are ventilated, 2) they are ventilated to ductsas uniformly as possible, 3) the alveoli are reinforced at their mouths forstructural integrity. In this model the basic unit of the acinus is a truncated14-hedron, henceforth termed as an order 1 polyhedron (Fig. 2.9a). A single14-hedron, surrounded by identical polyhedra on each of its faces forms asingle solid unit in space (Fig. 2.9b), henceforth termed as an order 2, 14-hedron. This forms a basic unit of the acinar-structure. A combinationof order-2 and order-1 polyhedra can then be assembled into a ductal treeforming a space-filling structure (Fig. 2.9b.). The structures formed byclose-packing alveoli of this shape resemble honeycombs and have featuredin recent works which have tried to accurately depict the space-filling natureof the alveoli [8, 11]. The fact that the geometrical structure of the alveoliis significant enough to influence properties such as deposition of particles[40, 41] requires that the choice of any geometry for a particular simulationbe accompanied by proper justification.2.2 Numerical simulations of alveolar airflows2.2.1 OverviewThe airways of the lung present a complex environment in where lowvelocity flow transports gases through both convective and diffusive means.While the geometries may be approximated to some degree as cylindricalstructures, the dichotomous branching of the upper lung and the presence of202.2. Numerical simulations of alveolar airflowsFigure 2.9: Close packed acinar geometry developed by Fung (Adapted withpermission from [7])alveolar openings in the respiratory zone make for a complex environment.Consequently, computational fluid dynamics (CFD) simulations have beenused to gain additional insights into the flow structures during inhalationand exhalation. This section reviews the mathematical background of CFDsimulations and focuses specifically on its use for simulating alveolar flow.2.2.2 Governing equationsRespiratory air-flows are generally considered to be examples of Newto-nian and incompressible flows[1, 42]. They are therefore, described by theNavier-Stokes equations in a dimensionless form as given in [43].(Wo)2∂u∗∂t∗+Re(u∗ · O∗)u∗ = −O∗p∗ + O∗2u∗ (2.1)where:O∗ = DO, t∗ = ωt, u∗ = uU, p∗ =pµU/D212.2. Numerical simulations of alveolar airflowswith u being the velocity field, p is the pressure, t is the time, ρ isthe density, D is the airway diameter, µ is the molecular viscosity, andω = (2pi/T ) is the angular frequency of the breathing time period T .Two particular dimensionless numbers are important in the context ofthe given equations. The Reynolds number given byRe =ρUDµrepresents the relative influence of the inertial forces and viscous forces onthe fluid flow. In pulmonary airflows during quiet sedentary breathing, Reis O(1) or smaller [43], therefore the inertial forces are generally negligibleand the flow in the acini is laminar. Due to the cyclic nature of the flowin the pulmonary airways, a second dimensionless number is required tocharacterize this unsteadiness. This is the Womersley number,Wo = D(ωUµ)0.5which compares unsteady acceleration to viscous effects. For typical acinarflows, Wo is O(0.1) [43]. This implies the frequency of airflow cycles issufficiently low that a parabolic velocity profile has time to develop duringeach cycle. The dichotomous branching nature of the airway tree lends itselfto a straightforward computation of the airflow rate in the zth generation ofthe tree asU(z) = 4Q/(2zpiD(z)) (2.2)where Q is the flow rate at the trachea (z = 0) and D(z) is the diameter of222.2. Numerical simulations of alveolar airflowsthe airways in the zth generation. However, it is important to note that thisrelation is only valid for ideal cylindrical airways and identical diameters atany given generation. In the conducting region, the diameters of the airwaysconform to theoretical predictions [44], and thus the velocities conform toEqn. 2.2. In the acinus, however, the effective diameters are larger [32], andvelocities tend to decay more steeply. It should also be noted that Eqn. 2.2 isonly valid if the longitudinal paths to all airways in a particular generationare assumed to be identical. Even though the transpulmonary pressureacross the airways may be assumed to be homogeneous throughout the lung,the airways themselves may be of different longitudinal lengths. As a result,for a constant pressure drop ∆P, an airway with a longer longitudinal pathto the terminal alveolar sacs will have higher viscous resistance, and thuslead to an airway with a lower flow rate despite sharing the same diameteras its neighbor of the same generation with a shorter longitudinal path [43].2.2.3 The nature of alveolar flowTo the best of the author’s knowledge, the earliest numerical analysisof flow in the acinus was an investigation of creeping flows by the methodof solving the bi-harmonic equation for the stream-function [45]. Studiesof flow in packed acini followed thereafter [46, 47]. Simulation of flows inthe alveolar domain has been primarily motivated by the need to studyhow aerosols and other fine particles behave in this regime [40, 48, 49]. Tothis end, investigations have shown that the orientation of the branchingalveolar ducts with respect to each other as well as to the gravity vector areimportant in determining deposition sites of micro-particles [49, 50].232.2. Numerical simulations of alveolar airflowsThe alveolar topology makes alveolar flows resemble a ductal shear layerpassing over a cavity formed by the opening of the alveolar mouth. Thisis illustrated in Fig. 2.10. Such an arrangement can give rise to flow sepa-ration and formation of a recirculation zone, depending on geometrical andflow parameters. In the flow regime encountered in the acinar region, therelative arrangement of the duct and cavity determines the occurrence andnature of flow-separation, particularly the location of the separation surfaceor ‘seperatrix’ [8, 51].Incorporating wall-motion into the simulation mimics the distentionsof the lung parenchyma, and is crucial in capturing the three-dimensionalrecirculating flow structures that are characteristic of acinar flows [39].Figure 2.10: ’Seperatrix’ in alveolar flow (Adapted with permission from [8])Since alveolar flows are a result of pressure gradients arising from dis-tension of the lung airways, the effect of airway motion on the flow is alsocritical. While there is evidence that rigid walled alveolar models give riseto recirculatory flows, [49, 52], wall motion kinematics into numerical sim-ulations of acinar airflows result in flow patterns that are irreversible andchaotic in nature [9]. This phenomena is termed as chaotic advection andis frequently observed in 2D time dependent and 3D cavity flows [53, 54].242.2. Numerical simulations of alveolar airflowsThis makes the alveolar flows inherently complex. Recirculating flows cre-Figure 2.11: Recirculatory flow in the pulmonary alveolus (Adapted withpermission from [9])ate stagnation points that result in irreversible paths of massless particles[9, 38]. The first simulations to incorporate moving alveolar walls indicatedthat in realistic alveolar geometries have the stagnation point form at theproximal end of the alveolar cavity, as shown in in Fig. 2.11. Subsequentsimulations incorporating 14-hedron geometries showed that the positionand strength of the recirculatory centre is dependent on the location of thealveolar generation. The proximal generations have a strong recirculationcentre develops at nearly the centre of the cavity. Distal generation showa weakening of the centre accompanied by a shift to the proximal end ofthe alveolus while the distal generation show no recirculation as shown inFig. 2.12. The velocities within the alveolar cavities are on average one ortwo orders less in magnitude than in the duct [8, 39].252.2. Numerical simulations of alveolar airflowsFigure 2.12: Relative streamlines showing the flow- field differences in prox-imal, medial and distal generations of the acinus (Adapted with permissionfrom [10])Under healthy conditions the nature of flow within the alveolus is heavilydependent on its location within the acinar generation. In the proximalgenerations a strong recirculatory zone is observed with the vortex centreclose to the centre of the alveolus. In the medial generations the centremoves to the proximal end of the alveoli while in the distal generations therecirculatory centre disappears and is replaced by radial flows.[10]2.2.4 Numerical simulations of particulate transport inalveoliGoverning equationsThe transport of particles in the pulmonary acinus is described by thefollowing equation [43]mdupdt= FD + Fg + F (t) (2.3)262.2. Numerical simulations of alveolar airflowsWhere the particle acceleration on the left hand side is a result of the threeforces on the right. Fg is the gravitational force on the particle with Fg =mg. FD is the drag force on the particle, which in the instance of a low Reflow, is given by the Stokes Law:FD =−3pidpµ(up − uf )Cc(2.4)The drag force is therefore dependent on the particle diameter dp, the fluidviscosity µ and the relative velocity of the particle up with respect to the en-training fluid uf . Cc is the Cunningham slip factor [55]. The third term onthe right F (t), is a stochastic force representing the Brownian motion due tothe collisions of the fluid molecules with the entrained particles. These threeforces are representative of the three competing mechanisms that determinethe fate of an inhaled particle in the acinus. Namely, convection, sedimen-tation and diffusion. The extent to which each of these forces dominate theparticle dynamics in the pulmonary airways is a function of particle size andwill be discussed in the following section.Micron sized particle dynamics in the pulmonary acinusFor micron sized particles, where dp > 0.5µm, the diffusive force dueto the Brownian motion ceases to have a significant impact on the particledynamics [56]. Eqn. 2.3 may then be non-dimensionalized by some charac-teristic velocity u and length L, asStkdu∗pdt∗= −(u∗p − u∗) +ugg∗u(2.5)272.2. Numerical simulations of alveolar airflowswhereu∗p =upu; u∗f =ufu; t∗ =tuL; g∗ =g|g|Stk =ρpd2pu18µLis the Stokes’ number representing the ratio of the inertial forces to theviscous drag force,ugu=ρpd2p|g|18µ|u| ,and compares the relative effects of gravity to the viscous drag force with theterminal settling velocity given by ug. While the inertial term is significantlymore important in the conducting airways, in the acinar airways, due to thelow values of Re, it is the gravity term that dominates.It is important to note that the relative strength of the two mechanismsis highly dependent on particle sizes. Sznitman et. al. [11] investigated theparticle dynamics of 1 µm and 3 µm particles in two separate geometries;namely, a simple isolated alveoli on an airway, consisting of one half-sphericalalveolus and a space filling geometry of the entire acinus derived from themodel developed by Fung [7].It was observed, that while the larger 3 µm particles (as shown inFig. 2.13) remained relatively insensitive to the flow and were governedalmost entirely by the orientation of the gravitational field, the smaller par-ticles (1 µm) exhibited complex kinematics due to the coupling of the un-282.2. Numerical simulations of alveolar airflowssteady flow field, local convective effects due to the complex geometry ofthe acinus and the sedimentation effects of the gravitational field. For bothparticle sizes, the deposition patterns were largely non-uniform, determinedprimary by orientation of the gravity vector relative to the duct. Particle en-try into the alveoli itself was governed largely by the position of the alveolaropenings relative to the direction of sedimentation, which was the principaldeterminant of the particle kinematics. This also implies that in spite ofthe low flow velocities, under appropriate circumstances, particles may beswept deep into the acinus.Figure 2.13: Deposition patterns of 3µm particle trajectories under gravity,coloured by velocity magnitude (m/s) (Adapted with permission from [11])The deposition patterns of the two particle sizes was studied using thedeposition efficiency, defined as the ratio of the particles inhaled to thosedeposited over a specified time period. It was observed that while the heav-ier particles (3 µm) were deposited within a single breath cycle, the finerparticles (1 µm) would exhibit much longer residence times (3 times longer292.2. Numerical simulations of alveolar airflowsthan the heavier cases). As shown in Fig. 2.14, the heavier particles reach adeposition efficiency of approximately 1, within halfway of the first breathcycle (t/T = 0.5). In contrast the lighter particles in one case show a de-position efficiency of approximately 0.9 at (t/T = 1.5), while in the othercase the efficiency remains slightly above 0.8 even after t/T = 3. This factimplies that the finer particles remain longer in suspension within the flow-stream and concurs directly with their more complex coupled dynamics asdiscussed earlier.Figure 2.14: Particle deposition efficiencies under differing gravitational ori-entations (Adapted with permission from [11])302.3. CFD simulation of diseased acini2.3 CFD simulation of diseased aciniThe majority of the numerical studies detailed to date have dealt withundamaged alveolar geometries with flow and boundary conditions reflectinghealthy physiologies. The effect of emphysematous destruction associatedwith COPD on the airflow patterns has not received as much attention.Oakes et. al. [57] and Berg [58] have conducted experiments of airflowin emphysematous alveolar models, identifying the influence of emphyse-matous destruction on particulate deposition and flow patterns formed inrat airways. However, the effect of emphysematous destruction on the closepacked geometry on the acinus was not considered in the geometrical model.Recent work by Aghasafari et. al. [59] simulated the effects of emphysemaon the partial closing of the terminal alveolar sacs through CFD, but thesimulation neglected the influence of wall motion during the breath cycle.However, this model differed markedly from the close packed geometry thatis being considered in this study. As such it only studied airflow in a clusterof alveoli at the end of a duct, without any branching thereby overlook-ing the variation in airflow caused in multi-generational geometries due toemphysema. Also, the primary focus of this investigation was the collapseand closure of the alveolar sacs due to emphysema and not the spetal de-struction of the walls themselves. The most recent work to include a fullnumerical simulation of emphysema with multi-lobular geometry involves astudy of aerosol deposition in healthy and emphysematous rat lungs [60].In this study emphysema is simulated by the removal of inter-alveolar septaas well as changes in the lung motion and airflow rate, drawn from experi-312.4. Summarymental data. However, the conclusions drawn were based on the acinus asa whole and not at the level of individual ducts. Furthermore, the effectsof emphysema on oxygen transport were not studied in detail and whilesome conclusions were drawn regarding particle depositions in emphysema-tous acini, they effects on particle-deposition patterns with progressivelyincreasing levels of emphysema were not quantified.CFD modelling is therefore potentially a powerful tool for studying com-plex phenomena such as emphysema, where multiple factors affect the air-flow. CFD analysis not only allows the modeller to investigate the effects ofindividual symptoms of the disease, but also allows the study of the effectsof the disease as it progresses in severity. This allows the investigation andprediction of airflow characteristics in a controlled manner that is difficultto replicate in experimental studies.2.4 SummaryThe human airway tree is a complex multi-scale structure which is pri-marily responsible for the transport of air from the atmosphere to the blood-stream. The smallest units of these airways are the pulmonary acini whichare composed of alveolar structures arranged in a honeycomb like fashionaround the terminal airways. The airways are surrounded by a complexnetwork of elastic tissues which are responsible for providing mechanicalsupport and the means for lung inflation.Patients with emphysema have enlargement of their airspaces due to air-trapping. The progress of emphysema in the lung is heralded by the loss322.4. Summaryof small airways in the distal ends of the pulmonary airway tree. The air-ways also lose their natural elastic recoil which contributes to the difficultyof the patient breathing out, resulting in the afore-mentioned air-trapping.Emphysema also causes significant changes in breathing patterns and quan-tifiable variables based on spirometry are currently used to diagnose COPD.However such techniques are not sensitive enough to detect emphysema inits earliest stages.Computational Fluid Dynamics (CFD) is a powerful tool that has beenwidely utilized for studying human pulmonary airflow, diffusive transportof gases in the lungs and particle and aerosol transport and deposition inthe airways. While there has been some initial research in employing CFDfor the study of emphysematous alveoli, studies involving a complete modelutilizing the mathematically derived 14-hedron geometry of the acinus andincorporating the effects of both parenchymal destruction and airway wallmotion is lacking in the literature.33Chapter 3Scope and ObjectivesAs explained in Chapter 2, the study of emphysema, especially in itsearliest stages is hampered by the lack of suitable methods of diagnosisand the impossibility of performing in vivo study at the level of terminalbronchioles. CFD analysis, which has already been used for the study ofhealthy airways in the pulmonary acinus is uniquely suited to this task andtherefore, has been chosen as the principal method of investigation here.Therefore, the objectives of the study are to analyze the effects of em-physema on:(a) The flow-field in the ducts and alveoli of the acinus.(b) The transport of oxygen into the acinar airspaces during inspirationunder the competing influences of diffusion and convection.(c) The dynamics and deposition patterns of micron sized particles (dp =1 µm, 3 µm) in the acinus.To achieve these objectives, a computational analysis of airflow is per-formed on a section representing the fourth and fifth generations of the pul-monary acinus which was geometrically modelled using the close-packed 14hedron acinar geometry detailed in Section 2.1.2. The effects of emphysema34Chapter 3. Scope and Objectivesare incorporated into the overall computational model using a combina-tion of three factors, (i) destruction of intra-alveolar septa, consistent withdamage observed in emphysematous patients; (ii) reduction in airway wallmotion, due to loss of parenchymal tissue matrix; (iii) reduction in airflowrate during the breath cycle which is a result of the loss of elastic nature ofthe airways. These conditions are meant to simulate the conditions in theairways in the earliest stages of emphysema.35Chapter 4Methodology4.1 Computational domainTo accurately capture the flow through the airways in a CFD study, phys-iologically realistic geometries, boundary conditions and wall deformationsmust be used. The closely packed alveolar structures in the acinar regionare modeled using a geometry composed of fourteen sided polyhedrons (14-hedrons), following the work of Fung [7]. There are several examples ofthis geometry being used for CFD studies of alveolar airflow in literature[8, 10, 11, 41, 43]. The computational geometry for the present study isshown in Fig. 4.1. It consists of two generations of alveolar ducts that func-tionally correspond to the fourth and fifth generations of the pulmonaryacinus after the start of the respiratory zone which approximately corre-sponds to the twenty-second and twenty-third generations from the trachea.Each duct models the alveolar spaces shaped as a 14-hedron, with 33 of thesealveolar spaces in total. The collection of alveolar spaces that are ventilatedby a single airway are herein referred to as a duct. As seen in Fig. 4.1, thedomain consists of three ducts, one occurring proximal to the dichotomousbranching and associated with the fourth generation of the pulmonary ac-inus (Duct 1), and two occurring distal to the dichotomous branching and364.1. Computational domainFigure 4.1: Computational domain of section of pulmonary acinus forpresent study. The blue and yellow arrows indicate the flow direction duringinspiration and expiration, respectively.associated with the fifth generation of the pulmonary acinus (Ducts 2 and3).The destruction of the septal walls between the alveolar spaces due toemphysematous progression is modelled in the present study by progres-sively removing septa in the three ducts of the computational geometry. Asseen in Fig. 4.2, septal destruction progresses spherically from the centreof the duct towards the alveoli. The septa in each duct were degraded tothe maximum extent possible without affecting the exterior of the geome-try. This corresponds to the early stages of the disease in which large-scaledestruction of the alveolar spaces has not yet occurred. Four cases were sim-ulated to represent various stages in the progression of emphysema: Healthy,Stage I, Stage II, and Stage III. The healthy geometry has all of its septa374.2. Boundary conditions and wall motionFigure 4.2: Progressive destruction of alveolar septa in the computationaldomain. (a) Healthy case (b) Case I (c) Case II (d) Case III.intact, while Stages I, II and III have septal destruction in one, two, andthree ducts, respectively.4.2 Boundary conditions and wall motionBoundary conditions for the computational domain consist of no-slipwalls for all surfaces except the inflow and outflow boundaries. The bound-aries defined at the entrances and exits to the ducts, labelled with arrowsin Fig. 4.1, are defined as either inlet or outlet boundaries depending onthe time of the simulation within the breath cycle. During the inspiratoryphase, one inlet is defined at the entrance to Duct 1 and two outlets are de-fined for each of Ducts 2 and 3. During the expiratory phase, the situationreverses; two inlets are defined for each of Ducts 2 and 3, and one outlet is384.2. Boundary conditions and wall motiondefined for Duct 1.To ensure accuracy of the present study, the boundary information de-fined at the inlet and outlet boundaries must reflect the actual physiologyof the pulmonary acinus. Because in vivo measurements at this depth ofthe lung are impossible, the volume/flow rate relationships measured at themouth must be used to infer flow rates in the pulmonary acinus. A typicalvolume/flow rate relationship at the mouth was already presented in Fig. 2.5and can be approximated as follows. The flow rate during inspiration isQin = Qin, max sin(ωpV ) for 0 ≤ V < Vmax (4.1)and during expiration isQQexp = γV for Vmax, exp ≤ V < Vmax, (4.2)Qexp = −βV + κ for 0 ≤ V < Vmax, exp (4.3)where Qin, max is the maximum inspiratory flow rate, Vmax is the maxi-mum volume excursion of the lung/alveolus, ωp is the breathing frequency(reciprocal of the breath period), Vmax, exp is the volume corresponding tomaximum expiratory flow rate, and β, γ, and κ are constants. The effect ofCOPD on the volume/flow-rate relationship was also illustrated in Fig. 2.5.Although, the effects of COPD on the lung volume/air-flow rate relationshipis not exclusively due to emphysema (Small Airways Disease (SAD) may alsobe a component) - due to the lack of available in vivo data from patientssuffering exclusively from early stage emphysema, the boundary conditions394.2. Boundary conditions and wall motionTable 4.1: Variation of parameters governing flow-rate and airflow motiondue to emphysemaCase (Lm/2) β κ γ Reinsp ReexpHealthy 0.0265 -0.3816 0.0202 1.362 0.20 0.34Case I 0.0236 -0.39672 0.0187 1.3032 0.18 0.30Case II 0.0208 -0.41181 0.0171 1.2440 0.16 0.25Case III 0.0179 -0.4269 0.0154 1.186 0.14 0.23imposed on the simulation have been obtained from the interpolating thevalues between healthy and COPD observations in Fig. 2.5. This assump-tion might also justified due to the fact that SAD affects airways which arepresent higher in the airway tree than the acinus. Thus inflow conditionson the boundary of the diseased acinus should account for the combinedeffects of all the components of COPD. The effects of emphysema on thevolume/airflow rate relationship may be summarised as follows. Emphyse-matous degradation leads to the loss of elastic recoil of the lung parenchyma.As a result, airways do not inflate as effectively and the maximum inspira-tory flow rate and volume excursion at the end of inspiration is reduced.Due to the decreased elastic recoil, the expiratory flow rate also reduced.The reduction in flow rate can be modelled by altering the parameters β , κ,γ, Vmax and Qin,max in Eqns. 4.1-4.3 according to Table 4.1. The extent towhich the parameters are varied depends on the extent of emphysematousprogression. As is clearly shown in Table 4.1, the progression of the diseaseis accompanied by a reduction in mesh motion and fall in the inlet Reynoldsnumber. As mentioned before the values are obtained by interpolating be-tween the healthy and COPD observations from Fig. 2.5The airway wall-motion which drives the airflow in the airways during404.2. Boundary conditions and wall motionbreathing may be approximated as self-similar [39], meaning that the tem-poral variation of any length scale L(t) with respect to a reference lengthL0 is given byL(t) = L0[1 +Lm2+Lm2sin(ωp t− pi2)]. (4.4)By defining a parameter λ equal to the bracketed expression on the righthand side of Eqn. 4.4, the above expression can be simplified asL(t) = L0λ.The length scale Lm in Eqn. 4.4 is defined asLm = (C + 1)13 − 1 (4.5)where C is defined as the ratio of tidal lung volume to the functional residualcapacity (FRC),C =VtFRC. (4.6)The variation of the alveolar volume with time is therefore given byV (t) = V0λ3. (4.7)The variation in the alveolar volume is implemented in the simulationsthrough a user-defined function (UDF), reproduced in Appendix A, whichexpands and contracts the cells within the computational mesh in a mannergiven by Eqn. 4.7. Figure 4.3(a) plots the resulting relative displacement of414.2. Boundary conditions and wall motionFigure 4.3: (a) Displacement of the airway wall from rest during inspiriationand expiration for the healthy and diseased cases. (b) Temporal variationof flow rate into the pulmonary acinus during inspiration and expiration forthe healthy and diseased cases.the airway wall for the four simulated cases during inspiration and expira-tion. It should be noted that the values on the y-axis represent the relativechange from the resting position of the lung. The values are obtained bysubstituting the constants from Table 4.1 into Eqn. 4.7. The motion of thealveolar walls affects the flow rates into and out of the computational do-main derived in Eqns. 4.1-4.3. In order to capture this effect, Eqn. 4.7 needsto be substituted into Eqn. 4.1-4.3. The resulting inspiratory flow rate isQin = Qin, max sin(ωpλ3)for 0 ≤ t < tinsp (4.8)and the expiratory flow rate isQexp = γλ3 for tinsp ≤ t < texp, max (4.9)Qexp = −βλ3 + κ for texp, max ≤ t < texp. (4.10)These equations were imposed at the inlet boundaries of the computational424.2. Boundary conditions and wall motiondomain as temporally-varying inflow velocity conditions. At the outletboundary, a zero gauge static pressure condition is defined. These flowrate relationships are shown in Fig. 4.3(b) for the four simulated cases. Thevalues of the velocity at the inlet have been reported as scalars and it shouldbe noted that during expiration the air flows out of the domain. It is alsonotable that the largest deviation from the healthy case occurs for Case 3,as it has the largest degradation of the septa. The breath-cycle time periodτ for all simulations was chosen to be 3 seconds in accordance with similarstudies available in literature [11, 39]. This corresponded to a breathing-cycle period of ωp = 2.094. This value of breathing-period corresponds tofairly quick breathing which corresponds to light exercise.In healthy human lungs, the alveolar walls are richly supplied by bloodcapillaries, and function as sites of gas exchange between the inspired air andthe bloodstream. A concentration gradient exists at the alveolar-capillaryinterface which allows the exchange of oxygen and carbon-dioxide by dif-fusion. Computationally, Sapoval et. al. demonstrated that this may bethought of as a gradient dependent on oxygen and carbon dioxide concen-tration close to the alveolar-capillary boundary [61]. This gradient whenimposed on the walls of computational domain,enhances the diffusive trans-port of oxygen towards the walls during inspiration and carbon dioxide awayfrom the walls during expiration, which mimics the physiological process ofgas-exchange in actual lungs. Such an approach was used by Hofemeieret. al. [62] in their latest work simulating airflow changes in acinar struc-tures from infancy to adulthood where a constant gradient was assumed forthe wall-flux of oxygen. However, such an approach is imperfect, as this434.3. Spatial meshgradient is actually highly time-dependent and changes rapidly during thebreath-cycle [5]. Furthermore, the flow throughout the entirety of the acinusis not diffusive but significant convective transport may occur in the moreproximal generations of the alveolus during inspiration. In fact, availableliterature [63] suggests that the transition point during sedentary breathingfor convective and diffusive transport is the fourth and fifth generation ofthe acinus. For heavier breathing, which is simulated in the present study,this transition point may move deeper into the acinus. In the present study,the computational domain was chosen to be the fourth and fifth generationsof the acinus. This fact coupled with a lower breathing time period meansthat during a majority of the inspiratory cycle, the transport of the oxygenthroughout the domain is primarily carried out by convection rather thandiffusion. In addition to this the probe points for recording changes in oxy-gen concentration are located near the alveolar entrances, well away fromthe walls. These two factors justify the decision to not model the diffusivetransport of oxygen and carbon dioxide through the walls. Additionally,the precise effects of emphysematous destruction on the alveolar-capillarydiffusion are not known and thus could not have been correctly modeled asa factor in simulating progressive emphysema.4.3 Spatial meshThe geometry was meshed using ANSYS R© Meshing (Release 15.0) com-mercial software and the resulting mesh is shown in Fig. 4.4. An unstruc-tured meshing algorithm using tetrahedral elements was adopted. Care was444.3. Spatial meshFigure 4.4: Computational mesh composed of tetrahedral elements. Noterefined areas at inlet and outletTable 4.2: Mesh independence study, with the observed variable being thestatic pressure drop across Duct IMesh Mesh elements Static pressure (Pa) % change1 500,000 0.00325 —2 600,000 0.00332 2.153 700,000 0.00335 0.09taken to refine the mesh at the inlet(s) and outlet(s) as shown in Fig. 4.4. Amesh-independence study was carried out with sucessively refined meshes asshown in Table 4.2, with the static pressure in Duct 1 being the monitoredparameter to establish grid independence The finest mesh was composed of700,000 elements. It was observed that the monitored quantity changed byonly 0.09 % thereby showing that the mesh with 700,000 elements achievessuitable accuracy while ensuring reasonable solution times.454.4. Solution approach4.4 Solution approach4.4.1 Continuity and momentum equationsThe determination of a flow field requires the solution of the continuity(Eqn. 4.11) and momentum transport equations (Eqn. 4.12) along with thesolutions of the transport equations of any other quantity of interest. Theseequations are∂ρ∂t+ O.(ρu) = 0 (4.11)∂ρu∂t+ O(ρuu) = −Op+ O.τ + ρg (4.12)where τ represents the stress tensor. In the present study the commercialsolver ANSYS R© Fluent (Release 15.0) was used to solve the above equationsusing the finite volume method, which is described in detail in the subsequentsections.4.4.2 Finite volume methodIntegral formulationIn the finite volume method the integral forms of the continuity andmomentum transport equations are applied over each control volume orelement that makes up the mesh of the computational domain.∮V∂(ρu)∂t.dV +∮A(ρu).dA = 0 (4.13)464.4. Solution approach∮V∂(ρu)∂t.dV +∮A(ρuu).dA = −∮Ap.I.dA+∮Aτ .dA+∮Vρg.dV, (4.14)where the subscripts V and A indicate integrals over the volume and surfaceof the domain in question and I being the identity matrix.Spatial discretizationThe integral forms of the governing equations may then be discretized inorder to transform them into a system of linear equations. For the momen-tum and continuity equations this leads to a total of 4 equations, one foreach component of the velocity (Eqns. 4.15a, 4.15b, 4.15c) and the continuityequation (Eqn. 4.16):aPu =∑nbanbunb +∑pfA.ˆi+ Su (4.15a)aP v =∑nbanbvnb +∑pfA.jˆ + Sv (4.15b)aPw =∑nbanbwnb +∑pfA.kˆ + Sw (4.15c)∑NfacesJfAf = 0 (4.16)The coeffcients ap and anb represent known values at the cell centre in ques-tion and that of its immediate neigbours.The nature of the discretization dictates that the pressure and massflux at the cell faces be known in order to solve the equations. ANSYS474.4. Solution approachFluent uses a co-located scheme where both the pressure and cell velocitiesare both stored at the cell centres. However such a scheme is prone tounphysical “checkerboarding” of pressure and velocities. This is overcome bythe method outlined by Rhie and Chow [64]. Therefore, in order to calculatethe face flux Jf = ρunf where unf is the normal velocity through the face f ,instead of using linear interpolation a momentum weighted average is used:Jf = ρfap,c0un,c0 + ap,c1un,c1ap,c0 + ap,c1+df ((pc0 +Opc0.r0)− (pc1 +Opc1.r1)) (4.17)where un,c0 , un,c1 , pc0, pc1 are the pressures and normal velocities at the cellcentres on either side of the face and df is a function of the average of thetwo momentum coefficients ap,c0 and ap,c1 (from Eqn. 4.15) on either sideof the face. The pressure values at the cell faces are interpolated in thefollowing manner:pf =pc0ap.c0+ pc1ap.c11ap.c0+ 1ap.c1(4.18)Pressure velocity couplingWhile the discretization of the continuity and momentum equations pro-vide four equations (Eqns. 4.15- 4.16), the continuity equation does not haveany pressure terms and thus the pressure field needs to be coupled to thevelocity field in order to solve for the pressure field. In the present case, thesemi-implicit method for pressure-linked equations (SIMPLE) scheme orig-inally described by Patankar [65] is used for this. The SIMPLE algorithmuses a relationship between velocity and pressure corrections to enforce massconservation and to obtain the pressure field. The momentum equation is484.4. Solution approachfirst solved with a guessed pressure field p∗, which results in a face-flux J∗fthat does not satisfy the continuity equation. The corrected flux is thenobtained by adding a correction factor J′f defined asJf = J∗f + J′f (4.19)The SIMPLE algorithm then postulates that the correction factor J′f beexpressed asJ′f = df (p′c0 − p′c1) (4.20)where p′is the cell pressure correction. When Eqns. 4.20 and 4.19 aresubstituted into 4.16 we get a discretized equation for the pressure correctionterms as followsapp′=∑nbanbp′nb +∑NfacesJ∗fAf (4.21)All the terms in Eqn 4.21 except for p′are known. Once the lone unknownis determined the corrected pressure is determined asp = p∗ + αpp′(4.22)where αp is the under-relaxation factor. The corrected face flux which sat-isfies the continuity equation is then obtained asJf = J∗f + df (p′c0 − p′c1) (4.23)494.4. Solution approachTemporal discretizationTemporal discretization involves the integration of every term in thedifferential equations over a time step ∆t. For any generic variable φ and afirst order time evolution given by∂φ∂t= F (φ) (4.24)The equation can be discretized in time to second order (as is used inthe current model) as follows3φn+1 − 4φn + φn−12∆t= F (φ) (4.25)where the superscripts n−1, n, n+1 denote the previous,current and futuretime levels respectively. The simulation was integrated in time with a timestep size of ∆t = 10−4 seconds, which satisfies the CFL condition for alltimes within the breath period. At each time step, the algebraic system ofequations are solved iteratively until the root-mean-square residuals for con-tinuity, momentum conservation, and species mass fraction converge to lessthan 10−6, which is achieved in approximately 20 coefficient loop iterationsper timestep on average. The temporal algorithm for the entire solution isbest described by the means of the flowchart shown in Fig. Scalar transportA scalar transport equation was also solved in order to track the trans-port of respiratory gases into and out of the computational domain. The504.4. Solution approachFigure 4.5: Solution Algorithm514.4. Solution approachpassive transport of gas concentration is given by∂(φi)∂t+∇(~v.φi) = −1ρ∇ · ~Ji (4.26)where φi is the concentration for each species in the domain and Ji is thecorresponding flux of the species given by Fick’s law as shown in Equation4.27~Ji = −ρDi∇Yi (4.27)The diffusive coefficient of oxygen within air for this simulation (Di) has avalue of 0.213 cm2/s.In order to determine the values of the scalar φ at the cell faces, a secondorder upwinding (SOU) technique is used:φf,SOU = φ+ Oφ.r (4.28)where φ and Oφ are the values of the scalar and its gradient at the cellupstream of it with r is the displacement vector from the upstream cellcentroid to the face centroid.The incoming air was assumed to be composed solely of nitrogen, oxygen,carbon dioxide and water vapor with volumetric fractions of 0.78,0.21,0.0004and 0.0096 respectively. The exhaled air had reduced oxygen and increasedcarbon dioxide and water vapor with volumetric fractions of 0.15,0.04 and0.03 respectively with the nitrogen fraction remaining unchanged. The ini-tial conditions in the computational domain were assumed to be identicalto that of the expelled air.524.4. Solution approach4.4.4 Particle transportThe particle dynamics in the computational domain are governed by abalance of convective and gravitational forces given by Eqn. 2.5 in Section2.1. These equations were integrated over time in order to determine thetrajectories of inhaled particles. Eqn. 2.5 may be restated in a simplifiedform (for one velocity component) as follows:dupdt=(up − uf )τ+ a (4.29)where τ represents the drag force coefficients and a represents all otherexternal forces , which in this thesis consists of the gravitational force onlybut in general may include other forces related to Brownian motion if theparticle is very small. The position of a particle is given bydxpdt= up (4.30)Equations 4.30 and 4.29 are therefore coupled and need to be simultaneouslyintegrated in order to determine the particle trajectories. An explicit Eulertemporal discretization is applied to Eqn 4.29 to yieldun+1p =unp + ∆t(a+ unf/τ)1 + (∆t/τ)(4.31)while a trapezoidal discretization is applied to solve for xpxn+1p = xnp +∆t(unp + un+1p )2(4.32)534.4. Solution approachDuring the simulation, particles were introduced into the domain fromthe inlet from τ = 0 to τ = 0.5, with τ defined as the instantaneous timenormalised over one breath-period i.e. τ = t/T .Since the orientation of the gravity vector has been shown to have con-siderable influence on the particle trajectories and hence deposition, simu-lations were conducted for two cases where the gravity vector was orientednormal and tangential to the entrance of the solution domain respectively.4.4.5 Initial conditionsThe solutions were initialized with a zero velocity field and at atmo-spheric pressure mimicking conditions immediately preceding the start ofinspiration. The initial mass fraction of oxygen inside the alveolar ducts is0.15, which corresponds to the oxygen-poor condition that exists at the com-pletion of expiration. Each breath period was 3 seconds long correspondingto a breathing frequency of 20 breaths per minute.4.4.6 Solution hardwareThis solution approach was computed in parallel using a message-passinginterface (MPI) using 6 virtual cores of an Intel Xeon CPU with a clock-speed of 2.50 GHz such that the average simulation time for one completebreath period was approximately 4.5 hours.54Chapter 5Results and Discussion5.1 OverviewThis chapter details the results obtained from the simulations describedin Chapter 4. The results and the discussions that accompany them havebeen presented in a manner that meets the objectives detailed in Chapter 3.They have been divided into two sections each of which discusses the natureof the results obtained and then presents the effects of advancing emphysemaon those results. The first section discusses the oxygen transport and generalnature of the flow-field into the acinus, while the second details the natureof particle deposition dynamics in the respiratory zone and the effects ofemphysema on the observed deposition parameters.5.2 Flow field and oxygen transport5.2.1 Oxygen transportThe mass fraction of oxygen inside the computational domain increasesduring inspiration due to convective and diffusive transport of oxygen-richair from the inlet into the alveoli. To monitor the average oxygen concentra-tion within the three ducts of the computational domain, twelve sampling555.2. Flow field and oxygen transportpoints are arranged with four points within each duct. Each sampling pointis equidistant from the inlet of the domain and arranged in a radially sym-metrical fashion along the axis of the flow. As a result, the average of thereadings from the four sample points within a duct for each time step rep-resents a measure of oxygen transported to the alveolar spaces of that ductup to that time step. From these sampling points, the percentage increaseof oxygen mass fraction averaged over the four sampling points within eachduct, denoted C, can be defined asC =144∑n=1(Ct, n − CoCo)× 100% (5.1)where Ct, n is the mass-fraction of oxygen at time t at sampling point nand Co is the initial mass-fraction of oxygen in the domain at the startof inspiration. The temporal variation of C during inspiration for eachduct is plotted in Fig. 5.1. The trends obtained in Fig. 5.1 clearly showthat oxygen transport through all three ducts is severely reduced with theprogress of emphysematous destruction and the associated reduction in air-way wall motion and flow rate. This provides a qualitative explanation ofventilation-perfusion mismatch (conditions where there is an imbalance be-tween the oxygen transported into the lung and blood circulated throughit) difficulties that are commonly reported in patients with emphysema [66].Figure 5.1 also shows that there is a phse-shift in the time-signal of oxygenmass fraction between the proximal and distal generations of the compu-tational domain. The proximal generation (Duct 1) reaches a maximumC value at t/T = 0.3 while the distal generation (Ducts 2 and 3) have an565.2. Flow field and oxygen transportincreasing value of C until the end of inspiration. This is explained by thefact that Ducts 2 and 3 are distal to Duct 1 and thus receive oxygen laterin the inspiratory cycle.In order to determine the individual effects of septal destruction andthe decrease in airway wall motion and accompanying loss in inlet air-flowrate, parametric studies are performed in which these two effects are turnedon or off independently. Figure 5.2 plots the temporal variation of C forsimulations in which only septal degradation is present, and the airway wallmotion and inlet flow rates equal those in the healthy case. For, each duct,very little difference in the transport of oxygen into the alveoli is notedbetween the healthy and diseased cases despite the significant differencesin the degree of septal degradation between the four cases. Fig. 5.3 plotsthe temporal variation of C for simulations in which no septal degradationis present and only the emphysematous degradation of airway wall motionand flow rates are included. The trends noted in Fig. 5.3 mirror those inFig. 5.1 for which all features of the disease are included. This suggeststhat the factors primarily responsible for the decreasing oxygen transportinto emphysematous acini are the reduction in airway wall motion and flowrate; as septal destruction has a negligible impact on the reduced oxygentransport.5.2.2 Hydraulic lossesIn order to characterize the effect of emphysema on the flow-relatedpressure drop, termed the hydraulic loss, the average drop in static pressurethrough each alveolar duct is computed. While the total pressure is a more575.2. Flow field and oxygen transport0.0 0.1 0.2 0.3 0.4 0.5t/ T024681012CDuct IHealthyCase ICase IICase III0.0 0.1 0.2 0.3 0.4 0.5t/ T024681012CDuct II0.0 0.1 0.2 0.3 0.4 0.5t/ T024681012CDuct IIIFigure 5.1: Temporal variation of the percentage increase in oxygen con-centration within a duct (C) for cases with both septal destruction andemphysematous wall motion/inlet flow rates.585.2. Flow field and oxygen transport0.0 0.1 0.2 0.3 0.4 0.5t/ T024681012CDuct IHealthyCase ICase IICase III0.0 0.1 0.2 0.3 0.4 0.5t/ T024681012CDuct II0.0 0.1 0.2 0.3 0.4 0.5t/ T024681012CDuct IIIFigure 5.2: Temporal variation of the percentage increase in oxygen con-centration within a duct (C) for cases with only emphysematous septal de-struction.595.2. Flow field and oxygen transport0.0 0.1 0.2 0.3 0.4 0.5t/ T024681012CDuct IHealthyCase ICase IICase III0.0 0.1 0.2 0.3 0.4 0.5t/ T024681012CDuct II0.0 0.1 0.2 0.3 0.4 0.5t/ T024681012CDuct IIIFigure 5.3: Temporal variation of the percentage increase in oxygen con-centration within a duct (C) for cases with only emphysematous wall mo-tion/inlet flow rates.605.2. Flow field and oxygen transportappropriate metric for assessing hydraulic loss, the extremely low velocitieswithin the pulmonary acini generate dynamic pressures that are smallerthan 10−3 times the static pressure, and thus the static pressure is a veryclose approximation of the total pressure. The average pressure drop foreach duct is computed as the difference in static pressure between the oneinlet and the average of the two outlets during inspiration, and the averageof the two inlets and one outlet during expiration. The pressure values areobtained from probes located at the spatial center of the inlets and outletsof each duct, which were shown in Fig. 4.1.The temporal variation in the pressure drop across the three ducts duringinspiration is shown in Figs. 5.4-5.6. As in Figs. 5.1-5.3, Fig. 5.4 plots theresults for simulations in which both septal degradation and emphysematouswall motion and inlet flow rates are present, and Figs. 5.5 and 5.6 plot theresults for simulations with only septal destruction and only emphysematousdegradation in wall motion and flow rates, respectively. Figure 5.4 showsthat the pressure drop through all ducts varies almost sinusoidally duringinspiration with a maximum at the point of maximum inspiratory flow rateand approaching zero at the start and end of inspiration. This is consistentwith losses varying directly with the flow rate, i.e. velocity-squared losses.Comparing Duct 1 and Ducts 2 and 3 in Fig. 5.4, the lower pressure dropthrough Ducts 2 and 3 is because the flow rate reduces past the dichotomousbranching. It is notable that the pressure drops through Ducts 2 and 3are equal for the Healthy, Case I, and Case III cases, but for Case II, alarger pressure drop occurs across Duct 2 than Duct 3. Because Case II hasseptal destruction in Duct 2 but not in Duct 3, the flow resistance through615.2. Flow field and oxygen transportFigure 5.4: Temporal variation in the static pressure drop across the alve-olar ducts for cases with both septal destruction and emphysematous wallmotion/inlet flow rates.625.2. Flow field and oxygen transportFigure 5.5: Temporal variation in the static pressure drop across the alveolarducts for cases with only emphysematous septal destruction.635.2. Flow field and oxygen transportFigure 5.6: Temporal variation in the static pressure drop across the alveolarducts for cases with only emphysematous wall motion/inlet flow rates.645.2. Flow field and oxygen transportDuct 2 is lower than through Duct 3, resulting in the flow preferentiallytending towards Duct 2. The increased velocity through Duct 2 results inthe observed increased static pressure drop relative to Duct 3. The velocitiesthrough the ducts are discussed in more detail in section 5.1.4.Further insight into the mechanism behind hydraulic losses through thepulmonary acini is found by considering the independent effects of septaldestruction and emphysematous wall motion and flow rates. In Fig. 5.5, onlyseptal destruction is considered; the wall motion and flow rates equal those inthe healthy case. Comparing the healthy and diseased cases in Duct 1 showsthat the removal of septa reduces the static pressure drop through the duct.This result is somewhat counter intuitive, as emphysema is associated withincreased difficulty in breathing. It indicates that destruction of the intra-alveolar septa makes the duct less resistive to the flow and thus produceslower hydraulic loss through the duct. In Ducts 2 and 3, there is no septaldestruction for Case I, hence it exactly follows the healthy case, while thedifference between Duct 2 and 3 in Case II is caused by the asymmetryin the flow resistance and velocity between the ducts, noted earlier. Incontrast, when only the emphysematous wall motion and inlet flow rates areconsidered (Fig. 5.6), the pressure drop through the alveolar ducts is muchhigher than in the cases where only septal destruction occurs. Within eachduct, a steady reduction in the pressure drop is noted between the healthyand diseased cases, consistent with the reduction in inlet flow rates for eachcase. This suggests that hydraulic losses through the emphysematous aciniare influenced primarily by the loss of intra-alveolar septa. The loss in septareduces the hydraulic losses through the alveolar ducts, and the reduced655.2. Flow field and oxygen transportflow rates and wall motion due to emphysema has a secondary importanceon the hydraulic loss. Similar results were also obtained for the expiratorybreathing phase.5.2.3 Flow field visualizationThe complexity of the spatial geometry and the temporally-varying na-ture of the wall motion and inlet flow rates makes the complete charac-terization of the flow field through the acinus very challenging. Therefore,attention is given primarily to the patterns that develop in the ducts andalveolar spaces at peak inspiration and expiration.Figures 5.7 and 5.8 plots streamlines coloured by the velocity magnitude(U) in planes that bisect the computational domain at peak inspiration andexpiration, respectively.It may be noted that in Case II (Figs. 5.7(c) and 5.8(c)), the velocities inDuct 2 with emphysematous destruction are much higher than Duct 3 withhealthy septa. This is because the removal of the septa decreases overallresistance to the flow, which in turn causes the higher velocities in Duct2. It should not be concluded that emphysematous destruction leads to anincrease in flow rate through the acinus. As previously described, the in-spiratory and expiratory flow-rates are coupled to the negative and positivepleural pressures generated in the acini during breathing by the parenchy-mal mesh that surrounds the acinus. The degeneration of the parenchyma inemphysema results in reduced pleural pressures and correspondingly lowerflow rates. As a result, emphysema reduces flow resistance through the air-way via septal destruction while also simultaneously reducing the overall665.2. Flow field and oxygen transportFigure 5.7: Velocity contours at peak inspiration (a) Healthy case (b) Case I(c) Case II (d) Case III. Arrows indicate general flow direction in the acinus.675.2. Flow field and oxygen transportFigure 5.8: Velocity contours at peak expiration (a) Healthy case (b) Case I(c) Case II (d) Case III. Arrows indicate general flow direction in the acinus.685.2. Flow field and oxygen transportFigure 5.9: Streamlines in an individual alveolus from Duct 1 near the endof inspiration for the healthy and diseased cases. Streamlines are colouredaccording to velocity magnitude.flow rates via degeneration of the elastic parenchymal tissue. From Figs. 5.7and 5.8 it is apparent that velocities within the cavities are on average oneor two orders lower in magnitude than the main flow in the duct, whichalso accords well with previously published results [8, 39]. As such, due totheir low inertia, the flow patterns in the alveolus are highly susceptible tochange based on the conditions in the duct. This is illustrated in Fig. 5.9for the alveolus labelled in Fig. 5.7, by plotting the streamlines near the endof inspiration coloured by the velocity magnitude. This particular alveolusis highlighted because it is oriented such that the air flows tangentially. Arecirculatory centre occurs in the diseased cases but is absent in the healthy695.2. Flow field and oxygen transportcase. As recirculating flows possess stagnation points that result in irre-versible paths of massless particles [9, 38], the onset of recirculatory flowhas a significant influence on the transport into the alveoli. In literature,the onset of recirculatory flow in the alveoli is determined based on the ra-tio of the flow rate into the alveolus (Qa) and the flow rate through theduct (Qd). (Qd) is termed the shear flow [9, 43]. Quasi-steady recirculatorycentres like those noted in the diseased cases in Fig. 5.9 occur for relativelysmall Qa/Qd ratios, which typically occur in the more proximal and medialgenerations of the acinar tree [10]. In the more distal generations, this ratioincreases and the flow patterns gradually become more radial, resemblingthe healthy case in Fig. 5.9. In addition, literature has found that includ-ing wall motion into the numerical simulations increases the relative radialvelocity component, resulting in more radial flow patterns [9].From the above considerations, the following deductions can be maderegarding Fig. 5.9. In the healthy case, the radial flow magnitude is relativelystrong due to the presence of the relatively larger airway wall motion andthus recirculation is not observed. In the diseased cases, two simultaneouseffects are present. The removal of alveolar septa results in a decrease inflow resistance within the duct and a concomitant increase in the shear flowcomponent within the alveolus. Additionally, the weakening of wall motioncauses a reduction in radial airflow into the alveoli. The combination ofthese two effects results in the noted formation of recirculatory centres inthe alveolus of the diseased cases. The acinar flow field is an importantfactor in describing the transport of aerosols and fine particulate matter.Therefore, these changes in the flow-field will have a considerable impact on705.3. Particle transportthe fates of particles inhaled into diseased lungs, considered in the followingsection.5.3 Particle transport5.3.1 Simulation conditions and boundary conditionsIn order to characterize the deposition of particles in healthy and em-physematous geometries, the deposition fraction (f) quantity is defined asfollowsf =(total number of particles deposited)(total number of particles released into the duct)(5.2)In the context of this thesis a particle is assumed to have been depositedwhen its path intersects the wall of the computational domain. In orderto characterize the spatial variation of particle deposition within the duct,the deposition fraction is computed separately for each duct. The resultingvalues are plotted against the normalised breath cycle τ for the cases ofadvancing emphysema, described in the previous section. The simulationsare conducted for two particle sizes with diameters of 1 and 3 µm. Thesesizes represent the diameter range where particles in acinar flows are affectedby a combination of the flow-field and gravity. Particles smaller than 1 µmare primarily affected by diffusive forces within the flow (e.g. Brownianmotion) while particles larger than 3 µm in the acinus have their fates whollydetermined by the orientation of the gravitational vector [11].The simulations are also conducted for two separate orientations of the715.3. Particle transportTable 5.1: Simulation matrix for particle deposition study1 µm 3 µm 1 µm 3 µmHealthy Normal Normal Tangential TangentialCase I Normal Normal Tangential TangentialCase II Normal Normal Tangential TangentialCase III Normal Normal Tangential Tangentialgravitational vector relative to the incoming airflow into the computationaldomain. This is done to study the relative importance of the orientationof gravitational vector versus the airflow on the particle trajectories. Oneorientation is normally directed co-axially with the inflow boundary of Duct1, while the second is directed tangentially to the same space. The parti-cle transport simulations are conducted for the four cases (Healthy, Case I,Case II and Case III) presented in earlier chapters. For presentation pur-poses, the results are grouped into four categories according to particle sizeand orientation of the gravity vector: 1 µm Normal, 3 µm Normal, 1 µmTangential and 3 µm Tangential. This is shown in a graphical format inTable 5.1 as a simulation matrix. The simulations are conducted for fivesuccessive breath cycles while the particle trajectories and the depositionfractions were monitored at the end of each cycle.5.3.2 Temporal variation of deposition fractionFigures 5.11-5.14 plot the deposition fraction for Ducts 1, 2 and 3 againstthe normalised respiratory breath period τ . In general, these plots of f vsτ follow a similar pattern in which the majority of the deposition occursin the first few breath-cycles, after which f achieves a steady state value.Beyond this point, most of the particles in the domain have deposited, and725.3. Particle transportFigure 5.10: Orientation of the gravity vector with respect to the computa-tional geometry for the normal and tangential cases respectivelythe balance having been convected out of the computational domain.From Figures 5.11-5.14 it is observed that for all cases, the greatest val-ues of f were in Duct I, indicating that a majority of particles deposit inproximal generations of the acinus. These values of f in Duct I increasewith advancing emphysema across all particle sizes and gravity vector ori-entations. Although Identical values of f are observed in Ducts II and IIIin a majority of the cases, a notable exception to this exists in scenar-ios which have emphysematous destruction corresponding to Case II. ThisCase is unique in that it is the only arrangement where the geometry is non-symmetrical about the principal axis of the model. Thus, while Duct II andIII are in the same generation, only Duct II has emphysematous destructionwhile Duct III does not. It is observed that for all scenarios, at the Case IIlevel, the values of f are consistently greater in Duct II than in Duct III.As expected, particle size plays an important role in determining de-position patterns. By comparing Fig. 5.11 and 5.12 it is clear that the 3µm particles deposit completely within one breath cycle, while 1 µm parti-cles remain in suspension longer with their f values stabilizing around 2-2.5735.3. Particle transport0 1 2 3 4 FractionDuct IHealthyCase ICase IICase III0 1 2 3 4 FractionDuct II0 1 2 3 4 FractionDuct IIIτττFigure 5.11: Deposition fraction of 1 µm particles, with the gravity vectororiented in the normal direction745.3. Particle transport0 1 2 3 4 FractionDuct I0 1 2 3 4 Fraction Duct IIHealthyCase ICase IICase III0 1 2 3 4 FractionDuct IIIτττFigure 5.12: Deposition fraction of 3 µm particles, with the gravity vectororiented in the normal direction755.3. Particle transport0 1 2 3 4 5τ0. FractionDuct IHealthyCase ICase IICase III0 1 2 3 4 5τ0. FractionDuct II0 1 2 3 4 5τ0. FractionDuct IIIFigure 5.13: Deposition fraction of 1 µm particles, with the gravity vectororiented in the tangential direction765.3. Particle transport0 1 2 3 4 5τ0. FractionDuct IHealthyCase ICase IICase III0 1 2 3 4 5τ0. FractionDuct II0 1 2 3 4 5τ0. FractionDuct IIIFigure 5.14: Deposition fraction of 3 µm particles, with the gravity vectororiented in the tangential direction775.3. Particle transportFigure 5.15: Particle deposition efficiencies under differing gravitational ori-entations (Adapted with permission from [11])breath cycles. Cases with 3 µm particles are also observed to have a highersteady-state f value than their corresponding 1 µm particle cases.The role of the orientation of the gravity vector is also significant. Thecases with the gravitational vector oriented normally to the incoming flow(Figs. 5.11 and 5.12) show a greater value of f in Ducts 2 and 3 than forthe corresponding tangential cases (Figs. 5.13 and 5.14).Figures 5.11 - 5.14 also conform closely to previous simulations conductedby Sznitman et. al. [11] which have been reproduced here in Fig 5.15. Whilethe deposition fractions for both particle sizes obtained by Sznitman et. al.are much higher, this can be explained bue to differences in geometry andboundary conditions. The Sznitman et. al. simulations tracked the deposi-785.3. Particle transporttion fraction in the third generation of the acinus and with a higher initialvelocity in the boundary conditions. In the presently conducted simulationthe deposition fraction was tracked in ducts lower in the acinar generationaltree with differing geometry and boundary conditions.5.3.3 Variation of steady state deposition fractionPlots showing the variation in steady-state deposition fraction in thethree ducts, for all of the simulated cases are shown in Figs. 5.16-5.19. Theseplots may be interpreted by keeping in mind the two competing forces ofconvective flow and gravity as given by Equation 2.3. The convective forceand gravity reinforce each other when they are aligned as in the normalcases. In the tangential cases the gravity vector is perpendicular to thebulk incoming flow, thus directing part of it towards the acinar walls. Thiseffect is magnified in the emphysematous cases where the convective flow isreduced.1 µm NormalFrom Fig 5.16, it is immediately apparent that for the 1 µm particleswith a normal gravity vector the steady-state deposition fraction in Duct Iincreases with increasing levels of emphysema. While the normal case showsa f value of only 0.15, this increases in Case III to 0.32. The depositionfraction values also increase uniformly in Ducts II and III for the healthy,Case I and Case III simulations. In Case II it is observed that Duct II has ahigher steady state deposition fraction (0.1) than Duct III (0.05). Since inCase II, Duct II has septal degradation while Duct III does not, this leads795.3. Particle transportHealthy Case I Case II Case III0. FractionIII IIIIII IIIIIIIIIIII IIIFigure 5.16: Plot of final deposition fraction after five breath cycles for 1µm particles, gravity vector in normal direction. The filled circles indicatehealthy acinar duct while empty circles indicate emphysematous destruction.The labels indicate the duct number.to the conclusion that septal degradation preferentially enhances particulatedeposition between two ducts of the same generation.3 µm NormalFig 5.17 shows steady state deposition fraction values for 3 µm particleswith the gravity vector oriented in the normal direction. While the steadystate f in Duct I increase with advancing emphysema as with the 1 µmparticle cases, it should be noted that they are much higher than the cor-responding cases. The Healthy case has a steady state deposition fractionvalue of 0.3 while the corresponding value for Case III is 0.6. This may beattributed to the greater effect of gravity on the larger and hence heavierparticles. It should also be noted that a greater proportion of particles de-posit in Duct I than Ducts II and III for the diseased cases, than for the805.3. Particle transportcorresponding simulations with 1 µm particles. In fact, the deposition inthe proximal duct is so greatly enhanced, that the more distal Ducts II andIII have lower deposition fraction values than the ones corresponding to theHealthy case under identical simulation conditions. It should also be notedthat the increase in deposition fraction for Duct I is not uniform as in the1 µm case. In particular Case I and Case II have very similar steady statevalues of f for Duct I (0.51 and 0.52 respectively). This may be attributedto the fact that in Case II. while the deposition is increased due to emphy-sematous degradation, the normal direction of the gravity vector -collinearto the incoming flow- and its increased effects due to the larger particle sizedirect it away from the walls and eventually deeper into the acinus. Thuslowering the deposition fraction in Duct I making it almost comparable toCase I. Case III does show an increase in deposition fraction for Duct I ascompared to Case I and II. This may be attributed to a further loweringof the flow rate, which prevents the particles from being convected into themore distal generations and deposits them at the branching point betweenDucts I, II and III. It should also be noted that Duct II for Case II showsgreater deposition than Duct III for the same case as was observed in theprevious subsection.815.3. Particle transportHealthy Case I Case II Case III0. FractionIII IIIIII IIIIIIIIIIII IIIFigure 5.17: Plot of final deposition fraction after five breath cycles for 3µm particles, gravity vector in normal direction. The filled circles indicatehealthy acinar duct while empty circles indicate emphysematous destruction.labels indicate the duct number.1 µm TangentialFig. 5.18 shows the steady state f values for 1 µm particles with thegravity vector oriented in the tangential direction. On comparing with thesimulations with an identical particle size but a normal gravity vector, it isseen that the deposition fraction values in Duct I are consistently greater forthe present cases. This is explained by the fact that the tangential gravityvector directs the particles towards the walls of the acinus instead of alongthe bulk flow, thereby enhancing the effects of emphysematous destruction.This also leads to a decrease in deposition fraction in the more distal ductsas compared to the healthy case. The pattern of increased deposition inDuct II as compared to Duct III in Case II is also observed here.825.3. Particle transportHealthy Case I Case II Case III0. FractionIII IIIIII IIIIIIIIIIII IIIFigure 5.18: Plot of final deposition fraction after five breath cycles for 1 µmparticles, gravity vector in tangential direction. The filled circles indicatehealthy acinar duct while empty circles indicate emphysematous destruction.labels indicate the duct number.3 µm TangentialFig. 5.19 shows the steady state f values for 3 µm particles with thegravity vector oriented in the tangential direction. These cases show thehighest deposition fraction of all the cases due to the twin effects of a largeparticle size increasing the effect of the gravitational vector and the tan-gential direction of the vector that enhances deposition on the walls awayfrom the bulk flow. There’s a marked increase in deposition faction valuesbetween the Healthy case and Case III (0.5 and 0.75 respectively). Howeverthere is very little increase in deposition fraction values for Cases I, II andIII. This is because the gravitational effects completely overwhelm the con-vective effects and since the geometrical conditions in Duct I are the samefor Case I, II and III, they have identical deposition fraction values. The835.3. Particle transportHealthy Case I Case II Case III0. FractionIII IIIIII IIIIIIIIIIII IIIFigure 5.19: Plot of final deposition fraction after five breath cycles for 3 µmparticles, gravity vector in tangential direction. The filled circles indicatehealthy acinar duct while empty circles indicate emphysematous destruction.labels indicate the duct number.pattern of increased deposition in Duct II as compared to Duct III in CaseII is also observed here.The results presented in this section particularly with reference to theincreased deposition of particles within the emphysematous ducts followclosely with the results presented by Oakes et. al. [67]. However, there isconsiderable debate on this point as a more recent study [60] shows thatthe deposition in emphysematous geometries in enhanced in the healthyportions of a diseased model. This difference between experimental andnumerical simulations requires further investigation and it is proposed thatlarger models of the complete acinus incorporating the progressive nature ofemphysema as has been presented here might help to discover the mecha-nism responsible for enhancement or diminishment of particle deposition inemphysematous models.84Chapter 6Summary and Conclusions6.1 Model developmentA CFD model has been developed, to simulate the effects of COPD,specifically early-stage emphysema in the airways of the pulmonary acinus.The airway structure was based on a mathematical model developed by Fung[7], using a 14-hedrons geometry to provide a close packed alveoli structuresurrounding a central duct was achieved by the arrangement of 14-hedrons.The computational geometry consisted of two generations - the fourth andfifth- of the acinus. The fourth generation consisted of a single duct (Duct1) which branched into two identical ducts (Ducts 2 and 3) that formed thefifth generation. During inspiration, the incoming air entered via a singleentryway into Duct 1 and exited through Ducts 2 and 3 each of which hadtwo exits. During expiration this process was reversed. The computationaldomain underwent a self-similar time-dependent sinusoidal motion of ex-pansion and contraction, which was implemented in order to simulate thereal-world motion of airways that provides the pressure differential that isresponsible for the inspiratory and expiratory airflow. The boundary con-ditions in the model were adjusted so as to mimic the appropriate air-flowrates. The effects of emphysema were incorporated by two methods. The856.2. Oxygen transport and flow-fieldemphysematous destruction of parenchymal tissue was modelled as septaldegradation in the ducts. The progressive nature of the disease was alsoincorporated with the inclusion of successive cases where the septal destruc-tion was present in one, two and three of the ducts in the domain. Theloss of the parenchymal mesh that normally provides the elastic recoil nec-essary for the proper functioning of the airways was also incorporated intothe disease model by changing the amplitude of sinusoidal motion of theairway domain. This measure is also complemented by a corresponding lossin flow-rate that is observed under actual physiological conditions.6.2 Oxygen transport and flow-fieldThe oxygen transport through the computational domain was tracked bysolving a scalar transport diffusion equation. A parameter C was defined inorder to quantify the change in concentration of oxygen over a single inhala-tion. Tracking this parameter over the different, progressive cases of emphy-sema, it was observed that increasing levels of emphysema corresponded toa decrease in the parameter across all three ducts. It was further observedthat this behavior was mainly caused by a decrease in airway wall motionand the corresponding loss in flow rate rather than the septal degradation.This is readily explained by the fact that while the majority of the gaseoustransport occurs by diffusion rather than advection, the loss in airflow ratedue to disease is sufficient enough to decrease the oxygen levels by almost 5% in the most extreme cases.The effect of the disease on the flow-field was quantified through the866.3. Particle transportstatic pressure-drop across each of the ducts for the different cases. It wasobserved that both the septal destruction and the decrease in airway wallmotion and its attendant loss in flow-rate both led to an overall static pres-sure drop. However, the airflow was preferentially increased through thoseregions which had septal destruction. While this may seem to indicate thatacinar generations with septal destruction might be better ventilated thantheir healthy counterparts, it should be remembered that septal destructionis usually accompanied by a loss in static pressure recoil which decreases theavailable pressure driving the breath-cycle resulting in an overall decreasein flow-rate and consequently poorer ventilation.The diseased conditions were also observed to have a direct impact onthe temporal dynamics of the flow-field over a single breath cycle. Whilere-circulatory flow motions are usually observed in acinar generations prox-imal to those included in this simulation (up to the 3rd generation of theacinus), it was noted that with the inclusion of the effects of emphysema,a recirculatory centre developed at the proximal end of the alveolus duringinspiration, which increased in intensity in proportion to the extent of thedisease.6.3 Particle transportIn order to study the effects of the disease on particle transport into thelung, simulations were conducted using two different particle sizes (1 µmand 3 µm). At these sizes the principal forces affecting the particles areconvective forces due to the flow and gravitational force. As a result, two876.3. Particle transportsets of simulations were conducted with the gravity vector being oriented ina direction normal and tangential to the flow at the domain entrance respec-tively. Particle deposition was quantified by the particle deposition fractionwhich was the number of particles deposited divided by the total numberintroduced into the domain during the initial inspiration. The particle de-position fraction was tracked in the computational domain for a total of fiveconsecutive breath cycles.The results show that the deposition fraction was increased in emphyse-matous ducts due to both a decrease in airway wall motion, and emphyse-matous destruction, thereby indicating that the earliest stages of the diseaseenhanced depositions in the ducts. This results in in an increased risk for dis-ease propagation. The increase in deposition fraction in ducts distal to thefirst emphysematous ducts depended on the particle size and orientation ofthe gravitational vector. Smaller particles and a gravitational vector alignedwith the the flow tend to penetrate and deposit deeper into the acinus, whilethe larger particles and a tangential gravitational vector tended to promotethe deposition of particles in the proximal ducts of the acinus. It was alsoobserved that for two ducts in the same generation, particles would pref-erentially deposit in the one with emphysematous destruction as comparedto one without. Thus, while emphysema notably enhanced the depositionin the first duct encountered by particles, the deposition patterns in thesubsequent ducts depended on particle size, gravitational vector orientationand the presence of emphysematous destruction in the distal airways.886.4. Summary of conclusions6.4 Summary of conclusions• The onset of emphysema causes a significant drop in oxygen transportwithin the acinus.• This loss of oxygen transported to the acinus during inhalation is pri-marily due to the loss in elastic recoil which causes a related drop inincoming airflow rate.• Emphysema causes a drop in hydraulic losses within the acinus withthe flow resistance dropping appreciably in ducts affected by septaldestruction.• The effects of emphysema changes the flow field within the alveolarcavity with the appearance of a recirculatory centre with the alveoliin the diseased cases indicating an increase in shear flow and decreasein flow into the alveolar cavity.• While particle size and the orientation of the gravitational vector playan important part in deposition within the acinus, emphysematousdestruction increases the particle deposition in any duct across allparticle sizes and orientations. This is due to both the decrease inairway wall motion and emphysematous septal destruction.• For two ducts in the same alveolar generation, the duct with septaldestruction shows greater deposition regardless of particle size or grav-itational vector orientation.• Smaller particles with a gravity vector oriented coaxially to the incom-896.5. Future working flow allow particles to penetrate and deposit deep in the acinus,while larger particles with other gravitational orientations tend to de-posit in more proximal generations6.5 Future workWhile this study was successful in implementing a model of emphysemaand applying to an isolated section of the pulmonary acinus, it should benoted that the complete structure of the airway tree is extremely heteroge-neous. As a result, reliable predictions of the effects of emphysema on thefactors investigated in this study over multiple airway generations will re-quire the development of a larger computational domain spanning the entirepulmonary acinus or even perhaps multiple acini. Incorporation of realisticpatient data to set appropriate flow rate boundary conditions for varyingstages of emphysema is also an additional refinement that should be con-sidered. In short, therefore, future works along these lines must incorporatemultiple airway generations and a larger number of cases with varying lev-els of emphysema, in order to make accurate predictions of the effects of adisease that affects the extremely heterogeneous human pulmonary airwaysystem.90Bibliography[1] C Kleinstreuer and Z Zhang. Airflow and Particle Transport in theHuman Respiratory System. Annual Review of Fluid Mechanics, 42(1):301–334, 2010. → pages xi, 3, 11, 21[2] W Mitzner. Emphysema A Disease of Small Airways or LungParenchyma? New England Journal of Medicine, 365(17):1637–1639,2011. → pages xi, 7[3] ER Weibel, B Sapoval, and M Filoche. Design of Peripheral Airwaysfor Efficient Gas Exchange. Respiratory physiology & neurobiology, 148(1):3–21, 2005. → pages xi, xii, 12, 19[4] DL Fry and RE Hyatt. Pulmonary Mechanics* A Unified Analysis ofthe Relationship Between Pressure, Volume and Gas Flow in the Lungsof Normal and Diseased Human Subjects. The American Journal ofMedicine, 29, 1960. → pages xi, 13[5] JB West. Pulmonary Pathophysiology : The Essentials. LippincottWilliams and Wilkins, 2005. → pages xi, 7, 14, 16, 44[6] ER Weibel, AF Cournand, and DW Richards. Morphometry of the91BIBLIOGRAPHY BibliographyHuman Lung With a Foreword. Springer Berlin Heidelberg, 1963. →pages xii, 10, 11, 17, 19[7] Y C Fung. A Model of the Lung Structure and its Validation. Journalof Applied Physiology, 1988. → pages xii, 20, 21, 28, 36, 85[8] H Kumar, MH Tawhai, EA Hoffman, and CL Lin. The Effects ofGeometry on Airflow in the Acinar Region of the Human Lung. Journalof Biomechanics, 42:1635–1642, 2009. → pages xii, 20, 24, 25, 36, 69[9] A Tsuda, FS Henry, and JP Butler. Chaotic Mixing of AlveolatedDuct Flow in Rhythmically Expanding Pulmonary Acinus. Journal ofApplied Physiology, 1995. → pages xii, 19, 24, 25, 70[10] P Hofemeier, R Fishler, and J Sznitman. The role of respiratory flowasynchrony on convective mixing in the pulmonary acinus. Fluid Dy-namics Research, 46(4):041407, 2014. → pages xii, 26, 36, 70[11] J Sznitman, T Heimsch, JH Wildhaber, A Tsuda, and T Ro¨sgen. Res-piratory Flow Phenomena and Gravitational Deposition in a Three-Dimensional Space-Filling Model of the Pulmonary Acinar Tree. Jour-nal of Biomechanical Engineering, 2009. → pages xii, xv, 20, 28, 29,30, 36, 43, 71, 78[12] J Vestbo, SS. Hurd, AG. Agusti, PW. Jones, C Vogelmeier, A Anzueto,PJ Barnes, LM Fabbri, FJ Martinez, M Nishimura, et al. Global strat-egy for the diagnosis, management, and prevention of chronic obstruc-tive pulmonary disease: Gold executive summary. American Journal92BIBLIOGRAPHY Bibliographyof Respiratory and Critical Care Medicine, 187(4):347–365, 2013. →pages 1, 2, 3, 4, 5[13] N Mittmann, L Kuramoto, SJ Seung, JM Haddon, C Bradley-Kennedy,and JM Fitzgerald. The cost of moderate and severe copd exacerbationsto the canadian healthcare system. Respiratory Medicine, 102(3):413–421, 2008. → pages 1[14] SD Shapiro. The pathogenesis of emphysema: The elastase antielastasehypothesis 30 years later. Proceedings of the Association of AmericanPhysicians, 107(3):346–352, 1995. → pages 3[15] JL Wright and A Churg. Advances in the Pathology of COPD.Histopathology, 49(1):1–9, jul 2006. → pages 3[16] DM Mannino, AS Buist, TL Petty, PL Enright, and SC Redd. Lungfunction and mortality in the united states: Data from the first nationalhealth and nutrition examination survey follow up study. Thorax, 58(5):388–393, 2003. → pages 4[17] JC Hogg and S Van Eeden. Pulmonary and systemic response to at-mospheric pollution. Respirology, 14(3):336–346, 2009. → pages 4[18] P Lee, TR Gildea, and JK Stoller. Emphysema in nonsmokers: Alpha1-antitrypsin deficiency and other causes. Cleveland Clinic Journal ofMedicine, 69(12):928–9, 2002. → pages 4[19] NA Zwar, GB Marks, O Hermiz, S Middleton, EJ Comino, I Hasan,S Vagholkar, and SF Wilson. Predictors of accuracy of diagnosis of93BIBLIOGRAPHY Bibliographychronic obstructive pulmonary disease in general practice. MedicalJournal of Australia, 195(4):168–71, 2011. → pages 4[20] TMA Wilkinson, GC Donaldson, JR Hurst, TAR Seemungal, andJA Wedzicha. Early therapy improves outcomes of exacerbations ofchronic obstructive pulmonary disease. American Journal of Respira-tory and Critical Care Medicine, 169(12):1298–1303, 2004. → pages5[21] DB Coultas, D Mapel, R Gagnon, and EVA Lydick. The Health Impactof Undiagnosed Airflow Obstruction in a National Sample of UnitedStates Adults. American Journal of Respiratory and Critical CareMedicine, 164(3):372–377, aug 2001. → pages 5[22] S Leung, AE Gallagher, V Kvetan, and LA Eisen. Barriers to Ultra-sonography Training in Critical Care Medicine Fellowships: A Surveyof Program Directors. Chest, 136(4 MeetingAbstracts):49S–50S, 2009.doi: 10.1378/chest.09-0553. → pages 6[23] R Pellegrino, G Viegi, V Brusasco, RO Crapo, F Burgos, R Casaburi,A Coates, CPM Van Der Grinten, P Gustafsson, J Hankinson,R Jensen, DC Johnson, N MacIntyre, R McKay, MR Miller, D Nava-jas, OF Pedersen, and J Wanger. Interpretative Strategies for LungFunction Tests. European Respiratory Journal, 26(5):948–968, 2005. →pages 6[24] WM Thurlbeck. Chronic Airflow Obstruction in Lung Disease. WBSaunders Company, 1976. → pages 694BIBLIOGRAPHY Bibliography[25] GL Snider. Emphysema: The First Two Centuries-and Beyond. A His-torical Overview, With Suggestions of Future Research: Part 2. Amer-ican Review of Respiratory Disease, 146(6):1615–22, 1992. → pages6[26] GL Snider, J Kleinerman, WM Thurlbeck, and ZH Bengali. The Def-inition of Emphysema. American Review of Respiratory Disease, 132(1):182–185, jul 1985. → pages 6[27] Stephen Milne and Gregory G King. Advanced imaging in copd: in-sights into pulmonary pathophysiology. Journal of thoracic disease, 6(11):1570, 2014. → pages 8[28] Julia Ley-Zaporozhan, Sebastian Ley, and Hans-Ulrich Kauczor. Mor-phological and functional imaging in copd with ct and mri: present andfuture. European Radiology, 18(3):510–521, 2008. → pages 8[29] D Aykac, EA Hoffman, G Mclennan, and JM Reinhardt. Segmentationand Analysis of the Human Airway Tree from Three-Dimensional X-ray CT Images. IEEE Transactions on Medical Imaging, 22(8):940–950,Aug 2003. → pages 10[30] K Horsfield and G Cumming. Morphology of the bronchial tree in man.Journal of Applied Physiology, 24(3):373–383, 1968. → pages 10[31] Design of Peripheral Airways for Efficient Gas Exchange. In RespiratoryPhysiology and Neurobiology, 2005. → pages 10[32] B Haefeli-Bleuer and ER Weibel. Morphometry of the Human Pul-95BIBLIOGRAPHY Bibliographymonary Acinus. The Anatomical Record, 220(4):401–414, 1988. →pages 11, 17, 23[33] ER Weibel. The Pathway for Oxygen. Harvard University Press, 1984.→ pages 12[34] JH Mateika and J Duffin. A Review of the Control of Breathing DuringExercise. European Journal Applied Physiology, 71:1–27, 1995. → pages12, 13[35] DE O’donnell and P Laveneziana. Physiology and consequences of lunghyperinflation in copd. European Respiratory Review, 15(100):61–67,2006. → pages 14[36] RE Hyatt. Expiratory flow limitation. Journal of Applied Physiology,55(1):1–7, 1983. → pages 14[37] SV Dawson and EA Elliott. Wave-speed limitation on expiratory flow-aunifying concept. Journal of Applied Physiology, 43(3):498–515, 1977.→ pages 14[38] FS Henry, JP Butler, and A Tsuda. Kinematically Irreversible AcinarFlow: A Departure from Classical Dispersive Aerosol Transport Theo-ries. Journal of Applied Physiology, 92(2):835–845, 2002. → pages 19,25, 70[39] J Sznitman, F Heimsch, T Heimsch, D Rusch, and T Ro¨sgen. Three-Dimensional Convective Alveolar Flow Induced by Rhythmic Breathing96BIBLIOGRAPHY BibliographyMotion of the Pulmonary Acinus. Journal of Biomechanical Engineer-ing, 2007. → pages 19, 24, 25, 41, 43, 69[40] A Tsuda, JP Butler, and JJ Fredberg. Effects of Alveolated Duct Struc-ture on Aerosol Kinetics II. Gravitational Sedimentation and InertialImpaction. Journal of Applied Physiology, 76(6):2510–2516, 1994. →pages 20, 23[41] P Hofemeier and J Sznitman. Role of alveolar topology on acinar flowsand convective mixing. Journal of Biomechanical Engineering, 136(6):061007, 2014. → pages 20, 36[42] TJ Pedley. Pulmonary Fluid Dynamics. Annual Review of Fluid Me-chanics, 9(1):229–274, 1977. → pages 21[43] J Sznitman. Respiratory Microflows in the Pulmonary Acinus. Journalof Biomechanics, 46:284–298, 2013. → pages 21, 22, 23, 26, 36, 70[44] CD Murray. The Physiological Principle of Minimum Work: I. TheVascular System and the Cost of Blood Volume. Proceedings of theNational Academy of Sciences of the United States of America, 12(3):207–214, 1926. → pages 23[45] MR Davidson and JM Fitz-Gerald. Flow Patterns in Models of SmallAirway Units of the Lung. Journal of Fluid Mechanics, 52(01):161–177,1972. → pages 23[46] HK Chang, RT Cheng, and LE Farhi. A Model Study of Gas Diffusion97BIBLIOGRAPHY Bibliographyin Alveolar Sacs. Respiration Physiology, (18):386–397, 1973. → pages23[47] WJ Federspiel and JJ Fredberg. Axial Dispersion in Respiratory Bron-chioles and Alveolar Ducts. Journal of Applied Physiology, 64(6):2614–2621, 1988. → pages 23[48] CP Yu and S Rajaram. Diffusional Deposition of Particles in a ModelAlveolus. Journal of Aerosol Science, 9(6):521–525, 1978. → pages 23[49] C Darquenne and M Paiva. Two- and Three-Dimensional Simulationsof Aerosol Transport and Deposition in Alveolar Zone of Human Lung.Journal of Applied Physiology, 80(4):1401–1414, apr 1996. → pages 23,24[50] C Darquenne. A Realistic Two-Dimensional Model of Aerosol Transportand Deposition in the Alveolar Zone of the Human Lung. Journal ofAerosol Science, 32(10):1161–1174, 2001. → pages 23[51] M Horner, G Metcalfe, S Wiggins, and JM Ottino. Transport Enhance-ment Mechanisms in Open Cavities. Journal of Fluid Mechanics, 452:199–229, 2002. → pages 24[52] B Ma, V Ruwet, P Corieri, R Theunissen, M Riethmuller, and C Dar-quenne. CFD Simulation and Experimental Validation of Fluid Flowand Particle Transport in a Model of Alveolated Airways. Journal ofAerosol Science, 40(5):403–414, 2009. → pages 2498BIBLIOGRAPHY Bibliography[53] JM Ottino. The Kinematics of Mixing: Stretching, Chaos, and Trans-port, volume 3. Cambridge university press, 1989. → pages 24[54] PD Anderson, OS Galaktionov, GWM Peters, FN Van de Vosse, andHEH Meijer. Analysis of Mixing in Three-Dimensional Time-PeriodicCavity Flows. Journal of Fluid Mechanics, 386:149–166, 1999. → pages24[55] WH Finlay. The Mechanics of Inhaled Pharmaceutical Aerosols: anIntroduction. Academic Press, 2001. → pages 27[56] A Tsuda, FS Henry, and JP Butler. Particle Transport and Deposition:Basic Physics of Particle Kinetics. Comprehensive Physiology, 3(4):1437–1471, 2013. → pages 27[57] JM Oakes, AL Marsden, C Grandmont, S C Shadden, C Darquenne,and IE Vignon-Clementel. Airflow and particle deposition simulationsin health and emphysema: From in vivo to in silico animal experiments.Annals of Biomedical Engineering, 42(4):899–914, 2014. → pages 31[58] EJ Berg. Stereoscopic particle image velocimetry analysis of healthyand emphysemic acinus models. Journal of Biomechanical Engineering,2010. → pages 31[59] P Aghasafari, IBM Ibrahim, and RM Pidaparti. Investigation of theeffects of emphysema and influenza on alveolar sacs closure through cfdsimulation. Journal of Biomedical Science and Engineering, 9(06):287,2016. → pages 3199BIBLIOGRAPHY Bibliography[60] Jessica M. Oakes, Philipp Hofemeier, Irene E. Vignon-Clementel, andJosue´ Sznitman. Aerosols in healthy and emphysematous in silico pul-monary acinar rat models. Journal of Biomechanics, 49(11):2213–2220,2015. → pages 31, 84[61] Bernard Sapoval, M Filoche, and ER Weibel. Smaller is betterbut nottoo small: a physical scale for the design of the mammalian pulmonaryacinus. Proceedings of the National Academy of Sciences, 99(16):10411–10416, 2002. → pages 43[62] Philipp Hofemeier, Lihi Shachar-Berman, Janna Tenenbaum-Katan,Marcel Filoche, and Josue´ Sznitman. Unsteady diffusional screeningin 3d pulmonary acinar structures: from infancy to adulthood. Journalof biomechanics, 2015. → pages 43[63] J Sznitman. Convective gas transport in the acinus: revisiting the roleof effective diffusivity. In 6th World Congress of Biomechanics (WCB2010). August 1-6, 2010 Singapore, pages 370–373. Springer, 2010. →pages 44[64] CM Rhie and WL Chow. Numerical study of the turbulent flow past anairfoil with trailing edge separation. AIAA Journal, 21(11):1525–1532,1983. → pages 48[65] S Patankar. Numerical Heat Transfer and Fluid Flow. CRC press, 1980.→ pages 48[66] PD Wagner, DR Dantzker, R Dueck, JL Clausen, and JB West.Ventilation-Perfusion Inequality in Chronic Obstructive Pulmonary100BIBLIOGRAPHY BibliographyDisease. The Journal of Clinical Investigation, 59:203–216, 1977. →pages 56[67] Jessica M. Oakes, Alison L. Marsden, Cline Grandmont, Chantal Dar-quenne, and Irene E. Vignon-Clementel. Distribution of aerosolizedparticles in healthy and emphysematous rat lungs: Comparison be-tween experimental and numerical studies. Journal of Biomechanics,48(6):1147 – 1157, 2015. → pages 84101Appendix102Fluent Mesh Motion UDFThe mesh motion that simulates the rhythmic movement of pulmonaryairways is implemented in Fluent by the use of a User Defined Function(UDF). The function is written in a heavily templatted version of C++ andimplements Eqtn. 4.4 over all points in the computational domain.#include "udf.h"#include "unsteady.h"#include "mem.h"FILE *fout;/* &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& */DEFINE_ON_DEMAND(save_original_grid_to_nodes){#if !RP_HOSTsave_original_grid_to_nodes_func();Message0("\n\n Done! \n You can check this through plotting thecontours of user defined node memory (e.g. print onwalls)\n");#endif /*!RP_HOST */}void save_original_grid_to_nodes_func(){Domain *domain;cell_t c;Thread *t;Node *v;int n;domain=Get_Domain(1);/*Store the mesh node coordinates in user-defined node memory,this data is used in the dynamic mesh*/thread_loop_c (t,domain){103begin_c_loop (c,t){c_node_loop (c,t,n){v = C_NODE(c,t,n);N_UDMI(v,0) = NODE_X(v);N_UDMI(v,1) = NODE_Y(v);N_UDMI(v,2) = NODE_Z(v);}}end_c_loop (c,t)}}DEFINE_GRID_MOTION(Parenchyma1,domain,dt,time,dtime){#if !RP_HOSTThread *tf= DT_THREAD(dt);face_t f;Node *v;real NV_VEC(A);real NV_VEC(dx);real previous_time;int n;int v1;int i;float beta;/* calculate displacement vector at the node normal to the surface*/i=0;begin_f_loop(f,tf){previous_time=PREVIOUS_TIME;F_AREA(A,f,tf);beta = 0.0416;amplitude=beta/2;dx[0]=amplitude*(A[0]/NV_MAG(A)*waveform);dx[1]=amplitude*(A[1]/NV_MAG(A)*waveform);dx[2]=amplitude*(A[2]/NV_MAG(A)*waveform);f_node_loop(f,tf,n){v = F_NODE(f,tf,n);/* update node if the current node has not been previously *//* visited when looping through previous faces */104if ( NODE_POS_NEED_UPDATE (v)){/* indicate that node position has been update *//*so that it’s not updated more than once */NODE_POS_UPDATED(v);NODE_COORD(v)[0]=N_UDMI(v,0)*(1.0+ amplitude +amplitude*sin(1.62885*CURRENT_TIME - (3.141/2)));NODE_COORD(v)[1]=N_UDMI(v,1)*(1.0+ amplitude+amplitude*sin(1.62885*CURRENT_TIME - (3.141/2)));NODE_COORD(v)[2]=N_UDMI(v,2)*(1.0+ amplitude +amplitude*sin(1.62885*CURRENT_TIME - (3.141/2)));i=i+1;}}}end_f_loop(f,tf);#endif /*!RP_HOST*/}/* &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& */105


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items