基于剪切稀化模型的三通減阻數值模擬

張志，張浩，2，?，李茂林，石若冉

(1.山東建筑大學熱能工程學院， 山東 濟南 250101；2.山東省綠色建筑協同創新中心，山東 濟南 250101)

0 引言

在當今環境污染與能源危機的大背景下，節能減排已經成為了我國實現“雙碳”目標的重要方式，近年來有學者認為節約能源已經成為了繼煤炭、石油、天然氣和電氣之后的第五大能源[1]，從側面反映了節約能源對于解決能源危機的重要性。 我國建筑能耗在總能耗中占比約為45%，其中在供熱管路中，管道輸運阻力約占供熱量的30%。 為了減少碳排放、降低能源消耗，有必要減少供熱管路的輸運阻力[2]。 三通作為供暖管路主要的分流和匯流部件，對于增加管道阻力有著重要的影響。 當流體流經三通處時，由于不同速度流體的沖擊混合、流量變化以及壁面影響，會形成彎曲的流線和渦旋區，進而導致流體流動過程中的壓力和能量損失，增加了輸運過程中的能耗[3]。

供熱管道中的減阻方式，可以分為非光滑表面減阻(如溝槽減阻、凹坑減阻、仿生表面減阻)、疏水表面減阻以及添加減阻劑減阻[4]。 其中，以添加減阻劑減阻最為方便、減阻效率最高。 Toms 發現在管道中添加少量的高分子聚合物可以起到減阻效果(稱之為Toms 效應[5])，后來有研究人員發現在管道中添加表面活性劑同樣可以起到減阻的效果。 高分子聚合物具有添加微量便可減阻的特點，對環境影響小，但由于其易降解、不耐剪切，在空調及供暖管路中難以發揮作用，而表面活性劑抗剪切能力強，即使失效再次滿足條件仍可恢復減阻效果，可以更好地適應空調/供熱管路的減阻。 Kotenko 等[6]開展了減阻表面活性劑在區域供冷的實驗研究，選用了針對不同溫度范圍的水溶液開發的兩種減阻產品進行實證，發現兩種減阻率工作溫度范圍在5 ～55 ℃時，其最大減阻率在60%～80%之間。

為了探究表面活性劑溶液減阻機理，不少學者針對其流變特性展開研究。 謝程程等[7]測量分析了混合比例為1 ∶1 的十六烷基三甲基氯化銨(Cetyltrimethyl Ammonium Chloride，CTAC)和水楊酸鈉(Sodium Salicylate， NaSal)溶液流變特性，發現減阻溶液黏度變化與質量濃度及剪切率大小相關。 莫偉南等[8]對CTAC 溶液進行了流變實驗測量，發現減阻溶液剪切黏度隨剪切率變化過程大致分為剪切稀化、剪切增稠和二次剪切稀化3 個階段，剪切稀化模型Carreau－Brid 和黏彈性模型Giesekus 均能較好地擬合形成剪切誘導結構之后的流變特性曲線。

文章以匯流三通為研究對象，采用非牛頓流體模型中的剪切稀化Carreau－Bird 模型，模擬添加表面活性劑CTAC/NaSal 溶液后三通管內的流動現象，對比純水與CTAC 溶液不同雷諾數下模擬結果，通過對減阻率、速度、湍動能變化及渦旋結構變化探究管內減阻特性。

1 數值計算方法

1.1 物理模型

匯流三通幾何模型如圖1 所示，D 為管徑；Inlet1 和Inlet2 分別為兩個進口，Outlet 為出口；a－a、b－b、c－c 分別為進口1、2 及出口上的截面；x、z 為坐標軸方向。 管徑D 設置為20 mm，為了使流動進入充分發展階段，并且出口不影響匯流后的流動，設置Inlet1 進口段距離為6D，Inlet2 進口段距離為4D，Outlet 出口段距離為6D。

圖1 匯流三通幾何模型圖

1.2 數學模型

流體流動的連續性方程、動量方程[9－10]分別由式(1)和(2)表示為

式中ρ 為流體密度，kg/m3；t 為流動時間，s；ui、uj為不同方向的速度，m/s；p 為壓力，N/m2； xi、xj為不同方向的位移，m；μ 為流體黏度，N·s/m2；fi為流體質點的單位質量力，N。

匯流三通中的流動特性比較復雜，存在流線彎曲、渦旋及二次流現象。 相比k 模型以及k－ε 模型，雷諾應力模型(Reynolds Stress Models， RSM)更細致地考慮了流動過程中流線彎曲、旋轉、剪切應力快速變化及局部渦旋現象，可以更加精準地預測復雜流動。

RSM 湍流方程由式(3)表示為

湍動能k 的方程和耗散率ε 的方程分別由式(4)和(5)表示為

式中Pij為剪應力產生項；Gij為浮力產生項，對于不可壓縮流體，其值為零；μt為湍動黏度，N·s/m2；δij為克羅內克函數符號，當i 和j 兩個指標相同時，δij＝1，當指標不同時，δij＝0；經驗常數取默認值，C1＝1.8、C2＝0.6、C1ε＝1.44、C2ε＝1.92、C3ε＝0.09、σk＝0.82、σε＝1。

表面活性劑CTAC/NaSal 溶液黏度模型采用非牛頓流體模型中的Carreau－Bird 模型。 表面活性劑的剪切稀化特性是其減阻的重要特性之一，李恩田[11]、張紅霞[12]分別通過驗擬合出CTAC 溶液部分濃度及溫度下Carreau－Bird 模型參數，并驗證了數值模擬與實驗結果。

Carreau－Bird 模型由式(6)表示為

1.3 參數設置及邊界條件

模擬選擇穩態求解器，采用壓力修正(Pressurebased)的求解方法，通過壓力耦合方程組的半隱式方法(Semi － Implicit Method for Pressure Linked Equations， SIMPLE)耦合壓力及速度，對于收斂殘差值設置，連續性方程設為10－5，其耦合方程設為10－6，并監測截面a－a、b－b、c－c 上的總壓值，當總壓值穩定不變時認為模擬結果滿足要求。

邊界入口條件為速度進口(Inlet)；出口為壓力出口(pressure-outlet)，其值取默認值0；壁面采用無滑移壁面(wall)條件。

計算域中的流體分別采用采用牛頓流體純水與表面活性劑CTAC/NaSal 溶液，其中純水的物性參數為ρ＝998.2 kg/m3、μ ＝1×10－3N·s/m2。 溶液中CTAC 的濃度較小，除黏度變化外，其他物性參照水的物性，黏度大小通過Carreau－Bird 表征，文章假定表面活性劑流動已經進入穩定剪切黏度狀態下，模型相關參數根據文獻[11]得到，取η0＝7.392 2 mPa·s、η∞＝0.586 1 mPa·s、λ＝0.0147 1、n＝0.312 3，描述質量濃度為200 mg/L 的CTAC 溶液。

1.4 網格無關性與模型驗證

網格質量和數量直接決定了數值計算的時間和精度，理論上網格越密，計算結果越精準，但所需計算資源、時間也越多，所以在進行網格劃分時不可能無限制加密網格。 為了確定對計算結果影響最小的網格密度，文章選取了4 種不同尺寸的網格方案，網格數量分別為65 萬、207 萬、309 萬、404 萬(分別對應方案1、2、3、4)。

選取雷諾數為20 001.2 時，CTAC/NaSal 溶液在以上4 種網格方案中的模擬結果，如圖2 所示，網格數較少的方案1 與其余方案相比，在中軸線上靜壓值變化相差較大，方案2、3、4 基本吻合，方案2(207 萬)已滿足計算精度要求。 選取網格方案2，劃分三通的結構化網格，通過O－Block實現對圓柱結構的更好適應，加密管壁邊界層處的網格，邊界層第一層網格厚度大小設置為0.01 mm，網格增長率為1.15，得到的網格數的207 萬，而網格質量則在0.86～1之間。 網格劃分如圖3 所示。

圖2 網格無關性驗證圖

圖3 方案2 局部網格圖

為了驗證所采用的Carreau－Bird 模型參數是否合適，選取管徑D 為20 mm、長度為13D 的直管模型，模擬雷諾數為20 000 時純水和CTAC 溶液的管內流動，通過比較純水與CTAC 溶液在湍流近壁面流場參數變化來判斷是否發生減阻效應。

采用無量綱離壁距離與無量綱平均軸向速度表達方法，其中無量綱離壁距離y＋的計算公式由式(7)和(8)表示為

式中y 為距離壁面長度，m；uτ為摩阻流速，m/s；τw為壁面平均剪切應力，Pa； Δp 為進出口壓差，Pa；L為直管長度，m。

無量綱平均軸向速度U＋由式(9)表示為

式中U 為直管內流速，m/s。

無量綱平均軸向速度與無量綱離壁距離關系曲線如圖4 所示。 對于牛頓流體，迪恩(Dean)平均軸向速度分布曲線計算式由式(10)和(11)表示為

圖4 無量綱平均軸向速度分布曲線圖

對于減阻溶液，高分子聚合物在對數區的Virk平均軸向速度漸近線計算式[13]由式(12)表示為

由圖4 可知，純水的平均軸向速度在邊界層上的變化，比較符合Dean 曲線，在黏性底層與對流區軌跡基本吻合，CTAC 溶液在黏性底層區域的變化和純水變化相同，隨著y＋的增加，CTAC 溶液的無量綱平均軸向速度上升變緩，略小于Virk 曲線，位于Virk 曲線與Dean 對數區曲線之間。 在黏性底層和對數區，CTAC 溶液無量綱平均軸向速度曲線斜率大于純水，在對數區曲線向上偏移，與添加劑減阻規律一致，因此可以認為該模型參數能夠模擬出表面活性劑在管道內的減阻效應。

1.5 關鍵參數計算

根據伯努利方程，可以得到不同管路之間的局部阻力系數計算公式，分別由式(13)和(14)表示為

式中ξ01、ξ02分別為a、c 管段之間與b、c 管段之間的局部阻力系數；p1、p2、p3為在截面a、b、c 的平均靜壓，Pa；u1、u2、u3為a、b、c 管段的平均流速，m/s；λ1、λ2、λ3分別為a、b、c 管段的摩擦阻力系數[14]；L1、L2、L3分別為截面a、b、c 到原點的距離，mm。

牛頓流體純水的主管雷諾數Re 由式(15)表示為

對于表面活性劑溶液，黏度取決于剪切速率，無法直接應用牛頓流體的計算方法，因此采用文獻[15]的方法，通過特征剪切速率下的黏度計算，由式(16)和(17)表示為

減阻率RD的計算由式(18)表示為

式中ξ、ξ′分別為純水和CTAC 溶液的局部阻力系數。

2 結果與分析

2.1 減阻率規律

三通管路與直管路在結構上相差較大，在加入CTAC 減阻溶液之后可能存在不同的減阻規律，比如兩者減阻變化過程不同、減阻最值的不同以及減阻最值發生的雷諾數不同。 通過模擬同等管徑及長度下的三通與圓直管，得到減阻率隨Re 的變化規律，如圖5 所示。 可以發現三通與直管的減阻率變化趨勢是一致的，呈現隨雷諾數先增大后減小的規律。 但兩者的最大減阻率以及發生最大減阻的雷諾數范圍存在差異，直管最大減阻雷諾數約在40 000，最大減阻率＞45%，這個結果與實驗數據[9]是一致的；三通減阻最大值雷諾數約在10 000，并且僅能達到約30%。 相對于直管，三通中可能存在部分區域，而添加CTAC 溶液可能無法完全發揮其減阻效果，因而導致減阻效果不佳。

圖5 CTAC 溶液減阻率隨Re 數變化圖

根據減阻率的變化情況，可以將減阻過程分為3 個區域，分別為不完全減阻區Ⅰ(Re＜第一臨界雷諾數Re1，減阻率＜0)、完全減阻區Ⅱ(Re1≤Re≤第二臨界雷諾數Re2，減阻率＞0)以及過減阻區Ⅲ(Re＞Re2，減阻率＜0)。 圓直管在雷諾數＜12 000，減阻溶液即處于不完全減阻區Ⅰ，此時減阻率＜0，對于整個管路來說，添加減阻劑增加了阻力，隨著雷諾數的減小，管道增阻明顯。 這主要是由于處于此雷諾數區間，剪切應力很小，沒有形成剪切誘導結構，無法發揮其減阻效應。 三通管不存在不完全減阻區Ⅰ，由于三通流體間的沖擊與壁面分離，即使是在較小的雷諾數下，仍能形成剪切誘導結構。 在稍微增加雷諾數的情況下，便達到最大的減阻效果，隨雷諾數增加，達到過減阻區域Ⅲ，溶液中的剪切誘導結構被破壞，導致不再減阻。

2.2 流速變化規律

2.2.1 軸向速度分布

雷諾數為20 001.2 時，主管中軸線上的純水和CTAC 溶液軸向速度變化如圖6 所示，其中橫坐標x通過管道直徑D 做無量綱化，縱坐標軸向速度U 通過斷面平均流速Ub做無量綱化。

圖6 Re＝20 001.2 時軸向速度變化圖

隨著進口距離不斷增加，速度值逐漸增加，直至三通支管處，主、支管發生匯流，流量增加，導致區域內速度急劇上升；同時由于兩股流體合并，伴隨著劇烈的湍流運動，能量耗散，流動阻力增加，在匯流之后形成渦旋區域，導致速度不斷降低，速度下降到一定值后穩定在原始速度的2 倍。 比較純水和CTAC溶液可以發現，CTAC 溶液的無量綱化速度始終處于純水之上，這表明CTAC 溶液起到了提高流動速度的作用；在出口段的無量綱速度值，純水處于上下波動，表現為較強的湍流狀態，CTAC 溶液較為恒定，表現為穩定的層流狀態。 出口段速度分布，CTAC 溶液相比于純水的出口速度分布較為均勻，并未出現速度波動分布現象，表明添加CTAC 溶液可以促使出口段流動從湍流狀態向層流狀態轉變。

2.2.2 流速等值線分布

Re＝6 000.2 時，三通管路對稱中心面(y ＝0)上純水和CTAC 溶液的速度分布云圖如圖7 所示。 匯流三通主要阻力來源于兩股不同流速的流體混合紊流沖擊產生的損失，以及匯流之后在內管壁分離形成的回流區，該回流區通常伴隨著強烈的渦旋現象。在進口1 和2，流體以0.5 m/s 的速度進入管路，在匯合處左側靠近壁面區域，產生一個低速區，這是流體偏轉的向心力導致的。 在匯合處，兩股流體相遇，速度增加，分層明顯。 在匯合后管道內側壁面(以靠近支管側為內側)附近產生了一個渦旋區，在外側壁面處產生一個高速流動區，這主要是由于流體偏轉脫離壁面所引起的，渦旋區的大小在一定程度上反映了三通的阻力。

圖7 Re＝6 000.2 時，y＝0 平面上的速度分布云圖

從整體的流速分布來看，CTAC 溶液減小了低速區的范圍，增大了管道外側高速流動區的范圍，改變了三通管內流速分布，減少了能量耗散，實現了減阻，但也增大了渦旋中心區的速度，提高了該區域的湍流運動，可能會使阻力增加。

2.3 湍動能變化規律

2.3.1 沿程湍動能變化規律

為了研究三通內沿程湍動能變化規律，分析不同雷諾數下純水與CTAC 溶液在主管軸線上各斷面的沿程湍動能變化，結果如圖8 所示。 可以看到，在相同雷諾數下，多數截面CTAC 溶液的湍動能更小，并且在減阻率較大的雷諾數下(20 000～40 000)，兩者湍動能差值更大。 在進口段(x ＝－0.15 ～0 m)CTAC 溶液相對于純水湍動能降低約為50%～75%，說明減阻溶液實現了對湍流的較好的抑制作用，處于更好減阻雷諾數的流動中，對湍流運動的抑制越明顯。 而在圖8(a)較小雷諾數和圖8(d)較大雷諾數下，純水和CTAC 溶液的湍動能差別較小，在進口段湍動能僅降低了10%～25%，這正好與CTAC 溶液在三通管路的減阻率隨雷諾數變化相對應，較大或較小的剪切力都會降低減阻劑的減阻效果。

圖8 不同雷諾數下純水和CTAC 溶液沿程湍動能分布圖

由圖8(b)可以看出，在x ＝0 m 處湍動能急劇增加，此時CTAC 溶液湍動能相對于純水湍動能提高約8.5%，由于匯流與壁面分離效應，流動的剪切應力較大，超出CTAC 溶液減阻范圍，添加減阻劑導致湍動能增加，隨著流動的發展，遠離渦旋區后，減阻劑逐漸恢復減阻效果。

2.3.2 徑向湍動能分布規律

由于在軸向湍動能分析時，x ＝0 m 處(流動渦旋區)CTAC 溶液的湍動能值大于純水湍動能值，因此選取Re＝20 001.2 時純水和CTAC 溶液分別在截面x ＝0.02、0.03、0.04、0.05 m 處徑向湍動能分布圖，分析徑向湍動能變化及在不同截面上徑向湍動能發展，如圖9 所示，其中d 為徑向方向上的坐標。 由圖9(c)可以發現湍動能最大值時發生在管道中心靠內側的位置(d/D ＝－0.1 處)。 外側壁面(d/D ＝0.5 處)純水的湍動能大于CTAC 溶液；隨著離壁面距離增加，湍動能不斷降低，此時純水和CTAC 溶液的湍動能相差較小；接近渦旋區(d/D ＝0 附近)時，湍動能急劇增加，此時CTAC 溶液的湍動能是大于純水的；離開渦旋中心區后(d/D＜－0.2)，湍動能大小表現為CTAC 溶液小于純水。 由此可知，在高速流動區(d/D≥0.2)，CTAC 溶液起到了減阻的效果，并且減阻主要發生在近壁面附近；在渦旋中心區(－0.2≤d/D＜0)，表現為增阻，此時加入CTAC 導致流動的湍動能增加，渦旋區的剪切應力過大，CTAC分子無法形成剪切誘導的網狀結構，而是無序的分布在溶液中起到阻礙作用。

圖9 Re＝20 001.2 時純水和CTAC 溶液徑向湍動能分布圖

在內側壁面處(d/D＜－0.2)，CTAC 溶液與純水的湍動能處于發展階段，隨著流動的進行，CTAC 溶液的湍動能逐漸小于純水，減阻劑恢復減阻效應。由此可知，添加CTAC 時，減阻主要發生在壁面和遠離渦旋中心區的近壁面，而渦旋中心區增加了阻力，由于整個流動中渦旋中心區的范圍較小，且增阻效果不明顯(湍動能增值僅為8.5%)，因此添加CTAC整體表現為減阻效果，但達不到直管段中的減阻率。

2.4 渦旋結構分析

2.4.1 渦量云圖

流體流動過程中，渦旋結構的產生會增加流動阻力和能量耗散，向純水中添加表面活性劑溶液進行減阻，同樣會對流動過程中的渦旋結構產生影響。通過分析渦旋結構變化，可以對添加CTAC 減阻機理有更深的理解。 因此，文章采用第二代渦識別方法中的Ω 準則[15]，計算了當Re＝20 001.2 的純水和CTAC 溶液的渦量，得到的渦量云圖如圖10 所示。

圖10 純水和CTAC 溶液局部渦量云圖

比較分析軸向及徑向湍動能大小，可以發現位于匯流之后的渦旋區，CTAC 溶液減阻效果并不好，因此渦量云圖選取的是三通管內流動的下半段區域。 通過對比渦量分布，可以發現相比于純水，CTAC溶液的較大值渦量在外側(上)壁面與內側(下)壁面區域范圍變小，在渦旋區與高速流動區交界處附近變化不大，而在渦旋中心區域范圍增大，這種變化與之前所分析的徑向湍動能變化是一致的。

2.4.2 Q 準則渦旋結構

湍流結構通常通過具有旋轉運動的渦旋結構表示。 渦旋結構的旋轉大小可以通過Q 準則表示[15]。Q 準則能較好的表達渦大小，但容易受到壁面剪切層影響，同時Q 閾值的選取也會影響結構的顯示。根據Q 值特征及三通管內流動特性，選擇了Q 分別為200、500、1 000 的等值面，用速度云圖體現，得到的Re＝20 001.2 下的純水與CTAC 溶液局部渦旋結構如圖11 所示。

圖11 Q＝200 等值面下流體速度分布圖

兩種溶液的速度分布都比較符合匯流之后的流動軌跡，選擇Q＝200 的等值面來代表渦旋結構是可行的。 從圖11 可以得到，在三通管內流動的下半段，渦旋結構主要出現在內側壁面(靠近支管側)。外側壁面流動主要以剪切為主，旋轉為次，因此在管道外側壁面處并未出現或少量出現渦旋。 大的渦旋結構出現于流動匯流后，尤其是靠近內側壁面附近。匯流初始渦旋結構較大，隨著流體流動發展，大尺度渦旋結構逐漸破碎分離，形成小尺度渦旋結構，最終在靠近出口處渦旋結構幾乎消失。 通過對應速度可以發現，渦旋結構較大的區域有著較大的速度。 與CTAC 溶液相比，純水整體渦旋結構范圍更大，同時發展長度更長，渦破碎過程中，純水中渦破碎分離距離更遠，并且小尺度渦的數量下降顯著。 在渦旋中心區，兩者的渦旋結構差別不大，管內加入CTAC 溶液在此處并未起到抑制渦旋大小的作用，因而在較大剪切力的渦旋中心區，CTAC 溶液表現的特性與牛頓流體純水并無差異。

圖12(a)是當Q 值為500 時的流體速度分布，其整體渦旋結構范圍小于Q ＝200 的結構(如圖11所示)，這說明的旋轉強度有所下降。 在渦旋中心區，圖12(a)中純水及CTAC 溶液渦旋結構大小與圖11 差異不大，這是由于流動匯流區整體剪切力較大，在強剪切力下，Q 值選取具有更大的范圍，其對閾值變化敏感度降低。 而在渦破碎分離區，圖12(a)中純水及CTAC 溶液渦旋結構大小與圖11 呈現較大差異，比較純水渦旋結構變化，可以發現其密集程度明顯降低，流向渦與展向渦數量在一定程度上減少。 對于低剪切力的小尺度渦旋，Q 準則捕捉能力較弱，對閾值的依賴性更大。 同樣的，圖12(b)在選取Q 值為1 000 時，渦旋結構更小且更加分散。

圖12 Q＝500 及Q＝1 000 等值面下流體速度分布圖

盡管對于小尺度渦旋顯示的范圍不同，但通過對比相同Q 值下純水和CTAC 溶液渦旋結構大小，3 種Q 值等值面圖也呈現了相同的規律。 由此可知，加入CTAC 溶液影響了管內渦旋結構，降低了小尺度渦旋的產生和渦發展長度，尤其對于出口段的渦旋結構，有效地起到了抑制作用，使其朝層流方向轉變。

3 結論

文章通過CFD 軟件模擬了CTAC 溶液的剪切稀化特性在三通管內產生的減阻效應，通過將CTAC 溶液與同工況純水模擬結果相比較，得到以下結論:

(1) 與直管相比較，三通最大減阻率較低，在低雷諾數以下不存在減阻率＜0 的區域，隨雷諾數增加減阻率在達到峰值后，很快降至接近零。

(2) 添加CTAC 會改變三通管內流速分布。CTAC 溶液中的低速區范圍變小，高速流動區范圍變大，同時渦旋中心區存在速度增加的現象。 相同雷諾數下CTAC 溶液的主管軸向無量綱速度大于純水，在流動出口段CTAC 溶液的速度更加穩定，出口段湍流流動趨近于層流化。

(3) 湍動能降低并非發生在整體區域，對于渦旋中心區，CTAC 溶液湍動能是增加的，在遠離渦旋區及近壁面處小于純水的，在渦旋區內側不同截面處呈現不同的變化。

(4) CTAC 溶液渦旋結構發生改變。 溶液中整體渦旋結構變小，小尺度渦數量下降，壁面以及近壁面處渦量減少，渦旋中心區的渦量增加，其減阻效應主要發生在壁面附近和遠離渦旋中心區的近壁面。