Sage Journals: Discover world-class research

Abstract

In the process of solving differential algebraic equations of motion for constrained multi-body dynamic system, singular problems such as singular mass matrices, singular configurations, and redundant constraints may appear, which bring great difficulties to numerical simulation. The article presents a new equation of motion for multi-body dynamic system based on UK equations (explicit equations of motion proposed by Udwadia and Kalaba) and least square method, with which the dynamics problems are transformed into unconstrained optimization problems. Due to its pretty explicit and simple form, the new equation is very suitable for programming. Compared with the Lagrange method, the new equation does not require the introduction of Lagrange multipliers. Using this new equation, both the traditional optimization and modern optimization algorithms can be used in dynamics simulation for multi-body system. With the help of different optimization methods, approximate solutions can be found near the location of singular configurations to avoid instability of solution and numerical diffusion. In the end, several examples are given to illustrate that the new equation set is not only a practical modeling method for the above singular problems but also has faster solving speed and higher accuracy than the traditional methods.

Keywords

UK equations holonomic and nonholonomic singular problems redundant constraints multi-body system

Introduction

Multi-body system dynamics provides a powerful tool for the dynamic performance evaluation and optimization design for many engineering mechanical systems in the fields of weapon, aviation, spaceflight, vehicle, robot, and precision machinery. There is an urgent need for modern engineering on the research of dynamic modeling and the control design of complex multi-body systems. In recent years, the study on dynamic modeling, design, control solving strategy, and experiment of multi-body systems has been carried out in depth by many scholars at home and abroad and more and more papers based on this topic are published.^1–5

Modeling and solving are the core problems of multi-body systems dynamics. By now, different styles of multi-body system dynamics modeling methods have greatly promoted the development of modern technology and provide a variety of effective means for solving the problems of dynamics calculation of complex mechanical systems. Generally speaking, the existing modeling methods of the multi-body system dynamics are mainly based on the following basic principles and methods, such as New-Euler vector mechanics method, the mechanical analysis method derived from D’Alembert principle, the minimum extremum principle based on Gauss principle, and the transfer matrix method.^6–10

The dynamic equations gained by different modeling methods are not the same, and in most cases, they can be expressed as differential algebraic equations, nonlinear differential equations,⁵ and transfer equation of system boundary state.^6,11Although the dynamic equations of the multi-body system modeled by different methods differ from each other, they are theoretically equivalent. Scholars at home and abroad have made outstanding contributions to solving dynamic problems, and many numerical methods come into being. Q Wang and Q Lu,¹² TM Wasfy and AK Noor,¹ W Schiehlen,² and G Wang⁵ summarized the various computational methods of multi-body systems. The numerical solution of two-order differential equations of multi-body system dynamics based on Lagrange method is relatively mature. But the numerical solution of differential algebraic equation established based on the Newton–Euler method is still in the exploration.

In the process of solving differential algebraic equations of constrained systems, the following singular problems will be encountered such as constraint violation of speed and displacement, singular mass matrix, singular Jacobian matrix of constraints, and the redundant constraints which all may arise during numerical integration procedure. The singular mass matrix may appear when the complex mechanical system is described by redundant generalized coordinates. Generally speaking, when the system moves, it may approach some unusual locations, so-called singular configuration, which can result in acceleration mutation and simulation distortion. In these singular positions, the degree of freedom of the system has a sudden change, which leads to the singularity or redundancy of the constrained equations. The redundant problems can arise when one eagerly combines the sub-systems’ constraints together and ignores the facts whether the constraints are singular or not, although every sub-system’s constraints are independent.

In the process of finding an automatic modeling and solving method in multi-body dynamics, judging whether the constraint is singular or not brings a huge amount of work. Many theories assume that the constraint equations are independent from each other when dealing with the problems of redundant constraints.^13–15 However, this premise is often difficult to meet. First, it is hard to judge whether the constraint is independent or not during numerical integration procedure for a complex multi-body system. Then, the ultimate goal of multi-body dynamics is automation of modeling and simulation. It is almost inevitable to introduce redundant constraints in the process of numerical simulation. Finally, if the system has some singular configurations, the degree of freedom of the system changes with the motion of the system, resulting in the redundancy of the constraints. In a word, it is hard to handle the dynamic problems of the multi-body system on the premise of independent constraints.

It is challenging to develop a general method for solving differential algebraic equations, and more and more scholars are attracted to the work. One of the important methods to control the constraints is to modify the velocity constraint equations and the displacement constraint equations, so that the violation degree of the constraints is within the controllable precision range. Then, given that the number of the constraint equations is smaller than the variable of the system, we should find some ways to solve a set of underdetermined equations. The methods of solving the underdetermined equations mainly include generalized inverse method,^16,17 least-square method,¹⁸ constraint method, space correction method,^19,20 and so on. We can also make the constraint equation compatible with differential equations^21,22 by changing the integral form. If the constraint equations of the system are redundant, the independent constraint should be selected from the redundancy constraints.^19,20 There are many methods to identify independent constraints, such as the linear independence between the row vector of the Jacobi matrix of constraints or the constraint force.^15,23 Augmented Lagrange formulation claims to be able to solve dynamic problems of multi-body systems with redundant constraints and even achieves good results for the system with singular configurations.²⁴ Z Qi et al.²⁵ proposed a method of violation correction to guarantee the stability of solution.

For the system with singular configurations, many methods of solving differential algebraic equations, including augmented Lagrange method, are found to fail near the locations of singular configurations.²⁶ Because the generalized mass matrix of the system becomes singular in singular positions, the violation degree of the system is increased.

The matrix of the system is usually full rank in unconstrained dynamic systems, but often singular in the complex multi-body dynamic systems with constraints or redundant coordinates, as generalized coordinates exceed the minimum number required to describe the configuration of the system. In general, the higher the flexibility of the modeling method, the much more likely the mass matrix of the system to be singular.^27–30

UK equations were proposed by Udwadia and Kalaba to handle systems whether or not their mass matrices are singular, and they are general and explicit equations of the system motion. The UK equations are applicable to not only the holonomic systems but also the nonholonomic systems whether or not their mass matrices are symmetric and positive definite.²⁹ Furthermore, the generalized inverse matrix is introduced into UK equations, and the explicit expressions of the system acceleration are given without solving a set of overdeterminate equations. In this article, the application of UK equations in the multi-body mechanical systems with singular configurations or redundant constraints is discussed by several concrete examples.

Recently, some experts and scholars have studied the UK equations and have made some achievements, but there are some problems remained to be solved. Most of their concern is focused on the application of the UK equations, because the UK equations do not require the use of Lagrange multipliers or care about whether the matrices are singular or not. Some scholars applied the explicit analytic equations of motion to the flexible constraint systems, musical instruments, and dynamic control of redundant robot.^31–33 De Falco et al.³⁴ gave a report about UK equations on numerical efficiency in multi-body dynamics simulations and on teaching effectiveness.

But this article will discuss the applications of the UK equations to singular problems of the multi-body dynamic mechanical systems, such as singular matrix, redundant constraints, and singular constraints, by comparing with other methods from the aspects of the reliability, accuracy of results, and the speed of computation.

UK equations

In this section, the article mainly introduces derivation of UK equations proposed by Udwadia and Kalaba.^27–30

The motion of an unconstrained multi-body mechanical system at any time $t$ can be described, using Lagrange equations, by

M (q, \overset{\cdot}{q}, t) \overset{\cdot\cdot}{q} = Q (q, \overset{\cdot}{q}, t)

(1)

with the initial conditions

q (0) = q_{0}, \overset{\cdot}{q} (0) = {\overset{\cdot}{q}}_{0}

(2)

where $q$ is the generalized coordinate n-vector; $Q$ is n-vector which is the generalized force of the system; $M$ is the symmetric $n \times n$ matrix which is generalized mass matrix of the system.

Provided the system now subject to $m$ constraints

φ_{i} (q, \overset{\cdot}{q}, t) = 0, i = 1, 2, \dots, m

(3)

where $k \leq m$ , equations in equation (3) are functionally independent, and $k$ is the rank of the constraint set (3). Constraint set (3) describes the constraints including all holonomic and nonholonomic constraints.

Assuming that the constraint set is smooth enough, the differential relation of equation (3) with respect to time $t$ can be obtained

A (q, \overset{\cdot}{q}, t) \overset{\cdot\cdot}{q} = b (q, \overset{\cdot}{q}, t)

(4)

where $A \in R^{m \times n}$ and $b$ is a m-vector.

$Q^{c}$ will appear on the right side of equation (1) to ensure equation (4) is satisfied, when the system gets sudden constraints. The equations of motion of the system become the following expression

M \overset{\cdot\cdot}{q} = Q (q, \overset{\cdot}{q}, t) + Q^{c} (q, \overset{\cdot}{q}, t)

(5)

The virtual work by the generalized forces of constraints can be expressed as³⁵

w^{T} Q^{c} (q, \overset{\cdot}{q}, t) = w^{T} C (q, \overset{\cdot}{q}, t)

(6)

where $C (q, \overset{\cdot}{q}, t)$ is an n-vector constant determined by the mechanician or by experimentation and observation.³⁵ And when the system is ideal and the constraint forces do not work under virtual displacements, $C (q, \overset{\cdot}{q}, t)$ is a zero vector. The virtual displacement vector, $w (t)$ , is any non-zero n-vector which satisfies^36,37

A (q, \overset{\cdot}{q}, t) w = 0

(7)

Solving equation (7), $w (t)$ can be written as follows

w = (I - A^{+} A) r

(8)

where $r$ is an arbitrary n-vector and $A^{+}$ is the generalized inverse matrix $A$ .

Since each component of $r$ can be chosen to be independent, substituting equation (8) in equation (6), we can get

(I - A^{+} A) Q^{c} (q, \overset{\cdot}{q}, t) = (I - A^{+} A) C (q, \overset{\cdot}{q}, t)

(9)

By equations (5) and (9), we can get

(I - A^{+} A) M \overset{\cdot\cdot}{q} (q, \overset{\cdot}{q}, t) = (I - A^{+} A) (Q (q, \overset{\cdot}{q}, t) + C (q, \overset{\cdot}{q}, t))

(10)

By equations (4) and (10), we can get the following equations

[\begin{matrix} (I - A^{+} A) M \\ A \end{matrix}] \overset{\cdot\cdot}{q} = [\begin{matrix} (I - A^{+} A) (Q + C) \\ b \end{matrix}]

(11)

Defining

\bar{M} = [\begin{matrix} (I - A^{+} A) M \\ A \end{matrix}]

(12)

we can solve equation (11) and get

\overset{\cdot\cdot}{q} = {\bar{M}}^{+} [\begin{matrix} (I - A^{+} A) (Q + C) \\ b \end{matrix}] + (I - {\bar{M}}^{+} \bar{M}) η

(13)

where $η$ is an arbitrary n-vector.

When matrix ${\bar{M}}^{+}$ has full rank, $I - {\bar{M}}^{+} \bar{M} = 0$ and equation (11) can be changed into the form

\overset{\cdot\cdot}{q} = {\bar{M}}^{+} [\begin{matrix} (I - A^{+} A) (Q + C) \\ b \end{matrix}]

(14)

It is proved that equation (14) is applicable to these systems with holonomic constraints, nonholonomic constraints, or a combination of the two types.³⁵ However, UK equations also have some limitations. First, if $I - {\bar{M}}^{+} \bar{M} \neq 0$ , an arbitrary n-vector $η$ appears on the right side of equation (14), which will bring a lot of problems to our dynamic simulation. Second, as we all know that the computation of solving a generalized inverse is much bigger than solving a set of equations. The UK equations are not fit for complex systems because a generalized inverse item arises in the right side of equation (14). In the following section, the article just introduces an effective method to overcome the difficulties brought by the random item $η$ .

Numerical integration algorithm

Equation (14) which describes the acceleration of a complex dynamic system is general, simple, and understandable, but in the process of integrating with respect to time $t$ , there may be some problems, such as the constraints violation of velocity and displacement. In the process of numerical simulation, the constraint violation degree of velocity and displacement increases as operating time goes on.

In this section, the article discusses how to integrate equation (14) with respect to time $t$ and reduce the constraints violation degree of velocity and displacement.

After modeling the system with equation (14), it is necessary to choose an integral method to simulate the multi-body system. There are two types of numerical algorithms for solving dynamic equations: explicit methods and implicit methods. Many explicit algorithms are available, such as Runge–Kutta algorithm, central difference algorithm, and predictor-corrector algorithm. Although the explicit algorithm has a relatively high computational efficiency, its step size is generally relatively small in order to ensure the accuracy of the solution.

There are many differences between implicit and explicit algorithms. The implicit algorithms have larger step size and better stability compared with explicit algorithms under the condition of satisfying the accuracy of the system requirement. And they are suitable for solving stiff problems in differential equations which cannot be solved by explicit algorithms.

In this article, an explicit algorithm with two-order accuracy, the Newmark method, is introduced. And a new modeling method based on the UK equations and the Newmark method is proposed.

According to Newmark method proposed in 1959 by Professor Newmark,³⁸ we make the following assumptions in the calculation

{\overset{\cdot}{q}}_{t + Δ t} = {\overset{\cdot}{q}}_{t} + [(1 - δ) {\overset{\cdot\cdot}{q}}_{t} + δ {\overset{\cdot\cdot}{q}}_{t + Δ t}] Δ t

(15)

q_{t + Δ t} = q_{t} + {\overset{\cdot}{q}}_{t} Δ t + [(\frac{1}{2} - α) {\overset{\cdot\cdot}{q}}_{t} + α {\overset{\cdot\cdot}{q}}_{t + Δ t}] Δ t^{2}

(16)

where $Δ t$ is integral step size, and $α$ and $δ$ are integral constants. $α$ and $δ$ can be adjusted according to the stability and accuracy of the system requirement. It is proved by Newmark,³⁸ when $δ = 0.5$ and $α = 0.25$ , the Newmark method turns into constant average acceleration method which is unconditionally stable.

By equations (15) and (16), the following expression can be obtained

{\overset{\cdot\cdot}{q}}_{t + Δ t} = \frac{1}{α Δ t^{2}} (q_{t + Δ t} - q_{t}) - \frac{1}{α Δ t} {\overset{\cdot}{q}}_{t} - (\frac{1}{2 α} - 1) {\overset{\cdot\cdot}{q}}_{t}

(17)

and

{\overset{\cdot}{q}}_{t + Δ t} = \frac{δ}{α Δ t} (q_{t + Δ t} - q_{t}) + (1 - \frac{δ}{α}) {\overset{\cdot}{q}}_{t} + (1 - \frac{δ}{2 α}) {\overset{\cdot\cdot}{q}}_{t}

(18)

Substituting equations (17) and (18) in equation (13), we can obtain

q_{t + Δ t} = α Δ t^{2} ({\bar{M}}^{+} [\begin{array}{l} (I - A^{+} A) (Q + C) \\ b \end{array}] + (I - {\bar{M}}^{+} \bar{M}) η) + Δ t {\overset{\cdot}{q}}_{t} + (\frac{Δ t^{2}}{2} - α Δ t^{2}) {\overset{\cdot\cdot}{q}}_{t} + q_{t}

(19)

Provided the matrix $\bar{M}$ has full rank, then

I - {\bar{M}}^{+} \bar{M} = 0

(20)

And equation (19) can be transformed into the following form

[\begin{array}{l} (I - A^{+} A) M \\ A \end{array}] (\frac{1}{α Δ t^{2}} (q_{t + Δ t} - q_{t}) - \frac{1}{α Δ t} {\overset{\cdot}{q}}_{t} - (\frac{1}{2 α} - 1) {\overset{\cdot\cdot}{q}}_{t}) = [\begin{array}{l} (I - A^{+} A) (Q + C) \\ b \end{array}]

(21)

Since equation (21) is overdetermined, we choose least square method and unconstrained optimization method to get the approximate solution.

Defining

f = [\begin{array}{l} (I - A^{+} A) M \\ A \end{array}] (\frac{1}{α Δ t^{2}} (q_{t + Δ t} - q_{t}) - \frac{1}{α Δ t} {\overset{\cdot}{q}}_{t} - (\frac{1}{2 α} - 1) {\overset{\cdot\cdot}{q}}_{t}) - [\begin{array}{l} (I - A^{+} A) (Q + C) \\ b \end{array}]

(22)

then, the 2-norm of equation (23) can be written as

g = Minimize {‖ f ‖_{2}^{2}}

(23)

But if equation (20) is not satisfied, the solution of equation (23) may diverge. In order to ensure the convergence and stability of solution of the optimization, here we add a regularization item $‖ q ‖_{2}^{2}$ on the right side of equation (23). And the new equation has the following form

g = Minimize {‖ f ‖_{2}^{2} + μ ‖ q ‖_{2}^{2}}

(24)

where $μ \geq 0$ which is a small constant.

Equation (24) is called UK optimization equation in this article which can be applied to the multi-body systems with singular configuration. Besides, the error of the constraint equations is given as follows

error = ‖ A \overset{\cdot\cdot}{q} - b ‖_{2}

(25)

The constraint problems in complex, multi-body mechanical system are transformed into optimization problems by UK optimization equation which holds the advantages of both UK equations and Newmark method. And also the UK optimization equation has some significant advantages compared with classical modeling methods as follows:

The generalized coordinates used to describe the configuration of the system are no longer independent any more. Instead, they can exceed the minimum number required by the system, which greatly simplify the form of the equations of motion of the complex system.

Compared with Lagrange’s method, the UK optimization equation is explicit and general for constrained discrete dynamic systems. Besides, it avoids to calculate the Lagrange multipliers which are always too difficult to obtain for those complex systems with multiple degrees of freedom.

Whether the constraints are ideal or not, the UK optimization method can be used in such systems with holonomic, nonholonomic constraints, or a combination of the two types.

As the UK optimization equation is based on the Newmark method and UK equations, it is unconditionally convergent and also has a second-order accuracy. Because the inertial matrix of the system is semi-positive and equation (24) is a convex two-order programming, the global minimum point of equation (24) must exist.

The traditional optimization methods, the modern optimization algorithms, or even both kinds of methods can be used to carry out the numerical simulation for complex multi-body dynamics system modeled by the UK optimization equations, such as Lagrange multiplier method, penalty function method, genetic algorithm, and ant colony algorithm.

Baumgarte method³⁹ is one of the famous violation correction methods which has the following form

\overset{\cdot\cdot}{φ} + 2 χ \overset{\cdot}{φ} + γ^{2} φ = 0

(26)

where $χ$ and $γ$ are regular parameters that have no law to follow. Compared with traditional violation correction method like Baumgarte method, the UK optimization method does not need to determine the values of irregular parameters. Furthermore, the UK optimization method proposed by the article has the form of differential equation rather than a differential algebraic equation which can simplify the calculation process. So it can be seen that the UK optimization has a higher computational efficiency than the traditional violation correction method.

The flowchart of the program is shown in Figure 1 where $t$ is the simulation time and $Δ t$ is the simulation time step.

Figure 1.

Flowchart of the program.

In this article, we use the traditional unconstrained optimization algorithm as the main method to stimulate the dynamic behavior of the multi-body system.

Numerical simulation

In this section, the article provides four different examples to discuss how UK optimization equation handles the singular problems of the multi-body dynamics. The four examples are, respectively, a skate on an inclined plane with nonholonomic constraints, a four-bar linkage with a singular mass matrix, a crank slide with singular configurations, and a rolling disk in a vertical direction under redundant constraints.

The UK optimization equations, the Routh equations, and the Lagrange method are applied to the four examples. And comparison results are given in the form of figures.

Example 1

Consider a skate on an inclined plane with Cartesian coordinates $x, y$ , as shown in Figure 2. We assume the y-axis to be horizontal, and the x-axis directed downward. $(x, y)$ is the coordinate of the contact point of the balanced skate with the plane, and $φ$ is the system’s rotation angle measured from y-axis. The equation of the nonholonomic constraint is

$\overset{\cdot}{x} \sin φ - \overset{\cdot}{y} \cos φ = 0$
(27)

Figure 2.
Skate on an inclined plane.

The integral step $Δ t = 0.005$ . Choosing the appropriate units of mass, length, and time, we can represent the Lagrangian function in the form

$L = ({\overset{\cdot}{x}}^{2} + {\overset{\cdot}{y}}^{2} + {\overset{\cdot}{φ}}^{2}) / 2 + x$
(28)

When $t = 0$ , the initial conditions of the system compatible with the nonholonomic constraint are $x (0) = y (0) = φ (0) = \overset{\cdot}{x} (0) = \overset{\cdot}{y} (0) = 0$ and $\overset{\cdot}{φ} = 1$ .

The Routh equation is a product of the combination of ordinary indefinite multipliers and Lagrange’s equation of the second kind. It can be used in nonholonomic mechanical systems.^40,41

In this example, the UK optimization equations and Routh equations are, respectively, used to simulate the dynamic problem taking $x$ , $y$ , and $φ$ as the generalized coordinates. Specially, Routh equations use Newmark algorithm as the integration method in the numerical simulation.

Discussion 1

In example 1, the article uses genetic algorithm^42,43 as the optimization method to get the global minimum point of equation (24). Figures 3 and 4 show the curves of generalized coordinates $x$ , $y$ , and $φ$ versus time $t$ . Figures 5 and 6 show the violation degree of displacement constraint equation during the numerical computations. Here, the error of the constraint equation is defined as equation (25), which should theoretically be zero. It can be seen from the Figures 5 and 6 that the constraint equation has an error of 10^–5 within the error tolerance simulated by both methods. But UK optimization equations take less CPU time than Routh equations and the simulation time spent by the two methods is, respectively, 34 and 58 s. This example proves that the UK optimization equations not only can be applied to holonomic systems but also nonholonomic systems.

Figure 3.
Time histories of the displacement of $x$ , $y$ , and $φ$ obtained by UK optimization equations.

Figure 4.
Time histories of the displacement of $x$ , $y$ , and $φ$ obtained by Routh equations.

Figure 5.
The error curve of the constraint equation obtained by UK optimization equations.

Figure 6.
The error of the constraint equation obtained by Routh method.

Example 2

The second example deals with a four-bar linkage with three rotation angles $θ_{1}$ , $θ_{2}$ , and $θ_{3}$ , as shown in Figure 7. The three rods, which have the same quality $m$ and length $l$ , do periodic motion under the general force $ψ$ in the $xy$ plane. We assume the y-axis to be downward, and the x-axis directed horizontal. The system has 1 degree of freedom and we take $θ_{1}$ , $θ_{2}$ , and $θ_{3}$ as the system’s generalized coordinates along with the holonomic constraints given by

$Φ = [\begin{matrix} θ_{1} - θ_{3} = 0 \\ θ_{1} + θ_{2} = π \end{matrix}]$
(29)

Figure 7.
Four-bar linkage.

The system’s kinetic energy is

$T = \frac{1}{2} m {\overset{\cdot}{θ}}_{1}^{2} + \frac{1}{2} m {\overset{\cdot}{θ}}_{2}^{2} + \frac{1}{2} m ({\overset{\cdot}{θ}}_{1} l)^{2}$
(30)

And the system’s potential energy can be simply expressed as

$V = \frac{1}{2} mgl \sin θ_{1} + \frac{1}{2} mgl \sin θ_{2} + mgl \sin θ_{1}$
(31)

We can use Lagrangian function $L = T - V$ to get the system’s mass matrix which is singular

$M = [\begin{matrix} \frac{4}{3} m l^{2} & 0 & 0 \\ 0 & \frac{4}{3} m l^{2} & 0 \\ 0 & 0 & 0 \end{matrix}]$
(32)

The parameters describing the system and its constraints are $m = 12$ , $l = 2$ , $g = 10$ , $Δ t = 0.005$ , and $ψ = 100$ . The initial conditions are chosen to be $θ_{1} (0) = θ_{3} (0) = π / 3$ , ${\overset{\cdot}{θ}}_{2} (0) = (2 / 3) π$ , and ${\overset{\cdot}{θ}}_{1} (0) = {\overset{\cdot}{θ}}_{2} (0) = {\overset{\cdot}{θ}}_{3} (0) = 0$ so that the constraints are compatible with each other.

Discussion 2

The simulation results modeled by UK optimization equations can be described as Figures 8 and 9, which respectively show the time histories of the displacement and the error in satisfying the constraint. The above results can be seen to be equivalent to Figures 10 and 11 simulated by Lagrange’s method of the first kind (to within the error tolerance). Both of the modeling methods have the same local tolerance level of $10^{- 12}$ . The time that it takes for the simulation is, respectively, 12 and 25 s. This section takes a four-bar linkage as an example to show that the UK optimization equations not only can be used in singular mass matrix problems but also has a much more efficient simulation speed than Lagrange multiplier method.

Figure 8.
Time histories of the displacement of $θ_{1}$ , $θ_{2}$ , and $θ_{3}$ obtained by UK optimization equations.

Figure 9.
The error in constraint equations vs time t obtained by UK optimization equations.

Figure 10.
Time histories of the displacement of $θ_{1}$ , $θ_{2}$ , and $θ_{3}$ obtained by Lagrange’s equations of the first kind.

Figure 11.
The error curve of constraint equation obtained by Lagrange’s equations of the first kind.

Example 3

Consider a crank slider included two rods with the same quality and length, moving in the $xy$ plane, free of any given forces, as shown in Figure 12, with the gravitational acceleration $g$ downward. We assume that the x-axis is horizontal, and the y-axis is vertically downward. The slider can only slide along the x-axis and its Cartesian coordinates is $(x, 0)$ . The two rods have the same quality and length, $m$ and $l$ . $θ_{1}$ and $θ_{2}$ are both rotation angle measured from the x-axis. The system has 1 degree of freedom and we take $θ_{1}$ , $θ_{2}$ , and $x$ as the system’s generalized coordinates along with two holonomical constraints given by

$Φ = [\begin{matrix} l \cos θ_{1} + l \cos θ_{2} - x \\ l \sin θ_{1} - l \sin θ_{2} \end{matrix}]$
(33)

Figure 12.
The crank slider.

Differentiating equation (30) twice, we obtain

$A = [\begin{matrix} - l \sin θ_{1} - l \sin θ_{2} - 1 \\ \cos θ_{1} - \cos θ_{2} 0 \end{matrix}]$
(34)

When $θ_{1} = θ_{2} = (π / 2) + k π (k \in Z)$ , the constraint equation is singular and the rank of its coefficient matrix $A = 1$ . The parameters used to describe the configuration of the system are $m = 6$ , $l = 1$ , $g = 10$ , and $Δ t = 0.05$ . The initial conditions of the system compatible with the holonomic constraints are

$θ_{1} (0) = θ_{2} (0) = \overset{\cdot}{x} (0) = {\overset{\cdot}{θ}}_{1} (0) = {\overset{\cdot}{θ}}_{2} (0) = 0, x (0) = 2$

Discussion 3

Figures 13 –16, respectively, show the curves of generalized coordinates and the curves of constraints error versus time $t$ modeled by UK optimization equations and the first kind Lagrange’s method from the 2630th s to 2660th s. Figures 15 and 16 modeled by Lagrange’s method of the first kind show that the simulation results diverge with the reason that the coefficient matrix $A$ of the constraint equations becomes singular when the system approaches some singular configurations. And the violation of the constraint equations gradually becomes larger and larger in the process of numerical simulation which badly exceeds the error tolerance. But Figures 13 and 14 modeled by UK optimization equation smoothly pass through the singular configurations and the violation degree does not exceed permitted range. The simulation time modeled by UK optimization equations and Lagrange’s method is, respectively, 589 and 916 s, which proves the efficiency of the UK optimization equations method.

Figure 13.
Time histories of the displacement of $θ_{1}$ , $θ_{2}$ , and $x$ obtained by UK optimization equations.

Figure 14.
The error curve of constraints vs time t modeled by UK optimization equations.

Figure 15.
Time histories of the displacement of $θ_{1}$ , $θ_{2}$ , and $x$ obtained by Lagrange’s equations of the first kind.

Figure 16.
The error in constraint equations vs time t modeled by Lagrange's equations of the first kind.

Compared with the traditional modeling methods, the UK optimization equations can get an approximate solution at the location of singular figurations with optimization methods and avoid instability of solution and numerical diffusion. A conclusion comes into being that the UK optimization equations can be seen as a practical simulation model in solving singular configuration problems of multi-body system from the view of accuracy and efficiency of calculation.

Example 4

Consider a rolling disk of mass $m$ and radius $R$ in plumb plane, free of any given forces, as shown in Figure 17. We assume the y-axis to be horizontal, and the x-axis directed downward. Let $(x, y)$ be the Cartesian coordinates of the wheel center. $ψ$ and $φ$ are, respectively, the angle of the wheel itself and the angle between the wheel axis and the x-axis. The system has 2 degrees of freedom and we take $ψ$ , $φ$ , $x$ , and $y$ as the system’s generalized coordinates along with three nonholonomical constraints given by

$Φ_{1} = [\begin{matrix} \overset{\cdot}{x} + R \overset{\cdot}{φ} \cos ψ \\ \overset{\cdot}{y} + R \overset{\cdot}{φ} \sin ψ \\ \overset{\cdot}{x} \sin ψ - \overset{\cdot}{y} \cos ψ \end{matrix}]$
(35)

$Φ_{2} = [\begin{matrix} \overset{\cdot}{x} + R \overset{\cdot}{φ} \cos ψ \\ \overset{\cdot}{y} + R \overset{\cdot}{φ} \sin ψ \end{matrix}]$
(36)

Figure 17.
The rolling disk in vertical direction.

The parameters describing the system and its constraints are $m = 1$ , $R = 2$ , $g = 10$ , and $Δ t = 0.005$ . The initial conditions of the system which are compatible with the nonholonomic constraint are

$x (0) = y (0) = ψ (0) = φ (0) = \overset{\cdot}{y} (0) = 0, \overset{\cdot}{x} (0) = - 2, \overset{\cdot}{ψ} = 1, φ = 1$

The system’s constraints $Φ_{1}$ are redundant, and the rank of the nonholonomical constraints is two. Some modeling methods need to filter independent constraints from the redundant constraints, but UK optimization equations can simulate dynamics of nonholonomic systems with redundant constraints directly.

Considering the Routh method cannot solve the dynamics problems with redundant constraints, in this example, a comparison of simulation results is carried out between UK optimization equation with redundant constraints and Routh method with independent constraints.

Discussion 4

Figures 18 and 19, respectively, describe the curves of the system generalized coordinates and the error of the constraint equation versus time $t$ simulated by the UK optimization equations. Due to the instability of the Routh equation in the redundancy problems, Figures 20 and 21 show the simulation results with independent constraints $Φ_{2}$ . It can be easily seen that both modeling methods have the same local tolerance level of $10^{- 12}$ and also the displacement curves do not differ greatly seen from Figure 22. The time spent for simulation is, respectively, 155 and 153 s. The UK optimization equations do not have to filter independent constraints from the redundant ones and they nearly have a same solving speed as the Routh method, so there comes the conclusion that the UK optimization equations are much more suitable to the redundant constraint problems compared with Routh method.

Figure 18.
Time histories of the displacement of $ψ$ , $φ$ , $x$ , and $y$ with redundant constraints obtained by UK optimization equations.

Figure 19.
The error curve of redundant constraints obtained by UK optimization equations.

Figure 20.
Time histories of the displacement of $ψ$ , $φ$ , $x$ , and $y$ obtained by Routh equations with independent constraints.

Figure 21.
The error curve of the independent constraints obtained by Routh equations.

Figure 22.
The difference curve of the displacement $ψ$ obtained by two different equations.

Conclusion

In order to solve the differential algebraic equations of constrained systems, researchers have to be faced with singular problems such as singular mass matrix, singular configurations, and redundant constraints, which bring great difficulties to the numerical simulation work. In this article, a new UK optimization equation based on UK equations and the Newmark method is proposed. After comparing the different simulation results obtained by UK optimization equation and traditional methods, the following conclusions are drawn:

The constraint problems in complex, multi-body mechanical systems are transformed into unconstrained optimization problems by UK optimization equations which hold the advantages of both UK equations and Newmark equations. This change just opens up new ways for solving the dynamic problems of multi-body systems. Mature traditional optimization methods and modern optimization algorithm or even the combinations of both kinds of methods can be applied to the dynamic solving problems of complex multi-body systems. The UK optimization equations can get an approximate solution in the singular figuration with optimization methods to avoid the solution divergence.

The UK optimization equations can be applicable to holonomic, nonholonomic systems, or a combination of both types. And the equations do not require, nor use, the notion of Lagrange multipliers which simplify the calculation.

The UK optimization equations are explicit and general to handle the complex multi-body systems whether or not their mass matrices are singular, whether or not their constraint equations are redundant or singular, and whether or not their constraints are ideal. The example of a crank slider demonstrates that the system simulated by UK optimization equations can smoothly approach the singular configurations, and the violation degree does not exceed the permitted range without considering constraint violation correction. Finally, the article takes a rolling disk as an example to verify that the UK optimization equation can be a practical modeling method to simulate the multi-body system with redundant constraints.

In terms of computational efficiency, from the article’s examples, we can see that the UK optimization equations have the advantages of faster calculating speed and higher accuracy than traditional methods such as Lagrange’s method and Routh method when dealing with singular problems in multi-body dynamic system.

The UK optimization equation used to describe the motion of constrained discrete dynamic systems in this article is pretty explicit and general without constraint violation correction. Considering its higher computational efficiency and simpler expression than the classical modeling methods, the UK optimization equation is much suitable for programming.

Footnotes

Handling Editor: Chuanzeng Zhang

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This paper was supported by the National Natural Science Foundation of China (11272167).

References

Wasfy

Noor

AK.

Computational strategies for flexible multibody system. Appl Mech Rev 2003; 56: 553–613.

Schiehlen

Computational dynamics: theory and applications of multibody systems. Eur J Mech A: Solid 2006; 25: 566–594.

Yoo

Kim

et al . Developments of multi-body system dynamics: computer simulations and experiments. Multibody Syst Dyn 2007; 18: 35–58.

Liu

Hong

Research and prospect of flexible multi-body systems dynamics. Chin J Comput Mech 2008; 25: 411–416 (in Chinese).

Wang

Numerical algorithms of multibody system dynamics. Comput Simul 2006; 23: 86–89 (in Chinese).

Rui

Yuan

et al . Transfer matrix method of multibody system and its application. Beijing, China: Science Press 2008 (in Chinese).

Dynamics of flexible multibody system. Beijing, China: Higher Education Press, 1996 (in Chinese).

Schiehlen

Nils

Robert

Multibody dynamics in computational mechanics and engineering applications. Comput Method Appl M 2006; 195: 5509–5522.

Schiehlen

Research trends in multibody system dynamics. Multibody Syst Dyn 2007; 18: 3–13.

10.

Schiehlen

Multibody system dynamics: roots and perspectives. Multibody Syst Dyn 1997; 1: 149–188.

11.

Rui

et al . Discrete time transfer matrix method for multibody system dynamics. Multibody Syst Dyn 2005; 14: 317–344.

12.

Wang

Advances in the numerical methods for Lagrange’s equations of multibody systems. Adv Mech 2001; 31: 9–17 (in Chinese).

13.

Shabana

AA.

Dynamics of multi-body systems. 3rd ed.Cambridge: Cambridge University Press, 2005.

14.

Hong

Computational dynamics of multi-body systems. Beijing, China: Higher Education Press, 1999 (in Chinese).

15.

Dynamics of multi-body systems. Beijing, China: Science Press, 2008 (in Chinese).

16.

Hong

A new violation correction method for constrained multibody systems. Acta Mech Sinica 1998; 30: 300–306 (in Chinese).

17.

Zhao

Pan

Chen

Constraint violation stabilization of Euler-Lagrange equations with non-holonomic constraints. Acta Mech Solida Sin 2004; 17: 45–51.

18.

Zhao

Pan

Least square algorithms and constraint stabilization for Euler-Lagrange equations of multibody system dynamics. Acta Mech Sinica 2002; 34: 594–603 (in Chinese).

19.

Zhang

Approximation of displacement and model reforming in multi-body systems. J Dalian Univ Technol 2001; 4: 149–151 (in Chinese).

20.

A modified numerical method for dynamical equations of the multi-body systems with singularity position. J Dyn Control 2006; 4: 109–113 (in Chinese).

21.

Pogorelov

Differential-algebraic equations in multibody system modeling. Numer Algorithms 1998; 19: 183–194.

22.

Arnold

Fuchs

Fiihrer

Efficient corrector iteration for DAE time integration in multi-body dynamics. Comput Method Appl M 2006; 195: 6958–6973.

23.

Wang

Liu

Dong

Studies on dynamic solution of redundant multi-body system. J Tianjin Univ 1998; 31: 11–16 (in Chinese).

24.

Jalon De

Bayo

Kinematic and dynamic simulation of multibody systems: the real time challenge. New York: Springer-Verlag, 1994.

25.

Fang

Study on redundant constraints in multibody systems. Acta Mech Sinica 2011; 43: 390–399 (in Chinese).

26.

Neto

Ambrosio

Stabilization methods for the integration of DAE in the presence of redundant constraints. Multibody Syst Dyn 2003; 10: 81–105.

27.

Udwadia

Kalaba

RE.

Equations of motion for mechanical systems. J Aerospace Eng 1996; 9: 64–69.

28.

Udwadia

Kalaba

RE.

What is the general form of the explicit equations of motion for constrained mechanical system. J Appl Mech 2002; 69: 335–339.

29.

Udwadia

Phohomsiri

Explicit equations of motion for constrained mechanical systems with singular mass matrices and applications to multi-body dynamics. P Roy Soc A: Math Phy 2006; 462: 2097–2117.

30.

Udwadia

Kalaba

Phohomsiri

Mechanical systems with nonideal constraints: explicit equations without the use of generalized inverses. J Appl Mech 2004; 71: 615–621.

31.

Pennestrì

Valentini

Falco

DD.

An application of the Udwadia–Kalaba dynamic formulation to flexible multibody systems. J Frankl Inst 2010; 347: 173–194.

32.

Han

Zhen

Chen

YH.

Dynamic modeling and simulation of multi-body systems using the Udwadia-Kalaba theory. Chin J Mech Eng 2013; 26: 839–850.

33.

Antunes

Debut

Dynamical computation of constrained flexible systems using a modal Udwadia-Kalaba formulation: Application to musical instruments. J Acoust Soc Am 2017; 141: 764–778.

34.

De Falco

Pennestrì

Vita

. The Udwadia-Kalaba Formulation: A Report on its Numerical Efficiency in Multibody Dynamics Simulations and on its Teaching Effectiveness. In: Multibody Dynamics 2005 - ECCOMAS Thematic Conference, Madrid, Spain, 21–24 June 2005.

35.

Udwadia

Kalaba

RE.

Explicit equations of motion for systems with non-ideal constraints. J Appl Mech 2001; 68: 462–467.

36.

Udwadia

Kalaba

RE.

Analytical dynamics: a new approach. Cambridge: Cambridge University Press, 1996.

37.

Udwadia

Kalaba

Hee-Chang

Equations of motion for mechanical systems and the extended D’Alembert’s principle. Q Appl Math 1997; 55: 321–331.

38.

Newmark

NM.

A method of computation for structural dynamics. J Eng Mech Div: ASCE 1959; 85: 67–94.

39.

Baumgarte

A new method of stabilization for holomic constraints. J Appl Mech 1983; 50: 869–870.

40.

Thowsen

The Routh–Hurwitz method for stability determination of linear differential–difference systems. Int J Control 1981; 33: 991–995.

41.

Mei

Liu

Luo

Advanced analytical mechanics. Beijing, China: Science Press, 1991 (in Chinese).

42.

Holland

JH.

Adaptation in natural and artificial systems. Cambridge, MA: MIT Press, 1992.

43.

Goldberg

DE.

Genetic algorithms in search, optimization and machine. Reading, MA: Addison-Wesley, 1989.

UK method for solving multi-body system dynamics with singular problems

Abstract

Keywords

Introduction

UK equations

Numerical integration algorithm

Numerical simulation

Example 1

Discussion 1

Example 2

Discussion 2

Example 3

Discussion 3

Example 4

Discussion 4

Conclusion

Footnotes

Declaration of conflicting interests

Funding

References