摘要：在工農業生產中，存在大量非牛頓流體在管道中的流動現象，比如石油鑽井採集液的集輸、高分子聚合物塑膠製品的生產加工等。相比傳統的牛頓流體，關於非牛頓流體在管道內的流動的研究還有很大發展空間，因此對其機理的研究具有重要的現實意義。藉助多物理場模擬軟體COMSOLMultiphysics，對非牛頓流體在管道內的流動現象進行了數值模擬研究。主要研究了非牛頓流體在不同截面直管中的流動，分析了流體在管道內的速度場分佈情況。結果表明，在同一圓管道內，冪律指數的增大會使冪律流體表現出剪下增稠效應，降低流動性，進而縮短流體的速度入口段長度；對於相同的邊界條件，不同的管道橫截面形狀會影響冪律非牛頓流體的速度分佈；管內流動的Carreau非牛頓流體的剪下應變速率與黏度關係符合其本構方程特性，相比牛頓流體表現出了明顯較差的流動性。
關鍵詞： 一維剪下 冪律流體 管道 非牛頓流體
1687年牛頓提出，作一維剪下流動的水，其剪下應變速率與剪下應力的大小成正比，這個規律就是後來著名的牛頓內摩擦定律。在流變學中，流變性符合這一規律的流體被稱作牛頓流體；反之，則為非牛頓流體。相比於牛頓流體，非牛頓流體在工農業生產乃至醫學研究中出現得更為廣泛，比如石油鑽井採出液的集輸處理、聚合物塑膠製品加工、人體血液在血管中的流動等。上述這些情形都涉及非牛頓流體在管道內的流動問題，因此非牛頓流體在各種管道環境下的流動機理具有充分的研究價值。1867年，J.C.麥克斯韋提出線性黏彈性方程，開始了非牛頓流體力學的研究。但由於黏彈性流體問題具有的複雜性，一直到20世紀50年代後相關研究領域才得到較為迅速的發展，並逐漸成為一門獨立的學科[1]。相比牛頓流體在管道流動領域已較為成熟的研究成果，關於非牛頓流體的研究還有很大的發展空間。針對非牛頓流體的研究方法主要包括實驗法、解析解法與數值解法。實驗法最為直接，可檢驗其他方法的正確性，但是成本較高且實驗結果的普遍性不佳。解析解法是理論上最為準確理想的研究方法，透過建立合適的微分方程組，使用純數學方法得出方程的精確解；但對於非牛頓流體複雜的流動情況來說，求解的難度過大。數值解法則是應用計算機將物理場離散化，之後將流體微分方程組轉化為代數方程並求出各個節點上的引數值，屬於一種近似解法[2]。由於數值解法容易獲得且能保證足夠的求解精度，其已經成為研究非牛頓流體問題最為常用的方法。
國外學者對非牛頓流體的研究開始最早。在圓管道方面，CHEBBI[3]使用積分邊界層方法，分析了冪律流體層流在圓管入口段與充分發展區的流動特性，獲得了與實驗資料吻合度較高的結果。對於非圓管道，GUCKES[4]使用有限差分法求解冪律流體和賓漢流體，給出了2種流體在偏心圓環管中流動的平均體積流率。MITSUISHI等人[5]使用高聚合物非牛頓流體作為實驗材料，給出了偏心圓環管中非牛頓流體流速與壓降之間的關係。
國內對於非牛頓流體管內流動的研究起步較晚，成果也較國外少。姜篤志[6]分析了非牛頓流體在管道內流變性的一般規律並建立模型。楊旭等人[7]分析了非牛頓流體的本構及流動規律，基於空間分數階微積分方法，建立了分數階非牛頓流體本構模型，使用分數階的階數大小來反映非牛頓流體流動的空間記憶性強弱。
經過以上分析可以看出，中國對非牛頓流體管內流動這一研究領域的研究不太成熟，大多數停留在對非牛頓流體管內流動進行基本的數值模擬這一階段，而缺少對於模擬條件的多樣化拓展以及對結果資料之間關係更加詳盡的討論。關於非牛頓奈米流體管內流動這一新興領域更是空白。又考慮到非牛頓流體管內流動在工業上的廣泛應用，進一步進行相關方面的數值研究是很有必要的。
本文使用COMSOLMultiphysics這一軟體作為計算平臺，採用有限元分析法對非牛頓流體流經不同管道時體現出的流動特性進行數值模擬研究。
1、數學模型與數值模擬方法
1.1控制方程組
一切的流體流動過程，都以以下3個基本的物理學原理為基礎：質量守恆定律，牛頓第二定律與能量守恆定律。將這些物理學用於構建流動模型，將會匯出一組方程，即連續性方程、動量方程與能量方程。這些方程是上述物理學原理的數學描述，本文不討論傳熱，所以不引入能量守恆方程，質量守恆定律（連續性方程）：
（公式）
式（1）表示的是瞬態三維可壓流體流動的連續性方程，本文所分析的流體流動處於穩態且不可壓縮，密度ρ不會隨著時間的變化而改變，所以，流體流動的數學描述為：
（公式）
動量守恆方程：
（公式）
式（3)(4)(5）是對於任何流體都成立的動量守恆方程，是微元體內流體動量對於事件的變化率等於外界作用於該微元體上的各種力的和。簡稱動量方程，也稱納維斯托克斯（N-S）方程。
1.2本構方程
本構方程是反映物料宏觀性質的數學模型，又被稱為流變狀態方程或是流變方程。在流變學中，本構方程是在某些假定條件下，對流體或彈性體的材料力學行為的數學描述，可以用來區分流體型別。前述的牛頓內摩擦定律即最簡單的流體本構方程。本構方程與連續性方程、運動方程一起構成封閉的方程組，用於求解流體的流動特性。
在本文模擬中使用的冪律流體的本構方程式為：
（公式）
式（6）中：m為黏稠係數，Pa·sn，表示物料的黏稠程度；n為冪律流變指數（簡稱冪律指數），為無量綱量，表示非牛頓流體的流動特性偏離牛頓流體的程度（n=1時為牛頓流體）。
1.3非牛頓流體管道流動模型
1.3.1流體切應力
對於圓形管道內的流體流動，非管壁處流體剪應力ô與管壁處剪應力ôw的取值滿足下列均勻流動方程式：
（公式）
式（7)(8）中：∆p為管道壓降；r為圓柱座標系中的r方向座標位置；L為管道長度；R為管道半徑。
上述均勻流動方程式的推導並沒有涉及流體的性質與流動狀態，所以該方程式適用於所有的流體與流動狀態。
1.3.2剪下應變速率
剪下應變速率描述的是流體的剪下流動，定義為單位時間的剪下應變變化：
（公式）
式（9）中：γ為剪下應變。
值得注意的是，剪下應變速率常與速度梯度混淆。實際上二者是不同的概念。速度梯度是流體的速度對空間座標的導數，用du/dy來表示。在數學上，二者的數值有時相等，這是因為一般速度梯度符合：
（公式）
但二者的物理意義並不相同，且數值上有時並不相等（如流體在同軸圓筒之間的流動，此時有角速度的影響）。
1.3.3廣義雷諾數
對於圓管內非牛頓流體流動的雷諾數Re計算，目前多是仿照牛頓流體，近似按照黏度或者對比牛頓流體壓降的公式計算其廣義雷諾數Re´：
（公式）
式（11）中：D為圓管直徑；n為流變指數，不同流體二者的取值不同；k為流變係數。
對於牛頓流體來說，n´=1,k´=μ。所以，牛頓流體雷諾數為：
（公式）
在本文所設定的模擬環境中，非牛頓流體流動的廣義雷諾數均足夠小，整個管道流場處於層流狀態。
1.4求解器介紹
本次模擬使用的COMSOLMultiphysics是一款功能強大的多物理場模擬軟體，可以對多個領域的物理過程進行模擬計算。本文研究的是非牛頓流體在管道內的流動特性，因此選用的是COMSOLMultiphysics中的Laminarflow求解器模組。這個模組內建了許多流動的基本物理引數，可以在這些引數的基礎上結合前文所述的數學模型定義更多模擬所需要的變數。
1.5數值計算方法
本文涉及的數值計算使用的是有限元分析法（Finiteelementmethod,FEM）。它的基本思路為：一個物體或系統被分解為由多個相互聯結的、簡單、獨立的點組成的幾何模型。在這種方法中這些獨立的點的數量是有限的，因此被稱為有限元。由實際的物理模型中推匯出來的平衡方程式被使用到每個點上，由此產生了一個方程組，這個方程組可以用線性代數的方法求解。
1.6網格區域性加密與獨立性驗證
在數值模擬過程中，網格的數量與大小與模擬結果的精度密切相關。更多的網格數目雖然能夠提高結果的精度，但是也意味著要耗費更多的計算機算力資源。經衡量，使用COMSOLMultiphysics中的Meshrefine功能對模型的關鍵計算區域進行區域性網格加密，模型的其他部分則保持相對較大的網格密度。
區域性加密效果如圖1所示。
圖1網格區域性加密效果圖
此外，為了保證求解結果的精度，每組模擬計算時都需要多次調整網格的數量，檢查計算結果是否因為網格密度不同而出現較大的誤差。經過調整，確認每組模擬的網格數量保持在200萬～400萬之間為宜，這時求解精度已經足夠。
2、模擬結果及分析
2.1物理模型
2.1.1幾何模型
本組模擬選用了2種不同截面的長直管道，分別為圓形截面、矩形截面管道。這2種截面管道的形狀示意與幾何尺寸資料如圖2和表1所示，管道整體如圖3所示。
圖2管道截面形狀與尺寸示意圖
表1管道幾何尺寸資料
圖3管道整體示意圖
非牛頓流體從管道的一端流入，在管道內進行流動過程後，從另一端流出。
2.1.2流體性質
本組模擬選用的非牛頓流體為冪律流體（冪律指數n≠1），其黏度與剪下應變速率的關係遵循下式：
（公式）
其餘的流動特性引數列如表2所示。
表2不同截面管模擬中非牛頓流體的流動特性引數值
2.2邊界條件
選用如下邊界條件：(1)進口邊界條件。進口速度uin，取值固定為0.15m/s；進口溫度Tin為293.15K(20℃）。(2)出口邊界條件。出口壓力為0。(3)壁面邊界條件。無滑移邊界條件，即壁面處的流體速度為0；常壁溫邊界，壁面溫度TW為333.15K(60℃）。(4)其他。流體不可壓縮，所有管道區域內流體處於層流狀態。
2.3冪律指數對非牛頓流體管道流動的影響
對於冪律流體，當冪律指數n<1時，流體表現出剪下稀化效應；當冪律指數n>1時，流體表現出剪下增稠效應。為了探究冪律指數n對冪律流體在管道中流動特性的影響，控制進口速度uin=0.15m/s，管截面當量直徑de=0.1m，管長l=0.8m，引數相同，冪律流體n取0.6～1.4，對圓形截面管道內的非牛頓冪律流體進行流動特性的模擬。不同冪律指數下非牛頓冪律流體的速度場分佈如圖4所示。
圖4展示了不同冪律指數條件下冪律流體在管道各個橫截面上的速度場分佈。觀察靠近管道入口段的速度場橫截面，可以發現隨著n的增大，流體的速度入口段長度越來越短。這是因為當冪律指數n>1時，流體表現出剪下增稠效應，而且該效應會隨著n的增大而得到強化（剪下稀化效應同理）。剪下增稠效應下的流體流動性降低，速度梯度小，更容易達到充分發展狀態；反之，剪下稀化效應下的流體流動性提高，速度梯度大，速度入口段長。
此外，冪律指數也對流體充分發展區最大速度有一定影響，兩者呈正相關關係。上述結論均可以從MUKHERJEE的研究中得到驗證。
圖4不同冪律指數下非牛頓冪律流體的速度場分佈（單位：m/s)
2.4截面形狀對非牛頓流體管道流動的影響
保持進口速度uin=0.15m/s，冪律指數n=0.6，管截面當量直徑de=0.1m，管長l=0.8m，引數相同，對不同截面管道內的非牛頓流體進行流動特性模擬。
2.4.1截面形狀對速度入口段長度的影響
在本組模擬中，在各管道的軸向中心截面上繪製速度雲圖，觀察非牛頓流體在不同截面管道中的速度入口段長度。不同截面形狀管道軸向中心截面速度雲圖如圖5所示，從圖5可以看出，流體在管道入口段形成速度邊界層，邊界層逐漸向截面中心流動彙集，最後消失形成穩定的流態。矩形形狀的管道內非牛頓流體的速度入口段長度比圓形管道的長。這說明速度入口段的長度與管道截面形狀有關。截面為圓形的管道內的非牛頓流體更早到達充分發展階段，流態更為穩定。
圖5不同截面形狀管道軸向中心截面速度雲圖（單位：m/s)
2.4.2截面形狀對速度場分佈的影響
為了更直觀地觀察出管道內的速度場分佈情況，重新選取數個沿管道軸向等距分佈的管道橫截面來繪製速度雲圖。透過觀察每個截面雲圖的變化過程，可以得出非牛頓流體的速度沿徑向與軸向發展的情況。兩種截面形狀管道的橫截面速度場分佈如圖6所示，從圖6可以看出，流體在管壁處速度為0，最大速度中心位於各管道的橫截面幾何中心區域，但並不是一開始就在中心區域達到最大速度。在2種管道靠近入口段的截面，可看到速度最大值區域分散在截面的幾何中心周圍。對於矩形截面管道來說這種現象更為明顯，雲圖中能夠看到入口段截面形成兩三個速度中心區域，由此可見，速度最大區域首先出現在靠近管道截面角區的位置，之後在向軸向發展的同時向徑向截面幾何中心合併，最終在中間達到流動速度峰值。
3、結論
本文透過COMSOLMultiphysics平臺，使用有限元法與區域性最佳化網格技術建立了非牛頓流體在不同橫截面管道中的流動模型。
圖6兩種截面形狀管道的橫截面速度場分佈（單位：m/s)
透過以上數值模型，分析了非牛頓流體在管道內的速度分佈，並討論了流體流動特性引數之間的關係，以及管道截面形狀和冪律指數對流體流動效能的影響。具體總結如下：對於冪律流體，更高的冪律指數n值將使流體更多表現出剪下增稠效應，進一步使得速度入口段長度逐步減短；反之，則表現出剪下稀化效應，增長速度入口段的長度。在相同的邊界條件下，橫截面為矩形的管道中流體的速度入口段長度最長，隨後為圓形截面管道；管道內流體流動的發展從接近截面角區的位置開始，之後再向管道截面幾何中心處發展。
