Computer-aided generation of equations and structural diagrams for simulation of linear stationary mechanical dynamic systems

V. K. Augustaitis*, V. Gičan**, N. Šešok***, I. Iljin**** *Vilnius Gediminas Technical University, Basanaviciaus 28, 03224 Vilnius, Lithuania, E-mail: vytautas.augustaitis@vgtu.lt **Vilnius Gediminas Technical University, Basanaviciaus 28, 03224 Vilnius, Lithuania, E-mail: vladimir.gican@vgtu.lt ***Vilnius Gediminas Technical University, Basanaviciaus 28, 03224 Vilnius, Lithuania, E-mail: pgses@vgtu.lt ****Vilnius Gediminas Technical University, Basanaviciaus 28, 03224 Vilnius, Lithuania, E-mail: pgilj@vgtu.lt

In majority of cases, the application of computeraided simulation in the environment of MATLAB / Simulink programs for such investigation is purposeful [13,14].For this purpose, the equations for description of the dynamic processes of the equipment under investigation, hereinafter referred to as the dynamic system or shortly system, and the structural diagram of the system formed on the base of the said equations according to the dynamic model of the system with discrete elements (widely used for examination of automatic control systems) should be available.Such a diagram formed according to the requirements of MATLAB/Simulink program package and "understandable" for a computer is referred to as Simulinkmodel [13,14].There are no substantial differences between it and the structural diagram: if any of them is available, the other is easily found on its base.When a mechatronic system is examined, its electrical part, such as electric drive, components of sensors and so on, usually are provided as already known structural diagrams [15].Cases of formation of structural diagrams of the mechanical part of mechanical or mechatronic systems appear to be more complicated.Because of the wide variety of structures of mechanical systems, their structural diagrams are not predictable in any specific case, so their formation in cases of complicated systems requires considerable attempts and errors are hardly avoidable if the process is not computerized.
The purpose of the paper is provision of a methodology of the application of MATLAB/Simulink program package to the available linearized mechanical system (or is a part of more complicated system) for computer-aided generation of equations describing movement of the system using Lagrange equations of the second type, transformation of the generated equations into the convenient structural diagram of the said system and Simulink-model, formation of literal and digital analytic expressions of transfer functions included in them.
Software package MATLAB Simulink program is not directly intended for creation of the said model.However, there are enough resources to solve this problem, which is devoted to this work.
The obtained structural diagram can be easily integrated in a structural diagram and Simulink-model of a more complicated system.In addition, the generated equations of the system are of independent value as well.

Generation of the equations
The equations describing movement of the linear mechanical dynamic system with lumped parameters required for the formation of structural diagram of the system are obtained from Lagrange equations of the second type ( ) where Π , T are kinetic and potential energy of the system, respectively; Φ dissipation function; i-th generalized Lagrangian coordinate, shortly referred to as the coordinate, and its time (t) derivative; the generalized outside force acting along the coordinate .
It can be stated that each of coordinates corresponds to one equation of (1).In traditional case, the equation (1) include n generalized coordinates x i x i ; their number n equals to the number m of degrees of freedom of the system (n = m) when a structural diagram for of selfcontained mechanical system is generated of transformed Eq. (1).However, on generation of structural diagrams in more common case when the stationary linear mechanical system under discussion is a part (component) of a larger mechanical system with nonlinear or nonstationary components, a mechatronic system and so on, it can be m > n, i.e. we'll have n equations with m coordinates included in them and the structural diagram for such a system.Hereinafter, we'll mark such coordinates by x s (s = 1, ..., i, ..., m).In such a case, the Eq. ( 2) will include k = m -n redundant coordinates.For the simulation of a system described by equations with redundant coordinates, the values of such redundant coordinates are set or found from the "rejected" part of the system by connecting its structural diagram to the structural diagram of the mechanical part under discussion.
It is notable that in both cases, i.e. when m = n and when m > n, the methodology of formation of structural diagrams and Simulink-models remains the same.Incorporating of auxiliary coordinates in the model let us assume analyzed system as a part of another more complex system, which may have nonlinear and nonstationary elements, electrical elements and so on.http://dx.doi.org/10.5755/j01.mech.17.3.500 In the case under discussion, kinetic energy T is a function of derivatives of coordinates and (more rarely) a function of the coordinates themselves; potential energy is a function of coordinates , and dissipation function Ф -a function of derivatives .In addition, the expressions of can be direct functions of time t, when, for example, known kinematic excitations ( ) where ; are the totalities of coordinates and their derivatives .
In many cases [16,17], striving to facilitate the generation of analytic expressions of the functions for complicated systems, it is purposeful to use not only the generalized coordinates included in the Eq. ( 1), hereinafter referred to as the principal generalized coordinates, or shortly -to as the principal coordinates, but also auxiliary coordinates s x j δ .They should be chosen in accordance with the below condition, i.e. they should be expressed unambiguously by the principal coordinates according to linear dependence in the following equa-  , , , , ; , , Before differentiation of the functions Φ Π , , T included in the Eq. ( 1), the auxiliary coordinates and their derivatives should be eliminated by using the Eqs.(3) and (4), i.e.
If the auxiliary coordinates are properly chosen according to the Eq.(3), the analytic Eq. ( 5) of the functions Φ Π , , T developed by the investigator will be considerably simpler, as compared to Eq. ( 2).Then elimination of the auxiliary coordinates and their derivatives as well as differentiation of the functions Eq. ( 2) shall be carried out in a computer-aided way.
When the final Eq. ( 2) of the functions F , , , T Φ Π after differentiation and other relevant procedures are available, linear differential equations with constant coefficients of a degree not higher than second (in some cases, linear algebraic equations) are obtained for description of the dynamic model of the system under discussion.As it was mentioned, each i-th generalized coordinate corresponds to i-th equation from the Eq.(1).
where ( ) Eqs. ( 2) and ( 5) but not assessed in the Eq. ( 3).The expressions of the generalized forces F i (t) are developed by the investigator, so they are not discussed in details herein.
For the formation of structural diagram, each Eq. ( 6) is solved in respect of the coordinate x i included in it [ .
According to the equations (11), the structural diagrams corresponding to them are formed; they are connected into the complete typical structural diagram of the whole system shown in Fig. 1 and on the base of the latter, Simulink-model for simulation of the system is developed.
Fig. 1 The typical structural diagram of the system under discussion formed on the base of the Eq.(11)

The example
Let's suppose that low frequency vibrations of a vehicle moving on an uneven road are examined (Fig. 2).The vehicle is considered as absolutely solid body standing on elastic supports with damping; the supports simulate units of its wheels.The masses of the wheels and the suspension are neglected.For generation of equations for de-scribing movement of the vehicle, two systems of coordinates are chosen.One of them, i.e. the system is fixedly connected to the vehicle and moves together with it.The point of the origin of coordinates of the said system of coordinates coincides with the center of stiffness of the vehicle located in the plane of fixing four elastic supports to the vehicle.In the plane, the points of fixation of the supports 1, 2, 3, and 4 are situated.All these points are situated in the same distance from the point of origin of coordinates and are symmetrical in respect of it (their positions are defined by the distances and , respectively).The other (immovable) system of coordinates , upon the rest of the vehicle when it is affected by the force of gravity only coincides with the system (in such a state, both systems of coordinates are shown in Fig. 2).It is conditionally supposed that the velocity of the vehicle's movement equals to zero and the uneven road is "moving" with its velocity It is considered that vibrations of the vehicle in respect of the system of coordinates are small.During the vibrations, it rotates by small angles about the axes passing the point and parallel to the axes X and Y ; in addition, it moves in vertical direction to the distance z along the axis of coordinates Z.For simplifications of the example, the shifts along the axes X, Y and rotation about the axis Z are neglected.It is considered that coordinates of the center S of the mass of the vehicle in the system of coordinates are .On simulation of the wheels and units of their suspensions, all four elastic supports of the vehicle are considered alike (their coefficients of stiffness and resistance are k and h, respectively).It was supposed that vibrations of the vehicle are excited in a kinematic way by the road's inequalities and ( ) that are known functions of time t (kinematic excitation).In addition, excitation of vibrations of the vehicle by the forces and moments of forces that impact the engine is assessed as well; after reduction to the origin of coordinates of the system the said forces and moments of forces are expressed by the vertical force ( ) Thus, the dynamic model discussed upon in the example has three degrees of freedom and its movement is defined by the principal generalized coordinates z , , y x ϕ ϕ .
The following specific values of the parameters of the dynamic model presented in the example were accepted: m, J x , J y , J x,y Fig. 2 The dynamic model of a vehicle moving on an uneven road Kinetic energy of the machine according to [18] without auxiliary coordinates For the simplification of the initial expression of potential energy and dissipation function Π Φ , auxiliary coordinates are introduced.Let's suppose that on vibrations of the vehicle, the elastic elements used for simulation of the wheels and their suspensions are deformed by the values that are considered auxiliary generalized coordinates.So, potential energy The auxiliary coordinates are expressed by the principal ones using the Eq.(3) After insertion of the values of the auxiliary coordinates j δ in the Eq. ( 13), potential energy expressed by principal coordinates only is found The expression of dissipative function is found in an analogous way

Structural diagram of the dynamic model of the example
The generalized outside forces ( ) Then, using Lagrange Eq. ( 1), the equations of the vehicle's movement in principal generalized coordinates z , , y x ϕ ϕ (there are no auxiliary coordinates in the example under discussion) are generated.The expressions of Eqs. ( 12), ( 15) and ( 16) included in the Eq. ( 1) are differentiated and if the generalized outside forces Eq. ( 17) are known, the differential equations of the vehicle's vibrations are found; from the latter, Eq. ( 11) type are obtained The structural diagram shown in the Fig. 3 is formed on the base of the solutions of ( 18) and the typical structural diagram (Fig. 1).

Computerized realization of the proposed methodology
The application of the computerized methodology of computation of Eq. ( 11) type equations describing mechanical dynamic systems and the parameters of Simulinkmodel formed on their base in literal and digital form can be divided to the following phases.
1. Familiarization with the system, formation of its dynamic model, computation of the values of its parameters, choosing the principal generalized coordinates, choosing the auxiliary generalized coordinates (if they are used) generation of their Eqs.(3) and ( 4) in the principal generalized coordinates, generation of expressions of kinetic energy T, potential energy П and dissipative function Ф included in the Eq.1), formation of expressions of the generalized outside forces F i , kinematic excitations and that affect the system.All materials are developed by the investigator.
Entering the initial data in to a computer.The data to be entered: data names; control matrix K; names of auxiliary coordinates Eqs.(3) and, (4), if they are used, and their time t derivatives; analytic expressions of the functions Т, П, Ф, numerical values of the parameters of the dynamic model.
First of all, data names are entered.They are written in lines; each line starts from the word "syms".The order of entering is chosen freely.The entered data names are distributed in a line, for example, as follows.
In the first line, all principal generalized coordinates , the kinematic excitations (if they are used) in the expressions of Т, П, Ф and nonstationary members In the second line, the derivatives of all values listed in the first line are provided; they are marked with the letter D before any of such value, for example, Dz, Dx and so on.
In the third line, the names of all parameters included in the dynamic model of the system under discussion (all constants included in the expressions of the functions Т, П, Ф, the connection Eq. (3), kinematic and outside excitations and generalized forces) are provided.
In the fourth line, the names of generalized outside forces F i (t), the differentiation operator p, and time t included in the equation (1) are provided.
For the example under discussion, the data names are written as follows syms phix phiy Z z1 z2 z3 z4 % syms Dphix Dphiy DZ Dz1 Dz2 Dz3 Dz4 syms m Jx Jy Jxy xs ys l1 l2 k h syms Mx My Pz p t (20) The control matrix K is entered in four lines.Each line should include the same number of elements.
In the first line, the names of the principal generalized coordinates , kinematic excitations and nonstationary members are entered in any order.
In the second line, time t derivatives of all the values listed in the first line are provided; they are marked with the letter D before any of such value.Under the name of each value listed in the first line, the name of its derivative should be specified in the second line, i.e. the name of a value mentioned in the first line and the name of its derivative should be in the same column of the matrix K.
In the third line, names of the generalized outside forces F i (t) acting along the generalized coordinates x i shall be written.The said names shall be written in those columns with coordinate's x i of the matrix K that correspond to the directions of acting of the forces.Other elements of the third line are equal to zero.
In the fourth line, it is specified that the functions Т, П, Ф included in Lagrange equations ( 1) should be differentiated according to the coordinates and their derivatives.In the elements of the line situated in the columns of the names of the mentioned coordinates and their derivatives, any positive whole numbers, such as 1, is written.Other elements of the line are equal to zero.
After the matrix K, the analytic (literal) expressions of auxiliary coordinates j δ and their derivatives j δ are entered.They are entered in lines according to the order specified in the Eqs.(3) and (4).
For the example under discussion Then the analytic expressions of the functions Т, П, Ф are entered (on entering, they are marked as TM, PM, FM, respectively).They consist of the sum of the sum-monds of the line [1].The summonds of each of said functions with the auxiliary coordinates (when they are used) are entered as summarized elements of the column of the vector.
For the example under discussion The Program Functions carries out the following actions: develops the analytic (literal) expressions of the functions Т, П, Ф (2) without auxiliary coordinates; differentiates them according to coordinates , according their derivatives and time t (Eq.( 1)); transforms them into the operator form and solves in respect of the coordinates , i.e. generates the Eq. ( 11); forms the matrices Mark, Den and Denc (Table ) with all the data required for formation of the structural diagram and Simulink-model of the system.The four-column matrix Mark consists of n groups of lines.Each i-th group of lines (i = 1, 2, ..., n) conforms to the x i -th generalized coordinate in respect of which the i-th Eq. ( 11) is solved.The number of lines of a group is where k i -the number of the coordinates in the right part of the i -th Eq. ( 11), i -the total number of generalized output forces F i and kinematic excitations N i,e in the same Eq.( 11) ).The total number of all lines of the matrix Mark is: . Groups of the lines are situated in this matrix in the same order as the coordinates x i in the first line of the matrix K, i.e. the group of lines that corresponds to the coordinate x 1 is in the beginning of the matrix, then the group of lines that corresponds to the coordinate x 2 follows and so on.
The first element of a line of any group of lines is the number of the line in the matrix (the lines are numbered consecutively starting from the number 1 and ending by the number ψ ).Other three elements of the line depend on the group the line belongs to and on its place in the group.
Let's suppose that we have the i-th group of lines.All third elements of this group of lines (the third column of the group) are the serial number of the coordinate x i in the matrix K, i.e. the number "i" ( i = 1, 2, ..., n ).The second element of the first line of the group is the name of the coordinate x i and the fourth element is symbolic inscription "den" showing that the line with such inscription is to be used as a denominator of the coordinate x i .
The other lines of the group are usable for marking the names of the values included in the numerator of the right part of the Eq.(11).The fourth element of all these lines is the symbol "num" showing that the value mentioned in the line is in the numerator of the right part of the Eq.(11).The second element of the second line of the group is the name of the generalized force F i .Starting from the third line, total k i lines are used for the names of the coordinates included in the right part of the i-th Eq. (11).The second elements of all said lines are names of coordinates; the lines that correspond to the said coordinates are situated in the same order as the names of the said coordinates in the matrix K.The lines of the last, i.e. ith group (the number of them is -Eq.( 9)) are used for listing the kinematic excitations N i v i,e that impact the system.Their names are the second elements of the said lines.
In matrix Den, the analytic (literal) expressions of the coefficients are provided; their dislocation in this matrix is coordinated with dislocation of the elements of the second column of the matrix Mark.The said coefficients are polynomials of the second degree in respect of the operator p according to the Eqs.( 7), ( 8) and (10) and are described by the following matrices lines  The name of the coordinate x s , force F i or kinematic excitation N i,e written in the second elements (the elements of the second column) of lines of the matrix Mark corresponds to the one of the lines Eqs. (25) or (26) of the matrix Den (that includes one of the said names) dislocated according to the same order as the lines of the matrix Mark.
Thus, both the matrix Den and the matrix Mark consist of n groups of lines.Each group includes r i lines and corresponds to one of the coordinates x i (i = 1, ..., n).The first line of the i-th group of the matrix Mark (its second element is the name of the coordinate x i ) corresponds to the line d i,i of the same group of the matrix Den.The second line includes the force F i and it corresponds to the line of the matrix Den.Analogously, one of lines Eq. (25) of the matrix Den corresponds to relevant other lines of the i-th group of the matrix Mark.When the data of the matrices Mark, Den, Denc are available, structural diagram of the system under discussion is developed in accordance with the instructions of MATLAB / Simulink set of programs and its simulation is carried out.In addition, the data provided in the matrices are sufficient for generation of the equations of the type Eq. ( 11) for the systems under discussion that is required for the application of various other methods of computation and research.
In addition to the data provided in the Table, pictograms of blocks of transfer functions that present a basis of the structure of Simulink-model of the object under discussion are developed.
For convenience of simulation, the Simulinkmodel of the whole object is divided to subsystems that correspond to relevant coordinates .In Fig. 4, a, the subsystem that corresponds to the coordinate The name Trans is automatically classified to the blocks of transfer functions; after it, a fraction follows; the numerator and the denominator of the fraction specify the numbers of the lines of the matrices Mark, Den, Denc used in the block (Table ).
Blocks of transfer functions are provided in a compact form.
The pictogram of block Trans 2/1 in Fig. 4, a is enlarged to show more clearly the mark Num(s)/Den(s) (inside its contour) showing that the numerator and the denominator of the transfer functions are lines of the matrix.
The fraction following the word Trans, for example, 2/1 (Fig. 4, a The connecting lines of the whole model and its subsystems in Simulink-models are connected by the investigator using the graphical editor of Simulink set.The fully connected system of the example is shown in Fig. 5, a.The simulated transitional process obtained on an abrupt change of the road unevenness is shown in Fig. 5, b. 3.In course of generation of the structural diagrams, data of independent value required for the generation of differential equations of the second degree for the system under discussion and development of the normal form of the said equations are obtained.

COMPUTER -AIDED GENERATION OF EQUATIONS AND STRUCTURAL DIAGRAMS FOR SIMULATION OF LINEAR STATIONARY MECHANICAL DYNAMIC SYSTEMS S u m m a r y
There are reasonable to use computer with MATLAB/Simulink software for modeling of dynamic systems of machine mechanic components, using lumpedparameter systems of them.In order to build structural schematics (Simulink blocks) of these dynamic models, there is a lot of work.This paper is intended to analyze linearized stationary mechanical system, which may be a part of more complex (nonlinear, nonstationary) system.There are presented methodology to build such Simulink model in the paper.Here are given analytical expression of differential equations and transfer functions for Simulink as well as numerical values.Simulink model is built in such way, that it is easy to attach necessary transfer functions of nonlinear or/and nonstationary connections.Example of such model is provided.

γ
included in the expressions of the auxiliary coordinates Eq. (3) are listed.
and ψ -line matrix Den consists of the lines (25); they are supplemented by the lines individual coefficients at forces F i in the Eq.(11); their analytic expressions are developed by the investigator.The elements 25) in the matrix Den are provided in the analytic (literal) form expressed in the parameters of dynamical model of the system.

Matrix
Denc differs from matrix Den only by numeral values of the coefficients provided in the Eq.(25) instead of their analytic (literal) expressions.

Fig. 4 Fig. 5
Fig. 4 The Simulink-model (a) and the window for entering the parameters of the block of the transfer function Trans 2/1 (b) for the equation that corresponds to the coordinate x ϕ of the example

1 .
The methodology of investigation of linear stationary dynamic models of the mechanical part of mechanical and mechatronic equipment in the environment MATLAB/Simulink programs required for computer-aided generation of structural diagrams (Simulink-models) is provided.2.On the base of the analytic (literal) expressions of kinetic energy, potential energy, dissipation function and generalized forces included in Lagrange equations of second type, the literal or numerical expressions of transfer functions included in the structural diagram are generated in the computer-aided way.

Table
The elements of the matrices Mark, Den and Denc found on solving the example