CN110705184A - 一种反应堆堆芯精细化数值求解的虚拟体积力动量源法 - Google Patents

Abstract

本发明涉及三维精细化CFD数值计算方法领域，特别涉及一种反应堆堆芯精细化数值求解的虚拟体积力动量源法，所述方法包括如下步骤：将棒束通道内的交混翼片去除，在翼片原来所占空间内对流体施加虚拟体积力，作为动量方程的源项，以考虑交混翼片对流场的影响；将各个方向的虚拟体积力加入到动量守恒方程中，将各分力在直角坐标系下分解，求得直角坐标系下x，y，z三个方向的分力，并进行数学描述。本发明利用动量方程不断迭代求解流场，从而完全实现在无交混翼片的简单通道中模拟交混翼流场的效果，极大地提高了三维精细化CFD数值计算的效率和计算精度。

Description

一种反应堆堆芯精细化数值求解的虚拟体积力动量源法

技术领域

本发明涉及三维精细化CFD数值计算方法领域，特别涉及一种反应堆堆芯精细化数值求解的虚拟体积力动量源法领域。

背景技术

能源需求的快速增长与能源短缺之间的矛盾已成为当今世界主要问题之一。有研究表明，到2030年世界范围内的能源消耗将比2010年增加35％到49％。新型清洁能源的开发利用既能缓解日益严重的能源危机，又能在一定程度上抑制全球变暖问题。上个世纪中期，核能利用在这样的背景下应运而生。目前各种类型的核反应堆已广泛应用于电力、供热和动力工程等领域。核反应堆中，发展最为成熟可靠且应用最为广泛的是轻水堆(LWR)，包括压水堆(PWR)和沸水堆(BWR)，以水或汽水混合物作为慢化冷却剂。典型压水堆堆芯由方形排列的燃料棒组成，同时燃料组件竖直方向上设置4到7个起固定和交混作用的交混翼定位格架。典型压水堆燃料组件、燃料棒及交混翼格架如附图1a、附图1b、附图1c、附图2a及附图2b所示。

随着核电发展，对核反应堆堆芯的热工水力及结构优化设计提出了更高的要求。棒束通道作为压水堆堆芯的核心部件，其热工水力性能直接影响核电厂的安全性、可靠性和经济性。棒束通道本身是一种狭长型的非圆管流动通道，其特点是结构紧凑，有效传热面积大，换热效率高，非常适合于反应堆内热能的输送。但是，如何有效地控制所允许的沸腾程度以及在出现事故后有效排放余热不致引起“烧毁”，是反应堆热工水力安全设计的主要任务之一。另外，起固定和交混作用的定位格架对棒束通道热工水力性能有两个完全相反的影响：增大压降和加强换热。由此可见，交混翼定位格架的热工水力设计还存在最优问题。

随着高可靠性和高经济性新型反应堆的出现，对核反应堆热工水力设计和安全分析提出了更高要求，一维系统程序由于其无法提供堆芯内各复杂结构处局部流动和传热细节，越来越难以胜任新的设计要求。近二十年来，随着计算机和计算流体力学(CFD)技术的大力发展，三维计算流体力学求解程序由于其精细程度高、可重复性好、成本要求低、实施周期短及无安全忧患等诸多优点越来越多地应用于反应堆和核能系统热工水力分析中。然而，对于反应堆堆芯这种结构复杂且跨尺度问题严重的热工水力设备，要对其实施全尺度甚至单个燃料组件的数值计算，采用常规数值计算方法在现有计算机硬件条件下仍然难以实现，即使在超大型并行计算集群上勉强能够进行相关计算，仍然存在前处理、数值求解及后处理效率过低且成本过高和普及性差等问题。

反应堆堆芯燃料棒区域带交混翼格架的大规模棒束通道数值计算的难点主要在于由不规则翼片引起的巨大数量的非结构化网格。尽管“分段求解技术”和“多面体网格技术”相比于传统计算方法能大大提高计算效率，能有效扩大数值计算中的棒束规模和解决轴向求解难题，但考虑到真实反应堆压力容器内构件众多且结构复杂，同时燃料组件本身棒束规模巨大且跨尺度问题严重，上述方法仍然难以解决在现有计算机硬件水平下实现全堆芯全尺度范围内的数值求解难题。事实上，棒束通道内交混翼格架对流动及传热特性的影响主要有：所制造的二次流场对冷却剂横向交混能力的加强、对流换热能力的强化及压降的增大。进一步分析可知，对流换热能力及压降的增大本质上均由二次流场的搅混效应导致，而二次流场又主要由格架带条及翼片的阻碍和导流作用引起。因此，如何从棒束通道内流动和传热现象的物理本质及数值求解策略等方面出发，开发新的数值模拟方法或技术，摒弃传统的数值计算方法，实现更为高效的数值计算方法从而解决堆芯全尺度的数值求解是当前CFD技术能否全面和高效应用到反应堆设计和安全等分析中去的关键。

发明内容

本发明的目的是提出一种反应堆堆芯精细化数值求解的虚拟体积力动量源法，在兼顾三维精细化CFD数值计算的准确性和快速性的基础上，本发明提出的虚拟体积力动量源法是将棒束通道内的交混翼片去除，在翼片原来所占空间内对流体施加虚拟体积力，作为动量方程的源项，以考虑交混翼片对流场的影响。虚拟体积力动量源法在保证虚拟体积力只是当地速度场函数的同时在网格划分时不再需要顾及翼片所占空间的形状，从而整个棒束通道均可采用高质量全结构化网格，完全实现在无交混翼片的简单通道中模拟交混翼流场的效果。

本发明是这样实现的：

一种反应堆堆芯精细化数值求解的虚拟体积力动量源法，所述方法包括如下步骤：

(1)将棒束通道内的交混翼片去除，在翼片原来所占空间内对流体施加虚拟体积力，作为动量方程的源项，以考虑交混翼片对流场的影响；

(2)将各个方向的虚拟体积力加入到动量守恒方程中，将各分力在直角坐标系下分解，求得直角坐标系下x，y，z三个方向的分力，并进行数学描述；

(3)虚拟体积力与翼片形状、排列及偏折角有关，是翼片的形状参数A_sp、排列方式参数A_ag及包含偏折方向的偏折角A_θ等固有几何属性参数的函数f(A_sp,A_ag,A_θ)，确定A_sp,A_ag,A_θ及其与虚拟体积力的关系f(A_sp,A_ag,A_θ)，对各个方向虚拟体积力进行数学描述；

(4)采用封闭空间点包含识别算法来判断流场中的网格节点哪些位于固体区域内，哪些位于固体区域的外，哪些在固体区域的表面上，再把动量源项加载到位于固体内和固体表面节点的动量方程上；

(5)采用封闭空间点包含识别算法对通道内的刚凸或弹簧部件采用虚拟体积力动量源法进行数学描述，通过在上述各虚拟体积力的求解公式中通过对各系数进行修正来考虑刚凸或弹簧的影响；

(6)采用虚拟体积力动量源法在进行CFD数值迭代求解时，各个力的大小及方向由初始速度场得到，各个方向的虚拟体积力反过来作用于速度场，由于各虚拟体积力均为当地速度场函数，因此整个数值计算是一个自动调节过程，最终获得收敛解。

所述步骤(1)中虚拟体积力F_v可以分解为与交混翼垂直的方向n_n的分力F_n，与交混翼平行的方向n_p的分力F_p，与交混翼相切的方向n_t的分力F_t，则虚拟体积力的合力可表示为：F_v＝F_n+F_p+F_t。

所述翼片法向虚拟体积力F_n可表示为流体质点沿翼片法向方向加速度的形式，其计算式为：

其中，ρ表示流体密度，V_n表示翼片法向的流体速度，V_t表示翼片切向的流体速度，V_p表示翼片平行方向的流体速度。

所述法向速度V_n沿各个方向的梯度表示成如下形式：

其中，C为乘数因子，用于修正数值计算结果，由于实际情况下流体不能穿过翼片，因此C的取值必须大于2；L_n为翼片厚度；a_t为切向厚度因子，考虑到切向速度梯度必然小于法向梯度，因此该值应大于1；a_p为平行方向厚度因子，该方向速度梯度同样小于法向梯度，因此其取值也应大于1，a_t和a_p同样可用来修正计算结果。

所述平行分力F_p与径向分力F_t的计算方法为：

其中，f_t及f_p分别表示各自方向上的摩擦阻力因子。

所述步骤(2)中，直角坐标系下x,y,z坐标轴的三个分力F_x,F_y,F_z所述分力的数学描述为：

其中，F_n,x表示F_n在直角坐标系下x方向的分力，F_p,x表示F_p在直角坐标系下x方向的分力，F_n,z表示F_n在直角坐标系下z方向的分力，F_p,z表示F_p在直角坐标系下z方向的分力。

所述步骤(4)中的封闭空间点包含识别算法，其步骤为：

1)将交混翼片等固体区域当作一个封闭空间，输入封闭边界

2)输入待查询点q；

3)计算向量

4)计算p点处的“综合法向向量

5)判断：如果

则“q点位于封闭区域边界上”；如果

则“q点位于封闭区域边界外”；否则“q点位于封闭区域边界内”；

6)若查询点q位于封闭区域边界上或封闭区域边界内，则节点q的动量方程需加载虚拟体积力动量源项；若查询点q位于封闭区域边界外，则节点q的动量方程不加载虚拟体积力动量源项。

所述虚拟体积力动量源法进行CFD数值迭代求解，其求解步骤为：

步骤一、在进行数值迭代求解时，先给出一个初始速度场；

步骤二、翼片或刚凸或弹簧的各虚拟体积力均为当地速度场的函数，将三个方向的虚拟体积力表示成了翼片法向速度V_n、切向速度V_t及平行速度V_p的函数；

步骤三、判断节点的空间位置，若在固体区域边界上或固体区域内，则该节点的动量方程需添加虚拟体积力动量源项，否则不添加虚拟体积力动量源项，翼片或刚凸或弹簧的法向虚拟体积力F_n、平行分力F_P和径向分力F_t又反过来作用于速度场V，利用动量方程计算得到新的速度场，获得直角坐标系下的V_x，V_y，V_z；

步骤四、不断迭代计算步骤一～步骤二，直到获得收敛解。

本发明的有益效果是：

本发明采用的虚拟体积力动量源法，将虚拟体积力用当地速度场函数的形式表达出来，无需借助完整几何模型的LES计算结果。在计算节点网格划分时又不再需要顾及翼片所占空间的形状，所以整个棒束通道均可采用高质量全结构化网格，降低了网格数量，提高了网格质量。再利用“封闭空间点包含识别算法”判断哪些计算节点需要在动量方程里添加虚拟体积力动量源项，然后再利用动量方程不断迭代求解流场，从而完全实现在无交混翼片的简单通道中模拟交混翼流场的效果，极大地提高了三维精细化CFD数值计算的效率和计算精度。

附图说明

图1a为压水堆燃料组件中格架与交混翼结构示意图；

图1b为压水堆组件中燃料组件与棒束通道结构示意图；

图1c为压水堆组件中燃料棒结构示意图；

图2a为压水堆棒束通道中交混翼片布置示意图；

图2b为压水堆棒束通道中交混翼片结构示意图；

图3为虚拟体积力分解示意图；

图4为封闭空间点包含识别算法流程图；

图5为采用虚拟体积力动量源法进行复杂流场CFD数值迭代求解的流程图。

具体实施方式

结合附图对本发明作进一步详细描述。

本发明提出的虚拟体积力动量源法是将棒束通道内的交混翼片去除，在翼片原来所占空间内对流体施加虚拟体积力，作为动量方程的源项，以考虑交混翼片对流场的影响。虚拟体积力动量源法在保证虚拟体积力只是当地速度场函数的同时在网格划分时不再需要顾及翼片所占空间的形状，从而整个棒束通道均可采用高质量全结构化网格，完全实现在无交混翼片的简单通道中模拟交混翼流场的效果。

为达到上述目的，本发明的技术方案包含以下步骤：

步骤1.将棒束通道内的交混翼片去除，在翼片原来所占空间内对流体施加虚拟体积力，作为动量方程的源项，以考虑交混翼片对流场的影响；

步骤2.将各个方向的虚拟体积力加入到动量守恒方程中，将各分力在直角坐标系下分解，求得直角坐标系下x，y，z三个方向的分力，并进行数学描述，虚拟体积力分解示意图如附图3所示；

步骤3.虚拟体积力与翼片形状、排列及偏折角有关，是翼片的形状参数(A_sp)、排列方式参数(A_ag)及偏折角(A_θ，包含偏折方向)等固有几何属性参数的函数f(A_sp,A_ag,A_θ)，本步骤要确定A_sp,A_ag,A_θ及其与虚拟体积力的关系f(A_sp,A_ag,A_θ)，对各个方向虚拟体积力进行数学描述；

步骤4.本发明采用“封闭空间点包含识别算法”来判断流场中的网格节点哪些位于固体区域内，哪些位于固体区域的外，哪些在固体区域的表面上，再把动量源项加载到那些位于固体内和固体表面节点的动量方程上，本步骤具体流程如附图4所示；

步骤5.同样利用上述方法，不仅可以完全实现在无交混翼片的燃料组件与棒束通道中模拟交混翼流场的效果，还可以对通道内的刚凸或弹簧等部件采用虚拟体积力动量源法进行数学描述，也可在上述各虚拟体积力的求解公式中通过对各系数进行修正来考虑刚凸或弹簧的影响；

步骤6.采用虚拟体积力动量源法在进行CFD数值迭代求解时，各个力的大小及方向由初始速度场得到，各个方向的虚拟体积力反过来作用于速度场，由于各虚拟体积力均为当地速度场函数，因此整个数值计算是一个自动调节过程，最终获得收敛解，迭代计算过程如附图5所示。

本发明进一步描述如下：

本发明是一种涉及解决反应堆堆芯范围内复杂模型的三维精细化数值求解难题的虚拟体积力动量源法，将压水堆燃料组件棒束通道内交混翼片(刚凸或弹簧等部件)对流场的影响用虚拟体积力加入到动量守恒方程中，再考虑翼片形状、排列及偏折角等因素，实现在无交混翼片(刚凸或弹簧等部件)的简单通道中模拟交混翼流场的效果，在CFD数值迭代求解过程中获得收敛解，关键步骤如下：

步骤1.将棒束通道内的交混翼片去除，在翼片原来所占空间内对流体施加虚拟体积力，作为动量方程的源项，以考虑交混翼片对流场的影响：

如附图3所示，对于任一交混翼片上的虚拟体积力F_v，可以分解为三个方向：与交混翼垂直的方向n_n，与交混翼平行的方向n_p，与交混翼相切的方向n_t。若上述三个方向的分力分别为F_n，F_p及F_t，则虚拟体积力的合力可以表示为：

F_v＝F_n+F_p+F_t (1)

步骤2.将各个方向的虚拟体积力加入到动量守恒方程中，将各分力在直角坐标系下分解，求得直角坐标系下x，y，z三个方向的分力，并进行数学描述，包括以下几个子步骤：

步骤2.1将各分力在直角坐标系下分解，虚拟体积力分解示意图如附图3所示：

式中，F_n,x表示F_n在直角坐标系下x方向的分力，其它力依此类推。

步骤2.2将各方向虚拟体积力表示成当地速度场的函数。从附图3所示的三个虚拟体积力分力来看，平行分力F_p和径向分力F_t的数学表达相对简单，可直接采用通用摩擦阻力计算公式，法向分力F_n则需基于当地流场质点加速度。稳态工况下，翼片法向虚拟体积力F_n可表示为流体质点沿翼片法向方向加速度(质点导数)的形式：

式中，V_n，V_t，V_p分别表示翼片法向、切向及平行方向的流体速度。另外，本发明将公式(3)中法向速度V_n沿各个方向的梯度表示成如下形式：

式中，C为乘数因子，可用于修正数值计算结果，由于实际情况下流体不能穿过翼片，因此C的取值必须大于2；L_n为翼片厚度；a_t为切向厚度因子，考虑到切向速度梯度必然小于法向梯度，因此该值应大于1；a_p为平行方向厚度因子，该方向速度梯度同样小于法向梯度，因此其取值也应大于1。a_t和a_p同样可用来修正计算结果。

平行分力F_p和径向分力F_t采用通用摩擦阻力公式：

式中，f_t及f_p分别表示各自方向上的摩擦阻力因子；a_tL_n及a_pL_n与公式(4)相同，用以将摩擦阻力标准化为单位长度上的力；

步骤3.虚拟体积力与翼片形状、排列及偏折角有关，是翼片的形状参数(A_sp)、排列方式参数(A_ag)及偏折角(A_θ，包含偏折方向)等固有几何属性参数的函数f(A_sp,A_ag,A_θ)，本步骤要确定A_sp,A_ag,A_θ及其与虚拟体积力的关系f(A_sp,A_ag,A_θ)，对各个方向虚拟体积力进行数学描述；

步骤4.本发明采用“封闭空间点包含识别算法”来判断流场中的网格节点哪些位于固体区域内，哪些位于固体区域的外，哪些在固体区域的表面上，再把动量源项加载到那些位于固体内和固体表面节点的动量方程上，具体包括以下几个子步骤：

步骤4.1将交混翼片等固体区域当作一个封闭空间，输入封闭边界

步骤4.2输入待查询点q；

步骤4.3计算向量

步骤4.4计算p点处的“综合法向向量

步骤4.5判断：如果

则“q点位于封闭区域边界上”；如果

则“q点位于封闭区域边界外”；否则“q点位于封闭区域边界内”；

步骤4.6若查询点q位于封闭区域边界上或封闭区域边界内，则节点q的动量方程需加载虚拟体积力动量源项；若查询点q位于封闭区域边界外，则节点q的动量方程不加载虚拟体积力动量源项；

步骤5.同样利用上述方法，不仅可以完全实现在无交混翼片的燃料组件与棒束通道中模拟交混翼流场的效果，还可以对通道内的刚凸或弹簧等部件采用同样的虚拟体积力动量源方法进行数学描述，也可在上述各虚拟体积力的求解公式中通过对各系数进行修正来考虑刚凸或弹簧的影响；

步骤6.采用虚拟体积力动量源法在进行CFD数值迭代求解时，各个力的大小及方向由初始速度场得到，各个方向的虚拟体积力反过来作用于速度场，由于各虚拟体积力均为当地速度场函数，因此整个数值计算是一个自动调节过程，最终获得收敛解，迭代计算过程如附图5所示，具体分为以下几个步骤：

步骤6.1.在进行数值迭代求解时，先给出一个初始速度场；

步骤6.2翼片(刚凸或弹簧等部件)的各虚拟体积力均为当地速度场的函数，将三个方向的虚拟体积力表示成了翼片法向速度V_n、切向速度V_t及平行速度V_p的函数，显然法向虚拟体积力F_n、平行分力F_p和径向分力F_t的大小及方向可由直角坐标系下的V_x，V_y，V_z根据公式(3)和(5)求解得到；

步骤6.3.判断节点的空间位置，若在固体区域边界上或固体区域内，则该节点的动量方程需添加虚拟体积力动量源项，否则不添加虚拟体积力动量源项，翼片(刚凸或弹簧等部件)的法向虚拟体积力F_n、平行分力F_P和径向分力F_t又反过来作用于速度场V，利用动量方程计算得到新的速度场，获得直角坐标系下的V_x，V_y，V_z；

步骤6.4.不断迭代计算步骤6.1～步骤6.2，直到获得收敛解。

Claims (8)

1.一种反应堆堆芯精细化数值求解的虚拟体积力动量源法，其特征是：所述方法包括如下步骤：

(1)将棒束通道内的交混翼片去除，在翼片原来所占空间内对流体施加虚拟体积力，作为动量方程的源项，以考虑交混翼片对流场的影响；

(2)将各个方向的虚拟体积力加入到动量守恒方程中，将各分力在直角坐标系下分解，求得直角坐标系下x，y，z三个方向的分力，并进行数学描述；

(3)虚拟体积力与翼片形状、排列及偏折角有关，是翼片的形状参数A_sp、排列方式参数A_ag及包含偏折方向的偏折角A_θ等固有几何属性参数的函数f(A_sp,A_ag,A_θ)，确定A_sp,A_ag,A_θ及其与虚拟体积力的关系f(A_sp,A_ag,A_θ)，对各个方向虚拟体积力进行数学描述；

(4)采用封闭空间点包含识别算法来判断流场中的网格节点哪些位于固体区域内，哪些位于固体区域的外，哪些在固体区域的表面上，再把动量源项加载到位于固体内和固体表面节点的动量方程上；

(5)采用封闭空间点包含识别算法对通道内的刚凸或弹簧部件采用虚拟体积力动量源法进行数学描述，通过在上述各虚拟体积力的求解公式中通过对各系数进行修正来考虑刚凸或弹簧的影响；

(6)采用虚拟体积力动量源法在进行CFD数值迭代求解时，各个力的大小及方向由初始速度场得到，各个方向的虚拟体积力反过来作用于速度场，由于各虚拟体积力均为当地速度场函数，因此整个数值计算是一个自动调节过程，最终获得收敛解。

2.根据权利要求1所述的一种反应堆堆芯精细化数值求解的虚拟体积力动量源法，其特征是：所述步骤(1)中虚拟体积力F_v可以分解为与交混翼垂直的方向n_n的分力F_n，与交混翼平行的方向n_p的分力F_p，与交混翼相切的方向n_t的分力F_t，则虚拟体积力的合力可表示为：F_v＝F_n+F_p+F_t。

3.根据权利要求2所述的一种反应堆堆芯精细化数值求解的虚拟体积力动量源法，其特征是：所述翼片法向虚拟体积力F_n可表示为流体质点沿翼片法向方向加速度的形式，其计算式为：

其中，ρ表示流体密度，V_n表示翼片法向的流体速度，V_t表示翼片切向的流体速度，V_p表示翼片平行方向的流体速度。

4.根据权利要求3所述的一种反应堆堆芯精细化数值求解的虚拟体积力动量源法，其特征是：所述法向速度V_n沿各个方向的梯度表示成如下形式：

其中，C为乘数因子，用于修正数值计算结果，由于实际情况下流体不能穿过翼片，因此C的取值必须大于2；L_n为翼片厚度；a_t为切向厚度因子，考虑到切向速度梯度必然小于法向梯度，因此该值应大于1；a_p为平行方向厚度因子，该方向速度梯度同样小于法向梯度，因此其取值也应大于1，a_t和a_p同样可用来修正计算结果。

5.根据权利要求2所述的一种反应堆堆芯精细化数值求解的虚拟体积力动量源法，其特征是：所述平行分力F_p与径向分力F_t的计算方法为：

其中，f_t及f_p分别表示各自方向上的摩擦阻力因子。

6.根据权利要求1所述的一种反应堆堆芯精细化数值求解的虚拟体积力动量源法，其特征是：所述步骤(2)中，直角坐标系下x,y,z坐标轴的三个分力F_x,F_y,F_z所述分力的数学描述为：

其中，F_n,x表示F_n在直角坐标系下x方向的分力，F_p,x表示F_p在直角坐标系下x方向的分力，F_n,z表示F_n在直角坐标系下z方向的分力，F_p,z表示F_p在直角坐标系下z方向的分力。

7.根据权利要求1所述的一种反应堆堆芯精细化数值求解的虚拟体积力动量源法，其特征是：所述步骤(4)中的封闭空间点包含识别算法，其步骤为：

1)将交混翼片等固体区域当作一个封闭空间，输入封闭边界

2)输入待查询点q；

3)计算向量

4)计算p点处的“综合法向向量”；

5)判断：如果

则“q点位于封闭区域边界上”；如果

则“q点位于封闭区域边界外”；否则“q点位于封闭区域边界内”；

6)若查询点q位于封闭区域边界上或封闭区域边界内，则节点q的动量方程需加载虚拟体积力动量源项；若查询点q位于封闭区域边界外，则节点q的动量方程不加载虚拟体积力动量源项。

8.根据权利要求1所述的一种反应堆堆芯精细化数值求解的虚拟体积力动量源法，其特征是：所述虚拟体积力动量源法进行CFD数值迭代求解，其求解步骤为：

步骤一、在进行数值迭代求解时，先给出一个初始速度场；

步骤二、翼片或刚凸或弹簧的各虚拟体积力均为当地速度场的函数，将三个方向的虚拟体积力表示成了翼片法向速度V_n、切向速度V_t及平行速度V_p的函数；

步骤三、判断节点的空间位置，若在固体区域边界上或固体区域内，则该节点的动量方程需添加虚拟体积力动量源项，否则不添加虚拟体积力动量源项，翼片或刚凸或弹簧的法向虚拟体积力F_n、平行分力F_P和径向分力F_t又反过来作用于速度场V，利用动量方程计算得到新的速度场，获得直角坐标系下的V_x，V_y，V_z；

步骤四、不断迭代计算步骤一～步骤二，直到获得收敛解。

Images