一維無源對流擴散方程的近場動力學計算格式

2022-02-16 11:46:00王慶李家寶鞠磊薛彥卓賈定睿

哈爾濱工程大學學報 2022年1期

關鍵詞：模型

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

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

對流和擴散是自然環境中常見的現象，包括河流、大氣等污染中污染物的分布、流體流動、熱運動、離子運動等[1]。擴散是微觀上物質分子從高濃度區域向低濃度區域轉移運動的現象，其擴散速率與物質的濃度梯度成正比；對流是宏觀上物質傳遞的現象。對流擴散方程是流體力學中一個基本的運動微分方程，對于很多領域具有重要的應用價值，但是除了部分簡單的情形，大多數問題無法得到精確解，只能利用數值方法進行模擬，因此對流擴散方程數值解法的研究具有重要意義。

從求解微分方程的方法上來說，目前主要的計算方法包括有限差分法[2-4](FDM)、有限元法[5-6](FEM)、有限體積法[7](FVM)，基本思想是對求解區域進行網格劃分，將微分方程離散化，把連續問題化為有限形式的線性代數方程組，求出原問題的近似解[8]。當對流占優時，傳統的數值方法會產生數值震蕩現象[9]，為此需要合適的數值算法消除震蕩。

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

除了擴散方程，采用近場動力學解決對流方程的研究也有了一定成果，Zhao等[22]采用中心差分格式和迎風格式的混合權重模型來計算對流擴散方程，但是原文中只是簡單的采用擴散系數和對流速度絕對值的比值來判斷對流和擴散的占優情況，并未考慮到網格大小的影響，其列舉的對流占優例子實質為擴散占優，因此不能就此說明原文提出的計算格式適用于對流占優情況。

本文的主要目的是探究近場動力學理論求解一維無源項穩態對流擴散方程的可行性及數值上的準確性。在求解無源項對流方程的基礎上，研究其中時間步長與空間步長的關系。同時引用無量綱數Pe考慮對流和擴散的主導地位，確定對流擴散方程的近場動力學求解格式，并與理論值進行比對分析。

1 近場動力學擴散模型

近場動力學研究擴散作用已經相當深入，以無源無匯的濃度擴散方程為例，采用鍵基近場動力學可以表示為：

(1)

式中：fh稱為響應函數，代表物質點x′與x間的相互作用；Hx是點x的近場域范圍；Vx′為點x′的體積；θ′和θ分別代表t時刻x′和x處的濃度。

響應函數通常有3種表示形式[20-21]：

f3=α3(θ′-θ)

(2)

式中α1、α2、α3為不同形式的微擴散系數。

本文選擇第1種響應函數f1為例，對應的近場動力學擴散方程可表示為：

(3)

式中α為x′與x之間的微擴散系數。

為體現鍵長對擴散作用的影響，引入權函數ω。如式(4)所示，權函數ω通常有2種形式，常數形式表示鍵長不影響點與點間的相互作用；線性形式表示點與點間的相互作用隨鍵長增加而減弱。如式(5)所示，將α表示成ω的函數。

(4)

(5)

式中：α0是宏觀的擴散系數；A表示離散后一維物質點的橫截面積；h為二維板的厚度；δ表示物質點的近場域半徑。

2 近場動力學對流模型

假設存在如圖1所示的圓柱形界面[19，22]，有離子沿圓柱體軸向發生對流運動，圓柱側面無能量耗散。上、下界面面積為S,分別表示成S1、S2，上、下界面對應的濃度為θ1、θ2，界面相距為d，則單位時間離子總量變化可以寫為[22]：

圖1 離子沿圓柱軸向對流運動[19，22]

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

(6)

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

假設θa是圓柱體內的離子平均濃度，則單位時間內的離子總量變化可以寫為：

(7)

在近場動力學中認為每個點只與其近場范圍內的點具有相互作用，同樣以近場動力學表示的對流方程，也認為每個點只與其近場范圍內的點發生對流。依照近場動力學鍵基模型，認為作用域內的點與點之間存在相互作用，即為“鍵”，對流作用即是通過“鍵”產生的。

假定在空間域中的一點為x，在該點的作用域Hx中存在一點x′，則兩點間的對流作用為：

(8)

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

通過對x點整個近場域內的物質點進行積分，可以得到：

(9)

假設x點的平均濃度與鍵的平均濃度滿足：

(10)

可以得到：

(11)

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

(12)

3 數值實現

3.1 擴散方程離散

一維擴散方程離散形式為：

(13)

式中：m=δ/Δx，m稱為鄰域系數;Δx為離散后的粒子點間距。

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

(14)

3.2 對流方程離散

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

圖2 一維擴散物質點分布(m=3)

圖3 一維對流物質點分布(對流速度正向，m=3)

(15)

對流方程的離散同樣需要討論xj→xi的情況，不同于擴散方程的中心差分形式，對流計算需要考慮對流方向的影響。在對流的數值計算中，對流項處理不當會導致數值發散。流動方向會影響流動信息，理論上上游和下游會對彼此產生影響，但是兩者的影響程度是不一樣的，在同時包括從上到下的對流和自由擴散時，上游對下游的影響包括同方向的擴散和對流，而下游對上游的影響則取決于擴散和對流兩者的比值，即Peclet數Pe，因為二者是反方向的。考慮到流動方向上的信息傳輸更為明顯，xj→xi可以直接采用迎風格式進行近似：

以右為正方向，對流從左向右時：

(16)

對流從右向左時：

(17)

3.3 時間步長計算

時間積分統一采用向前顯式積分，即：

(18)

因為涉及到對流和擴散2種模型，因此需要分別計算2種模型的收斂時間步長。

擴散的收斂時間步長計算公式[23]：

(19)

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

(20)

同理可以得到：

(21)

對擴散方程的離散格式和時間步長進行驗證，計算模型為一塊有限尺寸的厚板[23]，在其表面施加一個隨時間線性增長的溫度邊界條件。具體模擬參數如表1所示。

表1 溫度擴散計算參數

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

(22)

溫度T的理論解：

(23)

式中：n可以取一個較大的數值，在本文中，取n=100。

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

圖4 一維熱擴散PD計算值與理論值對比

下面以正向對流運動為例對一維對流方程的收斂時間步長進行推導。此處及后文的u和U均為標量，由擴散過程的數值穩定性推導可以類比得到一維對流方程收斂的時間步長估算公式：

(24)

(25)

為了明確PD計算值與理論值的對比結果，將二者的相對誤差定義為：

(26)

式中N為模型內部的物質點數目。

4 對流模型驗證

本節以三角函數為研究對象，驗證近場動力學格式的對流方程在數值應用上的準確性。

一維無源對流方程可以寫為：

(27)

式中U為對流速度，上述方程的理論解可以寫為三角函數的形式：

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

(28)

圖5 近場動力學對流離散模型

初始條件：θ(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。

該三角函數周期為1 s，選定計算時長為半個周期，即為0.5 s。

從圖6可以看出，PD計算出的函數值與理論值較為接近，由式(26)得到的計算誤差為2.60%，在可以接受的范圍內，證實了采用近場動力學求解對流方程的可行性。

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

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

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