新近发展的MEMS技术可大面积､长时间､高精度测量壁面剪应力,深入研究近壁面湍流的物理机制及减少摩擦阻力,为各类航行器燃料消耗成本的降低及航行器性能的提高提供了可能｡然而,MEMS 壁面剪应力传感器标定方法的准确度和可信度，是应用该传感器开展实验测量的基础和前提[1, 2, 3, 4]。日本国防研究院的K.Hayashi采用2D-LDV获得明渠式水槽底部的速度剖面及雷诺应力分布，估算得到的壁面剪应力对浮式MEMS传感器进行了标定[5]；英国剑桥大学的A.De Luca采用2.5m×50cm×2.5cm (长×高×宽)风洞，在湍流环境下用槽道流法标定热膜MEMS传感器[6]；中国空气动力研究与发展中心(CARDC)空气动力学国家重点实验室采用2m×20cm×1.5cm (长×宽×高)风洞，沿程布置了16个静压传感器，以槽道流法标定了不同环境温度、不同电阻过热比的热膜MEMS传感器[7]。
本研究在充分调研国内外MEMS剪应力传感器标定试验的基础上，建立了一种新型的基于充分发展的槽道流原理的水下壁面剪应力给定装置，即精密校验水槽，通过建立静压差与壁面剪应力的关系，实现水下MEMS壁面剪应力传感器的精确标定[8]。作为一种有效的辅助手段，本文针对标定装置槽道流动进行数值仿真，开展理论分析，并通过激光多普勒测速仪(Laser Doppler Anemometer,LDA)测速实验验证，得到了可信的试验段流场特性，继而设计标定方案，综合对比了压力变化率及壁面剪应力的实验结果与数值计算和理论分析，为水下MEMS壁面剪应力传感器标定方案确定提供技术支撑。
1 标定方法简介水下MEMS壁面剪应力传感器可采用微型热敏单元方式，通过敏感强制对流换热来反映壁面剪应力情况[9]，如西北工业大学微纳米实验室研制的热膜式水下MEMS壁面剪应力传感器，如图 1所示。该传感器阵列总尺寸为长×宽×高=70mm×50mm×75μm，包括10个镍传感器单元，温度电阻系数为2800±10ppm/℃；采用Parylene C作为其防水材料[8]。在恒流驱动工作的情况下，电流加热贴附壁面的MEMS传感器单元，热能在不同流速下通过热对流的形式转移到流场中，传感器单元输出不同的电压(信号)。随着自由来流速度的提高，加热表面产生更陡的速度和温度梯度，进而加强电阻的冷却使其阻值降低，输出电压将更低。
|
| 图 1 热膜式水下MEMS壁面剪应力传感器阵列 Fig. 1 MEMS hot film shear stress sensors array for underwater applications |
标定方法的选择取决于流场条件和所使用的传感器，已发展的一些标定方法(公式)适用于恒流驱动下热膜式剪应力传感器的标定。文中讨论的是基于时间平均的静态标定法。在固定水温、相同温度电阻系数及层流环境下，MEMS壁面剪应力传感器的标定公式可用下式表示：
式中:τw是理想的时间平均壁面剪应力，e2是平均输出电压的平方，A和B是标定常数。B代表在静态环境下MEMS传感器热量的耗散，这个过程与热线标定类似。然而，层流常常满足不了要求的湍流壁面剪应力范围。Ramaprian和Tu[10]提出改良的标定方法，即在湍流环境下应用电压的高阶多项式，式(1)可重写为：
式中:C0、C2、C4、C6是新的标定常数。为提高MEMS传感器测量精度，其测量的环境(如流场温度、成分及密度等)应与标定环境尽可能保持一致。
有多种方法可获得式(2)左项τw，即给定标定使用的壁面剪应力。一类方法是通过热线、LDA或PIV等测得光滑平板某位置的边界层内层速度剖面，速度u的分布遵守“标准速度律”：
式中:k是卡门常数，取0.4；D为常数，由实验确定；ν为水运动粘性系数；ρ为水密度；τw为待求参数。需注意的是，测点距平板首部的位置不同，得到τw值也不同。而且，微观的速度剖面的精确测量对实验人员及设备提出了很高的要求。因此，应用该方法标定具有一定难度，但值得深入研究。
本文工作针对另一类宏观标定方法(所谓槽道流法)。扁平矩形槽道如图 2所示,流体从一端流入扁平矩形槽道内，忽略矩形流道侧面的摩擦力，可得到上下壁面剪应力计算公式[11]：
|
| 图 2 充分发展槽道流受力示意图 Fig. 2 Schematic diagram of force for fully developed channel flow |
式中，Δp=p1-p2，为左右端面压力差；H是槽道高度；L是槽道长度。值得注意的是，此静态标定方法得到的τw是沿槽道的平均值。而且在一定流速下，由充分发展槽道流产生的壁面剪应力沿流向是不变的。因此，该方法可简单有效地标定MEMS阵列。
2 计算模型及数值计算方法 2.1 计算模型已研发的标定装置校验水槽主体为扁平矩形截面流道(试验段)，截面尺寸为250mm×20mm(宽×高)，长度为2350mm，配备有稳定段、收缩段等以保持流动稳定，如图 3所示。设计加速流体在试验段中水速范围为0.5～5m/s，可给定最大壁面剪应力达50Pa，标定环境与后续拟开展的水中模型壁面剪应力测量试验使用环境保持一致。为充分了解试验段流动特性，这里对含稳定段、收缩段及试验段的一体化标定装置内流场进行CFD仿真。
|
| 图 3 矩形截面流道装置示意图 Fig. 3 Schematic diagram of rectangulare channel flow section |
计算网格采用结构化网格，物体对称面处的网格划分如图 4所示。本文对于近壁面流动的模拟采用近壁面函数，所以y+取30，第一层网格的高度可以按照公式(5)计算，此公式为中国船舶科学研究中心在大量计算基础上得到的拟合结果：
|
| 图 4 计算模型的网格划分 Fig. 4 Meshing of CFD model |
式中:l为模型总长(m)；Re=ρvl/μ，ρ为密度(kg/m3)，v为入口平均流速(m/s)，μ为动力粘度(Pa·s)。根据不同入口速度，由式(5)确定第一层网格高度，网格数量总计为465 000。
2.3 控制方程为了模拟湍流运动，得到时均的连续方程：
和时均的N-S方程，即雷诺平均N-S方程(RANS)：
都为时均量，u′i是脉动速度。
2.4 湍流模型
为封闭控制方程，湍流模型采用RNG κ-ε模型。该模型通过大尺度运动的计算和粘性项修正来体现小尺度的影响，而且提供了一个可以更加有效处理壁面区域的解析函数。其湍动能与耗散率方程：
式中：Gk表示由于平均速度梯度引起的湍动能，Gb表示由于浮力影响引起的湍动能；YM是可压缩湍流脉动膨胀对总的耗散率的影响。αk和αε分别是湍动能k和耗散率ε的有效湍流普朗特数的倒数。对于高雷诺数，湍流黏性系数μt=ρCμk2/ε，Cμ=0.0845。
2.5 边界条件标定装置计算域包括稳定段、收缩段及试验段，采用1/4计算模型。边界包括入口、出口和壁面等。
实验时，采用流量计监测流量，进而得到稳定段入口平均速度v，以此作为速度入口。这里针对6种流量设置了6种入口速度大小，同时根据实验中测得的温度对水流的动力粘度μ进行了设置，具体如表 1所示。水温变化对热膜式MEMS传感器的标定影响较为显著，本文仅考虑在水温基本不变的环境下研究标定方法。
| v/(10-2m·s-1) | 温度/℃ | μ/(10-3Pa·s) |
| 3.3 | 26 | 0.8578 |
| 5.0 | 25 | 0.8703 |
| 6.6 | 26 | 0.8578 |
| 8.3 | 26 | 0.8578 |
| 10.1 | 26 | 0.8578 |
| 11.7 | 26.5 | 0.8516 |
出口边界设为压力出口边界条件；壁面设为无滑移条件；对称面设为垂直于对称面的速度分量为0，平行于对称面的速度分量的法向导数为0。
2.6 差分格式及修正方法动量方程采用二阶迎风差分格式；湍动能方程采用二阶迎风差分格式；湍流耗散率采用二阶迎风差分格式。利用代数多重网格方法(AMG)加速计算收敛。压力速度耦合采用SIMPLE方法；这里采用标准壁面函数对近壁面区域进行处理。
3 计算结果与实验对比 3.1 流场特性分析实际流体进入稳定段前，湍流度较高，涡结构尺度大。经过稳定段整流之后，流动品质改善。经过收缩段之后，流体流速得到提升，并基本稳定于试验段；而流场静压沿试验段均匀变小。为便于仿真，数值计算模型简化了稳定段。图 5给出入口流速为0.100 6m/s时矩形截面流道装置流场速度分布，对应静压分布如图 6所示。
|
| 图 5 矩形截面流道装置流场速度分布云图 Fig. 5 Velocity distribution contour of rectangulare channel flow section |
|
| 图 6 矩形截面流道装置流场静压分布云图 Fig. 6 Static pressure distribution contour of rectangulare channel flow section |
试验段入口处，由于湍流度或脉动速度的存在，湍流中出现大量的掺混，这导致与主要流动方向垂直的横向流动而引起的横向动量交换，因此入口后边界层增长速度快，流体快速充分发展[12]。充分发展所需距离Le近似公式为：
式中:d为试验段水力直径；Red=ρud/μ为试验段内流体雷诺数，u为试验段平均流速(可通过流体运动的连续方程及稳定段入口流量获得)。通过估算可知，在最大实验流速下，流场在试验段入口约1m后充分发展。
3.2 LDA测速验证激光多普勒测速仪(LDA)是一种利用激光多普勒效应来测量流体或固体运动速度的仪器，主要优点是非接触测量、精度高、较高的时间和空间分辨率。设试验段入口为起始点(L=0m)，图 7给出不同入口流速下，通过LDA测得试验段中心流速u的分布及对应CFD结果。对比验证了CFD计算的准确性，且流场充分发展于试验段入口1m以后(L＞1m)，确定该区域作为标定的实验区。本文中待标定的水下MEMS壁面剪应力传感器阵列布置在该实验区中部(L=1.675m)。
|
| 图 7 试验段中心速度分布 Fig. 7 Central velocity distribution of test section |
同时，对试验段截面流场压力与壁面剪切应力进行理论计算。当雷诺数小于3×106，实验符合普朗特阻力定律(湍流光滑管)。求解压力的计算公式如下：
式(11)为普朗特公式，其中f为湍流摩擦系数；式(12)中，ξ′为沿程损失系数；式(13)为伯努利方程，其中水线高度z1=z2，u1=u2。将上式联立可以得到压力变化率Δp/L，将其代入式(4)，进而求出壁面剪应力τw。
5 标定方案设计与验证为得到压力变化率Δp/L，选取合适的压差传感器测量充分发展段压力差。试验段采用有机玻璃(亚克力)制作而成，在其下表面布置一系列间隔为100mm的压力孔来测量流向压力分布，如图 8所示。以相邻静压孔两两测量的方式测量试验段全程压力分布。图 9给出了不同流速下试验段压力的变化趋势。可见在距试验段入口1m以后，数值计算结果与实验结果吻合得较好。图 10和11给出了不同流速下试验段充分发展段压力变化率Δp/L及壁面剪应力τw的理论、数值与实验测试值的比较。可见理论解与数值解基本一致，实验求得的壁面剪应力值介于二者之间，三者吻合良好，验证了标定方案设计合理可行，满足标定要求。
|
| 图 8 试验段测压孔布置 Fig. 8 Pressure taps distribution of test section |
|
| 图 9 试验段压力分布 Fig. 9 Pressure distribution of test section |
|
| 图 10 不同水速下流动充分发展段压力变化率Δp/L分布 Fig. 10 Δp/L distribution with water velocities in fully developed flow section |
|
| 图 11 不同水速下流动充分发展段壁面剪应力分布 Fig. 11 Wall shear stress distribution with water velocities in fully developed flow section |
为辅助精确标定水下MEMS壁面剪应力传感器阵列，本文针对基于充分发展槽道流的标定方案开展了仿真分析与试验验证工作。
(1) 针对基于槽道流法的标定装置内部流场，采用RANS控制方程、RNG κ-ε湍流模型及结构网格法进行数值计算，并结合理论分析及细致的LDA测速实验验证，得出不同流速下试验段水流流场特性及压力分布；得到最大实验流速下试验段湍流充分发展所需距离约1m，可给定的壁面剪应力约26Pa；
(2) 制定MEMS传感器标定方案，选取合适的差压传感器实测试验段流场静压并确定壁面剪应力，实验结果与数值计算和理论分析的综合对比结果较好，表明仿真分析可靠地模拟了标定所需壁面剪应力给定条件，标定方案设计合理可行。
| [1] | 姜次平, 邵世明. 船舶阻力[M]. 上海:上海交通大学出版社, 1985. Jiang C P, Shao S M. Ship resistance[M]. Shanghai: Shanghai Jiao Tong University Press, 1985. |
| [2] | Marc Perlin, Steven Ceccio. Mitigation of hydrodynamic resistance[M]. Singapore: World Scientific Publishing Co. Pte. Ltd, 2015: 1-4. |
| [3] | Lennart Lfdahl, Mohamed Gadel-Hak. MEMS applications in turbulence and flow control[J]. Progress in Aerospace Sciences, 1999: 103-104. |
| [4] | Naughton J W, Sheplak M. Modern developments in shear stress measurements[J]. AIAA Journal, 2002: 515-570. |
| [5] | Hayashi K. Shear stress acting on the bed with vertical circular cylinders in open-channel flow[C]//Proceedings of the 11th International Conference on Hydrodynamics, 2014: 204-207. |
| [6] | De Luca A, Haneef I. High sensitivity single thermopile SOI CMOS MEMS thermal wall shear stress sensor[J]. IEEE Sensors Journal, 2015: 156-159. |
| [7] | Liang Jinmin, Yang Dangguo. Calibration of a thermal MEMS shear stress sensor array[J]. Symposium on Piezoelectricity, Acoustic Waves, and Device Applications, 2012, 433-437. |
| [8] | 田于逵, 谢华, 黄欢, 等. MEMS壁面剪应力传感器阵列水下标实验研究[J]. 实验流体力学, 2015, 29(2): 8-12. Tian Y K, Xie H, Huang H, et al. Calibration of MEMS wall shear stress sensors array for underwater applications[J]. Journal of Experiments in Fluid Mechanics, 2015, 29(2): 8-12. |
| [9] | Xu Y, Lin Q, Lin G. Micromachined thermal shear-stress sensor for underwater aplications[J]. Journal of Microelectromechanical Systems, 2005, 14: 1023-1029. |
| [10] | Ramaprian B R, Tu S W. Calibration of a heat flux gauge for skin friction measurement[J]. Journal of Fluids Engneering, 1983, 105: 455-457. |
| [11] | 白亚磊, 明晓. 底层隔板法流体壁面剪切力的测量研究[J]. 实验流体力学, 2011, 25(5): 85-88. Bai Y L, Ming X. The research measuring skin friction with a sublayer fence[J]. Journal of Experiments in Fluid Mechanics, 2011, 25(5): 85-88. |
| [12] | White M. 流体力学[M]. 台北: 晓园出版社, 1992. White M. Fluid mechanics[M]. Taipei: Xiaoyuan Press, 1992. |