flapping，UVLM，matlab code

使用非定常涡格法(UVLM)计算一种简单扑翼机的升力。

%%%%%%%%%%%%%% FORCES.M %%%%%%%%%%%%%%%%%%%%
%% This is a matlab function that will compute the forces experienced by the airfoil 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
start=cputime;
V=1;
s=（0:ds:sf）‘;
% wa （Eq. 2.7）;
wav=-sin（alpha） - zdot.*cos（alpha）/V;
% wadot （Eq. 2.8）
wadotv=-alphadot.*cos（alpha） -...
（zddot.*cos （alpha）-zdot.*sin（alpha）.*alphadot） /V;
% CLs （Eq. 2.14）
CLs=-2*pi*wav;
% gammaw （Eq. 2.24）
gammaw=2*pi*sqrt（s./（2+s））.*wadotv;
% CFw （Eq. 2.27）
X=1./（s+2）;
for n=l:sf/ds+l
CFw（n）=2*pi*wadotv（l:n）‘*X（n:-1:1）*ds;
end
CFw=CFw‘;
% CFi （Eq. 2.21）
CFi=pi*wadotv;
% CPx （Eq. 2.41）
CPx=2*pi*（wav-CFw/（2*pi））.^2;
% CL （Eq. 2.2）
CL=CLs+CFw+CFi;
% CT （Eq. 2.2）
CT= CPx-（CLs+CFw+CFi）.*sin（alpha）;
CTperiod= CT（round（length（s）-2*T/ ds）:length （s））;
CLperiod= CL（round （length （s）-2*T/ ds）:length（s））;
period= [s（length（s））-2*T:ds:s（length（s））];
ctint= 1/2/T*trapz（period CTperiod）;
clint=1/2/T*trapz（period CLperiod）;
stop=cputime;
run=stop-start;
info=[‘CT = ‘ num2str（ctint） ‘; sf = ‘ num2str（sf） ‘; ds = ‘ num2str（ds） ‘;
cputime = ‘ num2str（run） ‘ seconds‘]


%%%%%%%%%%%% % TRAILING.M %%%%%%%%%%%/%%%%%/%%
%% This is a matlab function that will plot the airfoil motion in a movie style %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
start=cputime;
load naca009
s=0:ds:sf;
x= naca009（:1）/50-1;
xtr=cos（alpha）;
ztr= sin （-alpha） +z;
xl=-cos（alpha）;
zl= sin（alpha）+z;
zu=naca009（:2）/50;
xc=x*cos（alpha‘）;
for n=l:length（s）
for m=1:18
zc （mn）= z （n）-xc（m）*sin（alpha（n））;
end
end
M=moviein（21）;
for n=13:10:143
N=（n-12+9）/10;
hold off
plot（xc（:n） zc（:n）+zu）
hold on
plot（xc（:n） zc（:n）-zu）
axis（[-2 2 -2 2]）; grid
set（gca ‘Visible‘ ‘on‘）
M（:N）=getframe;
end
movie（M）;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       1767  2013-03-21 00:25  flapping.m

----------- ---------  ---------- -----  ----

                 1767                    1