浮動堆設備閘門瞬態熱固耦合數值模擬研究

(武漢第二船舶設計研究所，武漢 430064)

0 引言

浮動堆是國家海洋戰略的重要工程之一，為海洋開發提供永續可移動能源的保障，是我國經略海洋的必然選擇[1]。浮動堆在世界民用核能應用領域屬于難度極高的工程，其相關設備安全性研究尤為重要。浮動堆設備閘門在事故工況下，與安全殼一起包容放射性物質，是安全殼壓力邊界較薄弱的環節。浮動堆設備閘門與陸上核電站不同，需要承受海洋環境引起的沖擊、縱橫搖擺載荷，更為復雜的是LOCA事故下安全殼內溫度和壓力相應引起的熱應力載荷變化將更加劇烈[2]，在瞬態溫度載荷下導致結構破壞的一個重要原因是來自結構不同位置的溫度差[3]，而不單純是最大的溫度載荷，這種熱應力足以使設備閘門產生過量的塑性變形或斷裂[4]。

杜坤等[5]基于ANSYS有限元分析方案和法國RCC-M《壓水堆核電廠機械設備設計和建造規則》的理論，分析設備閘門在內壓、重力、地震等載荷組合的影響，開展了穩定性分析和法蘭、螺栓的強度分析。錢浩等[6]開展人員閘門在溫度、內壓、重力等載荷下的強度評定，其中溫度載荷采用了靜態分析方法。施勣等[7-8]結合土建強迫位移，采用了等效靜力的方法研究地震等載荷對人員閘門的影響。左樹春[9]通過有限元靜力分析法研究安全殼強迫位移造成法蘭密封面分離及相對錯動。采用熱固耦合分析方法針對現有核電站設備閘門結構及密封安全問題的研究較少，而且均未考慮海洋環境下搖擺等特殊工況。因此，對設備閘門在事故工況下瞬態溫度場以及相應的耦合熱應力進行研究是具有實際意義的研究方向。

本文利用有限元方法求得瞬態溫度場，結合事故工況下的壓力、搖擺、沖擊載荷，求得設備閘門耦合熱應力，采用的分析方法可為其他浮動堆設備閘門的設計與安全評估提供參考和理論指導。

1 熱應力問題基本原理

1.1 瞬態溫度場有限元求解方程

求解瞬態溫度場問題的核心是利用相應的數值方法求解線性常微分方程組，通過將經典熱力學中結構三維瞬態溫度場的控制方程進行推導，可以得出該問題以時間t為獨立變量的線性常微分方程組[10]。

(1)

矩陣C，K，P的元素由單元相應的矩陣元素集成[11]，為：

(2)

(3)

(4)

通過上式，可將時間域和空間域的偏微分方程問題在空間域內轉變成N個節點溫度φi(t)常微分方程的初值問題。

1.2 熱-固耦合分析理論

熱應力是溫度場和應力場共同作用的結果，溫度場和應力場互相作用的問題叫做熱固耦合問題[12]，耦合結果滿足物理方程、平衡方程和協調方程。熱-結構分析計算的有限元方程[13]為：

(5)

式中，M為質量矩陣；u，T分別為位移、溫度載荷；C為結構阻尼矩陣；Ct為比熱矩陣；K為結構剛度矩陣；Kt為熱傳導矩陣；F為總等效結點力列陣；Q為總等效結點熱流率向量。

2 耦合熱應力的有限元分析

本文采用順序耦合法進行分析，首先對設備閘門進行熱分析，在得到瞬態溫度場分布以后，將求得的節點溫度作為體載荷施加到結構中，并利用給定的載荷和位移條件對耦合應力進行分析。

2.1 工況說明

浮動堆在正常運行工況，受到船體搖擺、自重作用，極端事故工況時，需要考慮受沖擊影響下，LOCA事故的高溫高壓影響。LOCA是壓水堆核電站的設計基準事故，當船用核動力裝置發生LOCA事故時，大量高溫高壓流體從破口噴入安全殼，導致堆艙內溫度和壓力急劇上升[10]。由于前10 s出現溫度及壓力峰值，本文選取冷斷斷裂情況下前10 s的溫度及內壓進行分析。

正常運行工況下，安全殼內壁初始溫度為50 ℃，外表面溫度為22 ℃。冷斷斷裂情況下，安全殼內溫度及壓力變化如表1所示。其余設計載荷為搖擺載荷、沖擊載荷、自重、壓力、螺栓預緊力，如表2所示。

表1 安全殼內溫度及壓力

表2 載荷說明

在考慮船舶搖擺對換料蓋的影響時，假設換料蓋為一質點，搖擺的中心位于水線面、船舶中心線與船舯的交點處。將其運動函數簡化為角位移函數[14]：

(6)

(7)

本文計算時，將角位移函數求解二階導數[15]，將該載荷轉化為瞬時加速度施加。搖擺運動加速度控制方程如下：

(8)

(9)

設備閘門上下法蘭材料選用ASME boiler and pressure vessel code,sectionⅡ,Material D篇中的SA-738Gr.B級鋼板，其熱導率性能列于表3。

表3 SA-738Gr.B力學參數

設備閘門螺栓材料選用SA-193 B7。根據ASME boiler and pressure vessel code,sectionⅡ,Material D篇，當溫度為158 ℃時材料基本力學性能見表4。

表4 SA-193 B7力學參數

2.2 建模及約束條件

本文采用有限元軟件Ansys Workbench19.2進行模擬分析，分別建立設備閘門本體、主螺栓、筒體以及下法蘭模型，并裝配成為整體。計算模型中球面蓋板、法蘭和螺栓均使用三維實體單元Solid 186；法蘭之間的接觸及螺栓與法蘭之間的接觸采用接觸單元Targe 170和Conta 174；螺栓預緊力采用Prets 179單元[16-17]，各處的摩擦系數取為0.2。在建模過程中，本文對結構進行了如下簡化：

(1)結構中的密封墊片只作密封用，不是主要受力件，建模不予考慮，不作評估；

(2)設備閘門受搖擺和沖擊載荷，具有1/2對稱性，在保證計算精度的前提下，為減小計算量，采用1/2模型進行建模。

2.3 瞬態溫度場計算以及應力計算

本文結構熱固耦合分析分為3個階段：第1階段先分析設備閘門正常運行工況下熱穩態溫度場；第2階段分析結構升溫過程中的瞬時溫度場;第3階段把所得溫度場作為體載荷施加到結構整體熱應力分析中。

(1)第1階段計算設備閘門熱穩態分布。正常運行工況下安全殼內壁初始溫度為50 ℃，外表面溫度為22 ℃，與空氣對流換熱，換熱系數為5×10-6W/(mm2·℃)。此時，設備閘門溫度分布比較均勻，螺栓端部溫度最低，與設備閘門內壁最大溫差約為2 ℃，圖1示出初始熱穩態時設備閘門溫度場分布。

圖1 初始熱穩態時設備閘門溫度場分布示意

(2)第2階段計算整個升溫和降溫過程的各瞬態溫度場。冷斷斷裂情況下，設備閘門內表面經過3 s升溫到208.55 ℃，隨后逐漸降低至170.20 ℃。計算時，溫度載荷的施加分為10個載荷步，設置結束時間10 s，采用子部設置，每一子步計算10個瞬態熱結果，開啟時間積分設置。

圖2是第3 s瞬時溫度峰值時溫度場分布。此時設備閘門由內表面至外表面呈現較為陡峭的梯度分布，特別是在內表面區域溫度變化最大。可見,在安全殼內突然升溫時，由于熱量不能快速傳遞至外表面，形成熱量在設備閘門內表面堆積的情況。

圖2 第3 s瞬時溫度峰值時溫度場分布

表5示出設備閘門內外密封面、封頭外表面、螺桿及螺母的溫升情況。在10 s內，內密封面和球殼外表面溫度變化要大于其他區域，溫度差引起的熱應力影響也應大于其他區域。

表5 設備閘門各區域溫度

(3)第3階段將所得的節點溫度映射到結構中，并結合給定的結構應力和位移條件求解。圖3是設備閘門在第3 s時，設備閘門von Mises及周圍典型螺栓應力云圖。由于設備閘門需要承受搖擺、沖擊等不均勻載荷影響，因此，選取設備閘門周邊5個典型位置作為應力監測點。

圖3 第3 s時von Mises應力云圖

表6示出設備閘門球面蓋、下法蘭、螺栓A、螺栓B、螺栓C、螺栓D、螺栓E的von Mises應力時程。可以看出，設備閘門球面蓋、下法蘭應力變化隨溫度變化明顯，螺栓由于溫升量較小，應力變換隨溫度變化不明顯，出現峰值應力的時間點并不相同。球面蓋、下法蘭最大von Mises應力發生時間為3 s，最大值為452.87 MPa及429.35 MPa，螺栓B的von Mises應力較其他螺栓高，發生時間為8 s，最大值為214.54 MPa。

表6 設備閘門各部件von Mises應力時程

考慮瞬態溫度載荷后，設備閘門整體von Mises應力明顯增大。與采用穩態熱應力計算相比，采用瞬態熱固耦合方法可以更為準確地反映設備閘門內部結構溫度分布，分析出由于球面蓋及下法蘭和螺栓最大von Mises應力發生時刻，從而更為精確地分析設備閘門在復雜海洋環境中的安全性能。

3 結構耦合熱應力有限元結果分析

3.1 結構耦合熱應力結構影響分析

基于ASME boiler and pressure vessel code,section Ⅲ,Rules for construction of nuclear facility components中表NE-3221-1應力強度限制的要求，采用第三強度理論校核設備閘門結構強度進行校核。分別在球面蓋、法蘭及螺栓前10 s最大等效應力處開展應力分類評定，分離出一次總體薄膜應力Pm、彎曲應力PL+Pb、二次應力PL+Pb+Q，如圖4～7所示。

圖4 球面蓋封頭應力線性化路徑G01

圖5 球面蓋上法蘭應力線性化路徑G02

圖6 下法蘭應力線性化路徑F01

圖7 螺栓應力線性化路徑L01

采用瞬態熱應力載荷開展應力評定時，應考慮全時間周期內線性化路徑上的應力變化。球面蓋在3 s時von Mises應力最大，應予以重點考慮。

圖8示出10 s內球面蓋法蘭彎曲應力值，可以看出球面蓋法蘭G02兩端的彎曲應力變化較內側大，受溫度載荷影響，內側彎曲應力在0～3 s時快速升高，隨后變化趨于平穩，在10 s時達到峰值，完全應力最大值為93.82 MPa，整體變化趨勢與內壓變化相近。

圖8 0～10 s球面蓋法蘭彎曲應力值云圖

圖9為10 s內球面蓋封頭二次應力值，在0～2 s時封頭內表面二次應力快速升高，達到峰值148.89 MPa，而后逐漸降低。同一時刻封頭內表面受溫度載荷影響，內表面二次應力值要明顯大于其他區域。圖10,11示出10 s螺栓平均應力和最大應力值云圖，在3 s前螺栓平均應力和最大應力值增速明顯，在第8 s時，平均應力和最大應力值均達到最大值，整體變化趨勢與內壓變化相近。綜合考慮整體瞬態溫度分析，可以確定在設備閘門內側區域，由于結構溫度差產生的熱應力而引起的應力變化占主導地位，必須予以重視。

圖9 0～10 s球面蓋封頭二次應力值云圖

圖10 0～10 s螺栓平均應力值云圖

圖11 0～10 s螺栓最大應力值云圖

表7示出設備閘門應力分類評定結果，在第4 s時球面蓋封頭完全應力比值達到95.91%，可以認為在該時刻，設備閘門處于最危險狀態，應針對該區域進行加厚處理或采用隔熱涂層等措施，以降低熱應力的影響。

表7 設備閘門應力分類評定結果

(續表7)

3.2 結構耦合熱應力密封影響分析

設備閘門整體受非對稱載荷，選取變形量最大時刻的最大變形截面進行分析[18-20]。球面蓋與下法蘭呈分離的趨勢，密封面截面發生明顯的扭轉。圖12為該截面處法蘭密封面形變示意，表8為該截面內外密封面上下表面位移時程表。在第9 s時，內密封面最大分離值為0.4 mm，外密封面最大分離值為0.21 mm，第10 s變化趨于穩定。

圖12 法蘭密封面形變示意

表8 內外密封面上下表面位移變化

在考慮瞬態熱應力的情況下，內密封面存在約17 ℃的升溫，而且密封面分離值較大，密封性能下降，外密封圈升溫較低，密封面分離值較小，密封性能影響較小。因此，進行設備閘門密封性試驗時，應根據計算的最危險結果增加升溫狀態下密封圈試驗。

4 結論

本文計算浮動堆設備閘門正常運行工況下的熱穩態溫度場，并開展事故工況前10 s時刻的瞬態熱應力、熱應力最大值的變化趨勢以及整體結構和主螺栓上的應力狀態研究，得到接近真實的模擬結構工作狀態。

(1)采用瞬態熱固耦合方法計算可以較為準確地分析設備閘門各部件最大von Mises應力、完全應力、二次應力極值發生時刻，更準確地發現設計薄弱環節，并予以改進。

(2)在10 s內，設備閘門溫度場由內表面至外表面呈現較為陡峭的梯度分布，特別是在內表面區域溫度變化最大。可見在安全殼內突然升溫，熱量不能快速傳遞至外表面，形成熱能在設備閘門內表面堆積的情況。計算結果顯示，設備閘門內表面彎曲應力、二次應力受溫度載荷變化顯著。

(3)在考慮瞬態熱應力的情況下，在第9 s時，內密封面存在約17 ℃的升溫，而且密封面分離值約為0.4 mm，初步確定可能發生密封失效時的狀態，為改進密封性試驗提供數據支撐。

(4)本文采用等效靜力的方法來處理沖擊問題，取值也相對保守，但在工程實際中，船用核動力裝置的沖擊由船舶撞擊、直升機墜落等外部事件引起，還需開展沖擊響應分析，以便更為精確分析設備閘門的安全性。