陳柏宇 周照春 孫 燕 陳 爽 歐陽斌 盧 濤*

(1.北京化工大學 機電工程學院， 北京 100029; 2.中國核動力研究設計院， 成都 610200)

引 言

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

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

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

1 數(shù)學模型

1.1 LES模型

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

(1)

亞格子雷諾應力τij可定義為

(2)

式中，ρ為流體密度，u為流體流動速度。

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

(3)

式中，LS表示亞格子尺度的混合長度，計算式為

(4)

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

1.2 應力計算理論模型

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

σ=D×(ε-ε0)

(5)

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

ε=B·δe

(6)

式中，B為應變矩陣，δe為節(jié)點位移，可定義為

K·δe=QT

(7)

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

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

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

(8)

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

2 計算模型

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

圖1 T型管道結構模型Fig.1 Structural model of the T-junction

(9)

(10)

(11)

(12)

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

(13)

(14)

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

表1 不同主支管動量比的計算工況Table 1 Calculation conditions of the momentum ratio ofdifferent main and branch pipes

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

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

為保證計算的準確性，對管道模型進行結構化網(wǎng)格劃分。為滿足大渦模擬的計算要求，劃分網(wǎng)格時控制邊界層第一層的網(wǎng)格高度，保證y+值小于5，同時進行網(wǎng)格無關性驗證。網(wǎng)格劃分方案如表2所示。

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

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

圖3 網(wǎng)格無關性驗證結果Fig.3 Results of the grid independence test

3 數(shù)值計算及結果分析

3.1 CFD計算

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

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

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

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

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

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

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

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

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

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

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

3.2 CFD- FEM流固耦合計算

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

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

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

本文采用第三強度理論將管道溫度載荷與壓力載荷帶來的多軸應力轉化為單軸應力，第三強度理論的當量應力σ定義為

σ=σ1-σ3

(15)

式中，σ1為第一主應力，σ3為第三主應力。

(16)

(17)

式中，N為總的計算時間步數(shù)，σi為每個時間步計算得到的監(jiān)測點應力。

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

圖8 相貫區(qū)各監(jiān)測點的時均應力Fig.8 Time average stress of each monitoring point in the intersecting area

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

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

3.3 熱疲勞分析

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

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

(18)

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

(19)

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

(20)

式中，τ為總的計算時長，τ=5 s。

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

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

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

4 結論

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

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

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