《高超声速气动力及导数计算报告.doc》由会员分享，可在线阅读，更多相关《高超声速气动力及导数计算报告.doc（37页珍藏版）》请在三一办公上搜索。
1、高超声速气动力及导数计算报告撰写人： 学 号: 班 级： 2012年10月25日一 引言实验目的：熟练运用面元法中有关网格划分的方法；掌握高超音速气动力工程估算方法中面元法及牛顿法；同时，比较两者的计算结果，并分析差异产生原因。实验条件：编程计算该旋成体的升力、阻力、升阻比及俯仰力矩系数，还有导数作出曲线。图（1）计算条件：几何参数：飞行器运动状态：二 计算方法规定导弹的体坐标系为：轴沿着导弹纵轴向后，轴垂直于弹体向上，轴于其它两轴构成右手坐标系，即指向导弹左侧，原点位于导弹前缘点。取原点为参考点。牛顿法：牛顿法的假设如下：攻角小于物面倾角；假设流体由大量均匀分布的，彼此独立无相互作用的质点所
2、组成，它们排列整齐、平行地沿着直迹线流向物体。流体质点流与物面碰撞时，流体质点将失去与物面垂直的法向动量，而保持原有的切向动量沿物面向下流下去。由于法向动量的变化从而引起流体作用在物体上的力。流体对物面的压力只作用在物面能与流体质点相碰撞的表面上，而遮蔽区上压力为零。牛顿公式：其中为来流速度方向与物面切面的夹角。由课本可知：对于圆球部分，由于用于验证，我们假设,:对于圆锥部分，由于每一部分圆锥各自的一直为常值，因此直接带入积分进行计算，仅与有关，当时，.其余为零。由于用于验证，我们假设,因此可得对于圆锥其相关系数为：对于圆柱部分，由于，因此，仅计算其法向力系数，可得于是对于本实验：面元法：1）
3、确定平面面元小曲面元四个角点的坐标为。其中，两个对角点形成矢量，分别把它们标记为这里通过两个矢量的矢量积，可以得到一个新的矢量，把它记为：的单位矢量定义为：我们知道：由一个法向矢量和空间任一点可以确定一个平面。因此，把矢量作为法向矢量，并选取曲面的四个角点的平均值作为一个空间点,从而可以确定一个平面，这就是我们要选取的面元所在的平面。其中为面元的外法线方向单位矢量，由如下公式计算：现在，将四个角点沿矢量的方向投影到面元平面上，可以求得四个投影点，其坐标分别为式中，2）建立面元坐标系在面元形成后，需要建立一个面元坐标系。现在我们选取面元法向单位矢量作为一个单位矢量，再把的单位矢量作为一个单位矢量
4、：由和的矢量积可以决定第三的单位矢量：至此，已经建立了面元坐标系，原点取，轴平行于，轴平行于，轴平行于。易知轴和轴在面元平面上。3）面元坐标系与飞行器参考坐标系之间的变换在具体计算中，需要把飞行器坐标系中的坐标值变换到面元坐标系中去，也需要把面元系中的值转换到飞行器坐标系中去。此时，转换公式为：飞行器坐标转换到面元坐标从面元坐标系转换到飞行器坐标系4）面元的面积和质心坐标现在把投影点在飞行器坐标系中的值变换到面元坐标系中去，在面元坐标系中的值用来表示，则在面元平面上，由凸四边形面积公式可求得面元的面积为这里由于平行于轴，因此，面积公式可以化简为在面元坐标系中，面元质心的坐标为利用面元坐标系转换
5、到飞行器坐标系的转换关系，可求得质心在飞行器坐标中的质心坐标5）面元冲击角在气动力的计算中，我们必须要知道平面面元与自由流速度的夹角，即冲击角。而第个面元平面的冲击角：式中为第面元平面的单位法向向量。速度矢量式中，总角速度 向径 ：为来流速度，为物体滚转、偏航、俯仰角速度，为参考点的坐标，在编程过程中取（0，0，0）。6）计算面元压力系数根据冲击角是否大于零，可以把面元所处的流场分为迎风面和背风面。若冲击角大于零，则可利用修正牛顿理论：，计算面元的压力系数。其中为来流速度与面元内法线的夹角若冲击角小于零，则可直接令。7） 面元上气动力系数及动静导数计算将飞行器表面分成若干面元后，用积分的方法可
6、以得到它的气动力系数及其静、动导数。参考坐标系中（轴沿着导弹纵轴向后，轴垂直于弹体向上，轴于其它两轴构成右手坐标系，即指向导弹左侧，原点位于导弹前缘点。），来流速度矢量，假设面元外法向矢量，则来流速度在面元内法向上的投影为：分别对和求导，得物体滚转、俯仰、偏航角速度为，由三个角速度引起的物体与气流相对运动的三个速度为：其中表示面元质心到参考点的距离。令在内法向的投影为，若记对的导数为式中参考面积。总速度在面元法向的投影：，则修正的牛顿理论亦可表示为。根据上述计算，压力系数的静、动导数计算公式为：根据前面计算的面元上的几何参数及压力系数和静、动导数，可用下列公式计算面元上气动力、力矩及其静、动导
7、数。面元上所受的气动力：面元上气动力对参考点的力矩：面元上气动力的静导数：面元上气动力产生的力矩的静导数：面元上气动力的动导数:面元上力矩的动导数:求出每个面元上的气动力、力矩及其静动导数以后，通过求和或积分就可求出整体气动力系数。轴向、法向和侧向力系数:轴向、法向和侧向力的导数:俯仰、偏航和滚转力矩系数:俯仰、偏航和滚转阻尼导数:其中体轴系与风轴系之间的转换矩阵：三 计算程序程序：%注意：该程序中，语句后的%为解释用，语句前的%是为加快速度没必要计算的滚转方向%由于画图未使用函数调用而采用循环速度较慢大概在20秒左右%面元法%clear allclcsyms alf blt grm p q
8、r Ma m n Dblt=0; %侧滑角grm=0; %滚转角n=5;m=15;brb=1.4; %比热比p=0; %滚转角速度q=0; %俯仰角速度r=0; %偏转角速度D=3.912;midu=0.413;S=pi*(D+tan(16.5*pi/180)*4.2*D/4*2)2/4;cfylj=zeros(3,71);cfyznlj=zeros(3,71);cphjnlj=zeros(3,71);%cgzznlj=zeros(3,71);csl=zeros(3,71);czl=zeros(3,71);cld=zeros(3,71);for iii=1:3 for jjj=1:71 Ma=
9、4.5+3*(iii-1); %Ma在4.5到10.5之间画图 alf=(-5+0.5*(jjj-1)*pi/180; %攻角在-5到30度之间画图 xx=zeros(m,4*n); yy=zeros(m,4*n); zz=zeros(m,4*n);for i=1:m xx(i,1:n)=0.1/n*(1:n); xx(i,n+1:2*n)=0.1+(4.2*D/4-0.1)/n*(1:n); xx(i,2*n+1:3*n)=4.2*D/4+4.2*D/2/5*(1:n); xx(i,3*n+1:4*n)=3*4.2*D/4+4.2*D/4/5*(1:n);end %确定每个点的x坐标for
10、j=1:4*nif xx(1,j)0.1 zz(1,j)=sqrt(1.32-(1.3-xx(1,j)2);elseif xx(1,j)4.2*D/4 zz(1,j)=0.5+(xx(1,j)-0.1)*tan(20*pi/180);elseif xx(1,j)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/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
11、)+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; ybb=(y(rem(i,m)+1,ceil(i/m)+1)+y(rem(i,m),ceil(i/m)+y(rem(i,m),ceil(i/m)+1)
12、+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(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),ce
13、il(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)+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
14、/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,ceil(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(
15、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平均值 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,
16、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+N(i,2)2+N(i,3)2); %求出平面面元法向量N及其单位矢量n m1(i,:)=cross(n1(i,:),t1(i,:);end xyzb=xb yb zb;cha=ce
17、ll(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); endend %求点到面元平均点的差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,:); endend %求dkxyz=cell(4*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 endend %求面元四个
18、角点投影到面元平面上的坐标xyzxyzpk=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)*(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; endend
19、 %求面元坐标系中四个角点坐标xyzpkA=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: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)
20、-xyzpki,1(3,2);xyzp0(i,2)=-xyzpki,1(1,2)/3;xyzp0(i,3)=0;end %求质心坐标xyzp0xyz0=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)=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 %质心在飞行器坐标系中的坐标x
21、yz0v=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)*sin(alf)*cos(blt)-n1(i,3)*sin(blt); %计算v*n/v来流vcva(i,1)=n1(i,1)*sin(alf)*cos(blt)-n1(i,2)*cos(alf)*co
22、s(blt); %计算v*n/v来流对攻角的导数 vcvb(i,1)=n1(i,1)*cos(alf)*sin(blt)+n1(i,2)*sin(alf)*sin(blt)-n1(i,3)*cos(blt); %计算v*n/v来流对侧滑角的导数end vp=zeros(4*m*n,3);vq=zeros(4*m*n,3);vr=zeros(4*m*n,3);vpc=zeros(4*m*n,1);vqc=zeros(4*m*n,1);vrc=zeros(4*m*n,1);vdp=zeros(4*m*n,1);vdq=zeros(4*m*n,1);vdr=zeros(4*m*n,1);for i=
23、1:4*m*nvp(i,:)=0 p*xyz0(i,3) -p*xyz0(i,2);vq(i,:)=-q*xyz0(i,3) 0 q*xyz0(i,1);vr(i,:)=r*xyz0(i,2) -r*xyz0(i,1) 0; %计算vp,vq,vr 注意一下符号不太确定vpc(i,1)=dot(vp(i,:),n1(i,:);vqc(i,1)=dot(vq(i,:),n1(i,:);vrc(i,1)=dot(vr(i,:),n1(i,:); %计算vp,vq,vr垂直来流的速度 注意一下符号不太确定vdp(i,1)=(xyz0(i,3)*n1(i,2)-xyz0(i,2)*n1(i,3)/D;
24、vdq(i,1)=(xyz0(i,1)*n1(i,3)-xyz0(i,3)*n1(i,1)/D;vdr(i,1)=(xyz0(i,2)*n1(i,1)-xyz0(i,1)*n1(i,2)/D; %计算vd=vp+vq+vr对p、q、r的导数 注意一下符号不太确定end vfull=zeros(4*m*n,1);for i=1:4*m*nvfull(i,1)=vcv(i,1)+vpc(i,1)+vqc(i,1)+vrc(i,1); %计算总速度在物面的总法向速度endcp02=2/(brb*Ma2)*(brb+1)*Ma2/2)(brb/(brb-1)*(brb+1)/(2*brb*Ma2-(b
25、rb-1)(1/(brb-1)-1); %计算cp02cp=zeros(4*m*n,1);for i=1:4*m*n if vcv(i,1)0 cp(i,1)=cp02*vfull(i,1)2; else cp(i,1)=0; endend %计算每个平面面元的cpdFc=zeros(4*m*n,1);dFz=zeros(4*m*n,1);dFf=zeros(4*m*n,1);%dmx=zeros(4*m*n,1);dmy=zeros(4*m*n,1);dmz=zeros(4*m*n,1);for i=1:4*m*ndFc(i,1)=-cp(i,1)*n1(i,3)*A(1,i); %计算每个
26、平面面元侧向力 注意此处没乘动压dFz(i,1)=-cp(i,1)*n1(i,1)*A(1,i); %计算每个平面面元轴向力 注意此处没乘动压dFf(i,1)=-cp(i,1)*n1(i,2)*A(1,i); %计算每个平面面元法向力 注意此处没乘动压%dmx(i,1)=xyz0(i,2)*dFc(i,1)-xyz0(i,3)*dFf(i,1); %计算每个平面面元滚转力矩 注意此处没乘动压dmy(i,1)=xyz0(i,3)*dFz(i,1)-xyz0(i,1)*dFc(i,1); %计算每个平面面元偏航力矩 注意此处没乘动压dmz(i,1)=xyz0(i,1)*dFf(i,1)-xyz0(
27、i,2)*dFz(i,1); %计算每个平面面元俯仰力矩 注意此处没乘动压enddongya=0.5*midu*(Ma*340)2; %计算动压Fc=sum(dFc)*dongya; %计算侧向力Fz=sum(dFz)*dongya; %计算轴向力Ff=sum(dFf)*dongya; %计算法向力CFc=sum(dFc)/S; %计算侧向力系数CFz=sum(dFz)/S; %计算轴向力系数CFf=sum(dFf)/S; %计算法向力系数Cx=-(-CFz*cos(alf)*cos(blt)-CFf*sin(alf)*cos(blt)-CFc*sin(blt);%计算阻力系数Cy=-CFz*
28、sin(alf)+CFf*cos(alf); %计算升力系数X=-(-Fz*cos(alf)*cos(blt)-Ff*sin(alf)*cos(blt)-Fc*sin(blt); %计算阻力Y=-Fz*sin(alf)+Ff*cos(alf); %计算升力LD=Y/X; %计算升阻比%Mx=sum(-dmx)*dongya; %计算滚转力矩My=sum(dmy)*dongya; %计算偏航力矩Mz=sum(-dmz)*dongya; %计算俯仰力矩%mx=sum(-dmx)/S/4.2/D; %计算滚转力矩系数my=sum(dmy)/S/4.2/D; %计算偏航力矩系数mz=sum(-dmz)
29、/S/4.2/D; %计算俯仰力矩系数%计算俯仰、偏航、滚转阻尼导数Cpa=zeros(4*m*n,1);Cpb=zeros(4*m*n,1);%Cpp=zeros(4*m*n,1);Cpq=zeros(4*m*n,1);Cpr=zeros(4*m*n,1);dmya=zeros(4*m*n,1);dmzb=zeros(4*m*n,1);%dmxb=zeros(4*m*n,1);dmyq=zeros(4*m*n,1);dmzr=zeros(4*m*n,1);%dmxp=zeros(4*m*n,1);for i=1:4*m*n Cpa(i,1)=2*cp02*vcv(i,1)*vcva(i,1)
30、; Cpb(i,1)=2*cp02*vcv(i,1)*vcvb(i,1); %Cpp(i,1)=2*cp02*vcv(i,1)*vdp(i,1); Cpq(i,1)=2*cp02*vcv(i,1)*vdq(i,1); Cpr(i,1)=2*cp02*vcv(i,1)*vdr(i,1); %以上是压力系数的动静导数计算式 dmya(i,1)=-Cpa(i,1)*n1(i,1)*A(1,i)*xyz0(i,3)+Cpa(i,1)*n1(i,3)*A(1,i)*xyz0(i,1); dmzb(i,1)=-Cpb(i,1)*n1(i,2)*A(1,i)*xyz0(i,1)+Cpb(i,1)*n1(i,
31、1)*A(1,i)*xyz0(i,2); %dmxb(i,1)=-Cpb(i,1)*n1(i,3)*A(1,i)*xyz0(i,2)+Cpb(i,1)*n1(i,2)*A(1,i)*xyz0(i,3); %以上是力矩静导数 dmyq(i,1)=-Cpq(i,1)*n1(i,1)*A(1,i)*xyz0(i,3)+Cpq(i,1)*n1(i,3)*A(1,i)*xyz0(i,1);dmzr(i,1)=-Cpr(i,1)*n1(i,2)*A(1,i)*xyz0(i,1)+Cpr(i,1)*n1(i,1)*A(1,i)*xyz0(i,2); %dmxp(i,1)=-Cpp(i,1)*n1(i,3)*
32、A(1,i)*xyz0(i,2)+Cpp(i,1)*n1(i,2)*A(1,i)*xyz0(i,3);%以上是力矩动导数endcma=sum(dmya)/S/4.2/D;%偏航力矩的导数cnb=sum(dmzb)/S/4.2/D;%俯仰力矩的导数%clb=sum(dmxb)/S/4.2/D;%滚转力矩的导数cmqcma=(sum(-dmyq)+sum(-dmya)/S/4.2/D; %偏航阻尼导数cnrcnb=(sum(-dmzr)+sum(-dmzb)/S/4.2/D; %俯仰阻尼导数%clp=sum(dmxp)/S/4.2/D;%滚转阻尼导数cfylj(iii,jjj)=mz;csl(ii
33、i,jjj)=Y;czl(iii,jjj)=X;cld(iii,jjj)=LD;cfyznlj(iii,jjj)=cnrcnb;cphjnlj(iii,jjj)=cmqcma;%cgzznlj(iii,jjj)=clp; endendfigure(1);surf(x;x(1,:),y;y(1,:),z;z(1,:) %网格划分图alpha=-5:0.5:30;figure(2);plot(alpha,csl(1,:),g.-);grid on;hold on;plot(alpha,csl(2,:),b*-);plot(alpha,csl(3,:),r+-);xlabel();ylabel(Y)
34、;title(升力Y与迎角、马赫数Ma关系曲线);legend(Ma=4.5,Ma=7.5,Ma=10.5); %升力Y与迎角、马赫数Ma关系曲线 figure(3);plot(alpha,czl(1,:),g.-);grid on;hold on;plot(alpha,czl(2,:),b*-);plot(alpha,czl(3,:),r+-);xlabel();ylabel(X);title(阻力X与迎角、马赫数Ma关系曲线); %阻力X与迎角、马赫数Ma关系曲线legend(Ma=4.5,Ma=7.5,Ma=10.5);figure(4);plot(alpha,cld(1,:),g.-)
35、;grid on;hold on;plot(alpha,cld(2,:),b*-);plot(alpha,cld(3,:),r+-);xlabel();ylabel(Y/X);title(升阻比Y/X与迎角、马赫数Ma关系曲线); legend(Ma=4.5,Ma=7.5,Ma=10.5); %升阻比Y/X与迎角、马赫数Ma关系曲线figure(5);plot(alpha,cfylj(1,:),g.-);grid on;hold on;plot(alpha,cfylj(2,:),b*-);plot(alpha,cfylj(3,:),r+-);xlabel();ylabel(mz);title(
36、俯仰力矩系数mz与迎角、马赫数Ma关系曲线);legend(Ma=4.5,Ma=7.5,Ma=10.5); %俯仰力矩系数mz与迎角、马赫数Ma关系曲线figure(6);plot(alpha,cphjnlj(1,:),g.-);grid on;hold on;plot(alpha,cphjnlj(2,:),b*-);plot(alpha,cphjnlj(3,:),r+-);xlabel();ylabel(偏航阻尼导数);title(偏航阻尼导数与迎角、马赫数Ma关系曲线);legend(Ma=4.5,Ma=7.5,Ma=10.5); %偏航阻尼导数与迎角、马赫数Ma关系曲线figure(7)
37、;plot(alpha,cfyznlj(1,:),g.-);grid on;hold on;plot(alpha,cfyznlj(2,:),b*-);plot(alpha,cfyznlj(3,:),r+-);xlabel();ylabel(俯仰阻尼导数);title(俯仰阻尼导数与迎角、马赫数Ma关系曲线);legend(Ma=4.5,Ma=7.5,Ma=10.5); %俯仰阻尼导数与迎角、马赫数Ma关系曲线%牛顿法%D=3.912;Rm=D/2+4.2*D/4*tan(16.5*pi/180); %旋成体截面最大半径Ma=7.5; %马赫数r1=0.5;r2=D/2;b=0.5;L=4.2*
38、D;L1=0.1;L2=L/4-L1;L3=0.5*L;L4=L/4;bg=40; %变量个数dh=0.4; %间隔Cx=zeros(1,bg+1);Cy=zeros(1,bg+1);for h=1:(bg+1) alpha=-5+dh*(h-1); %攻角 %圆球部分 Rq=1.3; %圆球部分球半径st0=acos(b/Rq); %圆球部分球心角的一半 Cxtq=r12/(Rm2)*(cosd(alpha)2+0.5*(sind(alpha)2)*r12/(Rm2);Cytq=r12/Rm2*sind(alpha)*cosd(alpha); %前锥段部分（k=0） Cxtqzhi=1/Rm2*(r22-r12)*(2*(cosd(alpha)2*(sind(20)2+ (sind(alpha)2*(cosd(20)2); Cytqzhi=4/Rm2*(r1*sind(alpha)*cosd(alpha)*sind(20)* cosd(20)*L2*(
