氣液兩相分離的免方程多尺度模擬方法

楊晨，何航行

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

引 言

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

近年來，有學(xué)者提出了一種計算機(jī)輔助分析復(fù)雜系統(tǒng)的新型多尺度模擬框架——免方程方法(equation-free method，EFM)[6-7]。其關(guān)鍵思路是繞過從微細(xì)尺度模型中獲得宏觀模型方程這一傳統(tǒng)過程，為復(fù)雜系統(tǒng)多尺度模擬提供了一種新的思路。目前免方程方法的研究在國外發(fā)展迅速，但國內(nèi)只有少數(shù)文獻(xiàn)提及。早在20世紀(jì)末，Theodoropoulos等[6]已涉及到粗?；瘯r間步進(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)勢，提出了一種用于模擬氣液兩相分離過程的免方程多尺度方法。利用該方法模擬水在van der Waals（vdW）、Redlich-Kwong（R-K）和Peng-Robinson（P-R）狀態(tài)方程控制下的相變過程，發(fā)現(xiàn)多尺度模擬能夠從微細(xì)底層很好地模擬出氣液兩相分離過程中密度的分布，所得結(jié)果與LBM 單尺度模擬十分吻合；多尺度模擬計算產(chǎn)生的偽速度與LBM 單尺度模擬結(jié)果基本相同，不會導(dǎo)致偽速度增加。這說明該方法在提高了計算效率的同時保證了計算的穩(wěn)定性，驗證了該多尺度模擬方法的準(zhǔn)確性和高效性。通過免方程方法在氣液兩相分離過程中的應(yīng)用研究，可為進(jìn)一步從細(xì)微尺度出發(fā)對流動沸騰中氣泡行為等復(fù)雜多相流過程的多尺度模擬奠定基礎(chǔ)。

1 免方程多尺度方法

1.1 基本框架

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

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

1.2 粗?；队胺e分

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

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

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

圖3 CPI 模擬全過程Fig.3 Entire CPI simulation process diagram

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

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

2 基本模型

2.1 格子Boltzmann 模型

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

式中，x、t分別表示空間及時間坐標(biāo)；i=0,1,2,…,b表示每個格點的不同離散速度方向；ci表示各個離散方向上的速度矢量；Δt表示格子離散的時間步長；τ為量綱1 松弛時間；為平衡態(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)速度。密度和動量的宏觀參數(shù)分別為

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

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

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

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

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

2.2 提升算子和約束算子

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

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

約束算子則選取為LB 模型密度的宏觀參數(shù)形式

2.3 狀態(tài)方程

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

式中，cs為LB 模型的格子聲速，對于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。壓力對密度求一階導(dǎo)數(shù)和二階導(dǎo)數(shù)并令其等于零，則可得到臨界點的數(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ù)值計算區(qū)域的格子劃分為256×256，平衡態(tài)分布函數(shù)采用D2Q9 模型，松弛因子 1τ=，上下邊界條件和左右邊界條件均為周期性邊界條件，偽勢密度函數(shù)為任意常數(shù)。初始密度值ρ=ρc，為了使方程演化，在初始密度附近取0.01 的擾動。模擬結(jié)果均在以內(nèi)存4GB、CPU 主頻3.1GHz 的四核Xeon 處理器為計算平臺的同一臺計算機(jī)上運(yùn)行得到。

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

圖4給出了不同狀態(tài)方程控制下LB 模型模擬所得的密度曲線與實驗值的對比，圖中的實驗值根據(jù)水蒸氣性質(zhì)圖表繪制得出，離散點為數(shù)值模擬計算得出的結(jié)果。通過對比可以發(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é)果偏離實驗值較多，密度相對誤差較大，氣相密度比實驗值大，隨著溫度升高模擬結(jié)果與實驗值的相對誤差逐漸減小，而液相密度模擬結(jié)果比實驗值小，與實驗值的相對誤差基本維持在13%左右。R-K 和P-R 狀態(tài)方程模擬所得曲線與實驗值吻合較好，對于氣相密度R-K狀態(tài)方程得到的結(jié)果與實驗值最接近[圖4(b)]，對于液相密度P-R 狀態(tài)方程得到的結(jié)果較好。

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

圖4 不同狀態(tài)方程控制下的飽和密度曲線與實驗值的對比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 模擬的對比

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖15 不同投影步數(shù)下R-K 模型的飽和密度曲線對比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é)果?？梢钥闯觯瑤追N情況下的TPM 多尺度模擬與LBM 單尺度模擬所得飽和密度曲線在氣相區(qū)域基本重合。在液相區(qū)域，投影步數(shù)為400 和600 時所得結(jié)果比LBM 單尺度模擬更接近實驗值，TPM 與LBM 模擬的最大偽速度基本相同；投影數(shù)為700 時，模擬結(jié)果在低溫部分與實驗值的偏差增大，模型的最大偽速度稍有增大；投影數(shù)為1900 時，模擬結(jié)果在低溫部分與實驗值最為接近，模型的偽速度與其他幾種情況基本相同。

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

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

圖18 P-R 模型最大偽速度的變化曲線對比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 單尺度模擬所得曲線在低溫部分與實驗值的偏差增大。當(dāng)選取700 和1900 的投影步長時，TPM 模擬的溫度取值范圍減小，而最大偽速度與LBM 模擬結(jié)果基本保持一致，投影步長為1900 的偽速度稍大。

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

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

表1 LBM 和TPM 模擬與實驗值的誤差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)用二階伸縮式投影方式對介觀仿真演化得到的氣液相密度進(jìn)行有效的外推處理，能夠快速準(zhǔn)確地獲取后續(xù)演化步的信息，從而實現(xiàn)對氣液兩相分離過程的多尺度模擬研究。

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

符 號 說 明

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

c——粒子遷移速度

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

cs——LB 模型的格子聲速

F——變量的校正值

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

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

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

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

——格林函數(shù)

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

i——離散速度方向

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

p——壓力

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

T——溫度

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

Δt——格子離散的時間步長

δt——伸縮式投影方法的微觀演化步長

u——宏觀速度

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

——偽勢函數(shù)

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

Δx——格子離散的空間步長

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

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

ρ——宏觀密度

τ——量綱1 松弛時間

φ——有效密度

ω——偏心因子

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

[1]Guo Zhaoli (郭照立),Zheng Chuguang (鄭楚光),Li Qing (李青),Wang Nengchao ( 王能超).Lattice Boltzmann Method for Hydrodynamics (流體動力學(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é)報),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é)報),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é)報),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.