反应堆压力容器密封系统瞬态热弹塑性接触的高效解法

核 动 力 工 程 Nuclear Power Engineering 第 26 卷 第 1 期 2 0 0 5 年 2 月 Vol. 26. No.1 Feb. 2 0 0 5 文章编号：0258-0926(2005)01-0040-05 反应堆压力容器密封系统瞬态热弹塑性 接触的高效解法 徐铭宇 1，林腾蛟 1，李润方 1，杜雪松 1，李水安 2，杨 宇 2 （1. 重庆大学机械传动国家重点实验室，400044；2．中国核动力研究设计院核反应堆系统设计技术国家重点实验室，成都，610041） 摘要：本文在自行开发的三维瞬态密封分析程序的基础上，采用弹塑性小变形问题的非线性求解 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 、 接触问题的混合求解算法、接触传热有限元混合法、线性方程组求解器等方法，提高了计算效率，并对三维 瞬态密封分析程序中相应的计算和分析模块进行了修改和补充。应用改进后的分析程序对压力容器模拟体进 行密封分析，计算效率提高了一倍以上。 关键词：压力容器；弹塑性问题；接触问题；热传导；有限元法 中图分类号：TH123 文献标识码：A 1 引 言 近年来，国内外学者在压力容器密封分析方 面已作了广泛的研究，得出了很多有重要价值的 结论。H Sacher等考虑温度对螺栓及法兰的影响， 采用 5 种有限元模型分别用二维轴对称模型及三 维模型计算了反应堆压力容器的应力及变形，并 对各种计算结果进行了比较[1]。曲家棣等对压力 容器法兰进行了热态密封分析与试验，给出了热 态工况时螺栓载荷增量的计算方法，并对螺栓温 度滞后进行了细致的测量，分析结果与试验值相 符合[2]。李润方等人把法兰接触面的表面非线性 与弹塑性结构的材料非线性耦合，开发出二维和 局部三维瞬态热弹塑性接触分析程序，分析计算 了压力容器法兰的密封性能[3、4]。左卫东等对高 温气冷堆压力容器主法兰结构进行了弹塑性接触 计算，利用有限元分步加载技术模拟了主螺栓预 紧以及加压过程中主法兰的应力和位移情况[5]。 随着计算机技术的迅猛发展，使求解大规模 的压力容器有限元分析模型成为可能，求解效率 的提高也随之成为值得关注的问题。以前开发的 反应堆压力容器三维瞬态密封分析程序，由于采 用波前求解器求解线性方程组，且没有进行单元 编号优化，解题规模较大时波宽很大，不但计算 速度慢，有时还无法顺利求解。 本文采用弹塑性小变形问题的非线性求解方 法、接触问题的混合求解算法、接触传热有限元 混合法、线性方程组求解器等方法，提高了计算 效率，并对三维瞬态密封分析程序中相应的计算 和分析模块进行了修改和补充。 2 反应堆压力容器密封分析的计算方法 2.1 密封分析基本方程[3] 密封的基本方程可以表示为 )(T ufHRPuK ≡+++= R (1) 式中，K为系统刚度矩阵，N/m；u为节点位移向 量，m；P 为已知载荷向量，N；R为待定接触力 向量，N； TR 为温度载荷，N；H 为密封环反力 向量，N；f(u)为右端项向量，N。 2.2 非线性迭代格式 非线性方程常用的解法有初刚度法、变刚度 法、初应力法和初应变法。对压力容器(PV)进行 弹塑性分析时，由于塑性变形仅出现在密封区域， 可能进入塑性的单元极少，且塑性应变较小，同 时考虑到接触分析中形成柔度矩阵需花费大量的 机时，多次形成刚度矩阵是不经济的，求解时采 取初刚度法较为合理。因此在增量加载及塑性迭 代过程中，本文采用修正的 Newton-Raphson法求 收稿日期：2003-06-20；修回日期：2003-12-01 基金项目：国防科工委项目子课题(科工委[1999]149号) 徐铭宇等：反应堆压力容器密封系统瞬态热弹塑性接触的高效解法 41 解非线性方程式(1)，迭代格式为 )}()( )(dv{ )( )()( V )()(1- 0 )(-1 0 )( nn nnT nn uHuR ufBK uKu ∆∆−∆∆− −−= −=∆ ∫ σ ψ (2) )n()n()1n( uuu ∆+=+ (3) 式中， 0K 为初始切线刚度矩阵，N/m； )( )n(uψ 为 失衡力向量，N；B为应变矩阵，m-1；σ 为应力 向量，N/m2。 2.3 接触问题混合求解方法 由于接触边界存在非线性，接触问题的求解 也是一个反复迭代的过程。为了减少运算工作 量，采用混合求解方法。即首先把整个系统的 总刚度方程凝缩到接触边界，形成关于接触内 力的柔度方程，整个接触迭代过程中只需要修 改柔度矩阵。RPV 接触区较小，柔度矩阵的阶 次远比刚度矩阵的阶次低，明显提高了迭代效 率。引入整体平衡条件后的柔度方程为 ⎭⎬ ⎫ ⎩⎨ ⎧ ∆ −∆−= ⎭⎬ ⎫ ⎩⎨ ⎧ ∆⎥⎦ ⎤⎢⎣ ⎡ − − ei 1ii ei iC d 0 P S u R Q FF εP (4) 式中，F 为柔度矩阵，m/N；FC 为刚体位移变 换矩阵；Q为接触力变换矩阵； idR 为接触力增 量向量，N；Δuei为附加刚体位移向量，m； ΔSpi 为外载荷作用下接触点对的相对间距向 量，m； 1−iε 为初始间隙向量，m；ΔPei为外载 荷增量向量，N。 2.4 瞬态温度场计算方法[6] 在瞬态温度场计算时，空间量纲用有限元离 散，采用伽辽金加权残数法建立迭代格式；时间 量纲用有限差分离散，采用直接瞬时积分法建立 迭代格式。考虑接触传热后，瞬态温度场的基本 方程为 QqTCTK +=+ �T (5) 式中，KT 为导热矩阵，W/K；C 为变温矩阵， W·s/K；T为节点温度向量，K；T�为节点温度对 时间的导数向量；q为广义节点热流量向量，W； Q为接触点未知热流量向量，W。 为了提高计算效率，式(5)中的接触点未知热 流量采用如下的温度柔度方程求解 TT PQF ∆−= (6) 式中，FT为温度柔度矩阵，K/W；ΔPT为已知 节点热流量引起的接触点温度差向量，K。 当接触状态不连续时，根据接触传热状态 (连续传热、夹层传热、夹层绝热)的定解条件， 修改柔度矩阵，即可求得接触点未知热流量 Q。 3 线性方程组的求解方法[7] RPV密封分析过程中需要反复求解刚度方程 组；形成式(4)中的柔度矩阵和式(6)中的温度柔度 阵也必须反复求解有效温度刚度方程组。因此， 求解器的合理选择以及求解方法的改进将大大提 高求解效率。三维瞬态密封分析程序有带分块功 能的三角分解法(LDLT)和波前法两种求解器，具 有求解数十万个自由度的线性方程组的能力。 对于三角分解法，求解的效率很大程度上取 决于方程系数矩阵的带宽，需要合理地对单元或 节点进行编号以减小带宽。本文提出如下的方法： (1)考虑接触问题的特殊性，接触体之间的单 元和节点是相对独立的。通过依次分别对各接触 物体进行单元和节点编号，可以大大缩小系数 矩阵的带宽。 (2)若问题的求解规模较大，需要对各物体 分别进行节点带宽优化，将稀疏矩阵重新排序， 降低它的带宽和带廓，缩短求解时间。考虑到 接触问题，采用混合求解法需要反复求解刚度 方程组以得到接触点对上的位移差，形成柔度 矩阵。因此，回代过程可以作如下简化：分别 统计优化后每个物体接触点出现的最小编号， 回代求解时只需求解每个物体从接触点最小编 号自由度对应的方程开始到该物体最后一个自 由度对应的方程结束。这样，若节点优化后接 触点编号在方程中位置靠后时，计算效率会有 较大的提高。 (3)若问题求解规模不大，可以对每一个接 触物体先编非接触区域节点号，后编可能接触 区域的节点号。当形成柔度矩阵时，回代只在 接触区域内进行。由于接触区域仅占整个求解 区域的很小部分，回代过程将明显加快。 (4)当方程数较多或带宽较大时，则采用分 块消元法求解方程。分块数与分配给方程求解 的内存大小有关。当分块数不只一块时，分解 后的左端项必须写入文件，回代求解时再从文件 读出，将降低求解效率。因此应合理确定方程求 核 动 力 工 程 Vol. 26. No. 1. 2005 42 解内存，若能使位移求解或温度求解时方程求 解仅分一块，将大大减少求解时间。 4 RPV密封分析结果 RPV模拟体由圆柱形壳体、带平台半球形顶 盖、平底封头和座圈组成；顶盖和筒体由 30根主 螺栓相联结。根据结构的对称性，计算模型取容 器环向的 1/30；大部分区域采用 8节点六面体单 元；局部过渡区域采用 6 节点三棱柱单元。RPV 模拟体三维有限元网格如图 1a所示，计算模型的 单元总数为 3876，节点总数为 5442，接触点对 数为 88。 压力容器运行工况为：①升温 0.95h，同时加 内压至规定值，继续升温至 3.133h,水温至 300℃； ②恒温至 6.733h；③降温至 8.65h，水温至 100℃。 利用改进后的计算程序，对反应堆压力容器 模拟体的预紧和运行工况进行瞬态密封性能仿真 分析，并与改进设计程序前的计算结果进行了比 较。结果表明，改进后的分析程序比原分析程序 的计算效率提高了近两倍。 4.1 预紧阶段计算结果 预紧后容器及接触区的变形如图 1b所示。螺 栓及接触区附近变形最为显著。其中螺栓伸长了 0.123mm；顶盖法兰外缘顺时针偏转了 3.87′； 筒体法兰外缘逆时针偏转了 4.01′。此时法兰密 封面内部部分接触点已经分离，其它接触点处于 径向滑动状态。内环槽相对分离量增加 0.026mm， 外环槽相对分离量减少 0.019mm，说明内环已处 于回弹状态，而外环仍处于压扁状态。预紧结束 后，最大等效应力出现在上下法兰接触区域，其 值为 430MPa。 4.2 运行过程计算结果 升压过程中，上下法兰转角增加，法兰密 封面分离的接触点数继续增加。随后的升温过程 和恒温过程上下法兰转角小幅波动。降温过程， 上法兰的转角减小，下法兰转角增加。图 2给出 了运行过程中上下法兰外缘转角随时间变化的 曲线。 图 3 给出了运行过程中密封面台阶分离量 随时间变化曲线。升压过程台阶分离量迅速增 加，随后的升温和恒温过程分离量缓慢增大， 降温初期台阶分离量增加较大。说明压力容器 图 1 反应堆压力容器模拟体三维有限元网格 及预紧后的变形 Fig. 1 3D Finite Element Mesh of RPV Model and Its Deformation after Pre-Tension 图 2 上下法兰外缘转角随时间变化曲线 Fig. 2 Curve of Rotation Angle Upper and Lower Flanges Changed with Time 图 3 台阶分离量随时间变化曲线 Fig. 3 Separation Curve of Shoulder Changed with Time 的降温过程是最危险的工况，这是由于螺栓温 度滞后所引起的螺栓预紧力下降造成的。 升温过程整个容器在内压和水温作用下向 外膨胀，内外壁最大温差为 110℃，恒温过程温 度及变形趋于均匀，降温过程容器向内收缩， 内外壁最大温差为 105℃。整个温度循环过程， 顶盖法兰的位移落后于筒体法兰，这是由于前 者壁厚较大。由于座圈远离内壁热源，温度变 化较慢，其变形总是滞后一段时间，与试验得 徐铭宇等：反应堆压力容器密封系统瞬态热弹塑性接触的高效解法 43 图 4 顶盖外壁圆弧处温度随时间变化曲线 Fig. 4 Curve of Temperature at Cover Ectotheca Changed with Time 出的规律相符合。图 4 给出了顶盖外壁圆弧中 心处温度随时间变化曲线。 5 结 论 (1)反应堆压力容器运行过程的计算结果均 具有良好的规律性，符合实际情况。表明本文方 法是正确的，程序是可靠的。 (2)改进后的三维瞬态密封分析程序新增了计 算功能，计算效率较原分析程序提高了一倍以上。 该程序不仅可以用于反应堆压力容器的密封分 析，也可用于其它同类结构的密封分析，具有广 阔的应用前景。 参考文献： [1] Sacher H, Mayk M. Structure Analysis of a Reactor Pressure Vessel by Two- and Three- Dimensional Mod- els[J]. Int J Pres Ves & Piping. 1982, (10): 147~150. [2] 曲家棣，王佩珠，窦一康等. 压力容器热态密封实验 与分析[J]. 压力容器，1997，14(4): 1~7. [3] 李润方，林腾蛟．反应堆压力容器局部三维冷热态密 封分析[J]．机械 工程 路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理 学报，1993，29(6): 104～108. [4] 李润方，李水安，李健等. 核容器法兰密封系统密封 性能数值分析[J]. 核动力工程.1994，15(1): 12～18. [5] 左卫东，于溯源，刘俊杰等. 10MW 高温气冷堆压力 容器主法兰结构的有限元接触分析. 核动力工程， 2001，22(3): 226～231. [6] 孔祥谦．热应力有限单元法分析[M]．上海：上海交 通大学出版社，1999. [7] 王秉愚．有限元法程序设计[M]．北京：北京理工 大学出版社，1991. High Efficiency Algorithm for 3D Transient Thermo- Elasto-Plastic Contact Problem in Reactor Pressure Vessel Sealing System XU Ming-yu1，LIN Teng-jiao1，LI Run-fang1， DU Xue-song1，LI Shui-an2，YANG Yu2 （1. National Key Laboratory of Mechanical Transmission，Chongqing University，400044，China； 2. National Key Laboratory of Reactor System Design Technology，Nuclear Power Institute of China，Chengdu，610041，China） Abstract：There are some complex operating cases such as high temperature and high pressure during the operating process of nuclear reactor pressure vessel. It is necessary to carry out mechanical analysis and ex- perimental investigation for its sealing ability. On the basis of the self-developed program for 3-D transient sealing analysis for nuclear reactor pressure vessel, some specific measures are presented in this paper, to en- hance the calculation efficiency in several aspects such as the non-linear solution of elasto-plastic problem, the mixed solution algorithm for contact problem as well as contact heat transfer problem and linear equation set solver. The 3-D transient sealing analysis program is amended and complemented, with which the sealing analysis result of the pressure vessel model can be obtained. The calculation results have good regularity and the calculation efficiency is twice more than before. Key words：Reactor pressure vessel，Elasto-plastic problem，Contact problem，Heat transfer，Finite element method (下转第 69页) 彭 倩等：热轧温度对 N18新锆合金板材织构的影响 69 Effect of Hot-Rolling Temperature on the Texture of N18 Zirconium Alloy Plate PENG Qian，LIU Yan-zhang，ZHAO Wen-jin，SUN Chang-long （National Key Laboratory for Nuclear Fuel and Materials，Nuclear Power Institute of China，Chengdu，610041） Abstract：The texture evolution of plates after rolling at different temperatures was investigated. The result shows that after hot-rolling, cold-rolling and re-crystallized annealing, textures of the N18 plates varied, i.e., fn and fr increased，and ft decreased. The effect of the temperature of hot rolling on the texture can be summarized as follows: when hot rolling temperature changed from 750°C to 780°C, fn decreased, and ft and fr increased slightly. But when hot rolling temperature was higher than 780°C, there was little difference in texture with the increasing of temperature. β phase can be formed in N18 when temperature was higher than 780°C, and the volume fraction of β phase increased with the increasing of temperature. The deformation of N18 at 750°C followed the HCP crystal deformation mechanism. When rolling temperature was 780°C ~820°C, the textures of N18 were mainly controlled by the deformation of α phase. But when temperature was high enough, non-crystal mechanism and crystal boundary slide mechanism were also activated, besides the slip mechanism. Crystal boundary slide mechanism can be promoted by β phase distributing along α grain boundaries. Key words：N18 Zirconium alloy，Texture，Rolling temperature 作者简介： 彭 倩(1969—)，女，副研究员。2003年毕业于四川大学金属材料专业，获硕士学位。现从事核材料研究工作。 刘彦章(1975—)，男，博士研究生。2003年毕业于中国核动力研究设计院核燃料及材料专业，获硕士学位。现从事 核材料研究工作。 赵文金(1954—)，男，研究员，1977年毕业于清华大学反应堆材料专业。现从事核材料研究工作。 ( 责任 安全质量包保责任状安全管理目标责任状8安全事故责任追究制幼儿园安全责任状占有损害赔偿请求权 编辑：孙华平) (上接第 43页) 作者简介： 徐铭宇(1978—)，女，硕士研究生。2001年毕业于重庆交通学院汽车工程系交通运输专业。现从事机械设计及理论 方面的研究。 林腾蛟(1968—)，男，副教授，博士后。1989年毕业于重庆大学机械工程学院机械设计及制造专业。现从事机械计 算机辅助工程、机械系统动力学、压力容器密封分析等方向的研究。 李润方(1935—)，男，教授，博士生导师，重庆大学机械传动国家重点实验室学术委员会常务副主任。1957年毕业 于北京科技大学机械系冶金机械专业。现从事机械动力学、计算机辅助工程分析、机械设计信息化、集成化等方向的研 究。 (责任编辑：孙华平)