高超声速气动力及导数计算报告

高超声速气动力及导数计算报告撰写人：学号：_班级:2012年四月翌日引言实验目的: 熟练运用面元法中有关网格划分的方法； 掌握高超音速气动力工程估算方法中面元法及牛顿法；同时，比较两者的计算 结果，并分析差异产生原因。实验条件：编程计算该旋成体的升力、阻力、升阻比及俯仰力矩系数，还有导数作出曲计算条件：a = 5。30。,M尸 4.510.5, H = 10km,。= y = 0几何参数：a = 0.1m,b = 0.5m,D = 3.912m飞行器运动状态：P = Q = R = 0二计算方法规定导弹的体坐标系为：x轴沿着导弹纵轴向后，丫轴垂直于弹体向上，Z 轴于其它两轴构成右手坐标系，即指向导弹左侧，原点位于导弹前缘点。取原点 为参考点。牛顿法:牛顿法的假设如下： 攻角a小于物面倾角。； 假设流体由大量均匀分布的，彼此独立无相互作用的质点所组成，它们排 列整齐、平行地沿着直迹线流向物体。 流体质点流与物面碰撞时，流体质点将失去与物面垂直的法向动量，而保 持原有的切向动量沿物面向下流下去。由于法向动量的变化从而引起流体作用在 物体上的力。 流体对物面的压力只作用在物面能与流体质点相碰撞的表面上，而遮蔽区 上压力为零。牛顿公式：C = 2sin2 9其中9为来流速度方向与物面切面的夹角。由课本可知：21 .八、，、C =j 欠2(cos2 以 sin2 9 + -sin2 以 cos2 9)(兀 一丫 )max+4cos 以 sin以 cos9 sin9 siny - sin2 以 cos2 9 cosy siny rdr八-2C =j 技2cosa sina cos9 sin9 cosy sin y - 2sin y (cos2 a sm2 9max2+sin2 a cos2 9) 一 2cosa sina cos9 sin9 (兀一y ) + sin2 a cos2 9 sm3 y rdx对于圆球部分，由于用于验证，我们假设a。，y0 :/八、r2r2(C ). = cos2 a + isin2 amaxmaxr2(C ) . = cos a sin amax对于圆锥部分，由于每一部分圆锥各自的9、y广a一直为常值，因此直接带入积分 进行计算，仅与yk、a有关，当a 9时，y& = arccos(9/a).其余为零。由于用于验证， 我们假设a k k y kZi = Zi + nd k k z k -式中，dk=气 G -X，)+ n Y - y)+ n G-Z )k yk zkk = 1,2,3,42）建立面元坐标系在面元形成后，需要建立一个面元坐标系。现在我们选取面元法向单位矢量n作为一个单位矢量，再把T的单位矢量1作为一个单位矢量： 11t = t i +1 j +1 k11x 1y 1zt =、, t = Ty , t = T1x T 1 y T 1z TT =（T 2 + T 2 + T 2 1 1X 1Y 1Z由n和t的矢量积可以决定第三的单位矢量m :1m = n x t = m i + m j + m km = n t 一 n tx y 1z z 1 ym = n t n t y z 1x x 1zm = n t n tz x 1y y 1x至此，已经建立了面元坐标系，原点o取G,Yi,万）x轴平行于t，PP1轴平行于m， Z轴平行于n。易知Xp轴和Y轴在面元平面上。3）面元坐标系与飞行器参考坐标系之间的变换在具体计算中，需要把飞行器坐标系中的坐标值变换到面元坐标系中去，也 需要把面元系中的值转换到飞行器坐标系中去。此时，转换公式为： 飞行器坐标转换到面元坐标X tttX XP1x1 y1zY=mmmY YPxyzZnnnZ ZP1-XyzL从面元坐标系转换到飞行器坐标系X XtmnX 1xXxPY Y=tmnY1yyyJzZLtmnZL1zzzP4）面元的面积和质心坐标现在把投影点Xi, Y ,Zi ）在飞行器坐标系中的值变换到面元坐标系中去， k k k在面元坐标系中的值用Xi ,Yi ,Zi来表示，则pk pk pkXi= tpk1xQ一x）+t YY）+t G一z） _）1 y（ ）1z（_）yi = mpk xXi -x4m Y4m Vi -ZAky kz kzi = 0k = 1,2,3,4在面元平面上，由凸四边形面积公式可求得面元的面积为AAAA =这里由于r平行于，轴，因此，面积公式可以化简为AA 。2 p 3在面元坐标系中，面元质心的坐标为xi (yip 4 pl-xi )(yi - yi )i 1Xi =p03 yi - yi1p pyip o = - 3 ypi不0 = 0-yi )+ xi (yi - yip 2p 2 p 4pl利用面元坐标系转换到飞行器坐标系的转换关系，可求得质心在飞行器坐标中的 质心坐标x - X1 +1 xi + m yiyi Yi +1 xi + m yi 01 y p oy p 0z Zl +1 x + m yi01z p 0 z p 0 J5) 面元冲击角在气动力的计算中，我们必须要知道平面面元与自由流速度的夹角，即冲击 角。而第i个面元平面的冲击角：0 i - COS-1 童以2fnv式中n为第i面元平面的单位法向向量。速度矢量V=v -Qx r式中，Q总角速度Q - Pi - Qj - Rkr向径：r = 3 x )了 + (y y ) j + (z z ) j0 r0 r0 rV为来流速度，P, Q, R为物体滚转、偏航、俯仰角速度，x , y , z为参考点*r r r的坐标，在编程过程中取(0，0，0)。6) 计算面元压力系数根据冲击角是否大于零，可以把面元所处的流场分为迎风面和背风面。若冲击角大于零，则可利用修正牛顿理论：C = C o2 sin20，计算面元的压 力系数。其中0为来流速度与面元内法线的夹角Cp 02y M 2a(丫+1叫 E 丫+11122 M 2 (y 1)a若冲击角小于零，则可直接令Cp = o。7) 面元上气动力系数及动静导数计算将飞行器表面分成若干面元后，用积分的方法可以得到它的气动力系数及其 静、动导数。参考坐标系中(X轴沿着导弹纵轴向后，Y轴垂直于弹体向上，Z轴于其它 两轴构成右手坐标系，即指向导弹左侧，原点位于导弹前缘点。)，来流速度矢量 V* = V cosacosP cosPsina sin。，假设面元外法向矢量 n = n n n ，则来 流速度在面元内法向上的投影为：V V =ncosacos。一 n cosP sin a n sin。VJV分别对a和P求导，得(V V ) = n sinacos。 n cos。cosa上* axy(V V ) = n cosasin。+ n sin。sina n cos。 * Pxyz物体滚转、俯仰、偏航角速度为p、q、r，由p、q、r三个角速度引起的物体与气流相对运动的三个速度为:yp =Oi-p乙j + pyjVq -qz-i + O j -qx-kVr -ry-i + rx- j+ 0-其中、y、z表示面元质心到参考点的距离。令VVVr在内法向的投影为( /V加(V /V用(V /V y ,若记 1_00001_00V =(Vp+Vg+Vr)M p q、的导数为动 111(V /V ) =(n z-n y)/D动上 00 pyz f(V /V ) =(n x-n Z)/D动-L( oo q z x(V /V ) =(n y-n x)/D动_1_ oo r x y 式中D参考面积。总速度在面元法向的投影：V =U +U ,则修正的牛顿理论亦可表示为Lfull 1 动上(y YC -C sin2 0 r C 顷 oP P02冲V )00根据上述计算，压力系数的静、动导数计算公式为：/vw 虹. (VG云&动J动 7 7 7 7 7 /V8/V8 也/匕 叩叩叩/(v/OK)/oo/ 800;尹尹/根据前面计算的面元上的几何参数及压力系数和静、动导数，可用下列公式计算面元上气动力、力矩及其静、动导数。面元上所受的气动力：dF = -C -n -dSXp XdF = -C -n -dSyp ydF = -C -n -dSzp z面元上气动力对参考点的力矩:dm = dF - y - dF - z dm = dF - z - dF - x dm = dF - x - dF - y面元上气动力的静导数:dF = -C - n - dSXa.a XdF =-C - n - dSdF =-c - n - dSyaa ydF =-C - n - dSdF =-cL- dS面元上气动力产生的力矩的静导数:dm = dF - y - dF - z dm = dF - z - dF - x dm = dF - x - dF - y面元上气动力的动导数：zayaXadF =-C - n - dSdF =-C - n - dSXpp p xddSdF =-C - n - dSyrPr ydF =-C - n - dSxp xdF =-C - n - dSZqp z面元上力矩的动导数：dm = dF - y - dF - z，rrdm = dF z - dF - xq%dm = dF x - dF yZppXp求出每个面元上的气动力、力矩及其静动导数以后，通过求和或积分就可求出整体气动力系数。轴向、法向和侧向力系数：C =(1/S )L(dF )C =(1S S(dF)轴向、法向和侧向力的导数：/ 、C =(1/S) J) CyM:欤：俯仰、偏航和滚转力矩系数:CxCyCz= 1/（D S）（dm ） = 1/（D S）（dm ） = 1/（D S ）Z（dm ）俯仰、偏航和滚转阻尼导数：CxrC七Cz p= 1/（Q S）0m ） = 1/（D S）冬） = 1/（D S ）X（dmq）其中体轴系与风轴系之间的转换矩阵:X fYfZL f=cos a cos Psin a一 cos a sin P sin a一 sin a cos Pcosacos Psin P0cos PXtYtZt计算程序程序：%注意：该程序中，语句后的%为解释用，语句前的%是为加快速 度没必要计算的滚转方向%由于画图未使用函数调用而采用循环速度较慢大概在20秒左 右% 面元法%clear allclcsyms alf blt grm p q 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)*42*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:3for jjj=1:71Ma=4.5+3*(iii-1); %Ma 在 4.5 到 10.5 之间画图alf=(-5+05*(jjj-1)*pi/180;%攻角在-5到30度之间画图xx=zeros(m,4*n);yy=zeros(m,4*n);zz=zeros(m,4*n);for i=1:mxx(i,1:n)=01/n*(1:n);xx(i,n+1:2*n)=01+(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*42*D/4+42*D/4/5*(1:n);end%确定每个点的x坐标for j=1:4*nif xx(1,j)01zz(1,j)=sqrt(13”2-(13-xx(1,j)”2);elseif xx(1,j)4.2*D/4zz(1,j)=05+(xx(1,j)-01)*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)+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)+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),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);elseT1=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/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(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);endxb(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,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及其单位矢量nm1(i,:)=cross(n1(i,:),t1(i,:);endxyzb=xb yb zb;cha=cell(4*m*n,1);for i=1:4*n*mfor ii=1:4chai,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*mfor j=1:4dk(i,j)=dot(n1(i,:),chai,1(j,:);endend%求dkxyz=cell(4*m*n,1);for i=1:4*n*mfor ii=1:3for j=1:4xyzi,1(j,ii)=xyzki,1(j,ii)+n1(i,ii)*dk(i,j);endendend%求面元四个角点投影到面元平面上的坐标xyzxyzpk=cell(4*m*n,1);for i=1:4*n*mfor j=1:4xyzpki,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%求面元坐标系中四个角点坐标xyzpkA=zeros(1,4*m*n);for i=1:4*n*mA(1,i)=05*(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*nxyzp0(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%求质心坐标xyzp0xyz0=zeros(4*m*n,3);for i=1:4*m*nxyz0(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%质心在飞行器坐标系中的坐标xyz0v=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*nvcv(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)*cos(blt);%计算v*n/v来流对攻角的导数 vcvb(i,1)=n1(i,1)*cos(alf)*sin(blt)+n1(i,2)*sin(alf)*sin(b%计算v*n/v来流对侧滑角的导数lt)-n1(i,3)*cos(blt);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=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;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的导数注意一下符号不太确定endvfull=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-(brb-1)(1/(brb-1)-1); %计算 cp02 cp=zeros(4*m*n,1);for i=1:4*m*nif vcv(i,1)0cp(i,1)=cp02*vfull(i,1)”2;elsecp(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*n dFc(i,1)=-cp(i,1)*n1(i,3)*A(1,i);%计算每个平面面元侧向力注意此处没乘动压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(i,2)*dFz(i,1);%计算每个平面面元俯仰力矩 注意此处没乘动压enddongya=05*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*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)/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*nCpa(i,1)=2*cp02*vcv(i,1)*vcva(i,1);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 ,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)*A(1,i) *xyz0(i,2)+Cpp(i,1)*n1(i,2)*A(1,i)*xyz0(i,3);% 以上是力矩动导数end cma=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(iii,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(a)；ylabel(Y);title(升力Y与迎角a、马赫数Ma关系曲线)； legend(Ma=45,Ma=75,Ma=105);%升力Y与迎角a、马赫数Ma关系曲线 figure(3);plot(alpha,czl(1,:),g.-);grid on;hold on;plot(alpha,czl(2,:),b*-);plot(alpha,czl(3,:),r+-);xlabel( a);ylabel(X);title(阻力X与迎角a、马赫数Ma关系曲线)；%阻力X与迎角a、马赫数Ma关系曲线 legend(Ma=4.5,Ma=75,Ma=10.5);figure(4);plot(alpha,cld(1,:),g.-);grid on;hold on;plot(alpha,cld(2,:),b*-);plot(alpha,cld(3,:),r+-);xlabel(a);ylabel(Y/X);title(升阻比Y/X与迎角a、马赫数Ma关系曲线)； legend(Ma=4.5,Ma=75,Ma=10.5);%升阻比Y/X与迎角a、马赫数Ma关系曲线 figure(5);plot(alpha,cfylj(1,:),g-);grid on;hold on;plot(alpha,cfylj(2,:),b*-);plot(alpha,cfylj(3,:),r+-);xlabel(a);ylabel(mz);title(俯仰力矩系数mz与迎角a、马赫数Ma关系曲线)； legend(Ma=4.5,Ma=75,Ma=105);%俯仰力矩系数mz与迎角a、马赫数Ma关系曲线figure(6);plot(alpha,cphjnlj(1,:),g.-);grid on;hold on;plot(alpha,cphjnlj(2,:),b*-);plot(alpha,cphjnlj(3,:),r+-);xlabel( a );ylabel(偏航 阻尼导数)；title(偏航阻尼导数与迎角a、马赫数Ma关系曲线)； legend(Ma=45,Ma=75,Ma=105);%偏航阻尼导数与迎角a、马赫数Ma关系曲线 figure(7);plot(alpha,cfyznlj(1,:),g.-);grid on;hold on;plot(alpha,cfyznlj(2,:),b*-);plot(alpha,cfyznlj(3,:),r+-);xlabel( a );ylabel(俯仰 阻尼导数)；title(俯仰阻尼导数与迎角a、马赫数Ma关系曲线)； legend(Ma=4.5,Ma=75,Ma=10.5);%俯仰阻尼导数与迎角a、马赫数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*D;L1=0.1;L2=L/4-L1;L3=05*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);% 前 锥 段 部 分 ( Y k=0)Cxtqzhi=1/Rm2*(r2”2-r1”2)*(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*(1-tan(20*pi/180)+0.5*sind(alpha)*cosd(20)*sind(20)*cosd(alpha)*tan(20*pi/180)*(L1+L2)2-L2);%柱段部分Cxtzhu=0;Cytzhu=8/(3*pi)*(sind(alpha)2*L3*r2/Rm2;%后锥段部分Cxthzhi=1/Rm2*(Rm2-r22)*(2*(cosd(alpha)2*(sind(165) 2 +(sind(alpha)“2*(cosd(16.5)2);Cythzhi=4/Rm2*(r2*sind(alpha)*cosd(alpha)*sind(165)* cosd(16.5)*L4*(1-tan(16.5*pi/180)+0.5*sind(alpha)* cosd(16.5)*sind(165)*cosd(alpha)*tan(165*pi/180)*(L2-(0 .75*L)2);Cxt=Cxtq+Cxtqzhi+Cxtzhu+Cxthzhi;Cyt=Cytq+Cytqzhi+Cytzhu+Cythzhi;Cx(h)=Cyt*sind(alpha)+Cxt*cosd(alpha);Cy(h)=Cyt*cosd(alpha)-Cxt*sind(alpha);endS=pi*Rm2;%参考面积dongya=0.413*(Ma*340)2/2;%动压Yn=Cy.*(S*dongya);%升力Xn=Cx.*(S*dongya);%阻力alpha=-5:dh:11;%画出升力随攻角变化曲线figure(8);plot(alpha,Yn,r);hold on;grid on;xlabel(a)；ylabel(Y);title(Ma=7.5时由牛顿法与面元法得到升力Y随攻角变化曲线)；%画出阻力随攻角变化曲线figure(9);plot(alpha,Xn,r);hold on;grid on;xlabel(a);ylabel(X);title(Ma=7.5时由牛顿法与面元法得到阻力X随攻角变化曲线)；%画出升阻比系数随攻角变化曲线figure(10);plot(alpha,Yn./Xn,r);hold on;grid on;xlabel(a);ylabel(Y/D);title(Ma=7.5时由牛顿法与面元法得到升阻比Y/D随攻角变化曲线);alpha=-5:0.5:11;figure(8);plot(alpha,csl(2,1:33),g.-);hold on;legend(牛顿法，面元法)；figure(9);plot(alpha,czl(2,1:33),g.-);hold on;legend(牛顿法，面元法)；figure(10);plot(alpha,cld(2,1:33),g.-);hold on;legend(牛顿法，面元法)；四计算结果1）面元法计算结果20图（1）图（2）图（4）偏航阻尼导数与迎角口、马赫数&关系曲线3.1iiiI：* Ma4.5M 或 51 Ma=10.6051015202S30nn图（6）俯仰阻尼导数与迎甬口、马新数关系曲线0-掠呻回昌与毋图（7）2）牛顿法与面元法计算结果对比取马赫数Ma=7.5,分别计算攻角从-5度变化到11度过程中升力与阻力的变化情况。攻角在这个区间内时，旋成体只有柱段存在被风区，计算相对较简单。图（8）Ma=7.5时由牛顿，去与面元法得到汗阻比HD随攻角变化曲线图（10）由图（10）、（11）、（12）可以看出，用牛顿法与面元法算得的阻力较为近似， 但升力差别非常大。5.结果分析及结论 从面元法结果图图中来看在攻角为-530变化区间内，该旋成体升力与攻 角接近于线性正相关变化，俯仰力矩与攻角接近于线性正相关变化，虽然阻力 与马赫数、迎角成正相关，但与升力不同的是，阻力随攻角变化率并没有减小， 而是越来越大，由此可见高超声速飞行器的减阻设计将会有很大的意义。 从升阻比曲线来看，当攻角约为15时此时升阻比最大，此位置接近于物体 将要出现Yk的位置，说明升阻比的极限值与Yk有关。 俯仰力矩系数、偏航阻尼导数、俯仰阻尼导数随马赫数的变化几乎没有改变； 而随着攻角的增大，三者的绝对值均有所增加，这是由于攻角在一定范围内变 化时，气动力增加的缘故。 从升力系数以及阻力系数的变化情况来看，使用牛顿法直接积分和使用网格法 来计算二者变化趋势一致，但是具体数值之间存在着一定的误差，通过对轴向 力和法向力大小的对比可知，误差原因主要出现在轴向力上，考虑到只有最前 端圆球部分主要影响轴向力对法向力影响较小，误差原因可能是由于圆球部分 划分网格较少，通过对此部分网格增加分析，发现误差变化不大，因此具体的 误差原因可能是由于模型之间结果的误差造