Modeling of Dynamic Systems in Simulation Environment MATLAB/Simulink – SimMechanics

Veiw figure View Figure

Dynamic equation of motion which applies to the mechanical system Figure 1 in its vector form

(1)

Substituting the forces equation of the motion

(2)

where m - the mass of the body,

- acceleration of the mass,

- reaction force in the spring,

- damping force,

- excitation force.

The reaction force in the spring is

(3)

The damping force is

(4)

The excitation force is

(5)

Where k - spring constant,

b - viscous damping coefficient,

x – displacement of the spring,

ω – natural frequency of the exciting force.

After substitution of relation (3), (4), (5) to the equation (2) we receive the acceleration of mechanical system with one degree of freedom as follows:

(6)

The derived equation of motion (6) is further solved in the program MATLAB. As a result we gain kinematic values (displacement, velocity and acceleration) of the mechanical system in relation to time t.

3. Solution of the Mechanical System with One Degree of Freedom in MATLAB

In MATLAB the derived non-homogenous linear differential equation of second order with constant coefficients (6) was solved in relation to mechanical system with one degree of freedom by means of Runge-Kutta method. The following parameters were taken into consideration:

- mass m=1 [kg],

- viscous damping coefficient b=0 [N/(m/s)],

- spring constant k=(2 π)² [N/m],

- exciting force F(t)= F₀.sin(ωt), F₀=1 [N],

considering three different frequencies of exciting force:

- natural frequency of exciting force ω=1,

- natural frequency of exciting force ω=π,

- natural frequency of exciting force

ω=√(k/m)=√((2π)²/1)=2π.

The following conditions are initial conditions for the solution of mechanical system:

- x(0)=0 [m],

- v(0)=0 [m/s],

with solution period 0 ≤ t ≤ 20.

Prior to solution in MATLAB it is needed to create two m-files. The first is the m-function-file defining the language of MATLAB while the second m-file shall specify initial conditions and the integration period.

The m-files have the following form:

functionTDS_F.m

function dx=function TDS_F(t,x)

m=1; k=(2*pi)^2; b=0;

F0=1;

omega=pi;

F=F0*sin(omega*t);

dx=[x(2); -k * x(1) / m - b * x(2) / m + F/m];

vibration.m

[t,x]=ode45(@functionTDS_F,[0 20],[0.0;0.0]);

% [0 20] – parameters of time simulation (start time, end time), [0;0] – initial conditions (displacement, velocity)

% ode45 – numerical solution of differential equations systems, application of Runge-Kutta method of fourth order

Figure 4

plot(t,x)

title('Kinematic variables x(t),v(t) in MATLAB: F=sin(\omega*t), \omega = \pi ');

xlabel('t [s]');

ylabel('x [m], v [m/s]');

grid on;

legend('x=x(t)','v=v(t)');

Results of displacements and velocity of the mechanical system in relation to frequencies of exciting force ω=1, ω=π, ω=√(k/m)=√((2π)²/1)=2π in relation to time t are shown in Figure 2, Figure 3 and Figure 4 by means of function plot(t, x).

Figure 2. Simulation results x(t) and velocity v(t) in relation to frequency of exciting force ω=1

Download as

Veiw figure View Figure

Figure 3. Simulation results x(t) and velocity v(t) in relation to frequency of exciting force ω= π

Download as

Veiw figure View Figure

Figure 4. Simulation results x(t) and velocity v(t) in relation to frequency of exciting force ω= 2π

Download as

Veiw figure View Figure

4. Solution of the Mechanical System with One Degree of Freedom in SimMechanics

SimMechanics is an extension to Simulink which was designed for modeling of systems with more degrees of freedom. It includes a library of blocks which constitute real parts of the mechanism (e.g. bodies and connections between them, actuators, sensors etc.). These blocks are in later stages processed on physical level. The model of actual technical object is designed when combining these blocks. Connecting physical blocks together a mathematical model is automatically designed for solution in Simulink. This method enables us to solve complex tasks in case of which manual derivation of equations of motion would be time-consuming or would cost much effort or might not be possible at all.

The calculation of kinematical values of mechanical system with one degree of freedom shall be done in SimMechanics for the same parameters of mechanical system with one degree of freedom as in case of kinematical values in MATLAB. In SimMechanics we design a block diagram in Figure 5 based on body geometry and connections regarding the environment in question. The block diagram shall further be processed in Simulink.

Figure 5. Block diagram of mechanical system with one degree of freedom in SimMechanics

Download as

Veiw figure View Figure

Simulation results of kinematical parameters of mechanical system (displacement, velocity acceleration) are shown in Figure 6 to Figure 8 by means of plot function of m-file:

% variables simout in Workspace– results parameters x(t),v(t),a(t) after processing in SimMechanics

Figure 3

set(3,'Name','results x(t), v(t), a(t) in SimMechanics')

plot(simout.time, simout.signals.values(:,1:3), 'LineWidth', 2)

title('SimMechanics:x=x(t),v(t),a=a(t): m=1kg,F0=1,k=(2*pi)^2,\omega=\pi'),

legend('x=x(t)','v=v(t)','a=a(t)'),

xlabel('t [s]'),

ylabel('x [m], v [m/s], a [m/s^2]'),

grid on

Figure 6. Simulation results of displacement x(t),velocity v(t) and acceleration a(t) in relation to frequency of exciting force ω=1

Download as

Veiw figure View Figure

Figure 7. Simulation results of displacement x(t),velocity v(t) and acceleration a(t) in relation to natural frequency of exciting force ω= π

Download as

Veiw figure View Figure

Figure 8. Simulation results of displacement x(t),velocity v(t) and acceleration a(t) in relation to natural frequency of exciting force ω=2π

Download as

Veiw figure View Figure

Simulation results of kinematical values (displacement, velocity) are shown in Figure 9, Figure 10, Figure 11 by means of plot function of m-file:

% variables simout in Workspace – results parameters x(t),v(t) after processing in SimMechanics

set(2,'Name','results x(t), and v(t) in SimMechanics')

plot(simout.time, simout.signals.values(:,1:2), 'LineWidth', 2)

title('SimMechanics:x=x(t),v(t),m=1kg,F0=1,k=(2*pi)^2, \omega=\pi'),

legend('x=x(t)','v=v(t)'),

xlabel('t [s]'),

ylabel('x [m], v [m/s]'),

grid on

As expected, the results of solution in SimMechanics in Figure 9, Figure 10, Figure 11 are equal to solution in MATLAB Figure 2, Figure 3, Figure 4.

Figure 9. Simulation results of displacement x(t) and velocity v(t) in SimMechanics – natural frequency of exciting force ω=1

Download as

Veiw figure View Figure

Figure 10. Simulation results of displacement x(t) and velocity v(t) in SimMechanics with natural frequency of exciting force ω= π

Download as

Veiw figure View Figure

Figure 11. Simulation results of displacement x(t) and velocity v(t) in SimMechanics with natural frequency of exciting force ω= 2π

Download as

Veiw figure View Figure

In the following stage the courses of displacements x(t) for frequencies of exciting force are depicted in an analogous way: ω=1, ω=π, ω=2π in Figure 12, Figure 13, Figure 14.

Figure 12. Simulation results in SimMechanics displacement x(t) with natural frequency of exciting force ω=1

Download as

Veiw figure View Figure

Figure 13. Simulation results in SimMechanics displacement x(t) with natural frequency of exciting force ω= π

Download as

Veiw figure View Figure

Figure 14. Simulation results in SimMechanics displacement x(t) with natural frequency of exciting force ω=2π

Download as

Veiw figure View Figure

5. Mechanical System with Two Degrees of Freedom

In the following part of the paper we present a model representing mechanical system with two degrees of freedom as solved in MATLAB and SimMechanics. The model consists of two bodies with masses m₁ and m₂ which are fixed to a solid frame with a spring of spring constant k₁ and a damper with linear damping coefficient b₁. The bodies are bound to each other with a spring with spring constant k₂ and a damper with linear damping coefficient b₂. Masses of the springs shall be ignored. The mass m₁ is excited by external harmonic force F(t) = F₀sin(ωt) in the axis direction of the spring and the damper, where F₀=konst. The output of the solution the model of mechanical system in Figure 15 in MATLAB and SimMechanics are simulation results of displacements, velocities and accelerations of individual system members.

Figure 15. Mechanical system with two degrees of freedom

Download as

Veiw figure View Figure

Differential second order equations of motions related to the mechanical system depicted in Figure 15

(7)

(8)

shall be processed in order to be solved in MATLAB and shall gain the form of first order differential equations through the following substitution:

(9)

(10)

(11)

(12)

while transforming second order differential equations (7) and (8) to four first order differential equations as follows:

(13)

(14)

(15)

(16)

After the above-described processing the set of equations are rewritten into the m-file of MATLAB in the following form:

function dx=functionTDS_2(t,x)

m1=75; m2=150; k1=500; k2=250;

b1=10; b2=50; F=100*sin(2*t);

dx=[x(2); ...

(-(b1 + b2) * x(2) + b2 * x(4) - (k1 + k2)*x(1)+ ...

k2*x(3) + F)/m1; ...

x(4); ...

(b2*x(2)-b2*x(4)-k2*x(3)+k2*x(1))/m2];

MATLAB includes predefined functions, which differ in relevant integration method, to solve differential equations (7) and (8) numerically. The solution of required kinematic values shall be done after the following m-file is initiated:

[t,x]=ode45('functionTDS_2',[0 10],[0; 0; 0; 0]);

Figure 1

plot(t,x,'LineWidth',1.5)

title('Simulation in Matlab:x1(t),v1(t),x2(t),v2(t)');

xlabel('t [s]');

ylabel('x1 [m], dx1/dt [m/s], x2 [m], dx2/dt [m/s]');

legend('x1','dx1/dt','x2','dx2/dt');

grid on;

Graphical representation of kinematical parameters Figure 16 related to mechanical system in question Figure 15 shall be gained through standard function ode 45.

Figure 16. Simulation results of displacements and velocities of mechanical system with two degrees of freedom in MATLAB

Download as

Veiw figure View Figure

6. Solution of the Mechanical System with Two Degrees of Freedom in SimMechanics

The solution of kinematic values of mechanical system with two degrees of freedom shall be done in SimMechanics for the same parameters of mechanical system with two degrees of freedom as in case of kinematic values in MATLAB. In SimMechanics we design a block diagram Figure 17 based on body geometry and connections regarding the environment in question. The block diagram shall further be processed in Simulink.

Figure 17. Block diagram of mechanical system with two degrees of freedom in SimMechanics

Download as

Veiw figure View Figure

Solution results of kinematical parameters for individual members of mechanical system Figure 15 have been depicted in MATLAB by means of plot function of m-file:

% Depiction of kinematic values of simout9 and simout10

Figure 2

set(2,'Name','Results in Simulink – x1,v1,a1,x2,v2,a2)

subplot(2,1,1)

plot(simout9.time, simout9.signals.values(:,1:3), …

'LineWidth', 1.5)

title('Member 1');

legend('x1 = x1(t)','v1 = v1(t) ','a1 = a1(t)]');

xlabel('t [s]');

ylabel('x1 [m], v1 [m/s], a1 [m/s^2]');

grid on;

subplot(2,1,2)

plot(simout10.time, simout10.signals.values(:,1:3), …

'LineWidth', 1.5)

title('Member 2');

legend('x2 = x2(t)','v2 = v2(t) ','a2 = a2(t)]');

xlabel('t [s]');

ylabel('x2 [m], v2 [m/s], a2 [m/s^2]');

grid on

As expected, the results of solution in SimMechanics Figure 18 are equal to solution in MATLAB Figure 16.

Figure 18. Result parameters - displacements, velocities and accelerations of mechanical system in SimMechanics

Download as

Veiw figure View Figure

The results parameters x₁=x₁(t), v₁=v₁(t) of member1 and x₂=x₂(t), v₂=v₂(t) of member2 of solution in SimMechanics Figure 19 are from m-file in Matlab where:

plot(simout9.time, simout9.signals.values(:,1:2),'LineWidth', 1.5)

and

plot(simout10.time, simout10.signals.values(:,1:2),'LineWidth', 1.5)

Figure 19. Result parameters - displacements and velocities of mechanical system with two degree of freedom

Download as

Veiw figure View Figure