李庚旺 張 騏 陳 爽 歐陽斌 盧 濤*

(1.北京化工大學(xué) 機(jī)電工程學(xué)院， 北京 100029； 2.中國核動力研究設(shè)計(jì)院， 成都 610200)

引 言

T型管是石油、化工和核電管道系統(tǒng)中使用頻率較高的結(jié)構(gòu)，在主支管流體動量比很大的工況下，主管流體會侵入到支管，即湍流穿透現(xiàn)象。當(dāng)湍流穿透區(qū)域與冷熱流體混合區(qū)域重疊時(shí)，會引發(fā)流場內(nèi)強(qiáng)烈的溫度波動，管壁內(nèi)部隨之出現(xiàn)溫度波動，使管道的熱應(yīng)力也產(chǎn)生波動[1]，進(jìn)而可能誘發(fā)高周熱疲勞，使管道失效。

Kickhofel等[2]對帶泄漏流的T型管內(nèi)的湍流穿透現(xiàn)象進(jìn)行了實(shí)驗(yàn)研究，得出湍流穿透速度波動的劇烈程度與主支管動量比之間的關(guān)系。Hirota等[3]對T型(三通)管內(nèi)冷熱流體混合過程進(jìn)行了可視化實(shí)驗(yàn)，結(jié)果表明混合界面的垂直方向上會出現(xiàn)震蕩以及旋渦疊加。盧濤等[4]利用大渦模擬(LES)方法對T型管道內(nèi)冷熱流體混合過程進(jìn)行了數(shù)值計(jì)算，得出浮升力與混合過程波動情況的關(guān)系。Paffumi等[5]對材料為316L不銹鋼的管道進(jìn)行疲勞試驗(yàn)，研究了溫差、軸向力以及壁厚對熱疲勞的影響。Lee等[6]運(yùn)用大渦模擬獲得了近壁面流體的溫度波動，研究發(fā)現(xiàn)冷熱流體溫差和由湍流混合引起的強(qiáng)化對流換熱是T型管道疲勞失效的主要因素。張?jiān)絒7]使用瞬態(tài)界面耦合方法對管道進(jìn)行應(yīng)力分析，計(jì)算了管道疲勞壽命。綜上所述，目前對于T型管內(nèi)流動與熱現(xiàn)象的研究多集中于熱分層。此外，關(guān)于流固耦合的研究多為穩(wěn)態(tài)耦合或瞬態(tài)面耦合，體節(jié)點(diǎn)的數(shù)據(jù)需要在耦合計(jì)算過程重新插值，而流場分析軟件一般使用有限體積法，固體分析軟件一般使用有限單元法，重新插值計(jì)算得到的數(shù)據(jù)存在誤差。

本文使用Fluent軟件對T型管湍流穿透區(qū)域的流體速度與溫度場進(jìn)行分析，得出管內(nèi)各物理量的時(shí)空演變規(guī)律；使用用戶自定義函數(shù)(UDF)、批處理文件(Journal)以及CFD- Post軟件中的宏計(jì)算器(Macrocalculator)進(jìn)行數(shù)據(jù)提取，通過ANSYS Workbench平臺動態(tài)加載到管道內(nèi)部，確保每一時(shí)間步內(nèi)網(wǎng)格節(jié)點(diǎn)溫度完全耦合，得到固體應(yīng)力波動信息并對危險(xiǎn)點(diǎn)的疲勞壽命進(jìn)行評估。

1 計(jì)算模型

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

T型管內(nèi)湍流穿透現(xiàn)象伴隨著流場溫度與速度的劇烈波動，因此需要進(jìn)行瞬態(tài)分析。與使用時(shí)均方法的k-ε模型相比，大渦模擬方法可以較好地捕捉管內(nèi)流體溫度與速度的瞬時(shí)波動細(xì)節(jié)[8-11]，因此本文使用大渦模擬方法進(jìn)行流場數(shù)值計(jì)算。對于不同的大渦模擬湍流模型，其本質(zhì)區(qū)別在于亞格子(SGS)模型不同，亦即如何描述亞格子雷諾應(yīng)力。

亞格子雷諾應(yīng)力可以定義為

(1)

SGS模型通常使用渦- 黏模型計(jì)算，即

(2)

本文采用Smagorinsky-lily亞格子湍流黏度

(3)

其中Ls表示亞格子尺度的混合長度，即

Ls=min (kd,CSV1/3)

(4)

式中，k為馮卡門常數(shù)，k=0.41；d為到最近的壁面的距離；CS為Smagorinsky常數(shù)，取值0.1；V為計(jì)算單元體積。

T型管冷熱流體摻混過程中，重力與浮升力共同發(fā)生作用。在大渦模擬模型的動量方程中，針對浮升力項(xiàng)的求解需要對流體進(jìn)行變物性設(shè)置。根據(jù)水的物性表[7]對其密度隨溫度的變化進(jìn)行擬合，得到三階多項(xiàng)式

ρ=750.596+2.145 03×T-0.005 12×T2+0.000 002 312×T3

(5)

式中ρ為流體密度，kg/m3；T為流體溫度，K。

使用有限元方法對耦合熱應(yīng)力與應(yīng)變進(jìn)行求解，將流場計(jì)算的溫度場作為載荷對固體域進(jìn)行計(jì)算，具體理論模型為[12]

σ=D(ε-ε0)

(6)

式中，D為彈性矩陣，ε0為溫度變化引發(fā)的變形量

ε0=αΔT[1 1 1 0 0 0]T

(7)

式中，α為材料的線膨脹系數(shù)，ΔT為溫度變化。

ε可由式(8)求得

ε=Bδe

(8)

其中B為應(yīng)變矩陣，節(jié)點(diǎn)位移δe可由式(9)求出

Kδe=QT

(9)

式中，K為剛度矩陣,QT為受到的熱載荷。

1.2 物理模型

主管長530 mm，內(nèi)徑120 mm，壁厚3.5 mm;支管長1 100 mm，內(nèi)徑60 mm，壁厚2.5 mm。對流體域與固體域同時(shí)劃分全結(jié)構(gòu)化網(wǎng)格。幾何模型與網(wǎng)格模型如圖1所示，網(wǎng)格數(shù)量為962 339，第一層網(wǎng)格高度為0.01 mm，邊界層網(wǎng)格增長率為1.1，層數(shù)為15，y+值為4。網(wǎng)格無關(guān)性驗(yàn)證如圖2所示，分別對3種不同網(wǎng)格數(shù)量的模型(mesh1 1 825 205；mesh2 962 339；mesh3 405 416)進(jìn)行計(jì)算，得到管內(nèi)H=2.5截面，90°位置點(diǎn)溫度隨時(shí)間的變化曲線。從圖2可以看出，mesh1、mesh2計(jì)算結(jié)果相差較小，為兼顧計(jì)算精度與計(jì)算效率，下文采用mesh2作為流固耦合計(jì)算網(wǎng)格。

為便于定量探究管內(nèi)流體溫度與速度隨時(shí)間的變化情況，分別定義速度及溫度的統(tǒng)計(jì)量。

定義軸向時(shí)均速度Vz,mean

(10)

定義軸向均方根速度Vz,rms

(11)

絕對速度Vi為

(12)

式中，Vx,i、Vy,i、Vz,i分別為i時(shí)刻x、y、z方向的速度大小，對于支管，z方向速度即為軸向速度，單位為m/s。

時(shí)均絕對速度Vmean

(13)

均方根絕對速度Vrms

(14)

(15)

(16)

(17)

定義支管軸向無量綱高度為

(18)

式中，N為采樣時(shí)間點(diǎn)總數(shù)，i為采樣時(shí)刻序號，Tcold為主管入口溫度，Thot為支管入口溫度，z為軸向高度坐標(biāo)，Dm為主管直徑，t為主管厚度，Db為支管直徑。

2 參數(shù)設(shè)置

2.1 流場計(jì)算

流場計(jì)算的具體邊界條件見表1。計(jì)算過程中，按主管入口邊界條件進(jìn)行初始化，設(shè)置操作壓力為15.5 MPa。首先使用k-ε模型進(jìn)行穩(wěn)態(tài)計(jì)算，當(dāng)流動達(dá)到穩(wěn)定時(shí)開啟瞬態(tài)LES模型，使用couple耦合求解器與二階迎風(fēng)格式求解。瞬態(tài)求解總時(shí)間為15 s，時(shí)間步長為0.005 s，數(shù)據(jù)監(jiān)測頻率為200 Hz。沿支管軸向方向共設(shè)置20個(gè)監(jiān)測截面，每個(gè)截面設(shè)置多個(gè)近壁面監(jiān)測點(diǎn)，其中0°到180°的區(qū)域(正半?yún)^(qū))正對來流方向，波動較劇烈，因此該區(qū)域沿周向每隔15°設(shè)置一個(gè)監(jiān)測點(diǎn)，背對來流方向的位置(負(fù)半?yún)^(qū))每隔30°設(shè)置一個(gè)監(jiān)測點(diǎn)，監(jiān)測近壁面1 mm處溫度與速度波動信息。具體設(shè)置情況見圖3。

表1 流場計(jì)算邊界條件Table 1 Boundary conditions in the CFD simulation

2.2 流固耦合計(jì)算

瞬態(tài)流固耦合計(jì)算分為單向和雙向，在實(shí)際運(yùn)行過程中，管道的變形相對較小，對流場的影響可以忽略不計(jì)，因此選擇單向耦合方法進(jìn)行計(jì)算。將流場計(jì)算結(jié)果如管道內(nèi)壁面壓力、固體域溫度等作為載荷動態(tài)加載到固體域，加載時(shí)間步長為0.025 s，耦合計(jì)算時(shí)間共5 s，導(dǎo)入5～10 s內(nèi)的流場數(shù)據(jù)。

根據(jù)實(shí)際情況選擇管道材料為316L不銹鋼，具體邊界條件如表2所示。

表2 流固耦合計(jì)算邊界條件Table 2 Boundary conditions of the CFD- FEM simulation

3 結(jié)果與分析

3.1 流場計(jì)算結(jié)果

3.1.1速度場分析

在湍流穿透發(fā)生的區(qū)域，流場的流動情況十分復(fù)雜，在各個(gè)方向上均存在強(qiáng)烈的速度波動，本文主要探究速度波動的強(qiáng)度沿管道軸向以及周向的分布，因此對絕對速度V的時(shí)均值及均方根值進(jìn)行分析。圖4為Vmean與Vrms沿軸向的分布情況。

分析圖4(a)，當(dāng)H<2時(shí)，對比同一截面不同位置的時(shí)均速度，90°處的值遠(yuǎn)高于其他位置的值，這是由于主管流體在流經(jīng)交匯處時(shí)，由于流動截面積突增，流動方向壓力變大，致使流動速度逐漸減小至0，由負(fù)半?yún)^(qū)侵入支管；當(dāng)H>2.5時(shí)，90°位置的時(shí)均速度迅速下降至與其他位置幾乎相同，原因是主管流體侵入支管后，遇到速度相反的支管流體，向上的動量被迅速抵消；從H=4位置開始，曲線斜率幾乎為0，說明此區(qū)域內(nèi)占主導(dǎo)地位的是支管的低速流體。

圖4(b)中曲線呈現(xiàn)的規(guī)律與圖4(a)相似，90°位置的波動在H<2時(shí)遠(yuǎn)高于同截面的其他位置，但在H=2.5附近，0°曲線對應(yīng)值略高于90°，說明速度波動峰值的周向位置并不完全固定，由于流體的湍流特性，偶爾會在其他位置出現(xiàn)；在H=4處，曲線值接近0，說明在當(dāng)前工況下，湍流穿透的深度為H=4；注意到在10

進(jìn)一步就速度波動峰值的周向位置不固定的現(xiàn)象進(jìn)行分析。圖5為不同高度截面上Vrms的周向分布情況。

對比圖5中的4幅圖，可以發(fā)現(xiàn)隨著截面軸向高度的增加，均方根速度峰值出現(xiàn)的位置在不斷變化。當(dāng)H分別為1.5、2.0、2.5時(shí)，峰值出現(xiàn)的位置為75°、75°、30°，其中前兩個(gè)截面峰值位置相同，但波動幅度較大區(qū)域存在明顯的偏轉(zhuǎn)；當(dāng)H=4.0時(shí)，曲線的周向分布相對均勻，峰值與谷值之間相差很小，可以認(rèn)為在該位置速度波動已趨于平緩，與圖4分析結(jié)果一致；此外，還注意到H=1.5～2.5截面曲線峰值的偏轉(zhuǎn)方向?yàn)轫槙r(shí)針，這與流體進(jìn)入支管后受到的初始擾動有關(guān)。

3.1.2溫度場分析

圖6為無量綱溫度統(tǒng)計(jì)量的軸向分布。由圖6(a)可以看出，在湍流穿透工況下，波動管內(nèi)近壁面流場溫度的波動強(qiáng)度隨軸向高度的增加先增大后減小，在H=4位置達(dá)到峰值，且在該高度附近位置波動強(qiáng)度變化幅度很大，原因是該截面以下為速度湍流穿透區(qū)域，主管低溫流體在該區(qū)域占主導(dǎo)地位，支管流體還未流到該區(qū)域就被冷卻；當(dāng)主管流體向上侵入至H=4高度時(shí)，其動量被支管流體完全抵消，此位置溫度升高的速度急劇增大；當(dāng)H>4時(shí)，流體溫度達(dá)到穩(wěn)定，處于熱穩(wěn)定區(qū)，波動強(qiáng)度減弱，最終降低為0。

由圖6(b)可以看出，流體平均溫度隨軸向高度增加而升高，曲線斜率先增大后減小，當(dāng)H=4時(shí)為最大值，該位置對應(yīng)流體溫度波動最強(qiáng)烈的區(qū)域。

圖7為不同截面上無量綱均方根溫度的周向分布情況?？梢园l(fā)現(xiàn)隨著軸向高度的增加，不同截面的均方根溫度峰值在正負(fù)半?yún)^(qū)交替出現(xiàn)，當(dāng)H為1.5、2.5時(shí)，峰值出現(xiàn)在正半?yún)^(qū)，當(dāng)H=4.0時(shí)，峰值出現(xiàn)在負(fù)半?yún)^(qū)，此截面處溫度波動達(dá)到最大；當(dāng)H=5.5時(shí)，均方根溫度的周向分布趨于均勻，峰值與谷值之間相差很小，說明此區(qū)域已處于熱穩(wěn)定區(qū)。

3.2 流固耦合計(jì)算結(jié)果

3.2.1應(yīng)力分析

首先對耦合計(jì)算的應(yīng)力結(jié)果進(jìn)行分析。圖8為T型管不同軸向高度截面的應(yīng)力云圖，可以看出，在H=0截面，最大應(yīng)力出現(xiàn)的區(qū)域分別在管道0°及90°位置的內(nèi)壁面處，這是由于該區(qū)域?yàn)橄嘭瀰^(qū)，存在應(yīng)力集中效應(yīng)；當(dāng)H>2時(shí)，應(yīng)力最大位置出現(xiàn)在270°位置，原因是支管受到熱載荷產(chǎn)生熱應(yīng)力與熱變形，從而沿管道軸向方向膨脹，由于主支管入口皆為固定約束，該膨脹導(dǎo)致主管出口產(chǎn)生Z方向的變形；對比不同截面高度的應(yīng)力分布，從H=0截面開始，在270°位置應(yīng)力先增大后減小，這是由冷熱流體在該區(qū)域摻混產(chǎn)生的熱應(yīng)力引起。

3.2.2疲勞分析

對危險(xiǎn)點(diǎn)處的應(yīng)力波動信息進(jìn)行處理，使用疲勞分析理論進(jìn)行評估。疲勞問題分為高周疲勞與低周疲勞，后者主要針對單次應(yīng)力循環(huán)較大且循環(huán)次數(shù)較低的情況，本文研究的湍流穿透工況管道單次循環(huán)應(yīng)力較小，實(shí)際循環(huán)次數(shù)遠(yuǎn)高于105，因此屬于高周疲勞，需要分析名義應(yīng)力，并使用S-N曲線進(jìn)行疲勞評估。

對于等效多軸應(yīng)力的波動信息，使用雨流計(jì)數(shù)法[13]統(tǒng)計(jì)其全循環(huán)次數(shù)，得到雨流矩陣N,再依據(jù)Goodman曲線進(jìn)行等效，得到等效對稱循環(huán)應(yīng)力矩陣NG。

Goodman理論的具體公式為

(19)

式中，Sa(-1)為修正后的等效對稱應(yīng)力，Sa為應(yīng)力幅，Sm為平均應(yīng)力，Su為材料的極限強(qiáng)度。本文使用的316L不銹鋼的強(qiáng)度極限取520 MPa。

根據(jù)應(yīng)力分析，選定H=0截面的最大應(yīng)力點(diǎn)作為分析對象，該點(diǎn)的等效多軸應(yīng)力波動情況如圖9所示。

對應(yīng)力波動信息進(jìn)行頻域分析，由圖9(b)可看出，功率譜密度曲線最大值對應(yīng)的頻率為0.59 Hz，其周期為1.69 s。由于在耦合計(jì)算的前5 s及其后續(xù)的計(jì)算時(shí)間內(nèi)流場與固體場各物理量的波動狀態(tài)較為穩(wěn)定，且5 s的時(shí)間包含應(yīng)力波動的主周期，因此以5 s內(nèi)的數(shù)據(jù)作為代表，為便于雨流計(jì)數(shù)法分析，將其延拓至25 s。

圖10為雨流計(jì)數(shù)法得到的應(yīng)力雨流矩陣，圖11為對雨流矩陣進(jìn)行修正處理后得到的Goodman修正應(yīng)力矩陣。

(20)

式中Nfitting為Goodman修正應(yīng)力對應(yīng)的循環(huán)次數(shù)。

根據(jù)Miner疲勞累計(jì)損傷準(zhǔn)則得到累計(jì)疲勞損傷Df為

(21)

由式(21)得Df=1.657 7×10-8，疲勞壽命L=1/Df=6.032 4×107個(gè)周期，每個(gè)周期25 s，最終得到該位置的壽命為47.82年。

4 結(jié)論

(1)使用大渦模擬湍流模型預(yù)測了T型管交匯處流體的流動與傳熱情況，發(fā)現(xiàn)在支管軸向無量綱高度H=4截面處溫度波動強(qiáng)度最大。

(2)構(gòu)建了熱- 流- 固耦合方法，獲得了管道熱應(yīng)力的空間分布情況與瞬態(tài)波動信息，確定了管道疲勞失效的危險(xiǎn)點(diǎn)為軸向無量綱高度H=0截面上的應(yīng)力集中位置。

(3)對危險(xiǎn)點(diǎn)的熱應(yīng)力波動信息進(jìn)行處理，通過疲勞分析理論得到了該結(jié)構(gòu)的疲勞壽命為47.82年。