基于參數(shù)化-超單元法的HFETR乏燃料貯存格架抗震分析

劉金龍，章航洲，王 帥，孫志軍，隆 濤

(中國核動力研究設(shè)計院，四川 成都 610213)

乏燃料貯存格架(簡稱格架)是用于貯存乏燃料組件的關(guān)鍵設(shè)備，在地震中可能出現(xiàn)變形、碰撞、傾倒等現(xiàn)象，影響乏燃料組件的貯存安全。目前，格架抗震分析一般采用數(shù)值模擬，常用的分析模型有梁模型和殼模型。梁模型計算速度快，但精度較低[1-3]；殼模型精度較高，但計算時間長，占用資源多，很難實(shí)現(xiàn)多格架同時分析[4-6]。另外，格架結(jié)構(gòu)在滿足抗震安全的同時還須符合臨界、熱工等要求，設(shè)計中可能需要多次調(diào)整結(jié)構(gòu)尺寸參數(shù)、重復(fù)建模和抗震計算，而傳統(tǒng)的CAD設(shè)計、CAE分析建模效率較低，導(dǎo)致格架迭代設(shè)計周期大幅增加。

本文在對格架抗震分析模型及建模方式調(diào)研和分析的基礎(chǔ)上，提出基于超單元法建立高通量工程試驗堆(HFETR)格架分析模型，解決模型精度和計算效率不能兼顧的問題[7]；然后基于ANSYS APDL程序?qū)崿F(xiàn)格架分析模型的參數(shù)化，通過VB平臺對格架程序進(jìn)行封裝和可視化開發(fā)，建立格架分析界面，解決格架建模效率低的問題；最后，利用該方法對整池HFETR格架開展抗震分析。

1 超單元分析模型的建立

HFETR乏燃料組件貯存格架主要部件有貯存單元、柵板、圍板、加固梁、底板、底座等，最大外形尺寸為2.4 m×1.8 m×5.2 m(長×寬×高)，格架結(jié)構(gòu)如圖1a所示。貯存單元為3×5排列，內(nèi)部存放乏燃料組件，底部與底板焊接，中部和上部與柵板孔間隙配合。圍板與加固梁、底板、柵板之間通過焊接連接。共設(shè)計6臺乏燃料貯存格架，每臺格架單獨(dú)放置在水池中，相互間不連接，格架與池底間不固接。

1.1 單元選擇及網(wǎng)格劃分

格架貯存單元、圍板、柵板、底板以及底座底板采用殼單元SHELL181模擬；格架外圍加固梁、底座與凸臺之間的支撐柱以及組件采用梁單元BEAM188模擬；格架與底板接觸部位使用實(shí)體單元SOLID185模擬；格架圍板與底板、柵板與圍板、貯存單元與底板、凸臺與底板設(shè)為綁定接觸，接觸單元使用CONTA175、TARGE170，為避免接觸單元影響超單元生成，設(shè)KEYOPT(2)=1(罰函數(shù)法)[8-9]。網(wǎng)格劃分后的有限元模型(稱普通有限元模型)如圖1b所示。

1.2 超單元模型建立

根據(jù)格架結(jié)構(gòu)和接觸邊界情況對普通有限元模型進(jìn)行子結(jié)構(gòu)劃分和主節(jié)點(diǎn)選取，將底座上半部分及以上的底板、圍板、柵板、貯存單元作為超單元，底座部分支撐柱和底墊板部分保留為非超單元，格架外部輪廓節(jié)點(diǎn)作為主節(jié)點(diǎn)，生成的超單元模型如圖1c所示。

a——格架幾何模型；b——普通有限元模型；c——超單元模型

1.3 模型驗證

分別對超單元模型和普通有限元模型底座底面施加固定約束并進(jìn)行模態(tài)分析，兩種模型的模態(tài)頻率列于表1，振型示于圖2、3。由表1可見，超單元模型與普通有限元模型的模態(tài)頻率相差較小，最大相對偏差僅1.122%。由圖2、3可見，兩種模型的振型基本一致。表明本文所建超單元模型精度滿足要求，可用于格架抗震分析。

圖2 普通有限元模型前3階振型Fig.2 First three mode shapes of ordinary finite element model

表1 兩種模型的模態(tài)頻率比較Table 1 Comparison of modal frequencies between two models

2 基于APDL和VB的格架模型參數(shù)化

格架模型參數(shù)化主要包括幾何模型參數(shù)化、邊界條件參數(shù)化以及可視化分析界面設(shè)計。

2.1 幾何模型參數(shù)化

本文基于APDL程序建立格架幾何尺寸變量和各部件間尺寸關(guān)系，實(shí)現(xiàn)格架從內(nèi)到外、從上到下的尺寸驅(qū)動建模，如圖4所示[10-11]。為使各部件之間保持關(guān)聯(lián)，利用APDL關(guān)鍵點(diǎn)命令建立關(guān)鍵點(diǎn)坐標(biāo)，通過點(diǎn)驅(qū)動方式實(shí)現(xiàn)格架部件之間尺寸的驅(qū)動，部件連接部位共用關(guān)鍵點(diǎn)，在模型修改時可避免部件之間出現(xiàn)交叉或分離等情況，保證了格架整體結(jié)構(gòu)的連續(xù)性。組件的陣列主要借助APDL的AGEN、VGEN、LGEN等命令實(shí)現(xiàn)，通過輸入組件排列方式，自動陣列貯存單元；多格架的排列基于超單元模型，通過SETRN、SE、*DO等命令可實(shí)現(xiàn)超單元模型提取、擺放距離調(diào)節(jié)和擺放組合變化等。最終通過程序設(shè)計實(shí)現(xiàn)格架模型的建立、尺寸修改和驅(qū)動、格架內(nèi)組件裝載和多格架擺放方式的變化等，為格架結(jié)構(gòu)及布置迭代設(shè)計提供有力支持。

圖3 超單元模型前3階振型Fig.3 First three mode shapes of superelement model

圖4 格架各部件驅(qū)動示意圖Fig.4 Diagram of rack components drive

2.2 邊界條件參數(shù)化

格架抗震分析時所需考慮的邊界條件主要包括水環(huán)境、地震、阻尼、摩擦等。其中水環(huán)境的仿真最為復(fù)雜和困難。在地震中，格架周圍的水對格架運(yùn)動的影響十分顯著，水對格架的影響不能忽略。目前，模擬格架與水、組件與水之間流固耦合作用常用的方法是附加質(zhì)量法，該方法適用于小位移附加質(zhì)量的計算，水動力質(zhì)量矩陣形式[12-13]如下：

(1)

式中：MH為流固耦合水動力質(zhì)量，本格架x、y方向為水平方向，分別為MHx和MHy；M1為格架所處位置置換水的質(zhì)量；M2為邊界內(nèi)無格架情況下水的質(zhì)量。

格架水間隙示意圖如圖5所示。對于圖5所示的兩個水下矩形柱體，中間矩形柱體代表格架，外圍矩形柱體代表周圍池壁或格架，中間間隙充滿水。兩柱體在水平方向上相對運(yùn)動而產(chǎn)生的水動力質(zhì)量在數(shù)值上表示為：

圖5 格架水間隙示意圖Fig.5 Schematic diagram of water gap of rack

MHx=2ρhC2(C/3g1+C/3g3+B/g2+B/g4)

(2)

MHy=2ρhB2(B/3g2+B/3g4+C/g1+C/g3)

(3)

(4)

(5)

式中：ρ為水質(zhì)量密度；h為格架高度；C、B為考慮水間隙的標(biāo)稱尺寸；g1、g2、g3、g4為格架周圍間隙尺寸。各格架外部的附加質(zhì)量通過APDL編程后實(shí)現(xiàn)自動計算，并自動賦值給流固耦合單元FLUID38實(shí)現(xiàn)附加質(zhì)量的模擬。乏燃料組件內(nèi)含水以及與貯存單元之間的水間隙附加質(zhì)量，使用MASS21質(zhì)量點(diǎn)分別施加在組件和貯存單元的水平方向上[14]。

地震作用采用加速度時程模擬[15-16]，考慮運(yùn)行基準(zhǔn)地震(OBE)和安全停堆地震(SSE)兩種地震工況，地震加速度時程曲線如圖6所示，分別在x、y、z3個方向上同時施加地震載荷，持續(xù)時間20 s。地震數(shù)據(jù)通過*CREATE、*VREAD、ACEL等命令實(shí)現(xiàn)自動載入和施加。

圖6 地震加速度時程曲線Fig.6 Seismic acceleration time history curve

格架自由放置在乏燃料水池中，為重點(diǎn)關(guān)注格架滑移情況，格架底座與水池底面的摩擦系數(shù)保守設(shè)為0.2[17]。結(jié)構(gòu)阻尼采用瑞利阻尼，OBE時阻尼比取2%，SSE時阻尼比取4%，根據(jù)模態(tài)分析結(jié)果選擇感興趣的2.2～30.1 Hz進(jìn)行計算[5]。OBE時，α=0.515，β=1.97×10-4；SSE時，α=1.03，β=3.94×10-4。單格架模型及邊界條件加載情況如圖7所示，格架底座與水池底面采用面面接觸(CONTA175、TARGE170)，抗震計算采用完全法瞬態(tài)動力分析(ANTYPE、TRANS)。

圖7 單格架模型環(huán)境設(shè)置示意圖Fig.7 Schematic diagram of single frame model environment setting

2.3 可視化分析界面設(shè)計

格架模型參數(shù)化程序編寫完成后，利用VB對APDL程序封裝并建立參數(shù)關(guān)聯(lián)，通過開發(fā)可視化界面，使格架參數(shù)修改更加便捷，建模和分析時通過調(diào)用ANSYS軟件(Batch功能)在后臺運(yùn)行，將命令流提交給ANSYS進(jìn)行分析，結(jié)束后自動關(guān)閉ANSYS，最后使用VB命令提取計算數(shù)據(jù)，實(shí)現(xiàn)格架分析結(jié)果的查看，其流程及原理如圖8所示。

圖8 VB調(diào)用ANSYS原理圖Fig.8 VB calling ANSYS schematic diagram

基于VB可視化平臺開發(fā)的界面可實(shí)現(xiàn)格架模型建立和修改、邊界條件加載、分析計算以及結(jié)果查看等整個流程的操作，主界面如圖9所示。

圖9 格架分析主界面Fig.9 Main interface of rack analysis

3 乏燃料貯存格架抗震分析

本文分別選取滿載、半載、空載3種典型裝載形式對整池6臺格架開展抗震分析，格架裝載情況如圖10所示。格架之間按設(shè)計間距擺放，6臺格架超單元模型如圖11所示。根據(jù)計算結(jié)果對格架應(yīng)力、碰撞、傾倒情況進(jìn)行分析。

圖10 單格架裝載示意圖Fig.10 Schematic diagram of single rack loading

圖11 6臺格架超單元模型示意圖Fig.11 Schematic diagram of superelement model of six racks

3.1 格架應(yīng)力分析

6臺格架滿載、半載、空載時在OBE、SSE中底座對水池底的最大應(yīng)力如圖12所示，由于格架重心向-x、-y方向傾斜，格架左上角的底座對水池的應(yīng)力強(qiáng)度較大，格架1在半載、SSE情況下與水池發(fā)生最大相互作用，發(fā)生時刻為9.79 s。對格架1超單元模型進(jìn)行擴(kuò)展，得到格架1最大應(yīng)力發(fā)生在底座。底座結(jié)構(gòu)如圖13所示，格架小幅度傾斜時，底墊板不會傳遞豎直方向的扭矩，螺柱底部受到的作用力由螺柱傳遞給螺套，豎直方向的力由螺套肩部傳遞給底板，水平方向的力由螺套頸部傳遞給底板。格架底座最大壓應(yīng)力發(fā)生在螺柱與螺套連接處，如圖14所示，最大值為230 MPa，滿足應(yīng)力要求。

圖12 水池底應(yīng)力云圖Fig.12 Stress intensity cloud image at bottom of pool

圖13 底座示意圖Fig.13 Schematic diagram of base

圖14 底座應(yīng)力強(qiáng)度最大部位云圖Fig.14 Maximum stress intensity cloud image of support feet

3.2 格架滑移分析

格架在OBE、SSE工況下x、y方向的最大滑移位移曲線如圖15所示。在x方向，格架空載情況下在SSE工況時產(chǎn)生的滑移距離最大，為-15.4 mm；同理，y方向最大滑移距離發(fā)生在格架空載且處于SSE工況時，為-27.4 mm，均小于格架之間、格架與水池之間的間隙，不會發(fā)生碰撞。在OBE、SSE兩種地震工況下均顯示格架空載時對地震的響應(yīng)最劇烈，在同一時刻等加速度條件下，滑移幅度從小到大依次為滿載、半載、空載，即隨著組件裝載量的增加，格架整體質(zhì)量增加，格架滑移幅度有減小的趨勢，因此在格架滑移分析時空載情況應(yīng)作為重點(diǎn)工況進(jìn)行驗證。

圖15 地震工況下格架在x、y方向的最大滑移位移曲線Fig.15 Maximum slip displacement curves of rack in x and y directions under seismic condition

另外，格架在半載情況下累積的滑移位移與滿載、半載情況下有方向相反的趨勢，這是由于受格架結(jié)構(gòu)及裝載情況影響，質(zhì)量大小及分布不同，格架受地震影響產(chǎn)生的傾斜姿態(tài)、底座摩擦力分布等情況各異，使得格架滑移位移累積效果各不相同。

3.3 格架傾倒分析

為保守估計格架傾倒情況，將抗震計算得到的格架底部最大起跳高度與格架傾倒限值比較。格架傾斜示意圖如圖16所示。由圖16可知，格架傾斜時重心抬高，外側(cè)底座為主要支撐點(diǎn)，當(dāng)重心位于支撐點(diǎn)正上方時，重心抬升至最高值，超過該點(diǎn)即會發(fā)生傾倒。格架傾斜時重心抬高距離Δh按下式計算：

圖16 格架傾斜示意圖Fig.16 Diagram of rack tilt

(6)

式中：h為格架重心高度；l為支撐點(diǎn)與重心的垂直距離。格架重心抬高Δh時，支撐點(diǎn)另一側(cè)底座將抬高Δz(該高度視為格架傾倒限值)，l越小傾倒限值越小。根據(jù)以上分析，分別計算得到格架滿載、半載和空載時最小的傾倒限值分別為531、472、557 mm。

地震工況下格架在z方向的最大起跳高度曲線示于圖17?？梢姡贠BE、SSE工況中，格架均在半載情況下獲得最大起跳高度，格架半載時OBE工況下最大起跳高度為0.907 mm，SSE工況下最大起跳高度為1.21 mm，兩者均遠(yuǎn)小于格架傾倒限值，因此格架在OBE和SSE工況下均不會發(fā)生傾倒。分析表明，格架在半載時傾倒限值最小，而在地震中格架半載時傾斜幅度最大，因而在地震中格架半載時發(fā)生傾倒的可能性最大，在格架傾倒分析中半載情況應(yīng)作為重點(diǎn)工況進(jìn)行驗證。

圖17 地震工況下格架在z方向的最大起跳高度曲線Fig.17 Maximum jump height curve of rack in z direction under seismic condition

4 結(jié)論

本研究基于超單元法建立了HFETR乏燃料貯存格架抗震分析模型，并基于APDL程序和VB平臺對格架分析模型進(jìn)行了參數(shù)化處理，建立了格架專用分析界面，解決了多格架抗震分析困難以及結(jié)構(gòu)迭代設(shè)計流程復(fù)雜的問題，為格架結(jié)構(gòu)的迭代設(shè)計及抗震分析提供了有力的支撐。利用參數(shù)化-超單元法對整池HFETR格架進(jìn)行了抗震分析，對比分析了格架在空載、半載和滿載3種情況下的表現(xiàn)。結(jié)果表明：格架在3種情況下應(yīng)力均滿足要求，不會發(fā)生碰撞和傾倒；在滑移分析和傾倒分析中，空載和半載工況需重點(diǎn)關(guān)注。