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

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

(中國核動力研究設計院 核反應堆系統設計技術重點實驗室，四川 成都 610213)

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

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

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

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

1 格林函數方法

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

(1)

若考慮溫度沖擊td時間后分析位置的應力達到穩定狀態，則有：

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

(2)

對于一般結構，無法通過解析的方法獲取具體點的格林函數，可通過有限元軟件計算階躍熱沖擊下結構的響應獲取關心位置的格林函數[7]。獲取格林函數后，通過數值計算積分可獲取溫度場和應力場。

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

2 應力線性化

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

三維實體單元的應力線性化公式如下。

1) 薄膜應力

(3)

若考慮路徑上等距分布的n個點，則可簡化為：

(4)

2) 彎曲應力

(5)

圖1 應力線性化示意圖Fig.1 Schematic of stress linearization

同樣考慮路徑上等距分布的n個點，則可簡化為：

(6)

3) 峰值應力

(7)

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

3 應力極值點選取和應力配對

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

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

圖2 極值時刻點在圖表上的示意圖Fig.2 Schematic of stress extremes

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

獲取極值時刻點后應采用雨流計數法進行應力時間點配對，雨流計數法按照ISO-12110-2：2013[10]實現，配對時對累積疲勞使用系數無貢獻的波動進行過濾。

4 Fen計算

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

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

(8)

(9)

(10)

對于PWR環境中的水化學情況，有O*=0.29。對于不銹鋼，若應變幅值小于0.1%或交變應力小于195 MPa，認為壓水堆環境不會導致疲勞壽命降低，因此Fen=1。但當采用修正率方法計算轉換應變率時該規定不適用。

(11)

其中：Fen,i為ti-1至ti時刻的Fen；n為1個應力循環之間的總時間步數；Δεi為應變變化量，應變減少時Δεi=0。

5 算例和驗證

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

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

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

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

2) 熱格林函數和壓力應力傳遞函數生成

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

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

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

圖6 內表面節點y方向應力Sy隨時間的變化Fig.6 y component stress Sy of inner node versus time

本例中壓力應力傳遞函數較為簡單，只有1個壓力邊界條件，在管道內表面施加單位壓力并在主管和支管截斷面施加等效端頭力。

3) 熱應力和溫度場計算結果對比

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

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

內表面節點溫度格林函數法計算結果和有限元法結果對比如圖8所示。從圖8可見，格林函數法計算溫度場與有限元法結果一致。內表面節點y方向應力Sy的計算結果對比示于圖9。從圖9可見，熱應力計算精度較高，極值點的應力值相差在1%之內，但格林函數法出現極值點的時刻相比有限元法略有滯后，對疲勞計算的影響可忽略不計。

從計算時間看，格林函數法的計算時間主要受瞬態總時長影響，本例中瞬態總時長為336 000 s，且有2個熱邊界，需分別計算熱應力再進行疊加。由于采用了計算圓卷積代替直

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

圖8 內表面節點溫度對比Fig.8 Comparison of inner node temperature

圖9 內表面節點Sy對比Fig.9 Comparison of Sy of inner node

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

4) 應力線性化結果對比

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

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

表1 線性化應力結果對比Table 1 Comparison of linearization

圖10 線性化應力結果對比(t=23 s)Fig.10 Comparison of stress linearization results (t=23 s)

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

5) 應力極值點查找和應力配對

根據Rubberband算法和Signed tresca stress算法編寫了相應的FORTRAN程序，并根據ISO-12110-2：2013實現了雨流計數。ANSYS等商業軟件不具備相應的功能，無法進行對比，圖11為根據兩種算法找出的應力極值點示意圖，表2為找出的應力極值點，表3為根據執行雨流計數法程序后的應力配對結果。在本例中，兩種方法的結果一致。

表2 應力極值點查找結果Table 2 Result of stress extremes searching

表3 兩種算法應力配對結果對比Table 3 Comparison of stress state pairing by two algorithms

6)Fen計算和疲勞使用系數結果

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

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

6 結論和展望

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

圖12 計算Fen積分區域示意圖Fig.12 Schematic of integral zone for Fen

表4 疲勞使用系數計算結果Table 4 Calculation result of fatigue usage factor

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