%Berman and Wang/Boledner ; Plotting F Vs t/T
4/29/12
%%
c=4.02; % in mm mean chord length
R= 13.2; % in mm max chord length
r=(0:3.3:R);%radius
C = ((4*c)/pi)*sqrt(1-((r.^2)/(R)^2)); % C= c(r)the length of chord of wing
M=mean(C)
figure(1);
plot(r,C);
xlabel('r in (mm)')
ylabel('C(r) (mm)')
title('Chord Length as a Function of Radius')
%%
rhof= 1.29; % kg/m^3 density of surrounding fluid
M22 = (pi/4)*rhof*C.^2;
%%
thetam= .21468 %12.3 degrees
N=1
f=164 %frequency
thetadeg = 1.783726 %-102.2 degree
thetanot= 0.0319395 %1.83;
t=[0:.2:1]
theta=thetam*cos(2*pi*N*f*t+thetadeg)+thetanot;
figure(2);
plot(t,theta)
xlabel('Time in (Periods)')
ylabel('Theta(t) in (degrees)')
title('Optimized theta')
%%
etam= 1.518436 % 87; %deg
Ceta= 0.1
phieta= -pi/2 % -90; %degrees
etanot= pi% 180; %degrees
eta = (etam/tanh(Ceta))*tanh(Ceta*sin(2*pi*f*t+ phieta))+ etanot;
figure(3);
plot(t,eta)
xlabel('Time in (Periods)')
ylabel('Eta(t) in (radians)')
title('Optimized Eta')
%%
K=.1;
phim= 0.214675 %12.3
phi= (phim/(sin(K).^-1))*sin(K*sin(2*pi*f*t)).^(-1);
figure(4);
plot(t,phi)
xlabel('Time in (Periods)')
ylabel('Phi(t) in (degrees)')
title('Optimized Phi')
%%
theta=thetam*cos(2*(pi)*N*f*t+thetadeg)+thetanot;
thetadot=diff(theta);
%%
phi=(phim/(sin(K).^-1))*sin(K*sin(2*pi*f*t)).^(-1);
phidot =diff(phi); %just in foward flight
%%
eta = (etam/tanh(Ceta))*tanh(Ceta*sin(2*pi*t.*f + phieta))+ etanot;
etadot = diff(eta); %just in foward flight
%%
%Vyp = r.*(thetadot.*cos(eta)-phidot.*cos(theta)*sin(eta))
Vyp= [0 .1 .2 .3 .4 ];%note this numbers are fake and used just for programming purposes
Vxp= [0 .6 .7 .8 .9];%note this numbers are fake and used just for programming purposes
V=sqrt((Vxp).^2+(Vyp).^2);
%%
alpha= .4886921 %28 degrees
Crot= pi;
CT=1.341;
W= 0.9774 ;%2*alpha
Q= [0 3 .01463 .22702 .8926] %c.*V orgiginal was [0 3 .01463 .22702 .8926 0]
U= -0.6705 %-.5*CT
CC= [26.1983 24.5609 19.6487 11.4617 0] %C^.2
SS= 0.8290 %sin(W)
%Gamma= -0.5*CT*C.*V.*sin(2*alpha)+.5*Crot*C.^2.*etadot;
Gamma= U*Q*SS+.5*Crot*CC.*etadot
%%
b=5.12;
M11 = (pi/4)*rhof*b^2; % added mass term
%%
%axp= r.*((phiddot*cos(theta)+thetadot*(etadot-phidot*sin(theta))*cos(eta)+...
% (thetaddot-etadot*phidot*cos(theta))*sin(eta)));
axp=[ 0 .1 .2 .3 .4]; %note this numbers are fake and used just for programming purposes
%%
CD_0= 0;
CD_pi = 2.93;
dFvxp= 0.5*rhof*C.*CD_0*(cos(alpha))^2+CD_pi*(sin(alpha))^2*V.*(Vxp);
%%
Mwing= 0.4;
D=[ 0.0386 0.0374 0.0334 0.0255 0] % C./(c*R)*Mwing
O=[0 0.1700 -0.0000 -0.5101 -0.4191] %Vyp.*etadot
Z=[0 8.2479 -0.0021 -11.8971 -0.2560] %rhof*Gamma.*Vyp
E=[ 0 2.6559 5.3119 7.9678 10.6238] %M11*axp
%dFxp=((((C./(c*R))*Mwing+M22)*Vyp.*etadot-rhof*Gamma.*Vyp-M11*axp))dr-dFvxp;
((D+M22)*O-Z-E) -dFvxp
4/23/12
%%
c=4.02; % in mm
R= 13.2; % in mm
r=(0:.131:R);
C = ((4*c)/pi)*sqrt(1-((r.^2)/(R)^2)); % C= c(r)
M=mean(C)
figure(1);
plot(r,C);
xlabel('r in (mm)')
ylabel('C(r) (mm)')
title('Chord Length as a Function of Radius')
%%
rhof= 1.29; % kg/m^3 density of surrounding fluid
M22 = (pi/4)*rhof*C.^2;
%%
thetam=12.3;
N=1
f=164
thetadeg= -102.2
thetanot= 1.83;
t=[0:.2:1]
theta=thetam*cos(2*(pi)*N*f*t+thetadeg)+thetanot;
figure(2);
plot(t,theta)
xlabel('Time in (Periods)')
ylabel('Theta(t) in (degrees)')
title('Optimized theta')
%%
etam=87; %deg
Ceta=0.1;
phieta=-90; %degrees
etanot=180; %degrees
eta = (etam/tanh(Ceta))*tanh(Ceta*sin(2*pi*f*t+ phieta))+ etanot;
figure(3);
plot(t,eta)
xlabel('Time in (Periods)')
ylabel('Eta(t) in (degrees)')
title('Optimized Eta')
%%
K=.1;
phim=12.3;
phi= (phim/(sin(K).^-1))*sin(K*sin(2*pi*f*t)).^(-1);
figure(4);
plot(t,phi)
xlabel('Time in (Periods)')
ylabel('Phi(t) in (degrees)')
title('Optimized Phi')
____________________________________________________________________
%Parameters
c=4.02; % in mm
R= 13.2; % in mm
r=(0:.131:R);
b=5.12; %in mm c(r=0)= (4c/pi)*sqrt(1-(r^2/R^2))
t=[0:.1:10];
Mwing= 0.4; % BB01 mg
rhof= 1.29; % kg/m^3 density of surrounding fluid
N=1;
thetam=12.3; % degrees
f=122; % in HZ
thetadeg= -102.2; %degrees
thetanot= 1.83;
etam=87; %deg
Ceta=0.1;
phieta=-91.8; %degrees
etanot=-90; %degrees
tdependence=(0:.01:1);
phim=12.3;
K=1;
CT=1.341; % translational lift coefficient
alpha=28; %degrees mid-stroke
Crot= pi; % rotational lift coefficient
CD_0 = 0; % page 160 from Berman and Wang table 2
CD_pi = 2.93; % CD_.5pi page 160 from Berman and Wang table 2
phiddot=0;
%Equations of Motion
%C=[0:.1:5]
%T = 1/f;
H = t;
C = ((4*c)/pi)*sqrt(1-(r.^2/R.^2)); % C= c(r)
M22 = (pi/4)*rhof*C.^2;
M11 = (pi/4)*rhof*b^2;
theta=thetam*cos(2*(pi)*N*f*t+thetadeg)+thetanot;
%thetadot=f*N*thetam*thetanot*sin(2*pi*f*N*t+thetadeg)
thetadot=0 %just in foward flight
eta = (etam/tanh(Ceta))*tanh(Ceta*sin(2*pi*tdependence.*f + phieta))+ etanot;
etadot = 0; %just in foward flight
phi= (phim/(sin(K).^-1))*sin(K*sin(2*pi*f*t)).^(-1);
phidot = 0; %just in foward flight
%Vyp = r.*(thetadot*(cos(eta))-phidot*(cos(theta))*(sin(eta)));
Vyp = 0;
Vxp=r.*(phidot*cos(theta)*cos(eta)+thetadot*sin(eta));
%Vxp= 0;
V=sqrt((Vxp)^2+(Vyp)^2); %speed
Gamma= -0.5*CT*C*V*sin(2*alpha)+.5*Crot*C.^2*etadot;
%axp= r.*((phiddot*cos(theta)+thetadot*(etadot-phidot*sin(theta))*cos(eta)+...
% (thetaddot-etadot*phidot*cos(theta))*sin(eta)));
axp=0;
%ayp= r.*((thetadot*(etadot- phidot*sin(theta))*sin(eta)+(thetaddot-...
% etadot*phidot*cos(theta))*cos(eta)));
ayp=0;
dFvxp=0.5*rhof*C*(CD_0*(cos(alpha))^2+CD_pi*(sin(alpha))^2)*V*(Vxp);
dFvyp=0.5*rhof*C*(CD_0*(cos(alpha))^2+CD_pi*(sin(alpha))^2)*V*(Vyp);
dFxp=((((C/(c*R))*Mwing+M22)*Vyp*etadot-rhof*Gamma*Vyp-M11*axp))-dFvxp;
W=trapz(dFxp,r);
%Plots
%Figure plot(H,W)
plot(H,Vxp)