基于能量守恒耗散粒子動力學(xué)方法的自然對流模擬改進(jìn)研究*

魯維 陳碩? 于致遠(yuǎn) 趙嘉毅 張凱旋

1) (同濟(jì)大學(xué)航空航天與力學(xué)學(xué)院力學(xué)系,上海 200092)

2) (上海理工大學(xué)能源與動力工程學(xué)院,上海 200093)

3) (南開大學(xué)醫(yī)學(xué)院,天津 300071)

能量守恒耗散粒子動力學(xué)(eDPD)是一種研究熱輸運(yùn)過程的介觀尺度數(shù)值模擬方法,然而在eDPD 系統(tǒng)內(nèi)引入Boussinesq 假設(shè)以研究自然對流問題時,eDPD 系統(tǒng)自身的熱膨脹性對模擬結(jié)果的影響常常被忽略.首先研究了eDPD 系統(tǒng)的熱膨脹性,通過模擬獲得eDPD 系統(tǒng)的熱膨脹系數(shù)β;并由此模擬了不同瑞利數(shù)Ra、不同幾何結(jié)構(gòu)下的自然對流;利用eDPD 系統(tǒng)自身的熱膨脹性,在不引入Boussinesq 假設(shè)的前提下獲得了合理的溫度場和速度場,與相同Ra 數(shù)下有限體積法模擬結(jié)果相比,誤差明顯小于以往研究中相同條件下的對比誤差.研究表明在eDPD 系統(tǒng)中引入Boussinesq 假設(shè)時,需要考慮eDPD 系統(tǒng)自身熱膨脹性的影響,并且進(jìn)一步對Ra 數(shù)的計(jì)算進(jìn)行了修正.

1 引言

耗散粒子動力學(xué)(dissipative particle dynamics,DPD)是一種介觀尺度流體動力學(xué)數(shù)值模擬方法,由Hoogerbrugge 和Koelman[1]于1992 年提出,隨后Espa?ol 和Warren[2]構(gòu)筑了該方法的理論框架.DPD 被稱為“彌補(bǔ)宏觀尺度和微觀尺度之間空白的方法”[3].DPD 系統(tǒng)中以離散粒子作為動量載體,粒子在離散的時間和空間上運(yùn)動,這些粒子的運(yùn)動狀態(tài)等價于一個指定范圍內(nèi)大量分子的集中行為體現(xiàn),也即DPD 研究介觀尺度上粗粒化粒子的運(yùn)動狀態(tài).

傳統(tǒng)的DPD 模擬的是等溫系統(tǒng),不能模擬介觀尺度上粒子之間熱輸運(yùn)問題,Avalos 和Mackie[4]與Espa?ol[5]在DPD 的基礎(chǔ)上提出了能量守恒耗散粒子動力學(xué)(eDPD),可描述介觀尺度上的粒子間的能量傳遞.Ripoll 等[6,7]通過與一維、二維傳熱問題的理論對比驗(yàn)證了eDPD 模型對于模擬熱輸運(yùn)過程的可行性;Mackie 等[8,9]構(gòu)建了eDPD 參數(shù)與熱傳導(dǎo)系數(shù)以及熱擴(kuò)散率之間的關(guān)系;Lukes[10]通過量綱分析對eDPD 進(jìn)行了模擬參數(shù)與物理參數(shù)之間的映射研究.近年來Homman 等[11]和Stoltz[12]分別研究了eDPD 并行迭代算法以及算法穩(wěn)定性,以謀求更大的時間步長.隨著理論基礎(chǔ)的不斷完善,eDPD 被廣泛應(yīng)用到各個領(lǐng)域.Abu-Nada 完成了一系列基于eDPD 的研究,從固體導(dǎo)熱問題[13,14]到Rayleigh-Bénard(RB)以及封閉空腔中的自然對流問題等[15,16],顯示eDPD 能模擬復(fù)雜自然、受迫熱對流等過程[17,18].但是上述研究中當(dāng)瑞利數(shù)Ra較小時,Abu-Nada 使用eDPD 的模擬結(jié)果不理想,模擬方腔內(nèi)自然對流時,相同條件下的eDPD 模擬結(jié)果與有限體積法(finite volume method,FVM)模擬結(jié)果等溫線對比存在較大誤差.我們分析,誤差的主要原因是以往研究中沒有考慮eDPD 本身熱膨脹性,這部分熱膨脹性在適當(dāng)?shù)臈l件下也會產(chǎn)生相應(yīng)的自然對流,如果忽略其帶來的影響則會在Ra數(shù)量級較小時帶來較大的誤差.

在Abu-Nada 等的研究中,eDPD 系統(tǒng)被視為不可壓縮系統(tǒng),直接引入Boussinesq 假設(shè)來提供模擬中所需的浮升力.然而Mai-Duy 等[19]和Pan等[20]的研究顯示由于DPD 采用粒子間“軟”相互作用,DPD 系統(tǒng)本身具有一定的可壓縮性,在系統(tǒng)溫度較高或保守力足夠大的情況下才可認(rèn)為系統(tǒng)近似不可壓縮.因此通常情況下eDPD 研究自然對流問題時應(yīng)當(dāng)考慮系統(tǒng)自身熱膨脹性的作用,否則會使Ra數(shù)的計(jì)算結(jié)果與實(shí)際的Ra數(shù)出現(xiàn)較大偏差,體現(xiàn)在等溫線對比圖上就會使相同溫度的等溫線無法精確擬合.

本文首先通過模擬獲得eDPD 系統(tǒng)的熱膨脹系數(shù)β,隨后在不引入Boussinesq 假設(shè)的前提下,利用eDPD 系統(tǒng)自身的熱膨脹性,研究了不同Ra數(shù)、不同幾何結(jié)構(gòu)下的自然對流問題,獲得了合理的溫度場和速度場.在相同條件下與FVM 的模擬結(jié)果進(jìn)行對比,誤差要小于Abu-Nada 等[13-18]的模擬誤差.依靠eDPD 系統(tǒng)自身的熱膨脹性可以形成自然對流,表明在引入Boussinesq 假設(shè)時,需要考慮eDPD 系統(tǒng)本身熱膨脹性的影響.本文以Ra為15000 時同心圓環(huán)內(nèi)自然對流的模擬為例,分析eDPD 系統(tǒng)熱膨脹性對不可壓縮假設(shè)下采用Boussinesq 假設(shè)進(jìn)行自然對流模擬所產(chǎn)生的影響.進(jìn)一步,本文對Ra數(shù)的計(jì)算進(jìn)行修正,結(jié)合實(shí)際模擬的Ra數(shù)大小,對系統(tǒng)自身引入Boussinesq 假設(shè)的條件進(jìn)行評估.

2 能量守恒耗散粒子動力學(xué)模型

與分子動力學(xué)類似,eDPD 方法是一種在一定的截?cái)喟霃絩c內(nèi)基于粒子與相鄰粒子之間相互作用的粒子方法,粒子的運(yùn)動符合牛頓第二定律:

這里,vi表示粒子速度,Fi表示粒子受到的力,ri表示粒子的位置分別為保守力、耗散力和隨機(jī)力,分別表示為

式中,aij為粒子之間的保守力系數(shù)與溫度有關(guān),這里我們?nèi)ij=A·kB(Ti+Tj)/2 ,A=18.75,kB是玻爾茲曼常數(shù),Ti和Tj是粒子溫度;ζij是高斯白噪音,同時具有對稱性,即ζij=ζji,以此保證整體的動量守恒;γij和σij是耗散力系數(shù)和隨機(jī)力系數(shù);eij是粒子i和粒子j的單位 向量;vij是粒子i和粒子j之間的相對速度; ?t是時間步長;耗散力ωD權(quán)函數(shù)與保守力權(quán)函數(shù)ωC和隨機(jī)力權(quán)函數(shù)ωR之間的關(guān)系分別為

其中rij=|ri-rj| 是粒子i和粒子j之間的距離.

傳統(tǒng)的DPD 模型里能量不守恒,只能對恒溫的模型進(jìn)行模擬計(jì)算.因此eDPD 在DPD 模型的基礎(chǔ)上引入粒子內(nèi)能變量,使DPD 方程中包含能量項(xiàng),流體粒子之間的熱量傳遞通過可以下式表示[21]:

式 中,CV為粒子的定容比熱容;Ti為粒子溫度;分別為黏性熱通量、碰撞熱通量和隨機(jī)熱通量,分別為

其中,κij和αij分別表示碰撞熱通量系數(shù)和隨機(jī)熱通量系數(shù)為隨機(jī)熱流中的隨機(jī)變量,具有反對稱性,即保證了能量守恒,權(quán)函數(shù)ω隨著粒子之間的距離單調(diào)遞減,超過截?cái)喟霃街缶蜑?,如(6)式和(7)式所示.

耗散力系數(shù)γij和隨機(jī)力系數(shù)σij,它們之間滿足耗散漲落定理.當(dāng)此關(guān)系得以滿足時,系統(tǒng)的溫度將保持在給定值T附近,這一功能也被稱為恒溫器,由(12)–(14)式來實(shí)現(xiàn):

其中kB是玻爾茲曼常數(shù),而ko是用來調(diào)整傳導(dǎo)系數(shù)的參數(shù).在eDPD 模擬中,(7)式的s=2,s也可以取其他值用來調(diào)節(jié)流體的黏性.與DPD 一樣可以通過對于熱能的積分得到eDPD 中各參數(shù)與實(shí)際熱物性的理論公式.其中最重要的就是熱傳導(dǎo)系數(shù)λ和熱擴(kuò)散率αC.根據(jù) Ripoll[22]的研究推導(dǎo)可得

而熱傳導(dǎo)系數(shù)滿足λ=的關(guān)系,從而可得λ的理論值.

3 結(jié)果和討論

3.1 eDPD 系統(tǒng)物性參數(shù)設(shè)置與熱膨脹系數(shù)的計(jì)算

在自然對流問題的模擬中兩個重要的無量綱參數(shù)分別是瑞利數(shù)Ra和普朗特?cái)?shù)Pr,可以通過下式計(jì)算得到

其中g(shù)為重力,β為熱膨脹系數(shù),TH為高溫溫度,TC為低溫溫度,H為計(jì)算域的高度,νC為運(yùn)動黏度,αC為熱擴(kuò)散率,下標(biāo)C 表示該物理量為低溫溫度對應(yīng)的物理量.Ra的大小能判斷流體的自然對流是否能產(chǎn)生,以及自然對流的強(qiáng)烈程度;Pr則反映了流體中能量和動量遷移過程中的相互影響.

在以往研究中eDPD 系統(tǒng)被視為不可壓縮系統(tǒng),Ra中的熱膨脹系數(shù)β是基于Boussinesq 假設(shè)的輸入量.然而如前所述,根據(jù)Mai-Duy[19]和Pan等[20]的研究eDPD 系統(tǒng)不能直接視為不可壓縮系統(tǒng).本文研究了eDPD 系統(tǒng)的熱膨脹系數(shù)β,通過以下方式計(jì)算: 首先構(gòu)建3 個方向?yàn)橹芷谛赃吔鐥l件的正方體計(jì)算域,通過調(diào)整計(jì)算域體積,實(shí)現(xiàn)不同平均溫度下的等壓環(huán)境,等壓條件下正方體計(jì)算域初始密度ρ0=4,粒子質(zhì)量m=1,比熱容CV=27800,保守力系數(shù)a=18.75,隨機(jī)力系數(shù)σ=3,介觀導(dǎo)熱系數(shù)k0=0.000126,系統(tǒng)溫度T變化范圍為1–1.8,所有參數(shù)均為無量綱參數(shù).壓強(qiáng)設(shè)置為35,等壓條件下,密度ρ隨著溫度T變化的規(guī)律如圖1所示,擬合后得到ρ=0.59279T2-2.681T+6.2659 .圖1 也顯示了壓強(qiáng)分別為51 和68 時,密度與溫度之間的關(guān)系,表明壓強(qiáng)對上述eDPD 系統(tǒng)的熱膨脹性影響不大.

圖1 不同壓強(qiáng)下eDPD 密度隨溫度的變化Fig.1.Density variation of eDPD with temperature at different pressure level.

根據(jù)定義,熱膨脹系數(shù)β由下式計(jì)算得到

將ρ=0.59279T2-2.681T+6.2659 代入(18)式,得到eDPD 系統(tǒng)自身熱膨脹系數(shù)β隨溫度T變化的規(guī)律:

β在溫度T=1–1.8 的范圍內(nèi)變化趨勢如圖2所示.

圖2 eDPD 系統(tǒng)熱膨脹系數(shù)隨溫度的變化Fig.2.Variation of thermal expansion coefficient of eDPD system with temperature.

基于泊肅葉流動,在上述系統(tǒng)參數(shù)條件下可以得到系統(tǒng)的動力黏性系數(shù)μ=1.077,對應(yīng)的運(yùn)動黏性系數(shù)ν=μ/ρ=0.265,與相同條件下的理論解[23]一致.

3.2 eDPD 模型驗(yàn)證

基于eDPD 系統(tǒng)自身的熱膨脹性,模擬了方腔中RB 問題.如圖3 所示,方腔模擬區(qū)域?yàn)?30 ≤x≤30,-15 ≤y≤15,-2 ≤z≤2,z方向施加周期性邊界條件.上板為冷板,溫度TC=1.0;下板為熱板,溫度TH=1.8;左右兩板為互相平行的絕熱板,W/H始終保持為2,方腔內(nèi)流體粒子受到向下的重力,g=0.02.將絕熱板的α,σ,k0都設(shè)置為0 以實(shí)現(xiàn)絕熱條件,并施加麥克斯韋反射邊界條件從而防止流體粒子穿透.模擬時間步長為0.02,每5000 步對區(qū)域數(shù)據(jù)進(jìn)行平均統(tǒng)計(jì),總步數(shù)為100000 步.

圖3 eDPD 模擬方腔內(nèi)RB 問題的模型Fig.3.Model of eDPD to simulate RB problem in square cavity.

通過改變溫度和重力的大小以調(diào)整Ra,Pr始終保持為1.分別得到Ra為1600,1800,2000,4000時的自然對流溫度云圖和速度矢量圖,如圖4 和圖5 所示.值得注意的是,模擬中設(shè)置的無量綱重力應(yīng)小于0.1,無量綱溫差也須小于1,否則會使eDPD 模型的邊界產(chǎn)生較大的滑移,并且使系統(tǒng)熱膨脹性降低[19,20],從而導(dǎo)致模擬的速度場、溫度場與實(shí)際結(jié)果產(chǎn)生較大誤差.基于eDPD 模擬中溫度云圖與速度矢量圖隨Ra的變化情況,可以判斷產(chǎn)生自然對流的臨界瑞利數(shù)Ra在1600–1800 之間.當(dāng)Ra小于1600 時流體的傳熱方式以熱傳導(dǎo)為主,Ra大于1800 時流體的傳熱方式以熱對流為主.此時從能量守恒耗散粒子動力學(xué)的控制方程來分析,由于eDPD 粒子之間存在隨機(jī)力,具有隨機(jī)效應(yīng).當(dāng)傳熱方式發(fā)生改變時,從介觀尺度來看系統(tǒng)的自然對流產(chǎn)生并非對稱發(fā)生的,而是隨機(jī)地從某處開始產(chǎn)生微小的環(huán)流.隨著Ra不斷的增大,小的環(huán)流帶動整個流體最終形成穩(wěn)定較大的兩個環(huán)流,這也符合Ra可以反映自然對流的激烈程度的特點(diǎn).

圖4 eDPD 模擬不同瑞利數(shù)下自然對流的溫度云圖 (a) Ra=1600;(b) Ra=1800;(c) Ra=2000;(d) Ra=4000Fig.4.Temperature clouds of natural convection under different Rayleigh numbers simulated by eDPD: (a) Ra=1600;(b) Ra=1800;(c) Ra=2000;(d) Ra=4000.

圖5 eDPD 模擬不同瑞利數(shù)下自然對流的速度矢量圖 (a) Ra=1600;(b) Ra=1800;(c) Ra=2000;(d) Ra=4000Fig.5.Natural convection velocity vectors at different Rayleigh numbers by eDPD simulations: (a) Ra=1600;(b) Ra=1800;(c) Ra=2000;(d) Ra=4000.

本文得到的臨界Ra的范圍與Koschmieder和Pallas[24]實(shí)驗(yàn)得到的臨界瑞利數(shù)Ra=1675 相吻合,同時也符合Abu-Nada[15]得到的臨界瑞利數(shù)Ra=1707 的結(jié)果;此外Zhang 和?nskog[25]通過朗之萬方法基于熱通量變化也得到臨界瑞利數(shù)Ra為1708 的結(jié)論.顯示,利用eDPD 自身的熱膨脹性能夠有效描述自然對流的發(fā)生.

3.3 eDPD 模擬方腔內(nèi)RB 問題

基于eDPD 系統(tǒng)自身的熱膨脹性,對不同Ra條件下方腔內(nèi)自然對流RB 問題進(jìn)行了驗(yàn)證.在3.2 節(jié)參數(shù)設(shè)置基礎(chǔ)上通過調(diào)節(jié)重力加速度,將自然對流模型的Ra降低到3100,Pr始終保持為1,模擬結(jié)果歸一化后得到的溫度云圖如圖6 所示.

圖6 eDPD 模擬Ra 為3100 時方腔內(nèi)自然對流的溫度云圖Fig.6.eDPD simulation of the temperature cloud of natural convection in the square cavity when Ra is 3100.

FVM 針對偏微分方程進(jìn)行數(shù)值求解,它將求解區(qū)域劃分為有限數(shù)量的體積元,并在每個體積元內(nèi)計(jì)算它們與相鄰體積元之間的通量.本文使用FVM 模擬相同的自然對流做對照組,左右兩邊同樣設(shè)置為絕熱邊界條件,上下板溫差為0.8,下板設(shè)置為熱源,計(jì)算域的尺寸、液體的密度、黏性、比熱容、熱擴(kuò)散率、熱膨脹系數(shù)、重力加速度等參數(shù)都與eDPD 系統(tǒng)建模完成之后所得到的參數(shù)一致,Ra和Pr都是無量綱參數(shù),通過(17)式計(jì)算后使FVM 的Ra也為3100,Pr也為1.在此基礎(chǔ)上每隔0.04 個單位提取等溫線,eDPD 的等溫線與相同Ra的FVM 模擬結(jié)果對比如圖7 所示,可以看到eDPD 的模擬結(jié)果與FVM 模擬結(jié)果誤差很小,整體誤差小于0.7%,邊界最大誤差小于1.8%.

圖7 Ra 為3100 時eDPD 和FVM 模擬的方腔內(nèi)自然對流溫度等溫線對比圖Fig.7.Comparison of natural convection temperature isotherms in the square cavity simulated by eDPD and FVM at 3100 of Ra.

繼續(xù)調(diào)整重力大小,其他條件保持不變,將Ra提高到7300,Pr保持為1,如圖8 所示對比關(guān)系,仍獲得較好的對比結(jié)果,整體誤差小于1.1%,邊界處最大誤差小于2.4%.在絕熱邊界處,eDPD的溫度值略微偏離FVM 結(jié)果,主要是由于eDPD中無滑移邊界的不完全滿足所導(dǎo)致的.

圖8 Ra 為7300 時eDPD 和FVM 模擬的方腔內(nèi)自然對流溫度等溫線對比圖Fig.8.Comparison of natural convection temperature isotherms in the square cavity simulated by eDPD and FVM at 7300 of Ra.

在此基礎(chǔ)上再進(jìn)一步對比eDPD 和FVM 的速度場.將速度場沿著X＝0.5 和X＝1.5 取VX沿著Y方向分布,沿著X=0.3 和X=1.0 方向取VY沿著Y方向分布,得到VX,VY速度場度場歸一化之后與FVM 的模擬結(jié)果對比,如圖9 所示.

圖9 Ra 為7300 時eDPD 和FVM 模擬的方腔內(nèi)自然對流速度線對比圖 (a) X=0.5 時VX 沿Y 方向變化曲線;(b) X=1.5 時VX 沿Y 方向變化曲線;(c) X=0.3 時VY 沿Y 方向變化曲線;(d) X=1.0 時VY 沿Y 方向變化曲線Fig.9.Comparison of natural convection velocity profile in the square cavity simulated by eDPD and FVM when Ra is 7300:(a) Variation curve of VX along Y-direction at X=0.5;(b) variation curve of VX along Y-direction at X=1.5;(c) variation curve of VY along the Y-direction at X=0.3;(d) variation curve of VY along the Y-direction at X=1.0.

由圖9 可知,在不同X坐標(biāo)方位上,eDPD 與FVM 的速度場同樣吻合較好,由于無滑移邊界條件不完全滿足,從而導(dǎo)致速度場在邊界處存在一定的偏差.此外,由于eDPD 中的介觀熱波動效應(yīng),統(tǒng)計(jì)的速度存在一定波動,在速度極值處也有一定差異.該現(xiàn)象也存在于以往研究的eDPD 自然對流模擬之中[15,16].

本文中eDPD 模擬方腔內(nèi)自然對流整體誤差為1.1%,相較于Abu-Nada 在相同條件下基于Boussinesq 假設(shè)[15,16]的模擬誤差已有明顯改善.綜上所述,不需要引入Boussinesq 假設(shè),基于系統(tǒng)自身的熱膨脹性的eDPD 可有效模擬不同Ra數(shù)下的自然對流RB 問題.

3.4 eDPD 模擬同心圓環(huán)和偏心圓環(huán)內(nèi)自然對流

通過同心圓環(huán)自然對流問題的模擬,進(jìn)一步驗(yàn)證了基于系統(tǒng)自身的熱膨脹性的eDPD 對復(fù)雜幾何結(jié)構(gòu)中自然對流問題模擬的可行性.圖10 為同心圓環(huán)幾何模型示意圖,內(nèi)圓半徑Ri=10,外圓半徑Ro=30,內(nèi)圓溫度為高溫TH=1.8,外圓溫度為低溫Tc=1.0,其余模擬參數(shù)與3.3 節(jié)中自然對流RB 問題中的設(shè)置相同.同心圓環(huán)的Ra數(shù)和Pr數(shù)由下式得到

圖10 同心圓環(huán)模型示意圖Fig.10.Schematic diagram of the concentric ring model.

調(diào)整重力g的大小和改變Ro大小來獲得不同Ra下的同心圓環(huán)自然對流,Pr始終保持為1,可以得到Ra從1000–10400 時eDPD 模擬結(jié)果與FVM 模擬結(jié)果的等溫線對比圖,如圖11 所示,其中圖11(a)和圖11(b)中的Ro=30,圖11(c)和圖11(d)中Ro=40.

圖11 不同Ra 的 eDPD 和FVM 同心圓環(huán)自然對流等溫線對比圖 (a) Ra=1000;(b) Ra=4100;(c) Ra=7200;(d) Ra=10400 Fig.11.Comparison of natural convection isotherms of eDPD and FVM simulation for concentric rings at different Ra: (a) Ra=1000;(b) Ra=4100;(c) Ra=7200;(d) Ra=10400.

由等溫線對比結(jié)果可知,在同心圓環(huán)自然對流問題的模擬中,eDPD 與FVM 結(jié)果在Ra數(shù)較小的情況下,如圖11(a)和圖11(b)所示,對比結(jié)果誤差較小,整體誤差小于1.4%;當(dāng)Ra提高后,如圖11(c)和圖11(d)所示,溫度在圓環(huán)左右兩側(cè)產(chǎn)生了較明顯的偏差,不過整體誤差仍小于2.7%,低溫邊界處最大誤差小于3.9%,原因在于此處的流體速度較大,邊界條件的影響比較顯著,導(dǎo)致統(tǒng)計(jì)過程中溫度受到較大影響,不過整體模擬結(jié)果的精度相比相同條件下基于Boussinesq 假設(shè)的研究也有了明顯的提高[21,26].同時,通過改變內(nèi)圓位置模擬了偏心圓環(huán)的自然對流,如圖12 所示.類似于同心圓環(huán)模擬,溫度在速度變化較大處仍存在較大偏差,該現(xiàn)象證明為了實(shí)現(xiàn)eDPD 對自然問題的準(zhǔn)確模擬,不僅需要考慮熱膨脹性產(chǎn)生的影響,同時也應(yīng)引入更完善的無滑移邊界條件[27].

圖12 eDPD 和FVM 偏心圓環(huán)自然對流等溫線對比圖 (a) Ra=4100;(b) Ra=10400Fig.12.Comparison of eDPD and FVM eccentric circular natural convection isotherms: (a) Ra=4100;(b) Ra=10400.

上述模擬對比表明,在使用eDPD 進(jìn)行自然對流模擬研究的過程中,系統(tǒng)自身熱膨脹性的作用不應(yīng)該被忽略;考慮熱膨脹性可以減小以往研究中在Boussinesq 假設(shè)下eDPD 與FVM 之間的模擬誤差.最后,本文將對比探討熱膨脹性在傳統(tǒng)基于Boussinesq 假設(shè)的eDPD 模擬中的具體影響.

3.5 eDPD 引入Boussinesq 假設(shè)的自然對流問題模擬

在上述研究的基礎(chǔ)上,進(jìn)一步討論用eDPD模擬自然對流時引入Boussinesq 假設(shè)的適用條件.基于不可壓縮條件,通常引入Boussinesq 假設(shè)后粒子的運(yùn)動方程由下式控制:

其中熱膨脹系數(shù)為βB,重力向量為g,T為流體溫度,To為初始溫度.

由于認(rèn)為流體不可壓縮,流體自身熱膨脹性的影響在上述算法中被忽略.為提高eDPD 模擬自然對流的精度,本文根據(jù)所需模擬的Ra數(shù)量級來判斷是否引入Boussinesq 假設(shè): 當(dāng)按(18)式通過調(diào)整重力、溫差、計(jì)算域大小,結(jié)合eDPD 系統(tǒng)自身的熱膨脹系數(shù),就能得到所需模擬的自然對流的Ra數(shù)時,那么只需要充分考慮系統(tǒng)自身的熱膨脹性即可;該可壓縮條件下的eDPD 自身的熱膨脹性產(chǎn)生的浮升力已經(jīng)可以滿足相應(yīng)自然對流的模擬,3.3 和3.4 節(jié)中方腔內(nèi)RB 問題的模擬和同心圓環(huán)內(nèi)自然對流的模擬均證明模擬精度已足夠準(zhǔn)確;所模擬的自然對流Ra數(shù)量級在103以下時,不需要引入Boussinesq 假設(shè),僅僅依靠eDPD 自身熱膨脹性可完成相應(yīng)自然對流問題的模擬.

當(dāng)所需模擬的Ra數(shù)量級增大時,根據(jù)圖1 可知系統(tǒng)自身所能提供的熱膨脹系數(shù)有限,圖2 也表明溫度的取值范圍有限,無量綱溫差不能超過1,設(shè)置無量綱重力應(yīng)小于0.1,否則邊界會產(chǎn)生較大的滑移.所以由(18)式可知只能擴(kuò)大計(jì)算域的大小來得到所需模擬的自然對流的Ra,但這會大大增加計(jì)算成本,為提高計(jì)算效率、節(jié)約計(jì)算資源,可以在考慮可壓縮條件對eDPD 系統(tǒng)自身的熱膨脹性影響的同時引入Boussinesq 假設(shè),以此來實(shí)現(xiàn)更大數(shù)量級Ra數(shù)的自然對流模擬.

以同心圓環(huán)自然對流為例,圖13 為按照以往研究簡單把eDPD 視為不可壓縮系統(tǒng),直接引入Boussinesq 假設(shè),得到理論瑞利數(shù)Ra為11600 的eDPD 模擬結(jié)果與FVM 模擬結(jié)果的對比圖.由圖13可知,此時eDPD 結(jié)果即使在速度較小的區(qū)域也明顯偏離FVM 的結(jié)果,整體誤差為14.6%.由此可以判斷在該數(shù)量級的Ra下使用eDPD 模擬自然對流時,不考慮eDPD 自身熱膨脹性的影響,誤差會比較明顯,這是因?yàn)閑DPD 系統(tǒng)自身具有一定的熱膨脹性,當(dāng)通過自身熱膨脹性所得到的Ra大于自然對流產(chǎn)生的臨界點(diǎn)時就會有自然對流產(chǎn)生,此時引入Boussinesq 假設(shè)得到的Ra并非是當(dāng)前自然對流對應(yīng)的真實(shí)Ra,需要綜合考慮系統(tǒng)本身和引入Boussinesq 假設(shè)的影響,才能真實(shí)反映當(dāng)前自然對流所對應(yīng)的Ra,這也是直接引入Boussinesq 假設(shè)eDPD 模擬結(jié)果與FVM 誤差較大的原因.

圖13 不考慮eDPD 自身熱膨脹性時同心圓環(huán)自然對流等溫線對比圖Fig.13.Comparison of natural convection isotherms in concentric rings without considering the thermal expansion of the eDPD system.

隨后在引入Boussinesq 假設(shè)的基礎(chǔ)上,考慮系統(tǒng)熱膨脹性對模擬結(jié)果的影響.由粒子運(yùn)動方程(21)和浮升力公式

可知,Boussinesq 假設(shè)本質(zhì)上給eDPD 粒子施加一個體力以提供實(shí)現(xiàn)自然對流所需的浮升力,在eDPD 中這個力的施加方式和重力一致.要保持良好的邊界條件,浮升力取值范圍也應(yīng)該和重力一樣小于0.1,同時如要滿足無量綱溫差不能超過1,那么可以得到Boussinesq 假設(shè)中的βB取值范圍為

由此可知Boussinesq 假設(shè)極限條件下可以提供3–10 倍于系統(tǒng)自身的熱膨脹系數(shù).

基于上述分析可知,Boussinesq 假設(shè)產(chǎn)生的作用與eDPD 自身熱膨脹性的作用產(chǎn)生了疊加影響.通過計(jì)算可知,本節(jié)中eDPD 構(gòu)建的同心圓環(huán)模型自身熱膨脹性對應(yīng)的瑞利數(shù)Ra1為3400,在此基礎(chǔ)上按照以往研究中引入Boussinesq 假設(shè)對應(yīng)的瑞利數(shù)Ra2為11600,最終Ra總和應(yīng)為15000.在相同條件下考慮熱膨脹性的eDPD 模擬結(jié)果,與Ra為15000 的FVM 模擬結(jié)果對比如圖14 所示.結(jié)果顯示誤差顯著減小,整體誤差降低至3.2%.上述結(jié)果證明eDPD 此時模擬的自然對流所對應(yīng)的Ra更接近15000 而不是11600,說明模擬該數(shù)量級Ra的自然對流時,由系統(tǒng)自身熱膨脹性產(chǎn)生的浮升力在實(shí)際計(jì)算中應(yīng)該被充分考慮,否則模擬所得的自然對流所對應(yīng)的Ra與預(yù)設(shè)需要得到的Ra會有較大偏差,從而導(dǎo)致模擬準(zhǔn)確度受到較大影響.由結(jié)果還可以繼續(xù)推斷出當(dāng)采用eDPD 繼續(xù)模擬更大數(shù)量級Ra的自然對流時,可在考慮自身熱膨脹性的條件下,通過結(jié)合Boussinesq 假設(shè)實(shí)現(xiàn)大Ra的模擬;并且系統(tǒng)熱膨脹性對整體熱膨脹系數(shù)影響會隨著所需模擬的Ra增大而減小;這種情況下引入的Boussinesq 假設(shè)可以達(dá)到提高計(jì)算效率的目的.綜上所述,應(yīng)當(dāng)考慮eDPD 自身的熱膨脹性對整體自然對流的影響.

圖14 考慮eDPD 自身熱膨脹性的同時引入Boussinesq假設(shè)情況下,同心圓環(huán)自然對流等溫線對比圖Fig.14.Concentric circular natural convection isotherm comparison when considering the combined effect of thermal expansion of eDPD system and Boussinesq assuming.

4 結(jié)論

本文基于eDPD 系統(tǒng)的熱膨脹性,首先獲得了eDPD 系統(tǒng)自身的熱膨脹系數(shù)β,在此基礎(chǔ)上模擬了不同Ra數(shù)下,方腔內(nèi)RB 問題、同心圓環(huán)和偏心圓環(huán)的自然對流問題.通過對Ra從1600–4000之間的方腔內(nèi)RB 問題溫度場和速度場的變化進(jìn)行分析,顯示流體產(chǎn)生自然對流的臨界瑞利數(shù)Ra在1600–1800 之間,該結(jié)果與現(xiàn)有研究結(jié)果一致;之后定量分析了在eDPD 系統(tǒng)中Ra數(shù)為3100 和7300 時方腔內(nèi)RB 問題的溫度場和速度場,并與FVM 的模擬結(jié)果對比,最大誤差小于2.4%.

隨后模擬同心圓環(huán)自然對流問題,將Ra從1000–10400 之間同心圓環(huán)自然對流的溫度場與FVM 的模擬結(jié)果進(jìn)行對比,最大誤差小于3.9%,在小Ra數(shù)下,本文通過eDPD 自身熱膨脹性得到的自然對流誤差要遠(yuǎn)小于以往研究中的誤差.通過分析偏心圓環(huán)的自然對流問題模擬結(jié)果,顯示由于eDPD 中的介觀熱波動效應(yīng),統(tǒng)計(jì)的速度存在一定波動,加上邊界不能完全滿足無滑移條件,導(dǎo)致模擬結(jié)果在邊界處和極值處存在一定的偏差,未來需要進(jìn)一步完善.

最后本文探討了在eDPD 系統(tǒng)中自身熱膨脹性對于引入Boussinesq 假設(shè)的影響,并對相應(yīng)的Ra數(shù)進(jìn)行修正.按照以往的研究,以同心圓環(huán)內(nèi)自然對流模擬為例,eDPD 直接引入Boussinesq 假設(shè)得到的瑞利數(shù)Ra為11600,與Ra為11600 的FVM模擬所得溫度場對比,發(fā)現(xiàn)誤差為14.6%;而考慮自身熱膨脹性后引入Boussinesq 假設(shè),相同條件下的同心圓環(huán)自身熱膨脹性所對應(yīng)的Ra1為3400,引入Boussinesq 假設(shè)對應(yīng)的Ra2為11600,兩者影響相互疊加,整體的Ra達(dá)到15000,eDPD 模擬結(jié)果與Ra為15000 的FVM 溫度場模擬結(jié)果擬合誤差大大減小,誤差降低至3.2%.

進(jìn)一步,本文根據(jù)所需模擬的Ra數(shù)量級來判斷是否引入Boussinesq 假設(shè): 當(dāng)所需模擬的自然對流的Ra依靠eDPD 系統(tǒng)自身的熱膨脹系數(shù)就能達(dá)到,那么只需要充分考慮系統(tǒng)熱膨脹性導(dǎo)致的浮升力帶來的影響即可,此時模擬精度已足夠準(zhǔn)確,不需要引入Boussinesq 假設(shè);當(dāng)所需模擬的Ra數(shù)量級增大時,由于系統(tǒng)自身所能提供的熱膨脹系數(shù)有限,只能擴(kuò)大計(jì)算域的大小,但這會增加計(jì)算成本,可以在考慮系統(tǒng)自身熱膨脹性影響的同時,引入Boussinesq 假設(shè)實(shí)現(xiàn)較大Ra的自然對流模擬;隨著所需模擬的Ra增大,eDPD 自身的熱膨脹性對整體自然對流的影響會減弱,這種情況下引入的Boussinesq 假設(shè)可以達(dá)到提高計(jì)算效率的目的.綜上所述,應(yīng)當(dāng)充分考慮eDPD 自身的熱膨脹性對整體自然對流的影響,以提高自然對流模擬的有效性.本文的工作可為eDPD 方法研究自然對流問題提供新的見解.