空气动力学实验报告.docx

实验一 边界层流动测量实验摘要：边界层，又称为流动边界、附面层，它是流体流动过程中，紧贴壁面的粘性阻力不可忽略的一层薄薄的流体，它对主要流体运动的影响很大。自普朗特提出该概念起，边界层研究就一直是流体力学研究中一个焦点和难点课题。本实验通过热线风速仪测量距离凹口平板前缘不同位置点流体的速度分布情况，并对实验数据加以分析处理，从而确定出在不同工况中的边界层的厚度、位移厚度，以及避免粘性力等参数，最终分析边界层的特性。关键词：边界层，热线风速仪，粘性力，雷诺数，拟合，标定1. 实验简介此次实验是在一个开口式风洞中进行的，该风洞试验段截面尺寸为：500mm*500mm。设置风洞风机的运行频率为20Hz和30Hz、，利用热线风速仪测量凹槽分离点20mm的边界层上的速度分布。然后用两种不同的方法拟合热线风速仪实验前后标定曲线，得出标定误差值，从而分析比较这两种拟合方法的优缺点，并分析出实验中热线性能的稳定性。2. 实验步骤1) 将皮托管固定在风洞试验段，轴线和来流速度方向平行。记录皮托管标定系数k。皮托管静压连接到压力传感器负压接口，皮托管总压连接到压力传感器通道1；2) 热线风速仪探头安装在二位坐标架上，连接热线探头与恒温控制器输入、输出。此时热线恒温控制器切勿通电！将热线探头移至和皮托管同一高度；3) 热线输出连接到数据采集卡AI0，皮托管输出连接到数据采集卡AI1；4) 将热线恒温控制器通电，打开MATLAB热线风速仪标定程序“hw calibration.m”，改变文件名运行程序；5) 将热线移动至测量点（距离凹腔分离点X=20mm）上方自由来流中，调整风洞风速，风机运行频率f=30Hz, MATLAB运行热线速度分布测量程序“hw measurement.m”改变文件存储名称。改变风洞风速，风机运行频率f=20Hz,重复步骤4；6) 打开MATLAB热线风速仪标定程序hw calibration.m，改变标定参数存储文件名，重新运行标定程序。7) 按要求处理实验数据。3. 实验设备1) 开口式风洞，实验段截面500mm*500mm。实验段风洞壁面装有凹腔模型，测量模型分离点上游边界层速度分布。2) 程控二位坐标架：计算机软件可精确控制探针在边界层内移动。3) 压力传感器：将压力信号转化为电信号，灵敏度为0.0245 V/mm.4) 皮托管：测量来流静压Ps和总压Pt，用于测量来流速度U和标定热线风速仪。5) 热线风速仪：测量气体流场（尤其是湍流场）中速度大小，具有较高的频率响应（10kHz）。6) 数据采集卡：NI-USB 6212 BNC（16位，8差分通道）；7) 计算机：控制二维坐标架，控制和执行数据采集；其中，热线风速仪的工作原理如下：热线风速仪前段的热线探头作为后端控制电路中电桥的一个桥端电阻。电桥工作时把探头上的电阻丝加热，当流体以速度流过探头时带走一部分热量Q电桥设计是为了保持探头热丝温度恒定，探头热丝损失热量Q后，两端电压发生变化从而提供相应的热量以弥补损失掉的那部分热量，最终保证探头热丝的温度恒定。通过标准速度源可以将速度U与电桥电压E之间的关系拟合出来。这样在测量时，只要测量到电桥电压E，就可以知道当地风速U。热线风速仪的标定：用风洞自由来流速度作为标准速度，自由来流速度U由皮托管测量得到。拟合关系式用：1）Kings Law：；2）。4. 实验数据处理1) 分别用两种不同方法拟合热线风速仪实验前后标定曲线，给出标定误差值，分析拟合方法优缺点。将实验前后热线风速仪标定曲线画在同一张图上，将实验后热线标定加上温度矫正和前面两组曲线画在同一张图上，分析实验中热线性能是否稳定及温度矫正的作用。表4-1 实验前标定时采集的数据流速/U电压/E 温度/T3.08361.696221.2363.99461.742521.2944.94061.782421.3925.86891.817121.5216.7911.847521.6547.7121.876121.7678.62041.901521.8729.55781.925621.98910.4681.947322.173表4-2 实验后标定时采集的数据流速/U电压/E温度/T3.05831.692422.0293.99231.738822.0094.93641.778922.0335.84981.813322.0916.78111.844422.1167.69911.872622.1918.61191.89822.2829.54051.921622.402按照Kings Law：关系拟合程序：clc;clear;data=load(E:2014标定实验hw calibration-xxxxxx(U_E)_0527.txt);%载入数据v=data (: 1);w=v. 0.45;e=data (:, 2);f=inline (a(1)+a(2)*w,a,w);b,c=lsqcurvefit(f,1.35,0.205,w,e);%b是(A,B)向量组合，c是目标函数值；b,cfor i=1:x v(i)=data(i,1) h(i)=b*1;v(i).0.45;endfigure;plot(v,h,-r)Xlabel(流体流速U);Ylabel(电桥电压E);title(电压E与流速U的关系);按照关系拟合程序：clc;clear;data=load(E:2014标定实验hw calibration-xxxxxx(U_E)_0527.txt);%载入数据v=data(:,1);e=data(:,2);f=inline(a(1)+a(2)*e+a(3)*e.2+a(4)*e.3+a(5)*e.4,a,e);b,c=lsqcurvefit(f,1 1 1 1 1,e,v);%b是(A,B)向量组合，c是目标函数值；b,cfor i=1:x e(i)=data(i,2); h(i)=b*1 e(i) e(i).2 e(i).3 e(i).4;endfigure;plot(e,h,-.g)Xlabel(电桥电压E);Ylabel(流体流速U);title(电压E与流速U的关系);拟合曲线为：图1 Kings Law法拟合曲线由图1可以看出，未拟合前实验前后的曲线间距明显，这就表明，在实验结束后热线探头发生了一定距离的偏移。而且使用Kings Law拟合后，拟合前后的曲线还是存在一定的误差。根据，实验前曲线的标定误差为2.618%，实验后曲线的标定误差为2.433%。图2 拟合方法二标定曲线使用同样的误差公式计算得到：对实验前的曲线的标定误差为：0.598%，对实验后的曲线的标定误差为：0.803%。由上图可以看出，使用拟合方法二拟合后的曲线可以与未拟合的曲线比较好地重合，说明拟合误差很小。通过比较上面两图，可以得出，第二种拟合方法明显比第一种方法得到的拟合效果好，这就说，第二种拟合方法的标定误差较小。同时，从实验前后的两条曲线可以看出，实验前后热线探头发生了一定距离的偏移，这证明该实验使用的热线性能还不够稳定。2) 利用文献方法，预测测量点距离壁面精确距离，测量点摩擦速度，计算壁面粘性应力，边界层厚度，边界层位移厚度，边界层动量厚度，基于动量厚度的边界层雷诺数等参数。根据文献指导并参考徐进同学的程序，然后进行必要的校正完善，得到以下程序及图形：clc;clear;file=E:2014边界层试验x20f20 profile.txt;data=dlmread(file);junk=sortrows(data);yi=junk(:,1)./1000;u=junk(:,2);umean=mean(data(1:20,4);u99=umean*0.99;%自由来流风速；yd=-9.9:0.01:-4;for j=1:length(yd) %截止距离循环迭代； y=yi+yd(j)./1000; ut=0.20:0.01:1; lp=1.5e-5./ut;k=0.41;B=5.0; for i=1:length(ut) %壁面摩擦应力循环迭代； yp=y./lp(i);up=u./ut(i); ye=up+exp(-k*B)*(exp(k*up)-1-k*up-(k*up).2/2-(k*up).3/6); err(i,j)=sum(abs(yp(1:35)-ye(1:35)./yp(1:35)/35;%残余函数计算； endendm indexi=min(err);%行数；n indexj=min(m);%列数；yoff=yd(indexj)%截止距离；utf=ut(indexi(indexj)%壁面摩擦应力；lpf=1.5e-5./utf;yp=(yi+yoff./1000)./lpf;upp=1:0.5:25;up=junk(:,2)./utf;ue=(1./k).*log(yp)+B;ye=upp+exp(-k*B)*(exp(k*upp)-1-k*upp-(k*upp).2/2-(k*upp).3/6);figure(1);semilogx(yp,up,-ob,linewidth,1);hold on;semilogx(ye,upp,-or,linewidth,1);hold on;semilogx(yp,ue,-k,linewidth,1);xlabel(y+);ylabel(U+);legend(Experimental data,Spalding profile,Log-law);axis(1 1000 0 26); yabs=yi+yoff./1000;ul=1-u/umean;yabsp=trapz(yabs,ul).*1000%边界层位移厚度；uul=(u/umean).*(1-u./umean);theta=trapz(yabs,uul).*1000%边界层动量厚度；Re=1.2*umean*theta/1000./(1.5e-5)%基于动量厚度的边界层雷诺数；figure(2);Yreal=yi.*1000+yabsp;Ureal=junk(:,2);plot(Ureal,Yreal,-oc)Xlabel(Ureal),Ylabel(Yreal);axis(0 8 0 120);取f=20Hz时，其风速近似为6.26m/s，经循环迭代得到，测点距壁面的误差距离为9.34mm；测量点摩擦速度为0.34m/s；计算出壁面的剪切应力为0.149N/m2；边界层厚度从上图可以近似得到，为12mm；边界层位移厚度为2.3755mm；边界层动量厚度为2.0819mm；基于动量厚度的雷诺数为1051；取f=30Hz时，其风速近似为9.3m/s，经循环迭代得到，测点距壁面的误差距离为7.99mm；测量点摩擦速度为0.46m/s；计算出壁面的剪切应力为0.273N/m2；边界层厚度从上图可以近似得到，为11mm；边界层位移厚度为2.0217mm；边界层动量厚度为1.5098mm；基于动量厚度的雷诺数为1122。3) 画出边界层速度分布曲线（基于边界层厚度和自由来流无量纲化），边界层速度rms分布曲线。画出边界层速度分布曲线（基于边界层内层参数无量纲化）。（图例为20Hz）clc;clear;file=E:2014边界层试验x20f20 profile.txt;data=dlmread(file);junk=sortrows(data);yi=junk(:,1);ui=junk(:,2);y=yi./12;u=ui./6.26;plot(u,y,-ob);Xlabel(基于自由来流风速无量纲化),Ylabel(基于边界层厚度无量纲化);axis(0.45 1.1 0 10)clc;clear;file=E:2014边界层试验x20f20 profile.txt;data=dlmread(file);junk=sortrows(data);yi=junk(:,1);rms=junk(:,3);y=yi-9.34;plot(rms,y,-ob);Xlabel(rms),Ylabel(边界层);axis(0 0.8 0 100)clc;clear;file=E:2014边界层试验x20f20 profile.txt;data=dlmread(file);junk=sortrows(data);yi=junk(:,1);ui=junk(:,2);y=(yi-9.34)./12;u=ui./0.34;plot(u,y,-ob);Xlabel(基于边界层摩擦速度无量纲化),Ylabel(基于边界层厚度无量纲化);axis(8 20 0 10)以上图形均为在20Hz下的数据得到的，30Hz的相应图形用同样的方法即可得到，在此不再赘述。4) 对比分析讨论不同风洞运行速度下边界层速度分布曲线特性，边界层基本参数。项目20Hz30Hz自由来流速度m/s6.269.3测量点摩擦速度m/s0.340.46边界层厚度/mm1211边界层位移厚度/mm2.37552.0217边界层动量厚度/mm2.08191.5098边界层雷诺数10511122通过mat lab编程将两种频率下的三类图形放在一起进行比较，由数据文件获知，在风洞运行频率f=20Hz时，自由来流速度约为6.26m/s，而在风洞运行频率f=30Hz时，自由来流速度为9.3m/s，从图上可以看到，无量纲化后的图形几乎一致，边界层速度分布曲线其基本走势均为距离壁面较近的区域速度逐渐增加，到一定程度后就不再增加，这段距离即为边界层厚度。5) 对比分析讨论不同风洞运行速度下剪切层速度大小，rms分布曲线特性。clc;clear;file1=E:2014边界层试验x20f20 profile.txt;file2=E:2014边界层试验x20f30 profile.txt;data1=dlmread(file1);data2=dlmread(file2);junk1=sortrows(data1);junk2=sortrows(data2);yi1=junk1(:,1);ui1=junk1(:,2);yi2=junk2(:,1);ui2=junk2(:,2);y1=yi1-9.34;u1=ui1;y2=yi2-7.99;u2=ui2;figure(1);plot(u1,y1,-og,u2,y2,-*c);legend(20Hz,30Hz);Xlabel(V),Ylabel(Y);title(不同频率下剪切层速度大小比较);hold onyi1=junk1(:,1);rms1=junk1(:,3);yi2=junk2(:,1);rms2=junk2(:,3);y1=yi1-9.34;y2=yi2-7.99;figure(2);plot(rms1,y1,-og,rms2,y2,-*c);legend(20Hz,30Hz);Xlabel(rms),Ylabel(Y);title(不同频率下rms大小比较);从上面剪切层速度大小分布比较图和rms分布比较图两图可以看到，距壁面相同距离时，20Hz下的剪切层速度小于30Hz下的剪切层速度，从rms图上可以看到，基本趋势两种风速下是一致的，但是距壁面相同距离时，30Hz下的rms值要大于20Hz时的值，而且30Hz时的rms最大值为1.18，而20Hz时的最大值为0.704。6) 在Urms最大值处，速度的能量谱图，分析能量谱特征。从图上可以看到峰值出现的地方，对应的频率即为该信号主要能量来源，从低到高频率阶次依次为100Hz,160Hz,220Hz,280Hz，而且由图形趋势看到其主要能量集中于低频处，频率越高处能量越低。实验二 数据采集与处理1. 实验仪器1) 信号发生器及示波器：DSOX2012A；2) 数据采集卡：NI-USB 6212 BNC（16位，8差分通道，400 kS/s）；3) 计算机：控制和执行数据采集；2. 实验步骤1) 连接信号发生器，数据采集卡，电脑等设备；2) 设置信号发生器产生正弦波信号,信号频率为f=200Hz,信号幅值为A=0.05v；3) 信号发生器产生信号连接到数据采集卡AI0；4) 打开MATLAB数据采集程序data acquisition.m， 设置采样频率Fs=100Hz， 通道输入范围（InputRange）为-0.2 0.2。更改数据存储文件名，运行程序；5) 检查步骤4中数据是否成功保持。重复步骤4改变采样频率Fs=200Hz。更改数据存储文件名，运行程序；6) 检查步骤5中数据是否成功保持。重复步骤4改变采样频率Fs=400Hz。更改数据存储文件名，运行程序；7) 检查步骤6中数据是否成功保持。重复步骤4改变采样频率Fs=1000Hz。更改数据存储文件名，运行程序；8) 检查步骤7中数据是否成功保持。重复步骤7，改变通道输入范围（InputRange）为-10 10。更改数据存储文件名，运行程序。3. 实验数据处理1) 画出数据随时间变化曲线（带点线）x(t);不同参数下所得到的曲线画到同一张图上，用不同颜色标示。讨论区别，及采样参数影响。clc;clear;fip1=fopen(E:2014数据采集试验g1y100r02.bin,rt);fip2=fopen(E:2014数据采集试验g1y200r02.bin,rt);fip3=fopen(E:2014数据采集试验g1y400r02.bin,rt);fip4=fopen(E:2014数据采集试验g1y1000r02.bin,rt);fip5=fopen(E:2014数据采集试验g1y3000r02.bin,rt);fip6=fopen(E:2014数据采集试验g1y3000r10.bin,rt);A1,NUM1=fread(fip1,inf,double)A2,NUM2=fread(fip2,inf,double)A3,NUM3=fread(fip3,inf,double)A4,NUM4=fread(fip4,inf,double)A5,NUM5=fread(fip5,inf,double)A6,NUM6=fread(fip6,inf,double)t1=(0:1:9)*10-2;figure(1);plot(t1,A1(1:10),- Or);%第一hold ont2=(0:1:19)/2*10-2;plot(t2,A2(1:20),-+m);hold ont3=(0:1:39)/4*10-2;plot(t3,A3(1:40),-*g);hold ont4=(0:1:99)*10-3;plot(t4,A4(1:100),-c);hold ont5=(0:1:299)/3*10-3;plot(t5,A5(1:300),-b);Xlabel(时间t/s);Ylabel(幅值V/v);title(不同采样频率下的曲线比较)figure(2);t6=(0:1:299)/3*10-3;plot(t6,A5(1:300),-400Hz），当把频率取到1000Hz时，采样信号有了原始信号的轮廓，但还是失真比较大，当把频率取为3000Hz时，采样信号已经很大程度的逼近了原始信号，然后将采样幅值范围取为-10，10，在峰值处有了平台信号失真，所以幅值取为-0.2,0.2非常接近真实信号。因此在采集数据时，必须满足采样定理，同时信号的灵敏度要求幅值与真实信号幅值接近。2) 计算采集到数据的平均值xm，均方差xrms。比较不同采样参数下得到值的区别。讨论区别，及采样参数影响。clc;clear;fip1=fopen(E:2014数据采集试验g1y100r02.bin,rt);A1,NUM1=fread(fip1,inf,double);sum1=0;for i=1:NUM1 sum1=A1(i)+sum1;endx1=sum1/NUM1sum2=0;for i=1:NUM1 cha=A1(i)-x1; fang=cha.2; sum2=sum2+fang;endxrms1=sqrt(sum2/NUM1)依次算出6种情况下的均值和均方差（即标准差）。表2-1 不同情况下采集数据的均值和均方差列表比较项目100Hz0.2200Hz0.2400Hz0.21000Hz0.23000Hz0.23000Hz10平均值0.0110-0.0168-3.7972e-4-3.5217e-4-3.5032e-4-0.0010均方差0.00209.5665e-40.00320.01770.01770.0177编程计算求得各自的均值和均方差。从上表可以看出，100Hz和200Hz时均方差较小，而均值较大，这与实际情况一致，因为其图形近乎为一直线；400Hz时严重失真，图形呈锯齿波形，均值较小近乎为0；而后面三种高频采集的信号均值都很小，而且均方差值接近，较好地反映了原始信号，相比之下低的采样频率采集的信号均方差的误差较大。3) 用mat lab编程计算所得数据的傅里叶变化X（k），验证Parsevals Theorem。 画出|X(k)|随频率变化曲线。不同采样参数得到的曲线画到同一张图上，用不同颜色标示。讨论区别，及采样参数影响。clc;clear;fip1=fopen(E:2014数据采集试验g1y100r0.2.bin,rt);fip2=fopen(E:2014数据采集试验g1y200r0.2.bin,rt);fip3=fopen(E:2014数据采集试验g1y400r0.2.bin,rt);fip4=fopen(E:2014数据采集试验g1y1000r0.2.bin,rt);fip5=fopen(E:2014数据采集试验g1y10000r0.2.bin,rt);fip6=fopen(E:2014数据采集试验g1y10000r10.bin,rt);data1=fread(fip1,512,double)data2=fread(fip2,512,double)data3=fread(fip3,512,double)data4=fread(fip4,512,double)data5=fread(fip5,512,double)data6=fread(fip6,512,double)N=512; %512个采样点figure(1);Fs=100; %采样频率t=(0:N-1)/Fs; %时域Y=fft(data1,N); %快速傅里叶变换f=(0:N-1)*Fs/N; %频域的点数plot(f,abs(Y),-sr) %频谱图hold onFs=200; t=(0:N-1)/Fs;Y=fft(data2,N);f=(0:N-1)*Fs/N;plot(f,abs(Y),-dm)hold onFs=400; t=(0:N-1)/Fs;Y=fft(data3,N); f=(0:N-1)*Fs/N; plot(f,abs(Y),-hg)hold onFs=1000; t=(0:N-1)/Fs; Y=fft(data4,N); f=(0:N-1)*Fs/N;plot(f,abs(Y),-*c) hold onFs=10000;t=(0:N-1)/Fs;Y=fft(data5,N);f=(0:N-1)*Fs/N;plot(f,abs(Y),-oy)hold onFs=10000;t=(0:N-1)/Fs; Y=fft(data6,N); f=(0:N-1)*Fs/N; plot(f,abs(Y),-pk)800Hz200Hz从上图6种采样信号单独的傅立叶变换后的图形，可以很明显地看到，100Hz和200Hz两种情况下，只有0Hz处有幅值，属于直流分量，而其余部分近乎为零，没有能量，没有采集到信号，或者说采集到的信号完全失真。一般来说，满足采样定理的信号频谱图是对称的，后面的几种情况就是这样，明显看到在200Hz处出现峰值，表示能量集中在此频率下，说明原始信号的频率为200Hz，这与实际情况是相符的。在后两种情况下，采样频率是相同的（10000Hz），但是采样范围是不同的(-0.2,0.2,-10,10)，而从图上看到其傅立叶变换的频谱图是完全相同的，表明频谱图与采样频率密切相关的。4) 画出能量谱图，不同采样参数得到的曲线画到同一张图上，用不同颜色标示。讨论区别，及采样参数影响。clc;clear;fip1=fopen(E:2014数据采集试验g1y100r0.2.bin,rt);fip2=fopen(E:2014数据采集试验g1y200r0.2.bin,rt);fip3=fopen(E:2014数据采集试验g1y400r0.2.bin,rt);fip4=fopen(E:2014数据采集试验g1y1000r0.2.bin,rt);fip5=fopen(E:2014数据采集试验g1y10000r0.2.bin,rt);fip6=fopen(E:2014数据采集试验g1y10000r10.bin,rt);data1=fread(fip1,256,double)data2=fread(fip2,256,double)data3=fread(fip3,256,double)data4=fread(fip4,256,double)data5=fread(fip5,256,double)data6=fread(fip6,256,double)N=256; %256个采样点figure(1);Fs=100; %采样频率t=(0:N-1)/Fs; %时域Y=fft(data1,N); %快速傅里叶变换Z=Y.*conj(Y);f=(0:N-1)*Fs/N; %频域的点数plot(f,Z,-sr) %频谱图hold onFs=200; t=(0:N-1)/Fs;Y=fft(data2,N); Z=Y.*conj(Y);f=(0:N-1)*Fs/N;plot(f,Z,-dm) hold onFs=400;t=(0:N-1)/Fs; Y=fft(data3,N);Z=Y.*conj(Y);f=(0:N-1)*Fs/N; plot(f,Z,-hg) hold onFs=1000;t=(0:N-1)/Fs; Y=fft(data4,N); Z=Y.*conj(Y);f=(0:N-1)*Fs/N; plot(f,Z,-*c)hold onFs=10000;t=(0:N-1)/Fs; Y=fft(data5,N); Z=Y.*conj(Y);f=(0:N-1)*Fs/N; plot(f,Z,-oy) hold onFs=10000;t=(0:N-1)/Fs;Y=fft(data6,N);Z=Y.*conj(Y);f=(0:N-1)*Fs/N; plot(f,Z,-pk)从图上可以看出能量集中于200Hz处，表明原始信号为200Hz，单独分析每一种情况的能量谱，发现其表现跟频谱图是一致的，100Hz和200Hz失真，没有能量，以后的几种情况能量均集中于200Hz处，这与真实信号的频率是一致的。