1 引言

核反应堆安全壳是构成压水反应堆(pressurized water reactor, 简称PWR)最外围的建筑, 是指包容了核蒸汽供应系统的大部分系统和设备的外壳建筑, 用以容纳反应堆压力容器以及部分安全系统, 将其与外部环境完全隔离, 是保护核反应堆安全运转的最后屏障.PCCS(passive containment cooling system, 非能动安全壳冷却系统)^[1]是核电工业中保护大型先进压水堆的重要安全系统, 它不需要配置冗余的电源、设备冷却水、通风与空调等能动安全系统, 仅利用非能动系统与部件, 依靠自然物理规律, 如重力注射、压缩气体膨胀、冷却剂蒸发冷凝等来实现安全功能^[2].

科学计算可视化的目的是将抽象的科学数据转换为视觉信息, 使人们能够充分地洞察数据所蕴含的丰富内涵.用户通过对可视化的感知, 使用可视化交互工具进行数据分析, 获取知识, 并进一步提升为智慧^[3].将科学计算可视化技术应用于核反应堆安全壳冷却系统相关测量数据、仿真计算数据的分析, 可以帮助人们更好地了解安全壳所处的物理状态, 更早地发现可能潜在的各类安全问题, 确保安全壳冷却系统及整个核电站的正常运行.

安全壳冷却系统中包含多种需要分析的物理数据.从可视化的角度来看, 这些数据主要分为两类, 即标量数据:温度、压力、放射性浓度等; 矢量数据:气相速度等.对于三维标量场, 应用最广泛的可视化方法为等值面法, 其中, 等值面提取算法包括经典的Marching Cubes^[4]及之后的各种优化算法, 如Ciżnicki等人^[5]提出的并行化的Marching Tetrahedra算法、Scholz等人^[6]提出的视点相关的LOD算法.对于三维矢量场, 可视化方法主要包括几何图标法、矢量线法、矢量面法、矢量体法等分类, 相关研究包括如Spencer等人^[7]提出的基于图像的曲面均匀分布流线生成算法、Krishnan等人^[8]提出的时间面和脉面生成算法.

在对安全壳冷却系统进行数值模拟计算时, 网格划分是前处理过程中的关键, 所采用的划分算法直接影响了计算过程的效率与计算结果的精度.网格划分相关的研究工作包括Su等人^[9]提出的改进波前法和曲面展开的自动化划分方法及各类并行网格划分算法^{[10, 11]}.对于网格划分得到的单元数据, 通常使用改进的八叉树结构进行存储, 如Wilhelms等人^[12]提出的BONO及Sutton等人^[13]提出的T-BON.

然而, 当前的科学计算可视化方法通常基于直角坐标系, 将它们应用于如安全壳等类圆柱体的可视化时, 容易使人们对物理量在空间位置上分布情况的感知产生偏差.因此, 本文根据反应堆安全壳的外形特征, 提出了一种基于圆柱坐标系的前处理可视化网格划分方法及后处理可视化数据分析方法.在前处理网格划分中, 根据安全壳的外形特征, 提出了基于横截面和纵剖面的均匀与非均匀可视化划分方法.根据如PCCSAP-3D(passive containment cooling system analysis program, 非能动安全壳冷却系统分析程序)^[14]等安全壳冷却系统仿真计算程序对部分网格单元的特殊计算需求, 本文采用八叉树结构组织网格单元, 在提高存储空间效率的同时, 显著提高了网格单元的查找速度.在后处理数据分析中, 根据空间位置特征, 针对标量场数据提出了设定单变量的截面图方法与设定多变量的截线图方法.对于矢量场数据, 通过将颜色映射应用于几何图标, 部分地缓解了采样稠密时导致的几何图标遮挡所带来的致用户感知偏差的问题.

2 本文方法 2.1 网格划分方法

反应堆安全壳在空间上由垂直段、上封头、下封头这3部分构成, 其不包含下封头的地面以上部分形状, 如图 1(a)所示.其中, 垂直段形状为圆柱体, 上下封头形状均为半椭球.通过在半径方向、圆周方向及垂直高度方向对安全壳进行三维网格划分, 可以得到基于圆柱坐标系的网格划分方法, 其中任意网格单元均可通过半径方向索引iR、圆周方向索引jTh及垂直高度方向索引kZ来表示.安全壳划分所采用的圆柱坐标系及划分后对应的内部空间网格示意图分别如图 1(b)、图 1(c)所示.

定义半径方向划分数为nR, 圆周方向划分数为nTh, 垂直段高度方向划分数为nX, 上封头与下封头高度方向划分数均为nD, 可得安全壳高度方向总划分数为

$ nZ=nX+2\times nD $

(1)

对于上下封头, 在半径方向通过调整划分数可逼近椭球.以上封头为例, 令nD=nR-1, 则从下至上共nR-1层, 在半径方向最下层有nR个网格, 每上一层减少1个网格, 最上层为2个网格.因此, 安全壳高度方向总划分数可表示为

$ nZ=nX+2\times (nR-1) $

(2)

基于以上的划分规则, 可通过使用横截面与纵剖面可视化表示安全壳网格划分情况.通过设定圆周方向划分位置(即jTh), 可得对应于该圆周位置下的半纵剖面图; 通过设定垂直高度位置(即kZ), 可得对应于该高度位置下的横截面图.纵剖面可视化表示安全壳划分示意图如图 2所示, 其中深灰色区块表示垂直段处网格单元, 浅灰色区块表示上下封头处网格单元, 白色区块表示安全壳外部非有效网格单元区域.对应的半径方向划分数nR=5, 垂直段高度方向划分数nX=8, 上封头与下封头高度方向划分数nD=4, 安全壳高度方向总划分数nZ=16.横截面可视化表示安全壳划分示意图如图 3所示, 其中半径方向划分数nR=5, 圆周方向划分数nTh=4.

在对安全壳进行仿真计算分析时, 不同区域的重要性通常不同.若对重要的感兴趣区域增加划分数可以得到更精确的计算结果, 对于非重要区域减少划分数可以加快计算过程.基于不同区域重要性的不同, 通过在均匀划分的基础上, 根据具体计算需求对各单元在各维度上的尺寸进行自由调整, 可得非均匀的网格单元划分方法.即对于重要区域可调小单元尺寸以相对增加划分数, 对于非重要区域可调大单元尺寸以相对减少划分数.通过按网格单元尺寸等比例绘制纵剖面图与横截面图, 可直观地表示非均匀的网格划分.在图 2所示的纵剖面划分中, 左图为均匀划分, 右图为非均匀划分; 在图 3所示的横截面划分中, 左图为均匀划分, 右图为非均匀划分.

以横截面非均匀划分为例, 算法描述如下:

Step 1.  指定划分数nR、nTh、nX、nD.

Step 2.  设定将划分网格对应高度索引kZ.

Step 3.  根据nR、nTh与kZ的位置约束关系绘制均匀的横截面网格划分图.

Step 4.  根据不同区域的重要性, 调整各单元在半径方向、圆周方向的尺寸.

Step 5.  等比例地更新横截面网格划分图以直观表示此非均匀供用户分析.

Step 6.  重复Step 4与Step 5, 直至满足计算需求.

由图 2中的白色区块可知, 在圆柱坐标系下对安全壳进行三维网格划分后, 半椭球上下封头之外的单元为无效网格.在对网格划分后得到的各单元信息进行存储时, 可存储安全壳内外的全部网格或仅存储安全壳内的有效网格.按上述的三维网格划分, 可得全部网格数为

$ N=nR\times nTh\times nZ=nR\times nTh\times (nX+2\times (nR-1)) $

(3)

有效网格数为

$ {{N}_{1}}=nR\times nTh\times nX+2\times nTh\times \frac{(2+nR)\times (nR-1)}{2} $

(4)

无效网格数为

$ {{N}_{2}}=N-{{N}_{1}}=nTh\times (nR-2)\times (nR-1) $

(5)

通过直接计算白色区域, 同样可得无效网格数N₂.由公式(5)可知, 当nR或nTh较大时, 仅存储有效网格, 可获得较大的空间性能提升.

在对安全壳进行仿真计算时, 如PCCSAP-3D等仿真程序需要根据计算任务对某些特定的网格单元执行特殊处理, 为提高计算速度, 需要有高效的网格单元查找方法.需要支持的查找方式包括以3个维度索引获取特定的单个单元、以单个或两个维度索引获取满足条件的多个单元.若以Hash Table存储网格单元, 其虽能保证以O(1)时间复杂度查找单个单元, 但却不支持以单个或两个维度查找多个单元.若以三维数组形式存储全部网格单元, 则能以Random Access的方式获取任意单元, 即能保证查找单个或多个单元的时间复杂度均为O(1).但根据前述的网格划分方式讨论可知, 由于上下封头处部分单元的缺失, 存储包含非有效单元的全部单元将浪费较高比例的存储空间.本文以圆柱坐标系的3个维度进行八叉划分, 使用八叉树结构存储网格单元, 在仅需存储有效单元的同时, 获得了O(log n)的时间复杂度.

2.2 数据分析方法

根据前述的基于圆柱坐标系的网格划分方法, 对于安全壳相关的包括测量数据及仿真数据的三维标量数据, 本文提出了基于圆柱坐标系的物理数据可视化方法, 包括截面图方法与截线图方法.相比于等值面方法的基于物理量取值不变性的可视化, 本文的方法基于物理量的空间位置在单个或多个维度上的不变性, 能够更好地展示安全壳内各位置的物理量取值情况.

截面图方法包括:

● 水平截面图:特定高度z下, 物理量随角度theta及半径r变化构成的平面;

● 竖直截面图:特定角度theta下, 物理量随半径r及高度z变化构成的平面;

● 圆周截面图:特定半径r下, 物理量随角度theta及高度z变化构成的平面.

以水平截面图方法为例, 可视化算法描述如下.

Step 1.  从物理测量或仿真计算等数据源获取数据并进行规范化等预处理.

Step 2.  根据给定的高度z值, 在圆柱坐标系下从预处理后的数据中获取上下相邻的划分点高度z₁和z₂.

Step 3.  根据z与z₁、z₂的差值比, 确定上下相邻的划分点数据在后续插值计算高度z处物理量值时对应的权值.

Step 4.  沿半径方向遍历所有划分点对应的半径r:

 Step 4.1  沿圆周方向遍历所有划分点对应的角度theta.

    Step 4.1.1  在当前遍历点对应的r与theta下, 分别获取对应于z₁、z₂的物理量值v₁、v₂.

    Step 4.1.2根据Step 3中计算得到的权值, 使用v₁、v₂插值计算得到(r, theta, z)处的物理量值v.

    Step 4.1.3  根据颜色映射公式, 将物理量值v映射为对应的颜色c.

    Step 4.1.4  在高度z下, 获取在半径方向及圆周方向当前划分点对应的相邻点颜色c (多个值).

    Step 4.1.5  使用颜色c及c’在对应的(r, theta, z)下绘制三角形或四边形等几何图元.

为保证在半径r的圆弧上颜色平滑地变化, 可在圆弧上对角度theta进行插值, 即以多个三角形逼近扇形.

在将物理量值映射为RGB颜色值时, 首先将其变换到[0,4]区间内, 然后使用如下颜色映射公式计算得到对应的RGB值.

$ \langle R, G, B\rangle =\left\{ \begin{matrix} \langle 0, v, 1\rangle, \text{if} 0\le v<1 \\ \langle 0, 1, 2-v\rangle, \text{if} 1\le v<2 \\ \langle v-2, 1, 0\rangle, \text{if} 2\le v<3 \\ \langle 1, 4-v, 0\rangle, \text{if} 3\le v\le 4 \\ \end{matrix} \right. $

(6)

对于安全壳中的矢量场, 使用几何图标法, 以带箭头的线段表示矢量数据, 其中箭头方向表示该点矢量方向, 线段长度表示矢量大小.同时, 使用与标量场相同的颜色映射公式, 将矢量大小映射为箭头图标对应的颜色.

截线图方法包括:

● 单点曲线图:空间位置单个点单个物理量随时间变化曲线图, 用于分析物理量随时间的变化情况;

● 多点曲线图:空间位置多个点单个物理量随时间变化曲线图, 用于比较分析不同位置点的物理量取值;

● 单点双坐标曲线图:空间位置单个点两个物理量随时间变化曲线图, 用于比较分析两个物理量取值;

● 水平截线图:特定时刻点, theta和z固定不变, 物理量随r变化曲线图, 对应于空间位置为以安全壳竖直中轴为起点、水平向外延伸到安全壳壁面的线段;

● 竖直截线图:特定时刻点, theta和r固定不变, 物理量随z变化曲线图, 对应于空间位置为从安全壳底部延伸到顶部的线段;

● 圆周截线图:特定时刻点, r和z固定不变, 物理量随theta变化曲线图, 对应于空间位置为以r为半径的圆周;

● 经线图:特定时刻点, theta固定不变, r取最大值, 物理量随z变化曲线图, 对应于空间位置为从安全壳上封头顶点沿壁面延伸到下封头顶点的曲线, 类比于地球即为经线.

以水平截线图为例, 可视化算法描述如下.

Step 1.  从物理测量或仿真计算等数据源获取数据并进行规范化等预处理.

Step 2.  根据给定的角度theta和高度z分别获取相邻的划分点对应的坐标值theta₁、theta₂与z₁、z₂.

Step 3.  根据theta与theta₁、theta₂的差值比及z与z₁、z₂的差值比, 分别确定后续插值计算时的权值.

Step 4.  沿半径方向遍历所有划分点对应的半径r:

  Step 4.1  在当前遍历点对应的r及theta₁下, 分别获取对应于z₁、z₂的物理量值v_{1, 1}、v_{1, 2}, 并根据z₁、z₂的权值计算得到(r, theta₁, z)处物理量值v₁.

  Step 4.2  在当前遍历点对应的r及theta₂下, 分别获取对应于z₁、z₂的物理量值v_{2, 1}、v_{2, 2}, 并根据权值计算得到(r, theta₂, z)处物理量值v₂.

  Step 4.3 使用v₁、v₂并根据theta₁、theta₂的权值计算得到(r, theta, z)处的物理量值v.

  Step 4.4  在曲线图中增加从上一点(r’, v’)到当前点(r, v)的连线.

  Step 4.5 在安全壳三维模型示意图中绘制从上一点(r , theta, z)到当前点(r, theta, z)的连线.

即曲线图由两部分构成:

(1)  表示物理量值随时间或空间变化的二维曲线;

(2)  物理量值对应于空间位置点或线在三维安全壳模型中的位置示意图.

3 实验结果

基于本文提出的安全壳可视化方法, 设计与实现了包括参数管理、网格显示、计算监控、后处理显示及参数传递接口五大功能的安全壳可视化分析系统.其中, 参数管理提供参数输入、参数导入导出及参数合法性检查功能; 网格划分提供网格三维划分、网格单元选取及网格三维显示功能; 计算监控提供监控仿真计算程序运行状态、输出结果数据功能; 后处理显示提供物理量标量场、矢量场数据解析与三维绘制, 点曲线图、截线图数据解析与二维绘制功能; 参数传递接口用于将本文设计的可视化分析系统中的参数传递给相关安全壳仿真计算程序.系统中各功能模块间的关系如图 4所示, 系统界面如图 5所示.

在圆柱坐标系下的网格划分二维表示如图 6所示, 其中左侧为划分结果对应的横截面图, 右侧为划分结果对应的纵剖面图.其中, 绿色单元为被选中的在计算过程中需要特殊处理的单元.网格划分的三维表示如图 7所示.结果表明, 通过使用横截面图与纵剖面图, 可直观地表示安全壳在圆柱坐标系下的网格划分情况.

在配置为CPU:Intel Pentium G630 2.70 GHz、内存:4 GB、显卡:NVIDA GeForce 405(512 MB显存)、操作系统:Windows 8.1的机器上, 本文对比测试了采用GCC 4.9.2中std::vector、std::list及八叉树组织网格单元时, 执行查找操作的时间消耗, 结果见表 1, 其中, nFound表示进行随机查找的元素个数, 大约为总网格单元N₁的20%.实验结果表明, 当各维度网格划分数增大时, 采用八叉树组织存储网格单元, 仍能获得较高的查找效率, 保证了用户交互操作的实时性.

本文方法及Tecplot 360商用可视化软件在测量数据或仿真数据上的三维可视化显示效果如图 8所示, 其中, 图 8(a)~图 8(d)为本文方法的可视化效果, 图 8(e)与图 8(f)为Tecplot 360的可视化效果.具体为:图 8(a)为壳内温度对应的水平截面图, 图 8(b)为空气分压对应的竖直截面图, 图 8(c)为氧气浓度对应的圆周截面图, 图 8(d)为速度场, 图 8(e)为Tecplot 360对壳内温度的可视化显示, 图 8(f)为Tecplot 360对空气分压的可视化显示.从图 8中可以看出, 本文的截面图或速度场图能够较好地反映安全壳内相关物理量取值随空间位置的变化情况.

测量数据或仿真数据二维可视化显示效果如图 9所示, 其中, 图 9(a)为壳内温度多点曲线图, 图 9(b)为壳内温度与气相密度单点双坐标曲线图, 图 9(c)为壳内温度水平截线图, 图 9(d)为蒸汽分压经线图.各图的左上区域为数据项选择、坐标值选择等功能操作区, 左下区域为被选中的点或线在三维安全壳中对应空间位置的示意图.对于各绘制的曲线, 支持对线型、线宽及颜色等进行自定义设置.从图中可以看出, 截线图能够较好地反映出物理量取值随时间或空间的变化情况.

4 结束语

本文基于圆柱坐标系提出了针对于核反应堆安全壳的网格划分可视化方法与物理量数据的截线图、截面图可视化方法.其中, 网格划分方法可用于以直观的方式在圆柱坐标系下进行三维网格划分, 为相关的仿真计算程序提供参数输入.同时, 使用八叉树存储网格单元, 在保证获得较优存储空间的情况下, 获得了良好的单元查找时间性能.截线图、截面图方法能够可视化表示测量数据与仿真结果数据, 直观地展示安全壳内物理量随时间或空间的变化情况.实验结果表明, 本文的方法能够较好地展示反应堆安全壳的运转状况, 辅助研究人员更直观、高效地发现其中潜在的物理规律及反应堆运行过程中可能存在的安全问题.