欒一剛， 殷 越， 李永曦， 符 昊， 馬正偉

(哈爾濱工程大學(xué)動(dòng)力與能源工程學(xué)院∥中國(guó)船舶集團(tuán)有限公司第七○三研究所， 哈爾濱 150001)

氣溶膠是由懸浮在氣體中的固體或液態(tài)顆粒共同組成的系統(tǒng). 在生活中，在密閉空間或通風(fēng)條件較差的空間內(nèi)，病毒氣溶膠可能以較高顆粒殘留率(數(shù)目百分比，全文同)的形式存在，封閉環(huán)境可能造成病毒的傳播.

2005年，趙鈞等[1]采用計(jì)算流體力學(xué)的方法首次分析了SARS病毒顆粒在空氣中繞建筑物后的運(yùn)動(dòng)過(guò)程和分布情況. 同年，鄧偉鵬等[2]分析了SARS病毒的特性及傳播機(jī)理，并給出了“雙層百葉風(fēng)口單側(cè)頂送，并在異側(cè)病床兩邊下部回風(fēng)”的氣流分布方式作為隔離病房設(shè)計(jì)優(yōu)選方案，為隔離病房的通風(fēng)設(shè)計(jì)理論研究奠定了基礎(chǔ). 2011年，WANG等[3]以空氣傳播感染病房為對(duì)象進(jìn)行研究，采用數(shù)值模擬研究了人的行走速率對(duì)懸浮氣溶膠數(shù)量的影響. 2012年，高乃平等[4]通過(guò)實(shí)驗(yàn)研究，建立了全尺寸置換通風(fēng)的條件，研究了人體呼出的氣溶膠顆粒在室內(nèi)垂直和水平方向的分布，在距離發(fā)生源較近的區(qū)域中，顆粒在呼吸域中的高度存在自鎖現(xiàn)象. 同年，宗青松等[5]利用數(shù)值模擬以及實(shí)驗(yàn)測(cè)試相結(jié)合的方式，首次研究了在顆粒僅來(lái)自室外的情況下送風(fēng)速率、顆粒粒徑、壁面粗糙度等因素對(duì)室內(nèi)顆粒分布和沉積的影響. 在以往的實(shí)驗(yàn)中，多數(shù)使用示蹤氣體代替固體顆粒，以便研究呼出顆粒污染物在室內(nèi)環(huán)境中的分布與傳播[6-7]. 2013年，韓云龍等[8]采用數(shù)值模擬建立了自然通風(fēng)條件下室內(nèi)不同顆粒的擴(kuò)散模型，分析了2.5、10、50 μm粒徑的顆粒物在室內(nèi)的分布情況，結(jié)果表明：小顆粒(粒徑2.5 μm)具有很強(qiáng)的氣流跟隨性，并且其分布也比中大顆粒(粒徑10、50 μm)的均勻，大顆粒由于受重力作用主要沉積于地板. 2017年，康智強(qiáng)等[9]用數(shù)值模擬方法模擬了在會(huì)議室中3種通風(fēng)方案對(duì)飛沫氣溶膠運(yùn)動(dòng)軌跡的影響，引入顆粒物在室內(nèi)滯留時(shí)間和運(yùn)動(dòng)距離的概念評(píng)價(jià)指標(biāo)，結(jié)果表明：在同側(cè)上送上回的條件下滯留時(shí)間和運(yùn)動(dòng)距離最短，是避免氣溶膠在室內(nèi)傳播的最佳氣流組織方式. 2020年，劉荔等[10]首次提出了方艙醫(yī)院在冬季條件下可能發(fā)生的熱分層現(xiàn)象導(dǎo)致發(fā)生交叉感染的新概念，研究發(fā)現(xiàn)將人員活動(dòng)區(qū)域的垂直溫度梯度降低0.42 K/m，可以降低吸入飛沫核數(shù)的60%，并建議增大新風(fēng)量. 2020年，殷平[11]首次研究了病毒和集中空調(diào)系統(tǒng)的關(guān)系，并提出采取有效措施的集中空調(diào)系統(tǒng)可大幅降低室內(nèi)的病毒殘留率. 因此，了解室內(nèi)氣溶膠顆粒運(yùn)動(dòng)擴(kuò)散特性及通風(fēng)環(huán)境對(duì)病毒殘留率的影響，有利于防止交叉感染，提高空氣質(zhì)量[12-13].

本文結(jié)合實(shí)際，通過(guò)數(shù)值仿真首次研究相對(duì)密閉空間內(nèi)局部觸發(fā)氣溶膠的彌散特性，為了防止室內(nèi)高殘留率氣溶膠存在的可能性，提出了通風(fēng)方案，對(duì)氣溶膠傳播病毒風(fēng)險(xiǎn)的預(yù)防作用具有深刻意義. 針對(duì)密閉空間室內(nèi)，人體咳嗽時(shí)氣溶膠顆粒擴(kuò)散范圍及運(yùn)動(dòng)規(guī)律進(jìn)行數(shù)值模擬，并建立3種通風(fēng)方案，對(duì)比分析方案中顆粒的滯留時(shí)間、運(yùn)動(dòng)距離和殘留率，為設(shè)計(jì)可以降低室內(nèi)高殘留率氣溶膠傳播病毒可能性的通風(fēng)方案提供參考.

1 計(jì)算模型及控制方程

1.1 封閉空間模型

建立1∶1尺寸比例的教室模型(圖1A)，室內(nèi)空間尺寸為8 m×6 m×3 m，包括門、窗、通風(fēng)口. 人體模型采用簡(jiǎn)化的方塊模型，高于周圍他人的為模擬咳嗽時(shí)局部觸發(fā)的氣溶膠源，噴口位置如圖1B所示，噴口半徑為30 mm.

圖1 密閉公共空間計(jì)算模型示意圖

1.2 具有通風(fēng)方案的空間模型

圖2為通風(fēng)方案的3種形式. 在“對(duì)側(cè)通風(fēng)”方案中，單側(cè)窗戶為進(jìn)風(fēng)口，對(duì)側(cè)門為排風(fēng)口；在“上送兩側(cè)回風(fēng)”方案中，頂部通風(fēng)口進(jìn)風(fēng)，窗口與門為排風(fēng)口；在“單側(cè)通風(fēng)-對(duì)側(cè)/上側(cè)回風(fēng)”方案中，單側(cè)窗戶為進(jìn)風(fēng)口，頂部通風(fēng)與門為排風(fēng)口.

圖2 公共場(chǎng)所通風(fēng)方案

1.3 控制方程

采用連續(xù)相氣流流場(chǎng)和顆粒相耦合的方法進(jìn)行數(shù)值模擬，氣流場(chǎng)采用室內(nèi)氣流流場(chǎng)效果較好的Realizablek-ε湍流模型計(jì)算[14]，通用方程為：

(1)

其中，ρ為氣體密度，Φ為通用變量，可以代表速率分量u、v、w和熱力學(xué)溫度T等求解變量，t為時(shí)間，?為梯度，K代表廣義源項(xiàng).

對(duì)顆粒相的計(jì)算采用離散相顆粒模型(DPM)，該方法是在拉格朗日法的基礎(chǔ)上建立的，因此可以追蹤顆粒的運(yùn)動(dòng)軌跡，可以計(jì)算兩相耦合和單相耦合問(wèn)題，在忽略顆粒對(duì)氣體相的影響、只考慮氣體對(duì)顆粒的影響情況下，氣固相間作用力采用曳力模型描述，對(duì)于單個(gè)顆粒的運(yùn)動(dòng)方程，由牛頓第二運(yùn)動(dòng)定律可得在笛卡爾坐標(biāo)下的表達(dá)式：

(2)

其中，ρp為顆粒密度，ρg為氣相密度，g為重力加速率，up為顆粒相速率，ug為氣相速率，F(xiàn)D為流體對(duì)顆粒的曳力，F(xiàn)x為除重力以外的其他力.Fx包括：流場(chǎng)中由于速率梯度而產(chǎn)生的saffman力、顆粒在流體中加速或者減速而存在的Basset力、流場(chǎng)中顆粒的旋轉(zhuǎn)而產(chǎn)生的Magnus力(附加質(zhì)量)、熱泳力、布朗力以及顆粒之間和顆粒在與壁面碰撞的力.

在室內(nèi)顆粒物的運(yùn)動(dòng)受力中，壓力梯度力、虛擬質(zhì)量力、由于顆粒旋轉(zhuǎn)而產(chǎn)生的Basset力均比曳力小2個(gè)量級(jí)，因此考慮由于速率梯度產(chǎn)生的saffman力、熱泳力和布朗力的作用，方程(2)可簡(jiǎn)化為：

(3)

其中，F(xiàn)b為布朗力，F(xiàn)th為熱泳力，F(xiàn)s為 saffman力.

2 網(wǎng)格劃分和邊界條件

2.1 網(wǎng)格劃分

采用ICEM軟件對(duì)計(jì)算區(qū)域創(chuàng)建非結(jié)構(gòu)化網(wǎng)格，并在顆粒源噴射口附近進(jìn)行網(wǎng)格加密，計(jì)算網(wǎng)格如圖3所示. 由于數(shù)值模擬使用的非穩(wěn)態(tài)計(jì)算方法，在有限的計(jì)算資源前提下保證計(jì)算速率，建立了網(wǎng)格數(shù)分別為77萬(wàn)、151萬(wàn)、186萬(wàn)和246萬(wàn)的網(wǎng)格進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證. 在室內(nèi)空間取一條中心垂直線(X=0、Z=0)，得到1 s時(shí)刻的速率分布(圖3)．

圖3 計(jì)算域網(wǎng)格

從圖4可以看出，在室內(nèi)空間,中心線高度隨顆粒速率的變化率從0開(kāi)始先迅速增大，達(dá)到峰值之后減小，當(dāng)速率低于0.003 m/s后略有提升，隨后速率迅速降低到0. 在4種網(wǎng)格空間中，網(wǎng)格數(shù)為186萬(wàn)的中心線高度變化速率分布與77萬(wàn)的相比有較大差距，151萬(wàn)網(wǎng)格的次之，但網(wǎng)格數(shù)為186萬(wàn)的速率分布曲線與246萬(wàn)的分布結(jié)果并無(wú)明顯變化. 因此，網(wǎng)格數(shù)為186萬(wàn)的網(wǎng)格敏感性已經(jīng)達(dá)標(biāo)，最終選擇186萬(wàn)網(wǎng)格數(shù)作為全局網(wǎng)格劃分方案的參數(shù).

圖4 4種不同網(wǎng)格數(shù)空間的中心線高變化速率分布

2.2 邊界條件

數(shù)值模擬研究由局部觸發(fā)的氣溶膠彌散特性，獲得氣溶膠顆粒在特定時(shí)間的位置和殘留率等信息，因此采用Fluent軟件進(jìn)行非穩(wěn)態(tài)計(jì)算. 簡(jiǎn)化人物模型的呼吸方式，設(shè)置空間進(jìn)氣條件，模擬呼氣中咳嗽產(chǎn)生顆粒的氣流. 設(shè)定1 s內(nèi)由咳嗽觸發(fā)的氣溶膠噴射條件，顆粒物采用離散相模型(DPM)進(jìn)行計(jì)算，顆粒物設(shè)置為水顆粒. 湍流模型選擇Realizable模型計(jì)算[15-16]，采用Coupled算法，離散格式為二階迎風(fēng)格式，單時(shí)間步長(zhǎng)內(nèi)收斂殘差低于1×10-6數(shù)量級(jí)，非穩(wěn)態(tài)數(shù)值模擬設(shè)置時(shí)間步長(zhǎng)為0.2 s，步數(shù)共300步，即模擬1 min氣溶膠顆粒殘留率的變化情況.

2.2.1 密閉空間的邊界條件 氣溶膠擴(kuò)散中，由于氣體顆粒尺寸不同，大質(zhì)量顆粒受重力的影響較大，因此計(jì)算域設(shè)置必要的重力環(huán)境條件. 具體邊界條件如表1所示，人體正常體溫為310 K(37 ℃)左右，但模擬時(shí)人體正常外表溫度為304 K(31 ℃)，空間環(huán)境為室內(nèi)常溫293 K(20 ℃). 本文主要模擬人體咳嗽后氣溶膠顆粒的彌散特性，因此人體口腔設(shè)置速率邊界條件：速率為1.3 m/s，氣流溫度為308 K(35 ℃)[17]. 咳嗽噴出的氣溶膠為氣液兩相混合物，設(shè)置水顆粒物性參數(shù)：直徑為0.5～100 μm，平均粒徑為10 μm. 設(shè)定沿口腔法向方向噴出顆粒物：顆粒物初速率為10 m/s，溫度為310 K(37 ℃). 設(shè)置周圍人體為粒子捕捉邊界條件：墻壁為粒子逃逸邊界條件.

表1 密閉空間邊界條件參數(shù)Table 1 The boundary condition parameters of a confined space

2.2.2 通風(fēng)方案的參數(shù)設(shè)定 本文旨在研究3種通風(fēng)方案對(duì)密閉氣溶膠殘留率隨時(shí)間擴(kuò)散規(guī)律的影響，因此采用瞬態(tài)計(jì)算，得到特定時(shí)間點(diǎn)顆粒物的狀態(tài). 在密閉空間中，簡(jiǎn)化人體模型均具有一定頻率的呼吸，且模擬人體口腔咳嗽產(chǎn)生1 s時(shí)長(zhǎng)的顆粒物. 呼吸、房?jī)?nèi)空氣設(shè)定為不可壓縮空氣，顆粒物為液態(tài)水顆粒.

計(jì)算湍流模型采用Realizablek-ε模型和Coupled算法. 單步計(jì)算收斂殘差低于1×10-6. 設(shè)定瞬態(tài)求解參數(shù)，時(shí)間步長(zhǎng)為0.2 s，共計(jì)算300 步，即1 min內(nèi)在通風(fēng)條件下顆粒的彌散規(guī)律.

方案1：對(duì)側(cè)通風(fēng).

邊界條件與密閉空間邊界條件相同，包括：墻壁、人體溫度、呼吸氣流、顆粒物條件. 通風(fēng)下在保證通風(fēng)流量相同的情況下設(shè)置通風(fēng)窗口為速率入口，速率大小為0.2 m/s，溫度為290 K(17 ℃). 對(duì)側(cè)門為回風(fēng)口，設(shè)置為壓力出口，出口壓力為標(biāo)準(zhǔn)大氣壓(101.325 kPa). 通風(fēng)主要從窗口流入，從門口處流出，所有邊界條件見(jiàn)表2.

表2 對(duì)側(cè)通風(fēng)空間邊界條件參數(shù)

方案2：上送兩側(cè)回風(fēng).

除通風(fēng)外其他邊界條件與上述密閉空間的邊界條件相同，在保證通風(fēng)流量相同的情況下設(shè)置頂部送風(fēng)口為速率入口，速率大小為1.3 m/s，溫度為290 K(17 ℃). 4個(gè)窗口和2個(gè)門為回風(fēng)口設(shè)置為壓力出口,出口壓力為標(biāo)準(zhǔn)大氣壓. 因此通風(fēng)主要從頂部窗口流入，從兩側(cè)窗口和門口處流出，所有邊界條件見(jiàn)表3.

表3 上送兩側(cè)回風(fēng)通風(fēng)空間邊界條件參數(shù)

方案3：?jiǎn)蝹?cè)通風(fēng)-對(duì)側(cè)/上側(cè)回風(fēng).

除通風(fēng)條件外，邊界條件與方案(1)、方案(2)的條件相同. 在保證通風(fēng)流量相同的情況下設(shè)置4個(gè)窗口為速率入口，速率大小為0.2 m/s，溫度為290 K(17 ℃). 頂部窗口和2個(gè)對(duì)側(cè)門為回風(fēng)口設(shè)置成壓力出口：出口壓力為大氣壓. 因此通風(fēng)主要從4個(gè)窗口流入，從頂部窗口和對(duì)側(cè)2個(gè)門口處流出，所有邊界條件見(jiàn)表4.

表4 單側(cè)通風(fēng)-對(duì)側(cè)/上側(cè)回風(fēng)通風(fēng)空間邊界條件參數(shù)

3 結(jié)果與討論

3.1 密閉空間的顆粒彌散分布

在密閉空間中，顆粒從口腔噴出后60 s內(nèi)的顆粒彌散結(jié)果如圖5所示，僅有較大顆粒(粒徑>10 μm)落到前方，其他較小顆粒(粒徑<10 μm)隨空氣流動(dòng)，而由于室內(nèi)回流，將顆粒附著在自身上. 當(dāng)顆粒噴出5 s時(shí)，顆粒僅隨氣流運(yùn)動(dòng)到身體前上方，擴(kuò)散性不強(qiáng)，顆粒團(tuán)處于聚集狀態(tài). 當(dāng)噴出后15 s時(shí)，小部分顆粒滯留在原地，但絕大多數(shù)顆粒繼續(xù)向前上方運(yùn)動(dòng)，此時(shí)顆粒團(tuán)范圍擴(kuò)大. 當(dāng)噴出后25 s時(shí)，顆粒運(yùn)動(dòng)到空間頂部，由于頂部墻壁的物理邊界壓迫，顆粒由團(tuán)狀散開(kāi)成片狀，隨后在35 s時(shí)刻顆粒呈片狀向前飄動(dòng)，在45、55、60 s時(shí)刻，上方顆粒到達(dá)空間頂部前上方邊界，并向兩側(cè)和后方下落移動(dòng).

圖5 在60 s內(nèi)密閉空間中顆粒的彌散分布

氣溶膠顆粒傳播最遠(yuǎn)距離與殘留率隨時(shí)間的變化如圖6A所示,從口腔噴射而出的顆粒在最開(kāi)始噴出的時(shí)間段內(nèi)移動(dòng)的速率最快，隨著時(shí)間推移，移動(dòng)逐漸變緩，在噴出后55 s，顆粒到達(dá)空間中前方壁面并向后方空間擴(kuò)散，最大傳播距離為5.3 m，與顆粒彌散圖的顆粒運(yùn)動(dòng)現(xiàn)象相符. 由于人體和四周壁面DPM邊界條件分別設(shè)置為Trap和Escape，顆粒的殘留率不斷降低，剛開(kāi)始噴出顆粒存在于空氣中，并未到達(dá)墻壁及人體表面，僅大顆粒落于地面，因此在前15 s內(nèi)顆粒殘留率變化較小，此后顆粒接觸上方和前方壁面，殘留率下降的速率明顯，最終在60 s時(shí)刻，殘留率為70.86%.

在噴出后60 s時(shí)刻，密閉空間不同位置的溫度分布如圖6B所示，位置分別為X=0 m和Z=3.75 m處(距前壁0.5 m處). 由于人體表面附近溫度較周圍環(huán)境的高，氣流向上運(yùn)動(dòng)，產(chǎn)生了熱羽流現(xiàn)象，并在空間前方和后方壁面附近形成環(huán)流區(qū). 在空間兩側(cè)，空氣沿壁面向下流動(dòng)，形成室內(nèi)環(huán)流，顆粒隨氣流運(yùn)動(dòng)的影響大，造成顆粒在密閉空間內(nèi)無(wú)法排出.

3.2 不同通風(fēng)方案的仿真結(jié)果

3.2.1 對(duì)側(cè)通風(fēng)方案 圖7A為對(duì)側(cè)通風(fēng)顆粒傳播最遠(yuǎn)距離、殘留率隨時(shí)間的變化關(guān)系. 由于通風(fēng)的關(guān)系，最遠(yuǎn)距離在45 s時(shí)到達(dá)峰值，隨后顆粒向后方擴(kuò)散，相對(duì)密閉空間的提前了10 s，但是最大傳播距離為5.7 m. 對(duì)側(cè)通風(fēng)方案中室內(nèi)顆粒殘留率隨時(shí)間的變化趨勢(shì)與密閉空間的相同(圖6A)，但隨時(shí)間的不同，顆粒殘留率在25 s之前并未改變，說(shuō)明顆粒未被人體捕捉. 隨后顆粒接觸到壁面停滯或被通風(fēng)氣流帶走，60 s后室內(nèi)顆粒殘留率為65.4%，少于密閉空間的殘留率(圖6A).

圖7B為60 s時(shí)刻對(duì)側(cè)通風(fēng)條件下空間的溫度分布，由于通風(fēng)效果室內(nèi)溫度變得較均勻，并且可以看出室內(nèi)氣流朝門口流動(dòng)，顆粒更容易被氣流帶出室外.

3.2.2 上送兩側(cè)通風(fēng)方案 圖8A為上送兩側(cè)回風(fēng)顆粒物傳播最遠(yuǎn)距離與殘留率隨時(shí)間的變化，可以發(fā)現(xiàn)在60 s內(nèi)顆粒傳播距離越來(lái)越遠(yuǎn)，但是傳播速率變緩，最大傳播距離為5.23 m. 在前35 s，顆粒殘留率的下降速率變慢，在35 s后迅速下降，并且在60 s時(shí)刻，室內(nèi)顆粒殘留率為63.22%.

圖8B為60 s時(shí)刻上送兩側(cè)回風(fēng)條件下空間的溫度分布，由于上通風(fēng)方案中室內(nèi)溫度中間較低，上方送風(fēng)的冷氣流和人體表面附近的熱羽流作用，在室內(nèi)前后方形成環(huán)流，并且室內(nèi)氣流朝門口和窗口流動(dòng)，顆粒更容易被氣流帶出室外.

3.2.3 單側(cè)通風(fēng)-對(duì)側(cè)/上側(cè)回風(fēng)方案 圖9A為上送兩側(cè)回風(fēng)顆粒物傳播最遠(yuǎn)距離與殘留率隨時(shí)間的變化，可以發(fā)現(xiàn)從0～45 s內(nèi)顆粒傳播距離越來(lái)越遠(yuǎn)，但是傳播速率變緩，45～60 s的傳播距離下降，最大傳播距離為5.49 m. 從15 s之后顆粒物殘留率明顯下降，但在60 s時(shí)刻，室內(nèi)顆粒物殘留率依然較高(70.32%).

圖9 單側(cè)通風(fēng)-對(duì)側(cè)/上側(cè)回風(fēng)條件下顆粒的移動(dòng)特性與溫度分布

圖9B為單側(cè)通風(fēng)-對(duì)側(cè)/上側(cè)回風(fēng)條件下60 s時(shí)刻空間的溫度分布，由于對(duì)比對(duì)側(cè)通風(fēng)模式，多了上方出口，氣流一部分從上方流出，但還存在較強(qiáng)的室內(nèi)環(huán)流，導(dǎo)致在60 s時(shí)刻顆粒殘留率仍然較高.

3.3 方案測(cè)評(píng)

圖10為1 min內(nèi)3種不同通風(fēng)方案室內(nèi)顆粒殘留率的變化情況，可以看出在對(duì)側(cè)通風(fēng)方案顆粒殘留率穩(wěn)步下降，但是上送兩側(cè)回風(fēng)方案在60 s時(shí)刻顆粒殘留率降低到較低的水平，比對(duì)側(cè)通風(fēng)方案的低2.18%；單側(cè)通風(fēng)-對(duì)側(cè)/上側(cè)回風(fēng)方案顆粒殘留率下降較慢. 可見(jiàn)，上送兩側(cè)回風(fēng)方式可以有效降低室內(nèi)氣溶膠殘留率，降低氣溶膠的干擾風(fēng)險(xiǎn).

圖10 3種不同通風(fēng)方案室內(nèi)顆粒殘留率隨時(shí)間的變化

4 結(jié)論

通過(guò)數(shù)值仿真方法，開(kāi)展在密閉空間由局部觸發(fā)的氣溶膠彌散的仿真研究，結(jié)果表明：大顆粒受重力影響落于地面，小顆粒隨室內(nèi)環(huán)流擴(kuò)散. 在顆粒物噴出后60 s時(shí)，有70.86%數(shù)量的顆粒依然存在于室內(nèi)空氣中. 因此，在室內(nèi)存留的高殘留率氣溶膠，容易引發(fā)氣溶膠感染風(fēng)險(xiǎn)，特別是顆粒源的前上方殘留率最大的區(qū)域. 基于密閉空間提出了3種通風(fēng)方案，在保證通風(fēng)氣流流量相同的情況下，上送兩側(cè)回風(fēng)方案顆粒傳播距離短，同時(shí)室內(nèi)顆粒殘留率較低，是較理想的密閉空間氣溶膠通風(fēng)方案. 對(duì)側(cè)通風(fēng)是常用的通風(fēng)方案，可以保證氣溶膠殘留率在室內(nèi)穩(wěn)步降低，并且氣流紊流度低，通風(fēng)體感較好，可以在平時(shí)采用對(duì)側(cè)通風(fēng)方案，必要時(shí)采用上送兩側(cè)回風(fēng)的通風(fēng)方案以保證較低的顆粒物濃度.