Integrated Anti-Crack Finite Element Analysis of Reactor Pressure Vessel

ATLAS是一个适用于核电工程的大规模结构快速有限元分析系统，集成了核电设计中所常用的分析功能模块，如弹塑性分析功能、大开孔优化功能、热机耦合分析功能、瞬态分析功能、时程分析功能、楼层响应谱生成功能、反应谱分析功能、应力线性化功能、压力容器防断裂评定功能等。其有限元分析系统由三部分组成，包括前处理器ATLAS Modeler，后处理器ATLAS Viewer和求解器ATLAS Solver ^[14,15]。软件的设计采用了先进的面向对象的设计模式，使用多态、虚拟继承与模板技术，在保证程序性能的前提下最大限度降低了程序模块之间的耦合性，具有非常高的可扩展性并且易维护。它最突出的特点是，在普通计算机上就可以很迅速地生成、处理上亿自由度的有限元分析模型，可适用于解决核电工程中核电设备、管道系统、核岛厂房三维有限元模型庞大、处理速度慢的相关的问题。

ATLAS前处理器能直接读取CAD软件的几何模型进行有限元网格划分等操作，并进行材料属性和边界条件的设置，再提交给ATLAS求解器或ABAQUS求解器进行计算。其前处理器具备与UG、Catia、ProE等CAD软件的接口，并支持ABAQUS、MARC等求解器。

后处理器用于有限元计算得到的结果数据显示，其具有三维图形显示非常快等一些独特功能，也可以交互式的显示数据、进行数据捕捉、对后处理数据进行标准化设置，还可以根据客户需要查看动画结果。

ATLAS Solver求解器基于高速并行原理开发，其计算收敛速度较常规有限元软件提升相当明显。

求解器中使用的有限元法的基本思路是把分析的连续体假想分割成有限个单元所组成的组合体，这些单元仅在定点处相互连接，这些连接点称为节点，单元之间只能通过节点传递内力，这些内力称为节点力，作用在节点上的载荷称为节点载荷，连续体受工况作用传递到节点上的位移称为节点位移。

以线性四面体单元为例，单元有四个节点I、J、K、M，每个节点有三个位移u、v、w和三个节点力U、V、W。节点力和节点位移的关系可表示为：

式中：[K]^e为刚度矩阵。利用结构力的平衡条件和边界条件把各个单元按照原来的结构连接起来，形成整体的有限元方程，对各个单元进行分析，建立节点载荷与节点位移的关系，以解出节点位移。最终将所有单元组合起来得到整体的方程：

式中：[K]是整体刚度矩阵；{δ}为全部节点位移组成的阵列；{R}为全部节点载荷组成的阵列。在位移法中，只有{δ}是未知的，求解该线性方程组就可得到各节点的位移，将节点位移带入各相应方程就可得到单元的应力分量。

在有限元方法中，当在分析中采用连续单元时，可得到总应力分布。因此，为得出薄膜和弯曲应力，总应力分布应线性化处理为应力分量成分，且用于计算其当量应力。如采用壳体单元，则由壳体总应力可直接得到薄膜和弯曲应力。

薄膜和弯曲应力都是在通过元件厚度的横截面上导出的，这些截面称为应力分类面(SCP)，在二维结构中，由缩小两对应侧的分类面至无穷小长度而得到应力分类线(SCL)。分类面是通过元件截面所切出的平面，分类线是通过元件截面所切出的直线，当在轴对称或平面几何中观察时，分类线都是平面。

由有限元分析采用二维或三维连续单元所得的应力结果可以采用应力积分法进行处理，沿SCLs通过壁厚对应力分量进行积分以确定薄膜和弯曲应力分量，采用这一方法由总应力分布扣除薄膜和弯曲应力分布可以直接求得峰值应力分量。采用这些分量按下式计算当量应力：

导出应力分布的薄膜、弯曲和峰值应力的方法在下面分步骤给出，用于计算的应力分量应以由SCL定位所规定的局部坐标为基础：

1)第一步，计算薄膜应力张量，薄膜应力张量是由沿应力分类线每一应力张量的平均值组成的张量，或：

2)第二步，分三方面进行弯曲应力张量计算

(1)仅对局部周向和径向(正)应力分量计算得弯曲应力，并不对平行于SCL的局部应力分量或面内的剪切应力计算得弯曲应力。

(2)仅对导致SCL扭转的剪切应力分布考虑剪切应力的线性部分(在垂直于周向平面中的面外剪切应力)。

(3)弯曲应力分量是由沿应力分类线每一应力分量的线性变化部分所组成，或：

3)第三步，计算峰值应力张量，峰值应力张量是分量等于下值中的张量：

4)第四步，以薄膜和薄膜加弯曲应力分量为基础，计算在SCL端部处的三个主应力。

5)第五步，以薄膜和薄膜加弯曲应力分量为基础，用式(1)计算在SCL端部处的当量应力。

计算得到的应力强度因子为允许的临界最大值。

本文分两种情况，临近接管、法兰区和远离接管、法兰区。

1)对于远离接管、法兰区域，计算过程如下：

(1)薄膜应力对应的应力强度因子

(2)弯曲应力对应的应力强度因子

(3)径向温度梯度对应的应力强度因子

(4)评价依据

壁厚温差可由ASME规范中的参考图得到。

2)对于临近接管、法兰的区域，计算过程如下：

(1)一次薄膜应力对应的应力强度因子

(2)一次弯曲应力对应的应力强度因子

(3)二次薄膜应力对应的应力强度因子

(4)二次弯曲应力对应的应力强度因子

(5)评价依据

根据ASME BPVC第III卷NB-3220，确定各评定截面的应力限制。对于D级使用限制，适用ASME BPVC第III卷附录F的规则。

ATLAS核电有限元软件中完成了上述分析功能的前后处理界面开发。如模型前处理界面开发如图1(a)所示，有限元应力结果处理界面开发如图1(b)所示，线性化界面开发如图2所示，防断裂校核界面开发如图3所示。

Figure 1.  Model creation and model tree of displacement distribution

Figure 2.  Linearization of stress

Figure 3.  Sample of fracture analysis

如图1(a)所示，根据模型尺寸填写网格的最大尺寸75、平均尺寸50和最小尺寸10，点击确定按钮即可划分面网格。最后我们得到节点数29 067，面网格数27 160，体网格数133 310。

在后处理器中选择左侧模型树中的基本结果—Displacement-u，并双击，然后选择工具栏中的“光滑云图”即可显示X方向位移云图，对应的v为Y方向位移、w为Z方向位移，如图1(b)所示。应力结果和温度结果查看同位移的查看方法。

按照图2所示过程，选择应力线、输入等分数确定之后，选择应力线性化文件的保存位置，则会在所选择的文件夹生成应力线性化结果文件，其后缀名为．opt。

依照图3，点击即可调出防断裂分析的界面，按照ASME BPVC III D1附录G的要求，需要输入参考无延性温度和允许RT_NDT的温度。点击确定之后，消息窗会将参考应力强度因子K_IR和计算得到的应力强度因子K_I显示出来。

不失一般性，我们首先建立一个立方体模型，其尺寸为10 mm×10 mm×10 mm，弹性模量为9.8 GPa，泊松比为0.14。立方体模型中随机分布40个大小不同的球型孔洞，整个立方体体积为1 000 mm³，所有孔洞体积和为50 mm³，孔洞体积与整个立方体体积之比为1/20。载荷分布为顶板上表面施加垂直于表面的面载荷P＝5 MPa，方向Z轴负方向。约束分布为底板下表面固定X、Y、Z方向位移，顶板分别固定4个侧面水平位移。

Figure 5.  Model (left) , loads, and boundary conditions (right)

Figure 6.  Distribution of displacements of model (left: ATLAS; right: ABAQUS)

对模型中孔洞占比1%、2%、3%、4%、5%、10%、15%、20%分别进行了计算，并计算了弹性模量，弹性模量计算公式如下：

式中：E为弹性模量，本文中为计算出的对应孔隙率的弹性模量；σ为应力，本文中为施加的面载荷；L为立方体边长，本文中为10 mm；ΔL为竖向位移，本文中为施加面载荷的面上节点平均位移。

Figure 7.  Relationship between porosity of cavities and elastic modulus

无孔洞模型的材料参数中弹性模量值为9.8 GPa，由于孔洞的存在，使得整体模型的弹性模量降低。随着孔隙率的增加，弹性模量成近似线性递减的规律^[17]。

算例模型如图8所示，该模型仅受拉力(大小为30 MPa)的作用。其中材料参数为：密度2 700 kg/m³，弹性模量：70 GPa，泊松比：0.3。网格为4 966个节点，17 264个线性四面体单元。

Figure 8.  Meshing the model with cavities

我们点击左侧列表中的“位移”，分别勾选x、y、z方向位移约束，在模型显示区选择要施加约束的面(对称边界条件)，软件左侧的选择列表会列出已经选择的面及编号，点击创建荷载边界对话框中的确定按钮，完成施加载荷和约束条件的设置，如图9所示。位移结果如图10所示，应力云图如图11所示。

Figure 9.  Loadings and boundary conditions

Figure 10.  Distribution of displacements

Figure 11.  Contour of stress around the cavities

从图11中可明显看出，孔洞导致了严重的应力集中，尽管孔洞1比孔洞2的体积要大，然而在孔洞2周围的应力集中最大，达到了76.65 MPa(途中白点处即为应力最大值处)，由此可见，靠近边界的孔洞附近容易产生疲劳裂纹。由此可以看出，应力集中现象不能仅仅考虑单一孔洞的大小，更要对整体模型进行一体化结构有限元分析^[17]。

为完成压力容器防断裂分析，我们针对压力容器设计的规范，按照容标委的规范设定一个应力分类线，然后对应力分类线上的应力强度分布进行应力分类，也就是按照力平衡原则将应力分解为线性成分和非线性成分，然后根据这些规范数据对压力容器设计做校验。

选择应力分类线的两个端点的原则是确定位于应力强度最大部位壁厚方向的一个内壁点和一个外壁点的大概位置，使得分类线使得分类应力达到最大为止。

由于铸造过程中气体的卷入或使用过程中的腐蚀，压力容器中可能存在微观孔洞现象，就会导致实际弹性模量的下降，影响疲劳寿命。另外，孔洞的大小、位置，以及形状，都会对应力有一定的影响。基于高精度三维X射线断层扫描技术对压力容器进行扫描，得到压力容器的三维模型。对压力容器真实模型进行一体化有限元分析，得到压力容器三维应力分布。然后根据防断裂分析规范，进行应力线性化，计算出应力强度因子。最后根据计算得到的应力强度因子与参考应力强度因子的比较，即可判断该压力容器是否满足防断裂的要求。因此，压力容器防断裂一体化分析方法对正在运行的压力容器采取措施进行延寿或按长寿命设计的新压力容器都能降低建造成本，提高经济效益。

我们采用实际的反应堆压力容器工程数据进行建模分析。因压力容器是对称的，取其1/4模型来分析，该模型仅受内压(大小为5 Pa)的作用。材料参数为：弹性模量：200 GPa、泊松比：0.3。

在图12中，我们添加弹性模量等材料参数，并施加载荷和约束条件。有限元计算结果位移云图如图13所示。

Figure 12.  Adding loadings and constraint conditions

Figure 13.  Distribution of displacements

按照ASME BPVC III D1附录G的要求，需要输入参考无延性温度和允许RT_NDT的温度。

该分析严格按照ASME BPVC III D1附录G的规程，针对接管、封头等不连续区和远离不连续区采用不同的公式。

本例中采用RT-RT_NDT＝10，位置为远离不连续区，并指定到上一步的应力线性化结果文件：

根据本文第一章节中的计算原理，以及第二章节中软件的操作，本案例通过软件计算后得到的应力强度因子数值的比较就可以得出是否满足防断裂要求的结论。本算例中，K_I＝ 8.456 797，K_IR ＝ 41.166 740，满足K_I<K_IR，因此该压力容器满足防断裂的要求。

1)ATLAS可很好地完成含微观孔洞结构的结构有限元分析。ATLAS计算结果的准确性可以得到保证，与ABAQUS相同模型的计算结果近似。

2)由于微观孔洞的存在，使得整体模型的弹性模量降低。随着孔隙率的增加，弹性模量成近似线性递减的规律。靠近边界的微观孔洞导致更严重的应力集中，更容易产生疲劳裂纹，这与疲劳实验的结果相吻合。

3)无需借助其他软件，自主开发的结构有限元分析软件ATLAS能独立完成反应堆压力容器防断裂有限元分析。压力容器防断裂一体化分析方法有助于解决压力容器延寿或设计长寿命压力容器等实际工程问题。