朱紅玉，劉立勝，2，劉齊文，賴 欣

(1.武漢理工大學力學系，湖北 武漢 430070；2.武漢理工大學材料復合新技術國家重點實驗室，湖北 武漢 430070)

流體工程問題廣泛應用于原油勘探、水文地質(zhì)學、多孔介質(zhì)滲流[1]等領域。過去50年,流體力學隨著計算機的發(fā)展迅速崛起，實驗方法較難滿足流體真實的運動狀態(tài)，且效率低、成本高，相比于風洞實驗方法，數(shù)值模擬方法可以更高效地研究不同因素對流體外流場性能的影響[2]。因此，有效發(fā)展數(shù)值方法研究流體的運動狀態(tài)意義重大。

目前描述流體連續(xù)介質(zhì)運動有兩種方法[3]：歐拉方法和拉格朗日方法。其中，歐拉方法將物理量表示成與流體粒子空間位置、時間相關的函數(shù)；拉格朗日方法跟蹤流體粒子，記錄物理量與流體粒子空間位置、時間的關系。本文采用拉格朗日方法跟蹤流體粒子，對流體的運動狀態(tài)進行描述[4]。近場動力學方法是基于流體連續(xù)介質(zhì)力學的一種非局域方法，由Silling[5]提出，不同于擴展有限元方法等需額外添加判據(jù)，近場動力學方法的流體運動控制方程是積分形式的，沒有空間導數(shù)，克服了網(wǎng)格依賴性、非物質(zhì)奇異性、失效準則選取[7]等困難，而近場動力學鍵理論與經(jīng)典理論的區(qū)別在于有限距離上的鍵與接觸力[8]。

近場動力學方法過去十幾年發(fā)展迅速，已成功運用于裂紋擴展[9]、熱傳導[10]、燒蝕[11]等領域。如廖洋等[12]采用近場動力學強度折減法分析了邊坡的穩(wěn)定性問題;Oterkus等[13]采用近場動力學方法，在不引用其他失效準則情況下研究了滲流裂紋構型中的壓力場與裂紋的擴散；Katiyar等[14]采用常規(guī)態(tài)理論研究了微壓縮單相對流流體在均勻多孔介質(zhì)中的流動；Yang等[15]采用浸入邊界法研究了流體結構與多相黏性流作用;王瀚霖等[16]提出的預測模型考慮了摩擦阻力與流體的可壓縮性對天然氣管道壓力變化的影響，且通過對比驗證了該模型的預測結果較為準確;周學君等[17]將近場動力學中描述流體粒子之間相互作用的本構力引入光滑粒子流體動力學方法中用于處理固壁邊界，能夠體現(xiàn)真實的邊界作用，很好地解決了流體粒子非物理穿透邊界問題；Cundall等[18]在研究巖石變形時，考慮粒子流中顆粒間的相互作用，分別表示了顆粒上的法向力和切向力。

基于上述研究，本文在近場動力學方法中“法向鍵力”的基礎上引入“切向鍵力”，共同描述流體粒子的受力狀態(tài)，這種方法相比傳統(tǒng)的近場動力學鍵理論方法，能更有效地描述流體粒子間的“切向黏性力”作用，根據(jù)均質(zhì)不可壓縮牛頓流體的運動特性，考慮粒子鄰域范圍內(nèi)粒子間的相互作用力與壓力、黏性力的關系，推導出基于拉格朗日描述的近場動力學形式的流體運動控制方程。采用Fortran語言編程利用該方法計算相應的算例，并與有限體積法的模擬結果進行了對比，從而驗證了本文提出的改進近場動力學鍵理論方法在研究流體問題中的可行性和有效性，促進了近場動力學在流體問題領域中的應用研究。

1 改進近場動力學鍵理論的計算方法

1. 1 流體本構方程

流體中一點的應力狀態(tài)由9個分量構成，由矩陣表示流體應力狀態(tài)P(表示應力矩陣)的張量形式可寫為[19]

pij=-pδij+τij

(1)

式中:pij為應力張量(i、j分別取x、y、z);p為壓強;δij為克羅內(nèi)克δ(當i=j時，δij=1;當i≠j時，δij=0);τij為偏應力張量。

(2)

式中:4πa2表示半徑為a的球的表面積;pnn為壓力積分項;dA為面積微元;pxx、pyy、pzz分別表示x、y、z方向的壓應力。

圖1 流體元變形前、后壓力示意圖Fig.1 Diagrams of the pressure of undeformed and deformed fluid element

對于二維流動問題，上述計算公式可轉化為

(3)

式中：2πa表示半徑為a的二維流體元的周長;dS為弧長微元。

假設平均壓強(以下簡稱為壓強)與熱力學壓強相等，則有：

(4)

當i=j時，式(1)中τij為附加法向應力，表示在流體元上3個法向應力分量偏離平均壓強的部分，該應力與流體元的變形率相關，當流體靜止時該量消失。

根據(jù)Stokes提出的三個假設[19]：①流體應力與變形率成線性關系；②流體是各向同性的，應力與變形率的關系與坐標系的選取無關；③流體靜止時，法向應力等于靜壓強。

可得到牛頓流體的本構關系如下：

(5)

式中:vx、vy分別表示沿x、y方向的速度;μ為分子黏性系數(shù)；λ為第二黏性系數(shù)。

μ與λ的關系式為

(6)

則二維流體本構方程可表示為

(7)

式中:σ為應力矩陣;I為單位矩陣;D為變形率張量。

1. 2 法向鍵力與切向鍵力的推導

本文基于近場動力學鍵理論思想[5]，并考慮了粒子鄰域δ范圍內(nèi)粒子間的相互影響，考慮粒子間沿鍵長方向的法向鍵力作用和沿鍵切線方向的切向鍵力作用。在流體運動中，可認為壓強與沿鍵長方向速度分量影響法向鍵力，沿鍵切向速度分量影響切向鍵力[19]。粒子x與鄰域范圍內(nèi)粒子x'之間的法向鍵力和切向鍵力示意圖見圖2。

圖2 粒子x與領域范圍內(nèi)粒子x'之間的法向鍵力和切向鍵力示意圖Fig.2 Schematic diagram of normal bonding force and tangential bonding force in peridynamics between particle x and x' in the neighborhood

在粒子鄰域范圍內(nèi)，假設粒子x與粒子x'之間的鍵與x軸夾角為θ，沿鍵長方向單位矢量為e1，沿鍵切線方向單位矢量為e2，粒子鄰域半徑為δ，見圖3。其中：

e1=cosθi+sinθj

e2=sinθi-cosθj

(8)

式中：i、j分別表示x、y方向的單位矢量。

圖3 粒子鍵力作用鄰域Fig.3 Neighborhood of particles’ bonding force

在鄰域范圍內(nèi)粒子x與粒子x'之間的作用力有沿著鍵方向的法向鍵力f和垂直于鍵方向的切向鍵力τ作用，則基于鍵理論的近場動力學運動平衡方程[4]為

(9)

式中:ρ為粒子的材料密度;u為位移矢量;ü(x,t)為粒子x在t時刻位移的二階導數(shù)，表示加速度矢量;x為粒子x的初始坐標;x′為粒子x′的初始坐標；f為點對相互作用函數(shù);τ為切向作用力。

將公式(9)沿著坐標軸分別分解至x和y方向，則有：

(10)

式中：üx、üy分別表示粒子沿x、y方向的加速度；dAx表示粒子x鄰域范圍內(nèi)面積積分微元。

由均質(zhì)不可壓縮牛頓流體Navier-Stokes方程中的動量方程：

(11)

式中:x方向加速度方程中，壓力與速度沿x方向變化產(chǎn)生的加速度由法向鍵力引起，壓力與速度沿y方向變化產(chǎn)生的加速度由切向鍵力引起。

同理，可分析y方向加速度方程。為了區(qū)分方向，下面用粒子i代表粒子x，用粒子j代表粒子x′，根據(jù)近場動力學思想[10]，則粒子i與粒子j間的“法向鍵力”可表示為

(12)

式中:p為壓力;v為速度矢量，ζij為粒子間相對位置矢量。

同理，粒子i與粒子j間的“切向鍵力”可表示為

(13)

式中：c1、c2、c3為鍵參數(shù)。

將公式(12)和(13)代入公式(10)中，并沿x方向分解，則x方向的流體運動控制方程為

(14)

采用近場動力學鍵理論方法分析的某粒子的加速度應與傳統(tǒng)方法中粒子加速度的表達式相等[10]，則有：

(15)

推導過程中，法向鍵力取的速度鄰域[20]見圖4，切向鍵力取的速度鄰域見圖5。

圖4 法向鍵力速度鄰域Fig.4 Neighborhood of normal bonding force by velocity effects

圖5 切向鍵力速度鄰域Fig.5 Neighborhood of tangential bonding force by velocity effects

可求得鍵參數(shù)c1、c2、c3為

(16)

將公式(16)代入公式(14),即可得到改進的近場動力學鍵理論方法的流體運動控制方程。

1. 3 邊界條件和初始條件的設置

本文選取無限域圓柱繞流問題進行研究，因此邊界條件涉及到圓柱表面和壁面的處理。圓柱表面定義為無滑移壁面條件，邊界上流體粒子的速度等于固體粒子的速度，由于本文研究的是靜止圓柱，所以在圓柱表面邊界Γf上，有：

vxf=0,vyf=0

(17)

式中：vxf、vyf分別表示流體粒子在x、y方向的速度。

壁面邊界條件為

vxs=vxf,vys=vyf

(18)

式中：vxs、vys分別表示固體粒子在x、y方向的速度。

公式(18)則表示壁面固體粒子的速度與其相鄰流體粒子的速度相同。

速度入口邊界條件需要確定來流速度，可以通過分量形式或矢量形式來指定。本文選取流向為x方向，垂直于x方向的逆時針為y方向正方向，則速度入口邊界上的速度應為來流速度v0，即

vxf=v0,vyf=0

(19)

出流邊界上的零擴散通量條件是指出流邊界面上的流動變量均由流場內(nèi)部變量值外插得出。本文由于求解二維均質(zhì)不可壓縮流體控制方程，因此使用了與來流速度相同的均勻流場對流場做了初始化處理。

2 圓柱繞流問題的計算流程

2. 1 計算流程

根據(jù)上述改進的近場動力學鍵理論方法推導出的流體運動控制方程，可根據(jù)速度場計算出壓力場和粒子加速度，并對粒子速度、位置進行更新。圓柱繞流問題的計算流程見圖6。

圖6 圓柱繞流問題的計算流程圖Fig.6 Flow chart of flow around a cylinder

2. 2 空間離散與時間積分

改進的近場動力學鍵理論方法的流體運動控制方程計算的離散格式如下：

(20)

對于某粒子i，假設其鄰域某一粒子為j，則根據(jù)方程式(20),可計算得到粒子當前時刻的加速度，再根據(jù)Velocity-Verlet積分方法更新粒子速度、位移，從而完成流體狀態(tài)的更新。

3 算例應用與分析

3. 1 二維圓柱繞流問題的計算模型

本文選取雷諾數(shù)約為10，均勻來流的二維圓柱繞流模型，流場的尺寸和邊界條件見圖7。圓柱直徑D為0.01 m，選取計算域?qū)挾葹?.05 m，圓柱前方流場和后方流場分別為0.06 m和0.04 m，入口處流入速度vx= 0.001 m/s、vy=0.0 m/s，流場兩側壁處vy=0.0，圓柱壁面處vx=vy=0.0 m/s。流場計算域離散為27 000個粒子，計算時間步長取0.000 25 s，計算時間步為200 000步。計算模型采用Fortran語言編程實現(xiàn)。

圖7 流場尺寸和邊界條件Fig.7 Geometry and boundary conditions of flow field

3. 2 數(shù)值計算結果

本文利用改進的近場動力學鍵理論方法和有限體積法(商業(yè)軟件Fluent)模擬計算得到圓柱繞流算例瞬態(tài)1.25 s、50 s時在x、y方向的速度分布云圖，詳見圖8至圖11。

圖8 1.25 s時圓柱繞流算例在x方向的速度分布云圖Fig.8 X-direction velocity distribution of flow around a cylinder based on the improved bondbased peridynamics and the finite volume method at 1.25 s

圖9 50 s時圓柱繞流算例在x方向的速度分布云圖Fig.9 X-direction velocity distribution of flow around a cylinder based on the improved bondbased peridynamics and the finite volume method at 50 s

圖10 1.25 s時圓柱繞流算例在y方向的速度分布云圖Fig.10 Y-direction velocity distribution of flow around a cylinder based on the improved bondbased peridynamics and the finite volume method at 1.25 s

圖11 50 s時圓柱繞流算例在y方向的速度分布云圖Fig.11 Y-direction velocity distribution of flowaround a cylinder based on the improved bondbased peridynamics and the finite volume method at 50 s

由圖8至圖11可見，本文提出的改進近場動力學鍵理論方法可在相同的初始條件、邊界條件下初步計算流體的速度場，與有限體積法的模擬結果較為接近，但在某些位置存在一些差異，產(chǎn)生這種差異的原因可能是：圓柱壁面處流體粒子與圓柱的作用，有限體積法對網(wǎng)格的依賴性較強，而改進的近場動力學鍵理論方法直接考慮了鄰域內(nèi)粒子間的相互作用，且與鄰域大小、粒子間距的選取有很大的關系，因而影響了數(shù)值模擬計算結果的精度。

當雷諾數(shù)為10時，采用本文提出的改進近場動力學鍵理論方法計算無限域中圓柱繞流算例的耗時大于商業(yè)軟件Fluent，但本文提出的改進方法與有限體積法的計算結果匹配較好，主要表現(xiàn)為：本文提出的改進近場動力學鍵理論方法與商業(yè)軟件Fluent的模擬結果相比，在x方向速度值大小和分布位置相近，終止時刻速度最大值誤差約為6.4%，總體分布相似，且分布范圍一致；在y方向速度分布結果雖不如x方向的精度高，但總體分布仍相似，且分布范圍和趨勢一致，終止時刻速度最大值誤差約為22.6%。

4 結論與建議

本文在傳統(tǒng)近場動力學鍵理論方法的基礎上，提出了一種改進的近場動力學鍵理論方法，并對給定條件下圓柱繞流問題算例進行了數(shù)值模擬分析，得到如下主要結論：

(1) 研究均質(zhì)不可壓縮黏性流體平面運動時，在傳統(tǒng)近場動力學鍵理論方法基礎上，考慮了“法向鍵力”作用，并引入“切向鍵力”作用，通過鍵力與流體本構之間的聯(lián)系，分別推導了“法向鍵力”與“切向鍵力”的鍵參數(shù)，并推導出基于拉格朗日描述的近場動力學形式的流體運動控制方程。

(2) 相比傳統(tǒng)的近場動力學鍵理論方法，本文提出的改進方法能夠考慮流體粒子間黏性力的影響，初步反映了流場中黏性剪切的作用。采用Fortran語言編程，利用改進的近場動力學鍵理論方法模擬計算了低雷諾數(shù)時無限域圓柱繞流問題算例，并與有限體積法(商業(yè)軟件Fluent)的模擬結果進行了對比分析，結果表明：本文提出的改進方法與Fluent軟件的模擬結果在x、y方向速度場的分布、變化趨勢基本一致，且終止時刻速度最大值誤差分別為6.4%、22.6%，能夠初步實現(xiàn)對給定條件下流體運動的描述，從而驗證了本文提出的改進近場動力學鍵理論方法在流體運動狀態(tài)捕捉中的可行性，對近場動力學方法在流體方面的應用與探究作出了一定的貢獻。在后續(xù)研究中可考慮引入圓柱固體材料參數(shù)，著手處理流體粒子與固體粒子之間的強耦合作用，分析流場對固體部分的作用，從而為實際工程中圓柱部分材料的強度設計、流場中圓柱構件的排布等提供參考。