改進(jìn)近場(chǎng)動(dòng)力學(xué)方法在圓柱繞流問(wèn)題中的應(yīng)用

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

(1.武漢理工大學(xué)力學(xué)系，湖北 武漢 430070；2.武漢理工大學(xué)材料復(fù)合新技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430070)

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

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

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

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

1 改進(jìn)近場(chǎng)動(dòng)力學(xué)鍵理論的計(jì)算方法

1. 1 流體本構(gòu)方程

流體中一點(diǎn)的應(yīng)力狀態(tài)由9個(gè)分量構(gòu)成，由矩陣表示流體應(yīng)力狀態(tài)P(表示應(yīng)力矩陣)的張量形式可寫為[19]

pij=-pδij+τij

(1)

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

(2)

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

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

對(duì)于二維流動(dòng)問(wèn)題，上述計(jì)算公式可轉(zhuǎn)化為

(3)

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

假設(shè)平均壓強(qiáng)(以下簡(jiǎn)稱為壓強(qiáng))與熱力學(xué)壓強(qiáng)相等，則有：

(4)

當(dāng)i=j時(shí)，式(1)中τij為附加法向應(yīng)力，表示在流體元上3個(gè)法向應(yīng)力分量偏離平均壓強(qiáng)的部分，該應(yīng)力與流體元的變形率相關(guān)，當(dāng)流體靜止時(shí)該量消失。

根據(jù)Stokes提出的三個(gè)假設(shè)[19]：①流體應(yīng)力與變形率成線性關(guān)系；②流體是各向同性的，應(yīng)力與變形率的關(guān)系與坐標(biāo)系的選取無(wú)關(guān)；③流體靜止時(shí)，法向應(yīng)力等于靜壓強(qiáng)。

可得到牛頓流體的本構(gòu)關(guān)系如下：

(5)

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

μ與λ的關(guān)系式為

(6)

則二維流體本構(gòu)方程可表示為

(7)

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

1. 2 法向鍵力與切向鍵力的推導(dǎo)

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

圖2 粒子x與領(lǐng)域范圍內(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)，假設(shè)粒子x與粒子x'之間的鍵與x軸夾角為θ，沿鍵長(zhǎng)方向單位矢量為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和垂直于鍵方向的切向鍵力τ作用，則基于鍵理論的近場(chǎng)動(dòng)力學(xué)運(yùn)動(dòng)平衡方程[4]為

(9)

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

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

(10)

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

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

(11)

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

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

(12)

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

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

(13)

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

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

(14)

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

(15)

推導(dǎo)過(guò)程中，法向鍵力取的速度鄰域[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),即可得到改進(jìn)的近場(chǎng)動(dòng)力學(xué)鍵理論方法的流體運(yùn)動(dòng)控制方程。

1. 3 邊界條件和初始條件的設(shè)置

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

vxf=0,vyf=0

(17)

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

壁面邊界條件為

vxs=vxf,vys=vyf

(18)

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

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

速度入口邊界條件需要確定來(lái)流速度，可以通過(guò)分量形式或矢量形式來(lái)指定。本文選取流向?yàn)閤方向，垂直于x方向的逆時(shí)針為y方向正方向，則速度入口邊界上的速度應(yīng)為來(lái)流速度v0，即

vxf=v0,vyf=0

(19)

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

2 圓柱繞流問(wèn)題的計(jì)算流程

2. 1 計(jì)算流程

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

圖6 圓柱繞流問(wèn)題的計(jì)算流程圖Fig.6 Flow chart of flow around a cylinder

2. 2 空間離散與時(shí)間積分

改進(jìn)的近場(chǎng)動(dòng)力學(xué)鍵理論方法的流體運(yùn)動(dòng)控制方程計(jì)算的離散格式如下：

(20)

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

3 算例應(yīng)用與分析

3. 1 二維圓柱繞流問(wèn)題的計(jì)算模型

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

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

3. 2 數(shù)值計(jì)算結(jié)果

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

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

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

4 結(jié)論與建議

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

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

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