3 反应堆压力容器密封分析程序改进.doc

3 反应堆压力容器密封分析程序改进3 反应堆压力容器密封分析程序改进为了拓展程序功能，扩大计算规模，提高求解速度，我们对程序作了多处修改和扩充。如新增了分块求解器及程序重启动功能，添加了三维可变节点二次等参元，修改了分布载荷输入方式，给出了球面螺栓连接时预紧过程接触修正方法。本章主要叙述了程序重启动设置、二次等参元形函数、螺栓修正方法，带宽优化模块和分块求解器（参见第四章）也是密封分析程序的主要改进之处。3.1 现有程序的功能及缺点现有的反应堆压力容器三维密封分析程序涉及弹塑性小变形理论、弹塑性接触理论、稳态及瞬态温度场理论、热接触理论以及多种耦合分析方法。该程序兼顾了通用性和专用性两个方面，既可以分析诸如压力容器密封系统之类复杂的三维瞬态耦合热弹塑性接触问题，也不丧失求解一般有限元问题的效率。程序通过了多个算例的考证，对压力容器模拟体密封分析的结果与试验结果吻合良好，说明了该程序理论及方法的准确性。该程序可以处理多种性能随温度变化的材料类型；可以处理已知位移及载荷历程问题；可以处理具有处理斜支撑的约束问题；可以选择不同的屈服准则来适应各种类型的材料；可以选择不同的接触判定条件以加快接触状态的收敛；可以采用不同的非线性解法加速塑性迭代的收敛。该程序可以采用三种类型的三维线性可变节点等参元，包括四面体等参单元、五面体等参单元、六面体等参单元及三维退化单元，以模拟复杂的边界及过渡区域。然而，现有的三维密封分析程序还存在着一些缺点，主要包括： 求解刚度方程采用波前求解器，且没有考虑波宽优化。当所求问题规模较大时，计算效率较低。为了提高分析规模和计算效率，开发包含带宽优化的分块消元求解器十分必要。 程序没有重启动功能。当问题规模很大需要长时间计算时，无法检查中间结果并进行修改续算。开发重启动功能可以在工况步间任意设置断点以控制程序运行。 单元类型只有线性可变节点等参元，没有高次等参单元。新增二次等参单元有利于提高计算精度。 关于预紧过程螺栓修正问题，当螺栓与垫片接触面为球面时，修正公式还有待进一步改进。3.2 程序重启动功能设置考虑到核压力容器运行工况复杂，密封分析过程工况步多且存在多个非线性耦合迭代，当计算模型较大时，计算时间往往很长，有时需要几十个小时。因此，程序新增了重启动功能，可以通过断点的设置来控制程序的运行，这不仅方便地实现了工况步长的临时调整，减少了由于停电等引起的程序意外中断所浪费的计算时间，同时也节约了程序调试过程所花费的机时。为了实现重启动功能，程序新增了重启动控制文件ModelName.rst和重启动临时文件Fearpv8.tmp两个数据文件。重启动功能不改变原有的输入数据卡，分析过程是否启用重启动功能取决于重启动控制文件是否存在，由程序自动搜索。重启动控制文件包括每次重启动计算的起始和终止工况步以及对应的时间步长等工况步信息。重启动临时文件Fearpv8.tmp中存储了所有的输入数据和需用于后续计算的中间结果。包含重启动功能的反应堆压力容器法兰密封分析主程序框图如图3.1所示。由图3.1可知，重启动功能程序段的编写穿插在密封分析主程序模块中进行，代码编写及程序使用时需注意以下问题： 不考虑重启动时，计算结束后程序将自动删除所有临时文件。启用重启动功能时，在所有计算工况步完成之前必须保留临时文件以备后续所用，全部工况完成后才能将所有临时文件自动删除。 由于启用重启动功能时，通过多次启动程序运行才能完整生成计算结果。为了保证所有计算结果文件的完整性，重启动运行时，计算结果应在结果文件最后一个记录后续写。因此打开文件时，需在OPEN语句中添加POSITION=APPEND选项。 由于重启动运行时，需要重新动态分配数组并赋值。因此读重启动临时文件Fearpv8.tmp时，应首先读入与确定数组大小有关的变量，而后动态分配数组内存，最后将前次的计算结果作为初始值赋给动态数组。重启动临时文件中变量的读写次序必须保证完全一致。 程序计算及调试时，保留的中间计算结果应该包含所有的重启动文件、结果文件以及临时文件，才能保证重启动运行的顺利进行。由于接触问题的求解采用刚度方程和柔度方程相结合的混合求解方法，柔度矩阵元素需要通过反复回代求解刚度方程得到。当计算模型较大且接触点对很多时，形成柔度矩阵将花费大量的时间。为此，程序能实时将已生成的柔度矩阵元素自动保存在数据文件ModelName.fm中，当程序意外中断后，程序能从断点继续运行。图3.1 反应堆压力容器密封分析主程序框图Fig.3.1 The main Program Flowchart of Seal Analysis for RPV3.3 三维可变节点等参元3.3.1 等参单元及形状函数8485图3.2 三维等参数单元Fig.3.2 The 3-D Isoparametric Elements(a) 4节点四面体线性单元(a) The 4-Node Tetrahedron Linear Element(b) 10节点四面体二次单元(b) The 10-Node Tetrahedron Quadratic Element(c) 6节点五面体线性单元(c) The 6-Node Pentahedron Linear Element(e) 8节点六面体线性单元(e) The 8-Node Hexahedron Linear Element(d) 15节点五面体二次单元(d) The 15-Node Pentahedron Quadratic Element(f) 20节点六面体线性单元(f) The 20-Node Hexahedron Quadratic Element本文选用的三维等参单元包括4节点四面体线性单元、10节点四面体二次单元、6节点五面体线性单元、15节点五面体二次单元、8节点六面体线性单元、20节点六面体二次单元。以上各种三维等参单元如图3.2所示，单元的局部节点的编号必须按照图中给定的顺序。各单元形状函数，可以按通常的规则写出，即： 在局部坐标下建立有关参量的插值函数； 是插值函数同次的多项式函数；在节点上的值为1，而在其余节点其值为0。下面仅对10节点四面体单元和15节点五面体单元的形函数进行推导，其它单元类型的形函数同理可得。对于10节点四面体二次单元，如图3.2(b)，取插值函数为 (3.1)式中 待定系数，由节点上的函数值所唯一确定。将插值函数（3.1）写为如下形式 (3.2)式中 形函数。形函数的表达式为形如式（3.1）的多项式，且满足， (3.3)式中 节点的局部坐标，有根据形函数的性质可以方便地推导形函数，以和为例。在节点1处的值为1，在其它节点处的值都为零，考虑到平面和平面分别通过点4，9，2，6，3，10和点5，7，8。可以求得在节点5处的值为1，在其它节点处的值都为零，考虑到平面和平面分别通过点4，9，2，6，3，10和点1，7，8。可以求得同理可以求得其余点的形函数，整理可得10节点四面体二次单元的形函数为： (3.4)对于4节点四面体线性单元，如图3.2(a)，取插值函数为 (3.5)同理可得形函数为： (3.6)对于15节点五面体二次单元，如图3.2(d)，取插值函数为 (3.7)由图3.2(d)中各节点的坐标可以推出节点形函数，以节点1和节点7为例。在节点1处值为1，而在其余点值为零。容易找出面过点4、10、5、11、6、12，面过点5、14、2、8、3、15、6、11，面过点7、9、13，这样三个面就包含了除1点以外所有点。可以得出：在节点7处值为1，而在其余点值为零。容易找出面过点4、10、5、11、6、12，面过点5、14、2、8、3、15、6、11、5，面过点13、1、9、3、15、6、12、4，可见这三个面就包含了除7点以外所有点。可以得出：同理可以求得其余点的形函数，整理可得15节点五面体二次单元的形函数为： (3.8)对于6节点五面体线性单元，如图3.2(c)，取插值函数为 (3.9)同理可得形函数为： (3.10)对于8节点六面体线性单元，如图3.2(e)，取插值函数为 (3.11)同理可得形函数为： (3.12)对于20节点六面体二次单元，如图3.2(f)，取插值函数为 (3.13)同理可得形函数为： (3.14)3.3.2 高斯积分公式各种等参单元的数值积分均采用高斯积分法。对于8节点和20节点的六面体等参单元，其高斯积分公式为83 (3.15)式中、 加权系数； 一个方向上的高斯积分点数（通常取2或3），其坐标及加权系数如表3.1所示。对于六面体等参元，本程序可取高斯点为8个或27个。表3.1 8节点和20节点的六面体等参单元的积分点坐标及加权系数Table 3.1 The Coordinates and Weight Coefficients of the Integral Pointsfor the 8-Node and the 20-Node Hexahedron Isoparametric Elements高斯点加权系数21305/98/9对于4节点和10节点四面体等参单元，其高斯积分公式为86 (3.16)式中 加权系数； 四面体内高斯积分点数，本文取4，其坐标及加权系数如表3.2所示。表3.2 4节点和10节点的四面体等参单元的积分点坐标及加权系数Table 3.2 The Coordinates and Weight Coefficients of the Integral Pointsfor the 4-Node and the 10-Node Tetrahedron Isoparametric Elements高斯点加权系数I0.138196600.138196600.138196601/4II0.138196600.138196600.585410201/4III0.138196600.585410200.138196601/4IV0.585410200.138196600.138196601/4对于6节点和15节点五面体等参单元，其高斯积分公式为 (3.17)式中 加权系数； 五面体内高斯积分点数，本文取6，其坐标及加权系数如表3.3所示。表3.3 6节点和15节点的五面体等参单元的积分点坐标及加权系数Table 3.3 The Coordinates and Weight Coefficients of the Integral Pointsfor the 6-Node and the 15-Node Pentahedron Isoparametric Elements高斯点加权系数I、II1/61/61/3III、IV1/62/31/3V、VI2/31/61/3为了验证新增的三维可变节点二次等参元的形函数和高斯积分公式的正确性，本文对上述各种等参单元进行了算例验证。用改进后的压力容器密封分析程序NVSAS的计算结果与用ANSYS软件计算结果的比较，两者基本吻合，说明了新增二次等参单元是正确的。3.4 预紧过程螺栓修正压力容器三维瞬态耦合密封分析中，预紧工况是第一个计算工况，其分析的准确与否直接影响到以后运行工况分析的准确性。螺栓预紧过程中螺栓与螺母之间存在相对运动，要真实模拟该过程十分复杂。由于螺栓预紧方法是将螺栓分成若干组按对称方式轮流加载，加载过程为先将螺栓拉长再下旋螺母，因此本文参照螺栓预紧过程采用拉伸螺栓而后修正螺母的处理方法。加载方法是在螺母和上法兰垫片接触面上施加一对分布预紧压力。这样预紧结束后，螺栓被拉长，而法兰被压缩，螺栓与法兰之间便出现了间隙，该间隙即是螺母向下拧紧的距离。由于预紧后螺母与法兰几乎没有相对间隙且存在预紧压力，因此必须通过螺栓修正达到真实的模拟效果。由于垫片和上法兰之间摩擦力较大，几乎不产生相对滑动，因而可将垫片与上法兰考虑成一个整体。由于压力容器通常采用平面螺栓（螺母与垫片接触面为平面）和球面螺栓（螺母与垫片的接触面为球面）两种螺栓形式，预紧过程的螺栓修正也考虑平面螺栓修正和球面螺栓修正两种情况，如图3.3所示。图3.3 螺栓修正示意图Fig.3.3 The Sketch of the Bolt Correction(a) 平面螺栓(a) Plain Bolt(b) 球面螺栓(b) Spherical Bolt3.4.1 平面螺栓修正当螺母与垫片接触面为平面时，假设螺母与垫片在整个变形过程中相对滑动量很小，因此修正后的螺母接触面上的接触点与垫片接触面上对应的接触点重合。由于预紧结束后，上法兰顺时针转动，下法兰逆时针转动，这将使螺栓产生弯曲，如图3.3(a)所示。图3.4为平面螺栓修正示意图，修正过程可分为以下三步：图3.4 平面螺栓修正方法Fig.3.4 The Correction Method for Plain Bolt 将变形后整个螺母由初始位置沿螺栓轴线下移到位置1，以模拟螺栓的旋紧过程移动的轴向距离。此时螺母接触面上接触点的Z向平均位移与垫片接触面上接触点的Z向平均位移基本相等。 考虑到上法兰的刚性比螺栓的刚性大得多，因此仅对螺母接触面上的接触点进行修正，使其与垫片上对应的接触点（同一接触点对）相重合。 在螺母上施加已知位移计算螺栓的弯曲变形，并将预紧压力转化为螺母与垫片接触面上的接触力。3.4.2 球面螺栓修正当螺母与垫片接触面为球面时，螺母和垫片将出现相对滑动，因为螺栓与法兰之间有一定的间隙，且上下法兰转角很小，因