液橋重開過程的兩相格子Boltzmann模型及計算

王寧寧,劉海湖,張楚華

(西安交通大學能源與動力工程學院,710049,西安)

?

液橋重開過程的兩相格子Boltzmann模型及計算

王寧寧,劉海湖,張楚華

(西安交通大學能源與動力工程學院,710049,西安)

根據最近提出并發展的相場格子Boltzmann方法,建立了阻塞氣道重開過程的不混溶兩相流動模型和計算方法。基于自由能理論、引入指示函數對兩相流體界面進行描述,指示函數的演化遵循Cahn-Hilliard方程,具有堅實的物理基礎;通過壓力分布函數對流場信息進行求解,可有效降低密度梯度離散所誘發的數值不穩定性;引入勢能形式的界面張力項,與壓力形式的相比可有效抑制界面處的假擬速度。應用該模型對液橋重開的兩相流動過程進行了數值研究,并著重分析了毛細數對阻塞液橋演化過程的影響,結果表明,阻塞液橋的軸向厚度變化存在臨界毛細數現象,當毛細數大于臨界值時,阻塞液橋的軸向厚度隨時間逐漸減小,最終發生破裂,實現氣道的重開;利用該模型重現了無液滴和有液滴形成的兩種氣道重開現象,對比研究發現,在無液滴形成的氣道重開過程中,壁面經歷的壓力變化更大。研究工作可為深入認識人體肺氣道的生理病理機制、微通道內不混溶兩相流體的流動規律提供一定的理論依據。

格子Boltzmann方法;兩相流動;微通道;界面張力;液橋重開

航天醫學工程領域中,太空失重條件下人體肺呼吸道阻塞與重開現象的重新分配是引起航天員心肺功能異常、血氣功能交換障礙的主要原因,相關生理病理機制吸引了很多學者的關注[1]。肺呼吸道內的空氣、吸附液膜及管壁組成了復雜的氣液彈動力系統,阻塞液橋的破裂將嚴重影響兩相流動狀態及空間分布[2]。由于呼吸道解剖結構的復雜性、管壁及其吸附液膜多變的流變性質,所以與阻塞及重開現象相關的呼吸流動機理研究對理論、實驗及數值方法均提出了很大的挑戰。

Prisk等通過空間實驗室研究發現,微重力條件下呼吸道阻塞的臨界肺氣量與地面環境下的測量值相當,但阻塞發生后將遍及肺部其他區域[3]。Pedley等采用生物力學方法對可塌陷管、槽道內流體/彈性體動力學系統進行了理論及數值研究[4-5],但是未考慮界面張力等因素的影響。Grotberg等采用直接數值模擬方法,對定常界面張力作用下的氣道閉合現象進行了分析,發現該現象會導致壁面剪切應力劇烈波動[6]。許世雄等采用剛性直圓管對下呼吸道阻塞液橋重開過程進行了理論和實驗研究[7-8]。Malashenko等將肺呼吸道簡化建模為剛性直圓管道,并采用VOF(Volume of Fluid)方法研究了阻塞液橋的遷移及破裂過程[9]。Baudoin等實驗研究了直管道及多分岔管道中單個或多個液橋同時存在時的流動行為[2]。

格子玻爾茲曼方法(LBM)可以方便地對兩相界面進行描述,易于處理復雜幾何邊界,近年來在微尺度多相流動領域取得了顯著的進展[10-12]。本文通過相場LBM兩相流動模型[10-11]對阻塞氣道的重開過程進行了數學建模和數值計算,重點分析了毛細數對阻塞液橋軸向厚度的影響規律,以及液橋破裂過程中的兩相流動現象,以期為深入認識呼吸道阻塞與重開機理、管內兩相流體輸運技術提供理論參考。

1 液橋重開過程建模

1.1 簡化條件

阻塞氣道重開現象的本質是管道內氣液兩相流動及其流型的轉變過程,阻塞狀態時為彈狀流動,重開狀態時為環狀流動。為了簡化該流動物理過程,進行了如下假設[9]:管道壁面是剛性的;忽略液膜的分層及纖毛特性,將液膜簡化為牛頓流體,且初始厚度呈均勻分布;不考慮Marangoni應力的作用;初始阻塞液橋形狀為圓柱形;附著液膜運動及內部空氣流均為不可壓縮流動。此外,本文計算條件下Bond數(重力與界面張力之比)遠小于1,因此模擬中不考慮重力的作用。

1.2 數學模型及參數條件

阻塞細支氣管道二維幾何模型[9]如圖1所示。深色區域為壁面吸附液膜及阻塞液橋,無色區域為被液橋隔斷的空氣區域,圖中氣道為阻塞狀態,(xl,yl)為計算域的大小,U為平均氣體流速,R為氣道半徑,h0為通道中液膜的初始厚度,b0為液橋的初始厚度,F0為氣體水平方向上施加的外力并驅動兩相運動。上、下壁面為無滑移邊界條件,入口及出口為周期邊界條件。

圖1 阻塞氣道幾何模型圖

液橋軸向厚度演化主要與7個無量綱參數相關,即

(1)

式中:下標a表示氣相,m表示液膜;t為時間;η為動力學黏度;ρ為密度;Rea為肺氣道中空氣雷諾數。在人體細支氣管道內部Rea=O(1),因此可認為氣道中的流動基本為蠕變流,Rea對人體氣道流動的影響不大。此外,細支氣管內氣相、黏液的密度和黏度以及通道內的氣體流速通常變化不大,而界面張力系數受年齡、呼吸道疾病等因素的影響具有一定的變化范圍。毛細數Ca可以描述界面張力的作用,Ca=ηmU/σ為黏性力與界面張力之比。一般認為,人體細支氣管內空氣-黏液的界面張力系數σ=0.01～0.1 N·m,對應毛細數Ca=0.1～1。

2 相場LBM兩相流動模型

相場LBM兩相流動模型,基于自由能的概念、通過指示函數的演化方程自動捕捉兩相界面,模型物理意義清晰,同時具有LBM方法本身的各項優點,如易于處理復雜幾何邊界、程序簡單、并行效率高等,因此在模擬人體肺氣道阻塞與重開過程中的復雜兩相流動問題方面,具有良好的應用潛力。

(2)

(3)

[η(u+uT)]+FS+F0

(4)

μ=φ3-φ-ε22φ

(5)

式中:u為流體的宏觀速度;Mφ為移動性參數;p為壓力;FS為與界面張力有關的作用力項[12],當界面張力為常數時FS=3×21/2σμφ/4ε;μ為化學勢,是指示函數和界面厚度的函數。

(6)

(7)

(8)

(9)

(10)

(11)

(12)

式中:ξ為微觀粒子速度;I為二階單位張量;ωα為各離散速度的權系數(ω0=4/9,ω1,2,3,4=1/9,ω5,6,7,8=1/36);Γ為與Mφ有關的可調參數,Mφ=δtΓ(τg-1/2)。宏觀量的計算式為

(13)

(14)

(15)

為模擬兩種動力學黏度和密度不同的流體,將運動學黏度υ和密度定義為指示函數的線性函數

(16)

(17)

3 液橋重開過程的模擬結果與分析

3.1 毛細數對液橋演化過程的影響

在Ca=0.1～1.0的范圍內,對不同毛細數下的液橋演化規律進行了數值模擬,參數設置為R=50,h0/R=0.2,b0/R=2.0,Re=0.67,其他的與文獻[9]的參數設置相同。由計算結果發現,兩相演化過程可大致分為3類,如圖2所示。

(a) Ca=0.125, (b) Ca=0.182, (c) Ca=0.5, 軸向變厚 軸向不變 軸向變薄圖2 3種毛細數下的液橋軸向厚度演化過程

在演化過程中,液橋整體沿流動方向遷移,各毛細數下的初始柱形液橋均快速發展為更接近實際的彎月形液橋,初始階段受毛細數的影響不大。不同毛細數下阻塞液橋的軸向厚度演化呈現出不同的演化規律,即存在臨界毛細數現象。本文模擬獲得的臨界值為Cacr=0.182,在該臨界毛細數下阻塞液橋的軸向厚度在演化后期保持不變,當毛細數大于臨界值時,阻塞液橋的厚度隨時間減小,反之則增大。

(a)本文模擬結果

(b)文獻[9]模擬結果圖3 不同毛細數下液橋軸向厚度演化情況

不同毛細數下液橋軸向厚度b的變化趨勢如圖3所示。由圖3a模擬結果可見,各毛細數下,初始階段形成的液橋厚度幾乎相等,本文將該軸向厚度定義為b1,并用于其他時刻軸向厚度的無量綱化處理。為方便結果分析,引入新變量t′=tphy-0.008 s對本文的時間項進行處理,使得t′=0時刻與文獻[9]中tphy=0時刻相對應。對比可知:本文模擬所得的臨界毛細數與文獻[9]的模擬結果相同,均為0.182;在非臨界毛細數下,相場LBM獲得的液橋厚度演化速率略偏小,總體演化趨勢均呈近線性發展,與文獻[9]的模擬結果定性地吻合良好,從而驗證了相場LBM的有效性和準確性。

為了探究液橋厚度發生變化的原因,本文以毛細數大于臨界值的情況為例,給出了液橋及其臨近區域的壓力分布云圖,如圖4所示。由壓力分布云圖可得,以液橋長度在中軸線處的中點B為界,在毛細數大于臨界值時液橋上游的壓差PAB小于下游的壓差PBC,因此相同時間內流入液橋的流體體積小于流出液橋的流體體積,宏觀上表現為液橋厚度變薄。其他毛細數下可同理分析。

圖4 毛細數大于臨界值時的液橋壓力分布云圖

基于以上結論可知,毛細數高于臨界值時,阻塞液橋厚度逐漸減小,且毛細數越大,液橋厚度變化速率越快。因此,當液橋厚度減小至某一臨界狀態后,液橋的形狀將無法繼續維持,氣道發生重開。

3.2 無液滴形成的液橋重開過程

給定h0/R=0.26,b0/R=0.33,F0=1.5×10-8,對應初始階段Ca=0.265,Re=1.24。所得模擬結果與文獻[9]的模擬結果對比如圖5所示。

(a)本文模擬結果 (b)文獻[9]模擬結果圖5 無液滴形成的液橋重開過程

圖5顯示,在兩相流體沿流動方向遷移的過程中,柱形的阻塞液橋在初始階段迅速發展為彎月形,并且液橋左側界面的曲率明顯大于右側。液橋軸向厚度逐漸變薄,中心處變薄的速率最快且最先發生斷裂;隨著裂口的繼續擴大,斷裂的液橋徑向長度減小,并與上、下壁面附著的貼壁液膜逐漸融合,進而被隔斷的兩側流體重新連通,即氣道發生重開,重開過程中無液滴生成,該結果與文獻[9]吻合良好。

3.3 有液滴形成的液橋重開過程

由文獻[2,9,15]知,若液橋體積和毛細數足夠大,則液橋重開過程中可能伴隨有微液滴生成。為此,設置h0/R=0.2、b0/R=0.83、F0=1.5×10-8,對應初始階段Ca=1.865、Re=1.867進行模擬,并與文獻[9]模擬結果、文獻[2]實驗結果進行了比較,如圖6所示。

(a)本文模擬結果

(b)文獻[9]結果 (c)文獻[2]結果圖6 有液滴形成的液橋重開過程

圖6表明,在大初始阻塞液橋體積和毛細數下,本文模擬獲得了有液滴形成的液橋重開現象,所得演化過程與文獻[9]及文獻[2]結果吻合良好。同時看到,柱形液橋很快發展為彎月形,在界面張力及外力場的共同作用下,液橋左側的氣相流體尖端被“拉長”,液橋的軸向厚度逐漸變薄。其與無液滴形成的情況不同在于,左右兩側流體的最初接觸位置并非處于氣道中軸線,而是在中軸線的兩側各形成了一個裂口,兩裂口之間的殘余液膜形成了一個小液滴。

3.4 2種液橋重開過程的對比分析

比較以上2種液橋重開現象可知:在液橋破裂之前的短時間內,若液橋左側界面的曲率半徑小于右側界面的曲率半徑,則液橋斷裂直接在中軸線處發生,僅形成了一個裂口,不形成小液滴的概率很大;反之,液橋兩側流體在連通時存在2個或多個接觸點,接觸點之間的液膜被截斷后并未及時排出,從而形成了小液滴,該液滴在流動過程中進一步受到主流流動影響,可能破裂為更多的小液滴。此外,液滴生成也取決于液橋破裂前不穩定[9]等因素。

(a)無液滴形成的液橋重開過程

(b)有液滴形成的液橋重開過程圖7 2種液橋重開過程中壁面的壓力分布變化

為分析管壁在氣道重開過程中的受力變化,本文給出了2種液橋重開現象中液橋破裂前后氣道壁面上的壓力分布變化,如圖7所示。無液滴及有液滴形成的液橋重開現象對應的重開時間分別為tphy=0.008 5 s和tphy=0.017 3 s。在無液滴形成的重開過程中,液橋破裂前后其所在區域出現較大的壓力突變;在有液滴形成的氣道重開過程中,壓力幅值在所關注的位置區段內整體增長較小。因此,無液滴形成的氣道重開過程對管壁有更大的潛在危害。

4 結 論

(1)基于相場格子Boltzmann兩相流動模型,對液橋重開過程的兩相流動進行了數值研究,所得模擬結果與相關實驗及數值研究結果吻合良好,從而驗證了該數值方法的有效性、準確性,體現了該模型在涉及界面拓撲結構變化的復雜兩相流動行為方面的模擬能力。

(2)毛細數對阻塞液橋演化過程的影響存在臨界毛細數現象,本文獲得的臨界毛細數為0.182。毛細數小于臨界值時,阻塞液橋的厚度隨演化過程的進行逐漸增大,氣道仍為阻塞狀態;反之,阻塞液橋的厚度逐漸減小,實現重開。

(3)毛細數高于臨界值時獲得了無液滴和有液滴形成的2種液橋重開現象。在無液滴形成的液橋重開過程中,氣道壁面上的壓力變化更為劇烈,具有更大的潛在危險。

[1] PRISK G K. Invited review: microgravity and the lung [J]. Journal of Applied Physiology, 2000, 89(1): 385-396.

[2] BAUDOIN M, SONG Y, MANNEVILLE P, et al. Airway reopening through catastrophic events in a hierarchical network [J]. Proceedings of the National Academy of Sciences, 2013, 110(3): 859-864.

[3] PRISK G K, GUY H J, ELLIOTT A, et al. Ventilatory inhomogeneity determined from multiple-breath washouts during sustained microgravity on Spacelab SLS-1 [J]. Journal of Applied Physiology, 1995, 78(2): 597-607.

[4] CARPENTER P W, PEDLEY T J. Flow past highly compliant boundaries and in collapsible tubes [C]∥IUTAM Symposium 2001. Neidelberg, Germany: Springer, 2003: 26-30.

[5] LIU H, LUO X, CAI Z, et al. Sensitivity of unsteady collapsible channel flows to modeling assumptions [J]. Communications in Numerical Methods in Engineering, 2009, 25(5): 483-504.

[6] TAI C F, BIAN S, HALPERN D, et al. Numerical study of flow fields in an airway closure model [J]. Journal of Fluid Mechanics, 2011, 677: 483-502.

[7] 許世雄, CHEW Y T, LOW H T, 等. 下呼吸道重開的生物流體力學研究: 理論分析 [J]. 水動力學研究與進展: A輯, 2001, 16(1): 35-43. XU Shixiong, CHEW Y T, LOW H T, et al. Biofluid mechanics investigation on obstructed pulmonary lower airway reopening: theoretical analysis [J]. Journal of Hydrodynamics: A, 2001, 16(1): 35-43.

[8] 許世雄, CHEW Y T, LOW H T, 等. 下呼吸道重開的生物流體力學研究: 實驗模擬 [J]. 生物物理學報, 2000, 16(2): 380-386. XU Shixiong, CHEW Y T, LOW H T, et al. Biofluid mechanics investigation on obstructed pulmonary lower airway reopening: experiment simulation [J]. Acta Biophysica Sinica, 2000, 16(2): 380-386.

[9] MALASHENKO A, TSUDA A, HABER S. Propagation and breakup of liquid menisci and aerosol generation in small airways [J]. Journal of Aerosol Medicine and Pulmonary Drug Delivery, 2009, 22(4): 341-353.

[10]LIU H, KANG Q, LEONARDI C R, et al. Multiphase lattice Boltzmann simulations for porous media applications [J/OL]. Computational Geosciences, 2015: 1-29. [2016-01-15]. doi: 10.1007/S10596-015-9542-3.

[11]王寧寧, 劉海湖, 張楚華. 基于格子Boltzmann方法的Rayleigh-Taylor兩相不穩定流動研究 [J]. 工程熱物理學報, 2016, 37(1): 89-94. WANG Ningning, LIU Haihu, ZHANG Chuhua. Numerical simulation of Rayleigh-Taylor instability based on lattice Boltzmann method [J]. Journal of Engineering Thermophysics, 2016, 37(1): 89-94.

[12]安紅妍, 張楚華, 王寧寧. 基于格子Boltzmann方法的液液不混溶兩相流動數值模擬 [J]. 工程熱物理學報, 2014, 35(1): 78-81. AN Hongyan, ZHANG Chuhua, WANG Ningning. Numerical simulation of liquid-liquid two-phase flow based on lattice Boltzmann method [J]. Journal of Engineering Thermophysics, 2014, 35(1): 78-81.

[13]LIU H, VALOCCHI A J, ZHANG Y, et al. Lattice Boltzmann phase-field modeling of thermocapillary flows in a confined microchannel [J]. Journal of Computational Physics, 2014, 256: 334-356.

[14]QIAN Y, D’HUMIRES D, LALLEMAND P. Lattice BGK models for Navier-Stokes equation [J]. Europhysics Letters, 1992, 17(6): 479.

[15]HU Y, BIAN S, FILOCHE M, et al. Flow and sound generation in human lungs: models of wheezes and crackles [M]. Neidelberg, Germany: Springer, 2014: 301-317.

(編輯 苗凌)

Numerical Simulation for Liquid Bridge Reopening Process with Two-Phase Lattice Boltzmann Method

WANG Ningning,LIU Haihu,ZHANG Chuhua

(School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China)

Based on the recently proposed and developed phase-field lattice Boltzmann method, a computational model is established to simulate the reopening process of blocked human pulmonary airway. The two-phase computational model is established following the free energy theory, and an order parameter is introduced for the description of two-phase interface, which evolves according to the Cahn-Hilliard equation. The model thus has a solid physical foundation. A pressure distribution function is utilized for the hydrodynamic equations, which facilitates minimizing the discretization error of the density gradient to weaken the numerical instability. In addition, an interfacial force of potential form is adopted, which produces much smaller spurious velocities at the interface by comparing with its counterpart of pressure form. The model is used to simulate the liquid bridge reopening process, and the effect of capillary number is analyzed. There exists a critical capillary number, beyond which the axial thickness of the liquid bridge decreases, and finally ruptures, leading to the reopening of blocked human pulmonary airway. Two kinds of reopening processes are reproduced, and they are classified by whether droplets are formed. The reopening process without droplet formed undergoes a more severe pressure jump. This approach is expected be used to further investigate physiology and pathology of human respiratory system and fundamental immiscible two-phase flow phenomenon in microchannels.

lattice Boltzmann method; two-phase flow; microchannel; interfacial force; liquid bridge reopening

2016-03-12。 作者簡介:王寧寧(1989—),女,博士生;劉海湖(通信作者),男,特聘研究員。 基金項目:國家自然科學基金資助項目(51176146,51506168);中組部“千人計劃”青年人才項目;西安交通大學青年拔尖人才支持計劃。

時間:2016-06-14

10.7652/xjtuxb201609009

TK124

A

0253-987X(2016)09-0055-06

網絡出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20160614.1717.006.html