基于改進(jìn)3D-R樹的流固耦合模擬網(wǎng)格插值研究

苗 雪，王昭順，朱 迎，董玲玉，吳明宇，楊 文，胡長軍

(1.北京科技大學(xué) 計(jì)算機(jī)與通信工程學(xué)院，北京 100083；2.中國原子能科學(xué)研究院，北京 102413)

反應(yīng)堆流固耦合(FSI)過程復(fù)雜，涉及多學(xué)科、受多方因素影響，其重要特征是研究兩相介質(zhì)間的交互作用，即固體在流體載荷作用下會(huì)發(fā)生形變或運(yùn)動(dòng)，而固體的形變或運(yùn)動(dòng)會(huì)改變流體載荷的分布和大小[1]。流固耦合模擬計(jì)算中，一般采用不同的求解器分別對(duì)流體域(FD)和固體域(SD)進(jìn)行求解計(jì)算，所以兩者離散后的兩套網(wǎng)格在耦合界面上是不匹配的，因此需要采用網(wǎng)格插值技術(shù)對(duì)兩套網(wǎng)格的單元和節(jié)點(diǎn)建立對(duì)應(yīng)關(guān)系以進(jìn)行耦合界面上信息的傳遞[2]。

網(wǎng)格插值整體分為局部和全局插值兩類，局部插值是根據(jù)一定的匹配準(zhǔn)則選取另一套網(wǎng)格在待插值節(jié)點(diǎn)附近的一些單元或節(jié)點(diǎn)上的載荷進(jìn)行插值計(jì)算，求出待插值節(jié)點(diǎn)對(duì)應(yīng)的載荷值；全局插值是在另一套網(wǎng)格所有節(jié)點(diǎn)載荷的基礎(chǔ)上構(gòu)造一個(gè)整體的載荷與坐標(biāo)的插值函數(shù)，根據(jù)此插值函數(shù)可直接求得待插值節(jié)點(diǎn)對(duì)應(yīng)的載荷值[2-3]。局部插值需對(duì)網(wǎng)格進(jìn)行預(yù)先搜索以確定待插值節(jié)點(diǎn)附近的網(wǎng)格單元或節(jié)點(diǎn)，分為基于節(jié)點(diǎn)和基于單元兩類插值方式。基于節(jié)點(diǎn)的插值有最鄰近點(diǎn)插值法(NNI)[4-5]和反距離加權(quán)平均法(DWA)[6-8]，廣泛應(yīng)用的單元插值有六面體Lagrange、四面體體積插值法和等參元逆變換插值法[5,9]等。全局插值主要有無限平板樣條插值(IPS)[10]、有限表面插值(FPS)[11]、薄板樣條插值(TPS)[12]和徑向基函數(shù)插值(RBF)[13]，全局插值可避免網(wǎng)格單元或節(jié)點(diǎn)間的搜索過程。

局部插值要預(yù)先對(duì)網(wǎng)格進(jìn)行搜索以確定與待插值節(jié)點(diǎn)匹配的網(wǎng)格單元或節(jié)點(diǎn)。傳統(tǒng)網(wǎng)格搜索采用直接搜索方式，即不采用或采用簡(jiǎn)單的數(shù)據(jù)結(jié)構(gòu)索引并搜索網(wǎng)格，如Naranjo等[14]采用的簡(jiǎn)單分層映射方式。高精細(xì)流固耦合模擬中流體域網(wǎng)格和固體域網(wǎng)格數(shù)量巨大，且每次迭代都涉及大規(guī)模的網(wǎng)格搜索計(jì)算，所以網(wǎng)格搜索占據(jù)整個(gè)耦合計(jì)算過程的大部分時(shí)間。采用的網(wǎng)格搜索技術(shù)會(huì)直接影響流固耦合計(jì)算效率，所以直接的搜索方式會(huì)大幅降低耦合計(jì)算效率。另一種是采用合適的數(shù)據(jù)結(jié)構(gòu)組織大規(guī)模網(wǎng)格，廣泛應(yīng)用的3D空間索引有網(wǎng)格索引[15]、八叉樹[16]、3D-R樹[17]和一些簡(jiǎn)單混合樹[18]，其中3D-R樹具有自動(dòng)平衡、適用于外存存儲(chǔ)、空間利用率高和查詢效率高等優(yōu)點(diǎn)。

中國數(shù)值堆項(xiàng)目CVR1.0的熱工水力軟件PACA對(duì)流體域進(jìn)行求解計(jì)算后會(huì)輸出孤立的網(wǎng)格節(jié)點(diǎn)，通過這些節(jié)點(diǎn)不能直接還原出流體域網(wǎng)格的拓?fù)浣Y(jié)構(gòu)。因此，采用基于節(jié)點(diǎn)插值中的DWA方法[6-8]為結(jié)構(gòu)力學(xué)軟件HARSA的固體域網(wǎng)格節(jié)點(diǎn)(待插值節(jié)點(diǎn))進(jìn)行插值計(jì)算，得到其承受的壓力。本研究采用改進(jìn)3D-R樹索引大規(guī)模流體域網(wǎng)格節(jié)點(diǎn)，并利用高效搜索算法提高固體域網(wǎng)格節(jié)點(diǎn)的插值計(jì)算效率。進(jìn)一步利用HARSA對(duì)插值結(jié)果進(jìn)行固體流致振動(dòng)計(jì)算并用Archard模型對(duì)固體振動(dòng)進(jìn)行磨損評(píng)估。

1 軟件介紹

反應(yīng)堆流固耦合指熱工水力和結(jié)構(gòu)力學(xué)間的耦合，熱工水力部分使用中國數(shù)值堆項(xiàng)目CVR1.0研發(fā)的計(jì)算流體力學(xué)軟件PACA，結(jié)構(gòu)力學(xué)部分使用CVR1.0研發(fā)的結(jié)構(gòu)力學(xué)軟件HARSA。PACA和HARSA的耦合中，前者計(jì)算關(guān)注的是燃料棒周圍流體域的變化，而后者計(jì)算關(guān)注的是作用在燃料棒上的壓力及壓力對(duì)燃料棒造成的影響。

1.1 熱工水力模擬軟件PACA

單相精細(xì)化計(jì)算流體力學(xué)熱工水力模擬軟件PACA采用精確大渦模擬模型的高精度譜元方法，可對(duì)反應(yīng)堆堆芯進(jìn)行精確的數(shù)學(xué)物理建模，能精確求解反應(yīng)堆容器內(nèi)部的壓力場(chǎng)、流場(chǎng)和溫度場(chǎng)。PACA的計(jì)算過程、網(wǎng)格模型和輸出數(shù)據(jù)結(jié)構(gòu)示于圖1。由圖1a可見，PACA的求解計(jì)算過程包括預(yù)處理、計(jì)算求解、后處理3部分。

預(yù)處理部分的“網(wǎng)格分區(qū)”將流體域離散化為四面體網(wǎng)格，且每個(gè)四面體網(wǎng)格又被細(xì)化為4個(gè)六面體網(wǎng)格，部分流體域的六面體網(wǎng)格示于圖1b。預(yù)處理結(jié)束后，每個(gè)節(jié)點(diǎn)都有1個(gè)位置坐標(biāo)(x，y，z)。之后將流體域的六面體網(wǎng)格傳遞給計(jì)算求解部分進(jìn)行求解計(jì)算，計(jì)算結(jié)束后每個(gè)節(jié)點(diǎn)都有1個(gè)速度U、壓力p和溫度T。

后處理部分會(huì)輸出求解計(jì)算生成的流體域網(wǎng)格文件，即.fld文件。如圖1c所示，.fld文件由head和body兩部分構(gòu)成。在head部分，num為總網(wǎng)格單元數(shù)，inf為body中記錄的信息類型。body的每一行記錄1個(gè)節(jié)點(diǎn)的信息：(x，y，z)為節(jié)點(diǎn)的位置坐標(biāo)，Ux、Uy、Uz為節(jié)點(diǎn)對(duì)應(yīng)的矢量速度U分別在x、y、z軸方向的分量，p和T為節(jié)點(diǎn)上的壓力和溫度。

圖1 PACA的計(jì)算過程、網(wǎng)格模型和輸出數(shù)據(jù)結(jié)構(gòu)Fig.1 Calculation process, mesh model and output data structure of PACA

1.2 結(jié)構(gòu)力學(xué)模擬軟件HARSA

反應(yīng)堆結(jié)構(gòu)力學(xué)高保真有限元模擬軟件HARSA可對(duì)堆芯組件在高溫、輻照、流體等因素作用下的力學(xué)行為進(jìn)行模擬。HARSA的計(jì)算過程、網(wǎng)格模型和輸出數(shù)據(jù)結(jié)構(gòu)示于圖2。由圖2a可見，HARSA的計(jì)算流程大致分為預(yù)處理、計(jì)算求解和后處理3部分。

圖2 HARSA的計(jì)算過程、網(wǎng)格模型和輸出數(shù)據(jù)結(jié)構(gòu)Fig.2 Calculation process, mesh model and output data structure of HARSA

預(yù)處理部分的“網(wǎng)格劃分模塊”將燃料棒離散化為六面體網(wǎng)格，部分燃料棒的六面體網(wǎng)格示于圖2b，每個(gè)網(wǎng)格節(jié)點(diǎn)都具有1個(gè)位置坐標(biāo)(x′，y′，z′)。HARSA使用1個(gè).dat文件存儲(chǔ)燃料棒上需要進(jìn)行插值計(jì)算的網(wǎng)格節(jié)點(diǎn)，即待插值節(jié)點(diǎn)，文件結(jié)構(gòu)如圖2c所示。其中，head部分記錄待插值節(jié)點(diǎn)總數(shù)sum，body的每一行記錄1個(gè)待插值節(jié)點(diǎn)的位置信息。插值算法需根據(jù)待插值節(jié)點(diǎn)附近的流體域網(wǎng)格節(jié)點(diǎn)及各節(jié)點(diǎn)的壓力p，為每個(gè)待插值節(jié)點(diǎn)計(jì)算其承受的壓力p′。

1.3 插值過程

HARSA和PACA預(yù)處理計(jì)算后的燃料棒和流體域網(wǎng)格的密度和大小均不匹配，使得兩套網(wǎng)格不能在耦合界面上直接進(jìn)行信息交換。因此，需借助網(wǎng)格插值技術(shù)將流體域網(wǎng)格上的壓力信息進(jìn)行轉(zhuǎn)換并傳遞到燃料棒的網(wǎng)格上。PACA輸出的流體域網(wǎng)格節(jié)點(diǎn)數(shù)目非常大，傳統(tǒng)插值方式非常耗時(shí)。因此，研究利用3D-R樹索引大規(guī)模流體域網(wǎng)格節(jié)點(diǎn)，以加快對(duì)燃料棒上待插值節(jié)點(diǎn)壓力的求解計(jì)算，從而提高整體流固耦合效率。網(wǎng)格插值過程如圖3所示，整個(gè)插值過程包括2個(gè)階段，即索引構(gòu)建階段和查詢處理階段。

圖3 網(wǎng)格插值過程Fig.3 Mesh interpolation process

1) 索引構(gòu)建階段

將所有流體域網(wǎng)格節(jié)點(diǎn)采用逐一插入的方式構(gòu)建3D-R樹，任意網(wǎng)格節(jié)點(diǎn)oi由位置坐標(biāo)及其對(duì)應(yīng)的速度、壓力、溫度描述。oi的插入過程為：從根結(jié)點(diǎn)root開始自頂向下查找最佳葉子結(jié)點(diǎn)并將oi插入其中，之后判斷此葉子結(jié)點(diǎn)是否溢出，溢出則執(zhí)行結(jié)點(diǎn)分裂操作并更新樹，否則繼續(xù)插入新的網(wǎng)格節(jié)點(diǎn)直至構(gòu)成完整3D-R樹。最佳葉子結(jié)點(diǎn)指插入oi后周長增加最小的葉子結(jié)點(diǎn)。

2) 查詢處理階段

對(duì)于每個(gè)待插值節(jié)點(diǎn)qj，依據(jù)高效搜索算法查找距離它最近的3個(gè)流體域網(wǎng)格節(jié)點(diǎn)，之后利用式(1)計(jì)算qj對(duì)應(yīng)的壓力pj，直至為所有待插值節(jié)點(diǎn)計(jì)算出對(duì)應(yīng)的壓力。

2 基于改進(jìn)3D-R樹的網(wǎng)格插值

根據(jù)PACA輸出的流體域網(wǎng)格節(jié)點(diǎn)無法直接還原網(wǎng)格拓?fù)浣Y(jié)構(gòu)，本研究采用基于節(jié)點(diǎn)插值中的反距離加權(quán)平均法DWA[6-8]為燃料棒上的待插值節(jié)點(diǎn)計(jì)算對(duì)應(yīng)的壓力。DWA首先為待插值節(jié)點(diǎn)查找距離它最近的幾個(gè)流體域網(wǎng)格節(jié)點(diǎn)，之后分別以這些流體域網(wǎng)格節(jié)點(diǎn)到該待插值節(jié)點(diǎn)的距離為權(quán)重，計(jì)算待插值節(jié)點(diǎn)對(duì)應(yīng)的壓力。燃料棒-流體域網(wǎng)格模型的可視化示意圖如圖4所示，左側(cè)為1組燃料棒-流體域離散后的網(wǎng)格模型，右側(cè)為其局部放大。圖4中qi為燃料棒上某待插值節(jié)點(diǎn)，對(duì)流體域網(wǎng)格節(jié)點(diǎn)搜索獲得距離qi最近的3個(gè)節(jié)點(diǎn)，分別為os、om、on，根據(jù)式(1)計(jì)算qi的壓力pi。

圖4 燃料棒-流體域網(wǎng)格模型的可視化示意圖Fig.4 Visualization schematic diagram of mesh model in fuel rod-fluid domain

(1)

式中：j=s,m,n;ps、pm、pn分別為os、om、on的壓力；ds、dm、dn分別為os、om、on與qi的距離，且dm

2.1 改進(jìn)的3D-R樹

本研究利用3D-R樹索引堆芯空間中的大規(guī)模流體域網(wǎng)格節(jié)點(diǎn)，1棵3D-R樹在2D空間中的投影示意圖如圖5所示。3D-R樹是一棵高度平衡樹，除根結(jié)點(diǎn)外任意結(jié)點(diǎn)的孩子結(jié)點(diǎn)數(shù)不得大于預(yù)設(shè)扇出值M，若結(jié)點(diǎn)孩子數(shù)大于M，則結(jié)點(diǎn)溢出，需執(zhí)行結(jié)點(diǎn)分裂操作保證樹平衡。樹中結(jié)點(diǎn)使用最小包圍盒(MBB)描述，MBB由1個(gè)6元組(xmin，xmax，ymin，ymax，zmin，zmax)表示，MBB為恰好包圍空間點(diǎn)的最小盒。3D-R樹將位置靠近的空間點(diǎn)聚集構(gòu)成葉子結(jié)點(diǎn)，它們對(duì)應(yīng)于空間中最小的包圍盒，如圖5a中的實(shí)線矩形；繼續(xù)將空間中較小的包圍盒聚集構(gòu)成上層較大的包圍盒，即中間結(jié)點(diǎn)，如圖5a中虛線矩形；如此逐層聚集最終得到1個(gè)包含所有空間點(diǎn)的最大包圍盒，即根結(jié)點(diǎn)。

a——3D-R樹的2D平面投影；b——對(duì)應(yīng)樹型結(jié)構(gòu)圖5 3D-R樹在2D空間投影結(jié)構(gòu)示意圖Fig.5 Schematic diagram of 3D-R tree projection in 2D space

3D-R樹允許結(jié)點(diǎn)間存在重疊，這會(huì)使訪問空間點(diǎn)分布不均勻的區(qū)域時(shí)出現(xiàn)大量重復(fù)路徑，影響查詢效率。傳統(tǒng)3D-R樹通過大量復(fù)雜比較計(jì)算，將分裂得到的兩個(gè)新結(jié)點(diǎn)包圍盒周長之和最小作為評(píng)判依據(jù)進(jìn)行結(jié)點(diǎn)分裂，以此尋求較小的重疊區(qū)域。但PACA的流體域網(wǎng)格節(jié)點(diǎn)數(shù)目巨大、密度大且分布相對(duì)均勻，建樹時(shí)會(huì)執(zhí)行大量的結(jié)點(diǎn)分裂操作，所以傳統(tǒng)分裂策略會(huì)嚴(yán)重影響建樹效率。為此，基于流體域網(wǎng)格節(jié)點(diǎn)密度大且分布相對(duì)均勻的特點(diǎn)，設(shè)計(jì)了一種體積均分結(jié)點(diǎn)分裂方式，其通過簡(jiǎn)單計(jì)算即可得到包圍盒體積之和較小的兩個(gè)新結(jié)點(diǎn)，有效降低了索引構(gòu)建時(shí)間。

體積均分結(jié)點(diǎn)分裂過程如圖6所示，node為1個(gè)溢出結(jié)點(diǎn)，其包圍盒(xmin，xmax，ymin，ymax，zmin，zmax)由平行于坐標(biāo)軸的3類邊構(gòu)成，它們的范圍分別為xmin→xmax、ymin→ymax、zmin→zmax，構(gòu)成集合MBBSet={(xmin，xmax)，(ymin，ymax)，(zmin，zmax)}。對(duì)于MBBSet中的第j個(gè)元素tupj=(os，oe)，首先計(jì)算其垂直分割面split，之后遍歷node的每個(gè)孩子chi，若chi位于split左側(cè)，將chi并入新結(jié)點(diǎn)N1；若位于split右側(cè)，將chi并入N2；若與split相交，將其暫存入集合S。處理完node的所有孩子后，判斷S是否為空，為空則計(jì)算N1和N2體積和Vj，否則將S中的每個(gè)孩子插入N1或N2并計(jì)算N1和N2體積和Vj。處理完MBBSet中所有元素后輸出Vj最小的N1和N2，它們是分裂得到的兩個(gè)新結(jié)點(diǎn)。

圖6 新樹結(jié)點(diǎn)分裂策略流程圖Fig.6 Flow chart of new splitting strategy for tree node

2.2 網(wǎng)格搜索

傳統(tǒng)插值采用直接搜索方式為每個(gè)待插值節(jié)點(diǎn)查找距離它最近的3個(gè)流體域網(wǎng)格節(jié)點(diǎn)，時(shí)間復(fù)雜度為O(3mn)，其中n為流體域網(wǎng)格節(jié)點(diǎn)數(shù)，m為待插值節(jié)點(diǎn)數(shù)。由于流體域網(wǎng)格節(jié)點(diǎn)數(shù)巨大，直接搜索方式非常耗時(shí)。為此，使用基于(改進(jìn))3D-R樹的快速搜索方式。

基于(改進(jìn))3D-R樹的網(wǎng)格搜索流程示于圖7，使用優(yōu)先隊(duì)列PrioQue來加速搜索。PrioQue中存儲(chǔ)樹結(jié)點(diǎn)及其到待插值節(jié)點(diǎn)qi的最短距離，距離越小優(yōu)先級(jí)越高。初始時(shí)，PrioQue中存儲(chǔ)根結(jié)點(diǎn)root及其到qi的距離。首先從PrioQue中彈出優(yōu)先級(jí)最高的元素node，若node為結(jié)點(diǎn)類型，則遍歷node的所有孩子，對(duì)每個(gè)孩子chi，計(jì)算chi到qi的距離dist(chi，qi)并入隊(duì)。若node為非結(jié)點(diǎn)類型，則它是一個(gè)最近鄰，將node存入集合R。判斷近鄰數(shù)R(集合R中的元素個(gè)數(shù))是否為3，為3則計(jì)算待插值節(jié)點(diǎn)qi的壓力pi，否則繼續(xù)從PrioQue中彈出優(yōu)先級(jí)最高的元素并重復(fù)上述過程。

圖7 基于(改進(jìn))3D-R樹的網(wǎng)格搜索流程圖Fig.7 Mesh search flow chart based on (improved) 3D-R tree

3 實(shí)驗(yàn)驗(yàn)證及分析

利用表1所列6組數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，測(cè)試了燃料棒上待插值節(jié)點(diǎn)的插值計(jì)算效率。進(jìn)一步利用HARSA計(jì)算燃料棒在流體壓力作用下的流致振動(dòng)變化，并利用Archard模型評(píng)估棒束流致振動(dòng)變化下的磨損程度。PACA預(yù)處理計(jì)算后的流體域網(wǎng)格模型及參數(shù)示于圖8。

圖8 流體域網(wǎng)格模型展示Fig.8 Visualization of fluid domain mesh model

表1 實(shí)驗(yàn)用到的數(shù)據(jù)Table 1 Data used in experiment

3.1 網(wǎng)格插值效率驗(yàn)證

本研究測(cè)試了基于直接搜索、3D-R樹和改進(jìn)3D-R樹3種插值方式(分別用BaseLine、RTree和ImRTree表示)在不同條件下的性能。

燃料棒長變化時(shí)，構(gòu)建3D-R樹和改進(jìn)3D-R樹執(zhí)行的總結(jié)點(diǎn)分裂次數(shù)示于圖9a。棒長增大時(shí)，RTree和ImRTree執(zhí)行的結(jié)點(diǎn)分裂次數(shù)均增大；棒長相同時(shí)，兩者執(zhí)行的總結(jié)點(diǎn)分裂次數(shù)相近。燃料棒長變化時(shí)，RTree和ImRTree執(zhí)行結(jié)點(diǎn)分裂消耗的總時(shí)間示于圖9b。可見，為PACA輸出的流體域網(wǎng)格節(jié)點(diǎn)構(gòu)建索引時(shí)，改進(jìn)結(jié)點(diǎn)分裂效率遠(yuǎn)高于傳統(tǒng)結(jié)點(diǎn)分裂效率。

圖9 燃料棒長變化時(shí)結(jié)點(diǎn)分裂次數(shù)和結(jié)點(diǎn)分裂時(shí)間Fig.9 node splitting times and time vs fuel rod length

燃料棒數(shù)變化時(shí)，3種網(wǎng)格插值消耗的時(shí)間示于圖10a?？梢?，BaseLine消耗時(shí)間最多，ImRTree消耗時(shí)間最少，RTree消耗時(shí)間位于兩者之間。由于流體域網(wǎng)格節(jié)點(diǎn)數(shù)量和密度均較大，建樹時(shí)要進(jìn)行大量復(fù)雜的結(jié)點(diǎn)分裂計(jì)算，因此RTree消耗的時(shí)間較多。ImRTree采用簡(jiǎn)單的體積均分結(jié)點(diǎn)分裂方式，所以消耗的時(shí)間最少。燃料棒長變化時(shí)，3種網(wǎng)格插值消耗的時(shí)間示于圖10b。棒長增大時(shí)，BaseLine的時(shí)間消耗最大，ImRTree的時(shí)間消耗仍最少，而RTree的時(shí)間消耗仍位于前兩者之間。由圖10可見，ImRTree和RTree在高精細(xì)插值計(jì)算中更具優(yōu)勢(shì)。

圖10 燃料棒數(shù)或長度變化時(shí)的插值時(shí)間Fig.10 Interpolation time vs fuel rod number or length

當(dāng)樹中結(jié)點(diǎn)孩子數(shù)大于預(yù)設(shè)扇出值M時(shí)，會(huì)執(zhí)行結(jié)點(diǎn)分裂操作及樹的更新操作。當(dāng)M設(shè)置不合適時(shí)，可能會(huì)執(zhí)行頻繁的分裂及更新操作影響建樹效率，從而影響插值計(jì)算效率。本文選取表1中組4和組6兩組數(shù)據(jù)并選取M為25、50、75，測(cè)試基于RTree和ImRTree的插值效率。M變化時(shí)插值消耗的時(shí)間如圖11所示。由于ImRTree的結(jié)點(diǎn)分裂消耗時(shí)間較少，M變化時(shí)ImRTree的插值效率仍明顯高于RTree的插值效率。當(dāng)M=25時(shí)，由于RTree和ImRTree執(zhí)行的分裂及更新操作均最少，所以兩者插值效率均最高。當(dāng)M≥50時(shí)基于ImRTree的插值效率趨于平穩(wěn)，而M增大時(shí)基于RTree的插值效率降低明顯(圖11b)。

圖11 M變化時(shí)插值消耗的時(shí)間Fig.11 Time consumed by interpolation vs M

3.2 流致振動(dòng)磨損評(píng)估

燃料棒受到周圍流體域施加的交替相間的壓力后會(huì)產(chǎn)生往復(fù)運(yùn)動(dòng)，燃料棒的這種振動(dòng)會(huì)伴隨著棒束間的磨損。本文選取表1中組3的模型，利用上述網(wǎng)格插值方式將0～3 s流體域網(wǎng)格節(jié)點(diǎn)的壓力轉(zhuǎn)換為燃料棒上待插值節(jié)點(diǎn)對(duì)應(yīng)的壓力，并將插值結(jié)果傳遞給HARSA進(jìn)行結(jié)構(gòu)力學(xué)計(jì)算，得到0～3 s燃料棒上待插值節(jié)點(diǎn)的位移、速度、加速度變化。進(jìn)一步利用Archard模型對(duì)棒束振動(dòng)造成的磨損程度進(jìn)行評(píng)估。Archard模型是工程中廣泛應(yīng)用的磨損模型之一，其形式為：

V=SFL/3H

(2)

式中：V為磨損體積；S為磨損系數(shù)；F為與接觸面垂直的接觸力；L為總的滑動(dòng)距離；H為材料硬度；3為適用于半球形磨粒的形狀因子。

燃料棒束在0～3 s的位移變化示于圖12。0～1 s，棒束位移變化明顯，最大位移出現(xiàn)在左側(cè)燃料棒頂端，為27.83 mm，最小位移發(fā)生在棒束底端，為6.854×10-22mm。如圖12b、c、d所示，1～3 s棒束的最大位移未發(fā)生變化，1～2 s棒束的最小位移降至6.411×10-22mm，2～3 s棒束的最小位移增至7.967×10-22mm。由于1～3 s各時(shí)刻棒束對(duì)應(yīng)位置壓力相差過小，圖12b、c、d棒束上位移分布相似。

圖12 燃料棒束在不同時(shí)刻的位移響應(yīng)Fig.12 Displacement response of fuel bundle at different time

燃料棒束在0～3 s的速度和加速度變化分別示于圖13、14。由圖13a、b和圖14a、b可看出，初始0～1 s，棒束在周圍流體域的壓力作用下速度(加速度)變化明顯，最大速度(加速度)出現(xiàn)在棒束頂端，為56.2 mm/s(113.5 mm/s2)；最小速度(加速度)出現(xiàn)在棒束底端，為1.384×10-21mm/s(2.794×10-21mm/s2)?？梢姡?～1 s時(shí)棒束頂端承受的壓力最大，底端承受的壓力最小，且在流體域的壓力作用下棒束開始發(fā)生振動(dòng)。

由圖13b、c、d可知，1～2 s最小速度增至3.405×10-21mm/s，2～3 s最小速度降至2.387×10-21mm/s。1～3 s棒束最小速度呈先增后降的趨勢(shì)，這是由于棒束不同位置受到的流體域壓力不同導(dǎo)致的。由圖14b、c、d可知，1～2 s最大加速度增至338.3 mm/s2，最小加速度增至1.086×10-20mm/s2；2～3 s最大加速度增至563.1 mm/s2，最小加速度增至3.131×10-20mm/s2。1～3 s棒束的最大和最小加速度均呈增大趨勢(shì)，可見棒束頂端和底端的受力是一直增大的。

圖13 燃料棒束在不同時(shí)刻的速度響應(yīng)Fig.13 Velocity response of fuel bundle at different time

圖14 燃料棒束在不同時(shí)刻的加速度響應(yīng)Fig.14 Acceleration response of fuel bundle at different time

可見，在周圍流體域壓力作用下燃料棒束的不同位置受力程度不同，這導(dǎo)致棒束不同位置上的速度和加速度響應(yīng)不同，從而導(dǎo)致棒束不同位置上的位移響應(yīng)不同。燃料棒束在這種不同程度的交替壓力作用下會(huì)產(chǎn)生往復(fù)運(yùn)動(dòng)，即發(fā)生流致振動(dòng)。

根據(jù)燃料棒束0～3 s的振動(dòng)變化，進(jìn)一步對(duì)表1中組3的兩根燃料棒進(jìn)行磨損評(píng)估。提取燃料棒繞絲與相鄰燃料棒包殼管間的接觸力為980.523 N，滑移距離分別為11.009 51 mm和11.009 37 mm，材料硬度取20，磨損系數(shù)取2×10-9mm3/J。根據(jù)式(2)所示磨損評(píng)估公式，計(jì)算得到組3棒束中左側(cè)燃料棒的磨損體積為3.598 36×10-7mm3，磨損深度為2.451 31×10-8mm；組3棒束中右側(cè)燃料棒的磨損體積為3.598 31×10-7mm3，磨損深度為2.451 31×10-8mm。

4 總結(jié)

本研究使用3D-R樹索引PACA計(jì)算輸出的大規(guī)模流體域網(wǎng)格節(jié)點(diǎn)，實(shí)現(xiàn)了對(duì)燃料棒上待插值節(jié)點(diǎn)的高效插值計(jì)算。通過分析PACA流體域網(wǎng)格特點(diǎn)和3D-R樹中的結(jié)點(diǎn)分裂方式，設(shè)計(jì)了一種簡(jiǎn)單的體積均分結(jié)點(diǎn)分裂策略。實(shí)驗(yàn)表明，基于此改進(jìn)3D-R樹的插值效率明顯高于基于3D-R樹和傳統(tǒng)插值的效率。利用HARSA根據(jù)插值結(jié)果進(jìn)行了燃料棒的流致振動(dòng)計(jì)算，并用Archard模型計(jì)算了燃料棒的磨損程度。本研究對(duì)堆芯安全設(shè)計(jì)、延壽等具有重要意義。后續(xù)工作將深入研究燃料棒的振動(dòng)對(duì)流體域產(chǎn)生的影響以設(shè)計(jì)高效的流體域網(wǎng)格更新策略。