第３２卷第３期２０１８年６月南华大学学报(自然科学版)ＪｏｕｒｎａｌｏｆＵｎｉｖｅｒｓｉｔｙｏｆＳｏｕｔｈＣｈｉｎａ(ＳｃｉｅｎｃｅａｎｄＴｅｃｈｎｏｌｏｇｙ)Ｖｏｌ ３２Ｎｏ ３Ｊｕｎ ２０１８收稿日期:２０１８－０３－１８基金项目:国家科技重大专项项目:先进乏燃料贮存技术研究(２０１５ＺＸ０６００４００２)作者简介:朱庆福(１９７３－)ꎬ男ꎬ研究员ꎬ博士ꎬ主要从事反应堆物理与临界安全方面的研究.Ｅ￣ｍａｉｌ:ｑｆｚｈｕ＠ｃｉａｅ.ａｃ.ｃｎＤＯＩ:１０ １９４３１/ｊ ｃｎｋｉ １６７３－００６２ ２０１８ ０３ ００１核燃料循环系统临界事故源项计算程序ＧＥＴＡＣ￣２.０开发朱庆福ꎬ张㊀驰ꎬ夏兆东(中国原子能科学研究院反应堆工程研究设计所ꎬ北京１０２４１３)摘㊀要:针对核燃料循环系统中不同物理形态的核燃料ꎬ建立相应的中子动力学 热工水力耦合模型ꎬ开发了用于固体㊁溶液㊁粉末㊁核燃料系统临界事故源项计算的程序ＧＥＴＡＣ￣２.０.利用国际上公开的基准实验数据对程序进行了验证ꎬ程序对功率(裂变率)峰值的计算结果与基准实验数据的相对误差在１２％以内ꎬ验证了ＧＥＴＡＣ￣２.０程序的准确性.关键词:核燃料系统ꎻ源项计算ꎻＧＥＴＡＣ￣２.０中图分类号:ＴＬ３６４.４文献标志码:Ｂ文章编号:１６７３－００６２(２０１８)０３－０００１－０７ＤｅｖｅｌｏｐｍｅｎｔｏｆＳｏｕｒｃｅＴｅｒｍＣａｌｃｕｌａｔｉｏｎＣｏｄｅＧＥＴＡＣ￣２.０ｆｏｒＣｒｉｔｉｃａｌｉｔｙＡｃｃｉｄｅｎｔｓｉｎＮｕｃｌｅａｒＣｙｃｌｅＳｙｓｔｅｍＺＨＵＱｉｎｇｆｕꎬＺＨＡＮＧＣｈｉꎬＸＩＡＺｈａｏｄｏｎｇ(ＣｈｉｎａＩｎｓｔｉｔｕｔｅｏｆＡｔｏｍｉｃＥｎｅｒｇｙꎬＢｅｉｊｉｎｇ１０２４１３ꎬＣｈｉｎａ)Ａｂｓｔｒａｃｔ:Ｂａｓｅｄｏｎｉｎｄｉｖｉｄｕａｌｎｅｕｔｒｏｎｉｃ￣ｋｉｎｅｔｉｃｓａｎｄｔｈｅｒｍａｌ￣ｈｙｄｒａｕｌｉｃｓｃｏｕｐｌｉｎｇｍｏｄｅｌｏｆｎｕｃｌｅａｒｆｕｅｌｗｉｔｈｄｉｆｆｅｒｅｎｔｐｈｙｓｉｃａｌｆｏｒｍｓꎬＧＥＴＡＣ￣２.０ｃｏｄｅｗａｓｄｅｖｅｌｏｐｅｄｆｏｒｓｏｕｒｃｅｔｅｒｍｃａｌｃｕｌａｔｉｏｎｏｆｓｏｌｉｄꎬｓｏｌｕｔｉｏｎａｎｄｐｏｗｄｅｒｆｕｅｌｉｎｔｈｅｎｕｃｌｅａｒｆｕｅｌｃｙｃｌｅｓｙｓｔｅｍｗｈｅｎｃｒｉｔｉｃａｌ￣ｉｔｙａｃｃｉｄｅｎｔｓｈａｐｐｅｎｅｄ.Ｔｈｅｒｅｌａｔｉｖｅｅｒｒｏｒｏｆｔｈｅｐｅａｋｐｏｗｅｒ(ｆｉｓｓｉｏｎｒａｔｅ)ｃａｌｃｕｌａｔｅｄｂｙＧＥＴＡＣ￣２.０ｃｏｄｅｉｓ１２％ｃｏｍｐａｒｅｄｗｉｔｈｔｈｅｂｅｎｃｈｍａｒｋｅｘｐｅｒｉｍｅｎｔａｌｄａｔａꎬｗｈｉｃｈｖｅｒｉｆｉｅｓｔｈｅａｃｃｕｒａｃｙｏｆｔｈｅｃｏｄｅ.ｋｅｙｗｏｒｄｓ:ｔｈｅｎｕｃｌｅａｒｆｕｅｌｃｙｃｌｅｓｙｓｔｅｍꎻｓｏｕｒｃｅｔｅｒｍｃａｌｃｕｌａｔｉｏｎꎻＧＥＴＡＣ￣２.００㊀引㊀言在核燃料循环的主工艺中ꎬ核燃料根据物理形态的差别可以分为固体㊁溶液和粉末.在发生临界事故时ꎬ不同物理形态的核燃料系统的功率响应各有差异ꎬ而热工性质的差别引起的反应性反馈机制的差异是不同状态核燃料系统功率响应不同的主要原因.因此ꎬ针对不同形态的核燃料系统㊀㊀㊀南华大学学报(自然科学版)２０１８年６月建立不同的物理热工耦合模型ꎬ才可对超临界事故裂变源项进行模拟计算.本文针对固体㊁粉末和溶液三种形态的核燃料ꎬ研究其超临界事故的发生机理ꎬ并结合点堆动力学方程㊁非稳态传热方程和不同形态核燃料特有的热工模型ꎬ探究其事故工况下的系统功率响应.１㊀基本理论在发生临界事故时ꎬ由于不同形态的核燃料产生的热工现象不同ꎬ导致其事故功率响应有差别.当正反应性引入使核燃料系统进入超临界状态时ꎬ系统功率会快速上升ꎬ热量在系统内积累.此时ꎬ固体系统会产生热膨胀的现象ꎬ导致中子泄漏增加ꎬ引入反应性反馈ꎻ溶液系统内的慢化剂密度会减小ꎬ同时还会因为裂变产物与慢化剂的作用产生辐照裂解气泡ꎬ导致慢化效果变差ꎬ引入反应性反馈ꎻ粉末系统慢化剂密度会减小ꎬ当温度上升到一定程度时ꎬ系统内的慢化剂会发生沸腾相变的现象ꎬ进而引入反应性反馈.整个事故过程功率响应的模拟是一个中子物理和热工水力耦合的问题.其中ꎬ物理计算结果为热工水力计算提供热源分布ꎬ热工水力计算结果为物理计算提供反应性反馈的数值.编程求解物理热工耦合问题时ꎬ将要模拟的时间划分为很多时间步长ꎬ在每一个时间步长内做一次中子物理㊁热工水力计算ꎬ并计算反应性反馈ꎬ循环求解即可得到系统的功率㊁反应性和温度等物理量随时间的变化情况.２㊀固体核燃料系统源项计算当固体核燃料系统进入超临界状态时ꎬ在初始阶段ꎬ系统功率会快速上升ꎬ热量逐渐地积累ꎬ导致固体开始产生膨胀ꎬ系统体积增大ꎬ密度减小ꎬ中子泄露增加ꎬ进而产生负的反应性反馈ꎬ致使功率下降.模拟这样一个事故过程ꎬ需要求解点堆动力学㊁非稳态传热和线弹性运动微分方程[１￣２].２.１㊀点堆动力学模型忽略临界事故过程中通量形状因子的变化ꎬ对通量进行时空分离ꎬ采用点堆动力学方程进行中子动力学部分的计算ꎬ见式(１).点堆动力学方程的求解采用三次Ｈｅｒｍｉｔ插值法.ｄＰｄｔ＝ρ－βｅｆｆ－ΔρΛＰ＋ð６ｉ＝１λｉＣｉ＋ＳｄＣｉｄｔ＝βｉꎬｅｆｆΛＰ－λｉＣｉ㊀ｉ＝１ꎬ２ꎬ ꎬ６(１)其中:Ｐ为功率密度ꎻρ为体系外加反应性ꎻDρ为反应性反馈ꎻβｅｆｆ为总缓发中子有效份额ꎻΛ为中子代时间ꎻＣｉ为第ｉ组缓发中子先驱核密度ꎻβｉꎬｅｆｆ为第ｉ组缓发中子有效份额ꎻlｉ为第ｉ组缓发中子先驱核衰变常量ꎻＳ为中子源项.２.２㊀非稳态传热模型非稳态传热方程见式(２)ꎬ可采用有限容积差分法进行数值求解[３].∂Ｔ∂ｔ＝ｋρｃÑ２Ｔ＋Ｑ ｆｉｓｓｉｏｎρｃ(２)其中:Ｔ为系统内某点的温度ꎻｋ为材料热导率ꎻρ为材料密度ꎻｃ为材料体积比热容ꎻＱｆｉｓｓｉｏｎ为裂变产热率ꎬ可采用单群稳态中子扩散方程计算得到的系统内裂变释能的空间分布近似地作为动态情形下裂变释能的空间分布.对于圆柱体:Ｑ ｆｉｓｓｉｏｎ(ｒꎬｚ)＝ＣＪ０２.４０５Ｒｒæèçöø÷ｃｏｓπＨｚæèçöø÷(３)㊀㊀对于球体:Ｑ ｆｉｓｓｉｏｎ(ｒ)＝ＣｓｉｎπＲｒæèçöø÷ｒ(４)２.３㊀线弹性膨胀模型假定固体发生的膨胀在弹性力学范围内ꎬ用线弹性运动微分方程[４]描述系统内各个指点的形变情况ꎬ见式(５).ρ∂２ｕ∂ｔ２＝Ñ[μÑｕ＋μ(Ñｕ)Ｔ＋λＩｔｒ(Ñｕ)＋(３λ＋２μ)αＴ](５)其中:ｕ为系统内某点的位移ꎻa为材料的热膨胀系数ꎻμ㊁λ为Ｌａｍｅᶄ系数ꎬ表示如下:μ＝Ｅ２１＋ν()λ＝νＥ１＋ν()１－２ν()其中:Ｅ为材料的杨氏模量ꎻν为泊松比.２.４㊀动力学参数计算动力学参数(如βｅｆｆ㊁Λ等)计算结果的准确程度在很大程度上决定了瞬态分析的精确性.在ＧＥＴＡＣ￣２.０程序中ꎬ动力学参数均采用日本原子力所开发的ＳＲＡＣ程序包进行计算ꎬ计算流程见图１.２第３２卷第３期朱庆福等:核燃料循环系统临界事故源项计算程序ＧＥＴＡＣ￣２.０开发图１㊀动力学参数计算流程图Ｆｉｇ.１㊀Ｆｌｏｗｃｈａｒｔｏｆｋｉｎｅｔｉｃｐａｒａｍｅｔｅｒｓｃａｌｃｕｌａｔｉｏｎ２.５㊀固体基准实验验证为了验证ＧＥＴＡＣ￣２.０程序的正确性ꎬ采用美国Ｇｏｄｉｖａ￣Ｉ固体核燃料装置的四组瞬态基准实验结果进行验证.程序计算的裂变率和Ｇｏｄｉｖａ￣Ｉ装置１㊁２号实验数据的对比见图２~图５ꎬ四组实验的裂变率峰值对比见表１.表１㊀裂变率计算结果与基准题对比Ｔａｂｌｅ１㊀Ｆｉｓｓｉｏｎｒａｔｅｃａｌｃｕｌａｔｉｏｎｒｅｓｕｌｔｃｏｍｐａｒｅｄｗｉｔｈｂｅｎｃｈｍａｒｋ实验号系统周期/μｓ基准实验裂变率峰值ｆｉｓｓｉｏｎｓ/ｓ计算的裂变率峰值ｆｉｓｓｉｏｎｓ/ｓ相对误差/％１２９.５２.６７ˑ１０１９２.７１ˑ１０１９１.５０２９０２.５４ˑ１０１８２.６２ˑ１０１８３.１５３１６０１.０７ˑ１０１８０.９５ˑ１０１８－１１.２４３２０２.６３ˑ１０１７２.６２ˑ１０１７－０.３８从图２~图５可以看到ꎬ当正反应性引入时ꎬ功率迅速上升ꎬ随着热量的积累ꎬ温度上升ꎬ固体开始膨胀ꎬ引入负的反应性反馈ꎬ使功率上升速度减缓ꎬ产生功率峰.从表１的结果可以看到ꎬＧＥＴＡＣ￣２.０程序在模拟Ｇｏｄｉｖａ￣Ｉ四组基准实验时ꎬ裂变率峰值的最大误差在１２％以内ꎬ其正确性可以得到保证.３㊀溶液核燃料系统源项计算当溶液核燃料系统内因为外界环境因素的影响而进入超临界状态之后ꎬ根据法国ＳＩＬＥＮＥ溶液核燃料装置的瞬态实验结果ꎬ其系统功率的变化可以分为三个阶段:图２㊀１号实验裂变率㊁反应性随时间的变化Ｆｉｇ.２㊀Ｆｉｓｓｉｏｎｒａｔｅａｎｄｒｅａｃｔｉｖｉｔｙｃｈａｎｇｅｗｉｔｈｔｉｍｅｉｎｔｈｅ１ｓｔｅｘｐｅｒｉｍｅｎｔ图３㊀１号实验温度随时间的变化Ｆｉｇ.３㊀Ｔｅｍｐｅｒａｔｕｒｅｃｈａｎｇｅｗｉｔｈｔｉｍｅｉｎｔｈｅ１ｓｔｅｘｐｅｒｉｍｅｎｔ图４㊀２号实验裂变率㊁反应性随时间的变化Ｆｉｇ.４㊀Ｆｉｓｓｉｏｎｒａｔｅａｎｄｒｅａｃｔｉｖｉｔｙｃｈａｎｇｅｗｉｔｈｔｉｍｅｉｎｔｈｅ２ｎｄｅｘｐｅｒｉｍｅｎｔ３㊀㊀㊀南华大学学报(自然科学版)２０１８年６月图５㊀２号实验温度随时间的变化Ｆｉｇ.５㊀Ｔｅｍｐｅｒａｔｕｒｅｃｈａｎｇｅｗｉｔｈｔｉｍｅｉｎｔｈｅ２ｎｄｅｘｐｅｒｉｍｅｎｔ１)首次瞬爆阶段在引入反应性大于１Ɣ时ꎬ系统进入瞬发超临界状态.系统功率在瞬间急剧上升ꎬ而后因为温度上升和辐照裂解气体的产生引起负反应性反馈ꎬ系统功率逐渐下降ꎬ产生第一个瞬爆功率峰.当引入的反应性小于１Ɣꎬ即缓发超临界时ꎬ第一个功率峰的来临就会迟缓很多.２)功率震荡阶段在产生首次功率峰之后ꎬ由于产生的辐照裂解气体不断逸出水面ꎬ空泡引起的反应性反馈效应减弱ꎬ会使得系统进入重返临界的状态.随后功率继续增加ꎬ空泡反馈效应有增强ꎬ产生新的功率峰.如此循环往复ꎬ系统就进入了功率震荡阶段.３)功率下降阶段由于外界的人为干预引入负反应性或者燃料温度上升至沸点后ꎬ液体沸腾喷出容器或慢化剂蒸干ꎬ导致系统进入次临界状态ꎬ从而引起功率下降ꎬ事故终止.为了对整个事故过程进行模拟ꎬ需要建立相应的中子物理㊁传热和辐照裂解气泡模型.其中ꎬ中子物理和传热部分的计算采用与２.１㊁２.２节中相同的方程ꎬ辐照裂解气泡模型见３.１节.３.１㊀辐照裂解气泡模型临界事故过程中会产生大量的裂变碎片ꎬ加之中子㊁伽马与水之间的相互作用ꎬ使水发生裂解ꎬ进而产生辐照裂解气体ꎬ其主要成分为氢气和氧气.辐照裂解气体的浓度达到溶液的饱和浓度之后ꎬ将会在溶液内形成气泡ꎬ影响溶液的体积和密度.本文采用文献[５]中的辐照裂解气泡产生㊁迁移平衡方程进行空泡体积的求解ꎬ见式(６).∂ＶＢ(ｔꎬｚ)∂ｔ＝(１＋１ζ)Ｇ(Ｈ２)Ｐ(ｔꎬｚ)ＲＴ(Ｘｗ＋ρｇｈ(ｚ)＋２αｒ)－ｖ(ｚ)∂ＶＢ(ｔꎬｚ)∂ｚ(６)式中ꎬＶＢ表示空泡体积ꎻｔ表示时间ꎻｚ表示轴向位置ꎻ１/ζ表示其他气体与Ｈ２产量的比值ꎻＧ(Ｈ２)表示产氢额ꎻＰ(ｔꎬｚ)表示ｔ时刻ｚ位置的功率密度ꎻＲ表示气体常数８.３１４５ꎻＴ表示气泡温度ꎻＸｗ表示环境压强ꎻρ表示燃料密度ꎻｇ表示重力加速度ꎻｈ(ｚ)表示气泡距离溶液表面的高度ꎻα表示液体的表面张力系数ꎻｒ表示气泡半径ꎻｖ(ｚ)表示气泡的轴向运动速度.３.２㊀溶液基准实验验证采用法国ＳＩＬＥＮＥ溶液核燃料实验装置的Ｓ１￣３００~Ｓ３￣３００共三组瞬态实验结果对ＧＥＴＡＣ￣２.０进行验证.程序计算的裂变率和ＳＩＬＥＮＥ装置Ｓ２￣３００㊁Ｓ３￣３００实验数据的对比见图６㊁图７ꎬ温度㊁气泡份额随时间的变化见图８㊁图９ꎬ程序计算的裂变率峰值与三组实验数据之间的对比见表２.图６㊀Ｓ２￣３００实验裂变率㊁反应性随时间的变化Ｆｉｇ.６㊀ＦｉｓｓｉｏｎｒａｔｅａｎｄｒｅａｃｔｉｖｉｔｙｃｈａｎｇｅｗｉｔｈｔｉｍｅｉｎＳ２￣３００ｅｘｐｅｒｉｍｅｎｔ图７㊀Ｓ３￣３００实验裂变率㊁反应性随时间的变化Ｆｉｇ.７㊀ＦｉｓｓｉｏｎｒａｔｅａｎｄｒｅａｃｔｉｖｉｔｙｃｈａｎｇｅｗｉｔｈｔｉｍｅｉｎＳ３￣３００ｅｘｐｅｒｉｍｅｎｔ４第３２卷第３期朱庆福等:核燃料循环系统临界事故源项计算程序ＧＥＴＡＣ￣２.０开发图８㊀Ｓ２￣３００实验温度㊁气泡份额随时间的变化Ｆｉｇ.８㊀ＴｅｍｐｅｒａｔｕｒｅａｎｄｖｏｉｄｓｈａｒｅｃｈａｎｇｅｗｉｔｈｔｉｍｅｉｎＳ２￣３００ｅｘｐｅｒｉｍｅｎｔ图９㊀Ｓ３￣３００实验温度㊁气泡份额随时间的变化Ｆｉｇ.９㊀ＴｅｍｐｅｒａｔｕｒｅａｎｄｖｏｉｄｓｈａｒｅｃｈａｎｇｅｗｉｔｈｔｉｍｅｉｎＳ３￣３００ｅｘｐｅｒｉｍｅｎｔ表２㊀裂变率计算结果与基准题对比Ｔａｂｌｅ２㊀Ｆｉｓｓｉｏｎｒａｔｅｃａｌｃｕｌａｔｉｏｎｒｅｓｕｌｔｃｏｍｐａｒｅｄｗｉｔｈｂｅｎｃｈｍａｒｋ数据来源裂变率峰值/ｓＳ１￣３００Ｓ２￣３００Ｓ３￣３００实验１.３０ˑ１０１５１.７０ˑ１０１６１.１０ˑ１０１９ＧＥＴＡＣ￣２.０１.４４ˑ１０１５１.７９ˑ１０１６１.０９ˑ１０１９计算的相对偏差/％１１.１５.６－０.９由图６㊁图８可以看出ꎬ在Ｓ２￣３００实验初期(前６ｓ)ꎬ温度反应性反馈起主要作用ꎻ之后辐照裂解气体迅速产生ꎬ使反应性下降.在１０ｓ左右ꎬ由于此时功率处于下降状态ꎬ导致辐照裂解气体的产生量减小ꎬ系统的反应性开始上升ꎬ功率的下降速度进一步变缓.Ｓ３￣３００实验的变化结果与Ｓ２￣３００实验类似.由表２可以看到ꎬＧＥＴＡＣ￣２.０程序计算的裂变率峰值与实验数据之间的误差最大为１１.１％ꎬ其正确性可以得到保证.４㊀粉末核燃料系统源项计算粉末核燃料在没有慢化剂的情况下ꎬ其临界质量很大ꎬ不易发生临界事故.在特殊情况下ꎬ如果粉末系统受潮或有慢化剂渗入ꎬ则容易发生临界事故[６].由于慢化剂的进入ꎬ会在粉末系统内引入正反应性ꎬ引起系统功率的迅速增加ꎬ进而整个系统温度上升ꎬ体积膨胀ꎬ慢化剂的密度减少ꎬ慢化能力减弱.一般而言ꎬ湿粉末中的含水量较少ꎬ整个系统处于欠慢化的状态ꎬ故而当慢化剂的密度减少时ꎬ相当于向系统引入了负反应性反馈.负的温度反馈是抑制湿粉末系统瞬态事故工况下功率上升的第一种方式.由于粉末系统核燃料的水分量较少ꎬ且在处于粉末间隙的水难以产生稳定的辐照裂解气泡ꎬ即使产生气泡也会迅速弥散开ꎬ故而可以忽略辐照裂解气体的产生.随着温度的持续上升ꎬ粉末系统内会产生沸腾.系统中的水的温度将保持在其对应压力下的沸点不再变化.此后ꎬ沸腾传热导致系统内的液态慢化剂开始蒸发ꎬ引起系统慢化剂的干度增大ꎬ引入大量的负反应性反馈ꎬ并使得系统进入次临界状态.慢化剂的沸腾引起的相变反馈是抑制湿粉末系统瞬态事故工况下功率上升的第二种方式.假若之后再没有新的液态慢化剂流入ꎬ系统的功率将保持次临界状态ꎬ直至功率下降至极低水平.对粉末核燃料事故工况的模拟ꎬ中子物理部分的计算采用与２.１节中相同的方程ꎬ传热部分计算包括固㊁液两相传热和沸腾相变计算两个部分.４.１㊀固㊁液两项传热模型湿粉末系统内存在固㊁液两相ꎬ其精确的传热计算应该求解固液两相Ｎａｖｉｅｒ￣Ｓｔｏｋｅｓ方程ꎬ但因为含水量不大ꎬ固体粉末周围的水将迅速被加热ꎬ故在传热计算时ꎬ可以将系统看作单相均匀系统进而等效为单相非稳态传热.不同于一般单相传热的地方是ꎬ这个等效单相均匀系统的热工参数将必须单独计算.５㊀㊀㊀南华大学学报(自然科学版)２０１８年６月传热方程中的两个重要热工参数是比热容Ｃｐ和导热系数ｋꎬ这两个参数在很大程度上决定了传热方程的计算精度.对于湿粉末系统而言ꎬ混合体系的比热容和导热系数无法通过查表得到ꎬ故而这里采用化学工艺中经典的ＥＭＴ￣Ｌａｎｄａｕｅｒ模型进行求解[７￣８].ＥＭＴ￣Ｌａｎｄａｕｅｒ模型不同于传统简单热工中的热阻并联(热流从上到下流经不同的材料)㊁串联(热流从上到下经过每一层材料)模型ꎬ是一种适用于两种材料随机分布ꎬ每一相之间既不连续又不分散ꎬ每一中组分能否形成导热路径ꎬ取决于组分的量.ＥＭＴ￣Ｌａｎｄａｕｅｒ模型等效导热系数ｋｅ计算公式为:ðＮｉ＝１νｉｋｉ－ｋｅｋｉ＋２ｋｅ＝０(７)㊀㊀对于固液两相流ꎬ其等效导热系数可由式(８)给出:ｋｅ＝３ε－１()ｋｆ＋３１－ε()－１[]ｋｓ４＋３ε－１()ｋｆ＋３１－ε()－１()ｋｓ[]２＋８ｋｆｋｓ４(８)其中ꎬｋε为有效导热系数ꎬε为固相体积份额ꎬｋｆ㊁ｋｓ分别为液相㊁固相导热系数.等效比热容可由能量守恒方程得:Ｃｐｅ＝ｗｆＣｐｆ＋ｗｓＣｐｓ(９)其中ꎬＣｐｅ为等效比热容ꎬｗｆ㊁ｗｓ分别为液相㊁固相质量分数ꎬＣｐｆ㊁Ｃｐｓ分别为液相㊁固相比热容.４.２㊀沸腾相变传热模型当粉末系统内热量不断积累ꎬ使水达到其沸点后ꎬ系统的传热将由非稳态热传导转变为沸腾传热.此时ꎬ系统的温度将不再增加(不再遵守非稳态传热方程)ꎬ而是维持在水的沸点不变ꎬ热量由水的气化潜热带走ꎬ直至慢化剂蒸干ꎬ系统进入次临界状态.在粉末中的水进入沸腾之后ꎬ定压下的饱和水随着热量的继续增加而逐渐变为湿饱和蒸汽而从系统表面蒸发.由于沸腾状态下系统温度不再变化ꎬ故系统内的热量将全部由气化潜热带走ꎬ系统的干度决定了蒸发的水的量ꎬ也决定了沸腾相变引起的反应性反馈.干度可由式(１０)计算得到.ｘ＝ʏｔｘｔ０ʏＶＰｒңꎬｔ()ｄＶｄｔｈᵡ－ｈᶄ(１０)４.３㊀程序计算结果由于目前国际上尚无粉末核燃料的瞬态实验结果ꎬ故暂时无法进行程序的实验验证ꎬ这里仅展示程序的计算结果.０.２Ɣ㊁１.２Ɣ反应性引入情况下粉末系统的裂变率㊁温度随时间的变化见图１０㊁图１１ꎬ反应性的变化见图１２㊁图１３.从图１０㊁图１２中可以看到ꎬ由于０.２Ɣ阶跃正反应性的引入ꎬ系统裂变率开始增加ꎬ随后慢化剂温度上升ꎬ引入了负的反应性反馈ꎬ导致反应性的减小.当反应性减少到０时ꎬ系统的裂变率达到最大值.由于引入的反应性不大ꎬ故而系统并未产生泡核沸腾.当系统的产热与外界环境对流换热之间达到平衡时ꎬ系统的功率将稳定下来ꎬ裂变率不再变化.整个过程中ꎬ系统的反应性反馈仅由温度反馈提供.图１０㊀０.２Ɣ初始反应性引入下裂变率㊁温度随时间的变化Ｆｉｇ.１０㊀Ｆｉｓｓｉｏｎｒａｔｅａｎｄｔｅｍｐｅｒａｔｕｒｅｃｈａｎｇｅｗｉｔｈｔｉｍｅｗｉｔｈ０.２Ɣｉｎｉｔｉａｌｒｅａｃｔｉｖｉｔｙｉｎｓｅｒｔｉｏｎ图１１㊀１.２Ɣ初始反应性引入下裂变率㊁温度随时间的变化Ｆｉｇ.１１㊀Ｆｉｓｓｉｏｎｒａｔｅａｎｄｔｅｍｐｅｒａｔｕｒｅｃｈａｎｇｅｗｉｔｈｔｉｍｅｗｉｔｈ１.２Ɣｉｎｉｔｉａｌｒｅａｃｔｉｖｉｔｙｉｎｓｅｒｔｉｏｎ６第３２卷第３期朱庆福等:核燃料循环系统临界事故源项计算程序ＧＥＴＡＣ￣２.０开发图１２㊀０.２Ɣ初始反应性引入下总反应性随时间的变化Ｆｉｇ.１２㊀Ｒｅａｃｔｉｖｉｔｙｃｈａｎｇｅｏｖｅｒｔｉｍｅｗｉｔｈ０.２Ɣｉｎｉｔｉａｌｒｅａｃｔｉｖｉｔｙｉｎｓｅｒｔｉｏｎ图１３㊀１.２Ɣ初始反应性引入下总反应性随时间的变化Ｆｉｇ.１３㊀Ｒｅａｃｔｉｖｉｔｙｃｈａｎｇｅｏｖｅｒｔｉｍｅｗｉｔｈ１.２Ɣｉｎｉｔｉａｌｒｅａｃｔｉｖｉｔｙｉｎｓｅｒｔｉｏｎ从图１１㊁图１３中可以看到ꎬ１.２Ɣ阶跃正反应性的引入ꎬ使系统达到了瞬发临界状态ꎬ裂变率在０.６ｓ内就上涨了１０１２倍ꎬ达到了非常高的水平.大量的裂变释热导致系统温度也急速上涨ꎬ很快便达到了慢化剂沸点.沸腾的产生使得慢化剂开始大量蒸发ꎬ干度增加ꎬ进而引入很大的负反应性反馈ꎬ系统进入深次临界状态.㊀㊀事故过程中ꎬ系统的反应性反馈由慢化剂的温度反馈和达到慢化剂沸点后的沸腾反馈组成.在温度到达慢化剂沸点之前ꎬ反应性的抑制完全依赖慢化剂的温度反馈ꎻ当温度上升至慢化剂的沸点之后ꎬ温度反馈不再起作用ꎬ沸腾反馈起到抑制系统反应性的主导作用.５㊀结㊀论本文采用中子物理 热工水力耦合的方法ꎬ对固体㊁溶液㊁粉末态的核燃料建立了各自的耦合模型ꎬ开发了核燃料循环系统超临界事故源项计算程序ＧＥＴＡＣ￣２.０ꎬ并应用国际上公开的实验数据进行了验证(粉末系统除外).结果表明ＧＥＴＡＣ￣２.０程序计算的裂变率变化与实验结果符合良好ꎬ计算的裂变率峰值与实验数据之间的最大偏差在１２％以内ꎬ其正确性可以得到保证.参考文献:[１]张驰ꎬ周琦ꎬ朱庆福ꎬ等.金属核燃料系统瞬态特性分析研究[Ｊ].原子能科学技术ꎬ２０１６ꎬ５０(１２):２１７０￣２１７４.[２]贺仁辅.快中子临界装置和脉冲堆实验物理[Ｍ].北京:国防工业出版社ꎬ２０１２.[３]陶文铨.数值传热学[Ｍ].西安:西安交通大学出版社ꎬ１９８８.[４]ＡＵＦＩＥＲＯＭꎬＦＩＯＲＩＮＡＣꎬＬＡＵＲＥＡＵＡꎬｅｔａｌ.Ｓｅｒｐｅｎｔ￣￣ＯｐｅｎＦＯＡＭｃｏｕｐｌｉｎｇｉｎｔｒａｎｓｉｅｎｔｍｏｄｅ:ｏｆａＧｏｄｉｖａｐｒｏｍｐｔｃｒｉｔｉｃａｌｂｕｒｓｔ[Ｊ].Ｊｏｕｒｎａｌｏｆｐｅｒｉｏｄｏｎｔａｌｒｅｓｅａｒｃｈꎬ２０１５ꎬ３７(２):１５４￣１６０.[５]汪量子.溶液堆的蒙特卡罗方法物理计算模型及特性研究[Ｄ].北京:清华大学ꎬ２０１１.[６]李喆ꎬ刘开武.核临界安全手册[Ｍ].北京:原子能出版社ꎬ２００３.[７]付文强ꎬ高辉ꎬ薛征欣ꎬ等.多孔材料有效导热系数的实验和模型研究[Ｊ].中国测试ꎬ２０１６ꎬ４２(５):１２４￣１３０.[８]ＳＵＮＤÉＮＢꎬＹＵＡＮＪ.Ｅｖａｌｕａｔｉｏｎｏｆｍｏｄｅｌｓｏｆｔｈｅｅｆｆｅｃｔｉｖｅｔｈｅｒｍａｌｃｏｎｄｕｃｔｉｖｉｔｙｏｆｐｏｒｏｕｓｍａｔｅｒｉａｌｓｒｅｌｅｖａｎｔｔｏｆｕｅｌｃｅｌｌｅｌｅｃｔｒｏｄｅｓ[Ｊ].Ｉｎｔｅｒｎａｔｉｏｎａｌｊｏｕｒｎａｌｏｆｃｏｍｐｕｔａｔｉｏｎａｌｍｅｔｈｏｄｓａｎｄｅｘｐｅｒｉｍｅｎｔａｌｍｅａｓｕｒｅｍｅｎｔｓꎬ２０１３ꎬ１(４):４４０￣４５４.(责任编辑:扶文静)７。