《超音速翼型气动力特性研究汇总13399.pdf》由会员分享，可在线阅读，更多相关《超音速翼型气动力特性研究汇总13399.pdf（16页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、 超音速翼型气动力特性研究 摘要:本文研究方程为0.3(1)zxx 的轴对称超音速翼形在马赫数为 2，攻角分别为 0，2情形下的气动力特性，基于对翼型进行离散化处理得到该翼型的物理参数及气动力的近似解，并逐步减小空间步长x来提高解的精度。在步长数分别为5、20、50 及攻角为 0、2的条件下，计算求得翼型头部斜激波后的流动参数，并由此求解各分区相应参数，列出：表面压力 Cp 分布曲线 Cpx，及表面密度、温度分布曲线/-x、T/T-x。在不同条件下得出的轴向力 Ca、法向力 Cn、升力 Cl、阻力Cd 及绕头部顶点俯仰力矩 Cm 的表格。最终分析了编程计算的准确性与精度，分析了压力系数、温度、
2、密度沿该翼型的分布特性，并分析了不同攻角对该翼型气动特性的影响。问题描述 已知方程为0.3(1)zxx 的薄翼形，求该翼型在来流马赫数为 2，攻角分别为 0，2情形下的受力情况。对 x 范围（0，1）内分别按 5 等份、20 等份和50 等份进行离散计算，得到表面压力 Cp 分布曲线 Cpx，表面密度、温度分别曲线/、T/T。计算得出出轴向力 Ca、法向力 Cn、绕头部顶点俯仰力矩 Cm 及升力 Cl、阻力 Cd。计算方案:（一）计算思路：超音速来流以一定攻角遇到类似于楔形体的机翼前缘，在上下面都有可能产生附体斜激波，要是攻角过大也有可能不产生附体斜激波，这里首先需要根据斜激波的 关系曲线图来
3、作出判断。经判断，如果顶点处产生斜激波，即使用斜激波前后的马赫数、密度、温度、压强计算公式计算出顶点斜激波后的各项物理参数。接着，根据翼型的形状可知，气流在通过膨胀波之后会经过一系列的向外的转折角，根据普朗特-迈耶膨胀波理论，超音速气流经过每一个折角都会产生膨胀波。根据数值计算的基本原理，计算机不能处理连续曲线上随 x 值变化而连续变化的折角，所以在计算之前必须对翼型的几何结构进行离散化处理。离散化之后即可根据膨胀波前后马赫数的关系公式计算出每一个折角处膨胀波后的马赫数，然后根据膨胀波前后密度、温度、压强的计算公式计算出每一个膨胀波后的密度、温度、压强。得到以上基本的物理参数之后，即可用压强
4、P 的分布计算出压力系数 Cp 的分布进而计算出翼型所受的轴向力、法向力、升力、阻力及力矩系数。在进行以上计算之前的首要工作是编制 PM 表、等熵流动角度与马赫数的关系表。在具体操作中可使用已知的显示函数式进行编制，无需手动输入。（二）基本假设：实际上的翼型气动力会受到很多因素的干扰，用激波-膨胀波法计算时对实际的物理模型做了一些简化，假设：1、离散化之后任意两个离散点之间的物理参数是均匀分布的。2、不考虑斜激波、膨胀波的相交与反射。3、翼型端部的斜激波不会在上部削弱，即斜激波不会滑移。（三）坐标系建立 建立坐标系时，以机翼的前段点为原点O，弦线为 x 轴，垂直于弦线的直线为 y轴，x 轴正向
5、指向机翼尾缘，y 轴正向垂直于x 轴向上，如图所示：在进行力学参数计算时，也是以 X 轴正向为正，z 轴正向为正。（四）物面离散 物面离散的步骤如下（以分成 5 段为例）：首先在 X（01）上以 0.2 为步长取出 0、0.2、0.4、0.6、0.8、1 六个点，再根据0.3(1)zxx 计算出 x=0、0.2、0.4、0.6、0.8、1 时所对应的 z 值，对上述的 12 个点进行画图，即可得到以 5 段离散得到的翼型图。分 20 段、50 段进行离散的步骤与分 5 段的步骤相同。各个离散物面与水平面的夹角也可以求出，以方便后来的气流偏转角的计算，计算公式如下：180z(i+1)-z(i)(
6、i)=arctan()x(i+1)-x(i)（五）各区域物理参数的计算 （1）第一区流动参数的计算 按照分段数的要求将物面离散后，首先需要计算出各个区段的物理参数，特别是马赫数，因为斜激波、膨胀波前后的马赫数关系是求解其他物理参数的基础。由于存在攻角的缘故，超音速气流流经上下翼面的第一个偏转角并不等于半顶角1。超音速气流遇到机翼前缘时的偏转角 w 计算公式如下：对于上翼面：1=1-w 对于下翼面：1=1+w 然后，根据斜激波的-关系表来作出判断超音速来流遇到机翼前缘是否会产生附体斜激波。程序会查找来流马赫数M1 为 2 时的-关系表找出此对应的，如果没有找到与之对应的，即说明该情况下，在机翼前
7、缘不会产生斜激波，程序会报错。如果产生了斜激波，由已知的 M1 和计算出斜激波后的马赫数、温度、密度、压强。计算公式如下：211221222222112cos1()21sin1sin112MMMMM 22212211sin11sin11MM 2221122221122221121sin112sin11sin21sinpMpMMTTM （2）第 i 区气流参数的计算 在一区之后的每区前端点均产生膨胀波，第 i 区(1i)相对音速来流的折转角为 22i111arctan(1)arctan(1)()11iiwMM 通过对 P-M 表进行插值，根据iw查等熵流动函数表可以求出第i 区对应的马赫数iM。
8、由于膨胀波为等熵波，波前后流动参数满足等熵关系式：121 1211 121 012!112112112112iiiiiiMppMMM 2112112112iiiiMTTM 然后，即根据 P 的分布推算 Cp 的分布，计算公式如下：22222111111112.822222iiiiiiPPPPPPPPPPPPCpiVVVMVRTPa对于本题（六）离散元力/力矩计算公式、翼型合力/合力矩计算公式 离散化后，根据最初的假设每区流动参数均匀分布，压力作用在分区中点。轴向力为 x 方向上压力的分力，法向力为 y 方向的分力，总的轴向力与法向力计算公式如下 2211sintancos112.81 11 1
9、22NNxiiSPiFPi xCaVV 2211coscos112.81 11 122NNyiiSPiFPi xCnVV 升力、阻力 与轴向力、法向力间的差距仅仅由攻角造成，可以得到以下的计算公式：cossinsincosClCnCaCdCnCa 每一段上的压力对于机翼前缘端点都会产生一个俯仰力矩，取逆时针（抬头）为正方向，计算公式如下：上翼面：11tan()iiiiiiMP xyP xx 下翼面：22tan()iiiiiiMP xyP xx 总的俯仰力矩系数计算公式如下：1211 12NiiMCmV 程序流程图 结束 输入分段数 N,攻角 alpha 编制PM表、等熵流动函数表 离散化，计算
10、上下翼面各分段处的偏转角 1，2 计算考虑攻角后的气流遇到翼型前缘端点后流经上下翼面分别对应的偏转角 w 计算斜激波后的马赫数 M、密度、温度 T、压力 P 等物理参数 根据关系表判断在端点是否产生斜激波？计算出气流经过每一个折角需偏转的角度 w，根据等熵流动函数计算各段的马赫数 M 根据前后段马赫数的关系及等熵流动函数计算各段的物理参数 T,P,Cp 并计算气动力特性参数 Cm,Ca,Cn,Cl,Cd.报错 NO YES 计算结果及分析：（一）结果曲线及图表（1）表面压力 Cp 分布曲线 Cpx 图如下 Figure 1 0攻角 Cpx Figure 2 2攻角 Cpx（2）表面密度曲线/如
11、下：Figure 3 分段数为 5 时，2、0攻角下的/-x Figure 4 分段数为 20 时，2、0攻角下的/-x Figure 5 分段数为 50 时，2、0攻角下的/-x（3）表面温度曲线 T/T如下 Figure 6 分段数为 5 时，2、0攻角下的 T/T-x Figure 7 分段数为 20 时，2、0攻角下的 T/T-x Figure 8 分段数为 50 时，2、0攻角下的 T/T-x （4）轴向力 Ca、法向力 Cn、绕头部顶点俯仰力矩 Cm 及升力 Cl、阻力 Cd 表格 0攻角 2攻角 5 等分 20等分 50等分 5 等分 20等分 50等分 轴 向 力Ca 0.06
12、487 0.06534 0.06573 0.06543 0.06965 0.07121 表格 1 轴向力 Ca 0攻角 2攻角 5 等分 20等分 50等分 5 等分 20等分 50等分 法向力 Cn 0 0 0 0.075951 0.073747 0.077441 表格 2 法向力 Cn 0攻角 2攻角 5 等分 20等分 50等分 5 等分 20等分 50等分 俯 仰 力 矩Cm 0 0 0-0.394322-0.030473-0.030586 表格 3 俯仰力矩 Cm 0攻角 2攻角 5 等分 20等分 50等分 5 等分 20等分 50等分 Cl 0 0 0 0.073621 0.07
13、1272 0.074908 Cd 0.06487 0.06534 0.06573 0.06805 0.07218 0.07387 表格 4 升力系数 Cl、阻力系数 Cd （二）结果分析 对上面列写的表格及曲线图分析如下：（1）计算精度：在各攻角下，当取不同等分进行离散化时，计算所得曲线有一定差异，尤其是分段数为 5 的计算结果与另外两种分段数下的计算结 果会有一定误差，对于/、T/T、Cp 而言偏差在 4%左右，而分段数为 20与分段数为 50 相比，偏差都在 1%以内，因此可以认为分段数为 5 不能作为精确解，而分段数为 20 或是 50 在一定程度上能够代表精确解。（2）随着 x 值逐渐
14、增加，气流不断经膨胀波，气流马赫数逐渐增加，上下翼面升力系数逐渐减少，证明升力绝大部分由翼型前缘部分产生，与实验所得到的情况相符合。（3）分析曲线可知，从翼型前缘到后缘，由于产生膨胀波，翼面气流压强，密度，温度逐渐减小。这是符合超音速等熵流动的特点的。（4）攻角的影响：从攻 Figure1 到 Figure8 可以看出：相比于下翼面，上翼面的压力更小，压力系数 Cp 更小，温度与密度也都要小于下翼面的温度、压强。攻角大于 0 时，相对于来流方向，上翼面机翼前缘的半顶角更大，产生的斜激波更弱，而下翼面前缘产生的斜激波更强，而经过斜激波后的流动与攻角无关，所以上翼面的流速更高，下翼面的流速更小。从
15、而使上翼面的压力系数、温度、密度小于下翼面。编程计算所得的数据反映出这种规律。（5）从表 1 到表 4 可以看出：在 0攻角下，由于上下翼面的流动参数对称，翼型上的坐标系 x 轴与空气流动方向相重合，翼面法向（y 轴）与空气流动的法向相重合。受力在进行矢量相加后，得出法向力系数 Cn=0，此时升力系数 Cl=Cn=0。阻力系数 Cd=Ca0。当攻角为 2时，Cl，Cn，Cd，Ca0，为保证翼型平衡，此时压力对前缘的俯仰矩0，是使翼型低头的矩。可见轴对称翼型在有一定的攻角的条件下才会提供升力。优缺点分析及个人体会 优点：本文建立的轴对称翼型离散化模型在分段数较大（20100）时，能较精确地模拟计
16、算出绕该翼型 的空气流动参数及气动力特性。所编写程序能做到根据输入的平面翼型方程、分段数、来流马赫数、及攻角大小即能判断出是否在机翼前缘产生附体斜激波，并根据激波膨胀波法求解出该绕翼型的空气物理参数及该翼型的气动力特性。缺点：有一个问题没有得到解决，在离散步长趋于无限小时，空气在机翼尾缘的马赫数会无限增大。这是不符合物理规律的。个人体会：编制 P-M 表、及等熵流动关系表时，一定要保证有相当高的精度，刚开始我编制的 P-M 表精度仅仅达到了 0.1 度，0.1 度的误差会导致斜激波后马赫数出现 0.20.3 的绝对误差，导致温度、密度、压强的计算结果会出现较大的误差。将 P-M表精度调至 0.
17、01时，精度会有很大提高。附录：程序%*输入 N=?alpha=?*%*编制马赫数为 2 时的斜激波 BataCita关系表*%b=28:0.1:70;Bata2Cita=zeros(421,2);for i=1:421 Bata2Cita(i,1)=b(i,1);Bata2Cita(i,2)=180/pi*atan(2*(1/tan(b(i)*pi/180)*(4*sin(b(i)*pi/180)2-1)/(4*(1.4+cos(2*b(i)*pi/180)+1);end%*编制 PM 关系表*%m=1:0.001:4;x=zeros(3001,1);for n=1:3001 end PM=
18、zeros(3001,2);for n=1:3001 PM(n,2)=(180/pi)*(6)0.5)*atan(1/6)*(m(n,1)2-1)0.5)-atan(m(n,1)2)-1)0.5);PM(n,1)=m(n,1);end%-*考虑攻角*-%*第一步，离散化*%x=0:1/N:1;for i=1:N+1 y1(i)=x(i)*(1-x(i)*0.3;y2(i)=-y1(i);end for i=1:N w1(i)=(180/pi)*atan(y1(i+1)-y1(i)/(x(i+1)-x(i);w2(i)=-(180/pi)*atan(y2(i+1)-y2(i)/(x(i+1)-x
19、(i);end%*w2(i)任然是取正值，是为了后面各项物理参数的计算方便，*%*但是在求解力的作用时，就应该充分考虑清楚力的方向性*%bata1=1;bata2=1;w1(1)=w1(1)-alpha;w2(1)=w2(1)+alpha;M=zeros(2,N+1);M(1,1)=2;M(2,1)=2;%*分别计算上下面对应的斜激波角 Bata，*%*对于不符合产生斜激波条件的状况报错*%for i=1:421 if(abs(w1(1)-Bata2Cita(i,2)=0.05)bata1=Bata2Cita(i,1);end end for i=1:421 if(abs(w2(1)-Bata
20、2Cita(i,2)=0.05)bata2=Bata2Cita(i,1);end end if(bata1=1)|(bata2=1)disp error else disp correct end%*进行各段的物理参数计算*%M(1,2)=xiejiboMahe(2,bata1);M(2,2)=xiejiboMahe(2,bata2);w1(1)=w1(1)+alpha;w2(1)=w2(1)-alpha;for i=3:N+1%*循环次数*%q=M2cita(M(1,i-1)+w1(i-2)-w1(i-1);for n=1:3001 if(abs(q-PM(n,2)=0.04)M(1,i)=
21、PM(n,1);end end end for i=3:N+1 q=M2cita(M(2,i-1)+w1(i-2)-w1(i-1);for n=1:3001 if(abs(q-PM(n,2)=0.04)M(2,i)=PM(n,1);end end end T=zeros(2,N+1);P=zeros(2,N+1);mi=zeros(2,N+1);T(1,1)=1;P(1,1)=1;mi(1,1)=1;T(2,1)=1;P(2,1)=1;mi(2,1)=1;T(1,2)=xiejiboWendu(2,bata1);P(1,2)=xiejiboYaqiang(2,bata1);mi(1,2)=xi
22、ejiboMidu(2,bata1);T(2,2)=xiejiboWendu(2,bata2);P(2,2)=xiejiboYaqiang(2,bata2);mi(2,2)=xiejiboMidu(2,bata2);for i=2:N T(1,i+1)=T(1,i)*(1+(0.2*(M(1,i)2)/(1+(0.2*(M(1,i+1)2);T(2,i+1)=T(2,i)*(1+(0.2*(M(2,i)2)/(1+(0.2*(M(2,i+1)2);P(1,i+1)=P(1,i)*(1+(0.2*(M(1,i)2)/(1+(0.2*(M(1,i+1)2)(7/2);P(2,i+1)=P(2,i)
23、*(1+(0.2*(M(2,i)2)/(1+(0.2*(M(2,i+1)2)(7/2);mi(1,i+1)=mi(1,i)*(1+(0.2*(M(1,i)2)/(1+(0.2*(M(1,i+1)2)(5/2);mi(2,i+1)=mi(2,i)*(1+(0.2*(M(2,i)2)/(1+(0.2*(M(2,i+1)2)(5/2);end%*下面着力解决几个力学参数的计算*%*先求压力系数*%Cp=zeros(2,N+1);for i=1:N+1 Cp(1,i)=(P(1,i)-1)/(0.7*M(1)2);Cp(2,i)=(P(2,i)-1)/(0.7*M(1)2);end%*求机翼轴向力和法
24、向力*%*先求上下翼面的力，求合力时注意方向*%Fx=zeros(1,3);Fy=zeros(1,3);for i=1:N Fx(1,1)=Fx(1,1)-P(1,i+1)*(1/N)*tan(w1(i)*pi/180);Fy(1,1)=Fy(1,1)-P(1,i+1)*(1/N);Fx(1,2)=Fx(1,2)-P(2,i+1)*(1/N)*tan(w2(i)*pi/180);Fy(1,2)=Fy(1,2)-P(2,i+1)*(1/N);end Fx(1,3)=Fx(1,1)+Fx(1,2);%Fx(1,3)表示上下翼面的总轴向力；Fy(1,3)=-Fy(1,2)+Fy(1,1);%Fy(1
25、,3)表示上下翼面的总法向力；%*接着求解升力和阻力*%Fs=zeros(1,1);Fz=zeros(1,1);Fs=Fy(1,3)*cos(alpha*pi/180)+Fx(1,3)*sin(alpha*pi/180);Fz=Fx(1,3)*cos(alpha*pi/180)-Fy(1,3)*sin(alpha*pi/180);%*求俯仰力矩，要有每个面上的作用力的作用点，假设为每段中点*%Moment=zeros(2,N);for i=1:N X(i)=0.5*(x(i)+x(i+1);Y1(i)=0.5*(y1(i)+y1(i+1);Y2(i)=0.5*(y2(i)+y2(i+1);en
26、d for i=1:N Moment(1,i)=-P(1,i+1)*(1/N)*tan(w1(i)*pi/180)*Y1(i)-P(1,i+1)*(1/N)*X(i);Moment(2,i)=-P(2,i+1)*(1/N)*tan(w2(i)*pi/180)*Y2(i)+P(2,i+1)*(1/N)*X(i);end%*合力矩计算公式的加减符号*%summoment=0;for i=1:N summoment=summoment+Moment(1,i)+Moment(2,i);end%*最后算出要求的量 Ca,Cn,Ci,Cd,Cm*%Ca=Fx(1,3)/(0.5*M(1)2*1.4);Cn=Fy(1,3)/(0.5*M(1)2*1.4);Ci=Fs/(0.5*M(1)2*1.4);Cd=Fz/(0.5*M(1)2*1.4);Cm=summoment/(0.5*M(1)2*1.4);