王慶， 李家寶， 鞠磊， 薛彥卓， 賈定睿

(哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱 150001)

對(duì)流和擴(kuò)散是自然環(huán)境中常見(jiàn)的現(xiàn)象，包括河流、大氣等污染中污染物的分布、流體流動(dòng)、熱運(yùn)動(dòng)、離子運(yùn)動(dòng)等[1]。擴(kuò)散是微觀上物質(zhì)分子從高濃度區(qū)域向低濃度區(qū)域轉(zhuǎn)移運(yùn)動(dòng)的現(xiàn)象，其擴(kuò)散速率與物質(zhì)的濃度梯度成正比；對(duì)流是宏觀上物質(zhì)傳遞的現(xiàn)象。對(duì)流擴(kuò)散方程是流體力學(xué)中一個(gè)基本的運(yùn)動(dòng)微分方程，對(duì)于很多領(lǐng)域具有重要的應(yīng)用價(jià)值，但是除了部分簡(jiǎn)單的情形，大多數(shù)問(wèn)題無(wú)法得到精確解，只能利用數(shù)值方法進(jìn)行模擬，因此對(duì)流擴(kuò)散方程數(shù)值解法的研究具有重要意義。

從求解微分方程的方法上來(lái)說(shuō)，目前主要的計(jì)算方法包括有限差分法[2-4](FDM)、有限元法[5-6](FEM)、有限體積法[7](FVM)，基本思想是對(duì)求解區(qū)域進(jìn)行網(wǎng)格劃分，將微分方程離散化，把連續(xù)問(wèn)題化為有限形式的線性代數(shù)方程組，求出原問(wèn)題的近似解[8]。當(dāng)對(duì)流占優(yōu)時(shí)，傳統(tǒng)的數(shù)值方法會(huì)產(chǎn)生數(shù)值震蕩現(xiàn)象[9]，為此需要合適的數(shù)值算法消除震蕩。

近場(chǎng)動(dòng)力學(xué)理論(peridynamics, PD)是近些年來(lái)快速發(fā)展的一種非局部理論[10-11]，對(duì)于求解固體裂紋擴(kuò)展的問(wèn)題有非常高的適用性[12-14]。目前也已經(jīng)成功應(yīng)用到了包括熱傳導(dǎo)、電傳導(dǎo)、水分濃度擴(kuò)散等一系列常規(guī)非局部擴(kuò)散問(wèn)題上[15-19]。Bobaru等[18-19]構(gòu)建了鍵基非穩(wěn)態(tài)熱傳導(dǎo)方程，驗(yàn)證了方程的適用性，并對(duì)含裂紋板的熱傳導(dǎo)問(wèn)題進(jìn)行研究。Agwai等[20]推導(dǎo)出了態(tài)基近場(chǎng)動(dòng)力學(xué)熱傳導(dǎo)方程，進(jìn)而得到了簡(jiǎn)化的鍵基近場(chǎng)動(dòng)力學(xué)熱傳導(dǎo)方程，提出了3種核函數(shù)形式及計(jì)算方法，給出了熱傳導(dǎo)方程在近場(chǎng)動(dòng)力學(xué)格式下詳細(xì)的數(shù)值實(shí)現(xiàn)步驟。Chen等[21]研究了不同響應(yīng)函數(shù)下鍵基近場(chǎng)動(dòng)力學(xué)熱傳導(dǎo)方程的計(jì)算參數(shù)收斂性問(wèn)題。

除了擴(kuò)散方程，采用近場(chǎng)動(dòng)力學(xué)解決對(duì)流方程的研究也有了一定成果，Zhao等[22]采用中心差分格式和迎風(fēng)格式的混合權(quán)重模型來(lái)計(jì)算對(duì)流擴(kuò)散方程，但是原文中只是簡(jiǎn)單的采用擴(kuò)散系數(shù)和對(duì)流速度絕對(duì)值的比值來(lái)判斷對(duì)流和擴(kuò)散的占優(yōu)情況，并未考慮到網(wǎng)格大小的影響，其列舉的對(duì)流占優(yōu)例子實(shí)質(zhì)為擴(kuò)散占優(yōu)，因此不能就此說(shuō)明原文提出的計(jì)算格式適用于對(duì)流占優(yōu)情況。

本文的主要目的是探究近場(chǎng)動(dòng)力學(xué)理論求解一維無(wú)源項(xiàng)穩(wěn)態(tài)對(duì)流擴(kuò)散方程的可行性及數(shù)值上的準(zhǔn)確性。在求解無(wú)源項(xiàng)對(duì)流方程的基礎(chǔ)上，研究其中時(shí)間步長(zhǎng)與空間步長(zhǎng)的關(guān)系。同時(shí)引用無(wú)量綱數(shù)Pe考慮對(duì)流和擴(kuò)散的主導(dǎo)地位，確定對(duì)流擴(kuò)散方程的近場(chǎng)動(dòng)力學(xué)求解格式，并與理論值進(jìn)行比對(duì)分析。

1 近場(chǎng)動(dòng)力學(xué)擴(kuò)散模型

近場(chǎng)動(dòng)力學(xué)研究擴(kuò)散作用已經(jīng)相當(dāng)深入，以無(wú)源無(wú)匯的濃度擴(kuò)散方程為例，采用鍵基近場(chǎng)動(dòng)力學(xué)可以表示為：

(1)

式中：fh稱(chēng)為響應(yīng)函數(shù)，代表物質(zhì)點(diǎn)x′與x間的相互作用；Hx是點(diǎn)x的近場(chǎng)域范圍；Vx′為點(diǎn)x′的體積；θ′和θ分別代表t時(shí)刻x′和x處的濃度。

響應(yīng)函數(shù)通常有3種表示形式[20-21]：

f3=α3(θ′-θ)

(2)

式中α1、α2、α3為不同形式的微擴(kuò)散系數(shù)。

本文選擇第1種響應(yīng)函數(shù)f1為例，對(duì)應(yīng)的近場(chǎng)動(dòng)力學(xué)擴(kuò)散方程可表示為：

(3)

式中α為x′與x之間的微擴(kuò)散系數(shù)。

為體現(xiàn)鍵長(zhǎng)對(duì)擴(kuò)散作用的影響，引入權(quán)函數(shù)ω。如式(4)所示，權(quán)函數(shù)ω通常有2種形式，常數(shù)形式表示鍵長(zhǎng)不影響點(diǎn)與點(diǎn)間的相互作用；線性形式表示點(diǎn)與點(diǎn)間的相互作用隨鍵長(zhǎng)增加而減弱。如式(5)所示，將α表示成ω的函數(shù)。

(4)

(5)

式中：α0是宏觀的擴(kuò)散系數(shù)；A表示離散后一維物質(zhì)點(diǎn)的橫截面積；h為二維板的厚度；δ表示物質(zhì)點(diǎn)的近場(chǎng)域半徑。

2 近場(chǎng)動(dòng)力學(xué)對(duì)流模型

假設(shè)存在如圖1所示的圓柱形界面[19，22]，有離子沿圓柱體軸向發(fā)生對(duì)流運(yùn)動(dòng)，圓柱側(cè)面無(wú)能量耗散。上、下界面面積為S,分別表示成S1、S2，上、下界面對(duì)應(yīng)的濃度為θ1、θ2，界面相距為d，則單位時(shí)間離子總量變化可以寫(xiě)為[22]：

圖1 離子沿圓柱軸向?qū)α鬟\(yùn)動(dòng)[19，22]

S(θ1-θ2)U·e(x，x′)

(6)

式中：U是對(duì)流速度矢量；e(x，x′)是沿圓柱體縱軸的單位矢量。

假設(shè)θa是圓柱體內(nèi)的離子平均濃度，則單位時(shí)間內(nèi)的離子總量變化可以寫(xiě)為：

(7)

在近場(chǎng)動(dòng)力學(xué)中認(rèn)為每個(gè)點(diǎn)只與其近場(chǎng)范圍內(nèi)的點(diǎn)具有相互作用，同樣以近場(chǎng)動(dòng)力學(xué)表示的對(duì)流方程，也認(rèn)為每個(gè)點(diǎn)只與其近場(chǎng)范圍內(nèi)的點(diǎn)發(fā)生對(duì)流。依照近場(chǎng)動(dòng)力學(xué)鍵基模型，認(rèn)為作用域內(nèi)的點(diǎn)與點(diǎn)之間存在相互作用，即為“鍵”，對(duì)流作用即是通過(guò)“鍵”產(chǎn)生的。

假定在空間域中的一點(diǎn)為x，在該點(diǎn)的作用域Hx中存在一點(diǎn)x′，則兩點(diǎn)間的對(duì)流作用為：

(8)

式中：θ(x,t)和θ(x′,t)分別代表點(diǎn)x和x′在t時(shí)刻的離子濃度。

通過(guò)對(duì)x點(diǎn)整個(gè)近場(chǎng)域內(nèi)的物質(zhì)點(diǎn)進(jìn)行積分，可以得到：

(9)

假設(shè)x點(diǎn)的平均濃度與鍵的平均濃度滿(mǎn)足：

(10)

可以得到：

(11)

式中：u是點(diǎn)x和x′之間鍵的微速度矢量，滿(mǎn)足u=ωU/VHx[22];VHx表示點(diǎn)x的近場(chǎng)域體積，同樣以式(4)中的權(quán)函數(shù)ω表示鍵長(zhǎng)對(duì)對(duì)流作用的影響：

(12)

3 數(shù)值實(shí)現(xiàn)

3.1 擴(kuò)散方程離散

一維擴(kuò)散方程離散形式為：

(13)

式中：m=δ/Δx，m稱(chēng)為鄰域系數(shù);Δx為離散后的粒子點(diǎn)間距。

xj是xi近場(chǎng)域內(nèi)的一點(diǎn)，需要特別注意的是xj→xi的情況，無(wú)限趨近時(shí)2點(diǎn)可視為一點(diǎn)，雖然這種情況只存在于數(shù)學(xué)理論上，但該項(xiàng)表示的是點(diǎn)xi處濃度的空間變化率[18]，不能忽略。離散的物質(zhì)點(diǎn)具有體積，為了不使得物質(zhì)點(diǎn)相互滲透，可以采用多種高階精度格式來(lái)進(jìn)行近似該項(xiàng)，通常情況下，采用二階中心差分格式計(jì)算即可。由此，式(13)可以表示為：

(14)

3.2 對(duì)流方程離散

圖3中U為宏觀對(duì)流速度矢量。一維對(duì)流方程的離散形式與一維擴(kuò)散方程的離散形式基本一致：

圖2 一維擴(kuò)散物質(zhì)點(diǎn)分布(m=3)

圖3 一維對(duì)流物質(zhì)點(diǎn)分布(對(duì)流速度正向，m=3)

(15)

對(duì)流方程的離散同樣需要討論xj→xi的情況，不同于擴(kuò)散方程的中心差分形式，對(duì)流計(jì)算需要考慮對(duì)流方向的影響。在對(duì)流的數(shù)值計(jì)算中，對(duì)流項(xiàng)處理不當(dāng)會(huì)導(dǎo)致數(shù)值發(fā)散。流動(dòng)方向會(huì)影響流動(dòng)信息，理論上上游和下游會(huì)對(duì)彼此產(chǎn)生影響，但是兩者的影響程度是不一樣的，在同時(shí)包括從上到下的對(duì)流和自由擴(kuò)散時(shí)，上游對(duì)下游的影響包括同方向的擴(kuò)散和對(duì)流，而下游對(duì)上游的影響則取決于擴(kuò)散和對(duì)流兩者的比值，即Peclet數(shù)Pe，因?yàn)槎呤欠捶较虻??？紤]到流動(dòng)方向上的信息傳輸更為明顯，xj→xi可以直接采用迎風(fēng)格式進(jìn)行近似：

以右為正方向，對(duì)流從左向右時(shí)：

(16)

對(duì)流從右向左時(shí)：

(17)

3.3 時(shí)間步長(zhǎng)計(jì)算

時(shí)間積分統(tǒng)一采用向前顯式積分，即：

(18)

因?yàn)樯婕暗綄?duì)流和擴(kuò)散2種模型，因此需要分別計(jì)算2種模型的收斂時(shí)間步長(zhǎng)。

擴(kuò)散的收斂時(shí)間步長(zhǎng)計(jì)算公式[23]：

(19)

式中ξij表示xi和xj的位置矢量，顯然‖ξij‖最小為Δx。將一維擴(kuò)散系數(shù)的常數(shù)權(quán)函數(shù)形式代入式(19)，放縮后可取等號(hào)得到：

(20)

同理可以得到：

(21)

對(duì)擴(kuò)散方程的離散格式和時(shí)間步長(zhǎng)進(jìn)行驗(yàn)證，計(jì)算模型為一塊有限尺寸的厚板[23]，在其表面施加一個(gè)隨時(shí)間線性增長(zhǎng)的溫度邊界條件。具體模擬參數(shù)如表1所示。

表1 溫度擴(kuò)散計(jì)算參數(shù)

板的初始溫度為0 ℃，邊界條件為:

(22)

溫度T的理論解：

(23)

式中：n可以取一個(gè)較大的數(shù)值，在本文中，取n=100。

從圖4可以看出，理論值和PD計(jì)算值的對(duì)比結(jié)果吻合較好。在文獻(xiàn)[23]中采用的時(shí)間步長(zhǎng)為10-6s，由式(20)計(jì)算出的范圍為：Δt1<3.5×10-7s，為了方便計(jì)算，本文取為10-7s，可以看出由于經(jīng)過(guò)放縮過(guò)程，由式(20)得到的時(shí)間步長(zhǎng)取值范圍比式(19)計(jì)算出的要小，但是其結(jié)果也是必定滿(mǎn)足收斂條件的。

圖4 一維熱擴(kuò)散PD計(jì)算值與理論值對(duì)比

下面以正向?qū)α鬟\(yùn)動(dòng)為例對(duì)一維對(duì)流方程的收斂時(shí)間步長(zhǎng)進(jìn)行推導(dǎo)。此處及后文的u和U均為標(biāo)量，由擴(kuò)散過(guò)程的數(shù)值穩(wěn)定性推導(dǎo)可以類(lèi)比得到一維對(duì)流方程收斂的時(shí)間步長(zhǎng)估算公式：

(24)

(25)

為了明確PD計(jì)算值與理論值的對(duì)比結(jié)果，將二者的相對(duì)誤差定義為：

(26)

式中N為模型內(nèi)部的物質(zhì)點(diǎn)數(shù)目。

4 對(duì)流模型驗(yàn)證

本節(jié)以三角函數(shù)為研究對(duì)象，驗(yàn)證近場(chǎng)動(dòng)力學(xué)格式的對(duì)流方程在數(shù)值應(yīng)用上的準(zhǔn)確性。

一維無(wú)源對(duì)流方程可以寫(xiě)為：

(27)

式中U為對(duì)流速度，上述方程的理論解可以寫(xiě)為三角函數(shù)的形式：

θ(x,t)=0.5sin(2π(x-Ut))

(28)

圖5 近場(chǎng)動(dòng)力學(xué)對(duì)流離散模型

初始條件：θ(x,0)=0.5sin(2πx)；

左右虛擬邊界條件：θ(x,t)=0.5sin(2π(x-Ut))；

變量取值：x∈[0,2]m，t∈[0,0.5]s，U=1 m/s。

該三角函數(shù)周期為1 s，選定計(jì)算時(shí)長(zhǎng)為半個(gè)周期，即為0.5 s。

從圖6可以看出，PD計(jì)算出的函數(shù)值與理論值較為接近，由式(26)得到的計(jì)算誤差為2.60%，在可以接受的范圍內(nèi)，證實(shí)了采用近場(chǎng)動(dòng)力學(xué)求解對(duì)流方程的可行性。

圖6 函數(shù)初始圖像及t=0.5 s時(shí)的函數(shù)圖像(Δx=L/100， m=3)

近場(chǎng)域范圍系數(shù)通常取m=3[12]，但是由于引入了安全系數(shù)，因此為了得到安全系數(shù)取值對(duì)計(jì)算結(jié)果誤差的影響，以Δx=L/100，L/200，L/500為例，分別取m為3、4、5，綜合考慮安全系數(shù)、離散格式和近場(chǎng)域范圍系數(shù)等因素對(duì)計(jì)算結(jié)果的影響，得到以下結(jié)果。

從表2～4中能夠看出:1)隨著離散物質(zhì)點(diǎn)數(shù)的增加，同一鄰域系數(shù)m及安全系數(shù)對(duì)應(yīng)的計(jì)算誤差逐漸減小；2)相同鄰域系數(shù)m及離散格式下對(duì)應(yīng)的誤差呈現(xiàn)出先減小后增大的趨勢(shì)，當(dāng)鄰域系數(shù)m與安全系數(shù)的倒數(shù)滿(mǎn)足1/ζ=2m時(shí)，可以得到最小的計(jì)算誤差；3)整體趨勢(shì)上，隨著鄰域系數(shù)m的增大，相同條件下的計(jì)算誤差增大，因此在計(jì)算時(shí)可以直接取m=3，當(dāng)物質(zhì)點(diǎn)間隔足夠小時(shí)，如Δx=L/500，在誤差可接受的范圍內(nèi)，為了方便計(jì)算，m也可以取其他值。

表2 m=3時(shí)不同的安全系數(shù)對(duì)應(yīng)的計(jì)算誤差

表3 m=4時(shí)不同的安全系數(shù)對(duì)應(yīng)的計(jì)算誤差

表4 m=5時(shí)不同的安全系數(shù)對(duì)應(yīng)的計(jì)算誤差

5 對(duì)流擴(kuò)散模型驗(yàn)證

將采用近場(chǎng)動(dòng)力學(xué)理論描述的擴(kuò)散和對(duì)流方程合并，可以得到無(wú)源項(xiàng)的一維對(duì)流擴(kuò)散方程：

(29)

式中θ′和θ分別代表t時(shí)刻x′和x處的濃度。將一維對(duì)流方程的離散式(14)和一維擴(kuò)散方程的離散式(16)或式(17)結(jié)合起來(lái)即可得到式(29)的離散形式。

本節(jié)仍然以三角函數(shù)為例對(duì)基于近場(chǎng)動(dòng)力學(xué)的一維對(duì)流擴(kuò)散方程進(jìn)行驗(yàn)證：

(30)

其理論解為：

θ(x,t)=e-4α0π2t0.5sin(2π(x-Ut))

(31)

從理論解的形式可以看出:在一個(gè)周期的時(shí)間內(nèi)，只發(fā)生擴(kuò)散時(shí)函數(shù)幅值會(huì)減小；只發(fā)生對(duì)流時(shí)函數(shù)的峰值位置會(huì)向?qū)α鞯恼较蛞苿?dòng)；若二者同時(shí)發(fā)生，當(dāng)其中一者占優(yōu)時(shí)，另一者的現(xiàn)象會(huì)較弱；而二者均不占優(yōu)時(shí)，雙方的現(xiàn)象會(huì)較明顯。

5.1 擴(kuò)散占主導(dǎo)地位

從圖7可以看出，隨時(shí)間的變化，函數(shù)的波峰、波谷的位置沒(méi)有發(fā)生明顯的改變，而函數(shù)幅值有了明顯下降，這說(shuō)明計(jì)算中主要發(fā)生了擴(kuò)散現(xiàn)象，導(dǎo)致幅值迅速減小，這個(gè)現(xiàn)象與計(jì)算時(shí)的假設(shè)一致。而且從理論值與計(jì)算值的比對(duì)可以看出計(jì)算誤差比較小，計(jì)算結(jié)束時(shí)的相對(duì)誤差為2.43%。

圖7 擴(kuò)散占優(yōu)時(shí)不同時(shí)刻的函數(shù)值對(duì)比

5.2 對(duì)流占主導(dǎo)地位

從圖8可以看出隨時(shí)間變化，函數(shù)的幅值未發(fā)生明顯的變化，但是波峰、波谷的位置產(chǎn)生了顯著的位移，也即是對(duì)流現(xiàn)象明顯，與假設(shè)的情況一致。計(jì)算結(jié)束時(shí)，相對(duì)誤差為1.653%。

圖8 對(duì)流占優(yōu)時(shí)不同時(shí)刻的函數(shù)值對(duì)比

5.3 對(duì)流和擴(kuò)散均不占優(yōu)

在Pe=1時(shí)，從圖9看出，隨時(shí)間的變化，函數(shù)的波峰、波谷均產(chǎn)生了明顯的正向位移，同時(shí)幅值減小。這2種現(xiàn)象都可以從圖上明顯觀察出來(lái)，也就證明了對(duì)流和擴(kuò)散實(shí)際上均不占優(yōu)，計(jì)算結(jié)束時(shí)的相對(duì)誤差為1.478%。

圖9 對(duì)流擴(kuò)散均不占優(yōu)時(shí)不同時(shí)刻的函數(shù)值對(duì)比

通過(guò)以上算例，證實(shí)了無(wú)量綱數(shù)Pe可以用于判斷對(duì)流和擴(kuò)散現(xiàn)象是否占優(yōu)，同時(shí)也驗(yàn)證了本文構(gòu)建的近場(chǎng)動(dòng)力學(xué)對(duì)流擴(kuò)散模型的準(zhǔn)確性。

6 結(jié)論

1)本文成功構(gòu)建了一維無(wú)源近場(chǎng)動(dòng)力學(xué)對(duì)流擴(kuò)散模型，并在數(shù)值上驗(yàn)證了該模型的有效性和準(zhǔn)確性；

2)當(dāng)采用一階迎風(fēng)格式計(jì)算一維無(wú)源近場(chǎng)動(dòng)力學(xué)對(duì)流模型時(shí)，時(shí)間步長(zhǎng)的安全系數(shù)ζ與近場(chǎng)范圍系數(shù)m滿(mǎn)足1/ζ=2m時(shí)，模型的相對(duì)誤差可以取到最小值；

3)無(wú)量綱數(shù)Pe可以用于判斷對(duì)流和擴(kuò)散的占優(yōu)性，推導(dǎo)的對(duì)流和擴(kuò)散的收斂時(shí)間步長(zhǎng)公式可用于確定同時(shí)滿(mǎn)足二者要求的時(shí)間步長(zhǎng)，最終的數(shù)值結(jié)果滿(mǎn)足預(yù)期結(jié)果；

文中引入的無(wú)量綱數(shù)Pe和物質(zhì)點(diǎn)間距有關(guān)，因此，離散格式對(duì)于Pe的影響以及單元相關(guān)性還需要進(jìn)行進(jìn)一步的研究。