风洞试验是评估各类飞行器设计性能最主要、最直接的手段，即便数值模拟技术及模型飞行试验技术也在高速发展且日趋完善，进行必要的风洞试验仍是飞行器研发过程中不可或缺的环节。特别是跨超声速风洞，为多数飞行器巡航状态气动性能的试验研究提供了有效途径。

对风洞试验最重要的要求是正确模拟气流流过飞行器模型的状态并提供精确的试验数据，而优异的风洞流场品质是生产精确试验数据的前提^[1]。风洞流场湍流度作为一项重要的动态流场品质，其量值可以很大程度影响风洞试验结果及其精确度，国内外学者很早便意识到了这一点，并开始了相关研究。Jones对于湍流度对平板边界层转捩的影响进行了试验研究，试验中平板边界层的转捩通过光学方法进行捕捉，当风洞试验段流场湍流度从0.7%变化到35%时，转捩发生位置随湍流度的提高而不断向上游方向移动^[2]，类似的研究还可见于文献[3-6]中。Liu等对于喷管出口处流场湍流度对射流发展及噪声特性的影响进行了研究，发现湍流度会对喷管的气动、声学特性产生影响^[7]，国内也有学者对于湍流度对气动特性的影响进行了研究^[8-10]。

由于风洞流场湍流度会对边界层转捩及飞行器的气动、声学、热学等特性产生影响，其对于各类风洞试验的结果均较为重要。为了量化评估风洞流场湍流度，国内外学者开展了大量对于湍流度测量技术的研究。Dryden等利用热线风速仪对低速风洞流场的速度脉动、湍流度、湍流尺度等湍流相关量进行了测量研究，初步建立并完善了利用热线风速仪对低速风洞流场湍流度进行测量的方法^[11-13]。由于热线风速仪在低速不可压缩流中的响应关系式很明确(即King公式^[14])，因此低速风洞流场湍流度的热线测量方法已被研究较为透彻并广泛应用于低速风洞湍流度测试中^[15-16]。而在可压缩流中，热线风速仪的输出信号受当地流场速度、密度、温度的共同作用，其响应关系式尚不明确，因此其研究难度要远高于在不可压缩流中。以Stainback、Horstman等为代表的学者采用控制变量法对热线风速仪进行校准，通过大量的校准试验数据来求解速度、密度、总温的灵敏度系数，利用校准后的热线探针再对未知流场进行测量，进而求解湍流度^[17-20]。然而以上校准方法可能会遇到求解方程过程中系数矩阵近似奇异而难以求解的情况，且校准需大量试验对热线探针寿命不利，因此该方法并没有得到广泛应用。还有学者采用激光多普勒测速技术、瑞利散射测速技术、粒子图像测速技术等光学测量方法对可压缩流湍流度进行直接测量，取得了一定的成果^[21-24]，但由于光学测量方法普遍频响不高，难以捕捉到速度脉动中的高频成分，因此无法准确评估湍流度。

综上所述，风洞流场湍流度测量方法在低速范围内已经较为成熟，但在高速可压缩流范围内还存在较多问题。本文完善了在研的变热线过热比湍流度测量方法^[25]，引入了压力脉动项以从理论上优化湍流度求解方法，进而更加准确评估高速风洞可压缩流湍流度。在马赫数Ma=0.3~0.7进行了湍流度测量试验，对比了优化前后湍流度求解方法所得湍流度结果，并利用蒙特卡洛模拟方法对湍流度的不确定度进行了求解。结果表明优化后的湍流度求解方法所得湍流度结果与前期试验结果量值相符，随马赫数的变化趋势更加符合客观物理规律，且不确定度量值远小于湍流度量值，验证了优化后方法的可行性，为高速风洞湍流度评估提供了参考。

1 湍流度求解方法优化 1.1 湍流度求解方法

由文献[25]中的式(47)可知，恒温热线风速仪在可压缩流中的响应关系式为

$ \frac{{\Delta E}}{E} = {F_{{\rm{CTA}}}}\frac{{\Delta m}}{m} - {G_{{\rm{CTA}}}}\frac{{\Delta {T_0}}}{{{T_0}}} $	（1）

式中：E为热线风速仪输出电压；m为热线探针测量点气体质量流量；T₀为热线探针测量点气体总温；F_CTA、G_CTA分别为恒温热线风速仪质量流量、总温灵敏度系数。

式(1)左右同时除以G_CTA，定义$\theta=\frac{1}{G_{\mathrm{CTA}}} \cdot \frac{\Delta E}{E}$为响应函数，$r=\frac{F_{\mathrm{CTA}}}{G_{\mathrm{CTA}}}$为自变量，可得

$ \theta = \frac{{\Delta m}}{m}r - \frac{{\Delta {T_0}}}{{{T_0}}} $	（2）

对式(2)左右同时取均方值，可得

$ {\left( {\sqrt {\overline {{\theta ^2}} } } \right)^2} = \overline {{{\left( {\frac{{\Delta m}}{m}} \right)}^2}} {r^2} - 2\overline {\left( {\frac{{\Delta m}}{m}} \right)\left( {\frac{{\Delta {T_0}}}{{{T_0}}}} \right)} r + \overline {{{\left( {\frac{{\Delta {T_0}}}{{{T_0}}}} \right)}^2}} $	（3）

利用双曲线拟合方法对式(3)进行拟合可求解出质量流量脉动项与总温脉动项，进而求解出湍流度。但湍流度求解过程中由于忽略了压力脉动项以简化求解，因此存在一定的偏差。

1.2 湍流度求解方法优化的理论推导

为了更加准确求解流场湍流度，需对式(2)进行进一步处理。对于质量流量项，由质量流量定义可知：

$ \frac{{\Delta m}}{m} = \frac{{\Delta \left( {\rho U} \right)}}{{\rho U}} = \frac{{\Delta \rho }}{\rho } + \frac{{\Delta U}}{U} $	（4）

式中：ρ、U分别为热线探针测量点气体密度、速度。

对于总温项，由一维等熵关系式，有

$ {T_0} = T\left( {1 + \frac{{\gamma - 1}}{2}M{a^2}} \right) $	（5）

式中：T为热线探针测量点气体静温；γ为气体比热比。

对式(5)进行先取自然对数，再求偏导数的处理，可得

$ \begin{array}{*{20}{c}} {\frac{{\partial {T_0}}}{{{T_0}}} = \frac{{\partial T}}{T} + \frac{{\partial \left( {1 + \frac{{\gamma - 1}}{2}M{a^2}} \right)}}{{1 + \frac{{\gamma - 1}}{2}M{a^2}}} = }\\ {\frac{{\partial T}}{T} + \frac{{\left( {\gamma - 1} \right)Ma}}{{1 + \frac{{\gamma - 1}}{2}M{a^2}}}\partial Ma} \end{array} $	（6）

由马赫数定义，有

$ Ma = \frac{U}{{\sqrt {\gamma RT} }} $	（7）

式中：R为气体常数。

对式(7)求偏导数，可得

$ \begin{array}{l} \partial Ma = \partial \left( {\frac{U}{{\sqrt {\gamma RT} }}} \right) = \frac{{\sqrt {\gamma RT} \partial U - U\frac{{\sqrt {\gamma R} }}{{2\sqrt T }}\partial T}}{{\gamma RT}} = \\ \;\;\;\;\;\;\frac{U}{{\sqrt {\gamma RT} }}\left( {\frac{{\partial U}}{U} - \frac{1}{2} \cdot \frac{{\partial T}}{T}} \right) = Ma\left( {\frac{{\partial U}}{U} - \frac{1}{2} \cdot \frac{{\partial T}}{T}} \right) \end{array} $

（8）

将式(8)代入式(6)，可得

$ \begin{array}{l} \frac{{\partial {T_0}}}{{{T_0}}} = \frac{{\partial T}}{T} + \frac{{\left( {\gamma - 1} \right)M{a^2}}}{{1 + \frac{{\gamma - 1}}{2}M{a^2}}}\left( {\frac{{\partial U}}{U} - \frac{1}{2} \cdot \frac{{\partial T}}{T}} \right) = \\ \;\;\;\;\;\;\frac{1}{{1 + \frac{{\gamma - 1}}{2}M{a^2}}} \cdot \frac{{\partial T}}{T} + \frac{{\left( {\gamma - 1} \right)M{a^2}}}{{1 + \frac{{\gamma - 1}}{2}M{a^2}}} \cdot \frac{{\partial U}}{U} \end{array} $

（9）

用Δ代替$\partial $(均为小量)，整理可得

$ \frac{{\Delta {T_0}}}{{{T_0}}} = \alpha \frac{{\Delta T}}{T} + \beta \frac{{\Delta U}}{U} $	（10）

式中：

$ \alpha = \frac{1}{{1 + \frac{{\gamma - 1}}{2}M{a^2}}} $	（11）

$ \beta = \frac{{\left( {\gamma - 1} \right)M{a^2}}}{{1 + \frac{{\gamma - 1}}{2}M{a^2}}} $	（12）

将式(4)、式(10)代入式(2)，可得

$ \theta = \left( {\frac{{\Delta \rho }}{\rho } + \frac{{\Delta U}}{U}} \right)r - \left( {\alpha \frac{{\Delta T}}{T} + \beta \frac{{\Delta U}}{U}} \right) $	（13）

式(13)中变量过多，不利于湍流度的求解，对密度、静温项进行继续处理。对于密度项，有等熵过程压力与密度的关系式^[26]：

$ p = C{\rho ^\gamma } $	（14）

式中：p为气体静压；C为常数。

对式(14)进行先取自然对数，再求偏导数的处理，并用Δ代替$\partial $(均为小量)，整理可得

$ \frac{{\Delta \rho }}{\rho } = \frac{1}{\gamma } \cdot \frac{{\Delta p}}{p} $	（15）

对于静温项，有理想气体状态方程：

$ p = \rho RT $	（16）

对式(16)进行先取自然对数，再求偏导数的处理，并用Δ代替$\partial $(均为小量)，整理可得

$ \frac{{\Delta p}}{p} = \frac{{\Delta \rho }}{\rho } + \frac{{\Delta T}}{T} $	（17）

联立式(15)、式(17)，可得

$ \frac{{\Delta T}}{T} = \frac{{\gamma - 1}}{\gamma } \cdot \frac{{\Delta p}}{p} $	（18）

将式(15)、式(18)代入式(13)，整理可得

$ \begin{array}{l} \theta = \left( {\frac{1}{\gamma } \cdot \frac{{\Delta p}}{p} + \frac{{\Delta U}}{U}} \right)r - \left( {\alpha \frac{{\gamma - 1}}{\gamma } \cdot \frac{{\Delta p}}{p} + \beta \frac{{\Delta U}}{U}} \right) = \\ \;\;\;\;\;\left( {\frac{1}{\gamma } \cdot \frac{{\Delta p}}{p} + \frac{{\Delta U}}{U}} \right)r - \beta \left( {\frac{1}{{\gamma M{a^2}}} \cdot \frac{{\Delta p}}{p} + \frac{{\Delta U}}{U}} \right) \end{array} $

（19）

由式(19)可知，压力脉动项与速度脉动项系数的量值相近，因此这两项在响应关系式中对响应函数的贡献量相近，并不能通过简单的忽略掉压力脉动项来简化求解。

对式(19)左右同时取均方值，可得

$ \begin{array}{l} {\left( {\sqrt {\overline {{\theta ^2}} } } \right)^2} = \\ \;\;\;\;\;\left[ {\overline {{{\left( {\frac{{\Delta U}}{U}} \right)}^2}} + \frac{2}{\gamma }\overline {\frac{{\Delta U}}{U} \cdot \frac{{\Delta p}}{p}} + \frac{1}{{{\gamma ^2}}}\overline {{{\left( {\frac{{\Delta p}}{p}} \right)}^2}} } \right]{r^2} - \\ \;\;\;\;\;2\beta r\left[ {\overline {{{\left( {\frac{{\Delta U}}{U}} \right)}^2}} + \frac{1}{\gamma }\left( {1 + \frac{1}{{M{a^2}}}} \right)\overline {\frac{{\Delta U}}{U} \cdot \frac{{\Delta p}}{p}} + } \right.\\ \;\;\;\;\;\left. {\frac{1}{{{\gamma ^2}M{a^2}}}\overline {{{\left( {\frac{{\Delta p}}{p}} \right)}^2}} } \right] + {\beta ^2}\left[ {\overline {{{\left( {\frac{{\Delta U}}{U}} \right)}^2}} + } \right.\\ \;\;\;\;\;\left. {\frac{2}{{\gamma M{a^2}}}\overline {\frac{{\Delta U}}{U} \cdot \frac{{\Delta p}}{p}} + \frac{1}{{{\gamma ^2}M{a^4}}}{{\left( {\frac{{\Delta p}}{p}} \right)}^2}} \right] \end{array} $

（20）

对比式(3)、式(20)，可得

$ \overline {{{\left( {\frac{{\Delta m}}{m}} \right)}^2}} = \overline {{{\left( {\frac{{\Delta U}}{U}} \right)}^2}} + \frac{2}{\gamma }\overline {\frac{{\Delta U}}{U} \cdot \frac{{\Delta p}}{p}} + \frac{1}{{{\gamma ^2}}}\overline {{{\left( {\frac{{\Delta p}}{p}} \right)}^2}} $	（21）

$ \begin{array}{l} \overline {\frac{{\Delta m}}{m} \cdot \frac{{\Delta {T_0}}}{{{T_0}}}} = \beta \left[ {\overline {{{\left( {\frac{{\Delta U}}{U}} \right)}^2}} + } \right.\\ \;\;\;\;\;\;\left. {\frac{1}{\gamma }\left( {1 + \frac{1}{{M{a^2}}}} \right)\overline {\frac{{\Delta U}}{U} \cdot \frac{{\Delta p}}{p}} + \frac{1}{{{\gamma ^2}M{a^2}}}\overline {{{\left( {\frac{{\Delta p}}{p}} \right)}^2}} } \right] \end{array} $

（22）

$ \begin{array}{l} \overline {{{\left( {\frac{{\Delta {T_0}}}{{{T_0}}}} \right)}^2}} = {\beta ^2}\left[ {\overline {{{\left( {\frac{{\Delta U}}{U}} \right)}^2}} + } \right.\\ \;\;\;\;\;\left. {\frac{2}{{\gamma M{a^2}}}\overline {\frac{{\Delta U}}{U} \cdot \frac{{\Delta p}}{p}} + \frac{1}{{{\gamma ^2}M{a^4}}}\overline {{{\left( {\frac{{\Delta p}}{p}} \right)}^2}} } \right] \end{array} $	（23）

联立式(21)~式(23)，可得

$ \left[ {\begin{array}{*{20}{c}} {\overline {{{\left( {\frac{{\Delta U}}{U}} \right)}^2}} }\\ {\overline {\frac{{\Delta U}}{U} \cdot \frac{{\Delta p}}{p}} }\\ {\overline {{{\left( {\frac{{\Delta p}}{p}} \right)}^2}} } \end{array}} \right] = \mathit{\boldsymbol{D}} \cdot \mathit{\boldsymbol{X}} $	（24）

式中：系数矩阵D和矩阵X分别为

$ \mathit{\boldsymbol{D}} = {\left[ {\begin{array}{*{20}{c}} 1&{\frac{2}{\gamma }}&{\frac{1}{{{\gamma ^2}}}}\\ \beta &{\frac{\beta }{\gamma }\left( {1 + \frac{1}{{M{a^2}}}} \right)}&{\frac{\beta }{{{\gamma ^2}M{a^2}}}}\\ {{\beta ^2}}&{\frac{{2{\beta ^2}}}{{\gamma M{a^2}}}}&{\frac{{{\beta ^2}}}{{{\gamma ^2}M{a^4}}}} \end{array}} \right]^{ - 1}} $	（25）

$ \mathit{\boldsymbol{X}} = \left[ {\begin{array}{*{20}{c}} {\overline {{{\left( {\frac{{\Delta m}}{m}} \right)}^2}} }\\ {\overline {\frac{{\Delta m}}{m} \cdot \frac{{\Delta {T_0}}}{{{T_0}}}} }\\ {\overline {{{\left( {\frac{{\Delta {T_0}}}{{{T_0}}}} \right)}^2}} } \end{array}} \right] $	（26）

分析式(3)可知，等式右端的3个未知项(即为矩阵X中的3个未知元素)均可通过变热线过热比方法利用双曲线拟合方法进行求解^[25]；而当马赫数Ma确定时，式(24)中的系数矩阵D为常数矩阵。因此，可先利用双曲线拟合方法求解质量流量项均方值、总温项均方值及其交叉项，再根据式(24)求解某个确定的马赫数Ma情况下的流场湍流度：

$ Tu = \sqrt {\overline {{{\left( {\frac{{\Delta U}}{U}} \right)}^2}} } $	（27）

2 湍流度测量试验及结果分析 2.1 风洞及测量仪器

本次湍流度测量试验在中国空气动力研究与发展中心的探针校准风洞中进行，该风洞采用直吹射流式布局，其主要技术参数如表 1所示，结构示意图如图 1所示。所采用的主要测量仪器为IFA300型恒温热线风速仪及一支TSI单丝热线探针。

2.2 试验结果分析

在上述探针校准风洞中进行湍流度测量试验，试验马赫数范围约为Ma=0.3~0.7，每间隔约0.1取一个马赫数状态进行测量，每个状态下改变10个热线风速仪过热比，待系统稳定后进行数据采集。对热线风速仪采集到的电压信号进行以10 kHz为阈值的低通滤波处理，然后对式(3)中的自变量及函数进行求解，各个马赫数下数据如表 2所示。

利用表 2中的数据进行双曲线拟合，进而求解矩阵X中的3个未知元素，再利用式(24)~式(27)对湍流度进行求解。各个马赫数下双曲线拟合结果如图 2所示，各个马赫数下湍流度值、拟合优度值及利用优化前湍流度求解方法计算所得的湍流度值Tu如表 3所示。

由表 3中计算结果可知，拟合优度全都在0.99以上，说明利用双曲线拟合方法对优化后的湍流度求解方法所得的计算结果进行拟合的效果较好。湍流度求解方法优化前后湍流度计算结果对比如图 3所示。由图可知，优化前的湍流度求解方法所得的湍流度随马赫数的提高呈现下降趋势，这与客观物理规律及部分文献测试结果^[27-28]相悖；而优化后的湍流度求解方法所得的湍流度随马赫数的提高呈现先平稳、后上升趋势，这与客观物理规律及文献测试结果基本相符，说明优化后的湍流度求解方法能够较为准确地求得高速风洞可压缩流湍流度值。

3 湍流度不确定度的评估

为了合理评估单次湍流度测量试验结果是否能够较为准确代表风洞湍流度真实值，需要对湍流度测量试验结果的不确定度进行评估。由于湍流度是采用拟合方法进行求解得到的，因此采用传统的不确定度传递方法进行求解较为困难，为避免通过大量湍流度测量试验对湍流度测量不确定度进行评估，考虑采用蒙特卡洛模拟方法评估湍流度的不确定度^[29]。蒙特卡洛模拟方法的具体步骤如下：

1) 构造概率统计模型。在各个马赫数下，利用双曲线进行拟合时，可利用计算机进行仿真试验生成大量待拟合的散点数据，待拟合的散点数据到湍流度测量试验结果散点数据(后文简称“已知散点”)的纵向距离符合以已知散点的纵向位置为均值μ、所有已知散点到拟合曲线的纵向距离均值的1/3为标准差σ的正态分布N(μ, σ²)(即符合“3σ”原则)^[30]。

2) 模型的随机抽样。利用MATLAB软件生成符合正态分布N(μ, σ²)的1 000组随机数，将已知散点与生成的1 000组随机数叠加，即为待拟合的散点数据。

3) 确定评估值。利用双曲线拟合方法对生成的1 000组待拟合的散点数据进行拟合，得到1 000组湍流度拟合结果，并对湍流度的不确定度进行求解，即

$ u = \sqrt {\frac{1}{{n - 1}}\sum\limits_{i = 1}^n {{{\left( {{x_i} - \bar x} \right)}^2}} } $	（28）

式中：u为不确定度；i为计数变量；n为湍流度拟合结果总数，此处n=1 000；x_i为第i组湍流度拟合结果；x为1 000组湍流度拟合结果均值。

各个马赫数下湍流度求解的平均值及对应的不确定度如表 4所示。

由表 4中不确定度数据可知，在试验马赫数范围内，马赫数较低时，湍流度不确定度为0.001%量级，马赫数在0.7左右时，不确定度最高在0.01%左右，湍流度的不确定度量值远小于湍流度量值，说明单次湍流度测量试验结果即能较为准确代表风洞流场湍流度真实值。

4 结论

1) 针对在研的变热线过热比湍流度测量方法忽略压力脉动项的问题进行了完善，从理论上优化了利用恒温热线风速仪对可压缩流湍流度进行测量的方法。

2) 在马赫数Ma=0.3~0.7进行了湍流度测量试验，利用优化后的湍流度求解方法对湍流度进行了求解，并与优化前的求解结果进行了对比。结果表明优化后的湍流度求解方法所得的湍流度更加符合客观物理规律及部分文献测试结果，验证了优化后方法的有效性。

3) 利用蒙特卡洛方法对湍流度的不确定度进行了求解，不确定度量值远小于湍流度量值，验证了湍流度测量结果的稳定性。