陳柏宇 周照春 孫 燕 陳 爽 歐陽(yáng)斌 盧 濤*

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

引 言

在核電管路系統(tǒng)中，T型管結(jié)構(gòu)主支管流體的溫差會(huì)使得管道內(nèi)流體在摻混時(shí)出現(xiàn)劇烈的溫度波動(dòng)，引發(fā)管道產(chǎn)生熱應(yīng)力，導(dǎo)致疲勞破壞的發(fā)生。改變主支管流體的動(dòng)量比可以使得管道內(nèi)流體的摻混流型發(fā)生變化，進(jìn)而改變管道內(nèi)溫度的波動(dòng)區(qū)域，從而改善管道的熱疲勞情況。Hosseini等[1]通過(guò)實(shí)驗(yàn)探究了在T型管主管上游有無(wú)彎管的情況下，不同主支管流體動(dòng)量比MR對(duì)管道內(nèi)流體流動(dòng)的影響，并根據(jù)動(dòng)量比大小將流型分為壁面流(MR>4)、再附著流(1.35

對(duì)于流固耦合問(wèn)題及熱疲勞評(píng)估問(wèn)題，國(guó)內(nèi)外學(xué)者也進(jìn)行了深入的研究。胡麗娜等[4]為探究在啟、停堆工況下管道的熱應(yīng)力狀態(tài)，在耦合計(jì)算時(shí)對(duì)三通管道施加突變的溫度載荷，得到管道熱疲勞壽命與支管入射角度之間的關(guān)系式。Hannink等[5]通過(guò)依次導(dǎo)入不同時(shí)刻下的管道溫度場(chǎng)，計(jì)算出多個(gè)時(shí)刻下的管道熱應(yīng)力狀態(tài)，并對(duì)計(jì)算結(jié)果進(jìn)行了分析。Kumar等[6]研究了熱分層引起的熱應(yīng)力及其對(duì)管道疲勞壽命的影響。Guo等[7]對(duì)水平T型管進(jìn)行模擬試驗(yàn)和分析，設(shè)計(jì)了一個(gè)監(jiān)測(cè)管路壁溫的程序，能夠通過(guò)測(cè)量外壁溫度來(lái)監(jiān)測(cè)水平T型管冷熱流體混合時(shí)的壁溫波動(dòng)。Garrido等[8]通過(guò)生成隨機(jī)的管內(nèi)壁面溫度，使之與流場(chǎng)計(jì)算得到的隨機(jī)溫度載荷類(lèi)似，進(jìn)而探究在非穩(wěn)態(tài)載荷下T型管的熱應(yīng)力波動(dòng)。

在目前的研究中，評(píng)估不同動(dòng)量比對(duì)管道熱疲勞影響的工作較少，而T型管內(nèi)流體的攪混流型對(duì)管道的熱疲勞又有直接影響。本文通過(guò)LES模型對(duì)T型管內(nèi)流體的攪混狀態(tài)進(jìn)行瞬態(tài)模擬，在保證使用同一套計(jì)算網(wǎng)格模型的前提下，通過(guò)CFD-post軟件導(dǎo)出固體溫度載荷及內(nèi)壁面壓力載荷，基于workbench平臺(tái)進(jìn)行流固耦合計(jì)算，并評(píng)估管道的熱疲勞情況，為研究T型管內(nèi)動(dòng)量比對(duì)管道熱疲勞的影響提供了一定的數(shù)據(jù)支撐。

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

1.1 LES模型

本文將k-ε湍流模型計(jì)算得到的流場(chǎng)穩(wěn)態(tài)解作為初始場(chǎng)，利用LES模型進(jìn)行流場(chǎng)的瞬態(tài)計(jì)算。大渦湍流模型利用亞格子(SGS)尺度模型，以更好地捕捉瞬態(tài)流場(chǎng)中大尺度渦與小尺度渦之間的關(guān)系，SGS模型通常采用渦-黏度模型表示

(1)

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

(2)

式中，ρ為流體密度，u為流體流動(dòng)速度。

在亞格子Smagorinsky-Lily模型中，SGS湍流黏度μt可以定義為

(3)

式中，LS表示亞格子尺度的混合長(zhǎng)度，計(jì)算式為

(4)

其中，k為馮卡門(mén)常數(shù)，取0.41，d為第一層網(wǎng)格與壁面間距離，CS為Smagorinsky常數(shù)，取0.1，V為計(jì)算單元體積。

1.2 應(yīng)力計(jì)算理論模型

本文的應(yīng)力計(jì)算基于ANSYS軟件的有限元(FEM)分析理論進(jìn)行，理論模型為

σ=D×(ε-ε0)

(5)

式中，D為彈性矩陣，ε為節(jié)點(diǎn)形變，可定義為

ε=B·δe

(6)

式中，B為應(yīng)變矩陣，δe為節(jié)點(diǎn)位移，可定義為

K·δe=QT

(7)

式中，K為剛度矩陣，QT為載荷。

ε0為由溫度變化引起的變形，可定義為

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

(8)

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

2 計(jì)算模型

T型管的管道結(jié)構(gòu)如圖1所示。主、支管管徑分別為Dm與Db，主、支管流體流速分別為Vm與Vb，主管流體動(dòng)量為Mm，支管流體動(dòng)量為Mb，主、支管流體動(dòng)量比MR計(jì)算式如式(11)所示，其中，ρm、ρb分別為主、支管流體的密度。

圖1 T型管道結(jié)構(gòu)模型Fig.1 Structural model of the T-junction

(9)

(10)

(11)

(12)

式中，Ti為流體溫度，Thot與Tcold分別為熱流體與冷流體溫度。

(13)

(14)

在本文中，主支管流體溫差保持在100 ℃不變，通過(guò)改變主支管流體管徑比或速度比的方式來(lái)改變流體的動(dòng)量比，具體的計(jì)算工況如表1所示。表中，case1、case2及case3工況保持支管入口質(zhì)量流量一致，改變主支管管徑比與速度比，以實(shí)現(xiàn)主支管流體動(dòng)量比的改變，使得T型管內(nèi)流體的攪混流型由偏轉(zhuǎn)流型轉(zhuǎn)變?yōu)楸诿媪餍?；case4工況則通過(guò)進(jìn)一步增大主支管流體速度比，來(lái)研究動(dòng)量比對(duì)管道熱疲勞的影響。

表1 不同主支管動(dòng)量比的計(jì)算工況Table 1 Calculation conditions of the momentum ratio ofdifferent main and branch pipes

管道內(nèi)監(jiān)測(cè)點(diǎn)布置方案如圖2所示，共布置如圖所示的P1～P9共9個(gè)監(jiān)測(cè)面，面與面之間間隔100 mm，每個(gè)監(jiān)測(cè)面上有L1～L8共8個(gè)監(jiān)測(cè)方向，每個(gè)方向布有內(nèi)壁面監(jiān)測(cè)點(diǎn)以及據(jù)內(nèi)壁面1 mm處的流體監(jiān)測(cè)點(diǎn)。同時(shí)，管道主支管交匯處內(nèi)壁面與主管相交的相貫線為PA監(jiān)測(cè)線，其上布置PA1～PA4共4個(gè)應(yīng)力監(jiān)測(cè)點(diǎn)，與支管相交的相貫線為PB監(jiān)測(cè)線，其上布置PB1～PB4共4個(gè)應(yīng)力監(jiān)測(cè)點(diǎn)，用以監(jiān)測(cè)主支管交匯區(qū)熱應(yīng)力的波動(dòng)情況。為便于后續(xù)分析，以監(jiān)測(cè)點(diǎn)P3-L5-I/F為例，其表示在P3監(jiān)測(cè)面L5方向上，內(nèi)壁面/近壁面流體的監(jiān)測(cè)點(diǎn)。

圖2 監(jiān)測(cè)點(diǎn)布置示意圖Fig.2 Location of the monitoring points

為保證計(jì)算的準(zhǔn)確性，對(duì)管道模型進(jìn)行結(jié)構(gòu)化網(wǎng)格劃分。為滿足大渦模擬的計(jì)算要求，劃分網(wǎng)格時(shí)控制邊界層第一層的網(wǎng)格高度，保證y+值小于5，同時(shí)進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證。網(wǎng)格劃分方案如表2所示。

表2 網(wǎng)格劃分方案Table 2 Meshing scheme of the tee pipe

圖3為在case4工況條件下，P3-L5-I監(jiān)測(cè)點(diǎn)的網(wǎng)格無(wú)關(guān)性驗(yàn)證結(jié)果。如圖所示，相較于網(wǎng)格數(shù)量最少的mesh1網(wǎng)格的數(shù)值計(jì)算結(jié)果，網(wǎng)格數(shù)量更多的mesh2與mesh3網(wǎng)格的數(shù)值計(jì)算結(jié)果相差更小。因此mesh2網(wǎng)格能在滿足計(jì)算精度的同時(shí)提高計(jì)算效率，故選擇該網(wǎng)格作為本文數(shù)值模擬的網(wǎng)格劃分方案。

圖3 網(wǎng)格無(wú)關(guān)性驗(yàn)證結(jié)果Fig.3 Results of the grid independence test

3 數(shù)值計(jì)算及結(jié)果分析

3.1 CFD計(jì)算

本文利用Fluent軟件對(duì)T型管道進(jìn)行數(shù)值模擬計(jì)算，通過(guò)標(biāo)準(zhǔn)k-ε湍流模型計(jì)算得到收斂的穩(wěn)態(tài)結(jié)果，以該結(jié)果作為L(zhǎng)ES模型的初始場(chǎng)進(jìn)行瞬態(tài)計(jì)算。求解選用couple耦合求解器及二階迎風(fēng)格式，主支管入口均為速度入口(velocity inlet)，出口為自由出流(outflow)，管道外壁面為絕熱壁面，內(nèi)壁面為耦合無(wú)滑移邊界。計(jì)算總時(shí)長(zhǎng)為20 s，時(shí)間步長(zhǎng)為0.005 s，數(shù)據(jù)監(jiān)測(cè)頻率為200 Hz。由于瞬態(tài)計(jì)算結(jié)果在10 s后逐漸趨于穩(wěn)定，故取15～20 s時(shí)間內(nèi)的數(shù)據(jù)進(jìn)行對(duì)比分析。

圖4為T(mén)型管在不同動(dòng)量比工況下支管入口截面的溫度分布云圖。由圖可以看出，隨著主支管流體動(dòng)量比逐漸增大，冷熱流體混合區(qū)域逐漸向支管方向發(fā)展。在4例工況中，case1與case2工況管道內(nèi)的攪混流型為偏轉(zhuǎn)流，兩種工況管道壁面的溫度幾乎沒(méi)有波動(dòng)，case3與case4為壁面流工況，其中case3工況管道溫度波動(dòng)范圍最大。

圖4 支管流體流入截面的溫度分布云圖Fig.4 Temperature distribution contours of the branch pipe fluid inflow section

為表征冷熱流體摻混區(qū)域流體的流動(dòng)狀態(tài)，比較不同工況下P3截面位置的主管速度矢量分布圖，結(jié)果如圖5所示。由圖5可以看出，在攪混流型為偏轉(zhuǎn)流條件下，流體在主管內(nèi)的摻混區(qū)域隨著動(dòng)量比的增大，從case1工況的L6方向發(fā)展到case2工況的L5方向；而在攪混流型為壁面流的case3與case4工況條件下，流體摻混發(fā)生在主管內(nèi)部的L4方向上，該現(xiàn)象是由在主支管交匯處支管上游彎管產(chǎn)生的二次流造成的。對(duì)比近壁面流體速度矢量的大小可以發(fā)現(xiàn)，壁面流工況下近壁面流體的速度及摻混區(qū)域均明顯大于偏轉(zhuǎn)流工況。

圖5 P3截面速度矢量分布圖Fig.5 Velocity vector distribution of P3 section

由于支管上游彎管的影響，在彎管產(chǎn)生的二次流的作用下，摻混區(qū)域內(nèi)流體從L6方向向L5方向旋轉(zhuǎn)，故以L5方向?yàn)槔?，比較不同動(dòng)量比工況下，近壁面流體的無(wú)量綱時(shí)均溫度與無(wú)量綱均方根溫度的分布曲線，如圖6、圖7所示。

圖6 各工況L5方向近壁面流體的無(wú)量綱時(shí)均溫度曲線Fig.6 Dimensionless time average temperature curves of the fluid near the wall in the L5 direction under various conditions

圖7 各工況L5方向近壁面流體的無(wú)量綱均方根溫度曲線Fig.7 Dimensionless root mean square temperature curves of the fluid near the wall in the L5 direction under various conditions

圖6為不同動(dòng)量比工況下，L5方向近壁面流體在各監(jiān)測(cè)點(diǎn)上的無(wú)量綱時(shí)均溫度大小。如圖6所示，偏轉(zhuǎn)流工況(case1與case2)下在近壁面位置的無(wú)量綱時(shí)均溫度均近似為1，即冷熱流體摻混主要發(fā)生在流體中心，近壁面流體溫度保持不變；而壁面流工況(case3和case4)下近壁面流體的溫度沿主管軸向先下降再升高，在P3截面位置最低，其中case3工況P3截面流體監(jiān)測(cè)點(diǎn)的無(wú)量綱時(shí)均溫度最低，為0.42。

圖7為不同動(dòng)量比工況下，L5方向近壁面流體在各監(jiān)測(cè)點(diǎn)上的無(wú)量綱均方根溫度大小。如圖7所示，對(duì)于偏轉(zhuǎn)流工況，case1與case2工況近壁面流體的無(wú)量綱均方根溫度近似為0，即偏轉(zhuǎn)流工況近壁面位置流體的溫度波動(dòng)很小。而在壁面流工況中，case3與case4工況近壁面流體的無(wú)量綱均方根溫度明顯增大，且case3工況溫度的波動(dòng)程度要明顯大于case4工況，其主要原因?yàn)?，在壁面流工況中，隨著動(dòng)量比逐漸增大，冷熱流體摻混區(qū)域會(huì)逐漸靠近主管上壁面，使得近壁面流體溫度的波動(dòng)區(qū)域減小。無(wú)量綱均方根溫度最大位置位于case3工況P2截面的流體監(jiān)測(cè)點(diǎn)處，為0.18。

本文中，case1、case2及case3工況主支管流體溫差及流量均未發(fā)生改變，而case4工況相對(duì)于case3工況僅改變流體的速度比，根據(jù)式(11)中對(duì)動(dòng)量比的定義，影響這3例工況動(dòng)量比變化的主要因素為主支管流體的速度比，隨著速度比的增大，主支管流體間動(dòng)量比增大，冷熱流體的摻混區(qū)域逐漸靠近主管上壁面。

3.2 CFD- FEM流固耦合計(jì)算

在T型管結(jié)構(gòu)中，最易發(fā)生疲勞破壞的位置主要出現(xiàn)在主支管相貫區(qū)以及冷熱流體摻混區(qū)域。在本文中，冷熱流體發(fā)生摻混的位置位于T型管相貫處，故相貫區(qū)的PA1～PA4、PB1～PB4這8個(gè)監(jiān)測(cè)點(diǎn)是應(yīng)力計(jì)算時(shí)要重點(diǎn)關(guān)注的位置。

根據(jù)Fluent數(shù)值模擬的計(jì)算結(jié)果，將15～20 s時(shí)間內(nèi)管道的瞬時(shí)溫度載荷與內(nèi)壁面瞬時(shí)壓力載荷導(dǎo)出，通過(guò)workbench平臺(tái)加載到固體計(jì)算域中，固體計(jì)算域的約束條件及所施加的載荷如表3所示，管道材料選用316L不銹鋼。

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

本文采用第三強(qiáng)度理論將管道溫度載荷與壓力載荷帶來(lái)的多軸應(yīng)力轉(zhuǎn)化為單軸應(yīng)力，第三強(qiáng)度理論的當(dāng)量應(yīng)力σ定義為

σ=σ1-σ3

(15)

式中，σ1為第一主應(yīng)力，σ3為第三主應(yīng)力。

(16)

(17)

式中，N為總的計(jì)算時(shí)間步數(shù)，σi為每個(gè)時(shí)間步計(jì)算得到的監(jiān)測(cè)點(diǎn)應(yīng)力。

相貫區(qū)各監(jiān)測(cè)點(diǎn)的時(shí)均應(yīng)力大小如圖8所示。在相貫區(qū)與主管相交的PA監(jiān)測(cè)線上，除位于上游的PA2監(jiān)測(cè)點(diǎn)處時(shí)均應(yīng)力不隨T型管內(nèi)流體攪混流型的變化而出現(xiàn)明顯改變外，其余各監(jiān)測(cè)點(diǎn)在壁面流工況下的時(shí)均應(yīng)力均大于偏轉(zhuǎn)流工況下的時(shí)均應(yīng)力。其主要原因?yàn)楸诿媪鞴r下近壁面流體的溫度波動(dòng)大，使得其內(nèi)壁面溫度波動(dòng)增大，時(shí)均應(yīng)力大于偏轉(zhuǎn)流工況，其中，case4工況PA3監(jiān)測(cè)點(diǎn)的時(shí)均應(yīng)力最大，為215.82 MPa。由于這4例工況冷熱流體摻混區(qū)域主要集中在主管下游，位于上游的PA2監(jiān)測(cè)點(diǎn)的溫度波動(dòng)較小，使得該位置的應(yīng)力主要受主支管溫度梯度大小的影響，本文中各工況溫差均為100 ℃，溫度梯度分布相近，故在PA2點(diǎn)的時(shí)均應(yīng)力受流型的影響小。在相貫區(qū)與支管相交的PB監(jiān)測(cè)線上，由于T型管內(nèi)沒(méi)有湍流穿透現(xiàn)象發(fā)生，支管溫度波動(dòng)小，其時(shí)均應(yīng)力主要受溫度梯度的影響，故各工況監(jiān)測(cè)點(diǎn)的時(shí)均應(yīng)力大小基本一致。

圖8 相貫區(qū)各監(jiān)測(cè)點(diǎn)的時(shí)均應(yīng)力Fig.8 Time average stress of each monitoring point in the intersecting area

圖9為各工況相貫區(qū)監(jiān)測(cè)點(diǎn)的均方根應(yīng)力圖。如圖9所示，case1與case2工況在PB監(jiān)測(cè)線上的應(yīng)力波動(dòng)大于PA監(jiān)測(cè)線，case3工況下PA監(jiān)測(cè)線應(yīng)力波動(dòng)劇烈，其中PA1監(jiān)測(cè)點(diǎn)的應(yīng)力波動(dòng)最大，為31.66 MPa；case4工況應(yīng)力波動(dòng)程度相比于其他各工況波動(dòng)程度較小，僅在PA3位置應(yīng)力波動(dòng)達(dá)到6.77 MPa。即在case3工況下，相貫區(qū)靠近主管位置的應(yīng)力波動(dòng)最大，其原因是流體的溫度波動(dòng)主要發(fā)生在這一位置，從而使得壁面溫度的波動(dòng)更劇烈，應(yīng)力波動(dòng)也隨之增大。

圖9 相貫區(qū)各監(jiān)測(cè)點(diǎn)均方根應(yīng)力Fig.9 Root mean square stress of each monitoring point in the intersecting area

3.3 熱疲勞分析

根據(jù)有限元分析方法得到各點(diǎn)應(yīng)力的波動(dòng)信號(hào)后，對(duì)應(yīng)力信號(hào)進(jìn)行預(yù)處理，保留信號(hào)的峰值與谷值，并對(duì)預(yù)處理信號(hào)進(jìn)行雨流計(jì)數(shù)分析[9]，統(tǒng)計(jì)在15～20 s時(shí)間段內(nèi)，不同平均應(yīng)力下各交變應(yīng)力幅的循環(huán)次數(shù)，得到雨流矩陣N。

由于計(jì)算得到的應(yīng)力循環(huán)載荷為非對(duì)稱循環(huán)載荷，故本文利用Goodman平均應(yīng)力修正曲線對(duì)應(yīng)力載荷進(jìn)行修正，得到其等效的對(duì)稱循環(huán)應(yīng)力載荷，如式(18)所示

(18)

根據(jù)美國(guó)核管理委員會(huì)制訂的NUREG/CR-6909文件[10]中的奧氏體不銹鋼疲勞壽命S-N曲線，求得在不同交變應(yīng)力下管道的許用循環(huán)次數(shù)矩陣Nf。雨流矩陣與許用循環(huán)次數(shù)矩陣相除即可求得管道的疲勞損傷矩陣。管道監(jiān)測(cè)點(diǎn)上的疲勞損傷矩陣D定義為

(19)

根據(jù)線性累計(jì)疲勞損傷準(zhǔn)則，將D中各元素相加可求得管道的累計(jì)疲勞損傷D，當(dāng)D大于1時(shí)，即可認(rèn)為發(fā)生了疲勞損傷，故疲勞壽命L可定義為

(20)

式中，τ為總的計(jì)算時(shí)長(zhǎng)，τ=5 s。

管道的疲勞壽命由最先發(fā)生疲勞位置的壽命決定，通過(guò)對(duì)比各工況下所有應(yīng)力監(jiān)測(cè)點(diǎn)的壽命，找出最先發(fā)生疲勞的位置，以判斷不同動(dòng)量比工況下管道的疲勞壽命。各工況管道的疲勞壽命如表4所示。

表4 各工況下的疲勞壽命Table 4 Fatigue life for different calculation conditions

在偏轉(zhuǎn)流工況case1與case2中，管道的預(yù)測(cè)疲勞壽命分別為1.47×109年和5.40×1015年，在正常使用時(shí)間內(nèi)不會(huì)出現(xiàn)熱疲勞現(xiàn)象，而在壁面流工況case3與case4中，管道的預(yù)測(cè)壽命分別為8.81×10-4年及4.62×10-2年。由此可知，降低主支管流體動(dòng)量比，使攪混流型從壁面流轉(zhuǎn)變?yōu)槠D(zhuǎn)流能有效改善管道的熱疲勞現(xiàn)象，但在攪混流型相同時(shí)，通過(guò)降低動(dòng)量比的方式并不能改善管道的熱疲勞。

4 結(jié)論

(1)在偏轉(zhuǎn)流工況條件下，T型管內(nèi)流體溫度的波動(dòng)區(qū)域相比于壁面流工況條件更接近主管中心，通過(guò)降低動(dòng)量比，使T型管內(nèi)流體攪混流型從壁面流發(fā)展為偏轉(zhuǎn)流，能有效降低管道近壁面流體的溫度波動(dòng)程度。

(2)相較于壁面流工況，采用偏轉(zhuǎn)流工況能通過(guò)減小近壁面流體的溫度波動(dòng)程度來(lái)降低管壁熱應(yīng)力的波動(dòng)程度，從而改善T型管道的熱疲勞現(xiàn)象，提高管道的熱疲勞壽命。

(3)只有當(dāng)動(dòng)量比降低至T型管內(nèi)流體的攪混流型發(fā)生變化，即由壁面流型發(fā)展為偏轉(zhuǎn)流型后，管道的熱疲勞現(xiàn)象才能得到改善，流型未發(fā)生改變時(shí)，調(diào)整動(dòng)量比并不能改善管道的熱疲勞現(xiàn)象。