Question on Aircraft Drone Modeling - Simulation - Data Selection and Entry
9 views (last 30 days)
Show older comments
Okay, I'll be blunt. I'm trying to learn about aircraft aerodynamics by working through a textbook on drone design. I have no instructor, and I'm 30 years distant from my own college degree, so I've forgotten a step or two. That means when matlab has a segment of code that say that the term that is added is m[J K L] (vertically aligned) I'm at something of a loss regarding what I'm adding. Do I add m(J) to the topmost line, then add m(K) to the second line equation, and m(L) to the bottom line?
The first part assignment (Chapter 2) was to produce a model of the an aircraft in Matlab and get it to move on the screen. I did this. The code is below:
The pn, pe, and pd variables are apparently location variables. The dpsi is a change in the horizontal orientation of the aircraft (it turns right or left). The pd = pd - 1 threw me for a moment, until I realized what I'd done was to increment the direciton of motion with a negative one to make the plane move in the direction it is pointing. I may need to correct this. The following segment of code has no tie to any real estimate of forces acting on the drone. It is purely my data producing an aircraft that moves on the scree.
pn=0;
pe=0;
pd=0;
phi=0;
theta=0;
psi=0;
handle=[];
%Simulation parameters
dpn=0.5; %Change in position.
dpsi=0.01; %Change in yaw angle.
simlength=125;%Number of Sim Steps
%Draw and update the plane's position.
for k = 1:simlength
pn = 75+dpn;
pe = pe+dpn;
%THE CODE THAT PRODUCES A PROBLEM IS BELOW
pd = pd -1;
psi = psi+dpsi;
%Draw or update the plane.
handle = drawPlaneBody(pn,pe,pd,phi,theta,psi,handle);
pause(0.1); %Pause to visualize the movement
end
function handle = drawPlaneBody(pn,pe,pd,phi,theta,psi,handle)
%define points on plane in local NED coordintates
NED = airplanepoints;
%rotate plane by (phi; theta; psi)
NED = rotate(NED, phi, theta, psi);
%translate plane to [pn; pe; pd]
NED = translate(NED, pn, pe, pd);
%transform vertices from NED to XYZ
R = [...
0, 1, 0;...
1, 0 , 0;...
0, 0, -1;...
];
XYZ = R* NED';
%plot plane
if isempty(handle)
handle = plot3(XYZ(1,:),XYZ(2,:),XYZ(3,:),'b');
xlim([-100 150]);
ylim([-100 150]);
zlim([0 150]);
grid on;%show grid
xlabel('X');ylabel('Y');zlabel('Z');%label axes
else
set(handle, 'XData', XYZ(1,:), 'YData',XYZ(2,:), 'ZData',XYZ(3,:));
drawnow
end
end
function XYZ = airplanepoints
%define points on the aircraft in local NED
XYZ = [...
0 0 0;%point1
-2 1 1;%point2
-2 1 -1;%point3
0 0 0;%point1
-2 -1 1;%point4
-2 -1 -1;%point5
0 0 0;%point1
-2 1 1;%point2
-2 -1 1;%point4
0 0 0;%point1
-2 1 -1;%point3
-2 -1 -1;%point5
0 0 0;%point1
-2 1 1;%point2
-18 0 0;%point6
-2 1 -1;%point3
-18 0 0;%point6
-2 -1 -1;%point5
-18 0 0;%point6
-2 -1 1;%point4
-18 0 0;%point6
-2 1 1;%point2
0 0 0;%point1
-5 0 0;%point7
-5 -10 0;%point8
-8 -10 0;%point9
-8 10 0;%point10
-5 10 0;%point11
-5 0 0;%point7
-15.5 0 0;%point12
-15.5 2 0;%point13
-17.5 2 0;%point14
-17.5 -2 0;%point15
-15.5 -2 0;%point16
-15.5 0 0;%point12
-18 0 0;%point6
-18 0 -3;%point17
-15.5 0 0;%point15
-18 0 0;%point16
];
end
function XYZ=rotate(XYZ,phi,theta,psi)
%define rotation matrix
R_roll = [
1, 0, 0;
0, cos(phi), -sin(phi);
0, sin(phi), cos(phi)];
R_pitch = [
cos(theta), 0, sin(theta);
0, 1, 0;
-sin(theta), 0, cos(theta)];
R_yaw = [
cos(psi), -sin(psi), 0;
sin(psi), cos(psi), 0;
0, 0, 1];
R = R_roll*R_pitch*R_yaw;
%rotate vertices
XYZ =R*XYZ';
XYZ = XYZ';
end
function XYZ = translate(XYZ, pn, pe, pd)
XYZ = XYZ' + repmat([pn;pe;pd],1,size(XYZ,1));
XYZ=XYZ';
end
Alright, this is the end of the purely "draw a drone" segment ofthe code. The book then says to produce code that will permit the aircraft to respond to real world forces. That code follows. I haven't even tried to run it. I'm not sure that I even have the right variables transferred to the functions. The equations generally serve to translate forces into a reference frame that is relevant to the aircraft. I haven't check which reference frame is the starting one and which is the final one. I leave the movement code alone and maintain the same size box as in the prior simulation. I'm hoping someone wtho knows something about this book, Small, Unmanned Aircraft Theory and Practice will find this and assist. I need some instructor level guidance. Maybe theres an aircraft designer who understands this subject. I'll post this and hope. (The G1, G2, G3, etc. in the following code is a capital gamma with those numbers afterward in the text. There's also a Jy at the end of the line of data that is transferred. I don't know what that is. Does there need to be a Jx and a Jz transferred as well?)
pn=0;
pe=0;
pd=0;
phi=0;
theta=0;
psi=0;
handle=[];
%Simulation parameters
dpn=0.5; %Change in position.
dpsi=0.01; %Change in yaw angle.
%SIMULATION
u = 1;
v = 1;
w = 1;
theta = 3.1415/8;
psi = 3.1415/9;
phi = 3.1415/10;
l = 2;
m =3;
n=4;
p = 5;
q = 6;
r = 7;
G1 = 1;
G2 = 1.1;
G3 = 1.2;
G4 = 1.3;
G5 = 1.4;
G6 = 1.5;
G7 = 1.6;
G8 = 1.7;
Jy = 2;
Mavequations(u,v,w,theta,psi,phi,l,m,n,p,q,r,G1,G2,G3,G4,G5,G6,G7,G8,Jy);
%simlength=125;%Number of Sim Steps
%Draw and update the plane's position.
%for k = 1:simlength
% pn = 75+dpn;
% pe = pe+dpn;
%THE CODE THAT PRODUCES A PROBLEM IS BELOW
% pd = pd -1;
% psi = psi+dpsi;
% %Draw or update the plane.
% handle = drawPlaneBody(pn,pe,pd,phi,theta,psi,handle);
% pause(0.1); %Pause to visualize the movement
%
function MAVEquations(u,v,w,theta,psi,phi,l,m,n,p,q,r,G1,G2,G3,G4,G5,G6,G7,G8,Jy)
pn = cos(theta)*cos(psi)*sin(phi)*sin(theta)*cos(psi)*u-cos(phi)*sin(psi)*cos(phi)*sin(theta)*cos(psi)*v+sin(phi)*sin(psi)*w;
pe = cos(theta)*sin(psi)*sin(phi)*sin(theta)*sin(psi)*u +cos(phi)*cos(psi)*cos(phi)*sin(theta)*sin(psi)*v+sin(phi)*cos(psi)*w;
pd = -sin(theta)*u + sin(phi)*cos(theta)*v + cos(phi)*cos(theta)*w;
ud = r*v-q*w +(fx+fy+fz)/m;
vd = p*w-r*u +(fx+fy+fz)/m;
wd = q*u-p*v +(fx+fy+fz)/m;
phid =1*p + sin(phi)*tan(theta)*q + cos(phi)*cos(theta)*r;
thetad = 0*p + cos(phi)*q - sin(phi);
psid = 0 * p + sin(phi)/cos(theta) *q + cos(phi)/cos(theta)*r;
pd = G1*p*q-G2*q*r + G3*l+G4*n;
qd = G5*p*r-G6(p^2-r^2)+m/Jy;
rd = G7*p*q-G1*q*r+G4*l+G8*n;
handle = drawPlanBody(pn,pe,pd,phi,theta,psi,handle);
pause(0.1);
end
function handle = drawPlaneBody(pn,pe,pd,phi,theta,psi,handle)
%define points on plane in local NED coordintates
NED = airplanepoints;
%rotate plane by (phi; theta; psi)
NED = rotate(NED, phi, theta, psi);
%translate plane to [pn; pe; pd]
NED = translate(NED, pn, pe, pd);
%transform vertices from NED to XYZ
R = [...
0, 1, 0;...
1, 0 , 0;...
0, 0, -1;...
];
XYZ = R* NED';
%plot plane
if isempty(handle)
handle = plot3(XYZ(1,:),XYZ(2,:),XYZ(3,:),'b');
xlim([-100 150]);
ylim([-100 150]);
zlim([0 150]);
grid on;%show grid
xlabel('X');ylabel('Y');zlabel('Z');%label axes
else
set(handle, 'XData', XYZ(1,:), 'YData',XYZ(2,:), 'ZData',XYZ(3,:));
drawnow
end
end
function XYZ = airplanepoints
%define points on the aircraft in local NED
XYZ = [...
0 0 0;%point1
-2 1 1;%point2
-2 1 -1;%point3
0 0 0;%point1
-2 -1 1;%point4
-2 -1 -1;%point5
0 0 0;%point1
-2 1 1;%point2
-2 -1 1;%point4
0 0 0;%point1
-2 1 -1;%point3
-2 -1 -1;%point5
0 0 0;%point1
-2 1 1;%point2
-18 0 0;%point6
-2 1 -1;%point3
-18 0 0;%point6
-2 -1 -1;%point5
-18 0 0;%point6
-2 -1 1;%point4
-18 0 0;%point6
-2 1 1;%point2
0 0 0;%point1
-5 0 0;%point7
-5 -10 0;%point8
-8 -10 0;%point9
-8 10 0;%point10
-5 10 0;%point11
-5 0 0;%point7
-15.5 0 0;%point12
-15.5 2 0;%point13
-17.5 2 0;%point14
-17.5 -2 0;%point15
-15.5 -2 0;%point16
-15.5 0 0;%point12
-18 0 0;%point6
-18 0 -3;%point17
-15.5 0 0;%point15
-18 0 0;%point16
];
end
function XYZ=rotate(XYZ,phi,theta,psi)
%define rotation matrix
R_roll = [
1, 0, 0;
0, cos(phi), -sin(phi);
0, sin(phi), cos(phi)];
R_pitch = [
cos(theta), 0, sin(theta);
0, 1, 0;
-sin(theta), 0, cos(theta)];
R_yaw = [
cos(psi), -sin(psi), 0;
sin(psi), cos(psi), 0;
0, 0, 1];
R = R_roll*R_pitch*R_yaw;
%rotate vertices
XYZ =R*XYZ';
XYZ = XYZ';
end
function XYZ = translate(XYZ, pn, pe, pd)
XYZ = XYZ' + repmat([pn;pe;pd],1,size(XYZ,1));
XYZ=XYZ';
end
Yes, this is a Hail Mary appeal. I hope that the right person will find it and help me to develop an answer. Thank you very much.
0 Comments
Answers (0)
See Also
Categories
Find more on Axes Transformations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!