陳榮前 聶德明

(中國計量大學計量測試工程學院，杭州310018)

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

陳榮前 聶德明2)

(中國計量大學計量測試工程學院，杭州310018)

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

格子Boltzmann方法，剪切流，橢圓顆粒，旋轉

引言

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

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

1 模型

1.1 數(shù)值模型

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

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

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

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

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

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

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

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

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

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

1.2 物理模型

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

圖2 橢圓形顆粒在剪切流中的運動示意圖Fig.2 Schematic of an elliptical particle subjected to simple shear fl w

2 驗證

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

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

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

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

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

3 計算結果及討論

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

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

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

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

圖6 橢圓顆粒轉速最小時對應的轉角Fig.6 The orientation of the elliptical particle when its rotation velocity is smallest

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 橢圓顆粒長/短軸之比的影響

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

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

圖13 不同α對應的橢圓顆粒旋轉角度和角速度的變化Fig.13 The orientation and rotational velocity of elliptical particle for di ff erent α

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

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

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

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

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

圖15 Re=10時不同長短軸之比對應的流線結構及壓力分布Fig.15 The streamlines and pressure for di ff erent α atRe=10

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

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

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

4 結論

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

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

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

(3)當橢圓顆粒長短軸之比α小于臨界值時，顆粒旋轉的周期與α呈冪函數(shù)關系，且隨著α的增大而延長.當α超過臨界值時，顆粒最終以一定的傾角保持靜止狀態(tài)，且傾角與α也呈冪函數(shù)關系，α越大，顆粒最終的傾角越小.另外，雷諾數(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 費明龍,徐小蓉,孫其誠等.顆粒介質固--流態(tài)轉變的理論分析及實驗研究.力學學報,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 胡平,張興偉,牛小東等.三圓形顆粒在通道中沉降運動的數(shù)值研究.力學學報,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網絡版發(fā)表.

1)國家自然科學基金(11272302,11132008)和浙江省自然科學基金(LY15A020004)資助項目.

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

陳榮前,聶德明.橢圓顆粒在剪切流中旋轉特性的數(shù)值研究.力學學報,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