基于音速喷嘴固体内部温度分布的“热效应”分析

王 超 ，王孝通 ，丁红兵 ，安海骄

(1. 天津大学电气自动化与信息工程学院，天津 300072；2. 天津市过程检测与控制重点实验室，天津 300072；3. 天津市计量监督检测科学研究院，天津 300192)

音速喷嘴广泛应用于流量测量领域[1]，由于音速喷嘴自身特殊的结构，气体在喷嘴中不断加速膨胀，并伴有较大的温降，低温气体与喷嘴固体内部结构的传热过程会对喷嘴固体内部结构和热边界层造成一系列复杂影响，称之为“热效应”[2]．音速喷嘴“热效应”的研究对其在低雷诺数领域的应用具有较大的指导意义．

对于音速喷嘴“热效应”的研究，已有不少国内外专家进行了探索．2015年，Wright等[3]根据喷嘴固体内部温度抽样误差和热边界层效应，发现了音速喷嘴对环境温度的敏感性，并给出了热边界层相关的半理论关系．在此之前，许多学者如 Мagyari等[4]、Shevchuk[5]、Вognár 等[6]对层流边界层的喷嘴固体内部温度分布做了大量的研究工作，结果表明边界层对音速喷嘴固体内部温度分布有很大的影响．然而，在他们的研究中，将喷嘴固体内部温度简化为均匀场．近几年来，音速喷嘴固体内部温度分布引起了不少的关注．最早的研究可以追溯到 1996年，Thomas等[7]通过热电偶获得了喷嘴固体内部的瞬态和稳态温度．之后，Вignell等[8]观察到喷嘴近壁面处有大量的温降．Ünsal等[9]进行了一系列对比实验，并得出了音速喷嘴固体内部温度分布的数值结果．2015年，Wright等[3]和 Ünsal等[10]提出了他们最新研究成果，修正了“热效应”下的流出系数．作为“热效应”的基本特性，音速喷嘴固体内部温度分布值得深入研究．大部分研究仍然停留在实验观察阶段，而且只能适用于喷嘴固体内部温度一致的情况下，很少有针对喷嘴固体内部温度分布规律的研究．同时，大多数研究将喷嘴固体内部结构热膨胀简化为自由膨胀修正“热效应”下的流出系数．本文根据实验和仿真研究了音速喷嘴固体内部温度的分布情况，并分析了“热效应”对流量计量的影响，为进一步研究音速喷嘴“热效应”提供了数据支撑．

1 音速喷嘴“热效应”基本原理

1.1 音速喷嘴温降过程原理

根据音速喷嘴气流相对温度与马赫数的关系，喉部处的气流温度计算式为

式中：Tcore为中心位置气流温度；T0为喷嘴入口气流滞止温度；k为气流比热容比；Ma为马赫数．

对于空气 k＝1.4，当 T0＝26.85,℃(300,K)时，喉部处 Ma＝1，Tcore＝-23.15,℃(250,K)，即空气从喷嘴入口处到喉部的最大温降ΔTcore达到 50,℃．在喉部之后，气流达到超音速状态，Ma＞1，膨胀降温过程持续进行，带来的影响将会更大．由于气流膨胀降温，气流与喷嘴固体内部结构之间存在较大的温差，喷嘴固体内部结构的热量将向气流传递，其传热过程如图 1所示．由传热学相关理论可知，其热量传递主要通过两种方式：热辐射与对流换热，其中对流换热作用占据主导地位．对于对流换热过程，喷嘴固体内部结构与气流之间的温差主要发生在紧贴喷嘴内壁的薄流体层内，这个薄层叫做热边界层．

图1 喷嘴固体内部结构与气流的传热过程Fig.1 Heat transfer process between the nozzle body and the airflow

1.2 音速喷嘴“热效应”影响规律

在临界流条件下，流过音速喷嘴的质量流量计算公式[11]为

式中：qm为质量流量，kg/s；At为喉部内截面积，m2；Cd为流出系数；C*为临界流函数；M 为摩尔质量，kg/mol；R 为通用气体常数，约为8.314,41,J/(mol·K)；T0为滞止温度，K；p0为滞止压力，Pa．

考虑到“热效应”问题对于音速喷嘴流量测量的影响，其质量流量公式可以修正为

式中：Cα为热膨胀修正系数；CT为热边界层修正系数.

气流温降会导致喷嘴固体内部温度下降，进而会影响其喷嘴固体内部结构和热边界层，并最终带来一定的测量误差．基于此，本文通过实验和仿真对喷嘴固体内部温度分布进行细致研究．

2 音速喷嘴固体内部温度分布实验

2.1 实验装置与传感器

音速喷嘴固体内部温度分布实验装置结构如图2所示．组合干燥机、高效除油器和过滤器提高了空压机的气源质量，自力式调压阀能够自动调节气体压力，起稳定压力的作用．喷嘴上游的流体处于滞止状态，该状态下流体温度和大气温度一致，喷嘴下游的流体通过缓冲罐进入到大气中．

图2 实验装置结构Fig.2 Structure of experiment device

实验所用的 3支音速喷嘴的基本结构如图 3所示，材质为 AISI 304．喉部直径 d分别为 2.15,mm、3.71,mm、5.25,mm，扩散角均为 3.5°，长度均为50.8,mm，外径均为 22,mm，即喉部直径越大，喷嘴固体结构越薄．

喷嘴入口压力 p0为 300～800,kPa，出口压力均为大气压．采用厚度为 5,mm 的珍珠棉进行绝热处理，隔绝外部自然对流．

2.2 保温层性能分析

实验采用珍珠棉作为保温材料．以喷嘴的喉部处线热流量为例，分析保温层的性能．假设实验室的环境温度为t∞、保温层外侧温度为tw1、保温层内侧温度(喷嘴固体外壁面温度)为 tw2、喷嘴固体内壁面温度为 tf．在 tf和 t∞给定的情况下，稳态时喷嘴固体内壁面到保温层内侧的线热流量、保温层内侧到保温层外侧的线热流量以及保温层外侧与环境之间的线热流量守恒，即有

式中：d为喷嘴喉部外径；d'为内径；tλ为喷嘴固体材料热导率(AISI 304热导率为 16.27,W/(m·K))；λ为保温材料热导率(珍珠棉热导率为 0.031,32,W/(m·K))；δ为保温层厚度；h为对流传热系数，且有

式中：g为重力加速度；Pr为普朗特数；v 为运动黏度．

根据保温层内侧温度tw2与喷嘴固体内壁面温度tf的差值Δt可以判断出保温材料的保温效果．以喉径为 5.25,mm 喷嘴的喉部处温度为例，取极端情况 tf为-23.15,℃(250,K)、t∞为 26.85,℃(300,K)时，未经保温处理时Δt＝37.8,℃．经过 5,mm 厚的珍珠棉保温处理后Δt＝0.2,℃．由此可见，采用5,mm的珍珠棉保温层保温性能较好．

图3 温度测量与采集系统结构Fig.3 Diagram of temperature acquisition system

2.3 温度测量与采集系统

温度测量与采集系统如图 3所示，18支精度为0.2,℃的热电偶传感器分成3层，插入直径为1.5,mm的盲孔中用于探测各点的温度．1～6号为一组，放置在距离壁面 1,mm 处探测温度；7～12号为一组，放置在喷嘴固体内部中间层探测温度；1+～6+号为一组，放置在喷嘴外壁处探测温度．1+～6+号仅存在于喉径为 5.25,mm 的喷嘴，用于两层与 3层测温点数据对比与分析，分析后采用更为可靠的两层结构对喷嘴固体内部温度分布进行研究．每个传感器在不同喷嘴中的轴向位置分布如表 1所示．2号、2+号和 8号传感器放置在喷嘴喉部的垂直位置上．盲孔处填满了导热硅胶以满足外边界条件．

温度测量与采集系统能够将热电偶温度信号进行滤波与放大，热电偶冷端温度补偿采用精度 0.5,℃的 AD590集成温度传感器．为保证信号一致性，采用信号多路复用的方式，通过数据选择器 ADG406进行多路信号的选择，多支热电偶的信号复用同一个调理电路．单个热电偶两端分别连接两个数据选择器．选择器的地址通过程序按固定的频率进行切换．通过信号调理模块，热电偶和 AD590信号分别被放大了 2,000倍和 10倍，以满足数据采集卡的量程范围．

表1 不同喷嘴测温点距离喷嘴入口的轴向位置Tab.1 Axial positions from the inlet for different nozzles mm

通过NI数据采集卡USВ-6251将这两个信号传输至电脑上，传输速率达 1.25,Мb/s．通过 LabVIEW上位机实现了冷端温度补偿，并完成数据记录与实时控制．

3 实验数据与分析

3.1 多点温度数据采集

多支热电偶探测了音速喷嘴固体内部不同位置的温度，然而喷嘴固体内部温度不仅是空间位置的函数，也是时间的函数，因此，本实验分别对喷嘴固体内部稳态与瞬态温度进行采集．由传热学知识可知，喷嘴固体结构热流量最终达到平衡状态，此时喷嘴固体内部温度为稳态温度，而瞬态温度则是某一时刻的喷嘴固体内部温度．

音速喷嘴固体内部初始温度 Ts0(Ts为喷嘴固体内部瞬态温度，Ts0为喷嘴固体内部初始温度)为24.85,℃．图 4是不同喉径喷嘴在入口压力均为400,kPa时 1～6号探测点稳态温度与 4号探测点瞬态温度变化．

由式(1)可知，理想条件下，出口处的喷嘴固体内部温度最低．然而，从图 4(a)可知，喷嘴固体内部稳态温度最低点并非在出口，而是在中间，这是由于实际扩散段处存在气动激波和回流区造成．对比图4(a)和图 4(b)可知，喷嘴固体内部温降与入口压力和喉径成正相关．

根据探测点的温度值能简单地分析出喷嘴固体内部温度分布影响规律，但是并不足以描绘喷嘴固体内部的稳态和瞬态温度分布情况．为了更全面地分析喷嘴固体内部温度分布影响规律、深入研究音速喷嘴“热效应”，结合克里金空间插值算法，利用热电偶温度信号预测整个喷嘴固体内部的稳态、瞬态温度分布云图．

图4 p0＝400,kPa时不同喉径喷嘴固体内部温度Fig.4 Body inlet temperatures of different nozzles at the same inlet pressure p0＝400,kPa

3.2 克里金空间插值算法

克里金空间插值算法是一种空间随机过程的最优线性预测算法，基于有限数据点计算出整个空间数据[12]．区域化变量 Z(x)满足二阶平稳假设和本征假设，其数学期望为m，则其协方差函数c(h)为

其变异函数γ(h)为

假设待估计点 x的临域内共有 n个空间样本点点，即 x1，x2，…，xn，其样本值为 Z(xi)，则克里金空间插值算法的插值公式为

式中λi为权重系数，表示空间样本点 xi的样本值值Z(xi)对估计值Z*(x)的贡献程度[13]．

在满足无偏性、最优性条件下，估计方差为

为使估计方差最小，令

令F对λi和μ的偏导数分别为0，得到克里金方程组为

3.3 2层与3层测温点数据对比与分析

以喉径为5.25,mm的喷嘴在入口压力为300,kPa时的实验数据为例，分析2层和3层径向位置的探测点对于克里金空间插值的影响．

采用 3层探测点作为插值基点的稳态温度分布云图如图 5(a)所示，去除喷嘴外壁的 6个探测点数据，采用两层探测点作为插值基点的稳态温度分布云图如图 5(b)所示．2层和 3层探测点的稳态温度分布误差如图 6所示，误差最大值为 0.197,℃．热电偶传感器的精度为 0.2,℃，误差最大值低于传感器的精度，由此可见，喷嘴外壁的 6个测量点数据对插值结果影响较小．同时，对于小喉径喷嘴，喷嘴喉径越小，增加一层探测点对于传感器结构设计、测量精度的要求越高．

图5 不同测温点下的喷嘴固体内部稳态温度分布云图(单位：℃)Fig.5 Steady temperature distributions in isotherm maps of nozzle body at different detection point numbers(unit：℃)

图6 2层和 3层径向位置的探测点的喷嘴固体内部稳态温度分布误差Fig.6 Error between the steady temperature distributions in isotherm maps of nozzle body at two layers and three layers of detection points

上述分析说明在两层径向位置的测量点基础上，增加一层位于喷嘴外壁的探测点对插值结果影响较小．此外，在音速喷嘴“热效应”的相关研究中，Thomas等[7]采用了 2层径向位置的探测点作为插值基点，取得了较好的实验效果．综上所述，本文采用了更为可靠的 2层结构对喷嘴固体内部温度分布进行研究．

3.4 喷嘴固体内部温度分布云图

利用克里金空间插值算法对采集到的12个探测点的稳态、瞬态温度数据进行处理，得到了喷嘴固体内部的稳态、瞬态温度分布云图．各探测点稳态温度分布云图和最大温降如图 7和表 2所示．在喉径及入口温度(T0＝24.85,℃)相同时，入口压力越大，温度分布越明显，温降幅度越大．这是由于当入口压力增大时，喷嘴雷诺数逐渐增大，结合热边界层理论，喷嘴雷诺数增大，热边界层变薄，热阻减小，气流与喷嘴固体内部结构之间的对流换热作用增强，因此气流的温降对喷嘴固体内部温度的影响增大．另一方面，当入口温度及入口压力相同时，喷嘴喉径越大，温度分布越明显、温降幅度越大．这是由于喷嘴喉径越大，固体结构越薄，面积与体积之比越大，气流的温降对喷嘴固体内部温度的影响越大，此现象也可通过喷嘴固体结构热流量变化来解释[2]．

图 8为在 d＝5.25,mm、p0＝800,kPa和 T0＝24.85,℃的条件下喷嘴固体内部稳态温度分布云图与CFD结果对比，如图所示，实验与仿真得到的喷嘴固体内部稳态温度轴向分布均呈现两端高中间低特点，与理想条件下出口处喷嘴固体内部温度最低的特点有较大差异．上述分析结合稳态温度分布云图，给出了音速喷嘴固体内部温度分布特性．而瞬态温度分布云图则可以绘制出音速喷嘴固体内部温度分布的变化过程，以d＝2.15,mm、p0＝400,kPa和T0＝24.85,℃的条件为例，给出喷嘴固体内部瞬态温度分布时序图，如图 9所示，随着实验时间的推移，喷嘴固体内部温度整体上是下降的，且轴向分布始终呈现两端高中间低的特点．

图7 不同条件下的稳态温度分布(单位：℃)Fig.7 Steady temperature distributions under different conditions(unit：℃)

表2 不同条件下喷嘴固体内部温降Tab.2 Body temperature drop of each nozzle under different conditions

图8 瞬态温度分布实验与仿真对比(d＝5.25,mm，p0＝800,kPa)Fig.8 Comparison between experiment and CFD of transient temperature distribution(d＝5.25,mm，p0＝800,kPa)

图9 瞬态温度分布时序图(d＝2.15,mm，p0＝400,kPa，T0＝24.85,℃)(单位：,℃)Fig.9 Isotherm maps of transient temperature distribution(d＝2.15,mm，p0＝400,kPa，T0＝24.85,℃)(unit：℃)

4 “热效应”对流量计量的影响

4.1 喷嘴固体内部结构约束膨胀对流动特性影响

上述实验和分析得到了音速喷嘴固体内部稳态与瞬态温度分布云图，为进一步研究音速喷嘴“热效应”提供了更丰富的数据支撑．下面从喷嘴固体内部结构热膨胀和热边界层变化两方面分析音速喷嘴“热效应”对流量计量的影响．

由金属热膨胀原理可知，喷嘴固体内部结构热膨胀会使喷嘴喉部截面积改变，进而影响喷嘴的流动特性．Wright[3]和 Вignell等[8]的相关研究中，将喷嘴固体内部结构热膨胀简化为自由膨胀修正 “热效应”下的流出系数．若喷嘴固体内部结构热膨胀为自由膨胀，则喉部的流道半径相对变化量为αΔT，因此，式(3)中 Cα＝A/Aref＝1+2αΔT，α 为热膨胀系数，ΔT＝Ts-Ts0．

然而，喷嘴固体结构两侧会被夹装，其热膨胀为约束膨胀，特性与自由膨胀有较大差异．喷嘴喉径越小，热膨胀的影响越明显，因此，以 d＝2.15,mm喷嘴为例，利用有限元分析，研究喷嘴固体内部约束膨胀对流动特性的影响．

对于各向同性的喷嘴固体材料，无内部热源，瞬态热传导方程为

式中：热扩散系数 a2＝k/(ρcp.)，k为喷嘴固体热导率， ρ为喷嘴固体密度；cp为喷嘴固体比热容；Ts为喷嘴固体内部温度；t为时间．

假定热变形是线性可逆的，不考虑不可逆的塑性和蠕变变形，线性本构方程为

式中：ε 和σ 为喷嘴固体内部应变和应力；υ为喷嘴固体泊松比，δi,j为 Kronecker增量；Tref为喷嘴固体内部参考温度；E为喷嘴固体杨氏模量．

对于d＝2.15,mm的AISI 304音速喷嘴，假定三维瞬态喷嘴固体内部态温度场为轴对称，参考温度场各点温度均为 24.85,℃，其杨氏模量和屈服强度为1.9×1011,N/m2和 2.1×108,N/m2，比热容为 500J/(kg·K)，热导率为 16.27,W/(m·K)，热膨胀系数为 1.8×10-5，泊松比为 0.29．

结合图 9中不同时刻的瞬态温度分布云图对喷嘴固体内部结构进行有限元分析．图 10为 d＝2.15,mm 喷嘴在 t＝900,s时喷嘴固体内部约束膨胀的热应力与热应变有限元分析结果，喷嘴固体结构两侧设置为不发生扩展或收缩的固定端，其喷嘴固体内部瞬态温度分布云图对应图9(f)．

对于自由膨胀，在热应力方面，材料允许自由扩展或收缩，喷嘴固体内部没有应力存在；在热应变方面，当喷嘴固体内部温度低于参考温度时，喷嘴喉径是减小的．

对于约束膨胀，应力和应变将变得更加复杂．在热应力方面，如图 10(a)所示，喷嘴固体内部热应力不为零，最大热应力出现在喷嘴固体两侧夹装处．其最大热应力为 8.2×107,N/m2，小于 AISI 304的屈服强度，不存在不可逆的塑性和蠕变变形；在热应变方面，如图10(b)所示，喷嘴喉径是增大的．

图10 约束膨胀的热应力与热应变(d＝2.15,mm，t＝900,s)Fig.10 Thermal stress and strain of constrained thermal expansion(d＝2.15,mm，t＝900,s)

从上述定性分析可知，对于d＝2.15,mm的AISI 304音速喷嘴，自由膨胀与约束膨胀引起的喷嘴喉径变化是完全相反的．下面对 d＝2.15,mm喷嘴在不同时刻喷嘴固体内部约束膨胀进行有限元分析，以研究其对流动特性的影响．图11为d＝2.15,mm喷嘴不同时刻流道半径(r)变化，图中 t＝200,s、300,s和 900,s时的喷嘴固体内部瞬态温度分布云图分别对应图9(a)、(c)和(f)，其喉部半径变化量Δr分别 2.06×10-5mm、4.35×10-5,mm 和 6.21×10-5,mm．由于喷嘴固体内部温度分布不是均匀分布，因此，以喷嘴固体喉部壁面处的温度(Ts,th)作为特征参数，则ΔT＝Ts,th-Ts0．在 t＝200,s、300,s 和 900,s 时，Ts,th分别为 22.76,℃、20.38,℃和 18.64,℃，求出其约束膨胀修正系数Cα为 1-1.84×10-5ΔT，如表 3 所示(Ts0＝24.85,℃)，而不是自由热膨胀时的 1+3.60×10-5ΔT．由于实际工程应用中夹装方式各有不同，但约束膨胀与自由膨胀确实存在较大差异，需引起注意．

图11 d＝2.15,mm喷嘴不同时刻流道半径变化Fig.11 Radius changes of flow channel for 2.15 mm nozzle

表3 约束膨胀修正系数Tab.3 Thermal deformation of constrained thermal expansion

4.2 热边界层变化对喷嘴流动特性的影响

上述分析给出了喷嘴固体内部温度变化时喷嘴固体内部结构约束膨胀对流动特性的影响，而喷嘴固体内部温度变化时，热边界层的变化也会对喷嘴流动特性产生影响．

在本课题组之前的研究中[14]，对于非绝热壁，比热容比γ＝1.4、喷嘴固体内部温度 Ts＝(0.8～1.2)T0时，式(3)中热边界层修正系数CT表示为

式中：ΔT＝Ts-T0；Re为喉部雷诺数．

如表 2所示，实验中的温降在 8.5～14.2,℃之间，由于实验条件所限，实验所用喷嘴喉径较大、喉部雷诺数值在 1.0×105附近．对于 d＝2.15,mm 喷嘴，喷嘴固体内部温度每下降 5,℃，引起的流量误差为 0.019%,．音速喷嘴在低雷诺数领域的实际应用中，喉部雷诺数的量级[2]基本为 103和 104，当喉部雷诺数为 1.0×103时，由于热边界层的变化，喷嘴固体内部每 10,℃的温降引起的流出系数的变化达到0.408%,．本文实验采用厚度为 5,mm 的珍珠棉(热导率为 0.031,32,W/(m·K)作为保温材料，通过改进绝热条件能有效减少环境温度的干扰，进而可以进行标定，修正“热效应”对流出系数的影响．

5 结 语

本文首先结合音速喷嘴固体内部瞬态温度分布测量传感器及克里金空间插值算法，对比了2层与3层测温点数据，获得了固体内部瞬态温度分布云图．接着研究了音速喷嘴固体内部瞬态温度分布特性，结果发现喷嘴固体内部温度轴向分布呈现两端高中间低的特点，实验最大温降达14.2,℃．最后分析了“热效应”对流量计量的影响．对于 d＝2.15,mm喷嘴，其约束膨胀修正系数 Cα为 1-1.84×10-5ΔT，这与自由膨胀下特性有较大差异．此外，根据热边界层理论修正式计算发现，当温降为 10,℃、喉部雷诺数为1.0×103时，热边界层变化引起的流出系数误差达到 0.408%,．因此，研究结果有望用于评估“热效应”对流量计量的影响，拓宽音速喷嘴应用范围，并提高其计量精度．

参考文献：

［1］International Standard Organization. ISO 9300. Measurement of gas flow by means of critical flow venture nozzles[S]. Switzerland：The International Organization for Standardization，2005.

［2］王 超，王 刚，丁红兵. 小喉径音速喷嘴”热效应”对流量影响的热边界层分析[J]. 机械工程学报，2015，51(16)：164-170.Wang Chao，Wang Gang，Ding Hongbing. Analysis of thermal boundary layer of thermal effects on flow characteristics of small sonic nozzles[J]. Chinese Journal of Scientific Instrument，2015，51(16)：164-170(in Chinese).

［3］Wright J D，Kang W，Zhang L，et al. Thermal Effects on Critical Flow Venturis[R]. 2015.

［4］Magyari E，Keller B. Exact solutions for self-similar boundary-layer flows induced by permeable stretching walls[J]. European Journal of Mechanics-B/Fluids，2000，19(1)：109-122.

［5］Shevchuk I V. A new type of the boundary condition allowing analytical solution of the thermal boundary layer equation[J]. International Journal of Thermal Sciences，2005，44(4)：374-381.

［6］Bognár G，Hriczó K，Bognár G. Similarity solution to a thermal boundary layer model of a non-Newtonian fluid with a convective surface boundary condition[J]. Acta Polytechnica Hungarica，2011，8(6)：131-140.

［7］Thomas M K，Richard W C. Some effects of thermal phenomena on the accuracy of critical flow Venturi based flow rate measurements[C]//ASME Fluids Engineering Conference. San Diego，CA，USA，1996：7-11.

［8］Bignell N，Choi Y M. Thermal effects in small sonic nozzles[J]. Flow Measurement & Instrumentation，2002，13(1/2)：17-22.

［9］Ünsal B，Park K A，Kaykısızlı H. Invertigations onthe effect of thermal inertia on sonic nozzle discharge coefficient[C]//International Flow Measurement Conference.Paris，France，2013：259-262.

［10］Ünsal B，Çalışkan U. Numerical assessment of discharge coefficient and wall temperature dependence of discharge coefficient for critical-flow Venturi nozzles[C]//International Symposium on Fluid Flow Measurement. Arlington，USA，2016：11-18.

［11］王 超，丁红兵，刘 琴. 正压法音速喷嘴标准装置气源系统设计[J]. 仪器仪表学报，2012，33(6)：1364-1371.Wang Chao，Ding Hongbing，Liu Qin. Design of the gas source system for sonic nozzle gas flow standard device based on positive pressure method[J]. Chinese Journal of Scientific Instrument，2012，33(6)：1364-1371(in Chinese).

［12］费 兰，杜世昌，面向零件平面度误差估计的空间泛克里金插值建模[J]. 机械工程学报，2014，50(15)：127-135.Fei Lan，Du Shichang. Flatness error estimation based on universal kriging interpolation method[J]. Chinese Journal of Mechanical Engineering，2014，50(15)：127-135(in Chinese).

［13］Burrough P A. Development of intelligent geographical information systems[J]. International Journal of Geographical Information Science，1992，6(1)：1-11.

［14］Ding H B，Wang C，Wang G. Thermal effect on mass flow-rate of sonic nozzle[J]. Thermal Science，2018，22(1)：247-262.