6-DOF SIMULATION OF AIRCRAFT in Flight gear

版本 1.0.0 (223.6 KB) 作者: RAHUL N

Includes Simulink model and a live function, which animate the aircraft dynamic response using MATLAB animation and flight gear software.

关注

5.0

(2)

470.0 次下载

更新时间 2023/1/14

查看许可证

A Nonlinear, 6-DOF Dynamic Model of Aircraft: -

The Research Civil Aircarft Model. RAHUL N

Link for the Garteur web page GARTEUR

Link for RCAM Document RCAM Documentd

Note:- Here we are trying to simulate the similar sircraft like Boing 757-200

function[XDOT]=RCAM_model(X,U)

STATE VECTOR

% extract the state vector

x1=X(1); % Velocity in body x axis

x2=X(2); % Velocity in body y axis

x3=X(3); % Velocity in body z axis

x4=X(4); % angular rate in body x axis

x5=X(5); % angular rate in body x axis

x6=X(6); % angular rate in body x axis

x7=X(7); % phi (Bank angle)

x8=X(8); % theta (Pitch angle)

x9=X(9); % phi (Yaw angle)

CONTROL VECTOR

% extract the control vector

u1=U(1); %d_A (aileron)

u2=U(2); %d_T (stabilizer) here the whole horizontal stabilzer moves not the elevator

u3=U(3); %d_R (rudder)

u4=U(4); %d_th1 (throttle 1)

u5=U(5); %d_th2 (throttle 2)

CONSTANTS

Here we are trying to enter all the constant value by the RCAM research value.

Note all values entered are in the SI unit.

% Nominal vehicle constants

m=120000; % Aircraft total mass

% NOTE:- We will define the Ib and invIb later

cbar=6.6; % Mean aerodynamic chord

lt=24.8; % Distance by AC of tail and body (m)

S=260; % Wing planform area (m^2)

St=64; % Tail planform area (m^2)

Xcg=0.23*cbar; % X Position of cg in Fm (m)

Ycg=0; % Y Position of cg in Fm (m)

Zcg=0.10*cbar; % Z Position of cg in Fm (m)

Xac=0.19*cbar; % X position of aerodynamic center in Fm (m)

Yac=0; % Y position of aerodynamic center in Fm (m)

Zac=0; % Z position of aerodynamic center in Fm (m)

% Engine Constants

Xapt1=0; % X position of engine 1 force in Fm (m)

Yapt1=-7.94; % Y position of engine 1 force in Fm (m)

Zapt1=-1.9; % Z position of engine 1 force in Fm (m)

Xapt2=0; % X position of engine 2 force in Fm (m)

Yapt2=7.94; % Y position of engine 2 force in Fm (m)

Zapt2=-1.9; % Z position of engine 2 force in Fm (m)

% Other constants

rho=1.225; % Air density in (kg/m^3)

g=9.81; % Acceleration due to gravity (m/s^2)

depsda=0.25; % Change in downwash w.r.t alpha (rad/rad)

alpha_L0=-11.5*(pi/180); % Zero lift angle of attack (rad)

n=5.5; % Slope of linear region of lift slope

a3=-768.5; % coefficient of alpha^3

a2=609.2; % coefficient of alpha^2

a1=-155.2; % coefficient of alpha^1

a0=-15.212; % coefficient of alpha^0 (Different from RCAM Model)

alpha_switch=14.5*(pi/180); % alpha value lift slope goes from linear to non-linear

CONTROL SATURATION

These control saturation are defined in the Simulink Model. Thus the code is commented

% Note that these can alternatively be enforced in simulink

% u1min=-15*(pi/180);

% u1max=25*(pi/180);

% u2min=-25*(pi/180);

% u2max=10*(pi/180);

% u3min=-30*(pi/180);

% u3max=30*(pi/180);

% u4min=0.5*(pi/180);

% u4max=10*(pi/180);

% u5min=0.5*(pi/180);

% u5max=10*(pi/180);

% if(u1>u1max)

% u1=u1max;

% elseif(u1<u1min)

% u1=u1min;

% end

% if(u2>u2max)

% u2=u2max;

% elseif(u2<u2min)

% u2=u2min;

% end

% if(u3>u3max)

% u3=u3max;

% elseif(u3<u3min)

% u3=u3min;

% end

% if(u4>u4max)

% u4=u1max;

% elseif(u4<u4min)

% u4=u4min;

% end

% if(u5>u5max)

% u5=u5max;

% elseif(u5<u5min)

% u5=u5min;

% end

INTERMEDIATE VARIABLES

All the intermediate variables are expressed in terms of the control and state vector

which implies

These are the intermediate variables that are already described below

ASSUMPTION

The Density is assumed to be constant 1.225 kg/m^3.

There is no time lapse between the control surface transition, means the control surface movement is fast

The throttle position from 0 to1 is not dependent on time

Reverse stall is not considered

CL is only a function of

As we are creating a function which takes the input of control parametrs only, therefore whatever the input give to the system should be a function of states

Thus we are expressing all the Intermediate variable as a function of U

The actual equation are given as follows: -

% calculate the airspeed

Va=sqrt(x1^2+x2^2+x3^2);

%claculate the alph and beta

alpha=atan2(x3,x1);

beta=asin(x2/Va);

%claculation of the dynamic pressure

Q=0.5*rho*Va^2;

%Also define the vectors wbe_b and V_b

wbe_b=[x4;x5;x6];

V_b=[x1;x2;x3];

NON DIMENSIONAL AERODYNAMIC FORCES COEFFICIENTS IN THE STABILITY FRAME

Here we are trying to impliment the nonlinear part of Cl curve

where,

The is only a function of alpha

We are including the stall characteristics

Did not included reverse stall condition

is in the stability frame

Lift Generated by the tail

The general equation for the tail angle of attcak is given by

where, is downwash parameter

NOTE:- Here in the horizontal stabilizer the entire tail goes up and down (all movable tail)

here we are trying to model the downwash

=0.25

= Dynamic pitch response of the aircraft

Tail moment-arm of the aircarft

The total contribution from the tail is given by

Total Lift Coefficient of the aircraft

Total drag Coefficient is given by

Till now we can note that there is no beta dependence

Total Side Force

% Calculate the CL_wb

if alpha<=alpha_switch

CL_wb=n*(alpha-alpha_L0);

else

CL_wb=a3*alpha^3+a2*alpha^2+a1*alpha+a0;

end

%Claculate CL_t

epsilon=depsda*(alpha-alpha_L0);

alpha_t=alpha-epsilon+u2+1.3*x5*lt/Va;

CL_t=3.1*(St/S)*alpha_t;

% claculation of the Total lift force

CL=CL_wb+CL_t;

% Total drag force (Neglecting the tail)

CD=0.13+0.07*(5.5*alpha+0.654)^2;

% Calculate the side force

CY=-1.6*beta+0.24*u3;

DIMENSIONAL AERDOYNAMIC FORCES

note, here were are rotating the aerodynamic forces from the stability frame to body frame

% calculate the actual dimensional forces. These are in F_s (Stability axis)

FA_s=[-CD*Q*S;CY*Q*S;-CL*Q*S];

Rotating the forces from the stability frame to body frame

% Rotate these forces to F_b (Body axis)

C_bs=[cos(alpha) 0 -sin(alpha);0 1 0;sin(alpha) 0 cos(alpha)];

FA_b=C_bs*FA_s;

NON-DIMENSIOANL AERODYNAMIC MOMENT COEFFICIENT ABOUT THE AC IN THE BODY FRAME

Non dimensional Aeromoment coefficient cabout AC in the body frame is given by:-

where,

%Calculate the moments in Fb. Define eta,dCMdx and dCMdu

eta11=-1.4*beta;

eta21=-0.59-(3.1*(St*lt)/(S*cbar))*(alpha-epsilon);

eta31=(1-alpha*(180/(15*pi)))*beta;

eta=[eta11;eta21;eta31];

dCMdx=(cbar/Va)*[-11 0 5;0 (-4.03*(St*lt^2)/(S*cbar^2)) 0;1.7 0 -11.5]; % needed correcction

dCMdu=[-0.6 0 0.22;0 (-3.1*(St*lt)/(S*cbar)) 0;0 0 -0.63];

% Now calculate CM=[cl;cm;cn] about aerodynamic center in Fb

CMac_b=eta+dCMdx*wbe_b+dCMdu*[u1;u2;u3];

AERODYNAMIC MOMENT ABOUT AC IN THE BODY FRAME

%Normalize to an aerodynamic moments

MAac_b=CMac_b*Q*S*cbar;

AERODYNAMIC MOMENT ABOUT CG IN BODY FRAME

For moment transfer is given by

equation for the moment transfer

%Transfer moment to CG

rcg_b=[Xcg;Ycg;Zcg];

rac_b=[Xac;Yac;Zac];

MAcg_b=MAac_b+cross(FA_b,rcg_b-rac_b);

PROPULSION EFFECTS

Forces and Moments due to engine

thrust to weight ratio of the aircraft is given by

Total Propulsive force in the body frame is given by

where, is the distance between the CG of the aircraft and the application point of thrust.

Data for the above matrix

% First, calculate the thrust of each engine

F1=u4*m*g;

F2=u5*m*g;

%Assuming that the engine thrust is aligned with the Fb, we have

FE1_b=[F1;0;0];

FE2_b=[F2;0;0];

FE_b=FE1_b+FE2_b;

%Now the engine moment due to offset of engine thrust from cog

mew1=[Xcg-Xapt1;Yapt1-Ycg;Zcg-Zapt1];

mew2=[Xcg-Xapt2;Yapt2-Ycg;Zcg-Zapt2];

MEcg1_b=cross(mew1,FE1_b);

MEcg2_b=cross(mew2,FE2_b);

MEcg_b=MEcg1_b+MEcg2_b;

GRAVITY EFFECTS

Here we are considering the gravity effects to the aircraft

Let us consider earth frame Fe

We need to rotate to the body frame

as the gravity acts at the CG of the aircraft

% Calculate the gravitational forces in the body frame. This causes no moments abot the CG

g_b=[-g*sin(x8);g*cos(x8)*sin(x7);g*cos(x8)*cos(x7)];

Fg_b=m*g_b;

EXPLICIT FIRST ORDER FORM

Cosidering the FLAT EARTH EQUATIONS

Assuming that the mass is constant

where,

Euler Kinematical equation

%Inertia matrix

Ib=m*[40.07 0 -2.0923;0 64 0;-2.0923 0 99.92];

% Inverse of the inertia matrix

invIb=(1/m)*[0.0249836 0 0.000523151;0 0.015625 0;0.000523215 0 0.010019];

% From F_b ( all the forces in Fb ) and calculate the udot, vdot,wdot

F_b=Fg_b+FE_b+FA_b;

x1to3dot=(1/m)*F_b-cross(wbe_b,V_b);

%From Mcg_b (all the moments about the cog in Fb) and calculate pdot,qdpt,rdoot

Mcg_b=MAcg_b+MEcg_b;

x4to6dot=invIb*(Mcg_b-cross(wbe_b,Ib*wbe_b));

% Calculation of the phidot,thetadot,psidot

H_pi=[1 sin(x7)*tan(x8) cos(x7)*tan(x8);0 cos(x7) -sin(x7);0 sin(x7)/cos(x8) cos(x7)/cos(x8)];

x7to9dot=H_pi*wbe_b;

NAVIGATION EQUATION

To calculate the position of the aircraft in North, East and Down Direction

Direction cosine matrix

Navigational Eqaution

% Rotational Matrix from the 1 frame to 2 frame

C1v=[cos(x9) sin(x9) 0;

-sin(x9) cos(x9) 0;

0 0 1];

% Rotational Matrix from the 1 frame to 2 frame

C21=[cos(x8) 0 -sin(x8);

0 1 0;

sin(x8) 0 cos(x8)];

% Rotational Matrix from the 1 frame to 2 frame

Cb2=[1 0 0;

0 sin(x7) cos(x7);

0 -sin(x7) cos(x7)];

% Implimenting the Navigational Equation

Cbv=Cb2*C21*C1v;

Cvb=Cbv';

x10to12dot=Cvb*V_b;

EXPLICIT FIRST ORDER FORM

XDOT=[x1to3dot;

x4to6dot;

x7to9dot;

x10to12dot;

];

引用格式

RAHUL N (2024). 6-DOF SIMULATION OF AIRCRAFT in Flight gear (https://www.mathworks.com/matlabcentral/fileexchange/123285-6-dof-simulation-of-aircraft-in-flight-gear), MATLAB Central File Exchange. 检索来源 2024/4/28.