CN108763748A - 一种用于热管堆瞬态特性分析的方法 - Google Patents

一种用于热管堆瞬态特性分析的方法 Download PDF

Info

Publication number
CN108763748A
CN108763748A CN201810523851.2A CN201810523851A CN108763748A CN 108763748 A CN108763748 A CN 108763748A CN 201810523851 A CN201810523851 A CN 201810523851A CN 108763748 A CN108763748 A CN 108763748A
Authority
CN
China
Prior art keywords
temperature
safety rod
formula
area
active region
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201810523851.2A
Other languages
English (en)
Other versions
CN108763748B (zh
Inventor
曹良志
屈伸
周生诚
吴宏春
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN201810523851.2A priority Critical patent/CN108763748B/zh
Publication of CN108763748A publication Critical patent/CN108763748A/zh
Application granted granted Critical
Publication of CN108763748B publication Critical patent/CN108763748B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/06—Power analysis or power optimisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/08—Thermal analysis or thermal optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Operations Research (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Monitoring And Testing Of Nuclear Reactors (AREA)

Abstract

一种用于热管堆瞬态特性分析的方法，1、将热管堆沿径向由内而外依次划分为安全棒区、活性区以及径向反射层区；分别列出一维导热微分方程和边界条件；解得不同功率水平下各区的温度分布；2、基于体膨胀系数是线性膨胀系数3倍的假设，列出三区的直径、高度以及密度随温度的函数关系；3、对扇环进行三角网格划分研究。基于组合几何方法将1/6热管堆分割成18个扇环，从而对1/6热管堆进行三角网格划分；4、使用SARAX程序对其输运计算。将堆芯径向功率分布提供给步骤一重新进行温度分布计算；再进行步骤二的膨胀计算；再经过步骤三，将划分好的三角形网格提供给下一时间步的输运计算；依次循环；可计算热管堆的堆芯功率水平和反应性在瞬态过程中的变化。

Description

一种用于热管堆瞬态特性分析的方法
技术领域
本发明涉及热管堆动力学计算领域，具体涉及一种用于热管堆瞬态特性分析的方法。
背景技术
空间堆可为未来太阳系内行星间的旅行和探索提供长期稳定的电源，并可为航天器提供动力。目前国际上通常采用蒙卡程序进行空间堆的物理计算和堆芯设计。但该方法无法进行空间堆的动力学以及瞬态计算且计算成本偏高。
空间堆反射层温度和膨胀反应性系数对堆芯中子动力学特性影响显著。空间堆通常采用控制鼓作为反应性控制和调节系统，控制鼓的转动将导致中子通量密度空间分布的剧烈变化。空间堆的上述特点使得传统点堆中子动力学分析方法不再适用，显式考虑反射层反应性反馈系数及其内部的中子输运过程和控制鼓的转动，建立空间堆时空中子动力学分析模型势在必行。对空间堆建立详细的系统分析模型是空间堆设计和研发的关键基础之一。由于空间堆系统与地球上运行的反应堆系统显著不同，目前广泛采用的地球运行反应堆瞬态分析模型并不适用于空间堆瞬态分析。开展空间堆时空中子动力学及系统瞬态特性研究，耦合堆芯时空中子动力学和空间堆系统各部分瞬态热工水力学模型，建立高精度的空间堆系统瞬态分析模型，满足空间堆启动、停堆、重启和事故工况下系统瞬态分析需求，将为我国空间堆系统设计和研发提供坚实的技术基础。
发明内容
为解决以上问题，本发明提供了一种空间热管反应堆瞬态多物理分析计算方法，在现有的快堆时空中子动力学程序SARAX基础上添加动态网格程序、二维热膨胀模型以及热管热工瞬态程序；能够既考虑热管堆在瞬态过程中由热管提供的热工反馈，又能考虑由于热膨胀提供的反馈，获得热管堆在多物理耦合下的精确瞬态过程。
为了达到上述目的，本发明的技术方案概括如下：
一种用于热管堆瞬特性分析的方法，包括如下步骤：
步骤1：对于以六角形方式排列的具有234根燃料热管组件以及6根控制鼓的热管堆，沿径向将其由内而外依次划分为安全棒区、活性区以及径向反射层区；安全棒区为一圆柱体，直径为安全棒绝热层的外径；活性区为一圆环柱体，内径为安全棒绝热层外径，外径为内容器壁外径；径向反射层区仍为一圆环柱体，内径为内容器壁外径，外径为外容器壁外径；各区域的高均与热管堆的高相同；其中，安全棒区沿径向满足一维导热微分方程，如公式(1)所示；活性区沿径向满足一维含内热源导热微分方程，如公式(2)所示；径向反射层区沿径向满足一维导热微分方程，如公式(3)所示。
式中：
r——距离安全棒区中心线的距离；
λs——安全棒区导热系数；
λf——活性区导热系数；
λr——径向反射层区导热系数；
Ts(r)——安全棒区径向温度；
Tf(r)——燃料活性区径向温度；
Tr(r)——径向反射层区径向温度；
——活性区的体积释热率；
安全棒区的中心区满足绝热边界条件；安全棒区与活性区边界处温度连续并且热流密度连续；活性区与径向反射层区边界处温度连续并且热流密度连续；径向反射层的外边界满足辐射边界条件；因此可得如公式(4)到(9)所示的边界条件；
Ts(r1)＝Tf(r1) 公式(5)
Tf(r2)＝Tr(r2) 公式(7)
式中：
r——距离安全棒区中心线的距离；
r1——安全棒区外径；
r2——活性区外径；
r3——径向反射层区外径；
λs——安全棒区导热系数；
λf——活性区导热；
λr——径向反射层区导热系数；
Ts(r)——距离安全棒中心线距离为r处的安全棒温度；
Tf(r)——距离安全棒中心线距离为r处的活性区温度；
Tr(r)——距离安全棒中心线距离为r处的径向反射层温度；
Ts(r1)——安全棒区外边界的温度；
Tf(r1)——活性区内边界的温度；
Tf(r2)——活性区外边界的温度；
Tr(r2)——径向反射层区内边界的温度；
Tr(r3)——径向反射层区外边界的温度；
ε——热管堆外容器壁的发射率；
σ——斯忒潘-珀耳兹曼常量；
可得到三个区域的温度分布的解析解，如公式(10)到公式(12)所示：
Ts(r)＝c1lnr+c2公式(10)
Tr(r)＝c5lnr+c6 公式(12)
其中
c1＝0 公式(13)
式中：
r——距离安全棒区中心线的距离；
r1——安全棒区外径；
r2——活性区外径；
r3——径向反射层区外径；
λs——安全棒区导热系数；
λf——活性区导热；
λr——径向反射层区导热系数；
Ts(r)——距离安全棒中心线距离为r处的安全棒温度；
Tf(r)——距离安全棒中心线距离为r处的活性区温度；
Tr(r)——距离安全棒中心线距离为r处的径向反射层温度；
Ts(r1)——安全棒区外边界的温度；
Tf(r1)——活性区内边界的温度；
Tf(r2)——活性区外边界的温度；
Tr(r2)——径向反射层区内边界的温度；
Tr(r3)——径向反射层区外边界的温度；
ε——热管堆外容器壁的发射率；
σ——斯忒潘-珀耳兹曼常量；
热管堆在瞬态过程中，随着功率的改变(即体积释热率的改变)，三个区的温度分布会随之进行变化；定义三个区域在瞬态前后的温度变化如公式(19)到公式(21)所示：
式中：
r——距离安全棒区中心线的距离；
r1——安全棒区外径；
r2——活性区外径；
r3——径向反射层区外径；
Ts(r)——瞬态前距离安全棒中心线距离为r处的安全棒温度；
Tf(r)——瞬态前距离安全棒中心线距离为r处的活性区温度；
Tr(r)——瞬态前距离安全棒中心线距离为r处的径向反射层温度；
Ts'(r)——瞬态后距离安全棒中心线距离为r处的安全棒温度；
Tf'(r)——瞬态后距离安全棒中心线距离为r处的活性区温度；
Tr'(r)——瞬态后距离安全棒中心线距离为r处的径向反射层温度；
ΔTs——瞬态前后安全棒区的温度变化量；
ΔTf——瞬态前后活性区的温度变化量；
ΔTr——瞬态前后径向反射层区的温度变化量；
步骤2：调研得到三个区域的材料在不同温度下的线性膨胀系数。线性热膨胀系数γlte是指固体物质的温度每改变1℃时，其长度的变化和它在0℃时长度之比；相应的，体积热膨胀系数γvte是指固体物质的温度改变1℃时，其体积的变化和它在0℃时体积之比。通常，γvte是γlte的3倍；基于体膨胀系数是线性膨胀系数的3倍的假设，分别列出三个区域的直径、高度以及密度随温度的变化函数，如公式(22)到公式(30)所示：
式中：
r1——安全棒区外径；
r2——活性区外径；
r3——径向反射层区外径；
r1(T)——温度为T时，安全棒区外径；
r2(T)——温度为T时，活性区外径；
r3(T)——温度为T时，径向反射层区外径；
——安全棒区的线性膨胀系数；
——活性区的线性膨胀系数；
——径向反射层区的线性膨胀系数；
ΔTs——瞬态前后安全棒区的温度变化量；
ΔTf——瞬态前后活性区的温度变化量；
ΔTr——瞬态前后径向反射层区的温度变化量；
h——常温下热管堆的高度；
hs(T)——温度为T时安全棒区高度；
hf(T)——温度为T时活性区高度；
hr(T)——温度为T时径向反射层区高度；
ρs——常温下安全棒区材料的密度；
ρf——常温下活性区材料的密度；
ρr——常温下径向反射层区材料的密度；
ρs(T)——温度为T时安全棒区材料的密度；
ρf(T)——温度为T时活性区材料的密度；
ρr(T)——温度为T时径向反射层区材料的密度；
步骤3：对热管堆进行输运计算之前，需要对其进行三角形网格划分；几何上，热管堆可由若干个扇环拼接而成；因此，对扇环进行三角形网格划分研究，并基于组合几何的方式由若干个扇环拼成1/6热管堆，从而对1/6热管堆进行三角网格划分；
圆、圆环、扇形均可视为特殊的扇环。其中，扇形为内径为0的扇环；圆环为辐角等于2π的扇环；圆为内径为0，辐角等于2π的扇环。因此，仅需对任意扇环的三角形网格划分进行研究，便可运用于1/6热管堆的三角形网格划分；
三角形网格的划分过程先后分为离散点的生成和三角形网格的生成；离散点的生成即三角形网格节点的生成过程；三角网格的生成即在离散点生成的基础上，将其连线，生成三角形网格的过程；
用户需要输入第i个扇环的圆心坐标(xi,yi)、初始辐角αi,1与终止辐角αi,2、扇环的内径ri,1与外径ri,2以及三角形网格的名义边长δi；如公式(31)所示，第i个扇环的径向可均匀插入ni个节点；式(31)中表示向下取整，由于取整的存在，三角形网格的名义边长δi并不一定等于三角形网格的实际边长，但其大小能反应网格的疏密程度；
式中：
δi——第i个扇环中三角形网格的名义边长；
ni——第i个扇环径向均匀插入的节点数；
ri,1——第i个扇环的内径；
ri,2——第i个扇环的外径；
相应地，扇形含有n+2条同圆心角的弧，每条弧长可由式(32)表示：
式中：
ri,1——第i个扇环的内径；
ri,2——第i个扇环的外径；
αi,1——第i个扇环的初始辐角；
αi,2——第i个扇环的终止辐角；
ni——第i个扇环径向均匀插入的节点数；
li，j——第i个扇环中第j条弧的长度；
同理可根据式(33)计算得到每一条弧上均匀插入的节点数mi：
式中：
mi，j——第i个扇环中第j条弧上均匀插入的节点数；
δi——第i个扇环中三角形网格的名义边长；
li，j——第i个扇环中第j条弧的长度；
因此，可以得到各个节点的坐标如公式(34)与公式(35)所示：
式中：
xi,j,k——第i个扇环中第j条弧上第k个节点的横坐标；
yi,j,k——第i个扇环中第j条弧上第k个节点的纵坐标；
ri,1——第i个扇环的内径；
ri,2——第i个扇环的外径；
ni——第i个扇环径向均匀插入的节点数；
mi，j——第i个扇环中第j条弧上均匀插入的节点数；
αi,1——第i个扇环的初始辐角；
αi,2——第i个扇环的终止辐角；
xi——第i个扇环圆心的横坐标；
yi——第i个扇环圆心的纵坐标；
基于以上离散点生成方法，本发明采用三角网生长法实现三角形网格的生成，在此不予以赘述；最终实现三角形网格的划分；
基于对扇环三角形网格划分的研究，本发明采用组合几何的思想，将1/6热管堆视为由18个扇环拼接后的模型。基于步骤三得到的不同区域在温度T下的内径和外径，可计算得到18个扇环中每个扇环的圆心坐标、内径、外径、初始辐角以及终止辐角；可指定每个扇环的三角形网格名义边长并实现1/6堆芯的三角形网格剖分；
步骤4：使用快堆时空中子动力学程序SARAX对1/6堆芯进行输运计算。将输运计算得到的堆芯径向功率分布提供给步骤1重新进行三个区域的温度分布计算；再进行步骤2的膨胀计算；再经过步骤3，将划分好的三角形网格提供给下一时间步的输运计算；依次循环，可计算热管堆的堆芯功率水平和反应性在瞬态过程中的变化曲线。
与现有的快堆时空中子动力学程序SARAX相比，本发明有如下突出优点：
1.本发明考虑了热管堆瞬态过程的热膨胀效应，并采用每个时间步更新三角形网格的方法进行瞬态计算。可更加真实、准确地模拟热管堆的瞬态过程。
2.针对中子泄漏严重的热管堆，相比于SARAX程序中采用虚拟密度理论考虑膨胀效应的方法，采用三角形网格剖分的动态方法可更加精确地描述中子输运方程中泄漏项，提高输运计算的精度。
附图说明
图1是热管堆堆芯精细几何模型图径向截面图。
图2是热管堆堆芯精细几何模型图轴向截面图。
图3是热管堆径向区域划分图。
图4是不同三角形名义边长下圆的网格剖分图，其中：图4(a)名义尺寸为1，图4(b)名义尺寸为1，图4(c)名义尺寸为1。
图5是1/6热管堆堆芯网格剖分图。
图6是热管堆瞬态分析流程图。
具体实施方式
本发明基于现有的快堆时空中子动力学程序SARAX，考虑了热管堆在瞬态过程中的热膨胀效应。通过进行安全棒区、活性区和径向反射层区的温度分布计算、热膨胀引起的结构尺寸计算以及三角形网格划分处理，将SARAX程序扩充了进行热管堆瞬态计算的功能。该计算方法包括以下方面：
1)对于以六角形方式排列的具有234根燃料热管组件的热管堆，沿径向将其由内而外依次划分为安全棒区、活性区以及径向反射层区。分别在各区域列出一维导热微分方程和边界条件。解得在不同功率水平下热管堆各区域的温度分布。
2)调研得到三个区域的材料在不同温度下的线性膨胀系数。基于体膨胀系数是线性膨胀系数的3倍的假设，分别列出三个区域的直径、高度以及密度随温度的变化函数；
3)对扇环进行三角形网格划分研究。并基于组合几何的方式由若干个扇环拼成1/6热管堆，从而对1/6热管堆进行三角网格划分；
4)使用SARAX程序对1/6堆芯进行输运计算。将输运计算得到的堆芯径向功率分布提供给步骤1重新进行温度分布计算；再进行步骤2的膨胀计算；再经过步骤3，将划分好的三角形网格提供给下一时间步的输运计算。依次循环。可计算热管堆的堆芯功率水平和反应性在瞬态过程中的变化曲线。
步骤1：如图1和图2所示，对于以六角形方式排列的具有234根燃料热管组件以及6根控制鼓的热管堆，沿径向将其由内而外依次划分为安全棒区、活性区以及径向反射层区，如图3所示。安全棒区为一圆柱体，直径为安全棒绝热层的外径；活性区为一圆环柱体，内径为安全棒绝热层外径，外径为内容器壁外径；径向反射层区仍为一圆环柱体，内径为内容器壁外径，外径为外容器壁外径。各区域的高均与热管堆的高相同。其中，安全棒区沿径向满足一维导热微分方程，如公式(1)所示；活性区沿径向满足一维含内热源导热微分方程，如公式(2)所示；径向反射层区沿径向满足一维导热微分方程，如公式(3)所示。
式中：
r——距离安全棒区中心线的距离；
λs——安全棒区导热系数；
λf——活性区导热系数；
λr——径向反射层区导热系数；
Ts(r)——安全棒区径向温度；
Tf(r)——燃料活性区径向温度；
Tr(r)——径向反射层区径向温度；
——活性区的体积释热率；
安全棒区的中心区满足绝热边界条件；安全棒区与活性区边界处温度连续并且热流密度连续；活性区与径向反射层区边界处温度连续并且热流密度连续；径向反射层的外边界满足辐射边界条件。因此可得如公式(4)到(9)所示的边界条件。
Ts(r1)＝Tf(r1) 公式(5)
Tf(r2)＝Tr(r2)公式(7)
式中：
r——距离安全棒区中心线的距离；
r1——安全棒区外径；
r2——活性区外径；
r3——径向反射层区外径；
λs——安全棒区导热系数；
λf——活性区导热；
λr——径向反射层区导热系数；
Ts(r)——距离安全棒中心线距离为r处的安全棒温度；
Tf(r)——距离安全棒中心线距离为r处的活性区温度；
Tr(r)——距离安全棒中心线距离为r处的径向反射层温度；
Ts(r1)——安全棒区外边界的温度；
Tf(r1)——活性区内边界的温度；
Tf(r2)——活性区外边界的温度；
Tr(r2)——径向反射层区内边界的温度；
Tr(r3)——径向反射层区外边界的温度；
ε——热管堆外容器壁的发射率；
σ——斯忒潘-珀耳兹曼常量；
可得到三个区域的温度分布的解析解，如公式(10)到公式(12)所示：
Ts(r)＝c1lnr+c2 公式(10)
Tr(r)＝c5lnr+c6 公式(12)
其中
c1＝0 公式(13)
式中：
r——距离安全棒区中心线的距离；
r1——安全棒区外径；
r2——活性区外径；
r3——径向反射层区外径；
λs——安全棒区导热系数；
λf——活性区导热；
λr——径向反射层区导热系数；
Ts(r)——距离安全棒中心线距离为r处的安全棒温度；
Tf(r)——距离安全棒中心线距离为r处的活性区温度；
Tr(r)——距离安全棒中心线距离为r处的径向反射层温度；
Ts(r1)——安全棒区外边界的温度；
Tf(r1)——活性区内边界的温度；
Tf(r2)——活性区外边界的温度；
Tr(r2)——径向反射层区内边界的温度；
Tr(r3)——径向反射层区外边界的温度；
ε——热管堆外容器壁的发射率；
σ——斯忒潘-珀耳兹曼常量；
热管堆在瞬态过程中，随着功率的改变(即体积释热率的改变)，三个区的温度分布会随之进行变化。定义三个区域在瞬态前后的温度变化如公式(19)到公式(21)所示：
式中：
r——距离安全棒区中心线的距离；
r1——安全棒区外径；
r2——活性区外径；
r3——径向反射层区外径；
Ts(r)——瞬态前距离安全棒中心线距离为r处的安全棒温度；
Tf(r)——瞬态前距离安全棒中心线距离为r处的活性区温度；
Tr(r)——瞬态前距离安全棒中心线距离为r处的径向反射层温度；
Ts'(r)——瞬态后距离安全棒中心线距离为r处的安全棒温度；
Tf'(r)——瞬态后距离安全棒中心线距离为r处的活性区温度；
Tr'(r)——瞬态后距离安全棒中心线距离为r处的径向反射层温度；
ΔTs——瞬态前后安全棒区的温度变化量；
ΔTf——瞬态前后活性区的温度变化量；
ΔTr——瞬态前后径向反射层区的温度变化量；
步骤2：调研得到三个区域的材料在不同温度下的线性膨胀系数。线性热膨胀系数γlte是指固体物质的温度每改变1℃时，其长度的变化和它在0℃时长度之比；相应的，体积热膨胀系数γvte是指固体物质的温度改变1℃时，其体积的变化和它在0℃时体积之比。通常，γvte是γlte的3倍。基于体膨胀系数是线性膨胀系数的3倍的假设，分别列出三个区域的直径、高度以及密度随温度的变化函数，如公式(22)到公式(30)所示：
式中：
r1——安全棒区外径；
r2——活性区外径；
r3——径向反射层区外径；
r1(T)——温度为T时，安全棒区外径；
r2(T)——温度为T时，活性区外径；
r3(T)——温度为T时，径向反射层区外径；
——安全棒区的线性膨胀系数；
——活性区的线性膨胀系数；
——径向反射层区的线性膨胀系数；
ΔTs——瞬态前后安全棒区的温度变化量；
ΔTf——瞬态前后活性区的温度变化量；
ΔTr——瞬态前后径向反射层区的温度变化量；
h——常温下热管堆的高度；
hs(T)——温度为T时安全棒区高度；
hf(T)——温度为T时活性区高度；
hr(T)——温度为T时径向反射层区高度；
ρs——常温下安全棒区材料的密度；
ρf——常温下活性区材料的密度；
ρr——常温下径向反射层区材料的密度；
ρs(T)——温度为T时安全棒区材料的密度；
ρf(T)——温度为T时活性区材料的密度；
ρr(T)——温度为T时径向反射层区材料的密度；
步骤3：对热管堆进行输运计算之前，需要对其进行三角形网格划分。几何上，热管堆可由若干个扇环拼接而成。因此，对扇环进行三角形网格划分研究。并基于组合几何的方式由若干个扇环拼成1/6热管堆，从而对1/6热管堆进行三角网格划分。
圆、圆环、扇形均可视为特殊的扇环。其中，扇形为内径为0的扇环；圆环为辐角等于2π的扇环；圆为内径为0，辐角等于2π的扇环。因此，仅需对任意扇环的三角形网格划分方法进行研究，便可运用于1/6热管堆的三角形网格划分。
三角形网格的划分过程先后分为离散点的生成和三角形网格的生成。离散点的生成即三角形网格节点的生成过程；三角网格的生成即在离散点生成的基础上，将其连线，生成三角形网格的过程。
用户需要输入第i个扇环的圆心坐标(xi,yi)；初始辐角αi,1与终止辐角αi,2；扇环的内径ri,1与外径ri,2；三角形网格的名义边长δi。如公式(31)所示，第i个扇环的径向可均匀插入ni个节点。式(31)中表示向下取整。由于取整的存在，三角形网格的名义边长δi并不一定等于三角形网格的实际边长，但其大小能反应网格的疏密程度。
式中：
δi——第i个扇环中三角形网格的名义边长；
ni——第i个扇环径向均匀插入的节点数；
ri,1——第i个扇环的内径；
ri,2——第i个扇环的外径；
相应地，扇形含有n+2条同圆心角的弧。每条弧长可由式(32)表示：
式中：
ri,1——第i个扇环的内径；
ri,2——第i个扇环的外径；
αi,1——第i个扇环的初始辐角；
αi,2——第i个扇环的终止辐角；
ni——第i个扇环径向均匀插入的节点数；
li，j——第i个扇环中第j条弧的长度；
同理可根据式(33)计算得到每一条弧上均匀插入的节点数mi：
式中：
mi，j——第i个扇环中第j条弧上均匀插入的节点数；
δi——第i个扇环中三角形网格的名义边长；
li，j——第i个扇环中第j条弧的长度；
因此，可以得到各个节点的坐标如公式(34)与公式(35)所示：
式中：
xi,j,k——第i个扇环中第j条弧上第k个节点的横坐标；
yi,j,k——第i个扇环中第j条弧上第k个节点的纵坐标；
ri,1——第i个扇环的内径；
ri,2——第i个扇环的外径；
ni——第i个扇环径向均匀插入的节点数；
mi，j——第i个扇环中第j条弧上均匀插入的节点数；
αi,1——第i个扇环的初始辐角；
αi,2——第i个扇环的终止辐角；
xi——第i个扇环圆心的横坐标；
yi——第i个扇环圆心的纵坐标；
基于以上离散点生成方法，本发明采用三角网生长法实现三角形网格的生成，在此不予以赘述。最终实现三角形网格的划分。如图4所示，为三种情况下相同的圆的三角形网格划分图，依次为三角形名义边长为1，2，3的情况。
基于对扇环三角形网格划分方法的研究，本发明采用组合几何的思想，将1/6热管堆视为由18个扇环拼接后的模型。基于步骤3得到的不同区域在温度T下的内径和外径，可计算得到18个扇环中每个扇环的圆心坐标、内径、外径、初始辐角以及终止辐角。可指定每个扇环的三角形网格名义边长并实现1/6堆芯的三角形网格剖分，如图5所示。
步骤4：如图6所示，使用快堆时空中子动力学程序SARAX对1/6堆芯进行输运计算。将输运计算得到的堆芯径向功率分布提供给步骤一重新进行三个区域的温度分布计算；再进行步骤二的膨胀计算；再经过步骤三，将划分好的三角形网格提供给下一时间步的输运计算。依次循环。可计算热管堆的堆芯功率水平和反应性在瞬态过程中的变化曲线。

Claims (1)

1.一种用于热管堆瞬态特性分析的方法，包括如下步骤：
步骤1：对于以六角形方式排列的具有234根燃料热管组件以及6根控制鼓的热管堆，沿径向将其由内而外依次划分为安全棒区、活性区以及径向反射层区；安全棒区为一圆柱体，直径为安全棒绝热层的外径；活性区为一圆环柱体，内径为安全棒绝热层外径，外径为内容器壁外径；径向反射层区仍为一圆环柱体，内径为内容器壁外径，外径为外容器壁外径；各区域的高均与热管堆的高相同；其中，安全棒区沿径向满足一维导热微分方程，如公式(1)所示；活性区沿径向满足一维含内热源导热微分方程，如公式(2)所示；径向反射层区沿径向满足一维导热微分方程，如公式(3)所示。
式中：
r——距离安全棒区中心线的距离；
λs——安全棒区导热系数；
λf——活性区导热系数；
λr——径向反射层区导热系数；
Ts(r)——安全棒区径向温度；
Tf(r)——燃料活性区径向温度；
Tr(r)——径向反射层区径向温度；
——活性区的体积释热率；
安全棒区的中心区满足绝热边界条件；安全棒区与活性区边界处温度连续并且热流密度连续；活性区与径向反射层区边界处温度连续并且热流密度连续；径向反射层的外边界满足辐射边界条件；因此得如公式(4)到(9)所示的边界条件；
Ts(r1)＝Tf(r1) 公式(5)
Tf(r2)＝Tr(r2) 公式(7)
式中：
r——距离安全棒区中心线的距离；
r1——安全棒区外径；
r2——活性区外径；
r3——径向反射层区外径；
λs——安全棒区导热系数；
λf——活性区导热；
λr——径向反射层区导热系数；
Ts(r)——距离安全棒中心线距离为r处的安全棒温度；
Tf(r)——距离安全棒中心线距离为r处的活性区温度；
Tr(r)——距离安全棒中心线距离为r处的径向反射层温度；
Ts(r1)——安全棒区外边界的温度；
Tf(r1)——活性区内边界的温度；
Tf(r2)——活性区外边界的温度；
Tr(r2)——径向反射层区内边界的温度；
Tr(r3)——径向反射层区外边界的温度；
ε——热管堆外容器壁的发射率；
σ——斯忒潘-珀耳兹曼常量；
得到三个区域的温度分布的解析解，如公式(10)到公式(12)所示：
Ts(r)＝c1lnr+c2 公式(10)
Tr(r)＝c5lnr+c6 公式(12)
其中
c1＝0 公式(13)
式中：
r——距离安全棒区中心线的距离；
r1——安全棒区外径；
r2——活性区外径；
r3——径向反射层区外径；
λs——安全棒区导热系数；
λf——活性区导热；
λr——径向反射层区导热系数；
Ts(r)——距离安全棒中心线距离为r处的安全棒温度；
Tf(r)——距离安全棒中心线距离为r处的活性区温度；
Tr(r)——距离安全棒中心线距离为r处的径向反射层温度；
Ts(r1)——安全棒区外边界的温度；
Tf(r1)——活性区内边界的温度；
Tf(r2)——活性区外边界的温度；
Tr(r2)——径向反射层区内边界的温度；
Tr(r3)——径向反射层区外边界的温度；
ε——热管堆外容器壁的发射率；
σ——斯忒潘-珀耳兹曼常量；
热管堆在瞬态过程中，随着功率的改变，即体积释热率的改变，三个区的温度分布会随之进行变化；定义三个区域在瞬态前后的温度变化如公式(19)到公式(21)所示：
式中：
r——距离安全棒区中心线的距离；
r1——安全棒区外径；
r2——活性区外径；
r3——径向反射层区外径；
Ts(r)——瞬态前距离安全棒中心线距离为r处的安全棒温度；
Tf(r)——瞬态前距离安全棒中心线距离为r处的活性区温度；
Tr(r)——瞬态前距离安全棒中心线距离为r处的径向反射层温度；
Ts'(r)——瞬态后距离安全棒中心线距离为r处的安全棒温度；
Tf'(r)——瞬态后距离安全棒中心线距离为r处的活性区温度；
Tr'(r)——瞬态后距离安全棒中心线距离为r处的径向反射层温度；
ΔTs——瞬态前后安全棒区的温度变化量；
ΔTf——瞬态前后活性区的温度变化量；
ΔTr——瞬态前后径向反射层区的温度变化量；
步骤2：调研得到三个区域的材料在不同温度下的线性膨胀系数，线性热膨胀系数γlte是指固体物质的温度每改变1℃时，其长度的变化和它在0℃时长度之比；相应的，体积热膨胀系数γvte是指固体物质的温度改变1℃时，其体积的变化和它在0℃时体积之比；基于体膨胀系数是线性膨胀系数的3倍的假设，分别列出三个区域的直径、高度以及密度随温度的变化函数，如公式(22)到公式(30)所示：
式中：
r1——安全棒区外径；
r2——活性区外径；
r3——径向反射层区外径；
r1(T)——温度为T时，安全棒区外径；
r2(T)——温度为T时，活性区外径；
r3(T)——温度为T时，径向反射层区外径；
——安全棒区的线性膨胀系数；
——活性区的线性膨胀系数；
——径向反射层区的线性膨胀系数；
ΔTs——瞬态前后安全棒区的温度变化量；
ΔTf——瞬态前后活性区的温度变化量；
ΔTr——瞬态前后径向反射层区的温度变化量；
h——常温下热管堆的高度；
hs(T)——温度为T时安全棒区高度；
hf(T)——温度为T时活性区高度；
hr(T)——温度为T时径向反射层区高度；
ρs——常温下安全棒区材料的密度；
ρf——常温下活性区材料的密度；
ρr——常温下径向反射层区材料的密度；
ρs(T)——温度为T时安全棒区材料的密度；
ρf(T)——温度为T时活性区材料的密度；
ρr(T)——温度为T时径向反射层区材料的密度；
步骤3：对热管堆进行输运计算之前，需要对其进行三角形网格划分；几何上，热管堆由若干个扇环拼接而成；因此，对扇环进行三角形网格划分研究，并基于组合几何的方式由若干个扇环拼成1/6热管堆，从而对1/6热管堆进行三角网格划分；
圆、圆环、扇形均视为特殊的扇环；其中，扇形为内径为0的扇环；圆环为辐角等于2π的扇环；圆为内径为0，辐角等于2π的扇环；因此，仅需对任意扇环的三角形网格划分进行研究，便能运用于1/6热管堆的三角形网格划分；
三角形网格的划分过程先后分为离散点的生成和三角形网格的生成；离散点的生成即三角形网格节点的生成过程；三角网格的生成即在离散点生成的基础上，将其连线，生成三角形网格的过程；
用户需要输入第i个扇环的圆心坐标(xi,yi)、初始辐角αi,1与终止辐角αi,2、扇环的内径ri,1与外径ri,2以及三角形网格的名义边长δi；如公式(31)所示，第i个扇环的径向均匀插入ni个节点；式(31)中表示向下取整，由于取整的存在，三角形网格的名义边长δi并不一定等于三角形网格的实际边长，但其大小能反应网格的疏密程度；
式中：
δi——第i个扇环中三角形网格的名义边长；
ni——第i个扇环径向均匀插入的节点数；
ri,1——第i个扇环的内径；
ri,2——第i个扇环的外径；
相应地，扇形含有n+2条同圆心角的弧，每条弧长由式(32)表示：
式中：
ri,1——第i个扇环的内径；
ri,2——第i个扇环的外径；
αi,1——第i个扇环的初始辐角；
αi,2——第i个扇环的终止辐角；
ni——第i个扇环径向均匀插入的节点数；
li，j——第i个扇环中第j条弧的长度；
同理根据式(33)计算得到每一条弧上均匀插入的节点数mi：
式中：
mi，j——第i个扇环中第j条弧上均匀插入的节点数；
δi——第i个扇环中三角形网格的名义边长；
li，j——第i个扇环中第j条弧的长度；
因此，能够得到各个节点的坐标如公式(34)与公式(35)所示：
式中：
xi,j,k——第i个扇环中第j条弧上第k个节点的横坐标；
yi,j,k——第i个扇环中第j条弧上第k个节点的纵坐标；
ri,1——第i个扇环的内径；
ri,2——第i个扇环的外径；
ni——第i个扇环径向均匀插入的节点数；
mi，j——第i个扇环中第j条弧上均匀插入的节点数；
αi,1——第i个扇环的初始辐角；
αi,2——第i个扇环的终止辐角；
xi——第i个扇环圆心的横坐标；
yi——第i个扇环圆心的纵坐标；
基于以上离散点生成方法，采用三角网生长法实现三角形网格的生成，在此不予以赘述，最终实现三角形网格的划分；
基于对扇环三角形网格划分的研究，采用组合几何的思想，将1/6热管堆视为由18个扇环拼接后的模型；基于步骤3得到的不同区域在温度T下的内径和外径，计算得到18个扇环中每个扇环的圆心坐标、内径、外径、初始辐角以及终止辐角；指定每个扇环的三角形网格名义边长并实现1/6堆芯的三角形网格剖分；
步骤4：使用快堆时空中子动力学程序SARAX对1/6堆芯进行输运计算；将输运计算得到的堆芯径向功率分布提供给步骤1重新进行三个区域的温度分布计算；再进行步骤2的膨胀计算；再经过步骤3，将划分好的三角形网格提供给下一时间步的输运计算；依次循环，则能够计算热管堆的堆芯功率水平和反应性在瞬态过程中的变化曲线。
CN201810523851.2A 2018-05-28 2018-05-28 一种用于热管堆瞬态特性分析的方法 Active CN108763748B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810523851.2A CN108763748B (zh) 2018-05-28 2018-05-28 一种用于热管堆瞬态特性分析的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810523851.2A CN108763748B (zh) 2018-05-28 2018-05-28 一种用于热管堆瞬态特性分析的方法

Publications (2)

Publication Number Publication Date
CN108763748A true CN108763748A (zh) 2018-11-06
CN108763748B CN108763748B (zh) 2019-04-09

Family

ID=64003090

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810523851.2A Active CN108763748B (zh) 2018-05-28 2018-05-28 一种用于热管堆瞬态特性分析的方法

Country Status (1)

Country Link
CN (1) CN108763748B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113362971A (zh) * 2021-05-21 2021-09-07 西安交通大学 一种用于静态转换的紧凑型热管堆堆芯结构

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101546612A (zh) * 2009-04-21 2009-09-30 中国原子能科学研究院 钠冷快堆核电站冷却剂系统和部件的设计瞬态确定方法
JP6071404B2 (ja) * 2012-10-12 2017-02-01 株式会社東芝 原子力プラントおよび静的格納容器冷却系
CN106528985A (zh) * 2016-11-03 2017-03-22 哈尔滨工程大学 一种核动力装置冷凝器的分区化仿真方法
US20170162281A1 (en) * 2014-07-03 2017-06-08 Kabushiki Kaisha Toshiba Passive containment cooling and filtered venting system, and nuclear power plant
CN107478715A (zh) * 2017-07-03 2017-12-15 岭东核电有限公司 核电站热交换器传热管的无损检测分析方法、装置及系统

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101546612A (zh) * 2009-04-21 2009-09-30 中国原子能科学研究院 钠冷快堆核电站冷却剂系统和部件的设计瞬态确定方法
JP6071404B2 (ja) * 2012-10-12 2017-02-01 株式会社東芝 原子力プラントおよび静的格納容器冷却系
US20170162281A1 (en) * 2014-07-03 2017-06-08 Kabushiki Kaisha Toshiba Passive containment cooling and filtered venting system, and nuclear power plant
CN106528985A (zh) * 2016-11-03 2017-03-22 哈尔滨工程大学 一种核动力装置冷凝器的分区化仿真方法
CN107478715A (zh) * 2017-07-03 2017-12-15 岭东核电有限公司 核电站热交换器传热管的无损检测分析方法、装置及系统

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
WENWEN ZHANG.ETC.: ""Transient thermal-hydraulic analysis of a space thermionic reactor"", 《ANNALS OF NUCLEAR ENERGY》 *
屈伸等: ""热管堆高温数据库的制作及堆芯初步物理计算"", 《现代应用物理》 *
张文文等: ""热管改进型热离子反应堆瞬态分析程序开发"", 《原子能科学技术》 *
郑友琦等: ""快堆中子学计算程序NECP-SARAX的改进与验证"", 《强激光与粒子束》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113362971A (zh) * 2021-05-21 2021-09-07 西安交通大学 一种用于静态转换的紧凑型热管堆堆芯结构
CN113362971B (zh) * 2021-05-21 2022-10-28 西安交通大学 一种用于静态转换的紧凑型热管堆堆芯结构

Also Published As

Publication number Publication date
CN108763748B (zh) 2019-04-09

Similar Documents

Publication Publication Date Title
CN103150424B (zh) 一种获取反应堆堆芯三维中子通量密度精细分布的方法
Bouhacina et al. Numerical investigation of a novel tube design for the geothermal borehole heat exchanger
CN107066745B (zh) 获取快中子堆堆芯瞬态过程三维中子通量密度分布的方法
Beck et al. Geometric arrangement and operation mode adjustment in low-enthalpy geothermal borehole fields for heating
Bandos et al. Finite cylinder-source model for energy pile heat exchangers: Effects of thermal storage and vertical temperature variations
Li et al. Influence of ground surface boundary conditions on horizontal ground source heat pump systems
Liu et al. A novel 2-D ring-tubes model and numerical investigation of heat and moisture transfer around helix ground heat exchanger
Bouhal et al. Towards an energy efficiency optimization of solar horizontal storage tanks and circulation pipes integrating evacuated tube collectors through CFD parametric studies
Cui et al. 3D transient heat transfer numerical analysis of multiple energy piles
Wang et al. A novel composite-medium solution for pile geothermal heat exchangers with spiral coils
CN108763748B (zh) 一种用于热管堆瞬态特性分析的方法
Yang et al. Effects of groundwater pumping on ground surface temperature: A regional modeling study in the North China Plain
Li A new constant heat flux model for vertical U-tube ground heat exchangers
Jia et al. Analytical heat transfer model for coaxial heat exchangers based on varied heat flux with borehole depth
Xu et al. Shape optimization of composite porous vapor chamber with radial grooves: A study on the minimization of maximum pressure drop
Cao et al. A modified particle filter‐based data assimilation method for a high‐precision 2‐D hydrodynamic model considering spatial‐temporal variability of roughness: Simulation of dam‐break flood inundation
Mingzhi et al. Simplified heat transfer analysis method for large-scale boreholes ground heat exchangers
Yang et al. Optimal cooling channel layout in a hot enclosure subject to natural convection
Liao et al. A novel 2-D equivalent numerical model of helix energy pile based on heat transfer characteristics of internal heat convection
Gong et al. Analysis of thermal interaction coefficient for multiple borehole heat exchangers in layered soil considering groundwater seepage
Yamaguchi et al. Fast crustal deformation computing method for multiple computations accelerated by a graphics processing unit cluster
CN110489912A (zh) 一种太阳能跨季节土壤蓄热分层切片数值模拟的方法
Zhang et al. Calculation of particle volume fraction in computational fluid dynamics-discrete element method simulation of particulate flows with coarse particles
Filimonov et al. Optimal simulation of design and operation of geothermal systems
Sun et al. Enhanced geothermal system productivity analysis of a well-group in a limited area based on the flow field split method

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant