Nicolaou, L

Zaki, TA

English

Spiral - Imperial College Digital Repository

3rd International Conference on Computational and Mathematical Biomedical Engineering - CMBE201316-18 December 2013, Hong-KongP. Nithiarasu and R. Löhner (Eds.)An immersed boundary method for patient-specific modelling of flowand aerosol deposition in the respiratory airwaysLaura Nicolaou? and Tamer A. Zaki??Department of Mechanical Engineering, Imperial College London, Exhibition Road, LondonSW7 2AZ, UK, laura.nicolaou-fernandez03@imperial.ac.uk and t.zaki@imperial.ac.ukSUMMARYAn immersed boundary (IB) method for patient-specific modelling of flow in the respiratory air-ways is presented. No-slip conditions are enforced via momentum forcing and mass conservationat immersed boundaries is satisfied via a mass source term developed for moving walls. An it-erative scheme to compute the forcing term implicitly reduces the errors at the boundary andenhances stability. In addition, the mass source term accurately accounts for airway movementduring breathing, and reduces the spurious force oscillations which typically arise in IB simula-tions of moving body problems. The method is developed for use in a generalised curvilinearsystem, which is particularly useful in bioflow problems involving complex geometries. Flow inrealistic extrathoracic airways is examined with this method.Key Words: immersed boundary, moving boundaries, curvilinear, respiratory airways.1 INTRODUCTIONUnderstanding the flow dynamics in the respiratory system is important in determining the effi-ciency of aerosolized drug delivery and the toxicity of airborne pollutants. Accurate and efficientprediction poses a challenge due to the complexity of the flow in the airways and the complexityof the airway geometries. Due to the bifurcating nature of the lung and the effect of the breathingcycle on the flow, patterns in the respiratory airways are complex even during quiet breathing. Theflow remains turbulent in the mouth and throat and the larger airways, before extensive branchingincreases the total cross-sectional area and therefore reduces the gas velocity through each airway,causing laminarisation. In most of the tracheobronchial tree, the flow is intermittently turbulentwith eddies forming at the branching points. Accurate solution of the flow field therefore motivatesthe use of direct numerical simulations (DNS) in order to resolve all the scales in the flow. Due tothe high grid resolution requirements, it is desirable to adopt an efficient computational strategy.To this end, an immersed boundary method for curvilinear coordinates is developed, which allowsthe use of structured grids to model the complex patient-specific airways. Use of structured gridsgreatly simplifies the task of grid generation, particularly for moving airway boundaries and cir-cumvents the need for remeshing at every time step which leads to more efficient computationalalgorithms. Compared to a Cartesian grid, curvilinear grids can further improve efficiency byminimising the number of grid points outside the fluid domain, can have a more natural alignmentwith the streamlines which is desirable for the accuracy of the solution, and can provide a betterwall-normal resolution.2 NUMERICAL METHOD AND RESULTSThe no-slip condition at the immersed boundary (i.e. the airway walls) is applied via a directforcing approach, which consists in adding a momentum forcing term, f , on the boundary andoutside the lung geometry. The forcing ensures that the velocity at the airway walls satisfies theboundary conditions. A mass source/sink, q, is applied to cells containing the immersed boundaryin order to ensure mass conservation. The equations are discretized on a staggered curvilineargrid using a finite volume scheme, following the method described in [1]. Time integration isperformed using a second-order fractional step method. The diffusive terms are treated implicitly414using the Crank-Nicolson scheme and the non-linear convective terms are treated explicitly usingan Adams-Bashforth scheme. The discretized equations are given byû   un 1 t=    N(un 1) +  N(un 2)   rpn 1 + 1Re ↵L(û) +  L(un 1) + fn (1)r2 n = 1 t(r · û   qn) (2)un = û    tr n (3)pn = pn 1 +  n (4)where N(u) are the convective terms, L(u) are the implicit diffusive terms and (↵,  ,  , ✓) areweighting coefficients which depend on the numerical scheme adopted. In our case, ↵ = 3/2,  =  1/2 for the Adams-Bashforth scheme and   = ✓ = 1/2 for the Crank-Nicolson scheme.2.1 Momentum forcing termThe accuracy and stability of direct forcing methods used in conjunction with the fractional stepalgorithm depend on the computation of the forcing term. In the semi-implicit fractional stepmethod, the velocities and forcing terms are coupled through the implicit diffusive terms. Kim etal. [2] proposed a solution which consisted in provisionally advancing the velocity field explicitlyin order to compute the forcing term and then adding it to the semi-implicit momentum equations.However, this scheme introduces errors near the immersed boundary as the velocities are advancedimplicitly but the forcing term is calculated explicitly. The error due to the explicit treatment ofthe diffusive terms in the calculation of fn can render the scheme unstable in two scenarios: (i) thelow-Reynolds-number flow in the terminal bronchioles; (ii) the turbulent flow in the extrathoracicairways where near-wall resolution can render the viscous stability constraint more restrictive.The solution proposed herein consists in splitting the intermediate velocity equation into two steps,ũ   un 1 t=1Re ↵L(ũ) +  L(un 1)   Gpn 1    N(un 1)   ✓N(un 2) + fn 1 (5)û   ũ t=↵Re(L(û)   L(ũ)) +  fn, (6)which fully recovers Eq. 1. In the first step, an approximation, ũ, of the intermediate velocity isdetermined. This velocity field does not satisfy the exact boundary conditions at the immersedboundary since the forcing applied is evaluated from the previous time step. In the second step,the forcing is updated and the no-slip constraint at the immersed boundary is satisfied. The exactexpression for  fn is obtained from Eq. 6 by applying the no-slip constraint at the immersedboundary (û = unIB at IB points): fn =unIB   ũ t  ↵Re(L(û)   L(ũ)) . (7)Since the diffusive term, L(û), in Eq. 7 is unknown, an iterative method is required to solve forû and  fn implicitly. The first term in the forcing (Eq. 7) drives the velocity at the boundary tothe target value (i.e. it enforces the no-slip boundary constraint) while the second term ensuresthat the diffusive term at the boundary and, hence, the velocities û at surrounding non-IB pointsare correct. Compared to the explicit forcing method, the implicit approach remains stable up toapproximately 2.5 times larger  t. The second-order accuracy of the scheme in both space andtime is demonstrated in fig. 1 using the standard decaying vortex test case [2].2.2 Mass source term for moving boundariesApplication of IB methods to moving boundary problems can lead to spurious force oscillations(SFOs) due to the violation of local mass conservation [3,4]. In the present method, an extensionof the mass source term by Kim et al. [2] for use with moving boundaries is proposed in order to41510ï2 10ï1 10010ï810ï610ï410ï26 xErrorsecond orderfirst order10ï4 10ï3 10ï2 10ï110ï810ï610ï410ï26 tErrorfirst ordersecond order(a) (b) Δx Δt Figure 1: L2 ( ) and L1 ( ) norms for the u (⇤) and v ( ) velocities versus (a) gridspacing, and (b) time step for decaying vortex problem at Re = 500.satisfy mass conservation near the boundary and suppress SFOs. The present mass source term isdesigned to enforce local and global mass conservation and is given byqn =1 V 6Xi=1 sû · n Si + un  · n   !, (8)where  V is the cell volume,  Si is the area of each cell face, n is the unit normal vector outwardat each cell face and  s is the fraction of the cell face inside the solid. The quantities u , n  and   are the surface velocity, the surface outward unit normal vector (pointing towards the fluid),and the area of the boundary within the cell, respectively.The first term on the right-hand-side of Eq. 8 mimics the cut-cell/virtual cell-merging technique[3] without requiring construction of new polyhedral cells. It ensures that the continuity equationis satisfied for the fraction of the cell inside the fluid by excluding the contribution from outsidethe lung geometry. Since two neighbouring cells which share a common face have the same valueof the surface flux, but with opposite signs, the global contribution of this term vanishes [2]. Thesecond term in Eq. 8 is an addition to the original mass source term, and ensures that the continuityequation is satisfied when the boundary is moving by including the flux at the boundary. Figure 2demonstrates the ability of the proposed mass source term to suppress spurious force oscillationsin a moving body problem.0 1 2 3 4 5−2−1.5−1−0.500.511.52t Umax/DF DFigure 2: Temporal evolution of drag force on an oscillating cylinder at Remax = 100: , mo-mentum forcing without a mass source; , momentum forcing with a mass source.2.3 Flow in the extrathoracic airwaysThe new IB method has been applied to perform the first set of fully-resolved simulations of theflow in realistic extrathoracic airways [5]. Through realistic geometric representation of the air-ways and detailed flow fields, the effect of geometric variation on the mean flow as well as the416turbulent fluctuations has been studied. Comparison of the main flow features with PIV measure-ments in the same geometries [6] shows good agreement with our numerical predictions. Figure 3shows the grid, flow field and deposition pattern in one of these geometries. The velocity profilesin the mouth are highly skewed towards the inner wall due to the airway curvature. The flowaccelerates at the back of the mouth due to the restriction in cross-sectional area, developing apharyngeal jet which impinges onto the posterior wall. Due to the bend in the airway, the flowshifts towards the outer wall, separating from the inner wall and leading to a recirculation region(fig. 3b). The maximum kinetic energy occurs in the upper pharynx near the jet (fig. 3c). Theflow field has a large effect on the particle deposition as can be observed in fig. 3d. A deposi-tion hot-spot is observed at the posterior wall of the upper pharynx where particles deposit viaimpaction due to the pharyngeal jet. The method is currently being applied to the intrathoracicairways, taking into account the airway movement during the breathing cycle.2.70 0 0.27 0.54 0.81 1.08 1.35 1.62 1.89  2.16 2.43 |ū|Umean(b) 0.48 0 0.05 0.10 0.14 0.19 0.24 0.29 0.34  0.38 0.43 kU2mean(c) (a) (d) mouth pharynx larynx trachea ξ η Figure 3: Flow in realistic extrathoracic airways. (a) Airway geometry and curvilinear grid (everyeighth grid line in ⇠ and ⌘ is plotted); (b) mean velocity magnitude; (c) mean turbulent kineticenergy; (d) particle deposition.3 CONCLUSIONSAn immersed boundary method for patient-specific simulations of flow in the respiratory airwayshas been presented. The implicit iterative scheme decreases the errors at the boundary and en-hances stability. The method can also be applied in generalized curvilinear coordinates for effi-cient modelling of the complex airway geometries. Finally, the new mass source term improveslocal mass conservation at moving geometry walls and suppresses spurious force oscillations. Themethod has been applied for detailed simulations of the turbulent flow in realistic extrathoracicairways.REFERENCES[1] M. Rosenfeld, D. Kwak, M. Vinokur, A fractional step solution method for the unsteady in-compressible Navier-Stokes equation in generalized coordinate systems, J Comput Phys, 94,102–137, 1991.[2] J. Kim, D. Kim, H. Choi, An immersed-boundary finite-volume method for simulations of flowin complex geometries, J Comput Phys, 171, 132–150, 2001.[3] J. H. Seo, R. Mittal, A sharp-interface immersed boundary method with improved mass con-servation and reduced spurious pressure oscillations, J Comput Phys, 230, 7347–7363, 2011.[4] J. Lee, J. Kim, H. Choi, K.-S. Yang, Sources of spurious force oscillations from an immersedboundary method for moving-body problems, J Comput Phys, 230, 2677–2695, 2011.[5] L. Nicolaou, T. A. Zaki, Direct numerical simulations of flow in realistic mouth-throat geome-tries, J Aerosol Sci, 57, 71–87, 2013.[6] A. Heenan, W. Finlay, B. Grgic, A. Pollard, P. Burnell, An investigation of the relationshipbetween the flow field and regional deposition in realistic extra-thoracic airways, J Aerosol Sci,35, 1013–1023, 2004.417

An immersed boundary method for patient-specific modelling of flow and aerosol deposition in the respiratory airways

http://spiral.imperial.ac.uk/bitstream/10044/1/24896/2/Nicolaou_CMBE_2013.pdf

An immersed boundary method for patient-specific modelling of flow and aerosol deposition in the respiratory airways

Abstract

Similar works

Full text

Available Versions

Spiral - Imperial College Digital Repository