CC BY
217
UDC 531.3
Dynamic Equation of Constrained Mechanical System
A. W. Beshaw
Department of Mathematics Bahir Dar University Bahir Dar, Ethiopia
This paper modifies an explicit dynamic equation of constrained mechanical system. Kinematic position of the system is defined by generalized coordinates, which are imposed on constraints. The equations of motion in the form of the Lagrange equations with undetermined multipliers are constructed based on d'Alambert-Lagrange's principle. Dynamic equations are presented to the mind, resolved relative accelerations. Expressions for the undetermined multipliers are defined by considering the possible deviations from the constraints equations. For constraints stabilization additional variables used to estimate the deviations caused by errors in the initial conditions and the use of numerical methods. For approximation of ordinary differential equations solution, in particular, the nonlinear equations of first order, use explicit numerical methods. Linear equations of the constraints perturbation are constructed. The matrix of the coefficients of these equations is selected in the process of the dynamic equations numerical solution. Stability with respect to initial deviations from the constraints equations and stabilization of the numerical solution depend on the values of the elements of this matrix. As a result values for the matrix of coefficients corresponding to the solution of the dynamics equations by the method of Euler and fourth order Runge-Kutta method are defined. Suggested method for solving the problem of stabilization is used for modeling of the disk motion on a plane without slipping.
Key words and phrases: unconstrained system, holonomic constraints, nonholonomic constraints, stabilization, Taylor series, numerical solution.
1. Introduction
Consider a discrete mechanical system of n particles Pi, i = 1, 2,... ,k of masses TOi,TO2,... ,m,k. The position of a particle in a system can be denoted by an ordered duple of scalar coefficients (x, y, z) in an inertial reference frame. The total number of displacement components in the system will be denoted ^-vector [1]:
u(t) = (ui(t),u2 (t),...,uN (t)),
where u1 = x1, u2 = y1, u3 = z1, u4 = x2, ..., un = Zk, N = 3k. When the configuration coordinates (u1,u2,... ,un) are not all independent variables, a set of reduced-order variables q = (q1,q2,...,qn) exists, where n < N, that is sufficient to define a system configuration [2]. Reduced-order coordinates are related to the configuration coordinates through displacement transformation equations such that
Ui = Ui(q 1,q2,... ,qn,t), i = 1,...,N.
Now let us consider the system as unconstrained whose configuration is described by the n generalized coordinates q = [q1,q2,... ,qn]T, here
q(to) = q°, q(to) = q°. (1)
When we say the Mechanical system is unconstrained, we mean that the components (¡i of the velocity of the system can be assigned independently at any given initial time, say t = t0. The equation of motion of the system can be obtained, using Lagrange equation
Received 26th June, 2014.
dHS) - S - = »• ' = •■ (2)
where —^ = The equation of motion of the system (2) can be rewritten by a
relation Mq = f, where M is an n x n symmetric, positive-definite, generalized mass matrix and f = f (q,q,t) is an n x 1 column array of generalized applied forces and generalized inertia force terms (including the so-called "centrifugal" and "Coriolis" terms). The generalized acceleration of the unconstrained system, which we denote it by the n-vector au = au(q, q,t), is then given by
g = M-1f = au. (3)
2. Construction of Dynamic Equations of the System
Suppose the system is subjected to m constraint equations of the form
4>(q,q,t) = 0, (4)
where ^ = ... •^m) and the constraint equations (4) include all the usual
varieties of holonomic and nonholonomic cases. The constraint equations are assumed to satisfy the initial conditions q(t0) = q0,q(t0) = q°: ^(q°, q°,to) = 0.
Under the assumption of the differentiability of constraint equations, we can differentiate equations (4) with respect to time:
ijj = Vq q + ^q q + tt. (5)
To avoid the stability problems during numerical integrations of constraints let us add terms to compensate the deviations. So that equation (5) can be rewritten as
^q <i + ^q q + ^t + g = 0, (6)
where g = g(0, q, q, t) = 0. Rearranging the coefficients of (6) along with acceleration, velocity and other terms we get
Aq = Bq + C, (7)
here A = ^q, B = q, C = -(^t + g).
Equation (7) shows the kinematical relations in connection with the constraints. Now we consider the dynamical conditions. The presence of the constraints (4) imposes additional constraint forces on the system which change its acceleration. Using Lagrange's method of undetermined multipliers the equation of motion of the constrained system can be obtained from the relation [3]
d ( dT \ dT n + ™
dU%J - % = Qi + S ■ (8)
; I r)nA ) 8n; +
k=1
From equation (8) an equivalent form can be reset [4]:
\T
Mq = f + AT /J. (9)
This is solved using the method of Lagrange multipliers whereby an additional set of m variables ¡ik are introduced. Solving for acceleration from (9) will give us [5]
9 = M -1(f + A? ¿). (10)
If the combined mass and constraint matrix is nonsingular (i.e., A has full row rank, and all particles have nonzero mass), equation (10) can be solved as
< = < + iik, (11)
where ^ is the Lagrange multiplier obtained from (10) , a%u is unconstrained acceleration of the system which is obtained from (3) and a\ is the actual acceleration of the constrained system. The relations in equation (11) are second order differential equations, we may rewrite them as systems of first order ordinary differential equations:
dql i -— = v , di. ' d^
— = a\ + D% nk.
(12)
Substitute (10) in equation (7) and solve for ^ we get
» = (AM-1AT)-1 Bq + (AM-1AT)-1 (C - AM-1 f). (13)
The dynamic equations of motion of the system can be obtained by replacing the Lagrange multiplier (13) from relation (10)
q = M-1f + M-1AT(AM-1AT)-1[Bq + (C - AM-1f)]. (14)
Equation (14) further simplified and can be written in the form:
aa = au + H(Bv + C - Aau), (15)
where H = M' ~1AT (AM-1AT )-1 is an n by m matrix.
Again relation (15) may be put in the form
aa = Pau + Sv + R, (16)
where P = I — HA, S = HB and R = HC. To solve equation (16) numerically we should change the equation from second order to first order ordinary differential equations as:
( dqi i dt '
n ' n (17)
d^ n dvl „ i ?
dif = Z ^ "df + £ Sii* + ^
i=1 3=1
The terms added in (6) to correct deviation of constraints during numerical integrations may be represented as a multiple of constraints themselves:
m
gi = ^ kii^i, i = 1, 2,...,m. (18)
=1
Using relation (6) equation (18) can be expressed as
ip = K^,
(19)
where
K =
ku ki2 ^21 &22
m1
k
■ 2
klm
k2m
km,m.
As the result of this the value of C in (7) will be replaced by C = (K^ — i^t). In turn, as its consequence, the value of R in the equation of motion (16) will be changed.
3. Stabilization of Constraints During Numerical solutions
Let the initial values q0, v0 satisfy the condition 1| < a and the system (17) is solved by Euler method [6]
<1
i+1 _ ,.i
q + rq,
+1
vi + TVi
where ql = q(ti), t = ti+1 — ti, qz = vl and
(20)
n AÏ n
E+E+ri, i = o, i,2,
3 = 1
3 = 1
Taking t sufficiently small and the inequality || < a holds for t = to, and expanding the components of ^ in powers of t in relation with (16) we get
^+1 = ^ + t^ + (2), where ^ = ^qq + q +
1+1 = r + (Pau + Sv + R) + 4q v + W T + ^ V(2),
(21)
where ^ = ^(qz,ql,ti) and ^r^(2) is the remainder of Taylor expansion. From the relation (21) and the supposition (19) it follows that
r+1 = (i + tK% + y r(2),
(22)
where I is the identity matrix. The following statements can be proposed by estimating the right-hand side of (22).
Theorem 1. In the power series expansion of ^+1, if II ^ II ^ a, 11/ + tK1 II ^
Ö < 1 and
Ç r(2)
< (1 — S)a, then II^+1jj < a.
That is, +1| < ||j + tKi|| || + £ (2)|| < Sa + (1 — 5)a = a.
The proof of this follows from (21)-(22) as it is proved in [7]. If the system (16) or (17) is solved by the fourth order Runge-Kutta method the constraint equation can be expressed as (the detail mathematical formulation is given in [8])
V+1 = V + - (k1 + 2k2 + 2kz + k4), 6
where k1 = tpl, and let tp1 = f (U,ipl),
k2=f(u+2+2f(ti,v)), k3=f(u+2+2f&+2+2f(^))),
h=f (u+r,r+rf (u+2 ,r+2 f &+2 +27 (t%, .
Expanding k2, k3, k4 up to order three using Taylor series with respect to two variables and substituting the result in (23) we get
A^ = rf * + y F* + ^ (Gi + f^ F *) + 2^ (H * + FF* + ft Gi + (4 )2F + r 5R)5, (24)
where F* = ft + ^, G% = 4 + 2^ + , P = f^ + ^^, W = 4t +
3vlflt^ + + (^)3 4^ and vl = ipz. Consider the row decomposition of 4+1 =
f (r +1,U+i):
f+1 = f + fiAr+rfi +1 f(i2) +1 f(i 3) +1 f(i4)+r 5Rk5,
3!'
4P
(25)
where
f^2) = f^ AipiAipi + 2r fa Af + r2 4,
/(*3) = 4^ Ar Ar Ar + 3r ftw Ar Ar + 3r2 Ar + r3 fttt,
f(U) = 4w ArArArAr + 4r f^
+ 6r2 few Ar Ar + 4r3 fini, Ar + r4 fittt.
Theorem 2. If a solution of (25) used the fourth order of accuracy (24) and for all values of the variables ^ = r> t = U, (i = 0,1, 2,..., N), and the matrix K(ip, t), r and the remainder t5Ri5 in Taylor series expansion satisfy the conditions:
l^ll < e, t5 (Kl! < (1 - S)e,
. T2 . T3 . T4 .
i + rK * + + yp - + 24^
< 1,
then WrW < e for all n = 1,2,...,N, where t > 0, 5 > 0, L = K + K2 P = K + 3KK + K3, Q = K + 4KK + 3K2 + 6KK + K4.
Proof. Substitute (24) in (25) and rearranging terms gives:
T2 T3
fi+1 = f + tF* + —(Gi + /¿F*) + - {W + 3IiFi + f^G* + (Ii )2F>) +
+ 2^(( fi )2p 1 + ( fi )2°l + 1 + 41*0*+
(26)
+ 6 NiFi + 74F iF + 3 (Fi)2 + M ^ + r5Rf
where N = v2 f^+2v+ftt^, M = w4 +4w3+6v2 ftt^ +4vfttt^+ftttt. Let us consider the relation ijj = K^ and respective derivatives:
$ = (K + K2)tp = Uv + ft, if = (K + 3KK + K 3)%p = f2F + 4G + 3F/ + H,
= (# + 4KK + 3K2 + 6KK2 + K 4)tp =
= f$F + + UH + 4GI + 6^F + 7 f^FI + 3 f^F2 + M.
Therefore, equation (26) can be written in the form
1-2 • T3 ••
^+1 = ^ + rKip1 + —(K + K2)ipl + — (K + 3KK + K3)^+ 26
4
+ K + 4KK + 3K2 + 6KK2 + K4W + r5^5. (27) Taking into account the given conditions and assuming that
2 3
I + TK1 + —(K + K2)* + — (K + 3KÎK + K 3)*+ 26
r4 ....
+ —( K + 4 KK + 3K2 + 6KK2 + K4V <(5 < 1, 24
t5 M.Rfc5M < (1 - S)e,
so we obtain
№+1 <
. r2 . r3 . r4 .
I + TK1 + —Ll + —Pl + —Q* 2 6 24
|| e + (1 - S)e = e.
4. Example
A disk that rolls on a plane without gliding can be considered as a system with differential coordinates. The disk shall always stand perpendicular to the xy-plane (Fig. 1).
The center of the disk is exactly above the contact point ( x, ), and the velocity of the circumference Rp of the edge equals the velocity of the contact point in the xy-plane [9] (Fig. 2):
v = rp, X = rp sin <&, y = rp cos
Figure 1. Position of the disk Figure 2. Orientation of velocity
of the disk
These nonholonomic constraint equations can be put in the form [5] ^1 = x — rp sin § = 0, = y + rp cos § = 0.
We differentiate these (a) once with respect to time:
x — rp sin § — rpD cos § = 0, q + rp cos § — rp§ sin § = 0,
(a)
or in matrix form where
Ay = b,
(b)
A =
The kinetic energy is
10 — r sin ß 0 0 1 r cos ß 0
b= [r p ß cos ß rp) ß sin ß]
T = 2mx 2 + 1M y2 + 2 hp2 + 2 hè2,
(c)
where h is the moment of inertia of the disk about the axis perpendicular to the disk through the center, and I2 is the moment about the axis through the center and the contact point (x, y).
Applying Lagrange's equations on (c) we get
(Mx = Qx + Xi, M y = Qy + X2,
hp = Qv — X^ sinß + X2r cos ß,
[l2$ = Q#.
(d)
with Qx, Qy, Qv, Q# as possible external forces in respective directions. We consider the system without such forces and therefore let them equal to zero. This transforms (d) into
' mx = X1, my = X2,
hp = —X^ sin § + X2r cos {h§ = 0.
In matrix form where
M q = ATX,
(e)
M =
'm 0 0 0 xy
0 m 0 0 y
0 0 1 0 , q = ?
0 0 0 2
X =
X1 X2
Since M is positive definite matrix, we can solve for q in (e) and write it in the form
q = M -1ATX. (f)
Substituting this (f) in relation (b) we obtain:
q = (AM-1AT )-1 b.
(g)
Replacing A in (f) by the expression (g) we get:
q = M-1AT (AM-1AT )-1 b.
(h)
Note that it is very difficult to calculate these results without help of machines; symbolic math program in MATLAB produces such cases within fraction of seconds.
Regarding that I1 = 1/2mr2, I2 = 1/4mr2 and supposing m =1 kg, r = 1 unit let us try to solve equation (h) using numerical methods. As the result,
M =
1 0 0 0 1 0
0 1 0 0 , AT = 0 1
0 0 1/2 0 — sin è cos è
0 0 0 1/4 0 0
T
b = \ipè cosè pè sinè] . (i)
As we can see from equation (i), equation (h) can't be solved analytically. But the position and velocity of the disk can be approximated using numerical integration.
To solve equation (i) using numerical methods let us reduce it into systems of first order ordinary differential equations:
x = x(1), y = x(3), p = x(5), § = x(7), X(1) = x(2), x(3) = x(4), x(5) = x(6), x(7) = x(8), x(8) = 0,
x(2) = + 2cos(x(7))^/( 1 + 2cos(x(7))2 + 2 sin(x(7))2) ■ x(6) x(8) cos(x(7)) +
+ 2sin(x(7))2 cos(x(7))/( 1 + 2cos(x(7))2 + 2sin(x(7))2) -x(6)x(8), x(4) = 2sin(x(7)) ■ cos(x(7))2/( 1 + 2cos(x(7))2 + 2 sin(x(7))2) -x(6)x(8)+
+ (1 + 2sin(x(7))2)/( 1 + 2cos(x(7))2 + 2 sin(x(7))2) ■ x(6) x(8) sin(x(7)),
x(6) = 2sin(x(7)) ■ (1 + 2cos(x(7))^/( 1 + 2cos(x(7))2 + 2 sin(x(7))2) + + 4 ■ cos(x(7))2 sin(x(7))/( 1 + 2cos(x(7))2 + 2sin(x(7))^ x x(6)x(8) cos(x(7)) +
+ ^—4 ■ sin(x(7))2 cos(x(7))/( 1 + 2cos(x(7))2 + 2 sin(x(7))2) + +2 cos(x(7)) ■ ( 1+2 sin(x(7))2) / ( 1+2 cos(x(7))2+2 sin(x(7))2) ^ xx(6) x(8) sin(x(7)).
Now let us apply stabilization of the constraint equations as discussed above (Fig. 3). In doing so, we have
x — rp sin § — rp§ cos § + g1 = 0, y + rp cos § — rp§ sin § + g2 = 0, or in matrix form Aq = b + Dip, where g = D^ and D = D(t) is a two by two matrix. From this the equation of motion for the disk will be
q = M-1AT (AM-1AT )-1(b + Dip).
position X
0 ^-1-1-1-1-1-1-1-1-1-
0 1 2 3 4 5 6 7 8 9 10 time t
Figure 3. Absolute difference of stabilization of constraint equations
for position x
Finally taking the following initial values for position and velocity we get the corresponding values of position and velocity at any time t
x0 = ®°(1) = 0, x0 = ®°(2) = 1, y0 = ®°(3) = 0, j/0 = ®°(4) = -V3,
= £0(5) = 0, p0 = ^0(6) = 2, tf0 = ^0(7) = ^/6, tf0 = ^0(8) = 1.
5. Conclusion
In this paper a modified method of constructing dynamic equation of constrained mechanical system is presented. The equation is done applying the principle of Lagrange with stabilization of constraint equations. It contains unknown coefficient matrix which determines the stability of the solution. The choice of the matrix is performed experimentally using MATLAB program during numerical solution.
The stability of the initial value problem we considered depends on the size and sign of elements of this matrix.
Finally an example is given and investigation is made to determine the possible values of the elements of the proposed matrix. During experimental investigation of the elements of the coefficients, the results show better approximation for small simulation time.
References
1. P. P. Teodorescu, Mechanical Systems; Classical Models, Springer Sci-ence+Business Media B.V., 2009.
2. R. A. Layton, Principles of Analytical System Dynamics, Springer, 1998.
3. A. Arabyan, F. Wu, An Improved Formulation for Constrained Mechanical Systems, Multibody System Dynamics 213 (2) (1998) 49-69.
4. R. G. Mukharlyamov, Equation of Motion of Mechanical Systems, RUDN, Moscow, 2001.
5. J. G. de Jalon, E. Bayo, Kinematic and Dynamic Simulation of Multibody Systems, Springer-verlag, 1988.
6. R. G. Mukharlyamov, Differential-Algebraic Equations of Programmed Motions of Lagrangian Dynamical System, Mechanics of Solids 46 (4) (2011) 534-543.
7. R. G. Mukharlyamov, On the Numerical Solution of Equations of Extremals of Variational Problems with Constraints, Mathematics 479 (4) (2002) 36-43.
8. R. G. Mukharlyamov, A. W. Beshaw, Solving Differential Equations of Motion for Constrained Mechanical Systems, Bulletin of PFUR. Series "Mathematics. Information Sciences. Physics" (3) (2013) 81-91, in Russian.
9. M. D. Ardema, Analytical Dynamics: Theory and Applications, Kluwer Academic / Plenum Publishers, 2005.
УДК 531.3
Уравнения динамики несвободной механической системы
А. В. Бешау
Кафедра математики Бахрдарский университет Эфиопия, Бахрдар
Работа посвящена модификации уравнений динамики механической системы со связями. Кинематическое положение системы определяется обобщёнными координатами и скоростями, на которые наложены связи. На основе принципа Даламбера—Лагранжа составляются уравнения движения в форме уравнений Лагранжа с неопределёнными множителями. Уравнения динамики приводятся к виду, разрешённому относительно ускорений. Выражения для неопределённых множителей определяются с учётом возможных отклонений от уравнений связей. Для стабилизации связей вводятся дополнительные переменные, используемые для оценки отклонений, вызванных погрешностями задания начальных условий и использования численных методов. Для аппроксимации решений обыкновенных дифференциальных уравнений, в частности, нелинейных уравнений первого порядка, используются явные численные методы. Построены линейные уравнения возмущений связей, матрица коэффициентов которых выбирается в процессе численного решения уравнений динамики. Устойчивость по отношению к начальным отклонениям от уравнений связей и стабилизация численного решения зависят от значений элементов этой матрицы. В результате исследования определяются допустимые значения матрицы коэффициентов, соответствующие решению уравнений динамики методом Эйлера и методом Рунге—Кутта четвёртого порядка. Предложенный метод решения задачи стабилизации используется для моделирования движения диска по плоскости без проскальзывания.
Ключевые слова: свободная система, голономные связи, неголономные связи, стабилизация, ряд Тейлора, численное решение.
Литература
1. Teodorescu P. P. Mechanical Systems; Classical Models. — Springer Sci-ence+Business Media B.V., 2009.
2. Layton R. A. Principles of Analytical System Dynamics. — Springer, 1998.
3. Arabyan A., Wu F. An Improved Formulation for Constrained Mechanical Systems // Multibody System Dynamics. — 1998. — Vol. 213, No 2. — Pp. 49-69.
4. Мухарлямов Р. Г. Уравнения движения механических систем. — Москва: РУДН, 2001. — 99 с.
5. de Jalon J. G., Bayo E. Kinematic and Dynamic Simulation of Multibody Systems. — Springer-verlag, 1988.
6. Mukharlyamov R. G. Differential-Algebraic Equations of Programmed Motions of Lagrangian Dynamical System // Mechanics of Solids. — 2011. — Vol. 46, No 4. — Pp. 534-543.
7. Мухарлямов Р. Г. О численном решении уравнений экстремалей вариационной задачи с ограничениями // Математика. — 2002. — Т. 479, № 4. — С. 36-43.
8. Mukharlyamov R. G., Beshaw A. W. Solving Differential Equations of Motion for Constrained Mechanical Systems // Bulletin of PFUR. Series "Mathematics. Information Sciences. Physics". — 2013. — No 3. — Pp. 81-91.
9. Ardema M. D. Analytical Dynamics: Theory and Applications. — Kluwer Academic / Plenum Publishers, 2005.
