基于格林函數(shù)方法的核部件疲勞分析方法研究

謝 海，邵雪嬌，張毅雄，盧喜豐，艾紅雷，白曉明，高世卿

(中國核動(dòng)力研究設(shè)計(jì)院 核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610213)

疲勞失效是核電站設(shè)備的主要失效模式之一，在核電工程設(shè)計(jì)中一直備受關(guān)注。疲勞損傷可導(dǎo)致核電站設(shè)備或管道萌生裂紋并持續(xù)擴(kuò)展，嚴(yán)重情況下將導(dǎo)致冷卻劑泄漏等后果。核一級(jí)部件出現(xiàn)疲勞失效的原因主要有兩種，一種是機(jī)械振動(dòng)產(chǎn)生的機(jī)械疲勞，一般發(fā)生在回路系統(tǒng)上小支管等部件，另一種是由核電廠運(yùn)行過程中冷卻劑溫度和壓力變化導(dǎo)致的熱疲勞失效，如穩(wěn)壓器波動(dòng)管熱分層效應(yīng)導(dǎo)致的疲勞失效和一回路主管道上支管流導(dǎo)致的熱疲勞失效等。機(jī)械振動(dòng)通常導(dǎo)致高周疲勞，可通過改變結(jié)構(gòu)自身的固有頻率來降低出現(xiàn)疲勞失效的風(fēng)險(xiǎn)；而熱疲勞通常為低周疲勞，由熱疲勞導(dǎo)致的失效通常不易發(fā)現(xiàn)而難以防止，因此對(duì)熱疲勞失效敏感位置進(jìn)行監(jiān)測和評(píng)估是保證核電廠運(yùn)行安全的重要手段。

為防止核島壓力邊界出現(xiàn)的疲勞失效，設(shè)計(jì)階段即從保守的角度對(duì)設(shè)備和管道進(jìn)行疲勞計(jì)算分析。首先要根據(jù)熱工水力計(jì)算和核電廠運(yùn)行的經(jīng)驗(yàn)反饋提出各種設(shè)計(jì)瞬態(tài)，并給出這些瞬態(tài)發(fā)生的次數(shù)。事實(shí)上，大部分設(shè)計(jì)瞬態(tài)在波動(dòng)幅度和發(fā)生次數(shù)的假設(shè)上都具有相當(dāng)?shù)谋Ｊ匦浴Ｆ浯?，由于瞬態(tài)發(fā)生次序是未知的，因此疲勞分析中只能保守考慮最嚴(yán)厲的瞬態(tài)組合方式，進(jìn)一步增加了分析的保守性。從設(shè)計(jì)角度，一定的保守性是必要的，但合理挖掘疲勞分析的保守余量將給核電廠的經(jīng)濟(jì)性帶來好處，而采用實(shí)際監(jiān)測的數(shù)據(jù)進(jìn)行疲勞計(jì)算可達(dá)到這一目的。但采用實(shí)際檢測數(shù)據(jù)存在著一定的困難，主要的困難來自于熱應(yīng)力計(jì)算：若采用商業(yè)有限元軟件計(jì)算結(jié)構(gòu)熱應(yīng)力場，由于實(shí)際冷卻劑溫度和壓力時(shí)間歷程很長，計(jì)算將耗費(fèi)大量資源而失去了可行性。因此，能否快速計(jì)算結(jié)構(gòu)熱應(yīng)力是采用實(shí)際監(jiān)測數(shù)據(jù)進(jìn)行疲勞計(jì)算的關(guān)鍵所在。

近年來大量研究表明，相對(duì)于空氣環(huán)境中的疲勞壽命，反應(yīng)堆冷卻劑環(huán)境下金屬的疲勞壽命有顯著降低[1-4]，業(yè)界對(duì)此進(jìn)行了大量研究并提出兩種評(píng)價(jià)方法，本文采用其中應(yīng)用較多的環(huán)境疲勞修正因子(Fen)對(duì)考慮冷卻劑環(huán)境的疲勞壽命進(jìn)行評(píng)估。

本文將以核電廠一回路主管道支管模型為例，采用格林函數(shù)法計(jì)算實(shí)現(xiàn)快速計(jì)算熱應(yīng)力，并用快速傅里葉變換求解其中涉及的卷積過程，以進(jìn)一步加速熱應(yīng)力計(jì)算速度。闡述后續(xù)涉及的應(yīng)力線性化、應(yīng)力配對(duì)、應(yīng)力極值點(diǎn)選取和考慮環(huán)境影響的Fen求解過程，最后以算例驗(yàn)證計(jì)算方法的正確性。

1 格林函數(shù)方法

在數(shù)學(xué)理論中，格林函數(shù)是一種用于求解有初始條件或邊界條件的非齊次微分方程的函數(shù)，最早由英國人喬治·格林于1828年提出。經(jīng)過長期研究和發(fā)展，格林函數(shù)已成為物理學(xué)中的一個(gè)重要函數(shù)，有時(shí)也稱為源函數(shù)或影響函數(shù)。格林函數(shù)方法的原理可簡單理解為：任意一段時(shí)程溫度加載都可離散為一系列溫度變化，每段溫度變化產(chǎn)生一個(gè)向前發(fā)展的波動(dòng)，這個(gè)波動(dòng)可以是溫度或是應(yīng)力的波動(dòng)，而總的溫度或應(yīng)力效果為該時(shí)間點(diǎn)之前各段的線性疊加[5-7]：

(1)

若考慮溫度沖擊td時(shí)間后分析位置的應(yīng)力達(dá)到穩(wěn)定狀態(tài)，則有：

σT(P,t)=G0(P)(φ(t)-Tref)+

(2)

對(duì)于一般結(jié)構(gòu)，無法通過解析的方法獲取具體點(diǎn)的格林函數(shù)，可通過有限元軟件計(jì)算階躍熱沖擊下結(jié)構(gòu)的響應(yīng)獲取關(guān)心位置的格林函數(shù)[7]。獲取格林函數(shù)后，通過數(shù)值計(jì)算積分可獲取溫度場和應(yīng)力場。

格林函數(shù)法是基于線性疊加理論，而實(shí)際計(jì)算中材料的物理特性是隨溫度變化的，這導(dǎo)致了非線性的情況，因此嚴(yán)格來說格林函數(shù)法不適用于熱應(yīng)力計(jì)算。為從工程上解決這一問題，可保守地采用包絡(luò)的方法，在生成格林函數(shù)時(shí)采用最大的Eα(E為彈性模量，α為線膨脹系數(shù))。

2 應(yīng)力線性化

根據(jù)ASME第Ⅲ卷[8]NB-3228.5的要求，當(dāng)一次加二次應(yīng)力強(qiáng)度范圍Sn超過3Sm時(shí)，要對(duì)交變應(yīng)力Sa進(jìn)行適應(yīng)性修正，原因是基于彈性力學(xué)計(jì)算得出的交變應(yīng)力不能直接用于在S-N曲線上獲得許用循環(huán)次數(shù)，必須進(jìn)行彈塑性修正。在此過程中需計(jì)算評(píng)定路徑首尾節(jié)點(diǎn)的一次加二次應(yīng)力強(qiáng)度范圍Sn，即薄膜加彎曲強(qiáng)度范圍，因此需對(duì)分析路徑的應(yīng)力進(jìn)行線性化處理，將沿路徑的應(yīng)力分為薄膜應(yīng)力σm、彎曲應(yīng)力σb和峰值應(yīng)力σp。應(yīng)力線性化示意圖如圖1所示。

三維實(shí)體單元的應(yīng)力線性化公式如下。

1) 薄膜應(yīng)力

(3)

若考慮路徑上等距分布的n個(gè)點(diǎn)，則可簡化為：

(4)

2) 彎曲應(yīng)力

(5)

圖1 應(yīng)力線性化示意圖Fig.1 Schematic of stress linearization

同樣考慮路徑上等距分布的n個(gè)點(diǎn)，則可簡化為：

(6)

3) 峰值應(yīng)力

(7)

對(duì)于二維軸對(duì)稱情況，應(yīng)力分量線性化更為復(fù)雜，在本文中不進(jìn)行闡述。從以上公式可看出，為進(jìn)行應(yīng)力線性化，需要計(jì)算評(píng)定路徑上若干點(diǎn)的應(yīng)力分量，本文將計(jì)算路徑上等距6個(gè)點(diǎn)的應(yīng)力分量。從ANSYS軟件說明文檔[9]可知，ANSYS在路徑上取49個(gè)點(diǎn)(包括路徑起點(diǎn)和終點(diǎn))進(jìn)行應(yīng)力線性化，為與ANSYS保持一致，本文將6個(gè)點(diǎn)的應(yīng)力分量等距插值成49個(gè)點(diǎn)。插值方法可選擇分段線性插值或拉格朗日插值，也可選擇4次多項(xiàng)式最小二乘法擬合后將起點(diǎn)和終點(diǎn)強(qiáng)制指定為原值。本文將這幾種方法的結(jié)果與ANSYS進(jìn)行對(duì)比以確定合適的算法。

3 應(yīng)力極值點(diǎn)選取和應(yīng)力配對(duì)

要確定應(yīng)力循環(huán)應(yīng)先找出局部的應(yīng)力極值點(diǎn)，單應(yīng)力分量(標(biāo)量)很易找出，但對(duì)于多軸應(yīng)力狀態(tài)簡單的算法不具有可行性。本文擬分別采用兩種算法進(jìn)行極值點(diǎn)搜索，第一種是Rubberband方法(橡皮筋算法)，EPRI在開發(fā)FatiguePro4.0時(shí)使用了該方法，文獻(xiàn)[6]中對(duì)該方法已有詳細(xì)描述并提供了偽代碼；另一種是采用Signed tresca stress(帶正負(fù)號(hào)的tresca等效應(yīng)力)方法，即將tresca等效應(yīng)力乘以絕對(duì)值最大主應(yīng)力的正負(fù)號(hào)作為應(yīng)力時(shí)程，再尋找極值時(shí)刻點(diǎn)。

由于多軸應(yīng)力狀態(tài)的復(fù)雜性，選取的應(yīng)力極值點(diǎn)與其他時(shí)間點(diǎn)組合后可能無法得到最大的應(yīng)力范圍，因此文獻(xiàn)[6]中采用了Time Window方法(時(shí)間窗口方法)，即將極值點(diǎn)前后擴(kuò)大一定的范圍，取極值點(diǎn)附近的一系列時(shí)間點(diǎn)代替極值點(diǎn)進(jìn)行交變應(yīng)力范圍計(jì)算，如圖2、3所示，范圍大小可由材料的疲勞持久極限確定。

圖2 極值時(shí)刻點(diǎn)在圖表上的示意圖Fig.2 Schematic of stress extremes

圖3 Time Window示意圖Fig.3 Schematic of Time Window

獲取極值時(shí)刻點(diǎn)后應(yīng)采用雨流計(jì)數(shù)法進(jìn)行應(yīng)力時(shí)間點(diǎn)配對(duì)，雨流計(jì)數(shù)法按照ISO-12110-2：2013[10]實(shí)現(xiàn)，配對(duì)時(shí)對(duì)累積疲勞使用系數(shù)無貢獻(xiàn)的波動(dòng)進(jìn)行過濾。

4 Fen計(jì)算

環(huán)境疲勞修正因子Fen是室溫下空氣中的疲勞壽命Nair,RT與反應(yīng)堆運(yùn)行條件下水中的疲勞壽命NPWR之比，即Fen=Nair,RT/NPWR。

根據(jù)美國NRC研究文獻(xiàn)NUREG/CR6909第1版(2018)[4]，F(xiàn)en表達(dá)式中考慮了冷卻劑溫度、應(yīng)變率和水中溶解氧含量的影響，對(duì)于低合金鋼材料還應(yīng)考慮材料中硫含量的影響。本文只給出不銹鋼材料的Fen表達(dá)式：

(8)

(9)

(10)

對(duì)于PWR環(huán)境中的水化學(xué)情況，有O*=0.29。對(duì)于不銹鋼，若應(yīng)變幅值小于0.1%或交變應(yīng)力小于195 MPa，認(rèn)為壓水堆環(huán)境不會(huì)導(dǎo)致疲勞壽命降低，因此Fen=1。但當(dāng)采用修正率方法計(jì)算轉(zhuǎn)換應(yīng)變率時(shí)該規(guī)定不適用。

(11)

其中：Fen,i為ti-1至ti時(shí)刻的Fen；n為1個(gè)應(yīng)力循環(huán)之間的總時(shí)間步數(shù)；Δεi為應(yīng)變變化量，應(yīng)變減少時(shí)Δεi=0。

5 算例和驗(yàn)證

本文采用核電站典型部件-主管道安注接管作為分析模型，為驗(yàn)證疲勞分析過程涉及的算法，構(gòu)造了一個(gè)包括多個(gè)溫度變化歷程的瞬態(tài)，每個(gè)變化的幅度隨機(jī)產(chǎn)生。算例中不包括管道熱膨脹等機(jī)械外載。

1) 幾何、有限元模型和材料

有限元模型如圖4所示，模型包括主管道及安注接管，材料為X2CrNiMo18.12(控氮)，材料特性可從RCC-M規(guī)范[12]附錄ZI中獲得。生成格林函數(shù)時(shí)使用20 ℃的彈性模量和350 ℃下的熱導(dǎo)率和線膨脹系數(shù)。為使格林函數(shù)法計(jì)算出的熱應(yīng)力與有限元法結(jié)果存在可比性，用有限元法計(jì)算瞬態(tài)熱應(yīng)力時(shí)也采用與生成格林函數(shù)時(shí)相同的值而不采用隨時(shí)間變化的材料性能。

圖4 有限元模型示意圖Fig.4 Finite element model diagram

2) 熱格林函數(shù)和壓力應(yīng)力傳遞函數(shù)生成

發(fā)生安注時(shí)安注管內(nèi)壁面將承受安注水的溫度，安注水進(jìn)入主管道與冷卻劑攪混，嚴(yán)格來說主管道流動(dòng)方向下游將承受攪混后的冷卻劑溫度，該溫度可通過CFD軟件計(jì)算得出。但在工程上不會(huì)做如此詳細(xì)的分析，而直接認(rèn)為主管道內(nèi)壁仍承受主冷卻劑的溫度，而支管承受安注水的溫度。熱邊界條件劃分如圖5所示。

圖5 熱邊界條件Fig.5 Thermal boundary condition

對(duì)于存在兩個(gè)不同邊界條件的模型，計(jì)算格林函數(shù)時(shí)應(yīng)分別施加相應(yīng)的溫度沖擊，同時(shí)在另一邊界施加T=Tref的溫度沖擊。計(jì)算實(shí)際熱應(yīng)力場或溫度場時(shí)，分別計(jì)算各自邊界承受的瞬態(tài)下的值，然后疊加獲得實(shí)際值。在本例中，取Tref=0 ℃，令圖5中熱邊界H1和熱邊界H2的熱交換系數(shù)均為20 kW/(℃·m)，生成格林函數(shù)時(shí)分別施加0～50 ℃的瞬時(shí)溫度沖擊，同時(shí)在另一邊界保持T=0 ℃。圖6為分析路徑內(nèi)表面節(jié)點(diǎn)y方向應(yīng)力Sy隨時(shí)間變化的曲線。從圖6可見，分析點(diǎn)對(duì)兩個(gè)邊界的溫度沖擊響應(yīng)不同，由于分析點(diǎn)位于H1邊界，因此H1邊界對(duì)應(yīng)的格林函數(shù)對(duì)分析點(diǎn)的影響更大、更迅速。

圖6 內(nèi)表面節(jié)點(diǎn)y方向應(yīng)力Sy隨時(shí)間的變化Fig.6 y component stress Sy of inner node versus time

本例中壓力應(yīng)力傳遞函數(shù)較為簡單，只有1個(gè)壓力邊界條件，在管道內(nèi)表面施加單位壓力并在主管和支管截?cái)嗝媸┘拥刃Ф祟^力。

3) 熱應(yīng)力和溫度場計(jì)算結(jié)果對(duì)比

施加在主管道和安注接管內(nèi)壁的溫度和壓力變化示于圖7。為更好地驗(yàn)證算法，假設(shè)溫度的變化曲線較為極端和隨機(jī)。

本文采用FORTRAN 95編寫格林函數(shù)法程序，采用正反快速傅里葉變換求解涉及的圓卷積計(jì)算。為應(yīng)對(duì)超長瞬態(tài)的情況，采用分段卷積的方法，保證不出現(xiàn)數(shù)組大小超出編譯器支持的情況。

內(nèi)表面節(jié)點(diǎn)溫度格林函數(shù)法計(jì)算結(jié)果和有限元法結(jié)果對(duì)比如圖8所示。從圖8可見，格林函數(shù)法計(jì)算溫度場與有限元法結(jié)果一致。內(nèi)表面節(jié)點(diǎn)y方向應(yīng)力Sy的計(jì)算結(jié)果對(duì)比示于圖9。從圖9可見，熱應(yīng)力計(jì)算精度較高，極值點(diǎn)的應(yīng)力值相差在1%之內(nèi)，但格林函數(shù)法出現(xiàn)極值點(diǎn)的時(shí)刻相比有限元法略有滯后，對(duì)疲勞計(jì)算的影響可忽略不計(jì)。

從計(jì)算時(shí)間看，格林函數(shù)法的計(jì)算時(shí)間主要受瞬態(tài)總時(shí)長影響，本例中瞬態(tài)總時(shí)長為336 000 s，且有2個(gè)熱邊界，需分別計(jì)算熱應(yīng)力再進(jìn)行疊加。由于采用了計(jì)算圓卷積代替直

圖7 測試算例冷卻劑溫度和壓力時(shí)間歷程Fig.7 Coolant temperature and pressure history of test case

圖8 內(nèi)表面節(jié)點(diǎn)溫度對(duì)比Fig.8 Comparison of inner node temperature

圖9 內(nèi)表面節(jié)點(diǎn)Sy對(duì)比Fig.9 Comparison of Sy of inner node

接計(jì)算線卷積的方法，熱應(yīng)力計(jì)算小于10 s。采用有限元軟件計(jì)算熱應(yīng)力的時(shí)間長短主要由計(jì)算模型規(guī)模、總計(jì)算步長和模型復(fù)雜程度決定，本例中計(jì)算模型規(guī)模僅為不足20 000個(gè)8節(jié)點(diǎn)三維實(shí)體單元，采用順序耦合的方法，總計(jì)算步數(shù)為800(包括30個(gè)計(jì)算步，共800個(gè)計(jì)算子步)，采用惠普工作站HPZ840(Intel Xeon E5-2687W V3+96G內(nèi)存)，使用8核SMP方式和Sparse Direct Solver(INCORE)，最終計(jì)算時(shí)間約為147 min。

4) 應(yīng)力線性化結(jié)果對(duì)比

在本算例中，計(jì)算了沿路徑6個(gè)節(jié)點(diǎn)的應(yīng)力，并分別采用分段線性插值、拉格朗日插值或4次多項(xiàng)式最小二乘法擬合后將起點(diǎn)和終點(diǎn)強(qiáng)制指定為原值的方法進(jìn)行應(yīng)力線性化，將6個(gè)點(diǎn)的數(shù)據(jù)處理為49個(gè)點(diǎn)的數(shù)據(jù)。從本質(zhì)上來說，只有6個(gè)點(diǎn)的數(shù)據(jù)進(jìn)行線性化無法獲得與ANSYS完全一致的結(jié)果，只能在3種方法中找出誤差最小的算法。

按3種方法對(duì)沿路徑應(yīng)力進(jìn)行插值后對(duì)比，由于格林函數(shù)法計(jì)算的結(jié)果相比有限元法有所滯后，只有都選取極值時(shí)刻的值有比較意義，故有限元法選取t=20 s，而格林函數(shù)法取t=23 s的結(jié)果。線性化結(jié)果列于表1。線性化結(jié)果僅用于Sn(即σm+σb)的計(jì)算(圖10)，因此選擇拉格朗日插值更為合適。本文僅給出了1個(gè)時(shí)刻的線性化結(jié)果，其他時(shí)刻可能出現(xiàn)不同的結(jié)論。

表1 線性化應(yīng)力結(jié)果對(duì)比Table 1 Comparison of linearization

圖10 線性化應(yīng)力結(jié)果對(duì)比(t=23 s)Fig.10 Comparison of stress linearization results (t=23 s)

圖11 應(yīng)力極值點(diǎn)查找示意圖(第3和第4極值點(diǎn))Fig.11 Schematic of stress extremes searching (the third and fourth extreme point)

5) 應(yīng)力極值點(diǎn)查找和應(yīng)力配對(duì)

根據(jù)Rubberband算法和Signed tresca stress算法編寫了相應(yīng)的FORTRAN程序，并根據(jù)ISO-12110-2：2013實(shí)現(xiàn)了雨流計(jì)數(shù)。ANSYS等商業(yè)軟件不具備相應(yīng)的功能，無法進(jìn)行對(duì)比，圖11為根據(jù)兩種算法找出的應(yīng)力極值點(diǎn)示意圖，表2為找出的應(yīng)力極值點(diǎn)，表3為根據(jù)執(zhí)行雨流計(jì)數(shù)法程序后的應(yīng)力配對(duì)結(jié)果。在本例中，兩種方法的結(jié)果一致。

表2 應(yīng)力極值點(diǎn)查找結(jié)果Table 2 Result of stress extremes searching

表3 兩種算法應(yīng)力配對(duì)結(jié)果對(duì)比Table 3 Comparison of stress state pairing by two algorithms

6)Fen計(jì)算和疲勞使用系數(shù)結(jié)果

計(jì)算Fen時(shí)，由于已找出所有應(yīng)力極值點(diǎn)并進(jìn)行了配對(duì)，只需根據(jù)每對(duì)組合的時(shí)間歷程選取對(duì)應(yīng)時(shí)間范圍的Δε和Fen,iΔε，即可計(jì)算出該組合的Fen。圖12為正確選取所需計(jì)算范圍的示意圖。如圖12a所示，當(dāng)組合的兩個(gè)時(shí)間點(diǎn)中間相隔多個(gè)應(yīng)力極值點(diǎn)時(shí)，只需積分部分區(qū)域，如t1至t1下一個(gè)峰值點(diǎn)，t2前一個(gè)谷值點(diǎn)至t2時(shí)刻。如圖12b所示，當(dāng)t1至t2中間不存在其他極值點(diǎn)但應(yīng)變減少時(shí)，積分區(qū)域不包括下降區(qū)域。

表4為累積疲勞使用系數(shù)計(jì)算結(jié)果。最終結(jié)果不考慮EAF的累積疲勞使用系數(shù)時(shí)為CUF=∑PUi=0.268，考慮EAF時(shí)則為CUFen=∑(PUi·Fen,i)=0.534。

6 結(jié)論和展望

從以上對(duì)比結(jié)果可看出，在不考慮非線性因素的情況下，格林函數(shù)法用于計(jì)算熱應(yīng)力精確度較高，計(jì)算速度和所需的存儲(chǔ)空間相比有限元法優(yōu)勢巨大，且隨著模型的增加優(yōu)勢越明顯。本文采用FORTRAN語言編制格林函數(shù)法計(jì)算熱應(yīng)力程序，并利用正反快速傅里葉變換加速計(jì)算速度，極大減少了其中數(shù)值積分所需時(shí)間。疲勞分析后續(xù)涉及的應(yīng)力線性化、應(yīng)力極值點(diǎn)查找、雨流計(jì)數(shù)法配對(duì)和計(jì)算Fen流程等都已用FORTRAN程序?qū)崿F(xiàn)，與有限元法的對(duì)比結(jié)果很好地證明了算法和程序的正確性。

圖12 計(jì)算Fen積分區(qū)域示意圖Fig.12 Schematic of integral zone for Fen

表4 疲勞使用系數(shù)計(jì)算結(jié)果Table 4 Calculation result of fatigue usage factor

但同時(shí)也必須看到，由于無法處理熱應(yīng)力計(jì)算中的非線性因素，在實(shí)際工程中格林函數(shù)法具有一定的局限性。對(duì)于非線性因素不顯著的情況，格林函數(shù)法可快速、高效地計(jì)算出熱應(yīng)力，為疲勞分析提供數(shù)據(jù)。為解決格林函數(shù)法無法應(yīng)對(duì)非線性因素的問題，近年來已有部分學(xué)者通過權(quán)系數(shù)等方法對(duì)格林函數(shù)進(jìn)行修正來解決該問題[13-15]，但生成對(duì)應(yīng)權(quán)系數(shù)的方法并不簡單，后續(xù)的研究可從這方面著手，進(jìn)一步提高格林函數(shù)法在實(shí)際工程中的應(yīng)用價(jià)值。