周照春 陳柏宇 張琪琪 張經(jīng)豪 田 源 盧 濤

(北京化工大學(xué) 機電工程學(xué)院， 北京 100029)

引 言

T型管是核電管路系統(tǒng)中最常見的一種管道結(jié)構(gòu)形式。T型管道的熱疲勞是由冷熱流體交匯引起的，由于冷熱流體存在溫差，管道在溫度波動的作用下產(chǎn)生熱應(yīng)力，最終導(dǎo)致出現(xiàn)熱疲勞現(xiàn)象，對核電的安全運行造成嚴(yán)重威脅。

Lee等[1]運用大渦模擬(LES)對T型管道內(nèi)的冷熱流體摻混進(jìn)行數(shù)值模擬，并得出冷熱流體湍流混合所引起的強化對流換熱是管道熱疲勞失效的主要原因。Zhang等[2]基于計算流體動力學(xué)- 有限元分析方法(CFD- FEM)和雨流計數(shù)法對T型管內(nèi)的冷熱流體混合進(jìn)行了數(shù)值分析，構(gòu)建了T型管冷熱流體湍流混合的高周熱疲勞計算方法。Hu等[3]利用LES模型對T型管內(nèi)冷熱流體的摻混過程進(jìn)行了數(shù)值模擬，并根據(jù)主管與支管的動量比MR將混合流型分為壁面射流、偏轉(zhuǎn)射流和沖擊射流。Wegner等[4]通過改變T型支管角度提高了湍流的混合效果，并發(fā)現(xiàn)傾角的變化影響了渦結(jié)構(gòu)和二次流，從而產(chǎn)生了不同的湍流混合效果。張越[5]對壓器波動管進(jìn)行大渦模擬，使用瞬態(tài)界面單向流固耦合方法對管道進(jìn)行應(yīng)力分析，并對管道進(jìn)行了熱疲勞評估，計算了管道的疲勞壽命。Kamaya[6-7]針對直管、彎管、T型管中的流動與傳熱現(xiàn)象進(jìn)行簡化，實現(xiàn)了管道的熱- 流- 固耦合分析。

綜上所述，目前對T型管道內(nèi)冷熱流體摻混過程中的流動與傳熱以及熱疲勞評估已取得了很大的研究進(jìn)展，但關(guān)于支管上游帶有彎管的T型管結(jié)構(gòu)尚無文獻(xiàn)報道。本文改變了T型管的支管結(jié)構(gòu)，研究支管上游不同曲率的彎管結(jié)構(gòu)對管道的溫度波動及熱疲勞的影響。利用Fluent軟件對T型管內(nèi)冷熱流體的湍流混合過程進(jìn)行CFD數(shù)值模擬，對模擬結(jié)果進(jìn)行分析，得出湍流混合過程的流動與傳熱特性；利用CFD- Post軟件提取溫度場和壓力場的數(shù)據(jù)，并通過ANSYS Workbench平臺將其動態(tài)加載到固體域模型，開展CFD- FEM瞬態(tài)流固耦合計算，得到固體熱應(yīng)力的分布及波動信息，并對危險位置進(jìn)行了熱疲勞分析。

1 計算模型

本文所用的物理模型如圖1(a)所示，主管內(nèi)為高速高溫的流體，支管內(nèi)為低速低溫的流體。在支管位置增加了不同曲率(曲率k=0.006，0.01，0.02)的彎管，對3種曲率支管彎管的T型管進(jìn)行數(shù)值模擬及分析，研究不同曲率彎管對T型管冷熱流體混合熱疲勞的影響。T型管整體的結(jié)構(gòu)參數(shù)如表1所示，其中Lm為主管總長，L1為主管下游長，Lb為支管長，d1為管外徑，d2為管內(nèi)徑,e為壁厚。對整個物理模型的固體域和流體域同時進(jìn)行結(jié)構(gòu)化網(wǎng)格劃分，并對主支管相貫區(qū)及邊界層位置進(jìn)行加密處理。T型管結(jié)構(gòu)化網(wǎng)格模型如圖1(b)所示，網(wǎng)格數(shù)量為1 358 638，第一層網(wǎng)格高度為0.05 mm，邊界層數(shù)為10，邊界層網(wǎng)格增長率為1.1。

圖1 T型管計算模型Fig.1 Model of the simulation of T-junction

表1 T型管結(jié)構(gòu)參數(shù)Table 1 Structural parameters of T-junction

T型管內(nèi)冷熱流體的摻混是一個復(fù)雜、多變的過程，為了更好地捕捉管內(nèi)湍流混和瞬時波動及流動傳熱的特征細(xì)節(jié)，本文采用大渦模擬方法對T型管內(nèi)冷熱流體的摻混過程進(jìn)行數(shù)值模擬[8-9]，模擬過程考慮重力和浮升力的影響[10]。在大渦模擬動量方程的浮升力項中，對水的各項參數(shù)進(jìn)行變物性設(shè)置。利用NIST軟件對水的密度隨溫度的變化進(jìn)行二階多項式擬合，如式(1)所示。

ρ=-0.002 3T2+0.988 7T+909.5

(1)

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

2 數(shù)值模擬

2.1 流場模擬

流體計算的邊界條件如表2所示。計算過程中對水的密度、導(dǎo)熱系數(shù)、動力黏度等進(jìn)行變物性設(shè)置，并設(shè)置操作壓力為7 MPa，將主管入口邊界條件進(jìn)行初始化。首先選用k-ε模型將T型管內(nèi)冷熱流體混合過程計算至穩(wěn)態(tài)收斂，將穩(wěn)態(tài)收斂的計算結(jié)果作為LES瞬態(tài)計算的初始場。為了更好地捕捉冷熱流體摻混瞬時熱波動的特征細(xì)節(jié)，使用couple耦合求解器與二階迎風(fēng)式進(jìn)行求解。瞬態(tài)計算的時間步長為0.005 s，溫度監(jiān)測點的監(jiān)測頻率為200 Hz，計算總時長為20 s。

表2 流場計算邊界條件Table 2 Boundary conditions in the CFD simulation

如圖2所示，為定量分析流場的溫度波動情況，沿主管軸向每80 mm設(shè)置一個監(jiān)測面(P1～P8)，P0截面位于主支管相貫區(qū)，P1截面距離支管軸線50 mm；每個監(jiān)測面間隔30°設(shè)置一個監(jiān)測方向(N1～N12),每個監(jiān)測方向上設(shè)置2個監(jiān)測點，即內(nèi)壁面(I)和近壁面(F)1 mm處的流體監(jiān)測點，如P2- N4- F- T(V)為P2截面N4方向近壁面流體溫度(速度)監(jiān)測點。

圖2 監(jiān)測點的設(shè)置Fig.2 The arrangement of monitoring points

2.2 網(wǎng)格無關(guān)性驗證

為了排除網(wǎng)格大小帶來的計算誤差，對網(wǎng)格模型進(jìn)行無關(guān)性驗證，其結(jié)構(gòu)化網(wǎng)格劃分方案如表3所示。對4種網(wǎng)格進(jìn)行數(shù)值計算，得到管內(nèi)P7截面N6方向流體的平均溫度和P4截面N4方向流體的平均速度隨網(wǎng)格量增加所產(chǎn)生的變化如圖3所示。從圖3可以看出，Mesh1、Mesh2的計算結(jié)果與Mesh3、Mesh4相差較大，Mesh3與Mesh4兩者的計算結(jié)果相差較小，考慮到計算的精度和效率，采用Mesh3作為后續(xù)數(shù)值計算的網(wǎng)格劃分方案。

表3 網(wǎng)格無關(guān)性驗證劃分方案Table 3 Meshing scheme of the tee pipe for the gridindependence test

圖3 網(wǎng)格無關(guān)性驗證Fig.3 The grid independence test

2.3 流固耦合計算

流固耦合計算需要設(shè)置管道固體材料的屬性，本文T型管道所選擇材料為316L不銹鋼，在實際運行中管道的形變很小，可以忽略其對流場的影響，因此采用一種單向瞬態(tài)流固耦合的計算方法。將流場瞬態(tài)計算所得到的內(nèi)壁面的壓力場和固體域的溫度場通過批處理的方法按照時間序列進(jìn)行排列，并將其動態(tài)加載到固體計算域作為固體計算的邊界條件和載荷。流固耦合計算的時間步長為0.01 s，計算時長為5 s(15～20 s)。固體監(jiān)測點的設(shè)置如圖2所示，在流體監(jiān)測點的基礎(chǔ)上增加了主支管相貫區(qū)內(nèi)壁面位置的8個監(jiān)測點(A1、A2、B1、B2、C1、C2、D1、D2)。在流固耦合計算過程中，提取監(jiān)測點的應(yīng)力波動信號，用于后續(xù)的應(yīng)力分析和熱疲勞分析。固體計算的邊界條件如表4所示。

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

3 計算結(jié)果及分析

3.1 流場分析

3.1.1溫度場

對3個不同曲率的模型的軸向溫度進(jìn)行無量綱處理后，發(fā)現(xiàn)其分布和波動具有一致性，冷、熱流體在主支管相貫區(qū)位置匯合，在支管彎管的影響下，熱交換區(qū)域發(fā)生偏轉(zhuǎn)，熱交換區(qū)域主要位于主管下游N4～N8方向，其中溫度變化最為明顯的是N7方向。溫度波動是由冷熱流體摻混引起的，溫度波動主要發(fā)生在主管下游N3～N8方向，其中溫度波動最劇烈的是N4方向。3種不同曲率支管彎管T型管沿N7方向的無量綱時均溫度和沿N4方向的無量綱均方根溫度的軸向分布如圖4所示。

圖4 無量綱溫度軸向分布Fig.4 Axial distributions of

由圖4(a)可以看出，沿主管N7方向3種結(jié)構(gòu)的溫度均先降低后升高。在靠近相貫區(qū)的位置由于支管低溫流體的流入，主管溫度降低，而在遠(yuǎn)離相貫區(qū)的位置由于冷熱流體充分混合，主管溫度上升，其中在k=0.006的支管彎管結(jié)構(gòu)中流體的偏轉(zhuǎn)更為明顯，這是由于該彎管曲率小，低溫流體在主管發(fā)生了很大的周向偏轉(zhuǎn)，冷熱流體充分混合。圖4(b)對比了3種結(jié)構(gòu)溫度波動最劇烈位置(N4方向)的無量綱均方根溫度沿主管軸向的變化情況，可以看出沿主管軸向3種結(jié)構(gòu)的溫度波動均呈先增大后減小的趨勢。低溫流體在相貫區(qū)與高溫流體匯合，冷流體沖入主管，由于主管流速大，熱交換區(qū)從主管軸心位置向內(nèi)壁面偏轉(zhuǎn)，在P2截面處熱交換區(qū)偏轉(zhuǎn)到內(nèi)壁面位置。因此從P1到P2截面近壁面溫度波動逐漸增大。隨著冷熱流體的充分混合，管內(nèi)流體溫度逐漸趨于穩(wěn)定，溫度波動的劇烈程度也逐漸降低,溫度波動最劇烈的區(qū)域在P2N4處。當(dāng)支管彎管的曲率較大時，低溫流體在主管內(nèi)的周向偏轉(zhuǎn)小，熱交換區(qū)集中，導(dǎo)致溫度波動較大。因此隨著曲率的增大，溫度波動的劇烈程度增大。

3.1.2速度場

對不同曲率支管彎管T型管道內(nèi)冷熱流體的混合過程進(jìn)行研究。支管彎管能夠影響主管內(nèi)摻混流體的周向速度。20 s時3種結(jié)構(gòu)在P2截面的速度矢量圖如圖5所示，可以看出，周向速度主要位于支管入口區(qū)域(N4和N6)，低溫流體流經(jīng)曲率小的彎管會產(chǎn)生較大的周向速度，因此隨著曲率的增大周向速度減小。

圖5 P2截面周向速度矢量圖(20 s)Fig.5 Vector of circumferential velocity of P2 at 20 s

為了分析和研究主管內(nèi)流體速度分布和波動的特征規(guī)律，本文引入時均速度Vmean和均方根速度Vrms，時均速度用來表征速度的分布規(guī)律，均方根速度用來表征速度的波動程度。

圖6為主管下游N5方向的速度的軸向分布情況。N5為支管入口位置，在此方向的速度變化及波動最為明顯。由圖6(a)可以看出，3種結(jié)構(gòu)的速度分布基本一致，均沿主管軸向逐漸增大。由于支管的低速流體在靠近支管入口位置占主導(dǎo)地位，因此近壁面流速相對較小；在遠(yuǎn)離支管入口位置，近壁面流體逐漸恢復(fù)為高速流體。圖6(b)為3種結(jié)構(gòu)在N5方向沿主管軸向的均方根速度變化情況，可以看出3種結(jié)構(gòu)的速度大小波動趨勢一致，均隨著主管軸向先增大后減小，這與溫度波動具有一致性。支管低溫流體與主管高溫流體在相貫區(qū)混合，由于高溫流體速度較大，熱交換區(qū)由主管軸心位置向內(nèi)壁面偏移，在P2截面處熱交換區(qū)偏移至內(nèi)壁面，速度波動最為劇烈，隨著冷熱流體的充分混合，速度波動逐漸趨于平穩(wěn)。

圖6 速度的軸向分布(N5)Fig.6 Axial distributions of Vmean and Vrms(N5)

3.2 流固耦合分析

3.2.1應(yīng)力

對3種結(jié)構(gòu)經(jīng)流固耦合計算得到的應(yīng)力結(jié)果進(jìn)行分析。圖7為20 s時3種結(jié)構(gòu)在P0截面的等效應(yīng)力σ分布圖?？梢钥闯?，熱應(yīng)力集中區(qū)主要位于相貫區(qū)內(nèi)壁面。熱應(yīng)力是由溫度波動引起的，冷熱流體在主支管相貫區(qū)交匯，產(chǎn)生較大的熱應(yīng)力；隨著曲率的增大，熱應(yīng)力集中現(xiàn)象逐漸減弱,這與溫度及速度波動劇烈程度的變化趨勢是一致的。

圖7 P0截面的等效應(yīng)力分布(20 s)Fig.7 The equivalent stress distribution of P0 at 20 s

3.2.2熱疲勞

如圖8所示，提取流固耦合計算結(jié)果中監(jiān)測點D1(曲率k=0.006)的應(yīng)力波動信號，并對其進(jìn)行預(yù)處理(只保留拐點信息)。

圖8 D1監(jiān)測點的應(yīng)力波動信號Fig.8 Stress fluctuation information at the point D1

圖9為監(jiān)測點D1的熱疲勞損傷率評估結(jié)果。利用雨流計數(shù)法對預(yù)處理后的熱應(yīng)力波動信號進(jìn)行統(tǒng)計學(xué)分析，處理后得到如圖9(a)所示的雨流矩陣N及實際循環(huán)次數(shù)n，利用Goodman曲線進(jìn)行應(yīng)力修正，將等效熱應(yīng)力轉(zhuǎn)化為等效對稱循環(huán)應(yīng)力σ-1，并得到如圖9(b)所示的修正應(yīng)力矩陣NG。數(shù)值模擬的工況較為理想化，而在實際工況中存在很多復(fù)雜的因素(環(huán)境、焊縫等)，因此本文利用安全系數(shù)K修正對稱循環(huán)應(yīng)力，并放大載荷使結(jié)果更加保守，K取值為2。Goodman曲線為

圖9 D1監(jiān)測點的熱疲勞損傷率評估Fig.9 Estimation of the thermal fatigue damage rate at the point D1

(2)

利用應(yīng)力- 壽命曲線(S- N曲線)進(jìn)行熱疲勞評估。根據(jù)奧氏體不銹鋼的疲勞壽命S- N曲線插值獲得相關(guān)等效應(yīng)力下的許用循環(huán)次數(shù)，當(dāng)?shù)刃ΨQ循環(huán)應(yīng)力大于100 MPa時，滿足S- N曲線的變化規(guī)律；當(dāng)?shù)刃ΨQ循環(huán)應(yīng)力小于100 MPa時，許用循環(huán)次數(shù)無限增長，即不發(fā)生疲勞，但實際工況很難達(dá)到理論工況的標(biāo)準(zhǔn)，因此假定此種情況也滿足S- N曲線的變化規(guī)律。對等效對稱循環(huán)應(yīng)力小于100 MPa時的疲勞壽命S- N曲線進(jìn)行擬合，擬合公式為[5]

(3)

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

(4)

(5)

(6)

式中，Dτ,i,j為疲勞損傷率矩陣Dτ中的各元素，τ為熱應(yīng)力波動信號時長，τ=5 s。

對提取出的所有監(jiān)測點的熱應(yīng)力波動信號進(jìn)行熱疲勞損傷率評估，發(fā)現(xiàn)無量綱累計疲勞損傷率隨支管彎管曲率的增大呈減小的趨勢，k=0.02結(jié)構(gòu)彎管的無量綱累計疲勞損傷最小，并得出危險位置為k=0.006結(jié)構(gòu)的B1和D1監(jiān)測點處，其無量綱疲勞累計損傷率分別為52.98和52.27。

4 結(jié)論

(1)利用大渦模擬對曲率k為0.006、0.01、0.02這3種結(jié)構(gòu)的T型管冷熱流體的湍流混合過程進(jìn)行了數(shù)值預(yù)測，結(jié)果表明3種結(jié)構(gòu)的溫度分布和波動的趨勢相似，曲率大的支管彎管結(jié)構(gòu)產(chǎn)生了劇烈的溫度波動。

(2)低溫流體流經(jīng)彎管進(jìn)入主管后，在主管內(nèi)產(chǎn)生周向速度，且隨著支管彎管曲率的增大，周向速度逐漸減小。

(3)利用熱- 流- 固耦合計算方法，獲得了熱應(yīng)力的分布及波動情況，發(fā)現(xiàn)熱應(yīng)力集中主要位于主支管相貫區(qū)，熱應(yīng)力集中現(xiàn)象隨支管彎管曲率的增大而減弱。

(4)通過對熱應(yīng)力波動信號的統(tǒng)計學(xué)分析及處理完成了對管道的熱疲勞評估，發(fā)現(xiàn)無量綱累計疲勞損傷率隨支管彎管曲率的增大而減小，并確定了易疲勞點出現(xiàn)在k=0.006結(jié)構(gòu)彎管的B1位置處，其無量綱疲勞累計損傷率為52.98。