16.07 Dynamics Fall 2004 Version 1.0 Lecture D3- Equations of Motion in Cartesian Coordinates 2DOF Aircraft model Newton's second law F=ma is a vector equation that relates the magnitude and direction of the force vector, to the magnitude and direction of the acceleration vector. In the previous lecture we derived expressions for the acceleration vector expressed in cartesian coordinates. This expressions can now be used in Newtons second law, to produce the equations of motion expressed in cartesian coordinates Equations of Motion in Cartesian Coordinates In rectangular cartesian coordinates, ayz, we write F=Fri+ Fy+ Fik and a=a22+ayj+a,k. Thus, ve have, in component form F n. F This system of Ordinary Differential Equations(ODE's) must be augmented with appropriate initial condi- tions. Typically we will know the initial position ro and velocity vo. Thus, r(0)=r(0)+y(0)+z)k=ro andv(0)=x(0)+y(0)j+2(0)k=t0 In general, forces on a body depend in general on its position r, its velocity v and, sometimes, explicitly on e time t(e.g. the control inputs on an aircraft). In addition, they may also depend on the position and velocity of other neighboring bodies or particles. In this situation, we may need to solve for the motion of all the bodies involved simultaneously. Fortunately, one can often make the approximation that the effect of he surrounding bodies is small and hence simplify the problem considerably. Analytical Integration- Some General Cases There are only a few instances where the equations of motion can be solved analytically. In the general cases numerical integration is required. Among the special situations where analytical solution is possible, there are the problems where F is either constant or depends on t only. In such cases we can integrate the above
J. Peraire 16.07 Dynamics Fall 2004 Version 1.0 Lecture D3 - Equations of Motion in Cartesian Coordinates (2DOF Aircraft Model) Newton’s second law F = ma is a vector equation that relates the magnitude and direction of the force vector, to the magnitude and direction of the acceleration vector. In the previous lecture we derived expressions for the acceleration vector expressed in cartesian coordinates. This expressions can now be used in Newton’s second law, to produce the equations of motion expressed in cartesian coordinates. Equations of Motion in Cartesian Coordinates In rectangular cartesian coordinates, xyz, we write F = Fxi + Fyj + Fzk and a = axi + ayj + azk. Thus, we have, in component form Fx = max = mx¨ Fy = may = my¨ Fz = maz = mz . ¨ (1) This system of Ordinary Differential Equations (ODE’s) must be augmented with appropriate initial conditions. Typically we will know the initial position r0 and velocity v0. Thus, r(0) = x(0)i+y(0)j+z(0)k = r0, and v(0) = ˙x(0)i + ˙y(0)j + ˙z(0)k = v0. In general, forces on a body depend in general on its position r, its velocity v and, sometimes, explicitly on the time t (e.g. the control inputs on an aircraft). In addition, they may also depend on the position and velocity of other neighboring bodies or particles. In this situation, we may need to solve for the motion of all the bodies involved simultaneously. Fortunately, one can often make the approximation that the effect of the surrounding bodies is small and hence simplify the problem considerably. Analytical Integration – Some General Cases There are only a few instances where the equations of motion can be solved analytically. In the general cases numerical integration is required. Among the special situations where analytical solution is possible, there are the problems where F is either constant or depends on t only. In such cases we can integrate the above 1
uations independently to compute the components of the velocity vector. v(t)=v(0)+ F:()l,n()=n0)+1/F()d,n2()=n20)+/F(出 and the position ()=0+mb2(),y(t)=01 m)dt,2()=20)+m/n:()d There are also one-dimensional cases in which the equations of motion can be integrated analytically(see review on rectilinear motion). In some other situations, it is possible to partially solve these equations nalytically. For instance, if the forces can be derived from a potential, we might be able to write expressions e.g. conservation of total energy, that are satisfied by the motion even though we may not know exactly know what the motion is Numerical Integration The most general way of solving the equations of motion is by numerical integration. In this case we do not compute the solution of the problem but an approximation to it. It is typically useful to work only with first order equations, so we writ 1 -Fr(a, y, z, vx, Vy, U2, t) . Fy(a, y, 2,t) 1 F2(x,3,2,Ux,Uy,U2,t) This is a set of 6 first order odes, with initial conditions x(0)=r0,y(0)=30,z(0)=20,v(0)=tx0,vy(O)=vy0,v2(0) The above set of equations is sometimes written in shorthand notation as X=F(X, t), with X(O=Xo where X=(r, g, 2, Ux, Uy, U2), F=(ux, Uy, U2, F=/m, Fy/m, F:/m), and Xo =(ro, y0, 20, 020, Uy0, vxo) The array X is usually referred to as the vector of state variables, even though it is not truly a vector. It hould not be confused with at physical vector such as the velocity or acceleration. Note that the components of x do not even have the same units ortunately, the theory for numerically solving systems of ordinary differential equations, is very mature and good algorithms are available. Here we use a simple"Forward Euler Method" for illustration. Consider a time increment At(this At should be"small enough"so that the computed results are accurate, but on the
equations independently to compute the components of the velocity vector, vx(t) = vx(0) + 1 m Z t 0 Fx(t) dt , vy(t) = vy(0) + 1 m Z t 0 Fy(t) dt , vz(t) = vz(0) + 1 m Z t 0 Fz(t) dt , and the position, x(t) = x(0) + 1 m Z t 0 vx(t) dt , y(t) = y(0) + 1 m Z t 0 vy(t) dt , z(t) = z(0) + 1 m Z t 0 vz(t) dt . There are also one-dimensional cases in which the equations of motion can be integrated analytically (see review on rectilinear motion). In some other situations, it is possible to partially solve these equations analytically. For instance, if the forces can be derived from a potential, we might be able to write expressions, e.g. conservation of total energy, that are satisfied by the motion even though we may not know exactly know what the motion is. Numerical Integration The most general way of solving the equations of motion is by numerical integration. In this case we do not compute the solution of the problem but an approximation to it. It is typically useful to work only with first order equations, so we write, x˙ = vx y˙ = vy z˙ = vz v˙x = 1 m Fx(x, y, z, vx, vy, vz, t) v˙y = 1 m Fy(x, y, z, vx, vy, vz, t) v˙z = 1 m Fz(x, y, z, vx, vy, vz, t) This is a set of 6 first order ODE’s, with initial conditions x(0) = x0, y(0) = y0, z(0) = z0, vx(0) = vx0, vy(0) = vy0, vz(0) = vz0 . The above set of equations is sometimes written in shorthand notation as X˙ = F(X, t), with X(0) = X0 (2) where X = (x, y, z, vx, vy, vz) T , F = (vx, vy, vz, Fx/m, Fy/m, Fz/m) T , and X0 = (x0, y0, z0, vx0, vy0, vz0) T . The array X is usually referred to as the “vector” of state variables, even though it is not truly a vector. It should not be confused with at physical vector such as the velocity or acceleration. Note that the components of X do not even have the same units. Fortunately, the theory for numerically solving systems of ordinary differential equations, is very mature and good algorithms are available. Here we use a simple “Forward Euler Method” for illustration. Consider a time increment ∆t (this ∆t should be “small enough” so that the computed results are accurate, but on the 2
other hand, we will see that the smaller the timestep the more time steps we will have to compute). The idea is to start from the initial condition, and compute the solution at times At, 2At, 3At, etc. Supp that we have already progressed along the trajectory to time t, when the components of X(t) are all known (this could clearly be t=0 because of the initial condition). We then update all these quantities to the next discrete time t+△t,by X(t+△t)=X(t)+△F(X(t),t) Solving ODE's with MATLABE The task of integrating a systems of ODE's is particluarly simple if we use a computer package such as MATLAB. There exist a number of solution algorithms that one can choose from, but for our purposes almost any of them will work. The MATLAB syntax is typically of the form t, X]= solver(odefun, tspan, xo) The returned arrays t and X contain the vector of time points and the solution, X at these time points respectively. A function, odefun, that calculates F, from X, and t, must be provided, as well as a vector pecifying the interval of integration tspan =[to, tf], and the initial value of x, xo. In order to set the parameters for the solution algorithm one can use the MATLAB function odeset Seehttp://www.mathworks.com/access/helpdesk/help/techdoc/ref/ode45.htmlfordetailsaswellas the lab 1 handout 2 Degrees-of-Freedom(2DOF) Aircraft Model In this section we will develop a simple model to simulate the fight of an aircraft(or spacecraft). The aircraft is treated as a particle. Therefore, we assume that all its mass is concentrated at a point located at the aircrafts center of mass. In addition, we will assume that we have perfect attitude control. That is, we will decouple the motion of the center of mass with the rotation of the aircraft about the center of mass. For the later, we will assume that we have a control system that allows us to prescribe the angle of attach of the aircraft. In reality this is clearly not possible and any moment applied to the aircraft, through he control surfaces, to change its attitude will take a certain time to be effective. A more detailed model which the aircraft is treated as a rigid body, with rotational inertia, will be considered later on in the course. Although approximate, the current model is found to be useful for many applications. In addition ve will assume that the aircraft fies in a vertical plane. The extension to three dimensions, although not hard will make things too complicated at this point. Therefore, the two degrees of freedom of the aircraft are the horizontal and vertical positions, T and y. In order to be able to write the equations of motion we eed to look at the forces that act on the aircraft
other hand, we will see that the smaller the timestep the more time steps we will have to compute). The idea is to start from the initial condition, and compute the solution at times ∆t, 2∆t, 3∆t, etc. Suppose that we have already progressed along the trajectory to time t, when the components of X(t) are all known (this could clearly be t = 0 because of the initial condition). We then update all these quantities to the next discrete time t + ∆t, by X(t + ∆t) = X(t) + ∆tF(X(t), t) Note Solving ODE’s with MATLABr The task of integrating a systems of ODE’s is particluarly simple if we use a computer package such as MATLAB. There exist a number of solution algorithms that one can choose from, but for our purposes almost any of them will work. The MATLAB syntax is typically of the form [t,X] = solver(odefun,tspan,x0) . The returned arrays t and X contain the vector of time points and the solution, X at these time points, respectively. A function, odefun, that calculates F, from X, and t, must be provided, as well as a vector specifying the interval of integration tspan = [t0,tf], and the initial value of X, x0. In order to set the parameters for the solution algorithm one can use the MATLAB function odeset. See http://www.mathworks.com/access/helpdesk/help/techdoc/ref/ode45.html for details as well as the lab 1 handout. 2 Degrees-of-Freedom (2DOF) Aircraft Model In this section we will develop a simple model to simulate the flight of an aircraft (or spacecraft). The aircraft is treated as a particle. Therefore, we assume that all its mass is concentrated at a point located at the aircraft’s center of mass. In addition, we will assume that we have perfect attitude control. That is, we will decouple the motion of the center of mass with the rotation of the aircraft about the center of mass. For the later, we will assume that we have a control system that allows us to prescribe the angle of attach of the aircraft. In reality this is clearly not possible and any moment applied to the aircraft, through the control surfaces, to change its attitude will take a certain time to be effective. A more detailed model in which the aircraft is treated as a rigid body, with rotational inertia, will be considered later on in the course. Although approximate, the current model is found to be useful for many applications. In addition, we will assume that the aircraft flies in a vertical plane. The extension to three dimensions, although not hard will make things too complicated at this point. Therefore, the two degrees of freedom of the aircraft are the horizontal and vertical positions, x and y. In order to be able to write the equations of motion we need to look at the forces that act on the aircraft. 3
With reference to the figure below, we define the following angles. Let B, be the angle that the velocity vector(tangent to the fight path) makes with the x-axis. Thus, we will have tan B=vu/ur. Let a be the ngle of attack of the aircraft relative to velocity vector. As we will see, the angle of attack is relevant, as it determines the lift(and drag) for a given velocity and altitude. The angle of attack, is defined so that the lift generated when a= 0, is equal to zero. In addition, we define the angle ar as the angle that the thrust ector T makes with the zero-lift line. In most cases, this angle will be constant, but in some situations (e. g thrust vectoring in aircraft or rockets)it can be an input, or control, variable. All the angles are positive as thrust lin tero-lift line flight path X-axis The aerodynamic forces, are resolved into a lift component L, orthogonal to the velocity, and a drag compo- nent D, parallel to the velocity. In addition, we will have, the weight of the vehicle, w, and the thrust t, due to the propulsive system. When these forces are referred to the center of mass of the aircraft, we may have an additional pitching moment M. For this aircraft model, we will assume that the control surface guarantees that the pitching moment about the center of mass is always zero Forces on the aircraft For a standard wing-tail configuration aircraft we have the following Lift Most of the aircraft's lift is generated by the wing and horizontal stabilizer. For a given angle of attack a the lift is proportional to the relative wind s dynamic pressure q =(1/2)pu, and to the wings area S, L=louse
With reference to the figure below, we define the following angles. Let β, be the angle that the velocity vector (tangent to the flight path) makes with the x-axis. Thus, we will have tan β = vy/vx. Let α be the angle of attack of the aircraft relative to velocity vector. As we will see, the angle of attack is relevant, as it determines the lift (and drag) for a given velocity and altitude. The angle of attack, is defined so that the lift generated when α = 0, is equal to zero. In addition, we define the angle αT as the angle that the thrust vector T makes with the zero-lift line. In most cases, this angle will be constant, but in some situations (e.g. thrust vectoring in aircraft or rockets) it can be an input, or control, variable. All the angles are positive as shown in the diagram. CM x-axis flight path zero-lift line thrust line T β y x v The aerodynamic forces, are resolved into a lift component L, orthogonal to the velocity, and a drag component D, parallel to the velocity. In addition, we will have, the weight of the vehicle, W, and the thrust T, due to the propulsive system. When these forces are referred to the center of mass of the aircraft, we may have an additional pitching moment M. For this aircraft model, we will assume that the control surface guarantees that the pitching moment about the center of mass is always zero. Forces on the Aircraft For a standard wing-tail configuration aircraft we have the following: Lift Most of the aircraft’s lift is generated by the wing and horizontal stabilizer. For a given angle of attack α, the lift is proportional to the relative wind’s dynamic pressure, q = (1/2)ρv2 , and to the wing’s area S, L = 1 2 ρv2SCL. (3) 4
Here, CL, is the non-dimensional lift coefficient, which is a function of the shape, and the angle of attack a. The variation of Cl with the angle of attack is linear, up to a maximum value Clma, of the order of 1.5-2.0. Thus, we have where a= aCL/aa. For an idealized fat plate wing, a= 2T rad, and more generally for finite wings a=5 rad- to a=6 rad-I. Note that for a=0, CL=0, which implies that we are assuming that a=0 corresponds to the zero-lift line. Pitching Moment The lift force, L, is the result of integrating a distributed force over the surface of the wing. This distributed force is equivalent to L, acting through a particular point called the center of pressure(CP). For this model, re will assume that through proper control of the flaps and tail, the pitching moment about the center of mass is zero. Moreover, we shall assume that any attitude changes commanded by the control system are instantaneous The drag force can als exPress ed as a function of a non-dimensional drag coefficient, CD (5) The drag coefficient in general has a complicated dependence on a and Cl. But commonly it is assumed to have two components, one which is independent of the lift, and one which increases quadratically with the Cp=Cpo +KCi Typical values for CDo for the whole aircraft, are of the order of 0.003 to 0.02. An important ratio used to determine the aircraft's performance is L/D= CL/CD, which is typically of the order of 10-20 for most subsonic aircraft, and up to 60 for some sailplanes Thrust The thrust generated by the engine is a force directed at an angle ar measured with respect to the zero-lift line. In general, both the magnitude of T, and its direction ar, are inputs to the model. Weight Typically we will simple have W=mg and it will be directed along the negative y-axis. Note, however, that we are interested in spacecraft in orbit, the weight might be directed to a fixed point, e. g. the earth center and may be a function of the height
Here, CL, is the non-dimensional lift coefficient, which is a function of the shape, and the angle of attack α. The variation of CL with the angle of attack is linear, up to a maximum value CLmax , of the order of 1.5 − 2.0. Thus, we have CL = aα , (4) where a = ∂CL/∂α. For an idealized flat plate wing, a = 2π rad−1 , and more generally for finite wings a = 5 rad−1 to a = 6 rad−1 . Note that for α = 0, CL = 0, which implies that we are assuming that α = 0 corresponds to the zero-lift line. Pitching Moment The lift force, L, is the result of integrating a distributed force over the surface of the wing. This distributed force is equivalent to L, acting through a particular point called the center of pressure (CP). For this model, we will assume that through proper control of the flaps and tail, the pitching moment about the center of mass is zero. Moreover, we shall assume that any attitude changes commanded by the control system are instantaneous. Drag The drag force can also be expressed as a function of a non-dimensional drag coefficient, CD, as D = 1 2 ρv2 S CD. (5) The drag coefficient in general has a complicated dependence on α and CL. But commonly it is assumed to have two components, one which is independent of the lift, and one which increases quadratically with the lift, CD = CD0 + KC2 L . (6) Typical values for CD0 for the whole aircraft, are of the order of 0.003 to 0.02. An important ratio used to determine the aircraft’s performance is L/D = CL/CD, which is typically of the order of 10 − 20 for most subsonic aircraft, and up to 60 for some sailplanes. Thrust The thrust generated by the engine is a force directed at an angle αT measured with respect to the zero-lift line. In general, both the magnitude of T , and its direction αT , are inputs to the model. Weight Typically we will simple have W = mg and it will be directed along the negative y-axis. Note, however, that we are interested in spacecraft in orbit, the weight might be directed to a fixed point, e.g. the earth center, and may be a function of the height. 5
Equation of Motion Written in cartesian coordinates, the equations of motion are mar= Tcos(B +a+ar)-Dcos B-Lsin B may=Lcos B+Tsin(B +a+ar)-Dsin B-w State Space Form The equations of motion can be written in the form of equation(2)by introducing the velocities as variables (T/m)cos(B +a+aT)-(D/m)cosB-(L/m)sin B (L/m)cos B+(T/m)sin(B+a+ar)-(D/m)sin B-(w/m) where, T, a and ar are inputs to the model which will typically be a function of time and perhaps x through the control system. Recall that B= tan(ay/uz), and L is by expressions(3),(4)and D, by expressions (5),(6),and(4) ADDITIONAL READING J. L. Meriam and L.G. Kraige, Engineering Mechanics, DYNAMICS, 5th Edition 3/3, 3/4, 3/5(rectangular coordinates only)
Equation of Motion Written in cartesian coordinates, the equations of motion are max = T cos(β + α + αT ) − D cos β − L sin β may = L cos β + T sin(β + α + αT ) − D sin β − W (7) State Space Form The equations of motion can be written in the form of equation (2) by introducing the velocities as variables and taking X = x y vx vy , F = vx vy (T /m) cos(β + α + αT ) − (D/m) cos β − (L/m) sin β (L/m) cos β + (T /m) sin(β + α + αT ) − (D/m) sin β − (W/m) (8) where, T , α and αT are inputs to the model which will typically be a function of time and perhaps X through the control system. Recall that β = tan−1 (vy/vx), and L is by expressions (3), (4) and D, by expressions (5), (6), and (4). ADDITIONAL READING J.L. Meriam and L.G. Kraige, Engineering Mechanics, DYNAMICS, 5th Edition 3/3, 3/4, 3/5 (rectangular coordinates only) 6
