Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A novel SPH method for investigating the role of saliva in swallowing using 4D CT images Ho, Andrew Kenneth 2017

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

Item Metadata


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

Full Text

A Novel SPH Method forInvestigating the Role of Saliva inSwallowing using 4D CT imagesbyAndrew Kenneth HoB.E.Sc., The University of Western Ontario, 2008M.E.Sc., The University of Western Ontario, 2010A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Electrical & Computer Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)April 2017c© Andrew Kenneth Ho 2017AbstractThe thesis presents novel computer methods towards simulation of oropha-ryngeal swallowing. The anatomy and motion of the human upper airwaywas extracted from dynamic Computed Tomography (CT) data using a noveltool and workflow. A state-of-the-art SPH method is extended to accommo-date non-Newtonian materials in the extracted geometries. A preliminarynumerical experiment of six human oropharyngeal swallows using SmoothedParticle Hydrodynamics (SPH) demonstrates that the methods are robustand useful for simulation of oropharyngeal swallowing.The presence of saliva is well known to be important for mastication,swallowing, and overall oral health. However, clinical studies of patientswith hyposalivation are unable to isolate the effect of saliva from other con-founding factors. The simulation presented in this thesis examines fluidboluses under lubricated and non-lubricated boundary conditions. Uponcomparison with medical image data, the experiments suggest that salivadoes not provide a significant lubricative effect on the bolus transit times,but it may serve to reduce residue and therefore improve overall swallowingefficacy. Our findings, while preliminary, corroborate with existing clinicalresearch that finds that groups with hyposalivation do not have significantlydifferent transit times with control groups, but that residue may be increasediiAbstractin the hyposalivation group.Previous studies using computer simulation of fluid flow in the orophar-ynx typically make use of simplified geometries. Our work uses dynamic320-row Area Detector Computed Tomography (ADCT) images as the ba-sis for the simulations, and therefore does not require simplifying geometricassumptions. Since the data are dynamic, motion trajectories are all sup-plied by the ADCT data, and extrapolation from 2D sources such as bi-planevideofluoroscopy is not required. Processing the image data required the de-velopment of a novel workflow based on a new tool, which we call BlendSeg.We utilize and extend Unified Semi-Analytic Wall (USAW) SPH methodsso that orophrayngeal swallowing simulations may be performed. Theseextensions include the simulation of non-Newtonian boluses, and moving3D boundaries. Partial validation of the extended USAW SPH method isperformed using canonical flows.iiiPrefaceA swallowing proof-of-concept model using SPH to represent a bolus in abiomechanical model of a human head was published in [55].The extensions to USAW SPH described in Chapter 3 were performedby myself in collaboration with Prof. Sheldon Green and Prof. Sidney Fels.Prof. Green and Prof. Fels advised me on the suitable assumptions fororopharyngeal swallowing, and the extensions required in USAW SPH. Theyalso advised me on suitable experiments to run for the partial validation.I derived the extensions for the moving boundaries in 3D, and the non-Newtonian materials. I also performed the validation experiments and wrotethe chapter.The experiments in Chapter 4 were conducted at UBC by myself incollaboration with Dr. Yoko Inamoto, Dr. Eiichi Saitoh, Prof. RebeccaAffoo, Prof. Nicole Rogus-Pulia, Prof. Mark Nicosia, Prof. Sheldon Green,and Prof. Sidney Fels. The 320-row ADCT data was gathered by Dr.Inamoto, Dr. Saitoh, and their colleagues at Fujita Health University inJapan. They allowed us to use the data for our numerical simulations andprovided verification of the airway identification using BlendSeg. They alsogenerated the CT images and provided feedback on the chapter. Dr. Affooand Dr. Rogus-Pulia performed the bolus transit time measurements of theivPrefaceADCT data and the simulations. They also provided clinical analysis of theresults, and feedback on the chapter. Prof. Nicosia, Prof. Green and Prof.Fels provided support with the interpretation of fluid simulation results, aswell as feedback on the chapter. I extracted the airway from the ADCTimages using BlendSeg, and ran the extended USAW SPH simulation of theslip and no-slip simulations. I also wrote the bulk of the chapter. The workin Chapter 4 was approved by the UBC Clinical Research Ethics Boardunder certificate no. H16-01546.A version of section 4.3 describing a workflow for extracting dynamicboundaries from 4D data using BlendSeg, has been published in [54]. Dr.Yoko Inamoto and Dr. Eiichi Saitoh provided the 320-row ADCT data. Iwas responsible for the design of the workflow, the creation of the BlendSegtool, and also wrote the manuscript.The rheology measurements of the bolus, described in Section 4.4, wereperformed at UBC in Prof. Dana Grecov’s rheology lab, by Nick Yeh andmyself. Nick performed the rheology measurements using the KINEXUSrheometer. I prepared the samples, and fit a Cross model to the data.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xixDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xx1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Computer simulation of swallowing . . . . . . . . . . . . . . 31.2 Smoothed Particle Hydrodynamics . . . . . . . . . . . . . . . 41.3 Oropharyngeal bolus simulation with the extended USAWSPH method . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.4 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 72 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10viTable of Contents2.1 Anatomy and physiology of swallowing . . . . . . . . . . . . 102.1.1 Imaging of swallowing . . . . . . . . . . . . . . . . . . 132.1.2 Saliva . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.2 Fluid simulation . . . . . . . . . . . . . . . . . . . . . . . . . 202.2.1 Mesh-free methods . . . . . . . . . . . . . . . . . . . 232.2.2 Smoothed Particle Hydrodynamics . . . . . . . . . . 242.3 Fluid simulations of swallowing . . . . . . . . . . . . . . . . . 292.4 Segmentation and modelling from 4D imaging data . . . . . 312.5 ArtiSynth biomechanics toolkit . . . . . . . . . . . . . . . . . 333 Smoothed Particle Hydrodynamics for fluid flow . . . . . . 343.1 Classic Smoothed Particle Hydrodynamics formulation . . . 353.1.1 Spatial derivatives . . . . . . . . . . . . . . . . . . . . 373.1.2 SPH for the incompressible, viscous, Navier-Stokesequations . . . . . . . . . . . . . . . . . . . . . . . . . 393.1.3 Boundary handling in SPH . . . . . . . . . . . . . . . 413.1.4 Smoothing pressures with δ-SPH . . . . . . . . . . . . 423.2 Unified Semi-Analytic Wall SPH . . . . . . . . . . . . . . . . 433.2.1 Solid boundaries and insufficient domain representa-tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 443.2.2 Moving boundaries in 3D . . . . . . . . . . . . . . . . 513.2.3 Non-Newtonian fluids in USAW SPH . . . . . . . . . 523.2.4 Time integration . . . . . . . . . . . . . . . . . . . . . 543.3 Verification and validation . . . . . . . . . . . . . . . . . . . 543.3.1 Couette flow . . . . . . . . . . . . . . . . . . . . . . . 55viiTable of Contents3.3.2 Startup Hagen-Poiseuille flow for a Newtonian fluid . 573.3.3 Lid-driven cavity . . . . . . . . . . . . . . . . . . . . 603.3.4 Non-Newtonian Hagen-Poiseuille flow . . . . . . . . . 613.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 654 Fluid simulation of oropharyngeal swallowing . . . . . . . . 674.1 The role of saliva in swallowing . . . . . . . . . . . . . . . . . 684.2 320-row ADCT swallowing sequences . . . . . . . . . . . . . 694.3 Deriving geometric boundaries from 4D ADCT data . . . . . 704.3.1 BlendSeg . . . . . . . . . . . . . . . . . . . . . . . . . 734.3.2 Initial mesh generation using Amira . . . . . . . . . . 734.3.3 Visualizing 3D volumetric data in Blender . . . . . . 744.3.4 Generating intersection contours . . . . . . . . . . . . 754.3.5 Interpolating in time to create a moving boundary . . 774.3.6 Extracting the boundaries . . . . . . . . . . . . . . . 784.4 Viscosity of thickened boluses . . . . . . . . . . . . . . . . . 794.4.1 Slip and no-slip boundary conditions . . . . . . . . . 814.4.2 Bolus measurements . . . . . . . . . . . . . . . . . . . 824.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 834.5.1 Thickened boluses . . . . . . . . . . . . . . . . . . . . 854.5.2 Thin boluses . . . . . . . . . . . . . . . . . . . . . . . 924.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 964.6.1 Thickened boluses . . . . . . . . . . . . . . . . . . . . 964.6.2 Thin boluses . . . . . . . . . . . . . . . . . . . . . . . 98viiiTable of Contents4.6.3 Corroboration with existing simulation and clinical re-search . . . . . . . . . . . . . . . . . . . . . . . . . . . 994.6.4 Limitations and assumptions . . . . . . . . . . . . . . 1014.6.5 Future directions . . . . . . . . . . . . . . . . . . . . 1024.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1035 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1055.1 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1075.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . 110Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113AppendicesA Estimate of effect of saliva on a gravity driven bolus . . . 133B Using full-slip to approximate a lubricative saliva layer . . 138B.1 Problem description and assumptions . . . . . . . . . . . . . 140B.1.1 Steady-state solution . . . . . . . . . . . . . . . . . . 141B.1.2 Finite-Volume time-dependent solution . . . . . . . . 143B.1.3 Code verification . . . . . . . . . . . . . . . . . . . . . 146B.2 Time dependent flows . . . . . . . . . . . . . . . . . . . . . . 148B.2.1 Measured values of saliva thickness and viscosity . . . 151ixList of Tables3.1 L2 Error and apparent convergence for Hagen-Poiseuille flowsimulation. Samples were for a slice of particles near thecentre of the pipe perpendicular to the flow direction. . . . . 593.2 L2 Error and apparent convergence for a power-law Hagen-Poiseuille flow simulation. Samples were for a slice of particlesnear the centre of the pipe perpendicular to the flow direction. 644.1 Scans from three female subjects provided the data for thissimulation study. All subjects were healthy with no historyof dysphagia, and not taking any medication. . . . . . . . . . 704.2 Parameters of the Cross model for the Nectar- and Honey-thick simulations. . . . . . . . . . . . . . . . . . . . . . . . . . 844.3 Summary of transit times and residue for CT, and simulatedslip/no-slip for each swallow sequence. . . . . . . . . . . . . . 864.4 Summary of transit times and residue for thin boluses. CT,simulated slip, and no-slip for each swallow sequence. Somesimulated boluses were “aspirated”, and percentage volumeis indicated as well. . . . . . . . . . . . . . . . . . . . . . . . . 94xList of TablesB.1 L2 error convergence w.r.t. mesh size. Ra = Rs = Rb = 0.01,Ns = Nb, µb = mus = 1.0. . . . . . . . . . . . . . . . . . . . . 148B.2 L2 error convergence w.r.t. mesh size. Ra = 1., Rs = 0.001,Rb = 0.01, Ns = Nb, µb = 1, µs = 0.01. . . . . . . . . . . . . . 148xiList of Figures2.1 An illustrated lateral view of a cut-away human head showingthe oral, pharyngeal and laryngeal region with major struc-tures labeled [41]. . . . . . . . . . . . . . . . . . . . . . . . . . 112.2 Videofluoroscopic image of subject swallowing barium. Theposterior of the bolus, contained in the oral cavity, has lowimage intensity and is visible on the left side of the image.Image by Hellerhoff, distributed under a CC BY 3.0 license [50]. 142.3 Fibreoptic Endoscopic view of the larynx. The epiglottis,base of tongue, pharynx, and glottis are visible. Image byMed Chaos, distributed under a CC BY-SA 3.0 license [16]. . 162.4 320-row area detector CT mid-sagittal slice of swallowing.4 ml nectar-thick bolus is visible in the oral cavity with highimage intensity. Image courtesy of Fujita Health University [57]. 162.5 Realtime-MRI of mid-sagittal slice of a normal subject swal-lowing pineapple juice. Image by Martin Uecker, distributedunder a CC BY-SA 3.0 license [118]. . . . . . . . . . . . . . . 17xiiList of Figures2.6 Ultrasound image showing sagittal view of hyoid (low im-age intensity) and thyroglossal duct cyst. Image by LaughlinDawes, distributed under a CC BY-NC 2.5 license [47]. . . . . 183.1 Example 2-D particle distribution for an SPH fluid particle(black circle) with full support over its domain Ω (dottedarea). The domain boundary is designated by ∂Ω (dashedcircle). A typical 1-dimensional kernel function W (x) is over-laid. Notice that W (x) = 0 for x ∈ ∂Ω. . . . . . . . . . . . . 373.2 Example particle distribution for an SPH fluid particle (blackcircle) with truncated support over its domain Ω (dottedarea). The domain boundary is designated by ∂Ω (dashedcircle). A typical 1-dimensional kernel function W (x) is over-laid. Notice that unlike Figure 3.1, W (x) 6= 0 for ∀x ∈ ∂Ω. . 443.3 This figure illustrates the difficulty in getting an estimate for∇v · n for the differential surface element dS, using a fluidpoint (filled circle) perpendicular to the surface normal, n.This situation occurs for geometries where the fluid is con-tained in a non-convex boundary. There is inadequate infor-mation to estimate any shear components since ∆y betweendS and the fluid particle is nil. . . . . . . . . . . . . . . . . . 48xiiiList of Figures3.4 The vector∇v·n is assumed to be constant over ∂Ω for a fluidparticle a where it is truncated by the wall. To compute ∇v ·n, the point on the wall nearest to a, denoted g and indicatedwith a grey circle, is used in the computations. Additionally,it is assumed that the normal n points from g to a. Thisassumption is good near solid walls with low curvature, butgives the appearance of a “smoothed” wall near sharp, convexcorners. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.5 Simulation results vs. exact solution for Couette flow betweentwo parallel planes. Ten particles span the distance betweenthe top and bottom planes. . . . . . . . . . . . . . . . . . . . 563.6 Particle distribution for Hagen-Poiseuille in a plane perpen-dicular to the flow direction for initial spacing of approxi-mately 0.05. The tube wall is shown as a dashed circle. . . . 583.7 The rate of convergence for Hagen-Poiseuille flow appears tobe approximately first-order. . . . . . . . . . . . . . . . . . . 593.8 Hagen-Poiseuille flow simulation vs exact solution for themoderate (a) and high (b) resolution spacing. For the mod-erate resolution simulations, the particle spacing is approxi-mately 0.05 and 40 particles span the pipe diameter. For thehigh resolution simulations the particle spacing is approxi-mately 0.0125 and 160 particles span the diameter. Only halfof the pipe is shown due to symmetry. . . . . . . . . . . . . . 60xivList of Figures3.9 Velocity profiles for the lid-driven cavity flow for Re= 100(top row) and Re= 1000 (bottom row). Y-velocity throughthe horizontal centre of the cavity (left column) and x-velocitythrough the vertical centre of the cavity (right column). TheSPH results show good agreement with the results of Ghia(1982). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 623.10 The streamlines for the lid-driven cavity experiment for Re =100 (left) and Re = 1000 (right). A speed-of-sound of c = 5was used, and the simulation was run for 60 time units. Thepseudo-2D simulation was seeded with 100× 100 (Re = 100)and 200 × 200 (Re = 1000) particles. Running the Re =100 simulation with 200 × 200 resolution resulted in nearlyidentical results. . . . . . . . . . . . . . . . . . . . . . . . . . 633.11 The rate of convergence for power-law Hagen-Poiseuille flowalso appears to be approximately first-order. . . . . . . . . . . 653.12 Hagen-Poiseuille flow simulation for a power-law fluid plottedagainst the exact steady-state solution. The power-law fluidhas K = 1000 and n = 0.8. . . . . . . . . . . . . . . . . . . . 664.1 The segmented initial mesh (translucent blue) overlaid on themid-sagittal slice to illustrate the anatomical area being con-sidered. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 724.2 Generating an initial mesh from time t=0. Mid-sagittal slicesare shown. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74xvList of Figures4.3 Blender’s quad-view, when combined with BlendSeg, give theusual anterior-posterior, superior-inferior, and lateral views(bottom left, top left, and bottom right, respectively). Theupper right view gives an interactive 3D view which can beused to select which slice is being visualized, as well as to per-form sculpting on the mesh. Here, the mesh is being hiddenso that the intersections (orange) can be seen more clearly. . 754.4 Example mid-sagittal slice from time t=0.5s showing the in-tersection of the mesh with the slice. . . . . . . . . . . . . . . 784.5 Plot of apparent viscosity vs. shear-rate for nectar- and honey-thick boluses at a temperature of 10◦C. Note that these arelog-log axes. The materials have very similar curves, and areboth highly shear thinning. . . . . . . . . . . . . . . . . . . . 844.6 Timing chart for 4ml nectar-thick bolus, reclined position. . . 864.7 Timing chart for 10ml honey-thick bolus, semi-prone position 874.8 Timing chart for 10ml honey-thick bolus, reclined position . . 874.9 Nectar-thick bolus, semi-prone position. Columns (L-R): slip,3D CT, no-slip. . . . . . . . . . . . . . . . . . . . . . . . . . . 884.10 Honey-thick bolus, semi-prone position. Columns (L-R): slip,3D CT, no-slip. The shape of the air-bolus interface at thebolus head shows better agreement to the CT images in theno-slip condition than the slip condition in the oral cavity andpharynx. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91xviList of Figures4.11 Detail of bolus-air interface at a single mid-sagittal slice, atthe initiation of swallowing. The interface of the no-slip bolusshows better agreement with the original images than the slipbolus. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 924.12 Honey-thick bolus, reclined position. Columns (L-R): slip,3D CT, no-slip. The slip bolus escapes the oral cavity veryeasily as shown in the second and third rows. The no-slipbolus remains in the oral cavity, showing better agreementwith the CT images. . . . . . . . . . . . . . . . . . . . . . . . 934.13 Timing chart for 10ml thin bolus, reclined position. . . . . . . 944.14 Timing chart for 10ml thin bolus, semi-prone position . . . . 954.15 Timing chart for 10ml thin bolus, reclined position . . . . . . 95A.1 Simplified geometry of a gravity driven bolus with lubricatingsaliva layer. Wall is on the left, gray arrows represent thevelocity of the bolus. The g arrow shows the direction ofgravity. Ls and Lb are the widths of the saliva and bolus, resp.134A.2 Fluid element (dashed box) under consideration with heightH depth D into the page. τw is the shear stress at the wall,and τ(x) is the shear stress at position x. Mg is the force ofgravity on the element. . . . . . . . . . . . . . . . . . . . . . . 134xviiList of FiguresB.1 Simplified geometry of a gravity driven bolus with lubricatingsaliva layer. Wall is on the right, gray arrows represent thefluid velocity. The g arrow shows the direction of gravity.Ra, Rb and Rs are the widths of the air core, bolus, andsaliva, resp. Axial symmetry is assumed, with the centre ofthe air core at r = 0. The geometry has a total radius ofR = Ra +Rb +Rs. . . . . . . . . . . . . . . . . . . . . . . . 139B.2 Theoretical startup flow solution vs. FV solver for Newtonianflow in a pipe. Ra = 0, Rb = Rs = 0.02, Ns = Nb = 5,µb = µs = 1, ρ = 1000, g = 9.81. . . . . . . . . . . . . . . . . 147B.3 Time dependent solutions of a two-material start-up pipeflow at various times. Low viscosity saliva lubricates a high-viscosity bolus in the centre of the pipe. For this plot, theparameters were N = 20, Ra = 0, Rb = 10−2, Rs = 10−3,µb = 1, µs = 2e − 3. For high viscosity boluses, the no-slipshould give a much better approximation to the bolus flowthan a no-slip condition. . . . . . . . . . . . . . . . . . . . . . 149B.4 Contour plot of flow rates, as percentage of the full-slip flowrate, at 0.25s, as a function of both µs and Rs. . . . . . . . . 150B.5 Regions to the left of each contour represent choices of width,Rs, and viscosity, µS , for which a full-slip condition gives anapproximation of flow rate within 85% of the actual flow rate. 151xviiiAcknowledgementsFirst and foremost, I would like to thank Dr. Yoko Inamoto and Dr. EiichiSaitoh of Fujita Health University. Without the research they have con-ducted, as well as their willingness to collaborate, the bulk of this thesiswould have been impossible. I must thank Prof. Mark Nicosia for his guid-ance and wisdom in fluid simulation and swallowing. I would also like tothank Prof. Rebecca Affoo, who has spent an enormous amount of timegetting me up to speed on saliva. To Prof. Bill Pearson and Dr. JanaRieger, thank you for your invaluable support with anatomy and physiologyof swallowing, as well as the interpretation of the 320-row ADCT images.Dr. Nicole Rogus-Pulia, thank you for the support and interpretation of theoropharyngeal swallowing results. Nick Yeh and Dr. Dana Grecov, thankyou for your help with measuring the rheology of the bolus.Last and most importantly, I would like to thank my advisors for theirsupport and guidance. Thank you Sheldon, for all the chocolate, and Sid,for believing in me.xixI dedicate this thesis to my dad, who always believed in me, my mum, whopushed me to graduate and get a job, and Ronnie, who put up with methrough all these years.xxChapter 1IntroductionSwallowing disorders (dysphagia) occur for a multitude of reasons, includ-ing, but not limited to, muscle and tissue changes associated with aging,neurological disorder such as stroke or Parkinson’s disease, and oncologytreatment for tumours in the head and neck region. Dysphagia may leadto a significant decrease in quality-of-life and malnutrition. Aspiration offoreign materials into the lungs can lead to pneumonia and death. In aretrospective cross-sectional study, researchers found that over 50% of headand neck cancer patients treated with surgery combined with radiotherapyor chemotherapy suffered from oropharyngeal dysphagia, with malnutritionpresent in 20% of patients [37].One factor that is thought to contribute to dysphagia is reduced salivaproduction, or hyposalivation. For example, Pauloski et al. [101] found thathead and neck cancer subjects who underwent intensity-modulated radio-therapy, which was used to target the tumour while sparing the salivaryglands, performed better in swallowing tasks than a control group that un-derwent conventional radiotherapy. This suggests that sparing the salivaryglands improves swallowing function. However in a study by Rogus-Puliaet al. [107], the authors found that in subjects with Sjo¨gren’s syndrome,1Chapter 1. Introductionan autoimmune disorder that results in an extremely dry mouth, 96% ofswallows were judged to be functional, and no strong correlation was foundbetween temporal measures of swallowing and salivary flow rate. One of thechallenges in drawing strong conclusions from patient data is that often itis impossible to separate the variable of interest, for example salivary flowrate, from other confounding factors, such as reduced mobility of the tongue,hindering the analysis of the data.This thesis describes the extension of a state-of-the-art fluid simulationtechnique, Unified Semi-Analytic Wall (USAW) SPH [32], for the purposeof simulating oropharyngeal swallowing of non-Newtonian fluids. The sim-ulation method has been widely used for wave simulation, and possessesdesirable features such as an implicit representation of the air-liquid inter-face, and a natural handling of complex 3D geometry.320-row Area Detector Computed Tomography (ADCT) datasets give usan unprecedented view of the human head and neck during swallowing [57].The data also presents unique processing challenges when trying to extractgeometry that can be used in a computer simulation. A novel workflowcentred around a newly created tool, BlendSeg, was developed in order toextract the airway boundary correctly from the ADCT data.A computer simulation of oropharyngeal swallowing is performed us-ing the airways extracted from the ADCT data, and the extensions to theUSAW SPH method. These preliminary experiments of starch-thickenedliquid boluses in the oropharynx give us insight into the role of saliva innormal swallows. The results suggest that in healthy subjects saliva doesnot significantly increase the speed of starch-thickened boluses, contrary to21.1. Computer simulation of swallowingour initial expectations. However, saliva may help to decrease the quantityof residue after swallowing. Although these are preliminary experiments,they highlight the need for further study of this complex phenomenon.1.1 Computer simulation of swallowingComputer simulation of fluids can help us achieve a greater understandingof the physics in the human body because a simulation is able to isolateand change a single variable in order to understand its contribution to theoverall phenomena. It is also useful in situations where detailed measure-ments and observations are difficult or impossible to perform. For example,in-vivo measurements of the bolus flow in swallowing are usually limitedto 2D videofluoroscopic studies, fibreoptic endoscopic studies, or manome-try. Measuring quantities such as the saliva thickness on the surface of thetongue during swallowing, shear-rates and intra-bolus velocity is impossiblewith modern technology. Because of these difficulties, computer simulationpresents an attractive alternative that allows us to predict these quantitieswith a mathematical model.Swallowing is a complex phenomenon that most people take for granted.People swallow a large variety of materials that undergo significant struc-tural and chemical changes due to both mastication and mixture with saliva.As the aerodigestive tract is used for both breathing and swallowing, the in-terface between the bolus and the air must be considered. To further compli-cate matters, the tongue, soft-palate, pharynx, and larynx all undergo verylarge deformations during the oropharyngeal phase of swallowing in order31.2. Smoothed Particle Hydrodynamicsto propel the bolus into the esophagus while preventing aspiration into thelungs and nasal passage. Only a limited number of numerical swallowingsimulations in the oropharynx have been performed due to these challenges.This thesis limits its scope to the study of the fluid dynamics of liquid bo-luses in the oropharynx. It does not explicitly consider the structural andmuscular biomechanics of the oral, pharyngeal, and laryngeal complexes.1.2 Smoothed Particle HydrodynamicsSmoothed Particle Hydrodynamics (SPH) is a computational fluid simula-tion method that belongs to the “mesh-free” family of techniques [76]. Itexcels in simulating transient, free-surface flows because it does not requireexplicit handling of the fluid-air interface, making it particularly attractivefor an oropharyngeal swallowing simulation. One of the main challengeswith SPH is the application and enforcement of solid wall boundary con-ditions. Ghost particles, mirroring, and repulsive force particles are threeof the most common techniques for representing solid walls. Of the three,only the repulsive force particles are suitable for handling a complex mov-ing geometry. The repulsive particles line the surface of the solid walls andare effective at preventing the fluids from penetrating the solids. Howeverthey are not effective in enforcing the typical “no-slip” boundary conditionassumption in viscous fluid simulations.The Unified Semi-Analytic Wall (USAW) SPH method modifies the orig-inal SPH formulation to allow solid walls to be represented by arbitrarilyshaped surfaces [32]. The method can handle complex boundary shapes and41.2. Smoothed Particle Hydrodynamicsis capable of enforcing a “no-slip” condition. However, most of the develop-ment focuses on high Reynolds number flows of Newtonian materials suchas water. We extend the USAW SPH method to simulate non-Newtonianmaterials by deriving an equation for estimating shear-rate that includes anew surface integral term. Our implementation also presents a discretizedmethod of handling moving boundaries in 3D. This thesis describes the ex-tended USAW SPH method, and an implementation is tested using canoni-cal flows, such as Hagen-Poiseuille and the lid-driven cavity. We are able todemonstrate convergence of the method in the tested scenarios.These extensions are necessary for simulating oropharyngeal fluid flow,but are also applicable if USAW SPH is to be used in other applications.For example, in arterial blood flow, shear stress on the inner walls of arter-ies causes atherosclerosis, plaque formation on the inner walls of the bloodvessel, which may lead to a stroke or heart attack [23]. Numerical simu-lations of the blood flow can predict wall shear stress [13] and neglectingthe non-Newtonian properties of blood may significantly influence the sim-ulation results [59, 60]. Arterial walls are sometimes simulated with staticgeometry, however in reality they are elastic and numerical studies haveshown that the displacement of the walls has a significant effect on wallshear stress [82].51.3. Oropharyngeal bolus simulation with the extended USAW SPH method1.3 Oropharyngeal bolus simulation with theextended USAW SPH methodWe use the extended USAW SPH method to perform a numerical experi-ment using computer fluid simulation of liquid bolus flow in the oropharynxin order to study saliva’s role in swallowing. The experiments utilize geome-try derived from state-of-the-art dynamic, 320-row Area Detector ComputedTomography (ADCT) images of human swallowing. The ADCT creates dy-namic 3D images that sometimes contain artifacts due to subject or bolusmotion, and these images require significant human intervention to processcorrectly. A new workflow is described that allows for relatively efficientmanual extraction of moving boundaries from the ADCT images. The exe-cution of the workflow required the development of a purpose built softwaretool, BlendSeg, which allows an operator to deform a surface mesh whilemaintaining both the mesh and shape topology.A small series of experiments, limited in scope to a few human subjects,is performed by using the boundaries generated with BlendSeg to drive SPHparticles in both slip and no-slip simulations. The boluses swallowed by theparticipants in the ADCT images are starch-thickened fluids which exhibitnon-Newtonian behaviour. Their properties are measured in a rheometerand modelled as a Cross-type fluid using the extended USAW SPH method.We expected to find that the full-slip condition gives a good approximationto normal swallowing because of the presence of the saliva layer, which hasa much lower viscosity than the thickened boluses. However, the resultsindicate that the no-slip condition provides a much better approximation to61.4. Contributionsthe bolus flow, especially with respect to oral containment. On the otherhand, we observed an increase in oropharyngeal residue after the swallowwhen comparing the no-slip to the slip simulation. These findings suggestthat saliva does not have a large influence on the speed of the bolus, but itmight help to reduce oropharyngeal residue and therefore increase swallow-ing efficacy, consistent with findings of clinical researchers [106].The next section lists the contributions presented in this thesis. Chap-ter 2 describes the anatomy and physiology of swallowing, background influid simulation, bolus simulation, SPH, segmentation and modelling from4D imaging data, and the ArtiSynth biomechanics toolkit upon which thework is built. Chapter 3 describes the USAW SPH method and the ex-tensions required to simulate oropharyngeal bolus flow. We validate thesimulations by comparison with well-known Couette, Hagen-Poiseuille, andlid-driven cavity flows. Chapter 4 describes a computer fluid simulation oforopharyngeal swallowing using USAW SPH. The geometry motion is basedon dynamic 3D ADCT images of healthy participants swallowing starch-thickened boluses. The BlendSeg tool is also described in this chapter.1.4 ContributionsThe contributions of this thesis are:Extension of USAW SPH for oropharyngeal swallowing: The Uni-fied Semi-Analytic Wall (USAW) SPH method is a modern approach tocorrecting the “particle deficiency” problem inherent in classical SPH for-71.4. Contributionsmulations for fluid flow. We demonstrate that, with some modification,the method is suitable for fluid simulation of non-Newtonian boluses in theoropharynx. The method was extended to include a formulation for comput-ing shear-rates, which is in turn used to compute an instantaneous viscosityfor non-Newtonian fluids. We demonstrate numerical convergence using alid-driven cavity experiment and Hagen-Poiseuille flows.Workflow for creating dynamic 3D boundaries from 4D imagingdata: We describe a technique for creating dynamic, 3D moving bound-aries derived from 4D (3D and time) imaging data. Correctly extractinga suitable boundary from the images requires significant human interven-tion, and could not be fully automated. We created a tool, BlendSeg, thatgives a human operator a large degree of control in the creation of a movingboundary that closely follows dynamic 3D volumetric images.Numerical study examining saliva’s role in bolus transport andresidue in healthy individuals: These techniques are used to performa set of 3D oropharyngeal swallowing simulations. We test the performanceof both slip and no-slip boundary conditions in normal swallows of starch-thickened liquids. Boundary motion is derived from dynamic 3D imagesof healthy participants. The numerical predictions are compared to mea-surements gathered from the source images. The results suggest that innormal individuals saliva does not significantly increase the bolus speedfor thickened liquids. However, saliva may help to decrease the quantity oforopharyngeal residue, thereby improving swallowing efficacy. Further inves-81.4. Contributionstigation could expand the number of subjects and swallows under analysis,making the findings more clinically relevant.9Chapter 2BackgroundThis chapter begins with a brief overview of the anatomy and physiologyof human swallowing, followed by a review of literature in fluid simulationand Smoothed Particle Hydrodynamics. Previous research performing fluidsimulation of the bolus is then reviewed, followed by a review of segmentationof 4D data, and an introduction to the ArtiSynth simulation platform.2.1 Anatomy and physiology of swallowingSwallowing is a complex process consisting of four phases that are roughlydelineated: the preperatory and oral phases, which are voluntary, and thepharyngeal and esophageal, which are involuntary [25]. See Figure 2.1 for alabelling of anatomical structures in the head and neck.Preparatory and oral phaseDuring the oral preparatory phase, solids are broken down mechanically bycoordinated motion of the tongue, teeth, and cheeks. During this process,saliva is introduced to the food bolus, aiding in the break down of food andbeginning the digestion process. For liquids, the oral preparatory phase doesnot include mastication.102.1. Anatomy and physiology of swallowingFigure 2.1: An illustrated lateral view of a cut-away human head showingthe oral, pharyngeal and laryngeal region with major structures labeled [41].When the bolus is ready to be swallowed, it is placed in a trough formedon the superior aspect of the tongue. The tip of the tongue touches the backof the upper incisors. The soft palate forms a seal with the postero-superior112.1. Anatomy and physiology of swallowingaspect of the tongue that contains the bolus, preventing it from leaking pre-maturely into the pharynx. The swallowing sequence begins with the move-ment of the soft-palate postero-superiorly to close the nasopharynx. Almostsimultaneously, the tongue pushes up towards the palate, creating a travel-ling wave of pressure [96]. In healthy adults, the closure of the nasopharynxcoincincides with the beginning of the stripping wave that propels the bolusinto the esophagus [42]. As the tongue moves upwards, the prepared bolusmay fall into the pharynx due to gravity, or it may be pushed through dueto the motion of the tongue. The superior surface of the tongue continuespushing upwards, exerting pressure on the hard palate and squeezing thetail of the bolus out of the oral cavity and into the pharynx.Pharyngeal PhaseThe pharyngeal phase of swallowing is involuntary, and begins as the bolusenters the oropharynx. It is triggered by touch sensors on the pharynx. Asthe glossopalatal junction (GPJ) opens and the bolus passes through it, thepharynx is raised by the pharyngeal levator muscles to receive the bolus,while airway protection mechanisms are triggered in order to prevent aspi-ration of the bolus into the nasopharynx and the larynx. Nasopharyngealclosure is achieved by the soft palate and superior pharyngeal constrictors.Elevation of the tongue base is accompanied by elevation of the hyoid boneand lowering of the epiglottis. The protection of the airway is achieved pri-marily by closure of the vocal folds in the inferior to superior direction [78].As the bolus nears the larynx, the arytenoid cartilege begins to close andthe epiglottis lowers to a horizontal position. The epiglottis continues to122.1. Anatomy and physiology of swallowinglower and direct the bolus away from the vocal folds. Pharyngeal peristalsisis achieved through posterior motion of the base of the tongue and the ac-tivation of the middle and inferior pharyngeal constrictors. These motionssqueeze any remnants of the bolus down towards the esophagus. The pha-ryngeal phase ends when the bolus reaches the upper esophageal sphincter,which is open to receive the bolus.Esophageal PhaseLike the pharyngeal phase, the esophageal phase is also involuntary. Theesophagus contains rings of constrictor muscles along its length. The con-strictor muscles contract sequentially to form another peristaltic wave, trig-gered by touch sensors lining the esophageal wall, to push the bolus into thestomach. Peristalsis continues as long as there is food within the esophagus.2.1.1 Imaging of swallowingMedical imaging is a powerful tool in modern medicine for both clinical prac-tice and research. This section describes very briefly some of the most com-mon imaging methods used in swallowing. Videofluoroscopic (VFL) studiesand fibreoptic endoscopic evaluation of swallowing (FEES) are currently be-ing used in diagnostic imaging of dysphagia. In Japan, some researchers havebegun using dynamic computed tomography for diagnosis. Other imagingtechniques used for swallowing research include magnetic resonance imaging,and ultrasound.132.1. Anatomy and physiology of swallowingFigure 2.2: Videofluoroscopic image of subject swallowing barium. Theposterior of the bolus, contained in the oral cavity, has low image intensityand is visible on the left side of the image. Image by Hellerhoff, distributedunder a CC BY 3.0 license [50].142.1. Anatomy and physiology of swallowingVideofluoroscopic swallowing studyThe most common method for diagnostic imaging of swallowing is the vide-ofluorographic (VFL) study, also known as a Modified Barium Swallow(MBS) study. It is the current gold-standard method for clinical diagnosisof dysphagia. The VFL creates and records 2D projections of the subjectusing x-rays. There is good contrast between bones, soft tissues, and air(see Figure 2.2). Differentiating soft-tissues is challenging as they all appearwith similar intensities in the resulting images. Subjects swallow food andliquid mixed with barium, which provides contrast in the VFL images. Thework of Martin-Harris et al. [84] attempts to establish a standardized pro-tocol for VFL studies combined with a grading scheme, called the ModifiedBarium Swallow Impairment Profile (MBSImP).Fibreoptic Endoscopic Evaluation of SwallowingFibreoptic Endoscopic Evaluation of Swallowing (FEES) provides a videoimage of the pharynx, base of tongue, and larynx through the use of anendoscope inserted through the nasal passage (see Figure 2.3). AlthoughFEES is more invasive than VFL, it is more accessible to clinicians, and doesnot require the use of ionizing radiation. For swallowing, FEES has beenshown to be as sensitive as VFL in detecting premature spillage of the bolusfrom the oral cavity, residue after the swallow, and laryngeal penetrationand aspiration [53].152.1. Anatomy and physiology of swallowingFigure 2.3: Fibreoptic Endoscopic view of the larynx. The epiglottis, base oftongue, pharynx, and glottis are visible. Image by Med Chaos, distributedunder a CC BY-SA 3.0 license [16].Figure 2.4: 320-row area detector CT mid-sagittal slice of swallowing. 4 mlnectar-thick bolus is visible in the oral cavity with high image intensity.Image courtesy of Fujita Health University [57].162.1. Anatomy and physiology of swallowingFigure 2.5: Realtime-MRI of mid-sagittal slice of a normal subject swal-lowing pineapple juice. Image by Martin Uecker, distributed under a CCBY-SA 3.0 license [118].Dynamic Computed TomographyComputed Tomography (CT) uses ionizing radiation to create 3D volumetricimages. Recently, 320-row Area Detector Computed Tomography (ADCT)has been used to capture dynamic 3D images of swallowing at 10 framesper second [57]. The subject is seated in either a 45◦ supine or semi-proneposition during imaging, and swallows a bolus mixed with barium, similarto a VFL study (see Figure 2.4). Currently, ADCT is the only imagingmodality that captures fully 3D and dynamic images at rates sufficient fororopharyngeal swallowing.Real-time Magnetic Resonance ImagingMagnetic Resonance Imaging (MRI) uses a magnetic field to excite hydro-gen atoms in the body. In response to the magnetic field, the atoms emit172.1. Anatomy and physiology of swallowingFigure 2.6: Ultrasound image showing sagittal view of hyoid (low imageintensity) and thyroglossal duct cyst. Image by Laughlin Dawes, distributedunder a CC BY-NC 2.5 license [47].a signal that is received by the MRI and processed to create 3D volumetricimages. Real-time MRI has been used to image healthy subjects swallowingpineapple juice [118, 123] (see Figure 2.5). The sequence captures a singleslice at speeds of 24 frames per second. MRI does not involve the use ofionizing radiation, and is a very attractive alternative to VFL and CT. How-ever, it is not yet widely used, perhaps because of its relative inaccessibilitycompared to VFL and FEES. Some concerns are that it requires the subjectto lie in a supine position, which may be dangerous for individuals unableto maintain oral containment of a bolus.UltrasoundUltrasound (US) images have been used successfully in swallowing research,however is not currently used in diagnostic imaging (see Figure 2.6). Theyare an attractive alternative to VFL and FEES because it is non-invasiveand does not use ionizing radiation. US images have been used to image182.1. Anatomy and physiology of swallowingtongue propulsion of boluses, hyoid motion, pharyngeal motion, and theesophagus [17]. One interesting application of ultrasound is in imaging in-fants during breastfeeding [38], where it can be used to image a milk boluspassing from a breast or bottle, through the oral cavity, pharynx, and esoph-agus.2.1.2 SalivaSaliva plays an important role in oral health. Decreased salivary flow, or hy-posalivation, results in profound deterioration of oral homeostasis, increasedsusceptibility to dental caries and infections, decreased regulation and con-trol of the oral microflora [114]. Individuals with hyposalivation may alsoexperience impaired swallowing, or dysphagia [34, 100, 103].Salivary gland hypofunction can be caused by developmental or con-genital disorders [28], increased medication usage [113], systemic disorderssuch as Sjo¨grens syndrome [35], radiotherapy-induced damage to salivaryacinar tissue [51], anxiety [8], and aging [1]. Evidence examining the ef-fect of hyposalivation on swallowing, using an objective assessment methodsuch as videofluoroscopy, suggests that saliva is necessary for lubricationand moistening of the oral mucosa and bolus, and that hyposalivation mayaffect mastication, bolus formation and transport [43, 105]. A recent studyof chemoradiation patients finds a relationship between salivary flow rateand swallowing efficacy [106]. Despite the fact that some of these studiessuggest that there is a critical relationship between saliva and swallowing,other findings have failed to show a significant association [111]. For exam-ple, a study of patients suffering from Sjo¨gren’s syndrome [107], found no192.2. Fluid simulationstrong correlations between salivary flow rates and bolus transit times. Theexact nature of the relationship between saliva and swallowing is still notwell understood.2.2 Fluid simulationA fluid is any material in a state that continually deforms under shear stress,such as liquids, gasses, and plasmas. In fluid dynamics, the material is as-sumed to be a continuum even at infinitesimally small scales, ignoring in-dividual molecules and their interactions. Fluids are assumed to obey thelaws of conservation of mass, momentum, and energy. In this work, we as-sume that temperature remains constant throughout a simulation and doesnot have an impact on the momentum of the fluid, allowing us to ignorethe conservation of energy equation. Although all fluids are somewhat com-pressible, liquid flows at low speeds can be approximated as incompressible,meaning the density is a constant.When performing fluid simulation, one of the first choices that a modellerneeds to make is the method of representing the simulation domain. Theoptimal choice depends on the problem being studied. For flows with staticgeometry, for example air moving over an aeroplane wing or through a rigidpipe, static Eulerian discretizations are generally the best choice becausethey can be solved efficiently and accurately once a mesh is generated [27,64]. A mesh discretizes the simulation domain with a finite number of small“cells” or “elements” with non-zero volume. Quantities of interest, suchas the density, velocity, and pressure, are stored relative to the mesh. For202.2. Fluid simulationexample, pressure might be stored at the centre of each cell, while velocitiesmight be stored at the middle of a cell wall, or at a cell corner.An Eulerian discretization is one where the mesh that represents thedomain is not changed or moved for the duration of the simulation. Forexample, a network of weather stations on the ground that measures windspeed, pressure, and temperature tracks these quantities as a function oftime at fixed locations, forming an Eulerian representation of the geograph-ical area that they span. A Lagrangian discretization is one where each cellin the mesh always surrounds the same material. As the material moves,the cell walls moves as well so that the same volume of material remainswithin the cell from the start to the end of the simulation. Lagrangianmethods have limited use in fluid simulation because of the large amount ofdistortion that occurs for many fluid flows. As the mesh distortion increases,cells can become long and thin causing the volume of each cell to approachzero, or even become negative, causing the simulation to become unstable.At this point, the simulation stops and can only continue once the domainis remeshed, and quantities are transferred from the old to the new mesh.Generating a suitable mesh may be a non-trivial problem and, dependingon the type of elements used, may require significant manual intervention.Interpolation of variables from one mesh to another may result in undesir-able numerical diffusion, an artificial source of damping. If there is only amoderate amount of movement of solid geometry, the Arbitrary-Lagrangian-Eulerian (ALE) technique can accommodate the motion without the needto remesh, while maintaining most of the desirable properties of an Euleriansimulation [26].212.2. Fluid simulationFree-surface simulationSimulating a liquid bolus in the oropharynx is challenging because there areactually two fluids to consider, the bolus and the air, with an unknown in-terface between them, a so-called free-surface problem. Explicit simulationof the air phase is probably unnecessary for swallowing because of its lowdensity relative to most liquids. Meshed Lagrangian simulations are attrac-tive for free-surface problems because the cell walls of the mesh can be usedto represent the free surface. However this approach encounters issues if itreaches a point there is a large amount of movement in the fluid that causessignificant mesh distortion. For example, a thin streak of fluid is hard torepresent without using an intractable number of elements.Eulerian simulations have been modified to represent free-surfaces us-ing techniques such as Volume-of-Fluid (VOF) [52, 117], where an addi-tional scalar, between zero and one, representing the fraction of fluid oc-cupying each cell is tracked and advected in addition to the other physicalquantities. Level-set techniques [31, 108] work similarly to VOF, but thefree surface is represented implicitly as the zero level of a scalar field rep-resenting signed distance to the free-surface. Hybrid techniques such asParticle-In-Cell (PIC)/Marker-and-Cell (MAC) [46], Fluid-Implicit-Particle(FLIP) [11, 125] employ a background Eulerian grid and advect discreteparticles in order to track the fluid phase.222.2. Fluid simulationMoving boundary handlingAnother major difficulty of bolus simulation is the moving boundaries inthe aerodigestive tract which deform significantly in order to both propelthe bolus and protect the airway. This result in parts of the simulation do-main collapsing completely and then reopening, making meshing an extremechallenge.Eulerian methods can be modified to accommodate irregular and mov-ing boundaries using techniques such as Immersed Boundaries [88, 102]. Ofthese Immersed Boundary techniques, there are those designed for thin fila-ments [9, 124], which are two-way coupled and must include a model of thesolid mechanics, and those for irregular boundaries for which the solids aretreated as “rigid” (though they may change position) [29, 65, 90]. Hybridmethods such as PIC/FLIP have been reformulated to handle arbitrarilyaligned solids by using a variational formulation to derive a coupled pres-sure solve [5].2.2.1 Mesh-free methodsRecently, mesh-free methods have become popular for solving fluid free-surface problems because they can easily handle the fluid-air interface [6].Their most defining feature is that lack of the mesh-defined connectivitythat traditional Lagrangian and Eulerian methods rely on, resulting in abig advantage when there is a moving interface because mesh-free meth-ods do not require remeshing. Mesh-free methods are Lagrangian, as theircomputation points move with the material, meaning there is no need to232.2. Fluid simulationsolve for the advection term in the momentum equation, a non-linear termthat can be challenging to handle correctly. In general, however, they arecomputationally more expensive than the traditional meshed methods, andimposition of boundary conditions can be challenging. The increased com-putational cost is partially offset because mesh-free methods do not requirea mesh to be generated.The oldest and most popular mesh-free method for fluid simulation isSmoothed Particle Hydrodynamics (SPH), which was originally formulatedto solve problems in astronomy [40, 80, 91]. Many mesh-free methods havebeen developed for fluid simulation, such as the Moving Least Squares (MLS)method [7] and Moving Particle Semi-Implicit (MPS) [67, 68], but this workwill focus on the SPH method.2.2.2 Smoothed Particle HydrodynamicsSmoothed Particle Hydrodynamics (SPH) was first adapted for incompress-ible free-surface fluid simulation by Monaghan et al. [92]. Typical SPHderivations begin by considering that a smooth function in a domain χ maybe reproduced exactly under convolution with the Dirac delta function:f(x) =˚χf(x′) δ(x′ − x) dx′ (2.1)The representation is then transformed to a smoothed approximation byreplacing the delta function with a “kernel” function, usually a polynomial242.2. Fluid simulationapproximation to the Gauss function, with local support:〈f(x)〉 =˚χf(x′)W (x′ − x, h) dx′, (2.2)where W (x′ − x, h) is the kernel function, centered at 0 and scaled by h,a “smoothing length” which determines the extent of the interaction. Ash approaches zero, W (x, h) tends towards the Dirac delta function. Thekernel is defined to have local support, so that W (x′ − x, h) = 0 for anythingoutside of its support domain, Ω. This allows the integral in equation (2.2)to be evaluated over the truncated domain Ω:〈f(x)〉 =˚Ωf(x′)W (x′ − x, h) dΩ. (2.3)Monaghan’s first paper [92] applies the SPH equations to the inviscid,incompressible Navier-Stokes equations (the Euler equations). It uses weakcompressibility and an artificial speed-of-sound to solve for pressure. An ar-tificial damping term is included to damp pressure oscillations and maintainstability. Special repulsive-force particles are placed along boundaries toprevent SPH particles from passing through solid boundaries. A Newtonianviscosity term was first described by Morris et al. [95], who also describea surface-tension term [94]. The viscosity term was later validated in 3Dby Sigalotti et al. [110]. It has been used for diverse applications such aslarge wave modelling in dam breaks [71], mud-flows [109], and high-pressuredie-casting [18].Monaghan’s damping/stabilization term is capable of adding a signifi-252.2. Fluid simulationcant amount of viscous damping (energy loss) to a simulation. The workof Colagrossi et al. [19] use a moving least squares (MLS) filter to smoothdensity variations in the fluid, and also show that it is applicable for explicittwo-phase simulations with high density ratios between the liquid and gas.The MLS filter conserves energy and also reduces the amount of energy lossin Monaghan’s damping term. However, this filter was found to disrupt thetotal volume of fluid for long simulations by Antuono et al [3]. Marrone etal. [83] describe a numerical viscous damping term that, when applied toweakly compressible solvers, smooths the pressure field significantly with-out adding unwanted viscosity to the flow while also preventing numericalinstability.Non-Newtonian fluids in SPHGeneralized non-Newtonian fluid models are idealized representations offlows where the apparent viscosity, represented in this thesis by µeff , de-pends only on the shear rate, γ˙. Some of the common viscosity modelsinclude power-law, Herschel-Bulkley, Bingham, and Cross. Non-Newtonianfluids have been modelled with SPH, mostly explicitly. Shao and Lo [109]simulated a Cross fluid to model mud flow. Xu et al. [122] perform a 3Dsimulation of a Cross fluid to model droplet impacts, thin-plate injectionmold, and jet buckling. Hosseini et al. [56] simulate power law, Herschel-Bulkley, and Bingham viscosity models using an explicit solver. Explicitviscosity solvers place a severe time-step restriction on the simulation. Themaximum time step size typically depends on the order of ∆x2, while anexplicit pressure solve limits time steps to the order of ∆x. To make mat-262.2. Fluid simulationters worse, the viscosity time step restriction scales with the inverse of theviscosity. For a power-law fluid at rest, the effective viscosity is infinity. Fanet al. [30] describe an implicit formulation in 2D for the viscosity term inthe momentum equation which does not limit the time step size, howeverit requires solving a large system of equations. For high effective viscosi-ties, this can result in much larger time step sizes and faster simulations,although the cost of solving the each time step is high.Incompressible SPH formulationsThe first truly incompressible SPH method (PSPH) was described by Cum-mins and Rudman [22]. They use an approach analogous to a “fractional-step” method in meshed techniques, and describe an approximate Laplacianoperator, similar to Morris et al.’s viscous diffusion term, to perform an im-plicit pressure projection. They take a fractional step to advance velocityand then use the divergence of velocity as a source term to solve for thepressure. Boundary handling was performed using reflected particles, withspecial handling for 90◦ corners. Contrary to explicit SPH methods, the im-plicit pressure solve requires an explicit definition of the free-surface bound-ary so that the zero boundary condition can be applied correctly. Shao andLo [109] use a symmetric version of Cummins and Rudman’s Laplacian oper-ator but use a source term derived from a fractional computation of densityinstead of divergence of velocity. Their boundary handling used fixed ghostparticles instead of fictitious reflected particles. Both of these methods weredeveloped for 2D, and with walls containing 90◦ degree bends only. It is notobvious how these methods can be extended in general to arbitrary boundary272.2. Fluid simulationshapes. Incompressible solvers generally show much smoother pressure con-tours than weakly-compressible solvers [71], however setting up and solvingthe pressure projection equation can be time consuming, though this maybe offset by the ability to take larger time steps while still remaining stable.Additionally, care must be taken when defining the free-surface interfaceas this forms a boundary condition in the pressure solve and as such, thepressures are sensitive to this definition.Insufficient kernel support and the γ correction termNear solid walls, two problems arise with SPH: applying boundary condi-tions, and handling the truncated kernel support. The truncated kernelsupport arises near walls and causes the integral in equation 2.3 to underrepresent the true value of f(x). Bonet and Kulasagaram re-formulate SPHas an Euler-Lagrange problem and the resulting formulation, in the absenceof viscosity, conserves linear and angular momentum. During this process,they introduce a corrective term, γ, into the SPH formulation that correctsfor the insufficient kernel support problem [10, 69],γ(x) =˚ΩW (x′ − x, h) dΩ, (2.4)which evaluates to unity within the fluid and less than unity when close toa solid wall or free surface. The modified SPH interpolation equation is〈f(x)〉 = 1γ(x)˚Ωf(x′)W (x′ − x, h) dΩ. (2.5)282.3. Fluid simulations of swallowingOne of the main challenges with this approach is that γ and its gradient aredifficult to compute accurately. They address this problem by pre-computingγ and then approximating its value with a spline function based on distanceto walls. Special handling is required for corners and irregular geometries.The Unified Semi-Analytic Wall (USAW) SPH method extends the workof Bonet and Kulasegaram. Ferrand et al. describe an efficient and accuratemethod for computing γ by describing an analytic solution for ∇γ whenusing the 5th-order Wendland kernel [32], and solve for γ by integrating itin time along with other physical variables. Leroy et al. extend USAW SPHto solve for pressure implicitly [72], as well as to include a k −  turbulencemodel. The USAW methods mentioned so far have only been applicable for2D solvers. Mayrhofer et al. extend the method to 3D [85]. One of thelimitations of the USAW method is that the analytic calulation of ∇γ hasonly been formulated for the Wendland kernel. Mayrhofer et al. describea theoretical approach for accomodating moving boundaries in 3D, howeverthey do not describe a discretized approach and application. Additionally,they have not yet adapted their formulation for non-Newtonian materials.2.3 Fluid simulations of swallowingComputer simulations of swallowing in the oropharynx has been performedby other researchers in the past, however many of the models were cre-ated with simplified geometries of the airway. For example, Chang et al.,in [15], model a Newtonian fluid in a axisymmetrized pharynx, from theglossopalatal junction to upper esopahagel sphincter (UES). They use a292.3. Fluid simulations of swallowingfinite-element solver on a mesh with prescribed boundary motion, and as-sume there is no air. The boundary motion is interpolated from the work ofKahrilas et al. [62] and Cook et al. [21], who use bi-plane videofluoroscopyto determine opening diameters of the oropharynx and UES, respectively.Inlet conditions are prescribed by applying a pressure for the duration ofGPJ opening, and a velocity estimated from the literature. A pressure ofzero is prescribed at the outlet when the UES is open. The authors suggestthat this work could be used to predict pressure profiles in the pharynx.This work is extended to 3D in [14] using a K-SNAKE algorithm appliedto bi-plane fluoroscopy images, and assuming that pharyngeal cross sectionshave an elliptical shape. Meng et al. [87] extend the 2D work in [15] bymodelling a non-Newtonian power-law fluid instead of a Newtonian fluid.3D visualization of oropharyngeal swallowing from medical images werecreated in [61] and [73], as an aggregate of multiple swallows. Their workdid not simulate the bolus, but used image data to visualize it.Nicosia et al. [99] approximate bolus ejection from the oral cavity bysimulating two approaching parallel plates in 2D axisymmetric coordinatesby solving for the stream and vorticity functions. They used this model tostudy the pressure required to clear half of the bolus from the oral cavity asviscosity varies. The work was extended in [98] to study the effect of partial-slip boundaries on the shear rates and ejection times of a bolus. In anotherstudy by Nicosia [97], an ALE simulation was solved with finite-elements tostudy oral containment of a bolus on a tongue subject to “lingual pumping”gestures exhibited in some subjects prior to swallowing. The simulation wasin 2D and used geometry derived from videofluoroscopic images [62].302.4. Segmentation and modelling from 4D imaging dataMizunuma et al. [89] used a static CT scan as the basis for creating arealistic 3D geometry of the oral cavity and pharynx. Movement of solidswas derived from videofluoroscopic data The solids simulation used finite-elements, and a solid visco-elastic finite element model represented the bolus.Sonomura et al. [112] used similar geometry but replaced the volumetricFE organs with shell elements, and performed a fluid simulation of thebolus using volume-of-fluids to represent both water and a starch-thickenednon-Newtonian bolus. Kikuchi et al. [63] use a mesh-free MPS method tosimulate the organs and bolus, using structures derived from CT data andmotion derived from videofluoroscopic images.Outside of swallowing itself, mastication has been simulated using SPHby Harrison et al. [48, 49]. They use an elastobrittle SPH material model torepresent agar gel and its fracture during mastication.2.4 Segmentation and modelling from 4Dimaging dataRecent advances in dynamic Area Detector CT (ADCT) now allow for 3Dand time imaging of dynamic phenomena such as swallowing. The currentstate of the art uses a 320-row ADCT to capture dynamic 3D images [36, 57].The amount of detail in these images is stunning, and allows us to studydynamic oropharyngeal swallowing in 3D. A fluid simulation of the liquidbolus in the oropharynx would allow us to study physical phenomena relatedto swallowing, such as xerostomia (dry-mouth), and how it affects bolus flow.Commercial software, such as Amira or Mimics, are not currently ca-312.4. Segmentation and modelling from 4D imaging datapable of extracting a dynamic surface from 4D CT data that is suitablefor simulation. Much of the published literature on 4D segmentation hasbeen focused on the heart, and use model based approaches to identifyingstructures. McInerney et al. [86] use a model based approach that utilizesan expanding sphere in order to segment the left ventricle and track its mo-tion in 4D. A Finite Element model provides regularizing forces based onstretching and curvature. Montagnat et al. [93] define a 4D simplex-basedframework in order to segment the left ventricle. Their framework allowsfor curvature-based regularization, or a statistical-model-based regulariza-tion in both space and time. Some approaches use statistical models thatmake use of prior knowledge of the shape of the heart and its motion [58, 81],or take advantage of the cyclic nature of the heart beat and use a Fourieranalysis [24, 74].Because the ability to generate 4D images of the airway is so new, to thebest of the author’s knowledge there has been no work done yet in modellingthe airway from 4D CT images. Statistical models can not be leveraged untila sufficiently sized training set has been created. 4D images of oropharyngealswallowing contain challenges that are not present in images of, for example,the left ventricle. The airway is the void space defined by surrounding organssuch as the tongue, hard- and soft-palate, pharynx, and the esophagus. Theesophagus is initially closed and pressed flat, typically providing no contrastunless an air bubble is trapped or the bolus is passing through it. It runsparallel and very near to the trachea, which has very strong edges betweenthe airway and the air. During the propulsion phase of swallowing, the oralcavity and pharynx collapse completely, also resulting in very little edge322.5. ArtiSynth biomechanics toolkitcontrast over large parts of the region. The vocal folds and nasopharynx arealso squeezed shut in order to protect the airway. To add further difficulty,the bolus, when stationary, provides ample contrast with surrounding softtissues, however once it begins to move, its contrast lowers to the pointwhere it has similar intensity to the soft tissues, making edge identificationmore difficult. Blurring and ghosting artifacts in the images further edgedetection methods. Volumetric model based methods that take advantage ofphysical constraints such as incompressibility can not be used to regularizethe segmentation unless all of the organs are registered individually.2.5 ArtiSynth biomechanics toolkitArtiSynth is a Java-based toolkit for solving multibody physics problems [77].It includes a physics library that can perform coupled rigid-body and finite-element simulations, along with complex tasks such as the ability to resolvecollisions and solve inverse simulation problems. ArtiSynth’s main applica-tions are for simulation of biomechanics problems in the human body. Ithas been used to simulate many anatomical structures and functions suchas speech [45], tongue motion [115], face [33], foot [79], and mandibularmotion [2]. The SPH model and oropharyngeal swallowing simulations de-scribed in this thesis are implemented in ArtiSynth. They utilize its graphicsrendering capabilities, user interface, vector math library, hierarchical treerepresentations, and proximity detection methods.33Chapter 3Smoothed ParticleHydrodynamics for fluid flowSmoothed Particle Hydrodynamics (SPH) is now a popular method for sim-ulating free-surface flows of liquids. This chapter describes extensions to astate-of-the-art SPH method, the ability to simulate non-Newtonian materi-als in 3D moving walls, that are necessary for a swallowing simulation. SPHis a Lagrangian method, meaning simulation points remain “fixed” to a ma-terial point. The simulation points are referred to as “particles”. Movementof a point represents movement of a quanta of material. SPH is a mesh-freemethod, meaning there is no formally defined connectivity between a parti-cle and its neighbours. Instead, each particle can be thought of as exertingan influence on its neighbours, where the strength of this interaction scaleswith a “kernel” function.SPH is most commonly used to represent free-surface flows of high-Reynolds problems, such as dam breaking [83] and sloshing tanks [3]. Webelieve that with some modification, SPH can be adapted for use in anoropharyngeal swallowing simulation. Its advantages over traditional meth-ods are that it implicitly handles free surfaces and splitting and merging of343.1. Classic Smoothed Particle Hydrodynamics formulationfluids without any special handling. This chapter describes the basic for-mulation of SPH, and the Unified Semi-Analytic Wall (USAW) approach todealing with solid boundaries. One of the primary contributions of this the-sis are extensions to the USAW SPH method to be able to simulate movingwalls in 3D, and non-Newtonian fluids. They are described in sections 3.2.2and 3.2.3, resp. At the end of the chapter, Couette flow, startup Hagen-Poiseuille flow, a lid-driven cavity problem, and non-Newtonian Hagen-Poiseuille flow are simulated with our extended SPH method. The resultsare compared to theoretical solutions, or in the case of the lid-driven cavity, acanonical solution to the problem. The extended USAW SPH method showsgood agreement with accepted solutions, and approximately first-order con-vergence.3.1 Classic Smoothed Particle HydrodynamicsformulationTypical SPH derivations begin by considering that a smooth field in a do-main χ may be reproduced exactly under convolution with the Dirac deltafunction [75, 76]:f(x) =˚χf(x′) δ(x′ − x) dx′ (3.1)The representation is then transformed to a smoothed approximation byreplacing the delta function with a “kernel” function, usually a polynomial353.1. Classic Smoothed Particle Hydrodynamics formulationapproximation to the Gauss function.〈f(x)〉 =˚χf(x′)W (x′ − x, h) dx′ (3.2)where W (x′ − x, h) is the kernel function, centered at 0 and scaled by h, a“smoothing length” which determines the extent of the interaction.The kernel is defined to have local support, so that W (x′ − x, h) = 0for anything outside of its support domain, Ω. This allows the integral inequation (3.2) to be evaluated over the truncated domain Ω:〈f(x)〉 =˚Ωf(x′)W (x′ − x, h) dΩ (3.3)The kernel should also be normalized to unity,ˆΩW (x′ − x, h)dx′ = 1. (3.4)Spatial discretization is performed by seeding the domain with SPH par-ticles, each having mass and density. The continuous approximation (3.3)is then discretized:(f(x)) =N∑bf(xb)W (xb − x, h) ∆Vb (3.5)where the summation is over the N particles near x for which the kernelfunction W is not equal to 0, and xb is the position of neighbour particleb. See Figure 3.1 for a 2-D example. ∆Vb is the volume represented by theneighbour particle, and is defined as the mass divided by density ∆Vb =mbρb.363.1. Classic Smoothed Particle Hydrodynamics formulation∂ΩΩxW (x)Figure 3.1: Example 2-D particle distribution for an SPH fluid particle(black circle) with full support over its domain Ω (dotted area). The do-main boundary is designated by ∂Ω (dashed circle). A typical 1-dimensionalkernel function W (x) is overlaid. Notice that W (x) = 0 for x ∈ ∂Ω.3.1.1 Spatial derivativesSPH offers a simple method for estimating spatial derivatives. For a diver-gence,〈∇ · f(x)〉 =˚Ω[∇ · f(x′)]W (x′ − x, h) dΩ (3.6)and after integrating by parts,〈∇ · f(x)〉 =˚Ω∇ · [f(x′)W (x′ − x, h)]dΩ−˚Ωf(x′)[∇W (x′ − x, h)] dΩ(3.7)Application of Gauss’s (divergence) theorem to the first term of the373.1. Classic Smoothed Particle Hydrodynamics formulationright-hand side of this equation yields〈∇ · f(x)〉 =‹∂Ω[f(x′)W (x′ − x, h)]dS−˚Ωf(x′)[∇W (x′ − x, h)] dΩ(3.8)where S is the area element in the direction of the face normal.SinceW (x′ − x, h) is defined to be 0 (see Figure 3.1) on the boundary ∂Ω,the first term is usually dropped, leading to the final gradient approximation:〈∇ · f(x)〉 = −˚Ωf(x′)[∇W (x′ − x, h)] dΩ. (3.9)The derivative is now applied to the kernel function, defined analytically, andso it is very simple to evaluate. This equation is discretized and summed,as in equation 3.5.The above assumption holds as long as the fluid integration domain iscovered adequately with particles. However, if the fluid domain is truncatedby a solid wall, this assumption no longer holds. This is one of the sourcesof the “particle deficiency” problem.In practice, the following equivalent forms of the spatial derivatives arediscretized since they demonstrate improved stability.∇ · f(x) = 1ρ[∇ · (ρf(x))− f(x) · ∇ρ] , (3.10)and,∇ · f(x) = ρ[∇ ·(f(x)ρ)+f(x)ρ2· ∇ρ]. (3.11)After applying the SPH discretization method, the symmetrized forms383.1. Classic Smoothed Particle Hydrodynamics formulationof the derivative appear:(∇ · f(x))a =1ρaN∑bmb (f(xb)− f(xa)) · ∇W, (3.12)and,(∇ · f(x))a = ρaN∑bmb(f(xb)ρ2b+f(xa)ρ2a)· ∇W, (3.13)3.1.2 SPH for the incompressible, viscous, Navier-StokesequationsThe typical SPH approach to solving the incompressible Navier-Stokes wasfirst presented by Monaghan [92]. Incompressibility is weakly enforcedthrough a stiff equation of state utilizing an artificial speed of sound. Forthe viscosity term, we use the method described by Morris [95].For a Lagrangian system, the incompressible Navier-Stokes equationsmay be written asDvDt= −1ρ∇p+ 1ρ(∇ · µ∇) v + B (3.14)DρDt+ ρ∇ · v = 0 (3.15)where v is the velocity, DDt is the total derivative, ρ is the density, p ispressure, µ is the Newtonian viscosity, and B is the body force (e.g. gravity)per unit volume.We use the following equation of state to relate the density changes to393.1. Classic Smoothed Particle Hydrodynamics formulationthe pressure [94]:p = c2 (ρ− ρ0) , (3.16)where ρ is the density of the current particle, ρ0 is the reference densityof that particle, and c is an artificial speed of sound in the fluid, chosenlarge enough so that incompressibility is approximately obtained, but smallenough so as not to make the time step size too restrictive.Applying the SPH discretization operator (3.12) to the continuity equa-tion (3.15) yields the formula for updating each particle’s density:DρaDt= −N∑bmbvab · ∇W (3.17)where we use vab as short-hand for va − vb. The density of each particle isthen updated explicitly.For the momentum equation, the contribution from the pressure term iscomputed using equation (3.13):(1ρ∇p)a=N∑bmb(paρ2a+pbρ2b)∇W. (3.18)The diffusion term in the momentum equation is evaluated using the methoddescribed by Morris [95]:(1ρ(∇ · µ∇)v)a=N∑bmb(µa + µb)vabρaρb(1rab∂W∂ra)(3.19)where rab is the distance between particles a and b.403.1. Classic Smoothed Particle Hydrodynamics formulation3.1.3 Boundary handling in SPHWhen fluids interface with the air, the air pressure is assumed to be zero,and this is handled implicitly by the pressure formulation by the lack ofparticles representing an air phase [76]. Similarly, shearing forces from theair are assumed to be negligible. These assumptions are reasonable for manyapplications.Interfaces with solid walls are much more difficult to handle properly,and their representation in SPH is one of the primary differentiators betweendifferent methods. The traditional methods for dealing with solid walls arethose of repulsive forces [92], dummy particles [95], and mirror particles [22].Dummy particles require special placement around irregular boundaries inorder to maintain correct densities. Similarly, mirroring requires specialhandling of corners to maintain a correct particle density. For irregularmoving boundaries such as those found in the human body these methodswould be unable to correctly represent the solid wall. In this chapter we userepulsive force particles that prevent fluid particles from penetrating intosolid walls. They are easily placed along an irregular boundary have noissues with moving geometry. However they are unable to enforce a no-slipboundary condition correctly.413.1. Classic Smoothed Particle Hydrodynamics formulation3.1.4 Smoothing pressures with δ-SPHIn order to stabilize an explicit time integration scheme, Monaghan [92] usesan artificial viscosity term that damps vibrations.Πab =−αcµabρab, uab · rab < 00, otherwise(3.20)where µab = (µa+µb)/2, ρab = (ρa+ρb)/2, c is the speed of sound, and α isa parameter that can be tuned based on the simulation. This term is addedwithin the pressure summation term of the momentum update equation.In our experiments we found that this term adds some undesirable viscousdamping in addition to numerical damping.The work of [3] proposes a new damping term that is purely numerical.The method is given the name of δ-SPH [83]. In their work they use fixedboundary particles to represent boundaries. We find that it works well,without modification, with the USAW SPH scheme and produces relativelysmooth pressure distributions even with an explicit solver and a high speedof sound. In our experiments it stabilizes the simulation without addinga significant amount viscous damping, and therefore replaces Monaghan’sdamping term entirely.The δ-SPH method adds additional damping terms to the SPH updateequations,Dρa = ξhcN∑bmbρbψab · ∇W, (3.21)423.2. Unified Semi-Analytic Wall SPHDpa = αhcρaN∑bmbρbpiab∇W, (3.22)where the terms Dρa and Dpa are added to the density and momentum updatesfor particle a, respectively, c is the speed of sound, and h is the kernelsmoothing length. ξ and α are constants that that control the amount ofdamping. In our experiments, both are set to 0.02, and with some informaltesting we found that simulation results are not sensitive to this choice. ψaband piab are defined asψab = 2(ρb − ρa) rb − ra|rb − ra|2 , (3.23)piab =(ub − ua) · (rb − ra)|rb − ra|2 , (3.24)3.2 Unified Semi-Analytic Wall SPHThis section briefly describes the evolution from the classic SPH formula-tion to the Unified Semi-Analytic Wall (USAW) method of boundary han-dling [32, 72, 85]. We demonstrate how 3D moving boundaries, such as thosefound in the airway, can be simulated using the USAW SPH method. Next,we present a formulation for non-Newtonian viscosities using USAW SPH.433.2. Unified Semi-Analytic Wall SPH∂ΩΩW (x)W (x) 6= 0,∀x ∈ ∂ΩFigure 3.2: Example particle distribution for an SPH fluid particle (blackcircle) with truncated support over its domain Ω (dotted area). The domainboundary is designated by ∂Ω (dashed circle). A typical 1-dimensional ker-nel function W (x) is overlaid. Notice that unlike Figure 3.1, W (x) 6= 0 for∀x ∈ ∂Ω.3.2.1 Solid boundaries and insufficient domainrepresentationThe truncated support domain causes two problems in the SPH approxima-tion, the first is that the kernel is no longer normalized, and the second isthat the value of the kernel on the boundary is no longer zero, and so theboundary terms can not be dropped. Bonet and Kulasegaram address thefirst issue in [10, 69] by beginning from a variational formulation for a parti-cle system and introducing a scaling term to 3.2 the interpolation equationby γ, a renormalization factor.〈f(x)〉 = 1γ˚χf(x′)W (x′ − x, h) dx′, (3.25)443.2. Unified Semi-Analytic Wall SPHγ =ˆΩW (x′ − x, h)dx′ (3.26)However the computation of γ is non-trivial and an accurate method ofevaluation is required. Their variational formulation leads to a new term inthe pressure equation that accounts for boundary forces, which includes aderivative of γ term.Ferrand et al. [32] reformulate the boundary force term (∇γ) and theirgroup further expand the method to include turbulence, a 2D truly incom-pressible pressure solve [72], and into 3D [85]. Much of the challenge is infinding an accurate and efficient method for evaluating γ and ∇γ.The density and momentum updates, equations (3.17 - 3.19), have hadtheir boundary terms dropped due to the assumption that the kernel-functionW (x′−x, h) is nil on the boundaries. This assumption is true for fluid par-ticles that are far from a wall, but does not hold near solid boundaries (seeFigure 3.2). The equations (3.12 - 3.13) are re-written here with the surfaceterms and renormalization terms included:(∇ · f(x))a =1γaρaN∑bmb (f(xb)− f(xa)) · ∇W− 1γaρaK∑kρk (f(xk)− f(xa)) · nW∆S (3.27)453.2. Unified Semi-Analytic Wall SPH(∇ · f(x))a =ρaγaN∑bmb(f(xb)ρ2b+f(xa)ρ2a)· ∇W− ρaγaK∑k(f(xk)ρk+ ρkf(xa)ρ2a)· nW∆S (3.28)where the second term on the right-hand side of each equation is the bound-ary term, K is the set of boundary elements that intersect the sphere ofinfluence, Ω, of particle a, xk is the centroid of the surface element k, nis the unit-normal perpendicular to the surface element k, and ∆S is themagnitude of the area of the surface element.These modified boundary terms are applied to equations (3.17) and (3.18),the density update equation and pressure term of the momentum equation,respectively.DρaDt= − 1γaN∑bmbvab · ∇W + 1γaK∑kρkvak · nW∆S (3.29)(1ρ∇p)a=1γaN∑bmb(paρ2a+pbρ2b)∇W − 1γaK∑k(ρkpaρ2a+pkρk)· nW∆S(3.30)The same procedure can be applied to the diffusion term (3.19) and thefollowing equation is obtained:(1ρ(∇ · µ∇)v)a=1γaN∑bmb(µa + µb)vabρaρb(1rab∂W∂ra)− 1γaρaK∑k(µa + µk)W (∇v · n)k ∆S. (3.31)463.2. Unified Semi-Analytic Wall SPHEvaluating equations (3.29 - 3.31) requires for the velocity gradient dot-ted with the surface normal (∇v · n), for each surface element k. The nextsection discusses how this estimate is computed.Estimating velocity gradient at the wallThe computation of the velocity gradient is not trivial due to the fact thatparticle distributions may be sparse and irregular. Consider the case whereonly one particle is near enough to a boundary element to influence it;there is inadequate information to compute the full tensor ∇v. However,note that the full gradient is not required in equation (3.31), and that itonly contributes to the viscous term after being multiplied by the surfacenormal, ∇v · n.In the viscous formulation by Morris for fluid particles, the velocity gra-dient is only evaluated in the relative direction between the two particles.A finite-difference approximation gives a good estimate for this quantity.After transforming the integral using the divergence theorem, the term thatappears is ∇v · n. Using the finite-difference approach for the boundaryterms would work for flat and concave boundaries, such as plane Couetteflow and Hagen-Poiseuille flow. However it would be unstable for bound-aries with convex curvatures, such as over a sphere or a bluff body. Thisinstability arises because there is not enough information to give a goodestimate for ∇v · n when the surface normal and relative position vectorsare perpendicular or near-perpendicular (see Figure 3.3).To overcome this difficulty, it is assumed that for any fluid particle a,∇v · n is constant over the surface area under integration. It is defined by473.2. Unified Semi-Analytic Wall SPHyxdSn∆y = 0Figure 3.3: This figure illustrates the difficulty in getting an estimate for∇v·n for the differential surface element dS, using a fluid point (filled circle)perpendicular to the surface normal, n. This situation occurs for geometrieswhere the fluid is contained in a non-convex boundary. There is inadequateinformation to estimate any shear components since ∆y between dS and thefluid particle is nil.using the point on the solid surface which is closest to a, denoted rg, and byassuming the surface normal at rg is in the direction of rag (see Figure 3.4).The following formula is then used:(∇v · n)a =[1ragvag ⊗ ng]· ng (3.32)ng =ragrag. (3.33)Equation (3.32) simplifies to:(∇v · n)a =1ragvag, (3.34)which is simply the finite difference approximation of the velocity gradient inthe direction of particle a from point g using reflective boundary conditions.483.2. Unified Semi-Analytic Wall SPH∂ΩyxangFigure 3.4: The vector ∇v · n is assumed to be constant over ∂Ω for a fluidparticle a where it is truncated by the wall. To compute ∇v · n, the pointon the wall nearest to a, denoted g and indicated with a grey circle, is usedin the computations. Additionally, it is assumed that the normal n pointsfrom g to a. This assumption is good near solid walls with low curvature,but gives the appearance of a “smoothed” wall near sharp, convex corners.This approximation results in the following viscous term being used:(1ρ(∇ · µ∇)v)a=1γaN∑bmb(µa + µb)vabρaρb(1rab∂W∂ra)− 1γaρa(1ragvag)aK∑k(µa + µk)W∆S (3.35)This scheme has the added benefit of being computationally efficient, since(∇v · n)a is only evaluated once for each particle in the vicinity of a solidwall, as opposed to being evaluated for each differential surface element inthe vicinity of the particle.Surface discretizationWhen using a triangular mesh to represent a solid boundary, the trianglesmay not be of a similar size to the inter-particle spacing. Triangles that493.2. Unified Semi-Analytic Wall SPHare much smaller than the particle spacing do not pose an issue except thatevaluating surface integrals may be inefficient. If triangles are much largerthan the particle spacing, a pre-processing step can be performed checkingeach triangular element for size. If the length of the longest edge is largerthan the initial particle spacing, the triangle is subdivided into four smallertriangles. This process is repeated recursively until none of the elementshas any edge longer than the initial particle spacing. The centre of eachsub-element is used in distance computations to fluid particles; that is, atriangle’s area may be excluded from surface integration if the barycentriccentre of the triangle is too far from a particle, even if some part of thetriangle does lie within the fluid particles’ sphere of influence Ω. Conversely,a triangle’s full area may be included in the computation even though someof it may fall outside of Ω. This results in the method having some sensitivityto the quality of the surface discretization.Care must also be taken when defining the initial distribution of fluidparticles. Each SPH fluid particle can be considered as representing a volumeof fluid, with volume dV = m/ρ. One edge of this volume of fluid is typicallythe plane halfway between this particle and its nearest immediate neighbour.With this in mind, particles adjacent to a wall should initially be placed awayfrom the wall at a distance of half the initial spacing between particles. Ifone is not careful, an incorrect volume of fluid may be created. Placingparticles on the surface of the mesh can alleviate this issue, however forirregular meshes that deform, this method may not be suitable.503.2. Unified Semi-Analytic Wall SPH3.2.2 Moving boundaries in 3DThe work of Mayrhofer et al. [85] extends the USAW SPH method to threedimensions. The main challenge is the evaluation of the∇γ term in 3D. Theydemonstrate an analytical method for determining ∇γ, however it requiresthe use of the 5th order Wendland kernel [120]. In order to evaluate γ ata time t, γ can be initialized at time zero and then integrated numericallyas the particle moves through the ∇γ field. Beginning with the followingidentity, which can be derived from the chain rule,tn+1ˆtndγdtdt =tn+1ˆtn∇γ · u dt, (3.36)γ(tn+1)− γ(tn) =xn+1ˆxn∇γ · dx, (3.37)where γ follows the material point. This formulation assumes walls arestatic, e.g. for a static particle with a wall approaching it, dx is zero howeverγ(tn) 6= γ(tn+1). Mayrhofer et al. briefly describe a possible approach tothe solution of γ when walls are moving, and here we present our discretizedimplementation,xn+1/2 = xn +12(ua(tn)− um(tn))dt (3.38)γa(tn+1)− γa(tn) =K∑k∇γ(xn+1/2) · (ua(tn)− uk(tn)) dt, (3.39)513.2. Unified Semi-Analytic Wall SPHwhere a is the particle under consideration, k is a triangular surface elementin the set K that are close to particle a.This method has been implemented in ArtiSynth and simulations areable to run stably with moving, irregular boundaries while maintaining sta-bility and good approximations of γ. For time-integrated quantities, numer-ical drift can present an issue. In these cases, re-initialization of γ explicitlyis possible. In our implementation, particles that are far from any solid wallshave their γ value forced to either 0 or 1, depending on whether it is insideor outside a containing geometry. This step effectively “re-initializes” γ sothat numerical drift does not become an issue even for long running simula-tions. For instance, simulation times of 30–60 s were used for the lid-drivencavities presented later in this chapter and drift of γ is never an issue.3.2.3 Non-Newtonian fluids in USAW SPHUntil now, the USAW SPH method has only described flow for Newto-nian boluses. Previous SPH simulations for non-Newtonian flows use staticboundary particles to represent solid walls, for example [30, 56, 109, 122].However, static boundary particles are not easy to define for irregular bound-aries, and deforming boundaries may change the boundary particle distri-bution in unpredictable ways causing errors in the SPH integral.Here we present a USAW SPH method for explicitly integrated non-Newtonian flow, extending the work of [122]. The main challenge is toevaluate the shear-rate, γ˙, from which an effective viscosity µeff(γ˙) can becomputed using one of many models such as power-law, Cross, or Herschel-Bulkley.523.2. Unified Semi-Analytic Wall SPHApplying the SPH formulation to the velocity gradient for a particle iyields the following definition for its components:kαβa =(duαdxβ)a=N∑bmbρb(uαb − uαa )∂Wab∂xβ, (3.40)where α, β ∈ {0, 1, 2}, and xα, uβ are the components of position and velocityin the α and β directions, resp. Wab is shorthand for W (xb − xa, h). In 3Dthe shear rate is defined asγ˙ =√2 tr(dαβ · dβα), (3.41)dαβ =12(duαdxβ+duβdxα). (3.42)γ˙a =[2((k00a )2 + (k11a )2 + (k22a )2) + (k10a + k01a )2 + (k20a + k02a )2 + (k21a + k12a )2]1/2.(3.43)Extending the definition for USAW SPH requires including boundaryterms into kαβa . For the typical USAW SPH derivation, the following ex-pression is obtained,kαβa =(duαdxβ)a=1γN∑bmbρb(uαb − uαa )∂Wij∂xβ− 1γK∑k(uαk − uαa )nβWim∆S,(3.44)where nβ is the βth component of the wall normal. This definition can beused in equation 3.43 in order to produce the correct shear rate using USAW533.3. Verification and validationboundaries.3.2.4 Time integrationThe simulation is initialized at time zero by seeding the domain with parti-cles and assigning initial values to velocity, density, γ. We use a semi-implicitEuler time integration scheme to advance the simulation. In order to ad-vance a simulation from time tn to tn+1,while tn < tend doUn+1 = Un + ∆t · F (Xn,Rn,Γn,Un)Rn+1 = Rn + ∆t ·D(Xn,Rn,Γn,Un+1)Xn+1 = Xn + ∆t ·Un+1Γn+1 = Γn + ∆t ·G(Xn+1)tn+1 = tn + ∆tend whilewhere Un, Xn, Rn, and Γn, are the vectors of all particle velocities, posi-tions, densities, and γ, respectively, at time n. The function F computesthe rate of change of velocities, D the rate of change of densities, and G therate of change of γ. ∆t is the timestep tn+1 − tn.3.3 Verification and validationThe USAW SPH method described in this chapter was implemented en-tirely in Java, as part of the ArtiSynth multibody simulation platform [77].It is used to simulate four viscous flows with known solutions: Couetteflow between two parallel planes, Hagen-Poiseuille startup flow in a cylin-543.3. Verification and validationder, lid-driven cavity flow, and steady-state Hagen-Poiseuille flow for a non-Newtonian fluid. For the Couette and Hagen-Poiseuille flows, theoreticaltransient and steady-state solutions are known [121]. For the lid-driven cav-ity flows, our results are compared to those computed by Ghia et al. in [39].Fluid particles were initialized with zero velocity, and mass equal tothe fluid density (initially constant throughout the domain) multiplied bythe differential volume they occupy (e.g. dV = dx dy dz or in the case ofcylindrical coordinates, dV = r sinφdφ dr dz).3.3.1 Couette flowIn the first flow regime, two infinitely large, solid planes are parallel to eachother, with the bottom plane at z = 0 having velocity 0, and the top planeat z = 1 having velocity 1 ıˆ. The space between the planes is filled withviscous fluid with ν = µρ = 1. At steady state, the theoretical solution forvelocity in the fluid is simply proportional to the height above the stationaryplane. For x = (x, y, z) = xˆı + yˆ + zkˆ,v(x) = zˆı. (3.45)For this simulation, particles were arranged in a regular lattice, withspacing size of 0.1 between particles. The planes were actually size 2 × 2,and velocity samples were taken for a single column of particles near thecentre of the plane. The simulation was driven by the top plane’s assignedvelocity.For the Couette flow, the theoretical solution is recovered almost exactly553.3. Verification and validation0 0.2 0.4 0.6 0.8 heightv xvelocityExactSPHFigure 3.5: Simulation results vs. exact solution for Couette flow betweentwo parallel planes. Ten particles span the distance between the top andbottom planes.even at low resolution, with only 10 particles spanning the channel (seeFigure 3.5). This simulation tests only the viscous part of the update,since density updates and pressures are not used in finding the solution. Inderiving the viscous term over the boundary, the assumption was made that∇v·n is constant over the intersection of the surface with a particles supportsphere. This holds exactly for this simple case. The assumption involvingthe face normal direction is also true for this case.563.3. Verification and validation3.3.2 Startup Hagen-Poiseuille flow for a Newtonian fluidThe second flow, Hagen-Poiseuille flow in a circular pipe, is driven by auniform pressure differential inside the pipe. The pipe geometry is definedwith its axis aligned with kˆ, with radius R = 1, and a fluid with ν = 1driven by a pressure differential equivalent to an acceleration of 1kˆ.The exact time dependent solution of laminar, fully-developed Hagen-Poiseuille flow is given by the infinite series:vz(r, t) =|B|4ν(r2 −R2) +∞∑m=1|B|R2να2mJ2(αm)J21 (αm)J0(rαmR)exp[−να2mtR2],(3.46)where vz is the axial component of the velocity, r is the distance from theaxial centre of the pipe, t is the time, and the pipe radius is R. ν = µρ isthe dynamic viscosity. Jn are Bessel functions of the first kind with order n,and αm is the mth root of J0(kr). B is the pressure force driving the flow,a body force in our simulations. In the limit as t → ∞, the steady statesolution is reached,vz(r) =|B|4ν(r2 −R2) (3.47)Particles were distributed regularly in the kˆ direction, and radially inthe direction perpendicular to the flow (see Figure 3.6 for an example ofthe distribution used in a moderate-resolution simulation). Care was takento initialize the particle positions to the centroid of the differential volumeelement they are supposed to represent. The simulation was driven by abody force in the axial direction of the pipe and evolved in time until the573.3. Verification and validation−1 −0.5 0 0.5 1−1−0.500.51Figure 3.6: Particle distribution for Hagen-Poiseuille in a plane perpendicu-lar to the flow direction for initial spacing of approximately 0.05. The tubewall is shown as a dashed circle.simulation reached a steady state. Periodicity was enforced with a methodsimilar to [110], where master particles near the outlet of the pipe havetheir densities, velocities, and pressures copied to their corresponding slaveparticles near the inlet. A similar but opposite scheme is used to enforceperiodicity at the outlet.For this experiment, the velocities for half of a single row of particlesspanning the diameter of the pipe are reported. Figure 3.8 shows the simu-lation predictions plotted against the transient solutions for both moderate-and high-resolution simulations. The predicted velocities show excellentagreement with the transient solution, even at moderate resolutions.Table 3.1 shows the L2 Error for a slice of particles perpendicular to583.3. Verification and validationthe flow direction and in the centre of the simulated pipe. The simulationappears to be converging to the correct solution at an approximately first-order rate (see Figure 3.7).Approx. Spacing L2 Error Apparent Conv. Rate0.1 1.01978× 10−3 −0.05 3.00186× 10−4 1.764330.025 2.38488× 10−4 0.331940.0125 1.25535× 10−4 0.92583Table 3.1: L2 Error and apparent convergence for Hagen-Poiseuille flowsimulation. Samples were for a slice of particles near the centre of the pipeperpendicular to the flow direction.10−2 10−110−410−3Approximate spacingL2ErrorConvergence plot for Hagen-Poiseuille flowL2 error1st orderFigure 3.7: The rate of convergence for Hagen-Poiseuille flow appears to beapproximately first-order.593.3. Verification and validation0 0.5 = 0.05t = 0.1t = 0.2t = 0.4t = 1.5r radial positionv zvelocityexactSPH(a)0 0.5 = 0.05t = 0.1t = 0.2t = 0.4t = 1.5r radial positionv zvelocityexactSPH(b)Figure 3.8: Hagen-Poiseuille flow simulation vs exact solution for the mod-erate (a) and high (b) resolution spacing. For the moderate resolution sim-ulations, the particle spacing is approximately 0.05 and 40 particles spanthe pipe diameter. For the high resolution simulations the particle spacingis approximately 0.0125 and 160 particles span the diameter. Only half ofthe pipe is shown due to symmetry.3.3.3 Lid-driven cavityThe third flow is the classic lid-driven cavity experiment. We solve a pseudo-2D problem where all of the kernels and solid geometries are 3D, but thecomputational SPH particles are only distributed in a single plane. Out-of-plane particles are simply duplicates of the computation particles, offsetby a distance equal to the initial particle spacing. In total, 8 out-of-planeparticles were created for each computation particle: 4 on either side. Inorder to compensate for small, out of plane velocities and displacementsaccumulating for each particle, every computation particle had the out-of-plane component of its position and velocity clamped to zero every 100 timesteps.In our tests, a box with size 1 × 1 was created, with a velocity of 1 on603.3. Verification and validationthe top boundary. The density of the fluid was set to 1000, and Reynoldsnumbers of 100 and 1000 were simulated for 30 and 60 time units, respec-tively.Streamlines results are generated with Paraview using the pv-meshlessplugin for SPH flows. Figure 3.10 shows the streamlines predicted by SPHafter steady-state has been achieved. The steady-state velocity profilesthrough the horizontal and vertical centres of the cavity show good agree-ment with the results reported by Ghia et al. [39] (see Figure 3.9). Thepredicted location of the centre of the primary vortex was within one gridcell of the results reported in [39] (grid cell sizes are defined by Ghia et al.).For Re = 100, secondary and tertiary vortices appeared and disappeared fre-quently during the transient stage, but were not present after steady-statehad been reached. For Re = 1000 the bottom right vortex was present atsteady-state and the location also matched well with that reported by Ghiaet al. In order to test that the results were “grid-converged”, the Re = 100simulation was run with a 200 × 200 particle distribution with only smalldifferences in both the velocity profiles, centre of main vortex, and stream-line shapes. The L2-norm of the difference in velocity profiles is less than1%.3.3.4 Non-Newtonian Hagen-Poiseuille flowThe Hagen-Poiseuille flow can also be run for non-Newtonian boluses. Wetest the shear-rate estimate (equation 3.44) by simulating a power-law fluid,for which an analytical solution of the velocity profile at steady-state isknown for flow in a circular pipe.613.3. Verification and validation0 0.2 0.4 0.6 0.8 1−0.200.2x positionyvelocityGhia(1982)SPH(a)0 0.5 velocityypositionGhia(1982)SPH(b)0 0.2 0.4 0.6 0.8 1−0.50x positionyvelocityGhia(1982)SPH(c)−0.5 0 0.5 velocityypositionGhia(1982)SPH(d)Figure 3.9: Velocity profiles for the lid-driven cavity flow for Re= 100 (toprow) and Re= 1000 (bottom row). Y-velocity through the horizontal centreof the cavity (left column) and x-velocity through the vertical centre ofthe cavity (right column). The SPH results show good agreement with theresults of Ghia (1982).623.3. Verification and validation(a) (b)Figure 3.10: The streamlines for the lid-driven cavity experiment for Re =100 (left) and Re = 1000 (right). A speed-of-sound of c = 5 was used, andthe simulation was run for 60 time units. The pseudo-2D simulation wasseeded with 100 × 100 (Re = 100) and 200 × 200 (Re = 1000) particles.Running the Re = 100 simulation with 200 × 200 resolution resulted innearly identical results.A power-law fluid is characterized by the parameters n and K, wheren < 1 models a shear thinning fluid, n = 1 models a Newtonian fluid, andn > 1 models a shear thickening fluid. K is a “flow-consistency index”. Theshear-stress of a power-law fluid is defined in 1D asτ = K(∂u∂y)n, (3.48)which can be written as a generalized Newtonian fluid,τ = µeff(∂u∂y), (3.49)µeff = K(∂u∂y)n−1. (3.50)633.3. Verification and validationIn multiple dimensions, the shear rate (γ˙) is defined by equations 3.41i–3.42. The effective viscosity isµeff(γ˙) = Kγ˙n−1. (3.51)A Hagen-Poiseuille flow for a shear-thinning power-law fluid was set up ina similar manner to the Newtonian simulation of flow through a cylindricalpipe. The velocity profile for a power-law fluid at steady state can be foundanalytically, and the solution is given byvz(r) =nn+ 1(∂p∂z12K) 1n(Rn+1n − |r|n+1n ), (3.52)where vz is the velocity in the axial direction, and x is the distance from thecentre of the pipe.A Hagen-Poiseuille flow for a power-law fluid was simulated with threedifferent resolutions, resulting in 11, 21, and 41 particles spanning the di-ameter of the pipe. The fluid has a density ρ = 1000, K = 1000, n = 0.8,making it moderately shear-thinning. The pipe has a radius R = 1. Thesimulations were run for one second of simulation time.Approx. Spacing L2 Error Apparent Conv. Rate0.2 7.67764× 10−3 −0.1 3.40116× 10−3 1.174630.05 1.71635× 10−3 0.98668Table 3.2: L2 Error and apparent convergence for a power-law Hagen-Poiseuille flow simulation. Samples were for a slice of particles near thecentre of the pipe perpendicular to the flow direction.643.4. Conclusion10−1.6 10−1.4 10−1.2 10−1 10−0.8 10−0.610−310−2Approximate spacingL2ErrorConvergence plot for power-law Hagen-Poiseuille flowL2 error1st orderFigure 3.11: The rate of convergence for power-law Hagen-Poiseuille flowalso appears to be approximately first-order.3.4 ConclusionThis chapter summarizes USAW SPH, the SPH method that we chose touse for simulating fluid boluses. We present USAW SPH derivation and itsrelation to the original derivations. We then show a illustrate a method forhandling moving boundaries as well as a method for calculating shear-rates,necessary for simulating non-Newtonian boluses. We demonstrated that themethod gives the expected results when simulating some canonical New-tonian and non-Newtonian flows including Couette flow, Hagen-Poiseuilleflow, and a lid-driven cavity flow. The Hagen-Poiseuille flow was performedfor both Newtonian and power-law fluids. Both demonstrate approximatefirst-order convergence in the L2 error norm.653.4. Conclusion0 0.5 100.10.2r radial positionv zvelocityexactSPH(a) 11 particles span the pipe0 0.5 100.10.2r radial positionv zvelocityexactSPH(b) 41 particles span the pipeFigure 3.12: Hagen-Poiseuille flow simulation for a power-law fluid plottedagainst the exact steady-state solution. The power-law fluid has K = 1000and n = 0.8.66Chapter 4Fluid simulation oforopharyngeal swallowingA multitude of complexities make it difficult to simulate oropharyngeal swal-lowing. The flow is truly 3D, and axisymmetric simulations oversimplify theproblem. There is a high degree of irregular boundary motion that drivesthe flow. The material properties of fluid boluses can also be very complex,and there is an interface with air that needs to be handled appropriately.Collecting 3D images of human swallowing is difficult, and until recentlythere have been no full 3D captures of human swallowing.After applying the extensions described in Chapter 3, USAW SPH iscapable of simulating swallowing in the oropharynx. This chapter describesa practical test of these capabilities, a numerical fluid simulation of swallow-ing. Section 4.3 describes a novel workflow and purpose-built tool used tosegment the airway from 320-row ADCT sequences of swallowing. The nu-merical simulation tests the assumption that saliva provides lubrication forstarch-thickened liquid boluses. The assumption is tested by running bothslip and no-slip simulations in dynamic 3D geometries of the aerodigestivetract based on dynamic 3D images. The results suggest that saliva does not674.1. The role of saliva in swallowinghave a significant effect on clinical timing measures of the bolus, however itmay decrease oropharyngeal residue at the end of the swallow. The smallnumber of subjects (3) limits the scope of the findings, but the experimentsare a successful proof of concept.4.1 The role of saliva in swallowingSaliva plays important roles in mastication, swallowing, and maintainingoral and general health. Decreased salivary flow may result in dyspha-gia [34, 100, 103] and videofluorosopic studies suggest that saliva is requiredfor lubricating the airway and bolus [43, 105]. However, other work finds thatthe relationship between hyposalivation and dysphagia is unclear [107, 111].One of the difficulties in understanding this relationship is that the hypos-alivation may be confounded by other factors, such as reduced tongue-baseretraction. Using a numerical simulation, lubrication as a result of salivacan be simulated in isolation from other pathologies.This work uses state-of-the-art 320-row Area Detector CT (ADCT) im-ages as the basis for computer fluid simulations of the bolus [57]. 3D simula-tions of lubricated and non-lubricated aerodigestive tracts in normal subjectswere performed, and the results were compared to the bolus position in theimages. We assumed that the lubricated simulations would produce betterresults than their non-lubricated counterparts. Six different swallows werestudied; three swallows were completed with starch-thickened boluses, bothnectar- and honey-thick. The non-Newtonian viscosities of the nectar- andhoney-thick boluses were measured experimentally. The results indicate that684.2. 320-row ADCT swallowing sequencesfor liquid boluses in normal subjects, the no-slip condition produces moreaccurate bolus transit times than a slip condition, and that saliva does notlubricate the bolus through the oropharynx. However, the no-slip simula-tions show significant amounts of residue compared to the slip condition.4.2 320-row ADCT swallowing sequencesWe collaborated with swallowing researchers from Fujita Health Universitywho acquire dynamic 320-row ADCT images of human swallowing in theirclinic [36, 57]. Their imaging technique is the first with the ability to capturedynamic, fully 3D images of human swallowing in a single shot (that is,without requiring repetitions of the swallow to form a single image sequence).The following dataset of six 320-row ADCT swallows were used to cre-ate the simulations. Each swallow was of a healthy Japanese subject whoswallowed normally. The swallows were from two subjects, as described inTable 4.1. Each image sequence consisted of approximately 2–3 seconds ofdata sampled at 10 Hz resulting in 20–25 3D “frames”. Each frame has512 × 512 × 320 voxels, with resolution of 0.468 × 0.468 × 0.5 mm3. Thisstudy was approved by the UBC Clinical Research Ethics Board, certificateno. H16-01546.694.3. Deriving geometric boundaries from 4D ADCT dataSubject Sequence description34 y.o. female reclining 45◦ - nectar thick 4ml35 y.o. female reclining 45◦ - thin 10ml31 y.o. female reclining 45◦ - thin 10mlreclining 45◦ - honey thick 10mlsemi-prone - thin 10mlsemi-prone - honey thick 10mlTable 4.1: Scans from three female subjects provided the data for this sim-ulation study. All subjects were healthy with no history of dysphagia, andnot taking any medication.4.3 Deriving geometric boundaries from 4DADCT dataFor each ADCT sequence, a dynamic mesh needs to be created from theimages. A mesh is an interconnected set of triangles that approximates asmooth, 3D surface. The vertices of each triangle can be moved in time,causing the mesh surface to move and thus the 3D surface changes shape.Each moving mesh is used as a solid wall boundary condition for the fluidsimulation of the bolus.This section describes a workflow and newly created tools for extractingmoving boundaries from 4D images of swallowing. It is published in [54].Existing methods, both commercial and research, were found to be inade-quate for extracting a dynamic boundary of the upper airway from dynamicCT for fluid simulation (see Section 2.4).Working with dynamic ADCT data has its challenges. There does notappear to be a suitable automated technique for automatically labelling thedata. Each sequence has high spatial resolution but only moderate temporal704.3. Deriving geometric boundaries from 4D ADCT dataresolution (10 Hz). During the reconstruction of the images from the scannerdata, artifacts can appear as a result of metal fillings, subject motion, orbolus motion. The artifacts may appear as double images (for example, twohyoid bones in the same instance), blurring of edges, dark and bright radialspokes that emanate within an axial slice, and dark spots that appear next tobright moving objects. Additionally, openings such as the upper esophagealsphincter (UES) are normally closed in the CT with no discernible edges(for example, see Figure 4.2(a)), however for USAW SPH simulation, theboundary is required to be topologically constant in time, and thereforethese closed sections need to be represented at all time instances. Despitethese shortcomings, dynamic ADCT is the only imaging modality that hasdemonstrated its ability to capture human swallowing in full 3D.A triangular surface mesh would provide a suitable representation ofthe aerodigestive tract (the airway). By deriving the motion of the meshvertices from the ADCT images, a moving 3D boundary would be obtained.However, achieving temporal continuity is a challenge. For example, theUES forms the boundary between the pharynx and the esophagus. It mustbe segmented correctly in all time instances, even when it is closed. However,in the CT images the UES is squeezed flat, and thus is invisible for most ofthe sequence (see Figure 4.2(a)). The same issue applies to tongue/palatecontact and other anatomic regions that close during the swallow in orderto propel the bolus or protect the airway.We attempted to use existing 3D image-to-image registration algorithms(such as the Diffeomorphic Demons algorithm [119] as implemented in theInsight Toolkit) to generate deformation fields that could warp an exist-714.3. Deriving geometric boundaries from 4D ADCT dataFigure 4.1: The segmented initial mesh (translucent blue) overlaid on themid-sagittal slice to illustrate the anatomical area being considered.724.3. Deriving geometric boundaries from 4D ADCT dataing mesh at one time to a subsequent time frame, but this approach wasunsuccessful. A similar approach using the Elastix toolkit [66] was also un-successful. We soon concluded that even if a suitable automated methodexisted, it would require at least some manual intervention.4.3.1 BlendSegTo address the need for manual deformable registration in 3D, BlendSeg wascreated. It is a tool based on the open-source and extensible 3D modelingsoftware Blender (Blender Foundation, Amsterdam, Netherlands). BlendSegfacilitates the registration of a triangulated mesh to a 3D volumetric image.BlendSeg makes use of Blender’s “sculpting” feature that allows the operatorto morph a mesh using a variety of sculpting brushes while maintaining thesame mesh topology (connectivity). Keeping the topology constant allowslinear interpolation between vertices in the temporal direction, leading toa continuous deforming mesh. BlendSeg consists of two features missing inBlender : the first is a method to efficiently display volumetric data as slicesin the same space as the mesh, and the second is a method to compute andrender the intersection between the mesh and a plane.4.3.2 Initial mesh generation using AmiraThe initial triangular mesh can be created using any number of tradi-tional segmentation tools. Amira was chosen to create the initial segmen-tation manually. The surface was tesselated into triangles, simplified, andsmoothed using both Amira and MeshLab software. This initial mesh cannow be deformably registered to subsequent time frames (see Figure 4.2).734.3. Deriving geometric boundaries from 4D ADCT data(a) At this initial time, there is no motionand the airway is clearly rendered. How-ever there is closure of the esophagus, andthe soft palate seals the oral cavity fromthe pharynx.(b) The segmented initial mesh overlaidon the mid-sagittal slice to illustrate theanatomical area being considered.Figure 4.2: Generating an initial mesh from time t=0. Mid-sagittal slicesare shown.4.3.3 Visualizing 3D volumetric data in BlenderBlender’s primary function is to create and modify piecewise-planar surfacemeshes. Currently, there is no method for visualizing 3D imaging data inBlender. A powerful application programming interface (API) to Blender’sunderlying engine is provided through the Python programming language.BlendSeg is a Python plugin, and it facilitates the visualization of 3Dvolumetric data in the Blender sculpting window. 3D volumetric data (e.g.DICOM data) are converted to three “stacks” of images, each comprised ofan axial, coronal, or sagittal stack of images, and imported into Blender.Three orthogonal planes are created in Blender with each corresponding toone of the anatomical planes. The planes may be moved interactively bythe user, however the motion is constrained to translation along the normaldirection only. When the planes are translated in space, the displayed tex-ture is updated. In a 4-pane setup (Blender’s quad-view option), this gives744.3. Deriving geometric boundaries from 4D ADCT dataFigure 4.3: Blender’s quad-view, when combined with BlendSeg, give theusual anterior-posterior, superior-inferior, and lateral views (bottom left,top left, and bottom right, respectively). The upper right view gives an in-teractive 3D view which can be used to select which slice is being visualized,as well as to perform sculpting on the mesh. Here, the mesh is being hiddenso that the intersections (orange) can be seen more clearly.the typical anterior, superior, and lateral views (see Figure 4.3). The fourthpane is used for an interactive 3D view, where the planes can be translatedby the operator and sculpting of the mesh can be performed.4.3.4 Generating intersection contoursThe second component of BlendSeg is intersection contour rendering. Theintersection contour between two surfaces, S and G, is defined as the set ofpoints x ∈ (S ∩ G). BlendSeg computes the intersection contours betweenthe mesh and plane by searching for triangles that cross over the currentposition of each plane, using an algorithm similar to that described in [4].Each triangular mesh stores its connectivity using a “quad-edge” data754.3. Deriving geometric boundaries from 4D ADCT datastructure. For the general problem of finding intersections between twodiscretized surfaces, a tree of axis-aligned bounding boxes (AABBs) can beused to quickly minimize the amount of tests that need to be performed.From this set of potentially intersecting pairs of triangles, a set of MeshIntersection Points (MIPs) is created.Any intersection contour between two triangular meshes can be describedby an ordered set of MIPs, with straight lines connecting each MIP to thenext. Considering two triangles that intersect with one another, the MIPsare the points where an edge of one triangle intersects the face of the othertriangle. For a single pair of triangles that potentially intersect, zero ortwo MIPs exist. Once a set of MIPs are found, at least one contour ex-ists. Beginning at an arbitrary MIP, one can find neighbouring MIPs byutilizing the quad-edge connectivity data to “walk” around the intersectioncontour, removing MIPs from the set until reaching the beginning MIP.Open contours are possible as well, and are handled by reversing directionupon encountering an open edge (such as the edge of a plane).For each intersection point and edge representing the contour, new ver-tices and edges were defined in Blender and rendered. Computing inter-section contours and rendering them is performed at interactive rates ona laptop computer running Debian Linux with an Intel Core i5 proces-sor and 16 GB of memory. Moving the orthogonal, textured planes inBlender updates the image as well as the intersection of the mesh withthe plane, allowing the user to view how well the mesh lines up with the3D image data (see Figure 4.3). BlendSeg is freely available for downloadat Deriving geometric boundaries from 4D ADCT data4.3.5 Interpolating in time to create a moving boundaryBlender’s sculpting feature provides a number of “brushes” allowing theoperator to push, pull, grab, smooth, and otherwise modify the mesh locallyand interactively. The mesh will maintain the same connectivity as long asthe option for dynamic topology (DynoTopo) is disabled. This feature allowsthe operator to sculpt the mesh to match each 3D volume of the swallowsequence. Each vertex position can be interpolated between time points tocreate a continuous 4D mesh representation of the airway during the swallowsequence.Different interpolation schemes could be used, but simple linear inter-polation provides the best results. Linear interpolation guarantees that theinterpolated value will exactly match the specified value at knot points andthat there won’t be any overshoot. We found this to be important whendealing with opening and closing of the glossopalatal junction.Special mention should be made regarding tangential mesh motion, i.e.vertex motion constrained to the surface of the mesh. For a fluid simulationusing no-slip boundary conditions, shearing forces would be exerted on thefluid in the direction of the vertex motion. In reality, surfaces such as thehard-palate exhibit no tangential motion, but other surfaces like the tonguemay exhibit substantial tangential motion that is not captured in the CTimages. The CT data is unable to directly capture this type of tangentialboundary motion, however it may occur in the 4D mesh as an artifact ofthe sculpting process. To eliminate this effect on the fluid simulation, whencomputing the vertex velocities on the boundary, the tangential component774.3. Deriving geometric boundaries from 4D ADCT dataFigure 4.4: Example mid-sagittal slice from time t=0.5s showing the inter-section of the mesh with the slice.of the velocity is ignored.4.3.6 Extracting the boundariesThe process of using BlendSeg to register a mesh to a single 3D image cantake anywhere from a few minutes to a few hours depending on the amountof sculpting required. After creating meshes for four swallow sequences,each consisting of between 22 and 25 individual 3D images, the total timefor creating a moving airway is estimated to be around 30 hours. This doesnot include the time it takes to create an initial mesh. The first time a meshof an anatomic structure is created (for example the tongue, or airway), amesh needs to be created using traditional tools (such as Amira) which cantake a significant amount of time depending on the tools used. If a meshwith approximately similar geometry already exists, it can be registered toa new initial frame quickly using BlendSeg. Since temporal continuity is notrequired from the prior mesh to the new initial mesh, new vertices may beinserted and topology changes are permitted.To evaluate each registered mesh, two experts familiar with the anatomy784.4. Viscosity of thickened bolusesand physiology of swallowing were enlisted. One is a practising clinicianand the other is an anatomist who specializes in swallowing function. Theexperts verified the registrations by visual inspection of the intersectionsin BlendSeg. The method of handling airway closures (such as the UES,vocal folds, palatal-glottal contact, and nasopharynx) was explained to theexperts. To evaluate each registration, a first pass was made looking at thesagittal slices that covered the airway. The lateral videofluoroscopic imageis the most common method of evaluating dysphagia, and so the sagittalslices were chosen for the first pass (see Figure 4.4). Following the firstpass, closer examination was performed for areas of the image that includedmotion artifacts such as blurring or ghosting. In these difficult-to-interpretsections, context from nearby sagittal, axial, coronal, and temporal imageswere required.4.4 Viscosity of thickened bolusesViscosity is a material property that represents the amount of resistance toshearing of a fluid, here represented with the symbol µ. Imagine the tongueholding a liquid bolus against the hard-palate. A shearing motion is onewhere the tongue moves parallel to the hard palate. For a “thicker” fluidsuch as honey, there is significant resistance to that shearing motion, andthus the viscosity is higher. For a “thin” fluid such as water, there is littleresistance, and thus the viscosity is low.A Newtonian material, such as water, has a value of µ which is con-stant. In our simulations we assigned water a viscosity of 1 cP (centi-794.4. Viscosity of thickened bolusesPoise). Boluses thickened with starch exhibit non-Newtonian shear-thinningbehaviour, and the value of µ depends on the shear-rate. To go back to thetongue analogy, a shear-thinning bolus would show high resistance to shear-ing when the tongue is moving slowly, but as the tongue speeds up, theresistance would lower. The Cross model of viscosity computes an effectiveviscosity based on the shear rate according to the following equation:µeff(γ˙) =µ01 + (µ0τ∗ γ˙)n+1, (4.1)where µ0, τ∗, and n characterize the behaviour of the fluid. µ0 has units ofviscosity [Pa·s], τ∗ has units of stress [Pa], and n is dimensionless.The viscosities of the nectar- and honey-thick boluses were measuredin order to simulate them. A 100 ml sample of nectar-thick bolus consistsof 95 ml water, 5ml barium suspension, and 2.5 grams of thickener. Thehoney-thick bolus is the same as the nectar except 5 grams of thickener isused. The thickening agent is Top Balance III (Neo-hai-torome-ru; FoodcareInc., Sagamigahara, Japan) which comes in 2.5 gram packets. Before beingadministered to the subject for imaging, the bolus is chilled in a refrigeratorfor 4-5 hours.A calibrated KINEXUS rheometer (Malvern Instruments, Malvern UK)was used to measure the viscosity of the boluses. Samples of the nectar- andhoney-thickened boluses were prepared using distilled water. A cup-and-bobgeometry was used to measure the boluses. The boluses were measured at atemperature of 10◦C. Each bolus type underwent the following protocol: afive minute pre-shear was run in order to stabilize the sample and the tem-804.4. Viscosity of thickened bolusesperature. A test was performed by varying shear rate between 0.01 s−1 and1000 s−1. After the first ramp-up test was performed, a second test was per-formed on the same sample after a five minute delay. A new sample was thenloaded into the machine, underwent a five minute pre-shear/temperaturestabilization, and a third ramp-up test was performed. To check whetherthe refrigeration period of the sample makes any difference, a final pre-shearand ramp-up test was performed on a sample which had been refrigeratedfor 3 hours in order to make sure that refrigeration has no effect on themeasurements.4.4.1 Slip and no-slip boundary conditionsThe typical boundary condition used in fluid flow problems is no-slip. Fluids“stick” to solid boundaries and their velocity is equal to that of the boundary.For fluid flowing through a stationary rigid pipe, the velocity of the fluid atthe wall is zero. We use the no-slip boundary condition to approximate adry-mouth swallow, since there is no saliva lubricating the bolus.For healthy subjects, such as those in the images that are being used inthis study, a layer of saliva coats the solid boundaries. The saliva lowers thefriction between solid-solid contact (such as between the tongue and palate).As a boundary condition for fluid flow in a saliva coated airway, we use afull-slip boundary condition on the fluid. The full-slip condition constrainsthe fluid velocity to prevent penetration into the solid, but does not restrictits velocity in a direction tangent to the solid surface.As shown in Appendix B, the full-slip condition can be a reasonable ap-proximation for lubricated bolus flow, however the quality of the assumption814.4. Viscosity of thickened bolusesdepends on the viscosity and width (quantity) of saliva coating the airway,as well as the duration of the swallow. In reality, saliva viscosity and widthcan vary based on many factors. It can change between individuals, withinthe same individual based on time of day, whether the salivary flow is stim-ulated or not, and even depending on the location in the airway where themeasurement is made. We expect that the results of a full-slip and no-slipsimulation will bracket the actual behaviour of the bolus flow, however forthe time-scales involved in swallowing, and for some choices of saliva viscos-ity and quantity, the full-slip condition is a reasonable approximation to asaliva-lubricated flow.4.4.2 Bolus measurementsBolus based measurements were made on both the CT data and the sim-ulation data. The measurements chosen are entirely determined by bolusposition, and as such are suitable for measuring the simulation data. Theoral transit time (OTT), pharyngeal transit time (PTT), and residue weremeasured. Three time points were marked for each swallowing sequence:the time of the first backwards motion of the bolus head, the time whenthe bolus head cross the ramus of the mandible, and the time when thetail of the bolus clears the cricopharyngeal region. The OTT is defined asthe difference between the first two times, and the PTT is defined as thedifference between the second and third times.All images were presented from a lateral, orthographic view. The CTdata was transformed into a ”virtual videofluoroscope” that compressed each3D time frame into a 2D view that mimics a traditional VFL image. The824.5. Resultssimulations were presented from a lateral, orthographic view (see Figure 4.9,left and right columns). CT frames were 0.1 seconds apart, and the simu-lation time frames were 0.02 seconds apart. The measurements were madeby a Speech Language Pathologist (SLP). A second SLP also evaluated thedata to confirm the results of the first. The second SLP was blinded to theresults prior to the evaluation. In cases where their measures differed by0.05 seconds or more, the SLPs came to a joint consensus about what timeshould be chosen. The time interval 0.05 s corresponds to three or moreframes in simulation, and one or more frames of ADCT data.4.5 ResultsThe viscosity measurements for the starch-thickened boluses are presentedin a plot measuring the apparent viscosity as a function of strain-rate, inFigure 4.5. For each bolus type (nectar and honey), four tests were run on3 different samples, and there were almost no differences observed betweenthem, demonstrating high repeatability.The viscosity curves were modelled using the Cross viscosity model, asdescribed in equation 4.1. A power-law model was first fit to the experimen-tal data, and then the Cross model parameters were chosen to match thepower-law model. The parameters are given in the Table 4.2. The modelsof the viscosity curves are also plotted in Figure 4.5.All of the boluses were simulated with initial density of ρ0 = 1000kgm3,and an initial spacing of 1 mm between particles, with 1000 particles foreach millilitre of fluid. The artificial speed-of-sound was c = 10 ms for all834.5. Results10−2 10−1 100 101 102 10310−1100101102103Shear-rate (γ˙) [s−1]ApparentViscosity(µ)[Pa·s]Viscosity vs. Shear-RateExper. NectarCross NectarExper. HoneyCross HoneyFigure 4.5: Plot of apparent viscosity vs. shear-rate for nectar- and honey-thick boluses at a temperature of 10◦C. Note that these are log-log axes.The materials have very similar curves, and are both highly shear thinning.Nectar-thick Honey-thickn 0.2156 0.2156µ0 [Pa·s] 3000 4000τ∗ [Pa] 9 14.5Table 4.2: Parameters of the Cross model for the Nectar- and Honey-thicksimulations.844.5. Resultssimulations.The simulations were able to run mostly unsupervised and each took be-tween 1−21 days to complete using machines with multi-core Intel i7 CPUs.The reason for the long computation times were that starch-thickened bo-luses required smaller simulation time-steps to avoid numerical instabilitybecause of the viscous diffusion limit on the maximum step size, while thethinner boluses had time steps limited by the CFL condition.4.5.1 Thickened bolusesFor the starch-thickened bolus simulations, we expected to see a large dif-ference between the slip condition and no-slip condition simulations. For allthree of the thickened bolus sequences, we found that the no-slip conditionproduced a more accurate approximation of the bolus flow compared withthe slip condition. The position related timings of the bolus are in Table 4.3.As a test for resolution independence, the nectar-thick simulation was per-formed with an initial spacing of 0.8 mm, resulting in 7812 particles insteadof 4000. The results for both the slip and no-slip simulations showed nosignificant differences.OTT and PTT are useful measures when comparing both inter- andintra-patient swallows, however since we are looking at a single swallow andcomparing simulation to CT scans, we are actually performing intra-swallowcomparisons. In this case the absolute time measures become relevant. InFigures 4.6–4.8, the times chosen for “first backwards motion of the bo-lus”, “bolus head crosses the intersection of the ramus of the mandible andbase of tongue”, and “bolus tail leaves cricopharyngeal area” are indicated854.5. ResultsCT Slip No-slipNectar(reclined)OTT [s] 0.3 0.3 0.26PTT[s] 0.5 0.5 0.52Residue [%] 0 0.4 43.0Honey(semi-prone)OTT [s] 1.1 0.98 1.12PTT[s] 0.5 0.72 0.58Residue [%] 0 18.53 29.02Honey(reclined)OTT [s] 0.9 0.32 0.8PTT[s] 0.7 0.98 0.7Residue [%] 0 15.71 31.07Table 4.3: Summary of transit times and residue for CT, and simulatedslip/no-slip for each swallow sequence.0 0.5 1 1.5 2 2.5SlipCTNo-slipTime [s]4 ml Nectar-thick, reclinedOTT PTTFigure 4.6: Timing chart for 4ml nectar-thick bolus, reclined position.discretely. The red and green bars indicate the OTT and PTT, resp.Nectar-Thick - Reclined PositionFor the 4 ml nectar-thick (reclined) sequence, the simulated slip conditionbolus leaves the oral cavity before the no-slip bolus. Comparing the oraltransit times of the CT, slip, and no-slip conditions (Table 4.3), both of thesimulated OTTs were 0.04 s away from the CT measurement, however theylie on opposite sides of the CT measurement, with the no-slip simulation864.5. Results0 0.5 1 1.5 2 2.5SlipCTNo-slipTime [s]10 ml Honey-thick, semi-proneOTT PTTFigure 4.7: Timing chart for 10ml honey-thick bolus, semi-prone position0 0.5 1 1.5 2 2.5SlipCTNo-slipTime [s]10 ml Honey-thick, reclinedOTT PTTFigure 4.8: Timing chart for 10ml honey-thick bolus, reclined position874.5. ResultsSlip CT No-slipFigure 4.9: Nectar-thick bolus, semi-prone position. Columns (L-R): slip,3D CT, no-slip.884.5. Resultshaving a shorter OTT. The reason for this is that the first backwards move-ment of the slip simulation was much earlier due to poor oral containment.For PTT the slip and no-slip give very similar measurements. In this case,looking at the timing chart in Figure 4.6 illustrates how the OTT is betterapproximated by the no-slip condition. The bolus leakage from the oralcavity in the slip condition means that the PTT begins earlier. Therefore,the bolus tail leaving the cricopharyngeal region is also earlier for the slipcondition.These differences can be observed in Figure 4.9. In the CT images, itappears that the squeezing motion of the tongue, palate, and pharynx isthe dominant driver of bolus motion as evidenced by the lack of air nearthe bolus tail (see Figure 4.9, row 2). However, for the slip condition an airbubble is visible trailing the bolus tail.When quantifying the residue remaining in the oral cavity and pharynxafter the swallow, the no-slip simulation has a large quantity of residue inthe oral cavity and pharynx, however the slip has almost none. The CTimages show only trace amounts of residue after the swallow.Honey-thick - Semi-Prone PositionThis sequence consists of a normal subject lying semi-prone in the CT scan-ner swallowing a 10 ml honey-thick starch-thickened bolus. The measuredOTT and PTT of the no-slip bolus condition shows better agreement withthe CT data compared with the slip bolus condition (see Table 4.3).Visually, the slip and no-slip simulations both show similar results to theCT images (see Figure 4.10). Neither the slip nor no-slip conditions have an894.5. Resultsair bubble trailing the bolus tail. The most obvious visual difference is foundwhen looking at the interface between the bolus head and the air when itstarts to leave the oral cavity (see Figure 4.10, second row). An alternateview can be found in Figure 4.11 where slices of the simulations are showninstead of lateral projections.When considering residue, both of the simulations have significant residuein the oral cavity, valleculae, and piriforms after the swallow, with moreresidue observed in the no-slip condition simulation compared with the slipcondition simulation (see Table 4.3).Honey-thick reclinedIn this sequence the subject sat in a reclined position and swallowed a 10 mlhoney-thick starch-thickened bolus. The most striking visual difference be-tween the slip and no-slip simulations is the significant pre-swallow leakageof the bolus from the oral cavity during the slip condition simulation (seeFigure 4.12, second and third rows). Similar to the nectar thick swallows,for both the CT and no-slip images, the bolus tail is cleared from the oralcavity by the squeezing motion of the tongue against the palate, while theslip condition exhibits a bolus tail that is trailed by an air bubble.Looking at Table 4.3, as well as Figure 4.8, the no-slip condition simu-lation shows better agreement with the CT images while the slip conditionmeasure for OTT is much shorter due to the leakage of the bolus and thespeed with which the leakage travels to the intersection of the ramus of themandible and the tongue base.Both the slip and no-slip simulations show significant residue in the904.5. ResultsSlip CT No-slipFigure 4.10: Honey-thick bolus, semi-prone position. Columns (L-R): slip,3D CT, no-slip. The shape of the air-bolus interface at the bolus headshows better agreement to the CT images in the no-slip condition than theslip condition in the oral cavity and pharynx.914.5. ResultsSlip CT No-slipFigure 4.11: Detail of bolus-air interface at a single mid-sagittal slice, atthe initiation of swallowing. The interface of the no-slip bolus shows betteragreement with the original images than the slip bolus.oral cavity and pharynx after the swallow, with more residue observed inthe no-slip simulation. The quantities of residue are similar to the honey-thick semi-prone swallow simulations. For the no-slip simulation, numericalinstability occured before the end of the swallow, however, we have estimatedthe quantity and location of residue based on the last stable time.4.5.2 Thin bolusesThe thin boluses were simulated as Newtonian materials with a viscosity of1 cP. All three sequences used boluses of 10 ml. Because the viscosity was solow, we expected the slip and no-slip simulations to be very similar. In thefirst reclined simulation, the no-slip simulation showed had a small amountof residue (9 particles) and a single particle was aspirated. In the other twosimulations, both the slip and no-slip simulations showed small amounts ofaspiration (less than 3% of total volume) while the CT images showed none.The semi-prone swallow exhibited approximately 5% oropharyngeal residue,while none was detectable in the CT images.924.5. ResultsSlip CT No-slipFigure 4.12: Honey-thick bolus, reclined position. Columns (L-R): slip, 3DCT, no-slip. The slip bolus escapes the oral cavity very easily as shown inthe second and third rows. The no-slip bolus remains in the oral cavity,showing better agreement with the CT images.934.5. ResultsCT Slip No-slipThin(reclined)OTT [s] 0.2 0.18 0.14PTT[s] 0.3 0.54 0.56Residue [%] 0 0 0.09Aspiration [%] 0 0 0.01Thin(reclined)OTT [s] 0.3 0.24 0.22PTT[s] 0.1 0.3 0.38Residue [%] 0 0.69 1.36Aspiration [%] 0 2.92 2.47Thin(semi-prone)OTT [s] 0.2 0.14 0.14PTT[s] 0.5 0.4 0.44Residue [%] 0 4.42 5.65Aspiration [%] 0 1.75 2.55Table 4.4: Summary of transit times and residue for thin boluses. CT, sim-ulated slip, and no-slip for each swallow sequence. Some simulated boluseswere “aspirated”, and percentage volume is indicated as well.0 0.5 1 1.5 2 2.5SlipCTNo-slipTime [s]10 ml thin, reclinedOTT PTTFigure 4.13: Timing chart for 10ml thin bolus, reclined position.944.5. Results0 0.5 1 1.5 2 2.5SlipCTNo-slipTime [s]10 ml thin, reclinedOTT PTTFigure 4.14: Timing chart for 10ml thin bolus, semi-prone position0 0.5 1 1.5 2 2.5SlipCTNo-slipTime [s]10 ml thin, semi-proneOTT PTTFigure 4.15: Timing chart for 10ml thin bolus, reclined position954.6. Discussion4.6 DiscussionThe experiments in this work test the assumption that in healthy individuals,saliva lubricates starch-thickened boluses during swallowing. We performeda series of simulations of healthy swallows based on dynamic CT data and,for thickened boluses, found that the no-slip condition generally gives betterpredictions of bolus transport and transit times than a slip condition. How-ever, the no-slip condition also predicts a much higher quantity of residuethan the slip condition, and that are visible in the ADCT images.The results of the simulations seem to imply that in normal subjects,the presence of saliva does not have the effect of lubricating the bolus toshorten its transport times through the airway. Saliva has a much lowerviscosity than the thickened boluses in the sequences we studied. Witha moderately thick (or wide) layer of saliva coating the airway, we wouldexpect to see measurable shortening of bolus transport times and a quickertransport through the oral cavity and pharynx (see Appendix B). However,the no-slip simulations exhibit a higher quantity of oropharyngeal residueafter the swallow. This could be interpreted as evidence that saliva helps todecrease the residue after a swallow.4.6.1 Thickened bolusesOne of the biggest differences between the two boundary conditions werethat with the slip condition, the bolus appears to be influenced much moreby gravity. This is evidenced by the air near the tail of the slip bolus, whichis absent in the no-slip and the CT images. A number of physical phenomena964.6. Discussionmight explain the observations. One explanation is that the saliva is presentin a thin enough layer that it does not significantly affect the bolus. For thisto be realistic, consider a gravity driven bolus (such as the reclined positionwith leakage out of the oral cavity, Figure 4.12). For a saliva layer at least100 microns thick, we would expect the full-slip approximation to give agood approximation to the CT images (see Appendix B).Instead, the no-slip condition appears to give a much better approxima-tion. One possible explanation is that the saliva layer has a thickness on theorder of only a few microns or less. If we consider a flow-rate based analysis,as in Appendix A, thin portions of the bolus would be influenced more bya thin layer of saliva than thicker portions. However, a microns thick salivalayer is extremely thin and it is not obvious if this is the case in reality.Another explanation is that the saliva is quickly absorbed into the bolusso that it no longer provides lubrication, or similarly, that the saliva isquickly removed from the airway surface by the passing bolus so that itslubricating effect would not significantly speed up the transport of the bolus.It could also be a combination of the two, however the important takeaway isthat it does not appear that saliva lubricates liquid boluses in this manner.Solids do not follow the same rules as fluids with respect to friction, thereforesaliva might still play a very important role in lubricating solid and semi-solid foods during the swallow and lowering the coefficient of friction.ResidueThe thickened bolus simulations all demonstrated higher quantities of residuewith a no-slip condition than a slip condition. This might suggest that saliva974.6. Discussionmay play a role in decreasing residue, however there are other possible ex-planations for this behaviour. One is that the residue in the simulation is afunction of the moderate number of particles. The no-slip condition makesthe assumption that an infinitesimally small layer of bolus always coats thesurface of the aerodigestive tract. For the particle resolution used in thiswork (1 mm diameter), the small layer has a finite and significant volumethat contributes to the residue at the end of the swallow. Increasing theresolution of the simulation may decrease the amount of residue because theno-slip layer would then consist of smaller particles. However it is not clearif the volume of residue would approach zero as the particle size tended tozero, or if it would converge to some small but finite volume.The residue in the simulations might also be due to the difficulty infully closing the virtual pharynx (a small gap must be left in order for thesimulation to remain stable), or obtaining the correct timing of the pharynxclosure. It is also possible that saliva has a larger effect when lubricating thinstreaks of boluses, such as the bolus tail, and so they are more effectivelycleared by the squeezing motion. Another possibility is that the temperatureof the bolus tail is higher than 10◦C, and thus has a lower viscosity than theassigned viscosity.4.6.2 Thin bolusesBecause the thin boluses were simulated with a very low viscosity (similar towater), we did not expect the slip and no-slip to have very different results,and that is what we observed. The measured timings were all very similar,within 0.05 s or less, except in one case where the PTT was increased by984.6. Discussion0.08 s in the no-slip simulation. In addition, the differences in simulatedresidue and aspiration between slip and no-slip were within 1% of the totalbolus volume.We suspect that the simulations show aspiration because thin bolusesmove faster than thickened boluses, and as such, the 10 Hz may be insuffi-cient to adequately capture the speed of the epiglottic inversion and vocalfold closure. ‘Because the simulations showed aspiration while the originalCT images did not show any, and because the slip and no-slip simulationswere very similar, they did not give evidence towards saliva’s role in swal-lowing with respect to lubrication.4.6.3 Corroboration with existing simulation and clinicalresearchIn Sonomura et al. [112] the authors use a slip condition to model the pres-ence of saliva. However, they also found that they needed to decrease thestrength of gravity to one third of its actual value in order to obtain real-istic bolus velocities. Our findings suggest an explanation for their results,that the slip condition would have unrealistically sped up the bolus in theirsimulation. In the clinical research literature, Sonies et al. [111] found thatdespite the fact that healthy individuals demonstrate a wide range in sali-vary flow rates (salivary volumes), they are very similar with regard to thecharacteristics of their oropharyngeal swallowing of saliva and thin fluidboluses. Rogus-Pulia & Logemann [107] examined swallowing in patientswith Sjo¨gren’s Syndrome compared with healthy controls. They found thatalthough patients with Sjo¨gren’s syndrome perceive their swallows to be994.6. Discussionimpaired, few statistically significant differences were found between thepatients with Sjo¨gren’s syndrome and the normal age-matched controls ontemporal measures of swallowing of thin, pudding thick, and solid food bo-luses. These findings appear to be consistent with our simulation results.Hamlet et al. [43] performed a retrospective study of head and neck cancerpatients with postradiation xerostomia and normal controls swallowing forbarium liquid, paste, and shortbread cookies. They found that bolus transittimes were unaffected compared to healthy controls, but oral and pharyngealresidue was increased in the patient group, also consistent with our results.However, other studies have reported that hyposalivation is associatedwith significant lengthening of bolus transit times [12, 105], but this may bedue in part, to variation in study designs, assessment methods, and selectionof outcome measures. For instance, in [12], the authors perform an age-matched comparison between subjects with Sjo¨gren’s Syndrome and healthycontrols, and found that the Sjo¨gren’s subjects exhibited longer times for awater bolus than the healthy controls. However, they used ultrasound tomeasure hyoid bone motion and did not observe the bolus flow directly withVFL. In [105], the authors also compare patients suffering from Sjo¨gren’ssyndrome with healthy controls, and use VFL to measure timings for basalswallows (no bolus) and 10 ml water bolus swallows. However, their studydid not use barium as a contrast agent and therefore do not observe thebolus directly. Their timing measures are based on the motion of observableanatomical landmarks, such as the hyoid bone.1004.6. Discussion4.6.4 Limitations and assumptionsEvery computer model relies on making simplifying assumptions to keepsimulation times tractable. The main assumptions in this simulation are thefollowing: temperature is constant (isothermal), chemical effects are ignored(e.g. amylase in saliva is known to significantly decrease the viscosity ofstarch-thickened boluses), the effect of air is ignored, the fluid is allowed tobe slightly more compressible than in reality, and the simulation method onlyallowed for a moderate number of particles. The simulations also assumethat the saliva presence is uniform throughout the oropharynx, however inreality the thickness of saliva varies and depends on multiple factors.Of these assumptions, the chemical effects are the most likely to influencethe simulation. Air has a low density relative to liquids so that its effectis likely negligible. Temperature changes could affect the simulation in oneof two ways: by energy exchange and therefore having an effect on themomentum of the fluid, or by changing the material properties of the bolus.Although there is a significant temperature difference between the humanbody and the bolus, the short duration of the event makes it unlikely that thetemperature of the bolus changes significantly during the swallow. Artificialcompressibility is controlled by a speed of sound parameter, and it has beenreported that choosing a speed of sound 10 times greater than the maximumexpected velocity results in a less than 1% change in volume. Our chosenspeed of sound (10 m/s) is much higher than ten times the maximum bolusviscosity, and therefore should not have a significant effect on the simulationoutcome. Chemical effects, that is, the mixture of amylase with saliva,1014.6. Discussionhave been shown to have a significant effect on viscosity over short periodsof time [44]. However, their study found that without vigorous mixing,the effect of the amylase was restricted to a small local area, and for oursimulation we assume that the bolus was not vigorously mixed with salivabefore swallowing.The most significant source of error in the simulation is in the boundaryextraction from the CT images. Since the work was done manually, humanerrors are inevitable but hopefully small. Spatial resolution of the CT im-ages is very high and more than adequate for thin structures such as theepiglottis, however the temporal resolution, 10 Hz may not be enough tocapture sufficient temporal information. However, despite this shortcoming,the dynamic 320-row ADCT is the only imaging modality that allows forthis type of simulation study to be performed.4.6.5 Future directionsThis work can be extended in a number of directions. As a preliminaryexperiment, a number of issues were highlighted and these should be ad-dressed in subsequent work. Imaging dry-mouth, but otherwise normal sub-jects, would give a good test case for the no-slip boundaries. Increasing thenumber of subjects, as well as postures and bolus types, would improve thestrength of the findings. The study of solid and semi-solid boluses wouldrequire more complex material properties definitions, measurements, andsimulation techniques than those present in this study.Another possible extension is to modify the geometry and/or timings ofthe boundaries in a way that mimics physical pathologies, such as base-of-1024.7. Summarytongue tumour removal. This may help illuminate some of the underlyingcauses for dysphagia in patients who have had partial glossectomy with afreeflap and/or radiotherapy.The airway could be coupled to finite-element computer models of thehead and neck anatomy using software such as ArtiSynth [77], instead ofhaving its deformations kinematically controlled. This would allow us tostudy how changes to the organs may affect the bolus flow.4.7 SummaryThis chapter presents a numerical study of six healthy oropharyngeal swal-lows. We test the assumption that a slip boundary condition would providea better approximation to reality than a no-slip condition. The reason-ing behind the assumption is that in normal individuals, saliva coats theaerodigestive tract and has a much lower viscosity than thickened boluses.We found that the no-slip condition gave a better approximation to the CTimages when observing bolus transit times, however oropharyngeal residueafter the swallow was significantly increased.The results of the thin swallows were inconclusive. Because of the verylow viscosity of water, the slip and no-slip simulations were nearly identical,and did not always agree with the CT images. Nevertheless, because ofthe nearly identical behaviour between slip and no-slip simulations, theseresults would neither support or refute the assumption being tested. It isnot obvious what could have caused this result, but one possible explanationis that thin boluses move quicker than thickened boluses and the 10 Hz frame1034.7. Summaryrate of the CT images could not adequately capture the solid motion.In the existing literature, there is conflicting evidence with respect tothe effect of hyposalivation on dysphagia. Our observations corroboratewith modified barium studies that indicate bolus transit times are unaf-fected by hyposalivation, but residue increases with hyposalivation. Furtherinvestigation in the clinic, using simulation, and in the laboratory studyingsaliva at the microscopic level may shed light on this complex phenomenon.104Chapter 5ConclusionThis thesis presents a numerical study of fluid flow in the oropharynx, withthe goal of understanding the relationship between saliva and dysphagia. Inthe current literature, there are only a handful of publications that studyfluid flow in the oropharynx, and many of them used simplified geometriesof the airway due to a lack of imaging modalities that could support true 3Dfluid analysis. With the advent of dynamic full 3D imaging of human swal-lowing using 320-row ADCT, true dynamic 3D simulation is now possible.However, a number of challenges first needed to be overcome.The first challenge is deciding on an appropriate fluid simulation tech-nique for representing the bolus. 3D fluid simulation is a challenging prob-lem, and simulation of swallowing requires a transient, free-surface repre-sentation and proper handling of moving boundaries. The Unified Semi-Analytic Wall (USAW) SPH method offers an attractive solution to some ofthese problems. However, it required extension to be able to handle mov-ing 3D boundaries and non-Newtonian flows. With these modifications, itbecame suitable for the study of thickened liquids in the oropharynx. Weperform partial validation of the methods through the simulation of Couette,Hagen-Poiseuille, and lid-driven cavity flows.105Chapter 5. ConclusionExtracting the airway and its motion from the ADCT data required thedevelopment of a new tool, BlendSeg. Using BlendSeg, airways from a seriesof ADCT scans were extracted and used to drive the modified USAW SPHsimulations of the bolus. The nectar- and honey-thick boluses were thickenedwith starch. Their viscosity was measured using a KINEXUS rheometer.They are both shear thinning materials, and in the extended USAW SPHmethod, they were modelled as Cross-type fluids.BlendSeg was shown to be practical and effective as it was used to extractgeometry from six ADCT image sequences of healthy subjects swallowingthin, nectar-thick, and honey-thick boluses. We conducted a preliminary nu-merical investigation using these geometries and the extended USAW SPHmethod. For each of the six geometries, both no-slip and slip boundary con-ditions were simulated. Bolus transit times were measured and comparedto measurements of the ADCT data. The results for the thin bolus simu-lations were inconclusive. The thickened bolus simulations suggest that innormal subjects saliva does not speed up the bolus transit times. However,the no-slip simulations had higher amounts of residue after the swallow.Impaired saliva production is thought to cause dysphagia, however a re-view of the literature indicates that the relationship between them is unclear.A number of existing clinical studies have compared swallowing of subjectswith hyposalivation to normals, and found that bolus transit times wereunaffected in the patient group. Further, increased residue was observedin patient groups experiencing hyposalivation [43, 107], corroborating oursimulation findings. On the other hand, other studies find that swallowingdurations measured using anatomical landmarks increase for subjects with1065.1. Limitationshyposalivation [12, 105].To summarize, the contributions of this thesis are firstly presenting anextension to a state-of-the-art SPH simulation method, in order to run sim-ulations in the oropharynx with non-Newtonian liquids. A new method ofmanually extracting boundaries from 4D data was created and describedhere. This method was used to extract a boundary from dynamic ADCTdata and used to represent the aerodigestive tract in our swallowing simu-lations. These methods were used in a preliminary examination of saliva’srole in swallowing liquid boluses. The experiments suggest that saliva doesnot increase bolus speed, but may decrease oropharyngeal residue.5.1 LimitationsThe contributions presented in this thesis have a number of limitations.The tool developed to create moving boundaries from ADCT data, Blend-Seg, is a tool that allows an operator to make fine, detailed adjustmentsto a 3D mesh. Currently we use it as the primary tool in the creation ofdynamic 3D geometries. However, the process is time consuming and atrisk of error. The ADCT data has a relatively low temporal resolution of10 Hz. As a result of subject and bolus motion during the acquisition ofeach frame, artifacts such as blurring, ghosting, and streaking are present inthe images. These artifacts confound automated registration and segmen-tation techniques, and also manual extraction with BlendSeg as boundariesbecome difficult to delineate.Because the oropharyngeal bolus simulations are primarily driven by the1075.1. Limitationsboundary motion, any errors in airway identification because of operatormistakes, or due to artifacts in the images, culminate in errors of the fluidsimulation and may influence the results. We expect this to be the primarysource of error in the oropharyngeal bolus simulations. In fact, we suspectthis may contribute to the disagreement between thin bolus simulations andthe ADCT data. The biggest hurdle to improving the 3D bolus simulationsmay be obtaining 3D data with a higher frame rate than 10 Hz.The extended USAW SPH methods could be better validated. A strongervalidation would include moving boundaries, a free surface, and non-Newtonianfluid, with a size and fluid velocity comparable to that of oropharyngealswallowing. This type of experiment would be useful not only for validat-ing the extended USAW SPH method, but for measuring improvements toit, and also to validate other simulation methods as they become available.A significant issue with our extension is that the non-Newtonian solver isexplicit. For shear-thinning boluses, this puts a severe time step size restric-tion on the solver, limiting the resolution of the particles and increasing thesimulation time.Oropharyngeal swallowing is an extremely complex phenomenon. Ourtreatment of it was limited to the consideration of fluids. Additionally, onlythree healthy subjects with six swallows between them were studied. Anexpanded study with perhaps a dozen healthy subjects swallowing differentmaterials would result in stronger clinical impact of the findings.The inconsistent results of the thin bolus simulations were troubling.For the thickened boluses, a number of assumptions were made. The first isthat the air in the airway does not contribute significantly to the flow of the1085.1. Limitationsbolus. In terms of shearing, this is most likely a good assumption. However,without considering air explicitly, SPH can not predict the trapping of airbubbles and their effect on the flow.The temperature of the bolus was assumed to be constant at 10◦C. Thisis the temperature at which the rheology was measured. For all fluids, andparticularly for non-Newtonian fluids, the viscosity can be a strong functionof temperature. observation of the thickened boluses indicates that they aresensitive to temperature, but the sensitivity was not quantified.The presence of amylase in saliva could also have a large effect on thebolus viscosity. In a study of starch-thickened boluses, a 1 ml volume ofsaliva was introduced to 10 ml of starch-thickened bolus, the viscosity ofthe bolus was reduced by 90% after only 10 s of mixing [44]. However, thisfinding seems to be sensitive to mixing. If the subject mixed saliva withthe bolus vigorously before swallowing, it would likely have a large effecton the viscosity. In the ADCT study, this effect is hopefully small if thesubject holds the bolus stationary in their mouth before swallowing. It isalso possible that mixing with saliva occurs during the swallow, and thatthe bolus viscosity is decreasing during the oropharyngeal swallow. Oursimulations do not account for this factor.The resolution of the particles is another potential source of error inthe swallowing simulation. Each particle has a diameter of 1 mm, and eachmillilitre of bolus is represented by 1000 particles. The 4 ml nectar bolus wasre-run with twice the number of particles, each having a diameter of 0.8 mm,without any significant change in the bolus timing measures. However itremains a possibility that further increasing the number of particles could1095.2. Future workhave an effect on the results. For example, the upper esophageal sphincteris normally closed. This closure is represented in the model by opposingsurfaces overlapping. We observed that as the particles shrink in size, moredemand is placed on the airway model because smaller particles may slipthrough the closure if the regions do not overlap adequately. In addition, thesqueezing of the stripping wave is similarly modelled by the two opposingsurfaces (representing the tongue and the pharynx) pushing together. If thesqueezing does not proceed in the proper peristaltic direction, instabilitycould occur as fluid is trapped with nowhere to go. Therefore, increasingthe resolution of the simulation significantly might require a faster solver,increased temporal resolution to properly resolve the stripping wave, andimproved techniques for representing closure.5.2 Future workLooking forward, a number of future directions present themselves in thefields of fluid simulation, biomechanics, dysphagia research, and image pro-cessing. One of the most time consuming parts of this work was extractingmoving boundaries from ADCT data. Finding new automated techniquesthat could speed up this process would significantly reduce the time it takesto study the CT data. Further improvements in the temporal resolution ofADCT data could possibly make it easier for automated methods to work,as this might decrease the severity of motion artifacts.In our extension of USAW SPH, shear-thinning boluses integrated ex-plicitly require very small time steps and increases the computation time1105.2. Future workgreatly. This could be alleviated by performing an implicit time integrationfor the viscous term, however existing implicit viscosity solvers would needto be extended to the USAW SPH formulation.The current oropharyngeal swallowing study only included two healthysubjects swallowing liquids. More subjects of different ages and patholo-gies should be included in order to have a more complete picture of salivain swallowing. Another approach forward is to gather 3D scans of sub-jects swallowing semi-solid and solid boluses. A large gap in knowledgeexists here because of the complexity in the physics behind mastication,food breakdown, and swallowing of solid foods.Moving away from the bolus, existing biomechanics models of the humanhead and neck region can now leverage the dynamic 3D ADCT scans. Theimaging data provides a detailed look at the 3D motion in the oropharynxand can potentially be used to study the biomechanics of the head and neckduring swallowing. A biomechanics simulation coupled to a fluid simulationof the bolus has the capability of predicting swallowing outcomes that couldresult from surgical or neurological changes to the body. The work presentedin [116] illustrates a possible approach to achieving this coupling by using aunified skinning technique.In dysphagia research, the observations of the simulation should be fur-ther corroborated with clinical studies. These could be done by re-analyzingolder data, or by designing new studies that complement the simulation data.Saliva is a complex fluid and so far, most dysphagia research has focused onquantity of saliva. Studying the molecular changes in saliva due to medica-tions and radiotherapy, and its effect on swallowing, may be a fruitful area1115.2. Future workof research.112Bibliography[1] Rebecca H Affoo, Norine Foley, Rushlee Garrick, Walter L Siqueira,and Ruth E Martin. Meta-analysis of salivary flow rates in young andolder adults. Journal of the American Geriatrics Society, 63(10):2142–2151, 2015.[2] Sug-Joon Ahn, Ling Tsou, C Antonio Sa´nchez, Sidney Fels, and Ho-Beom Kwon. Analyzing center of rotation during opening and closingmovements of the mandible using computer simulations. Journal ofbiomechanics, 48(4):666–671, 2015.[3] Matteo Antuono, Andrea Colagrossi, Salvatore Marrone, and DiegoMolteni. Free-surface flows solved by means of sph schemes withnumerical diffusive terms. Computer Physics Communications,181(3):532–549, 2010.[4] David Baraff, Andrew Witkin, and Michael Kass. Untangling cloth.In ACM Transactions on Graphics (TOG), volume 22, pages 862–870.ACM, 2003.[5] Christopher Batty, Florence Bertails, and Robert Bridson. A fast vari-113Bibliographyational framework for accurate solid-fluid coupling. In ACM Transac-tions on Graphics (TOG), volume 26, page 100. ACM, 2007.[6] Ted Belytschko, Yury Krongauz, Daniel Organ, Mark Fleming, andPetr Krysl. Meshless methods: an overview and recent developments.Computer methods in applied mechanics and engineering, 139(1):3–47,1996.[7] Ted Belytschko, Yun Yun Lu, and Lei Gu. Element-free galerkinmethods. International journal for numerical methods in engineer-ing, 37(2):229–256, 1994.[8] M Bergdahl and J Bergdahl. Low unstimulated salivary flow and sub-jective oral dryness: association with medication, anxiety, depression,and stress. Journal of Dental Research, 79(9):1652–1658, 2000.[9] Richard P Beyer. A computational model of the cochlea using theimmersed boundary method. Journal of Computational Physics,98(1):145–162, 1992.[10] J Bonet, S Kulasegaram, MX Rodriguez-Paz, and M Profit. Varia-tional formulation for the smooth particle hydrodynamics (sph) sim-ulation of fluid and solid problems. Computer Methods in AppliedMechanics and Engineering, 193(12):1245–1256, 2004.[11] JU Brackbill and HM Ruppel. Flip: A method for adaptively zoned,particle-in-cell calculations of fluid flows in two dimensions. Journalof Computational Physics, 65(2):314–343, 1986.114Bibliography[12] Anthony J Caruso, Barbara C Sonies, Jane C Atkinson, and Philip CFox. Objective measures of swallowing in patients with primarysjo¨gren’s syndrome. Dysphagia, 4(2):101–105, 1989.[13] Thanapong Chaichana, Zhonghua Sun, and James Jewkes. Computa-tion of hemodynamics in the left coronary artery with variable angu-lations. Journal of biomechanics, 44(10):1869–1878, 2011.[14] Michael W. Chang, Eugene Lin, and Jenq-Neng Hwang. Contourtracking using a knowledge-based snake algorithm to construct three-dimensional pharyngeal bolus movement. Dysphagia, 14(4):219–227,1999.[15] Michael W. Chang, Brigette Rosendall, and Bruce A. Finlayson.Mathematical modeling of normal pharyngeal bolus transport; a pre-liminary study. Journal of rehabilitation research and development,35:327–334, 1998.[16] Med Chaos. Normal epiglottis, 2013.[17] Gloria Chi-Fishman. Quantitative lingual, pharyngeal and laryngealultrasonography in swallowing research: a technical review. Clinicallinguistics & phonetics, 19(6-7):589–604, 2005.[18] PW Cleary, J Ha, M Prakash, and T Nguyen. 3d sph flow predictionsand validation for high pressure die casting of automotive components.Applied Mathematical Modelling, 30(11):1406–1427, 2006.[19] Andrea Colagrossi and Maurizio Landrini. Numerical simulation of in-115Bibliographyterfacial flows by smoothed particle hydrodynamics. Journal of Com-putational Physics, Volume 191, Issue 2, p. 448-475., 191:448–475,nov 2003.[20] LMC Collins and C Dawes. The surface area of the adult humanmouth and thickness of the salivary film covering the teeth and oralmucosa. Journal of dental research, 66(8):1300–1302, 1987.[21] Ian J Cook, Wylie J Dodds, RO Dantas, Benson Massey, Mark K Kern,Ivan M Lang, James G Brasseur, and WJ Hogan. Opening mecha-nisms of the human upper esophageal sphincter. American Journal ofPhysiology-Gastrointestinal and Liver Physiology, 257(5):G748–G759,1989.[22] Sharen J. Cummins and Murray Rudman. An sph projection method.Journal of computational physics, 152(2):584–607, 1999.[23] Kristopher S Cunningham and Avrum I Gotlieb. The role of shearstress in the pathogenesis of atherosclerosis. Laboratory investigation,85(1):9–23, 2005.[24] Je´roˆme Declerck, Jacques Feldmar, and Nicholas Ayache. Definitionof a four-dimensional continuous planispheric transformation for thetracking and the analysis of left-ventricle motion. Medical Image Anal-ysis, 2(2):197–213, 1998.[25] W.J. Dodds. The physiology of swallowing. Dysphagia, 3(4):171–178,1989.116Bibliography[26] J Donea, Antonio Huerta, J-Ph Ponthot, and A Rodriguez-Ferran. En-cyclopedia of computational mechanics vol. 1: Fundamentals., chapter14: Arbitrary lagrangian-eulerian methods, 2004.[27] Jim Douglas and James E Gunn. A general formulation of alternatingdirection methods. Numerische Mathematik, 6(1):428–453, 1964.[28] John W Eveson. Xerostomia. Periodontology 2000, 48(1):85–91, 2008.[29] EA Fadlun, R Verzicco, Paolo Orlandi, and J Mohd-Yusof. Combinedimmersed-boundary finite-difference methods for three-dimensionalcomplex flow simulations. Journal of computational physics,161(1):35–60, 2000.[30] X-J Fan, RI Tanner, and R Zheng. Smoothed particle hydrody-namics simulation of non-newtonian moulding flow. Journal of Non-Newtonian Fluid Mechanics, 165(5):219–226, 2010.[31] Stanley Osher Ronald Fedkiw and S Osher. Level set methods anddynamic implicit surfaces. Surfaces, 44:77, 2002.[32] Martin Ferrand, D. R. Laurence, B. D. Rogers, Damien Violeau, andChristophe Kassiotis. Unified semi-analytical wall boundary condi-tions for inviscid, laminar or turbulent flows in the meshless sphmethod. International Journal for Numerical Methods in Fluids,71(4):446–472, 2013.[33] Cormac Flynn, Ian Stavness, John Lloyd, and Sidney Fels. A finiteelement model of the face including an orthotropic skin model under117Bibliographyin vivo tension. Computer Methods in Biomechanics and BiomedicalEngineering, 0(0):1–12, 0. PMID: 23919890.[34] Philip C Fox, Peter F van der Ven, Barbara C Sonies, James M Weif-fenbach, and Bruce J Baum. Xerostomia: evaluation of a symptomwith increasing significance. The Journal of the American Dental As-sociation, 110(4):519–525, 1985.[35] Robert I Fox, Michael Stern, and Paul Michelson. Update in sjo¨grensyndrome. Current opinion in rheumatology, 12(5):391–398, 2000.[36] Naoko Fujii, Yoko Inamoto, Eiichi Saitoh, Mikoto Baba, SumikoOkada, Satoshi Yoshioka, Toshiaki Nakai, Yoshihiro Ida, KazuhiroKatada, and Jeffrey B Palmer. Evaluation of swallowing using 320-detector-row multislice ct. part i: single- and multiphase volume scan-ning for three-dimensional morphological and kinematic analysis. Dys-phagia, 26(2):99–107, 2011.[37] P Garcia-Peris, L Paro´n, C Velasco, C De la Cuerda, M Camblor,I Breto´n, H Herencia, J Verdaguer, C Navarro, and P Clave. Long-term prevalence of oropharyngeal dysphagia in head and neck cancerpatients: impact on quality of life. Clinical Nutrition, 26(6):710–717,2007.[38] Donna T Geddes, Lynda M Chadwick, Jacqueline C Kent, Cather-ine P Garbin, and Peter E Hartmann. Ultrasound imaging of infantswallowing during breast-feeding. Dysphagia, 25(3):183–191, 2010.118Bibliography[39] U. K. N. G. Ghia, Kirti N. Ghia, and C. T. Shin. High-re solutions forincompressible flow using the navier-stokes equations and a multigridmethod. Journal of computational physics, 48(3):387–411, 1982.[40] Robert A Gingold and Joseph J Monaghan. Smoothed particle hy-drodynamics: theory and application to non-spherical stars. Monthlynotices of the royal astronomical society, 181(3):375–389, 1977.[41] Henry Gray. Anatomy of the human body. Lea & Febiger, 1918.[42] S. Hamlet, L. Jones, R. Mathog, M. Bolton, and R. Patterson. Boluspropulsive activity of the tongue in dysphagic cancer patients. Dys-phagia, 3(1):18–23, 1988.[43] Sandra Hamlet, Jennifer Faull, Barbara Klein, Amr Aref, JamesFontanesi, Robert Stachler, Falah Shamsa, Lewis Jones, and MarkSimpson. Mastication and swallowing in patients with postirradiationxerostomia. International Journal of Radiation Oncology* Biology*Physics, 37(4):789–796, 1997.[44] Ben Hanson, Mark T OLeary, and Christina H Smith. The effect ofsaliva on the viscosity of thickened drinks. Dysphagia, 27(1):10–19,2012.[45] Negar M Harandi, Jonghye Woo, Maureen Stone, Rafeef Abughar-bieh, and Sidney Fels. Subject-specific biomechanical modelling of thetongue: analysis of muscle activations during speech. Proc Int SeminSpeech Prod (ISSP), pages 174–177, 2014.119Bibliography[46] Francis H. Harlow and J. Eddie Welch. Numerical calculation oftime-dependent viscous incompressible flow of fluid with free surface.Physics of Fluids, Volume 8, Issue 12, p.2182-2189, 8:2182–2189, dec1965.[47] H Ric Harnsberger. Handbook of head and neck imaging. Mosby In-corporated, 1995.[48] Simon M Harrison, Paul W Cleary, Graham Eyres, Matthew D Sin-nott, and Leif Lundin. Challenges in computational modelling of foodbreakdown and flavour release. Food & function, 5(11):2792–2805,2014.[49] Simon M. Harrison, Graham Eyres, Paul W. Cleary, Matthew D. Sin-nott, Conor Delahunty, and Leif Lundin. Computational modeling offood oral breakdown using smoothed particle hydrodynamics. Journalof Texture Studies, 45(2):97–109, 2014.[50] Hellerhoff. Normaler schluck-11, 2011.[51] BS Henson, MR Inglehart, A Eisbruch, and JA Ship. Preserved sali-vary output and xerostomia-related quality of life in head and neckcancer patients receiving parotid-sparing radiotherapy. Oral oncology,37(1):84–93, 2001.[52] Cyril W. Hirt and Billy D. Nichols. Volume of fluid (vof) method forthe dynamics of free boundaries. Journal of computational physics,39(1):201–225, 1981.120Bibliography[53] Susan G Hiss and Gregory N Postma. Fiberoptic endoscopic evaluationof swallowing. The Laryngoscope, 113(8):1386–1393, 2003.[54] Andrew Kenneth Ho, Yoko Inamoto, Eiichi Saitoh, Sheldon Green,and Sidney Fels. Extracting moving boundaries from dynamic, multi-slice ct images for fluid simulation. Computer Methods in Biomechan-ics and Biomedical Engineering: Imaging & Visualization, pages 1–6,2016.[55] Andrew Kenneth Ho, Ling Tsou, Sheldon Green, and Sidney Fels. A 3dswallowing simulation using smoothed particle hydrodynamics. Com-puter Methods in Biomechanics and Biomedical Engineering: Imaging& Visualization, 2(4):237–244, 2014.[56] SM Hosseini, MT Manzari, and SK Hannani. A fully explicit three-stepsph algorithm for simulation of non-newtonian fluid flow. InternationalJournal of Numerical Methods for Heat & Fluid Flow, 17(7):715–735,2007.[57] Yoko Inamoto, Naoko Fujii, Eiichi Saitoh, Mikoto Baba, SumikoOkada, Kazuhiro Katada, Yasunori Ozeki, Daisuke Kanamori, andJeffrey B Palmer. Evaluation of swallowing using 320-detector-rowmultislice ct. part ii: kinematic analysis of laryngeal closure duringnormal swallowing. Dysphagia, 26(3):209–17, 2011.[58] Razvan Ioan Ionasec, Ingmar Voigt, Bogdan Georgescu, Yang Wang,Helene Houle, Fernando Vega-Higuera, Nassir Navab, and Dorin Co-maniciu. Patient-specific modeling and quantification of the aortic121Bibliographyand mitral valves from 4-d cardiac ct and tee. IEEE transactions onmedical imaging, 29(9):1636–1651, 2010.[59] Barbara M Johnston, Peter R Johnston, Stuart Corney, and DavidKilpatrick. Non-newtonian blood flow in human right coronary arter-ies: steady state simulations. Journal of biomechanics, 37(5):709–720,2004.[60] Barbara M Johnston, Peter R Johnston, Stuart Corney, and DavidKilpatrick. Non-newtonian blood flow in human right coronary arter-ies: transient simulations. Journal of biomechanics, 39(6):1116–1128,2006.[61] Peter J Kahrilas, Shezhang Lin, Jim Chen, and Jerilyn A Logemann.Three-dimensional modeling of the oropharynx during swallowing. Ra-diology, 194(2):575–579, 1995.[62] Peter J Kahrilas, Shezhang Lin, Jerilyn A Logemann, Gulchin A Er-gun, and Frank Facchini. Deglutitive tongue action: volume accommo-dation and bolus propulsion. Gastroenterology, 104(1):152–162, 1993.[63] Takahiro Kikuchi, Yukihiro Michiwaki, Tetsu Kamiya, Yoshio Toyama,Tasuku Tamai, and Seiichi Koshizuka. Human swallowing simulationbased on videofluorography images using hamiltonian mps method.Computational Particle Mechanics, 2(3):247–260, 2015.[64] John Kim and Parviz Moin. Application of a fractional-step methodto incompressible navier-stokes equations. Journal of computationalphysics, 59(2):308–323, 1985.122Bibliography[65] Jungwoo Kim, Dongjoo Kim, and Haecheon Choi. An immersed-boundary finite-volume method for simulations of flow in complex ge-ometries. Journal of Computational Physics, Volume 171, Issue 1, pp.132-150 (2001)., 171:132–150, jul 2001.[66] Stefan Klein, Marius Staring, Keelin Murphy, Max A Viergever, andJosien PW Pluim. Elastix: a toolbox for intensity-based medical imageregistration. Medical Imaging, IEEE Transactions on, 29(1):196–205,2010.[67] Seiichi Koshizuka, Atsushi Nobe, and Yoshiaki Oka. Numerical analy-sis of breaking waves using the moving particle semi-implicit method.International Journal for Numerical Methods in Fluids, 26(7):751–769,1998.[68] Seiichi Koshizuka and Y Oka. Moving-particle semi-implicit methodfor fragmentation of incompressible fluid. Nuclear science and engi-neering, 123(3):421–434, 1996.[69] Sivakumar Kulasegaram, Javier Bonet, RW Lewis, and M Profit.A variational formulation based contact algorithm for rigid bound-aries in two-dimensional sph applications. Computational Mechanics,33(4):316–325, 2004.[70] William E Langlois and Michel O Deville. Slow viscous flow. Springer,1964.[71] E.-S. Lee, C. Moulinec, R. Xu, D. Violeau, D. Laurence, andP. Stansby. Comparisons of weakly compressible and truly incom-123Bibliographypressible algorithms for the sph mesh free particle method. Jour-nal of Computational Physics, Volume 227, Issue 18, p. 8417-8436.,227:8417–8436, sep 2008.[72] Agnes Leroy, Damien Violeau, Martin Ferrand, and Christophe Kas-siotis. Unified semi-analytical wall boundary conditions applied to 2-dincompressible sph. Journal of Computational Physics, 261:106–129,2014.[73] Shezhang Lin, Jim Chen, Paul Hertz, and Peter J. Kahrilas. Dynamicreconstruction of the orophanryngeal swallow using computer basedanimation. Computerized medical imaging and graphics, 20(2):69–75,1996.[74] Xiang Lin, Brett R Cowan, and Alistair A Young. Automated de-tection of left ventricle in 4d mr images: experience from a largestudy. In International Conference on Medical Image Computing andComputer-Assisted Intervention, pages 728–735. Springer, 2006.[75] G.R. Liu and MB Liu. Smoothed particle hydrodynamics: a meshfreeparticle method. World Scientific Pub Co Inc, 2003.[76] M. B. Liu and G. R. Liu. Smoothed particle hydrodynamics (sph): anoverview and recent developments. Archives of computational methodsin engineering, 17(1):25–76, 2010.[77] John E. Lloyd, Ian Stavness, and Sidney Fels. Artisynth: a fast inter-active biomechanical modeling toolkit combining multibody and finite124Bibliographyelement simulation. In Soft tissue biomechanical modeling for com-puter assisted surgery, pages 355–394. Springer, 2012.[78] Jeri A Logemann, PETER J Kahrilas, JOAN Cheng, BR Pauloski,PATRICIA J Gibbons, ALFRED W Rademaker, and SHEZHANGLin. Closure mechanisms of laryngeal vestibule during swallow.American Journal of Physiology-Gastrointestinal and Liver Physiol-ogy, 262(2):G338–G344, 1992.[79] Vincent Luboz, Antoine Perrier, Ian Stavness, JE Lloyd, Marek Bucki,Francis Cannard, Bruno Diot, Nicolas Vuillerme, and Yohan Payan.Foot ulcer prevention using biomechanical modelling. Computer Meth-ods in Biomechanics and Biomedical Engineering: Imaging & Visual-ization, 2(4):189–196, 2014.[80] Leon B Lucy. A numerical approach to the testing of the fission hy-pothesis. The astronomical journal, 82:1013–1024, 1977.[81] Michael Lynch, Ovidiu Ghita, and Paul F Whelan. Segmentation ofthe left ventricle of the heart in 3-d+ t mri data using an optimizednonrigid temporal model. IEEE Transactions on Medical Imaging,27(2):195–203, 2008.[82] M Malve`, A Garcia, J Ohayon, and MA Martinez. Unsteady bloodflow and mass transfer of a human left coronary artery bifurcation:Fsi vs. cfd. International communications in heat and mass transfer,39(6):745–751, 2012.125Bibliography[83] S Marrone, M Antuono, A Colagrossi, G Colicchio, D Le Touze´, andG Graziani. δ-sph model for simulating violent impact flows. ComputerMethods in Applied Mechanics and Engineering, 200(13):1526–1542,2011.[84] Bonnie Martin-Harris, Martin B Brodsky, Yvonne Michel, Donald OCastell, Melanie Schleicher, John Sandidge, Rebekah Maxwell, andJulie Blair. Mbs measurement tool for swallow impairment–mbsimp:establishing a standard. Dysphagia, 23(4):392–405, 2008.[85] Arno Mayrhofer, Martin Ferrand, Christophe Kassiotis, Damien Vio-leau, and Franc¸ois-Xavier Morel. Unified semi-analytical wall bound-ary conditions in sph: analytical extension to 3-d. Numerical Algo-rithms, 68(1):15–34, 2015.[86] Tim McInerney and Demetri Terzopoulos. A dynamic finite elementsurface model for segmentation and tracking in multidimensional med-ical images with application to cardiac 4d image analysis. Computer-ized Medical Imaging and Graphics, 19(1):69–83, 1995.[87] Y. Meng, M. A. Rao, and A. K. Datta. Computer simulation of thepharyngeal bolus transport of newtonian and non-newtonian fluids.Food and bioproducts processing, 83(4):297–305, 2005.[88] Rajat Mittal and Gianluca Iaccarino. Immersed boundary methods.Annual Review of Fluid Mechanics, vol. 37, Issue 01, pp.239-261,37:239–261, jan 2005.126Bibliography[89] H. Mizunuma, M. Sonomura, K. Shimokasa, H. Ogoshi, S. Nakamura,and N. Tayama. Numerical modeling and simulation on the swallowingof jelly. Journal of texture studies, 40(4):406–426, 2009.[90] J Mohd-Yusof. Combined immersed-boundary/b-spline methods forsimulations of ow in complex geometries. Annual Research Briefs.NASA Ames Research Center= Stanford University Center of Turbu-lence Research: Stanford, pages 317–327, 1997.[91] Joe J. Monaghan. Smoothed particle hydrodynamics. Annual reviewof astronomy and astrophysics, 30:543–574, 1992.[92] Joe J. Monaghan. Simulating free surface flows with sph. Journal ofcomputational physics, 110(2):399–406, 1994.[93] Johan Montagnat and Herve´ Delingette. 4d deformable models withtemporal constraints: application to 4d cardiac image segmentation.Med Image Anal, 9(1):87–100, 2005.[94] Joseph P. Morris. Simulating surface tension with smoothed particlehydrodynamics. International journal for numerical methods in fluids,33(3):333–353, 2000.[95] Joseph P. Morris, Patrick J. Fox, and Yi Zhu. Modeling low reynoldsnumber incompressible flows using sph. Journal of computationalphysics, 136(1):214–226, 1997.[96] M.A. Nicosia, J.A. Hind, E.B. Roecker, M. Carnes, J. Doyle, G.A.Dengel, and J.A. Robbins. Age effects on the temporal evolution of127Bibliographyisometric and swallowing pressure. The Journals of Gerontology SeriesA: Biological Sciences and Medical Sciences, 55(11):M634, 2000.[97] Mark A Nicosia. A planar finite element model of bolus containmentin the oral cavity. Comput. Biol. Med., 37(10):1472–8, 2007.[98] Mark A. Nicosia. Theoretical estimation of shear rate during the oralphase of swallowing: effect of partial slip. Journal of Texture Studies,44(2):132–139, 2013.[99] Mark A. Nicosia and JoAnne Robbins. The fluid mechanics of bolusejection from the oral cavity. Journal of Biomechanics, 34(12):1537–1544, 2001.[100] Kerstin E O¨hrn, Ylva-Britt Wahlin, and Per-Olow Sjo¨de´n. Oral statusduring radiotherapy and chemotherapy: a descriptive study of patientexperiences and the occurrence of oral complications. Supportive carein cancer, 9(4):247–257, 2001.[101] Barbara Roa Pauloski, Alfred W Rademaker, Jerilyn A Logemann,Muveddet Discekici-Harris, and Bharat B Mittal. Comparison ofswallowing function after intensity-modulated radiation therapy andconventional radiotherapy for head and neck cancer. Head & neck,37(11):1575–1582, 2015.[102] Charles S Peskin. The immersed boundary method. Acta numerica,11:479–517, 2002.[103] Philippe Poisson, Thibault Laffond, Sandra Campos, Veronique128BibliographyDupuis, and Isabelle Bourdel-Marchasson. Relationships between oralhealth, dysphagia and undernutrition in hospitalised elderly patients.Gerodontology, 2014.[104] Panu JF Rantonen and Jukka H Meurman. Viscosity of whole saliva.Acta Odontologica Scandinavica, 56(4):210–214, 1998.[105] Nelson L Rhodus, Stephen Colby, Karlind Moller, and Janna Bereuter.Quantitative assessment of dysphagia in patients with primary andsecondary sjo¨gren’s syndrome. Oral Surgery, Oral Medicine, OralPathology, Oral Radiology, and Endodontology, 79(3):305–310, 1995.[106] Nicole M Rogus-Pulia, Charles Larson, Bharat B Mittal, Marge Pierce,Steven Zecker, Korey Kennelty, Amy Kind, and Nadine P Connor. Ef-fects of change in tongue pressure and salivary flow rate on swallowefficiency following chemoradiation treatment for head and neck can-cer. Dysphagia, pages 1–10, 2016.[107] Nicole M Rogus-Pulia and Jeri A Logemann. Effects of reduced salivaproduction on swallowing in patients with sjo¨grens syndrome. Dys-phagia, 26(3):295–303, 2011.[108] James Albert Sethian. Level set methods and fast marching methods:evolving interfaces in computational geometry, fluid mechanics, com-puter vision, and materials science, volume 3. Cambridge universitypress, 1999.[109] Songdong Shao and Edmond Y. M. Lo. Incompressible sph methodfor simulating newtonian and non-newtonian flows with a free sur-129Bibliographyface. Advances in Water Resources, Volume 26, Issue 7, p. 787-800.,26:787–800, n/a 2003.[110] Leonardo Di G. Sigalotti, Jaime Klapp, Eloy Sira, Yasmin Melea´n,and Anwar Hasmy. Sph simulations of time-dependent poiseuille flowat low reynolds numbers. Journal of Computational Physics, Volume191, Issue 2, p. 622-638., 191:622–638, nov 2003.[111] BC Sonies, JA Ship, and BJ Baum. Relationship between saliva pro-duction and oropharyngeal swallow in healthy, different-aged adults.Dysphagia, 4(2):85–89, 1989.[112] M. Sonomura, H. Mizunuma, T. Numamori, H. Michiwaki, andK. Nishinari. Numerical simulation of the swallowing of liquid bo-lus. Journal of Texture Studies, 42(3):203–211, 2011.[113] Leo M Sreebny and Steven S Schwartz. A reference guide to drugsand dry mouth–2nd edition. Gerodontology, 14(1):33–47, 1997.[114] Leo M Sreebny and Arjan Vissink. Dry mouth, the malevolent symp-tom: a clinical guide. John Wiley & Sons, 2010.[115] Ian Stavness, John E Lloyd, and Sidney Fels. Automatic predictionof tongue muscle activations using a finite element model. Journal ofbiomechanics, 45(16):2841–2848, 2012.[116] Ian Stavness, C Antonio Sa´nchez, John Lloyd, Andrew Ho, JohntyWang, Sidney Fels, and Danny Huang. Unified skinning of rigid and130Bibliographydeformable models for anatomical simulations. In SIGGRAPH Asia2014 Technical Briefs, page 9. ACM, 2014.[117] H. S. Udaykumar, R. Mittal, P. Rampunggoon, and A. Khanna. Asharp interface cartesian grid method for simulating flows with com-plex moving boundaries. Journal of Computational Physics, Volume174, Issue 1, pp. 345-380 (2001)., 174:345–380, nov 2001.[118] Martin Uecker, Shuo Zhang, Dirk Voit, Alexander Karaus, Klaus-Dietmar Merboldt, and Jens Frahm. Real-time mri at a resolutionof 20 ms. NMR Biomed, 23(8):986–94, 2010.[119] Tom Vercauteren, Xavier Pennec, Aymeric Perchant, and NicholasAyache. Diffeomorphic demons: Efficient non-parametric image regis-tration. NeuroImage, 45(1):S61–S72, 2009.[120] Holger Wendland. Piecewise polynomial, positive definite and com-pactly supported radial functions of minimal degree. Advances incomputational Mathematics, 4(1):389–396, 1995.[121] Frank M White and Isla Corfield. Viscous fluid flow, volume 3.McGraw-Hill New York, 2006.[122] Xiaoyang Xu, Jie Ouyang, Binxin Yang, and Zhijun Liu. Sph simula-tions of three-dimensional non-newtonian free surface flows. ComputerMethods in Applied Mechanics and Engineering, 256:101–116, 2013.[123] Shuo Zhang, Arno Olthoff, and Jens Frahm. Real-time magnetic reso-131nance imaging of normal swallowing. Journal of Magnetic ResonanceImaging, 35(6):1372–1379, 2012.[124] Luoding Zhu and Charles S Peskin. Simulation of a flapping flexiblefilament in a flowing soap film by the immersed boundary method.Journal of Computational Physics, 179(2):452–468, 2002.[125] Yongning Zhu and Robert Bridson. Animating sand as a fluid. InACM Transactions on Graphics (TOG), volume 24, pages 965–972.ACM, 2005.132Appendix AEstimate of effect of saliva ona gravity driven bolusTo estimate the influence of a thin layer of lubrication on a gravity drivenflow, we assume a geometry of a flat wall, perpendicular to the ground. Twolayers of fluid exist between the wall and the air: a higher viscosity fluid thatis closer to the air, and a lower viscosity lubricating layer that is betweenthe first fluid and the wall (see Figure A.1).The lubricating layer has a viscosity µs, and thickness Ls, while thethicker fluid viscosity µb and thickness Lb. Gravity drives the flow down-wards in the negative y direction. The density of the two fluids is assumedto be equal, ρ = ρs = ρb.For a fluid element as shown in Figure A.2, τw is the stress from thewall, and τ(x) is the stress at distance x from the wall. Setting up the forcebalance equation, ∑Fy = 0, (A.1)HDτw −HDτ(x) = Mg = ρgHDx, (A.2)133Appendix A. Estimate of effect of saliva on a gravity driven bolusyxgsalivaLsbolusLbairFigure A.1: Simplified geometry of a gravity driven bolus with lubricatingsaliva layer. Wall is on the left, gray arrows represent the velocity of thebolus. The g arrow shows the direction of gravity. Ls and Lb are the widthsof the saliva and bolus, resp.xτw τ(x)MgHFigure A.2: Fluid element (dashed box) under consideration with height Hdepth D into the page. τw is the shear stress at the wall, and τ(x) is theshear stress at position x. Mg is the force of gravity on the element.134Appendix A. Estimate of effect of saliva on a gravity driven boluswith boundary conditionsHDτw = ρgHD(Ls + Lb), (A.3)where τ(Ls + Lb) = 0 because the shear stress due to air is negligible.ThereforeρgHD(Ls + Lb)−HDτ(x) = ρgHDx, (A.4)τ(x) = ρg(Ls + Lb − x). (A.5)Inside the saliva, for a velocity V (x),τ(x) = µsdVdx, (A.6)dVdx=ρgµs(Ls + Lb − x), (A.7)V (x) =ρgµs((Ls + Lb)x− x22)+ C1, (A.8)But V (0) = 0 due to the no-slip condition, therefore C1 = 0V (x) =ρgµs((Ls + Lb)x− x22), x ≤ Ls, (A.9)135Appendix A. Estimate of effect of saliva on a gravity driven bolusandV (Ls) =ρgµs(L2s2+ LsLb). (A.10)In the bolus,dVdx=ρgµb(Ls + Lb − x), (A.11)V (x) =ρgµb((Ls + Lb)x− x22)+ C2, (A.12)and using the condition that the velocity at Ls is continuous,C2 = (ρgµs− ρgµb)(L2s2+ LsLb)=ρgµsµb(µb − µs)(L2s2+ LsLb). (A.13)The volumetric flow rate Q is the amount of fluid flowing through a fixedplane, and is defined for the bolus asQ =ˆ Ls+LbLsV (x) ·Ddx. (A.14)Combining equations A.12−A.14,Q = Dˆ Ls+LbLsρgµb[(Ls + Lb)x− x22]+ρgµsµb(µb − µs)(L2s2+ LsLb)dx,(A.15)Q =ρgDµb[(Ls + Lb)x22− x36]Ls+LbLs+ρgD(µb − µs)Lbµsµb(L2s2+ LsLb),(A.16)136Appendix A. Estimate of effect of saliva on a gravity driven bolusQ =ρgD6µb(3L2sLb+6LsL2b +2L3b)+ρgD(µb − µs)Lbµsµb(L2s2+ LsLb). (A.17)Notice that the first term for Q in equation A.17 is independent of theviscosity of saliva, while the second term does depend on it. Taking theratio of the second term to the first term and then simplifying gives us anestimate, which we call Rs, of how much saliva contributes to the overallbolus flow rate, assuming steady state.Rs =3(µb − µs)µsL2s + 2LsLb3L2s + 6LsLb + 2L2b. (A.18)Using equation A.18 and assuming the viscosity µb = 200µs, for Rs tohave a value less than 0.1 (a 10% increase in bolus flow rate due to saliva),Lb ≈ 6000 Ls.137Appendix BUsing full-slip toapproximate a lubricativesaliva layerIn this appendix, we show that a full-slip boundary condition can be used toapproximate viscous bolus flow with a lubricative saliva layer. Dependingon the choices of saliva viscosity and thickness, a full-slip condition gives areasonable approximation for time durations similar to those encountered inswallowing.Oropharyngeal swallowing occurs in less than one second, and the timerequired for the bolus to leave the oral cavity is much shorter. A time-dependent finite volume (FV) solution of the velocity and flow-rates fora two material, gravity driven flow in a pipe (Figure B.1) is described andanalyzed. The solver is verified by comparison with the theoretical transientand steady-state solutions, and shown to have second-order convergence inspace.The geometry is axisymmetric with two Newtonian materials represent-ing the bolus and the saliva. An air core can be included in the model and138Appendix B. Using full-slip to approximate a lubricative saliva layerzrgsalivaRsbolusRbairRaFigure B.1: Simplified geometry of a gravity driven bolus with lubricatingsaliva layer. Wall is on the right, gray arrows represent the fluid velocity.The g arrow shows the direction of gravity. Ra, Rb and Rs are the widthsof the air core, bolus, and saliva, resp. Axial symmetry is assumed, withthe centre of the air core at r = 0. The geometry has a total radius ofR = Ra +Rb +Rs.139B.1. Problem description and assumptionsdoes not exert shear forces on the bolus. All fluids are assumed to be ini-tially stationary and driven by gravity. A transient analysis shows that forthe time durations associated with swallowing, and for realistic choices ofsaliva radius and viscosity, the full-slip boundary condition gives a reason-able approximation to the expected flow rate.B.1 Problem description and assumptionsA layer of saliva is sandwiched between a solid stationary wall and a layerof bolus. The saliva and bolus are assumed to have Newtonian viscosities µsand µb, with radii Rs and Rb. They are assumed to have equal density, ρ.An air core is at the centre of the pipe, with radius Ra, and does not exertany shear forces on the bolus. The total pipe radius is R = Ra+Rb+Rs. Aconstant pressure gradient drives the flow with acceleration g. At time zerothe velocity is zero everywhere.We assume the flow is axisymmetric and infinite in the z direction. Theseassumptions simplify the Navier-Stokes equations in fluid m (where m couldbe bolus or saliva) to∂uy∂t=µmρm1r∂∂r(r∂uy∂r)+ g. (B.1)The velocity between saliva and wall is = 0. (B.2)140B.1. Problem description and assumptionsThe air exerts no shear stress on the fluid∂uy∂r(Ra) = 0. (B.3)The velocity of the saliva and bolus are equal at the saliva-bolus interface,uy(Ra +Rb)− = uy(Ra +Rb)+, (B.4)as are their stresses:τs(Ra +Rb) = τb(Ra +Rb) (B.5)µs∂uy∂r(Ra +Rb)− = µb∂uy∂r(Ra +Rb)+ (B.6)The + and − superscripts indicate a small distance to the left and right ofRx, resp.B.1.1 Steady-state solutionThe theoretical steady-state solution can be derived as follows,∂uy∂t=µmρm1r∂∂r(r∂uy∂r)+ g. (B.7)and at steady state,0 =µmρm1r∂∂r(r∂uy∂r)+ g. (B.8)141B.1. Problem description and assumptionsAt this point we drop the subscript y from the velocity. Inside the saliva,Ra +Rb ≤ r ≤ Ra +Rb +Rs,∂u∂r= −ρsgµsr2+1rC0, (B.9)u(r) = −ρsg2µsr24+ C0ln|r|+ C1. (B.10)In the bolus Ra ≤ r ≤ Ra +Rb,∂u∂r= −ρbgµbr2+1rC2, (B.11)u(r) = −ρbg2µbr24+ C2ln|r|+ C3. (B.12)The integration constants Ci can be found by applying the boundary con-ditions, leading to the following equations for velocity at steady state,u(r) = − gρ4µb(r2 − 2R2aln|r|)+ C1 Ra ≤ r ≤ Ra +Rb, (B.13)u(r) = − gρ4µsr2C2ln|r|+ C3 Ra +Rb ≤ r ≤ R, (B.14)C1 = −gρ4(1µb− 1µs)(Ra+Rb)2 +(C2 − gρ2µbR2a)ln|Ra+Rb|+C3 (B.15)142B.1. Problem description and assumptionsC2 =gρ2µsR2a, (B.16)C3 =gρ4µs(R2 − 2R2aln|R|). (B.17)B.1.2 Finite-Volume time-dependent solutionA finite-volume (FV) scheme can give us the time-dependent solution ifsolved numerically. We use a second-order space discretization, and a second-order implicit Crank-Nicholson time-stepping method. Special care is re-quired for correct treatment at the saliva-bolus interface.We discretize the one-dimensional problem into two sections, saliva andbolus. There are Ns and Nb saliva and bolus elements, respectively, and Nelements total (N = Ns + Nb). Each cell, i, has a width ∆ri and a height∆z. Velocities ui are computed at cell centres.For an element i, the FV formulation is derived by integrating equa-tion B.1 over each control volume,ˆdVρi∂ui∂tdV =ˆdVµi1r∂∂r(r∂ui∂r)dV + g. (B.18)Integrating and applying the divergence theorem,ρi∂ui∂t∆V =˛∂Vµi∂ui∂r∆S + g. (B.19)∂ui∂t=µiρiri∆ri(∂ui+ 12∂rri+ 12−∂ui− 12∂rri− 12)+ g, (B.20)143B.1. Problem description and assumptionswhere ri+ 12= ri +∆ri2 and ri− 12 = ri −∆ri2 .Within either fluid, the gradients can be discretized using a second-ordercentred scheme,∂ui+ 12∂r=ui+1 − ui∆ri+ 12, (B.21)where∆ri+ 12=∆ri + ∆ri+12, (B.22)so that variable spacing can be handled correctly.Discrete boundary conditionsFor the wall boundary condition at r = R, i = (Nb + Ns) − 12 = N − 12 , weuse u(R) = 0,∂uN− 12∂r=2uN−1∆rN−1. (B.23)For the boundary condition at r = Ra, i = −12 ,∂u− 12∂r= 0. (B.24)At the saliva-bolus interface, x = Ra+Rb, i = Nb− 12 , we can not assumethe left and right gradients are equal. The velocity gradient can be definedon either side to be a function of uNb− 12 .µb∂u−Nb− 12∂r= µs∂u+Nb− 12∂r, (B.25)144B.1. Problem description and assumptionsµbuNb− 12 − uNb−112∆rNb−1= µsuNb − uNb− 1212∆rNb, (B.26)leading to the following equation for uNb− 12 ,uNb− 12 =11 + auNb +a1 + auNb−1, (B.27)wherea =µbµs∆rNb∆rNb−1. (B.28)This gives us the discrete forms of the velocity gradient adjacent to thesaliva-bolus boundary for both sides,∂u−Nb− 12∂r=2(1 + a)∆rNb−1(uNb − uNb−1), (B.29)∂u+Nb− 12∂r=2a(1 + a)∆rNb−1(uNb − uNb−1). (B.30)Time integrationThe time integral can be discretized using the second-order implicit Crank-Nicholson method,∂ui∂t=un+1i − uni∆t=12(Fn+1i− 12+ Fni− 12)+12(Fn+1i+ 12+ Fni+ 12)+ g, (B.31)145B.1. Problem description and assumptionswhere the left and right fluxes for a cell i at time n are given byFni− 12=µiri− 12ρiri∆ri(∂ui− 12∂r)n, (B.32)Fni+ 12=µiri+ 12ρiri∆ri(∂ui+ 12∂r)n. (B.33)Performing some basic algebra results in a tri-diagonal system of equationsthat is closed by the boundary conditions given previously.B.1.3 Code verificationThe numerical solver is implemented in Java. We assume the density in boththe saliva and bolus are constant ρ = 1000. An acceleration of g = 9.81 isused to drive the flow. The velocity is initialized to zero everywhere and thesystem is evolved in time until the following convergence criteria is met,|un+10 − un0 | < 10−13. (B.34)At steady state, for a variety of choices of µb, µs, Rs, Rb, Ra, we areable to demonstrate second order convergence in space by comparing to thetheoretical steady-state solution of Section B.1.1.For a single material, µb = µs, and Ra = 0, the flow is the NewtonianHagen-Poiseuille flow. A theoretical time-dependent solution can be cal-culated from an infinite series of Bessel functions [70]. In these cases, thesolver shows good agreement with the theoretical results.146B.1. Problem description and assumptions0 0.5 1 1.5 2·10− = 0.05t = 0.10t = 1.75r radial positionv zvelocityexactFVFigure B.2: Theoretical startup flow solution vs. FV solver for Newtonianflow in a pipe. Ra = 0, Rb = Rs = 0.02, Ns = Nb = 5, µb = µs = 1,ρ = 1000, g = 9.81.147B.2. Time dependent flowsMesh size(Nb +Ns) L2 Error Order10 0.00290976 -20 7.27920675E-4 1.9990492440 1.82010841E-4 1.9997568680 4.55046340E-5 1.99993900160 1.13762741E-5 1.99998534Table B.1: L2 error convergence w.r.t. mesh size. Ra = Rs = Rb = 0.01,Ns = Nb, µb = mus = 1.0.Mesh size(Nb +Ns) L2 Error Order10 0.00485712 -20 0.00121428 1.9999998140 3.03569988E-4 2.0000001280 7.58924447E-5 2.00000099160 1.89730527E-5 2.00000445Table B.2: L2 error convergence w.r.t. mesh size. Ra = 1., Rs = 0.001,Rb = 0.01, Ns = Nb, µb = 1, µs = 0.01.B.2 Time dependent flowsThe time-dependent startup velocity profile for a full-slip condition is:u(t) = gt, (B.35)where g is acceleration due to gravity and t is the elapsed time. Consider ashort-duration two-material startup flow where the bolus viscosity is muchgreater than the saliva viscosity. This section shows that the full-slip solutiongiven above is a reasonable approximation for the bolus velocity profile forcertain choices of saliva viscosity, µs, and the thickness, Rs. These choices ofµs and Rs fall within the range of values that have been published in the lit-148B.2. Time dependent flows0 0.5 1·10−2024xuyvelocityfull-slip0.5s0.4s0.3s0.2s0.1s0.5s unlubricatedFigure B.3: Time dependent solutions of a two-material start-up pipe flowat various times. Low viscosity saliva lubricates a high-viscosity bolus inthe centre of the pipe. For this plot, the parameters were N = 20, Ra = 0,Rb = 10−2, Rs = 10−3, µb = 1, µs = 2e− 3. For high viscosity boluses, theno-slip should give a much better approximation to the bolus flow than ano-slip condition.erature, as long as the time duration is low enough. For example, Figure B.3shows the velocity profiles at 0.1s intervals up to 0.5s. The bolus has µb = 1,Rb = 0.01, while the saliva has µs = 0.002 and thickness Rs = 0.001. Forshort time durations, the full-slip velocity gives a good approximation to thesimulated bolus velocity profile. As the simulation progresses, the full-slipapproximation begins to overshoot the simulated bolus velocity, however itstill gives a much better approximation than an unlubricated (no-slip) ve-locity profile, also shown in Figure B.3. The unlubricated velocity profile at0.5s was obtained by setting µs = µb = 1.The quality of the full-slip approximation is most sensitive to three fac-tors, µs, Rs, and the time duration. For a specific time, for example 0.25safter the start of the flow, we can plot the bolus flow rate, Q as a function149B.2. Time dependent flowsSaliva viscosity [cP]Saliva Width [microns]Flow rate at 0.25s as percentage of full−slip flow rate0. 1 1.5 2 2.5 3 3.5 4 4.5 550100150200250300350400450500Figure B.4: Contour plot of flow rates, as percentage of the full-slip flowrate, at 0.25s, as a function of both µs and Rs.of µs and Rs. Figure B.4 shows a contour plot where the lines correspondto the flow rate Q as a percentage of the theoretical full-slip flow-rate. Onthis plot, a no-slip (non-lubricated) flow can be found at the bottom edgeof the plot, where Rs approaches zero. At 0.25s, the full-slip approximationmay be reasonable for many choices of Rs and µs.We can plot the constant contour line representing Q = 85% of the full-slip flow rate as a function of Rs and µs. Four such contour lines are shownin Figure B.5, representing the 85% cutoff line for times of 0.1s, 0.25s, 0.5s,and 1.0s. Areas to the left of the contour lines show the possible choices ofRs and µs that would be well approximated by a full-slip condition. These150B.2. Time dependent flows0.5 1 1.5 2 2.5 3 3.5 4 4.5 550100150200250300350400450500Flow rate cuts within 85 % of full−slipSaliva viscosity [cP]Saliva Width [microns]t = 0.10st = 0.25st = 0.50st = 1.00sFigure B.5: Regions to the left of each contour represent choices of width,Rs, and viscosity, µS , for which a full-slip condition gives an approximationof flow rate within 85% of the actual flow rate.experiments used a fixed time step ∆t = 10−3s, and bolus viscosity of µb = 1.Decreasing the time step to 10−4 did not change the results significantly. Theresults were not sensitive to changes in bolus viscosity.B.2.1 Measured values of saliva thickness and viscosityRantonen and Meurman [104] measured the mean viscosity of saliva at 90s−1and 37◦C of stimulated saliva to have a viscosity around 2–3 cP. The meanviscosity of unstimulated saliva was measured at around 4–8 cP. The thick-ness (Rs) of salivary film in the oral cavity has been estimated [20] to have151B.2. Time dependent flowsan average thickness of 50–100 microns, however this was for unstimulatedsaliva. This might underestimate the thickness of saliva during swallowing,since for example, the flow rate of stimulated saliva was measured to beabout 5 times greater than for unstimulated saliva [104]. However it is notclear how this could be extrapolated to thickness of saliva during swallowing.152


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