Re: matlab solver
Does the code that you gave me have a solution? I want to see if there's an error in my code.
On thing to think about is if you really want to search for the answer over the entire space - phi is probably some angle (correct) - would you expect to get an answer that had |phi| > 80?
This code does what is called a grid search. It steps over the space at a large step size (.5), finds the best value at that step size. Then, it looks at a range near that optimum at a smaller step size(.1), and finds the best value at that step size. This is repeated at (.01).
There is an error in the printf statement (at least on my version of matlab), but it is working.
PHP Code:
%%Constants%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rho0 = 1.225;
labda = -0.0065;
T0 = 288.15;
g0 = 9.80665;
R = 287.05;
%%Aircraft 1 parameters (index 1)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S_1 = 154;
W_1 = [0 0 -17463.49];
Cd0_1 = 0.01232;
AR_1 = 23;
e_1 = 0.8765;
%%Aircraft 2 parameters (index 2)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S_2 = 154;
W_2 = [0 0 -17463.49];
Cd0_2 = 0.01232;
AR_2 = 23;
e_2 = 0.8765;
%%Cable parameters%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigma_ultimate = 1.1E9;
j = 2;
rho_c = 700;
%%Inputs%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Jet stream inputs
H_delta = 1978.50;
W_delta = 16.48;
windratio = 0.7;
%Aircraft 1 inputs
AOA_1 = 9;
Vy_1 = [0 29 0];
Vz_1 = [0 0 0];
H_cruise_1 = 11000;
%Kite inputs
Vz_2 = [0 0 0];
%Cable inputs
d_c = 5;
%%Atmosphere calculations%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rho_1 = rho0*(1+labda*H_cruise_1/T0)^(-(g0/(R*labda)+1));
H_cruise_2 = H_cruise_1 - H_delta;
rho_2 = rho0*(1+labda*H_cruise_2/T0)^(-(g0/(R*labda)+1));
%%Aircraft calculations%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Lift and drag coefficient
Cl_1 = AOA_1*pi/180*2*pi*AR_1/(2+sqrt(4+AR_1^2));
Cd_1 = Cd0_1 + Cl_1^2/(pi*AR_1*e_1);
%Speed
Vx_1 = [W_delta*windratio 0 0];
V_1 = Vx_1 + Vy_1 + Vz_1;
V_1_size = sqrt(sum(V_1.^2));
%%LOOP: Change phi_1 until Fy_1 = 0%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
phi_1 = 0;
T(2) = 1;
while(abs(T(2)) > 0.01)
phi_1 = phi_1+0.0001;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Lift and drag
L_1_size = 0.5*rho_1*V_1_size^2*Cl_1*S_1;
L_1 = [-L_1_size*sin(phi_1*pi/180)*Vy_1(2)/V_1_size L_1_size*sin(phi_1*pi/180)*Vx_1(1)/V_1_size L_1_size*cos(phi_1*pi/180)];
D_1_size = 0.5*rho_1*V_1_size^2*Cd_1*S_1;
D_1 = [-D_1_size*Vx_1(1)/V_1_size -D_1_size*Vy_1(2)/V_1_size 0];
%%Cable calculations 1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Cable height 1
H_cable_1 = H_delta*windratio;
%Maximum cable tension
T_max = sigma_ultimate/j*(d_c/1000)^2*pi/4;
%Cable drag 1
Dx_c1 = -1/3*1/2*rho_1*Vx_1.^2*d_c/1000*H_cable_1;
Dy_c1 = -1/3*1/2*rho_1*Vy_1.^2*d_c/1000*H_cable_1;
Dz_c1 = -1/3*1/2*rho_1*Vz_1.^2*d_c/1000*H_cable_1;
D_c1 = Dx_c1 + Dy_c1 + Dz_c1;
%%Tension force equilibrium%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fx_1 = [L_1(1)+D_1(1)+D_c1(1) 0 0];
Fy_1 = [0 L_1(2)+D_1(2)+D_c1(2) 0];
Fz_1 = [0 0 L_1(3)+ D_1(3)+ D_c1(3)+W_1(3)];
T = Fx_1 + Fy_1 + Fz_1;
T_size = sqrt(sum(T.^2));
%cable angle
beta_1 = atan(Fz_1(3)/Fx_1(1))*180/pi;
end
%kite calculations%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%LOOP: vary AOA_2, Vy_2, phi_2 so that T1 = T2%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Vy_2 = [0 0 0];
minA=0;
minI=0;
minPhi=0;
minVal=inf;
ARange =[-2 14];
iRange =[-40 40];
phiRange = [-90 90]
for step = [.5 .1 .01]
for AOA_2 = ARange(1):step:ARange(2)
for i = iRange(1):step:iRange(2)
Vy_2(2) = i;
for phi_2 = phiRange(1):step:phiRange(2)
%Lift and drag coefficient
Cl_2 = AOA_2*pi/180*2*pi*AR_2/(2+sqrt(4+AR_2^2));
Cd_2 = Cd0_2 + Cl_2^2/(pi*AR_2*e_2);
%Speed
Vx_2 = [W_delta-Vx_1(1) 0 0];
V_2 = Vx_2 + Vy_2 + Vz_2;
V_2_size = sqrt(sum(V_2.^2));
%Lift and drag
L_2_size = 0.5*rho_2*V_2_size^2*Cl_2*S_2;
L_2 = [-L_2_size*sin(phi_2*pi/180)*Vy_2(2)/V_2_size L_2_size*sin(phi_2*pi/180)*Vx_2(1)/V_2_size L_2_size*cos(phi_2*pi/180)];
D_2_size = 0.5*rho_2*V_2_size^2*Cd_2*S_2*sign(V_2_size);
D_2 = [-D_2_size*Vx_2(1)/V_2_size -D_2_size*Vy_2(2)/V_2_size 0];
%%Cable calculations 1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Cable height 1
H_cable_2 = H_delta - H_cable_1;
%Cable drag 1
Dx_c2 = -1/3*1/2*rho_2*Vx_2.^2*d_c/1000*H_cable_2;
Dy_c2 = -1/3*1/2*rho_2*Vy_2.^2*d_c/1000*H_cable_2;
Dz_c2 = -1/3*1/2*rho_2*Vz_2.^2*d_c/1000*H_cable_2;
D_c2 = Dx_c2 + Dy_c2 + Dz_c2;
%%Tension force equilibrium%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fx_2 = [L_2(1)+D_2(1)+D_c2(1) 0 0];
Fy_2 = [0 L_2(2)+D_2(2)+D_c2(2) 0];
Fz_2 = [0 0 L_2(3)+ D_2(3)+ D_c2(3)+W_2(3)];
T_2 = Fx_2 + Fy_2 + Fz_2;
T_size_2 = sqrt(sum(T.^2));
%Cable angle
beta_2 = atan(Fz_2(3)/Fx_2(1))*180/pi;
dFx = Fx_1(1) - Fx_2(1);
dFy = Fy_1(2)^2 + Fy_2(2)^2;
dFz = Fz_1(3) + Fz_2(3);
dF = dFx^2 + dFy^2 +dFz^2;
if dF < 1
printf('done');
break;
elseif dF < minVal
minA=AOA_2;
minI=i;
minPhi=phi_2;
minVal = dF
end
end
% if dF < 1
% printf('done');
% break;
% end
end
% if dF < 1
% printf('done');
% break;
% end
end
ARange = [minA-4*step minA+4*step]
iRange = [minI-4*step minI+4*step]
phiRange = [minPhi-4*step minPhi+4*step]
end
Similar Threads
