李勇霖，陳榮華，田文喜，秋穗正，蘇光輝

(1.西安交通大學(xué) 動力工程多相流國家重點實驗室，陜西 西安 710049；2.西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

20世紀后期，隨著計算機運算能力的提高，數(shù)值模擬技術(shù)得到了飛速發(fā)展。數(shù)值模擬技術(shù)從對物理問題的數(shù)學(xué)描寫出發(fā)，使用各種數(shù)值方法求得物理問題的數(shù)值解，分析數(shù)值解進而探求物理問題的本質(zhì)規(guī)律。數(shù)值模擬的實現(xiàn)關(guān)鍵在于有效算法的開發(fā)，在計算流體力學(xué)(CFD)領(lǐng)域尤為如此。現(xiàn)有傳統(tǒng)的數(shù)值方法如有限差分法(FDM)、有限元法(FEM)及有限體積法(FVM)等都是在對求解區(qū)域劃分網(wǎng)格后進行數(shù)值求解。傳統(tǒng)的數(shù)值算法雖已日趨成熟，但在處理復(fù)雜幾何形狀、高速撞擊、流固耦合及自由面追蹤等特殊問題時，還面臨著網(wǎng)格扭曲等挑戰(zhàn)性問題[1]。

相比于網(wǎng)格法，使用移動粒子方法來進行數(shù)值模擬可對上述問題進行很好的分析。尤其是流體的拓撲變形不能通過移動及固定網(wǎng)格進行分析，卻能通過移動粒子方法進行分析。移動粒子方法的另一個優(yōu)點是，在計算對流時，直接通過粒子的運動計算而不需考慮數(shù)值擴散問題[2]。因為這些優(yōu)勢，近幾十年來涌現(xiàn)出了粒子數(shù)值分析方法。這些粒子數(shù)值分析方法大致能分成兩種類型：第1類是建立在概率論模型上，如分子動力學(xué)方法[3]、蒙特卡羅方法[4]等；第2類是建立在確定論模型上，如光滑粒子動力學(xué)(SPH)方法[5]和移動粒子半隱式(MPS)方法[2]等。第1類方法需對大量粒子的行為進行長時間跟蹤，直到獲得一精確的平均值，相對而言，第2類方法不需如此大的內(nèi)存和計算機時間。

在此之后，一些粒子法與網(wǎng)格法耦合的方法被提出，如MPS-MAFL[6]及有限體積粒子(FVP)[7]方法。MPS-MAFL方法的液相通過粒子離散，氣相不需要離散而只是為液相提供壓力邊界條件，氣相的壓力通過氣體狀態(tài)方程確定；FVP方法中，氣相使用歐拉網(wǎng)格表示，液相通過拉格朗日粒子進行離散。FVP方法實際是將MPS方法計算中的粒子定義為1個正方體，其體積即為原粒子直徑的立方。在實際計算中，如果按固定的有限體積進行劃分，當(dāng)粒子發(fā)生移動后，粒子排布不再如初始狀態(tài)那般規(guī)則，有限體積在計算域內(nèi)會出現(xiàn)重疊，而且在相鄰有限體積粒子中間會出現(xiàn)間隙。

由于拉格朗日粒子法的本質(zhì)，MPS方法中的固壁邊界條件很難像網(wǎng)格法一樣嚴格實施，這是因為當(dāng)計算的粒子位于計算域內(nèi)部時，粒子的支持域沒有被邊界截斷，計算結(jié)果能達到較高精度，當(dāng)粒子位于邊界附近時，支持域會被邊界截斷，此時，普通的核近似方法在邊界的計算上不再適用。這個問題是限制粒子法發(fā)展的關(guān)鍵性問題。對于此類問題，在SPH方法的改進上已做了較多嘗試[8-15]，包括邊界施加方法和虛粒子法。邊界施加方法是對邊界的粒子施加一個額外的作用力，實現(xiàn)起來較為簡單，但物理意義不明確且數(shù)值穩(wěn)定性較差；虛粒子法則在邊界外側(cè)額外布置一些粒子以修復(fù)邊界附近粒子缺失現(xiàn)象，此法理論依據(jù)較完善、精度較高，但對于復(fù)雜形狀邊界的處理比較困難。在計算存在自由表面的流體流動問題時，常通過某種條件確定自由表面粒子，并認為自由表面粒子的壓力為0，通過此處理，能使MPS方法在計算自由表面流動問題中獲得較精確的解。

在計算流體力學(xué)領(lǐng)域，傳熱問題的研究是普遍存在的，當(dāng)使用MPS等粒子方法進行傳熱計算時，也會出現(xiàn)計算精度的問題。對于傳熱問題而言，如果計算邊界不是恒溫，那么將粒子近似法直接應(yīng)用于求解傳熱問題是不正確的，即粒子模型的流動邊界處理方式不能直接應(yīng)用于傳熱問題的計算，即便是FVP方法在求解傳熱問題時得到的結(jié)果也不盡人意。

為克服MPS方法存在的不足，本文在MPS方法的基礎(chǔ)上，提出一種粒子-網(wǎng)格混合求解方法，即使用MPS-FVM耦合的方法對對流傳熱問題進行求解，既能保留粒子法在處理復(fù)雜幾何形狀、高速撞擊、流固耦合和自由面追蹤問題的優(yōu)勢，又能在處理傳熱問題時，消除求解域邊界的計算誤差以及整體計算獲得精確的數(shù)值解。

1 數(shù)學(xué)物理模型

1.1 控制方程

粒子-網(wǎng)格混合方法中的連續(xù)性方程、動量守恒方程和能量守恒方程分別為：

(1)

(2)

Ein+Eg=Est

(3)

其中：ρ為密度；u為速度；p為壓力；ν為運動黏性系數(shù)；f為外力；t為時間；Ein為單位時間控制體積內(nèi)能量凈入量；Eg為單位時間控制體積內(nèi)能量產(chǎn)生量；Est為控制體及內(nèi)能量變化率。

在粒子-網(wǎng)格混合方法中，MPS方法用于求解動量守恒方程，F(xiàn)VM用于求解能量守恒方程，則前兩個方程以微分形式給出，能量方程以積分形式給出。

1.2 粒子相互作用模型

在MPS方法中，粒子之間的相互作用通過核函數(shù)實現(xiàn)，每個粒子僅與其作用半徑內(nèi)的有限個粒子發(fā)生相互作用。本文算例中使用的核函數(shù)為：

w(|rj-ri|)=

(4)

其中：ri、rj分別為i粒子和j粒子的坐標(biāo)；re為粒子作用半徑。

i粒子的粒子數(shù)密度ni定義為：

(5)

對于不可壓縮流體，流體密度保持定常。由于流體密度與粒子數(shù)密度呈正比，其粒子數(shù)密度維持為1個常數(shù)。

某物理量在i粒子處的梯度矢量可通過梯度模型求解，梯度模型在MPS方法中被用來求解壓力的梯度項。在動量守恒方程中，擴散項以Laplace算子表示，使用Laplace模型求解。梯度模型和Laplace模型分別為：

(6)

(7)

其中：φ表示物理量；d為計算的空間維度；n0為初始粒子數(shù)密度；λ為鄰居粒子和中心粒子間距的方差，它的出現(xiàn)是為了考慮粒子擴散范圍的時間效應(yīng)。

1.3 能量守恒方程的求解

在對粒子運動進行求解后，將獲得每個粒子的坐標(biāo)位置，然后根據(jù)粒子坐標(biāo)，對計算域進行網(wǎng)格劃分，每個粒子中心將對周圍有1個控制體積。而這個控制體積是根據(jù)該時刻的粒子間位置決定的。實際上，在計算中，求解獲得粒子在該時層的坐標(biāo)后，可對i粒子的鄰居粒子進行分析，每個鄰居粒子對i粒子有1個作用面積，在二維計算中，作用面積Aj定義為：

Aj=ljl0

(8)

其中：lj為j粒子與i粒子連線的垂直面作用線；l0為初始粒子間距。

那么對于i粒子的控制體積則為：

(9)

對于無內(nèi)熱源的導(dǎo)熱問題，式(3)的能量方程可寫成：

(10)

其中：κ為導(dǎo)熱系數(shù)；Ti和Tj分別為i、j粒子溫度；cp為比定壓熱容；上標(biāo)n表示時層。

1.4 算法流程

粒子-網(wǎng)格混合方法求解對流傳熱問題的算法流程如圖1所示，在求解實際問題時，先根據(jù)具體問題設(shè)定粒子初始布置，然后使用MPS方法求解流動問題，求解得到粒子坐標(biāo)位置后對其進行網(wǎng)格劃分，再使用FVM對傳熱進行求解。

2 導(dǎo)熱問題的驗證

圖2示出一維半無限平板導(dǎo)熱問題算例模型，其計算域尺寸為41×20，左端紅色粒子設(shè)定為恒溫，即Tl=1 K，右邊綠色粒子初始溫度為0 K，且溫度可變。圖2b是按照粒子位置劃分的FVM網(wǎng)格，可看出，在平板線邊緣，網(wǎng)格的體積為原始體積的一半；而在平板角邊緣，網(wǎng)格的體積則為原體積的1/4。平板的物性分別取為密度ρ=1 000 kg/m3、比熱容c=1 J/(kg·K)、導(dǎo)熱系數(shù)κ=1 W/(m·K)。

圖1 粒子-網(wǎng)格混合方法的流程Fig.1 Flow chart of particle-grid hybrid method

算例的初始條件如下。t=0 s：T=Tl=1 K，x=0；T=0 K，0

在這個算例中，使用MPS方法及粒子-網(wǎng)格混合方法分析了不同時刻平板溫度的分布，結(jié)果如圖3所示，并將計算結(jié)果與精確解進行了比較，結(jié)果如圖4所示?？砂l(fā)現(xiàn)，在導(dǎo)熱前期，使用MPS方法求得的結(jié)果與解析解之間有較大差距，特別是在靠近熱源(即恒溫)處的粒子，溫度誤差相當(dāng)大，但隨時間的變化，這個誤差逐漸減小。使用粒子-網(wǎng)格混合方法進行計算時，前期誤差雖然存在，但誤差較小，且隨時間的變化，誤差逐漸展平至可忽略不計。圖5示出距離恒溫源最近(x=0.025)的一列粒子的縱向溫度分布。實際計算中，一維半無限大平板的上、下端處設(shè)置為絕熱，等溫線應(yīng)呈直線分布，使用MPS方法計算時出現(xiàn)了等溫線彎曲現(xiàn)象，這是由于粒子相互作用模型，即Laplace模型不適合求解溫度Laplace算子。使用粒子-網(wǎng)格混合方法計算的結(jié)果，前期雖存在誤差，但誤差很小，且誤差隨時間的變化而逐漸變得更小，后期誤差幾乎可忽略不計，且在計算時間范圍內(nèi)，等溫線一直呈直線分布，符合理論分析的結(jié)果。

a——MPS粒子布置；b——FVM網(wǎng)格劃分圖2 半無限平板導(dǎo)熱算例Fig.2 Example of semi-infinite plate thermal conductivity

a——前期；b——后期圖3 半無限平板溫度分布Fig.3 Temperature distribution in semi-infinite plate

a——前期；b——后期圖4 數(shù)值計算溫度與精確解間的誤差Fig.4 Error between numerical temperature and exact solution

a——前期；b——后期圖5 x=0.025處的溫度分布Fig.5 Temperature distribution at x=0.025

3 方腔內(nèi)自然對流

不同溫度的兩個垂直側(cè)邊方腔內(nèi)的自然對流是用來測試不同程序計算傳熱問題有效性的典型算例。本文選取了方腔自然對流的標(biāo)準題[16]進行計算，所述問題如圖6所示，方腔左側(cè)壁面為高溫恒定溫度Tl=1 K，右側(cè)壁面為低溫恒定溫度Tr=0 K，上、下壁面絕熱。腔體為正方形，邊長為L，流體Pr為0.701，四周為固體壁面，取無滑移條件，流體符合Boussinesq假設(shè)。

本文方腔內(nèi)自然對流的初始條件如下。t=0 s：u=v=0，T=Tl=1 K，x=0；u=v=0，T=Tr=0 K，0

圖6 物理模型Fig.6 Physical model

在模擬計算方腔自然對流時，設(shè)定了如圖6的初始條件及邊界條件，隨計算的進行，逐漸達到穩(wěn)態(tài)過程。圖7示出瑞利數(shù)Ra=106的瞬態(tài)及穩(wěn)態(tài)結(jié)果。

a～d——瞬態(tài)；e——穩(wěn)態(tài)圖7 Ra=106的瞬態(tài)及穩(wěn)態(tài)結(jié)果Fig.7 Transient and steady-state results with Ra=106

a——溫度等溫線圖；b——流線圖圖9 Ra=105的計算結(jié)果Fig.9 Calculation result of Ra=105

本文還計算了Ra分別為103及105時方腔自然對流的結(jié)果，如圖8、9所示，當(dāng)Ra較小時，方腔內(nèi)部的等溫線幾乎沿豎直方向，流動的典型特性是在方腔中間出現(xiàn)了1個大渦。當(dāng)Ra增大時，方腔內(nèi)部等溫線逐漸變?yōu)樗?，并且在熱冷壁面附近的薄邊界層?nèi)仍保持垂直，這說明，當(dāng)Ra較小時，換熱主要是熱傳導(dǎo)為主，隨著Ra的增加，換熱逐漸變?yōu)橛蓪α鲹Q熱為主導(dǎo)地位。Ra增大后，方腔內(nèi)部的大渦分裂成兩個。將計算所得的平均Nu與基準題解進行對比，結(jié)果如圖10所示，可發(fā)現(xiàn)粒子-網(wǎng)格混合方法在計算方腔內(nèi)自然對流時具有較高的精度。本文選擇了方腔對流達到穩(wěn)定狀態(tài)后的平均Nu進行對比。為了對計算所得的平均Nu進行不確定度分析，本文對穩(wěn)定狀態(tài)后20個時刻的平均Nu進行分析，可得到Ra為103～106的方腔對流平均Nu的標(biāo)準偏差uNu為：

圖10 平均Nu結(jié)果對比 Fig.10 Comparison of average Nu

平均Nu的標(biāo)準偏差如圖11所示。由圖11可見，粒子-網(wǎng)格混合方法計算得到的平均Nu的標(biāo)準偏差在2%以內(nèi)，說明了粒子-網(wǎng)格混合方法在計算對流換熱問題時所得到的結(jié)果是比較可靠的。

為研究粒子-網(wǎng)格耦合方法的工程適用性，對不同計算量的方腔對流進行了計算時間對比，因為粒子-網(wǎng)格混合方法已實現(xiàn)OpenMP的單節(jié)點多核并行算法計算，圖12示出不同計算量的方腔對流問題在使用計算機不同核數(shù)的單個時間步長計算耗時，計算機型號為RR CPU E5-2697 v3 2.6 GHz。由圖12可見，計算方腔內(nèi)對流問題時，使用的粒子數(shù)量越多，單個時間步長耗時越長，當(dāng)使用更多的CPU參與計算時，能在一定程度上縮小單個時間步長耗時，總地來說使用并行算法的粒子-網(wǎng)格混合方法在一定程度上提高了計算效率。目前在并行效率上，還沒有完成更高并行效率的MPI并行及GPU并行算法的搭建。若完成更高并行效率的算法搭建，基于粒子-網(wǎng)格耦合方法有可能應(yīng)用于大規(guī)模的工程計算。

圖11 平均Nu的標(biāo)準偏差Fig.11 Standard deviation of average Nu

圖12 不同粒子數(shù)量、CPU數(shù)量的單個時間步長耗時 Fig.12 Usage time of single step under different particle numbers and different numbers of CPU

4 結(jié)論

本文針對MPS方法在求解導(dǎo)熱問題時存在的邊界不精確情況，對MPS方法進行了改進，提出了使用MPS方法求解動量守恒方程、使用FVM求解能量守恒方程的粒子-網(wǎng)格混合方法。該方法使用MPS方法求解流動問題，保留了MPS方法在求解大變形問題及相變問題方面所具有的優(yōu)勢，且還解決了原方法求解導(dǎo)熱方程不精確的問題。本文研究的粒子-網(wǎng)格混合方法僅適用于二維情況，對于三維情況還需做進一步研究。