热量表流量计量稳定性的流场动态数值分析

官燕玲，江 超，肖 洋,2，张宇昊,3

（ 1. 长安大学环境科学与工程学院，陕西 西安 710054； 2. 成都基准方中建筑设计有限公司西安分公司，陕西 西安 710054；3. 信息产业电子第十一设计研究院科技工程股份有限公司西安分院，陕西 西安 710100 ）

由于供暖分户计量的要求，户用热量表成为集中供暖系统用户端常用的计量设备．在运行中，由于用户端流量可自行调节，热量表会经常处于变流量工作状态；同时，流体在基表中要绕过立柱反射体，呈现非定常性；再有，热量表前、后接管的直管段长度也会对流场产生一定的影响．热量表在非稳定状态下计量的可靠性会直接影响热计量收费的合理性，最终影响集中供暖系统的运行节能效果．因此本文提出对一款常用的户式热量表，进行流量计量可靠性的流场动态数值研究．

热量表主要由流量传感器，配对温度传感器和计算器三部分组成．流量传感器是热量表最主要的部件，按其测量原理可以分为机械式、超声波式和电磁式三类．由于超声波的测量优势，被广泛应用，并以时差法最为普遍[1-2]．立柱式时差法超声波热量表在基表上装有一对立柱反射体，其柱体上端相对斜面作为反射面，这是目前用的较多的一种形式，见图1．

图1 热量表基表几何图Fig.1 Geometric figure of heat meter base table

从该表的工作原理可以看到，其流量计量是基于超声波传播所经过的流场轴线上的线速度，通过修正得到流通断面上的平均速度，再由流通断面面积计算得到流量．在热量表的生产中，由线速度转换为面速度的修正值K，是通过实验测量标定的．这个修正值K的可靠性与基表的水流特性有直接关系．近年来，为了进一步提高热量表的计量精度，在超声波反射装置、超声波声路的改进方面人们作了大量的工作[3-4]．但是，针对热量表的非稳定工作条件以及安装条件的影响，从流场特性方面对其流量计量的可靠性进行研究的还很少．

本文提出对一款DN25的立柱反射体时差超声波热量表，应用大涡模拟数值方法，基于基表内的水流特性，围绕流量计量的可靠性，分析了流场从初场到充分发展的动态过程，以及流场的非定常性、流量的改变以及安装直管段长度对K值的影响．本文旨在为热量表运行计量的可靠性给予评价，为热量表的改进、管理及安装提供帮助．

1 数值建模

1.1 计算模型的选择

该问题为三维非定常的湍流问题．湍流数值模拟方法有3种，即直接数值模拟(DNS)、大涡模拟(LES)和雷诺时均方程模拟(RAMS)．直接数值模拟可以给出所有湍流脉动量，但对计算机的性能要求很高，很难达到；雷诺平均数值模拟方法只能给出统计平均量；大涡数值模拟方法可以给出大于惯性区尺度的脉动信息，特别是大尺度脉动信息．

该基表几何形状见图1．基表流道前后有两个圆柱体，两圆柱体上端相对有45°的斜面为超声波的反射面，两柱中间为超声波传播区．基表内流体依次流过这两个柱体，形成受限空间两个圆柱体的绕流流场．这是一个复杂的非定常流场，其流场特性受到很多因素的影响．关于受限空间绕流问题，文献[5-6]进行过实验研究，认为用LES湍流数值方法能更可靠地对该问题进行仿真．

根据以上分析，综合考虑所研究的流场的特殊性，决定选用LES湍流数值方法进行问题的研究．

1.2 几何模型的建立

1.2.1 超声波热量表的几何模型

本文研究的时差法超声波热量表基表几何形状如图1所示．在基表过流直管段内，装有一对直径为12 mm的圆柱反射体；圆柱上端相对45°斜面为超声波反射面，反射面中心在过流直管的轴线上；基表总长为160 mm，入口及出口内径为28 mm；基表缩径断面通道，总长为58 mm，内径为D为17 mm；图中a、b所示为一对换能器安装位置．

1.2.2 计算几何模型

按照图1的内廓几何形状，绘制计算几何模型，图2为其外形图．三维坐标原点设在入口的中心点上；X轴为轴向坐标，与过流通道轴线重合，方向由入口指向出口(见图1)，计算范围为0～160 mm；Y轴为竖向坐标，指向上，计算范围为-14～25 mm；Z轴为径向坐标，指向前，计算范围为-14～14 mm．

图2 计算几何外形图Fig.2 Computational geometry outline drawing

1.3 网格划分

将整个基表作为计算区域，由于该基表形状复杂，因此选用四面体网格对全场进行网格划分，网格数约为5×105个，离壁面第一排网格的无量纲距离y＋小于1，可以满足计算精度要求[7]．图3为通过轴线的竖向断面和水平断面网格分布图．

图3 网格划分断面图Fig.3 Cross-section diagramof meshing

1.4 数值计算方法

本文选用LES湍流数值模型．采用Deardorff的Box的滤波方法；亚格子应力模型为涡粘模型，其中亚格子湍流粘性系数选择Smagorinsky-Lily模型；Smagorinsky常数Cs为0.1．对流项采用二阶中心差分格式；非稳态项采用二阶隐式格式；算法采用SIIVIPLE算法．边界条件为接管入口为均匀速度边界条件；接管出口为自由出流(outflow)条件；壁面均为无滑移边界条件．

2 基表内有限空间绕流流场的特性分析

2.1 水流从初场到充分发展的过程

本节通过水流从初场到充分发展的过程的分析，来研究实际运行中由于流量发生变化可能对流场稳定性带来的影响．

设定时间步长为0.000 05 s，约为超声波单射程所需时间(超声波在水中的传播速度约为1 500 m/s[8])．入口流速为1.5 m/s(流量为3.32 m3/h)，流体密度为998.2 kg/m3，雷诺数Re为6.8×104．计算的问题是，以水为流动介质，在基表内的流道内，初始各点流速均为1.5 m/s，突然出现前后两个立柱及周围受限空间的流动，观察流场从最初到充分发展的过程．

图4 0.000 05 s时间步长瞬时线平均速度计算曲线Fig.4 The instantaneous average linear velocity of 0.00005 s time step

图4 所示为计算过程曲线．图中的横坐标为计算时间步，纵坐标为流速，图中曲线为一对反射面中心之间的瞬时线平均速度(即超声波传播路径上流体的平均速度)．图中看到，从第1 000步开始，曲线基本趋于平稳．可以认为流场从初态到充分发展只用了1 000步，即用了0.05 s的时间．由于流场从初始状态到充分发展的时间很短，这期间流场对流量测量精度的影响就很小，因此，在评价流量计量的准确性时，可以不考虑因为流量的瞬间变化所带来的测量误差．

图5给出了从开始到充分发展这段时间的流场速度分布图．从图中可以看到流场的形成过程和两个柱体绕流所形成的流场特征．图5共取了5个时刻的计算结果，每个时刻截取y为-5 mm和0.0 mm的两个水平断面；y=0.0 mm即为轴心高度，是超声波传播路径所在的断面．

图5 充分发展流场形成过程流场断面流速等值分布图Fig. 5 The contour maps of cross-section velocity in fully developed flow forming process

图5 可以看到充分发展流场的形成过程，y=0 mm断面中心主流区随着时间的推移由第一个柱体后面向第二个柱体靠近，到1 000步时达到第二个绕柱，即流场达到稳定．此时，在超声波传播通道中，形成紊流的速度分布特征，即管道断面速度分布较均匀．图中还可以看到，到1 000步流场达到充分发展时，在小断面通道前后之间(两个绕柱之间)，流场断面分布有所不同，向第二个柱体方向，主流区逐渐缩小，很明显，这样的变化会影响流量计量的精度．

2.2 流场的非定常性

从图5中的1 000步（流场已达到充分发展）可以看到，第一个柱体绕流和第二个柱体绕流完全不同．由于第一个绕流之后要进入小断面通道，所以流体在柱体后只有小的分离，在y=0 mm断面可以看到有小的对称蜗旋，在y=-5 mm断面基本就不存在蜗旋了，在进入缩径通道时，就已形成管内紊流断面分布状态，因此第一个绕流对后面的测量区域影响不大；但在第二个绕流后有明显的分离，通过大量的计算，可以看到类似涡街现象，但这已在测量区域之外了，应不会影响流量的测量．

3 流量计量的可靠性分析

3.1 计算流量及时间步长的选择

2007年建设部颁布了行业标准《热量表》(CJ 128-2007)(替代CJ 128-2001)[9]，2001年国家质量技术监督局颁布了国家计量检定规程《热能表》(JJG225-2001)[10]．标准[9]规定了热表的常用流量qp及常用流量与最小流量 qmin之比的最小限值．公称直径25 mm的常用流量 qp为3.5 m3/h，常用流量与最小流量之比为50，从而得到最小流量 qmin为0.07 m3/h．最大流量 qmax为水流经热量表时在短时间内正常运行的极限流量，由厂家提供，该热量表的最大流量maxq 为7.0 m3/h．

文献[9-10]给出了基本相同的流量检测点，文献[9]规定的出厂检验的三个测量点为： qmin≤q≤1.1qmin、0.1qp≤q≤0.11qp、0.9qp≤q≤1.0qp；型式检验规定的五个测量点为： qmin≤q≤1.1qmin、0.1qp≤q≤0.11qp、0.3 qp≤q≤0.33qp、0.9qp≤q≤1.0qp、0.9qmax≤q≤1.0qmax．

根据以上规定，本文选取了对应型式检验的5个流量进行了计算分析(其中包括了出厂检验的三个流量)，流量分别为：6.65、3.325、1.05、0.376 8和0.075 m3/h，详见表1．通过多工况的计算比较，0.000 05 s步长与0.000 5 s 步长的计算结果基本相同，考虑机时限制，确定计算时间步长为0.000 5 s．

3.2 计算结果及分析

3.2.1 不同流量的线速度的稳定性分析

由于超声波传播路径是在反射面中心两点之间，这之间的线速度的稳定性是影响热量表计量精度的一个重要因素．图6为反射点之间的线速度的计算曲线．

表1 不同流量的线速度稳定性分析数据Tab.1 The stability analysis data of linear velocity in different flow

8.6

图6 0.000 5 s时间步长线速度计算曲线（纵坐标为速度m/s，横坐标为时间步）Fig.6 Linear velocity calculation curve of 0.000 5 s time step (ordinate for velocity m/s，abscissa for time step)

图6(a)为瞬时线平均速度(两个反射面之间轴心速度平均值)计算曲线，可以看到，对应不同流量，其流场达到稳定的步数有不同，流量由大到小其稳定步数分别为100、200、300、400、2 200步,即流量越小，流场从初始状态达到稳定所需要的时间就越长，但对应最小流量，稳定时间也只需要1.1 s．

图6(b)为各流量工况流场稳定后300步的瞬时线平均速度的计算曲线，图中的水平线为瞬时线平均速度的时均值．通过这个曲线可以看到流场的非定常性引起的线平均速度的波动．下面通过误差值的分析，比较不同流量大小瞬时线速度平均值的稳定性．

计算数据分析见表1．表中瞬时线平均速度均方差的计算公式为：

式中： vi为瞬时线平均速度，m/s；v为瞬时线平均速度时均值，m/s；n为瞬时值个数(同计算步数)．

表1中给出了瞬时线平均速度时均值v、相对v的最大绝对误差和相对误差、瞬时线平均速度均方差σ．表中看到，对于相对v值最大绝对误差，随着流量的减小而减小；对于相对v值最大相对误差，各流量下的都很小，没有明显的变化规律；对于瞬时线平均速度的均方差σ，随着流量的减小而减小，所有流量，最大为0.006 59，最小为0.000 04．

通过以上数据分析看到，无论流量大小，该流道的线速度均比较稳定．

3.2.2 流量修正系数K值的稳定性分析

时差式超声波热量表，测出的是超声波传播路径上流体的平均流速，并非计算流量所需要的管道横断面上的平均流速，因而在计算流量时要采用K系数进行修正．K系数定义式为[1-3]

式中：v为超声波传播路径上的流体平均流速，m/s；u为管道横断面上流体平均流速，m/s ．

K的稳定性是衡量不同超声波反射装置及超声波声路优劣的主要依据．要求K值在不同的流量下能够较稳定地保持为某个定值．

如上计算的瞬时线平均速度即为超声波传播路径上的流体平均流速，再根据流量可以得到流道断面流速，即可得到各流量下的K值．K值的稳定性分析见表2．

表2 不同流量下的K值的稳定性分析数据Tab.2 The stability analysis data of the value of K in different flow

同表1，表2中瞬时线平均速度时均值v是对应各流量流场稳定后的300步的平均值，表中流量对应型式检验的5个流量，其中带星号“☆”的流量同为出厂检验的3个流量．从表2看到，对应型式检验5个流量下的K的平均值为1.076 8，各流量下的K值与之相比，最大绝对误差为0.127，最大的相对误差为11.8%，都出现在最小流量工况下；其它流量的K值，其误差则随流量由大到小而依次减小，次小流量工况绝对误差为0.004，相对误差为0.4%，最大流量工况的绝对误差为0.054，相对误差为5%．

若以出厂检验的3个流量的K值的平均值作为比较，如表2，K的平均值为1.104，此时，最大绝对误差、相对误差仍然是出现在最小流量工况下，分别为0.10和9.0%；其他两个流量工况，流量大，则绝对误差和相对误差大，次小流量工况绝对误差为0.023，相对误差为2.1%，最大流量工况的绝对误差为0.077，相对误差为7%．

以上分析表明，无论是型式检验的5个流量还是出厂检验的3个流量工况，K值的不稳定性都是存在的，最大误差均出现在最小流量工况下，其次是出现在最大流量工况下．

4 热量表前后直管段长度对流量计量的影响

热量表在系统中的安装应考虑前后一定的直管段长度，一般产品说明书上会给出最小直管段安装长度，这是考虑热量表前后流场的不稳定会造成流量计量的误差而提出的基本要求．为了说明直管长度对流量计量的影响，本节在前文的基础上，对同一款热量表进行了计算分析．

4.1 计算几何模型

本节在图1所示几何模型的基础上(缩径管直径为17 mm)，加长热量表前后直管段尺寸,并带弯头，见图7所示．

图7 热量表实际安装距离示意图Fig.7 The schematic diagram of the actual installation distance of the heat meter

图中，L1表示热量表上游直管段长度、L2表示热量表下游直管段长度．接口直径D为28mm，根据直管段不同长度设三种工况：工况一，L1=10D (280 mm)、L2=5D (140 mm)；工况二， L1=10D (280 mm)、L2=2.5D (70 mm)；工况三，L1=5D (140 mm)、L2=5D(140 mm)．

将整个基表以及弯头区域作为计算区域，由于该基表形状复杂，因此选用四面体网格对全场进行网格划分，网格数约为2×106个．

4.2 定解条件及模型

入口边界设置为均匀速度入口(velocity inlet)，速度取值大小依次为0.2、0.4和1.5 m/s．出口边界设置为自由出流(outflow)．其他条件与前文一致．同样采用大涡模拟方法．

4.3 计算结果及分析

4.3.1 不同直管段长度流场特征

图8所示为0.4 m/s入口流速下的，不同直管段长度，在流场稳定条件下的，Y=0和Z=0断面的速度等值线图．

图8 0.4 m/s流速下热量表管道内流速等值分布图Fig.8 The contour maps of heat meter pipe flow velocity distribution in 0.4 m/s velocity

工况一和工况二，热量表出口直管段长度不同，分别为5D和2.5D，但入口直管段长度相同，均为10D．从图中可以看到，入口直管中，前段有因为弯头引起的涡流，但接近热量表时涡流已消失，后面接出管长不同，对前面接入管中的流场没有明显的影响；出口接管中，无论是5D还是2.5D管长，流场都不稳定．

工况一和工况三，热量表入口直管段长度不同，分别为10D和5D，但出口直管段长度相同，均为5D．从图中可以看到，在入口直管段中，对因弯头引起的涡流，工况一在热量表入口处已消失，但工况三其涡流在热量表入口处没有消失．

下面将通过数据分析，说明接管直管段长度对流量计量稳定性的影响．

4.3.2 热量表前后直管段长度对K值稳定性的影响

分析方法同上，由数值计算得到每个工况基表过流断面轴心的瞬时线平均速度的时均值，以及该断面上流体平均流速，根据式(2)得到每个工况的K值．对应每个安装工况分别有三个流速，计算其K的平均值，从而得到相对误差和均方差，由此来评价各安装工况流量计量的稳定性．

计算分析结果见表3．由表3可以得到以下结论．

表3 K 值稳定性分析Tab.3 The stability analysis of the value of K

对这三个安装工况，K值均随着流速的增加而减小，1.5 m/s流速下，工况一至工况三K值分别为1.026 7、1.033 5、0.974 4．

工况一和工况二相比较，K值的均方差，工况一为0.018 1、工况二为0.017 8，两者差别不大，说明当热量表出口的直管段长度由5D减小到2.5D时，对流量计量的稳定性影响不大．

工况一和工况三比较，K值的均方差，工况一为0.018 1，工况三为0.050 1，很明显工况三远大于工况一，说明热量表上游的直管段长度由10D减小到5D时，会明显地影响流量计量的稳定性．

5 结论

(1) 对于DN25的立柱反射体超声波热量表，基表水流从初始状态到稳定状态所用时间只有0.05 s，说明因流量的改变要达到一个新的稳定流场所需时间很短．因此，在评价流量计量的准确性时，可以不考虑因为流量的瞬间变化所带来的测量误差．

(2) 在超声波传播区域内，无论流量大小，超声波传播路径的瞬时线平均速度随时间均比较稳定．

(3) 对应不同的流量，该热表的K值的不稳定性是存在的．对应型式检验流量和出厂检验流量两种情况，与各平均K值相比，两者最大误差均出现在最小流量工况下，其次是出现在最大流量工况下．对应型式检验流量，K值的最大误差为11.8%，其次为5.0%；对应出厂检验流量，最大误差为9.0%，其次为7.0%．

(4) 热量表进、出口接管的直管段长度对流量计量稳定性有影响，特别是进口直管段长度对其影响较大．通过对三种流量工况计算，进口直管段分别为10D和5D时，其K值的均方差分别为0.018 1和0.050 1．在实际工程中应保证足够的直管段长度．

References

[1] SANDERSONM L, YEUNG H. Guidelines for the use of ultrasonic non-invasive metering techniques[J]. Flow Measurement and Instrumentation. 2002, 13( 4): 125-142.

[2] 强发红, 毛协柱. 时差法超声流量计的应用技术[J].石油化工自动化, 2001(1): 60-64.QIANG Fahong, MAO Xiezhu. Application technology of ultrasonic flowmeter under time-difference method[J].Automation inpetro-chemical industry, 2001(1): 60-64.

[3] 贺胜, 彭黎辉, 仲里敏. 基于CFD的超声波流量计最优声道位置研究[J]. 仪器仪表学报. 2009, 30(4): 852-856.HE Sheng, PENG Lihui, ZHONG Limin. Computational fluid dynamics based sound path optimization for ultrasonic flow meter[J]. Chinese Journal of Scientific Instrument, 2009, 30(4): 852-856.

[4] 沈芳. DN25户用热量表基表性能的研究[D]. 济南: 山东大学, 2011.SHEN Fang. Study on the performance of DN25 household heat meter[D]. Jinan: Shandong University, 2011.

[5] 贾晓荷. 单圆柱及双圆柱绕流的大涡模拟[D]. 上海:上海交通大学, 2008.JIA Xiaohe. Large eddy simulation of Single cylinder and double cylinder of streaming[D]. Shanghai: Shanghai Jiao Tong University, 2008.

[6] 陈晓春. 基于并行计算的大涡模拟方法及其工程应用基础研究[D]. 西安: 西安建筑科技大学, 2004.CHEN Xiaochun. Large eddy simulation method based on parallel computing and its engineering application in basic research[D]. Xi'an: Xi'an University Of Architecture And Technology, 2004.

[7] 张兆顺, 崔桂香, 许春晓. 湍流大涡数值模拟的理论与应用[M]. 北京: 清华大学出版社, 2008.ZHANG Zhaoshun, CUI Guixiang, XU Chunxiao. Theory and application of numerical simulation of turbulent large eddy[M]. Beijing: Tsinghua University press, , 2008.

[8] 岑敏锐. 超声波在液体中的传播速度与温度的关系[J].物理实验. 2008, 28(5): 39-41.CEN Minrui. The relationship between velocity and temperature in liquid of ultrasonic[J]. Physical experiment.2008, 28(5): 39-41.

[9] 中华人民共和国建设部. CJ 128-2007 热量表[S]. 北京: 中国标准出版社, 2008.The ministry of construction of the People's Republic of China. CJ 128-2007 Heat meter[S]. China Standard Press,Beijing, 2008.

[10] 国家质量技术监督局. JJG225-2001 热能表[S]. 北京:中国标准出版社, 2001.National Bureau of quality and technical supervision.JJG225-2001 Heat meter[S]. Beijing: China Standard Press, 2001.