《高超声速气动力及导数计算报告》由会员分享，可在线阅读，更多相关《高超声速气动力及导数计算报告（38页珍藏版）》请在金锄头文库上搜索。
1、高超声速气动力及导数计算报告撰写人：学 号 : 班 级：2012 年 10 月 25 日一 引言实验目的：熟练运用面元法中有关网格划分的方法；掌握高超音速气动力工程估算方法中面元法及牛顿法； 同时， 比较两者的计算结果，并分析差异产生原因。实验条件：编程计算该旋成体的升力、 阻力、 升阻比及俯仰力矩系数， 还有导数作出曲线。图（ 1）计算条件： 5 30 , 4.5 10.5, 10 , 0M H km几何参数： 0.1 , 0.5 , 3.912a m b m D m飞行器运动状态： 0P Q R二 计算方法规定导弹的体坐标系为： X 轴沿着导弹纵轴向后， Y 轴垂直于弹体向上， Z轴于其它
2、两轴构成右手坐标系， 即指向导弹左侧， 原点位于导弹前缘点。 取原点为参考点。牛顿法：牛顿法的假设如下：攻角 小于物面倾角 ；假设流体由大量均匀分布的， 彼此独立无相互作用的质点所组成， 它们排列整齐、平行地沿着直迹线流向物体。流体质点流与物面碰撞时， 流体质点将失去与物面垂直的法向动量， 而保持原有的切向动量沿物面向下流下去。 由于法向动量的变化从而引起流体作用在物体上的力。流体对物面的压力只作用在物面能与流体质点相碰撞的表面上， 而遮蔽区上压力为零。牛顿公式：22 sinpC其中 为来流速度方向与物面切面的夹角。由课本可知：2 2 2 22 0m ax2 22 12(cos sin sin
3、 cos )( )24 cos sin cos sin sin sin cos cos sin krxt kk k kCrrdr2 22 0max2 2 2 2 32 2 cos sin cos sin cos sin 2sin (cos sin2sin cos ) 2cos sin cos sin ( ) sin cos sin 3kxyt k k kk kCrrdx对于圆球部分，由于用于验证，我们假设 0 T1=x(rem(i,m)+1,ceil(i/m)+1)-x(rem(i,m),ceil(i/m) y(rem(i,m)+1,ceil(i/m)+1)-y(rem(i,m),ceil(i
4、/m) z(rem(i,m)+1,ceil(i/m)+1)-z(rem(i,m),ceil(i/m); T2=x(rem(i,m),ceil(i/m)+1)-x(rem(i,m)+1,ceil(i/m) y(rem(i,m),ceil(i/m)+1)-y(rem(i,m)+1,ceil(i/m) z(rem(i,m),ceil(i/m)+1)-z(rem(i,m)+1,ceil(i/m); xbb=(x(rem(i,m)+1,ceil(i/m)+1)+x(rem(i,m),ceil(i/m)+x(rem(i,m),ceil(i/m)+1)+x(rem(i,m)+1,ceil(i/m)/4; y
5、bb=(y(rem(i,m)+1,ceil(i/m)+1)+y(rem(i,m),ceil(i/m)+y(rem(i,m),ceil(i/m)+1)+y(rem(i,m)+1,ceil(i/m)/4; zbb=(z(rem(i,m)+1,ceil(i/m)+1)+z(rem(i,m),ceil(i/m)+z(rem(i,m),ceil(i/m)+1)+z(rem(i,m)+1,ceil(i/m)/4; %记录 i 不等于 m时面元 xyz 平均值xyzk(i,1)=x(rem(i,m),ceil(i/m) y(rem(i,m),ceil(i/m) z(rem(i,m),ceil(i/m) x(
6、rem(i,m)+1,ceil(i/m) y(rem(i,m)+1,ceil(i/m) z(rem(i,m)+1,ceil(i/m) x(rem(i,m),ceil(i/m)+1) y(rem(i,m),ceil(i/m)+1) z(rem(i,m),ceil(i/m)+1) x(rem(i,m)+1,ceil(i/m)+1) y(rem(i,m)+1,ceil(i/m)+1) z(rem(i,m)+1,ceil(i/m)+1); else T1=x(1,ceil(i/m)+1)-x(m,ceil(i/m) y(1,ceil(i/m)+1)-y(m,ceil(i/m) z(1,ceil(i/m
7、)+1)-z(m,ceil(i/m); T2=x(m,ceil(i/m)+1)-x(1,ceil(i/m) y(m,ceil(i/m)+1)-y(1,ceil(i/m) z(m,ceil(i/m)+1)-z(1,ceil(i/m); xbb=(x(1,ceil(i/m)+1)+x(m,ceil(i/m)+x(m,ceil(i/m)+1)+x(1,ceil(i/m)/4; ybb=(y(1,ceil(i/m)+1)+y(m,ceil(i/m)+y(m,ceil(i/m)+1)+y(1,ceil(i/m)/4; zbb=(z(1,ceil(i/m)+1)+z(m,ceil(i/m)+z(m,cei
8、l(i/m)+1)+z(1,ceil(i/m)/4; % 记录 i=m 面元 xyz 平均值xyzk(i,1)=x(m,ceil(i/m) y(m,ceil(i/m) z(m,ceil(i/m) x(1,ceil(i/m) y(1,ceil(i/m) z(1,ceil(i/m) x(m,ceil(i/m)+1) y(m,ceil(i/m)+1) z(m,ceil(i/m)+1) x(1,ceil(i/m)+1) y(1,ceil(i/m)+1) z(1,ceil(i/m)+1); end xb(i,1)=xbb; yb(i,1)=ybb; zb(i,1)=zbb; % 记录面 元xyz 平均值
9、 t(i,:)=T1; t1(i,1)=t(i,1)/sqrt(t(i,1)2+t(i,2)2+t(i,3)2); t1(i,2)=t(i,2)/sqrt(t(i,1)2+t(i,2)2+t(i,3)2); t1(i,3)=t(i,3)/sqrt(t(i,1)2+t(i,2)2+t(i,3)2); %求 出t1 矢量N(i,:)=cross(T2,T1); n1(i,1)=N(i,1)/sqrt(N(i,1)2+N(i,2)2+N(i,3)2); n1(i,2)=N(i,2)/sqrt(N(i,1)2+N(i,2)2+N(i,3)2); n1(i,3)=N(i,3)/sqrt(N(i,1)2+
10、N(i,2)2+N(i,3)2); %求 出平面面元法向量 N 及其单位矢量 n m1(i,:)=cross(n1(i,:),t1(i,:); end xyzb=xb yb zb; cha=cell(4*m*n,1); for i=1:4*n*m for ii=1:4 chai,1(ii,1:3)=xyzb(i,1:3)-xyzki,1(ii,1:3); end end % 求点到面元平均点的差dk=zeros(4*n*m,4); for i=1:4*n*m for j=1:4 dk(i,j)=dot(n1(i,:),chai,1(j,:); end end % 求 dk xyz=cell(4
11、*m*n,1); for i=1:4*n*m for ii=1:3 for j=1:4 xyzi,1(j,ii)=xyzki,1(j,ii)+n1(i,ii)*dk(i,j); end end end %求面元四个角点投影到面元平面上的坐标 xyz xyzpk=cell(4*m*n,1); for i=1:4*n*m for j=1:4 xyzpki,1(j,1)=t1(i,1)*(xyzi,1(j,1)-xyzb(i,1)+t1(i,2)*(xyzi,1(j,2)-xyzb(i,2)+t1(i,3)*(xyzi,1(j,3)-xyzb(i,3); xyzpki,1(j,2)=m1(i,1)*
12、(xyzi,1(j,1)-xyzb(i,1)+m1(i,2)*(xyzi,1(j,2)-xyzb(i,2)+m1(i,3)*(xyzi,1(j,3)-xyzb(i,3); xyzpki,1(j,3)=0; end end %求面元坐标系中四个角点坐标 xyzpk A=zeros(1,4*m*n); for i=1:4*n*m A(1,i)=0.5*(xyzpki,1(4,1)-xyzpki,1(1,1)* (xyzpki,1(2,2)-xyzpki,1(3,2); %由于标注问题标号为 3 和 4 的点与课文相反end % 求面元面积xyzp0=zeros(4*m*n,3); for i=1:
13、4*m*n xyzp0(i,1)=(xyzpki,1(1,2)-xyzpki,1(2,2)*xyzpki,1(3,1)+(xyzpki,1(3,2)-xyzpki,1(1,2)*xyzpki,1(2,1)/3/(xyzpki,1(2,2)-xyzpki,1(3,2); xyzp0(i,2)=-xyzpki,1(1,2)/3; xyzp0(i,3)=0; end % 求质心坐标 xyzp0 xyz0=zeros(4*m*n,3); for i=1:4*m*n xyz0(i,1)=xyzb(i,1)+xyzp0(i,1)*t1(i,1)+xyzp0(i,2)*m1(i,1); xyz0(i,2)=
14、xyzb(i,2)+xyzp0(i,1)*t1(i,2)+xyzp0(i,2)*m1(i,2); xyz0(i,3)=xyzb(i,3)+xyzp0(i,1)*t1(i,3)+xyzp0(i,2)*m1(i,3); end % 质心在飞行器坐标系中的坐标 xyz0 v=Ma*340*cos(alf)*cos(blt) Ma*340*sin(alf)*cos(blt) Ma*340*sin(blt); % 无穷远处的来流速度vcv=zeros(4*m*n,1); vcva=zeros(4*m*n,1); vcvb=zeros(4*m*n,1); for i=1:4*m*n vcv(i,1)=-n1(i,1)*cos(alf)*cos(blt)-n1(i,2
