搖擺條件下水平圓管內(nèi)湍流壓力脈動(dòng)特性分析

朱仕斌，艾華寧

（中廣核研究院有限公司，深圳 518026）

核能作為一種高效、清潔能源，在使用安全性、穩(wěn)定性以及對(duì)環(huán)境的保護(hù)性上具有明顯優(yōu)勢(shì)以及極大的經(jīng)濟(jì)效應(yīng)［1］。而隨著我國(guó)“海洋強(qiáng)國(guó)”戰(zhàn)略的提出，海上小型反應(yīng)堆對(duì)島嶼的供電、供熱以及海洋油氣開(kāi)發(fā)等具有獨(dú)特的優(yōu)勢(shì)。其作為一種靈活、安全、經(jīng)濟(jì)的核電新形式，能夠彌補(bǔ)海上能源短缺的問(wèn)題。因此，發(fā)展海上小堆具有重大意義。

反應(yīng)堆冷卻劑流動(dòng)會(huì)誘發(fā)堆內(nèi)構(gòu)件的振動(dòng)，即流致振動(dòng)。長(zhǎng)期流致振動(dòng)可能使堆內(nèi)構(gòu)件結(jié)構(gòu)產(chǎn)生疲勞損傷或連接件發(fā)生松動(dòng)、磨損等問(wèn)題，可能造成事故。因此，堆內(nèi)構(gòu)件流致振動(dòng)理論分析和計(jì)算是核電設(shè)備安全分析不可缺少的重要內(nèi)容［2］。

工程經(jīng)驗(yàn)和大量研究表明，反應(yīng)堆堆內(nèi)構(gòu)件的流致振動(dòng)主要由隨機(jī)湍流激勵(lì)引起［3］。海上小堆在海洋上進(jìn)行作業(yè)時(shí)，在海浪、海風(fēng)等的作用下將引起平臺(tái)及核電站內(nèi)的管線產(chǎn)生起伏、傾斜、搖擺等不穩(wěn)定的運(yùn)動(dòng)狀態(tài)，從而改變湍流運(yùn)動(dòng)特性與空間分布，導(dǎo)致壓力脈動(dòng)的頻率、幅值發(fā)生變化，進(jìn)而影響結(jié)構(gòu)所受載荷的大小與分布，因此我們需要對(duì)海上環(huán)境因素重點(diǎn)考慮。目前，國(guó)內(nèi)外針對(duì)海洋條件下的流動(dòng)與傳熱研究主要集中探討流速、壓降、傳熱特性等宏觀量的變化規(guī)律［4-9］，對(duì)海洋條件引起的壓力脈動(dòng)的研究較少。

壓力脈動(dòng)主要由大尺度旋渦及小尺度的湍流產(chǎn)生，在一回路流道中，壓力脈動(dòng)主要由湍流引起。本文依據(jù)文獻(xiàn)［4］實(shí)驗(yàn)?zāi)Ｐ?，開(kāi)展數(shù)值模擬，選用該模型的原因在于：（1）有實(shí)驗(yàn)數(shù)據(jù)，便于數(shù)學(xué)模型的驗(yàn)證；（2）模型簡(jiǎn)單，壓力脈動(dòng)量均為湍流引起，從而剔除了大尺度旋渦對(duì)壓力脈動(dòng)的影響，為完善海洋條件下堆內(nèi)構(gòu)件流致振動(dòng)分析理論提供參考。

1 數(shù)理模型

1.1 數(shù)學(xué)模型

雷諾平均N-S方程方法（Reynolds-averaged Navier-Stokes equations，RANS），大渦數(shù)值模擬方法（Large eddy simulation，LES）及湍流直接數(shù)值模擬方法（Direct Numerical Simulation，DNS）是目前研究湍流和工程應(yīng)用的3種主要數(shù)值模擬方法。RANS無(wú)法捕捉脈動(dòng)量，而DNS的計(jì)算成本太高，因此本文采用LES。LES的基本思想［10］：在湍流數(shù)值模擬中只計(jì)算大尺度的脈動(dòng)，將小尺度脈動(dòng)對(duì)大尺度運(yùn)動(dòng)的作用建立模型。由于放棄了直接計(jì)算小尺度的脈動(dòng)，數(shù)值模擬的時(shí)間和空間步長(zhǎng)可以放大，從而緩解對(duì)計(jì)算資源的苛刻要求。

1.1.1 控制方程

守恒形式的納維-斯托克斯方程（Navier-Stokes equation，N-S方程）如下［11］：

式中，τ —粘性應(yīng)力張量；

ρ—密度（kg/m3）；

u —速度矢量（m/s）；

t—時(shí)間（s）；

P—壓力（Pa）；

g —重力加速度。

式（2）是通過(guò)分子粘性τ將能量耗散。只有當(dāng)渦足夠?。ㄐ〉終olmogorov尺度時(shí)）時(shí)，分子粘性才算足夠大，能將渦耗散掉。

LES中，通過(guò)增加一個(gè)額外的應(yīng)力項(xiàng)（亞格子應(yīng)力）來(lái)增大耗散，且能夠達(dá)到恰好耗散掉大于網(wǎng)格的渦的要求，此時(shí)N-S方程變?yōu)椋?/p>

式中，τsgs—亞格子應(yīng)力張量。

在亞格子模型中，普遍采用渦粘模型進(jìn)行計(jì)算：

式中，vsgs—亞格子粘性；

ksgs—亞格子湍動(dòng)能（J）；

—變形率張量；

δij—狄克拉函數(shù)。

采用該模型后，還存在一個(gè)未知量（亞格子粘性，vsgs），各亞格子模型的區(qū)別在于對(duì)亞格子粘性的計(jì)算方法不同，主要有Smagorinsky-Lilly子模型，WALE子模型，WMLES子模型。

本文中采用WALE子模型，該模型是在原始的Smagorinsky-Lilly模型基礎(chǔ)上，通過(guò)修改速度尺度來(lái)解決粘性子層中亞格子應(yīng)力不為零的情況。

1.1.2 慣性力模型

流體與實(shí)驗(yàn)回路壁面的相對(duì)速度是影響流動(dòng)與傳熱特性的決定因素，而不是流體相對(duì)于地面的絕對(duì)速度［12］，因此我們需要確定非慣性系下的動(dòng)量方程，慣性坐標(biāo)系（O0X0Y0Z0）與非慣性坐標(biāo)系（oxyz）如圖1所示。

圖1 慣性坐標(biāo)系與非慣性坐標(biāo)系Fig.1 Inertial & non-inertail coordinates

通過(guò)建立慣性坐標(biāo)系與非慣性坐標(biāo)系的關(guān)系，并將慣性坐標(biāo)系下的動(dòng)量方程代入，獲得非慣性系下的動(dòng)量方程：

式中，ao—平動(dòng)加速度（m/s）；

β —轉(zhuǎn)動(dòng)角加速度（rad/s2）；

ω —轉(zhuǎn)動(dòng)加速度（rad/s）；

r —控制體在非慣性坐標(biāo)系的位置（m）；

ur—相對(duì)速度（m/s）。

對(duì)比式（4）、式（7），并結(jié)合等效原理［13］可知，非慣性坐標(biāo)系下的慣性力F為：

通過(guò)自定義函數(shù)（User Defined Function，UDF）將重力及慣性力添加至Fluent的動(dòng)量源項(xiàng)中。

1.2 物理模型

海洋條件下的實(shí)際運(yùn)動(dòng)復(fù)雜，為六自由度耦合作用的結(jié)果。為清晰了解海洋條件下的運(yùn)動(dòng)對(duì)湍流壓力脈動(dòng)的影響與關(guān)系，我們通常研究單個(gè)自由度的影響。本文以文獻(xiàn)［2］的實(shí)驗(yàn)?zāi)Ｐ妥鳛閹缀文Ｐ烷_(kāi)展搖擺運(yùn)動(dòng)對(duì)湍流壓力脈動(dòng)的數(shù)值研究，其搖擺臺(tái)以及搖擺方式見(jiàn)圖2，其中ab段為實(shí)驗(yàn)段，bc-cd段為出口段，管路繞x軸進(jìn)行搖擺運(yùn)動(dòng)。管長(zhǎng)為4.5 m，管徑為0.0345 m。為消除入口段效應(yīng)，我們?nèi)?shí)驗(yàn)段入口1.6 m后的數(shù)據(jù)進(jìn)行測(cè)量，得到管內(nèi)MN段（2 m）壓力隨時(shí)間變化的實(shí)驗(yàn)數(shù)據(jù)。

圖2 搖擺臺(tái)及搖擺方式示意圖Fig.2 Rolling bench and rolling mode

根據(jù)實(shí)驗(yàn)段模型，建立如圖3所示的數(shù)值計(jì)算模型。在截面D壁面處均布12個(gè)監(jiān)測(cè)點(diǎn)，用于監(jiān)測(cè)壓力隨時(shí)間變化。

圖3 數(shù)值計(jì)算檢測(cè)面布置圖Fig.3 The scheme of measuring sections for numerical calculation

2 數(shù)值計(jì)算

2.1 計(jì)算工況

為獲得搖擺周期及雷諾數(shù)對(duì)湍流壓力脈動(dòng)的影響，我們計(jì)算了如表1所示的7個(gè)工況，其中前三個(gè)工況無(wú)搖擺，后四個(gè)工況有搖擺。

表1 水平管道壓力脈動(dòng)數(shù)值計(jì)算工況表Table 1 The case table of horizontal pipeline pulsation numerical simulation

2.2 計(jì)算網(wǎng)格

水平圓管模型簡(jiǎn)單，因此本文采用ICEM軟件生成結(jié)構(gòu)化網(wǎng)格，如圖4所示。為進(jìn)行網(wǎng)格敏感性測(cè)試，以工況1為例，本文計(jì)算無(wú)不同網(wǎng)格總數(shù)下的壓力脈動(dòng)均方根（Root Mean Square，RMS）值，見(jiàn)表2?？梢钥吹剑S著網(wǎng)格總數(shù)的增加，RMS值逐漸增大，但增長(zhǎng)趨勢(shì)逐漸減小，在綜合考慮計(jì)算資源與計(jì)算效率的前提下，我們最終選定604萬(wàn)的網(wǎng)格進(jìn)行后續(xù)分析。該網(wǎng)格的MN段（關(guān)注區(qū)）的網(wǎng)格尺寸為1 mm，第一層網(wǎng)格高度設(shè)置為9×10-3mm，保證y plus約為1。

表2 網(wǎng)格無(wú)關(guān)性測(cè)試對(duì)比表Table 2 The comparison table of grid independence test

圖4 結(jié)構(gòu)網(wǎng)格示意圖Fig.4 The scheme of structured grid

2.3 邊界條件

管道入口設(shè)置為速度入口，管道出口設(shè)置為壓力出口，壁面為無(wú)滑移邊界，操作壓力設(shè)置為一個(gè)大氣壓。

2.4 求解器設(shè)置

考慮到流場(chǎng)介質(zhì)為不可壓液體，因此本文采用壓力基求解器（Pressure-based solver）。速度壓力耦合算法為SIMPLE算法，空間與時(shí)間的離散格式均為2階。

為快速獲取瞬態(tài)計(jì)算結(jié)果，本文開(kāi)始時(shí)采用RANS模型進(jìn)行穩(wěn)態(tài)計(jì)算，得到穩(wěn)態(tài)場(chǎng)后生成脈動(dòng)速度場(chǎng)，以此作為L(zhǎng)ES的初始值。LES的時(shí)間步長(zhǎng)應(yīng)滿(mǎn)足每一時(shí)間步流體所流過(guò)的距離小于網(wǎng)格尺寸，以此得到時(shí)間步長(zhǎng)為0.5 ms。待計(jì)算消除初始效應(yīng)后，我們將統(tǒng)計(jì)量歸零，并開(kāi)啟采樣，以獲取壓力、速度等參數(shù)的統(tǒng)計(jì)量。

3 計(jì)算結(jié)果分析

3.1 數(shù)學(xué)模型的驗(yàn)證

圖5給出了工況4條件下的MN段壓降隨時(shí)間的變化曲線，其中，參考曲線［8］為數(shù)值計(jì)算所得，并且采用與本文相同的湍流模型?？梢钥吹剑?jì)算曲線、試驗(yàn)曲線［4］與參考曲線相近，因此加載至Fluent的慣性力模型是可信的。數(shù)值計(jì)算所得結(jié)果均表現(xiàn)為后半周期更接近實(shí)驗(yàn)值，并且本文計(jì)算結(jié)果優(yōu)于參考文獻(xiàn)所得結(jié)果。

圖5 MN段壓降隨時(shí)間變化曲線Fig.5 The pressure drop curve of MN segment with time

3.2 速度場(chǎng)定性分析

圖6為工況1、工況4和工況5在不同搖擺相位角下y-z截面的速度分布云圖，為便于觀察，在后處理時(shí)將軸向長(zhǎng)度與徑向長(zhǎng)度的比例按20∶1進(jìn)行了縮減顯示。工況1無(wú)搖擺運(yùn)動(dòng)，因此本文取1、2、3、4秒時(shí)刻的結(jié)果進(jìn)行對(duì)比。由于慣性力及重力的影響，我們可以在工況4和工況5的y-z截面上觀察到速度場(chǎng)隨相位的上下擺動(dòng)，而工況1則表現(xiàn)出各個(gè)時(shí)刻下速度場(chǎng)相近。

圖6 y-z截面速度分布云圖Fig.6 The contour of velocity distribution of y-z section

3.3 時(shí)域分析

從圖6可以看出流體從入口段的均勻來(lái)流向湍流轉(zhuǎn)捩的過(guò)程。在D截面處，湍流已充分發(fā)展。圖7給出了工況5在D截面的12個(gè)監(jiān)測(cè)點(diǎn)的時(shí)程曲線，可以看到，每個(gè)節(jié)點(diǎn)的壓力值在各個(gè)時(shí)刻均相近，因此本文以D截面的1節(jié)點(diǎn)為例，給出不同工況的對(duì)比曲線，如圖8所示?？梢钥吹?，湍流的壓力脈動(dòng)量相對(duì)于平均值而言是小量。從圖8（a）可知，節(jié)點(diǎn)壓力的峰值隨搖擺頻率的增加而增加。圖8（a）、圖8（c）表明，監(jiān)測(cè)點(diǎn)處壓力的響應(yīng)與搖擺周期是同步變化的。

圖7 監(jiān)測(cè)點(diǎn)壓力隨時(shí)間變化曲線（工況5）Fig.7 Monitoring point pressure curve with time（case5）

圖8 監(jiān)測(cè)點(diǎn)壓力隨時(shí)間變化曲線Fig.8 Monitoring point pressure curve with time

圖8 （續(xù)）

受搖擺運(yùn)動(dòng)的影響，工況4—工況7存在一個(gè)與搖擺周期相同的壓力波動(dòng)，該波動(dòng)是對(duì)平均場(chǎng)的貢獻(xiàn)，在計(jì)算RMS值時(shí)需要將平均場(chǎng)的貢獻(xiàn)剔除。本文采用高通濾波的方法，將截止頻率設(shè)置為0.5 Hz（略大于搖擺頻率），對(duì)所有工況的所有監(jiān)測(cè)點(diǎn)進(jìn)行濾波，并計(jì)算濾波后的壓力脈動(dòng)RMS值，RMS計(jì)算公式如下：

式中，pi—預(yù)測(cè)值，為采樣結(jié)果；

ai—真實(shí)值，被高通過(guò)濾器濾掉的波動(dòng)值；

n—樣本總數(shù)。

7組工況下12個(gè)監(jiān)測(cè)點(diǎn)的RMS值的對(duì)比曲線如圖9所示。圖9（a）為同一雷諾數(shù)下不同搖擺頻率的RMS值對(duì)比圖，可以看到，F(xiàn)luent采樣所得結(jié)果與濾波后的結(jié)果存在偏差，這是對(duì)數(shù)據(jù)進(jìn)行了濾波所致，該雷諾數(shù)下所得RMS均為小量，搖擺頻率對(duì)壓力脈動(dòng)的影響極小。圖9（b）為無(wú)搖擺條件下雷諾數(shù)倍增后的RMS值對(duì)比圖，可以看到，隨著雷諾數(shù)的倍增，工況2的RMS值比工況1增加了約45 Pa，而工況3比工況2增加了約10 Pa，這說(shuō)明隨著雷諾數(shù)的增加，其對(duì)壓力脈動(dòng)的影響逐漸減弱，從圖9（c）也可以得出相同的結(jié)論，這與參考文獻(xiàn)［14］所述一致。

圖9 監(jiān)測(cè)點(diǎn)的RMS值對(duì)比圖Fig.9 Comparison of RMS values at monitoring point

圖10為搖擺和無(wú)搖擺條件下，RMS值隨雷諾數(shù)的變化曲線，可以看到，搖擺條件下的斜率比無(wú)搖擺條件斜率大，即搖擺條件的引入增大了雷諾數(shù)對(duì)RMS值的影響。

圖10 RMS值隨雷諾數(shù)變化曲線Fig.10 The value of RMS curve with Reynolds number

3.4 頻域分析

觀察圖7可以發(fā)現(xiàn)，初始效應(yīng)在1 s后消失，因此本文取1 s后的樣本進(jìn)行頻域分析，為控制變量，均取1~5 s共4 s的采樣數(shù)據(jù)進(jìn)行高通濾波，后對(duì)該數(shù)據(jù)做快速傅立葉變換，最終可得到各工況下各監(jiān)測(cè)點(diǎn)的功率譜密度（Power Spectral Density，PSD）曲線。從時(shí)域分析可知，在同一工況下相同截面上的各個(gè)監(jiān)測(cè)點(diǎn)的時(shí)程曲線及RMS值相近，其PSD曲線也應(yīng)相近，因此本文取1節(jié)點(diǎn)進(jìn)行對(duì)比，如圖11所示?？梢钥吹?，壓力脈動(dòng)表現(xiàn)為寬頻脈動(dòng)。對(duì)比圖9和圖11可知，RMS值大的工況，雖然局部PSD值會(huì)更小，但總趨勢(shì)呈現(xiàn)出全頻率段的幅值更大，這是因?yàn)殡S著搖擺頻率/雷諾數(shù)的增加，系統(tǒng)的平均場(chǎng)獲得了更多的能量，而該能量又通過(guò)能量串級(jí)的方式將低頻能量傳遞到高頻能量處，直至通過(guò)粘性耗散將動(dòng)能轉(zhuǎn)化為熱能，而平均場(chǎng)能量越大，傳遞給高頻的能量就越多，即表現(xiàn)為PSD曲線的幅值越高。

4 結(jié)論

本文采用ANSYS Fluent求解器，以水平圓管為例，通過(guò)UDF將慣性力模型添加至動(dòng)量方程的源項(xiàng)中，開(kāi)展LES數(shù)值模擬，并對(duì)監(jiān)測(cè)點(diǎn)處的數(shù)據(jù)進(jìn)行時(shí)域和頻域分析，得出以下結(jié)論：

（1）相對(duì)于平均壓力而言，湍流壓力脈動(dòng)值占比極小，表現(xiàn)為監(jiān)測(cè)點(diǎn)處壓力時(shí)程曲線上出現(xiàn)小的“毛刺”；

（2）觀察有搖擺條件的工況發(fā)現(xiàn)，速度場(chǎng)與監(jiān)測(cè)點(diǎn)處的壓力曲線的波動(dòng)均呈現(xiàn)與搖擺同步的正弦擺動(dòng)；

（3）對(duì)比搖擺頻率的變化發(fā)現(xiàn)，RMS值隨著搖擺頻率的增加而增加，而無(wú)論有無(wú)搖擺，RMS值都隨雷諾數(shù)的增加而增加；

（4）對(duì)比搖擺及不搖擺工況，隨著雷諾數(shù)的增加，其對(duì)壓力脈動(dòng)的影響將逐漸減弱，表現(xiàn)為斜率逐漸減小，但搖擺條件增大了雷諾數(shù)對(duì)RMS值的影響。