一种六角形组件几何堆芯的时空动力学求解方法
【专利摘要】本发明公开了一种六角形组件几何堆芯的时空动力学求解方法，所述方法包括：采用对角线隐式龙格库塔方法分离时间项，将时空动力学方程转化为固定源问题，基于保角变换思想的格林函数函数节块法求解所述固定源问题，实现了六角形组件几何堆芯时空动力学求解，本方法能够用于六角形组件几何堆芯瞬态事故物理分析，为六角形组件几何堆芯的相关的理论研究和工程设计奠定了技术基础。
【专利说明】
-种六角形组件几何堆巧的时空动力学求解方法
技术领域
[0001] 本发明设及核反应堆堆忍设计领域，尤其设及一种六角形组件几何堆忍的时空动 力学求解方法。
【背景技术】
[0002] 在核能利用中，核安全一直是一个头等重要的问题，尤其是日本福岛事故之后，国 内外对反应堆安全重视程度上升到了一个前所未有的高度;在核反应堆安全分析和负荷跟 踪过程中，堆忍功率或中子通量随时间的瞬态变化至关重要，运就需要求解反应堆的中子 动力学方程；W往的点堆或一维动力学模型无法精确模拟功率分布，事故分析的精度已经 不能满足工程要求，必须求解=维时空动力学方程。
[0003] 另一方面，六角形组件几何稠密堆忍是一种先进的反应堆堆忍设计技术，具有高 功率密度和长寿期的特点；目前广泛应用于俄罗斯、乌克兰、保加利亚和匈牙利等多个国 家；中国田湾核电站采用的俄AES-91也是在VVER-1000/V320基础上改进的，除PWR夕h六角 形组件几何稠密堆忍在高溫气冷堆、液态金属冷却快中子增殖堆W及超临界水冷反应堆和 新概念的行波堆中也将得到应用；同时六角形排列的堆忍特别适合于军用核动力装置；因 此，六角形组件几何堆忍的设计研发、事故分析是非常有价值的;但是由于六角形几何特殊 性，矩形的堆忍物理时空动力学软件不能直接应用到六角形堆忍，六角形堆忍物理时空动 力学软件研发是必须和急迫的。
[0004] 目前国内外处理时空动力学空间离散方法仍采用粗网节块法，有多项式展开法、 解析节块法、半解析节块法等；1979年，R.D丄awrence提出的格林函数节块，采用第S类边 界条件格林函数作为积分核，求解出射流响应矩阵;90年代，清华大学胡永明、赵险峰采用 第二类边界条件格林函数作为积分核；由于净流响应矩阵为=对角，可直接求解无需内迭 代过程，计算效率较高;格林函数节块法能够带不连续因子计算，具有模型推导近似少，不 受能群限制，精度高等优点;但是只能应用于矩形组件几何堆忍计算;时空动力学时间离散 方法有向后欧拉格式、刚性限制、龙格库塔方法等;清华大学2012年在向后欧拉格式基础上 提出了对角线隐式龙格库塔方法，具有较高计算精度和效率，但也只适用于矩形堆忍时空 动力学求解。
[0005] 综上所述，本申请发明人在实现本申请实施例中发明技术方案的过程中，发现上 述技术至少存在如下技术问题：
[0006] 在现有技术中，现有的时空动力学求解方法存在无法适用于六角形组件几何堆忍 时空动力学求解的技术问题。
【发明内容】
[0007] 本发明提供了一种六角形组件几何堆忍的时空动力学求解方法，解决了现有的时 空动力学求解方法存在无法适用于六角形组件几何堆忍时空动力学求解的技术问题，实现 了六角形组件几何堆忍时空动力学求解，本方法能够用于六角形组件几何堆忍瞬态事故物 理分析，为六角形组件几何堆忍的相关的理论研究和工程设计奠定了技术基础。
[0008] 为解决上述技术问题，本申请实施例提供了一种六角形组件几何堆忍的时空动力 学求解方法，所述方法包括：
[0009] 采用对角线隐式龙格库塔方法分离时间项，将时空动力学方程转化为固定源问 题，基于保角变换思想的六角形格林函数函数节块法求解所述固定源问题。
[0010] 进一步的，所述采用对角线隐式龙格库塔方法分离时间项，将时空动力学方程转 化为固定源问题，具体包括：
[0011] 六角形中子时空动力学方程组如下式，
[0012]
(1)
[001引式中，下标g代表能群，能群数为G;下标i代表缓发先驱核组，组数为ND;參(/V)为g 群中子通量密度，单位为n/(cm2Xs);c,护，O为第i组先驱核浓度，单位为l/cm3;vg为g群中子 速度，单位为cm/s; Dg为g群扩散系数，单位为cm; X r,g和X f,g分别为g群移出截面和裂变截 面，5： 气为g/群到g群的散射截面，单位为cnfi; V为每次裂变释放的中子数;Xg为瞬发中子 裂变谱;Xg, 1为缓发中子谱分额;、为先驱核衰变常数，单位为; 01为缓发中子份额；
[0014] 记
it为当前待求解时刻，to为上一时刻， A t = t-to，yi = y(t)，yo = y(to)，f (t，y)代表式（1)中非时间偏导数项的总和，将式（1)简记 为如下形式：
[0015]
(ii
[0016] 其中，S为龙格库塔级数，由aij组成的S阶方阵A称为龙格库塔矩阵，由bj构成的向 量b称为龙格库塔权，由Cj构成的向量C称为龙格库塔节点，用服(s，p)表示一个S级P阶精度 的龙格库塔格式；
[0017] 对式(2)含时微分方程，对角线隐式龙格库塔方法的离散形式为：
[001引 獄
[0019]
[0020] (4)
[0021] 设?g,-仍为第i级第g群通量分布，为第i级第n组先驱核浓度分布；
[0022] (5)
[002；3] 脚
[0024] 将DIRK格式具体应用于式（I)，首先离散先驱核浓度方程，第i级先驱核浓度方程 转化为：
[0025]
(7)
[0026] 中子通量的时间偏导数采用相同的格式，代入第i级通量方程得到最终求解的固 定源问题，如下式：
[0027]
[002引
[0029]
[0030]
[0031] 固定源方程采用六角形节块格林函数方法求解，解得第i级Cl后，由式(3)得到f(to +Cl A t点)用于后续级九"1的计算；
[0032] 综合考虑计算精度、稳定性及计算量，采用2级2阶精度的DI服格式：
[00削
。巧
[0034] 进一步的，所述基于保角变换思想的格林函数函数节块法求解所述固定源问题， 具体包括：
[0035] 六角形节块k内=维分群扩散方程为：
节块的长和宽，hk为当前节块的轴向高度，g(U，V'）为保角变换因子；
[0036； (12)
[0037；
[0038； (13)
[0039； ^角形节块保角变换后矩形
[0041] (14)
[0040] 对式(13)沿V'，z两个方向积分，得到横向积分方程：
[0042]
[0043] (15) 闺 (16)
[0045] (17)
[0046] (峭
[0047] (19)
[004引 （20)
[0049]引入第二类边界条件格林函数巧>為)，经过格林函数节块求解公式推导得到偏 中子通量方程的积分解为：
[0化3] 其中在相邻节块的界面上，满足非均匀通量连续，引入通量不连续因子，/:<+为k节块U方向正端点的g群不连续因子;指k节块的U方向+1的节块的负端点不连续因子；
[00 加； (2:1)
[0化1;
[0052； 泛）
[i ；0化41
[i
[0057]在节块内体积分得到六角形节块的中子平衡方程：
(26)
[0化引
[0059] 其中與、巧。为六角形节块k第g群，S方向左（右)面的平均净流，S = U, v，w;g = 1，--G;惡为六角形节块k第g群的平均通量，Rk，A%六角形节块的径向面积及体积；g(〇 和巧輪)分别为节块体积平均源项和固定源项；
[0060] (27)
[0061]
[0062] 本申请实施例中提供的一个或多个技术方案，至少具有如下技术效果或优点：
[0063] 提出了一种六角形组件几何堆忍的时空动力学求解方法具体方法包括:采用对角 线隐式龙格库塔方法分离时间项，将时空动力学方程转化为固定源问题，基于保角变换思 想的六角形格林函数函数节块法求解所述固定源问题的技术方案，即空间离散采用格林函 数节块法，时间离散采用对角线隐式龙格库塔方法，且通过基准题的计算验证了方法的正 确性，所W，有效解决了现有的时空动力学求解方法存在无法适用于六角形组件几何堆忍 时空动力学求解的技术问题，进而实现了本方法六角形组件几何堆忍时空动力学求解，本 方法能够用于六角形组件几何堆忍瞬态事故物理分析，为六角形组件几何堆忍的相关的理 论研究和工程设计奠定了技术基础的技术效果。
【附图说明】
[0064] 图1是本申请实施例一中龙格库塔格式示意图；
[0065] 图2是本申请实施例一中国际参考程序计算的堆忍相对功率随时间变化趋势示意 图；
[0066] 图3是本申请实施例一中基于本方法研制程序计算的堆忍相对功率随时间变化趋 势示意图；
[0067] 图4是本申请实施例一中15 %扰动下中屯、节块热群通量随时间变化示意图；
[0068] 图5是本申请实施例一中15.1 %扰动下中屯、节块热群通量随时间变化示意图。
【具体实施方式】
[0069] 本发明提供了一种六角形组件几何堆忍的时空动力学求解方法，解决了现有的时 空动力学求解方法存在无法适用于六角形组件几何堆忍时空动力学求解的技术问题，实现 了六角形组件几何堆忍时空动力学求解，能够用于六角形组件几何堆忍瞬态事故物理分 析，为六角形组件几何堆忍的相关的理论研究和工程设计奠定了技术基础。
[0070] 本申请实施中的技术方案为解决上述技术问题。总体思路如下：
[0071] 提出了一种六角形组件几何堆忍的时空动力学求解方法，具体方法包括:采用对 角线隐式龙格库塔方法分离时间项，将时空动力学方程转化为固定源问题，基于保角变换 思想的六角形格林函数函数节块法求解所述固定源问题的技术方案，即空间离散采用格林 函数节块法，时间离散采用对角线隐式龙格库塔方法，且通过基准题的计算验证了方法的 正确性，所W，有效解决了现有的时空动力学求解方法存在无法适用于六角形组件几何堆 忍时空动力学求解的技术问题，进而实现了六角形组件几何堆忍时空动力学求解，本方法 能够用于六角形组件几何堆忍瞬态事故物理分析，为六角形组件几何堆忍的相关的理论研 究和工程设计奠定了技术基础。
[0072] 为了更好的理解上述技术方案，下面将结合说明书附图W及具体的实施方式对上 述技术方案进行详细的说明。
[0073] 实施例一：
[0074] 在实施例一中，提供了一种六角形组件几何堆忍的时空动力学求解方法，请参考 图1-图5,所述方法包括：
[0075] 采用对角线隐式龙格库塔方法分离时间项，将时空动力学方程转化为固定源问 题，基于保角变换思想的六角形格林函数函数节块法求解所述固定源问题。
[0076] 本发明提出了适用于一种新型六角形组件几何堆忍的时空动力学数值计算方法， 包括空间离散和时间离散两方面：
[0077] 1.时间离散
[0078] 六角形中子时空动力学方程组如下式，
[0079]
[0080] 式中，下标g代表能群，能群数为G;下标i代表缓发先驱核组，组数为ND;在护，O为g 群中子通量密度，n/(cm2Xs);c,巧。为第i组先驱核浓度，l/cm3;vg为g群中子速度，cm/s;Dg 为邑群扩散系数，cm; Sr,g和Xf,g分别为g群移出截面和裂变截面，Xs,-g为g/群到g群的散 射截面，cnfi;v为每次裂变释放的中子数;Xg为瞬发中子裂变谱;Xg,1为缓发中子谱分额;Al 为先驱核衰变常数，S^l; 01为缓发中子份额。
[0081] 为便于讨论，巧
，t为当前待求解时刻，to 为上一时刻，A t = t-to，yi = y(t)，yo = y(to)，f(t，y)代表方程组（1)中非时间偏导数项的 总和，将方程组（1)简记为如下形式：
[0082]
(30)
[0083] 龙格库塔方法如图1所示，S为龙格库塔级数(stage)，由曰^组成的S阶方阵A称为龙 格库塔矩阵，由bj构成的向量b称为龙格库塔权，由Cj构成的向量C称为龙格库塔节点，通常 用RK(s，p)表示一个S级P阶精度的龙格库塔格式。
[0084] 对形如式(2)的含时微分方程，对角线隐式龙格库塔化I服:矩阵A的上S角元素为 零)方法的离散形式为：
[0085] (31)
[0086]
[0087] (32)
[0088] 设护)为第i级第g群通量分布，C。,,护)为第i级第n组先驱核浓度分布（因为第i 级已经代表时间变量，所W变量中只写空间分布量F，不写时间变量）。
[0089] (33)
[0090] (34)
[0091] 将DIRK格式具体应用于瞬态方程组（1)。首先离散先驱核浓度方程，第i级先驱核 浓度方程转化为：
[0092]
(35)
[0093] 中子通量的时间偏导数采用相同的格式，代入第i级通量方程得到最终求解的固 定源问题，如下式。
[0098] 固定源方程采用节块格林函数方法求解。解得第i级Cl后，由式(3)得到f(to+Ci A t，Ci)用于后续级而.,,1的计算。
[0094]
[0095]
[0096]
[0097]
[0099] 综合考虑计算精度、稳定性及计算量，采用2级2阶精度的的DI服格式：
[0100] 削:
[0101]
[0102]
[0103] (40)
[0104]
[0105] (4U
[0106] 六角形节块保角变换后矩形 节块的长和宽，hk为当前节块的轴向高度，g(U，V'）为保角变换因子。
[0107] 对上式沿V'，z两个方向积分，得到横向积分方程
[010引 （侣）
[0109]
[0110] (43)
[0111] (44)
[0112] (巧
[0113] (46)
[0114] (47)
[0115] (48)
[0116] 引入第二类边界条件格林函数(?, (M，u。），经过格林函数节块求解公式推导得到偏 中子通量方程的积分解为：
[0117] (49)
[011 引
[0119]
[0120] 其中在相邻节块的界面上，满足非均匀通量连续，引入通量不连续因子，为k节 块U方向正端点的g群不连续因子。/gt"指k节块的U方向+1的节块的负端点不连续因子。 「01911
[0124] 在节块内体积分得到六角形节块的中子平衡方括：
[0126] 其中巧，，巧C为六角形节块k第g群，S方向左（右)面的平均净流，S = U, v，w;g = 1，…G;g为六角形节块k第g群的平均通量，Rk，Ak为六角形节块的径向面积及体积。)
[0125]
[0127] W) 和巧的)分别为节块体积平均源项和固定源项。
[012 引
[0129] 偏中子通量求解公式，净中子流禪合方程W及节块平衡方程构成完备求解方程 组，采用剩余权重法展开，并通过源迭代法进行求解。
[0130] 本发明提出了一种新型六角形时空动力学求解方法，下面通过计算基准题，与国 际相应程序计算结果对比验证方法的正确性。
[0131] 1.二维基准题
[0132] 该基准题类似于一个重水堆的试验装置，该堆忍的装载模式使得中子通量的峰值 位置出现在堆忍燃料组件的边缘，瞬态过程模拟了 1号材料的热群吸收截面在0.2秒内减少 4.5%，该扰动相当于向堆内引入了75.1分的正反应性。
[0133] 图2和图3分别给出了在瞬态开始3秒内，国际相应程序和基于本方法编制程序计 算的堆忍相对功率随时间的变化趋势图，DIF3D-K是采用目方法离散时间项的节块法程序， FX2-TH是细网差分程序，FEM肥X-K是使用全隐式向后差分格式和时间积分方法离散时间导 数项的解析节块法；由该图可W看出，本方法的数值结果和各程序符合较好。
[0134] 2. S维基准题
[0135]该基准题有20圈组件，组件对边距17.78cm，轴向高500cm，瞬态过程模拟了最中屯、 组件及其相邻的一圈组件的热群吸收截面在to时刻阶跃变化了 15%和15.1%，引起了堆忍 中屯、热群中子通量显著增长；图4、图5给出了瞬态开始5秒内，基于本方法编制的程序和 肥NKO程序计算结果比较。肥NKO程序是犯M时间离散格式的六角形节块展开法瞬态程序；由 该图可W看出，本方法的数值结果和肥NKO程序符合较好。
[0136] 六角形时空动力学瞬态方程采用剩余权重法求解，将偏中子通量密度和横向泄漏 项在节块内用二阶Legerulre正交多项式Pn-I(U) ,n = l ,2,3展开，有：
[0137]
(57)
[013引其中儀，4?如=1，2,3)为节块4第邑群偏中子通量密度和横向泄露项的展开系数。
[0139] 对于=维问题，径向横向积分对应的横向泄漏分为径向横向泄漏I置(U)和轴向横 向泄漏 <二(U)两部分，即：
[0140]
[0141]
[0142]
[0143] 。
[0144] 将中子源项和泄漏项展开式代入式(23)中得到[(义)±t的求解公式，将偏通量、源 项和泄漏项展开式代入第二类边界条件格林函数积分方程(式(21))得到偏通量展开系数 求解公式。偏通量展开系数方程(式(21))、净中子流禪合方程(式(22)) W及中子平衡方程 (式(26))构成一组完备的求解公式，通过源迭代方法求解，得出节块界面平均净中子流、节 块内中子通量展开系数、节块平均通量Wkeff等物理量。
[0145] 上述本申请实施例中的技术方案，至少具有如下的技术效果或优点：
[0146] 提出了一种六角形组件几何堆忍的时空动力学求解方法具体方法包括:采用对角 线隐式龙格库塔方法分离时间项，将时空动力学方程转化为固定源问题，基于保角变换思 想的六角形格林函数函数节块法求解所述固定源问题的技术方案，即空间离散采用六角形 格林函数节块法，时间离散采用对角线隐式龙格库塔方法，且通过基准题的计算验证了方 法的正确性，所W，有效解决了现有的时空动力学求解方法存在无法适用于六角形组件几 何堆忍时空动力学求解的技术问题，进而实现了六角形组件几何堆忍时空动力学求解，本 方法能够用于六角形组件几何堆忍瞬态事故物理分析，为六角形组件几何堆忍的相关的理 论研究和工程设计奠定了技术基础。
[0147] 尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造 性概念，则可对运些实施例作出另外的变更和修改。所W，所附权利要求意欲解释为包括优 选实施例W及落入本发明范围的所有变更和修改。
[0148] 显然，本领域的技术人员可W对本发明进行各种改动和变型而不脱离本发明的精 神和范围。运样，倘若本发明的运些修改和变型属于本发明权利要求及其等同技术的范围 之内，则本发明也意图包含运些改动和变型在内。
【主权项】
1. 一种六角形组件几何堆芯的时空动力学求解方法，其特征在于，所述方法包括： 采用对角线隐式龙格库塔方法分离时间项，将时空动力学方程转化为固定源问题，基 于保角变换思想的六角形格林函数函数节块法求解所述固定源问题。2. 根据权利要求1所述的方法，其特征在于，所述采用对角线隐式龙格库塔方法分离时 间项，将六角形时空动力学方程转化为固定源问题，具体包括： 六角形中子时空动力学方程组如下式：⑴ 式中，下标g代表能群，能群数为G;下标i代表缓发先驱核组，组数为ND;么OV)为g群中 子通量密度，单位为n/(cm2Xs) ;?(?,?)为第i组先驱核浓度，单位为l/cm3;vAg群中子速 度，单位为cm/s;D gSg群扩散系数，单位为cm;Xr,g和Xfig分别为g群移出截面和裂变截面， 群至Ijg群的散射截面，单位为cnfSv为每次裂变释放的中子数;X g为瞬发中子裂 变谱;xg,:为缓发中子谱分额;M为先驱核衰变常数，单位为^T1 为缓发中子份额； 记j =[為(/V),…,九(Λ0,q(尸，.，t为当前待求解时刻，t〇为上一时亥1J，Δ t = t_to，yi = y(t)，yo = y(to)，f (t，y)代表式（1)中非时间偏导数项的总和，将式（1)简记为如 下形式：(2) 其中，s为龙格库塔级数，由aij组成的s阶方阵A称为龙格库塔矩阵，由bj构成的向量b称 为龙格库塔权，由Cj构成的向量c称为龙格库塔节点，用RK(s，p)表示一个s级p阶精度的龙 格库塔格式； 对式(2)含时微分方程，对角线隐式龙格库塔方法的离散形式为：设Φρ·Ρ)为第i级第g群通量分布，为第i级第η组先驱核浓度分布；(5) (3) (4)(6) 将DIRK格式具体应用于式（1)，首先离散先驱核浓度方程，第i级先驱核浓度方程转化 为： (7) 中子通量的时间偏导数采用相同的格式，代入第i级通量方程得到最终求解的固定源 问题，如下式：固定源方程采用六角形节块格林函数方法求解，解得第i级1后，由式(3)得到f(to+Cl AtJ1)用于后续级凡用的计算； 综合考虑计算精度、稳定性及计算量，采用2级2阶精度的DIRK格式：111)3.根据权利要求2所述的方法，其特征在于，所述基于保角变换思想的六角形格林函数 函数节块法求解所述固定源问题，具体包括： 六角形节±#k内=维分群扩散方稈为：块的长和宽，hk为当前节块的轴向高度，g(u，v'）为保角变换因子；对式(13)沿ν'，z两个方向积分，得到横向积分方程： (12) (13) <角形节块保角变换后矩形节弓丨入第二类边界条件格林函数<(?，《。），经过格林函数节块求解公式推导得到偏中子 通畺玄耜的葙分艇为.其中在相邻节块的界面上，满足非均匀通量连续，引入通量不连续因子，/i+sk节块u 方向正端点的g群不连续因子;./^1指k节块的u方向+1的节块的负端点不连续因子；在节块内体积分得到六角形节块的中子平衡方程：其中、.为六角形节块k第g群，S方向左(右)面的平均净流，s = u，v，w;g=l,…G; g为六角形节块k第g群的平均通量，Rk，Ak为六角形节块的径向面积及体积;g(〇和$仏） 分别为节块体积平均源项和固定源项；
【文档编号】G06F17/11GK106021184SQ201610298272
【公开日】2016年10月12日
【申请日】2016年5月6日
【发明人】安萍, 刘东彬, 芦韡, 李庆, 姚栋
【申请人】中国核动力研究设计院