利用 矩阵指数计算方法 求解 燃耗方程.pdfVIP

CORPHY-2014——第15 届反应堆数值计算与粒子输运学术会议暨2014 年反应堆物理会议 四川成都，2014 年8 月25 日至27 日 利用矩阵指数计算方法求解燃耗方程 马纪敏 宋娟 王少华 （中国工程物理研究院核物理与化学研究所, 四川 绵阳,621900 ） 摘要：矩阵指数是求解大规模核素系统燃耗方程的重要方法之一，但是燃耗矩阵的大型、稀疏和刚性的特 点限制了矩阵指数计算方法的应用。本研究梳理了级数展开、有理近似和正交多项式逼近三种典型的矩阵 指数计算方法的原理，针对燃耗方程给出了各方法的求解公式，并给出了各方法在应用中需要注意的方面。 最后，针对铀堆和钍堆算例，以线性子链法的结果为参考结果，对几种方法进行了数值验证。 关键词：燃耗计算；矩阵指数；级数展开；有理近似；正交多项式逼近 1 引言 燃耗问题研究各种核素随时间的变化规律，求解燃耗方程对于反应堆而言至关重要。燃 耗计算方法可分为三类： ① 核素链方法。利用燃耗问题的特点，将核素转换链网分解为一系列的线性子链后解 析求解。单链的解由Bateman 形式或者一般形式给出。该方法简单直观、容易理解， 主要涉及链的构造、解析解形式和链的截断，缺点是截断会引入误差，且对于复杂 链网计算时间显著增长。 ② 差分方法。将燃耗方程写成一组方程组，利用差分方法迭代求解，如显式差分、隐 式差分及龙格-库塔方法等。该方法实现较简单，但对核素系统规模和燃耗步长限制 较多，且易发生数值不稳定性，通常在确定论燃耗计算程序（如WIMS 、CASMO ） 应用较多（只针对简单几条反应链计算）。 ③ 矩阵指数计算方法。将燃耗方程写成矩阵格式，其解可直接写为矩阵指数的形式， 求解过程转化为矩阵指数的计算。该方法对大规模核素系统较为适用，但需要根据 燃耗矩阵大型、稀疏、刚性的特点选择合适的方法。 燃耗计算方法始于十九世纪初期，早期由Harry Bateman 提出描述衰变反应线性链的方 程（即 Bateman 方程）并给出初步解析解[1]。上个世纪六七十年代，针对燃耗计算提出了 矩阵指数法，并在Origen 程序中得到应用[2-3] 。为解决燃耗矩阵的刚性问题，Origen 采用 瞬时衰变和长期平衡近似来处理短寿命核素，然后对简化过渡矩阵采用泰勒展开法计算矩阵 指数。尽管Origen 结果较为准确并且计算速度快，但是Origen 在某些特殊情况下核素密度 会发生突变，导致出现大的误差，其对短周期元素的处理会在一些特定情况下出现明显的物 理性错误[4]。矩阵指数的计算方法有很多[5]，但是燃耗矩阵的大型、稀疏和刚性的特点极 大地限制了这些方法的应用，这导致在后来的很长一段时间内基于矩阵指数的燃耗计算发展 缓慢。 在2010 年芬兰的M. Pusa 研究了燃耗矩阵的特点，发现燃耗矩阵谱点分布在复平面靠 近负实轴区域内，提出了有理函数近似方法，包括求积组近似和切比雪夫有理近似方法 （CRAM ），并应用于蒙卡堆物理程序Serpent 内[6-7] 。之后A.E. Isotalo 等人对 CRAM、 ORIGEN、线性子链法（TTA ）等进行了比较，发现综合考虑精度和速度时 CRAM 最有优 势[8]。2012 年清华的佘顶提出了拉盖尔多项式展开方法，并开发了燃耗程序Depth[9-10]。 基于矩阵指数的燃耗方程求解方法在近期快速发展，并开始得到广泛应用。 本文研究了利用矩阵指数的计算方法求解燃耗方程的理论，讨论了各方法的特点，并利 用典型反应堆算例对这些方法进行了数值验证。 基金项目：中物院核物理与化学研究所科技创新基金项目(2012GL06) ；中物院中子物理学重点实验室重点 课题(2012AC01) 作者简介：马纪敏 （1985- ），男，陕西延安人，助理研究员，博士，核科学与技术专业，majm03@yeah.net CORPHY-2014——第15 届反应堆数值计算与粒子输运学术会议暨2014 年反应堆物理会议 四川成都，2014 年8 月25 日至27 日 2 理论模型 2.1 燃耗矩阵 在一个均匀化材料区域内，核素总量是时间的函数，可简单描述为“核素浓度变化率= 产生率-消失率”。核素浓度可通过非齐次