曾凱磊 陳燕 譚堰琴 朱清源

東華大學(xué)環(huán)境科學(xué)與工程學(xué)院

室內(nèi)環(huán)境中的呼吸道傳染病一直是困擾人們的公共衛(wèi)生問題，比較典型的呼吸道傳染病有——從發(fā)現(xiàn)之日起一直與人類共存的流感和肺結(jié)核，曾造成巨大危害但逐步被消滅的SARS 以及2019 年年末爆發(fā)的至今仍在世界各地蔓延肆虐的新型冠狀肺炎疫情。呼吸道傳染病的病原體主要以人體呼出的飛沫為載體實現(xiàn)人群之間的傳播，而人體飛沫在室內(nèi)的傳播輸運(yùn)特性受唾液滴的釋放位置，釋放速度，粒徑大小以及室內(nèi)通風(fēng)方式的影響[1-3]，Chen 等[4]人在分析了影響飛沫液滴擴(kuò)散的眾多因素后指出室內(nèi)通風(fēng)方式和初始呼氣速度的影響占主導(dǎo)地位。當(dāng)前所進(jìn)行的飛沫液滴等顆粒污染物的擴(kuò)散研究大多基于機(jī)械通風(fēng)背景，而機(jī)械通風(fēng)背景下的飛沫等顆粒污染物的擴(kuò)散分析并不一定都適用于自然通風(fēng)，因此有必要對自然通風(fēng)背景下的飛沫液滴進(jìn)行擴(kuò)散分析。學(xué)校教室內(nèi)人員相對密集，且大多數(shù)不具備機(jī)械通風(fēng)條件，本文以自然通風(fēng)教室為背景，數(shù)值分析教室進(jìn)風(fēng)風(fēng)速對人員咳嗽飛沫液滴顆粒擴(kuò)散的影響。

1 物理數(shù)學(xué)模型

1.1 物理模型

原型為東華大學(xué)第一教學(xué)樓1107 教室，桌椅為組合式桌椅，未固定可靈活擺放，具體的擺放方式為9排5 列。根據(jù)實地測量結(jié)果對教室及教室內(nèi)桌椅、門窗進(jìn)行尺寸取整，桌椅均簡化為實體方塊，人體模型則參考相關(guān)文獻(xiàn)[5]進(jìn)行簡化，僅考慮人體模型的身軀和頭部。實際情形下教室開窗一側(cè)的外墻由磚墻和玻璃外墻組成，而開門一側(cè)的內(nèi)墻上方位置設(shè)有常年關(guān)閉的廊窗，本文暫時不考慮廊窗和外墻上懸窗的開啟。模型及相關(guān)尺寸數(shù)據(jù)如圖1 和表1 所示。

圖1 教室模型示意圖

表1 人體及教室內(nèi)相關(guān)布置物個數(shù)與尺寸

1.2 數(shù)學(xué)模型的選擇

模擬飛沫液滴擴(kuò)散的CFD 方法有歐拉-歐拉法和歐拉-拉格朗日法兩種，前者將飛沫顆粒相看為連續(xù)介質(zhì)，可以預(yù)測飛沫顆粒濃度分布，而后者將飛沫顆粒相看作離散相，可以對每一個粒子的運(yùn)動軌跡進(jìn)行跟蹤[6]。何其斌等人[7]指出由于歐拉法忽略了室內(nèi)顆粒物的顆粒特性，故只適用于模擬微小的能夠較好跟隨氣流流動的中性浮力顆粒物，而拉格朗日法則主要用于預(yù)測粒子的運(yùn)動軌跡和瞬態(tài)運(yùn)動。本文主要研究的時咳嗽飛沫液滴的瞬態(tài)擴(kuò)散過程，因此選用歐拉-拉格朗日法對咳嗽飛沫顆粒進(jìn)行軌跡跟蹤研究。對于室內(nèi)連續(xù)場，綜合計算速度和流場模擬精度要求，選擇計算速度適中且能很好地預(yù)測自然通風(fēng)空間流場的基于RANS 方程的標(biāo)準(zhǔn)k-ε 模型進(jìn)行模擬[8]。

2 邊界條件設(shè)置

本文研究的是夏季呼吸道傳染病疫情爆發(fā)期間教室前后門、4 扇平推外窗完全開啟時，自習(xí)狀態(tài)下的咳嗽飛沫液滴的瞬態(tài)擴(kuò)散。設(shè)定自習(xí)狀態(tài)下教室內(nèi)人數(shù)為七人，位置編號從1 到7，具體的人員分布如圖1所示。七人中僅有一人為潛在患者（潛在患者標(biāo)紅），忽略其他人員咳嗽或呼吸時的液滴擴(kuò)散影響，在設(shè)定的不同進(jìn)風(fēng)風(fēng)速下分別對潛在患者處于位置1、位置3時的咳嗽液滴擴(kuò)散進(jìn)行分析。正常呼吸過程中產(chǎn)生的顆粒物液滴較少對室內(nèi)影響較小，而咳嗽會產(chǎn)生更多的微?；蛞旱危@些微?；蛞旱螘鞑ズ苓h(yuǎn)的距離[9]，所以本文忽略潛在患者正常呼吸時產(chǎn)生的飛沫液滴擴(kuò)散。

2.1 連續(xù)場邊界條件設(shè)置

自然通風(fēng)空氣流動的動力為室內(nèi)外壓力差（熱壓、風(fēng)壓或二者都有），空氣流動物理過程的復(fù)雜性以及室外環(huán)境的多變性導(dǎo)致通風(fēng)率難以預(yù)測，因此，控制自然通風(fēng)以獲得所需的室內(nèi)環(huán)境條件具有挑戰(zhàn)性，并且自然通風(fēng)無法控制通過門窗口的氣流方向，室內(nèi)自然通風(fēng)流場的模擬難度較大[10-11]。而有關(guān)自然通風(fēng)的模擬方法主要有兩種，一種是建立建筑的實際模型，根據(jù)當(dāng)?shù)氐臍庀筚Y料直接給定通風(fēng)口的風(fēng)度大小和溫度。另一種是在建筑周圍加一個外場，擴(kuò)大計算模型，在外場的入口給定當(dāng)?shù)貧庀髼l件下對應(yīng)的速度大小和溫度[12]。如果對進(jìn)口風(fēng)速、風(fēng)向、溫度進(jìn)行實時模擬，不僅計算成本巨大而且也很難實現(xiàn)，因此本文采用第一種自然通風(fēng)模擬方法對進(jìn)風(fēng)邊界條件進(jìn)行簡化。根據(jù)當(dāng)?shù)貧庀髼l件假設(shè)入口進(jìn)風(fēng)平均速度（方向垂直于進(jìn)風(fēng)口），并假定進(jìn)風(fēng)參數(shù)在模擬過程中不發(fā)生變化，即在特定的通風(fēng)模式下進(jìn)行穩(wěn)態(tài)模擬。

具體數(shù)值設(shè)置參照《民用建筑供暖通風(fēng)與空氣調(diào)節(jié)設(shè)計規(guī)范》[13]選取，上海地區(qū)夏季室外平均風(fēng)速為3.1 m/s，本文分別選取1 m/s，2 m/s 的進(jìn)風(fēng)風(fēng)速進(jìn)行研究，進(jìn)風(fēng)溫度和相對濕度選用夏季通風(fēng)室外計算溫濕度，數(shù)值分別為31.2 和69%。進(jìn)風(fēng)口選用速度進(jìn)口（velocity-inlet），水蒸氣質(zhì)量分?jǐn)?shù)計算約為0.0194，出風(fēng)口（門）采用出流邊界條件（outflow），認(rèn)為出口流動達(dá)到完全發(fā)展?fàn)顟B(tài)。教室墻壁、天花板、地面、桌椅表面設(shè)置為絕熱壁面邊界條件，僅考慮人員顯熱忽略潛熱和輻射，人員皮膚發(fā)熱率取40 W/m2[14]。

2.2 離散相邊界條件設(shè)置

關(guān)于咳嗽過程的飛沫液滴，已有大量文獻(xiàn)使用不同的方法進(jìn)行了研究，加之人員個體存在差異使得研究結(jié)果略有不同。孫煒等人[15]發(fā)現(xiàn)初始直徑大于150的飛沫液滴釋放之后很快會沉降到了地板上，Chao 等人[16]通過粒子圖像測速技術(shù)（PIV）發(fā)現(xiàn)健康志愿者咳嗽時平均呼氣流速為11.7 m/s，而Zhu 等人[1]所測量的咳嗽平均氣流數(shù)據(jù)為11.2 m/s。實際咳嗽過程中的呼氣速度是變化的，本文對人體以咳嗽方式產(chǎn)生飛沫的過程進(jìn)行簡化——飛沫僅從嘴巴以水平方向咳出，飛沫速度與咳嗽呼氣速度相同，取平均流速11.7 m/s，一次咳嗽持續(xù)時間設(shè)定為0.5 s[17]，潛在患者進(jìn)行間隔時間為0.5 s 的兩次咳嗽，咳嗽噴出的氣流相對濕度設(shè)定為100%，氣流與飛沫顆粒溫度均取體溫309 K，呼出氣流中的水蒸氣質(zhì)量分?jǐn)?shù)計算約為0.0372。

本文選取粒徑范圍為2～150 μm 的飛沫液滴為研究對象，此范圍的飛沫液滴短時間內(nèi)不會迅速沉降到地面上[15]，一旦攜帶有病原體微生物則潛在危害極大。數(shù)值模擬時將飛沫液滴的物性簡化為水，密度稍微夸大設(shè)定為1003 kg/m3，蒸發(fā)溫度設(shè)定為273 K，蒸發(fā)后的液滴核粒徑設(shè)為初始液滴粒徑的32%[18]，通過計算得出蒸發(fā)性組分（水）的質(zhì)量分?jǐn)?shù)約為96.6%。結(jié)合Chao 等人[16]總結(jié)的咳嗽產(chǎn)生的不同粒徑范圍的液滴數(shù)量，依據(jù)密度估算出各個液滴組的質(zhì)量，再根據(jù)飛沫液滴對應(yīng)的質(zhì)量百分?jǐn)?shù)采用python 軟件對在2～150 μm 范圍內(nèi)非均勻分布的飛沫液滴進(jìn)行Rosin-Rammler 粒徑分布擬合，質(zhì)量分?jǐn)?shù)表和擬合圖像見表2 和圖2。

表2 一次咳嗽飛沫液滴質(zhì)量分?jǐn)?shù)表

圖2 咳嗽飛沫液滴粒徑分布擬合

此外，設(shè)定飛沫顆粒在門窗開口，嘴巴開口處的離散相邊界條件為“escape”；在室內(nèi)人體表面、桌椅表面、地板的離散相邊界條件為“trap”，即沉降到接觸表面；而各個壁面則設(shè)為彈性壁面條件“reflect”，即顆粒到達(dá)壁面會反彈回室內(nèi)空間，考慮到顆粒在壁面切向和法向上的動量損失，恢復(fù)系數(shù)均減小為0.5。

3 結(jié)果與分析

本文模擬了自然通風(fēng)教室內(nèi)兩種進(jìn)風(fēng)風(fēng)速下，潛在呼吸道疾病患者分別在位置1，位置3 進(jìn)行連續(xù)兩次咳嗽后釋放的特定粒徑范圍飛沫液滴的分布擴(kuò)散情況。由于計算機(jī)性能和內(nèi)存有限，從第一次咳嗽起共模擬了100s 的咳嗽液滴擴(kuò)散情形。數(shù)值模擬首先利用SIMPLEC 算法求解穩(wěn)態(tài)的自然通風(fēng)室內(nèi)連續(xù)相流場，之后流場中加入飛沫液滴離散相并采用PISO 算法進(jìn)行非穩(wěn)態(tài)的兩相相間耦合計算[19]，直至收斂。不同進(jìn)風(fēng)風(fēng)速下的兩種情形下的咳嗽飛沫液滴分布圖如圖3 所示。

圖3 潛在患者位于1、3 位置時不同進(jìn)風(fēng)風(fēng)速下的飛沫液滴分布圖

根據(jù)飛沫液滴分布圖可以看出，在位置1、3 處咳嗽釋放的飛沫液滴伴隨著進(jìn)風(fēng)風(fēng)速的增大，空間殘留的飛沫液滴數(shù)量明顯減少，位置1 右前方為飛沫液滴的主要影響區(qū)域，而在位置3 釋放的飛沫液滴，由于中間兩扇外窗的間距較大，位置3 沒有直接位于氣流流動路徑上，位置3 附近存在氣流回流，飛沫液滴向位置3 左右兩側(cè)空間擴(kuò)散，液滴彌漫于整個教室空間，風(fēng)速一定程度上加快了飛沫向左右側(cè)的擴(kuò)散。

當(dāng)人員咳嗽位置確定后，進(jìn)風(fēng)風(fēng)速是唯一變量，而進(jìn)風(fēng)風(fēng)速對教室空間殘存液滴數(shù)量的影響，可以用特定時刻的室內(nèi)空間液滴包（parcels）個數(shù)與第二次咳嗽結(jié)束后的室內(nèi)空間液滴包個數(shù)的比值的變化來表征，詳見圖4。每個液滴包可以包含多個但不限于整數(shù)個飛沫液滴顆粒。從圖4 可以看出，進(jìn)風(fēng)風(fēng)速增大可以加速1、3 位置處釋放飛沫液滴的空間數(shù)量減少速度，通過對飛沫液滴的沉降，捕獲，逃逸作用的促進(jìn)來實現(xiàn)這一過程。與1 位置咳嗽情形不同，3 位置處釋放的飛沫液滴在100 s 后的殘存數(shù)量隨著風(fēng)速增大明顯減少。由于教室空間較大，為了進(jìn)一步了解咳出飛沫液滴的擴(kuò)散情況，對潛在患者位于1、3 位置時不同進(jìn)風(fēng)風(fēng)速下的坐姿呼吸平面高度（y=1.2 m）上的飛沫液滴濃度云圖進(jìn)行分析，液滴濃度云圖如圖5 所示。本文對濃度云圖中的最大液滴濃度進(jìn)行了縮放（由實際值調(diào)整縮小為1.00e-10kg/m3），以便于更好的顯現(xiàn)出飛沫影響范圍隨時間和進(jìn)口風(fēng)速的變化情況。

圖4 不同進(jìn)風(fēng)風(fēng)速下的教室空間飛沫液滴數(shù)量變化情況

圖5 潛在患者位于1、3 位置時不同進(jìn)風(fēng)風(fēng)速下坐姿呼吸高度平面上的飛沫濃度云圖

從飛沫液滴濃度云圖可以看出，當(dāng)進(jìn)風(fēng)風(fēng)速為1 m/s 時，1 位置處釋放的飛沫在呼吸平面上的濃度影響大約60 s 后很輕微，而增大風(fēng)速為2 m/s 后該時間約提前到了30 s，飛沫影響區(qū)域局限于位置1 右前方。對于3 位置咳嗽情形，飛沫液滴的影響范圍更大，影響時間更長，并且飛沫對于3 位置左側(cè)附近空間呼吸高度平面區(qū)域的影響明顯強(qiáng)于右側(cè)，增大進(jìn)風(fēng)風(fēng)速同樣會加快飛沫液滴濃度向3 位置左側(cè)的擴(kuò)散。對于咳嗽位置周圍飛沫液滴的具體變化情況，本文在潛在患者左右方坐姿呼吸平面高度上各截取一個1200 mm×800 mm 長方形平面，分析截取平面上飛沫液滴平均濃度隨時間和進(jìn)風(fēng)風(fēng)速的變化情況。位置3 釋放的飛沫液滴向該位置左右兩側(cè)擴(kuò)散，并有可能對位置2，位置4 上的人員產(chǎn)生影響，在位置2，位置4 前方坐姿呼吸平面高度截取相同的長方形平面，依據(jù)是否出現(xiàn)平均濃度非零情形來評估飛沫影響。坐姿呼吸高度小截面上的飛沫液滴平均濃度變化折線圖及相關(guān)散點(diǎn)圖如圖6 所示。

圖6 坐姿呼吸高度(y=1.2 m)小截面上的飛沫液滴平均濃度變化情況

折線圖a 為潛在患者在位置1 咳嗽時兩種進(jìn)風(fēng)風(fēng)速下的右方呼吸高度小平面上的平均液滴濃度隨模擬時間的變化情況。折線圖b，c 分別為潛在患者在位置3 咳嗽時兩種進(jìn)風(fēng)風(fēng)速下的位置3 右方和左方呼吸高度小平面上的平均液滴濃度的變化情況，折線圖d，e 分別為1 m/s 和2 m/s 進(jìn)風(fēng)風(fēng)速下的位置3 左方和右方呼吸高度小平面平均飛沫液滴濃度的變化情況。散點(diǎn)圖f 為潛在患者在3 位置咳嗽時100 s 模擬時間內(nèi)飛沫液滴對位置2，位置4 上同室人員影響情況分析。

依據(jù)圖像可以看出各個呼吸高度小平面上的平均飛沫液滴濃度在1.5 s 左右到達(dá)峰值，并且峰值濃度隨進(jìn)風(fēng)風(fēng)速增大而增大。對于在位置1 咳嗽的情形，2 m/s 進(jìn)風(fēng)風(fēng)速下的平均液滴濃度變化較為劇烈，大約在4 s 時兩種風(fēng)速下位置1 右方的平均液滴濃度開始平齊。對于在位置3 咳嗽的情形，當(dāng)進(jìn)風(fēng)風(fēng)速由1 m/s增大至2 m/s 后，初始時刻位置3 右側(cè)平均液滴濃度值總體上稍有增加，而位置3 左側(cè)平均液滴濃度則在8s 之前有大幅度提升。此外，1 m/s 進(jìn)風(fēng)風(fēng)速下，初始時刻位置3 左右兩側(cè)平均飛沫液滴濃度差值較小，大約5s 以后左側(cè)液滴平均濃度開始一直高于右側(cè)，直至兩方濃度平齊趨向于0。2 m/s 進(jìn)風(fēng)風(fēng)速時，初始時刻位置3 左右側(cè)平均飛沫液滴濃度差較大，左側(cè)濃度一直高于右側(cè)，大約8s 以后兩側(cè)液滴濃度開始平齊并趨向于0。位置2 和位置4 分別位于距位置3 較遠(yuǎn)的左右兩側(cè)，根據(jù)散點(diǎn)圖f 可知：位置3 釋放的飛沫液滴會對位于位置2 和位置4 的同室人員產(chǎn)生影響，盡管影響可能很輕微。并且相較于位置2，飛沫液滴會很快擴(kuò)散至位置4 呼吸高度平面，但隨著進(jìn)風(fēng)風(fēng)速的增加，大風(fēng)速下的位置2 前方的最大平均液滴濃度要高于位置4 前方的最大平均濃度，飛沫對位置2 的影響要強(qiáng)于位置4。

4 結(jié)論

本文在模擬過程中對進(jìn)風(fēng)狀況，教室熱邊界條件，人員的呼吸狀態(tài)，飛沫液滴的物性參數(shù)均進(jìn)行了一定簡化。得出的主要結(jié)論如下：

1）進(jìn)風(fēng)風(fēng)速增大可以加快1、3 位置處釋放的飛沫液滴數(shù)量的減少。但相比3 位置處釋放的飛沫液滴，模擬時間內(nèi)增大進(jìn)風(fēng)風(fēng)速對1 位置處釋放的飛沫液滴的最終殘存數(shù)量影響不大。

2）對于位置1 情形，大風(fēng)速時位置1 右方呼吸高度小平面上的平均液滴濃度在4 s 之前始終比小風(fēng)速下的高。對于位置3 情形，隨著進(jìn)風(fēng)風(fēng)速增加飛沫液滴向左側(cè)的擴(kuò)散效應(yīng)加劇，左側(cè)平均液滴濃度在8 s之前均有大幅度提升，且同一進(jìn)風(fēng)風(fēng)速下位置3 的左側(cè)平均濃度在初始時刻一直高于右側(cè)，大進(jìn)風(fēng)風(fēng)速時左右兩側(cè)的平均飛沫液滴濃度差較大。

3）位置1 釋放的飛沫液滴僅對該位置右前方區(qū)域產(chǎn)生影響，而位置3 釋放的飛沫液滴影響范圍較廣，會對位于位置2 和位置4 的同室人員產(chǎn)生影響。進(jìn)風(fēng)風(fēng)速的增加會加快位置3 釋放的飛沫液滴向位置2，位置4 呼吸平面的擴(kuò)散，飛沫向位置4 呼吸平面的擴(kuò)散速度相對要快些，但大進(jìn)風(fēng)風(fēng)速下位置2 前方飛沫液滴最大平均濃度要高于位置4 前方的最大平均濃度，飛沫在大風(fēng)速下對位置2 的影響較強(qiáng)。