管道不可壓縮流動壓力脈動(偽聲)的數(shù)值研究

蔡建程， 謝新俊， MOCHALIN Ievgen， BRAZHENKO Volodymyr， 鄂世舉

(1. 浙江師范大學 工學院，浙江 金華 321004；2. 浙江省城市軌道交通智能運維技術與裝備重點實驗室，浙江 金華 321005)

管道系統(tǒng)出現(xiàn)在眾多工程應用中，管道振動噪聲造成管道系統(tǒng)的安全隱患，廣受人們關注。引發(fā)管道振動的一個因素是主機動平衡差以及基礎設計不當，而實際生產(chǎn)情況表明，引發(fā)管道振動另一個主要因素是流體壓力脈動[1]。由于間歇加壓，管道內(nèi)壓力在平均值上下脈動，而在管道彎曲部位，壓力脈動會產(chǎn)生相應的隨時間變化的激振力，引發(fā)振動與噪聲。

在經(jīng)典流體力學中，學者們在直管及彎頭的水力損失方面進行了廣泛的研究。對于彎管，彎頭處存在二次流、射流-尾跡等現(xiàn)象的研究，具有較高的學術研究價值。近年來，管道氣流脈動研究越來越熱，主要集中在基于管道聲學的傳遞矩陣法和基于計算流體動力學(computational fluid dynamics，CFD)的數(shù)值研究。

在管道聲學壓力的傳遞矩陣法研究方面，賀尚紅等[2]利用傳遞矩陣法建立液壓管網(wǎng)壓力-流量傳遞函數(shù)模型，研究液壓系統(tǒng)受泵源周期性吸排油機制影響產(chǎn)生流量脈動，流量脈動遇到系統(tǒng)負載阻抗后形成壓力脈動。Cyklis等[3-4]也進行了類似的工作，用以研究壓縮機管道系統(tǒng)氣流脈動。聲學傳遞矩陣法的優(yōu)點是計算量小，通過把長管道系統(tǒng)分解成四端參數(shù)聲學模型，進行快速預測。但當壓力脈動幅度超過8%時，傳遞矩陣法將產(chǎn)生較大誤差[5]，所以有必要進行管道內(nèi)壓力脈動的CFD數(shù)值研究。

CFD方法能更為詳細地得到管道內(nèi)流場，特別是在彎管流動方面。對于彎管流動，彎頭轉彎半徑和管道直徑之比Rc/D是重要的無量綱參數(shù)。Tanaka等[6]利用大渦模擬(large eddy simulation，LES)結合壁面函數(shù)求解Rc/D=2的彎管流動，分析彎頭處的非定常流動。Tan等[7]基于特征線分離算法的大渦模擬，計算了Rc/D=1和Rc/D=2兩種情況的彎管流動，結果顯示彎頭曲率對于彎頭橫截面上的壓力、速度以及二次流分布有很大影響。Dutta等[8]則用k-ε湍流模型求解Rc/D=2的彎管流動，研究分離點、再附點隨雷諾數(shù)Re的變化情況。陳敏等[9]利用FLUENT軟件的RNGk-ε湍流模型計算了不同雷諾數(shù)下90°彎管內(nèi)的二次流。李靜等[10]也利用相同方法數(shù)值研究了雷諾數(shù)和半徑比Rc/D對90°彎管流動的局部阻力的影響。魏志等[11]利用LES方法對閥體后90°彎管內(nèi)部流場進行了三維數(shù)值模擬，獲得了蝶閥在全開狀態(tài)下彎管內(nèi)不同截面的流場特性和不同徑向截面的渦旋結構，并與粒子圖像測速試驗結果進行了對比。梁開洪等[12]則研究了入流角對90°圓弧彎管流內(nèi)的二次流、分離、偏流及流動阻力的影響。

在管道壓力脈動研究方面，韓文龍等[13]利用ANSYS FLUENT軟件建立了管道系統(tǒng)流動的CFD模型，提出了合理的邊界條件，分析了往復式壓縮機管道系統(tǒng)氣流脈動的數(shù)值，并與試驗比較，表明標準k-ε湍流模型在計算管道系統(tǒng)氣流脈動時最為準確。馬斕擎[14]利用ANSYS CFX軟件模擬了進口質(zhì)量流量脈動時彎管和閥門下游的壓力脈動特性。

學者們利用CFD模擬管內(nèi)流動時經(jīng)常使用不可壓縮流動假設，這是因為依據(jù)流體力學理論，對于Ma<0.3的流動通常看成不可壓縮(流體密度ρ為常數(shù))。對于不可壓縮流動，水動力壓力脈動為聲學理論的偽聲，管內(nèi)水動力壓力脈動特性的研究較少。本文先理論分析不可壓縮流場壓力脈動(偽聲)與聲波方程的區(qū)別，再利用CFD數(shù)值研究直管與彎管入口處設置簡諧水動力壓力脈動條件下，管內(nèi)非定常流動的特性。

1 管道流動的數(shù)學模型

1.1 不可壓縮流動數(shù)學模型

不可壓縮流動Navier-Stokes方程

?·v=0

(1)

(2)

式中，v，ρ，p，μ為流體速度、密度、壓力及動力黏度。該連續(xù)方程和動量方程可以組成封閉的方程組，與相應的邊界條件構成數(shù)學上的適定問題，而不需要求解能量方程。不可壓縮流動壓力p沒有以與時間相關的形式出現(xiàn)，這歸因于連續(xù)方程的非發(fā)展特點。

式(1)對時間求偏導后減去式(2)的散度，得到

?2p=-ρ?·[(v·?)v]

(3)

可見，給定速度場時其壓力場的控制方程為Poisson方程，壓力與速度以橢圓型方程耦合。壓力場的Poisson方程意味著擾動的傳播速度為無窮大，即整個壓力場同步變化，不同位置點擾動無時間差。對于一維模型，式(3)簡化為

(4)

式中，u為x方向的流體速度。

1.2 一維管內(nèi)可壓縮流動模型(小擾動傳播)

此處推導管內(nèi)可壓縮壓力脈動的波動方程。一維的連續(xù)方程和動量方程分別為

(5)

(6)

根據(jù)Stokes假設有

λ=-2μ/3

(7)

假設波動為等熵過程，有

(8)

式中：κ為絕熱指數(shù)；c為聲速。連續(xù)方程式(5)轉化為

(9)

對于管道里流動，把物理量分解成平均量和脈動量之和φ=φ0+φ′，認為脈動量與平均量相比是小量φ′?φ0，把脈動量及其對自變量x，t的導數(shù)視為一階微小量。式(6)對x求導，并且忽略流體黏性，得到

(10)

式(9)按φ=φ0+φ′變成

(11)

然后分別對x和t求導，得到

(12)

(13)

將式(12)、式(13)代入式(10)，得到

(14)

考慮到u0相比于聲速c很小，忽略u0項，得到經(jīng)典聲波方程

(15)

其行波解形式為f(x-ct)，由于不考慮黏性，聲壓不隨時間和距離衰減。周紅等[15]把黏性力項(式(6)右端第二項)進行簡化建模：把速度u的二階導數(shù)項簡化為阻尼系數(shù)α和速度u相乘，并且用分離變量法得到解析解，其中包括隨時間衰減項exp(-αt/2ρ0)。

式(15)的行波解f(x-ct)表明擾動傳播速度為c，當c→∞擾動對時間的求導項可以忽略，此時變成不可壓縮流動的壓力脈動情況。

2 管道不可壓縮流動的數(shù)值計算

2.1 CFD模型

本文以管道內(nèi)壓力脈動為研究對象，分析其不可壓縮流動的壓力脈動特性。流體介質(zhì)為不可壓縮空氣，密度為1.225 kg/m3，動力黏度為1.789 4×10-5Pa·s。所研究的直管直徑D=0.190 m、長度為7.0 m。彎管由轉彎半徑為Rc=0.11 m的90°彎頭連接兩根直徑D=0.190 m、長度為3.5 m直管構成，Rc/D=0.579。

直管和彎管的示意圖，如圖1所示。對于直管，沿管道中心線以及上部表面附近(離壁面0.2 mm處的邊界層內(nèi))，以2.0 m為間距設置流場計算的壓力及速度監(jiān)測點，從進口側開始編號1～34(進出口處因為設置邊界條件，所以不設置監(jiān)測點)。彎管的前后兩段直管內(nèi)也類似地設置監(jiān)測點，編號與直管一致，另外在彎頭中心線以及上側、內(nèi)側、外側表面附近的邊線上各設置7個監(jiān)測點，參見圖1(c)的局部放大圖，以研究彎頭處的流場脈動。對于彎管設置流向中間剖面和縱向中截面，用于分析內(nèi)部流動。對于直管，管內(nèi)流動基本上呈軸對稱，所以僅設置流向中間剖面，參見圖1(d)。

圖1 所研究的直管和彎管

表1 計算模型網(wǎng)格統(tǒng)計

圖2 圓管截面計算網(wǎng)格分布

SSTk-ω湍流模型的湍動能k和耗散率ω的對流-擴散方程為

(16)

(17)

式中：Gk，Gω為k，ω的產(chǎn)生項；Yk，Yω為耗散項；Sk，Sω為用戶自定義源項；Dω為交叉-擴散項；擴散系數(shù)Γk=μ+μt/σk；Γω=μ+μt/σω；σk和σω為k，ω湍流Prandtl數(shù)；μt為湍流黏性系數(shù)，定義為

(18)

式中：S為應變率張量； 系數(shù)α1=0.31；α*系數(shù)用于衰減湍流黏度；F2為第二混合系數(shù)，它在邊界層內(nèi)取值為1而在自由剪切層內(nèi)取值為0。

2.2 邊界條件及計算方法

進口處施加隨時間正弦變化的總壓力邊界條件：p=131+23.6sin(3 644.24t)，基頻為580 Hz。該邊界條件模擬某臺具有12個葉片轉速為2 900 r/min的風機出口壓力脈動。風機內(nèi)部旋轉葉片產(chǎn)生的流場脈動，風機出口處的壓力脈動包括水動力壓力脈動和聲學壓力脈動，而前者為主要分量。管道進口的湍流條件設為湍流強度5%，水力直徑為0.19 m。出口設置為壓力出口條件，平均靜壓為0。管道壁面使用無滑移壁面邊界條件。

由于彎管流動的壓力損失比直管大，所以彎管流量比直管流量小。計算得到直管內(nèi)平均流速為11.3 m/s，彎管為9.267 m/s，對應的Re分別為1.471×105和1.206×105，為湍流流動。

壓力與速度耦合采用SIMPLE算法，動量、湍動能k以及耗散率ω方程均采用二階迎風格式。時間步長設為2.020 47×10-5，即一個周期分為2 048個時間步。迭代斂殘差值設置為10-6。經(jīng)過大約8個周期，計算趨于穩(wěn)定后進行數(shù)據(jù)采集及保存，同時打開FLUENT瞬態(tài)計算中的Statistics選項以得到各時間步的速度、壓力的時均值和均方根值，共采集了16個周期的數(shù)據(jù)。

基于彎管流動更為復雜這一事實，對網(wǎng)格無關性、時間步長無關性以及數(shù)值模擬的準確性進行驗證。分別對單元數(shù)151.7萬、267.7萬、419.2萬三套網(wǎng)格進行計算，計算穩(wěn)定后進行時間步的平均，得到時均流場。用網(wǎng)格無關性驗證說明的彎頭處平均壓力場，如圖3所示。由圖3可知，三套網(wǎng)格的計算結果基本相同，267.7萬、419.2萬網(wǎng)格的結果非常接近，說明267.7萬網(wǎng)格的流場可視為網(wǎng)格無關。同時對267.7萬網(wǎng)格，進行時間步長為1.01×10-5s的計算，比較圖3(c)和圖3(d)可知流場幾乎無區(qū)別，表明時間步長已經(jīng)足夠小。

圖3 網(wǎng)格無關性驗證彎頭處平均壓力場

3 結果分析

3.1 平均流場

沿管道中心的平均靜壓分布，如圖4所示。由于兩個管道入口處總壓相等(都為131 Pa)，而直管內(nèi)的流動速度大，所以靜壓較低。對于直管，中心線上的靜壓分布與周圍壁面分布基本一致(橫截面上壓力分布均勻)，所以圖中僅顯示中心線處的數(shù)值，可以看出平均靜壓基本上呈線性下降。對于帶90°彎頭的彎管，顯示了中心線處及以上側、內(nèi)側、外側邊界層中的靜壓分布，可以看出上游直管靜壓基本上呈線性下降直到靠近彎頭進口截面。在彎頭進口附近，管壁內(nèi)側面壓力有所下降，外側面則上升，中心線和上側面基本一致。在彎頭下游的靜壓先下降然后恢復，且不同部位的壓力特性區(qū)別較大(橫截面上壓力分布不均勻)，再往下游管道橫截面上壓力分布又恢復均勻，表現(xiàn)出與直管相同的分布特性。

圖4 沿管道流向的靜壓分布

為研究流動沿管道的發(fā)展情況，在管道設置了若干橫截面，截面位置參見圖1。這些橫截面上的時均壓力分布，數(shù)據(jù)取自流向中間剖面，如圖5所示。由圖5可知，直管橫截面上的平均靜壓基本為常數(shù)；對于彎管，彎頭附近橫截面上壓力分布不均勻，但遠離彎頭區(qū)域的橫截面上壓力逐漸恢復均勻。

圖5 管道截面上的平均壓力分布

圖5中對應橫截面上的時均速分布，如圖6所示。由圖6可知，直管平均速度場為典型的管內(nèi)湍流分布形狀；對于彎管，90°彎頭對下游直管內(nèi)的速度分布影響非常明顯。

圖6 管道橫截面上的平均速度分布

彎頭處的流向及縱向中截面上的時均流場，如圖7所示，截面位置參見圖1(d)。從平均靜壓圖上可以看出，在彎頭上游直管壓力分布均勻，基本為平面波，在彎頭附近開始復雜化。結合壓力和速度場，可以看到在彎頭下游管道內(nèi)側部分存在明顯的低壓、低速的尾跡區(qū)，該尾跡區(qū)的存在導致外側流體流速增加，從而形成一定的射流現(xiàn)象，之后沿管道下游流動逐漸恢復均勻。

圖7 90°彎頭處平均流場

90°彎頭中間橫截面上的某一時刻速度投影矢量圖，如圖8所示。由圖8可以明顯看出二次流現(xiàn)象：彎頭內(nèi)側流體(圖中左邊)因轉彎離心力作用向外側運動，而在壁面附近由于速度低離心效果不明顯，流體在壓力差驅(qū)動下從外側擠向內(nèi)側，兩者效果疊加最終在橫截面上形成一對旋渦流。

圖8 90°彎頭中間橫截面的二次流

3.2 脈動流場

以脈動時域信號的標準差作為脈動強度的評價標準，沿管道流動方向的壓力脈動強度分布，如圖9和圖10所示，參考點位置參見圖1。

對于直管，其壁面邊界層的壓力脈動與管道中心的壓力脈動幾乎一樣，所以圖中省去直管壁面附近壓力脈動情況。由圖9可知，直管內(nèi)不可壓縮流場水動力壓力脈動幅度基本呈線性衰減。在出口處脈動強度為0，因為那里設置了靜壓為0的壓力邊界條件。動提高了局部壓力脈動強度。彎頭后直管壓力脈動呈線性下降。由圖10可知，盡管彎頭處平均流場紊亂，但管道中心與上側、外側面附近壓力脈動強度幾乎一致，且沿管道方向線性衰減。由于彎頭內(nèi)側曲率半徑小，其測點上的壓力脈動強度分布稍有不同，脈動強度線性衰減斜率小，這歸因于彎頭內(nèi)側出口處流場更為紊亂。

圖9 沿管道壓力脈動強度分布

圖10 90°彎頭處壓力脈動強度分布

對于彎管內(nèi)流動，由于90°彎頭的影響，使上半截直管內(nèi)的壓力脈動降低速度比直管慢。彎頭處流場擾

沿管道流動方向的脈動流場，如圖11和圖12所示。可以看出對于直管，壓力脈動在橫截面上基本一致，呈一維形式，而速度脈動在橫截面上并不均勻。對于彎管，脈動流場也有類似結論：盡管有90°彎頭存在，其內(nèi)部的壓力脈動仍呈一維形式；速度脈動相當紊亂，并且彎頭后直管內(nèi)速度脈動增加。

圖11 直管內(nèi)流向中間截面上的脈動流場

圖12 彎管內(nèi)的脈動流場

管內(nèi)中心線上3個測點的壓力脈動時間信號，如圖13所示。其距離指離管道入口的直管距離，彎頭拐角處指彎頭中間點(見圖 1中圓圈內(nèi)的點)?？梢钥闯鲋惫軆?nèi)各處的壓力脈動幾乎為正弦波，沿管道下游脈動幅度逐漸下降。各點處脈動波形無相位差，這是因為不可壓縮流動壓力脈動(偽聲)的微分方程為Poisson方程，聲速為無窮大，所以聲場內(nèi)各點的脈動具有同步性。對于彎管，管內(nèi)壓力脈動大體上仍為正弦波形狀，但由于彎頭對流動影響，波形中夾帶其他頻率分量。

圖13 管道內(nèi)中心線上監(jiān)測點壓力脈動

上述直管和彎管中心處3個測點的壓力脈動頻譜，如圖14所示。由圖14可知，580 Hz基頻最為突出，其低次諧波離散分量也較顯著，高于9次諧波的離散分量已經(jīng)不明顯。彎管內(nèi)流動的寬頻壓力脈動分量較直管大，這歸因于彎頭對流動的影響。整體上，沿管道壓力脈動強度逐漸下降。壁面附近邊界層內(nèi)壓力頻譜與管道中心處基本一致，說明水動力壓力脈動形態(tài)沿管道呈準一維方式，盡管速度場分布呈復雜三維流動。

圖14 管道中心線上監(jiān)測點壓力脈動頻譜

直管和彎管彎頭中間截面處測點的速度脈動頻譜，如圖15所示。由圖15可知， 580 Hz基頻及其低次諧波離散分量較明顯，這是因為流體黏性壁面附近的速度脈動略低于管道中心處的速度脈動。

圖15 管道內(nèi)監(jiān)測點速度脈動頻譜

4 結 論

本文分析了管內(nèi)脈動水動壓力(偽聲)的橢圓型微分方程，并與聲波方程進行了比較?；谟嬎懔黧w動力學，利用SSTk-ω湍流模型，數(shù)值研究管道入口處施加隨時間正弦變化壓力邊界條件下管內(nèi)水動壓力脈動特性，得到的結論如下：

(1) 直管道內(nèi)水動壓力波基本上呈一維形式，其橫截面上平均壓力及脈動壓力分布均勻。沿管道下游，壓力脈動幅度線性衰減。平均速度分布為典型的管道湍流速度分布。管道內(nèi)速度脈動十分紊亂，在管道橫截面上分布不均勻，這與壓力脈動的特性不同。

(2) 彎管內(nèi)流動在90°彎頭后面有明顯的尾跡區(qū)，而在橫截面上存在二次流。彎頭對彎內(nèi)流動的影響導致在彎頭附近流動呈三維復雜化，局部流動擾動明顯增強，使壓力脈動和速度脈動幅度都有所提高。盡管彎頭處流動呈復雜的三維流動，但壓力脈動的傳播整體上仍為一維形式。

(3) 壓力脈動頻譜中基頻及低次諧波明顯，9次諧波以上的離散分量已不明顯。90°彎頭使其附近流場的壓力脈動和速度脈動寬頻分量有所增強。由于壓力脈動微分方程為Poisson方程，擾動傳播速度為無窮大，所以管道各點的壓力脈動信號無相位差。