面元法计算机身的气动参数的matlab程序代码

面元法计算超音速机身的气动参数的编程实验机身参数如下图所示：面元法计算机身的气动参数的matlab程序代码（仅供参考）function[]=surface(a,b,D,n)P=0;Q=0;R=0;L=D*4.2;A=16.5*pi/180;B=20*pi/180;S=pi*(L/4*tan(A)+D/2)^2fori=1:11X(i)=-a*(i-1)/10;r1=sqrt(1-(X(i)+a)*(X(i)+a)/(a*a))*b;forj=1:(n+1)Y(i,j)=r1*cos(2*pi*(j-1)/n);Z(i,j)=r1*sin(2*pi*(j-1)/n);endendfori=12:21X(i)=-a-(i-11)*(L/4-a)/10;r2=-(X(i)+a)*tan(B)+b;forj=1:(n+1)Y(i,j)=r2*cos(2*pi*(j-1)/n);Z(i,j)=r2*sin(2*pi*(j-1)/n);endendfori=22:31X(i)=-L/4-(i-21)*(L/2)/10;r3=D/2;forj=1:(n+1)Y(i,j)=r3*cos(2*pi*(j-1)/n);Z(i,j)=r3*sin(2*pi*(j-1)/n);endendfori=32:41X(i)=-L*3/4-(i-31)*(L/4)/10;r4=(-X(i)-L*3/4)*tan(A)+D/2;forj=1:(n+1)Y(i,j)=r4*cos(2*pi*(j-1)/n);Z(i,j)=r4*sin(2*pi*(j-1)/n);endendfori=1:40forj=1:nk=(i-1)*n+j;p1(k,:)=[X(i)Y(i,j)Z(i,j)];p2(k,:)=[X(i)Y(i,(j+1))Z(i,(j+1))];p3(k,:)=[X(i+1)Y((i+1),(j+1))Z((i+1),(j+1))];p4(k,:)=[X(i+1)Y((i+1),j)Z((i+1),j)];T1(k,:)=p3(k,:)-p1(k,:);t0(k)=sqrt(dot(T1(k,:),T1(k,:)));t(k,:)=T1(k,:)/t0(k);T2(k,:)=p4(k,:)-p2(k,:);N(k,:)=cross(T2(k,:),T1(k,:));n0(k)=sqrt(dot(N(k,:),N(k,:)));n1(k,:)=N(k,:)/n0(k);s0(k,:)=(p1(k,:)+p2(k,:)+p3(k,:)+p4(k,:))/4;d1(k)=dot(n1(k,:),(s0(k,:)-p1(k,:)));d2(k)=dot(n1(k,:),(s0(k,:)-p2(k,:)));d3(k)=dot(n1(k,:),(s0(k,:)-p3(k,:)));d4(k)=dot(n1(k,:),(s0(k,:)-p4(k,:)));t1(k,:)=p1(k,:)+n1(k,:)*d1(k);t2(k,:)=p2(k,:)+n1(k,:)*d2(k);t3(k,:)=p3(k,:)+n1(k,:)*d3(k);t4(k,:)=p4(k,:)+n1(k,:)*d4(k);m(k,:)=cross(n1(k,:),t(k,:));x1(k)=dot(t(k,:),(t1(k,:)-s0(k,:)));y1(k)=dot(m(k,:),(t1(k,:)-s0(k,:)));z1(k)=0;x2(k)=dot(t(k,:),(t2(k,:)-s0(k,:)));y2(k)=dot(m(k,:),(t2(k,:)-s0(k,:)));z2(k)=0;x3(k)=dot(t(k,:),(t3(k,:)-s0(k,:)));y3(k)=dot(m(k,:),(t3(k,:)-s0(k,:)));z3(k)=0;x4(k)=dot(t(k,:),(t4(k,:)-s0(k,:)));y4(k)=dot(m(k,:),(t4(k,:)-s0(k,:)));z4(k)=0;s(k)=(x3(k)-x1(k))*(y2(k)-y4(k))/2;m0(k)=(x4(k)*(y1(k)-y2(k))+x2(k)*(y4(k)-y1(k)))/(3*(y2(k)-y4(k)));m1(k)=-y1(k)/3;m2(k)=0;p(k,:)=s0(k,:)+(t(k,:)*m0(k))+(m(k,:)*m1(k));endendforq=-5:30i=q+6;V=[-cos(q*pi/180)0sin(q*pi/180)];cx=0;cy=0;mz=0;fork=1:(40*n)r(k)=pi/2-acos(-dot(n1(k,:),V)/norm(n1(k,:)));d=2.38+0.03792*r(k)-0.002521*r(k)^2+0.00004583*r(k)^3+2.917*10^(-7)*r(k)^4;ifr(k)>=0cp(k)=d*sin(r(k))*sin(r(k));elsecp(k)=0;endcx=cx+cp(k)*n1(k,1)*s(k)/S;cy=cy-cp(k)*n1(k,3)*s(k)/S;mz=mz+(-n1(k,3)*p(k,1)+n1(k,1)*p(k,3))*cp(k)*s(k)/(L*S);ende=cx;f=cy;Cy(i)=-cx*sin(q*pi/180)+cy*cos(q*pi/180);Cx(i)=cx*cos(q*pi/180)+cy*sin(q*pi/180);Mz(i)=mz;szb(i)=Cy(i)/Cx(i);endfori=1:36;x(i)=-6+i;endsubplot(2,3,1);plot(x,Cy,'r.-');xlabel('α');ylabel('Cy');holdon;subplot(2,3,2);plot(x,Cx,'b.-');xlabel('α');ylabel('Cz');holdon;subplot(2,3,3),plot(x,Mz,'c.-');xlabel('α');ylabel('Mz');holdon;subplot(2,3,4);plot(x,szb,'g.-');xlabel('α’);ylabel('L/D');holdon;fori=1:36ma(i)=4.5+(i-1)*6/35;endforq=-5:5:30i=q+6;forj=1:36Y0(j)=Cy(i)*(ma(j)*299.46)^2*0.4127/2*S;Z0(j)=Cx(i)*(ma(j)*299.46)^2*0.4127/2*S;endsubplot(2,3,5);plot(ma,Y0,'m.-');xlabel('Ma');ylabel('Y');text(10,(-2