AA numerical cough machine
Cesar Pairetti
Sorbonne Universit´e and CNRS, Institut Jean Le Rond d’Alembert, UMR 7190, Paris, FranceCentro Internacional Mec´anica Computacional (CONICET - UNL), Santa Fe, Argentina andFacultad de Ciencias Exactas, Ingenier´ıa y Agrimensura (UNR), Rosario, Argentina
Rapha¨el Villiers
Sorbonne Universit´e and CNRS, Institut Jean Le Rond d’Alembert, UMR 7190, Paris, France
St´ephane Zaleski
Sorbonne Universit´e and CNRS, Institut Jean Le Rond d’Alembert, UMR 7190, Paris, France andInstitut Universitaire de France, IUF, Institut Jean Le Rond d’Alembert, UMR 7190, Paris, France ∗ (Dated: January 19, 2021)We introduce a simplified model of physiological coughing or sneezing, in the form of a thin liquidlayer subject to a rapid (30 m/s) air stream. The setup is simulated using the Volume-Of-Fluidmethod with octree mesh adaptation, the latter allowing grid sizes small enough to capture theKolmogorov length scale. The results confirm the trend to an intermediate distribution between aLog-Normal and a Pareto distribution P ( d ) ∝ d − . for the distribution of droplet sizes in agree-ment with a previous re-analysis of experimental results by one of the authors. The mechanism ofatomisation does not differ qualitatively from the multiphase mixing layer experiments and simu-lations. No mechanism for a bimodal distribution, also sometimes observed, is evidenced in thesesimulations. Coughing and sneezing are two processes by which alarge number of droplets of muco-salivary fluid are ex-haled and subsequently travel large distances in the envi-ronment [1, 2]. These phenomena have acquired an acuteinterest in the context of the so-called aerosol transmis-sion route of the Covid-19 pandemic [3], but have beenstudied for near a century in the context of respiratorydiseases in general [4]. Early investigations by Duguid[5] and Loudon & Roberts [6] of the number and size ofdroplets emitted in coughing and sneezing events haveyielded rich data, incorporating very large numbers ofdroplets. A comprehensive analysis of the characteris-tics of these droplets, both in size and velocity, wouldbe of immense interest, as it is a prerequisite for themodeling of the droplet cloud propagating further down-stream from the mouth. Rapid photographic imaging [7]has revealed features similar to those observed in other atomization processes, including thin liquid sheets, liga-ments and droplets. In this context, the statistical dis-tribution of droplet sizes is of particular interest. Al-though the log-normal distribution has been frequentlymentionned in connection with exhalations [8, 9] as wellas other atomizing flows [10] many other distributions N ( d ) of the diameter d have been put forward such ascompound gamma distributions [11] and many others. Are-analysis of the data of Duguid and Loudon & Robertsfor sneezes has however revealed [3] a N ( d ) ∼ d − scalingover an impressive three orders of magnitude. This scal-ing allows determining the proportion of millimeter-sizeddroplets that travel short distances and the proportion ofmuch smaller droplets that can be incorporated in a tur-bulent puff or particle-laden jet and travel long distancesas discussed in ref. [1, 3]. The fraction of the exhaled muco-salivary fluid in each class of droplet sizes may alsobe determined in this way and the probability of havingviral loads in each class can be inferred under adequatehypotheses, such as a homogeneous distribution of thevirus in volume or surface.In order to better understand the fluid mechanics ofexhalations, King, Brock & Lundell [12] have designed aphysical model of violent exhalation that may be nick-named a “cough machine”. Air is flowing at high speed(from 10 to 30 m/s) in a rectangular-section duct, witha flow rate analog of the observed human cough. Athin layer of muco-salivary fluid is deposited at the bot-tom of the duct. While King, Brock & Lundell use thecough machine to observe non-atomizing waves on thethin layer, it can also be used to simulate droplet forma-tion at higher speeds and/or lower velocities. The devicewould then achieve atomization through a process sim-ilar to that of shear flow atomization of planar sheets[10, 11, 13–16]. (See also the recent review of numericalapproaches in [17].) However there are important dif-ferences since the liquid layer is initially at rest and theairflow is impulsive. These differences could lead to dif-ferent distributions of droplet sizes and velocities, whichare the focus of the current investigation.We model the muco-salivary fluid and the surround-ing gas as a Newtonian fluid (The muco salivary fluidis non-newtonian but these effects only kick-in at verysmall scales). Thus the flow is described by the Navier-Stokes equations, which we solve by a simple finite vol-ume discretization in the one-fluid numerical approachto two-phase flow[18], using a Volume-Of-Fluid methodfor tracking fluid interfaces[18], an octree grid and theBell-Collela-Glaz advection scheme on staggered grids as a r X i v : . [ phy s i c s . f l u - dyn ] J a n L (cid:96) x (cid:96) y (cid:96) z e h ρ l µ l ρ a µ a U σ −
30 0.03TABLE I: Dimensional values of the fluid and geometricalparameters described in [19], a Crank-Nicholson method for viscos-ity [20], and a height-function method for curvature andsurface tension [21]. The octree grid is refined using awavelet estimate of the local error, as described in theself-documented code Basilisk ( http://basilisk.fr ).We use thirteen levels of refinement at the maximum,which results in the equivalent of 2 (cid:39) or sixteenbillion grid cells. However, thanks to the octree natureof the simulation, only 10-100 million grid cells are usedin practice. The physical domain subject to the simula-tion is a cube of dimensions L with flow inlet at x = 0and flow outlet at x = L . The flow is channeled througha tube of length (cid:96) x , and rectangular cross section withheight (cid:96) y and width (cid:96) z , with walls or plates of thickness e centered on the planes x = L/ − (cid:96) z / y = ± (cid:96) y / h , density ρ l and vis-cosity µ l . The gas phase, modelling the exhaled air, hasdensity ρ a , and viscosity µ a . It enters the flow from theleft wall at x = 0 with uniform and constant velocity U through an inflow boundary condition in the region L/ − (cid:96) y / < y < L/ (cid:96) y / L/ − (cid:96) z / h 2. This condition ensures the inflow of airblows just above the initial position of the liquid layer.The surface tension of the liquid is noted σ . Gravityis neglected as it is small compared to inertial effects g(cid:96) z /U (cid:28) 1. The main dimensionless parameters are theReynolds number of the air based on the channel heightRe a = ρ a U (cid:96) z /µ a = 36 , a = ρ a U (cid:96) z /σ = 720, andthe Reynolds number of the liquid based on the heightof the liquid Re l = ρ l U h/µ l = 30 , Direct Numerical Simulations (DNS). Indeed the assessment of the DNS character ofthe simulation may be inferred from the value of the dis-sipation (cid:15) in the similar setup of refs. [10, 16]. Therethe kinetic energy dissipation per unit volume (cid:15) [22] wasmeasured and it was observed that the maximum value of (cid:15)/ ( ρ a U /(cid:96) z ) was about 0.01 which yields an estimate of (cid:15) . The Kolmogorov length scale is then η = ( ρ a ν a /(cid:15) ) / With the value of (cid:15) estimated above, we have η (cid:39) . − L (cid:39) . /η ≤ . /η (cid:39) . 75 Even if one uses the moreconservative estimate η = (cid:96) z Re − / a as in [17] one gets∆ /η (cid:39) . 4. The perhaps surprising result that the simu-lations may be qualified as DNS can be explained by thefact that while the less extreme simulations of [10, 16]were limited to (cid:96) z / ∆ = 256 grid points in the gas jetthickness, our simulations using octree refinement go upto the equivalent of (cid:96) z / ∆ = 1092.Simulations are initialized with zero velocity in the airand liquid although the incompressibility condition re-sults in a non-zero velocity field everywhere in the gasimmediately after time zero. The simulation is continuedfor a total time T=9 ms. The liquid surface is quicklysignificantly perturbed with waves present over the en-tire length (cid:96) x of the tube, with a much larger wave nearthe inlet and some secondary waves near the outlet (Fig.1). The wave stretches into a thin liquid layer that fillswith air under the combined effect of pressure and vor-ticity. Eventually, two mechanisms lead to the formationof droplets: the formation of fingers or ligaments at theend of the sheets and the puncturing of the sheets. Thelatter causes the formation of holes in the sheets andthe subsequent expansion of those holes, leading to theformation of ligaments that eventually pinch and breakinto droplets. This process has been described in otherexperimental [7] and numerical [10] investigations. Inthe current “closed channel” configuration similar mech-anisms of droplet formation are observed. It is seen onFig. 1 that the small dark red droplets near the channelexit have moved quickly since the drag law scaling for asphere F d ∼ ρ a d U / d is the droplet radius) im-plies that the small droplets accelerate much faster thanthe large ones.A clear phenomenon that couples with droplet forma-tion is growth of a large sheet or wave near the inlet andits progression downstream. The velocity of such wavesis typically the Dimotakis velocity [23] and it is expressedas U D = U √ ρ a / ( √ ρ a + √ ρ l ). A similar kind of solitarywave progressing at the Dimotakis velocity has been ob-served in the simpler setting of an infinite vortex sheetbetween air and liquid [24]. The liquid surface is excitedby a local perturbation of the flat vortex sheet, and inthis case the inhomogeneity of the flow at the entranceplays the role of the localised perturbation, while less lo-calised waves are seen further downstream. In our casethe Dimotakis velocity is U D ∼ . 004 m/s which may becompared to a rough measurement from the simulations(Fig. 2) of U D, num (cid:39) . 25 m/s. The identification ofdroplets or connected fluid components, and the compu-tation of their volumes V d allows to define an equivalentdroplet diameter d = (6 V d /π ) / . Because of the pres-ence of some non-spherical droplets in the sample, wecompared the filtered counts, incorporating only near-spherical droplets using a sphericity index. We find thatthis filtering removes droplets with d > X [m/s]-30 30U X [m/s] FIG. 1: A view of the air-liquid interface at t=9 ms. The interface is colored by the velocity of the droplets, with green thesmallest velocity and dark red the largest.FIG. 2: In a cross section of the channel, the liquid phaseand the solid are shown in black and the gas phase in white.Three different snapshots at regularly spaced time intervalsare shown. The oblique line connects the positions, thus giv-ing a graphical display of the wave velocity U D , agreing withthe Dimotakis velocity discussed in the text. frequency is shown in Fig. 3 together with the statisticalerror bars, the error being defined as one standard devi- d [ m]10 N L11L12L13Paretto, =-3.32 FIG. 3: Droplet counts in each bin, proportional to the Prob-ability Distribution Function of droplet diameters, for t = 9ms and for two grid resolutions L13: ∆ x = 18 µ m, L12:∆ x = 36 µ m and L11: ∆ x = 72 µ m. ation of the binomial. It is seen that the distribution isclose to the power law N ( d ) ∼ d − . for small sizes andthen inches down. On Fig. 3 we do not show dropletcounts for d < 100 microns and d > 400 microns, sincethe small sizes are plagued by grid resolution errors andthe larger sizes by statistical errors.An often considered number frequency (NF)model is the log-normal, which reads N ( d ) =( B/d ) exp (cid:2) − (ln d − ˆ µ ) / (2ˆ σ ) (cid:3) where B is a nor-malization constant, ˆ µ is the expected value of ln d ,also called the geometric mean , and σ is the standarddeviation of ln d , also called the geometric standarddeviation (GSD, see [8]). If we plot y = dP ( d ) versus x = ln d the Log-Normal frequency distribution appearsas a parabola. This is done in Fig. 4. That figure shows d [ m]10 N ( d ) x d log-normal, = 3.73, =0.97L11L12L13 FIG. 4: Droplet counts in “log-normal coordinates” in whichthe log-normal NF appears as a parabola, with labels L11-13as in Fig. 3. a fit with a GSD ˆ σ = 0 . 98. The dimensionless GSD isof the same order of magnitude as the GSD (ˆ σ ∼ 1) insimilar experiments and simulations [10, 25, 26], butin the authors’ opinion, this may not capture universalunderlying physics but the similar range of scales that iswithin numerical or experimental reach in the literaturecited. Comparing the NF at various resolution in [10, 25]it is seen that as the grid size is reduced the NF shiftsto the left, with its geometric mean decreasing andits GSD increasing. This also seen partially in Fig. 3where the L11 resolution peaks “before” the two others.It is thus possible that as resolution is increased thecurved NF seen on Figs 3 and 4 would progressivelyasymptote to a power law, that is a Pareto, instead ofa converged Log-Normal. Such a Pareto NF necessarilyhas a lower bound, if only the molecular size. This lowerbound is not attainable numerically or perhaps evenexperimentally. A tempting hypothesis is to associate itto the thickness of the sheets as they break or perforate,which involves mechanisms that are not modelled in thisstudy.To conclude, we have demonstrated a simple physi-cal analog of the physiological mechanism of coughingand sneezing, that is strongly reminiscent of planar sheetatomization processes. A Pareto d − distribution wasnot found but is not excluded at very small diameters.Perspectives include higher resolutions simulations and laboratory experiments in the regime of this numericalexperiment, and using the characteristics of the numeri-cally estimated droplet sizes and velocity to predict thefurther evolution of the droplet cloud using Lagrangianparticle methods such as those of Chong et al. [27].We acknowledge the funding of the Bec.ar programme,the ERC-ADV grant TRUFLOW, the ANR NAN-ODROP grant funded by the Fondation de France andthe PRACE Covid grant of computer time. We thank Ly-dia Bourouiba, Pallav Kant, Detlef Lohse and St´ephanePopinet for stimulating and fruitful discussions, and forsharing useful material, visualisations and codes. ∗ [email protected][1] L. Bourouiba, E. Dehandschoewercker, and J. W. Bush,J. Fluid Mech. , 537 (2014).[2] L. Bourouiba, JAMA , 1837 (2020).[3] S. Balachandar, S. Zaleski, A. Soldati, G. Ahmadi, andL. Bourouiba, International Journal of Multiphase Flow , 103439 (2020).[4] L. Bourouiba, Annual Review of Fluid Mechanics (2020).[5] J. P. Duguid, Epidemiology & Infection , 471 (1946).[6] R. Loudon and R. Roberts, Nature , 95 (1967).[7] B. Scharfman, A. Techet, J. Bush, and L. Bourouiba,Experiments in Fluids , 24 (2016).[8] M. Nicas, W. W. Nazaroff, and A. Hubbard, Journal ofOccupational and Environmental Hygiene , 143 (2017).[9] W. F. Wells et al., Airborne Contagion and Air Hy-giene. An Ecological Study of Droplet Infections. (Cam-bridge: Harvard University Press (for The Common-wealth Fund), Mass., USA., 1955).[10] Y. Ling, D. Fuster, S. Zaleski, and G. Tryggvason, Phys.Rev. Fluids , 014005 (2017).[11] E. Villermaux and B. Bossa, J. Fluid Mech. , 412(2011).[12] M. King, G. Brock, and C. Lundell, physiology.org ,1776 (1985).[13] M. Gorokhovski and M. Herrmann, Annu. Rev. FluidMech. , 343 (2008).[14] G. Agbaglah, R. Chiodi, and O. Desjardins, J. FluidMech. , 1024 (2017).[15] A. Zandian, W. A. Sirignano, and F. Hussain, Interna-tional Journal of Multiphase Flow , 117 (2019).[16] Y. Ling, D. Fuster, G. Tryggvason, and S. Zaleski, J.Fluid Mech. , 268 (2019).[17] C. R. Constante-Amores, L. Kahouadji, A. Batchvarov,S. Shin, J. Chergui, D. Juric, and O. K. Matar, arXivpreprint arXiv:2012.01887 (2020).[18] G. Tryggvason, R. Scardovelli, and S. Zaleski, DirectNumerical Simulations of Gas-Liquid Multiphase Flows (Cambridge University Press, 2011).[19] S. Popinet, J. Comput. Phys. , 572 (2003).[20] P.-Y. Lagr´ee, L. Staron, and S. Popinet, J. Fluid Mech. , 378 (2011).[21] S. Popinet, J. Comput. Phys. , 5838 (2009).[22] S. B. Pope, Turbulent flows (2001).[23] P. Dimotakis, AIAA J , 1791 (1986).[24] J. Hoepffner, R. Blumenthal, and S. Zaleski, Phys. Rev. Lett , 104502 (2011).[25] M. Herrmann, AAS21