李亞飛，周懿，胡鉞，高政

非牛頓流體在不規(guī)則管道中的流動特性研究＊

李亞飛，周懿，胡鉞，高政

（船舶動力工程技術(shù)交通行業(yè)重點實驗室，湖北 武漢 430063；武漢理工大學 能源與動力工程學院，湖北 武漢 430063）

非牛頓流體被大量應用于工業(yè)生產(chǎn)中，如利用剪切稀化流體應用的醬類食品生產(chǎn)運輸、高分子聚合物在工業(yè)生產(chǎn)中的管道運輸及鉆井用洗井液完井液等。相比傳統(tǒng)的牛頓流體，關(guān)于非牛頓流體在管道內(nèi)的流動的研究還有很大發(fā)展空間，因此對其機理的研究具有重要的現(xiàn)實意義。借助多物理場仿真軟件COMSOL Multiphysics，對非牛頓流體在不規(guī)則管道內(nèi)的流動現(xiàn)象進行了數(shù)值模擬研究，分析了非牛頓流體在管道內(nèi)的速度場分布、剪切應變速率分布以及黏度分布情況。結(jié)果表明，在同一不規(guī)則管道中，管內(nèi)流動的非牛頓Carreau流體的剪切應變速率與黏度關(guān)系符合其本構(gòu)方程特性，相比牛頓流體表現(xiàn)出了明顯較差的流動性，且速度場的分布一直保持一致；實驗組非牛頓流體的剪切應變速率量級比對照組牛頓流體小很多，解釋了兩種流體在速度場分布的差異性；驗證了Carreau流體構(gòu)型中剪切應變速率與黏度的負相關(guān)關(guān)系。

非牛頓流體；管道；流動特性；冪律流體

非牛頓流體是剪切力與剪切應變率之間不是線性關(guān)系的流體。非牛頓流體的應用十分普遍，涉及采用剪切增稠液體的裝甲、食品工廠設(shè)計果醬類食品的運輸加工及高分子聚合物在化工生產(chǎn)中的運輸?shù)戎T多方面。這類非牛頓流體在工業(yè)生產(chǎn)方面的應用研究，直接影響工業(yè)生產(chǎn)效率及國家經(jīng)濟，因此對非牛頓流體在管道中的流動機理具有較大的研究價值。

對非牛頓流體的研究方法包括實驗法、解析解法與數(shù)值解法。實驗法直接，但成本高且實驗普遍性不高。解析解法理論上為最理想的研究方法，通過建立合適的微分方程，使用純數(shù)學方法得出方程的精確解；但對于非牛頓流體復雜的流動情況來說，求解的難度過大。數(shù)值解法則是應用計算機將物理場離散化，之后將流體微分方程組轉(zhuǎn)化為代數(shù)方程并求出各個節(jié)點上的參數(shù)值，屬于一種近似解法。由于數(shù)值解法容易獲得且能保證足夠的求解精度，其已經(jīng)成為研究非牛頓流體問題最為常用的方法。

本文選取一種注塑噴嘴模型作為研究對象，以COMSOL Multiphysics這一軟件作為計算平臺，采用有限元法對非牛頓流體流經(jīng)不規(guī)則管道時表現(xiàn)出的流動特性進行數(shù)值模擬研究。在不規(guī)則腔體內(nèi)，對照研究非牛頓Carreau流體與冪律牛頓流體的流動特性，討論非牛頓流體的剪切應變速率與黏度對于速度大小與速度場分布的影響。

1 管內(nèi)流動數(shù)學模型

1.1 控制方程組

一切的流體流動過程，都以三個基本的物理學原理為基礎(chǔ)：質(zhì)量守恒定律，牛頓第二定律與能量守恒定律。將這些物理學用于構(gòu)建流動模型，將會導出一組方程，即連續(xù)性方程、動量方程與能量方程。這些方程是上述物理學原理的數(shù)學描述，本文不討論傳熱，因此不引入能量守恒方程：

質(zhì)量守恒定律（連續(xù)性方程）：

式（1）表示瞬態(tài)三維可壓流體流動的連續(xù)性方程，本文所分析的流體流動處于穩(wěn)態(tài)且不可壓縮，密度不會隨著時間的變化而改變，因此流體流動的數(shù)學描述為：

動量守恒方程：

式（3）（4）（5）是對于任何流體都成立的動量守恒方程，是微元體內(nèi)流體動量對于事件的變化率等于外界作用于該微元體上的各種力的和，簡稱動量方程，也稱納維斯托克斯（N-S）方程。

1.2 本構(gòu)方程

本構(gòu)方程是反映物料宏觀性質(zhì)的數(shù)學模型，又被稱為流變狀態(tài)方程或是流變方程。在流變學中，本構(gòu)方程是在某些假定條件下，對流體或彈性體的材料力學行為的數(shù)學描述，可用來區(qū)分流體類型。本構(gòu)方程與連續(xù)性方程、運動方程一起構(gòu)成封閉的方程組，用于求解流體的流動特性。

在本文仿真中使用的冪律流體的本構(gòu)方程式如下：

式（6）中：為黏稠系數(shù)，表示物料的黏稠程度，Pa·sn；為冪律流變指數(shù)（簡稱冪律指數(shù)），為無量綱量，表示非牛頓流體的流動特性偏離牛頓流體的程度（=1時為牛頓 流體）。

1.3 非牛頓流體管道流動模型

剪切應變速率描述的是流體的剪切流動，定義為單位時間的剪切應變變化：

式（7）中：為剪切應變。

值得注意的是，剪切應變速率常與速度梯度混淆。實際上二者是不同的概念。速度梯度是流體的速度對空間坐標的導數(shù)，用/來表示。在數(shù)學上，二者的數(shù)值有時相等，這是因為一般速度梯度符合：

但二者的物理意義并不相同，且數(shù)值上有時并不相等（如流體在同軸圓筒之間的流動，此時有角速度的影響）。

2 物理模型及計算方法

在規(guī)則管道內(nèi)，流體的速度場分布較為均勻，各流動特性參數(shù)的變化梯度較小。本文在這一部分選取一不規(guī)則管道腔體進行管道流動特性仿真，以更清楚地展示流體的流動特性。另選取適當?shù)姆桥ｎD流體與牛頓流體進行變量控制，構(gòu)成實驗組與對照組。通過對比兩種流體的仿真結(jié)果，說明非牛頓流體不同于牛頓流體的流動特性。

2.1 物理模型

2.1.1 幾何模型

本組仿真的計算區(qū)域如圖1所示，為一個呈圓柱對稱的不規(guī)則彎曲管道腔體。流體由腔體的圓形入口處流入，中途流經(jīng)一截面當量直徑較小的喉口，最后由腔體的環(huán)形出口處流出。腔體的總長度為42 mm，圓形入口直徑為18 mm，環(huán)形出口寬度為12 mm，其余具體幾何尺寸如圖2所示。A—B為管道進口段；F—G為管道出口段；B—C—D—E—F 與 G—H—I—J—A為管道壁面。使用此不規(guī)則腔體進行數(shù)值模擬的原因：在腔體喉口處由于其的截面當量直徑相對較小，流體經(jīng)過時會出現(xiàn)較大的壓力梯度，進而顯著地影響速度、剪切率與動力黏度。所以在喉口當量直徑最窄處選取一個截面，計算分析兩種流體各項流動參數(shù)在這個截面上的分布，以更精確地說明流體的流動特性。截面具體位置如圖3所示。喉口處截面的位置與圖2中R12與R3圓心連線重合。

圖1 不規(guī)則管道腔體示意圖

圖2 不規(guī)則管道腔體的具體尺寸（單位：mm）

圖3 喉口處截面的位置（單位：mm）

2.1.2 流體性質(zhì)

2.1.2.1 非牛頓流體性質(zhì)（實驗組）

本組仿真計算所使用的非牛頓流體為線性聚苯乙烯在1-氯苯中的溶液（A solution of linear polystyrene in 1-chloronaphthalene），其被廣泛用于塑料制品制造中。該溶液在圖1中不規(guī)則腔體內(nèi)的流動過程，可視為對塑料制品生產(chǎn)中注塑工序的溶液流經(jīng)注塑噴嘴工段過程的模擬。該聚苯乙烯溶液屬于非牛頓Carreau流體，其黏度eff與流體剪切應變速率相關(guān)，具體關(guān)系式如下：

式（9）中：0為零剪切應變速率時的黏度；inf為無限大剪切應變速率時的黏度；為松弛時間，s；為模型指數(shù)。

本溶液的相關(guān)非牛頓Carreau流體模型參數(shù)值與仿真需要的其他物性如表1所示。對照組牛頓流體的物性參數(shù)值如表2所示。

2.1.2.2 牛頓流體性質(zhì)（對照組）

根據(jù)控制單一變量的原則，本組仿真所使用的對照組牛頓流體只進行黏度相關(guān)系數(shù)的調(diào)整，其他物性參數(shù)與實驗組非牛頓流體保持一致。對照組牛頓流體使用冪律流體模型，冪律指數(shù)=1。當冪律指數(shù)=1時，流體的黏度eff實際上固定為1 Pa·s。

表1 線性聚苯乙烯在1-氯苯中的溶液的物性參數(shù)值

參數(shù)名參數(shù)值 零剪切應變速率時黏度μ00 無限大剪切應變速率時黏度μinf /（Pa·s）166 松弛時間λ/s1.73×10-2 Carreau模型指數(shù)n 0.583 密度ρ/（kg/m3）450

表2 對照組牛頓流體的物性參數(shù)值

參數(shù)名參數(shù)值 稠度指數(shù)m1 冪律指數(shù)n1 最低剪切應變速率/s-10.01 密度ρ/（kg/m3）450

2.2 邊界條件

進口邊界條件：進口為圖2中A—B段。設(shè)定進口壓力in，取值從10 kPa開始，以40 kPa的梯度遞增至210 kPa，模擬不同進口壓力下實驗組與對照組流體的流動特性。

出口邊界條件：出口為圖2中F—G段，設(shè)定出口壓力為0。

壁面邊界條件：壁面為圖2中B—C—D—E—F與G—H—I—J—A段。設(shè)定無滑移邊界條件，即壁面處的流體速度為0。

無黏性應力條件：·[（▽+（▽）T）]。

其他：流體不可壓縮；所有管道區(qū)域內(nèi)實驗組非牛頓流體處于層流狀態(tài)，進口壓力較大時，對照組牛頓流體會出現(xiàn)湍流。

2.3 求解器介紹

本次模擬使用的COMSOL Multiphysics是一款功能強大的多物理場仿真軟件，可以對多個領(lǐng)域的物理過程進行模擬計算。

本文研究的是非牛頓流體在管道內(nèi)的流動特性，因此選用的是COMSOLMultiphysics中的Laminar flow求解器模塊。這個模塊內(nèi)置了許多流動的基本物理參數(shù)，可以在這些參數(shù)的基礎(chǔ)上結(jié)合前文所述的數(shù)學模型定義更多仿真所需要的變量。

2.4 數(shù)值計算方法

本文涉及的數(shù)值計算使用的是有限元分析法（Finite element method，F(xiàn)EM）。它的基本思路為：一個物體或系統(tǒng)被分解為由多個相互聯(lián)結(jié)、簡單、獨立的點組成的幾何模型。在這種方法中這些獨立的點的數(shù)量是有限的，因此被稱為有限元。由實際的物理模型中推導出來的平衡方程式被使用到每個點上，由此產(chǎn)生了一個方程組，這個方程組可以用線性代數(shù)的方法來求解。

3 仿真結(jié)果分析

3.1 兩種流體的速度場對比

當進口壓力分別為10 kPa、90 kPa、170 kPa時，從實驗數(shù)據(jù)分析可知，隨著進口壓力的增大，牛頓流體的速度場分布情況發(fā)生了較為明顯的變化。進口壓力為10 kPa時，流體速度在腔體入口段達到最大值，喉口處的中心速度也達到了局部峰值，喉口段之后的流體速度漸緩。進口壓力增大時，腔體入口段的速度場分布情況差異不大，但喉口處的局部速度峰值與入口段的速度最大值之間的梯度差異會逐漸減小，最終在喉口處達到全局速度最大值。在進口壓力為90 kPa時，速度入口段的速度最大值約為 57.669 m/s，喉口處速度峰值則在52.111 m/s附近；進口壓力為170 kPa時，速度入口段的速度最大值約為82.369 m/s，喉口處速度峰值則已經(jīng)達到85.698 m/s。此外，喉口后部的擴散段速度場分布也隨著進口壓力增大而出現(xiàn)較為明顯的差異。當進口壓力為90 kPa時，喉口處的速度峰值區(qū)域已經(jīng)開始變得狹長，峰值點更接近下方壁面。峰值區(qū)域狹長化和峰值點向下方壁面移動的程度隨著進口壓力進一步上升而強化，這一特點從進口壓力達到170 kPa之后就非常明顯。另一方面，流場內(nèi)流體速度的絕對值也隨著入口壓力的增加出現(xiàn)明顯的上升。速度峰值由10 kPa入口壓力時的2.746 m/s上升至90 kPa入口壓力時的57.669 m/s，在170 kPa時則是82.369 m/s。

不同進口壓力下非牛頓流體的速度場分布如圖4所示。從圖4可以看出，在進口壓力為10 kPa時，非牛頓流體的速度場分布與牛頓流體的分布情況相似；但隨著進口壓力的增加，非牛頓流體的速度場分布并不會出現(xiàn)明顯變化，這與牛頓流體不同。此外，非牛頓流體的速度絕對值也遠小于牛頓流體，進口壓力為10 kPa時速度峰值約為0.015 m/s，進口壓力為90 kPa時速度峰值約為0.157 m/s，進口壓力為 170 kPa時速度峰值約為0.346 m/s。很明顯，實驗組的非牛頓流體流動性相比對照組的牛頓流體要差得多。

造成這種流動性差異的原因可以通過比較兩種流體的剪切應變速率和黏性來得出。非牛頓流體剪切應變速率的分布情況如圖5所示。

從實驗數(shù)據(jù)中可得牛頓流體的剪切應變速率值遠大于非牛頓流體的（3個數(shù)量級的差異），說明單位時間內(nèi)牛頓流體的剪切應變變化量更多，表現(xiàn)為較高的流速；非牛頓流體則相反。

另外，如前文所述，對照組牛頓流體由于冪律指數(shù)=1，所以流體的黏度大小固定為1 Pa·s；而實驗組非牛頓流體的黏度由Carreau模型決定，為一個與剪切應變速率呈負相關(guān)的變量，其具體大小與分布如圖6所示?？梢郧宄氐贸?，實驗組非牛頓流體的黏度遠遠大于對照組牛頓流體的黏度，大部分區(qū)域黏度值達140 Pa·s以上，這也可以解釋非牛頓流體的流動性為何遠不如牛頓流體。

圖4 不同進口壓力下非牛頓流體的速度場分布

圖 5 不同進口壓力下非牛頓流體的剪切應變速率分布

圖6 不同進口壓力下非牛頓流體的黏度分布

3.2 喉口截面處的流動特性對比分析

從之前對兩種流體速度場的分析可以得知，速度、黏度等變量在流動區(qū)域的喉口處均會出現(xiàn)比較明顯的梯度變化。利用流動腔體圓柱對稱的特性，通過二維函數(shù)圖像即可確定各項流動參數(shù)在喉口處截面上的分布。

兩種流體流速沿喉口截面徑向的分布曲線如圖7和圖8所示。

圖7與圖8所反映的速度大小變化情況，與前文圖5反映的喉口處速度場分布相符。隨著進口壓力的增大，喉口處的流體流速也相應增加。圖7中各條速度線的峰值，對應的是實驗中牛頓流體的喉口處局部速度峰值；峰值位置的橫坐標隨進口壓力增大而向左橫移，對應的是實驗時牛頓流體局部速度峰值點區(qū)域隨進口壓力增大的下移。圖8中各速度線峰值橫坐標則幾乎沒有出現(xiàn)橫移的現(xiàn)象，也符合圖4中非牛頓流體不同進口壓力下速度場分布大致相同的情況。

從式（9）與式（10）可以看出，實驗組非牛頓流體的速度分布與剪切應變速率分布呈正相關(guān)，剪切應變速率的分布又與黏度分布呈負相關(guān)。為了驗證這兩組關(guān)系，繪制實驗組流體的剪切應變速率與黏度沿截面徑向的分布曲線，如圖9、圖10所示。

從圖9、圖10可以看出，隨著進口壓力增大，實驗組非牛頓流體在截面上的剪切應變速率曲線逐漸上凹，黏度曲線逐漸下凸。速度曲線與剪切應變速率曲線的變化方向相同，剪切應變速率曲線又與黏度曲線變化的趨勢相反，這就驗證了由式（9）與表1數(shù)據(jù)所表示的剪切應變速率與黏度之間的負相關(guān)關(guān)系。

圖7 牛頓流體在截面上的速度分布

圖8 非牛頓流體在截面上的速度分布

圖9 非牛頓流體在截面上的剪切應變速率分布

圖10 非牛頓流體在截面上的黏度分布

4 結(jié)論

本文通過 COMSOL Multiphysics平臺，使用有限元法與局部優(yōu)化網(wǎng)格技術(shù)，建立了非牛頓流體在不規(guī)則管道腔體內(nèi)的流動模型。通過以上數(shù)值模型，分析了非牛頓流體在管道內(nèi)的速度分布、剪切應變速率分布及黏度分布，并討論了流體流動特性參數(shù)之間的關(guān)系。

具體總結(jié)如下：①隨著進口壓力的增大，對照組牛頓流體的流速不斷提高，速度場的分布中出現(xiàn)狹長狀局部峰值區(qū)域；實驗組非牛頓流體流速一直非常小，表現(xiàn)出很差的流動性，且速度場的分布情況一直保持一致。②在高速管道流動中，非牛頓Carreau流體的剪切應變速率遠小于牛頓流體，進而導致其黏度比牛頓流體的黏度大得多，表現(xiàn)為管道內(nèi)極差的流動性與不隨進口壓力變化而改變的速度場分布。③通過對喉口處非牛頓流體的剪切應變速率與黏度分布曲線的分析，驗證了Carreau流體構(gòu)型中剪切應變速率與黏度的負相關(guān)關(guān)系。

O373

A

10.15913/j.cnki.kjycx.2020.08.009

2095－6835（2020）08－0024－05

國家自然科學基金青年科學基金項目（編號：51806157）

李亞飛（1998—），男，湖北武漢人，本科在讀。

〔編輯：嚴麗琴〕