超音速翼型气动力特性研究汇总.docx
- 文档编号:28445553
- 上传时间:2023-07-13
- 格式:DOCX
- 页数:18
- 大小:679.75KB
超音速翼型气动力特性研究汇总.docx
《超音速翼型气动力特性研究汇总.docx》由会员分享,可在线阅读,更多相关《超音速翼型气动力特性研究汇总.docx(18页珍藏版)》请在冰豆网上搜索。
超音速翼型气动力特性研究汇总
超音速翼型气动力特性研究
摘要:
本文研究方程为
的轴对称超音速翼形在马赫数为2,攻角分别为0°,2°情形下的气动力特性,基于对翼型进行离散化处理得到该翼型的物理参数及气动力的近似解,并逐步减小空间步长
来提高解的精度。
在步长数分别为5、20、50及攻角为0°、2°的条件下,计算求得翼型头部斜激波后的流动参数,并由此求解各分区相应参数,列出:
表面压力Cp分布曲线Cp-x,及表面密度、温度分布曲线/-x、T/T-x。
在不同条件下得出的轴向力Ca、法向力Cn、升力Cl、阻力Cd及绕头部顶点俯仰力矩Cm的表格。
最终分析了编程计算的准确性与精度,分析了压力系数、温度、密度沿该翼型的分布特性,并分析了不同攻角对该翼型气动特性的影响。
问题描述
已知方程为
的薄翼形,求该翼型在来流马赫数为2,攻角分别为0°,2°情形下的受力情况。
对x范围(0,1)内分别按5等份、20等份和50等份进行离散计算,得到表面压力Cp分布曲线Cp-x,表面密度、温度分别曲线/、T/T。
计算得出出轴向力Ca、法向力Cn、绕头部顶点俯仰力矩Cm及升力Cl、阻力Cd。
计算方案:
(一)计算思路:
超音速来流以一定攻角遇到类似于楔形体的机翼前缘,在上下面都有可能产生附体斜激波,要是攻角过大也有可能不产生附体斜激波,这里首先需要根据斜激波的
关系曲线图来作出判断。
经判断,如果顶点处产生斜激波,即使用斜激波前后的马赫数、密度、温度、压强计算公式计算出顶点斜激波后的各项物理参数。
接着,根据翼型的形状可知,气流在通过膨胀波之后会经过一系列的向外的转折角,根据普朗特-迈耶膨胀波理论,超音速气流经过每一个折角都会产生膨胀波。
根据数值计算的基本原理,计算机不能处理连续曲线上随x值变化而连续变化的折角,所以在计算之前必须对翼型的几何结构进行离散化处理。
离散化之后即可根据膨胀波前后马赫数的关系公式计算出每一个折角处膨胀波后的马赫数,然后根据膨胀波前后密度、温度、压强的计算公式计算出每一个膨胀波后的密度、温度、压强。
得到以上基本的物理参数之后,即可用压强P的分布计算出压力系数Cp的分布进而计算出翼型所受的轴向力、法向力、升力、阻力及力矩系数。
在进行以上计算之前的首要工作是编制P-M表、等熵流动角度与马赫数的关系表。
在具体操作中可使用已知的显示函数式进行编制,无需手动输入。
(二)基本假设:
实际上的翼型气动力会受到很多因素的干扰,用激波-膨胀波法计算时对实际的物理模型做了一些简化,假设:
1、离散化之后任意两个离散点之间的物理参数是均匀分布的。
2、不考虑斜激波、膨胀波的相交与反射。
3、翼型端部的斜激波不会在上部削弱,即斜激波不会滑移。
(三)坐标系建立
建立坐标系时,以机翼的前段点为原点O,弦线为x轴,垂直于弦线的直线为y轴,x轴正向指向机翼尾缘,y轴正向垂直于x轴向上,如图所示:
在进行力学参数计算时,也是以X轴正向为正,z轴正向为正。
(四)物面离散
物面离散的步骤如下(以分成5段为例):
首先在X(0~1)上以0.2为步长取出0、0.2、0.4、0.6、0.8、1六个点,再根据
计算出x=0、0.2、0.4、0.6、0.8、1时所对应的z值,对上述的12个点进行画图,即可得到以5段离散得到的翼型图。
分20段、50段进行离散的步骤与分5段的步骤相同。
各个离散物面与水平面的夹角θ也可以求出,以方便后来的气流偏转角的计算,计算公式如下:
(五)各区域物理参数的计算
(1)第一区流动参数的计算
按照分段数的要求将物面离散后,首先需要计算出各个区段的物理参数,特别是马赫数,因为斜激波、膨胀波前后的马赫数关系是求解其他物理参数的基础。
由于存在攻角的缘故,超音速气流流经上下翼面的第一个偏转角并不等于半顶角θ1。
超音速气流遇到机翼前缘时的偏转角w计算公式如下:
对于上翼面:
对于下翼面:
然后,根据斜激波的θ-β关系表来作出判断超音速来流遇到机翼前缘是否会产生附体斜激波。
程序会查找来流马赫数M1为2时的θ-β关系表找出此θ对应的β,如果没有找到与之对应的β,即说明该情况下,在机翼前缘不会产生斜激波,程序会报错。
如果产生了斜激波,由已知的M1和β计算出斜激波后的马赫数、温度、密度、压强。
计算公式如下:
(2)第i区气流参数的计算
在一区之后的每区前端点均产生膨胀波,第i区(
)相对音速来流的折转角为
通过对P-M表进行插值,根据
查等熵流动函数表可以求出第i区对应的马赫数
。
由于膨胀波为等熵波,波前后流动参数满足等熵关系式:
然后,即根据P的分布推算Cp的分布,计算公式如下:
(六)离散元力/力矩计算公式、翼型合力/合力矩计算公式
离散化后,根据最初的假设每区流动参数均匀分布,压力作用在分区中点。
轴向力为x方向上压力的分力,法向力为y方向的分力,总的轴向力与法向力计算公式如下
升力、阻力与轴向力、法向力间的差距仅仅由攻角造成,可以得到以下的计算公式:
每一段上的压力对于机翼前缘端点都会产生一个俯仰力矩,取逆时针(抬头)为正方向,计算公式如下:
上翼面:
下翼面:
总的俯仰力矩系数计算公式如下:
输入分段数N,攻角alpha
攻角alpha
编制P-M表、等熵流动函数表
离散化,计算上下翼面各分段处的偏转角θ1,θ2
计算考虑攻角后的气流遇到翼型前缘端点后流经上下翼面分别对应的偏转角w
计算斜激波后的马赫数M、密度ρ、温度T、压力P等物理参数
根据𝜃−𝛽关系表判断在端点是否产生斜激波?
计算出气流经过每一个折角需偏转的角度w,
根据等熵流动函数计算各段的马赫数M
根据前后段马赫数的关系及等熵流动函数计算各段的物理参数T,ρ,P,Cp并计算气动力特性参数Cm,Ca,Cn,Cl,Cd.
报错
NO
NO
YES
程序流程图
结束
计算结果及分析:
(一)结果曲线及图表
(1)表面压力Cp分布曲线Cp-x图如下
Figure10°攻角Cp-x
Figure22°攻角Cp-x
(2)表面密度曲线/如下:
Figure3分段数为5时,2°、0°攻角下的/-x
Figure4分段数为20时,2°、0°攻角下的/-x
Figure5分段数为50时,2°、0°攻角下的/-x
(3)表面温度曲线T/T如下
Figure6分段数为5时,2°、0°攻角下的T/T-x
Figure7分段数为20时,2°、0°攻角下的T/T-x
Figure8分段数为50时,2°、0°攻角下的T/T-x
(4)轴向力Ca、法向力Cn、绕头部顶点俯仰力矩Cm及升力Cl、阻力Cd表格
0°攻角
2°攻角
5等分
20等分
50等分
5等分
20等分
50等分
轴向力Ca
0.06487
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.071272
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值逐渐增加,气流不断经膨胀波,气流马赫数逐渐增加,上下翼面升力系数逐渐减少,证明升力绝大部分由翼型前缘部分产生,与实验所得到的情况相符合。
(3)分析曲线可知,从翼型前缘到后缘,由于产生膨胀波,翼面气流压强,密度,温度逐渐减小。
这是符合超音速等熵流动的特点的。
(4)攻角的影响:
从攻Figure1到Figure8可以看出:
相比于下翼面,上翼面的压力更小,压力系数Cp更小,温度与密度也都要小于下翼面的温度、压强。
攻角大于0时,相对于来流方向,上翼面机翼前缘的半顶角更大,产生的斜激波更弱,而下翼面前缘产生的斜激波更强,而经过斜激波后的流动与攻角无关,所以上翼面的流速更高,下翼面的流速更小。
从而使上翼面的压力系数、温度、密度小于下翼面。
编程计算所得的数据反映出这种规律。
(5)从表1到表4可以看出:
在0°攻角下,由于上下翼面的流动参数对称,翼型上的坐标系x轴与空气流动方向相重合,翼面法向(y轴)与空气流动的法向相重合。
受力在进行矢量相加后,得出法向力系数Cn=0,此时升力系数Cl=Cn=0。
阻力系数Cd=Ca>0。
当攻角为2°时,Cl,Cn,Cd,Ca>0,为保证翼型平衡,此时压力对前缘的俯仰矩<0,是使翼型低头的矩。
可见轴对称翼型在有一定的攻角的条件下才会提供升力。
优缺点分析及个人体会
优点:
本文建立的轴对称翼型离散化模型在分段数较大(20~100)时,能较精确地模拟计算出绕该翼型的空气流动参数及气动力特性。
所编写程序能做到根据输入的平面翼型方程、分段数、来流马赫数、及攻角大小即能判断出是否在机翼前缘产生附体斜激波,并根据激波膨胀波法求解出该绕翼型的空气物理参数及该翼型的气动力特性。
缺点:
有一个问题没有得到解决,在离散步长趋于无限小时,空气在机翼尾缘的马赫数会无限增大。
这是不符合物理规律的。
个人体会:
编制P-M表、及等熵流动关系表时,一定要保证有相当高的精度,刚开始我编制的P-M表精度仅仅达到了0.1度,0.1度的误差会导致斜激波后马赫数出现0.2~0.3的绝对误差,导致温度、密度、压强的计算结果会出现较大的误差。
将P-M表精度调至0.01°时,精度会有很大提高。
附录:
程序
%%*****输入N=?
alpha=?
*****%%
%%*****编制马赫数为2时的斜激波Bata~Cita关系表*****%
b=[28:
0.1:
70]';
Bata2Cita=zeros(421,2);
fori=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);
forn=1:
3001
end
PM=zeros(3001,2);
forn=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];
fori=1:
N+1
y1(i)=x(i)*(1-x(i))*0.3;
y2(i)=-y1(i);
end
fori=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(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,*****%%
%%****对于不符合产生斜激波条件的状况报错*****%%
fori=1:
421
if(abs(w1
(1)-Bata2Cita(i,2))<=0.05)
bata1=Bata2Cita(i,1);
end
end
fori=1:
421
if(abs(w2
(1)-Bata2Cita(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;
fori=3:
N+1%%***循环次数***%%
q=M2cita(M(1,i-1))+w1(i-2)-w1(i-1);
forn=1:
3001
if(abs(q-PM(n,2))<=0.04)
M(1,i)=PM(n,1);
end
end
end
fori=3:
N+1
q=M2cita(M(2,i-1))+w1(i-2)-w1(i-1);
forn=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)=xiejiboMidu(2,bata1);
T(2,2)=xiejiboWendu(2,bata2);
P(2,2)=xiejiboYaqiang(2,bata2);
mi(2,2)=xiejiboMidu(2,bata2);
fori=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)*((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);
fori=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
%%*****求机翼轴向力和法向力*****%%
%%*****先求上下翼面的力,求合力时注意方向*****%%
Fx=zeros(1,3);
Fy=zeros(1,3);
fori=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,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);
fori=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));
end
fori=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;
fori=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);
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 超音速 气动力 特性 研究 汇总