【2017年整理】3 反应堆压力容器密封分析程序改进

【2017年整理】3 反应堆压力容器密封分析程序改进 3 反应堆压力容器密封分析程序改进 3 反应堆压力容器密封分析程序改进 为了拓展程序功能，扩大计算规模，提高求解速度，我们对程序作了多处修改和扩充。如新增了分块求解器及程序重启动功能，添加了三维可变节点二次等参元，修改了分布载荷输入方式，给出了球面螺栓连接时预紧过程接触修正。本章主要叙述了程序重启动设置、二次等参元形函数、螺栓修正方法，带宽优化模块和分块求解器(参见第四章)也是密封分析程序的主要改进之处。 3.1 现有程序的功能及缺点 现有的反应堆压力容器三维密封分析程序涉及弹塑性小变形理论、弹塑性接触理论、稳态及瞬态温度场理论、热接触理论以及多种耦合分析方法。该程序兼顾了通用性和专用性两个方面，既可以分析诸如压力容器密封系统之类复杂的三维瞬态耦合热弹塑性接触问题，也不丧失求解一般有限元问题的效率。程序通过了多个算例的考证，对压力容器模拟体密封分析的结果与试验结果吻合良好，说明了该程序理论及方法的准确性。 该程序可以处理多种性能随温度变化的材料类型;可以处理已知位移及载荷历程问题;可以处理具有处理斜支撑的约束问题;可以选择不同的屈服准则来适应各种类型的材料;可以选择不同的接触判定条件以加快接触状态的收敛;可以采用不同的非线性解法加速塑性迭代的收敛。 该程序可以采用三种类型的三维线性可变节点等参元，包括四面体等参单元、五面体等参单元、六面体等参单元及三维退化单元，以模拟复杂的边界及过渡区域。 然而，现有的三维密封分析程序还存在着一些缺点，主要包括: ? 求解刚度方程采用波前求解器，且没有考虑波宽优化。当所求问题规模较大时，计算效率较低。为了提高分析规模和计算效率，开发包含带宽优化的分块消元求解器十分必要。 ? 程序没有重启动功能。当问题规模很大需要长时间计算时，无法检查中间结果并进行修改续算。开发重启动功能可以在工况步间任意设置断点以控制程序运行。 ? 单元类型只有线性可变节点等参元，没有高次等参单元。新增二次等参单元有利于提高计算精度。 ? 关于预紧过程螺栓修正问题，当螺栓与垫片接触面为球面时，修正公式还有待进一步改进。 19 重庆大学硕士学位论文 3.2 程序重启动功能设置 考虑到核压力容器运行工况复杂，密封分析过程工况步多且存在多个非线性耦合迭代，当计算模型较大时，计算时间往往很长，有时需要几十个小时。因此，程序新增了重启动功能，可以通过断点的设置来控制程序的运行，这不仅方便地实现了工况步长的临时调整，减少了由于停电等引起的程序意外中断所浪费的计算时间，同时也节约了程序调试过程所花费的机时。 为了实现重启动功能，程序新增了重启动控制文件ModelName.rst和重启动临时文件Fearpv8.tmp两个数据文件。重启动功能不改变原有的输入数据卡，分析过程是否启用重启动功能取决于重启动控制文件是否存在，由程序自动搜索。重启动控制文件包括每次重启动计算的起始和终止工况步以及对应的时间步长等工况步信息。重启动临时文件Fearpv8.tmp中存储了所有的输入数据和需用于后续计算的中间结果。包含重启动功能的反应堆压力容器法兰密封分析主程序框图如图3.1所示。 由图3.1可知，重启动功能程序段的编写穿插在密封分析主程序模块中进行，代码编写及程序使用时需注意以下问题: ? 不考虑重启动时，计算结束后程序将自动删除所有临时文件。启用重启动功能时，在所有计算工况步完成之前必须保留临时文件以备后续所用，全部工况完成后才能将所有临时文件自动删除。 ? 由于启用重启动功能时，通过多次启动程序运行才能完整生成计算结果。为了保证所有计算结果文件的完整性，重启动运行时，计算结果应在结果文件最后一个记录后续写。因此打开文件时，需在OPEN语句中添加POSITION='APPEND' 选项。 ? 由于重启动运行时，需要重新动态分配数组并赋值。因此读重启动临时文件Fearpv8.tmp时，应首先读入与确定数组大小有关的变量，而后动态分配数组内存，最后将前次的计算结果作为初始值赋给动态数组。重启动临时文件中变量的读写次序必须保证完全一致。 ? 程序计算及调试时，保留的中间计算结果应该包含所有的重启动文件、结果文件以及临时文件，才能保证重启动运行的顺利进行。 由于接触问题的求解采用刚度方程和柔度方程相结合的混合求解方法，柔度矩阵元素需要通过反复回代求解刚度方程得到。当计算模型较大且接触点对很多时，形成柔度矩阵将花费大量的时间。为此，程序能实时将已生成的柔度矩阵元素自动保存在数据文件ModelName.fm中，当程序意外中断后，程序能从断点继续运行。 20 3 反应堆压力容器密封分析程序改进 开始 是是重启动运行否 否 输入弹塑性信息输入有限元模型信息 是输入载荷信息是否计算弹塑性问题 读重启动数据文件 否确定动态数组大小 是输入温度场信息是否计算温度场 否分配动态数组是输入接触信息是否计算接触问题 否是读重启动数据文件获取前输入螺栓信息是否计算螺栓载荷 次计算所有中间数据否是输入密封环信息是否计算密封环反力 读入时间循环步数否 带宽优化 分配动态数组、初始化变量 输入时间循环步数 输入本循环步信息 瞬态温度场分析是本时间步是否进行 接触状态是瞬态温度场计算 否需要修改时否间是(热)弹塑性接增本时间步是否进行量触问题分析弹塑性接触耦合计算循否环是修改计算模型计本时间步是否为预紧结束算螺栓接触力 否否 达到规定时间步数否 是否将本次计算所有中间数据所有时间步结束否写入重启动数据文件 是 释放动态数组 结束 图3.1 反应堆压力容器密封分析主程序框图 Fig.3.1 The main Program Flowchart of Seal Analysis for RPV 21 重庆大学硕士学位论文 3.3 三维可变节点等参元 [84][85]3.3.1 等参单元及形状函数 本文选用的三维等参单元包括4节点四面体线性单元、10节点四面体二次单 元、6节点五面体线性单元、15节点五面体二次单元、8节点六面体线性单元、20 节点六面体二次单元。以上各种三维等参单元如图3.2所示，单元的局部节点的编 号必须按照图中给定的顺序。 ζζ 44 810 97ηη 51133 6 ξξ22 (a) 4节点四面体线性单元 (b) 10节点四面体二次单元 (a) The 4-Node Tetrahedron Linear Element (b) The 10-Node Tetrahedron Quadratic Element ζζ41266410 11ηη551513ξξ 14933117 822 (c) 6节点五面体线性单元 (d) 15节点五面体二次单元 (c) The 6-Node Pentahedron Linear Element (d) The 15-Node Pentahedron Quadratic Element 855168 ζζ7η14η 1918ξ41124ξ9111 3210 (e) 8节点六面体线性单元 (f) 20节点六面体线性单元 32(e) The 8-Node Hexahedron Linear Element (f) The 20-Node Hexahedron Quadratic Element 图3.2 三维等参数单元 Fig.3.2 The 3-D Isoparametric Elements 22 3 反应堆压力容器密封分析程序改进 各单元形状函数，可以按通常的规则写出，即: N(,,,,,)i,1,2,？,ni ? 在局部坐标下建立有关参量的插值函数; (,,,,,) ? 是插值函数同次的多项式函数; N(,,,,,)i ?i在节点上的值为1，而在其余节点其值为0。 N(,,,,,)j(j,i)i 下面仅对10节点四面体单元和15节点五面体单元的形函数进行推导，其它 单元类型的形函数同理可得。 对于10节点四面体二次单元，如图3.2(b)，取插值函数为 222 (3.1) u,a，a,，a,，a,，a,，a,，a,，a,,，a,,，a,,12345678910式中 — 待定系数，由节点上的函数值所唯一确定。 au(i,1,2,？,10)ii 将插值函数(3.1)写为如下形式 10 (3.2) u,N(,,,,,)u,ii,i1式中 — 形函数。 N(,,,,,)(i,1,2,？,10)i 形函数的达式为形如式(3.1)的多项式，且满足 ， (3.3) N(,,,,,),0N(,,,,,),1ijjjiiii 式中 i — 节点的局部坐标，有 (,,,,,)iii (,,,,,),(0,0,0)(,,,,,),(1,0,0)111222 (,,,,,),(0,1,0) (,,,,,),(0,0,1)333444 (,,,,,),(0.5,0,0) (,,,,,),(0.5,0.5,0) 555666 (,,,,,),(0,0.5,0) (,,,,,),(0,0,0.5) 777888 (,,,,,),(0.5,0,0.5)(,,,,,),(0,0.5,0.5) 999101010 N(,,,,,)根据形函数的性质可以方便地推导形函数，以和N(,,,,,)i1N(,,,,,)为例。 5 N(,,,,,)在节点1处的值为1，在其它节点处的值都为零，考虑到平面1 和平面分别通过点4，9，2，6，3，10和点5，7，1,,,,,,,00.5,,,,,,,0 8。可以求得 ,,,,,,(1,,,)(0.5,,,)N,,2(1,,,,,,)(0.5,,,,,,) 1(1,,,)(0.5,,,),,,,,,(0,0,0) N(,,,,,)在节点5处的值为1，在其它节点处的值都为零，考虑到平面5 1,,,,,,,0,,0和平面分别通过点4，9，2，6，3，10和点1，7，8。可以求 得 23 重庆大学硕士学位论文 ,,,,(1,,,)N,,4,(1,,,,,,) 5(1,,,),,,,(0.5,0,0) 同理可以求得其余点的形函数，整理可得10节点四面体二次单元的形函数为: N,2(1,,,,,,)(0.5,,,,,,)N,4,(2,,1)12 N,4,(2,,1)N,4,(2,,1)34 (3.4) N,4,(1,,,,,,)N,4,,56 N,4,(1,,,,,,)N,4,(1,,,,,,)78 N,4,,N,4,,910 对于4节点四面体线性单元，如图3.2(a)，取插值函数为 (3.5) u,a，a,，a,，a,1234 同理可得形函数为: N,1,,,,N,,12 (3.6) N,,N,,34 对于15节点五面体二次单元，如图3.2(d)，取插值函数为 222,,,,,,,,u,a，a，a，a，a，a，a，a，12345678 (3.7) 2222 ,,a ，a,,，a,,,，a,,，a,,，a,,，a,,9101112131415 由图3.2(d)中各节点的坐标可以推出节点形函数，以节点1和节点7为例。 在节点1处值为1，而在其余点值为零。容易找出面过点4、N(,,,,,)1,,,01 10、5、11、6、12，面过点5、14、2、8、3、15、6、11，面1,,,,,0,，,，0.5,,0过点7、9、13，这样三个面就包含了除1点以外所有点。可以得出: ,,,,,,(1,,)(1,)(，，0.5)N,,,(1,,,,)(1,,)(,，,，0.5,) 1[(1,,)(1,)(，，0.5)],,,,,,(0,0,0) N(,,,,,)在节点7处值为1，而在其余点值为零。容易找出面过点1,,,07 4、10、5、11、6、12，面过点5、14、2、8、3、15、6、11、5，面1,,,,,0,,0过点13、1、9、3、15、6、12、4，可见这三个面就包含了除7点以外所有点。可以得出: ,,,,(1,,)(1,)N,,2,(1,,)(1,,,,) 7[(1,,)(1,)],,,,(1,0,,1) 同理可以求得其余点的形函数，整理可得15节点五面体二次单元的形函数为: N,(1,,,,)(1,,)(,，,，0.5,)N,,,(1,,)(1,,，0.5,) 12 N,,,(1,,)(1,,，0.5,)N,,(1,,,,)(1，,)(,，,,0.5,) 34 N,,,(1，,)(1,,,0.5,)N,,,(1，,)(1,,,0.5,) 56 N,2,(1,,)(1,,,,)N,2,,(1,,) 78 N,2,(1,,)(1,,,,)N,2,(1，,)(1,,,,) (3.8) 910 N,2,,(1，,)N,2,(1，,)(1,,,,) 1112 24 3 反应堆压力容器密封分析程序改进 N,(1，,)(1,,)(1,,,,)N,,(1，,)(1,,)1314 N,,(1，,)(1,,)15 对于6节点五面体线性单元，如图3.2(c)，取插值函数为 (3.9) u,a，a,，a,，a,，a,,，a,,123456 同理可得形函数为: N,0.5(1,,)(1,,,,)N,0.5,(1,,)12 (3.10) N,0.5,(1,,)N,0.5(1，,)(1,,,,)34 N,0.5,(1，,)N,0.5,(1，,)56 对于8节点六面体线性单元，如图3.2(e)，取插值函数为 (3.11) u,a，a,，a,，a,，a,,，a,, ，a,,，a,,,12345678 同理可得形函数为: N,(1,,)(1,,)(1,,)8N,(1，,)(1,,)(1,,)812 N,(1，,)(1，,)(1,,)8N,(1,,)(1，,)(1,,)834 (3.12) N,(1,,)(1,,)(1，,)8N,(1，,)(1,,)(1，,)856 N,(1，,)(1，,)(1，,)8N,(1,,)(1，,)(1，,)878 对于20节点六面体二次单元，如图3.2(f)，取插值函数为 222,,,,,,,,,,,,u,a，a，a，a，a，a，a，a，a ，a，12345678910 222222 (3.13),,,,,,,,,,,,,,, a，a，a，a，a，a，a，11121314151617 222 a,,,，a,,,，a,,,181920 同理可得形函数为: N,(1，,)(1,,)(1,,)(,,,,,,2)8N,(1,,)(1,,)(1,,)(,,,,,,,2)821 N,(1，,)(1，,)(1,,)(,，,,,,2)8N,(1,,)(1，,)(1,,)(,,，,,,,2)834 N,(1,,)(1,,)(1，,)(,,,,，,,2)8N,(1，,)(1,,)(1，,)(,,,，,,2)856 N,(1，,)(1，,)(1，,)(,，,，,,2)8N,(1,,)(1，,)(1，,)(,,，,，,,2)878 22 N,(1,,)(1,,)(1,,)4N,(1，,)(1,,)(1,,)4109 22 (3.14) N,(1,,)(1，,)(1，,)4N,(1,,)(1,,)(1,,)41112 22 N,(1,,)(1,,)(1，,)4N,(1，,)(1,,)(1，,)41413 22 N,(1,,)(1，,)(1，,)4N,(1,,)(1,,)(1，,)41516 22 N,(1,,)(1,,)(1,,)4N,(1，,)(1,,)(1,,)41718 22 N,(1，,)(1，,)(1,,)4N,(1,,)(1，,)(1,,)41920 3.3.2 高斯积分公式 各种等参单元的数值积分均采用高斯积分法。对于8节点和20节点的六面体 [83]等参单元，其高斯积分公式为 25 重庆大学硕士学位论文 nnn 1 1 1 (3.15) f,(,,,,)d,d,d,,WWWf(,,,,,),,,ijkijk,,, ,,,1 1 1,,,111ijk 式中 、、 — 加权系数; WWWjik — 一个方向上的高斯积分点数(通常取2或3)，其坐标及加权系数如n 表3.1所示。对于六面体等参元，本程序可取高斯点为8个或27个。 表3.1 8节点和20节点的六面体等参单元的积分点坐标及加权系数 Table 3.1 The Coordinates and Weight Coefficients of the Integral Points for the 8-Node and the 20-Node Hexahedron Isoparametric Elements 高斯点 加权系数 ,Wnkk 2 1 ,13 5/9 ,353 8/9 0 [86]对于4节点和10节点四面体等参单元，其高斯积分公式为 n 1 1,, 1,,,,1 (3.16) f(,,,,,)d,d,d,,Wf(,,,,,),kkkk,,, 0 0 06,1k 式中 — 加权系数; Wk — 四面体内高斯积分点数，本文取4，其坐标及加权系数如表3.2所示。 n 表3.2 4节点和10节点的四面体等参单元的积分点坐标及加权系数 Table 3.2 The Coordinates and Weight Coefficients of the Integral Points for the 4-Node and the 10-Node Tetrahedron Isoparametric Elements 高斯点 加权系数 n,4 ,,,W kkkk I 0.13819660 0.13819660 0.13819660 1/4 II 0.13819660 0.13819660 0.58541020 1/4 III 0.13819660 0.58541020 0.13819660 1/4 IV 0.58541020 0.13819660 0.13819660 1/4 对于6节点和15节点五面体等参单元，其高斯积分公式为 n 1 1 1,,1 (3.17) f(,,,,,)d,d,d,,Wf(,,,,,),kkkk,,, ,1 0 02,1k 26 3 反应堆压力容器密封分析程序改进 式中 — 加权系数; Wk — 五面体内高斯积分点数，本文取6，其坐标及加权系数如表3.3所示。 n 表3.3 6节点和15节点的五面体等参单元的积分点坐标及加权系数 Table 3.3 The Coordinates and Weight Coefficients of the Integral Points for the 6-Node and the 15-Node Pentahedron Isoparametric Elements 高斯点 加权系数 n,6 ,,,Wkkkk I、II 1/6 1/6 1/3 ,13 III、IV 1/6 2/3 1/3 ,13 V、VI 2/3 1/6 1/3 ,13 为了验证新增的三维可变节点二次等参元的形函数和高斯积分公式的正确性，本文对上述各种等参单元进行了算例验证。用改进后的压力容器密封分析程序NVSAS的计算结果与用ANSYS软件计算结果的比较，两者基本吻合，说明了新增二次等参单元是正确的。 3.4 预紧过程螺栓修正 压力容器三维瞬态耦合密封分析中，预紧工况是第一个计算工况，其分析的准确与否直接影响到以后运行工况分析的准确性。螺栓预紧过程中螺栓与螺母之间存在相对运动，要真实模拟该过程十分复杂。由于螺栓预紧方法是将螺栓分成若干组按对称方式轮流加载，加载过程为先将螺栓拉长再下旋螺母，因此本文参照螺栓预紧过程采用拉伸螺栓而后修正螺母的处理方法。加载方法是在螺母和上法兰垫片接触面上施加一对分布预紧压力。这样预紧结束后，螺栓被拉长，而法兰被压缩，螺栓与法兰之间便出现了间隙，该间隙即是螺母向下拧紧的距离。由于预紧后螺母与法兰几乎没有相对间隙且存在预紧压力，因此必须通过螺栓修正达到真实的模拟效果。 由于垫片和上法兰之间摩擦力较大，几乎不产生相对滑动，因而可将垫片与上法兰考虑成一个整体。由于压力容器通常采用平面螺栓(螺母与垫片接触面为平面)和球面螺栓(螺母与垫片的接触面为球面)两种螺栓形式，预紧过程的螺栓修正也考虑平面螺栓修正和球面螺栓修正两种情况，如图3.3所示。 27 重庆大学硕士学位论文 (a) 平面螺栓 (b) 球面螺栓 (a) Plain Bolt (b) Spherical Bolt 图3.3 螺栓修正示意图 Fig.3.3 The Sketch of the Bolt Correction 3.4.1 平面螺栓修正 当螺母与垫片接触面为平面时，假设螺母与垫片在整个变形过程中相对滑动量很小，因此修正后的螺母接触面上的接触点与垫片接触面上对应的接触点重合。由于预紧结束后，上法兰顺时针转动，下法兰逆时针转动，这将使螺栓产生弯曲，如图3.3(a)所示。图3.4为平面螺栓修正示意图，修正过程可分为以下三步: 螺栓初始轴线 螺母初始位置 位置1 垫片接触平面 修正后螺栓轴线 z R下法兰 图3.4 平面螺栓修正方法 Fig.3.4 The Correction Method for Plain Bolt ? 将变形后整个螺母由初始位置沿螺栓轴线下移到位置1，以模拟螺栓的旋紧过程移动的轴向距离。此时螺母接触面上接触点的Z向平均位移与垫片接触面上接触点的Z向平均位移基本相等。 ? 考虑到上法兰的刚性比螺栓的刚性大得多，因此仅对螺母接触面上的接触 28 3 反应堆压力容器密封分析程序改进 点进行修正，使其与垫片上对应的接触点(同一接触点对)相重合。 ? 在螺母上施加已知位移计算螺栓的弯曲变形，并将预紧压力转化为螺母与垫片接触面上的接触力。 3.4.2 球面螺栓修正 当螺母与垫片接触面为球面时，螺母和垫片将出现相对滑动，因为螺栓与法兰之间有一定的间隙，且上下法兰转角很小，因而在预紧结束后螺栓往往不会发生弯曲。施加了预紧力后螺母接触球面与垫片接触球面变形后的相对位置如图3.3(b)所示。图3.5为球面螺栓修正示意图，垫片接触球面最低点位于位置，螺O0母接触球面最低点位于位置。修正过程可分为以下四步: O1 螺栓轴线 o螺母初始位置1 垫片初始位置位置1 o位置2o02 o位置33 z 下法兰R 图3.5 球面螺栓修正方法 Fig.3.5 The Correction Method for Spherical Bolt ? 将变形后整个螺母由初始位置沿螺栓轴线移动到位置1，球面最低点由O1到。此时螺母接触面上接触点的Z向平均位移与垫片接触面上接触点的Z向平O2 均位移基本相等，二者接触球面最低点在同一水平面上。 ? 在保证螺栓不弯曲的前提下，对螺母接触面上点进行修正。将接触面沿轴 O线下移到位置2，螺母接触球面最低点在垫片接触球面上位置。 3 O? 将螺母接触面上的点修正到与中心点为的部分垫片接触球面重合。即3 O使螺母接触球面在位置2时以球面最低点为支点，球面旋转到位置3与垫片接3 触球面在同一球面上。 ? 修改螺母与垫片接触面的接触状态及初始间隙，将预紧压力转化为螺母与垫片接触面上的接触力。 29 重庆大学硕士学位论文 3.5 本章小结 本章叙述了现有压力容器三维法兰密封分析程序的功能及有待改进的主要内容，给出了具有重启动功能的密封分析程序框图及重启动程序模块的编写要点，推导了三维可变节点二次等参单元的形函数，并通过了算例验证，给出了平面螺栓和球面螺栓连接时预紧过程接触修正方法。通过以上新增功能的开发和程序模块编制，使三维法兰密封分析程序得到了进一步的完善。 30 3 反应堆压力容器密封分析程序改进 31