超声波体积流量计（ultrasonic volumetric flowmeter，UFM）因其高精度、无压降、响应快速等特点而被广泛应用于计量领域. 在诸多类型的UFM中，基于渡越时间差原理设计的UFM应用最广泛，其原理是通过测量分别位于上游和下游的超声换能器激发的声波的传播时间并计算它们之间的时间差来预测体积流量 ^[1]. 声波的传播时间主要通过分析接收的声信号获得，因此，声信号的质量会在较大程度上影响时差法UFM的计量精度. 由Lynnworth等^[2-3]对UFM的总结可以看出，测量过程主要涉及声波传播与流场耦合的过程. 在流动介质存在的工作区域中，由超声换能器激发的声波会在传播过程中受到流动介质的调制. 由此可见，在实际工况下UFM管段内的流速分布会对其计量的准确性和可重复性产生极大的影响. 因此，能够预测UFM内部与流场耦合的声场对于提高UFM的计量性能具有重要意义，相关内容已有诸多学者进行了广泛研究.

流声耦合问题的研究主要基于几何（射线）声学和波动声学. 前者被用于研究移动介质中声波传播，因为在流动介质存在的工作区域中，声波传播的实际射线轨迹长度不再等于2个对射换能器之间的直线距离，会在基于时间差法原理获取的测量结果中引入非线性误差. Lygre等^[4]将简化的射线追踪应用于具有均匀流速分布的UFM仿真之中. Boone等^[5]开发新的射线追踪算法，该算法直接通过对适合任意流速分布的射线追踪式进行有限差分来实现. Iooss等^[6]使用射线追踪方法对时差法体积流量计在湍流工况下测量结果的不确定性进行研究. Kupnik等^[7]在适当的射线追踪系统的基础上，研究高温下时差法UFM的性能. Zheng等^[8-9]分析流动介质引起的测量非线性和不确定性，以提高气体超声体积流量计的计量性能.

基于射线声学的研究方法虽然能够计算出声波的传播轨迹，却无法获取一些关键的声学信息，例如声波的相位信息以及射线之间的声压，并且在计算时须将声源理想化为点声源. 这些问题可以通过波动声学的计算方法来补偿，特别是随着计算声学的发展，使用基于有限元法（finite element method，FEM）、有限差分法、边界元法等数值计算方法对复杂声场问题进行求解是较好的选择. 其中，FEM作为标准的数值计算方法得到了广泛应用^[10]. 例如，Lerch等^[11-12]使用先进有限元技术对压电换能器进行仿真. Eccardt等^[13]使用FEM模拟流动介质中的声波传播，使用的控制式是基于Pierce等^[14]推导的近似精确且适用于非均匀流体的波动式. Luca等^[15]使用线性欧拉式描述声波在UFM中的传播. 如上所述，完全依靠FEM可以分析绝大多数关于UFM的问题，但也面临一些困境，比如FEM难以对具有复杂内部结构的传感器进行精确建模. 例如，对于时差法UFM中常用的压电式超声波换能器，在实际应用中，即使在换能器设计方案已知的情况下，换能器内部压电材料的极化特性、匹配层的声学、机械特性以及黏合剂的材料属性等均难以在FEM模型中进行确切定义. 因此，若仅依靠传统的数值法构建计算模型对UFM中的声场分布进行预测会引入较大的不确定性.

为了克服这一缺点，本研究采用将实验测量数据和数值法相结合的混合计算方法，该方法能够将计算模型中具有最大建模不确定性的换能器部分通过实验数据进行替代. 通过测量的实验数据构建出发射端换能器表面的振动速度分布并作为计算模型的边界条件来预测流动介质中的声场分布，从而避免建立换能器数值模型以及求解其中的机-电场的复杂过程. 特别是在换能器内部结构未知的情况下仍可对其激发的声场分布进行预测.

湿式气体超声波体积流量计通常可以由如图1所示的计算模型进行描述. 研究对象是由课题组研发的型号为DN50的时差法超声体积流量计，观测平面是管段的中心轴截面. DN50的口径D=50 mm，由安装角度α=50°的1对超声换能器构成测量系统. 可以利用有限元软件COMSOL 5.4对该计算域进行建模，具体可以参照COMSOL自带的气动声学接口. 在计算之前须先获取发射端换能器表面的振动速度和计算域中的流场分布. 测量换能器表面上有限的测量点上的振动数据再进行数据拟合确定边界条件，利用计算流体动力学（computational fluid dynamics，CFD）模块求解充分发展的湍流得到流场分布. 最后，求解计算域中的控制式对UFM内部的声场分布进行预测，详细步骤如下.

所研究的超声体积流量计的空间尺度远大于热黏滞声学模型须考虑的狭小空间尺度^[16]，因此不考虑温度变化带来的影响并且假定理想的流动介质（例如理想空气），即不考虑介质黏度的影响，假定在一个稳定的背景流场中给定一个声速度势 $\varPhi $，介质的密度为 ${\rho _0}$，声速为 ${c_0}$. 须注意的是声粒子速度 ${{u}}$相和声速度势的关系^[17]为

(1) ${{u}} = - \nabla \varPhi .$

在非均匀流速分布区域，流速 ${{v}}$ 随着空间位置而变化，适用于该区域的线性声学方程组^[18]如下：

(2) ${{\partial \rho }}/{{\partial t}} + \nabla \cdot \left( {{\rho _0}\nabla \varPhi + \rho {{v}}} \right) = 0,$

(3) $\rho = - \frac{{{\rho _0}}}{{{c_0^2}}}\left( {\frac{{\partial \varPhi }}{{\partial t}} + {{v}} \cdot \nabla \varPhi } \right).$

将式（3）代入式（2），可以得到：

(4) $\begin{split} & \nabla \cdot \left[ {{\rho _0}\nabla \varPhi - \frac{{{\rho _0}}}{{c_0^2}}\left( {\frac{{\partial \varPhi }}{{\partial t}} + {{v}} \cdot \nabla \varPhi } \right){{v}}} \right] - \\ & \;\;\;\;\;\;\;\;\;\frac{{{\rho _0}}}{{c_0^2}}\frac{\partial }{{\partial t}}\left( {\frac{{\partial \varPhi }}{{\partial t}} + {{v}} \cdot \nabla \varPhi } \right) = 0. \end{split} $

式中： $\rho $为声扰动引起的介质密度相对于平均密度 ${\rho _0}$的变化量， $t$为时间变量. 式（4）即为如图1所示计算域的控制式，考虑到主要研究内容是换能器辐射声场在管段内的分布情况，因此没有必要使用时域上的控制式进行求解. 可以消去时间相关的计算参数，假定 $\varPhi = \phi {{\rm{e}}^{{\rm{i}}\omega t}}$， $\omega $为角频率， $\phi $为 $\varPhi $幅值，代入式（4）得到频域形式：

(5) $ \begin{split} & \nabla \cdot \left( {{\rho _0}\nabla \phi - \frac{{{\rho _0}}}{{c_0^2}}({\rm{i}}\omega \phi + {{v}} \cdot \nabla \phi ){{v}}} \right) - \\ & \;\;\;\;\;\;\;\;\;\;\frac{{{\rho _0}}}{{c_0^2}}{\rm{i}}\omega ({\rm{i}}\omega \phi + {{v}} \cdot \nabla \phi ) = 0. \end{split} $

由于切向速度对流-固耦合计算模型的影响有限，仅考虑等熵声学情况，此时流场对换能器的振动影响是可以忽略的，认为换能器表面附近为无滑移流体边界，此时振动边界附近的声粒子速度可以表示为

(6) ${{u}} \cdot {{n}} = {v_{n}}.$

式中： ${v_{n}}$为换能器表面的法向振动速度， ${{n}}$为换能器表面的外法向（指向流体介质）. 对于连续介质并且在满足无滑移边界的情况下，式（6）已被证明是有效的^[19]，因此计算域中振动边界条件可以表示为

(7) $ \begin{split} & - {{n}} \cdot \left( {{\rho _0}\nabla \phi - \frac{{{\rho _0}}}{{c_0^2}}({\rm{i}}\omega \phi + {{v}} \cdot \nabla \phi ){{v}}} \right){\rm{ = }}\\& \quad\quad{\rho _0}\left( {{v_{n}} + \frac{1}{{{\rm{i}}\omega }}{{v}} \cdot \nabla {v_{n}}} \right). \end{split} $

最终，声压可以表示为

(8) $p = {\rho _0}({\rm i}\omega \phi + {{v}} \cdot \nabla \phi ).$

使用德国Polytec公司型号为PSV-500的激光扫描测振仪对发射端超声换能器表面的振动速度进行测量，测量方式如图2所示. 图中还显示了由图形动画系统生成的换能器表面的振动模态，正常工作的换能器的振动模式是类弯曲振动，这也是圆柱形超声换能器在应用中所需的表面振动模式. 由于测量使用的激光光斑尺寸的限制，在超声换能器表面沿直径（12 mm）方向上只能容下7个测量点.

选择如下线性函数来描述换能器表面的振动速度幅值：

(9) ${v_{{n}}}(r) = {v_{\rm{o}}}\sum\limits_{m,n} {{A_{mn}}} {\left[ {1 - {{(r/a)}^{2m}}} \right]^n}.$

式中： ${v_{\rm{o}}}$为换能器中心的速度，r为距离换能器中心的距离，a为换能器的半径， ${A_{mn}}$为合适的拟合系数. 该函数仅与径向坐标r相关. 在应用该函数时需要换能器表面振动幅值在其中心处达到最大值，且随着r的增加而单调递减^[20]. 具有辊支撑约束的辐射器表面上的振动速度幅度分布可以用 $m = 1$和 $n = 2$来描述，且表面振动模式是类弯曲振动. 最终拟合结果如图3所示，拟合系数 ${A_{11}}$、 ${A_{12}}$分别为0.6565、0.3460，以此通过插值法获得的拟合数据和测量数据之间的均方根误差为0.0145，拟合程度较高. 如图3（a）所示为超声换能器表面某一直径上的振动速度拟合曲线，如图3（b）所示为整个换能器表面3D振动速度的拟合结果. 另外须注意的是，式（9）作为插值函数使用有一定的限制条件，并不适合描述所有的表面振动模式. 目前研究仅针对在良好工况下工作的超声换能器，因为在超声体积流量计的实际应用中，为了获得较好的声学信号，在谐振频率附近工作的换能器要尽可能的满足一些条件：近似弯曲振动；振动幅度是轴对称的；振幅最大位于换能器表面中心处.

使用COMSOL中的CFD流体计算模块对DN50在标定中使用的20、50、80、160、200 m³/h这些工作体积流量下的稳态的背景流场进行计算，对于口径为50 mm的管段而言，这些标定体积流量对应的平均流速分别为2.83、7.08、11.32、22.64、28.30 m/s. 虽然UFM在实际工作中通常要求安装前置直管段来保证最终的流场环境为充分发展的湍流，但出于研究的目的这里考虑安装效应对UFM计量的影响^[21]，并选择国际法制计量组织建议的测试UFM性能时使用的前置安装管段作为上游管段^[22]. 如图4所示，分别选取典型的带同心扩管的空间90°双弯管、带同心扩管的90°单弯管和直管作为研究对象.

圆管的临界雷诺数为2320，依据雷诺数的计算公式，可知介质为空气（293 K，0.1014 MPa）的50 mm圆管内的临界速度为0.69 m/s，该速度远小于研究中涉及的最小平均流速2.83 m/s. 因此，采用SST湍流模型来计算管段内的流速分布. 如图5（a）所示为图1中观测平面上的流速v分布的仿真结果，这里仅以200 m³/h体积流量且体积流量计上游为带同心扩管的双弯管的仿真结果为例. 图中，在 $ x $、 $ y $、 $ z $方向上，流速v的分量分别为 $ {v}_{x} $、 $ {v}_{y} $、 $ {v}_{z} $；l为声道中心线 l_AB的长度. 由图5（b）可以看出， $ {v}_{y} $、 $ {v}_{z} $起伏较大的区域主要位于超声波换能器的安装空腔中，即对应A1、A2这2个区域，而在区域B则趋于平稳. 即使在最复杂的流场情况下，作为主要速度分量的 $ {v}_{x} $依然远大于 $ {v}_{y} $、 $ {v}_{z} $这2个方向的速度分量，如果是充分发展的湍流，在区域B这2个方向的速度分量几乎趋于0. 从控制式推导的原理上来看，该式不适合处理有旋流场中的声场问题，但气体超声体积流量计的应用环境是低马赫数，此时可以忽略由位于换能器安装空腔处的涡流引发的声噪音问题，只考虑涡流对声波传播产生的影响.

在进行声场计算前，首先将2.1节计算获得的不同平均体积流量下的流速分布v导入计算域中作为背景流场，并设定如下的边界条件. 将1.3节中通过测量数据拟合的发射端换能器表面的振动速度作为计算域发射端换能器的振动边界. 设置模型中的管段的管壁为硬边界条件，在计算域中管段的两端（对应UFM的入口和出口）添加完美匹配层（perfect match layer，PML）来吸收声波，契合UFM安装在沿管段方向可认为是无限长的管网中的应用环境. 另外，在接收端换能器表面的后端添加PML来吸收声波，这样的设置会与接收端换能器表面实际对应的阻抗边界条件有所差异. 与PML相比，阻抗边界会透射部分声波同时反射一部分声波，这些反射出去的声波后续会以回波的形式影响接收端换能器表面所处空间的声压分布，影响程度与其传播的距离和反射声波的比率有关. 考虑到气体超声体积流量计在使用电压信号计算声波传播时间过程中一般不会用到接收端反射的声波形成的回波贡献的电信号，因此选择PML作为边界条件对声场进行计算. 虽然获取的结果不能完全反映实际管段内的声场情况，但却是研究发射端换能器在管段中激发的辐射声场特性需要的分布结果，也与后续采用的实验验证的方案相契合，另一个优点是在不知道换能器具体制作工艺和材料时，通过该方法依然可以获取进行有效分析的声场分布.

在完成以上的边界条件设置后，通过求解1.2节介绍的控制式即可对UFM内部声场分布进行预测. 与实际UFM使用的超声换能器的工作频率200 kHz相对应，在解过程中设定角频率为200 kHz的连续工作模式. 须注意的是在使用FEM求解声场时，网格的最大尺寸须小于1/4波长以保证计算结果的准确性. 所获得的不同平均体积流量下的声场分布预测结果如图6所示. 图中仅显示20、80、200 m³/h这3个体积流量下UFM上游换能器作为发射端激发的声场在观测平面上的分布. 由图6的预测结果来看，上游换能器激发的声场的主瓣会沿流场流动方向偏移，且偏移量会随着平均体积流量的增加而变大. 主瓣的偏移程度会影响作用在接收换能器表面的力的大小. 当主瓣区域逐渐偏离接收换能器表面时，换能器通过压电效应产生的输出电压信号的幅值也会逐渐减少. 与实际测量中，体积流量越大对应的输出电压信号幅值越小这一情况一致. 在获取了上游换能器作为发射端激发的声场分布后，可以对下游的接收端换能器表面所处空间的压力分布进行研究，如图7所示为前置双弯管下的预测结果. 图7（a）、（b）、（c）分别对应平均体积流量为20、80、200 m³/h时的声压分布结果. 作用于接收端换能器表面所处空间上的平均声辐射力表达式可以写作：

(10) ${F_S} = \iint_S p{ \rm{d}}S.$

式中：S为接收端换能器的表面区域. S在不同平均体积流量 ${\overline Q _{{V}}}$下的结果如表1所示. 表中， ${F_{{\rm{sp}}}}$、 ${F_{{\rm{sb}}}}$、 ${F_{{\rm{db}}}}$分别表示前置安装管段为直管、单弯管、双弯管时下游（顺流）换能器作为接收端所获得的结果， ${F_{{\rm{\_sp}}}}$为前置直管段时上游（逆流）换能器作为接收端所获得的结果. 接收端换能器表面受到的辐射力 ${F_S}(\omega )$在压电效应的作用下会转化成电能，最终转化为接收端换能器的输出电压信号 ${V_{\rm{R}}}(\omega )$. 接收过程的传递函数^[23]如下：

(11) ${t_{\rm{R}}}(\omega ) = {V_{\rm{R}}}(\omega )/{F_S}(\omega ).$

由上述可知，接收换能器表面受到的作用力大小和最终产生的系统电压信号的幅值在相同的工作频率 $\omega $下成线性关系，因此作用于接收换能器表面的平均声压可以使用由压电效应产生的电信号的幅值来表示. 因此，搭建能够测量超声体积流量计接收端电压信号幅值的实验台，将测量得到的电压信号幅值和估计的换能器表面所在空间的平均作用力的数值进行归一化处理后再进行对比，以验证所提出的声场混合计算方法.

使用如图8所示的实验测量装置来测量接收端换能器输出的电压信号幅值. 实验装置由音速喷嘴标定装置、气体超声体积流量计、示波器和由NI的 PXIe-1082&PXIe-5122构成的数据采集系统组成. 数据采集系统用于收集和记录不同平均体积流量下的下游换能器输出的电压信号幅值信息. 示波器主要用于显示输出电压信号的波形，可以直观地进行观察. 以前置管段为直管为例，分别测量不同平均体积流量时下游和上游的换能器作为接收端所输出的电压信号幅值U，其结果分别如图9（a）、（b）所示，图中显示了100个采样点的数据. 由图中的测量数据可以看出，每次测量得到的电压信号的幅值均不同且在一定的范围之内波动.

使用这100次测量的平均电压幅值作为该体积流量点下的有效电压幅值，并定义 ${V_{{\rm{sp}}}}$、 ${V_{{\rm{\_sp}}}}$分别为顺流和逆流下的有效测量电压幅值，这样可以与仿真结果相对应. 同样，定义前置管段为单弯管和双弯管安装条件下获得的下游换能器输出的有效电压幅值分别为 ${V_{{\rm{sb}}}}$、 ${V_{{\rm{db}}}}$，结果如表2所示.

对表1、2这2组数据的数值归一化后再进行比较，最终结果如图10所示. 可以看出，仿真预测和实验测量的数据整体上较一致，其中当前置管段为直管时一致性最好，双弯管时的表现最不佳，并且最大偏差均出现在200 m³/h的情况下. 双弯管作为前置管段时得到的结果不佳，是由于CFD预测的流场与实际流场之间的差异引起的，数值计算模型有一定简化，对于复杂流场环境的模拟有一定的不足. 此外，声场在200 m³/h时的偏移是最大的，结合接收端换能器指向性的影响，实际测量的电压信号幅值会小于接收到的声能量，这也是测量数据和预测数据之间的最大偏差出现在200 m³/h的原因. 通过这些数据比较和分析，验证了所提出的声场混合计算方法. 实际上在UFM应用中基本处于直管作为前置管段的工况下，在这一工况下，预测的声场结果与实验结果最为吻合.

（1）针对超声体积流量计声场研究中，因换能器数值建模而引入的较高不确定性问题，提出混合使用测量边界条件和数值计算方法建立半物理场的计算模型，对超声体积流量计内与流场耦合的声场进行求解，以减少计算中的不确定性因素并提高计算效率和预测结果的准确性. 通过该方法对不同工况下的声场进行预测，能够可视化地分析换能器激发的声场在测量空间内与流场作用后的结果，通过分析这些预测的结果可以对超声体积流量计的整体设计优化进行指导.

（2）所提出的半物理场的计算方法还可以扩展到分析其他涉及换能器的问题，不局限于超声体积流量计方面的研究. 此外，所研究的超声换能器的振动模式是类弯曲振动，因此对于非谐振频率换能器或者振动模式较复杂的情况，须进一步研究. 另外，考虑时间参数，可以进一步研究时域上的声波传播.