-
Notifications
You must be signed in to change notification settings - Fork 0
/
Hover.m
143 lines (119 loc) · 4.46 KB
/
Hover.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
%Free Wake Method - Hover Version 1.3.2
% - Aim: Self induced velocity
% - Comment: only self induced velocity tip vortex
% - Comment: check if works: 2 running variables, reason to prevent ra==rp
% or rb==rp -> leads to NaN with Biot Savart
%initialization
close all;
%input variables for discretisation
R=10;
Nr=10;
Nb=1; %Note: Functions only with one Nb atm, TODO: make adjustments to others
t_rev=2; %time it takes for 1 revolution
omega=-(2*pi/t_rev);
t_final=1;
dt=0.05;
t_iter=ceil(t_final/dt)+1;
%wake input
p_wake=zeros(3,t_iter,Nr); %(x,y,z), progresses with time, discretisation Nr
Gamma_blade=zeros(Nr-1,1); %only one Blade for now
Gamma_wake=zeros(t_iter,Nr); %Gamma variable for Nr with time
V=zeros(3,t_iter,Nr); %Velocity components for each point p_wake
rc=0.5;
%rotational Matrix
Rot_M=[cos(omega*dt) sin(omega*dt) 0; ...
-sin(omega*dt) cos(omega*dt) 0; 0 0 1];
r=linspace(0.1*R,R,Nr); %seperation along r
p=zeros(3,Nr,Nb); %p_blade decleration
%blade configuration, initialization
for k=1:Nb
p(1,:,k)=r*cos((2*pi/Nb)*(k-1));
p(2,:,k)=r*sin((2*pi/Nb)*(k-1));
end
%Gamma_blade (simplification: Gamma=-4*pi*sin(alpha)*q_inf
sin_alpha=sin(10*pi/180); %assumption 10 deg constant along r, TODO: subprogram variable alpha(r)
q_inf=zeros(Nr,1);
q_inf=r*abs(omega); %flow velocity component (only tangential to blade atm)
Gamma_temp=-4*pi*sin_alpha*q_inf;
for ii=1:Nr-1
Gamma_blade(ii)=0.5*(Gamma_temp(ii)+Gamma_temp(ii+1)); %constant gamma distribution along Blade
end
% Simulation
for ii=1:1:t_iter
%shed wake
for tt=2:Nr-1 %for each timestep the wake is shed, based on previous gamma_blade along r
Gamma_wake(ii,1)=-Gamma_blade(1);
Gamma_wake(ii,tt)=Gamma_blade(tt-1)-Gamma_blade(tt);
Gamma_wake(ii,Nr)=Gamma_blade(Nr-1);
end
for jj=1:Nr
%position, coordinates
p(:,jj,k)=Rot_M*p(:,jj,k); %new blade position
p_wake(:,ii,jj)=p(:,jj,k); %new wake position
end
% ----------------------------------------------------------------%
%BiotSavart - Euler
%function [ V ] = BiotSavartLaw3( ra,rb,rp,Gamma )
%p_wake=zeros(3,t_iter,Nr);
%Gamma_wake=zeros(t_iter,Nr);
%V=zeros(3,t_iter,Nr);
if (ii>1) %only calculate when timestep is bigger than 1
V(:,:,:)=0; %clear V before each new iteration
for kk=2:ii-1
for ll=1:kk-1
V(:,ll,Nr)=V(:,ll,Nr)+(BiotSavartLaw3(p_wake(:,kk,Nr),p_wake(:,kk+1,Nr),...
p_wake(:,ll,Nr),Gamma_wake(kk,end)));
if ll>2
for mm=1:ll-2
V(:,ll,Nr)=V(:,ll,Nr)+(BiotSavartLaw3(p_wake(:,mm,Nr),p_wake(:,mm+1,Nr),...
p_wake(:,ll,Nr),Gamma_wake(kk,end)));
end
end
%implementation notes
%V=(u,v,w), till current time point ii NOTE: should it be ii-1? flow
%condition, along Nr (ignoring wake tip)
%ra=outer p_wake (kk=ii)
%rb=outer p_wake (kk=ii+1)
%rp=every point besides on the outer wake itself
%Gamma=only tip vorticies
end
end
% ----------------------------------------------------------------%
%PROBLEM AREA: Velocity self inducing tip vortex
% for kk=1:ii-1
% V(:,1,Nr)=sum(BiotSavartLaw3(p_wake(:,kk,Nr),p_wake(:,kk+1,Nr),...
% p_wake(:,1,Nr),Gamma_wake(kk,end)));
% end
%problem wenn p_wake and p_gamma are on the same point -> NaN
%reason being:
% ----------------------------------------------------------------%
%Euler Equation
%p_wake()=p_wake()+dt*V();
%Note: vi is only used for displaying purposes, from Momentum Theory
%assumption: vi is added to every component, as a constant
%in hover, vi=sqrt(Thrust/(2*rho*pi*R^2)), Thrust=m*g
vi=zeros(3,1);
vi(3)=0.1;
for oo=1:Nr
for pp=1:ii-1
p_wake(:,pp,oo)=p_wake(:,pp,oo)+dt*V(:,pp,oo)-vi(:);
end
end
% ----------------------------------------------------------------%
end %if ii>1
iteration=ii
end
%plot end result
figure(1)
%plot near wake using plot3
%for ii=1:Nr
plot3(p_wake(1,:,Nr),p_wake(2,:,Nr),p_wake(3,:,Nr));
hold on;
%end
%plot current blade position
plot3(p(1,:,:),p(2,:,:),p(3,:,:),'k','LineWidth',3)
axis([-2*R 2*R -2*R R*2 -10 0]);
ylabel('y');
xlabel('x');
zlabel('z');
title('Wake for hover flight, one blade')