快堆堆芯單排組件抗震計(jì)算關(guān)鍵算法研究與驗(yàn)證

周世安,黑寶平,高付海

(中國(guó)原子能科學(xué)研究院,北京 102413)

對(duì)于反應(yīng)堆及核設(shè)施,在強(qiáng)烈的外部載荷如地震作用下確保安全運(yùn)行是最基本的要求之一。液態(tài)金屬快堆(LMFR)中包含低壓運(yùn)行的系統(tǒng)以及薄壁柔性組件,這使其在地震作用下受到的影響極其顯著。對(duì)于快堆堆芯,其地震下的運(yùn)動(dòng)主體十分龐大,全堆芯的六角形套管組件多達(dá)成百上千根,每根組件與相鄰的6根組件均可能發(fā)生碰撞,且每時(shí)每刻的運(yùn)動(dòng)狀態(tài)都在發(fā)生變化,從而導(dǎo)致碰撞的部位也隨之變化,同時(shí)會(huì)受到液體金屬冷卻劑與組件流固耦合的影響,這使得堆芯抗震問題變得尤為復(fù)雜,因此建立針對(duì)快堆堆芯的抗震分析方法很有必要。

快堆堆芯的抗震分析方法通常分為設(shè)計(jì)試驗(yàn)和程序分析,試驗(yàn)數(shù)據(jù)直觀可靠,但開展全堆芯組件抗震試驗(yàn)成本高昂,相比之下,研發(fā)專用程序以準(zhǔn)確高效地進(jìn)行全堆芯抗震計(jì)算和評(píng)價(jià)是一種經(jīng)濟(jì)可行的方式。國(guó)外許多國(guó)家已研發(fā)出了成熟的快堆堆芯抗震分析軟件[1-5],有的軟件是專用于堆芯抗震的計(jì)算程序,如日本研發(fā)的FINAS及REVIAN-3D、印度研發(fā)的CORE-SEIS、韓國(guó)研發(fā)的SAC-CORE和意大利研發(fā)的CORALIE及CLASH;有的則是在原有的通用成熟軟件中進(jìn)行了更新,使其具備快堆堆芯抗震計(jì)算能力,如法國(guó)的CASTEM,美國(guó)的COSMOS,俄羅斯的DINARA以及日本的FINDS、SAFA和SALCON。

程序設(shè)計(jì)采用的計(jì)算方法通常分為2種:模態(tài)疊加法和直接積分法。模態(tài)疊加法可提升計(jì)算速度,但忽略了高階模態(tài)的影響;直接積分法則計(jì)算速度較慢。使用模態(tài)疊加法的程序有COSMOS、FINDS和SAFA等,SAFA還配置了專用于單自由度系統(tǒng)的Nigam時(shí)間積分法,更大程度地提高了計(jì)算效率。使用直接積分法的程序有CORE-SEIS、SACCORE、FINAS和SALCON等,SALCON在直接積分法的基礎(chǔ)上采用GUYAN自由度凝聚方法減少矩陣維度,并通過(guò)變時(shí)間步長(zhǎng)進(jìn)一步提升計(jì)算效率。程序中堆芯組件管腳約束方式通常有3種:固支約束+等效彈簧、簡(jiǎn)支約束+間隙彈簧,以及固支約束+等效旋轉(zhuǎn)彈簧元。大多數(shù)程序都可對(duì)管腳底端施加固支或簡(jiǎn)支約束,并在管腳球形密封處添加等效彈簧以施加間隙等效剛度,而FINDS與FINAS除對(duì)底端施加固支約束或簡(jiǎn)支約束外,還能使用旋轉(zhuǎn)彈簧元來(lái)模擬管腳處的運(yùn)動(dòng)。碰撞動(dòng)力方程的常用求解方法有中心差分法、Newmark法和基于Newmark法衍生的Wilson-θ法,其中Newmark法的參數(shù)γ和β變化時(shí),算法的顯隱式、計(jì)算速度和穩(wěn)定性條件隨之改變。本程序取γ=1/2、β=1/4,使Newmark法為平均加速度法,這種參數(shù)取法使算法為隱式且無(wú)條件穩(wěn)定,通過(guò)這種取法保證求解的穩(wěn)定性。冷卻劑與堆芯組件的流固耦合計(jì)算方式有3種:集中附加質(zhì)量法、均質(zhì)化方法[6]和有限元計(jì)算法[5,7]。大部分程序常用集中附加質(zhì)量法,該方法可調(diào)整組件的1階固有頻率以體現(xiàn)冷卻劑對(duì)組件的附加作用。

中國(guó)原子能科學(xué)研究院曾使用ANSYS、FINAS以及CASTEM對(duì)快堆堆芯組件進(jìn)行了初步的單根[8]以及單排組件抗震分析計(jì)算[9],但一直沒有形成自主開發(fā)的堆芯抗震分析專用軟件程序。

針對(duì)我國(guó)在快堆堆芯抗震計(jì)算自研專用程序方面的空白,本文研究建立適用于大規(guī)模快堆堆芯抗震分析的理論模型,研發(fā)能解決地震載荷作用下堆芯組件復(fù)雜非線性多點(diǎn)碰撞問題的時(shí)間歷程關(guān)鍵算法,并從簡(jiǎn)單的懸臂梁振動(dòng)理論解到國(guó)際原子能機(jī)構(gòu)(IAEA)組織的法國(guó)RAPSODIE堆單排19根組件試驗(yàn)結(jié)果進(jìn)行從下到上的逐步驗(yàn)證,完成單排組件情況下地震分析理論的建立與驗(yàn)證,證明算法的正確性和可行性,為后續(xù)多排組件乃至全堆芯抗震分析奠定基礎(chǔ)。

1 理論模型

快堆堆芯組件在地震動(dòng)時(shí)的動(dòng)力學(xué)方程為:

(1)

鑒于組件在地震過(guò)程中主要表現(xiàn)的行為是水平碰撞且在碰撞中基本沒有軸向變形和剪切變形,故將組件簡(jiǎn)化為歐拉梁,并以此為根據(jù)構(gòu)建初始的質(zhì)量矩陣與剛度矩陣。

1.1 質(zhì)量矩陣和剛度矩陣

僅考慮兩端受載的情況,則歐拉梁的梁平衡方程為:

(2)

其中:u為梁的撓度,即梁?jiǎn)卧?jié)點(diǎn)處垂直于軸線方向的位移;E為梁的彈性模量;I為梁的慣性矩。

通過(guò)Hermite插值方程表示的梁形函數(shù)與虛位移原理,可推導(dǎo)出梁的單元?jiǎng)偠染仃嚭蛦卧|(zhì)量矩陣,通過(guò)對(duì)不同參數(shù)梁的單元矩陣進(jìn)行裝配最終獲得總的質(zhì)量矩陣和剛度矩陣。為保證計(jì)算精度足夠高并大幅提升計(jì)算效率,本文采用GUYAN縮減法[10]進(jìn)行相應(yīng)的自由度縮減。首先需要選定主從自由度,其中主自由度是施加了力、計(jì)算邊界條件和需要輸出結(jié)果的所在自由度,其余的則是從自由度。本質(zhì)上,縮減矩陣是使用主自由度來(lái)表示從自由度。對(duì)于包含阻尼和外力作用的系統(tǒng),其動(dòng)力學(xué)方程為:

(3)

其中,{F}為受到的外力向量。

(4)

根據(jù)靜平衡近似,可得到:

(5)

[K][T]{Xb}={F}

(6)

(7)

(8)

1.2 結(jié)構(gòu)阻尼矩陣

對(duì)于地震作用下的堆芯組件,除流固耦合產(chǎn)生的附加阻尼及碰撞時(shí)產(chǎn)生的碰撞阻尼,組件結(jié)構(gòu)本身存在結(jié)構(gòu)阻尼,工程上常用Rayleigh阻尼[11]來(lái)表示。Rayleigh阻尼可表示為:

[C]=a0[M]+a1[K]

(9)

(10)

對(duì)于新型結(jié)構(gòu),通常情況下需通過(guò)試驗(yàn)測(cè)量獲取各階精確阻尼比和頻率?？於讯研窘M件材料一般為不銹鋼,但燃料組件和屏蔽組件內(nèi)部結(jié)構(gòu)材料和工藝設(shè)計(jì)有所差異,因此二者阻尼比不同。由于快堆堆芯組件為細(xì)長(zhǎng)梁結(jié)構(gòu),第1階和第2階振型對(duì)相鄰組件的碰撞行為起主要作用,故在求解堆芯抗震總動(dòng)力方程的結(jié)構(gòu)阻尼矩陣時(shí),分別選取對(duì)應(yīng)上述振型的第1階和第2階固有頻率和阻尼比。在確定ω1和ω2及ξ1和ξ2的基礎(chǔ)上,若ξ1=ξ2=ξ,則a0和a1可表示為:

(11)

1.3 碰撞剛度及阻尼

通常情況下,快堆堆芯組件的碰撞剛度可通過(guò)以下兩種方法獲取:1) 直接法,通過(guò)對(duì)真實(shí)組件進(jìn)行碰撞試驗(yàn)直接獲取程序所需碰撞剛度;2) 間接法,通過(guò)組件壁面的擠壓剛度間接計(jì)算求得程序所需碰撞剛度。本文采用間接法計(jì)算碰撞剛度Kgap,可通過(guò)式(12)求得[1]。

(12)

其中,Ket1和Ket2分別為發(fā)生碰撞的兩根組件的擠壓剛度。

擠壓剛度可通過(guò)對(duì)組件施加如圖1所示的壓力使其發(fā)生形變來(lái)求得。通過(guò)有限元對(duì)組件進(jìn)行建模,兩端施加力F,計(jì)算出ΔD,便可求出組件的擠壓剛度Ket:

(13)

圖1 擠壓剛度計(jì)算示意圖Fig.1 Schematic diagram of extrusion stiffness calculation

當(dāng)發(fā)生碰撞的兩根組件的結(jié)構(gòu)材料和截面形狀一致時(shí),Ket1=Ket2=Ket。根據(jù)碰撞剛度Kgap計(jì)算得到?jīng)_擊阻尼Cgap:

(14)

其中,e為材料的彈性系數(shù),對(duì)于不銹鋼,e一般為0.55。

在單排組件的情況下,組件墊塊處節(jié)點(diǎn)的位移無(wú)需根據(jù)幾何位置關(guān)系進(jìn)行分解,此時(shí)剛度的設(shè)置得到相應(yīng)的簡(jiǎn)化。當(dāng)?shù)趈根組件受到前后兩根組件的沖擊時(shí),總剛度矩陣和阻尼矩陣變?yōu)槿缦滦问?

(15)

aj-1=cj-1+Pgap

(16)

bj-1=-Pgap

(17)

其中:c為原始總剛度矩陣和阻尼矩陣在碰撞組件處的對(duì)應(yīng)項(xiàng);Pgap為碰撞產(chǎn)生的沖擊剛度或沖擊阻尼;a為碰撞對(duì)剛度矩陣和阻尼矩陣在對(duì)角項(xiàng)上產(chǎn)生的影響;b為碰撞對(duì)剛度矩陣和阻尼矩陣在非對(duì)角項(xiàng)上產(chǎn)生的影響。

2 基準(zhǔn)例題驗(yàn)證

2.1 等截面懸臂梁固有頻率驗(yàn)證

對(duì)于簡(jiǎn)單的等截面懸臂梁(圖2),其存在固有頻率理論解。將梁截面參數(shù)輸入到研發(fā)的自研程序內(nèi),計(jì)算固有頻率并與理論解進(jìn)行對(duì)比,以驗(yàn)證梁?jiǎn)卧臉?gòu)建與裝配,以及自由度縮減法則的正確性。

圖2 簡(jiǎn)單等截面懸臂梁示意圖Fig.2 Simple uniform cantilever beam

圖2所示等截面矩形梁,長(zhǎng)為0.02 m、寬為0.01 m、高為1 m。其楊氏模量為2.1×105MPa、密度為7.85×103kg/m3、泊松比為0.3。將該梁分別在ANSYS與自研程序中建模并與理論解進(jìn)行對(duì)比,結(jié)果列于表1。對(duì)比發(fā)現(xiàn):自研程序計(jì)算的前3階固有頻率與理論解和ANSYS的計(jì)算結(jié)果完全吻合,隨著階數(shù)的升高,誤差逐漸變大,在前5階,頻率相對(duì)誤差均小于1%。

2.2 中國(guó)實(shí)驗(yàn)快堆堆芯燃料組件固有頻率驗(yàn)證

真實(shí)組件管腳和柵板插座為間隙配合裝配。對(duì)中國(guó)實(shí)驗(yàn)快堆(CEFR)堆芯燃料組件進(jìn)行質(zhì)量等效與剛度等效,等效參數(shù)列于表2。

表1 簡(jiǎn)單等截面懸臂梁固有頻率計(jì)算結(jié)果對(duì)比Table 1 Comparison between theoretical and calculation results of nature frequency in simple uniform cantilever beam example

表2 CEFR堆芯燃料組件等效后的梁參數(shù)Table 2 CEFR core fuel assembly equivalent beam parameter

組件可等效為變截面梁,在管腳底端設(shè)置為簡(jiǎn)支約束,且在管腳球形密封處設(shè)置剛度為8.8×105N·m的等效彈簧[12],將模型參數(shù)輸入自研程序和ANSYS中計(jì)算固有頻率并進(jìn)行對(duì)比,結(jié)果列于表3。從表3可知,二者吻合較好,最大誤差在第3階處,僅為0.49%。2.1和2.2節(jié)中梁的約束條件不同、截面參數(shù)不同,因此梁?jiǎn)卧臉?gòu)建和矩陣自由度縮減方式均不同,通過(guò)這兩個(gè)案例可驗(yàn)證組件等效方法的合理性,同時(shí)也可驗(yàn)證程序梁?jiǎn)卧獦?gòu)建以及GUYAN縮減法的正確性。

表3 CEFR堆芯燃料組件ANSYS與自研程序固有頻率計(jì)算結(jié)果對(duì)比Table 3 Comparison between nature frequenciesin CEFR core fuel assembly by ANSYS and code

2.3 CEFR單根組件空氣中時(shí)程分析驗(yàn)證

在堆芯組件大規(guī)模非線性碰撞情況下,ANSYS計(jì)算精度低且計(jì)算時(shí)間長(zhǎng),但單根組件不存在碰撞時(shí)其計(jì)算結(jié)果是具有可信度的。通過(guò)單根組件在自研程序和ANSYS中計(jì)算結(jié)果的對(duì)比,可驗(yàn)證時(shí)程分析法的正確性。

圖3 KOBE波地震輸入Fig.3 Kobe earthquake seismic input

針對(duì)2.2節(jié)中提到的CEFR燃料組件模型,采用ANSYS和自研程序分別進(jìn)行時(shí)程分析計(jì)算,其中KOBE波地震輸入如圖3所示,計(jì)算結(jié)果對(duì)比如圖4所示。由圖4可知,二者基本吻合,其中ANSYS計(jì)算的組件頂點(diǎn)處位移最大值為5.263×10-4m,自研程序計(jì)算的組件頂點(diǎn)處位移最大值為5.320×10-4m,二者相對(duì)誤差為1.08%,由此證明單組件的時(shí)程分析法是正確的。

圖4 ANSYS與自研程序?qū)EFR單組件時(shí)程分析結(jié)果對(duì)比Fig.4 Comparison between time historical calculation results calculated by ANSYS and code

2.4 RAPSODIE堆單排19組件空氣中時(shí)程分析驗(yàn)證

IAEA于1987年到1995年以法國(guó)、日本和意大利快堆堆芯抗震試驗(yàn)為基準(zhǔn)組織并開展了各成員國(guó)快堆堆芯抗震程序?qū)Ρ闰?yàn)證的大型會(huì)議活動(dòng),試驗(yàn)的組織情況列于表4。

表4 IAEA試驗(yàn)組織情況Table 4 Model reactor and seismic testconducted by IAEA

為驗(yàn)證相鄰組件碰撞時(shí)碰撞剛度和阻尼設(shè)置的正確性,采用表4中法國(guó)RAPSODIE堆空氣中單排19根組件碰撞試驗(yàn)數(shù)據(jù)對(duì)自研程序進(jìn)行驗(yàn)證。該試驗(yàn)包含燃料組件(F/A)和屏蔽組件(N/S),排列方式如圖5所示。組件上下兩個(gè)碰撞墊塊分別位于1.28 m和0.71 m處,燃料組件的1階、2階阻尼比為3%,屏蔽組件的1階、2階阻尼比為1%,梁截面參數(shù)與試驗(yàn)測(cè)得的碰撞剛度和碰撞阻尼參考日本FINDS程序的設(shè)置[2]。

在自研程序中,燃料組件與屏蔽組件底端均設(shè)置為固定約束,而燃料組件在管腳球形密封處設(shè)置一個(gè)剛度為6×106N/m的彈簧以等效間隙帶來(lái)的影響,計(jì)算固有頻率并與試驗(yàn)結(jié)果對(duì)比,如表5所列。可發(fā)現(xiàn),第1階固有頻率程序計(jì)算結(jié)果和試驗(yàn)結(jié)果符合良好。試驗(yàn)未給出第2階和第3階固有頻率。

圖5 RAPSODIE單排試驗(yàn)排列簡(jiǎn)圖Fig.5 Layout of RAPSODIE reactor in single row assemblies test

表5 程序計(jì)算和試驗(yàn)所測(cè)固有頻率對(duì)比Table 5 Comparison between natural frequency results calculated by code and measured from test

圖6 程序計(jì)算和試驗(yàn)所測(cè)組件最大位移對(duì)比Fig.6 Comparison between max assembly displacement calculated by code and measured from test

將水平地震波RAP087[2]輸入自研程序,在考慮組件間碰撞的情況下進(jìn)行時(shí)程計(jì)算,提取每根組件頂部的最大位移并和試驗(yàn)所測(cè)結(jié)果作對(duì)比,如圖6所示?？煽闯龆哂?jì)算結(jié)果接近,其中燃料組件的最大位移基本吻合,誤差主要出現(xiàn)在屏蔽組件上,但堆芯抗震更加關(guān)注燃料組件的沖擊與變形,因此證明程序時(shí)程分析法的計(jì)算結(jié)果是可靠的。

3 結(jié)論

本文在國(guó)內(nèi)外快堆抗震程序開發(fā)和試驗(yàn)驗(yàn)證調(diào)研的基礎(chǔ)上,建立了快堆堆芯抗震理論模型,闡述了堆芯組件等效簡(jiǎn)化以及碰撞剛度和阻尼處理方法,解決了堆芯組件多點(diǎn)復(fù)雜非線性碰撞的關(guān)鍵性問題,選擇適用于快堆堆芯組件在地震載荷作用下的隱式Newmark法進(jìn)行了時(shí)程分析計(jì)算,提出了快堆單排組件抗震分析理論算法,并將該理論算法進(jìn)行代碼開發(fā)。通過(guò)從簡(jiǎn)單懸臂梁到復(fù)雜變截面簡(jiǎn)支梁的固有頻率計(jì)算,以及單組件和單排多組件的時(shí)程分析計(jì)算,從下而上逐步驗(yàn)證了本文所建立的抗震模型和算法的正確性。后續(xù)將進(jìn)一步考慮液體對(duì)組件運(yùn)動(dòng)的流固耦合影響并將二維算法擴(kuò)展到三維,拓展程序適用范圍的同時(shí)提高計(jì)算精度。