核動(dòng)力系統(tǒng)蒸汽發(fā)生器U 形管流致振動(dòng)數(shù)值模擬

王瑩杰，王明軍,*，王樹強(qiáng)，田文喜，秋穗正，蘇光輝

（1.西安交通大學(xué)核科學(xué)與技術(shù)學(xué)院，陜西 西安 710049；2.大亞灣核電運(yùn)營管理有限責(zé)任公司運(yùn)行部，廣東 深圳 518124）

蒸汽發(fā)生器是壓水堆連接一、二回路的核心設(shè)備，其有效性關(guān)乎核反應(yīng)堆的安全穩(wěn)定運(yùn)行，是核動(dòng)力系統(tǒng)的重點(diǎn)研究設(shè)備[1]。傳熱管作為蒸汽發(fā)生器中的關(guān)鍵部件，是核動(dòng)力系統(tǒng)的一回路壓力邊界和傳熱邊界。壓水堆核電站的運(yùn)行經(jīng)驗(yàn)表明，蒸汽發(fā)生器傳熱管斷裂事故在核電廠事故中居首要地位[2]。在高溫高壓的環(huán)境中，二次側(cè)流體沖刷傳熱管束，會引起傳熱管振動(dòng)，即流致振動(dòng)現(xiàn)象，這一現(xiàn)象可能進(jìn)一步導(dǎo)致?lián)Q熱管的破裂損傷[3]，造成冷卻劑泄漏，破壞蒸汽發(fā)生器的安全性和完整性。因此，對蒸汽發(fā)生器傳熱管的流致振動(dòng)現(xiàn)象進(jìn)行研究具有重要的意義。隨著計(jì)算機(jī)技術(shù)的發(fā)展，流固耦合方法也被廣泛應(yīng)用于核電廠結(jié)構(gòu)完整性及流致振動(dòng)等現(xiàn)象研究中[4,5]，關(guān)于蒸汽發(fā)生器的雙向流固耦合研究[6,7]基于WORKBENCH 平臺，計(jì)算過程數(shù)據(jù)交換緩慢，且只分析了傳熱管的穩(wěn)態(tài)熱應(yīng)力。馮志鵬[8]利用CFD 計(jì)算軟件用戶自定義函數(shù),通過Newmark-β 方法計(jì)算圓管結(jié)構(gòu)動(dòng)力學(xué)方程求得渦激振動(dòng)圓管的振動(dòng)特性。傳熱管潛在的激勵(lì)源包括存在于管子內(nèi)側(cè)的一次側(cè)流體、機(jī)械誘發(fā)振動(dòng)，以及存在于管子外側(cè)的二次側(cè)流體作用，分析表明：引起傳熱管惡化的主要振動(dòng)原因來自于二次側(cè)流體的激勵(lì)作用[9]。本文將基于Newmark-β 方法，開發(fā)基于FLUENT 的先進(jìn)壓水堆蒸汽發(fā)生器U 形管流致振動(dòng)模型，開展蒸汽發(fā)生器傳熱管在二次側(cè)流體沖刷作用下流致振動(dòng)現(xiàn)象的數(shù)值模擬。

1 數(shù)值計(jì)算方法

1.1 控制方程及求解流程

流體域控制方程為二維不可壓縮粘性流體的連續(xù)性方程和Navier-Stokes 方程，在慣性坐標(biāo)系下，二者表達(dá)式如下：

式中：ρ——流體密度；

v——運(yùn)動(dòng)粘度系數(shù)；

——速度矢量；

p——流場壓力。

傳熱管在流場中的振動(dòng)可簡化為質(zhì)量 -彈簧 -阻尼系統(tǒng)，其兩自由度模型結(jié)構(gòu)振動(dòng)方程為：

其中：U∞——均勻來流速度；

ζ、ω0——系統(tǒng)的阻尼比和圓頻率；

M——單位長度管的質(zhì)量；

Fl和Fd——單位長度管體受到的升力和阻力；

Cl和Cd——升力系數(shù)和阻力系數(shù)；

Y、、——管橫向位移、速度和加速度；

X、、——管流向位移、速度和加速度。

本文基于松耦合分區(qū)方法求解傳熱管的流固耦合振動(dòng)，流體域計(jì)算采用基于有限體積法的CFD 通用求解程序FLUENT 完成；結(jié)構(gòu)域的求解通過在C 語言環(huán)境下自編Newmark-β 算法完成[10]，通過Fluent 用戶自定義函數(shù)（UDF）與主程序鏈接求解管的振動(dòng)響應(yīng)；網(wǎng)格域根據(jù)振動(dòng)響應(yīng)，采用動(dòng)網(wǎng)格模型不斷更新，本文數(shù)值方法實(shí)現(xiàn)流程圖如圖1 所示[11]。

圖1 流固耦合計(jì)算流程Fig.1 Fluid-structure coupling calculation process

1.2 計(jì)算模型及邊界條件

本文以典型壓水堆蒸汽發(fā)生器單根U 形管半圓形彎頭為研究對象[12]，計(jì)算模型如圖2 所示，U 形管的材料為Inconel690，其結(jié)構(gòu)和物理參數(shù)如表1 所示，當(dāng)流體從下至上流動(dòng)時(shí)，流體掠過管外壁，則從側(cè)方看，計(jì)算幾何模型可簡化為圓形截面的管外繞流，如圖3 所示。矩形計(jì)算區(qū)域均為35D×20D，模型上游來流區(qū)域?yàn)?0D，下游尾流區(qū)域?yàn)?5D，離上下邊界各為10D，D為與來流方向垂直的特征長度，即管外徑。流場采用結(jié)構(gòu)化網(wǎng)格進(jìn)行離散，對管壁附近和尾流區(qū)域等參數(shù)梯度變化較大計(jì)算敏感區(qū)進(jìn)行局部網(wǎng)格加密，網(wǎng)格劃分示意圖如圖4 所示。蒸汽發(fā)生器U 形管彎頭區(qū)域二次側(cè)流體為空泡份額0.9 左右的高溫高壓水蒸氣，計(jì)算所得流體計(jì)算域的平均物性參數(shù)如表2 所示。流體從左至右流動(dòng)（相當(dāng)于二次側(cè)流體從下至上流動(dòng)），左側(cè)邊界為速度入口，來流為均勻速度，右側(cè)邊界為壓力出口，上下邊界為自由滑移邊界，振動(dòng)管外形所形成的邊界為運(yùn)動(dòng)邊界，該邊界為數(shù)據(jù)傳遞的耦合面，設(shè)定為無滑移壁面。

表1 U 形傳熱管幾何及物理參數(shù)Table 1 Geometric and physical parameters of U-tube

續(xù)表

表2 流體物性參數(shù)Table 2 Physical parameters of fluid

圖2 單根U 形管半圓形彎頭模型Fig.2 Semicircular elbow model of single U-tube

圖3 計(jì)算模型Fig.3 Computational model

圖4 網(wǎng)格劃分結(jié)果Fig.4 Meshing results

2 數(shù)值方法驗(yàn)證

根據(jù)蒸汽發(fā)生器二次側(cè)流體速度范圍，求得雷諾數(shù)大小在3×104到3×105范圍內(nèi)，這一范圍屬于圓管繞流雷諾數(shù)的亞臨界區(qū)，管壁上分離的剪切層開始向不規(guī)則的湍流狀態(tài)轉(zhuǎn)變，一部分動(dòng)能由湍流所攜帶，出現(xiàn)周期性的漩渦脫落。CFD 中的DES 湍流模型可以有效的描述亞臨界區(qū)圓管繞流的特性，并且較于大渦模擬和雷諾時(shí)均方法在計(jì)算效率和準(zhǔn)確性上有一定優(yōu)勢[12]。本文采取DES 湍流模型開展計(jì)算，計(jì)算可得截面的升力及阻力系數(shù)，對升力及阻力系數(shù)進(jìn)行快速傅里葉變換（FFT）即可得到升力、阻力系數(shù)無量綱主頻，反應(yīng)漩渦脫落頻率的Strouhal 數(shù)St 可由式（7）求得：

式中：fs——升力系數(shù)無量綱主頻，即漩渦脫落頻率；

D——特征長度，即管外徑；

u——入口平均速度。

在計(jì)算前首先對網(wǎng)格無關(guān)性進(jìn)行驗(yàn)證，圖5所示為計(jì)算得到的阻力系數(shù)均值以及Strouhal 數(shù)隨網(wǎng)格數(shù)量的變化情況。由圖可知，當(dāng)網(wǎng)格數(shù)量達(dá)到56 742 時(shí)（第5 套網(wǎng)格），阻力系數(shù)Cd和St數(shù)幾乎不再隨著網(wǎng)格數(shù)量的增加而變化。表3為56 742 網(wǎng)格數(shù)量下的計(jì)算結(jié)果和文獻(xiàn)中（相同雷諾數(shù)下）的結(jié)果對比。由圖5 和表2 可知，第五套網(wǎng)格可以認(rèn)為是本模型的網(wǎng)格獨(dú)立性解，且驗(yàn)證了文中數(shù)值模擬的正確性。圖6 給出了阻力系數(shù)和升力系數(shù)的時(shí)間歷程曲線，可以看到阻力系數(shù)和升力系數(shù)均呈現(xiàn)規(guī)則的周期性震蕩，由于是二維模擬，沒有體現(xiàn)出繞流的三維特性，所以和文獻(xiàn)[13]的三維模擬的震蕩效果稍有差異，但圖中阻力系數(shù)和升力系數(shù)的波動(dòng)趨勢和文獻(xiàn)[13]是一樣的，同樣沒有相位差異。圖8為流體域的速度云圖，可以看到流體流過傳熱管后，在尾流區(qū)域形成周期性脫落的對稱漩渦。

圖5 網(wǎng)格無關(guān)性驗(yàn)證Fig.5 Grid independence analysis

表3 計(jì)算結(jié)果和文獻(xiàn)結(jié)果比較Table 3 Comparison of calculation results and literature results

圖6 驗(yàn)證模型升力系數(shù)和阻力系數(shù)時(shí)程圖Fig.6 Time history diagram of lift coefficient and drag coefficient for verification model

圖7 文獻(xiàn)[13]升力系數(shù)和阻力系數(shù)時(shí)程圖Fig.7 Time history diagram of lift coefficient and drag coefficient for reference [13］.

圖8 流體域速度云圖Fig.8 Velocity contour of fluid domain

3 模態(tài)分析

模態(tài)分析是研究結(jié)構(gòu)動(dòng)力學(xué)動(dòng)力特性的一種方法，一般應(yīng)用在工程振動(dòng)領(lǐng)域。模態(tài)是指機(jī)械結(jié)構(gòu)的固有振動(dòng)特性，每一個(gè)模態(tài)都有特定的固有頻率、阻尼比和模態(tài)振型[17]。為求解兩自由度結(jié)構(gòu)振動(dòng)方程，首先需要對U 形管彎管段進(jìn)行模態(tài)分析，以獲得系統(tǒng)的阻尼比和圓頻率。系統(tǒng)的模態(tài)分析只與系統(tǒng)的質(zhì)量、轉(zhuǎn)動(dòng)慣量、剛度、邊界條件有關(guān)，而與外加載荷沒有關(guān)系，由于系統(tǒng)的低階模態(tài)對振動(dòng)響應(yīng)影響較大，而高階模態(tài)可以忽略不記[18]，一般只進(jìn)行較低階次的模態(tài)分析。圖9 給出了A、B、C、D 四種不同約束方式的U 形管模型，彎管兩底端采用固支約束，圓弧上約束為簡支約束，分別模擬具有不同數(shù)目防振條的U 形管約束。表4 給出了四種模型模態(tài)分析得到的U 形管前四階振動(dòng)頻率。由表可知，隨著防振條數(shù)目增加，U 形管一階振動(dòng)頻率增大，當(dāng)U 形管沒有防振條時(shí)（A 模型），一階振動(dòng)頻率和二階振動(dòng)頻率差別較大；當(dāng)U 形管增加防振條固定時(shí)，一階振動(dòng)頻率和二階振動(dòng)頻率差別較小，且均出現(xiàn)兩兩相鄰階數(shù)頻率相近的現(xiàn)象。

圖9 不同約束方式的U 形管模型Fig.9 Model of U-tube with different constraints

表4 模態(tài)分析結(jié)果Table 4 Results of modal analysis Hz

4 結(jié)果分析

一般情況下，模型一階振動(dòng)模態(tài)對結(jié)構(gòu)影響最大，因此本文取U 形管一階振動(dòng)模態(tài)分析結(jié)果求解兩自由度振動(dòng)方程。根據(jù)第三部分模態(tài)分析結(jié)果以及傳熱管物性及結(jié)構(gòu)參數(shù)確定公式（3）、式（4）中的系數(shù)，再通過Newmark-β數(shù)值方法，求解兩自由度的振動(dòng)方程，實(shí)現(xiàn)傳熱管渦激振動(dòng)現(xiàn)象的模擬。圖10 所示為A 模型在雷諾數(shù)等于1×105時(shí)，傳熱管渦激振動(dòng)的升力系數(shù)時(shí)程圖，圖11（a）、圖11（b）分別為渦激振動(dòng)阻力系數(shù)、升力系數(shù)經(jīng)過傅里葉變換得到的頻譜圖，由圖可知，阻力系數(shù)和升力系數(shù)都存在振動(dòng)的主頻率和一個(gè)附加頻率，這說明渦激振動(dòng)傳熱管漩渦脫落頻率同時(shí)受到傳熱管自身固有頻率和激振力頻率的影響，阻力系數(shù)振動(dòng)變化主頻率為升力系數(shù)振動(dòng)變化主頻率的2 倍。圖12 為傳熱管流向、橫向的無量綱位移時(shí)程曲線，由圖可知，在渦激振蕩作用下，傳熱管在流向方向上產(chǎn)生位移，且在一定位移范圍內(nèi)來回振動(dòng)，而在橫向方向上下振動(dòng)，位移均值為零。圖13 為傳熱管橫向位移頻譜圖，其主要頻率大小約為 20 Hz，這與渦激振動(dòng)傳熱管的渦脫頻率、固定傳熱管的渦脫頻率接近；橫向振動(dòng)幅值較小，說明傳熱管渦激振動(dòng)對流動(dòng)擾動(dòng)不大。圖14 為不同約束條件下，傳熱管的橫向振動(dòng)無量綱位移時(shí)程曲線圖，由圖可知，是否安裝防振條對傳熱管的橫向位移影響較大，B 型約束傳熱管橫向振動(dòng)位移為A 型約束傳熱管橫向振動(dòng)位移的1/100，而防振條的數(shù)目對橫向振動(dòng)位移的影響不大，這是由于本文只研究了U 形管半彎頭的最頂部圓截面，以及只考慮了一階振動(dòng)模態(tài)。圖15 展示了不同雷諾數(shù)條件對A 型約束傳熱管渦激振動(dòng)的阻力系數(shù)、升力系數(shù)、傳熱管橫向位移、流向位移影響時(shí)程圖。由圖15（a）、圖15（b）可知，流動(dòng)的升力系數(shù)和阻力系數(shù)都呈現(xiàn)周期性變化現(xiàn)象，隨著雷諾數(shù)的增加，阻力系數(shù)和升力系數(shù)的周期都逐漸變小，同時(shí)阻力系數(shù)振動(dòng)的幅值越來越小，升力系數(shù)振動(dòng)幅值基本不變。由圖15（c）、圖（d）可知，隨著雷諾數(shù)的增加，傳熱管的流向位移逐漸增大，橫向位移也逐漸增大，說明傳熱管的振動(dòng)現(xiàn)象受流體湍流激振力的影響較大。

圖10 A 模型傳熱管升力系數(shù)時(shí)程圖Fig.10 T ime history diagram of lift coefficient for model A tube

圖11 頻譜圖Fig.11 Spectrum of lift and drag coefficients

圖12 U 形管無量綱位移時(shí)程Fig.12 Time history of dimensionless displacement of U-tube

圖13 U 形管橫向位移頻譜圖Fig.13 Transverse displacement spectrum of U-tube

圖14 不同約束條件下U 形管無量綱位移時(shí)程圖Fig.14 Time-history diagram of dimensionless displacement of U-pipe under different constraints

圖15 不同雷諾數(shù)下渦激振動(dòng)計(jì)算結(jié)果時(shí)程圖Fig.15 Time history diagram of vortex-induced vibration of U-tueb at different Reynolds numbers

5 結(jié)論

本文通過在FLUENT 程序中求解結(jié)構(gòu)振動(dòng)Newmark-β 方程，基于動(dòng)網(wǎng)格模型開發(fā)了適用于壓水堆蒸汽發(fā)生器傳熱管流致振動(dòng)數(shù)值模的雙向流固耦合方法，開展了典型壓水堆蒸汽發(fā)生器U 形傳熱管彎管的流致振動(dòng)數(shù)值模擬。結(jié)論主要如下：

（1）隨著U 形管防振條數(shù)目的增加，U 形管的固有頻率增大，當(dāng)U 形管沒有防振條時(shí)，一階振動(dòng)頻率和二階振動(dòng)頻率差別較大；當(dāng)U形管增加防振條固定時(shí)，一階振動(dòng)頻率和二階振動(dòng)頻率差別較小，且均出現(xiàn)兩兩相鄰階數(shù)頻率相近的現(xiàn)象；傳熱管的流向振動(dòng)會對橫向振動(dòng)產(chǎn)生影響，因此在計(jì)算過程中應(yīng)該求解雙自由度振動(dòng)方程。

（2）在傳熱管渦激振動(dòng)中，傳熱管在流向方向上產(chǎn)生位移，且在一定位移處來回振動(dòng)，而在橫向方向上圍繞中心點(diǎn)上下振動(dòng)，位移均值為零；傳熱管渦激振動(dòng)對流體的擾動(dòng)不大，單傳熱管漩渦脫落頻率同時(shí)受到自身固有頻率和激振力頻率頻率的影響，阻力系數(shù)傅里葉變化得到的主頻率是升力系數(shù)傅里葉變換得到的主頻率的2 倍。

（3）傳熱管在渦激振動(dòng)影響下，是否安裝防振條對傳熱管的橫向位移影響較大，B 型約束傳熱管橫向振動(dòng)位移為A 型約束傳熱管橫向振動(dòng)位移的1/100；隨著雷諾數(shù)的增加，阻力系數(shù)和升力系數(shù)的周期都逐漸變小，但阻力系數(shù)振動(dòng)的幅值越來越小，升力系數(shù)振動(dòng)幅值基本不變。隨著雷諾數(shù)的增加，管的流向位移、橫向振動(dòng)幅值都逐漸增大，傳熱管的振動(dòng)現(xiàn)象受流體湍流激振力的影響較大。