Matlab平板（2D）空中运动轨迹模拟代码

%--------------------------------------------------------% % Trajectory of 2D plate based on Holmes's model %--------------------------------------------------------% close all clear clc % properties of plate m=0.196; % mass A=0.133; % area l=sqrt(A); % chord/span I=0.001; % inertia moment % properites of air rou=1.225; % density of air % time step T=10; % Total time dt=0.01; % time step nt=T/dt; % number of time steps t=(0:dt:T-dt); % constant g=9.81; % initialization ax=zeros(nt,1); % horizontal acceleration ay=zeros(nt,1); % vertical acceleration aw=zeros(nt,1); % angular acceleration u=zeros(nt,1); % horizontal velocity v=zeros(nt,1); % vertical velocity w=zeros(nt,1); % angular velocity V=zeros(nt,1); % total velocity x=zeros(nt,1); % horizontal displacement y=zeros(nt,1); % veritical displancement theta=zeros(nt,1); % the angle of the plate to the horizontal, rad % initial condition U=35; % wind speed or relative horizontal velocity v(1)=0; theta(1)=deg2rad(0.0); y(1)=75; % calculation for i=1:nt-1 V(i)=sqrt((U-u(i))^2+v(i)^2); % square of total velocity beta=atan(-v(i)/abs(U-u(i))); % rad alpha=rad2deg(theta(i)+beta); % angle of attack, deg, [0,180] alpha=mod(alpha,180); ax(i)=sign(U-u(i))*0.5*rou*(Cd(alpha)*cos(beta)-(Cl(alpha)+C_Lr(w(i),U,l))*sin(beta))*A*V(i)^2/m; ay(i)=0.5*rou*(Cd(alpha)*sin(beta)+(Cl(alpha)+C_Lr(w(i),U,l))*cos(beta))*A*V(i)^2/m-g; aw(i)=0.5*rou*Cm(alpha)*A*l*V(i)^2/I; u(i+1)=u(i)+ax(i)*dt; v(i+1)=v(i)+ay(i)*dt; w(i+1)=w(i)+aw(i)*dt; x(i+1)=x(i)+(U-u(i))*dt; y(i+1)=y(i)+v(i)*dt; theta(i+1)=theta(i)+w(i)*dt; % rad if (y(i+1)<=0) disp('plate dropped onto the ground'); iEnd=i; break end end % plot figure plot(x,y,'.') grid on title('x-y') figure plot(t,x,'.') grid on title('t-x') figure plot(t,y,'.') grid on title('t-y') figure; plot(t,theta,'.') grid on title('t-theta') figure; plot(t,V,'.') grid on title('t-V')