楊晨，何航行

（重慶大學(xué)動(dòng)力工程學(xué)院低品位能源利用技術(shù)及系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室，重慶400044）

引 言

在現(xiàn)實(shí)問(wèn)題中，流體或多或少總會(huì)含有其他物質(zhì)，絕大多數(shù)流動(dòng)屬于多相流。氣液兩相流動(dòng)現(xiàn)象廣泛存在于能源、動(dòng)力、石油、航天等多個(gè)技術(shù)領(lǐng)域，如核反應(yīng)堆內(nèi)的流動(dòng)沸騰現(xiàn)象等。過(guò)去幾十年中，各國(guó)學(xué)者對(duì)于一些典型的兩相流動(dòng)傳熱過(guò)程進(jìn)行了大量的數(shù)值模擬研究。傳統(tǒng)的數(shù)值方法很難模擬相界面動(dòng)力學(xué)行為，無(wú)法刻畫出兩相流動(dòng)的細(xì) 節(jié)[1-3]。同時(shí)，隨著研究的深入，僅著眼于宏觀尺度的模擬研究已經(jīng)無(wú)法滿足研究和發(fā)展的需求[4]。目前，基于微細(xì)尺度（微觀尺度、介觀尺度等）的模擬方法如動(dòng)力學(xué)Monte Carlo 方法（KMC）、分子動(dòng)力學(xué)方法（MD）和格子Boltzmann 方法（LBM）等是深入研究這些復(fù)雜過(guò)程的有效手段。直接運(yùn)用這些方法在微細(xì)尺度上能夠建立反映過(guò)程內(nèi)在機(jī)理的微觀模型和介觀模型，但其模擬過(guò)程需要完全依賴計(jì)算量巨大的微細(xì)仿真器（fine scale simulator，F(xiàn)SS），經(jīng)常面臨計(jì)算機(jī)內(nèi)存開銷巨大、超出當(dāng)前 計(jì)算能力的問(wèn)題，以致無(wú)法實(shí)現(xiàn)對(duì)某些復(fù)雜實(shí)際 問(wèn)題的模擬，導(dǎo)致微細(xì)尺度模擬研究面臨極大的 困難[5]。

近年來(lái)，有學(xué)者提出了一種計(jì)算機(jī)輔助分析復(fù)雜系統(tǒng)的新型多尺度模擬框架——免方程方法(equation-free method，EFM)[6-7]。其關(guān)鍵思路是繞過(guò)從微細(xì)尺度模型中獲得宏觀模型方程這一傳統(tǒng)過(guò)程，為復(fù)雜系統(tǒng)多尺度模擬提供了一種新的思路。目前免方程方法的研究在國(guó)外發(fā)展迅速，但國(guó)內(nèi)只有少數(shù)文獻(xiàn)提及。早在20世紀(jì)末，Theodoropoulos等[6]已涉及到粗粒化時(shí)間步進(jìn)（coarse time-stepper，CTS）這種跨尺度耦合思想；2001年Gear 等[7-9]開始從事粗?；队胺e分（coarse projective integration，CPI）的研究；之后Kevrekidis 和Gear 等[10-13]運(yùn)用免方程方法進(jìn)行了一系列研究。由此，該多尺度模擬框架逐漸受到復(fù)雜系統(tǒng)模擬研究者的廣泛關(guān)注，其應(yīng)用研究也逐漸增多。

LBM 方法由于具有微觀粒子特性，可以方便地描述不同相之間的相互作用，清晰地捕捉氣液相界面[14]。因此，本研究充分運(yùn)用這一優(yōu)勢(shì)，提出了一種用于模擬氣液兩相分離過(guò)程的免方程多尺度方法。利用該方法模擬水在van der Waals（vdW）、Redlich-Kwong（R-K）和Peng-Robinson（P-R）狀態(tài)方程控制下的相變過(guò)程，發(fā)現(xiàn)多尺度模擬能夠從微細(xì)底層很好地模擬出氣液兩相分離過(guò)程中密度的分布，所得結(jié)果與LBM 單尺度模擬十分吻合；多尺度模擬計(jì)算產(chǎn)生的偽速度與LBM 單尺度模擬結(jié)果基本相同，不會(huì)導(dǎo)致偽速度增加。這說(shuō)明該方法在提高了計(jì)算效率的同時(shí)保證了計(jì)算的穩(wěn)定性，驗(yàn)證了該多尺度模擬方法的準(zhǔn)確性和高效性。通過(guò)免方程方法在氣液兩相分離過(guò)程中的應(yīng)用研究，可為進(jìn)一步從細(xì)微尺度出發(fā)對(duì)流動(dòng)沸騰中氣泡行為等復(fù)雜多相流過(guò)程的多尺度模擬奠定基礎(chǔ)。

1 免方程多尺度方法

1.1 基本框架

免方程方法（EFM）建立了一種基于閉包（closure）的系統(tǒng)辨識(shí)，同時(shí)構(gòu)建了微觀/隨機(jī)仿真和傳統(tǒng)連續(xù)科學(xué)計(jì)算、數(shù)值分析的橋梁。目前免方程方法中常見的應(yīng)用是時(shí)間尺度上的多尺度模擬，因此一般以時(shí)間尺度上的多尺度模擬為其基本框架。該基本框架由粗粒化時(shí)間步進(jìn)（CTS）和相應(yīng)的數(shù)值技術(shù)構(gòu)成，其中CTS 相當(dāng)于整個(gè)框架中的“內(nèi)核”。在EFM 中，粗?；瘯r(shí)間步進(jìn)包括“提升”（lifting，借助于適當(dāng)?shù)姆植己瘮?shù)從宏觀集中獲得分子狀態(tài)）、演化和“約束”（restriction，通過(guò)平均、篩選、光滑處理從微觀或分子模型獲得宏觀狀態(tài)集）3 個(gè)過(guò)程[6]。其中，提升過(guò)程是該組合中首先需實(shí)現(xiàn)的過(guò)程，它是一種跨尺度過(guò)程，能夠把宏觀系統(tǒng)量離散轉(zhuǎn)化成與之相匹配的細(xì)小量離散分布，如在氣體動(dòng)態(tài)模擬中提升過(guò)程把某氣體的密度離散成大量氣體粒子的隨機(jī)離散分布。演化過(guò)程把提升過(guò)程得到的細(xì)小量離散分布投入到細(xì)小尺度模型中進(jìn)行運(yùn)算，經(jīng)過(guò)一些短的演化時(shí)間后從中得到新的細(xì)小量離散分布。在CTS 中演化不需要尺度跨越，整個(gè)過(guò)程都在微細(xì)尺度上完成。最后執(zhí)行約束過(guò)程，其目的恰好與提升過(guò)程相反，即把新得到的細(xì)小量離散分布反饋到宏觀水平。接下來(lái)把從連續(xù)數(shù)次CTS中得到的一系列宏觀系統(tǒng)值抽離出來(lái)，運(yùn)用合理的數(shù)值技術(shù)如Euler 法、Runge-Kutta 法、Newton- Raphson 法等對(duì)它們進(jìn)行相應(yīng)的數(shù)值處理。重復(fù)進(jìn)行上述3 個(gè)過(guò)程，達(dá)到對(duì)系統(tǒng)分析研究的目的。其基本框架如圖1所示。需要強(qiáng)調(diào)的是，免方程方法并不是一種具體的模擬方法，而是一種多尺度模擬框架。

圖1 基于EFM 的多尺度框架Fig.1 EFM-based multi-scale framework

1.2 粗?；队胺e分

粗?；队胺e分（CPI）是免方程方法中一種具有代表性的方法。其基本思想是在復(fù)雜系統(tǒng)的微細(xì)仿真器上運(yùn)用傳統(tǒng)的數(shù)值外推方法，達(dá)到對(duì)復(fù)雜系統(tǒng)進(jìn)行分析研究的目的，即將粗?；瘯r(shí)間步進(jìn)獲取的宏觀量實(shí)施外推投影以預(yù)測(cè)和捕獲系統(tǒng)的特 性[15-17]。

這里簡(jiǎn)要說(shuō)明CPI 的基本原理（圖2）：假設(shè)模擬過(guò)程中的某一小段總共有n個(gè)連續(xù)演化步，當(dāng)?shù)? 步的初始宏觀物理值提升至細(xì)微仿真器后，必須通過(guò)解析相對(duì)困難的細(xì)微仿真器依次計(jì)算這n個(gè)演化步，才能得到第n步的宏觀物理值。這樣雖然能計(jì)算出第n步的宏觀值，但其運(yùn)算效率低。而采用CPI 進(jìn)行多尺度計(jì)算，則只需通過(guò)該細(xì)微仿真器計(jì)算出前k個(gè)演化步的宏觀物理值，然后對(duì)這k個(gè)已獲得的宏觀物理值進(jìn)行有效的外推處理，就能快速準(zhǔn)確地預(yù)測(cè)出第n步的宏觀值。由于中間省去了h（即n k- ）個(gè)演化步的計(jì)算量以及運(yùn)用了較為準(zhǔn)確的外推處理，CPI 能在保證系統(tǒng)模擬精度的同時(shí)有效地降低內(nèi)存消耗，提高CPU 的計(jì)算效率。CPI 模擬全過(guò)程如圖3所示。

圖2 某一小段CPI 模擬過(guò)程Fig.2 CPI simulation diagram of a short part

圖3 CPI 模擬全過(guò)程Fig.3 Entire CPI simulation process diagram

由于CPI 無(wú)須獲取系統(tǒng)的宏觀模型方程，模擬研究是在現(xiàn)象及過(guò)程的細(xì)微仿真器上直接運(yùn)用傳統(tǒng)外推數(shù)值方法快速準(zhǔn)確地預(yù)測(cè)系統(tǒng)的宏觀性能，著 眼于過(guò)程細(xì)微尺度模型的多尺度模擬方法能更真實(shí)地反映實(shí)際過(guò)程的內(nèi)在特性。本文應(yīng)用伸縮式投影方法（TPM），選用具有二階精度的Adams-Bashforth預(yù)測(cè)校正數(shù)值方法進(jìn)行外推處理，其格式可表示為

式中，F(xiàn)為變量的預(yù)報(bào)值，F(xiàn)為變量的校正值，k為演化數(shù)，h為投影數(shù)，α、β、χ、φ為Adams- Bashforth 方法中的二級(jí)系數(shù),以下仿真計(jì)算中α取0.4，β取0.6，χ取0.3，φ取0.7。

2 基本模型

2.1 格子Boltzmann 模型

單組分多相LB 方程為[18-19]

式中，x、t分別表示空間及時(shí)間坐標(biāo)；i=0,1,2,…,b表示每個(gè)格點(diǎn)的不同離散速度方向；ci表示各個(gè)離散方向上的速度矢量；Δt表示格子離散的時(shí)間步長(zhǎng)；τ為量綱1 松弛時(shí)間；為平衡態(tài)粒子分布函數(shù)。

由于D2Q9 模型比D2Q7 模型的數(shù)值穩(wěn)定性好，本研究選取D2Q9 模型，因此Maxwell 平衡態(tài)分布函數(shù)為

式中，權(quán)重系數(shù)ω0=4/9，ωi=1/9（i=1,2,3,4），1/36（i=5,6,7,8）；c=Δx/Δt為粒子遷移速度；equ為平衡態(tài)速度。密度和動(dòng)量的宏觀參數(shù)分別為

根據(jù)文獻(xiàn)，要使氣液兩相分離，必須考慮相間相互作用力。Shan 和Chen 提出了一種能夠直接刻畫粒子間相互作用的LB 模型，即利用一個(gè)偽勢(shì)函數(shù)來(lái)反映不同相之間的作用力

式中，?k(x)、分別是x處第k相和x′處第相的有效密度。若只考慮臨近格點(diǎn)間的相互作用力，格林函數(shù)為

參數(shù)的絕對(duì)值決定了相間相互作用的強(qiáng)度，其符號(hào)決定二者之間是相互吸引還是相互排斥，Δx為格子長(zhǎng)度。在由Shan 和Chen 提出的單組分多相LB模型中，位置x和位置x′的粒子之間的相互作用力為[13]

式中，c0為常數(shù)，對(duì)于D2Q9 格子模型c0=6.0。

Shan-Chen 偽勢(shì)模型通過(guò)平衡態(tài)速度來(lái)體現(xiàn)相間作用力的影響，新的平衡速度為

2.2 提升算子和約束算子

結(jié)合LB 模型氣液兩相流動(dòng)過(guò)程的演變規(guī)律，經(jīng)過(guò)多次試驗(yàn)分析得到提升算子的具體形式

式中，1/48＜λ＜1/20。

約束算子則選取為L(zhǎng)B 模型密度的宏觀參數(shù)形式

2.3 狀態(tài)方程

粒子間相互作用力由式(8)給定，系統(tǒng)的狀態(tài)方程可表示為[20]

式中，cs為L(zhǎng)B 模型的格子聲速，對(duì)于D2Q9模型cs=1/。

本研究選取了以下3 種狀態(tài)方程（equation of state，EOS）：van der Waals（vdW）狀態(tài)方程，Redlich- Kwong（R-K）狀態(tài)方程以及Peng-Robinson（P-R）狀態(tài)方程。

vdW 狀態(tài)方程

式中，R為通用氣體常數(shù)，格子Boltzmann 方法中取值為1.0。壓力對(duì)密度求一階導(dǎo)數(shù)和二階導(dǎo)數(shù)并令其等于零，則可得到臨界點(diǎn)的數(shù)值：ρc=1/(3b)，Tc=8a/(27bR)。a和b為常量，它們分別取值a=9/49、b=2/21。

R-K 狀態(tài)方程

式中，a和b的值由臨界參數(shù)表達(dá)：a=

P-R 狀態(tài)方程

式中，a和b的值由臨界參數(shù)表達(dá)：a=為偏心因子，與工質(zhì)本身有關(guān)，水的偏心因子ω=0.344。

3 仿真結(jié)果與分析

本研究數(shù)值計(jì)算區(qū)域的格子劃分為256×256，平衡態(tài)分布函數(shù)采用D2Q9 模型，松弛因子 1τ=，上下邊界條件和左右邊界條件均為周期性邊界條件，偽勢(shì)密度函數(shù)為任意常數(shù)。初始密度值ρ=ρc，為了使方程演化，在初始密度附近取0.01 的擾動(dòng)。模擬結(jié)果均在以內(nèi)存4GB、CPU 主頻3.1GHz 的四核Xeon 處理器為計(jì)算平臺(tái)的同一臺(tái)計(jì)算機(jī)上運(yùn)行得到。

3.1 不同狀態(tài)方程的結(jié)果對(duì)比

圖4給出了不同狀態(tài)方程控制下LB 模型模擬所得的密度曲線與實(shí)驗(yàn)值的對(duì)比，圖中的實(shí)驗(yàn)值根據(jù)水蒸氣性質(zhì)圖表繪制得出，離散點(diǎn)為數(shù)值模擬計(jì)算得出的結(jié)果。通過(guò)對(duì)比可以發(fā)現(xiàn)，vdW 狀態(tài)方程的溫度取值范圍最小，T/Tc取值范圍約為[0.86，1.0]；P-R 狀態(tài)方程的溫度取值范圍最大，約為[0.76，1.0]；R-K 狀態(tài)方程中T/Tc取值范圍約為[0.8，1.0]。在可取的溫度范圍內(nèi)，vdW 狀態(tài)方程控制下的模擬結(jié)果偏離實(shí)驗(yàn)值較多，密度相對(duì)誤差較大，氣相密度比實(shí)驗(yàn)值大，隨著溫度升高模擬結(jié)果與實(shí)驗(yàn)值的相對(duì)誤差逐漸減小，而液相密度模擬結(jié)果比實(shí)驗(yàn)值小，與實(shí)驗(yàn)值的相對(duì)誤差基本維持在13%左右。R-K 和P-R 狀態(tài)方程模擬所得曲線與實(shí)驗(yàn)值吻合較好，對(duì)于氣相密度R-K狀態(tài)方程得到的結(jié)果與實(shí)驗(yàn)值最接近[圖4(b)]，對(duì)于液相密度P-R 狀態(tài)方程得到的結(jié)果較好。

圖5和圖6分別給出了vdW 和R-K、R-K 和P-R 狀態(tài)方程控制下最大偽速度隨密度比的變化曲線。通過(guò)對(duì)比可以看出，隨著密度比的增大，各模型的最大偽速度均會(huì)增大。對(duì)于某一給定的密度比，R-K 模型的最大偽速度明顯小于vdW 模型，而P-R 模型的最大偽速度小于R-K 模型；同時(shí)，P-R 模型密度比的取值范圍最大，可以處理較高密度比的問(wèn)題。綜上所述，P-R 狀態(tài)方程能較好地模擬水的相變過(guò)程。

圖4 不同狀態(tài)方程控制下的飽和密度曲線與實(shí)驗(yàn)值的對(duì)比Fig.4 Comparison of coexistence curves obtained from simulations with experimental values for different EOS

圖5 vdW 和R-K 狀態(tài)方程最大偽速度的變化曲線Fig.5 Maximum magnitude of spurious current changes with density ratio for vdW and R-K EOS

圖6 R-K 和P-R 狀態(tài)方程最大偽速度的變化曲線Fig.6 Maximum magnitude of spurious current changes with density ratio for R-K and P-R EOS

3.2 LBM 與TPM 模擬的對(duì)比

圖7～圖12給出了LBM 單尺度模擬與TPM 多尺度模擬在不同狀態(tài)方程控制下所得結(jié)果的對(duì)比及其與實(shí)驗(yàn)值的對(duì)比?？傆?jì)算步數(shù)為20000，此部分TPM 多尺度模擬選取的演化數(shù)和投影數(shù)均為400。

圖7 vdW 模型的飽和密度曲線對(duì)比Fig.7 Comparison of coexistence curves for vdW EOS

圖8 vdW 模型最大偽速度對(duì)比Fig.8 Comparison of spurious current for vdW EOS

圖9 R-K 模型的飽和密度曲線對(duì)比Fig.9 Comparison of coexistence curves for R-K EOS

圖10 R-K 模型最大偽速度對(duì)比Fig.10 Comparison of spurious current for R-K EOS

圖11 P-R 模型的飽和密度曲線對(duì)比Fig.11 Comparison of coexistence curves for P-R EOS

圖12 P-R 模型最大偽速度對(duì)比Fig.12 Comparison of spurious current for P-R EOS

由于vdW 狀態(tài)方程本身形式簡(jiǎn)單，從圖7可以看出vdW 模型的飽和密度曲線偏離實(shí)驗(yàn)值較大，但 其變化趨勢(shì)與實(shí)驗(yàn)值保持一致。TPM 多尺度模擬結(jié)果與LBM 單尺度模擬結(jié)果十分接近，而且TPM 多尺度模擬結(jié)果比LBM 模擬結(jié)果更接近實(shí)驗(yàn)值。

圖8給出了vdW模型最大偽速度隨密度比的變化情況?？梢钥闯?，TPM 模擬產(chǎn)生的偽速度與LBM模擬產(chǎn)生的偽速度基本相同，TPM 模擬產(chǎn)生的偽速度除了最后一個(gè)點(diǎn)之外都比LBM 模擬小，說(shuō)明采用此種情況下的TPM 多尺度模擬不會(huì)導(dǎo)致模型偽速度增大，保證了模型計(jì)算的穩(wěn)定性。

圖9給出了R-K 模型經(jīng)過(guò)20000 步迭代后的飽和密度曲線?？梢钥闯觯琓PM 多尺度模擬的結(jié)果與LBM 單尺度模擬的結(jié)果在氣相部分吻合良好，而且與實(shí)驗(yàn)值十分接近；液相部分TPM 多尺度模擬結(jié)果比LBM 單尺度模擬結(jié)果更接近實(shí)驗(yàn)值。模擬結(jié)果比vdW 模型更加準(zhǔn)確，模擬計(jì)算出的密度值準(zhǔn)確性較好。在低溫處，液相密度的相對(duì)誤差小于2%；氣相密度絕對(duì)誤差很小，但相對(duì)誤差較大。這是由狀態(tài)方程的性質(zhì)造成的，因?yàn)楫?dāng)溫度遠(yuǎn)低于臨界溫度時(shí)密度比會(huì)變得很大。隨著溫度逐漸接近飽和溫度，氣相密度的相對(duì)誤差從27%減小到7%。

圖10給出了R-K 模型最大偽速度隨密度比的變化趨勢(shì)?？梢钥闯?，TPM 多尺度模擬結(jié)果與LBM單尺度模擬結(jié)果基本保持一致，TPM 產(chǎn)生的偽速度平均值小于LBM 單尺度模擬產(chǎn)生的偽速度平均值，說(shuō)明采用此種情況下的TPM 多尺度模擬降低了模型的偽速度，提高了模型的穩(wěn)定性。

圖11給出了P-R 模型經(jīng)過(guò)20000 步迭代后的飽和密度曲線?？梢钥闯?，TPM 多尺度模擬結(jié)果與LBM 單尺度模擬結(jié)果十分吻合，氣液兩相的密度變化趨勢(shì)與實(shí)驗(yàn)值較為接近，但P-R 模型產(chǎn)生的誤差比R-K 模型產(chǎn)生的誤差大。在低溫處，液相密度的相對(duì)誤差小于10%；但氣相密度相對(duì)誤差還是較大，盡管絕對(duì)誤差很小。

圖12給出了P-R 模型最大偽速度隨密度比的變化趨勢(shì)。可以看出，P-R 模型的密度比變化范圍比R-K 模型大，能處理高密度比情況下的問(wèn)題。TPM 多尺度模擬的偽速度小于LBM 單尺度模擬的偽速度，說(shuō)明此種情況下采用TPM 多尺度模擬不僅沒(méi)有導(dǎo)致模型的偽速度增大，反而降低了模型的偽速度，提高了模型的穩(wěn)定性，充分體現(xiàn)了該多尺度模擬方法的有效性。

3.3 不同投影步數(shù)的模擬結(jié)果對(duì)比

改變投影和演化的步數(shù)，分別選取演化400 步、投影400 步（TPM1），演化200 步、投影600 步（TPM2），演化100 步、投影700 步（TPM3），演化100 步、投影1900 步（TPM4），其他計(jì)算條件不變，比較TPM 多尺度模擬在不同投影情況下的模擬結(jié)果。

圖13～圖18給出了選取不同的投影步長(zhǎng)時(shí)3種狀態(tài)方程控制下取得的模擬結(jié)果。

圖13和圖14為vdW 狀態(tài)方程控制下的模擬結(jié)果?？梢钥闯?，當(dāng)投影步數(shù)分別取值400、600、700和1900 時(shí)，TPM 多尺度模擬與LBM 單尺度模擬所得飽和密度曲線基本重合。投影數(shù)為400、600 和700 時(shí)，模型的最大偽速度保持一致；投影數(shù)為1900時(shí)，模型的最大偽速度有所增大。

圖13 不同投影步數(shù)下vdW 模型的飽和密度曲線對(duì)比Fig.13 Comparison of coexistence curves obtained from simulations using different extrapolation steps with experimental values for vdW EOS

圖14 vdW 模型最大偽速度的變化曲線對(duì)比Fig.14 Comparison of spurious current obtained from simulations using different extrapolation steps for vdW EOS

圖15 不同投影步數(shù)下R-K 模型的飽和密度曲線對(duì)比Fig.15 Comparison of coexistence curves obtained from simulations using different extrapolation steps with experimental values for R-K EOS

圖15和圖16為R-K 狀態(tài)方程控制下的模擬結(jié)果?？梢钥闯?，幾種情況下的TPM 多尺度模擬與LBM 單尺度模擬所得飽和密度曲線在氣相區(qū)域基本重合。在液相區(qū)域，投影步數(shù)為400 和600 時(shí)所得結(jié)果比LBM 單尺度模擬更接近實(shí)驗(yàn)值，TPM 與LBM 模擬的最大偽速度基本相同；投影數(shù)為700 時(shí)，模擬結(jié)果在低溫部分與實(shí)驗(yàn)值的偏差增大，模型的最大偽速度稍有增大；投影數(shù)為1900 時(shí)，模擬結(jié)果在低溫部分與實(shí)驗(yàn)值最為接近，模型的偽速度與其他幾種情況基本相同。

圖16 R-K 模型最大偽速度的變化曲線對(duì)比Fig.16 Comparison of spurious current obtained from simulations using different extrapolation steps for R-K EOS

圖17 不同投影步數(shù)下P-R 模型的飽和密度曲線對(duì)比Fig.17 Comparison of coexistence curves obtained from simulations using different extrapolation steps with experimental values for P-R EOS

圖18 P-R 模型最大偽速度的變化曲線對(duì)比Fig.18 Comparison of spurious current obtained from simulations using different extrapolation steps for P-R EOS

圖17和圖18為P-R 狀態(tài)方程控制下的模擬結(jié)果?？梢钥闯?，TPM 多尺度模擬與LBM 單尺度模擬所得飽和密度曲線在氣相區(qū)域基本重合。在液相區(qū)域，TPM 多尺度模擬比LBM 單尺度模擬所得曲線在低溫部分與實(shí)驗(yàn)值的偏差增大。當(dāng)選取700 和1900 的投影步長(zhǎng)時(shí)，TPM 模擬的溫度取值范圍減小，而最大偽速度與LBM 模擬結(jié)果基本保持一致，投影步長(zhǎng)為1900 的偽速度稍大。

在模擬精度和運(yùn)算效率方面，以R-K 模型為例（見表1和表2）。由表1可以看出，對(duì)于氣相密度，TPM 多尺度模擬的平均誤差和均方根誤差都大于LBM 模擬的誤差。當(dāng)投影步數(shù)為400、600 和700時(shí)，氣相密度的平均誤差隨投影步長(zhǎng)增大而增加；繼續(xù)增大投影步數(shù)至1900 時(shí)，相對(duì)誤差的平均值并沒(méi)有增大，而是保持為15%。液相密度方面，除投影步數(shù)為400 的情況之外，TPM 模擬的相對(duì)誤差均比LBM 單尺度模擬的相對(duì)誤差??；平均誤差最小為4.0%，分別是投影步數(shù)為600 和1900 時(shí)取得的。而從液相的均方根誤差可以看出，TPM 多尺度模擬的誤差（10%左右）均小于LBM 單尺度模擬的誤差（12.2%），說(shuō)明TPM 多尺度模擬在這一部分有效地提高了模型的準(zhǔn)確性。由表2可以看出，TPM多尺度模擬的CPU 耗時(shí)明顯小于LBM 單尺度模擬的CPU 耗時(shí)，可見TPM 多尺度模擬能夠大大提高模型的運(yùn)算效率；同時(shí)，由總體誤差可以看出TPM多尺度模擬的誤差均小于LBM 單尺度模擬的誤差，充分說(shuō)明了該方法能夠在提高運(yùn)算效率的同時(shí)提高模型的精確度，體現(xiàn)了其在復(fù)雜計(jì)算中的優(yōu)勢(shì)。當(dāng)投影步數(shù)從700（TPM3）增大至1900（TPM4）時(shí),誤差保持在8%左右，體現(xiàn)了該方法在相當(dāng)范圍內(nèi)的穩(wěn)定性。

綜上所述，在合理的投影范圍內(nèi)運(yùn)用TPM 多尺度模擬方法對(duì)氣液兩相分離過(guò)程進(jìn)行模擬研究能夠大幅降低對(duì)細(xì)微仿真器的依賴性，說(shuō)明所提多尺度模擬方法對(duì)現(xiàn)階段細(xì)微單尺度模擬中出現(xiàn)的巨大內(nèi)存開銷和CPU 耗時(shí)問(wèn)題實(shí)現(xiàn)了實(shí)質(zhì)性的解決，體現(xiàn)了其有效性和高效性。當(dāng)在細(xì)微尺度上模擬更復(fù)雜的氣液兩相流動(dòng)過(guò)程或更復(fù)雜的多尺度系統(tǒng)時(shí)，采用該多尺度方法能有效地解決細(xì)微單尺度模擬難以實(shí)現(xiàn)的難題，從而為復(fù)雜兩相流動(dòng)傳熱問(wèn)題的模擬研究提供新思路。

表1 LBM 和TPM 模擬與實(shí)驗(yàn)值的誤差Table 1 Error of simulation results obtained from LBM and TPM model verified by experimental values

表2 模擬運(yùn)算效率Table 2 Operational efficiency of simulations

4 結(jié) 論

本研究提出了一種用于模擬氣液兩相分離的免方程多尺度模擬方法，該方法以氣液兩相的格子Boltzmann（LB）模型作為介觀仿真器，LB 模型進(jìn)行少量演化步后，運(yùn)用二階伸縮式投影方式對(duì)介觀仿真演化得到的氣液相密度進(jìn)行有效的外推處理，能夠快速準(zhǔn)確地獲取后續(xù)演化步的信息，從而實(shí)現(xiàn)對(duì)氣液兩相分離過(guò)程的多尺度模擬研究。

對(duì)氣液兩相分離過(guò)程中密度分布和模型產(chǎn)生的偽速度的仿真分析表明，本研究所提出的多尺度模擬方法從介觀尺度出發(fā)真實(shí)地反映了氣液兩相分離的過(guò)程，在確保多尺度模擬的準(zhǔn)確性和真實(shí)性的基礎(chǔ)上有效地解決了細(xì)微尺度模擬中出現(xiàn)的CPU耗時(shí)過(guò)長(zhǎng)、內(nèi)存開銷過(guò)大的問(wèn)題。與以往的LBM單尺度模擬相比，本研究所提出的多尺度模擬方法具有明顯的計(jì)算優(yōu)勢(shì)以及多尺度特點(diǎn)，說(shuō)明所提方法對(duì)于研究更為復(fù)雜的流動(dòng)傳熱問(wèn)題非常具有潛力。通過(guò)EFM 方法的改進(jìn)以及LBM 模型本身的改進(jìn)，可為下一步工作實(shí)現(xiàn)流動(dòng)沸騰中氣泡動(dòng)力學(xué)行為等更復(fù)雜和更深層次的兩相流動(dòng)傳熱過(guò)程的數(shù)值模擬研究打下堅(jiān)實(shí)的基礎(chǔ)。另外，將該方法的理念運(yùn)用于其他能源動(dòng)力系統(tǒng)的模擬也是重要的發(fā)展方向，以進(jìn)一步推動(dòng)多尺度模擬的研究和發(fā)展。

符 號(hào) 說(shuō) 明

a，b——狀態(tài)方程中的系數(shù)

c——粒子遷移速度

ci——離散方向上的速度矢量

cs——LB 模型的格子聲速

F——變量的校正值

F——變量的預(yù)報(bào)值

F?——粒子間相互作用力

fi——粒子分布函數(shù)

——平衡態(tài)粒子分布函數(shù)

——格林函數(shù)

h——伸縮式投影方法的微觀投影數(shù)

i——離散速度方向

k——伸縮式投影方法的微觀演化數(shù)

p——壓力

R——通用氣體常數(shù)

T——溫度

t——時(shí)間坐標(biāo)

Δt——格子離散的時(shí)間步長(zhǎng)

δt——伸縮式投影方法的微觀演化步長(zhǎng)

u——宏觀速度

——平衡態(tài)速度

——偽勢(shì)函數(shù)

x——空間坐標(biāo)

Δx——格子離散的空間步長(zhǎng)

α，β，χ，?——Adams-Bashforth 方法中的系數(shù)

λ——提升算子中的系數(shù)

ρ——宏觀密度

τ——量綱1 松弛時(shí)間

φ——有效密度

ω——偏心因子

ωi——權(quán)重系數(shù)

[1]Guo Zhaoli (郭照立),Zheng Chuguang (鄭楚光),Li Qing (李青),Wang Nengchao ( 王能超).Lattice Boltzmann Method for Hydrodynamics (流體動(dòng)力學(xué)的格子 Boltzmann 方法)[M].Wuhan:Hubei Science & Technology Press,2002:158.

[2]He Yaling (何雅玲),Wang Yong (王勇),Li Qing (李慶).Lattice Boltzmann Method:Theory and Applications (格子Boltzmann 方法的理論及應(yīng)用) [M].Beijing:Science Press,2008:156.

[3]Zeng Jianbang (曾建邦),Li Longjian (李隆鍵),Liao Quan (廖全),Chen Qinghua (陳清華),Cui Wenzhi (崔文智),Pan Liangming (潘良明).Application of lattice Boltzmann method to phase transition process [J].Acta Physica Sinica(物理學(xué)報(bào)),2010,59 (1):178-185.

[4]Deen N G,van Sint Annaland M,Kuipers J A M.Multi-scale modeling of dispersed gas-liquid two-phase flow [J].Chemical Engineering Science,2004,59 (8):1853-1861.

[5]Armaou A,Kevrekidis I G,Theodoropoulos C.Equation-free gaptooth-based controller design for distributed complex/multiscale processes [J].Computers & Chemical Engineering,2005,29 (4):731-740.

[6]Kevrekidis I G,Gear C W,Hyman J M,Kevrekidis P G,Runborg O,Theodoropoulos C.Equation-free coarse-grained multiscale computation:enabling microscopic simulators to perform system-level analysis [J].Communications in Mathematical Sciences,2003,1 (4):715-762.

[7]Kevrekidis I G,Gear C W,Hummer G.Equation-free:the computer-aided analysis of complex multiscale systems [J].AⅠChE Journal,2004,50 (7):1346-1355.

[8]Gear C W.Projective integration methods for distributions [J].NEC Trans.,2001,130:1-9.

[9]Gear C W,Kevrekidis I G.Projective methods for stiff differential equations:problems with gaps in their eigenvalue spectrum [J].SⅠAM Journal on Scientific Computing,2003,24 (4):1091-1106.

[10]Cisternas J,Gear C W,Levin S,et al.Equation-free modelling of evolving diseases:coarse-grained computations with individual-based models [J].Proceedings of the Royal Society of London.Series A:Mathematical,Physical and Engineering Sciences,2004,460 (2050):2761-2779.

[11]Makeev A G,Kevrekidis I G.Equation-free multiscale computations for a lattice-gas model:coarse-grained bifurcation analysis of the NO + CO reaction on Pt (100) [J].Chemical Engineering Science,2004,59 (8):1733-1743.

[12]Erban R,Kevrekidis I G,Adalsteinsson D,et al.Gene regulatory networks:a coarse-grained,equation-free approach to multiscale computation [J].Journal of Chemical Physics,2006,124 (8):084106.

[13]Kevrekidis I G,Samaey G.Equation-free multiscale computation:algorithms and applications [J].Annual Review of Physical Chemistry,2009,60:321-344.

[14]Zeng Jianbang (曾建邦),Li Longjian (李隆鍵),Liao Quan (廖全),Cui Wenzhi (崔文智),Chen Qinghua (陳清華),Pan Liangming (潘良明).Simulation of boiling process with lattice Boltzmann method [J].Journal of Xi’an Jiaotong University(西安交通大學(xué)學(xué)報(bào)),2009,43 (7):25-29.

[15]Maluckov A M,Ishiguro S,?koric M M.A macro projective integration method in 2D microscopic system applied to nonlinear ion acoustic waves in a plasma [J].Communications in Computational Physics,2008,4 (3):556-574.

[16]Gear C W,Li J,Kevrekidis I G.The gap-tooth method in particle simulations [J].Physics Letters A,2003,316 (3):190-195.

[17]Yang Chen (楊晨),Peng Wei (彭偉).Telescoping projective methods for the multiscale simulation of reaction-diffusion processes [J].Journal of Engineering Thermophysics(工程熱物理學(xué)報(bào)),2011,32 (3):377-381.

[18]Shan X,Chen H.Lattice Boltzmann model for simulating flows with multiple phases and components [J].Physical Review E,1993,47 (3):1815.

[19]Shan X,Chen H.Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation [J].Physical Review E,1994,49 (4):2941.

[20]Yuan P,Schaefer L.Equations of state in a lattice Boltzmann model [J].Physics of Fluids(1994-present),2006,18 (4):42101.