§§ WMINIWIWWIWIHHI'HHM’WI'IIJHI

 

 

 

THESIS ‘

200%

This is to certify that the
thesis entitled

NUMERICAL SIMULATION OF AIRFLOW IN THE EQUINE UPPER
AIRWAY TO STUDY THE PHENOMENON OF DORSAL
DISPLACEMENT OF SOFT PALATE

 

presented by

f.-

SRIKANTH SRIDHAR

.LIBRARY
Michigan State
Universrty

has been accepted towards fulfillment
of the requirements for the

 

 

 

MS. degree in MECHANICAL ENGINEERING

 

 

W

Major Professor’s Signatdﬁ
2/015’ /M
I I v

Date

MSU is an afﬁrmative-action, equal-opportunity employer

—--.-L- ‘l-l—‘-I-n-l-“J-n—‘*

*-'4~‘-'-l‘&w-‘ -

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/07 p:/CIRC/DaleDue.indd-p.1

NUMERICAL SIMULATION OF AIRFLOW IN THE EQUINE UPPER AIRWAY To
STUDY THE PHENOMENON OF DORSAL DISPLACEMENT OF SOFT PALATE
By

Srikanth Sridhar

A THESIS

Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of

MASTER OF SCIENCE
Department of Mechanical Engineering

2008

ABSTRACT

NUMERICAL SIMULATION OF AIRF LOW IN THE EQUINE UPPER AIRWAY TO
STUDY THE PHENOMENON OF DORSAL DISPLACEMENT OF SOFT PALATE

By
Srikanth Sridhar
Dorsal Displacement of Soft Palate (DDSP) is a commonly observed upper—airway
respiratory disorder in horses, especially while racing. Either solely or in combination
with other conditions, DDSP is one of the most common causes for performance
retardation and abnormal respiratory noise in equine athletes. However, the exact causes
for the occurrence of DDSP are not yet understood clearly. Besides, little is known about
the dynamic pressure and velocity distributions that may contribute to dynamic collapse

of the soft palate into the upper airway.

In this work, the airﬂow through the equine upper airway was studied using
Computational Fluid Dynamics in order to investigate if there is correlation between a
particular ﬂow regime and the occurrence of DDSP. Post-mortem CT scan images of the
upper airway were used to develop a computational domain for the upper airway
geometry. Numerical computations were performed in order to solve the ﬂow ﬁeld for
different geometries and breathing conditions. The results of the numerical solutions were
validated by comparing with experimental results recorded in literature. Pressure and
velocity distributions were studied throughout the domain and throughout the breathing
cycle in order to ﬁnd out if there is a correlation between the ﬂow ﬁeld and the

occurrence of DDSP. The various assumptions applied to the modeling were discussed.

ACKNOWLEDGMENTS

I take this opportunity to thank my thesis advisor Dr. Mei Zhuang for her guidance with
this work and more importantly, for her patience, kindness and generosity of spirit. I
appreciate her knowledge and skill in the ﬁeld of ﬂuid mechanics. She has been a
constant source of encouragement and support, even at those times when I had fallen

short of her expectations.

I’d also like to thank Dr. N. Edward Robinson (College of Veterinary Medicine, MSU).
In addition to providing partial funding for this work and co-guiding this thesis, he has
also generously shared with me, gems from his vast treasure of knowledge and
experience in the a wide spectrum of ﬁelds ranging from biology to music! I hold great

appreciation for his fervor for knowledge and passion for science and science-teaching.

I thank Dr. Frederik Derksen (College of Veterinary Medicine, MSU) for his guidance
with this work and for not getting annoyed with my endless questions to him on equine

respiratory physiology, but answering them with patience.
The airway geometries used in this work were from Dr. Rob Malinowski (College of

Veterinary Medicine, MSU). I thank him for spending his time in developing the airway

models and sharing his valuable insights. I greatly admire his skill with computers!

iii

The assistance provided by Ms. Diana Rosenstein and Ms. Maggie Hofmann (College of

Veterinary Medicine) in developing the airway models is greatly appreciated.

I’d like thank Aida Montalvo (Graduate Secretary, Department of Mechanical
Engineering) for all her help throughout my stay at MSU. Her secretarial skills,

experience and communication ability have saved me lot of time, effort and energy.

Dr. Farhad Jaberi (Department of Mechanical Engineering, MSU), and Dr. Charles Petty
(Department of Chemical Engineering, MSU) are greatly appreciated for sharing their
valuable insights on turbulence modeling. Dr. Jaberi was one of the members of the thesis

committee.

I thank my family, friends and my lab-mates at the CFD lab, MSU for their moral
support. In particular, I’m very grateful to my ﬁ'iends Deebika and Saravan for having
helped me complete my thesis in several ways, including proof-reading and offering their

comments.

1v

TABLE OF CONTENTS

LIST OF TABLES ............................................................................................................. vi
LIST OF FIGURES .......................................................................................................... vii
Chapter 1: Introduction ....................................................................................................... 1
1.1. Dorsal Displacement of Soft Palate ......................................................................... 1
1.2. CFD and its application to modeling respiratory ﬂows ........................................... 3
1.3. Scope of the work .................................................................................................... 3
1.4. Organization of the thesis ........................................................................................ 4
Chapter 2: Problem Description .......................................................................................... 5
2.1. Description of DDSP ............................................................................................... 5
2.2. Nature of flow within the upper airway ................................................................... 6
2.3. Problem statement .................................................................................................... 7
Chapter 3: Method .............................................................................................................. 8
3.1. The upper airway model .......................................................................................... 8
3.1.1. Computational domain corresponding to the upper airway .............................. 8
3.1.2. Manipulating the upper airway geometry ....................................................... 10
3.1.3. Pre-processing and Meshing ........................................................................... 13
3.2. Mathematical formulation ...................................................................................... 17
3.2.1. Governing equations ....................................................................................... 17
3.2.2. Boundary conditions ....................................................................................... 19
3.2.3. Initial condition ............................................................................................... 25
3.3. Numerical solution methods .................................................................................. 29
3.4. Choice of turbulence model ................................................................................... 30
3.5. Assumptions in modeling ...................................................................................... 33
Chapter 4: Results and Discussion .................................................................................... 35
4.1. Validation ............................................................................................................... 35
4.2. Velocity and pressure distributions: ....................................................................... 37
4.2.1 Velocity and Pressure distributions through the tracheal pressure cycle ......... 37
4.2.1. Variation in velocity distribution with respiratory rate ................................... 50
4.2.2. Variation of velocity distribution with ﬂexure in upper airway ..................... 52
4.3. Wall pressure along the dorsal surface of lower wall ............................................ 55
4.3.1. Variation of wall pressure distribution with respiratory rate .......................... 57
4.3.2. Variation of wall pressure distribution with induced ﬂexure in the upper
airway ........................................................................................................................ 59
Chapter 5: Conclusions and Future Work ......................................................................... 61
References ......................................................................................................................... 65

LIST OF TABLES

Table 3.1: Unsteady pressure boundary condition corresponding to different rates of
respiration ....................................................................................................... 24

Table 4.1: Comparison between experimental and numerical results ............................. 35

vi

LIST OF FIGURES

IMAGES IN THIS THESIS ARE PRESENTED IN COLOR

Figure 2.1. Sagittal sectional view of equine upper airway ................................................ 5
Figure 3.1. Image from CT scan of upper airway ............................................................... 9
Figure 3.2. Two-dimensional model developed from CT scan ........................................ 10
Figure 3.3. Model representing airway with no induced ﬂexure ...................................... 11
Figure 3.4. Model representing airway with mild induced ﬂexure ................................... 12
Figure 3.5. Model representing airway with extreme induced ﬂexure ............................. 12
Figure 3.6. Boundaries and zones during inhalation ......................................................... 13
Figure 3.7. Boundaries and zones during exhalation ........................................................ 14
Figure 3.8. Velocity gradient based adaptive reﬁnement of mesh ................................... 16
Figure 3.9. Generic boundary conditions for respiratory ﬂow ......................................... 20
Figure 3.10. Tracheal pressure cycle during rest .............................................................. 21
Figure 3.11. Tracheal pressure cycle during mild exercise .............................................. 22
Figure 3.12. Tracheal pressure cycle during intense exercise .......................................... 23
Figure 3.13. Locations of test points within the domain for monitoring solutions ........... 26

Figure 3.14. Variation of pressure with respect to time at the test point ‘a’ indicated in
Figure 3.13 ................................................................................................... 27

Figure 3.15. Variation of x-velocity with respect to time at the test point ‘b’ indicated in
Figure 3.13 ................................................................................................ 28

Figure 3.16. Locations of test surfaces within the computational domain for comparing
solutions obtained using different models ................................................... 31

vii

Figure 3.17. Comparison between k-epsilon, k-omega and RSM model based solutions —

velocity proﬁle at cross-section B1-B2 ........................................................ 32
Figure 3.18. Orientation of the coordinate axes within the domain .................................. 33
Figure 4.1. Velocity distribution (m/s) at peak inspiration ............................................... 38
Figure 4.2. Occurrence of peak inspiration relative to the tracheal pressure cycle .......... 38
Figure 4.3. Pressure contours (Pa) at peak inspiration ...................................................... 40
Figure 4.4. Occurrence of peak inspiration relative to the tracheal pressure cycle .......... 40

Figure 4.5. Velocity distribution (m/s) at end of inspiration and just before expiration .. 42

Figure 4.6. Occurrence of end of inspiration relative to the tracheal pressure cycle ........ 42
Figure 4.7. Pressure contours (Pa) at end of inspiration and just before expiration ......... 43
Figure 4.8. Occurrence of end of inspiration relative to the tracheal pressure cycle ........ 43
Figure 4.9. Velocity distribution (m/s) at peak expiration ................................................ 45
Figure 4.10. Occurrence of peak expiration relative to the tracheal pressure cycle ......... 45
Figure 4.11. Pressure contours (Pa) at peak expiration .................................................... 46
Figure 4.12. Occurrence of peak expiration relative to the tracheal pressure cycle ......... 46

Figure 4.13. Velocity distribution (m/s) at end of expiration and just before inspiration 48
Figure 4.14. Occurrence of end of expiration relative to the tracheal pressure cycle ....... 48
Figure 4.15. Pressure contours (Pa) at end of expiration and just before inspiration ....... 49
Figure 4.16. Occurrence of end of expiration relative to the tracheal pressure cycle ....... 49
Figure 4.17. Velocity distribution (m/S) at end of inspiration — rest condition ................. 51
Figure 4.18. Velocity distribution (m/s) at end of inspiration - mild exercise condition . 51

Figure 4.19. Velocity distribution (m/s) at end of inspiration - intense exercise condition
...................................................................................................................... 52

Figure 4.20. Velocity distribution (m/s) at end of inspiration — no induced ﬂexure ........ 53

viii

Figure 4.21. Velocity distribution (m/s) at end of inspiration — mild induced ﬂexure ..... 54
Figure 4.22. Velocity distribution (m/s) at end of inspiration — extreme induced flexure 54

Figure 4.23. Wall pressure along dorsal surface - variation through the tracheal pressure

cycle ............................................................................................................. 56
Figure 4.24. Wall pressure along dorsal surface — variation with rate of respiration ....... 58
Figure 4.25. Wall pressure along dorsal surface — variation with induced ﬂexure .......... 59

ix

CHAPTER 1: INTRODUCTION

In this work, Computational Fluid Dynamics (CF D) was applied to predict the unsteady
pressure and velocity distributions within the equine upper airway, corresponding to a set
of different respiratory conditions. The motivation behind this work was to determine in
what way the pattern of airﬂow established within the upper airway inﬂuences the
development of Dorsal Displacement of Soft Palate (DDSP) which is a common

performance-retarding respiratory disorder in equine athletes.

This chapter ﬁrst introduces the reader to the condition of Dorsal Displacement of Soft
Palate in horses. This is followed by a very brief introduction to Computational Fluid
Dynamics and a crisp overview, with examples, of the application of CPD to modeling

biological ﬂows. Finally, the organization of this thesis is concisely described.

1.1. Dorsal Displacement of Soft Palate

Dorsal Displacement of the Soft Palate which is commonly referred as ‘DDSP’ is a
common upper airway respiratory disorder that manifests in horses, most often, during
intense exercise [1]. In a horse exhibiting this condition the soﬁ palate within the upper
airway whose caudal end (or the end in the direction toward the animal’s tail), which is
usually held in position by the epiglottis muscles is dynamically displaced dorsally
(toward the animal’s back) into the upper airway, typically during intense exercise.
DDSP results in a dynamic obstruction to airﬂow, particularly during expiration [2], and

1

is often associated with exercise intolerance and/or abnormal respiratory noise [3].
Intermittent DDSP, either solely or in combination with other disorders is one of the most
common causes for poor performance and exercise intolerance in equine athletes [4].
Hence it is important to study the problem and investigate its causes. However, the fact
that this condition, quite often, manifests itself only during exercise while without a
clinical sign while resting, poses a special challenge for diagnosis and treatment. It is
evident from literature that multiple, often interrelated factors, may contribute to DDSP

including dysfunction of the epiglottis [5], Skeletal muscles [6], and nerves [7].

A few researchers had hypothesized that the ﬂuid mechanics in the equine upper airway
may be contributing to DDSP. For instance, Robinson et al. [8] suggested negative
inspiratory pressures in the pharynx may be a cause for DDSP. However, Rehder et al.
[9] measured dynamic upper airway pressures and correlated the data with results from
endoscopic examination in exercising horses. Based on this, they disagreed that there is a
link between the magnitude of negative inspiratory pharyngeal pressures and occurrence
of DDSP. Currently, little is known about the distributions of unsteady pressures and
velocities inside the upper airway which might be contributing to this condition. Hence it
is worthwhile to study the mechanics of ﬂow behind this problem in order to investigate

if there is a link.

1.2. CFD and its application to modeling respiratory flows

Computational Fluid Dynamics (CF D) is a technique wherein the complex partial
differential equations governing the parameters of a ﬂow ﬁeld are approximated to sets of
algebraic equations that are solved using specialized numerical techniques to solve for the
unknown ﬂow parameters. There have been several instances where researchers have
used CFD to successfully predict the unknown respiratory ﬂow parameters. For instance,
Kepler et al. [10] successfully predicted and validated with experimental results, the nasal
airﬂow and gas uptake patterns in rhesus monkey by using CF D simulations, while
Kimbell et al. [11], predicted ﬂow-dependent nasal uptake patterns of formaldehyde in
rat, rhesus monkey and human nasal passages. CFD simulations can not only facilitate
analysis by experimentation, but can also be a vital tool in situations where
experimentation is not feasible or when the accuracy of experimental results is not within
acceptable limits. Moreover, with a computer model to represent the physical problem,
changes in the analysis can be easily incorporated by altering the corresponding features

of the computer model, obviating lengthy experimental set ups.

1.3. Scope of the work

The objective of this work was, to solve the ﬂow ﬁeld inside the equine upper airway in
order to study the distribution of ﬂow parameters such as velocity and pressure on order

to investigate whether there is a possible correlation to the occurrence of DDSP. In

particular, it was desired to study the effect of the following on the distribution of ﬂow
parameters:

1. Intensity/ Rate of respiration

2. Flexure in the airway (which is frequently linked with increase in the respiratory

rate)

1.4. Organization of the thesis

A detailed description of the problem is provided in Chapter 2 along with a description of
the nature of the ﬂow within the upper airway. Chapter 3 elaborates the methods
followed in this work, beginning with creating the upper airway model. The mathematical
formulation, the governing equations solved, turbulence model applied and the numerical
solution techniques used for solving them are presented. The assumptions applied to the
approach are also described. In chapter 4, the results obtained and their validation details

are discussed in detail. Finally, conclusions based on our results are presented.

CHAPTER 2: PROBLEM DESCRIPTION

This chapter describes the problem of this thesis in detail.

2.1. Description of DDSP

 

’/

Ill/W

mu... \ - ’\’////// .
//"“"'\‘ ‘ //////
ﬂf////////%y/////Imml,,,....n.OI:/.//////,

7 I 0
ll// Hard

Palate
Palate Trachea

    
 
 

 
 
  

 

  

I

Epiglottis Pharynx 0.

 

 

 

 

Figure 2.1. Sagittal sectional view of equine upper airway

Figure 2.1 depicts upper airway comprising of the nasal cavity, hard palate, soft palate,
epiglottis and the pharynx. A portion of the lower airway (trachea) is also shown in the
Figure. The locations of the hard and the soft palates, the epiglottis and the trachea are to

be noted. The soft palate, as the name suggests is pliant and, except during swallowing is

kept beneath the epiglottis by muscular activity. Unlike humans, horses cannot breathe
through their mouth. The epiglottis acts as a control valve opening the rear of the mouth
either solely into the trachea (as in the Figure) while breathing or solely into the
esophagus (while swallowing). When DDSP occurs during breathing, the soft palate
displaces above the epiglottis and into the airway, obstructing expiratory air coming from
the trachea, allowing a portion of the expiratory air into the mouth, thereby vibrating the

palate resulting in a loud expiratory noise [12].

2.2. Nature of ﬂow within the upper airway

The airﬂow is driven by the pressure difference between atmospheric pressure at the
nares and the tracheal pressure. The latter is governed by the oscillating pressure induced
by the diaphragm which expands and contracts alternately, causing positive tracheal
pressure for one half-cycle (thereby causing exhalation) and negative tracheal pressure

for the next half-cycle (thereby causing inhalation).

Experiments recorded in literature [13] indicate that the ﬂow, during intense exercise,
ranges from transition-to-turbulent to fully turbulent. The ﬂow is three dimensional and
unsteady since it’s driven by an unsteady pressure difference. Since air is a highly
compressible ﬂuid, and since the ﬂow is achieved by forcing oscillatory pressures, the
ﬂow is expected to exhibit some compressibility effects even though the velocities are

expected to be far lesser when compared to that of sound.

The frequency of these time-periodic ﬂows is governed by the forcing pressure’s
frequency. Low frequencies and amplitudes characterize resting conditions and higher

ﬁequencies and amplitudes characterize exercising conditions.

Due to the pulling of reins on the horse, typically while racing, there is expected to be an
artiﬁcially induced ﬂexure in the upper airway near the pharynx. This additional ﬂexure
at the pharyngeal cavity could have a bearing on the ﬂow parameters within the whole of

the domain.

2.3. Problem statement

The statement of the problem was:
1. To solve the ﬂow ﬁeld within the domain of the upper airway for various
respiratory rates and ﬂexure angles.
2. To study and understand the distribution of the ﬂow parameters for these different
cases.
3. To investigate if the results obtained could indicate that the mechanics of the ﬂow

contributes to Dorsal Displacement of Soft Palate.

CHAPTER 3: METHOD

In this chapter, the steps involved in the work, beginning with obtaining the geometric
model representing the upper airway to solving the ﬂow within the domain, applying the
appropriate boundary conditions are described in detail. A detailed description of the
mathematical modeling applied to this problem, the numerical techniques used and the

simplifying assumptions incorporated are also included.

3.1. The upper airway model

The ﬁrst part of this section describes the development of the computational domain
representing the upper airway. The next part elaborates the manipulation of the geometry
to obtain models corresponding to different degrees of induced ﬂexures. A description of

the schemes used for meshing these domains is described in the last section.

3.1.1. Computational domain corresponding to the upper airway

In order to obtain a representative geometry for the upper airway, a post-mortem CT scan
was performed on a 15 year old, mixed breed horse. The horse was euthanized for other
reasons. The CT scan generated images of the transverse section of the horse’s upper
airway. These images were stored in the ‘dicom’ format, which one of the standard
formats for medical images from CT. The images were imported into the MPR (Multi-

Planar Reconstruction) module of the program called ANALYZE 6.0 (Mayo Clinic). This
8

program generated two-dimensional images of the sagittal and rostral sections ﬁ'om the
transverse section images generated by CT scan. In Figure 3.1 one such CT scan image

which is used to obtain the sagittal section of the upper airway is shown.

 

    
  

INostril

 

   

 

 

 

Figure 3.1. Image from CT scan of upper airway

The sagittal section images were saved in ‘tiff’ format and exported to Adobe Illustrator.
This program traced the sagittal section images and generated the vector format,

consisting of points and lines and saved the data in ‘.ai’ format.

Finally, the program 3DS Max (AutoDesk) converted the vector format data to IGES
format and scaled to appropriate dimensions. This data can now be imported into
GAMBIT (Fluent Inc.) for meshing and other pre-processing activities. Figure 3.2
shows the image of sagittal section of the upper airway reconstructed from the CT scan
data using the procedure mentioned above. The location of the different anatomical

features — nares, nasal cavity, turbinates, pharynx and larynx are indicated.

Nasal Turbinates

 

Figure 3.2. Two-dimensional model developed from CT scan

3.1.2. Manipulating the upper airway geometry

As mentioned earlier, one of the objectives of this study was to study the effects of the
artiﬁcial ﬂexure induced in the upper airway. In order to solve the ﬂow ﬁeld
corresponding to these conditions, it is needed to develop the corresponding
computational domains. In order to achieve this, the ‘3DS Max’ program was used to

manipulate the geometry to create models representing different degrees of ﬂexure in the

10

upper airway, keeping in mind to not alter the relative positions of the anatomical
features. Hence from the original upper airway 2D model described in Figure 3.2 and
Figure 3.3, two additional models were thus created representing slightly ﬂexed and
extremely ﬂexed airways. The models were ﬂexed around the pharynx as shown in
Figure 3.4 and Figure 3.5. These geometries represent the ﬂexing of the upper airway
which is induced by pulling of reins on a horse, while racing. The ﬂow ﬁeld would be
solved later on, within each of these domains, corresponding to identical sets of boundary

conditions.

/L

No Artiﬁcial Flexure

Figure 3.3. Model representing airway with no induced flexure

ll

Mild Artiﬁcial Flexure Induced

Figure 3.4. Model representing airway with mild induced ﬂexure

 

 

 

/

Extreme Artiﬁcial Flexure Induced

Figure 3.5. Model representing airway with extreme induced flexure

12

3.1.3. Pre—Processing and Meshing

The pre-processing activities that are needed to be performed in order to prepare for the
solving the ﬂow ﬁeld, are described in this section. Firstly, the IGES data ﬁles described
in the previous section were imported into GAMBIT 2.3.16 (Fluent Inc.). The edges in
the two-dimensional model corresponding to different zones, viz. inlet, outlet, walls were
deﬁned. This is needed to be done so that the appropriate boundary conditions can be
applied to the different zones. The upper and lower boundaries of the geometry were set
as no slip wall zones. The left (nostril) and the right (tracheal opening) boundaries were
deﬁned as inlet/outlet zones. During inhalation, the nostril boundary would be an inlet

and the tracheal opening boundary would be an outlet as shown in Figure 3.6.

Upper wall

   

  

Lower wall

Nostn'l \

(Inlet) Tracheal
Opening
(Outlet)

Figure 3.6. Boundaries and zones during inhalation

13

And during exhalation, the tracheal opening boundary would be an inlet and the nostril
boundary would be an outlet as shown in Figure 3.7 . The pressure difference between the

nostril boundary and the tracheal opening boundary would govern the ﬂow direction.

Upper wall

/%Lower:a\
Nostril \,

(Outlet) Tracheal
Operung
(Inlet)

Figure 3.7. Boundaries and zones during exhalation

Once the different zones are deﬁned it’s needed to divide the domain into a ﬁnite number
of tiny cells. This operation, known as meshing and is needed to be done so that the
partial differential equations governing the ﬂow of the ﬂuid continuum can be
approximated into sets of algebraic equations solved at the discrete nodes introduced by
the meshing. The physical domain’s extents are of the order 70 cm x 20 cm. (Since the
anatomical structures do not ﬁt regular geometric shapes, it is not possible to specify the
exact extents in each coordinate direction). The area of the domain corresponding to the
non ﬂexed model was about 4.5 x 10'2 m2. The domain was meshed with an unstructured
meshing scheme by paving quadrilateral elements. Due to the irregular shape of the
anatomical structures, structured meshing was not appropriate for this situation.

14

The size of the quadrilateral elements chosen determined the mesh density. An initial
mesh size of 0.005 m was selected. The solution obtained for this mesh density was
recorded. (A steady-state problem was solved for this purpose, with inlet pressure at
atmospheric pressure and the outlet pressure at the maximal diaphragm pressure during
extreme exercise.) Then, the mesh was progressively reﬁned and the solutions obtained
thereby were compared with the earlier solution(s) in order to check for grid
independence. At a mesh size of 0.0005 m, the solution was more or less independent of
the grid. All the models were meshed with this same density. And for this mesh density,

the number of quadrilateral elements throughout the domain was of the order of 100,000.

In order to adequately resolve the regions with high velocity gradients, a steady state
solution was obtained, corresponding to the maximal inspiratory ﬂow rate during intense
exercise and grid adaptation, based on velocity gradients was at this particular solution. It
was assumed that the maximal velocity gradients would exist in during peak inspiration
and hence would require the maximal spatial resolution of grid, representing a “worst-
case scenario”. The mesh was progressively reﬁned in the same way until a gradual
change in mesh density from the walls to the center of the streamline was observed.
Figure 3.8 illustrates this process of adaptive grid reﬁnement that resulted in a ﬁner mesh
density near the walls and other regions of high velocity gradients than those of gradual
variations in velocity. In addition, after each reﬁnement, the mesh quality was examined
by checking the variation in diagonal ratio between adjacent elements throughout the

domain and was thus reﬁned.

15

 

Regions with high velocity gradients resolved ﬁnely

O
O . v.‘
U Q Q . O 0
. .9009 ’9’.
5‘0“.“ ...:..O:¢:O:9:O: I
man-a. ' w: "-
' ”-5.33.” "
I“...:.‘ O
l:.'. n

.::. ...

o. O .‘
swam
6.0. .o“'l ‘“
. . . «ma.»-

0 0.0... . n

so. 0..

00. O I

«swat
‘.|

 

9 e
. .o‘o’o’o’o 0.. 0:0 0 o o .. b.
$9....on o’o o’o’o‘v a...
"3:3:‘9‘9’0’9 0’0?0'0’0’6’:°:0’0
> “1'90‘9’9’9’Po’o 0‘0”.” ' 0'
... .QQQQOOOOQQJOO
i"='9o’em§::’:6,o“"0:o“ o 6609 .eeo‘
ﬁzz-9“ Acme... 3.990. 0.049;

 

Figure 3.8. Velocity gradient based adaptive reﬁnement of mesh

The mesh thus obtained was saved and exported to the F LUENT program which was

used to solve the governing equations for the problem.

3.2. Mathematical formulation

3.2.1. Governing equations

The governing equations solved for this problem are the Reynolds Averaged Navier
Stokes Equations (RANS) in two dimensions. The generalized RANS equations are

given as follows [14, 15]:

Continuity equation:

500) . = (3.1)
at +V (<pV>) 0

 

Momentum equation:

 

 

6<v> , __ . _ v '
6t +V (<pvv>)— V<p>+V (<;> (pv v >> (3.2)
Energy equation:
8C T
(p >+V-((pC1,,Tv>)=9333+—1—(Irw>-V<p>+<v'-Vp>+<<1>>
at at (p)
(3.3)

Where

= 9:. (3'4)
<¢> (Tuax->

17

The k-epsilon model, which belongs to the broader class of Two-equation models, is used
to introduce to achieve the closure for the ( p v ' v ') term in the momentum equation,

which is also known as the Reynolds Stress term..

According to this method, in addition to the RANS equations, two other equations are
solved. Each of these equations conserve one additional scalar introduced in order to

close the Reynolds stress term. The additional scalar quantities that are introduced are

turbulent kinetic energy (k) and Turbulent Kinetic Energy Dissipation Rate (8). The

following are the two additional conservation equations introduced because of the k-s

turbulence model [15, 16]:

k-equation:

,u
Lé[)tl(l+V-(pkv)=V- [(,u+-O%)Vk]+ Gk+Gb—pa—YM +Sk

(3.5)
s-equation:
2

6(p8) 1“: ‘ a a
——+V- v =V- +— V + — G +C G —C —+S

(3.6)
where,
C” = 0.09
Cg1 =1.44

18

‘92
0k =1.0
05:13

Gk and Ch represent the changes in turbulent kinetic energy due to mean velocity
gradients and due to buoyancy respectively. YM represents the contribution of the
ﬂuctuating dilatation to the overall dissipation rate [16]. 3., and SE represent source terms.

K and epsilon are related by the following equation:

k2
#t = PC). —g- (3.7)

Finally, the closure is obtained as follows by applying the turbulent-viscosity hypothesis:

—pv.v. +—p ..= pp —+—
l j 3 l_] t ax}. 5x1. (3.8)

3.2.2. Boundary conditions

The ﬂow within the domain is driven by the pressure difference between the tracheal
opening and the nostril. The nostril is at atmospheric pressure and the pressure at the
tracheal opening is an unsteady pressure induced by the cyclic expanding and contracting
action of the lungs which in turn is achieved by the work done by the muscles of the
diaphragm. At the nostril, the pressure was set as atmospheric pressure at sea level and
room temperature. At the tracheal end of the upper airway, a sinusoidal oscillating
pressure was applied. Although the functional form of the unsteady tracheal pressure is

not known, the sinusoidal function with amplitude equal to the peak tracheal pressure and

19

frequency equal to breath frequency was assumed to be a suitable approximation. At all
the wall zones, the no-slip boundary condition was applied. A generic illustration of the

boundary conditions applied to this problem is shown in Figure 3.9.

bk)Shp

 
 

Atmospheric \

Pressure /“\

 

Presarre

 

 

 

Tune

Figure 3.9. Generic boundary conditions for respiratory flow

The frequency and the amplitude of the sinusoidal function applied at the tracheal
opening were based on the tracheal pressures and number of breaths per minute for
different resting and peak exercising conditions recorded vastly in literature. (e. g.,
Lumsden et al.). An intermediate state — ‘mild exercise’ was also introduced. An
unsteady pressure function whose amplitude and frequency are about midway between
rest and peak exercise was formulated, corresponding to this condition. The sinusoidal
unsteady pressure functions applied at the tracheal opening are illustrated in Figure 3.10

(rest), Figure 3.11 (mild exercise) and Figure 3.12 (intense exercise).

20

Tracheal Pressure (Pa)
- Rest

 

3000 r . . r . .

2000

T

a

I

-1000

-2000 r

 

 

'30000 0.5 1 1.5 2 2.5 3

Time (s)

Figure 3.10. Tracheal pressure cycle during rest

21

3.5

 

Tracheal Pressure (Pa)
- Mild Exercise

'5

 

3000 I I I I I I

§

3

i

 

 

 

Time (s)

Figure 3.11. Tracheal pressure cycle during mild exercise

22

3.5

 

3000

 

2000 ~ _

 

 

1000

T

Tracheal Pressure (Pa)
-Intense Exercise
O

4000* a

-2000 - *

-3...U .U ,U ,U .U .U U I

0 0.5 1 1.5 2 2.5 3.5 4
Time (s)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

mr

Figure 3.12. Tracheal pressure cycle during intense exercise

As a summary, Table 3.1 lists the amplitude and the frequency of these time-periodic

functions and the breathing rates for the different intensities of exercise.

23

Table 3.1. Unsteady pressure boundary condition corresponding to different rates of

 

 

 

 

respiration
Exercise Maximum Tracheal Breaths per min. Frequency
S. No. Condition Pressure ‘ I Pm, I ’ (Pa) ‘1" (Hz)
1. Rest 490 15 0.25
2. Mild 1500 60 1
Exercise
3. Intense 2940 1 20 2
Exercise

 

 

 

 

 

 

 

The tracheal pressure boundary condition was incorporated in such a way that the
inhalation half-cycle began ﬁrst and was followed by the exhalation half-cycle. Hg. for

peak exercise, the function applied was:

p = — 2940 sin (4m)

In order to solve the ﬂow corresponding to the each exercise condition (rest, mild
exercise and intense exercise), the corresponding unsteady pressure function was applied.
As mentioned earlier, one of the objectives of this study was to study the effects of the
artiﬁcial ﬂexure introduced into the upper airway on the ﬂow ﬁeld. Hence, a set of nine
numerical computations were performed, each corresponding to a particular degree of
artiﬁcially induced ﬂexure (none, mild, extreme) and a particular rate of breathing (rest,

mild exercise and intense exercise).

24

The boundary conditions at the top and the bottom walls were standard no-slip wall
boundary condition, wherein the ﬂuid velocity at the wall is equal to the velocity of the

wall, which in this case equals zero.

Since the energy equation was also solved, it was required to supply appropriate
temperature boundary conditions. At the nostril, a standard room temperature of 25
degree centigrade was imposed. At the trachea, temperatures of 30, 35 and 40 degree
centigrade were imposed for rest, mild exercise and intense exercise conditions
respectively. In order to solve the equations describing the turbulent quantities, it was
needed to specify the intensity or turbulence at the inlet and the outlet. This was speciﬁed

as 1% for resting, 5% for mild exercise and 10% for intense exercise.

3.2.3. Initial condition

Since the ﬂow was initialized to almost no ﬂow throughout the domain (which is an
unphysical condition), the solution was monitored for a large number of time steps,
covering several time-cycles in order to make sure that the solution attained a time-
periodic state. In order to ensure that the time-periodic solution has been achieved,
pressures and velocities were monitored at randomly chosen points in the domain and the
unsteady state solution was calculated out until time-periodic state was established at all
of these points. Once this state was attained, this converged time-periodic solution was
saved and that solution was used as a starting point for further simulations. Figure 3.13

illustrates the location of two points ‘a’ and ‘b ‘that were chosen within the domain.

25

Figure 3.14 shows the variation of pressure with respect to time at point ‘a’. And Figure
3.15 shows the variation of the x component of velocity with respect to time at point ‘b’.
It is seen that these quantities exhibit a time periodic behavior alter a few complete

cycles.

 

Figure 3.13. Locations of test points within the domain for monitoring solutions

26

 

300

3

mm

 

100~

 

...............

U U U W

I I I l

L l l
0 20 40 60 80 100 120 I40 160 180 200
Time Step

 

 

 

 

 

 

 

 

 

 

 

Pressure (Pa) at Point ’a'

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 3.14. Variation of pressure with respect to time at the test point ‘a’ indicated
in Figure 3.13

27

x-velocity(m/s) at Point 'b’

10 I I I I I I I I I

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

-10 I 4 I I I I I

J I
0 20 40 60 80 100 I20 I40 160 180
Time Step

Figure 3.15. Variation of x-velocity with respect to time at the test point ‘b’
indicated in Figure 3.13

28

200

3.3. Numerical solution methods

A commercial ﬁnite volume based CFD solver (F LUENT 6.2.16, 2 D version) was used
to perform the computations. Unsteady, segregated, ﬁrst-order implicit solver with time
steps of 0.025 sec (for exercise), 0.125 (for mild exercise) 0.2 sec (for rest) were used.
These time steps were found to be suitable in achieving convergence in a reasonable
number of iterations. An unsteady pressure was applied at the outlet, as described earlier.
K—epsilon turbulence model was applied and the turbulent quantities were speciﬁed by
setting the turbulent intensity (1% for rest condition, 5% for mild exercise condition and
10% for intense exercise condition) and hydraulic diameter (0.05 m at the inlet and outlet
for all models and respiratory conditions). Convergence of the solution was observed by
monitoring the residuals of the ﬁve conservation equations. Convergence criteria were set
to 1e-3 for the scaled residuals of continuity and momentum equations and 1e-6 for those
of k, e and energy equations. First order upwind scheme was used for the convective

terms.

For each case, in order to compute the solution over one complete breathing cycle, it took
about 3 hours on a Pentium 4, 2.4 GHz CPU. The solution was computed for 10 complete
breathing cycles, starting with zero velocity initial condition at all nodes, in order to
achieve a time-periodic state. The solution at the end of 10 complete breathing cycles was
used as an initial condition for computations from that point onwards. Once a time-
periodic state was achieved, the computations were performed for one complete breathing

cycle and the solutions at each time step was recorded in the form of a data ﬁle.

29

 

 

 

 

3.4. Choice of turbulence model

This section describes the reasoning behind the choice of the k-epsilon turbulence model

for this problem.

In order to close the Reynolds stress term introduced by the averaging operator applied in
the RANS equations, various models have been developed. The k-e and k-(I) are models
belong to the broader class of ‘two-equation models’, which are based on the turbulent
viscosity hypothesis. Yet another model, known as the Reynolds Stress Model has been
developed in order to provide closure to the RANS equations. This model is widely
regarded as very appropriate when resolving the boundary layers, near wall regions are
important and when the mainstream velocity exhibits appreciable variation over
macroscopic time-scales. For this reason, Reynolds Stress Model would be the most
appropriate choice for this problem. However, the computational effort required is
signiﬁcantly large in this case because of the additional equations introduced, each
corresponding to a particular component of the Reynolds stress (presenting these
equations here is beyond the scope of this work). For the sake of a benchmark
comparison, the ﬂow was solved for mild exercising conditions for the mild induced
ﬂexure geometry using RSM, k-epsilon and k-omega models. The solutions were
compared by observing the velocity proﬁles established at different cross-sections
throughout the domain. Figure 3.16 shows two such cross sections — Al-A2 and Bl-B2
that were chosen to compare the velocity proﬁles established using each of the turbulence
models mentioned above. Figure 3.17 is a plot of the velocity proﬁles across cross

section BI-BZ.
30

From such comparisons, it was observed that there was no signiﬁcant difference in the
average velocity across the cross-section between the k-epsilon and the RSM models,
whereas the k-omega model differed slightly. Due to greater ease of use and greater
efﬁciency in achieving convergence faster, the k-epsilon model was used to model

turbulence in this problem.

Al Bl

Figure 3.16. Locations of test surfaces within the computational domain for
comparing solutions obtained using different models

31

 

I-‘N
0°C

 

 

 

—-
O‘

 

 

g—e
.5

.—|
0

Velocity Magnitude (m/s)
co 5

 

 

 

1 J

0 20 40 60 80 100
Percentage Curvelength along Cross-Section

Figure 3.17. Comparison between k-epsilon, k-omega and RSM model based
solutions - velocity proﬁle at cross-section 81-82

A similar comparison across section Al-A2 yielded similar inferences.

32

3.5. Assumptions in modeling

The following are the important assumptions of this work:

The model that was developed from a dead horse adequately represents all the
features of the living horse, especially with respect to collapsible tissues like
nostrils.

Variation of ﬂow parameters in the z direction (indicated in Figure 3.18) is

negligible.

L).

Figure 3.18. Orientation of the coordinate axes within the domain

The models corresponding to the different degrees of ﬂexure closely correspond
to the ﬂexure induced by pulling the reins.

The model is closely comparable to horses of different sub-species and genders.
The ﬂow is single phase. i.e., the effects of moisture during exhalation half-cycle
are not important.

The tracheal pressure is a time periodic function given by a sinusoidal function of

time.

33

The dynamic changes in geometry due to the horse’s gait and gallop are
negligible. I

The dynamic changes in the geometry due to collapse of tissues such as soft
palate, nostril etc. into the airway are negligible.

The air entering the horse’s nostril has zero velocity, i.e. the effect of air entering
the nostrils with a certain speed when the horse is racing or due to environmental
factors like wind etc. are not important.

Inspired and expire air have same properties. i.e., Changes in properties like
speciﬁc heat, density etc., are negligible.

The weight of the jockey riding the horse has negligible effect in the driving
pressures in the trachea.

The walls of the model are rigid. While in reality, the tissues that make up the
wall of the domain are pliant and change the cross sectional areas throughout the

breathing cycle.

34

 

CHAPTER 4: RESULTS AND DISCUSSION

Validation details, results and discussions are presented in this chapter.

4.1. Validation

Experimental (recorded in Lumsden et. al [16]) and numerically calculated values of ﬂow

rates at peak inspiration and at peak expiration were compared. The comparison is shown

in Table 4.1.

Table 4.1. Comparison between experimental and numerical results

 

Exercise condition

Peak inspiratory flow rate
(L/s) - (experimental)

Peak expiratory flow rate
(L/s) - (numerical)

 

 

 

 

 

Rest 3.50 5.92
Mild exercise 42.98 46.92
Peak exercise 78.21 81.00

 

The experimental values presented in the table are interpolated values from the actual
experimental values. The interpolation was done in order to match the ﬂow rates with the

exact breathing frequencies corresponding to which the numerical calculations were

performed.

35

 

From the unsteady state computation, the time steps corresponding to peak inspiration
and peak expiration were noted. At these time instants, ﬂow rates across the nostril inlet
boundary were calculated. The two-dimensional numerical solutions delivered the ﬂow
rates based on a unit channel width, which in this case equals one meter. These values
were scaled down to an equivalent channel width for the nostril opening corresponding to
rest (1 cm), mild exercise (3 cm) and intense exercise (6 cm) in order to account for the
expansion of the nostrils during exercise. Numerically calculated values for peak ﬂow

rates were found to be comparable with experimental values.

36

4.2. Velocity and pressure distributions:

It is worthwhile to study the velocity and pressure distribution throughout the domain at
important time instants in an unsteady state solution in order to be able understand the
solution. Hence, the pressure contours and mean velocity vectors were drawn at the time
instants corresponding to the occurrences of peak inhalation, end of inspiration and
beginning of expiration, peak expiration and end of expiration and beginning of
inspiration. These were the time instants at which imique patterns of distributions of
velocity and pressure were expected. In order to know the exact time steps at which these
events occurred, the velocity vector distributions were studied from the solutions at each
time step from the start to the end of the respiratory cycle. For instance, the instant when
the maximum magnitudes of velocities were observed at the inlet in the inspiratory ﬂow

direction was considered as the instant corresponding to ‘peak inspiration’.

4.2.1 Velocity and Pressure distributions through the tracheal pressure
cycle

Figure 4.1 shows the velocity ﬁeld during peak inspiration and during intense exercise
for a non-ﬂexed model. The velocity vectors are colored and scaled according to their
magnitudes. The graph in Figure 4.2 indicates the portion of the tracheal pressure cycle
that is completed at the instant when peak inspiration occurs. The beginning of the
tracheal pressure cycle is considered here as the instant at which the pressure at the
trachea is zero and is about to become negative, as dictated by the boundary condition
that was supplied. The shaded area in the Figure represents the part of the tracheal

pressure cycle completed just when peak inspiration occurs.

37

 

 

Tracheal Pressure (Pa)

 

l

. M I I I I L 4r 4 I
30000 10 20 30 40 50 60 70 80 90 100

 

 

Percentage of Completion of Tracheal Pressure Cycle

Figure 4.2. Occurrence of peak inspiration relative to the tracheal pressure cycle

38

The peak velocities were observed near the inlet and the outlet (viz., the nostril and the
tracheal opening). Although one would expect the peak velocities (and hence peak
inspiratory ﬂow) to occur at a time instant when the driving pressure is maximum, this
was not what was observed. The tracheal pressure is at its negative peak at the time
instant corresponding to the completion of 25% of the tracheal pressure cycle. After this
instant, the pressure decreases in magnitude, while still remaining negative until 50% of
the cycle is completed. From Figure 4.2, it can be seen that the peak inspiration occurred
when 40% of the tracheal pressure cycle is completed, when the driving pressure is lower
than that at 25% completion of pressure cycle. This lag between the expected occurrence
of peak inspiration and the actual occurrence of peak inspiration indicates that the effect
of compressibility of the ﬂuid elements is pronounced in time-periodic oscillatory ﬂows,
especially at higher pulsating frequencies. If the ﬂuid elements were totally
incompressible, then a pressure change at the end of the airway should be ‘sensed’

instantly by the ﬂuid elements, thus making the velocities in phase with the pressure.

Yet another interesting observation was that the streamlines followed a curved path
during inspiration in order to maneuver around the swirl created near the nostril area due
to the presence of the turbinates where almost zero velocity was observed. The ﬂow
through the remaining part of the nasal cavity and the nasopharyngeal region was
observed to be devoid of any eddies. Figure 4.3 shows the pressure contours during the
same time instant (peak inspiration). An important observation from this ﬁgure is that the
pressure in the pharynx and the mid-nasal cavity is fairly constant at -1000 Pa (gauge)

and about -2000 Pa at certain pockets around the soft palate and the pharynx.

39

0 00
,1 I. 475 00

-350 00

-525 00

-700 00

~875 00

.105000
4225 00
-1400 00
4575 00
4750 00
.1925 00
-2100 00
-2275 00
4450 00
-2625 00
4800.00
4975 00
-3150.00
-3325 00
4500.00

 

   

Figure 4.3. Pressure contours (Pa) at peak inspiration

 

2000~

  
  
 

E

I
I
a
E

'3

-2000 -

 

 

 

.3mo—L—L__—l¥ I l M m I l
I

10 20 30 40 50 60 70 80 90 100
Percentage of Completion of Tracheal Pressure Cycle

Figure 4.4. Occurrence of peak inspiration relative to the tracheal pressure cycle

40

Figure 4.5 shows the velocity distribution during the end of the inspiration half-cycle and
just at the beginning of the expiratory half-cycle. At this instant, the ﬂuid elements are
about to reverse their directions and go from the trachea to the nostril due to positive
gauge pressure at the trachea. In Figure 4.6, the occurrence of the end of inspiration
relative to the tracheal pressure cycle is shown. However, as is evident from the Figure,
all the ﬂuid elements do not reverse their directions simultaneously. From observing the
ﬁgure, it can be hypothesized that the ﬂuid elements near the trachea start ‘sense’ the
change in tracheal pressure from negative to positive gauge earlier than those ﬂuid
elements near the nostril. Hence, the ﬂuid elements near the pharynx and the tracheal
opening start following an expiratory ﬂow direction whereas those ﬂuid elements near
the nostril still follow an inspiratory ﬂow direction. As a result, swirls are formed due to
the mixing of a set of ﬂuid elements moving in one direction with those moving in the

opposite direction.

Due the lag effect due to compressibility, as discussed earlier, ﬂow reversal from
inspiration to expiration does not occur at the instant at which zero pressure difference
occurs between the nostril and the trachea (which is at 50% of completion of the tracheal
pressure cycle). Instead, this occurs at about 60% into the pressure cycle which is way

past the instant at which the tracheal pressure becomes zero from negative gauge.

Figure 4.7 shows the pressure contours at end of inspiration and beginning of expiration
and Figure 4.8 indicates the occurrence of this even relative to tracheal pressure cycle.

Variation in pressure from the nasal cavity to the tracheal opening appears to be gradual.

41

40 00
:- :1 38 00
. 36 00
34.00
32 00
30 00
28.00
26 00
24.00
22.00
20.00
18.00
16 00
14.00
12.00
10 00
8 00

6 00

4 00

2 00

0 00

rd

 

 

 

Figure 4.5. Velocity distribution (m/s) at end of inspiration and just before
expiration

1AM
J‘MI

 

Tracheal Pressure (Pa)

 

 

l L I J I I

I
10 20 30 40 50 60 70 80 90 100
Percentage of Completion of Tracheal Pressure Cycle

 

 

Figure 4.6. Occurrence of end of inspiration relative to the tracheal pressure cycle

42

n
It
_.
co
~I
or
8

 

 

 

 

 

_ g l 1 I I 4.: I I I I
30000 10 20 30 40 50 60 70 80 90 100

Percentage of Completion of Tracheal Pressure Cycle
Figure 4.8. Occurrence of end of inspiration relative to the tracheal pressure cycle

43

Figure 4.9 shows the distribution of velocity distribution during peak inspiration and
Figure 4.10 and show the time instant in the tracheal pressure cycle during which this
event occurs. Here again, the lag between the expected occurrence and the actual
occurrence of the event is observed. The maximal expiratory driving pressure would
occur at the instant of 75% completion of the tracheal pressure cycle. Hence peak
expiration is found to occur at this time instant. Instead, peak expiration is found to occur
at about 90% completion of the unsteady pressure cycle, when the driving pressure is
much lower compared to that at 75% completion. The streamlines near the nostril take a
dorsal (towards the horse’s backbone) curved path before exit as opposed to the ventral
(towards the horse’s belly) curved path taken by the streamlines during peak inspiration.
The ﬂuid elements appear to take this path because of negotiating a swirl established in
the nasal cavity as seen in Figure 4.9. At this instant, peak velocities were observed near
the inlet and the outlet (viz. the tracheal opening and the nostril). At the pharynx and the
inner nasal cavity, velocity magnitudes are smaller and vary more gradually. This could
be attributed to the increase in cross-sectional area as a ﬂuid element moves from the

trachea to the pharynx during expiration.

Figure 4.11 describes the pressure contours during peak expiration and Figure 4.12
indicates the occurrence of peak expiration relative to the tracheal pressure cycle. The
maximum pressure occurs at the pharyngeal cavity which is at about 1700 Pa (gauge). It
is to be noted that this pressure is higher than the driving pressure at the tracheal opening

which is at about 1200 Pa (gauge).

44

 

Tracheal Pressure (Pa)

Figure 4.9. Velocity distribution (m/s) at peak expiration

 

l I | I

     

     

I I I I l

E

 

       

 

_ ___L___J— I J I 1 I I I
30000 10 20 30 40 50 60 70 80 90 100

Percentage of Completion of Tracheal Pressure Cycle

Figure 4.10. Occurrence of peak expiration relative to the tracheal pressure cycle

45

 

 

Figure 4.11. Pressure contours (Pa) at peak expiration

 

2000L

   

E

Tracheal Pressure (Pa)

 

 

l l L l I

I
10 20 30 40 50 6O 70 80 90 100
Percentage of Coupletion of Tracheal Pressure Cycle

 

Figure 4.12. Occurrence of peak expiration relative to the tracheal pressure cycle

46

Figure 4.13 shows the velocity distribution at the end of the expiration-cycle and at the
beginning of the inspiration-cycle. Figure 4.14 shows the percentage of the tracheal
pressure cycle completed when this event occurs. The ﬂow reverses from exhalation to
inhalation at this instant. Most of the domain has close to zero velocity, except near the
nostril and at the pharynx where there are swirls with velocities ranging from 0 to 16 m/s.
Both of these regions have a geometric similarity — both of them involve rapid expansion
in area. The swirl that is seen at the pharynx is just above the epiglottis muscle. The
velocity in the vicinity of the posterior region of the soft palate and the epiglottis appears

to be greater than the surrounding regions due to the effect of the swirl.

Here again, due to inertial effects there is a lag between the time instant at which this
event is to occur and the time instant at which this actually occurs. Rather than at the start
of the tracheal pressure cycle, beginning of inspiration occurs at about 8% into the

tracheal pressure cycle.

Figure 4.15 shows the pressure contours at end of expiration and beginning of inhalation
and Figure 4.16 indicates the portion of the tracheal pressure cycle completed when this
event occurs. The pressure distribution varies gradually from 100 Pa at the nostril end
and about -1600 Pa at the tracheal end, as indicated by the striations in the contours along

the upper airway.

47

1‘rachealPressure(Pa)

 

10,00
15 20
14 40
13 60
12.80
12.00
11.20
10.40
9 60
8.80
8 00
7.20
6.40
5.60
4.80
4 00
3.20
2.40
1 60
0 80
0 00

 

 

Figure 4.13. Velocity distribution (m/s) at end of expiration and just before
inspiration

annn
JWV I T I I I I I I

 

2000

1000

-1000

-2000

 

 

 

 

_3m l 1 l M A l l l .L
0 10 20 30 40 50 60 70 80 90 100

Percentage ochrrpIetion ofTracheal Pressure Cycle
Figure 4.14. Occurrence of end of expiration relative to the tracheal pressure cycle

48

 

Tracheal Pressure (Pa)

 

‘1 Mn
JWV

 

2000

1000

-2000

 

 

 

I

l 1
0 10 20 30 40 50 60 70 80 90 100

GMA I J J I I I
JWU

 

Percentage of Completion of Tracheal Pressure Cycle

Figure 4.16. Occurrence of end of expiration relative to the tracheal pressure cycle

49

4.2.1. Variation in velocity distribution with respiratory rate

The ﬂow ﬁeld in this problem is dependent upon two broad parameters: Rate of
respiration and degree of artiﬁcially induced ﬂexure. In order to study the sensitivity of
the solution to the rate of respiration alone, the velocity distribution given by solutions
corresponding to rest, mild exercise and intense exercise for a given degree of induced
ﬂexure and at different time instants were plotted and studied. Such velocity vector plots
for no induced ﬂexure model and at end of inspiration are shown in this section. In
section 4.2 it was observed that the time instant corresponding to the end of inspiration
and beginning of expiration had the strongest swirls. Hence, this time instant was chosen

to compare solutions corresponding to the different breathing rates.

Figure 4.17, Figure 4.18 and Figure 4.19 show the velocity distribution corresponding
to rest, mild exercise and intense exercise respectively at end of inspiration (beginning of
expiration) for the non ﬂexed geometry. Each of these Figures has been plotted to a

common scale so that they can be compared with ease.

50

 

Figure 4.17. Velocity distribution (m/s) at end of inspiration - rest condition

 

Figure 4.18. Velocity distribution (m/s) at end of inspiration - mild exercise
condition

51

 

Figure 4.19. Velocity distribution (m/s) at end of inspiration - intense exercise
condition

From the above comparisons, one can observe that at rest, the swirl formed in the anterior
nasal cavity has velocities ranging from 0 — 18 m/s. Whereas, during mild and intense
exercise conditions, the swirl(s) in the vicinity of the soft palate is (are) not only larger in
terms of physical dimensions, but are also more signiﬁcant in terms of the velocity
gradients (0-30 m/s for mild exercise and 0-40 m/s for intense exercise). It’s also evident
from this comparison that the velocity at any point increases with increase in the
respiration rate, which is quite obviously due to greater driving forces and higher

frequencies.

4.2.2. Variation of velocity distribution with flexure in upper airway

In order to study the sensitivity of the solution to the degree of induced ﬂexure alone, the
velocity distribution given by solutions corresponding to non-ﬂexed, mildly ﬂexed and

extremely ﬂexed geometries for a given rate of respiration at different time instants are

52

plotted and studied. These velocity plots during intense exercise and at the instant
corresponding to end of inspiration are shown in this section. Figure 4.20, Figure 4.21
and Figure 4.22 show the velocity distribution corresponding to non-ﬂexed model,
model with mild induced ﬂexure and model with extreme induced ﬂexure respectively at
end of inspiration (begimring of expiration) for intense exercising condition. Each of
these Figures has been plotted to a common scale so that they can be compared with ease.
Unlike what was observed in the previous section, the intensity of the swirls, in terms of
velocity range was not found to vary signiﬁcantly in a deﬁnite way between non-ﬂexed,
mildly ﬂexed and extremely ﬂexed geometries. This might indicate that the solution is

more sensitive to the rate of respiration than it is to the degree of induced ﬂexure.

 

Figure 4.20. Velocity distribution (m/s) at end of inspiration — no induced ﬂexure

53

 

Figure 4.22. Velocity distribution (m/s) at end of inspiration — extreme induced
ﬂexure

54

4.3. Wall pressure along the dorsal surface of lower wall

As mentioned in chapter 1, in a horse with the DDSP condition, the soft palate which is
usually held in position by the epiglottis muscles is displaced upwards into the upper
airway. Hence, it expected that there will be a net upward force on the soft palate of the
lower wall of the geometry in general, at least during brief periods of time in a respiratory
cycle. In order to understand if there is a lift created on the lower wall of the geometry
and to know the time instants at which this situation occurs, wall pressures along the
lower wall were plotted from the nostril end to the tracheal end at different time instants.
One such plot is shown in Figure 4.23 where the wall pressures from the solution
corresponding to extremely ﬂexed geometry under intense exercise condition are plotted
at different time instants throughout the breathing cycle. If T is the time period of the
tracheal pressure function, the lower wall pressures are plotted at the following time
instants: when t = T, 0.1 T, 0.5 T, and 0.6 T. The relative positions of the various
anatomical features relevant to the current discussion are also marked along the length of

the lower wall.

55

 

 

 

 

 

 

 

 

 

2000» I l I I I I I I I ,5, ---Att 2 T
At t = 0.1T
- '"Att = 0.3T
—AI 1 = 0.5T
IOOOI """"""""""""""" - i d """"" Atr=0.6T
,,,, ("~""---—--. ‘ 3". 3 L: E
’/J ................. N”, : g
I, .I 51“¢~‘§ :
A 0L j i" ‘ _
a _ g .........
a I . ............... 5
g ..... ‘ ' i
E P. \ A. II
I
2 1m '1. \\ I,
o - .1. .. _‘
I': -"‘"‘ - 'I
=. I‘? A _ ﬂ _ - “ .r ..... ~ _ , - w“ . ,1 ,.
3 I ’I.,¢-"' ‘1‘“; ‘19 l
'2m0r .\ i: “-.‘./s- I, E I? ||' _
X =; I "I I I I II I.
\’ I I ’ ‘ t
\ ‘ I I ' l‘ _i
. .a. | ‘ ”I“ (I \J
I l.‘ I!
i I!
.3000. . 3 d
I
4000 1 1 I l I 1 1 1 1
0 IO 20 30 4O 50 60 70 80 90 100
Percentage Curvelength
Narcs Hard Palate Soft Epiglottis Trachea
Palate

Figure 4.23. Wall pressure along dorsal surface - variation through the tracheal
pressure cycle

From the above Figure, it can be observed that the lowest pressure along the bottom wall
occurred when at about 30% of completion of the pressure cycle. At this instant, the ﬂow
was observed to be close to peak inspiration (from section 4.1). The pressure at the dorsal
surface of the bottom wall exhibited time-periodic variation. Near the nostril at time
instants between 0.3 T and 0.6 T, high gradients of wall pressure were observed. This
could be due to the entrance effects contributing to increased turbulence when the ﬂow
rates are close to being peak inspiratory ﬂow rates. Over the most of the hard palate, the

wall pressure was observed to remain either near-constant or vary gradually, except at the

56

end connecting to the soft palate. In that region and the one that continues further into
about 30% of the soft palate, high wall pressure gradients were observed. This could be
attributed to the rapid change in the cross-sectional area and also due to the ﬂexure in the
geometry (even the non—ﬂexed model has a mild natural ﬂexure) resulting in curvature of
the streamlines and hence in the formation of eddies contributing to greater turbulence.
High magnitudes of wall pressure gradients were found to occur at the portion of the

trachea connecting to the epiglottis as well.

In general, lower pressures were found to occur at times close to those corresponding to
the peak inspiration times (0.3 T). At this time instant, pressures were found to reach up
to maximum negative values of -3000 Pa in regions adjacent to the soft palate and the
nostrils, and up to -4000 Pa in the regions near the epiglottis-tracheal interface. For these

times, the wall pressure gradients were also found to be greater.

4.3.1. Variation of wall pressure distribution with respiratory rate

To understand the effect of the rate of respiration on the lift created on the lower wall and
the soft palate in particular, the lower wall pressures at time instant t = 0.3 T were plotted
(Figure 4.24) from the solutions corresponding to non-ﬂexed geometry for rest, mild
exercise and intense exercise conditions. The time instant corresponding to t = 0.3 T was
chosen because from the earlier section, it was observed that the least wall pressures

occurred at this time step.

57

 

 

 

 
  
 

 

 

 

 

 

\\ -------------- v ...................... "'11“
I \‘\o"-‘\ I """""""""""""""""""" ,4 "Mild
"I v .
-500 i‘ ' \ _ ,, \ Exercise
I ,_,---..\,-‘ ” "'V'V‘ \\ [3m
\ I I-" \ -\ I" \ — .
\\ ,. "’1 If II \\. Exercme
.1000 r ‘1! ‘ \\ II I” IIIIIIIII I! i \I\\\_ ,"' I
9,:
1? r
E-ISOOP i a
3
g 4000+ q
'3
3
-2500~ A
-3000~ _
_ I I I I I I I I I
35000 10 20 30 40 50 60 70 80 90 100
Percentage Curvelength
Nares -1 Hard Palate Soft Epiglottis“ Trachea A
Palate

Figure 4.24. Wall pressure along dorsal surface - variation with rate of respiration

It can be inferred from Figure 4.24 that for a given ﬂexure, the wall pressure during peak
inspiration decreases with intensity of respiration. Moreover, wall pressure gradients are
also found to be more pronounced at intense exercise conditions. This might explain
could possibly explain why DDSP should, quite often, manifest only during intense

exercise.

58

4.3.2. Variation of wall pressure distribution with induced flexure in the
upper airway

To understand the effect of artiﬁcially induced ﬂexure on the lift created on the lower
wall and the soft palate in particular, the lower wall pressures at time instant t = 0.3 T
were plotted (Figure 4.25) from the solutions corresponding to intense exercise condition
for the non-ﬂexed and the extremely ﬂexed models. Here again, the time instant
corresponding to t = 0.3 T was chosen because from the section 4.2, it was observed that

the least wall pressures occurred at this time step.

 

 

 

 

 

 

 

 

 

 

 

0 I I I I I I r I
__No
Flerau'e
_Extreme
Heme
53’
d:
*3
3
4m 1 l l l I 1 L l l
0 10 20 30 40 50 60 70 80 90 100
Percentage Curvelengh
A Nares U HardPalste -4 Soft Epiglottis“ Trachea -

 

Palate

Figure 4.25. Wall pressure along dorsal surface - variation with induced flexure

59

Figure 4.25 clearly shows that with a greater ﬂexure in the airway, the pressure gradients
along the wall, especially around the soft palate and the epiglottis are greater. This could
be attributed to the additional curvature of the streamlines. The soft palate region, the
epiglottis region and the entrance to the trachea region of the lower are observed to have

high pressure gradients (gradients of wall pressure along the surface of the lower wall).

60

CHAPTER 5: CONCLUSIONS AND FUTURE WORK

The equations governing the respiratory ﬂow of air within the equine upper airway were
solved numerically. Solutions were obtained for airﬂow under different rates of
respiration viz. rest, mild exercise and intense exercise and corresponding to different
airway geometries viz. non-ﬂexed, mildly ﬂexed and extremely ﬂexed. Each of these
solutions was analyzed by plotting the velocity and pressure distributions at various time
instants throughout the respiratory cycle. The pattern of airﬂow established was found to

be uniquely determined by the upper airway geometry and by the rate of respiration.

At peak inspiration during intense exercise, the pressures in the pharynx, the posterior
and mid-nasal cavity were fairly constant at a value of -1000 Pa with pockets of even
lower pressures reaching up to about -2000 Pa at certain regions along the dorsal surface
of the lower wall that are adjacent to the soft palate, the epiglottis and anterior pharynx.
Wall pressures along the dorsal surface of the lower wall were found to reach up to
maximum negative values of -3000 Pa adjacent to the soﬁ palate and the nostrils, and up
to -4000 Pa near the epiglottis-tracheal interface. The lift forces set up on the lower wall
due of these highly negative wall pressures and luminal pressures could reach such high
values that might not be balanced by the downward forces applied on the pliant structures
like soft palate by the restraining muscles of the upper airway. Hence the soft palate and
other such pliant tissues of the upper airway may have a natural tendency to get displaced

into the airway.

61

In addition to being subjected to highly negative wall pressures and luminal pressures, the
regions along the dorsal wall adjacent to the soft palate and the epiglottis-tracheal
interface are also subjected to high wall pressure-gradients. These sharp variations in wall
pressure may be instrumental in destabilizing the muscles of the soft palate and/or the

epiglottis, facilitating the dorsal displacement of the soft palate.

The velocity plots drawn at various instants through the respiratory cycle indicate that in
the nostril, the anterior nasal cavity and the posterior pharyngeal cavity, the ﬂow is
complex with swirls observed frequently. The intensity of these swirls in terms of ranges
of velocity magnitudes increase with increase in respiration rate and driving pressures.
Since increase in swirling motion is often linked with increase in the intensity of
turbulence, it is possible that the increased magnitudes of ﬂuctuating components of
velocities during intense exercise may be inducing to the soft palate muscles, an unsteady
ﬂuttering motion, or a tendency to ﬂutter, thereby contributing to dorsal displacement of

soft palate.

The phenomena mentioned above viz. swirling ﬂows in the posterior pharyngeal cavity
and highly negative dorsal wall pressures and dorsal wall pressure-gradients adjacent to
the soft palate-epiglottis region were observed to be more pronounced under conditions
of intense exercise. This might explain why the condition of DDSP is not diagnosable in
horses at rest. Similarly, the highly negative dorsal wall pressures and dorsal wall
pressure-gradients adjacent to the soft palate-epiglottis region were more pronounced

with increase in degree of induced ﬂexure in the upper airway. This might explain why

62

the condition of DDSP is at times not diagnosable through tests at high speeds on a
treadmill in a horse with a history of the condition while racing, since this testing
situation is not representative of an actual racing situation where there is a ﬂexure
introduced in the airway. This may also explain why this condition in many cases is
manifested almost exclusively while racing as under this situation, not only do rapid
breathing rates prevail, but also an additional ﬂexure to the airway is introduced by the

pulling of the reins.

It is possible that the condition of DDSP is initiated by strong negative pressures along
the dorsal surface of the soft palate and that the displacement of the soft palate into the
airway causes a reduction in the cross-sectional area of the upper airway which in turn
decreases the pressure inside the airway in the vicinity of the soft palate, facilitating a
dynamic collapse of the soft palate into the airway. Whether DDSP will occur or not,
given a prevalence of large negative dorsal wall pressures along the soft palate would
depend on other factors like — stiffness of the soft palate and its holding muscles
including the epiglottis (which in turn may be inﬂuenced by other factors including the
age of the horse), existence of other upper airway respiratory disorders such a pharyngeal

collapse, existence of certain neuromuscular disorders, etc.

From this study, it is possible to conclude that the pattern of airﬂow established is at least

partly responsible for the dorsal displacement of soft palate. Intense exercise and high

degree ﬂexure in the upper airway may be favorable factors for DDSP to occur. A

63

qualitative estimate of the individual contribution of these factors in causing DDSP was

beyond the scope of this study.

Our analysis assumed the ﬂow to be varying only along the plane, deﬁned by the sagittal
section of the geometry. A three dimensional representative geometry would provide
better insights into the ﬂow simulations. The soft palate portion of the lower wall could
be modeled to represent a pliant tissue instead of a rigid wall. This might help one
visualize the occurrence of DDSP. For the sake of simplicity, a sinusoidal function of
time was applied to describe the unsteady pressure boundary condition at the tracheal
boundary. An analytical function of time, constructed from curve-ﬁtting time history data
of the actual outlet pressure values recorded from experiments, to describe the variation
in the driving pressure. It might be worthwhile to study the effects of the entry
conditions of the air entering the nostril on the distribution of pressure and velocity
throughout the domain. That study might throw light on the dependence of DDSP on the
racing speed and the wind conditions. If more computational resources are available, one
might apply a more sophisticated turbulence model in order to resolve the near-wall
regions adequately. A physical model of the upper airway could be constructed with
which, the ﬂow could be studied in detail using ﬂow visualization techniques. With a
complete understanding of how the ﬂuid mechanics within the upper airway causes the

condition of DDSP, it might be possible to develop a treatment for this condition.

64

10.

REFERENCES

. Parente, E.J., Martin, B.B., Tulleners, E.P., Ross, M.W., 2002, “Dorsal

Displacement of the Soft Palate in 92 Horses During High-Speed Treadmill
Examination (1993-1998),” Vet. Surg., 31, pp. 507-512.

Ducharme, N.G., Hackett, R.P., Woodie, J.B., Dykes, N., Erb, H.N., Mitchell,
L.M., Solderholm, L.V., 2003, “Investigations into the role of the thyrohyoid

muscles in the pathogenesis of dorsal displacement of the soft palate in horses,”
Equine Vet. J ., 35, pp. 258-263.

Derksen, F.J., Holcombe, S.J., Hartmann, W., Robinson, N.E., Stick, J.A., 2001,
“Spectrum Analysis of respiratory sounds in exercising horses with

experimentally induced laryngeal hemiplegia or dorsal displacement of the soft
palate,” Am. J. Vet. Res., 62, pp. 659-664.

Martin, B.B., Reef, V.B., Parente, E.J., Sage, AD, 2000, “Causes of poor
performance of horses during training, racing, or showing: 348 cases (1992-
1996),” J. Am. Vet. Med. Assoc., 216, pp. 554-558.

Duggan, V.B., MacAllister, C.G., Davis, M.S., 2002, “Xylazine-induced
attenuation of dorsal displacement of the soft palate associated with epiglottic
dysfunction in a horse,” J. Am. Vet. Med. Assoc., 221, pp. 399-402.

Holcombe, S.J., Comellisse, C.J., Berney, C., Robinson, NE, 2002,
“Electromyographic activity of the hyoepiglotticus muscle and control of
epiglottis position in horses,” Am. J. Vet. Res., 63, pp. 1617-1621.

Llewellyn, H.R., Petrowitz, A.B., 1997, “Stemothyroideus Myotomy for the
Treatment of Dorsal Displacement of the Soft Palate,” Proc. Am. Assoc. Equine
Pract. , 43, pp. 239-243.

Robinson, N.E., Sorenson, P.R., Goble, DO, 1975, “Patterns of airﬂow in normal
horses and horses with respiratory disease,” Proc. Am. Assoc. Equine Pract., 21,
11-21

Rehder, R.S., Ducharme N.G., Hackett R.P., Nielan G.J., 1995, “Measurement of
upper airway pressures in exercising horses with dorsal displacement of the soft
palate,” Am. J. Vet. Res., 56, pp. 269-274.

Kepler, G.M., Richardson, R.B., Morgan, K.T., Kimbell, 1.8., 1998, “Computer
Simulation of Inspiratory Nasal Airﬂow and Inhaled Gas Uptake in a Rhesus
Monkey,” Toxicology and App. Pharmacology, 150, pp. 1-11.

65

11.

12.

13.

14.

15.

16.

Kimbell, J.S., Subramaniam, R.P., Gross, E.A., Schlosser, P.M., Morgan, K.T.,
2001, “Dosimetry Modeling of Inhaled Formaldehyde: Comparison of Local Flux

Predictions in the Rat, Monkey and Human Nasal Passages,” Toxicological
Sciences, 64, pp 100-110.

Anderson, J.D., Tulleners, E.P., Johnston, J.K., Reeves, M.J., 1995,
“Sternothyrohyoideus myectomy or staphylectomy for treatment of intermittent

dorsal displacement of the soft palate in racehorses: 209 cases (1986-1991),” J.
Am. Vet. Med. Assoc., 206, pp. 1909-1912.

Lumsden, J .M., Derksen, F.J., Stick, J.A., Robinson, N.E., 1993, “Use of ﬂow-
volume loops to evaluate upper airway obstruction in exercising Standardbreds”,

54, 5, pp 766-775.

Pope, B., 2000, Turbulent Flows, 1St ed., Cambridge University Press, Cambridge,
UK.

Lecture notes, “Turbulence Modeling and Simulation”, fall 2007, Mechanical
Engineering, MSU.

Fluent 6.2. User’s Guide, 2005, Fluent Inc., Lebanon, New Hampshire, USA.

66

.

 

 

 

 

..
n
m
m
u

LIBRARIES

SIY’

R
E
N
d
,L

ICHIGAH S”: I

1IIIIIIII

I

III

1293 02956 4873