王 博，陳一鳴，杜勝男，王衛(wèi)強

（遼寧石油化工大學(xué)石油天然氣工程學(xué)院，遼寧撫順113001）

21世紀(jì)，隨著計算機技術(shù)的出現(xiàn)和發(fā)展，數(shù)值分析逐漸成為和理論分析、實驗研究相并列的三大重要科學(xué)研究手段之一。數(shù)值方法可以對理論分析難以求解的復(fù)雜問題進行模擬求解，擴大了研究范圍。同時，數(shù)值方法可以節(jié)約實驗所帶來的成本，加快科學(xué)發(fā)展的速度。從方法論的角度，流體流動與傳熱的描述可以分為宏觀、介觀和微觀三個層次[1]。傳統(tǒng)計算流體力學(xué)（Computational Fluid Dynamics，CFD）主要從宏觀角度出發(fā)，建立連續(xù)介質(zhì)模型，將非線性微分方程轉(zhuǎn)換成非線性代數(shù)方程組，通過反復(fù)迭代進行求解。其中，具有代表性的是有限體積法（Finite volume method，F(xiàn)VM）和有限差分法（Finite Difference Method，F(xiàn)DM）[2]。

N-S方程是從復(fù)雜的流體運動中簡化出來的一個重要模型，該方程是一個非線性偏微分方程組，運用傳統(tǒng)CFD方法求解該方程非常困難，到目前為止，只有極少數(shù)非常簡單的流動問題才能求其精確解，大多數(shù)還是要通過離散方法求其數(shù)值解，尤其對于不可壓縮流動N-S方程的求解，壓力方程具有橢圓形，無法推進求解，收斂性較差，因此傳統(tǒng)CFD方法具有一定的局限性。國內(nèi)學(xué)者閻超等[3]、國外學(xué)者 H.Charles[4]、F.H.Moukalled等[5]對此進行了詳細(xì)的分析。為了克服傳統(tǒng)CFD方法求解困難及復(fù)雜邊界難以處理的問題，基于介觀層次的格子玻爾茲曼方法（Lattice Boltzmann Method，LBM）應(yīng)運而生，該方法是將流體離散成流體粒子，物理區(qū)域離散成格子，時間變量離散成時間步長，按照簡化的動力學(xué)模型，求解格線上宏觀特征值的平均值，進而求解流場[6]。相比之下，LBM方法具有物理意義清晰、邊界處理簡單、程序易于實施、模型健壯性高等優(yōu)點。因此，LBM方法是目前最具前途的數(shù)值模擬方法之一，國內(nèi)學(xué)者何雅玲等[7]、郭照立等[8]，國外學(xué)者S.Chen等[9]、穆罕默德·阿卜杜勒馬吉德[10]對該方法的應(yīng)用及發(fā)展都做出了較為詳盡的詮釋。李江濤等[11]將LBM方法與FDM方法相結(jié)合對頁巖氣的升尺度滲流進行了模擬；李校等[12]基于LBM理論，運用標(biāo)準(zhǔn)k-ε模型對二維室內(nèi)湍流空氣進行了數(shù)值模擬；劉海湖[13]運用LBM-FDM方法對不相溶的兩種表面活性劑混合問題進行了數(shù)值模擬，對混合后的表面張力做出了詳細(xì)的分析。

通過上述分析可知，LBM方法對于求解流體流動與傳熱等物理問題表現(xiàn)出較強的適用性，但目前對LBM方法的研究主要以應(yīng)用技術(shù)為主，并且眾多研究是將該方法與傳統(tǒng)的CFD方法相結(jié)合，而對LBM方法與傳統(tǒng)CFD方法相異之處的對比研究卻鮮有報道。因此，為了對比分析傳統(tǒng)CFD方法與LBM方法的差異，以無限大平板熱擴散問題為例，運用3種數(shù)值模擬方法對平板間的溫度分布進行研究，并將LBM方法的求解結(jié)果與FDM方法和FEA方法的求解結(jié)果進行對比，進而從直觀角度明確LBM方法的計算可行性及優(yōu)越性。

1 物理模型

1.1 問題描述

假設(shè)平板初始溫度T=0℃，當(dāng)時間t≥0時，平板左側(cè)施加高溫T≈100℃，平板長度L=100 m，平板厚度相較于平板長度極小，可忽略不計。平板幾何模型如圖1所示。

1.2 擴散方程

由于平板厚度遠(yuǎn)小于平板長度，且平板理論上為無限大平板。因此，可將平板間的熱擴散問題看成一維熱擴散問題。一維熱擴散方程可以表示為：

式中，T為溫度，℃；ρ為介質(zhì)的密度，kg/m3；c為比熱容，J/(kg?℃)；k為熱傳導(dǎo)系數(shù)，W/(m?℃)；α為熱擴散系數(shù)，取 0.25[14]。

圖1 平板幾何模型Fig.1 Flat geometry model

2 數(shù)值方法

2.1 有限差分法

2.1.1 基本思路 FDM方法的基本思想是把連續(xù)的定解區(qū)域用有限個離散點所構(gòu)成的網(wǎng)格來代替，離散點稱作網(wǎng)格節(jié)點；用網(wǎng)格節(jié)點上定義的離散變量函數(shù)近似代替連續(xù)定解變量的函數(shù)，用有限差分方程組代替原微分代數(shù)方程組，用離散差商解近似代替連續(xù)微商解，然后再利用插值方法將離散解還原為計算域的近似解[15]。

FDM方法的求解步驟為：（1）區(qū)域離散化，即將偏微分方程的求解區(qū)域細(xì)分成由有限個格點組成的網(wǎng)格。（2）近似替代，即用有限差分計算值替代每一個格點的導(dǎo)數(shù)值。（3）逼近求解，即用一個插值多項式及其微分來代替原偏微分方程進行求解。

2.1.2 區(qū)域離散化 針對求解區(qū)域離散化的節(jié)點個數(shù)及節(jié)點間距沒有固定的要求，為了增加求解效率，對計算域進行均勻離散化處理，結(jié)果如圖2所示。

圖2 離散化節(jié)點分布Fig.2 Discretized node distribution

2.1.3 近似替代 按照不同的劃分標(biāo)準(zhǔn)，有限差分的近似方法可分為：（1）按照差分精度可分為一階格式、二階格式和高階格式。（2）按照差分的空間形式可分為中心格式和逆風(fēng)格式。（3）按照時間因子的影響可分為顯格式、隱格式、顯隱交替格式等。針對平板間的熱擴散問題，選擇顯式有限差分方法，具體形式為：

時間導(dǎo)數(shù)采用一階前向差分，即

空間導(dǎo)數(shù)采用二階中心差分，即

式中，β為相關(guān)系數(shù)，取 0.01；Δx為步長；Δt為時間間隔。

因此，節(jié)點i處的擴散方程可進一步表示為：

式中，ω為松弛頻率；α為擴散系數(shù)，m2/s；t為時間，s。

2.2 格子玻爾茲曼法

2.2.1 基本思路 LBM方法是一種基于介觀尺度的計算流體力學(xué)方法。該方法相比于傳統(tǒng)模擬方法，具有介于微觀分子動力學(xué)模型和宏觀連續(xù)模型的介觀模型特點。因此，該方法具備流體相互作用描述簡單、復(fù)雜邊界易于設(shè)置、易于并行計算、程序易于實施等優(yōu)勢[16]。玻爾茲曼理論及LBM方法已經(jīng)廣泛地被認(rèn)為是描述流體運動與處理工程問題的有效手段。LBM方法的求解流程如圖3所示。2.2.2 求解設(shè)置 針對恒溫?zé)o限大平板的熱擴散問題，基于LBM理論，選用D2Q4模型，平板軸向劃分100個格子區(qū)間，平板左側(cè)端面初始化溫度為100℃，熱擴散系數(shù)α為 0.25。圖 4、5為二維（D2Q4）格子排布及熱擴散格子分布示意。

根據(jù)Mohamad提出的，溫度分布函數(shù)的動力學(xué)方程可寫成[17]：

式中，ck為格子運動速度，m/s；Ωk為碰撞算子，根據(jù)BGK近似求解，即

式中，τ為達到平衡態(tài)分布函數(shù)的松弛時間，s；Teqk為

平衡態(tài)分布函數(shù)，定義式如下：

式中，Tk(x,t)為溫度分布函數(shù)；wk為k方向權(quán)重因子。

采用BGK近似后，溫度分布函數(shù)離散化表達式為：

圖3 LBM方法求解流程Fig.3 LBM method solution flow chart

圖4 一維格子排布示意Fig.4 Schematic diagram of one-dimensional grid arrangement

圖5 一維熱擴散格子分布示意Fig.5 Schematic diagram of one-dimensional thermal diffusion lattice distribution

將式（10）代入式（11）化簡得一維空間擴散方程為：

2.3 有限體積法

2.3.1 基本思路 FVM方法作為目前CFD領(lǐng)域最成熟的算法。該算法是將流體的Euler控制方程在單元控制體內(nèi)進行積分后離散求解[18]。就離散方法而言，有限體積法可視作有限單元法和有限差分法的中間物。該方法求解流程如圖6所示。

圖6 FVM法求解流程Fig.6 FVM method solution flow chart

目前常用的CFD軟件，例如Fluent、CFX等都是基于該方法。因此，采用Fluent作為FVM方法的代表對平板間的熱擴散問題進行求解。恒溫?zé)o限大平板的上、下壁面為絕熱邊界，因此平板法向方向無溫度變化現(xiàn)象，由于計算域內(nèi)介質(zhì)的溫度僅一個方向有變化且傳熱系數(shù)不變，每一個物質(zhì)元流入與流出的熱量均相等，平板間溫度的變化情況與模型的維度無關(guān)，可將計算域簡化為一維溫度場，為了更好地展現(xiàn)流場細(xì)節(jié)，利用FVM方法對問題進行求解時，以三維的模型進行展示。

2.3.2 網(wǎng)格劃分 采用結(jié)構(gòu)化網(wǎng)格對平板進行網(wǎng)格劃分，平板長度遠(yuǎn)遠(yuǎn)大于寬度且計算域為規(guī)則對稱區(qū)域，因此選取平板上截面的二維模型進行模擬分析，網(wǎng)格具體劃分如圖7所示。

圖7 網(wǎng)格劃分Fig.7 Mesh generation

2.3.3 邊界條件 平板左側(cè)端面設(shè)定初始溫度為100℃，右側(cè)端面無限延展，較遠(yuǎn)端面處設(shè)置為0℃，上下設(shè)置為壁面；計算域材料的擴散系數(shù)為0.25。

3 結(jié)果分析

3.1 FDM方法結(jié)果分析

基于FDM方法，運用Fortran編譯器進行程序編寫，設(shè)置空間域離散節(jié)點之間的距離Δx=1.0，時間域離散步長Δt=0.5，軸線方向離散格子數(shù)為100，平板左側(cè)端面初始化溫度為1℃，平板間介質(zhì)的熱擴散系數(shù)α=0.25，通過歸一化求解[19]可得時間步數(shù)為400，求解得平板間溫度變化情況如圖8所示。

由圖8可知，F(xiàn)DM方法求解得到的平板間溫度變化呈現(xiàn)遞減的趨勢，且變化曲線近似服從對數(shù)分布函數(shù)，即溫度遞減的速率在逐漸降低；恒溫條件下，平板左側(cè)端面溫度取最大值，沿著軸線方向溫度逐漸減小，軸向20 m附近溫度降至0℃附近，30 m后至無窮遠(yuǎn)處溫度始終保持在0℃。

3.2 LBM方法結(jié)果分析

基于LBM方法，運用Fortran編譯器進行程序編寫，沿平板軸線方向均勻劃分100個離散格子，設(shè)置時間步長Δt=1.0，格子間距Δx=1.0，平板左側(cè)端面初始化溫度為100℃，平板間介質(zhì)的熱擴散系數(shù)α=0.25，計算當(dāng)t=200 s時，平板間的溫度分布情況，結(jié)果如圖9所示。由圖9可知，LBM方法求解得到的平板間溫度變化規(guī)律與FDM方法求解的結(jié)果近似相同，兩種方法下的溫度變化曲線吻合度較高。

圖9 基于LBM方法的溫度變化Fig.9 Temperature change based on LBM method

為了驗證兩種求解方法得到的溫度變化分布情況，選取軸向變動數(shù)據(jù)進行擬合處理，結(jié)果如圖10所示。由圖10可知，LBM方法與FDM方法求解得到的平板間溫度變化曲線近似服從對數(shù)分布，其擬合曲線為y=-35.89lnx+117.47，曲線擬合優(yōu)度R2為96.79%。

圖10 溫度變化擬合Fig.10 Temperature change fitting

3.3 FVM方法結(jié)果分析

針對恒溫條件下的平板間熱擴散問題，運用Ansys Workbench中穩(wěn)態(tài)熱分析模型對其進行求解，由于平板軸線方向尺寸為無限大，平板厚度及寬度遠(yuǎn)遠(yuǎn)小于其長度，可選取平板上表面作為分析系統(tǒng)的代表計算域；計算域左側(cè)邊界設(shè)置高溫T=100℃，平板間介質(zhì)的熱擴散系數(shù)α=0.25，求解得平板間溫度變化情況如圖11所示。

圖11 溫度云圖Fig.11 Temperature cloud map

由圖11可知，基于FVM方法得到的平板間溫度變化云圖沿軸線方向逐漸減小，且溫度條帶近似均勻分布；平板左側(cè)端面溫度取最大值，平板右側(cè)區(qū)域溫度逐漸降低至0℃，即平板無限遠(yuǎn)處溫度降到最低；為直觀地分析溫度沿軸線方向的變化規(guī)律，繪制溫度變化曲線如圖12所示。由圖12可知，F(xiàn)VM方法求解得到的平板間溫度呈減小的趨勢，且減小的速率基本保持不變，即平板間的溫度變化曲線近似呈現(xiàn)線性分布，該結(jié)論與上述兩種方法相比略有差異。

3.4 對比分析

圖13為3種求解方法得到的溫度變化。由圖13可知，3種數(shù)值方法皆可以很好地模擬恒溫邊界條件下的熱擴散問題，且平板間的溫度變化皆呈減小的趨勢；根據(jù)溫度變化曲線的切線變化率可知，LBM方法和FDM方法的溫度變化曲線表現(xiàn)非線性趨勢，且溫度變化情況較為相近。相比之下，F(xiàn)VM方法求解結(jié)果表現(xiàn)出較大的差異，其溫度變化曲線近似線性趨勢；根據(jù)熱擴散理論[20-21]可知，介質(zhì)熱擴散速度與溫度呈正比例變化，即溫度越高，分子無規(guī)則運動越劇烈，分子間碰撞頻率及強度會加大，單個分子動能增加，最終會導(dǎo)致平板間的分子總動能增加；由LBM方法及FDM方法的求解結(jié)果可知，平板左側(cè)端面溫度達到最大值，因此擴散初期速度較快，即溫度變化曲線切線斜率達到最大值，隨著溫度的逐漸降低，平板間熱擴散速度逐漸減小，即溫度曲線切線斜率越來越?。煌ㄟ^對比可知，針對恒溫條件的熱擴散問題，LBM方法和FDM方法求解得到的溫度變化曲線熱擴散速度在逐漸減小，符合熱擴散機理，而FVM方法求解得到的溫度變化曲線熱擴散速度保持不變，說明該方法對于熱擴散的求解是可行的，但是求解精度要低于另外兩種方法。

圖12 基于FVM方法的溫度變化Fig.12 Temperature change based on FVM method

圖13 平板軸向溫度變化Fig.13 Plate axial temperature curve

4 結(jié) 論

針對無限大平板的熱擴散問題，運用3種CFD方法對溫度分布進行了求解并得到：

（1）相同物理模型條件下，針對恒溫?zé)徇吔鐥l件下的熱擴散問題，LBM方法和FDM方法的求解準(zhǔn)確性要略優(yōu)于FVM方法。

（2）運用歸一化和尺度準(zhǔn)則，計算相同條件下的熱擴散問題，F(xiàn)DM的求解步數(shù)要大于LBM方法，即針對擴散問題，LBM方法求解時間比較短，計算效率較高。

（3）基于LBM方法的一維熱擴散問題，該方法具有較好的準(zhǔn)確性和適用性，為后期研究二維、三維的熱擴散問題提供了一定的參考價值。