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

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

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

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

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

1 物理模型

1.1 問(wèn)題描述

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

1.2 擴(kuò)散方程

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

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

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

2 數(shù)值方法

2.1 有限差分法

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

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

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

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

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

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

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

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

因此，節(jié)點(diǎn)i處的擴(kuò)散方程可進(jìn)一步表示為：

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

2.2 格子玻爾茲曼法

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

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

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

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

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

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

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

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

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

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

將式（10）代入式（11）化簡(jiǎn)得一維空間擴(kuò)散方程為：

2.3 有限體積法

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

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

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

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

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

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

3 結(jié)果分析

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

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

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

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

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

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

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

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

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

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

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

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

3.4 對(duì)比分析

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

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

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

4 結(jié) 論

針對(duì)無(wú)限大平板的熱擴(kuò)散問(wèn)題，運(yùn)用3種CFD方法對(duì)溫度分布進(jìn)行了求解并得到：

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

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

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