魯維 陳碩? 于致遠 趙嘉毅 張凱旋

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

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

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

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

1 引言

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

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

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

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

2 能量守恒耗散粒子動力學模型

與分子動力學類似,eDPD 方法是一種在一定的截斷半徑rc內基于粒子與相鄰粒子之間相互作用的粒子方法,粒子的運動符合牛頓第二定律:

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

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

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

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

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

其中,κij和αij分別表示碰撞熱通量系數(shù)和隨機熱通量系數(shù)為隨機熱流中的隨機變量,具有反對稱性,即保證了能量守恒,權函數(shù)ω隨著粒子之間的距離單調遞減,超過截斷半徑之后就為0,如(6)式和(7)式所示.

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

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

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

3 結果和討論

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

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

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

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

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

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

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

β在溫度T=1–1.8 的范圍內變化趨勢如圖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,對應的運動黏性系數(shù)ν=μ/ρ=0.265,與相同條件下的理論解[23]一致.

3.2 eDPD 模型驗證

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

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

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

圖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ù)Ra=1675 相吻合,同時也符合Abu-Nada[15]得到的臨界瑞利數(shù)Ra=1707 的結果;此外Zhang 和?nskog[25]通過朗之萬方法基于熱通量變化也得到臨界瑞利數(shù)Ra為1708 的結論.顯示,利用eDPD 自身的熱膨脹性能夠有效描述自然對流的發(fā)生.

3.3 eDPD 模擬方腔內RB 問題

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

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

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

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

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

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

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

圖9 Ra 為7300 時eDPD 和FVM 模擬的方腔內自然對流速度線對比圖 (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坐標方位上,eDPD 與FVM 的速度場同樣吻合較好,由于無滑移邊界條件不完全滿足,從而導致速度場在邊界處存在一定的偏差.此外,由于eDPD 中的介觀熱波動效應,統(tǒng)計的速度存在一定波動,在速度極值處也有一定差異.該現(xiàn)象也存在于以往研究的eDPD 自然對流模擬之中[15,16].

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

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

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

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

調整重力g的大小和改變Ro大小來獲得不同Ra下的同心圓環(huán)自然對流,Pr始終保持為1,可以得到Ra從1000–10400 時eDPD 模擬結果與FVM 模擬結果的等溫線對比圖,如圖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.

由等溫線對比結果可知,在同心圓環(huán)自然對流問題的模擬中,eDPD 與FVM 結果在Ra數(shù)較小的情況下,如圖11(a)和圖11(b)所示,對比結果誤差較小,整體誤差小于1.4%;當Ra提高后,如圖11(c)和圖11(d)所示,溫度在圓環(huán)左右兩側產(chǎn)生了較明顯的偏差,不過整體誤差仍小于2.7%,低溫邊界處最大誤差小于3.9%,原因在于此處的流體速度較大,邊界條件的影響比較顯著,導致統(tǒng)計過程中溫度受到較大影響,不過整體模擬結果的精度相比相同條件下基于Boussinesq 假設的研究也有了明顯的提高[21,26].同時,通過改變內圓位置模擬了偏心圓環(huán)的自然對流,如圖12 所示.類似于同心圓環(huán)模擬,溫度在速度變化較大處仍存在較大偏差,該現(xiàn)象證明為了實現(xiàn)eDPD 對自然問題的準確模擬,不僅需要考慮熱膨脹性產(chǎn)生的影響,同時也應引入更完善的無滑移邊界條件[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 進行自然對流模擬研究的過程中,系統(tǒng)自身熱膨脹性的作用不應該被忽略;考慮熱膨脹性可以減小以往研究中在Boussinesq 假設下eDPD 與FVM 之間的模擬誤差.最后,本文將對比探討熱膨脹性在傳統(tǒng)基于Boussinesq 假設的eDPD 模擬中的具體影響.

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

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

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

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

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

以同心圓環(huán)自然對流為例,圖13 為按照以往研究簡單把eDPD 視為不可壓縮系統(tǒng),直接引入Boussinesq 假設,得到理論瑞利數(shù)Ra為11600 的eDPD 模擬結果與FVM 模擬結果的對比圖.由圖13可知,此時eDPD 結果即使在速度較小的區(qū)域也明顯偏離FVM 的結果,整體誤差為14.6%.由此可以判斷在該數(shù)量級的Ra下使用eDPD 模擬自然對流時,不考慮eDPD 自身熱膨脹性的影響,誤差會比較明顯,這是因為eDPD 系統(tǒng)自身具有一定的熱膨脹性,當通過自身熱膨脹性所得到的Ra大于自然對流產(chǎn)生的臨界點時就會有自然對流產(chǎn)生,此時引入Boussinesq 假設得到的Ra并非是當前自然對流對應的真實Ra,需要綜合考慮系統(tǒng)本身和引入Boussinesq 假設的影響,才能真實反映當前自然對流所對應的Ra,這也是直接引入Boussinesq 假設eDPD 模擬結果與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 假設的基礎上,考慮系統(tǒng)熱膨脹性對模擬結果的影響.由粒子運動方程(21)和浮升力公式

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

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

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

圖14 考慮eDPD 自身熱膨脹性的同時引入Boussinesq假設情況下,同心圓環(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 結論

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

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

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

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