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

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

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

（船舶動(dòng)力工程技術(shù)交通行業(yè)重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430063；武漢理工大學(xué) 能源與動(dòng)力工程學(xué)院，湖北 武漢 430063）

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

非牛頓流體；管道；流動(dòng)特性；冪律流體

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

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

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

1 管內(nèi)流動(dòng)數(shù)學(xué)模型

1.1 控制方程組

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

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

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

動(dòng)量守恒方程：

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

1.2 本構(gòu)方程

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

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

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

1.3 非牛頓流體管道流動(dòng)模型

剪切應(yīng)變速率描述的是流體的剪切流動(dòng)，定義為單位時(shí)間的剪切應(yīng)變變化：

式（7）中：為剪切應(yīng)變。

值得注意的是，剪切應(yīng)變速率常與速度梯度混淆。實(shí)際上二者是不同的概念。速度梯度是流體的速度對(duì)空間坐標(biāo)的導(dǎo)數(shù)，用/來(lái)表示。在數(shù)學(xué)上，二者的數(shù)值有時(shí)相等，這是因?yàn)橐话闼俣忍荻确希?/p>

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

2 物理模型及計(jì)算方法

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

2.1 物理模型

2.1.1 幾何模型

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

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

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

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

2.1.2 流體性質(zhì)

2.1.2.1 非牛頓流體性質(zhì)（實(shí)驗(yàn)組）

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

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

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

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

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

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

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

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

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

2.2 邊界條件

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

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

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

無(wú)黏性應(yīng)力條件：·[（▽+（▽?zhuān)㏕）]。

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

2.3 求解器介紹

本次模擬使用的COMSOL Multiphysics是一款功能強(qiáng)大的多物理場(chǎng)仿真軟件，可以對(duì)多個(gè)領(lǐng)域的物理過(guò)程進(jìn)行模擬計(jì)算。

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

2.4 數(shù)值計(jì)算方法

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

3 仿真結(jié)果分析

3.1 兩種流體的速度場(chǎng)對(duì)比

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

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

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

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

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

圖4 不同進(jìn)口壓力下非牛頓流體的速度場(chǎng)分布

圖 5 不同進(jìn)口壓力下非牛頓流體的剪切應(yīng)變速率分布

圖6 不同進(jìn)口壓力下非牛頓流體的黏度分布

3.2 喉口截面處的流動(dòng)特性對(duì)比分析

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

O373

A

10.15913/j.cnki.kjycx.2020.08.009

2095－6835（2020）08－0024－05

國(guó)家自然科學(xué)基金青年科學(xué)基金項(xiàng)目（編號(hào)：51806157）

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

〔編輯：嚴(yán)麗琴〕