飞行器机翼非定常气动力计算及颤振计算程序.7z

clear clf s=7.5; c=2; m=100; kappa_freq=5; theta_freq=10; xcm=0.5*c; xf=0.48*c; e=xf/c-0.25; velstart=1; velend=180; velinc=0.1; a=2*pi; rho=1.225 Mthetadot=-1.2; M=(m*c^2-2*m*c*xcm)/(2*xcm); damping_Y_N=1; if damping_Y_N==1; z1=0.0; z2=0.0; w1=2*2*pi; w2=14*2*pi; alpha=2*w1*w2*(-z2*w1+z1*w2)/(w1*w1*w2*w2); beta=2*(z2*w2-z1*w1)/(w2*w2-w1*w1); end a11=(m*s^3*c)/3+M*s^3/3; a22=m*s*(c^3/3-c*c*xf+xf*xf*c)+M*(xf^2*s); a12=m*s*s/2*(c*c/2-c*xf)-M*xf*s^2/s; a21=a12; A=[a11,a12;a21,a22]; k1=(kappa_freq*pi*2)^2*a11; k2=(kappa_freq*pi*2)^2*a22; E=[k1 0;0 k2]; icount=0; for V=velstart:velinc:velend icount=icount+1; if damping_Y_N==0; C=[0,0;0,0]; else C=rho*V*[c*s^3*a/6,0;-c^2*s^2*e*a/4,-c^3*s*Mthetadot/8]+alpha*A+beta*E; end K=(rho*V^2*[0,c*s^2*a/4;0,-c^2*s*e*a/2])+[k1,0;0,k2]; Mat=[[0,0;0,0],eye(2);-A\K,-A\C]; lambda=eig(Mat); for jj=1:4 im(jj)=imag(lambda(jj)); re(jj)=real(lambda(jj)); freq(jj,icount)=sqrt(re(jj)^2+im(jj)^2); damp(jj,icount)=-100*re(jj)/freq(jj,icount); freq(jj,icount)=freq(jj,icount)/(2*pi); end Vel(icount)=V; end figure(1) subplot(2,1,1); plot(Vel,freq,'k'); vaxis=axis;xlim=([0 vaxis(2)]); xlabel('Air speed(m/s)'); ylabel('Freq(Hz)'); grid subplot(2,1,2); plot(Vel,damp,'k'); xlim=([0 vaxis(2)]); axis([xlim ylim]); xlabel('Air speed(m/s)'); ylabel('Damping Ratio(%)'); grid