CN104697705A - 一种油膜干涉表面摩擦应力测量方法 - Google Patents

Abstract

本发明公开了一种油膜干涉表面摩擦应力测量方法，其具体包括以下的步骤：步骤一、在t1时刻和t2时刻之间，获取油膜粘性系数 随时间的变化关系，并测量得到t1时刻和t2时刻的油膜厚度；步骤二、设定初始摩擦应力，并根据步骤一得到的油膜粘性系数随时间的变化关系，t1时刻的油膜厚度，数值计算薄油膜方程，并采用显式时间步推进到t2时刻，从而计算得到t2时刻的油膜厚度；步骤三、比较计算得到的t2时刻的油膜厚度与测量得到的t2时刻的油膜厚度，并根据比较结果对摩擦应力进行修正，重新从t1时刻开始显式时间步推进到t2时刻，当计算所得的油膜厚度与测量得到的油膜厚度一致时循环结束，此时的即为求解的摩擦应力。通过上述方法得到油膜干涉表面摩擦应力。

Description

一种油膜干涉表面摩擦应力测量方法

技术领域

本发明涉及试验流体力学技术领域，本发明公开了一种油膜干涉表面摩擦应力测量方法。

背景技术

摩擦阻力在飞行器的设计中占有非常重要的地位。对于商业飞机而言，在巡航状态时摩擦阻力几乎占到总阻的40％～50％，有评估认为减少1％的阻力大约可以增加10名乘客[1]，无论是从经济效益还是节能环保方面来说，减阻都是非常重要的；对于高超声速带动力飞行器而言，摩擦阻力占到了总阻的30％，直接关系到飞行器是否可以获得正推力以及设计方案是否成功；对于超高声速再入滑翔类武器而言，摩擦阻力的大小直接影响其射程，关系到武器设计的成败。

传统的模型表面摩擦应力测量技术有摩阻天平、Preston管等方法，这些方法需要在实验前进行仔细的校准，在实验后对测量数据进行修正，因此整个实验非常繁复，测量误差较大。尤其是在高超声速风洞实验中，由于来流速压较小，模型表面的摩擦应力只有几帕斯卡到几十帕斯卡，基于应变原理的传感器很难测量准确。而Preston管等一些方法需要用探针测量靠近壁面的压力分布，探针直接插入流场，对探针附近的流动不可避免的产生干扰。因此需要发展非插入、更准确的测量方法。近二十年来得到大力发展的测量方法主要有MEMS、剪切液晶和油膜干涉(Oil Film Interferometry)。

其中油膜干涉方法能适应较大速度范围的流动情况。它是一种非插入、不需要校准、能应用于曲面的摩擦应力测量方法，分辨率较高。除此之外，它对模型的加工要求不高。油膜干涉方法的主要原理为，覆盖在模型表面的油膜，受边界层的剪切作用会形成楔形形状，其厚度随时间的变化率主要与表面摩擦应力有关。油膜的厚度通常在微米量级，对流场基本没有干扰，可以通过光学干涉的方法准确测量。

当入射光(波长为λ)照射模型表面时，分别在空气/油膜界面和油膜/模型界面上发生反射。两道反射光在无穷远或相机透镜的聚焦下相互干涉，光强相长或相消，表现为黑白相间的条纹。

1962年Squire提出的薄油膜公式给出了油膜厚度与摩擦应力的关系(见公式1)。Tanner等人首次将其应用于干涉法摩擦应力测量。Zilliac提出用油滴代替油线，将之应用于三维流动，用前缘附近两三道条纹求解摩擦应力大小，用油滴流动方向指示摩擦应力的方向。

$\frac{\∂h}{\∂t}+\frac{\∂}{\∂x}\left(\frac{{\mathrm{\τ}}_{w,x}{h}^2}{2\mathrm{\μ}}\right)+\frac{\∂}{\∂z}\left(\frac{{\mathrm{\τ}}_{w,z}{h}^2}{2\mathrm{\μ}}\right)=0---\left(1\right)$

传统的油膜干涉摩擦应力测量技术主要应用于低速和亚跨声速范围，在整个吹风过程中油膜的粘性系数变化较小，通常被当成一个常数进行处理。只需要精确测量模型表面目标区域内一点在某时刻的温度，即可获得用于求解的粘性系数。有些解算方法虽然考虑了油膜的粘性系数随时间的变化情况，但在实际应用中需要连续测量油膜的温度并进行积分。这一过程非常繁琐，且实现起来很困难。如Garrison和Ackman提出的迭代求解方法，计算公式如下：

$C_{f,i+1}^{1/2}=\frac{\underset{0}{\overset{x}{\∫}}{(n/C_{f,i})}^{1/2}\mathrm{dx}}{{\mathrm{hn}}^{1/2}\underset{0}{\overset{t}{\∫}}\frac{q}{\mathrm{\μ}}\mathrm{dt}}---\left(2\right)$

其中C_f是摩擦应力系数，n是流线间距离，q是来流动压，i是迭代步数，迭代的初始值使用下式。

$C_{f,1}=\frac{\mathrm{\μx}}{\mathrm{qht}}---\left(3\right)$

(2)式使用的是单时间步求解方法。对于常规直加热式高超声速风洞而言，风洞启动时间较长，难以准确确定油膜变形的开始时间，且由于启动和关车阶段气流冲击载荷较大，使用单时间步求解方法误差较大。因此，常规高超声速风洞宜采用双时间步求解方法，即在吹风过程中，在t1和t2时刻分别记录油膜干涉图像。另外使用红外热像仪记录t1到t2时刻油膜温度的变化情况。

由于油膜的厚度在微米量级，可以认为模型表面的温度就是油膜的温度。在高超声速风洞实验的吹风过程中，油膜的温度随空间和时间均是变化的。对于使用油滴法的求解区域(5mm×5mm)，其随时间的变化更为明显。例如，来流M_∞＝7.94的平板模型，在实验时间内，油膜温度变化约为0.55℃/s。显然，当双时间步的时间间隔较长时，由于温度的变化引入的误差将非常明显。高超声速风洞实验的摩擦应力系数较小，以平板模型为例，距前缘100mm处约为0.00036，在较短的时间内(<1s)作用在油膜上的效果(油膜厚度的变化)不明显，由相应的干涉图像计算油膜厚度时相对误差较大。减小油膜的粘性系数似乎是一个可行的方法，当油膜粘性系数足够小时，短时间内油膜的厚度也会发生明显变化。但在实际的风洞实验中，当油膜粘性系数较小时，很容易受各种扰动的影响，导致干涉条纹扭曲变形，因而不能获得准确的油膜厚度。

具体到数据处理方面，油膜干涉摩擦应力测量技术主要由两部分构成，一是由光学干涉图像获得油膜的厚度，二是由油膜的厚度随时间的变化率获得摩擦应力。前者显然与温度的变化无关，其处理方法可参考任何一本有关光学干涉方面的书籍。后者在求解薄油膜方程时需要考虑油膜粘性系数随时间的变化。以一维情况为例，见公式(4)。

$\frac{\&PartialD;h}{\&PartialD;t}+\frac{\&PartialD;}{\&PartialD;x}\left(\frac{{\mathrm{\&tau;}}_x{h}^2}{2\mathrm{\&mu;}}\right)=0---\left(4\right)$

双时间步解法一种是将(4)式在t1到t2时间段积分，积分过程中将μ当作常数；一种是做定常摩擦应力假设(如(3)式)，分别测得t1时刻的条纹宽度为d1，t2时刻的条纹宽度为d2，条纹宽度对应的油膜厚度差为

$\mathrm{\&Delta;h}=\frac{\mathrm{\&lambda;}}{2}\left(\frac{1}{\sqrt{n_f^2-n_a^2{\sin }^2{\mathrm{\&theta;}}_i}}\right)---\left(5\right)$

其中λ为入射光波长，n_f、n_a分别为硅油和空气的折射率，θ_i为入射角，则摩擦应力为

$\mathrm{\&tau;}=\frac{0.5*(\mathrm{\&mu;}2+\mathrm{\&mu;}1)*(d2-d1)}{\mathrm{\&Delta;h}*(t2-t1)}---\left(6\right)$

上述两种方法要么根本没有考虑μ的变化，要么只是简单的取μ的平均值来代替粘性系数的变化。这样的处理带来的误差比较大，尤其是在高超声速实验中。以双时间步解法的第二种方法为例，t1时刻的不同、t1与t2之间间隔的不同会导致不同的结果。

发明内容

针对上述的摩擦应力测量方法存在没有考虑粘性系数变化或者对粘性系数变化采用简单平均的方式导致测量结果误差较大的技术问题，本发明提出了一种新的油膜干涉表面摩擦应力测量方法。

本发明的目的通过下述技术方案来实现：

一种油膜干涉表面摩擦应力测量方法，其具体包括以下的步骤：步骤一、在t1时刻和t2时刻之间，获取油膜粘性系数μ随时间的变化关系，并测量得到t1时刻和t2时刻的油膜厚度；步骤二、设定初始摩擦应力τ₀，并根据步骤一得到的油膜粘性系数μ随时间的变化关系，t1时刻的油膜厚度，数值计算薄油膜方程，并采用显式时间步推进到t2时刻，从而计算得到t2时刻的油膜厚度；步骤三、比较计算得到的t2时刻的油膜厚度与测量得到的t2时刻的油膜厚度，并根据比较结果对摩擦应力τ₀进行修正，重新从t1时刻开始显式时间步推进到t2时刻，当计算所得的油膜厚度与测量得到的油膜厚度一致时循环结束，此时的τ₀即为求解的摩擦应力τ。通过上述方法得到油膜干涉表面摩擦应力，这样的方法充分结合粘性系数的变化，同时又不需要连续测量油膜的温度，就能够得到准确的油膜干涉表面摩擦应力。

更进一步地，上述采用显式时间步推进到t2时刻，从而计算得到t2时刻的油膜厚度具体包括以下的步骤：将一维情况下的薄油膜公式进行离散，时间上采用向前差分，空间上采用中心差分，得到

其中n为时间步数，i为空间步数，μⁿ为第n个时间步时油膜的粘性系数，h表示油膜厚度，t表示时间，μ表示硅油的粘性系数，n＝(t2-t1)/Δt，x是物面坐标。

更进一步地，上述获取油膜粘性系数μ随时间的变化关系具体步骤为：首先多次测量油膜温度，获得此时间段内硅油油膜温度随时间的变化关系，并根据预先标定的硅油粘性系数和温度曲线获得油膜粘性系数随时间的变化关系。

更进一步地，上述步骤二中初始摩擦应力τ₀为：其中：μ1、μ2分别为t1和t2时刻油膜的粘性系数，d1、d2分别为t1和t2时刻油膜干涉条纹的间距(即相邻亮条纹中心线间距离)，Δh为一倍条纹间距对应的油膜厚度差。

更进一步地，上述步骤三中比较计算得到的t2时刻的油膜厚度与测量得到的t2时刻的油膜厚度具体为：当油膜干涉图像质量较高时，直接比较计算得到的t2时刻的油膜厚度值与测量得到的t2时刻的油膜度值。

更进一步地，上述步骤三中比较计算得到的t2时刻的油膜厚度与测量得到的t2时刻的油膜厚度具体为：当油膜干涉图像质量较差时，将计算得到的油膜厚度值进行线性拟合，计算得到油膜厚度斜率，并与测量得到的油膜厚度斜率进行比较，并判断计算与测量的油膜厚度斜率是否一致。

更进一步地，上述油膜粘性系数随时间的变化关系通过拟合得到，所述拟合为线性拟合或者多项式拟合。

通过采用以上的技术方案，本发明具有以下的有益效果：采用本发明的方法充分考虑了油膜的粘性系数随时间的变化情况，同时不需要连续测量油膜的温度并进行积分的复杂计算，就能够准确得到油膜干涉表面摩擦应力。

附图说明

图1为油膜干涉表面摩擦应力测量方法的流程图。

图2为测量系统示意图。

图3为t1时刻油膜中线灰度分布曲线。

图4为t2时刻油膜中线灰度分布曲线。

图5为t1油膜厚度分布。

图6为t2时刻油膜厚度分布

图7为计算与测量比较。

具体实施方式

下面结合说明书附图，详细说明本发明的具体实施方式。

本发明公开了一种油膜干涉表面摩擦应力测量方法，其具体包括以下的步骤：

步骤一、在t1时刻和t2时刻之间，获取油膜粘性系数μ随时间的变化关系，并测量得到t1时刻和t2时刻的油膜厚度。其中获取油膜粘性系数μ随时间的变化关系具体步骤为：首先多次测量油膜温度，获得此时间段内硅油油膜温度随时间的变化关系，并根据预先标定的硅油粘性系数和温度曲线获得油膜粘性系数随时间的变化关系。以某测量点为例，t1时刻到t2时刻测得的油膜温度如下表所示。

t'/s	0	0.4	1.2	2	2.8	3.6	4
								T/℃	69.80	70.03	70.45	70.87	71.28	71.74	71.90
t'/s	4.4	5.2	6	6.8	7.6	8
								T/℃	72.16	72.61	73.07	73.50	73.98	74.20

其中t'＝t-t1。根据实验得到的硅油的粘性系数和温度曲线标定结果，以及上表中检测得到的油膜温度随时间的变化表，拟合得到粘性系数与时间关系为：

μ＝-0.0007*t'+0.0894   (7)

其中拟合可以采用线性拟合，也可以是多项式拟合，多项式拟合时一般取2次多项式就够了。油膜厚度的测量与温度的变化无关，其处理方法可参考任何一本有关光学干涉方面的书籍就可以实现，其不属于本发明的重点，在此不相袭赘述其原理和实现方法。

步骤二、设定初始摩擦应力τ₀，并根据步骤一得到的油膜粘性系数μ随时间的变化关系，t1时刻的油膜厚度，数值计算薄油膜方程，并采用显式时间步推进到t2时刻，从而计算得到t2时刻的油膜厚度。采用显式时间步推进到t2时刻，从而计算得到t2时刻的油膜厚度具体过程为：

将(4)式进行离散，时间上采用向前差分，空间上采用中心差分，可得

$\frac{h_i^{n+1}-h_i^n}{\mathrm{\&Delta;t}}+\frac{\mathrm{\&tau;}}{{2\mathrm{\&mu;}}^n}*\left(\frac{{\left(h_{i+1}^n\right)}^2-{\left(h_{i-1}^n\right)}^2}{2\mathrm{\&Delta;x}}\right)=0---\left(8\right)$

其中n为时间步数，i为空间步数，μⁿ为第n个时间步时油膜的粘性系数，由(7)式给出

μⁿ＝-0.0007*n*Δt+0.0894   (9)

将(8)式整理得

$h_i^{n+1}=h_i^n-\frac{\mathrm{\&tau;}*\mathrm{\&Delta;t}}{2{\mathrm{\&mu;}}^n}*\left(\frac{{\left(h_{i+1}^n\right)}^2-{\left(h_{i-1}^n\right)}^2}{2\mathrm{\&Delta;x}}\right)---\left(10\right)$

其中初始摩擦应力τ₀可以根据(6)式预估，将初始摩擦应力τ₀和(9)式代入(10)式，取Δt＝0.005s，Δx＝1PIXEL，显式推进求解油膜厚度h，当n＝(t2-t1)/Δt时即得t2时刻的计算油膜厚度。(4)式表示一维情况，是(1)式二维情况的特例，为了简便起见用了一维情况。x，z是物面坐标，y是法向坐标，因为摩擦应力是平行于物面的，所以没有y下标。(1)式是薄油膜方程，是在一定假设下理论推导出来的。

步骤三、比较计算得到的t2时刻的油膜厚度与测量得到的t2时刻的油膜厚度，并根据比较结果对摩擦应力τ₀进行修正，重新从t1时刻开始显式时间步推进到t2时刻，当计算所得的油膜厚度与测量得到的油膜厚度一致时循环结束，此时的τ₀即为求解的摩擦应力τ。在t2时刻，可以直接比较计算得到的t2时刻的油膜厚度值与测量得到的t2时刻的油膜厚度值，也可以比较油膜厚度的斜率。因为在油膜厚度的计算中有时因为图像质量不是很好，每个像素点处的油膜厚度的相对误差较大，这时如果用厚度判据可能导致摩擦应力被过分修正，远远偏离正确结果。用斜率作为判据是因为在一小块的测量区域，如果物面的曲率不是很大时摩擦应力在空间的变化很小，表现的结果就是油膜为楔形，在t时刻油膜的厚度总体呈现出随x线性变化，所以对油膜厚度进行线性拟合，然后比较斜率。一般在油膜干涉图像质量较高时优先厚度判据。

油滴法测量摩擦应力时因为目标区域较小，温度在此区域的空间变化量较小，而CFD结果表明温度随时间在小范围(<10℃)的变化对摩擦应力影响很小。当目标区域较为平缓时(如平板、机翼上下表面等)，摩擦应力的空间梯度较小。通过对薄油膜方程的量纲分析可以发现，当摩擦应力的空间梯度较小时可以作定常摩擦应力处理。综合可知，获得τ精确解的最有效的方法是考虑μ随时间变化的影响。

上述的实施例中所给出的系数和参数，是提供给本领域的技术人员来实现或使用本发明的，本发明并不限定仅取前述公开的数值，在不脱离本发明的发明思想的情况下，本领域的技术人员可以对上述实施例作出种种修改或调整，因而本发明的保护范围并不被上述实施例所限，而应该是符合权利要求书提到的创新性特征的最大范围。

Claims (7)

1.一种油膜干涉表面摩擦应力测量方法，其具体包括以下的步骤：步骤一、在t1时刻和t2时刻之间，获取油膜粘性系数                                               随时间的变化关系，并测量得到t1时刻和t2时刻的油膜厚度；步骤二、设定初始摩擦应力，并根据步骤一得到的油膜粘性系数随时间的变化关系，t1时刻的油膜厚度，数值计算薄油膜方程，并采用显式时间步推进到t2时刻，从而计算得到t2时刻的油膜厚度；步骤三、比较计算得到的t2时刻的油膜厚度与测量得到的t2时刻的油膜厚度，并根据比较结果对摩擦应力进行修正，重新从t1时刻开始显式时间步推进到t2时刻，当计算所得的油膜厚度与测量得到的油膜厚度一致时循环结束，此时的即为求解的摩擦应力。

2.如权利要求1所述的种油膜干涉表面摩擦应力测量方法，其特征在于所述采用显式时间步推进到t2时刻，从而计算得到t2时刻的油膜厚度具体包括以下的步骤：将一维情况下的薄油膜公式进行离散，时间上采用向前差分，空间上采用中心差分，得到

，其中n为时间步数，i为空间步数，为第n个时间步时油膜的粘性系数，h表示油膜厚度，t表示时间，μ表示硅油的粘性系数，n=（t2-t1）/，x是物面坐标。

3.如权利要求1或者2所述的种油膜干涉表面摩擦应力测量方法，其特征在于所述获取油膜粘性系数随时间的变化关系具体步骤为：首先多次测量油膜温度，获得此时间段内硅油油膜温度随时间的变化关系，并根据预先标定的硅油粘性系数和温度曲线获得油膜粘性系数随时间的变化关系。

4.如权利要求1或者2所述的种油膜干涉表面摩擦应力测量方法，其特征在于所述步骤二中初始摩擦应力为：，其中：、分别为t1和t2时刻油膜的粘性系数，、分别为t1和t2时刻油膜干涉条纹的间距，为一倍条纹间距对应的油膜厚度差。

5.如权利要求1或者2所述的种油膜干涉表面摩擦应力测量方法，其特征在于所述步骤三中比较计算得到的t2时刻的油膜厚度与测量得到的t2时刻的油膜厚度具体为：当油膜干涉图像质量较高时，直接比较计算得到的t2时刻的油膜厚度值与测量得到的t2时刻的油膜度值。

6.如权利要求1或者2所述的种油膜干涉表面摩擦应力测量方法，其特征在于所述步骤三中比较计算得到的t2时刻的油膜厚度与测量得到的t2时刻的油膜厚度具体为：当油膜干涉图像质量较差时，将计算得到的油膜厚度值进行线性拟合，计算得到油膜厚度斜率，并与测量得到的油膜厚度斜率进行比较，并判断计算与测量的油膜厚度斜率是否一致。

7.如权利要求1或者2所述的种油膜干涉表面摩擦应力测量方法，其特征在于所述油膜粘性系数随时间的变化关系通过拟合得到，所述拟合为线性拟合或者多项式拟合。