Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

An assessment of lattice Boltzmann method for swallowing simulations Mazhari, Seyed Babak 2016

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

Item Metadata

Download

Media
24-ubc_2016_september_mazhari_babak.pdf [ 2.7MB ]
Metadata
JSON: 24-1.0305020.json
JSON-LD: 24-1.0305020-ld.json
RDF/XML (Pretty): 24-1.0305020-rdf.xml
RDF/JSON: 24-1.0305020-rdf.json
Turtle: 24-1.0305020-turtle.txt
N-Triples: 24-1.0305020-rdf-ntriples.txt
Original Record: 24-1.0305020-source.json
Full Text
24-1.0305020-fulltext.txt
Citation
24-1.0305020.ris

Full Text

An Assessment of Lattice-Boltzmann Method for Swallowing Simulations   by  Seyed Babak Mazhari  B.Sc. Mechanical Engineering, University of Tehran, 2012   A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  MASTER OF APPLIED SCIENCE  in  THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Mechanical Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)   June 2016   © Seyed Babak Mazhari, 2016  ii  Abstract Lattice Boltzmann is a fixed grid particle based method originated from molecular dynamics which uses a kinetic-based approach to simulate fluid flows. The fixed grid nature and simplicity of lattice Boltzmann algorithm makes it an appealing approach for preliminary swallowing simulations. However, the issues of compressibility effect and boundary/initial condition implementation can be the source of instability and inaccuracy especially at high Reynolds simulations. The current work is an assessment of the lattice Boltzmann method with respect to high Reynolds number flow simulations, compressibility effect of the method, and the issue of boundary and initial condition implementation. Here we investigate the stability range of the lattice Boltzmann single relaxation and multi relaxation time models as well as the issue of consistent boundary/initial condition implementation. The superior stability of multi relaxation time (MRT) model is shown on the lid-driven cavity flow benchmark as a function of Reynolds number. The computational time required for the SRT model to simulate the li-driven cavity flow at Re=3200 is about 14 times higher than the MRT model and it’s shown that computational time is related to the third power of lattice resolution. It is suggested that single relaxation time model is inefficient for simulations with moderately high Reynolds number Re>1000 and the use of multi relaxation time model becomes necessary. Compressibility effect is the next topic of study where the incompressible lattice Boltzmann method is introduced. The compressibility error of the method surpasses the spatial discretization error and becomes the dominant source of error as the flow Reynolds number increases. It is shown on a 2D Womersley flow benchmark that the physical time step required for LBM is about 300 times larger than the physical time step of the finite volume implicit solver while generating results with the same order of accuracy at Re=2000. Due to the compressibility error inherent to the method, lattice Boltzmann is not recommended for preliminary swallowing simulations with high Reynolds number, since implicit time advancement methods can generate results with the same order of accuracy in noticeably less computational time.             iii  Preface  This dissertation is original, unpublished, independent work by the author, S. B. Mazhari.                            iv  Table of Contents Abstract ................................................................................................................................................................. ii Preface .................................................................................................................................................................. iii Table of Contents ..................................................................................................................................................iv List of Tables .......................................................................................................................................................... v List of Figures .......................................................................................................................................................vi Nomenclature ..................................................................................................................................................... viii Acknowledgments .................................................................................................................................................ix Chapter 1. Introduction ........................................................................................................................................... 1 1.1 Swallowing simulation characteristics ......................................................................................................... 1 1.2 General review of the lattice Boltzmann method ......................................................................................... 2 1.3 Summary of the thesis .................................................................................................................................. 4 Chapter 2. Theory of the lattice Boltzmann method ............................................................................................... 6 2.1 Statistical mechanics and Boltzmann equation ............................................................................................ 7 2.1.1 Lattice-Boltzmann BGK Formulation .................................................................................................. 8 2.1.2 Multi relaxation time model................................................................................................................ 11 2.1.3 Chapman-Enskog expansion ............................................................................................................... 13 2.2 Physical and lattice systems ....................................................................................................................... 14 2.2.1 Sources of error in LBM ..................................................................................................................... 16 2.3 Swallowing simulation characteristics and LBM accuracy ........................................................................ 16 Chapter 3. Boundary conditions ........................................................................................................................... 18 3.1 Boundary conditions in a swallowing simulation ...................................................................................... 19 3.2 Solid boundary treatments .......................................................................................................................... 20 3.2.1 Manipulation of corner nodes ............................................................................................................. 21 3.2.2 Bounce-back scheme .......................................................................................................................... 21 3.3 Open boundary conditions ......................................................................................................................... 32 3.3.1 Zou/He boundary condition ................................................................................................................ 32 3.3.2 Extrapolation scheme .......................................................................................................................... 40 3.4 Curved boundary conditions ...................................................................................................................... 43 3.4.1 Filippova-Hanel boundary dynamics .................................................................................................. 44 3.5 Transient simulations ................................................................................................................................. 51 Chapter 4. Conclusion .......................................................................................................................................... 57 References ............................................................................................................................................................ 60 v  List of Tables Table 3.1- Comparison of the center of vortex ‘y’ coordinates for cavity flow between Ghia and LBM at different Reynolds and Mach numbers ............................................................................................. 23 Table 3.2- Computational time comparisons for cavity flow simulations at Re=100 and Re=400 ...... 25 Table 3.3- Center of vortex coordinate comparison between Ghia and MRT lattice Boltzmann for cavity flow at Re=10000 ................................................................................................................................. 31 Table 3.4- Minimum stable collision frequency for 2D Poiseuille flow with respect to Reynolds number and lattice resolution ............................................................................................................................ 35 Table 3.5- Center of vortex coordinate comparison for 2:1 contraction flow between lattice Boltzmann and Ansys Fluent finite volume solver ................................................................................................. 42 Table 3.6- Choice of relaxation rate 𝑠𝜗 as a function of Reynolds number ......................................... 56                  vi  List of Figures Figure 2.1- D2Q9 lattice structure .......................................................................................................... 8 Figure 3.1- Configuration of D2Q9 lattice structure with solid and fluid nodes sketched as solid and hollow points ........................................................................................................................................ 20 Figure 3.2- Schematic view of the halfway bounce back boundary treatment during time step “t” .... 21 Figure 3.3- Cavity streamline configuration (a) Re=100, collision frequency=1.66. (b) Re=3200, collision frequency=1.66 ...................................................................................................................... 22 Figure 3.4- Non-dimensional velocity along the vertical and horizontal centerlines (a) horizontal velocity at Re=100 and Re=1000 (b) horizontal velocity at Re=400 and Re=3200 (c) vertical velocity at Re=100 and Re=1000 (d) vertical velocity at Re=400 and Re=3200. .............................................. 24 Figure 3.5- CPU time as a function of lattice resolution (a) Re=100 (b) Re=400 ................................ 25 Figure 3.6- Convergence behavior of the average lattice Mach number for cavity flow (a) Re=100 (b) Re=400 ................................................................................................................................................. 26 Figure 3.7- Stability comparison between SRT and MRT models on the cavity flow ......................... 27 Figure 3.8- Centerline velocity comparison between MRT lattice Boltzmann and Ghia at Re=100 (a) vertical centerline (b) horizontal centerline .......................................................................................... 28 Figure 3.9- Unknown distribution functions configuration for boundary node ‘b’ .............................. 29 Figure 3.10- Streamlines for cavity flow at Re=10000 using MRT lattice Boltzmann model (𝜔 =1.994) ................................................................................................................................................... 30 Figure 3.11- Non-dimensional centerline velocity profile comparison between MRT lattice Boltzmann and Ghia at Re=10000 (a) horizontal centerline (b) vertical centerline ............................................... 30 Figure 3.12- Steady state velocity and pressure contours for 2D Poiseuille flow at Re=10 (a) non-dimensional velocity contour (b) lattice density contour ..................................................................... 34 Figure 3.13- Fully developed 2D Poiseuille flow velocity profile for lattice resolutions Ny=50, 100, and 200 compared to analytical solution ..................................................................................................... 34 Figure 3.14- L2 norm error analysis for Poiseuille flow ...................................................................... 35 Figure 3.15- Convergence history for Poiseuille flow with collision frequency set to 1.96 (𝜔 = 1.96) .............................................................................................................................................................. 36 Figure 3.16- Effect of relaxation parameter on the convergence rate and accuracy (a) history of convergence for average Mach number for 𝜏=0.6, 1, 2 (b) fully developed velocity profile for 𝜏=0.6, 1, 2 ............................................................................................................................................................ 37 Figure 3.17- Lattice density contour for the Poiseuille flow (a) iteration: 100 (b) iteration: 400 (c) iteration: 1000 (d) iteration: 15000 ....................................................................................................... 38 vii  Figure 3.18- Centerline lattice density profile for Poiseuille flow at different iterations ..................... 39 Figure 3.19- Velocity contour and streamline for 2:1 contracting channel flow using extrapolation scheme at Re=10 .................................................................................................................................. 41 Figure 3.20- Centerline non-dimensional velocity profile for 2:1 contracting channel flow at Re=10 42 Figure 3.21- Concave curved wall condition with 5 unknown distribution functions streaming from solid into fluid ...................................................................................................................................... 43 Figure 3.22- The illustrative geometry for Filippova-Hanel scheme ................................................... 44 Figure 3.23- Geometry of the flow over a bounded cylinder ............................................................... 45 Figure 3.24- Velocity contours for flow over a bounded cylinder at Re=10 (a) x-velocity (b) y-velocity .............................................................................................................................................................. 46 Figure 3.25- Velocity profile comparison between lattice Boltzmann and Ansys Fluent finite volume solver .................................................................................................................................................... 46 Figure 3.26- Pressure variations at point ‘M’ due to pressure wave interactions at the inlet boundary .............................................................................................................................................................. 47 Figure 3.27- Drag coefficient as a function of lattice nodes/cylinder diameter (a) Re=1 (b) Re=10 (c) Re=40 ................................................................................................................................................... 50 Figure 3.28- Non-dimensional velocity profile comparison between incompressible lattice Boltzmann and Ansys Fluent finite volume solver at Re=10 and 𝜔 = 2𝜋 ............................................................. 54 Figure 3.29- Non-dimensional velocity profile comparison between incompressible lattice Boltzmann and Ansys Fluent finite volume solver at Re=10 and 𝜔 = 20𝜋 ........................................................... 55          viii  Nomenclature Re  Reynolds number      CFD Computational fluid dynamics f   Particle distribution function (PDF)  𝜏   Single relaxation time rate 𝑓𝑒𝑞   Equilibrium PDF              𝑒𝑖      Lattice velocity in ith direction 𝑓𝑛𝑒𝑞 Non-equilibrium PDF     Ma   Mach number 𝑐𝑠   Lattice speed of sound       𝜗   Kinematic viscosity 𝜇   Dynamic viscosity        𝜌  Density SRT  Single relaxation time      MRT  Multi relaxation time 𝑥𝑖  Position of the ith particle    𝑃𝑖   Linear momentum of the ith particle   H   Hamiltonian operator      𝜀   Particle velocity vector c  Lattice unit velocity      PDE  Partial differential equation p   Pressure         𝑤𝑖   Equilibrium PDF weighting factor e   Energy        j  Momentum density q  Energy flux        T   Temperature Ω𝑖  Collision operator      L   Length 𝜔  Collision frequency      FHM  Filippova-Hanel boundary scheme 𝑓(∗)   Superficial PDF in FHM method    ∆  Length ratio in FHM method X  Fitting weight in FHM method   A  Area 𝜏𝑖𝑗  Shear stress         𝐶𝐷   Drag coefficient F   Force        D  Diameter 𝜉   Bulk viscosity       𝛼  Womersley number               ix  Acknowledgments  I would like to express my great appreciation to my supervisors Dr. Sidney Fels and Dr. Sheldon Green who guided me through this program with their insightful instructions. This work couldn’t have been done without their continuous support and patience.   I would like to thank Andrew Ho and Peter Anderson for their suggestions during this research. I’m grateful for their kindness and assistance. I would also like to thank Pouyan Jahangiri for his friendly suggestions and scientific discussions during this work.  I would like to thank the Natural Science and Engineering Research Council of Canada (NSERC) and members of the Parametric Modeling of Human Anatomy (PMHA) group for their help and support.  Lastly, I would like to thank my family and friends for their unconditional help and support during this work.    Seyed Babak Mazhari                 1  Chapter 1. Introduction The process of swallowing starts with the voluntarily squeezing action of the tongue against hard palate in order to drive the food bolus into the upper pharynx region where the involuntary phase of swallowing begins. The involuntary phase, by which an ingested food bolus is transported through the pharynx and the upper esophageal sphincter, involves very rapid motion of the anatomical structures. The quality and duration of bolus transport is dominantly controlled by these changes in geometry and the rheology of the bolus. The main characteristics of bolus transport during oral and pharyngeal phases of swallowing are (1) the complex boundary movements of the flow domain which are controlled by neural stimulations (2) the rheology of the bolus and (3) the free surface problem when dealing with multiphase simulations. Modeling the flow with such characteristics is computationally a very challenging problem and requires the use of numerical schemes capable of simulating complex moving boundary conditions and multi-phase flow. Although conventional CFD methods are capable of solving a vast range of fluid dynamics problems, they have many limitations when dealing with complex boundaries. Practically, one of the most challenging tasks of the traditional CFD approach is high quality mesh generation based on the geometry and the flow characteristics of the problem. It is very time consuming to generate a high quality 3D mesh in the complex geometry of the pharynx during swallowing. Furthermore, Body fitted grid methods require an automated remeshing algorithm at each time step to adopt the mesh elements based on the deformation of the pharynx during bolus transport. The high sensitivity of these techniques to the quality of mesh elements makes one not to consider them for modeling the bolus transport problem in the upper airway geometry. The central theme of this thesis is to investigate the capabilities and shortcomings of the lattice Boltzmann method for simulating the fluid flows showing similar characteristics. Being a fixed grid method, lattice Boltzmann can be considered as a capable numerical technique for flow simulations with complex curved and moving boundary conditions. On the other hand, there are many short comings to the method such as instability at high Reynolds number and compressibility effect that have to be investigated before one can apply the method to swallowing simulations. This work begins with a basic introduction to the lattice Boltzmann method and later focuses in detail on the characteristics of the method and its limitations when dealing with high Reynolds flows and boundary condition implementation in preliminary swallowing simulations. 1.1 Swallowing simulation characteristics As mentioned before, swallowing is the process of several interrelated events involving muscle contraction/relaxation by which the geometry of the pharynx changes and drives the bolus pass the upper esophageal sphincter and into the esophagus. The time needed for a bolus to pass through the pharynx and into the esophagus is usually within a second [1] once a swallow has been initiated. The velocity of the bolus is dependent on the rheology and the initial volume of the bolus. Range of velocity could be around 10-70 cm/s at the head of the bolus and 10-15 cm/s at the tail [2] depending on the rheology and the initial volume of the bolus. The Reynolds number of the bolus transport depends mainly on the bolus viscosity and can be in the range of laminar flows for highly viscous bolus to completely turbulent flows for very thin boluses.  2  As an example take the effective diameter of the pharynx to be D=2 cm and the characteristic velocity to be U0=10 cm/s for the flow of water inside the pharynx. This will result in a Reynolds number of Re≈2000 for the flow of water with kinematic viscosity 𝜗 ≈ 10−6𝑚2. 𝑠−1. It is obvious that Reynolds number increases for fluids with less viscosity or higher velocity compared to water. On the other hand, the Reynolds number of the flow can decrease significantly for thick fluids such as honey (𝑅𝑒 ≈ 1) with kinematic viscosity 𝜗 ≈ 0.0022 𝑚2. 𝑠−1. In the range of low Reynolds flows, the lubricating effect of saliva becomes significant and it is necessary to consider this effect in the simulations using multiphase models or using slip boundary conditions. Various classes of multiphase models are used to capture the gross hydrodynamic features of the multiphase flows. The most popular class is to capture the interface directly on a regular stationary grid. The best examples are the volume of fluid (VOF) and the marker and cell methods. In the marker and cell method, the marker particles are advected for each fluid, whereas in volume of fluid method, marker functions are advected. Recent applications of this class include the method of Yabe [3] and the phase field method of Jacqmin [4]. The main difficulty with this class is the maintenance of sharp boundaries between different phases and the calculation of surface tension. The second class which is potentially the most accurate class between the continuum multiphase models uses separate boundary fitted grids for each phase. This class of methods is mostly used for the relatively simple geometries and applications to complex 3D geometries with unsteady deforming boundaries, such as pharyngeal swallowing, are very rare. The unsteady simulation of a 3D bubble is the most impressive example [5]. The third class is specified to Lagrangian methods in which the grid follows the fluid. While this class is fairly hard to implement for complicated unsteady flows, Tezdayur [6] have published accurate results for the 3D unsteady motion of many spherical particles. Finally the fourth class is specified to front tracking where a separate front marks the interface on a fixed grid. The grid around the interface is then modified to make the grid line in alignment with the interface. The method of lattice Boltzmann uses an intermolecular interaction model proposed by Shan and Chen [7] to simulate two phase flows. This model does not require the explicit tracking of the interface which makes it straight forward to implement, but it is hard to simulate flows with density ratio higher than 50. Recent developments in lattice Boltzmann multiphase simulations include the work by Li [8] that allows density ratios around 500 at moderately high Reynolds (40<Re<1000) droplet splashing flows.    The issue of boundary condition implementation has significant importance in a swallowing simulation since the bolus is mainly driven by the tongue movement and pharyngeal wall contractions. In addition, the fluid structure interactions between the bolus and the anatomical structure of the pharynx is a problem of interest during swallowing simulations. Different techniques are used for flow simulations with arbitrarily shaped moving boundaries which can be categorized into two main classes of Lagrangian and Eulerian methods. In a Lagrangian method, the grid is reconfigured at each time step to conform to the shape of the boundary, therefore, boundary conditions can be applied precisely on the boundary location. On the other hand, Eulerian methods use a fixed grid through which the boundary moves, therefore, the boundary condition is usually applied using extrapolation schemes near the boundary. Although the Lagrangian methods yield more precise solutions, the extra computational effort of remeshing at each time step can offset the higher accuracy.  1.2 General review of the lattice Boltzmann method During the past two decades, lattice Boltzmann method has been introduced from the heart of statistical mechanics as an alternative numerical scheme for fluid flow simulations. In contrast to 3  conventional numerical techniques such as finite difference, finite volume, and finite element methods that are rooted in discretization of the macroscopic governing equations for flow motion, lattice Boltzmann is based on microscopic models and mesoscopic kinetic equations [9], thus lattice Boltzmann can be considered as a microscope for fluid mechanics and a telescope for molecular dynamics. It originated from the lattice gas automata (LGA) that was introduced in 1986 as a discrete particle kinetics based on discrete space and time. Frisch et al [10] recovered the Navier-Stokes equation from the LGA kinetic model and realized the importance of symmetry of the lattice for recovery of macroscopic conservation equations. The LGA model introduced by Frisch et al is based on the idea that particles can move between the sites of regular lattice and may collide with each other only on the lattice grid points. The evolution equation of lattice gas automata is written as:    𝑛𝑖(𝑥 + 𝑒𝑖∆𝑡, 𝑡 + ∆𝑡) = 𝑛𝑖(𝑥, 𝑡) + Ω𝑖(𝑥, 𝑡)  (1.1.1)  Where 𝑛𝑖 is a Boolean variable describing the presence of a particle at a specific lattice grid point and moving at a specific lattice direction prescribed by the grid structure. The discretized lattice velocity in the ith direction is denoted by 𝑒𝑖 and ∆𝑡 is the time step. The collision operator Ω𝑖(𝑛(𝑥, 𝑡)) is responsible for the recovery of the corresponding macroscopic equation and is defined by a set of collision rules that conserve mass and momentum on a symmetric lattice grid. By averaging the Boolean variable 𝑛𝑖 on sub regions of the lattice we can construct the macroscopic properties of the flow.   Lattice Boltzmann obeys the same underlying principles as LGA, however, the Boolean particle occupation variable 𝑛𝑖 is replaced with a continuous distribution function 𝑓𝑖(𝑥, 𝑡) in order to get rid of the statistical noise. The evolution equation for lattice Boltzmann method is written as:  𝑓𝑖(𝑥 + 𝑒𝑖∆𝑡, 𝑡 + ∆𝑡) = 𝑓𝑖(𝑥, 𝑡) + Ω𝑖(𝑥, 𝑡)  (1.1.2)  Using the Bhatnagar, Gross and Krook (BGK) model [11] for the collision operator Ω𝑖, we can recover the incompressible Navier-Stokes equations up to second order of Mach number using the Chapman-Enskog analysis. The lattice Boltzmann BGK (LBGK) equation is given by:  𝑓𝑖(𝑥 + 𝑒𝑖∆𝑡, 𝑡 + ∆𝑡) = 𝑓𝑖(𝑥, 𝑡) −1𝜏(𝑓𝑖(𝑥, 𝑡) − 𝑓𝑖𝑒𝑞(𝑥, 𝑡)) (1.1.3)  Where 𝜏 is the single relaxation time parameter and 𝑓𝑖𝑒𝑞 is the distribution function when the flow reaches the equilibrium state defined by the Maxwellian equilibrium function in order to recover Navier-Stokes equations for mass and momentum conservation [12]. Based on the evolution equation (1.1.3), the distribution function 𝑓𝑖 relaxes towards the equilibrium distribution function 𝑓𝑖𝑒𝑞 with a relaxation rate of (1/𝜏). The relaxation parameter 𝜏 in the model is related to the kinematic viscosity of the fluid by the equation 𝜗 = 𝑐𝑠2𝛿𝑡(2𝜏 − 1)/𝑑 where 𝑐𝑠 is the speed of sound, 𝛿𝑡 is the time step, and 𝑑 is the dimension of the space with 𝑑=2 for two dimensional space and d=3 for three dimensional space discretization. The macroscopic characteristics such as density and momentum are then calculated as the hydrodynamic moments of the distribution functions at each lattice node. In order to obtain a positive kinematic viscosity, the relaxation parameter has to be chosen in the range of 𝜏 > 0.5 with 𝜏 → 0.5 approaching the inviscid flow limit and 𝜏 → ∞ approaching the creeping flow limit. It’s worth mentioning that approaching the inviscid flow limit 𝜏 → 0.5 can be the source of instability in lattice Boltzmann simulations due to the high energy dissipation rate for very low viscous fluids. Generally, the lattice time step is chosen as 𝛿𝑡 = 1 during lattice Boltzmann simulations for simplicity 4  and the time mapping between the Lattice and physical simulations are performed using non-dimensional parameters. This issue is further explained in section (2.2). Two approaches can be used to overcome the instability issue when relaxation parameter is approaching the inviscid limit. First, we can use finer lattice resolutions in order to capture the high gradient flow characteristics, however, this can lead to extremely long computational time for complex flows such as pharyngeal bolus transport. The second approach is to use different relaxation rates for different hydrodynamics moments (multi relaxation time model) of distribution functions. In the latter, different moments such as density, momentum, and energy are relaxed towards their equilibrium state using different relaxation rates. This will allow for the energy to dissipate with a different speed which results in more stable solutions when 𝜏 → 0.5.  Since the birth of lattice Boltzmann method, it has been used in a variety of flow problems including suspension flows, porous media, and magnetohydrodynamics. The flexibility, intrinsic parallelism, and simplicity of the method have made lattice Boltzmann a popular numerical technique for flow simulations. Also, the fixed grid nature of the method has made it an appealing approach for moving boundary problems. Despite the mentioned advantages inherent to the lattice Boltzmann method, this novel technique suffers from many shortcomings such as instability at high Reynolds flows (low viscosity fluids), large computational error for flows with moderately high lattice Mach number (0.1<Ma<0.5), and extensive computational time when modeling flows with high gradient regions that requires very fine lattice resolution.  1.3 Summary of the thesis  The present work aims at investigating the capabilities and shortcomings of lattice Boltzmann method for preliminary simulations of swallowing. The main focus of the work can be divided into five aspects. First, the theoretical background of the lattice Boltzmann is introduced in Chapter (2). The single relaxation time and multi relaxation time models are explained and sources of error in a lattice Boltzmann simulations are discussed. Section (2.3) concludes the chapter with possible shortcomings of the method mainly in terms of low accuracy and instability at high Reynolds number flows. This issue will become critical when the method is applied to preliminary swallowing simulations at moderately high Reynolds number (Re>1000). The sensitivity of the method to the choice of relaxation parameter 𝜏 and lattice resolution is the topic of discussion later in Chapter (3) where different flow benchmarks are investigated and boundary condition implementation is explained for open and curved boundaries in preliminary simulations of swallowing.     Secondly, the issue of straight wall boundary condition implementation is studied in section (3.2). The solid boundaries such as the hard palate can be modeled as straight walls in preliminary swallowing simulations, so it is insightful to study the straight wall boundary condition implementation for lattice Boltzmann method. We have looked into the stability and accuracy of the SRT and MRT models coupled with the halfway bounce back boundary treatment when simulating the cavity flow benchmark. The halfway bounce back boundary treatment is the most popular scheme for implementation of no-slip condition at straight solid walls and achieves second order of spatial accuracy when the wall is placed exactly halfway between the boundary nodes. The main purpose of the section is to demonstrate the shortcomings of the SRT model when dealing with high Reynolds flow simulations. It is observed later in section (3.2.2) that computational time needed for the SRT model to simulate cavity flow with desirable accuracy at Re>3000 is large (≈3 days) since very fine lattice resolution (800*800) has to be used. The use of multi relaxation time model will allow to use relaxation rates closer to the inviscid 5  limit 𝜏 → 0.5 which leads to significantly less compressibility error (𝑂(𝑀𝑎2)), as a result, lower lattice resolution can yield more stable and accurate solutions. The superior stability of the multi relaxation time model is the main topic of study in section (3.2.2) and the effect of boundary condition implementation on the solution stability is further studied. Section (3.3) focuses on the open inlet and outlet boundary conditions in lattice Boltzmann simulations. The inlet/outlet conditions in a swallowing simulation are not very well defined. The constant pressure condition (Zou/He) is investigated as the choice of inlet boundary and unidirectional flow with zero velocity gradient (Extrapolation scheme) is investigated as the choice of outlet boundary for preliminary swallowing simulations. The main focus of the section is to observe the compressibility effect inherent to the lattice Boltzmann method and to measure the stability of the LBE solver as a function of Reynolds number and the choice of relaxation parameter on a 2D Poiseuille flow. The compressibility effect of the LBE method is the source of instability and inaccuracy during time dependent flows. The initial discontinuity between the imposed boundary values for pressure and velocity and the initial condition defined in the bulk will result in variable oscillations in the flow domain with an acoustic nature, therefore, consistent initial condition implementation becomes as important as boundary implementation in time dependent problems. The choice of relaxation parameter also affects the rate by which the initial continuity damps out inside the domain and the flow reaches equilibrium state. In section (3.4), we have looked into the implementation of boundary conditions as well as force evaluation on curved surfaces during lattice Boltzmann simulations. It is possible to break the curved boundary into stair wise steps and use the bounce back boundary scheme to simulate no-slip condition, however, bounce back only achieves first order accuracy since the boundary wall cannot be placed halfway between the lattice nodes for all boundary nodes. Additionally, very fine lattice resolution, i.e. high computational cost, is required to reduce the error introduced by the geometry breakdown. The Filippova-Hanel boundary treatment is introduced in section (3.4) which uses direct fitting distribution functions on the curved boundary and achieves second order accuracy while conserving the integrity of the geometry. The flow over a bounded cylinder is studied as the benchmark flow to investigate the accuracy of the lattice Boltzmann method coupled with Filippova-Hanel boundary treatment and momentum exchange force evaluation scheme. We have looked into the acoustic nature of pressure oscillations in lattice Boltzmann simulations due to compressibility effect. These oscillations are the source of instability at high lattice Mach number flow simulations and care should be taken when implementing the initial and boundary conditions and choosing the relaxation parameter. Additionally, the accuracy is investigated using the value of drag coefficient as a function of lattice resolution (lattice nodes/diameter of the cylinder) for three different Reynolds numbers Re=1, 10, 40.   Transient flow simulation is the subject of study in section (3.5). Accurate initial and boundary condition implementation is highly important in lattice Boltzmann simulations since the initial pressure or velocity discontinuity on the boundary is the source of density oscillations in lattice Boltzmann simulations. These adverse oscillations can cause instability when dealing with complex flows such as pharyngeal bolus transport. In section (3.5), the iterative method of Mei [13] is first introduced for consistent initial condition implementation, then the compressibility effect of lattice Boltzmann method is investigated on the channel flow with periodic inlet pressure boundary condition. Finally, the thesis is concluded in Chapter (4) with an assessment of the shortcomings and capabilities of lattice Boltzmann method for preliminary simulations of swallowing.   6  Chapter 2. Theory of the lattice Boltzmann method  Historically, lattice Boltzmann method was developed from its forerunner lattice gas automata (LGA) [14]. The fundamental idea behind LGA is to model the molecular dynamics of the fluid using artificial particles placed on a structured lattice grid. The microscopic interactions between particles can lead to the macroscopic equations to describe the same flow. The corresponding macroscopic equations will result in measurements of macroscopic variables such as velocity and pressure inside the flow domain. The microscopic interactions between particles are modeled using two consecutive steps called collision and streaming which will be explained later in the chapter. During the interactions, lattice symmetry plays a key role in conserving mass, momentum, and angular momentum [15]. The first lattice structure was developed by Hardy, de Pazzis and Pomeau [16]. In their model they used a square lattice with 4-fold rotational symmetry which was only able to achieve mass and momentum conservation. Later Frisch, Hasslacher and Pomeau took a step further and used a hexagonal lattice structure which was able to conserve the angular momentum as well as mass and momentum. The hexagonal lattice model made it possible not only to retrieve diffusion equation but also Navier-Stokes equations [10].  Lattice gas automata uses an integer value to represent the presence or absence of a particle at each lattice node moving at specific directions depending on the lattice structure. This model suffered from many shortcomings such as statistical noise, instability at high Reynolds, and lack of Galilean invariance. In order to overcome the statistical noise, McNamara and Zanetti [17] used the idea of real-valued particle distribution functions instead of integer values. The particle distribution function represents the probability of a particle being present at a specific lattice node and moving along a specific direction based on the lattice structure. Later, lattice Boltzmann was established as an independent method derived from direct discretization of the kinetic Boltzmann BGK equation [18]. The Boltzmann equation is the cornerstone of the kinetic molecular dynamics and with linearization of the collision operator in Boltzmann equation [19], the lattice Boltzmann BGK model was born as a successful method to retrieve Navier-Stokes equations.  After lattice Boltzmann BGK model was introduced, the popularity of lattice Boltzmann method has been increased as an alternative approach to conventional numerical methods such as finite volume, finite difference and finite element methods. Despite the increasing popularity, lattice Boltzmann method possesses several limitations [20] in terms of accuracy and stability especially when dealing with high Reynolds and high Mach number flows. In general, it is possible but very difficult for lattice Boltzmann method to achieve higher than second order accuracy in both spatial and temporal domain [21].  In section (2.1), we will introduce the Boltzmann BGK equation in the frame of statistical mechanics and will discuss various differences between Boltzmann BGK and Navier-Stokes equations with respect to numerical implementation. Section (2.1.1) explains the implementation of lattice Boltzmann BGK algorithm with single relaxation time in two consecutive steps called collision and streaming. Due to the shortcomings of single relaxation time BGK model for high Reynolds and low viscous flow simulations, multi relaxation time (MRT) model is introduced in section (2.1.2). In section (2.1.3), the Chapman-Enskog expansion is introduced and a few mathematical considerations for derivation of Navier-Stokes equations from Boltzmann BGK equation are discussed. The dimensional 7  considerations and sources of error in a lattice Boltzmann simulation are the topic of study for section (2.2) and the chapter is concluded in section (2.3) with a general discussion about possible stability and accuracy shortcomings that can be observed in a preliminary swallowing simulation using lattice Boltzmann method.  2.1 Statistical mechanics and Boltzmann equation  Statistical mechanics provides a framework for prediction of macroscopic bulk properties of materials based on the microscopic properties and interactions of atoms and molecules in the material. The motion of molecules in a material can be described by the deterministic Newton equations or Hamiltonian system and the thermodynamics conservation equations for mass, momentum and energy [22]. Ludwig Boltzmann is a pioneer in this field that used a probabilistic approach to provide a bridge between the thermodynamic equations and the statistical mechanics. The Boltzmann equation established in 1872 was the starting point for rapid development in statistical mechanics [23].  There are about 2.7 × 109 molecules in 1 cm3 of air in atmospheric conditions. The dynamics of this system can be described by solving the Newton equations for every molecule in the system, but it is obvious that the number of equations would be too large and consequently the computation time is extremely large. This problem is taken into account by considering a system of 𝑁 particles with distribution functions 𝑓𝑁(𝑥1, 𝑃1, 𝑥2, 𝑃2, … , 𝑥𝑁 , 𝑃𝑁 , 𝑡), where N is a much smaller number compared to actual number of molecules. 𝑥𝑖 is the position vector and 𝑃𝑖 is the linear momentum vector of the 𝑖𝑡ℎ molecule. The evolution of distribution function  𝑓𝑁 is governed by the Liouville equation. In three dimensional space Liouville equation is written as follow:  𝜕𝑓𝑁𝜕𝑡−∑(𝜕𝐻𝑁𝜕𝑥𝑖𝜕𝑓𝑁𝜕𝑃𝑖3𝐷𝑖=1−𝜕𝐻𝑁𝜕𝑃𝑖𝜕𝑓𝑁𝜕𝑥𝑖) = 0  (2.1.1)  Where 𝐻𝑁 is the Hamiltonian of the system. By integration over part of the phase space the reduced density is measured as:  𝐹𝑠(𝑥1, 𝑃1, 𝑥2, 𝑃2, … , 𝑥𝑠, 𝑃𝑠, 𝑡)= 𝑉𝑠∫𝑓𝑁(𝑥1, 𝑃1, 𝑥2, 𝑃2, … , 𝑥𝑁 , 𝑃𝑁 , 𝑡)𝑑𝑥1𝑑𝑃1…𝑑𝑥𝑁𝑑𝑃𝑁  (2.1.2) Where 𝑉𝑠 is a normalization coefficient. It is theoretically shown that a coupled system of differential equations for 𝐹𝑠(1 ≤ 𝑠 ≤ 𝑁) with the name of BBGKY [24] is equivalent to the Liouville equation. The Boltzmann equation has been derived from BBGKY by the assumption of two-partial local collisions with uncorrelated velocities before collision and free of external forces [25]. The continuous Boltzmann equation in partial differential form is written as follow:  𝜕𝑓𝜕𝑡+ 𝜀. ∇𝑓 = Ω(𝑓)  (2.1.3) 8  Where 𝑓 is the particle distribution function, 𝜀 is the particle velocity vector and Ω is the collision integral. The most popular model that converges to the solution of Navier-Stokes equations is the BGK model proposed by Bhatnagar, Gross and Krook [18]. In the BGK model collision integral is defined as:  Ω(𝑓) = −𝑓 − 𝑓𝑒𝑞𝜏  (2.1.4) In equation (2.1.4), 𝜏 is the single relaxation time required for the molecules to relax towards local equilibrium and 𝑓𝑒𝑞 defines the particle distribution functions at Maxwellian equilibrium state [12]. The Boltzmann equation coupled with the BGK collision integral will result in a set of first order PDEs to solve for distribution functions at different lattice directions. This differs from the Navier-Stokes based solvers in various aspects such as: 1. Navier-Stokes equations are inherently second order PDEs due to the viscous stress terms, while the Boltzmann BGK equation includes a set of first order PDEs. 2. The Boltzmann BGK equation avoids any non-linear term, while NS solvers inevitably have to take care of the non-linear convective term 𝑢. ∇𝑢.  3. NS based solvers need to solve the Poisson equation to obtain the pressure field which involves global data communication in the problem domain. Lattice Boltzmann measures the pressure locally through an equation of state 𝑝 = 𝜌𝑐𝑠2. 4. Since Boltzmann equation has arrived from the heart of molecular dynamics, the LBE model can easily be applied to micro-scale flow problems 2.1.1 Lattice-Boltzmann BGK formulation  In this section, we have looked into the implementation of lattice Boltzmann algorithm. First, we will discretize the Boltzmann BGK equation in phase space. The D2Q9 lattice structure is used for phase space discretization throughout this work. As shown in Figure (2.1), D2Q9 lattice structure consists of 9 discrete velocities in 2 dimensional space with c0 representing a stationary particle.  Figure 2.1- D2Q9 lattice structure 9  By discretizing the Boltzmann equation in phase space and using the BGK collision operator, the Boltzmann equation for the particle distribution function is given as:  𝜕𝑓𝛼(𝑥, 𝑡)𝜕𝑡+ 𝑒𝛼∇𝑓𝛼(𝑥, 𝑡) = −1𝜏(𝑓𝛼(𝑥, 𝑡) − 𝑓𝛼𝑒𝑞(𝑥, 𝑡))  (2.1.5) Where 𝑒𝛼 represents the αth discrete velocity and 𝑓𝛼  and 𝑓𝛼𝑒𝑞 are distribution function and equilibrium distribution function of the corresponding particle in αth direction. The single relaxation time 𝜏 is related to the lattice kinematic viscosity as:  𝜗 = (𝜏 −12)𝑐𝑠2𝛿𝑡   (2.1.6)  Where the speed of sound 𝑐𝑠 is a lattice dependent velocity equal to 1√3 for the D2Q9 lattice structure and 𝛿𝑡 is the lattice time step. Equation (2.1.6) provides a straightforward method for changing the fluid viscosity in the model. It is obvious that in order to have a positive viscosity, the relaxation time should satisfy the condition 𝜏 > 0.5. The limit 𝜏 → 0.5 corresponds to the flow with no viscosity effects, while 𝜏 → ∞ represents the Stokes (creeping) flow. The first limit for inviscid flows can be problematic since stability issues will appear if the lattice resolution is not high enough, especially when dealing with complex geometries and high velocity gradients. This is due to the fact that model cannot dissipate energy because of the very low viscosity. The equilibrium distribution function 𝑓𝛼𝑒𝑞 in equation (2.1.5) is calculated as:  𝑓𝛼𝑒𝑞 = 𝜌𝑤𝛼[1 +3𝑐2𝑒𝛼 . 𝑢 +92𝑐4(𝑒𝛼 . 𝑢)2 −32𝑐2𝑢. 𝑢]  (2.1.7)  Where 𝑐 = 𝛿𝑥/𝛿𝑡 represents the ratio of lattice space step over lattice time step during simulations. For simplicity, usually both space and time step sizes are chosen as 𝛿𝑥 = 𝛿𝑡 = 1 during LB simulations and the physical time step of the simulation is adjusted using the single relaxation time parameter 𝜏. This will be further explained later in section (2.2).      𝑤𝛼 is the weighting factor given as:   𝑤𝛼 ={    49                  𝛼 = 019        𝛼 = 1,2,3,4136        𝛼 = 5,6,7,8 (2.1.8) A set of nine velocity vectors for the D2Q9 lattice structure are described as follow:  𝑒𝛼 ={    0                                                                                       𝛼 = 0𝑐(cos ((𝛼 − 1)𝜋4) , sin ((𝛼 − 1)𝜋4) )            𝛼 = 1,2,3,4√2𝑐(cos ((𝛼 − 1)𝜋4) , sin ((𝛼 − 1)𝜋4))        𝛼 = 5,6,7,8 (2.1.9) 10   Discretization of equation (2.1.5) in time and space will result in a linear system of equations that can be solved for the distribution functions 𝑓𝛼 using a two-step process called collision and streaming.  𝐶𝑜𝑙𝑙𝑖𝑠𝑖𝑜𝑛 ∶  𝑓𝛼(𝑥, 𝑡 + 𝛿𝑡) = 𝑓𝛼(𝑥, 𝑡) −1𝜏[𝑓𝛼(𝑥, 𝑡) − 𝑓𝛼𝑒𝑞(𝑥, 𝑡)] (2.1.10)  𝑆𝑡𝑟𝑒𝑎𝑚𝑖𝑛𝑔 ∶  𝑓𝛼(𝑥 + 𝑒𝛼𝛿𝑡, 𝑡 + 𝛿𝑡) = 𝑓𝛼(𝑥, 𝑡 + 𝛿𝑡) (2.1.11) After the distribution functions are evaluated at each time step, density and momentum are calculated as hydrodynamic moments of distribution functions:  𝜌 = ∑ 𝑓𝛼8𝛼=0= ∑ 𝑓𝛼𝑒𝑞8𝛼=0 (2.1.12)  𝜌𝑢 = ∑ 𝑒𝛼𝑓𝛼8𝛼=0= ∑ 𝑒𝛼𝑓𝛼𝑒𝑞8𝛼=0 (2.1.13) Pressure is then obtained locally at each lattice node using the equation of state:  𝑝 = 𝜌𝑐𝑠2 = 𝜌/3  (2.1.14) The general algorithm for the lattice Boltzmann single relaxation time BGK model can be implemented in 6 steps: 1. First step in any LB simulation is to rescale the problem from physical system to the dimensionless system using characteristic variables for length and velocity. The mapping between the physical and the lattice system of units is done using the non-dimensional parameters which should be identical in both systems. (This will be further discussed in section (2.2)). 2. After rescaling, the problem is discretized on a uniform grid in Cartesian coordinates and the lattice spatial step 𝛿𝑥 is determined based on the geometry and flow properties. Time step 𝛿𝑡 is determined according to the choice of lattice resolution in order to have limited Mach number in the flow. 3. The next step is to initialize the distribution functions in the geometry of the problem. The difficulty is to construct the distribution functions based on the macroscopic initial velocity and pressure. However, initial pressure is not always given in the problem and we have to solve Poisson problem to find the initial pressure field.   4. After the initial condition is specified, the iterative loop for evolution of distribution functions starts. Similar to initial condition implementation, we need to construct the unknown distribution functions streaming from the solid into the bulk at the boundaries based on the imposed macroscopic boundary conditions. Different boundary dynamics are discussed in sections (3.2), (3.3), and (3.4) to construct the unknown distribution functions when dealing with Dirichlet and Neumann boundary conditions for velocity and pressure. 5. After boundary implementation, collision is carried out to find the post collision distribution functions at all lattice nodes. Collision is applied to all the lattice nodes in the bulk or sitting on the boundaries. 11  6. The final step inside the evolution loop is to apply the streaming. Streaming is performed in all directions for the nodes inside the fluid region and directions pointing into the solid for boundary nodes. If the convergence criteria inside the loop is met, the solution is converged and data can be extracted for post processing.  2.1.2 Multi relaxation time model  The fundamental shortcoming of the single relaxation time model is that it’s difficult to damp the acoustic modes in the transient pressure field when dealing with high Reynolds number flows. This is due to the fact that bulk and shear viscosities are considered identical in the SRT BGK model [26]. This problem can be solved by the use of multi relaxation parameters. The superior stability of the MRT compared to SRT method is investigated on the cavity flow benchmark in section (3.2). The lattice Boltzmann MRT collision model is written as:  𝐹(?⃗? + 𝑒𝛿𝑡, 𝑡 + 𝛿𝑡) − 𝐹(?⃗?, 𝑡) = −𝑀−1. ?̂?. [𝑅(?⃗?, 𝑡) − 𝑅𝑒𝑞(?⃗?, 𝑡)] (2.1.15)  In MRT model, a set of new variables 𝑅 = (𝜌, 𝑒, 𝜀, 𝑗𝑥, 𝑞𝑥 , 𝑗𝑦, 𝑞𝑦, 𝑝𝑥𝑥, 𝑝𝑥𝑦)𝑇 is introduced which is related to the set of distribution functions 𝐹 = (𝑓0, 𝑓1, 𝑓2, 𝑓3, 𝑓4, 𝑓5, 𝑓6, 𝑓7, 𝑓8)𝑇 using the linear transformation matrix M (equation (2.1.16)). ?̂? = 𝑑𝑖𝑎𝑔(𝑠0, 𝑠1, 𝑠2, 𝑠3, 𝑠4, 𝑠5, 𝑠6, 𝑠7, 𝑠8) is a non-negative 9×9 diagonal relaxation matrix with relaxation rates as diagonal components.   𝑅 =(      𝜌𝑒𝜀𝑗𝑥𝑞𝑥𝑗𝑦𝑞𝑦𝑝𝑥𝑥𝑝𝑥𝑦)      =(      1−44000000    1−1−21−20010    1−1−2001−2−10    1−1−2−120010    1−1−200−12−10    121111101    121−1−1110−1    121−1−1−1−101    12111−1−10−1 )      (       𝑓0𝑓1𝑓2𝑓3𝑓4𝑓5𝑓6𝑓7𝑓8)        (2.1.16)  In vector R, 𝜌 is the fluid density, 𝑒 is the energy, 𝜀 is related to the square of energy, 𝑗𝑥 and 𝑗𝑦 represent the mass flux or momentum density, 𝑞𝑥 and 𝑞𝑦 correspond to the energy flux, and 𝑝𝑥𝑥 and 𝑝𝑥𝑦 are related to the diagonal and off diagonal components of the viscous stress tensor. An inherent advantage of the MRT model is that different macroscopic variables are relaxed with different rates. In MRT simulations, density and momentum are the only conserved moments in the system (𝜌 = 𝜌𝑒𝑞 , 𝑗 =𝑗𝑒𝑞) and the non-conserved moments at equilibrium state are given as functions of the conserved moments [27]:  𝑒𝑒𝑞 = −2𝜌 + 3(𝑢2 + 𝑣2) (2.1.17)  𝜀𝑒𝑞 = 𝜌 − 3(𝑢2 + 𝑣2) (2.1.18)  𝑞𝑥𝑒𝑞 = −𝑢 (2.1.19)  𝑞𝑥𝑒𝑞 = −𝑢 (2.1.20) 12   𝑝𝑥𝑥𝑒𝑞 = 𝑢2 − 𝑣2 (2.1.21)  𝑝𝑥𝑦𝑒𝑞 = 𝑢𝑣 (2.1.22) The collision procedure is carried out on the macroscopic variables during MRT simulations as follow:  ?̃? = 𝑒 − 𝑠1(𝑒 − 𝑒𝑒𝑞) (2.1.23)  𝜀̃ = 𝜀 − 𝑠2(𝜀 − 𝜀𝑒𝑞) (2.1.24)  ?̃?𝑥 = 𝑞𝑥 − 𝑠4(𝑞𝑥 − 𝑞𝑥𝑒𝑞) (2.1.25)  ?̃?𝑦 = 𝑞𝑦 − 𝑠6(𝑞𝑦 − 𝑞𝑦𝑒𝑞) (2.1.26)  ?̃?𝑥𝑥 = 𝑝𝑥𝑥 − 𝑠7(𝑝𝑥𝑥 − 𝑝𝑥𝑥𝑒𝑞) (2.1.27)  ?̃?𝑥𝑦 = 𝑝𝑥𝑦 − 𝑠8(𝑝𝑥𝑦 − 𝑝𝑥𝑦𝑒𝑞) (2.1.28) Where the symbol ~ represents the post collision variable. Different collision frequency rates 𝑠𝑖 are used for different moments that help damping the acoustic modes in a simulation more effectively than single relaxation time model. Before the streaming step happens, we need to construct the post collision distribution function matrix ?̃? = (𝑓0, 𝑓1, 𝑓2, 𝑓3, 𝑓4, 𝑓5, 𝑓6, 𝑓7, 𝑓8)𝑇based on the post collision moment matrix ?̃?.  ?̃? = 𝑀−1?̃? (2.1.29) In practice, the vector form of the collision procedure on macroscopic variables is written as:  ?̃? = 𝑅 − ?̂?(𝑅 − 𝑅𝑒𝑞) (2.1.30) Combining the two equations (2.1.29) and (2.1.30) we have:  ?̃? = 𝐹 −𝑀−1?̂?(𝑅 − 𝑅𝑒𝑞) (2.1.31) After the collision step is performed based on equation (2.1.31), the streaming step is carried out in the same way as in SRT model. It is theoretically shown by Lallemand and Luo [28] that one can set 𝑠7 = 𝑠8 = 1/𝜏 in order to obtain the same value of shear viscosity in MRT and SRT models.  In general, there is more flexibility in choosing other relaxation parameters in MRT model. The value of s0, s3, s5 does not affect the convergence since density (𝜌) and momentum (𝑗) are conserved moments. A good rule of thumb for the relaxation parameters s1, s2, s4, s6 is to choose values slightly greater than 1 [27]. We can completely recover the SRT BGK model for incompressible flow by setting:  𝑠1 = 𝑠2 = 𝑠4 = 𝑠6 = 𝑠7 = 𝑠8 =1𝜏  (2.1.32) When simulating high Reynolds flows using lattice Boltzmann SRT model, the solution field for velocity and pressure can show spatial oscillations in the regions around the stagnation point and sharp 13  corners where large gradients are observed. Specifically near a sharp convex corner, there are local points of singularity for velocity and pressure, so large gradients can be seen in the pressure field. It is necessary to use sufficient lattice resolution in these regions, otherwise spatial oscillations are observed which, depending on the geometry of the flow, can propagate further into the bulk and cause instability. These spatial oscillations also adversely affect the convergence rate speed. Multi relaxation time model shows better computational stability since the relaxation of different kinetic moments are separated. Lallemand and Luo [28] proved using the linearized analysis that MRT model will also result in second order spatial accuracy as the SRT model does.    It is worth mentioning that SRT and MRT models behave similarly in the long wavelength limit for macroscopic variables. However, when there are high gradient regions in the flow, i.e., substantial short wavelength components, it is expected that MRT model performs with better stability due to separation of relaxation modes. The superior stability of MRT model when simulating high Reynolds cavity flow is investigated in section (3.2). It is also observed that choice of boundary condition can noticeably improve the numerical stability when coupled with MRT model. The MRT model has significant advantages compared to SRT model when handling the geometric singularities, since one can independently adjust the bulk viscosity in MRT model. In general, the MRT model can find the solution field for (u, v, p) with much less oscillations near a singularity. 2.1.3 Chapman-Enskog expansion  Chapman-Enskog expansion is applied to the lattice Boltzmann evolution equation in order to derive a set of partial differential equations in terms of variables 𝜌 and 𝜌?⃗⃗? that approximate Navier-Stokes equations to certain orders of lattice Mach number. The evolution equation of lattice Boltzmann BGK model is written as:  𝑓𝑖(𝑥 + 𝑒𝑖⃗⃗⃗ ⃗∆𝑡, 𝑡 + ∆𝑡) − 𝑓𝑖(𝑥, 𝑡) = −1𝜏(𝑓𝑖(𝑥, 𝑡) − 𝑓𝑖𝑒𝑞(𝑥, 𝑡)) (2.1.33) If we expand the variable 𝑓𝑖(𝑥 + 𝑒𝑖∆𝑡, 𝑡 + ∆𝑡) around 𝑓𝑖(𝑥, 𝑡) using Taylor expansion series, the equation (2.1.33) can be written as below to the first order in time and space:       ∆𝑡 (𝜕𝑓𝑖(𝑥, 𝑡)𝜕𝑡+ 𝑒𝑖⃗⃗⃗ ⃗. ∇𝑓𝑖(𝑥, 𝑡)) + ⋯ = −1𝜏(𝑓𝑖(𝑥, 𝑡) − 𝑓𝑖𝑒𝑞(𝑥, 𝑡)) (2.1.34) Summing equation (2.1.34) over all lattice directions and using equations (2.1.12) and (2.1.13) we can retrieve the mass conservation equation in the limit of ∆𝑥 → 0 and ∆𝑡 → 0. The right hand side is set to be zero for mass conservation (∑𝑓𝑖(𝑥, 𝑡) = ∑𝑓𝑖𝑒𝑞(𝑥, 𝑡) at all lattice nodes).  (𝜕∑𝑓𝑖(𝑥, 𝑡)𝜕𝑡+ ∇.∑𝑒𝑖⃗⃗⃗ ⃗𝑓𝑖(𝑥, 𝑡)) +⋯ = 0 → (𝜕𝜌𝜕𝑡+ ∇. 𝜌?⃗⃗?) + ⋯ = 0 (2.1.35) In order to retrieve the Navier-Stokes momentum conservation equations, it is necessary to write the Taylor expansion up to second order terms in space in order to retrieve the shear stress terms in the Navier-Stokes equations [29]. The coefficients used in the equilibrium distribution function 𝑓𝑖𝑒𝑞(𝑥, 𝑡) 14  and the relaxation parameter 𝜏 are then calculated to recover the momentum conservation equations. For the BGK collision operator and using the D2Q9 lattice structure the equilibrium distribution function 𝑓𝑖𝑒𝑞(𝑥, 𝑡) is calculated as equation (2.1.7) and the relaxation parameter 𝜏 is derived as equation (2.1.6). From a kinetic point of view, we can decompose the distribution functions into two equilibrium and non-equilibrium parts. The equilibrium part represents the physical phenomenon that all molecules tend to approach a local equilibrium state as defined by the collision operator. The non-equilibrium part describes the local oscillations due to molecular forces such as van der Waals, as well as other diverging motions from the equilibrium state.   𝑓𝛼(𝑥, 𝑡) = 𝑓𝛼𝑒𝑞(𝑥, 𝑡) + 𝑓𝛼𝑛𝑒𝑞(𝑥, 𝑡) (2.1.36) Looking from a macroscopic point of view, these diverging motions appear in the form of temporal and special derivative of macroscopic variables such as density and momentum. The approximation of distribution functions is the third step during Chapman-Enskog analysis. Approximating the distribution function 𝑓𝛼(𝑥, 𝑡) in the left hand side of equation (2.1.34) with the equilibrium distribution function 𝑓𝛼𝑒𝑞(𝑥, 𝑡) to zero order will result in:   𝑓𝛼𝑛𝑒𝑞(𝑥, 𝑡) = −𝜏∆𝑡 (𝜕𝜕𝑡+ 𝑒𝛼 . ∇) 𝑓𝛼𝑒𝑞(𝑥, 𝑡) + ⋯  (2.1.37) Asymptotic analysis shows that non-equilibrium distribution functions can be considered as a small perturbation from molecular equilibrium state [30]. The perturbation is with the scale of Knudsen number O(𝜀), defined as the mean free path length for molecular interactions over the physical characteristic length of the problem. Approaching the value of 1 for Knudsen number implies that mean free path is comparable to the characteristic physical length, therefore the continuum assumption of fluid mechanics is not valid anymore. In order to hold the continuum assumption, small Knudsen number values have to be chosen for simulation. Using equation (2.1.37) we can estimate the distribution function 𝑓𝛼(𝑥, 𝑡) = 𝑓𝛼𝑒𝑞(𝑥, 𝑡) + 𝑓𝛼𝑛𝑒𝑞(𝑥, 𝑡) up to first order when recovering the Navier-Stokes momentum equations.   In section (2.2), we have looked into two critical issues that has to be taken care of during a lattice Boltzmann simulation. First, we will discuss how to convert the physical variables with standard units (SI) to lattice units for problem definition in lattice system. Then, we briefly mention the sources of error in a lattice Boltzmann simulation and the chapter is concluded in section (2.3) with a general overview of lattice Boltzmann and possible shortcomings with respect to swallowing simulations. 2.2 Physical and lattice systems  Physical flow problems can be described using physical units of length, time and mass. These physical units are different from units internally used inside lattice Boltzmann algorithm. In order to have a practical lattice Boltzmann simulation, we need to perform a mapping between the lattice simulation units and the physical units using characteristic variables L0, T0, and M0 for length, time, and mass. 15   𝑥𝐿0 → 𝑥  ,   ?̂?𝑇0 → 𝑡   ,   ?̂?𝑀0 → 𝑚  (2.2.1)  In equation (2.2.1), ?̂? is the lattice time and 𝑡 is the physical time. The mapping is required to convert the macroscopic variables such as velocity from lattice spacing per time step in LB simulation to physical units such as meters per second. The characteristic length variable L0 is measured based on the choice of lattice resolution for the simulation. Take the cavity flow geometry as an example with a side length of 1 mm and a lattice resolution of 20 nodes across the side. The physical spatial step in the simulation would be ∆𝑥 = 5 × 10−5m and assuming that choice of lattice spacing ∆?̂? = 1 is chosen for the simulations, the characteristic length L0 would be equal to 𝐿0 = ∆𝑥/∆?̂? = 5 × 10−5 meters.     The characteristic mass M0 is measured based on the density of the fluid. Assume that water with a density of 𝜌 = 998.3 𝑘𝑔.𝑚−3 at T=20o C is used for the simulation. The choice of initial lattice density is dictated by the initial pressure boundary condition through the equation of state (equation (2.1.14)). However, for steady state simulations the initial density of ?̂? = 1 is chosen for all lattice nodes. The characteristic mass M0 is then measured by the equation 𝜌 = 𝑀0𝐿0−3?̂? and is equal to 𝑀0 = 1.25 ×10−10 kg.  The characteristic time parameter is determined based on the choice of user for relaxation parameter 𝜏. The kinematic viscosity of water is 𝜗 = 1.004 × 10−6 𝑚2𝑠−1 at T=20o C. In the lattice BGK model, kinematic viscosity is measured by the choice of relaxation parameter and lattice sound speed ?̂? = (𝜏 − 0.5)?̂?𝑠2∆?̂?. By definition, there is a restriction on the choice of relaxation parameter 𝜏 >0.5 in order to obtain positive lattice viscosity values. For the sake of simplicity, usually spatial and time step sizes of ∆𝑥 = ∆?̂? = 1 are chosen for lattice Boltzmann simulations. Through a dimensional analysis on the viscosity we can see that 𝜗 = ?̂?𝐿02𝑇0−1. Combining the two previous equations with the assumption of ∆?̂? = 1 will result in 𝜗 = (𝜏 − 0.5)?̂?𝑠2𝐿02𝑇0−1. The BGK model dictates that lattice sound speed is equal to ?̂?𝑠 = 1/√3, therefore the characteristic time parameter T0 is dependent on the choice of relaxation parameter 𝜏, for example a value of 𝜏 = 1 will result in 𝑇0 = 4.15 × 10−4s. Therefore the physical time step considered in the simulation is equal to ∆𝑡 = 4.15 × 10−4 seconds. We can decrease the time step size using smaller values for relaxation time. For example, A value of 𝜏 = 0.51 will result in ∆𝑡 = 9.30 × 10−6 seconds.  The dependency between time step and the choice of relaxation time place practical restrictions in LB simulations. In the cavity example, physical time step is measured using the choice of lattice resolution and relaxation parameter. When simulating transient flows, it is very important to choose a time step much smaller than the time scale of the flow evolution. Suppose that the cavity example described before had a side of 10 cm instead of 1mm. it is not practically possible to maintain the same spatial resolution ∆𝑥 = 5 × 10−5m since it is computationally very expensive, so we choose a new value of spatial resolution ∆𝑥=0.05 cm. If we use the same value for relaxation time 𝜏=1 and assuming all other variables are the same, the new physical time step would be equal to ∆𝑡=415 s. It is obvious that this new value of physical time step is extremely large for simulating transient flows. In order to have smaller time step size, we can either use a much finer lattice resolution which will result in large computational time, or use a different value for the relaxation parameter in a way that (𝜏 − 0.5) is small enough. However very small values for (𝜏 − 0.5) implies inviscid flow and very long simulation runtime is expected due to slow rate of energy dissipation. On the other hand, the simulation accuracy 16  decreases for large values of relaxation parameter. A good rule of thumb is to use 0.5 < 𝜏 < 3 for lattice Boltzmann BGK simulations [31].      One should be careful about the low Mach number assumption when choosing the relaxation parameter. The weak compressibility of the lattice Boltzmann BGK model leads to small density gradients in the flow domain which is in contradiction with the assumption of incompressible fluid. Using the Chapman-Enskog expansion, we can obtain the incompressible Navier-Stokes equations accurate up to 𝑂(𝑀𝑎2) for continuity equation and up to 𝑂(𝑀𝑎3) for the momentum equations written as follow:  {𝜕𝑡𝑢 + (𝑢. ∇)𝑢 = −∇𝑝 + 𝜗∆𝑢 + 𝑔 + 𝑂(𝑀𝑎3)∇. 𝑢 = 0 + 𝑂(𝑀𝑎2)  (2.2.2) The above system of equations is biased from incompressible Navier-Stokes equations due to the compressibility effect of lattice Boltzmann BGK model. The lattice time step has to be small with respect to the lattice space step to reduce this error. The time step is usually chosen as ∆𝑡 = 𝑂(∆𝑥2). Practically, we can decrease the lattice Mach number of a simulation using higher lattice resolution or smaller relaxation parameter.  2.2.1 Sources of error in LBM  The main sources of error in a lattice Boltzmann simulation are the spatial discretization error, temporal discretization error, and the compressibility error.   𝐸𝑡𝑜𝑡𝑎𝑙 = 𝐸∆𝑥 + 𝐸∆𝑡 + 𝐸𝑀𝑎 = 𝑂(∆𝑥2) + 𝑂(∆𝑡2) + 𝑂((∆𝑡∆𝑥)2)  (2.2.3) Similar to any other numerical scheme, the spatial and temporal discretization error exist in lattice Boltzmann simulations. The resolution of lattice Boltzmann simulations are usually chosen based on the geometry and characteristic length of the flow domain. The time step is then measured based on the choice of relaxation parameter 𝜏. Assuming that ∆𝑡 = 𝑂(∆𝑥2) is chosen for the LB simulations, the spatial discretization error is the dominant source of error according to equation (2.2.3), while for ∆𝑡 =𝑂(∆𝑥), the discretization error remains dominant for large values of ∆𝑥 and ∆𝑡 and switches to the compressibility error as the main contributing factor for small choices of space and time step. When ∆𝑡 = 𝑂(∆𝑥0), the compressibility effect is merely the dominant source of error and increases as we increase the lattice resolution.       2.3 Swallowing simulation characteristics and LBM accuracy  Generally speaking, accuracy and stability of the lattice Boltzmann BGK scheme is very sensitive to Reynolds number of the flow and the choice of relaxation parameter. Increasing the relaxation parameter will result in higher Mach number flows, i.e. more compressibility error. On the other hand, approaching to lower limit of 𝜏 → 1/2, i.e. inviscid flow, also causes instability in the high gradient regions of the flow due to low speed energy dissipation. With regards to swallowing simulations and their complexity, it is expected that extremely fine lattice resolutions are required to have a stable 17  solution. The food bolus will experience high gradients in velocity when moving through the pharynx. Very fine lattice resolution with an optimum choice of relaxation time is required for simulations of this complexity, therefore it is insightful to study the stability and accuracy of lattice Boltzmann method on benchmark flows.  Test case studies have to be done to investigate the stability of lattice Boltzmann method at different Reynolds number.  In a swallowing simulation, Reynolds number can rise to very high values for low viscosity materials. Due to shortcomings of SRT model for high Reynolds flow simulations, it is important to study multi relaxation time model and measure its stability at higher Reynolds simulations. To preserve the desired accuracy, one might have to use extremely fine lattice resolution to simulate high Reynolds flows which will result in huge computational cost. It is expected that multigrid methods can significantly decrease the computational cost while ensuring high accuracy in the simulations. However the random boundary movements of the pharynx during swallowing makes it very challenging to implement the multigrid method coupled with the moving boundary. Also, the high gradient region in a swallowing simulation is constantly changing as the bolus transports throughout the pharynx. Therefore, the fixed grid nature of the lattice Boltzmann solver makes it extremely hard to define a fixed region of high divergence for multigrid method implementation.                 18  Chapter 3. Boundary conditions Boundary conditions play a very important role for conventional numerical methods, such as finite element method, to solve any differential system. In a similar way, imposing boundary conditions accurately is crucial for lattice Boltzmann method. In this chapter we will investigate the implementation of different boundary conditions for LBM and will assess the stability and accuracy of these schemes and their limitations when dealing with flow characteristics of bolus transport during swallowing.     Lattice Boltzmann is quite different from the conventional methods in terms of boundary condition implementation, since the leading role playing on the stage is not macroscopic variables but rather microscopic distribution functions. For example when solving Navier-Stokes equations using macroscopic methods, the velocity or pressure conditions imposed on the boundaries are directly used to find the same functions on the boundary cells and then the discretized system of equations is solved for pressure and velocity in the complete domain of the problem. When the boundary conditions are approximated on the boundary the same way as that within the domain by a certain conventional numerical method, the problem ends up with solving a linear or nonlinear algebraic system. On the other hand, in lattice Boltzmann method first we have to reconstruct the microscopic distribution functions on the boundaries in a way that mass and momentum are conserved on lattice grid points and then implement the collision and streaming steps to propagate information into the domain. Reconstructing the unknown distribution functions on the boundary based on the macroscopic imposed conditions is one of the significant difficulties of lattice Boltzmann method. In this chapter, we will focus on how to impose Dirichlet conditions which specify the value of the function on the boundary, e.g. solid walls, and Neumann conditions which specify the normal derivative of the function on the boundary for example outflow condition with zero normal velocity gradient. We have looked into four different boundary conditions for the purpose of modeling straight solid walls, curved solid walls, constant pressure and unidirectional flow, and unidirectional flow with zero normal velocity gradient. The two later choices can be good approximations for inlet and outlet boundaries in preliminary swallowing simulations. In section (3.2) we have looked into the so called bounce back scheme which is the most commonly used solid wall boundary treatment in lattice Boltzmann simulations. The bounce back scheme is based on the idea that the boundary nodes send back what streams into the solid in the opposite direction and back into the bulk during the collision step. This scheme is mostly used to simulate straight wall conditions and will give first order accuracy when the lattice nodes are exactly on the boundary wall (full-way bounce back) and will result in second order accuracy when the wall is placed halfway between the lattice grid points at the boundary (half-way bounce back). When modeling curved boundaries using bounce back scheme, the wall is treated as a series of stair wise steps which will result in numerical instability when modeling high Reynolds flows. We will examine the stability of the scheme on lid driven cavity flow with respect to the Reynolds and Mach number of the flow. This scheme can be used to model simplified models of swallowing when the hard palate is taken as a straight wall. We will also look into the implementation of bounce back boundary treatment with both MRT and SRT models and will discuss the superior stability of MRT compared to SRT. MRT model is shown to be more stable than the SRT model theoretically by Lallemand and Luo analysis [28], but the 19  influence of boundary condition implementation is not considered in their work. When dealing with complex high Reynolds flows such as low viscous bolus transport during swallowing, the choice of a stable boundary condition which is compatible with MRT model becomes important. To investigate this issue, The two dimensional lid driven cavity benchmark is modeled using both bounce back and non-equilibrium extrapolation boundary treatments and the effect of boundary condition on the stability of the solution at higher Reynolds numbers are discussed.  In section (3.3), we will look into open boundary choices as inlet and outlet conditions. Zou/He scheme is introduced and studied for implementation of Dirichlet conditions on pressure and velocity. A 2D pressure driven channel flow is chosen as the benchmark for stability and accuracy measurements. The extrapolation scheme is introduced in section (3.3) in order to model outlet condition as unidirectional flow with zero normal velocity gradient. A 2D 2:1 contracting channel flow is used as the benchmark with constant velocity inlet and zero velocity gradient outlet conditions. Section (3.4) focuses on the issue of curved boundary condition implementation and force evaluation during lattice Boltzmann simulations. The improved Filippova-Hanel boundary scheme is introduced and the stability and accuracy of the force evaluation scheme and the flow solver are tested on the flow over a bounded cylinder placed in a 2D channel. Finally Section (3.5), looks into the transient channel flow benchmark with a periodic pressure inlet condition. First, an iterative initialization procedure is introduced that solves the position equation for pressure whit a divergent free initial velocity field. Then, the compressibility effect of lattice Boltzmann method is compared for LBE and incompressible LBE solutions. 3.1 Boundary conditions in a swallowing simulation Implementing the boundary conditions accurately in a swallowing simulation is crucial since the flow is mostly boundary driven. The mass and momentum should be accurately conserved at boundary lattice nodes as well as the interior nodes in order to have a stable solution. The complex curved boundary of the pharynx and the high Reynolds number in a swallowing simulation dictate the need for curved boundary schemes that preserve the integrity of the flow geometry. We will investigate the capabilities and limitations of Filippova and Hanel boundary treatment to model flows with such characteristics. We have used flow over a cylinder symmetrically placed in a 2D channel as our benchmark. Limitations on lattice resolution per diameter of the cylinder are investigated and the stable range of Reynolds number is measured.  The inflow and outflow boundary conditions in a swallowing simulation are not very well defined. In a simplified model the inflow condition can be approximated with a Dirichlet constant atmospheric pressure condition. To model the outflow condition, we can extend the pharynx with a straight tube so that the flow is approximately unidirectional and fully developed at the outlet of the tube [29]. In this case, we can use unidirectional (𝑢𝑥 = 0) and zero axial velocity gradient condition (𝜕𝑢𝑦𝜕𝑦= 0) at the outlet. We investigate the “Zou/He boundary treatment” for implementation of the inflow condition and the “extrapolation scheme” to implement the zero velocity gradient at the outlet. 20  3.2 Solid boundary treatments In this work, we have used the D2Q9 lattice structure In order to illustrate the idea of different kinds of boundary conditions and to assess their accuracy and stability. Consider the node Q in Figure (3.1). Since there is no streaming from the solid into the bulk, the distribution functions f3, f6, f7 are unknown and have to be reconstructed. The principle for tackling boundary conditions is that mass and momentum or even energy have to be conserved on the boundary so as to be consistent with that in the bulk region. In microscopic point of view, all the distribution functions on the boundary, which can be separated as equilibrium part related to density and velocity and non-equilibrium part represented also by velocity derivatives, are required to satisfy the conservation law.    Figure 3.1- Configuration of D2Q9 lattice structure with solid and fluid nodes sketched as solid and hollow points  Using equations (2.1.12) and (2.1.13) we have that:  {𝜌 = 𝑓0 + 𝑓1 + 𝑓2 + 𝑓3 + 𝑓4 + 𝑓5 + 𝑓6 + 𝑓7 + 𝑓8𝜌𝑢𝑥 = 𝑓1 − 𝑓3 + 𝑓5 − 𝑓6 − 𝑓7 + 𝑓8𝜌𝑢𝑦 = 𝑓2 − 𝑓4 + 𝑓5 + 𝑓6 − 𝑓7 − 𝑓8  (3.2.1) Assuming that both velocity components are given on the boundary node Q, there are 4 unknowns in the above system (𝜌, 𝑓3, 𝑓6, 𝑓7) and three equations to solve for these unknowns, so another equation is needed to solve the linear system for all unknown variables. The fourth equation will be chosen in a way that is consistent with the physics of the problem and is the topic of discussion in the following sections. Density on the node Q can be found using the consistency of the first two equations in (3.2.1). Adding the first two equations of the system and solving for 𝜌 will result in:  𝑢𝑥 = −1 +1𝜌(𝑓0 + 𝑓2 + 𝑓4 + 2(𝑓1 + 𝑓5 + 𝑓8))  (3.2.2) In general, it is always possible to find the third unknown macroscopic variable 𝜌, 𝑢𝑥, or 𝑢𝑦 on a straight boundary when two of them are given in the form of Dirichlet boundary conditions. Keep in mind that density can be written in terms of pressure as 𝑝 = 𝜌𝑐𝑠2 by the results for ideal gas dynamics. 21  3.2.1 Manipulation of corner nodes Special attention should be given to the corner nodes when implementing boundary conditions in LBM simulations since the number of unknown distribution functions will increase to five instead of three at these points. Take the upper right corner node in Figure (3.1) as an example. In addition to 𝑓3, 𝑓6, 𝑓7 two other distribution functions 𝑓4, 𝑓8 are also streaming from the solid into the bulk, therefore two more equations are needed to be able to solve the system of equations for the corner nodes. The macroscopic approach to this problem is to extrapolate the unknown density or velocities on the corner nodes based on their values on the neighboring nodes. In a similar way, the microscopic approach also uses extrapolation of the extra unknown distribution functions from neighboring nodes. 3.2.2 Bounce-back scheme Bounce back is the simplest boundary treatment that we investigate in this thesis and it is mostly used to simulate no slip condition on straight walls. The basic idea behind this scheme is to send back the distribution functions that stream from the flow domain into the solid at the opposite direction. Figure (3.2) demonstrates the implementation of bounce back scheme on the boundary node “P” placed above the bottom solid wall at time step “t”:  Figure 3.2- Schematic view of the halfway bounce back boundary treatment during time step “t” Since the velocity is equal to zero on the wall, conservation of momentum requires that:  {𝜌𝑢𝑥 = 𝑓1 − 𝑓3 + 𝑓5 − 𝑓6 − 𝑓7 + 𝑓8 = 0 → 𝑓1 + 𝑓5 + 𝑓8 = 𝑓3 + 𝑓7 + 𝑓6𝜌𝑢𝑦 = 𝑓2 − 𝑓4 + 𝑓5 + 𝑓6 − 𝑓7 − 𝑓8 = 0 → 𝑓2 + 𝑓5 + 𝑓6 = 𝑓4 + 𝑓7 + 𝑓8 (3.2.3) Bounce back is an obvious answer to equation (3.2.3). Conservation of energy is related to the second order moment of lattice velocity and is not taken into account by the bounce back scheme, so this method cannot be used for thermal flow simulations. It is shown that bounce back achieves second order accuracy when the wall is placed halfway between the two lattice grid points on the solid and fluid regions [32]. 22  We have looked into the stability and accuracy of this scheme on a 2D lid-driven cavity benchmark. The configuration of the benchmark consists of a two-dimensional square cavity whose top plate moves from left to right, while the other three walls are fixed. The most important observation is the limitation of this scheme on Reynolds number when used with both SRT and MRT models. Different Reynolds numbers ranging from 100 to 3200 are chosen and the stability of the scheme is investigated as a function of Reynolds number and the collision frequency. The Reynolds number in the lattice system is defined as 𝑅𝑒 = (𝑢0𝑁𝑦)/𝜗 where 𝑢0 is the lattice velocity of the moving lid, 𝑁𝑦 is the number of lattice unit lengths or the lattice resolution on the wall boundaries, and 𝜗 is the lattice fluid viscosity which is related to the collision frequency by the following equations: 𝑐𝑜𝑙𝑙𝑖𝑠𝑖𝑜𝑛 𝑓𝑟𝑒𝑞𝑢𝑒𝑛𝑐𝑦 ∶  𝜔 =1𝜏  𝜗 = (1ω−12) 𝑐𝑠2𝛿𝑡 (3.2.4) The lattice sound speed, 𝑐𝑠, is equal to 1√3 on a D2Q9 lattice structure and the most common choice of lattice time step is equal to unity in LBM simulations. Stability of lattice Boltzmann method is sensitive to the choice of collision frequency, therefore it is best to keep this parameter constant while increasing the lattice resolution in order to decrease the Mach number and have a more accurate and stable simulation. It is observed experimentally that in order to have a stable solution for the cavity flow at Re=100 or Re=3200 a minimum lattice resolution of 𝑁𝑦 = 10 or 𝑁𝑦 = 240 is required respectively when collision frequency is set to 𝜔=1.66. As a result, the computational cost of simulation increases greatly with increasing Reynolds number. To overcome this limitation, we will use MRT models and non-equilibrium extrapolation boundary conditions later in this chapter which will allow us to raise the Reynolds number significantly. Figure (3.3) shows the streamline configuration in a cavity flow for Reynolds number Re=100 and Re=3200. Lattice Boltzmann BGK model is used with a single relaxation time parameter coupled with bounce back boundary treatment for the simulations. The data for center of vortices ‘y’ coordinates are provided in Table (3.1) and are compared with the ones obtained by Ghia [33]. The good agreement suggests that lattice Boltzmann BGK model is capable of simulating rotating flows, which is a characteristic of pharyngeal flows, but is highly limited by the Reynolds number.    (a) (b) Figure 3.3- Cavity streamline configuration (a) Re=100, collision frequency=1.66. (b) Re=3200, collision frequency=1.66 23  Table 3.1- Comparison of the center of vortex ‘y’ coordinates for cavity flow between Ghia and LBM at different Reynolds and Mach numbers Re & Mach Center of Cavities (y) coordinates Major Vortex LBM GHIA Re=100 , Mach=0.014 0.7460 0.7344 Re=100 , Mach=0.058 0.7387 0.7344 Re=100 , Mach=0.173 0.7375 0.7344 Re=100 , Mach=0.289 0.7368 0.7344 Re=400 , Mach=0.058 0.6064 0.6055 Re=1000 , Mach=0.144 0.5656 0.5625 Re=3200 , Mach=0.231 0.5405 0.5469  Bottom Left Vortex Re=100 , Mach=0.014 0.0378 0.0391 Re=100 , Mach=0.058 0.0369 0.0391 Re=100 , Mach=0.173 0.0271 0.0391 Re=100 , Mach=0.289 No vortex 0.0391 Re=400 , Mach=0.058 0.0486 0.0469 Re=1000 , Mach=0.144 0.0802 0.0781 Re=3200 , Mach=0.231 0.1211 0.1094  Bottom Right Vortex Re=100 , Mach=0.014 0.0604 0.0625 Re=100 , Mach=0.058 0.0739 0.0625 Re=100 , Mach=0.173 0.0799 0.0625 Re=100 , Mach=0.289 0.0894 0.0625 Re=400 , Mach=0.058 0.1248 0.125 Re=1000 , Mach=0.144 0.1137 0.1094 Re=3200 , Mach=0.231 0.0849 0.0859   24  The two results show good agreement in Table (3.1). It is worth to mention that error increases for higher Mach numbers. The bottom left vortex cannot be detected with Re=100 & Ma=2.89 due to low lattice resolution. Figure (3.4) demonstrates the ‘x’ and ‘y’ centerline velocity profiles for Reynolds numbers 100, 400, 1000, and 3200. LBM results are compared with the Results obtained by Ghia. The solid circle and square symbols represent the data obtained by Ghia.   (a) (b)   (c) (d) Figure 3.4- Non-dimensional velocity along the vertical and horizontal centerlines (a) horizontal velocity at Re=100 and Re=1000 (b) horizontal velocity at Re=400 and Re=3200 (c) vertical velocity at Re=100 and Re=1000 (d) vertical velocity at Re=400 and Re=3200. 25   The data in Figure (3.3), (3.4), and Table (3.1) validates the single relaxation time BGK model capability to accurately solve the cavity flow using halfway bounce back boundary treatment, if the Mach number is much smaller than unity, i.e. high lattice resolution is used. Computationally speaking, the CPU time needed to simulate cavity flow at Re=3200 with a 800*800 lattice resolution, minimum required resolution when collision frequency is set to be 1.66, is 245221 seconds (about 70 hours) on a desktop computer with Intel(R) core™ i7 CPU at 3.30 GHz. Therefore, simulating high Reynolds flows such as swallowing low viscosity fluids requires extremely fine lattice resolution and huge computational time. It is worth mentioning that one of the great advantages of LBM is that it’s easy to parallelize the implementation of the streaming step which will reduce the computational time significantly. The other approach to simulate high Reynolds flows is to use Multi Relaxation Time models. We have compared the stability of SRT and MRT models on the cavity flow later in this section. Figure (3.5) shows the computational time as a function of lattice resolution. The collision frequency is set to be 1.66 in all simulations.     (a) (b) Figure 3.5- CPU time as a function of lattice resolution (a) Re=100 (b) Re=400  As can be seen in Figure (3.5), Three different lattice resolutions N=100, N=200, and N=400 are used for both Reynolds numbers 100 and 400. Both curve equations suggest that CPU time is proportional to the lattice resolution by a power of three, 𝐶𝑃𝑈𝑡𝑖𝑚𝑒 ∝ 𝑁3, which is consistent with results obtained in [29]. On the other hand, the flow Reynolds number does not affect the computational time noticeably compared to the choice of lattice resolution and collision frequency as suggested by the data represented in Table (3.2). Table 3.2- Computational time comparisons for cavity flow simulations at Re=100 and Re=400 Reynolds Lattice resolution 100*100 200*200 400*400 100 t=263 s t=1996 s t=19061 s 400 t=280 s t=2247 s t=23549 s  y = 0.0002x3.08971101001000100001000001 10 100 1000CPU time (s)Lattice Resolution (N)Re=100y = 0.0001x3.1971101001000100001000001 10 100 1000CPU time (s)Lattice Resolution (N)Re=40026  The convergence criteria is chosen based on the oscillation level of the average Mach number during simulations. The flow will reach steady state solution when the Mach number oscillations are completely damped and it reaches the average value of Mach number for the steady state solution. Figure (3.6) shows the average Mach number oscillations for different lattice resolutions N=100, 200, and 400 and different Reynolds numbers Re=100 and Re=400. Figure (3.6.a) suggests that flow with Re=100 reaches the steady state solution approximately at iterations N=30000, N=90000, and N=200000 with lattice resolution of 100*100, 200*200, and 400*400 respectively.   (a) (b) Figure 3.6- Convergence behavior of the average lattice Mach number for cavity flow (a) Re=100 (b) Re=400 According to the computational cost and stability shortcomings of SRT model, the need for more stable models when dealing with high Reynolds flows become apparent. Due to simplicity, the BGK equation has become the most popular lattice Boltzmann equation despite its well-known deficiencies such as flow simulation at high Reynolds numbers. It is theoretically shown by Lallemand and Luo analysis that using multi relaxation time parameters will increase the stability of lattice Boltzmann at higher Reynolds flows. Figure (3.7) compares the stability of SRT and MRT models when used with bounce back boundary condition on the lid driven cavity benchmark. High Reynolds number, low lattice resolution, and the choice of collision frequency can be the source of instability during simulations. The predicted minimum allowed lattice resolutions for a stable solution at different collision frequencies are presented on the vertical axis, as a function of Reynolds number on the horizontal axis. It is observed that stability of MRT model is only slightly higher than SRT model when used with bounce back boundary condition. Figure (3.8) performs the same centerline velocity comparison as done in Figure (3.4) to validate our MRT solver. Later we will investigate the stability of MRT model when used with non-equilibrium extrapolation boundary treatment on the cavity benchmark. The superior stability of MRT model is obvious when used with non-equilibrium extrapolation boundary treatment.  As described in section (2.1.2), Multi relaxation time model uses different relaxation parameters for different moments of distribution functions. Similar to the implementation of bounce back for SRT 27  model, distribution functions streaming from the fluid into the solid region are bounced back in the opposite direction before the collision step occurs. The macroscopic variables vector and the equilibrium vector in equations (2.1.15) are then reconstructed using the new distribution functions during collision step. Finally updated distribution functions are passed to the neighbor nodes during streaming step. We have looked into the stability of MRT model as a function of relaxation parameters 𝑆7 = 𝑆8 and Reynolds number. Equation (2.1.16) is used and values of relaxation parameters 𝑆0, 𝑆1, 𝑆2, 𝑆3, 𝑆4, 𝑆5 , 𝑆6 are taken from [34].   Figure 3.7- Stability comparison between SRT and MRT models on the cavity flow Figure (3.7) demonstrates the stable and unstable regions for LB SRT and MRT models when modeling the lid driven cavity flow using halfway bounce back boundary treatment. The minimum stable lattice resolution is shown on the vertical axis and has been measured for four different Reynolds numbers ranging from 100 to 3200. The effect of collision frequency is also studied for three different values of 1.66, 1, and 0.5. The upper region of each line will result in a stable solution with macroscopic variables converging to steady values. The MRT model is only slightly more stable that the SRT model for all collision frequency values. Also, notice the linear behavior between lattice resolution and Reynolds number at constant collision frequencies. This behavior is consistent with the results obtained 28  by E. Aslan [35]. The slope of the line fitted to the SRT model data points are 0.074, 0.238, and 1.313 for collision frequencies of 1.66, 1.0, and 0.5 respectively. This sheds light on the fact that minimum stable lattice resolution is much more sensitive to Reynolds number at higher collision frequencies. It’s worth mentioning that computational cost of the MRT model is slightly higher than SRT (<10%) when same collision frequency and lattice resolution are used. Figure (3.7) suggests that MRT model is only slightly more stable than SRT model when used with halfway bounce back boundary condition.    (a) (b) Figure 3.8- Centerline velocity comparison between MRT lattice Boltzmann and Ghia at Re=100 (a) vertical centerline (b) horizontal centerline  By coupling the multi-relaxation-time lattice Boltzmann method with the non-equilibrium extrapolation boundary treatment, the stability of the model is observed to be improved when dealing with high Reynolds number flows. The stability of an MRT lattice Boltzmann simulation is determined based on the choice of collision frequency, lattice resolution, and boundary condition implementation. The choice of non-equilibrium boundary treatment will allow to use collision frequency values very close to the inviscid limit 𝜔 → 2. In order to illustrate the superior stability of the non-equilibrium extrapolation boundary scheme compared to halfway bounce back scheme, we study the cavity flow benchmark with Reynolds number as high as Re=10000.   The non-equilibrium extrapolation scheme is based on the idea of decomposing the unknown distribution functions at the boundary node into their equilibrium and non-equilibrium parts. Suppose that node ‘b’ is placed on the boundary and node ‘f’ is the closest neighbor node in the fluid region as shown in Figure (3.9). The distribution functions f3, f6, and f7 are unknown during simulation, since they’re streaming from the solid into the fluid region. 29   Figure 3.9- Unknown distribution functions configuration for boundary node ‘b’   The idea is to decompose the unknown distribution functions (for example f3 in Figure (3.9)) into their equilibrium and non-equilibrium parts as follow:  𝑓3(𝑥𝑏 , 𝑡) = 𝑓3𝑒𝑞(𝑥𝑏, 𝑡) + 𝑓3𝑛𝑒𝑞(𝑥𝑏 , 𝑡)  (3.2.5)  In two dimensional space, generally two of the macroscopic variables u, v, or P will be specified on the boundaries in the form of Dirichlet or Neumann boundary conditions and one will remain unknown. In a cavity flow, pressure is not specified on the boundary walls. Non-equilibrium extrapolation scheme uses the lattice density of the nearest fluid neighbor node ‘f’ to calculate the equilibrium distribution function 𝑓3𝑒𝑞(𝑥𝑏 , 𝑡) at the boundary node. The non-equilibrium part 𝑓3𝑛𝑒𝑞(𝑥𝑏 , 𝑡) is also approximated by 𝑓3𝑛𝑒𝑞(𝑥𝑓 , 𝑡).   𝑓3𝑒𝑞(𝑥𝑏 , 𝑡) = 𝜌(𝑥𝑓 , 𝑡)𝑤3[3𝑐2𝑒3. 𝑢(𝑥𝑏 , 𝑡) +92𝑐4(𝑒3. 𝑢(𝑥𝑏 , 𝑡))2 −32𝑐2𝑢(𝑥𝑏 , 𝑡). 𝑢(𝑥𝑏 , 𝑡)] (3.2.6)  𝑓3𝑛𝑒𝑞(𝑥𝑏, 𝑡) = 𝑓3𝑛𝑒𝑞(𝑥𝑓 , 𝑡) (3.2.7)   Using equations (3.2.6) and (3.2.7), one can calculate the unknown distribution function f3 at the boundary node ‘b’. The same procedure is applied to the other two unknown distribution functions f6 and f7.  The relaxation rates for the MRT model are chosen as s0=s3=s5=0, s1=1.1, s2=1.0, s4=s6=1.2, s7=s8=𝜔 based on the reference work [34] and the collision frequency is set to be 𝜔=1.994 in order to get a Mach number of Ma=0.05 during simulations with lattice resolution of N=200*200. Figure (3.10) shows the streamlines for cavity flow at Re=10000.  30   Figure 3.10- Streamlines for cavity flow at Re=10000 using MRT lattice Boltzmann model (𝜔 = 1.994)  Figure (3.11) compares the horizontal and vertical centerline velocity profile between data obtained by lattice Boltzmann MRT model and data obtained by GHIA. The solid black symbols represent data points presented by GHIA and the dashed velocity profile is obtained using MRT simulations. Good agreement between the two velocity data sheds light on the validity of the implemented MRT simulations. Furthermore, center of vortex coordinates are compared in Table (3.3) and good agreement is seen between the lattice Boltzmann and GHIA results. As mentioned before, the choice of collision frequency 𝑠7 = 𝑠8 = 𝜔 = 1.994 will result in a flow Mach number Ma=0.05 which is 10 times larger than the spatial step 𝛿𝑥 = 1/200 = 0.005, therefore, the compressibility error is the dominant source of error in the simulations. The flow behavior reaches the steady condition after about 300000 iterations and a computational time t=17826s (≈5 hours). The long computational time is due to the choice of collision frequency which is approaching the limit of inviscid flow 𝜔 → 2.   (a) (b) Figure 3.11- Non-dimensional centerline velocity profile comparison between MRT lattice Boltzmann and Ghia at Re=10000 (a) horizontal centerline (b) vertical centerline  31   Table 3.3- Center of vortex coordinate comparison between Ghia and MRT lattice Boltzmann for cavity flow at Re=10000 Vortices LBM(xc) GHIA(xc) LBM(yc) GHIA(yc) Primary 0.5126 0.5117 0.5297 0.5333 BR1 0.7776 0.7656 0.0590 0.0586 BR2 0.9346 0.9336 0.0713 0.0625 BL1 0.0558 0.0586 0.1695 0.1641 BL2 0.0195 0.0156 0.0239 0.0195 T 0.0671 0.0703 0.9144 0.9141  In this section, we have looked into the stability and accuracy of the lattice Boltzmann SRT and MRT models coupled with the half-way bounce back and non-equilibrium extrapolation boundary schemes for straight walls. The lid driven cavity flow is chosen as the benchmark and the stability of the solver is investigated as a function of lattice resolution, Reynolds number, and choice of collision frequency. The validity of the SRT model is shown on Figures (3.3) and (3.4) and Table (3.1). As mentioned before in section (2.2.1), SRT model uses a single relaxation parameter 𝜏 in order to relax all of the hydrodynamic moments towards the Maxwellian equilibrium state. This will result in instability during simulations with choices of relaxation time 𝜏 → 0.5 or equivalently 𝜔 = 1/𝜏 → 2 due to the low rate of energy dissipation in flows with high gradient regions such as pharyngeal bolus transport. In addition, choices of relaxation time 𝜏 > 3 will result in high Mach number in the simulation, i.e. high compressibility error, therefore a range of  0.55 < 𝜏 < 3 is suggested for lattice Boltzmann SRT model. The computational time of lattice Boltzmann simulations is dominantly affected by the choice of lattice resolution and collision frequency. Figure (3.5) suggests that computational time is proportional to the third power of lattice resolution (𝐶𝑃𝑈𝑡𝑖𝑚𝑒 ∝ 𝑁3) which result in very long computational time for complex flow regimes such as pharyngeal swallowing that requires very fine lattice resolution in order to accurately capture the high gradient regions. Using parallel computing (GPU) and multigrid methods are two effective approaches and can be further investigated to overcome this issue, however, implementation of multigrid method is very challenging since the high gradient flow region is moving as the bolus transports through pharynx into the esophagus. In order to overcome the issue of low energy dissipation, we have used the multi relaxation time model that relaxes different hydrodynamic moments towards the equilibrium state using different relaxation rates. This will allow to use collision frequencies very close to the inviscid limit 𝑠7 = 𝑠8 = 𝜔 → 2. We have looked into the stability of the MRT model with two different boundary treatments. It is observed that the superior stability of the MRT model is dependent on the choice of boundary treatment as well as lattice resolution and collision frequency. Figure (3.7) compares the stability of SRT and MRT models coupled with half-way bounce back scheme and it is seen that MRT is only slightly more stable than the SRT model when coupled with bounce back boundary scheme. On the other hand, we were able to simulate cavity flow with Reynolds number as high as Re=10000 when MRT model is coupled with non-equilibrium boundary scheme. It is suggested that use of MRT model is mandatory when dealing with complex flow regimes such as pharyngeal bolus transport with high gradient regions and care has to be taken for the choice of boundary scheme. 32  3.3 Open boundary conditions As mentioned before the inlet and outlet boundary conditions in a swallowing simulation are not very well defined. The choice of inlet atmospheric pressure boundary condition can be a good approximation for the preliminary simulations since the bolus is at atmospheric pressure when the process of swallowing begins. The outlet condition can be approximated with an extended tube with unidirectional flow and zero normal velocity gradient at the outlet [29].   In this section we will look into the characteristics of Zou/He boundary conditions [36] when applied to the 2D pressure driven channel flow and will discuss the accuracy and stability of this boundary treatment when used with lattice Boltzmann BGK SRT model. Later we will investigate the extrapolation boundary scheme as the choice of outlet condition in a swallowing simulation. The extrapolation scheme is used when the flow is unidirectional with zero normal velocity gradient. 3.3.1 Zou/He boundary condition It is shown in section (2.1.2) that distribution functions are made of equilibrium and non-equilibrium parts. The equilibrium part is a function of macroscopic velocity and density values and the non-equilibrium part is a function of velocity, density, and gradients of velocity. Zou/He is based on the idea of bounce back for the non-equilibrium parts of the distribution functions that are normal to the boundary line. In other words, the fourth equation needed to solve equation (3.2.1) for the boundary node ‘Q’ in Figure (3.1) can be written as:  𝑓3𝑛𝑒𝑞(𝜌, 𝑢, ∇𝑢) = 𝑓1𝑛𝑒𝑞(𝜌, 𝑢, ∇𝑢)  (3.3.1) Separating 𝑓3 into equilibrium and non-equilibrium parts and using equation (3.2.1) we have:  𝑓3(𝑥, 𝑡) = 𝑓3𝑒𝑞(𝜌, 𝑢) + 𝑓3𝑛𝑒𝑞(𝜌, 𝑢, ∇𝑢) = 𝑓3𝑒𝑞(𝜌, 𝑢) + 𝑓1𝑛𝑒𝑞(𝜌, 𝑢, ∇𝑢) → 𝑓3(𝑥, 𝑡) = 𝑓3𝑒𝑞(𝜌, 𝑢) + 𝑓1(𝑥, 𝑡) − 𝑓1𝑒𝑞(𝜌, 𝑢) → 𝑓3(𝑥, 𝑡) = 𝑓1(𝑥, 𝑡) −23𝜌𝑢𝑥(𝑥, 𝑡) (3.3.2) With 𝑓3 solved in terms of the known parameters 𝑓1, 𝜌, and 𝑢𝑥 at the boundary node ‘Q’, the other two unknown distribution functions 𝑓6, 𝑓7 can be derived from equation (3.2.1):  {    𝑓3 = 𝑓1 −23𝜌𝑢𝑥𝑓6 = 𝑓8 −12(𝑓4 − 𝑓2) −12𝜌𝑢𝑦 −16𝜌𝑢𝑥𝑓7 = 𝑓5 +12(𝑓4 − 𝑓2) +12𝜌𝑢𝑦 −16𝜌𝑢𝑥  (3.3.3) The bounce back rule is also applied to the non-equilibrium parts of the two extra unknown distribution functions at the corner nodes. For example, 𝑓2𝑛𝑒𝑞(𝜌, 𝑢, ∇𝑢) = 𝑓4𝑛𝑒𝑞(𝜌, 𝑢, ∇𝑢) and 𝑓7𝑛𝑒𝑞(𝜌, 𝑢, ∇𝑢) = 𝑓5𝑛𝑒𝑞(𝜌, 𝑢, ∇𝑢) at the top right corner node in Figure (3.1). 33  We’ve analyzed the 2D steady state incompressible Poiseuille flow between two parallel plates with constant pressure inlet/outlet conditions. The inlet/outlet conditions are implemented using Zou/He boundary treatment and halfway bounce back is used for the fixed top and bottom wall conditions. Since the flow is pressure driven, accurate implementation of inlet/outlet conditions are of crucial importance. Theoretically, correct implementation of boundary conditions should result in second order accuracy, since both Zou/He and bounce back schemes are inherently second order accurate.  The flow geometry is a 2D channel with 𝐿𝑥 = 5 and 𝐿𝑦 = 1. The boundary conditions are as follow: {    𝑃𝑜𝑢𝑡 = 0 ,  𝑢𝑦 = 0  𝑎𝑡 𝑜𝑢𝑡𝑙𝑒𝑡𝑃𝑖𝑛 = 𝑃𝑜𝑢𝑡 +12𝑅𝑒𝐿𝑥𝜌𝐿𝑦3  ,  𝑢𝑦 = 0   𝑎𝑡 𝑖𝑛𝑙𝑒𝑡𝑢𝑥 = 𝑢𝑦 = 0   𝑜𝑛 𝑡𝑜𝑝 𝑎𝑛𝑑 𝑏𝑜𝑡𝑡𝑜𝑚 𝑤𝑎𝑙𝑙𝑠 The inlet pressure is chosen based on the approximated Reynolds number of the flow. Reynolds number of the flow is defined as 𝑅𝑒 = (𝑢𝑎𝑣𝑒𝐿𝑦)/𝜗 in which 𝑢𝑎𝑣𝑒 is calculated based on the analytical solution of the Navier-Stokes equations for 2D pressure driven Poiseuille flow. Since the pressure gradient in ‘x’ direction is a constant value, we have:  𝜕𝑝𝜕𝑥= 𝜇𝜕2𝑢𝜕𝑦2→ 𝑢 =12𝜇(𝜕𝑝𝜕𝑥) (𝑦2−(𝐿𝑦2)2)    𝑤ℎ𝑒𝑟𝑒  𝑦= 0 𝑎𝑡 𝑚𝑖𝑑 𝑝𝑙𝑎𝑛𝑒 → 𝑢𝑎𝑣𝑒 =1𝐿𝑦∫ 𝑢𝑑𝑦𝐿𝑦2−𝐿𝑦2= −13𝜇(𝜕𝑝𝜕𝑥) (𝐿𝑦2)2 → 𝑅𝑒 = −𝐿𝑦𝜌3𝜇2(𝜕𝑝𝜕𝑥) (𝐿𝑦2)2 (𝜕𝑝𝜕𝑥) =𝑝𝑜𝑢𝑡 − 𝑝𝑖𝑛𝐿𝑥→ 𝑃𝑖𝑛 = 𝑃𝑜𝑢𝑡 +12𝜇2𝑅𝑒𝐿𝑥𝜌𝐿𝑦3          (3.3.4) This analysis will allow us to investigate the stability range of the solver with respect to Reynolds number and to measure the critical Reynolds number as a function of collision frequency. Figure (3.12) provides the steady state lattice velocity and lattice density contours for the Poiseuille flow at Re=10. Figure (3.13) compares the fully developed velocity profile of the LBM solver with the profile obtained by the analytical solution. The good agreement between the two profiles sheds light on the validity of the implemented LBM solver. 34    (a) (b) Figure 3.12- Steady state velocity and pressure contours for 2D Poiseuille flow at Re=10 (a) non-dimensional velocity contour (b) lattice density contour  Figure 3.13- Fully developed 2D Poiseuille flow velocity profile for lattice resolutions Ny=50, 100, and 200 compared to analytical solution As can be seen in Figure (3.13), the numerical solution approaches the analytical solution as the lattice resolution is increased. The maximum value of the velocity ratio derived from analytical solution is 1.5 at the centerline of the channel which is also obtained by the LBM solver. 35   As mentioned before, lattice Boltzmann BGK model is inherently second order accurate in space, therefore we should obtain second order accurate results when boundary conditions with the same order of accuracy are implemented. To investigate this issue, L2 norm error analysis is performed on three different lattice resolutions Ny=50, 100, and 200 and at two different Reynolds numbers Re=10 and Re=40. The L2 norm is defined using the numerical and analytical velocity field as follow:  𝐿2 𝑛𝑜𝑟𝑚 𝑒𝑟𝑟𝑜𝑟 =√∑ ∑ ((𝑢𝑥(𝑖, 𝑗) − (𝑢𝑥𝑎(𝑖, 𝑗))2 + ((𝑢𝑦(𝑖, 𝑗) − (𝑢𝑦𝑎(𝑖, 𝑗))2)𝑁𝑦𝑖=0𝑁𝑥𝑖=0𝑁𝑥 × 𝑁𝑦 (3.3.5) In the above equation, 𝑢𝑥 and 𝑢𝑦 are numerical velocities while 𝑢𝑥𝑎 and 𝑢𝑦𝑎 denote the analytical velocities. Figure (3.) provides the L2 norm error as a function of lattice resolution. The slope of the line presents the order of accuracy for the solver and proves that second order accuracy is obtained using Zou/He and halfway bounce back boundary conditions.  Figure 3.14- L2 norm error analysis for Poiseuille flow  Stability of the Zou/He scheme has been measured for different choices of lattice resolution. Table (3.4) provides the minimum value of collision frequency for a stable solution as a function of lattice resolution and Reynolds number.  Table 3.4- Minimum stable collision frequency for 2D Poiseuille flow with respect to Reynolds number and lattice resolution Ny=10 Re 10 100 400 1000 minimum stable 𝝎 0.28 1.35 unstable  unstable Ny=20 Re 10 100 400 1000 minimum stable 𝝎 0.14 1.09 1.58 unstable Ny=40 Re 10 100 400 1000 minimum stable 𝝎 0.07 0.47 1.56 unstable  -5-4.5-4-3.5-3-2.5-2-1.5-1-0.500 1 2 3Log(L2)Log(Ny)Re=10Re=40Slope -236  As mentioned before, the relaxation time should satisfy the condition 𝜏 > 0.5 in order to have a physical flow with positive kinematic viscosity. Due to the reciprocal relationship between the relaxation time and the collision frequency, the same condition applies to the collision frequency in the form of 𝜔 < 2. Approaching the value of 𝜔 = 2 will cause instability in the high gradient regions of the flow. Figure (3.15) shows the convergence history of the Poiseuille flow with Re=1000 when the collision frequency is approaching the limit of 𝜔 → 2. The effect of low energy dissipation becomes apparent after around 10000 iterations and the residual starts raising for all variables. In order to overcome this issue, one should use finer resolutions to capture the high velocity gradients when simulating with 𝜔 → 2.  As can be seen in Table (3.44), the simulation is unstable at all choices of resolution for Re=1000. Therefore, we need to use higher lattice resolutions to decrease the Mach number and have a stable solution. The computational time needed for the mesh resolution 800*160 and collision frequency of 𝜔=1.8 (Mach=0.16) to converge on a Desktop PC with Intel(R) core™ i7 CPU at 3.30 GHz is 150125 seconds (about 42 hours). It is obvious that a complex pharyngeal flow with very high gradient regions and high Reynolds number (Re>1000) would require extremely fine lattice resolutions and as a result extremely large computational time.   Figure 3.15- Convergence history for Poiseuille flow with collision frequency set to 1.96 (𝜔 = 1.96) The compressibility effect of lattice Boltzmann method is shown in Figures (3.17) and (3.18).  The initial condition is no velocity everywhere and the pressure is constant and equal to outlet pressure 37  condition. As shown in Figure (3.1717.a), after 100 iterations the inlet pressure information is propagating into the first quarter of the domain while on the right boundary pressure remains unchanged until around iteration=400. The speed of convergence in a steady state problem and for a specific lattice resolution is dictated by the choice of collision frequency and the initial condition implementation. Increasing the single relaxation time 𝜏 or decreasing the collision frequency 𝜔 = 1/𝜏 will result in faster convergence. On the other hand, the Mach number of the flow will also increase as we raise the value of single relaxation time while keeping the lattice resolution and Reynolds number fixed, therefore higher compressibility error is observed. Figure (3.16) shows the effect of relaxation time parameter 𝜏 on the convergence rate and the accuracy of results for Poiseuille flow with a lattice resolution of 20*100 at Re=1.     (a) (b) Figure 3.16- Effect of relaxation parameter on the convergence rate and accuracy (a) history of convergence for average Mach number for 𝜏=0.6, 1, 2 (b) fully developed velocity profile for 𝜏=0.6, 1, 2  As can be seen in Figure (3.16.a), increasing the value of relaxation parameter will result in higher average flow Mach numbers at a fixed lattice resolution and Reynolds number. It is also clear that larger relaxation parameter will result in faster convergence rate as shown in the figure. The accuracy, on the other hand, is inversely related to the value of relaxation parameter due to the compressibility effect. Figure (3.16.b) compares the fully developed non dimensional velocity profile for different value of relaxation parameter 𝜏=0.6, 1, and 2. It is observed that compressibility error increases with larger relaxation values. The maximum non dimensional velocity on the channel centerline is equal to 𝑢/𝑢𝑖𝑛=1.54, 1.58, and 1.66 for 𝜏=0.6, 1, and 2 respectively. The error of the centerline velocity from the analytical solution with (𝑢𝑢𝑖𝑛)𝑎𝑛𝑎𝑙𝑦𝑡𝑖𝑐𝑎𝑙=1.5 is equal to 2.7%, 5.3%, and 10.7% for 𝜏=0.6, 1, and 2 respectively. In order to have a stable and accurate solution, an optimum value for the relaxation parameter has to be chosen based on the time and length characteristics of the flow.   38    (a) Iteration : 100 (b) Iteration : 400   (c) Iteration : 1000 (d) Iteration : 15000 Figure 3.17- Lattice density contour for the Poiseuille flow (a) iteration: 100 (b) iteration: 400 (c) iteration: 1000 (d) iteration: 15000    Figure (3.1717) proves that lattice Boltzmann BGK algorithm is able to achieve accurate results in the long run, but is very sensitive to the choice of initial conditions in the sense of time consumption. Inappropriate choice of initial condition can lead to very large computation cost in steady state lattice Boltzmann simulations, let alone the unsteady problems that depend on the time. Therefore precise implementation of initial velocity and pressure conditions become as important as implementation of boundary conditions in LB simulations.   39    (a) (b) Figure 3.18- Centerline lattice density profile for Poiseuille flow at different iterations  Figure (3.18) shows the lattice density profile along the channel centerline. The compressibility effect of lattice Boltzmann is obvious in Figure (3.1818.a). The pressure oscillations at every lattice node is very high at the first 1000 iterations and almost completely damped after 10000 iterations. In the Zou/He boundary scheme, the pressure waves in the interior of the computational domain interact with the inlet boundary in a specific manner. Recall that Zou/He is based on the idea of separation of distribution function into its equilibrium and non-equilibrium parts. Assuming that ‘A’ is a boundary lattice node with unknown distribution functions f1, f5, and f8 streaming from the solid into the bulk, implementing the Zou/He condition on the boundary node ‘A’ will result in:  𝑓1(𝑥𝐴) = 𝑓1𝑒𝑞(𝑥𝐴) + 𝑓1𝑛𝑒𝑞(𝑥𝐴) 𝑓1𝑛𝑒𝑞(𝑥𝐴) = 𝑓3𝑛𝑒𝑞(𝑥𝐴) 𝑓1𝑒𝑞(𝑥𝐴) = 𝐹𝑢𝑛𝑐𝑡𝑖𝑜𝑛 𝑜𝑓 [𝜌𝑖𝑛𝑙𝑒𝑡  , 𝑢𝑖𝑛𝑙𝑒𝑡] 𝑓1𝑛𝑒𝑞(𝑥𝐴) = 𝐹𝑢𝑛𝑐𝑡𝑖𝑜𝑛 𝑜𝑓 [𝜌𝑖𝑛𝑡𝑒𝑟𝑖𝑜𝑟 , 𝑢𝑖𝑛𝑡𝑒𝑟𝑖𝑜𝑟 , ∇𝑢𝑖𝑛𝑡𝑒𝑟𝑖𝑜𝑟]  (3.3.6)  Equation (3.3.6) suggests that interior pressure and velocity information streaming from the fluid region into the solid will reflect back into the bulk using the non-equilibrium distribution functions at boundary nodes. The equilibrium part, on the other hand, is only a function of macroscopic variables imposed on the boundary. This is clearly different from bounce back scheme when both equilibrium and non-equilibrium parts are reflected back into the flow domain. Generally, the value of non-equilibrium part is much smaller than the distribution part, therefore use of Zou/He scheme is preferred compared to the bounce back scheme in flows with singularity and high pressure gradients. This problem is further investigated on the flow over a cylinder benchmark in Figure (3.26).  Finding a solution for the compressibility effect observed in Figures (3.1717) and (3.1818) becomes crucial when dealing with time dependent flows. Another way of reducing the compressibility 40  effect is to use the incompressible lattice Boltzmann method with a slight modification to the weakly compressible method [25]. The equilibrium distribution function is modified as:  𝑓𝛼𝑒𝑞 = 𝜌𝑤𝛼 + ?̅?𝑤𝛼[3𝑐2𝑒𝛼 . 𝑢 +92𝑐4(𝑒𝛼 . 𝑢)2 −32𝑐2𝑢. 𝑢]  (3.3.7)  Where the constant ?̅? corresponds to the density of fluid. The lattice density is defined as 𝜌 = ?̅? +𝛿𝜌  with 𝛿𝜌 being the density fluctuations and the velocity is recovered using the value of ?̅? instead of 𝜌.  ?̅?𝑢 =∑𝑒𝑖𝑓𝑖8𝑖=0  (3.3.8) The boundary condition implementation is also slightly changed and equation (3.3.3) is modified to:  {    𝑓3 = 𝑓1 −23?̅?𝑢𝑥𝑓6 = 𝑓8 −12(𝑓4 − 𝑓2) −12?̅?𝑢𝑦 −16?̅?𝑢𝑥𝑓7 = 𝑓5 +12(𝑓4 − 𝑓2) +12?̅?𝑢𝑦 −16?̅?𝑢𝑥  (3.3.9)   The incompressible lattice Boltzmann scheme will be further investigated in section (3.5) when discussing the time dependent flow benchmark.  In this section we have analyzed the accuracy and stability of the Zou/He boundary treatment on a 2D pressure driven channel flow. The accuracy of the solver is measured using L-2 norm analysis in Figure (3.) and it’s shown that second order accuracy is obtained using lattice Boltzmann BGK model coupled with halfway bounce back and Zou/He boundary conditions. The stability of the solver is tested on three different lattice resolutions N=10, 20, and 40. The maximum Reynolds number that can be stably simulated, although with high compressibility error, using a value of 1.8 for collision frequency is measured to be Re=200,600, and 900 for the lattice resolutions N=10, 20, and 40 respectively. The value of 𝜔=1.8 is used for the collision frequency which will result in a flow Mach number of Ma=0.61 for the finest lattice resolution N=40 at Re=900. Large values for Mach number will result in high compressibility error, therefore collision frequency should be chosen with respect to the low Mach number rule. On the other hand, approaching the limit of 𝜔 → 2 for the collision frequency will result in instability around high gradient regions due to low energy dissipation in the model. A minimum lattice resolution of 160×800 is required to model Poiseuille flow at Re=1000 with Mach=0.16 (high Mach number) which will result in almost 42 hours of computational cost on an Intel(R) core™ i7 Desktop PC.  3.3.2 Extrapolation scheme As mentioned in section (3.1), the outlet boundary flow in a swallowing simulation can be approximated with a unidirectional flow with zero normal velocity gradient. In lattice Boltzmann 41  simulations the Neumann boundary condition is applied using second order differencing. Take the node ‘b’ in Figure (3.9) as an example, the second order backward differencing scheme is written for the zero normal velocity gradient as follow:  𝑢(𝑥𝑏 , 𝑡) = 2𝑢(𝑥𝑓 , 𝑡) − 𝑢(𝑥𝑓𝑓 , 𝑡)  (3.3.10) We have studied the extrapolation scheme introduced by [37] to apply the outflow boundary condition. Using equation (3.3.10) one can find the macroscopic velocity on the boundary node ‘b’. With the velocity known on the boundary, the problem turns into a Dirichlet boundary problem and will be solved in the same way as non-equilibrium extrapolation scheme described in section (3.2.2).  A 2D 2:1 contracting channel flow benchmark is used to investigate the validity of the implemented boundary condition. The uniform velocity condition is applied using Zou-He scheme for the inlet boundary and halfway bounce back is used to simulate no-slip condition at the walls. The described extrapolation scheme is implemented to simulate zero velocity gradient at the outlet. Figure (3.1919.a) shows the velocity contour at Reynolds Re=10. Reynolds is defined based on the inlet velocity and width of the outlet. A lattice resolution of 200×400 is used which means there are 200 lattice nodes on the inlet and 100 lattice nodes at the outlet. The relaxation time is chosen as 𝜏=0.6 which will result in a lattice Mach number of Ma=0.0058 during simulation. In order to compare the vortex dimensions with the lattice resolution, Figure (3.1919.b) focuses on the top vortex streamline configuration with the lattice grid shown on the background. The lattice spatial step is 𝛿𝑥 = 1/200 =0.005 which is in the same order as Mach number of the flow, therefore, both spatial discretization and compressibility error are the dominant sources of error during simulation.    (a) (b) Figure 3.19- Velocity contour and streamline for 2:1 contracting channel flow using extrapolation scheme at Re=10  Table (3.) compares the vortex center coordinate values obtained using LBM and Ansys Fluent finite volume solver. The error is defined as the difference between the values of the two methods over the characteristic length (width of the outlet). The data represented below suggests that center of vortex coordinate error is about 5 lattice nodes which can be due to spatial discretization error of both LBM 42  and FVM solvers and the compressibility error of lattice Boltzmann method. Higher choices of collision frequency will decrease the compressibility error while increasing the computational time as suggested by Figure (3.), therefore, care should be taken when choosing the collision frequency.  Table 3.5- Center of vortex coordinate comparison for 2:1 contraction flow between lattice Boltzmann and Ansys Fluent finite volume solver  LBM Fluent Error (%) Top Vortex (x) 0.945 0.970 5 Top Vortex (y) 0.947 0.968 4.2 Bottom Vortex (x) 0.945 0.970 5 Bottom Vortex (y) 0.053 0.032 4.2  We can see good agreement (Error<5%) between the results obtained by lattice Boltzmann and finite volume methods. The symmetry of the problem about the channel centerline is obvious in the results shown for top and bottom vortices in Table (3.). In order to further validate the implementation of the extrapolation boundary scheme, non-dimensional velocity profile along the channel centerline is shown in Figure (3.2020).   Figure 3.20- Centerline non-dimensional velocity profile for 2:1 contracting channel flow at Re=10 As shown in Figure (3.20), the zero velocity gradient condition is clearly satisfied at the outlet. The non-dimensional velocity is almost constant after x=1.5 and the flow becomes fully developed after this point. The value of the centerline non-dimensional velocity at the outlet is equal to 3 which is theoretically correct after the flow reaches its fully developed state and the velocity profile takes a parabolic shape. The maximum Reynolds number for the LB simulation to be stable, with all other parameters fixed, is measured to be Re=300. Changing the outlet boundary condition to constant pressure condition using Zou/He scheme will result in the same criteria for the stability and the 43  instability effects can be seen for Re>300. As a result, the extrapolation scheme can be considered as a more stable boundary scheme with respect to high Reynolds number.      3.4 Curved boundary conditions As mentioned before, the simplest lattice Boltzmann boundary treatment of the no-slip condition is the bounce back scheme. This method is able to model the no-slip condition on curved geometries by breaking them into stair-wise steps. It is proven that bounce back provides second order accuracy when the wall is placed exactly halfway between the two boundary nodes on the fluid and solid regions, but it’s not possible to achieve second order of accuracy when dealing with curved boundaries using this scheme. Also, breaking the geometry into stair wise steps would require very high lattice resolutions in order to achieve a stable solution. To overcome these issues, Filippova and Hanel [38] provided a second order boundary condition that preserves the integrity of geometry. Later, Mei et al, presented a modified condition based on Filippova and Hanel’s work to improve the numerical stability [39]. The main challenge with curved boundary is that they are arbitrarily located on the lattice domain and not exactly on the lattice nodes or halfway in between. Also, we cannot use the conservation law (see equation (3.2.1)) to derive the unknown macroscopic density or velocity because there are more unknown distribution functions at the boundary when dealing with concave curved walls compared to straight walls. As shown in Figure (3.21) for D2Q9 model, the boundary node ‘Q’ is surrounded by 5 nodes in the solid region shown as hollow circles and 3 nodes in the fluid region shown as solid black circles. Therefore, 5 distribution functions streaming from the solid to the node ‘Q’ are unknown and have to be reconstructed during simulations. The approach to model curved boundaries is based on the fitting of the distribution functions at the boundary nodes.   Figure 3.21- Concave curved wall condition with 5 unknown distribution functions streaming from solid into fluid   In section (3.4.1), the Filippova and Hanel boundary treatment (FH) for curved boundary condition as well as the improved version by Mei et al are studied and the stability and accuracy limit on the Reynolds number and lattice resolution is measured. We have tested the accuracy of the solver as a function of (lattice nodes/diameter of the cylinder) at different Reynolds numbers Re=1, 10, 40. The accuracy is measured using the drag coefficient value comparisons between LBM and Fluent results. Two different force evaluation schemes, Momentum exchange and Stress integration, are later introduced for drag coefficient calculations in lattice Boltzmann method. It is shown in Figure (3.27) that accuracy drops down significantly for low lattice resolutions due to the high compressibility error. 44  The effect of interior pressure waves and their interaction with the inlet boundary is also studied by monitoring the pressure fluctuations at point ‘M’ shown in Figure (3.23). This effect is shown and discussed in Figure (3.26). 3.4.1 Filippova-Hanel boundary dynamics The illustrative geometry for the method of direct fitting by Filippova and Hanel is shown in Figure (3.22). This figure is extracted from [23] with author’s permission.    Figure 3.22- The illustrative geometry for Filippova-Hanel scheme   The boundary wall in Figure (3.22) is denoted by the bold solid line and the solid black points represent the interstation of the boundary with lattice links. The streaming direction from fluid node ‘f’ to the solid node ‘b’ is denoted by 𝑒𝛼 and the opposite direction from solid to the fluid is denoted by 𝑒?̅?. The distance ratio ∆ is defined as:  ∆=|𝑥𝑓 − 𝑥𝑤||𝑥𝑓 − 𝑥𝑏|  ,   0 ≤ ∆≤ 1   (3.4.1) To construct the post collision distribution function streaming from the solid node ‘b’ to the fluid node ‘f’ based upon known information in the surrounding fluid nodes, Filippova and Hanel proposed the following linear interpolation: 45   𝑓?̅?(𝑥𝑏, 𝑡) = (1 − X)𝑓?̅?(𝑥𝑓 , 𝑡) + 𝑋. 𝑓?̅?(∗)(𝑥𝑏 , 𝑡) (3.4.2)  In equation (3.4.2), X is the fitting weight for the interpolation and 𝑓?̅?(∗)(𝑥𝑏 , 𝑡) is a fictitious distribution function, similar to the equilibrium function and defined as:  𝑓𝛼(∗)(𝑥𝑏 , 𝑡) = 𝑤𝛼𝜌(𝑥𝑓 , 𝑡)[1 +3𝑐2𝑒𝛼𝑢𝑏𝑓 +92𝑐4(𝑒𝛼 . 𝑢𝑓)2𝑢𝑓 . 𝑢𝑓] (3.4.3) Where 𝑢𝑓 = 𝑢(𝑥𝑓 , 𝑡) is the fluid velocity on the lattice boundary node ‘f’ and 𝑢𝑏𝑓 is to be chosen based on the distance ratio ∆. The fitting weight X is also chosen with the same criteria as 𝑢𝑏𝑓:  {𝑢𝑏𝑓 =(∆ − 1)𝑢𝑓∆+𝑢𝑤∆𝑋 =2∆ − 1𝜏     𝑤ℎ𝑒𝑛    ∆≥ 0.5 (3.4.4)  {𝑢𝑏𝑓 = 𝑢𝑓𝑋 =2∆ − 1𝜏 − 1     𝑤ℎ𝑒𝑛    ∆< 0.5 (3.4.5) Later, Mei et al used different nodes to measure 𝑓?̅?(𝑥𝑏, 𝑡) when ∆< 0.5. This would result in better numerical stability when the relaxation time 𝜏 is approaching the value of 1.  {𝑢𝑏𝑓 = 𝑢𝑓𝑓𝑋 =2∆ − 1𝜏 − 2     𝑤ℎ𝑒𝑛    ∆< 0.5 (3.4.6) Where 𝑢𝑓𝑓 = 𝑢(𝑥𝑓 + 𝑒?̅?𝛿𝑡, 𝑡) is the lattice velocity of the node ‘ff’. In this section, we have determined the stability and accuracy of the improved Filippova Hanel scheme (FH Mei) on a 2D flow over a bounded cylinder symmetrically placed in a channel. The geometry of the problem is shown in Figure (3.23).   Figure 3.23- Geometry of the flow over a bounded cylinder  Uniform velocity inlet boundary condition is applied using Zou/He scheme and zero axial velocity gradient condition is implemented for the outlet using the extrapolation scheme discussed in section (3.3.2) and no-slip condition is implemented using halfway bounce back algorithm for the top and bottom walls. 46  Figure (3.24) shows the horizontal and vertical velocity contours for flow at Re=10. Reynolds number is defined based on the inlet velocity as characteristic velocity and the cylinder diameter as characteristic length.    (a) (b) Figure 3.24- Velocity contours for flow over a bounded cylinder at Re=10 (a) x-velocity (b) y-velocity  To ensure the validity of the LB solver, velocity profile is obtained at specifically 10 lattice nodes after the last cylinder boundary node. The simulation is done with a node density of 20 lattice nodes per diameter of the cylinder. This means that a lattice resolution of 100×500 is chosen for the whole computational domain. Figure (3.25) compares the results obtained by LBM with the velocity profile obtained by Ansys Fluent finite volume solver at the same x-location of the flow domain x=1.2. The error of interpolation also exists in the comparison, since the lattice node coordinates does not coincide exactly with the cell coordinates of the mesh used by Fluent. The velocity profiles obtained by the two methods are in good agreement as suggested by the figure.        Figure 3.25- Velocity profile comparison between lattice Boltzmann and Ansys Fluent finite volume solver 47  Figure (3.26) investigates the pressure wave interactions between the interior and the inlet boundary implemented using Zou/He scheme as discussed in section (3.3.1). Point ‘M’ in Figure (3.23) is placed 10 lattice nodes away from the inlet boundary on the channel centerline. It is also 79 lattice nodes away from the front stagnation point of the cylinder at x=0.9. The initial velocity is set to uniform flow with inlet velocity everywhere in the flow domain, therefore initial discontinuity is observed on solid boundaries. The inlet width is increased to 2W so that the initial discontinuity on the front stagnation point of the cylinder reaches the point ‘M’ faster than the discontinuity on top and bottom walls. The variation of pressure at point ‘M’ is continuously monitored at all lattice time steps. Uniform flow initial condition is used at t=0, therefore discontinuity is introduced in the macroscopic variables pressure and velocity on the surface of the cylinder, as a result spurious pressure wave starts propagating into the domain. The earliest time that this wave associated with the initial discontinuity reaches to the point ‘M’ is dependent on the distance from the front stagnation point of the cylinder:  𝑡 =𝑙𝑎𝑡𝑡𝑖𝑐𝑒 𝑑𝑖𝑠𝑡𝑎𝑛𝑐𝑒𝑠𝑝𝑒𝑒𝑑 𝑜𝑓 𝑠𝑜𝑢𝑛𝑑= 79√3 ≈ 137 𝑙𝑎𝑡𝑡𝑖𝑐𝑒 𝑡𝑖𝑚𝑒 𝑠𝑡𝑒𝑝𝑠 (3.4.7) Figure (3.26) shows that first impulse of the pressure wave at point ‘M’ occurs at around t=137 lattice times steps, therefore the pressure wave propagation is actually following an acoustic behavior due to the compressibility effect. The inlet boundary partially reflects this pressure wave back into the interior domain using the non-equilibrium part of the distribution functions discussed in section (3.3.1). The second peak of the pressure variations represents the second interaction of this pressure wave at point ‘M’. The expected time for the second interaction using the forgoing analysis is:  𝑡 =𝑙𝑎𝑡𝑡𝑖𝑐𝑒 𝑑𝑖𝑠𝑡𝑎𝑛𝑐𝑒𝑠𝑝𝑒𝑒𝑑 𝑜𝑓 𝑠𝑜𝑢𝑛𝑑= (79 + 20)√3 ≈ 171 𝑙𝑎𝑡𝑡𝑖𝑐𝑒 𝑡𝑖𝑚𝑒 𝑠𝑡𝑒𝑝𝑠 (3.4.8) The second pressure peak in Figure (3.26) is seen at 𝑡 ≈ 171 lattice time steps which validates the acoustic propagation of the pressure waves. Also note the pressure (Lattice density) difference between the two peaks at times t=137 and t=171. The second peak is damped since the inlet boundary only partially reflects back the disturbance.    Figure 3.26- Pressure variations at point ‘M’ due to pressure wave interactions at the inlet boundary  48  In order to investigate the accuracy of the Filippova-Hanel scheme and the force evaluation scheme with respect to lattice resolution (lattice nodes/cylinder diameter), the drag coefficient is compared between the LBM solver and Ansys Fluent finite volume solver. Force evaluation on a curved boundary is definitely a quantity of interest during swallowing simulations. There are complex fluid-structure interactions happening between anatomical members and the food bolus during swallowing and an accurate force evaluation scheme is required to calculate these forces on the anatomical bodies. Two different force evaluation schemes, Momentum exchange and Stress Integration, are commonly used in lattice Boltzmann simulations. The stress Integration scheme uses great details about the surface geometry and extensive amount of extrapolation from neighboring nodes. It could be time consuming and challenging to program for 3D geometries. On the other hand, the momentum exchange scheme uses simpler calculations, but the accuracy of the results should be examined especially at high Reynolds flows [40]. The stress integration scheme is widely used by the NS equation based solvers, while the momentum exchange scheme is uniquely used by lattice Boltzmann method. Stress integration scheme, as the name mentions, integrates the normal and shear stress on the surface.     𝐹 = ∫[𝑝 + 𝜌𝜗(∇𝑢 + ∇𝑢𝑇)]𝑑𝐴 (3.4.9) In order to find an accurate value for the velocity gradient in equation (3.4.9), a fine lattice resolution should be used around the cylinder surface. In general, since velocity is not the primary variable in LBM calculations and its evaluation is based on the distribution functions, it is hard to find the gradient of velocity accurately on the curved boundary. The normal stress or pressure at boundary nodes can be easily calculated using the equation 𝑝 = 𝜌𝑐𝑠2 = 𝜌/3 associated with D2Q9 or D3Q15 lattice structures. Mei R [41] evaluated the shear stress using the non-equilibrium part of the particle distribution functions.  𝜏𝑖𝑗 = 𝜌𝜗(∇𝑖𝑢𝑗 + ∇𝑗𝑢𝑖) (3.4.10)  𝜏𝑖𝑗 = (1 −12𝜏)∑[𝑓𝛼(𝑥, 𝑡) − 𝑓𝛼𝑒𝑞(𝑥, 𝑡)]8𝛼=0× (𝑒𝛼𝑖𝑒𝛼𝑗 −12𝑒𝛼 . 𝑒𝛼𝛿𝑖𝑗) (3.4.11) In equation (3.4.11),[𝑓𝛼(𝑥, 𝑡) − 𝑓𝛼𝑒𝑞(𝑥, 𝑡)] is the non-equilibrium part of the distribution function in 𝛼th lattice direction and 𝑒𝛼𝑖 and 𝑒𝛼𝑗 are the ith and jth components of the discrete lattice velocity 𝑒𝛼 respectively.  The stress vector has to be computed on the curved body surface through an extrapolation procedure, since the Cartesian grid is used for lattice Boltzmann simulations. This will result in further loss of accuracy when the resolution is not sufficient to capture the shear layer effects near the boundary.  The momentum exchange method effectively breaks the geometry surface into stair wise steps and computes the force on each boundary lattice link placed halfway between the fluid and the solid nodes. Take node ‘b’ in the Figure (3.22) as an example, the force results from the momentum exchange rate [𝑒?̅?𝑓?̅?(𝑥𝑏 , 𝑡) − 𝑒𝛼𝑓𝛼(𝑥𝑓 , 𝑡)]𝛿𝑥/𝛿𝑡 between the two opposing directions of the neighboring lattices. The total force acting on a solid surface by fluid is calculated as: 49   𝐹 = ∑ ∑𝑒𝛼[𝑓𝛼(𝑥𝑏 , 𝑡) + 𝑓?̅?(𝑥𝑓 , 𝑡)] × [1 − 𝑤(𝑥𝑏 + 𝑒?̅?𝛿𝑡)]𝛿𝑥/𝛿𝑡8𝛼=1𝑎𝑙𝑙 𝑥𝑏 (3.4.12) Where 𝑤(𝑥𝑏 + 𝑒?̅?𝛿𝑡) is an indicator which is equal to zero at the fluid nodes and equal to 1 at solid nodes. The inner summation computes the momentum exchange of node ‘b’ with all it’s possible neighboring fluid nodes. The outer summation adds up the forces calculated for all of the boundary nodes like ‘b’. It is relatively easy to implement the momentum exchange scheme both for 2D and 3D simulations, however its accuracy has to be measured for curved boundaries. The performance of the two force evaluation schemes are compared by [40]. In summary, the momentum exchange method is simple to implement and can give reliable results when a sufficient lattice resolution is used for the simulation. The stress integration scheme generates similar results when sufficient lattice resolution is used, but there is a much larger uncertainty in the results when the resolution is limited. The computational effort of stress integration scheme in implementing the extrapolation and integration on the surface is extremely high when compared with the method of momentum exchange.  We have considered the accuracy of momentum exchange method at different Reynolds numbers Re=1, 10, 40 as a function of (lattice nodes/cylinder diameter). The drag coefficient is defined as:  𝐶𝐷 =𝐹𝑥12𝜌0𝑈02𝐴 (3.4.13) Where Fx is the net axial force on the cylinder and A is the effective area. Figure (3.27) shows the drag coefficient value on the vertical axis as a function of lattice resolution (lattice nodes/cylinder diameter) for a constant diameter of D=W/5. Different lattice resolutions are used ranging from Ny=30 to Ny=120 and the relaxation time is set to 𝜏=0.6. In Figure (3.27), the solid delta symbols represent the drag coefficient obtained by LB momentum exchange force evaluation scheme and the black line represents the value of drag coefficient obtained by Ansys Fluent finite volume solver with double precision accuracy using a fine structured mesh (≈25000 elements). The Mach number of the flow with Re=1 at inlet with a lattice resolution of Ny=30 (6 lattice nodes/Diameter) is about Ma=0.01 which is much smaller than the spatial step ∆𝑥=1/6≈0.17. As we increase the lattice resolution with a constant relaxation time (𝜏=0.6), the Mach number decreases and reaches the value of Ma=0.002 with a lattice resolution of Ny=120 (∆𝑥=1/24≈0.04), therefore the error of spatial resolution O(∆𝑥2) is dominant during all simulations at Re=1. The error of spatial discretization remains the same order for all of the three different Reynolds numbers Re=1, 10, 40. On the other hand, the inlet Mach number is directly proportional to Reynolds number and the value of Ma will change from Ma=0.1 to Mach=0.02 at Re=10 which means that error is dominant by both spatial resolution O(∆𝑥2) and compressibility effect O(𝑀𝑎2). The error of Mach number will further increase at Re=40 and will become the dominant error during the simulations. We can see from the figures that the minimum lattice resolution required to reach a reliable value for drag coefficient increases noticeably at Re=40 as the compressibility error becomes the dominant source of error. In order to have a drag coefficient with relative error (𝐶𝐷𝐿𝐵𝑀−𝐶𝐷𝐹𝑙𝑢𝑒𝑛𝑡𝐶𝐷𝐹𝑙𝑢𝑒𝑛𝑡 ) < 1% for the flow at Re=40, we need to have at 50  least a resolution of 24 lattice nodes/cylinder diameter which leads to (64592 s) 18 hours of computational time on an Intel(R) core™ i7 desktop PC. It is worth mentioning that the error of extrapolation for momentum exchange approach is also higher at low lattice resolutions and will further degrade the accuracy of the results obtained for drag coefficient.   (a) (b)  (c)  Figure 3.27- Drag coefficient as a function of lattice nodes/cylinder diameter (a) Re=1 (b) Re=10 (c) Re=40 In this section we have looked into the implementation and accuracy of improved Filippova Hanel curved boundary scheme as well as the accuracy of the momentum exchange force evaluation scheme. The validity of the lattice Boltzmann solver coupled with the FH boundary scheme is tested in Figures (3.25) and (3.27). The compressibility effect of lattice Boltzmann method and the acoustic nature of the pressure wave propagation into the computational domain is shown in Figure (3.26). This figure suggests that consistent initial condition implementation is of crucial importance in a transient swallowing simulation using lattice Boltzmann method. Initial discontinuity in velocity and pressure is 51  the source of acoustic wave propagations in the domain that can degrade the stability of the solver. An effective approach for consistent initial condition implementation will be introduced in section (3.5) to address this issue. We have also looked into the implementation and accuracy of the momentum exchange force evaluation scheme for lattice Boltzmann simulations, since the force on anatomical bodies such as the epiglottis is definitely a quantity of interest during swallowing simulations. It is observed in Figure (3.27) that in order to reach a reliable value (error < 2%) for drag coefficient of a bounded cylinder in a 2D channel, we need to have at least 20 lattice nodes/Diameter of the cylinder at Re=40. In order to have a lattice Mach number of Ma=0.05 when the Reynolds number is as high as Re=1000, i.e. low compressibility error, a lattice resolution of 500*2500 is required with relaxation rate set to 𝜏=0.51.  It is mandatory to use the multi relaxation time model when the relaxation rate is approaching the inviscid flow limit of 𝜏 → 0.5. A multi relaxation time lattice Boltzmann simulation with a lattice resolution of 500*2500 approximately requires 1 week of computational time on a desktop PC with coreTM i7 CPU at 3.30 GHz. It seems unfeasible to model a swallowing simulation at Re>1000 with complex curved boundaries using lattice Boltzmann method considering the extremely high computational cost.   3.5 Transient simulations  As seen in section (3.3.1) and (3.4.1), the initial condition implementation is of crucial importance in LB simulations. The spurious pressure waves due to initial discontinuity can cause instability and low accuracy during time dependent simulations. In this section, we will focus on the important issue of consistent initial condition implementation for distribution functions in lattice Boltzmann MRT simulations. Using the equilibrium distribution function 𝑓𝛼𝑒𝑞(𝜌, 𝑗) is a common choice to initialize the distribution function 𝑓𝛼 at all lattice nodes. In that case both initial velocity and pressure field u0, p0 are required in order to find the equilibrium function 𝑓𝛼𝑒𝑞(𝜌, 𝑗), but the initial pressure condition is not given in many problems. This problem is usually taken care of by solving the Poisson equation to obtain the pressure field and the density field via the equation of state. This approach of initialization can be laborious and not easy to implement. An iterative initialization procedure is proposed by Mei [13] which is based on the idea that density (𝜌) is the only conserved variable and the momentum (j) is relaxed to the state imposed by the initial conditions for velocity. The implementation of the proposed scheme is relatively easy and very consistent with multi-relaxation-time model implementation and it can be shown using a Chapman-Enskog analysis that the proposed initialization procedure will solve the Poisson equation for the pressure. In this section, the implementation of the iterative initialization procedure proposed by Mei is explained and the compressibility effect of the lattice Boltzmann model is investigated on the time dependent channel flow simulation with an oscillating inlet pressure condition in time.   The algorithm proposed by Mei uses an iterative process to initialize the distribution functions. Recall equation (2.1.15) for the evolution equation of distribution functions in lattice Boltzmann MRT model. As described in section (2.1.2), the vector of hydrodynamic moments 𝑅 =(𝜌, 𝑒, 𝜀, 𝑗𝑥 , 𝑞𝑥 , 𝑗𝑦, 𝑞𝑦, 𝑝𝑥𝑥 , 𝑝𝑥𝑦)𝑇is linearly related to the vector of distribution functions 𝐹 =(𝑓0, 𝑓1, 𝑓2, 𝑓3, 𝑓4, 𝑓5, 𝑓6, 𝑓7, 𝑓8)𝑇 by the linear transformation matrix M in equation (2.1.16). Density 𝜌 and momentum 𝑗𝑥, 𝑗𝑦 are the only conserved moments in the system. The other moments are non-conserved moments that are functions of the conserved moments in the equilibrium state (equations (2.1.17) to (2.1.22)). The relaxation matrix ?̂? = 𝑑𝑖𝑎𝑔(𝑠𝜌, 𝑠𝑒 , 𝑠𝜀 , 𝑠𝑗, 𝑠𝑞 , 𝑠𝑗, 𝑠𝑞 , 𝑠𝜗, 𝑠𝜗) is a diagonal 9×9 52  matrix including the relaxation rates for different hydrodynamic moments. Since the density and momentum are conserved variables, the corresponding relaxation rates 𝑠𝜌 and 𝑠𝑗 have no effect on the system and can take any value, however, the value of 𝑠𝑗 would be effective if an external force is applied to the system in terms of a body force in the evolution equation. The shear viscosity and the bulk viscosity of the flow in lattice simulations are obtained using the relaxation rates 𝑠𝜗 and 𝑠𝜀 [42].  We will use the incompressible lattice Boltzmann method to reduce the compressibility effect in transient simulations as described in section (3.3.1). The mean density is assumed to be ?̅? = 1 for simplicity and the density fluctuations are represented by 𝛿𝜌, thus the conserved moments will be written as:  𝛿𝜌𝑒𝑞 = 𝛿𝜌 =∑𝑓𝛼 (3.5.1)  𝑗𝑒𝑞 = 𝑗 = ?̅?𝑢 =∑𝑒𝛼𝑓𝛼 (3.5.2)  In an athermal MRT lattice Boltzmann simulation (internal energy is not a conserved quantity) the non-conserved moments at equilibrium state are given by equations (2.1.17) to (2.1.22). We can retrieve the incompressible Navier-Stokes equations using equations (2.1.17) to (2.1.22) in the limit of small Mach numbers through a Chapman-Enskog procedure. We can reduce the MRT equation (equation (2.1.15)) to single relaxation BGK equation (equation (2.1.5)) by setting all of the relaxation rates in Matrix ?̂? to a single value of 1/𝜏. The speed of sound in the model is 𝑐𝑠 = 1/√3 and the shear viscosity 𝜗 and the bulk viscosity 𝜉 are determined by:  𝜗 =13(1𝑠𝜗−12) (3.5.3)  𝜉 =16(1𝑠𝑒−12) (3.5.4)  As mentioned before, the initial pressure field p0 is not always given in flow simulations and we have to solve Poisson equation to find p0. A common initialization technique in lattice Boltzmann simulations is to set the initial lattice density 𝜌0 (related to the initial pressure field p0 through the equation of state) to a constant value while the initial velocity fields u0 and v0 are given by flow initial conditions. The equilibrium moment equations can then be calculated using the initial values of the conserved moments 𝜌0 and 𝑗0. Finally, the equilibrium distribution functions are determined by the linear mapping 𝑓𝑒𝑞 = 𝑀−1.𝑚𝑒𝑞 and used as the initial value for distribution functions 𝑓. This method of initialization is inaccurate, since the initialization error often persists throughout the entire simulation in the form of compressibility effect.    In a new initialization approach by Mei, the original LBE scheme with three conserved moments 𝜌, 𝑗𝑥, and 𝑗𝑦 is replaced by  the corresponding generalized diffusion–advection LBE with one conserved moment 𝜌 [43]. In this case, the flow momentum 𝑗 is not a conserved quantity anymore and relaxes towards the state imposed by the initial condition u0 and v0. The initial density fluctuations 𝛿𝜌0 is initialized to zero for simplicity, therefore the initial momentum can be written as 𝑗0 = ?̅?𝑢0. The relaxation is done through an iterative process until a steady state for distribution functions is obtained. 53  The flow momentum 𝑗 is iteratively relaxed towards the equilibrium state of flow momentum 𝑗𝑒𝑞 which is set to be equal to the initial flow momentum 𝑗0 = ?̅?𝑢0.  𝑗∗ = 𝑗 − 𝑠𝑗(𝑗 − 𝑗𝑒𝑞) (3.5.5)  After initialization of density 𝛿𝜌0 = 0 and momentum 𝑗0 = ?̅?𝑢0, the equilibrium moments are constructed using equations (2.1.17) to (2.1.22). After construction of equilibrium moments, the iterative process starts with the following steps: 1. The non-conserved moments (all moments except density) are relaxed towards their equilibrium state using equation (2.1.15) 2. Streaming of distribution functions 3. Update the density fluctuation using equation (3.5.1) 4. Repeat the process until a steady value is reached for density (∑ |𝛿𝜌(𝑥𝑖, 𝑡 + 1) − 𝛿𝜌(𝑥𝑖 , 𝑡)|𝑖 <𝑡𝑜𝑙)    Density is the only conserved variable during the proposed iterative process and it satisfies the diffusion equation which can be derived using the Chapman-Enskog analysis up to second order [43]. The diffusion equation is written as:  𝜕𝑡(𝑐𝑠2𝜌) = 𝜒∇2(𝑐𝑠2𝜌) + 𝜒∇∇: 𝑗0𝑗0 − 𝑐𝑠2∇. 𝑗0 (3.5.6) In equation (3.5.6), pressure 𝑝 is replaced by 𝑐𝑠2𝜌 from the equation of state. The diffusion coefficient 𝜒 is derived as:  𝜒 = 𝑐𝑠2(1𝑠𝑗−12) (3.5.7) After the iterative procedure reaches steady state, the flow momentum 𝑗 approaches the prescribed initial velocity field given by 𝑗0 = ?̅?𝑢0 and the initial value obtained for density satisfies the Poisson equation if ∇. 𝑢0 = 0.  ∇2(𝑐𝑠2𝜌0) = −∇. (𝑢0. ∇𝑢0) (3.5.8) Equation (3.5.8) is satisfied when the iterative process reaches steady condition and the value obtained for initial pressure (equivalently initial density) is independent of the choice of relaxation rates for MRT model. The athermal LBE and the diffusion-advection LBE follow the same evolution equation (equation (2.1.15)) for initialization process with the exception that flow momentum 𝑗 is identical to 𝑗𝑒𝑞 for athermal LBE while it relaxes towards 𝑗𝑒𝑞 for diffusion-advection initialization procedure. As a result, the distribution functions are the same for the two schemes when the steady state condition is achieved, i.e., when 𝑗 = 𝑗𝑒𝑞.   We have tested the accuracy of the incompressible lattice Boltzmann method on the transient 2D channel flow with dimensions Lx=2 and Ly=1. The boundary conditions for the inlet/outlet and the solid walls are defined as follow: {𝑃𝑜𝑢𝑡 = 0 ,  𝑢𝑦 = 0  𝑎𝑡 𝑜𝑢𝑡𝑙𝑒𝑡𝑃𝑖𝑛 = 𝑃𝑜𝑢𝑡 + 𝐴𝐿𝑥cos (𝜔𝑡) ,  𝑢𝑦 = 0   𝑎𝑡 𝑖𝑛𝑙𝑒𝑡𝑢𝑥 = 𝑢𝑦 = 0   𝑜𝑛 𝑡𝑜𝑝 𝑎𝑛𝑑 𝑏𝑜𝑡𝑡𝑜𝑚 𝑤𝑎𝑙𝑙𝑠 54  Where 𝜔 is the oscillation frequency and 𝐴 is the amplitude. The non-dimensional oscillation parameter (Womersley number) 𝛼 is defined as:  𝛼 =𝐿𝑥2√𝜔𝜗 (3.5.9)  In section (3.3.1) we displayed the compressibility effect of lattice Boltzmann method in Figures (3.17) and (3.18). The solution for this problem can be achieved by letting the simulation run for a long time until steady state is reached. However, for time dependent problems we need to use the incompressible lattice Boltzmann method described in section (3.3.1), otherwise the compressibility effect will keep the adverse variable oscillations persisting in the flow domain for a while causing instability and inaccuracy at later time steps. The main idea behind the incompressible lattice Boltzmann model is to eliminate the terms of 𝑂(𝑀𝑎2) that result in density fluctuations in the weakly compressible lattice Boltzmann model [44].   Figures (3.28) and (3.29) show the non-dimensional velocity profile at x=1 along the channel obtained using incompressible multi relaxation time lattice Boltzmann method compared with results obtained by Ansys Fluent finite volume solver for two different oscillation frequencies 𝜔 = 2𝜋 and 𝜔 = 20𝜋 and with initial velocity set to zero everywhere in the flow domain.   Figure 3.28- Non-dimensional velocity profile comparison between incompressible lattice Boltzmann and Ansys Fluent finite volume solver at Re=10 and 𝜔 = 2𝜋 55   Figure (3.28) shows the velocity profile comparison at Re=10 with the oscillation frequency 𝜔 =2𝜋 which results in 𝛼 = 3.96. The spatial resolution NY=50 is used for both lattice Boltzmann and finite volume solver and the lattice time step is chosen based on the small Mach number criteria (𝛿𝑡 ≈𝛿𝑥2). The characteristic velocity is defined using the analytical solution for the steady state poiseuille flow with constant pressure gradient (𝑃𝑖𝑛 − 𝑃𝑜𝑢𝑡)/𝐿𝑥 = 𝐴.        𝑈0 = 𝑈𝑎𝑣𝑒𝑟𝑎𝑔𝑒 = −13𝜇(𝜕𝑝𝜕𝑥) (𝐿𝑦2)2=𝐴3𝜇(𝐿𝑦2)2 (3.5.10)  Figure (3.29) shows the same velocity profile comparison at Re=10 for flow with inlet pressure oscillation frequency 𝜔 = 20𝜋. The same spatial resolution is used for both lattice Boltzmann and finite volume simulations, however, smaller lattice time step has to be chosen due to the change in physical time scale (𝑇 = 1/𝜔) of the flow. The velocity profile comparisons in Figures (3.28) and (3.29) suggest that incompressible lattice Boltzmann method is good at catching the transient flow if the Mach number limit 𝑀𝑎 → 0 is satisfied. In order to satisfy the low Mach number limit one can use very high lattice resolution which results in extremely high computational time or very low relaxation rate 𝜏 → 0.5 which dictates the use of multi relaxation time model.  Figure 3.29- Non-dimensional velocity profile comparison between incompressible lattice Boltzmann and Ansys Fluent finite volume solver at Re=10 and 𝜔 = 20𝜋 56  In order to have a meaningful comparison between the computational time of the lattice Boltzmann method and the finite volume solver, one has to make sure that both methods are generating results with the same order of numerical error. In the following analysis, we’ve chosen the grid resolution to be fixed Ny=100 during all LBM and FVM simulations, so the spatial discretization error is constant 𝑂(∆𝑥2) = 1𝑒−4. As mentioned before, the compressibility error of the lattice Boltzmann method can become the dominant source of error during simulations with high Reynolds number. Table (3.6) provides the value of relaxation rate 𝑠𝜗 required to keep the compressibility error 𝑂(𝑀𝑎2) in the same order of spatial discretization error 𝑂(∆𝑥2) as we increase the Reynolds number, i.e. 𝑀𝑎 ≈ ∆𝑥 ≈0.01. The oscillation frequency is set to 𝜔 = 2𝜋 and the physical time step of LBM simulation is determined based on the choice of relaxation rate. The fourth column provides the physical time step of the lattice Boltzmann simulations while keeping the lattice Mach number constant 𝑀𝑎 ≈ ∆𝑥 ≈ 0.01. Since multi relaxation time model is used for simulations, the choice of relaxation rate can accept values very close to the inviscid flow limit 𝑠𝜗 → 0.5.  Table 3.6- Choice of relaxation rate 𝑠𝜗 as a function of Reynolds number   𝒔𝝑 Mach ∆𝒕 (s) Re=10 0.70 0.0115 6.67e-04 Re=100 0.52 0.0115 6.67e-05 Re=500 0.504 0.0115 1.33e-05 Re=2000 0.501 0.0115 3.33e-06  As suggested by the table above, we need to decrease the value of relaxation rate  noticeably in order to keep the compressibility error in the same order as spatial discretization error as we increase the flow Reynolds number. It’s worth mentioning that the finite volume solver can accept much larger time steps, especially when using the implicit time advancement scheme. Since the compressibility error is specific to lattice Boltzmann method, the numerical accuracy of the finite volume solver does not change noticeably as we increase the flow Reynolds number. A time step of ∆𝑡=0.001 will result in a stable solution for the finite volume implicit solver at Re=2000. Therefore, the lattice Boltzmann method requires almost 300 times more time steps  to simulate the flow with Re=2000 at a specific physical time. As mentioned before, the range of Reynolds number for swallowing simulations can change from low Reynolds numbers Re<100 for thick materials, e.g. honey, to high Reynolds number Re>10000 for thin materials such as air; therefore, it becomes crucial to look after the compressibility error during swallowing simulations. As estimated in section 1.1, the Reynolds number of the flow can reach values as high as Re=2000 when dealing with single phase continuous flow of water in the pharynx; therefore, the CPU time required for the lattice Boltzmann MRT scheme to simulate the continuous flow of water inside the pharynx will be much larger compared to finite volume implicit solver due to the restrictions that compressibility effect will put on the choice of relaxation rate, i.e. time step.   57  Chapter 4. Conclusion  Lattice Boltzmann method has been widely used during the past decade in the land of computational fluid dynamics due to its simplicity and intrinsic parallelism. In addition, the fixed grid nature of this method eliminates the burden of remeshing at each time step for flow problems with moving boundaries. Aside from these advantages to the method, there are many shortcomings that has to be investigated with respect to the sensitivity of the method to the choice of lattice resolution and relaxation rate, especially for high Reynolds number flows. In this work we have focused mainly on three topics of compressibility effect, stability range for SRT and MRT models, and the important issue of boundary condition implementation. These appear to be three of the main issues that has to be investigated before applying the method to complicated 3D swallowing simulations. Generally speaking, increasing the lattice resolution or approaching the value of inviscid flow limit for the relaxation rate 𝜏 → 0.5 are the two ways of decreasing the lattice Mach number and observing less compressibility error during simulations. Using multigrid lattice Boltzmann method can decrease the computational time significantly, however, it is computationally very challenging to apply the multigrid scheme to swallowing simulations since the high gradient region of the flow domain moves with the bolus.      Looking at the lid-driven cavity flow benchmark in section (3.2.2), the minimum computational time required for the lattice Boltzmann SRT model to achieve a stable solution for flow at Re=3200 on a 800*800 lattice grid and with relaxation parameter 𝜏 = 0.6 (𝑀𝑎 = 0.23) is about 70 hours on a desktop computer with coreTM i7 CPU at 3.30 GHz. Note that flow Mach number 𝑀𝑎 = 0.23 will result in very high compressibility error and instability during time dependent simulations therefore higher lattice resolution, i.e. higher computational cost, is required for time dependent flow simulations using single relaxation time model. Take the continuous swallowing of water as a simple example of single phase swallowing problem. Assume that pharynx can be modeled as a tube with a diameter of D=2 cm and a typical length of L=11 cm [29]. Usually the flow of water reaches the upper esophageal sphincter within a time frame of 1 second once the pharyngeal phase of swallowing begins, therefore the average velocity of the flow can be estimated as V=11 cm/s. As a result the Reynolds number of continuous swallowing of water can reach values as high as Re=2200, therefore lattice Boltzmann single relaxation time model requires several days on a desktop PC to simulate this very simplistic 2D model of single phase continuous swallowing of water in the pharynx modeled as a perfect circular tube. Figure (3.5) shows that CPU time is related to lattice size (N) by a power of three, 𝐶𝑃𝑈𝑡𝑖𝑚𝑒 ∝ 𝑁3, which results in extremely huge computational time for high Reynolds flow simulations using SRT lattice Boltzmann model. To overcome this issue, the multi relaxation time model is introduced that allows the relaxation rate 𝜏 to accept values very close to the inviscid limit 𝜏 → 0.5 without compromising the solution stability. Using the MRT model coupled with the non-equilibrium extrapolation boundary scheme allows the relaxation rate to accept the value of 𝜏 = 0.5015 that results in a Mach number of 𝑀𝑎 =0.05 on a 200*200 lattice grid at Re=10000. Therefore, use of multi relaxation time model can significantly decrease the computational cost in preliminary swallowing simulations at high Reynolds. However, the stability range of multi relaxation time model has to be investigated when applied to problems with arbitrary moving boundary conditions such as pharyngeal swallowing.  The choice of inlet and outlet boundary conditions for lattice Boltzmann simulations can affect the stability of the simulations, especially when dealing with time dependent flows. The 58  information propagation in lattice Boltzmann scheme follows an acoustic nature, therefore the characteristic velocity of the problem should be much smaller than the speed of sound in the model. Larger density oscillations are observed in simulations with higher lattice Mach number due to the compressibility effect of the method. These adverse oscillations are initiated due to the initial discontinuity of variables on the boundaries and are reflected back into the domain as they interact with the flow boundaries at later time steps. The bounce back boundary treatment reflects the equilibrium and non-equilibrium parts of the distribution function back into the flow domain at the boundary nodes, whereas, boundary treatments such as Zou/He reflects back only the non-equilibrium part of the distribution function and the equilibrium distribution function is determined based on the imposed macroscopic variables on the boundaries. In the latter, the adverse oscillations are only partially reflected back into the flow domain by the non-equilibrium part of the distribution function. Initial condition implementation is also of crucial importance in time dependent lattice Boltzmann simulations. In order to eliminate the density oscillations we have used the incompressible lattice Boltzmann method to model the transient channel flow in section (3.5). The incompressible Navier-Stokes equations can be recovered on the low limits of Mach number 𝑀𝑎 → 0 using the incompressible lattice Boltzmann evolution equation through a Chapman-Enskog analysis.   The issue of force evaluation on anatomical bodies is of noticeable interest in swallowing simulations. In section (3.4.1), the accuracy of the momentum exchange force evaluation scheme is measured using the drag coefficient of a cylinder placed in a 2D channel. The measurements in Figure (3.27) show that the results obtained by the momentum exchange scheme for drag coefficient is very sensitive to the lattice resolution (lattice nodes/cylinder diameter). A minimum lattice resolution of 20 nodes/diameter is required to obtain reliable drag coefficient results for flow at Re=40 where the compressibility error is the dominant source of error in the simulation. Use of multigrid lattice Boltzmann method around the anatomical bodies with sharp edges will increase the accuracy of the momentum exchange scheme and requires further investigation. The force evaluation on the epiglottis is an example that has attracted scholars’ attention during the past few years. Assume that the tip of the epiglottis can be modeled as a sphere with a diameter equal to the thickness at the midline of the epiglottis. According to the literature on epiglottis dimensions, the average thickness of the tip of the epiglottis on the mid sagittal plane, i.e. diameter of the cylinder, is D=1.68 mm for a group of Japanese adult men and D=1.45 mm for Japanese adult women [47]. The average diameter of the pharynx during exhaling has measured to be 12.9 mm for men and 11.6 mm for women respectively. Assuming that flow Reynolds number is as low as Re=40 during swallowing, a minimum lattice resolution of 20 nodes/diameter for the tip of the epiglottis will result in a minimum 2D lattice resolution of 𝑁 = 154 ×1314 for the entire domain of an adult man’s pharynx with a length of L=11cm. The computational time required to simulate this 2D swallowing simulation using lattice Boltzmann MRT is estimated to be around 12 hours on a desktop PC with core™ i7 CPU at 3.30 GHz. Switching from 2D to 3D simulations will increase the computational time significantly and it is obvious that Reynolds number can be noticeably higher than Re=40, therefore finer lattice resolution is required which results in even more computational cost. In other words, a 3D lattice Boltzmann MRT simulation of swallowing including force evaluation on the epiglottis with Reynolds number as high as Re>1000 will require weeks of computational time.   The future steps for this work are aimed at deceasing the computational time of the lattice Boltzmann method. The first step is to look into the implicit formulation of the lattice Boltzmann method. This will allow to use much larger time steps for high Reynolds number simulations while sacrificing the simplicity and the intrinsic parallelism of the explicit formulation. Jing Liu introduces 59  first order, second order and fourth order implicit formulations of the Boltzmann equation in his work [48]. The next step is to look into multi grid lattice Boltzmann algorithms for the issue of force evaluation around boundaries with sharp edges such as epiglottis. Chen Peng [25] provides an overview on grid refinement algorithms for explicit lattice Boltzmann method and Jing Liu introduces an extension of the multi block method for the implicit formulation. The stability and accuracy of this scheme has to be investigated when coupled with moving boundary condition. The moving boundary condition is the next topic that has to be investigated. The simplest moving boundary condition is to use the bounce back scheme with interpolations at each time step [49]. However, the bounce back scheme reflects both the equilibrium and non-equilibrium parts of the distribution function back into the flow domain, therefore the adverse pressure oscillation due to initial discontinuity remains longer in the flow domain compared to more complicated boundary schemes. The stability of the bounce-back scheme has to be investigated when dealing with high Reynolds flows with sharp boundaries such as swallowing simulations. Also, further investigation is needed for lattice Boltzmann implicit immersed boundary methods. The immersed boundary method is an effective method for fluid–structure interaction problems. Jian Hao [50] proposed a lattice Boltzmann based implicit immersed boundary method which solves the fully nonlinear algebraic system resulting from discretization by an Inexact Newton–Krylov method. The stability of the immersed boundary scheme coupled with implicit lattice Boltzmann scheme requires further study for swallowing modeling problems.   To summarize, it is observed in this work that lattice Boltzmann method suffers from several shortcomings in handling high Reynolds number flows, sharp curved boundaries, and flows with high Mach number. The compressibility error which is inherent to the method of lattice Boltzmann increases with Reynolds number of the flow and causes high inaccuracy as Reynolds increases. The solution to this problem is to increase the lattice resolution in order to have lower compressibility error in the same order of spatial discretization error. However, the computational cost of a lattice Boltzmann simulation is proportional to the third power of lattice resolution 𝐶𝑃𝑈𝑡𝑖𝑚𝑒𝛼𝑁3 which restricts the choice of lattice resolution based on the available CPU power. The computational time measurements in this work suggest that lattice Boltzmann method is not an efficient choice to model swallowing simulations mainly due to the compressibility error inherent to the method.            60  References  [1]  McConnel, "Analysis of pressure generation and bolus transit during pharyngeal swallowing," The Laryngoscope, 1988.  [2]  M. A. F. Hendrix, "Relation between volume swallowed and velocity of the bolus ejected from the pharynx into the esophagus,"  Gastroenterology , 1978.  [3]  T. Yabe, "Unified solver CIP for solid. liquid, and gas," Computational Fluid Dynamics review, 1997.  [4]  D. Jacqmin, "An energy approach to the continuum surface tension method," Technical report AIAA, 1996.  [5]  S. Takagi, "Force acting on arising bubble in a quiescent liquid," ASME Fluids Engineering Division Conference, 1996.  [6]  T. Tezdayur, "Large scale fluid-particle interactions," Computational Fluid Dynamics review, 1996.  [7]  X. Shan, "Lattice Boltzmann model for simulating flows with multiple phases and components," Phys. Rev. E, 1993.  [8]  Q. Li, "Lattice Boltzmann modeling of multiphase flows at large density ratio with an improved pseudopotential model," Phys. Rev. E, 2013.  [9]  Y. Al-Jahmany, "Comparative study of lattice Boltzmann and finite volume methods for the simulation of non-Newtonian fluid flow through a planar contraction," PhD thesis dissertation, Majmaah university, 2004.  [10]  B. H. Y. P. U. Frisch, "Lattice-gas automata for the Navier Stokes equation," Phys. Rev, Apr 1986.  [11]  L. L.-S. He X, "A priori derivation of the lattice Boltzmann equation," Phys Rev, 1997.  [12]  D. d. P. L. Y. Qian, "Lattice BGK models for the navier-stokes equation," Europhys. Lett., 1992.  [13]  L.-S. L. P. L. D. d. `. Renwei Mei, "Consistent initial conditions for lattice Boltzmann simulations," Computers & Fluids , 2006.  [14]  J. J. F. J. Higuera, "Boltzmann approach to lattice gas simulations," Europhys. Lett., 1989.  [15]  C. S. J. S. M. D. Cao N, "Physicalsymmetry and lattice symmetry in lattice Boltzmann method," Phys Rev, 1997.  61  [16]  Y. P. O. d. P. J. Hardy, "Time evolution of a two dimensional model system with invariant states and time correlation functions," J. Math. Phys, 1973.  [17]  G. Z. G. R. McNamara, "Use of the Boltzmann equation to simulate lattice-gas automata," Phys. Rev, 1988.  [18]  L. L. X. He, "Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation," Phys. Rev, 1997.  [19]  J. Koelman, "A simple lattice Boltzmann scheme for Navier-Stokes fluid flow," Europhysics Letters, 1991.  [20]  M. P. K. S. W. A. V. Zecevic, "Stability and accuracy of various different schemes for the lattice Boltzmann method," Mathematical Soc, 2012.  [21]  P. L. F. Dubois, "Towards higher order lattice Boltzmann schemes," J. Stat. Mech., 2009.  [22]  D. Chandler, "Introduction to Modern Statistical Mechanics," Oxford University Press, 1987.  [23]  S. Succi, "The lattice Boltzmann equation for fluid dynamics and beyond," Oxford University Press, 2002.  [24]  N. N. Bogoliubov, "Kinetic Equations," Journal of Physics, 1946.  [25]  C. Peng, "The lattice Boltzmann Method for Fluid Dynamics: Theory and Applications," Master's thesis dissertation, EPFL, 2006.  [26]  M. Fink, "Simulation of Newtonian flows using lattice-BGK-Method," Ph.D. dissertation, Essen-Duisburg University, 2007.  [27]  R. M. L.-S. L. W. S. Dazhi Yu, "Viscous flow computations with the method of lattice Boltzmann equation," Progress in Aerospace Sciences, 2003.  [28]  L. L.-S. Lallemand P, "Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance, and stability," Phys Rev, 2000.  [29]  A. Pal, "Motility of the pharynx analyzed using lattice Boltzmann simulation," PhD Thesis, The Pennsilvania State University, 2000.  [30]  A. Caiasso, "Asymptotic Analysis of lattice Boltzmann method for Fluid-Structure interaction problems," PhD thesis, The university of Berlin, 2007.  [31]  E. Llewellin, "An extensible lattice Boltzmann framework for the simulation of geophysical flows. Part I: theory and implementation," Computers & Geosciences, 2010.  [32]  L. AJC, "Numerical simulation of particular suspensions via a discretized Boltzmann equation," Fluid Mech, 1994.  62  [33]  K. N. C. C. T. S. U. Ghia, "High Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method," Journal of Computational Physics, 1982.  [34]  M. P. A. N. D. M. Razzaghian, "Simulation of flow in lid driven cavity by MRT and SRT," International Conference on Mechanical and Robotics Engineering, 2012.  [35]  I. T. A. C. B. E. Aslan, "Investigation of the lattice Boltzmann SRT and MRT Stability for Lid Driven Cavity Flow," International Journal of Materials, Mechanics and Manufacturing, 2014.  [36]  Q. Z. a. X. He, "On pressure and velocity flow boundary conditions and bounceback for the lattice Boltzmann BGK model," Physics of fluids, 1996.  [37]  Z. Y. Michael Junk, "Outflow boundary conditions for the lattice Boltzmann method," Progress in Computational Fluid Dynamics, 2008.  [38]  H. Filippova O, "Grid refinement for lattice-BGK models," Journal of Computational Physics, 1998.  [39]  L. L.-S. S. W. Mei R, "An accurate curved boundary treatment in the lattice Boltzmann method," Journal of Computational Physics, 1999.  [40]  Q. C. Z. X. M. W. S. C. Yu Chen, "Momentum-exchange method in lattice Boltzmann simulations of particle-fluid interactions," Phys Rev, 2013.  [41]  Y. D. S. W. L. L.-S. Mei R, "Force evaluation in the lattice Boltzmann method involving curved geometry," Phys Rev, 2002.  [42]  L. L.-S. M. C. Pan C, "An evaluation of lattice Boltzmann schemes for porous medium flow simulation," Comput Fluids. in press, 2005.  [43]  G. I, "Equilibrium-type and link-type lattice Boltzmann models for generic advection and anisotropic-dispersion equation," Adv Water Res, 2005.  [44]  P. J. Dellar, "Incompressible limits of lattice Boltzmann equations using multiple relaxation times," Journal of Computational Physics, 2003.  [45]  D. Wolf-Gladrow, "Lattice-Gas Cellular Automata and Lattice Boltzmann Models - An Introduction," Springer, 2005.  [46]  S. PA, "Initial and boundary conditions for the lattice Boltzmann method," Phys Rev, 1993.  [47]  M. Kano, Y. Shimizu, "A Morphometric Study of Age-Related Changes in Adult Human Epiglottis Using Quantitative Digital Analysis of Cartilage Calcification," Cells Tissues Organs, 2005.  [48] J. Liu, "Time-implicit solution of the lattice Boltzmann equation," Master's thesis dissertation, University of Wyoming, 2008 63   [49] P. Lallemand, "Lattice Boltzmann method for moving boundaries," Journal of Computational Physics, 2003  [50] J. Liu, "A lattice Boltzmann based implicit immersed boundary method for fluid-structure interaction," Computers and mathematics with applications, 2010 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0305020/manifest

Comment

Related Items