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

(1.西安交通大學(xué) 動(dòng)力工程多相流國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710049；2.西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049)

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

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

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

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

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

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

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

1.1 控制方程

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

(1)

(2)

Ein+Eg=Est

(3)

其中：ρ為密度；u為速度；p為壓力；ν為運(yùn)動(dòng)黏性系數(shù)；f為外力；t為時(shí)間；Ein為單位時(shí)間控制體積內(nèi)能量?jī)羧肓浚籈g為單位時(shí)間控制體積內(nèi)能量產(chǎn)生量；Est為控制體及內(nèi)能量變化率。

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

1.2 粒子相互作用模型

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

w(|rj-ri|)=

(4)

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

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

(5)

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

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

(6)

(7)

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

1.3 能量守恒方程的求解

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

Aj=ljl0

(8)

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

那么對(duì)于i粒子的控制體積則為：

(9)

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

(10)

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

1.4 算法流程

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

2 導(dǎo)熱問(wèn)題的驗(yàn)證

圖2示出一維半無(wú)限平板導(dǎo)熱問(wèn)題算例模型，其計(jì)算域尺寸為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

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

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

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

a——前期；b——后期圖4 數(shù)值計(jì)算溫度與精確解間的誤差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)自然對(duì)流

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

本文方腔內(nèi)自然對(duì)流的初始條件如下。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

在模擬計(jì)算方腔自然對(duì)流時(shí)，設(shè)定了如圖6的初始條件及邊界條件，隨計(jì)算的進(jìn)行，逐漸達(dá)到穩(wěn)態(tài)過(guò)程。圖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的計(jì)算結(jié)果Fig.9 Calculation result of Ra=105

本文還計(jì)算了Ra分別為103及105時(shí)方腔自然對(duì)流的結(jié)果，如圖8、9所示，當(dāng)Ra較小時(shí)，方腔內(nèi)部的等溫線幾乎沿豎直方向，流動(dòng)的典型特性是在方腔中間出現(xiàn)了1個(gè)大渦。當(dāng)Ra增大時(shí)，方腔內(nèi)部等溫線逐漸變?yōu)樗剑⑶以跓崂浔诿娓浇谋∵吔鐚觾?nèi)仍保持垂直，這說(shuō)明，當(dāng)Ra較小時(shí)，換熱主要是熱傳導(dǎo)為主，隨著Ra的增加，換熱逐漸變?yōu)橛蓪?duì)流換熱為主導(dǎo)地位。Ra增大后，方腔內(nèi)部的大渦分裂成兩個(gè)。將計(jì)算所得的平均Nu與基準(zhǔn)題解進(jìn)行對(duì)比，結(jié)果如圖10所示，可發(fā)現(xiàn)粒子-網(wǎng)格混合方法在計(jì)算方腔內(nèi)自然對(duì)流時(shí)具有較高的精度。本文選擇了方腔對(duì)流達(dá)到穩(wěn)定狀態(tài)后的平均Nu進(jìn)行對(duì)比。為了對(duì)計(jì)算所得的平均Nu進(jìn)行不確定度分析，本文對(duì)穩(wěn)定狀態(tài)后20個(gè)時(shí)刻的平均Nu進(jìn)行分析，可得到Ra為103～106的方腔對(duì)流平均Nu的標(biāo)準(zhǔn)偏差uNu為：

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

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

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

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

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

4 結(jié)論

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