4. Linear Systems Modeling

Linear systems modeling is often used in engineering because this simplification makes it possible to easily analyze the properties of the system.
In particular, linear modeling is most suitable for:

1) analyzing the stability of a system
2) vibrations mechanics
3) developing control systems

> restart:  mylib :="C:\LIBS": libname := mylib,libname: with(MBsymba):

4.1 The Standard Procedure for Developing a Linear Model

In order to develop a linear model of a system, the following steps are usually followed:

1) developing the full mathematical model of the system, i.e. finding the set of first order differential equations which describe the system

(Note a second order diff. eqn. may be split into two first order diff. eqns.)

where it is supposed that only external actions depends on time

2) finding the equilibrium solution(s) in absence of external actions and with zero velocity (i.e. time dependency disappears)

3) linearizing the equations of motion as follows (Taylor series expansion)

Linear Models of Mechanical Systems

For a mechanical system the linear equations of motion may be written in the following form

where matrices are the well known mass, damping, and stiffness matrices.

Each of them is associated to a quadratic, respectively the kinetic energy, the Rayleigh's dissipation function, and the elastic energy, as follows:

the expressions above are widely used in modeling linear systems using the Lagrangian approach.

4.2 A Smart Approach to Linear Modeling

As seen above, the linearization of the equations of motion is the last step in mathematical modeling.
This implies that the derivation of the non-linear equations of motion is still required and it may be not easy.
Indeed, many difficulties arise when modeling mechanical systems due to the presence of non-linear rotational terms in kinematics and so on. Therefore, it would be very useful to be able to develop the linear equations of motion directly instead of the non-linear ones. In general this is not possible but there is a large class of problems which allow this operation. In particular, this is possible when the equilibrium solution is clear at the begin of modeling process and all variables represents small variation from the equilibrium configuration (It is important that all variables are small, otherwive there will be some trouble that will be discussed later).

Under these assumptions, it is possible to use linearized kinematic relations for computing the expression of kinetic and potential energy, and in conclusion for deriving the equations of motion without losing any information of the system properties.

Indeed, starting from the non-linear expression of the position, velocity, and acceleration of a body

the linearization leads to

and it should be noted that the linear expressions of velocity and accelerations may be obtained either by linearizing the corresponding non-linear expressions or by deriving the linear position and then differentiating with respect to time.

For quadratric expressions, such kinetic energy, the 2nd order Taylor expansion yelds:

and once again quadratic expressions may be calculate by direct use of the linear kinematic equations.

Therefore, the linear equations of motion may be derived using either the Newtonial or Lagrangian approach.

4.3 Example: the Pendulum

Definition of the the reference coordinate frame of the pendulum located a the center of mass of the pendulum bob.

PF := rotate('Y',alpha(t)) * translate(0,0,L);

Position of the center of mass.

> G := origin(PF): show(%);

Linearized for small angles of pendulum displacement.

GL := linearize(project(G,ground), [alpha(t)],t): show(%);

velocity: first derivation and then linearization

> VG := velocity(G): show(%);
linearize(VG, [alpha(t)],t): show(%);

velocity: first linearization and then derivation

> VGL := velocity(GL): show(%);

FOR DETERMINING THE VELOCITY IT CAN BE SHOWN THAT ANY ORDER OF LINEARIZATION AND DERIVATION YIELDS THE SAME RESULT!

acceleration: first derivation and then linearization

> AG := velocity(VG): show(%);
linearize(AG, [alpha(t)],t): show(%);

acceleration: first linearization and then derivation

> AGL := velocity(VGL): show(%);

FOR DETERMINING THE ACCELERATION IT CAN BE SHOWN THAT ANY ORDER OF LINEARIZATION AND DERIVATION YIELDS THE SAME RESULT!

kinetic energy: first derivation and then Taylor expansion

> KE := 1/2*m*dot_prod(VG,VG);
taylorF(KE,[alpha(t),diff(alpha(t),t)],3);

kinetic energy: first linearization and then calculation

1/2*m*dot_prod(VGL,VGL);

FOR DETERMINING THE KINETIC ENERGY IT CAN BE SHOWN THAT ANY ORDER OF LINEARIZATION AND DERIVATION YIELDS THE SAME RESULT!

4.4 Example: Sequence of Rotations

An important property in linear kinematics is that rotational transformations are commutative
Indeed:

> T1 := rotate('X',alpha) * rotate('Y',beta);

> taylorF(T1,[alpha,beta],2);

> T2 := rotate('Y',beta) * rotate('X',alpha);

> taylorF(T2,[alpha,beta],2);

4.5 Example: the Pendulum Mounted on a Rotating Support

kinematics

support rotating frame

> RF := rotate('Z',Omega*t);

pendulum frame

> alpha0 := 'alpha0'; equilibrium position

> sm_vars := [alpha(t)]; small displacements from equilibrium

> PF := RF * rotate('Y',alpha(t)+alpha0) * translate(0,0,L);

> G := project(origin(PF),ground): #show(%);
GL := linearize(G, sm_vars,t): show(%,'synthetic');

velocity:
first derivation and then linarization

> VG := velocity(G): #show(%);
VGL := linearize(VG, sm_vars,t): show(%);

position
of the pendulum Center of Mass

> G := project(origin(PF),RF): #show(%);
GL := linearize(G, sm_vars,t): show(%,'synthetic');

velocity:
first derivation and then linearization

> VG := project(velocity(G),RF): #show(%);
linearize(VG, sm_vars,t): show(%);

first linearization and then derivation

> VGL := project(velocity(GL),RF): show(%);

acceleration:
first derivation and then linearization

> AG := project(velocity(VG),RF): #show(%);
linearize(AG, sm_vars,t): show(%);

first linearization and then derivation

> AGL := project(velocity(VGL),RF): show(%);

Euler's equations

> pendulum := make_BODY(G,m): show(pendulum);
hub_torque := make_TORQUE(PF,MX,0,MZ,pendulum): show(hub_torque);
_gravity := make_VECTOR(ground,0,0,g): show(_gravity);

first derive equation, then linearize

> eqnE := euler_equations({pendulum,hub_torque},origin(ground),PF): #show(%);
eqnE := linearize(eqnE, sm_vars,t): show(eqnE,'synthetic');

first linearize then derive equation

> pendulumL := make_BODY(linearize(G,sm_vars,t),m): #show(%);

> hub_torqueL := make_TORQUE(linearize(PF,sm_vars,t),MX,0,MZ,pendulumL): #show(%);

> eqnEL := linearize(euler_equations({pendulumL,hub_torqueL},origin(ground),PF),sm_vars,t):
eqnEL[comps] := map(simplify,eqnEL[comps],trig):
show(eqnEL,'synthetic');

equilibrium solution

> Yeqn := collect(simplify(expand(Ycomp(eqnE)),trig),[alpha(t)]):

> subs(alpha(t)=0,expand(Yeqn/sin(alpha0)));

> equilibrium := alpha0 = solve(%,alpha0);

the equation of motion becomes:

> expand(subs(equilibrium, Yeqn));

4.6 Example: a Rigid Rotor Mounted on Flexible Bearings

> alpha(t), beta(t), rotations due to bearings elastiticy

> Omega; constant spin velocity

we assume that rotations due to elasticity are smalls
using the 'l inear modelind({small_vars}) ' command, linearization will be made automatically

> linear_modeling( {alpha(t),beta(t)} );

Warning, Linear Modeling option has been choosen for the following variables:

body frame

> BF := rotate('Y',alpha(t)) * rotate ('X',beta(t)) * rotate('Z',Omega*t);

> project(angular_velocity(BF),ground);

definition of rotor mass properties

> rotor := make_BODY(BF, m, It,It,Ia): show(%);

> k, linear stiffness of the bearing

> L; distance between beerigns

> MX := -k*(L/2)^2*alpha(t); elastic torque

> MY := -k*(L/2)^2*beta(t); elastic torque

reative torque due to the bearings

> elastic_bearings := make_TORQUE(ground, MX,MY,0, rotor): show(%);

equations of motion using Euler's Equation

> eqns := euler_equations({rotor,elastic_bearings},CoM(rotor),ground):
show(%);

>