陳榮前 聶德明

(中國(guó)計(jì)量大學(xué)計(jì)量測(cè)試工程學(xué)院，杭州310018)

橢圓顆粒在剪切流中旋轉(zhuǎn)特性的數(shù)值研究1)

陳榮前 聶德明2)

(中國(guó)計(jì)量大學(xué)計(jì)量測(cè)試工程學(xué)院，杭州310018)

研究顆粒在流體剪切作用下的運(yùn)動(dòng)特性是理解和預(yù)測(cè)顆粒懸浮流流動(dòng)行為的關(guān)鍵.當(dāng)流體的慣性不能忽略時(shí)，顆粒的運(yùn)動(dòng)往往變得非常復(fù)雜.本文采用格子Boltzmann方法對(duì)中等雷諾數(shù)下橢圓顆粒在剪切流中的旋轉(zhuǎn)運(yùn)動(dòng)進(jìn)行了模擬.首先,研究了雷諾數(shù)(0＜Re≤170)的影響,結(jié)果表明當(dāng)雷諾數(shù)低于臨界值時(shí)，顆粒以周期性的方式旋轉(zhuǎn)，角速度最小時(shí)對(duì)應(yīng)的長(zhǎng)軸方向隨著雷諾數(shù)的增大而逐漸遠(yuǎn)離水平方向，而且這一傾角與雷諾數(shù)呈分段線性關(guān)系；當(dāng)雷諾數(shù)大于臨界值時(shí)，橢圓形顆粒最終保持靜止?fàn)顟B(tài)，且靜止時(shí)的轉(zhuǎn)角與雷諾數(shù)呈冪函數(shù)關(guān)系，雷諾數(shù)越大，轉(zhuǎn)角越小，橢圓的長(zhǎng)軸越遠(yuǎn)離水平位置.其次,研究了橢圓顆粒的長(zhǎng)短軸之比α(1≤α≤10)的影響，結(jié)果表明顆粒旋轉(zhuǎn)的周期與α呈冪函數(shù)關(guān)系，α越大，顆粒旋轉(zhuǎn)周期越小.此外，當(dāng)α超過臨界值時(shí)，顆粒也在水平位置附近保持靜止?fàn)顟B(tài)，此時(shí)的轉(zhuǎn)角與α也呈冪函數(shù)關(guān)系，α越大，轉(zhuǎn)角越小.研究還發(fā)現(xiàn)，當(dāng)雷諾數(shù)較大時(shí)橢圓顆粒在旋轉(zhuǎn)過程中會(huì)產(chǎn)生過沖現(xiàn)象.

格子Boltzmann方法，剪切流，橢圓顆粒，旋轉(zhuǎn)

引言

由于目前對(duì)顆粒與流體相互作用的機(jī)理研究尚不完善，因此相關(guān)的基礎(chǔ)問題，如顆粒在剪切流體中的運(yùn)動(dòng)以及顆粒的沉降等兩相流問題，一直以來是學(xué)者的研究熱點(diǎn)[1-11].然而，目前許多針對(duì)懸浮顆粒運(yùn)動(dòng)的研究是在低雷諾數(shù)下進(jìn)行的.例如，Hwang等[12]在忽略慣性的前提下采用有限元方法模擬了懸浮粒子在黏彈性簡(jiǎn)單剪切流中的運(yùn)動(dòng)，發(fā)現(xiàn)了雙顆粒間的 KTT(kissing-tumbling-tumbling)現(xiàn)象.同樣，Choi等[13]也采用有限元方法研究了不同初始間距下雙顆粒在庫埃特流中的旋轉(zhuǎn)運(yùn)動(dòng).Lundell等[14]在蠕動(dòng)條件下模擬了橢球顆粒在剪切流中的運(yùn)動(dòng)，他們研究并討論了顆粒的運(yùn)動(dòng)軌跡與顆粒慣性的關(guān)系.Pasquino等[15]通過實(shí)驗(yàn)和直接數(shù)值模擬研究了顆粒在剪切黏彈性流中形成串列結(jié)構(gòu)的現(xiàn)象.

另一方面，如果考慮流體的慣性，則流體自身的運(yùn)動(dòng)及存在于流體中的顆粒的運(yùn)動(dòng)必然會(huì)產(chǎn)生變化.流場(chǎng)的非線性效應(yīng)使得顆粒的運(yùn)動(dòng)特性變得復(fù)雜.例如，Mikulencak等[16]研究了剪切流中圓形和球形顆粒的旋轉(zhuǎn)特性，他們發(fā)現(xiàn)當(dāng)雷諾數(shù)逐漸增大時(shí)，顆粒周圍閉合的流線結(jié)構(gòu)很快消失，而Subramanian等[17-18]進(jìn)一步詳細(xì)研究了這種流線結(jié)構(gòu)的改變對(duì)流場(chǎng)傳熱傳質(zhì)的影響.另外，Yacoubi等[19]采用浸入界面方法對(duì)多顆粒在牛頓流體中的沉降進(jìn)行直接數(shù)值模擬，其設(shè)定的雷諾數(shù)固定為200.研究發(fā)現(xiàn)當(dāng)流場(chǎng)中顆粒數(shù)目為奇數(shù)時(shí)顆粒整體分布呈 “凸”形，而當(dāng)顆粒數(shù)目為偶數(shù)時(shí)顆粒整體分布呈 “凹”形，且處于兩端的顆粒最容易出現(xiàn) DKT(draftingkissing-tumbling)現(xiàn)象.進(jìn)一步地，Nie等[20]采用格子Boltzmann--虛擬域方法模擬了雷諾數(shù)為248的類似問題，研究表明顆粒的初始間距對(duì)DKT有顯著的影響.當(dāng)顆粒間距減小時(shí)，不再是兩端的顆粒發(fā)生DKT現(xiàn)象，而是靠近中心的兩個(gè)顆粒.可見，即使對(duì)于最簡(jiǎn)單的流場(chǎng)條件，流體的慣性也會(huì)使得顆粒的運(yùn)動(dòng)變得復(fù)雜而豐富.此外，如果顆粒不是各向同性的圓形或球形，則情況可能更為復(fù)雜.最近，Nie等[21]對(duì)于沿軸向分布的顆粒沉降特性進(jìn)行了研究，他們發(fā)現(xiàn)了在不同雷諾數(shù)下顆粒的分組沉降的行為.Ding等[22]采用格子Boltzmann方法對(duì)橢圓形顆粒在簡(jiǎn)單剪切流中運(yùn)動(dòng)進(jìn)行直接數(shù)值模擬，研究發(fā)現(xiàn)橢圓形顆粒的旋轉(zhuǎn)具有周期性，也就是說顆粒的旋轉(zhuǎn)速度時(shí)快時(shí)慢，不是一個(gè)固定的值，這與圓形和球形的結(jié)果不同.然而，Ding等[22]沒有對(duì)橢圓的長(zhǎng)/短軸之比的影響進(jìn)行研究，很顯然，這個(gè)比值對(duì)顆粒的旋轉(zhuǎn)特性也同樣有顯著的影響.此外，研究表明無論是橢圓[22]還是橢球體[23]，當(dāng)雷諾數(shù)超過某一臨界數(shù)時(shí)，顆?？赡懿辉俎D(zhuǎn)動(dòng)，而是固定在剪切流場(chǎng)中，但此方面的研究缺乏定量的結(jié)果，如顆粒此時(shí)的方向與雷諾數(shù)及長(zhǎng)短軸之比的關(guān)系等.Je ff ery[1]在忽略流體慣性的前提下對(duì)此問題有理論解，表明當(dāng)橢圓處于水平位置時(shí)候角速度最小，但不為零.由于Je ff ery[1]從Stokes方程出發(fā)，因此無法預(yù)測(cè)雷諾數(shù)較大時(shí)橢圓顆粒的旋轉(zhuǎn)特性，也就無法預(yù)測(cè)橢圓靜止時(shí)所處的位置.本文的研究發(fā)現(xiàn)剪切流場(chǎng)中橢圓靜止時(shí)的傾角接近水平位置，但沒有到達(dá)水平位置，且傾角與雷諾數(shù)和橢圓的長(zhǎng)短軸之比均有關(guān)，存在著定量的關(guān)系.此外，橢圓顆粒周期性旋轉(zhuǎn)時(shí)角速度最小時(shí)對(duì)應(yīng)的位置與雷諾數(shù)有關(guān)，隨著雷諾數(shù)增大這一位置越來越偏離水平方向.這在以往的研究中未曾涉及過.最近，Huang等[24]細(xì)致地研究了雷諾數(shù)及受限比等對(duì)橢圓顆粒旋轉(zhuǎn)特性的影響，并給出了橢圓顆粒靜止時(shí)的臨界雷諾數(shù)，同時(shí)他們還分析了顆粒初始位置不在區(qū)域中心位置時(shí)的運(yùn)動(dòng)趨勢(shì)，但對(duì)于橢圓顆粒靜止時(shí)的狀態(tài)沒有進(jìn)行研究.基于以上考慮，本文將采用Lallemand等[25]提出的基于反彈格式的格子Boltzmann方法研究橢圓顆粒在剪切流中旋轉(zhuǎn)運(yùn)動(dòng)的特性，主要關(guān)注較大范圍內(nèi)雷諾數(shù)(0＜Re≤170)及橢圓顆粒的長(zhǎng)短軸之比(1≤α≤10)對(duì)顆粒旋轉(zhuǎn)特性的影響，同時(shí)研究在較大雷諾數(shù)下流場(chǎng)結(jié)構(gòu)的特性.基于反彈格式的格子Boltzmann方法最早由 Ladd[26-27]提出來，但其將顆粒的邊界點(diǎn)置于網(wǎng)格點(diǎn)的中心，因此在描述曲線邊界時(shí)具有較大的誤差，而且這種方法僅適用于顆粒密度遠(yuǎn)大于流體密度的情況.Aidun等[28]改進(jìn)了這一缺陷，但仍然將顆粒的邊界點(diǎn)置于網(wǎng)格點(diǎn)的中心.Lallemand等[25]采用插值的思想改進(jìn)了這一方法，使得格子Boltzmann方法可以更準(zhǔn)確地模擬具有曲線邊界的顆粒運(yùn)動(dòng).

1 模型

1.1 數(shù)值模型

本文采用帶有單個(gè)松弛因子的格子 Boltzmann方程模型(LBGK)，該模型是目前應(yīng)用最為廣泛的一種格子Boltzmann模型[29-31]，LBGK不僅在編程上較為精簡(jiǎn)，又能夠保證足夠的精度，因此適用于解決流體流動(dòng)問題.LBGK模型的離散方程如下

式中，fi(x,t)表示分子在i方向上的速度分布函數(shù)，表示i方向上的平衡態(tài)速度分布函數(shù)，ci為分子在i方向的速度，τ為無量綱的松弛時(shí)間.流體的宏觀速度u和密度ρ可以通過下式求解

采用Qian等[21]提出的D2Q9格式，離散的格子速度為

假定馬赫數(shù)很小，通過對(duì)式(1)進(jìn)行Chapman-Enskog展開，可以導(dǎo)出不可壓N-S方程組

式中，σ為應(yīng)力張量，σ=-pI+2μS，μ為流體動(dòng)力黏度，p為壓力，S為應(yīng)變率張量,I為單位矩陣.

以上方法可以解決流體的流動(dòng)問題，而對(duì)于顆粒與流體相互作用的耦合問題,則采用基于反彈格式的動(dòng)量交換法[25]來處理運(yùn)動(dòng)的邊界，該方法可以看作是對(duì)反彈格式的一種改進(jìn)，反彈格式是格子Boltzmann方法中處理固定平直邊界條件的一種常用方法，即假設(shè)粒子與固定壁面碰撞后速度反轉(zhuǎn).而對(duì)于運(yùn)動(dòng)的邊界，不僅要考慮邊界的速度對(duì)流體的作用，還要處理好計(jì)算域中固體節(jié)點(diǎn)與流體節(jié)點(diǎn)的關(guān)系.如圖1(a)所示，虛線為此時(shí)的邊界，邊界右側(cè)陰影部分表示固體區(qū)域，實(shí)心點(diǎn)xs表示固體節(jié)點(diǎn)，邊界左側(cè)的空心點(diǎn)表示流體節(jié)點(diǎn)，流體節(jié)點(diǎn)xft是為了插值構(gòu)造的虛擬節(jié)點(diǎn)，粒子在邊界處發(fā)生碰撞，利用流體節(jié)點(diǎn)進(jìn)行二次插值得到流體節(jié)點(diǎn)經(jīng)壁面反彈之后的速度分布函數(shù)，具體的插值格式如下

式中，q表示流體節(jié)點(diǎn)到界面的距離與固體節(jié)點(diǎn)到界面的距離之比，i表示該分方向指向固體區(qū)域，表示與i的方向相反，uw表示壁面的速度.隨著固體邊界的移動(dòng)，一些固體節(jié)點(diǎn)在下一時(shí)刻會(huì)變?yōu)榱黧w節(jié)點(diǎn)，如圖1(b)所示，虛線為上一時(shí)刻的固體邊界的位置，實(shí)線為這一時(shí)刻的位置，圖中陰影部分為此時(shí)的固體區(qū)域，圖中流體節(jié)點(diǎn)xf在前一時(shí)刻為固體節(jié)點(diǎn)，此時(shí)需要重新計(jì)算指向固體區(qū)域的速度分布函數(shù)，在模擬時(shí)可以用非平衡外推格式計(jì)算，在界面上發(fā)生的動(dòng)量交換可以用下式計(jì)算

圖1 反彈邊界條件Fig.1 Illustration of the bounced-back scheme

圖1 反彈邊界條件(續(xù))Fig.1 Illustration of the bounced-back scheme(continued)

利用式(8)得到的動(dòng)量可以計(jì)算出流體對(duì)顆粒的合作用力和合力矩，再由牛頓第二定律確定顆粒的運(yùn)動(dòng).

1.2 物理模型

建立物理模型如圖2所示，橢圓形顆粒被放置在穩(wěn)定的二維簡(jiǎn)單剪切流中，上下剪切面的速度固定為U0，計(jì)算區(qū)域固定為L(zhǎng)×H=2400×480，邊界條件采用非平衡外推格式.初始時(shí)刻橢圓形顆粒位于區(qū)域的中心處(L/2,H/2)，橢圓顆粒的長(zhǎng)半軸長(zhǎng)a=48，b為橢圓的短半軸長(zhǎng)，長(zhǎng)短軸的比α=a/b，橢圓顆粒長(zhǎng)軸與x軸正方向的夾角為θ，旋轉(zhuǎn)的角速度為ω.假定橢圓顆粒的密度與流體密度相等，因此橢圓顆粒懸浮在流體中.流場(chǎng)的剪切率為G=2U0/H，定義雷諾數(shù)為Re=4Ga2/ν，ν為流體的運(yùn)動(dòng)黏度.

圖2 橢圓形顆粒在剪切流中的運(yùn)動(dòng)示意圖Fig.2 Schematic of an elliptical particle subjected to simple shear fl w

2 驗(yàn)證

為了驗(yàn)證本文采用的算法和計(jì)算代碼，首先計(jì)算了在極低雷諾數(shù)(Re=0.08)下的橢圓旋轉(zhuǎn)角速度，此時(shí)設(shè)定橢圓顆粒長(zhǎng)短軸之比α=2.將結(jié)果與Je ff ery[1]的理論解進(jìn)行了對(duì)比.Je ff ery得到的橢圓形顆粒在剪切流中的角速度與角度的關(guān)系如下

如圖3所示(Gt為無量綱化時(shí)間)，可以看到，本文的模擬結(jié)果與精確解符合得很好.需要指出的是，Je ff ery[1]的理論解是在Re=0的前提下得到的，本文由于采用數(shù)值方法求解，因此只能選擇盡可能小的雷諾數(shù)進(jìn)行對(duì)比.從圖3可知，當(dāng)選擇Re=0.08時(shí)已經(jīng)足夠接近理論解了.

圖3 Je ff ery[1]理論解與本文模擬結(jié)果的對(duì)比Fig.3 Comparison of Je ff ery solution[1]and the present simulation result

另外，為了進(jìn)一步說明本文算法的可靠性，還計(jì)算了不同雷諾數(shù)下橢圓顆粒旋轉(zhuǎn)的周期隨雷諾數(shù)的變化情況，此時(shí)仍固定長(zhǎng)短軸之比α=2.將結(jié)果與Ding等[22]的結(jié)果進(jìn)行對(duì)比.如圖4所示，可以看到二者符合得較好.圖中T為實(shí)際周期，GT則為無量綱化的周期.

圖4 不同雷諾數(shù)下橢圓顆粒的旋轉(zhuǎn)周期的對(duì)比Fig.4 The period of the rotation of the elliptical particle at di ff erent Reynolds number

3 計(jì)算結(jié)果及討論

3.1 雷諾數(shù)變化的影響

首先考察雷諾數(shù)的變化對(duì)橢圓形顆粒在剪切流中運(yùn)動(dòng)的影響.通道的長(zhǎng)度L固定為2400，高度H固定為480，設(shè)定橢圓顆粒的長(zhǎng)半軸a=48，長(zhǎng)短軸之比α=2.初始時(shí)刻橢圓顆粒被放置在計(jì)算區(qū)域的中心 (L/2,H/2)，顆粒的轉(zhuǎn)角θ0=π/2，雷諾數(shù)Re分別定為 5,10,15,30.模擬結(jié)果如圖5所示，橢圓形顆粒在簡(jiǎn)單剪切流的作用下沿著逆時(shí)針做旋轉(zhuǎn)運(yùn)動(dòng)，從圖5(a)～圖5(c)可以看到，當(dāng)雷諾數(shù)Re＜30時(shí)橢圓形顆粒的旋轉(zhuǎn)的角度和角度變化都呈現(xiàn)周期性變化的特點(diǎn).對(duì)于一次完整的旋轉(zhuǎn)周期，顆粒角速度ω曲線呈現(xiàn)“雙駝峰”形，即當(dāng)轉(zhuǎn)角θ∈(0,π/2)和(π,3π/2)時(shí)，ω隨著轉(zhuǎn)角θ的增大而增大，當(dāng)θ∈(π/2,π)和(3π/2,2π)時(shí)，ω隨著θ的增大而減小，橢圓形顆粒的長(zhǎng)軸轉(zhuǎn)到平行于x軸位置(θ=0或π)附近時(shí)，剪切流對(duì)于顆粒的有效力臂較小，從而產(chǎn)生的力矩也較小，顆粒旋轉(zhuǎn)的角速度達(dá)到最小值，而當(dāng)橢圓形顆粒的長(zhǎng)軸轉(zhuǎn)到垂直于x軸位置(θ=π/2或3π/2)附近時(shí)，橢圓形顆粒的長(zhǎng)軸垂直于剪切流剪切方向，因此產(chǎn)生的力矩較大，顆粒旋轉(zhuǎn)的角速度達(dá)到最大值.

圖5 不同雷諾數(shù)下橢圓顆粒旋轉(zhuǎn)的角度和角速度Fig.5 The orientation and rotation velocity of the elliptical particle at di ff erent Reynolds numbers

此外，從圖5(b)和圖5(c)還可以看到，隨著雷諾數(shù)的增大，速度圖中的豎虛線(橢圓顆粒旋轉(zhuǎn)角θ=0對(duì)應(yīng)的位置)開始漸漸向右側(cè)偏移，說明橢圓角速度的最小值不再出現(xiàn)在傾斜角為0(或π)的位置，而是隨著雷諾數(shù)的增大逐漸提前，也就是說，當(dāng)橢圓的長(zhǎng)軸還沒有到水平位置時(shí)其角速度已經(jīng)為最小了，這與Je ff ery[1]的結(jié)果不同.為了定量說明這一問題，本文給出了上述傾角θm與雷諾數(shù)Re的關(guān)系，如圖6所示.很明顯，根據(jù)雷諾數(shù)的大小,上述關(guān)系可以分為兩個(gè)區(qū)域,即0＜Re≤2.875與2.875＜Re≤25，而且在這兩個(gè)區(qū)域中θm與Re近似呈線性關(guān)系.與此同時(shí)，還給出了兩個(gè)區(qū)域的擬合曲線，即式(10)，可以看到，在兩個(gè)區(qū)域中θm都隨Re的增大而增大，但Re較小時(shí)θm增大的速率快得多，這與Re較小時(shí)流場(chǎng)以黏性為主導(dǎo)有關(guān).

圖6 橢圓顆粒轉(zhuǎn)速最小時(shí)對(duì)應(yīng)的轉(zhuǎn)角Fig.6 The orientation of the elliptical particle when its rotation velocity is smallest

另外，從圖5還可以看到，隨著雷諾數(shù)的增大，橢圓顆粒的旋轉(zhuǎn)周期也逐漸變大，當(dāng)雷諾數(shù)Re≥30時(shí)，橢圓顆粒最終靜止在水平附近的位置(Re＞30的數(shù)據(jù)未列出，但趨勢(shì)相同)，這與Huang等[24]的研究結(jié)果相符合.

為了更進(jìn)一步分析雷諾數(shù)對(duì)橢圓顆粒旋轉(zhuǎn)周期的影響，在固定其他參數(shù)不變的情況下，本文還模擬了Re=0.05,0.08,1,2,3,4,5,10,15,20,24,26,28時(shí)橢圓顆粒在剪切流中的運(yùn)動(dòng)情況，結(jié)果如圖7所示.

圖7 不同雷諾數(shù)下橢圓顆粒在剪切流中的旋轉(zhuǎn)周期Fig.7 The period of the rotation of elliptical particle at di ff erent Reynolds number

通過計(jì)算可以獲得相應(yīng)的周期，對(duì)周期與雷諾數(shù)的關(guān)系進(jìn)行最小二乘法擬合，采用Ding等[22]提出的如下擬合函數(shù)

式中，Rec表示臨界雷諾數(shù)，此處選擇Rec=29，此時(shí)模擬的結(jié)果與擬合曲線符合得較好.從圖中可以看到，當(dāng)Re＜3時(shí)，橢圓顆粒的轉(zhuǎn)動(dòng)周期GT幾乎不變化，但當(dāng)Re＞3時(shí)，轉(zhuǎn)動(dòng)周期隨Re的增大而開始明顯延長(zhǎng)，該分界點(diǎn)與圖6所示的分界點(diǎn)相似.這應(yīng)該與Re不斷增大而引起的流體慣性作用有關(guān).從圖7可知，當(dāng)Re繼續(xù)增大時(shí)，周期趨于無窮大，說明此時(shí)顆粒不再轉(zhuǎn)動(dòng)，而是靜止在水平位置附近.

當(dāng)雷諾數(shù)大于臨界值時(shí)橢圓顆粒在剪切流中會(huì)呈靜止?fàn)顟B(tài)，為了說明這一點(diǎn)，首先給出Re=0.1, 1和10時(shí)的流線結(jié)構(gòu)及壓力分布，如圖8所示.圖中的橢圓轉(zhuǎn)角均為0.94π.從圖8可以看到，由于橢圓的存在，流線結(jié)構(gòu)分成兩部分，一部分是剪切層，處于顆粒的上下位置，對(duì)顆粒產(chǎn)生的力矩是正的，另一部分是回流層，處于顆粒的左右位置，產(chǎn)生的力矩是負(fù)的.很顯然，當(dāng)雷諾數(shù)很小時(shí),如Re=0.1和1時(shí)，回流層的區(qū)域很小且回流層的流體沒有接觸到橢圓顆粒，因此對(duì)顆粒產(chǎn)生的負(fù)力矩很小，此時(shí)以剪切層的流體作用為主導(dǎo)，橢圓顆粒會(huì)沿剪切方向進(jìn)行周期性地旋轉(zhuǎn).當(dāng)雷諾數(shù)增大時(shí)，如Re=10時(shí)，此時(shí)可以看到回流層的范圍明顯增大，且回流層的流體與顆粒直接接觸，此時(shí)對(duì)顆粒產(chǎn)生的力矩會(huì)增大，另一方面，從圖中還可以看到，橢圓顆粒周圍的壓力分布也隨雷諾數(shù)的增大而發(fā)生顯著改變，無論在顆粒的左端還是右端，其上下的壓差隨雷諾數(shù)的增大而增大，并且此壓差產(chǎn)生負(fù)力矩.綜合以上兩個(gè)因素可知，當(dāng)雷諾數(shù)增大到一定程度時(shí)，在剪切層流體、回流層流體以及壓力分布的共同作用下，橢圓所受力矩可能為零，此時(shí)橢圓會(huì)處于靜止?fàn)顟B(tài).

圖8 橢圓轉(zhuǎn)角θ=0.94π時(shí)不同雷諾數(shù)下的流線結(jié)構(gòu)及壓力分布Fig.8 The streamlines and pressure forθ=0.94π at di ff erent Reynolds numbers

為了觀察橢圓顆粒靜止時(shí)的狀態(tài)，給出了Re=30,90和150時(shí)的流線結(jié)構(gòu)以及壓力分布，如圖9所示.可以看到流動(dòng)為穩(wěn)定狀態(tài)，且橢圓轉(zhuǎn)角θ均不等于π，雷諾數(shù)越大，θ角越小，也就是說顆粒越來越遠(yuǎn)離水平位置.更進(jìn)一步地，圖10和圖11分別給出了雷諾數(shù)為30，40和150時(shí)橢圓傾角和角速度隨時(shí)間的變化曲線，由圖可知，當(dāng)雷諾數(shù)大于臨界雷諾數(shù)時(shí)，橢圓顆粒最終靜止在水平位置附近，但傾角小于π，且隨著雷諾數(shù)的增大，最終傾角逐漸減小.另外，對(duì)于較大的雷諾數(shù)，如Re=150時(shí),會(huì)造成橢圓顆粒旋轉(zhuǎn)的過沖現(xiàn)象，即橢圓的傾角先增大，達(dá)到最大值后逐漸減小，最后趨于穩(wěn)定.這在Ding等[22]的結(jié)果中并未出現(xiàn).從對(duì)應(yīng)的角速度變化圖中可以看到，當(dāng)顆粒發(fā)生過沖現(xiàn)象后，它的角速度為負(fù)值，隨后逐漸增大并最終趨于零.

圖9 橢圓顆粒靜止時(shí)不同雷諾數(shù)下的流線結(jié)構(gòu)和壓力分布Fig.9 The streamlines and pressure at di ff erent Reynolds numbers when the particle becomes stationary

圖10 不同雷諾數(shù)下橢圓顆粒傾角隨時(shí)間的變化Fig.10 Time history of orientation of elliptical particle for di ff erent Reynolds numbers

圖11 不同雷諾數(shù)下橢圓顆粒角速度隨時(shí)間的變化Fig.11 Time history of rotational velocity of elliptical for di ff erent Reynolds numbers

為了定量分析雷諾數(shù)對(duì)橢圓顆粒最終傾角的影響，考慮雷諾數(shù)30≤Re≤170，計(jì)算出橢圓顆粒達(dá)到穩(wěn)定狀態(tài)時(shí)的轉(zhuǎn)角.采用如下的擬合函數(shù)并通過最小二乘法擬合數(shù)據(jù)

結(jié)果如圖12所示，擬合曲線與計(jì)算的結(jié)果符合得很好，說明顆粒最終達(dá)到穩(wěn)定的靜止?fàn)顟B(tài)時(shí)，傾角與雷諾數(shù)也呈冪函數(shù)關(guān)系.另外可以看到，當(dāng)雷諾數(shù)大于臨界雷諾數(shù)時(shí)，隨著雷諾數(shù)的增大，最終的傾角逐漸減小，說明顆粒的長(zhǎng)軸逐漸遠(yuǎn)離水平位置.

3.2 橢圓顆粒長(zhǎng)/短軸之比的影響

圖12 不同雷諾數(shù)下橢圓顆粒的最終傾角Fig.12 Final orientation of elliptical particle at di ff erent Reynolds numbers

下面考察顆粒的長(zhǎng)短軸之比α對(duì)橢圓顆粒在剪切流中運(yùn)動(dòng)的影響.同時(shí)考慮兩組雷諾數(shù)的情況，即Re=5和10，模型中其他參數(shù)不變，固定顆粒的長(zhǎng)軸a=48,通過改變橢圓顆粒的短軸b來改變?chǔ)?由于兩組的結(jié)果比較相似，下文以雷諾數(shù)Re=5為例，給出α=2,3,4,5.5時(shí)的模擬結(jié)果，如圖13所示.從圖中可以看到結(jié)果與圖5的情況相似.當(dāng)長(zhǎng)/短軸之比α=5.5時(shí)，隨著α的增大，顆粒旋轉(zhuǎn)的周期越來越大.顆粒旋轉(zhuǎn)角速度到達(dá)小值時(shí)的角度也會(huì)提前，從圖13(c)中可以明顯地觀察到這一現(xiàn)象.而當(dāng)α≥5.5時(shí)，圖13(d)結(jié)果表明橢圓顆粒最終會(huì)停止在傾角θ=π附近.

圖13 不同α對(duì)應(yīng)的橢圓顆粒旋轉(zhuǎn)角度和角速度的變化Fig.13 The orientation and rotational velocity of elliptical particle for di ff erent α

圖13 不同α對(duì)應(yīng)的橢圓顆粒旋轉(zhuǎn)角度和角速度的變化(續(xù))Fig.13 The orientation and rotational velocity of elliptical particle for di ff erent α(continued)

為了定量研究長(zhǎng)短軸之比對(duì)橢圓顆粒旋轉(zhuǎn)周期的影響，對(duì)于兩組雷諾數(shù)Re=5和10，分別設(shè)置1≤α≤5和1≤α≤3.通過計(jì)算得到橢圓顆粒的旋轉(zhuǎn)周期，結(jié)果如圖14所示，從圖中可以看到，橢圓顆粒的旋轉(zhuǎn)周期隨著橢圓顆粒長(zhǎng)短軸之比的增大而延長(zhǎng).說明在相同雷諾數(shù)的情況下，越細(xì)長(zhǎng)的橢圓轉(zhuǎn)動(dòng)得越慢.當(dāng)α=1時(shí)，橢圓退化成圓，此時(shí)顆粒旋轉(zhuǎn)一周的時(shí)間最短.采用最小二乘法對(duì)離散數(shù)據(jù)進(jìn)行擬合，得到擬合函數(shù)如下

從圖14中可以看出，離散數(shù)據(jù)與擬合函數(shù)符合的很好，說明橢圓顆粒旋轉(zhuǎn)的周期與橢圓顆粒長(zhǎng)短軸之比α呈冪函數(shù)關(guān)系.對(duì)于Re=5，長(zhǎng)短軸之比的臨界值約為5.5，而對(duì)于Re=10，臨界值約為3.5.說明雷諾數(shù)越大，臨界的α值越小.

圖14 雷諾數(shù)分別為5和10對(duì)應(yīng)的周期與α的關(guān)系Fig.14 The period of the rotation of elliptical particle for di ff erent α atRe=5 and 10,respectively

當(dāng)橢圓顆粒的長(zhǎng)短軸之比大于臨界值時(shí)，橢圓顆粒不再以一定的周期旋轉(zhuǎn)，而是靜止在θ=π附近，我們給出雷諾數(shù)Re=10時(shí)α分別為4,6和8對(duì)應(yīng)的流線結(jié)構(gòu)，如圖15所示，可以看到此時(shí)流場(chǎng)處于穩(wěn)定的狀態(tài).然而，此時(shí)并沒有發(fā)現(xiàn)橢圓顆粒旋轉(zhuǎn)的過沖現(xiàn)象.

圖15 Re=10時(shí)不同長(zhǎng)短軸之比對(duì)應(yīng)的流線結(jié)構(gòu)及壓力分布Fig.15 The streamlines and pressure for di ff erent α atRe=10

下面研究最終傾角與α的關(guān)系，針對(duì)Re=5和Re=10，分別設(shè)定6≤α≤10和3.5≤α≤10,計(jì)算出顆粒靜止時(shí)刻的轉(zhuǎn)角，結(jié)果如圖16所示，可以看到橢圓顆粒的轉(zhuǎn)角隨著長(zhǎng)短軸之比的增大而減小.也就是說，越細(xì)長(zhǎng)的橢圓，靜止時(shí)越遠(yuǎn)離水平位置.采用最小二乘法擬合，得到如下擬合函數(shù)

從圖中可以看到，離散數(shù)據(jù)與擬合曲線符合得較好，與雷諾數(shù)的影響類似，橢圓顆粒靜止時(shí)的轉(zhuǎn)角與橢圓的長(zhǎng)短軸之比呈冪函數(shù)關(guān)系.

圖16 雷諾數(shù)分別為5和10對(duì)應(yīng)的最終傾角與α的關(guān)系Fig.16 The fina orientation of elliptical particle for di ff erent α atRe=5 and 10,respectively

4 結(jié)論

本文采用格子Boltzmann方法對(duì)橢圓顆粒在剪切流中的運(yùn)動(dòng)進(jìn)行來直接數(shù)值模擬.主要研究了雷諾數(shù)和橢圓顆粒的長(zhǎng)短軸之比對(duì)橢圓顆粒旋轉(zhuǎn)特性的影響，有以下結(jié)論：

(1)當(dāng)雷諾數(shù)小于臨界值時(shí)，橢圓顆粒的旋轉(zhuǎn)周期與雷諾數(shù)呈冪函數(shù)關(guān)系，雷諾數(shù)越大，旋轉(zhuǎn)周期越大；顆粒角速度最小時(shí)對(duì)應(yīng)的長(zhǎng)軸方向隨著雷諾數(shù)的增大而逐漸遠(yuǎn)離水平方向，且其傾角與雷諾數(shù)呈分段線性的關(guān)系.

(2)當(dāng)雷諾數(shù)超過臨界值時(shí)，橢圓顆粒最終靜止在剪切流場(chǎng)中，且傾角與雷諾數(shù)呈冪函數(shù)關(guān)系，雷諾數(shù)越大傾角越小.

(3)當(dāng)橢圓顆粒長(zhǎng)短軸之比α小于臨界值時(shí)，顆粒旋轉(zhuǎn)的周期與α呈冪函數(shù)關(guān)系，且隨著α的增大而延長(zhǎng).當(dāng)α超過臨界值時(shí)，顆粒最終以一定的傾角保持靜止?fàn)顟B(tài)，且傾角與α也呈冪函數(shù)關(guān)系，α越大，顆粒最終的傾角越小.另外，雷諾數(shù)越大，臨界的α值越小.

1 Je ff ery GB.The motion of ellipsoidal particles immersed in a viscous flui//Proceedings of the Royal Society A,1922,102(715): 161-179

2 Batchelor GK,Green JT.The hydrodynamicinteraction of two small freely-moving spheres in a linear fl w fieldJournal of Fluid Mechanics,1972,56(2):375-400

3 Feng J,Hu H,Joseph D.Direct simulation of initial value problems for the motion of solid bodies in a newtonian flui Part 1.Sedimentation.Journal of Fluid Mechanics,1994,261:95-134

4 Feng J,Joseph DD.The unsteady motion of solid bodies in creeping fl ws.Journal of Fluid Mechanics,1995,303:83-102

5 Ladd AJC.Sedimentation of homogeneous suspensions of non-Brownian spheres.Physics of Fluids,1997,9:491-499

6 Aidun CK,Ding EJ.Dynamics of particle sedimentation in a vertical channel:period doubling bifurcation and chaotic state.Physics of Fluids,2003,15(6):1612

7 費(fèi)明龍,徐小蓉,孫其誠(chéng)等.顆粒介質(zhì)固--流態(tài)轉(zhuǎn)變的理論分析及實(shí)驗(yàn)研究.力學(xué)學(xué)報(bào),2016,48(1):48-55(Fei Minglong,Xu Xiaorong,Sun Qicheng,et al.Studies on the transition between solidand fluid-li e states of granular materials.Chinese Journal of Theoretical and Applied Mechanics,2016,48(1):48-55(in Chinese))

8 Ardekani AM,Rangel RH.Unsteady motion of two solid spheres in Stokes fl w.Physics of Fluids,2006,18:103306

10 胡平,張興偉,牛小東等.三圓形顆粒在通道中沉降運(yùn)動(dòng)的數(shù)值研究.力學(xué)學(xué)報(bào),2014,46(5):673-684(Hu Ping,Zhang Xingwei,Niu Xiaodong,et al.Numerical study on the sedimented motion characteristics of three aligned circular particles in the inclined channels.ChineseJournalofTheoreticalandAppliedMechanics,2014,46(5): 673-684(in Chinese))

12 Hwang WR,Hulsen MA,Meijer HEH.Direct simulations of particle suspensions in a viscoelastic flui in sliding bi-periodic frames.Journal of Non-Newtonian Fluid Mechanics,2004,121(1):15-33

13 Choi YJ,Hulsen MA,Meijer HEH.An extended finitelement method for the simulation of particulate viscoelastic fl ws.Journal of Non-Newtonian Fluid Mechanics,2010,165(11):607-624

14 Lundell F,Carlsson A.Heavy ellipsoids in creeping shear fl w: Transitions of the particle rotation rate and orbit shape.Physical Review E Statistical Nonlinear&Soft Matter Physics,2010,81(1): 016323

15 Pasquino R,D’Avino G,Ma ff ettone PL,et al.Migration and chaining of noncolloidal spheres suspended in a sheared viscoelastic medium.Experiments and numerical simulations.Journal of Non-Newtonian Fluid Mechanics,2014,203(1):1-8

16 Mikulencak DR,Morris JF.Stationary shear fl w around fi ed and free bodies at finit Reynolds number.Journal of Fluid Mechanics, 2004,520:215-242

17 Subramanian G,Koch DL.Inertial e ff ects on the transfer of heat or mass from neutrally buoyant spheres in a steady linear velocity fieldPhysics of Fluids,2006,18:073302

18 Subramanian G,Koch DL.Centrifugal forces alter streamline topology and greatly enhance the rate of heat and mass transfer from neutrally buoyant particles to a shear fl w.Physical Review Letters, 2006,96:134503

19 Yacoubi AE,Xu S,Wang ZJ.Computational study of the interaction of freely moving particles at intermediate Reynolds numbers.Journal of Fluid Mechanics,2012,705(2):134-148

20 Nie D,Lin J,Zheng M.Direct numerical simulation of multiple particles sedimentation at an intermediate reynolds number.Communications in Computational Physics,2014,16(3):675-698

21 Nie D,Lin J,Chen R.Grouping behavior of coaxial settling particles in a narrow channel.Physical Review E Statistical Nonlinear&Soft Matter Physics,2016,93:013114

22 Ding E,Aidun CK.The dynamics and scaling law for particles suspended in shear fl w with inertia.Journal of Fluid Mechanics,2000, 423:317-344

23 Huang H,Yang X,Krafczyk M,et al.Rotation of spheroidal particles in Couette fl ws.Journal of Fluid Mechanics,2012,692:369-394

24 Huang SL,Chen SD,Pan TW,et al.The motion of a neutrally buoyant particle of an elliptic shape in two dimensional shear fl w:a numerical study.Physics of Fluids,2015,27(5):083303

25 Lallemand P,Luo LS.Lattice Boltzmann method for moving boundaries.Journal of Computational Physics,2003,184(2):406-421

26 Ladd AJC.Numerical simulations of particulate suspensions via a discretized Boltzmann equation.Part I.Theoretical foundation.Journal of Fluid Mechanics,1994,271:285-309

27 Ladd AJC.Numerical simulations of particulate suspensions via a discretized Boltzmann equation.Part II.Numerical results.Journal of Fluid Mechanics,1994,271:311-339

28 Aidun CK,Lu Y,Ding E.Direct analysis of particulate suspensions with inertia using the discrete Boltzmann equation.Journal of Fluid Mechanics,1998,373:287-311

29 Benzi R,Succi S,Vergassola MR.The lattice Boltzmann equation: theory and applications.Physics Reports,1992,222:145-197

30 Qian YH,D’Humires D,Lallemand P.Lattice BGK models for Navier-Stokes equation.Europhysics Letters,1992,17(6):479-484

31 Chen SY,Doolen GD.Lattice Boltzmann method for fluifl ws.Annual Review of Fluid Mechanics,1998,30:329-364

NUMERICAL STUDY ON THE ROTATION OF AN ELLIPTICAL PARTICLE IN SHEAR FLOW1)

Chen Rongqian Nie Deming2)
(College of Metrology and Measurement Engineering，China Jiliang University，Hangzhou310018，China)

A thorough understanding of the behavior of particles freely suspended in a shear fl w is fundamentally important for understanding and predicting fl w behavior of particle suspensions.The motion of particles is very complex when the flui inertia is taken into account.In the present study,the lattice Boltzmann method has been used to simulate the rotation of an elliptical particle in simple shear fl w at intermediate Reynolds numbers.Firstly,the e ff ect of the Reynolds number(0＜Re≤170)has been studied.Results show that the particle rotates periodically whenReis smaller than a critical value.The orientation of the particle at which the particle has its minimum angular velocity decreases asReincreases,which has a piecewise linear relationship withRe.Moreover,the rotation period has a power-law relationship withRe.The largerReis,the larger the rotation period is.However,whenReis greater than the critical value,the elliptical particle will reach a steady state.Results show that the fina orientation of the elliptical particle has a power-law relationship withRefor the steady state.The largerReis,the smaller the orientation is.Secondly,the e ff ect of the ratio of major axis/minor axis α(1≤α≤10)has also been studied.It shows that there is also a power-law relationship between the rotation period and α.The larger the value of α is,the smaller the rotation period is.Similarly,when αis greater than a critical value,the elliptical particle does not rotate.The fina orientation of the elliptical particle has a power-law relationship with α.The larger the value of α is,the smaller the orientation is.Furthermore,it also shows that the overshoot is observed when the elliptical particle is rotating ifReis larger enough.

lattice Boltzmann method,shear fl w,elliptical particle,rotation

O359

A

10.6052/0459-1879-16-002

2016–01–04收稿，2016–12–30錄用,2017–01–05網(wǎng)絡(luò)版發(fā)表.

1)國(guó)家自然科學(xué)基金(11272302,11132008)和浙江省自然科學(xué)基金(LY15A020004)資助項(xiàng)目.

2)聶德明,教授,主要研究方向：顆粒多相流、格子Boltzmann方法.E-mail:nieinhz@cjlu.edu.cn

陳榮前,聶德明.橢圓顆粒在剪切流中旋轉(zhuǎn)特性的數(shù)值研究.力學(xué)學(xué)報(bào),2017,49(2):257-267

Chen Rongqian,Nie Deming.Numerical study on the rotation of an elliptical particle in shear fl w.Chinese Journal of Theoretical and Applied Mechanics,2017,49(2):257-267