失水事故引起反應(yīng)堆堆內(nèi)結(jié)構(gòu)動響應(yīng)分析的聲學(xué)有限元方法

趙燮霖 馬英超 葉獻(xiàn)輝 姜乃斌 周進(jìn)雄

1（西安交通大學(xué)航天航空學(xué)院機(jī)械振動與強(qiáng)度國家重點(diǎn)實(shí)驗(yàn)室 西安 710049）

2（中國核動力研究設(shè)計院核反應(yīng)堆系統(tǒng)設(shè)計技術(shù)重點(diǎn)實(shí)驗(yàn)室 成都 610200）

大破口失水事故（Large-break Loss-of-coolant Accident，LB-LOCA）是核反應(yīng)堆結(jié)構(gòu)設(shè)計中必須要考慮的一種極限情況。失水事故引起從破口到反應(yīng)堆內(nèi)部泄壓波的傳播，引起吊籃等結(jié)構(gòu)的流固耦合振動。導(dǎo)致連接件松動、結(jié)構(gòu)動應(yīng)力過大而破壞，因此評估堆內(nèi)結(jié)構(gòu)由于失水事故引起的動響應(yīng)具有重要意義。

在核反應(yīng)堆結(jié)構(gòu)力學(xué)領(lǐng)域?qū)Υ笃瓶谑鹿室鸬亩褍?nèi)結(jié)構(gòu)動響應(yīng)分析始于20世紀(jì)60～70年代，特別是1979年美國三哩島核電站失水事故引起了人們對這一問題的進(jìn)一步關(guān)注。早期開展的研究如Au-Yang等［1-2］認(rèn)為流體的影響僅表現(xiàn)為一種附加水動力質(zhì)量，給出了一種純結(jié)構(gòu)力學(xué)的求解方程，由于忽略了流固耦合作用，此種求解方法的精度不高。同時期也有不少研究機(jī)構(gòu)，用內(nèi)部開發(fā)的流體力學(xué)和結(jié)構(gòu)有限元分析軟件，進(jìn)行單向或雙向流固耦合分析計算，但由于軟件缺乏通用性，分析結(jié)構(gòu)之間差異較大且高度依賴于研究機(jī)構(gòu)的工程設(shè)計經(jīng)驗(yàn)。近年來，隨著計算機(jī)軟件硬件水平的飛速發(fā)展，特別是大型通用CAE（Computer Aided Engineering）計算機(jī)輔助工程軟件在核工程領(lǐng)域的推廣和應(yīng)用，研究者開展了基于計算流體動力學(xué)（Computation Fluid Dynamics，CFD）的單向/雙向流固耦合LOCA事故分析。Hermansky等［3］采用MSC.Dytran有限元軟件對 LOCA 動態(tài)響應(yīng)進(jìn)行數(shù)值模擬；Nilsson［4］利用ADINA采用三種不同的方法對HDR V31.1實(shí)驗(yàn)進(jìn)行數(shù)值模擬，并對比了計算成本；Sheykhi等［5-6］編制相應(yīng)計算程序并結(jié)合ABAQUS對于不同反應(yīng)堆壓力容器在LOCA事故下的壓力極值與應(yīng)力分布進(jìn)行了預(yù)測；Brandt等［7］采用管理MpCCI軟件管理流固數(shù)據(jù)傳遞，基于Fluent和Abaqus對HDR（Heiss Dampf Reaktor）實(shí)驗(yàn)LOCA實(shí)驗(yàn)進(jìn)行了數(shù)值模擬。這些基于CFD的LOCA動響應(yīng)分析，可以使流體力計算更加精確，依賴商用軟件的耦合功能或第三方數(shù)據(jù)通訊交換軟件的功能，使LOCA動響應(yīng)的分析模型從早期的梁模型發(fā)展到三維復(fù)雜堆內(nèi)結(jié)構(gòu)建模和仿真，大大提高了處理復(fù)雜工程實(shí)際問題的能力，但同時也顯著增加了計算成本和設(shè)計分析周期。

在LOCA動響應(yīng)分析方法中，與基于CFD的分析方法并行的另一種方法是基于勢流體的聲學(xué)有限元方法。相比于精確的CFD計算，勢流體聲學(xué)方法忽略了流體的粘性，無須求解速度場，流體動量方程退化為簡單的聲波方程。Bathe等［8-10］提出了基于位移的聲學(xué)有限元方法求解流固耦合問題，基本變量為流體和固體的位移和流體壓力，非常適合用于有限元法求解，使得流固耦合數(shù)值計算得到簡化。此外 Luke［11］與 K?hk?nen［12］以 及 Sussman 與Sundqvist［13］也開展了類似的基于聲學(xué)有限元方法的LOCA動響應(yīng)研究。上述這些研究大都基于用戶自己編制的有限元程序，適用性和通用性有一定的局限性。本文基于ADINA軟件開展了基于聲學(xué)有限元的LOCA動響應(yīng)分析，通過與德國HDR實(shí)驗(yàn)反應(yīng)堆V32破口失水事故實(shí)驗(yàn)結(jié)果和基于ADINA軟件的CFD雙向流固耦合結(jié)果對比，表明聲學(xué)有限元方法用于模擬德國HDR實(shí)驗(yàn)反應(yīng)堆V32破口失水事故的動響應(yīng)是可行的。

1 基于位移的聲學(xué)有限元理論

一般的流體力學(xué)理論其基本的變量為速度、密度、壓力和溫度。采用有限差分法或有限體積法求解。而ADINA采用的是Bathe等［8-10］提出的基于位移的勢流體理論，采用有限元法求解，這樣就和以位移為基本變量的固體有限元方法實(shí)現(xiàn)求解策略的無縫銜接，這也是ADINA流固耦合分析技術(shù)的優(yōu)勢和特色之一。參考文獻(xiàn)［8-10］，基于位移的聲學(xué)有限元理論概述如下：

假設(shè)流體無粘、可壓縮、無旋，聲壓及位移、速度都是變量，則流體的控制方程如下：

其中：式（1）為連續(xù)性方程，式（2）為勢流體忽略重力下的動量方程，也稱為Euler方程。式中：ρ為流體密度；V為流體質(zhì)點(diǎn)速度；p為流體質(zhì)點(diǎn)的壓力。

式（1）、（2）為傳統(tǒng)的以速度、密度、壓力為基本未知量的流體力學(xué)方程。針對勢流體聲學(xué)有限元求解策略，可以將上述方程轉(zhuǎn)化為以位移和壓力為基本變量的方程。定義，其中：c為勢流體的聲速；U為流體質(zhì)點(diǎn)的位移，則連續(xù)性方程（1）兩邊同時對時間積分得到：

而動量方程（2）則簡單改寫為：

式（3）和（4）即為以位移、壓力為基本未知量的勢流體有限元基本方程。固體的平衡方程為：

式中：Us為結(jié)構(gòu)的位移；F為載荷矢量，包括流固耦合作用力；σ為結(jié)構(gòu)的應(yīng)力。式（3）～（5）并加上固體的本構(gòu)方程，就構(gòu)成了以流體/固體位移、流體壓力為基本未知量的流固耦合方程，可以很方便地采用有限元法進(jìn)行離散和求解，具體細(xì)節(jié)可參考Bathe等的文獻(xiàn)［14-16］。

2 計算模型與結(jié)果分析討論

圖1給出了德國HDR實(shí)驗(yàn)反應(yīng)堆的結(jié)構(gòu)圖。本文主要針對V32實(shí)驗(yàn)［17-18］，破口位于降壓管口位置，反應(yīng)堆其他管口在實(shí)驗(yàn)中都處于封閉狀態(tài)。HDR實(shí)驗(yàn)堆主要考慮失水事故引起的堆內(nèi)吊籃的動響應(yīng)，其余結(jié)構(gòu)都作了簡化處理，堆內(nèi)結(jié)構(gòu)對動響應(yīng)的影響通過吊籃下方的質(zhì)量環(huán)進(jìn)行模擬，吊籃下端自由，上端剛性連接于壓力容器內(nèi)壁上。

圖1 HDR壓力容器結(jié)構(gòu)圖Fig.1 Structural drawing of pressure vessel of HDR

圖2 給出了計算所采用的流體及固體域網(wǎng)格，網(wǎng)格由Hypermesh生成，經(jīng)過網(wǎng)格無關(guān)解檢驗(yàn)，最終確定流體域網(wǎng)格6.4萬，固體域網(wǎng)格4.5萬。固體域采用實(shí)體單元，流體域采用3D亞聲速流體單元。流固界面處網(wǎng)格共節(jié)點(diǎn)。破口壓力邊界條件施加簡化了HDR實(shí)驗(yàn)中傳感器實(shí)測的壓力值，初始時刻冷卻劑 壓 力 11 MPa，在 1 ms內(nèi) 驟 降 8.75 MPa，在1～1.5 ms內(nèi)壓力回升0.75 MPa，在1.5 ms后破口壓力維持在3 MPa。

圖2 HDR流體區(qū)域(a)和固體區(qū)域(b)模型網(wǎng)格Fig.2 HDR meshes of fluid region(a)and solid region(b)

聲學(xué)有限元計算所用參數(shù)如下：流體密度ρ0=730 kg·m-3，體積壓縮模量k=8.202×108Pa；固體材料密度ρ=7 800 kg·m-3，線彈性材料楊氏模量E=1.75×1011Pa，泊松比v=0.3。聲學(xué)有限元采用ADINA9.4中的ADINAstructures模塊，求解器采用隱式動力學(xué)求解，Newmark積分方法。為了與聲學(xué)有限元計算結(jié)構(gòu)對比，同時開展了基于CFD的雙向流固耦合模擬，流體采用可壓縮層流模型，流體動力粘度系數(shù)μ=0.000 01 kg·m-3，流固交界面定義流固耦合邊界條件。

圖3給出了不同時刻下流體壓力分布云圖，可以看出破口失水事故引起流體域內(nèi)壓力波的傳播過程。圖4給出了由于波傳播導(dǎo)致的反應(yīng)堆結(jié)構(gòu)位移云圖，可以看出由于泄壓波在流域內(nèi)的傳播，對浸沒在流體域的吊籃施加交變的壓力，導(dǎo)致吊籃產(chǎn)生搖擺振動。

圖3 不同時刻下流體區(qū)域壓力分布云圖Fig.3 Plots of pressure contours at selected instants of time

圖4 不同時刻下固體區(qū)域位移分布云圖Fig.4 Plots of displacement contours at selected instants of time

為了進(jìn)一步驗(yàn)證所采用的聲學(xué)有限元計算結(jié)果，圖5給出了破口下方1.7 m高度處吊籃與壓力容器之間的徑向位移差，圖6給出了破口高度處（8.85 m）吊籃外壁的軸向應(yīng)變。可以看出，基于聲學(xué)有限元方法和基于CFD雙向流固耦合計算方法所預(yù)測的徑向位移差與軸向應(yīng)變較實(shí)驗(yàn)測試結(jié)果吻合良好，相比于聲學(xué)有限元方法基于CFD雙向流固耦合計算方法給出的曲線更加光滑，這與已有文獻(xiàn)中的結(jié)果一致［13］。相比于CFD方法，聲學(xué)有限元方法計算效率較高，采用配置了4核i7CPU 3.6 Hz處理器的工作站進(jìn)行并行計算，聲學(xué)有限元方法計算耗時3 h 46 min，CFD流固耦合方法耗時5 h 16 min。

圖5 破口下方1.7 m處吊籃與RPV相對徑向位移之差Fig.5 Relative radial displacement between the core barrel and RPV under blowdown nozzle 1.7 m

圖6 吊籃外壁在破口高度（8.85 m）處軸向應(yīng)變Fig.6 Axial strain on the outer surface of the core barrel at height of blowdown nozzle(8.85 m)

3 結(jié)語

基于有限元軟件ADINA中的勢流流體模塊，對德國HDR實(shí)驗(yàn)反應(yīng)堆V32破口失水事故實(shí)驗(yàn)進(jìn)行了聲和結(jié)構(gòu)雙向耦合數(shù)值模擬，結(jié)果表明：聲學(xué)有限元模型計算結(jié)果與CFD雙向流固耦合以及實(shí)驗(yàn)結(jié)果吻合良好。由于聲學(xué)有限元計算方法只需要求解聲波方程而無需求解復(fù)雜的NS方程，因此較CFD流固耦合方法的計算效率略高。反應(yīng)堆破口失水事故過程中冷卻劑的壓力波包含由破口壓降導(dǎo)致的入射波和由于固體位移引起的誘發(fā)波，誘發(fā)波滿足小振幅假設(shè)，因此這類短時瞬態(tài)問題適合用聲學(xué)有限元求解。但勢流模型不包含流體粘性，因而用于長時程動響應(yīng)分析時要慎重。