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

(中國(guó)核動(dòng)力研究設(shè)計(jì)院，四川 成都 610213)

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

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

1 超單元分析模型的建立

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

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

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

1.2 超單元模型建立

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

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

1.3 模型驗(yàn)證

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

圖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è)計(jì)。

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

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

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

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

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

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

(1)

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

格架水間隙示意圖如圖5所示。對(duì)于圖5所示的兩個(gè)水下矩形柱體，中間矩形柱體代表格架，外圍矩形柱體代表周?chē)乇诨蚋窦?，中間間隙充滿(mǎn)水。兩柱體在水平方向上相對(duì)運(yùn)動(dòng)而產(chǎn)生的水動(dòng)力質(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)稱(chēng)尺寸；g1、g2、g3、g4為格架周?chē)g隙尺寸。各格架外部的附加質(zhì)量通過(guò)APDL編程后實(shí)現(xiàn)自動(dòng)計(jì)算，并自動(dòng)賦值給流固耦合單元FLUID38實(shí)現(xiàn)附加質(zhì)量的模擬。乏燃料組件內(nèi)含水以及與貯存單元之間的水間隙附加質(zhì)量，使用MASS21質(zhì)量點(diǎn)分別施加在組件和貯存單元的水平方向上[14]。

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

圖6 地震加速度時(shí)程曲線(xiàn)Fig.6 Seismic acceleration time history curve

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

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

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

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

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

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

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

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

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

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

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

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

6臺(tái)格架滿(mǎn)載、半載、空載時(shí)在OBE、SSE中底座對(duì)水池底的最大應(yīng)力如圖12所示，由于格架重心向-x、-y方向傾斜，格架左上角的底座對(duì)水池的應(yīng)力強(qiáng)度較大，格架1在半載、SSE情況下與水池發(fā)生最大相互作用，發(fā)生時(shí)刻為9.79 s。對(duì)格架1超單元模型進(jìn)行擴(kuò)展，得到格架1最大應(yīng)力發(fā)生在底座。底座結(jié)構(gòu)如圖13所示，格架小幅度傾斜時(shí)，底墊板不會(huì)傳遞豎直方向的扭矩，螺柱底部受到的作用力由螺柱傳遞給螺套，豎直方向的力由螺套肩部傳遞給底板，水平方向的力由螺套頸部傳遞給底板。格架底座最大壓應(yīng)力發(fā)生在螺柱與螺套連接處，如圖14所示，最大值為230 MPa，滿(mǎn)足應(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方向的最大滑移位移曲線(xiàn)如圖15所示。在x方向，格架空載情況下在SSE工況時(shí)產(chǎn)生的滑移距離最大，為-15.4 mm；同理，y方向最大滑移距離發(fā)生在格架空載且處于SSE工況時(shí)，為-27.4 mm，均小于格架之間、格架與水池之間的間隙，不會(huì)發(fā)生碰撞。在OBE、SSE兩種地震工況下均顯示格架空載時(shí)對(duì)地震的響應(yīng)最劇烈，在同一時(shí)刻等加速度條件下，滑移幅度從小到大依次為滿(mǎn)載、半載、空載，即隨著組件裝載量的增加，格架整體質(zhì)量增加，格架滑移幅度有減小的趨勢(shì)，因此在格架滑移分析時(shí)空載情況應(yīng)作為重點(diǎn)工況進(jìn)行驗(yàn)證。

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

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

3.3 格架傾倒分析

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

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

(6)

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

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

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

4 結(jié)論

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