STABILIZATION METHODS FOR SIMULATIONS OF C O N S T R A I N E D M U L T I B O D Y D Y N A M I C S By Hong Sheng Chin B. Sc.(Mathematics) Liaoning University, Shenyang, China, 1982. M. Sc.(Mathematics) Liaoning University, Shenyang, China, 1984. A T H E S I S S U B M I T T E D IN PARTIAL F U L F I L L M E N T O F THE REQUIREMENTS FOR THE DEGREE OF D O C T O R OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF MATHEMATICS We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA May 1995 © Hong Sheng Chin, 1995 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature) Department of f^cCtke^KMACS The University of British Columbia Vancouver, Canada Date DE-6 (2/88) Ap(vL <*/ , /3<?r Abstract The descriptor form of constrained multibody systems, and any general formulation of such systems with closed loops, yield index-3 differential-algebraic equations (DAEs). Generally, some index reduction techniques have to be used before numerical discretization can be safely applied. However, a direct differentiation of the constraints introduces mild instability. Hence one must consider index reduction with stabilization. One popular method for index reduction with stabilization is Baumgarte's technique. However the difficulty of choosing the parameters in practice makes this method's robustness unclear. Moreover our numerical experiments show that there are still large constraint drifts even with Baumgarte's stabilization. In the thesis, we employ concepts and techniques of dynamical systems in order to improve the situation. We first study the relationship between a DAE and its underlying vector fields. A general form of vector fields with stabilized invariant manifolds is given. We propose a new numerical stabilization method for semi-explicit index-2 and index-3 DAEs of Hessenberg form. Our stabilization method improves on Baumgarte's stabilization which is widely used in engineering as well as in simulations of multibody systems. We then develop a numerical code based on the new stabilization method with an adaptive step-size control for the descriptor form of the Euler-Lagrange equations in multibody systems. Numerical simulations have been conducted with our code on a variety of multibody systems including a spatial five-link-suspension model in a vehicle, Andrews squeezing mechanism, a two-arm manipulator with a prescribed motion and a mechanism with kinematic singularity. Our code is efficient, fast and therefore is more attractive for real-time simulations. 11 When high speed and light-weight substructures are involved in the multibody system, the rigid body model is usually no longer valid. In such a case we compute the elastic deformations and oscillations of the substructures using the finite element method. Satisfactory numerical results using our method are presented for deformable multibody systems as well. 111 Table of C o n t e n t s Abstract ii List of Tables vii List of F i g u r e s viii Aknowledgement 1 2 Introduction 1 1.1 Simulation of multibody dynamics 1 1.2 Objectives and contributions of the Thesis 4 I n d e x R e d u c t i o n s of D A E s and Reformulations 9 2.1 2.2 3 xi Differential Algebraic Equations and Index 16 2.1.1 Solvability and Index of DAEs 17 2.1.2 Index Reduction for Higher Index DAEs 19 Vector Field on a Manifold and its Stabilizations 23 2.2.1 DAEs and their Underlying Vector Fields 23 2.2.2 Stability of a DAE and its reformulations 26 2.2.3 Stabilizations of Underlying ODEs 28 2.3 Index reduction and stabilization for index-2 DAEs 30 2.4 Stabilization for Index-3 DAEs 39 S t a b i l i z a t i o n s of D i s c r e t e D y n a m i c s Formulations IV 44 4 5 3.1 Dynamics of Difference Equations 47 3.2 Stability of Invariant Manifolds in Difference Equations 50 3.3 Post-Stabilizations for ODEs with Invariant Manifolds 58 3.4 Stabilization for Multibody Dynamics with Holonomic Constraints . . . . 66 3.5 Stabilization for Nonholonomic Systems 71 3.6 Stabilizations and Projection Methods 73 S i m u l a t i o n of C o n s t r a i n e d Rigid M u l t i b o d y D y n a m i c s 74 4.1 Problem formulation 77 4.1.1 Absolute coordinates 77 4.1.2 Relative coordinates 79 4.2 Implementation with Runge-Kutta Methods 81 4.3 A Numerical Code: RKSTAB 88 4.4 Numerical simulations of holonomic constrained systems 91 4.4.1 Andrews squeezing mechanism 93 4.4.2 A vehicle suspension model 97 4.4.3 A planar two-link manipulator with prescribed path 106 4.4.4 Slider-Crank mechanism with a TSDA 112 4.5 RKSTAB on a singular slider-crank mechanism 116 4.6 A Nonholonomic constrained system 123 P o s t - S t a b i l i z a t i o n in a D e f o r m a b l e M u l t i b o d y S y s t e m 129 5.1 Introduction 129 5.2 Generalized coordinates in deformable multibody systems 131 5.3 Equations of motion for the deformable slider-crank mechanism 135 5.4 Simulations of the deformable slider-crank mechanism 142 v 6 Conclusions 150 A p p e n d i x : A Brief Review of Mechanics of Deformable Bodies 155 A A Brief Review of Mechanics of Deformable Bodies 155 VI List of Tables 1 Comparison of Baumgarte and two improved methods 35 2 Baumgarte and Baumgarte improved for an index-3 DAE 42 1 Maximum drifts of a two-arm robot 70 1 Comparison of B(20,100) with PStab on squeezing mechanism 94 2 Results of RKSTAB on Andrews squeezing mechanism 96 3 Comparison of MEXX and RKSTAB on squeezing mechanism 97 4 Results of RKSTAB on the suspension mechanism, case I 5 Comparison of B(40,400) and PStab on suspension model, case II 6 Comparison of B(12,70) with PStab on the manipulator in case I Ill 7 Comparison of Baumgarte's parameters on the manipulator in case II . . Ill 8 Comparison of RKSTAB and DASSL on crank-slider with a TSDA 116 1 Slider-crank parameters 102 . . . . ... 106 143 vii List of Figures 2.1 Drifts on position level of a 7-body system with rtol=l.d-4 12 2.2 Drifts on velocity level of a 7-body system with rtol=l.d-4 12 2.3 Drifts on position level of a 7-body system with rtol=l.d-7 13 2.4 Drifts on velocity level of a 7-body system with rtol=l.d-7 14 2.5 Drifts in position for different Baumgarte parameters 14 2.6 Drifts in velocity for different Baumgarte parameters 15 4.1 Andrews squeezing mechanism 93 4.2 Solution of Andrews squeezing mechanism 94 4.3 Drifts with B(20,100) on squeezing mechanism 95 4.4 Drifts with PStab on squeezing mechanism 95 4.5 A five-point vehicle suspension mechanism 98 4.6 Suspension graph with cut-joints indicated 99 4.7 Vertical position of the vehicle wheel, case I 101 4.8 Solution (<f>i, V>i,a, /3,7) of the suspension, case I 4.9 Drifts with B(10,30) for the suspension model, case I 101 102 4.10 Drifts with PStab for the suspension model, case I 103 4.11 Vertical position of the vehicle wheel, case II 104 4.12 Solution (c^x, ^ 1 , a,/?,7) of the suspension, case II 104 4.13 Drifts with B(40,400) for the suspension, case II 105 4.14 Drifts with PStab for the suspension, case II 105 4.15 A two-link planar robotic system 107 viii 4.16 Constraint path of the end-effector, Case I 108 4.17 Joint constraint torques of the manipulator, Case I 109 4.18 y-coordinate solution of the end-effector, Case I 109 4.19 Joint constraint torques of the manipulator for u) = 0.5 110 4.20 y-coordinate solution of the end-effector in Case II 110 4.21 Slider-crank mechanism with a TSDA 112 4.22 Solutions of slider-crank with a TSDA, £=1000 114 4.23 Velocities of slider-crank with a TSDA, £=1000 114 4.24 Solutions of slider-crank with a TSDA, £=10000 115 4.25 Velocities of slider-crank with a TSDA, £=10000 115 4.26 A singular slider-crank mechanism 117 4.27 Solutions of the singular system with rtol = l.d-5 118 4.28 Solutions of the singular system with rtol = l.d-8 118 4.29 Step sizes for the singular system without stabilization 119 4.30 Step sizes for the singular system with PStab 119 4.31 Solution of the singular system with B(20,100), rtol=l.d-10 120 4.32 Solution of the singular system with PStab, rtol=l.d-10 121 4.33 Solution of new formulation for singular slider-crank 122 4.34 A rolling disc on a horizontal plane 124 4.35 The orbit of the rolling coin center 124 4.36 Drifts of rolling disk with post stabilization 126 4.37 Drifts of rolling disk with Baumgarte's stabilization 126 4.38 The orbit of contact points in the plane 127 4.39 The height of the rolling disk center 127 5.1 132 A planar deformable body ix 5.2 A two-dimensional beam element 134 5.3 A slider-crank mechanism with deformable connecting rod 135 5.4 Connection condition between two beam elements 137 5.5 Transverse deflection at the midpoint for <^'(0) = 180 143 5.6 Transverse deflection at the midpoint for <^'(0) = 120 144 5.7 Transverse deflection at the midpoint for </>'(0) = 60 5.8 Drifts of deformable model with B(30,100) for c^'(0) = 180 145 5.9 Drifts of deformable model with PStab for ^'(0) = 180 146 144 5.10 Drifts of deformable model with B(16,90) for </>'(0) = 60 147 5.11 Drifts of deformable model with PStab for 0'(O) = 60 147 5.12 Drifts in the rigid model with B( 10,80) for <//(0) = 60 148 5.13 Drifts in the rigid model with PStab for ^'(0) = 60 148 A.l 156 Deformed and Undeformed State of a Flexible Body A.2 Notation for stress 158 x Aknowledgement The author would like to express his heartfelt appreciation and gratitude to his thesis advisor Dr. Uri Ascher for his continual encouragement and patient guidance throughout the course of this work. Many thanks are also extended to the thesis committee members: Drs. Wayne Nagata, Dinesh Pai, Jim Varah and Brian Wetton. Their helpful comments during the research and writing of this thesis are greatly appreciated. At last, but not least, the author would like to take this opportunity to thank his wife Dawn and his daughter Alina for all their help and understanding. Without the family's encouragement, love and support, this work would never have been accomplished. XI Chapter 1 Introduction 1.1 S i m u l a t i o n of m u l t i b o d y d y n a m i c s The dynamics of multibody systems, such as motion of robotic manipulators, mechanical chains and spacecrafts, is becoming increasingly important in engineering. A computer simulation of such multibody systems requires a concerted integration involving several computational aspects [47, 81, 82]. These include selection of a data structure for the system's configuration, computerized generation of governing equations of motion, incorporation of constraint conditions and implementation of suitable solution algorithms. Basic methods for multibody system simulations are provided by the disciplines of dynamics (the multibody formulations), numerical mathematics and computer science [26]. Algorithms for the simulation of complex mechanical systems require the understanding of the application in question. The configuration of a multibody system is identified by a set of variables called generalized coordinates that completely defines the location and orientation of each individual body and an approximation of deformation fields of deformable substructures in the system. Depending on the choice of generalized coordinates, the equations of motion for a multibody system can be obtained mainly by two kinds of formulations: the augmented and the recursive formulations. In an augmented formulation, the configuration of a (rigid) system is identified by using a set of Cartesian coordinates to indicate the locations and orientations of each individual component in the multibody system. This approach has the advantage that the 1 Chapter 1. Introduction 2 formulation of the equations is straightforward. It allows easy additions of complex force functions and constraint equations [19, 20]. Typically six or seven coordinates are used for each spatial rigid body in a system. The connectivity conditions between adjacent components in the system are introduced to the dynamic formulation by applying a set of nonlinear algebraic constraint equations. This set of constraint equations are adjoined to the system's equations of motion by using Lagrange multipliers [4, 46, 82, 81]. The dynamic equations are then index-3 DAEs for holonomic constrained systems and index2 DAEs for nonholonomic constrained systems. In such a formulation the mass matrix and the Jacobian matrix of the constraints have sparse structures. In a recursive formulation, relative joint variables are used to formulate a minimum set of differential equations [53, 4]. For multibody systems with a tree structure so-called O(N) formalisms have been found by various researchers independently, yielding the state space representation of the system dynamics in explicit form with a number of operations growing linearly with the number N of system bodies [53, 47, 26, 96]. The dynamics equations in a recursive formulation are written in terms of the system's degrees of freedom and are of lower dimensionality than those in an augmented formulation. However, the corresponding matrices are normally dense and the generalized force functions to be evaluated are often more cumbersome. The situation becomes more complicated when there is a deletion or addition of force or a change in the kinematic constraint structures in a multibody system [19, 20]. For a closed-loop multibody system, however, there are closure conditions which remain as algebraic equations. In this case a recursive formulation has to be combined with a cut-joint technique in order to derive the equations of motion. In general the equations of motion for a closed-loop system are DAEs with index-3 or index-2 as stated earlier. When high operating speeds and greater position accuracy are required in a simulation of a multibody system which contains lightweight elements, the rigidity assumption for Chapter 1. Introduction 3 those lightweight mechanical elements is no longer valid [19, 85]. The multibody system under investigation is then assumed to consist of a set of interconnected rigid and flexible substructures. The system configuration is identified by using a coupled set of reference and elastic coordinates [85]. Elastic coordinates can be the generalized coordinates in an assumed modes method or the displacements and the slopes of selected nodal points in a finite element method [95, 85]. The system differential and algebraic equations of motion which are developed by using the Lagrangian formulation are written in terms of a coupled set of reference and elastic modal coordinates [85]. Obviously the efficient integration of DAEs can play an important role in simulations of multibody dynamics. Such systems arise also in many other applications. There are various methods developed for the efficient integration of DAEs and the equations of motion for constrained multibody systems [3, 8, 12, 15, 25, 47, 65, 96]. A number of commercial codes have recently been developed for the generation and integration of the equations of motion of mechanical systems [46, 80, 47, 63, 65, 44, 5]. These codes perform well in many circumstances, but they are also known to be too slow in certain applications, especially in simulations involving person in-the-loop. This thesis concentrates on the efficient integration of the equations of motion for constrained multibody systems. To that end we investigate index reduction with stabilization for semi-explicit DAEs and stabilization for difference equations with invariant manifolds. We develop the theory and methods for numerically solving DAEs by viewing them as ODEs with invariants. Our study leads to post-stabilization methods for numerical solutions of vector fields with invariant manifolds. One of our goals is to eventually approach real-time simulations of multibody systems. Chapter 1. 1.2 Introduction 4 O b j e c t i v e s and contributions of t h e T h e s i s DAEs can be regarded as ODEs on certain invariant manifolds by index reduction. Thus the numerical solutions of the DAEs can be obtained by integration of their underlying ODEs. In certain circumstances difficulties may occur since the numerical solutions of the underlying ODE could drift away from the invariant manifold. Therefore index reduction with stabilization for higher index DAEs becomes important for numerical solutions of such DAEs. In the thesis, we employ concepts and techniques of dynamical systems to characterize the relationship between a semi-explicit DAE and its underlying vector fields. There are infinitely many underlying ODEs for a given DAE and the stability properties of the invariants can be greatly different in different underlying vector fields. Numerically even a mild unstability of the invariants could leads to an unacceptable large numerical drift when an underlying ODE is integrated [12, 8]. By constructing certain Lyapunov functions we propose a general form of underlying vector fields with stabilized invariant manifolds for a given semi-explicit DAE. Numerically those underlying ODEs provide better reformulations for the index reduction of the original DAE. For semiexplicit index-2 DAEs our general form of stabilization includes Baumgarte's method when a certain special form of stabilization matrix is chosen. For semi-explicit index3 DAEs, however, our stabilization is different from that of Baumgarte as the invariant manifold with the new stabilization has a monotonic property which Baumgarte's method does not share. It is important to notice that the stability property of an invariant manifold in an underlying vector field can perish after a discretization. With Baumgarte's stabilization the invariant manifold can be unstable in the difference equations even though the manifold is asymptotically stable in the underlying vector field. The stability of the manifold Chapter 1. Introduction 5 in the difference equations directly relates to the numerical drift of a numerical solution. This is partially the reason that large numerical drifts can result with Baumgarte's stabilization in certain cases. Inspired by the stabilization formulation we proposed for underlying vector fields, we consider stabilization of invariant manifolds for difference equations obtained from discretization of underlying vector fields. Since the stabilization is performed after the discretization we call it post-stabilization. The post-stabilization is efficient and it corresponds to one iteration of certain coordinate projection techniques [3, 11, 5, 8, 24]. Post-stabilization was mentioned also in [25]. Our consideration here is to reformulate the difference equations so that the invariant manifold is uniformly stable. Therefore the solutions of the post-stabilization algorithms have very small numerical drifts and retain high accuracy. Our novel contribution is in the mathematically rigorous development of this family of post-stabilization methods. In our consideration several post-stabilization algorithms for vector field with invariants and equations of motion in multibody systems are proposed. As an application of the post-stabilizations we derive a stabilization method for simulation of constrained multibody dynamics. Based on the latter post-stabilization method we have developed a numerical code with an automatic step-size control for the EulerLagrange equations for constrained multibody systems. The code, RKSTAB, is an adaptation of the ODE code DOPRI5(4) of Hairer, Norsett and Wanner [43]. The numerical scheme used in the code is an embedding Runge-Kutta 5(4) method, which originates from the formula due to Dormand and Prince [23] where a minor modification in the 4th order solution has been made. Certain comparisons of our numerical code with other codes are made for various multibody systems. Numerical results show that our code is reliable, efficient and fast, therefore it is more attractive for real-time simulations. We also apply the post-stabilization for deformable multibody systems. The equations of motion for a deformable slider-crank mechanism are set up by using the finite Chapter 1. Introduction 6 element method which yields a system of index-3 DAEs. Using the post-stabilization we proposed for simulation of constrained multibody dynamics, satisfactory numerical results are obtained for the deformable multibody systems. In general our formulations of index reduction with stabilization can be applied to any higher index ( > 3) semi-explicit DAEs. The family of stabilizations we proposed in the thesis includes some existing stabilization methods for index-2 DAEs when a special form of stabilization matrix is chosen. Improvements on existing stabilization methods are demonstrated by numerical examples. For difference equations resulting from a discretization of an underlying vector field we stabilize the invariant manifold in a similar way as we consider for a vector filed. With the post-stabilization methods the invariant manifold becomes uniformly stable in the stabilized difference equations. Therefore the drifts in the stabilized numerical solutions are small. With the analysis of the stability of invariant manifolds in difference equations we derive a general form of post-stabilization algorithms for numerical solutions of vector fields with invariants. For numerical solutions of Euler-Lagrange equations in multibody systems we develop poststabilization methods for holonomic and nonholonomic constrained systems. Numerical experiments based on our post-stabilization algorithms on a variety of multibody systems show that our algorithms are efficient and reliable. O u t l i n e of t h e t h e s i s In Chapter 2 we briefly review index reduction and stabilized index reduction methods for semi-explicit DAEs. We consider stabilizations for an ODE with an invariant manifold, because a DAE can be regarded as a vector field on a certain invariant manifold. A new kind of stabilized vector fields in which the invariant manifold is asymptotically stable is given by using Lyapunov's direct method. For semi-explicit index-2 DAEs this general form of stabilization includes Baumgarte's method as a special choice of the stabilization matrix. Improvements on Baumgarte's stabilization are demonstrated by Chapter 1. Introduction 7 numerical examples. Our stabilization methods can be used for any semi-explicit high index DAEs. We present index reduction with stabilization for semi-explicit index-2 DAEs using the general formulation. The relationship between our stabilizations and projected coordinates [12, 35] is also considered. We then turn to the differences in stabilization formulations for semi-explicit index-2 and index-3 DAEs. For semi-explicit index-3 DAEs the invariant manifold has a monotonic property in our stabilized vector fields and this is one of the differences between ours and Baumgarte's method. In Chapter 3 we consider stabilization methods in discrete dynamical systems. We present an example which shows that the invariant manifold becomes unstable after a discretization when Baumgarte's stabilization has been used. Therefore stabilization of the manifold in difference equations is needed as the drifts in numerical solutions have a direct link with the stability property of the manifold in difference equations. We propose stabilization methods in discrete dynamical systems which can be difference equations from a discretization of a vector field. We prove that with the post-stabilization the manifold is uniformly stable and therefore the numerical drifts are very small in the stabilized difference equations. truncation error 0(hp+1), stabilization algorithms. Furthermore for an actual discretization with a local we prove that numerical drifts are 0(h2(p+1)) in our post- For the Euler-Lagrange equations in multibody systems we derive our method of choice "PStab" which embodies a compromise between the economy of stabilization step and its effectivity for the drifts. In Chapter 4 we consider the issue of simulations of multibody dynamics. We first consider problem formulations in a multibody system by using absolute and relative coordinates. Then we come to discretization of the Euler-Lagrange equations. Combined with PStab we have to modify the 4th order solution in DOPRI(5,4) [23, 43] in order that the scheme remains essentially a 6-stage method. We prove that the modified solution has order 4 and therefore the local error estimation remains valid as that in Chapter 1. Introduction 8 DOPRI(5,4). For a general formulation of equations of motion in a multibody system the mass matrix can be singular [5, 65]. We therefore rewrite PStab into a new form. The new formulation of PStab also sheds some light on making a full use of certain linear algebra algorithms designed for large multibody systems. Therefore PStab can be expected to be implemented efficiently for those large systems as well. Based on PStab and the modified DOPRI(5,4) we develop a numerical code, RKSTAB, with an adaptive step-size control. With a small additional effort we include Baumgarte's stabilization in RKSTAB as an option left to the user. This also allows us to make a systematic comparison of PStab with Baumgarte's stabilization. Using the code we conduct a variety of multibody dynamics simulations. Numerical experiments show that our code is robust and fast. In Chapter 5 we first consider the issues in simulation and dynamic analysis of deformable multibody systems. Simulations of deformable systems can be more complicated as there are reference and elastic coordinate couplings in the dynamics. Furthermore, high positioning accuracy is required and the tolerance of numerical drifts can be very small in a deformable system. We therefore apply PStab in the simulation of the deformable multibody system. For a deformable slider-crank mechanism we set up the equations of motion using the finite element method considered in [85]. The motion of the mechanism is mainly caused by large initial angular velocity of the crank. We consider the different elastic deformations of the connecting rod for different initial velocities. We also compare a rigid model and deformable model for a relative small initial velocity. In the deformable model RKSTAB takes much smaller step-sizes than those in the rigid model and this is partially due to the fact that the equations of motion are essentially PDEs. W i t h Baumgarte's stabilization large numerical drifts are observed both in the rigid and deformable model. With PStab the numerical drifts are very small for the both models. Chapter 2 I n d e x R e d u c t i o n s of D A E s and Reformulations In many physical systems which are described by initial value or boundary value differential problems, certain physical quantities, such as the net charge or the total energy are conserved as the system evolves. Other facts like kinematic relations between adjacent substructures in multibody dynamics or the non-negativity of concentrations in chemical reactions may be implicitly implied on physical or mathematical grounds. Usually those physical systems are described by a mixed set of differential equations and algebraic equations. Since the last decade much interest has been shown in the numerical solutions of differential-algebraic equations(DAEs)[12, 17, 44, 76]. In general consider an implicit ordinary differential equation(ODE): F(t,x(*),x'(<)) = 0 where F ( t , x , z) is continuous on some domain in R2n+1. (2.1) If the Jacobian _F z (£,x,z) = g^-(i,x, z) is nonsingular at some point ( i , x , z) then by the implicit function theorem, (2.1) is equivalent to an explicit ODE near that point. If F z ( i , x , z) is singular for all (i, x, z) in some region, then (2.1) is a DAE. One special case for DAE (2.1) is the situation where the derivatives of certain components of x in (2.1) do not appear explicitly. DAEs are different from (explicit) ODEs. In general, the existence of the solution for such an initial value problem is not known. Not all initial value DAE problems admit solutions even for linear constant coefficient DAEs. In order to have a solution, the initial 9 Chapter 2. Index Reductions of DAEs and Reformulations 10 values must be consistent i.e. they must satisfy certain algebraic equations. Numerically a DAE can be an ill-posed problem and special care must be taken when some discretization scheme is applied to it. The difficulty level for a DAE can be indicated by its index. The higher the index the more difficult the DAE appears to be. For higher index DAEs, some index reduction technique typically has to be used before the discretization. But a well-observed phenomenon in solving such a reduced index DAE is the drift problem [12, 15, 5, 29], i.e. the numerical solutions of such reduced index DAE can end up far away from the manifold defined by the original constraints. This could be a serious problem as the numerical solution may quantatively lose its physical meaning or other important properties. To see the drift problem let us consider the simulation of constrained multibody dynamics. The descriptor form of constrained mechanical systems with holonomic constraints yields an index-3 DAE: q' = v (2.2a) M(q)v' = f(q,v)-GT(q)A (2.2b) 0 = g(q) (2.2c) where q G R n « is the generalized coordinate vector, A G R71-1' is the Lagrange multiplier vector, G(q) = gq(q) and GGT has a (bounded) inverse. In DAE (2.2), (2.2c) represents n\ constraints on the position coordinates which indicate the link or joint conditions between the adjacent bodies. (2.2c) may also represent a prescribed motion in robotic motion planning. In many methods the holonomic constraints are differentiated with respect to time t once or twice before the numerical discretization. With two differentiations of the position constraints, for example, we have 0 = G(q)v' + Jj(G(q)v)q' = G(q)v' + z(q,v) (2.3) Chapter 2. Index Reductions of DAEs and Reformulations 11 where z(q,v) = ^ ( G ( q ) v ) q ' = ^ g ( q ) ( v , v ) Using the equation (2.3) we can eliminate the Lagrange multiplier vector in (2.2b) and obtain an ODE for (q, v) which is called the state space form : where q' = v (2.4a) v' = M-1(q)f(q,v)-M-1(q)GT(q)A(q,v) (2.4b) 1 A = A(q,v) = (G(q)M-1(q)G'T)-1(G(q)M-1(q)f(q,v) + z(q,v)) (2.5) It is important to notice that there are differences between the original descriptor form (2.2) and the state space form (2.4). They are not equivalent to each other. The constraints on position level and those on velocity level which are obtained by one differentiation of the position constraints with respect to time t, together define an invariant manifold for the state space form (2.4). Only those solutions of (2.4) which lie on the invariant manifold are the solutions of the descriptor form. Also for the descriptor form (2.2), there are infinitely many different underlying ODEs. All of the underlying ODEs admit the solutions of the original DAE (2.2) but they can present different stability properties for the invariant manifold. In numerical integrations, even a mild instability of the invariant manifold may cause serious drift problems i.e. the numerical solutions of (2.4) violate the constraints (2.2c) and (2.3) badly. Figs. 2.1 and 2.2 show the drifts on the position and velocity levels for a multibody system (Andrew's squeezing mechanism, 1 In this Chapter and Chapter 3 we assume that the mass matrix M(q) is symmetric positive definite. A general consideration of M(q) is made in Chapter 4. Chapter 2. Index Reductions of DAEs and Reformulations drifts in position constraints 0.06 0.04 0.02 -0.02 -0.04 - -0.06. Figure 2.1: Drifts on position level of a 7-body system with rtol=l.d-4 drifts in velocity constraints -0.2 -0.4 Figure 2.2: Drifts on velocity level of a 7-body system with rtol=l.d-4 12 Chapter 2. Index Reductions x 41 of DAEs and -|Q-4 Reformulations 13 drifts in position constraints i i i 1 1 1 r 3- 2- 1 - 0- -1 - -2- -3.41 1 1 1 1 1 1 1 1 1 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 2.3: Drifts on position level of a 7-body system with rtol=l.d-7 considered in Chapter 4) when the state space form (2.4) is integrated. The integration is carried out by the code DOPRI 54 [43], which has an automatic step size control. For any choice of (the local error) tolerances the drift eventually will build up to an unacceptably large amount as t grows. Figs. 2.3 and 2.4 are the plots of drifts with a much smaller tolerance. The increasing tendency of the drifts remains. In 1972 Baumgarte [15] introduced a stabilization method aiming at reducing the drifts. In fact Baumgarte's stabilization technique is to reformulate the equations of descriptor form to an underlying ODE such that the invariant manifold is asymptotically stable. However, the difficulty of choosing Baumgarte's parameters in practice makes its robustness unclear [8, 7]. Figs. 2.5 and 2.6 are the drifts in position and velocity constraints in Andrew's squeezing mechanism for various choices of Baumgarte parameters. We chose the time interval to be [0,0.3] and the absolute error = l.d-5 and the relative error = l.d-4. The (usual) z-direction represents the drifts on the velocity Chapter 2. Index Reductions x of DAEs and Reformulations 1Q-3 14 drifts in velocity constraints 1.5 1 0.5 0 •0.5 -1 ' 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 2.4: Drifts on velocity level of a 7-body system with rtol=l.d-7 Figure 2.5: Drifts in position for different Baumgarte parameters Chapter 2. Index Reductions of DAEs and Reformulations 15 constraints. The ^-direction indicates the varying of Baumgarte's parameters and the y-direction represents the percentage of time t in the interval [0,0.3]. In this example, the different choices of Baumgarte parameters do not present significant differences and all the Baumgarte parameters (large or small) appear to be inefficient. The drifts will build up eventually for different Baumgarte parameters. In order to integrate the DAEs more efficiently we have to study other more efficient stabilization methods. To this end we begin our study with the asymptotic stability of invariant manifolds in a vector field. We derive different reformulations or underlying ODEs for a given DAE. Among the reformulations proposed in the index-2 case, Baumgarte stabilization is a special choice. In our stabilization reformulations the invariant manifolds have a decreasing property. With Baumgarte's stabilization the invariant manifolds do not have this decreasing property for index-3 DAEs. In this chapter various reformulations of index-2 and index-3 DAEs of Hessenberg form Chapter 2. Index Reductions of DAEs and Reformulations 16 [17] are studied. In order to study those reformulations we briefly present some concepts of DAEs and index reduction methods for DAEs in Section 1. In Section 2 we study the relationship between a DAE and its underlying ODEs. For a given DAE we present a general form of underlying ODEs in which the invariant manifold is asymptotically stable. In Section 3 we study index-reduction reformulations for semi-explicit index- 2 DAEs. Numerical examples show that the reformulations we proposed improve on Baumgarte stabilization. In Section 4 we consider the reformulations and stabilizations for index-3 DAEs, especially the differences in the reformulations for index-2 DAEs and for index-3 DAEs. The methods presented in this chapter are general and can be also used for even higher index(> 4) DAE reformulations. 2.1 Differential A l g e b r a i c E q u a t i o n s and I n d e x DAEs are not ODEs [72] although people often relate a given DAE to an ODE in some sense. For instance, a popular definition of the differential index of a DAE is the minimum number of constraint differentiations by which the DAE can be converted to an ODE for all original solution components. Analogous to those of ODEs, the existence and uniqueness of solution, the asymptotic behavior of the solutions and Lyapunov stability for DAEs have been considered [76, 67]. Those studies enrich the theory of DAEs and lead to a better understanding of DAEs. However, the classical numerical methods for stiff ODEs such as BDF and implicit Runge-Kutta methods can not be used efficiently for all DAEs [11, 17, 44]. This fact keeps us focusing on the key differences between DAEs and ODEs in order to make efficient integrations of DAEs. Chapter 2. Index Reductions 2.1.1 of DAEs and Reformulations 17 Solvability and I n d e x of D A E s Consider a nonlinear DAE: F(*,x,x') = 0 (2.6) where F ( i , x , z) is continuous in a domain Q C R 2 n + 1 . A solution x(t) of (2.6) on an interval I is denned as a continuously differentiable function on I that satisfies (2.6) for all t £ I. It is well known that every ODE with local Lipschitz continuity has a unique solution when an initial value is given, at least locally. In fact, the continuity of the right hand side function of an ODE is sufficient for existence of a solution for an initial value problem. For DAEs, however, it is not simply the matter of smoothness, but also the DAE's algebraic structure that determines the existence of solutions. To see this consider an initial value problem of a linear DAE with constant coefficients: Ax!(t) + Bx(t) where x = (xi,X2)T. + f (t) = 0, x(0) = 0 (2.7) If / = (0,0) T then (2.7) has infinitely many solutions given by x\{t) = — /J e T _ i ^ ( r ) J r , x%{t) = £(t) , where £(t) can be any continuously differentiable function satisfying {(0) = 0. When / = (0,1) T , (2.7) does not have any solution at all (just look at the second row in (2.7)). Therefore both the nonhomogenous term f and the structure of matrices A and B are crucial for the existence of solutions. In fact, for a constant linear DAE, it is the so-called regularity of the matrix pencil that determines the existence of the solution [39]. This leads to the concept of solvability[17] for general DAEs (2.6) . Definition 2.1 (solvable) Let I be an open subinterval of R, ft a connected open subset of R2n+1, and F a differentiable function from tt to Rn. Then DAE (2.6) is solvable on Chapter 2. Index Reductions of DAEs and I if there is an r-dimensional Reformulations 18 family of solutions (j)(t,c) defined on a connected open set I x fi, VI C Rr, such that: 1. <f)(t, c) is defined on I for each c £ fl 2. (t, <f>(t, c), <f>'(t, c)) £ 0 for (f,c)elx0 3. If ij>(t) is any other solution with (2, ?/>(£), ?/>'(£)) £ tt, then ip(t) = 4>(t,c) for some c £ ft 4- The graph of <f> as a function of(t,c) is an (r+1)-dimensional manifold. From the definition above DAE (2.7) is not solvable. For a fully implicit DAE (2.6), the criteria for solvability are not clear and many considerations [76] are made based on the constant-rank assumption: rank Z ) z F ( i , x , z) = constant < n For any function $ : R('+ 1 ) n + 1 —> R n , $ = $ ( t , ^ 0 , ^ , • • • ,£,), the pseudoderivative D $ of $ is a new function : R('+ 2 ) ri + 1 _>. R n defined by [17] as D$(t, &,&,•••, 6+i) = -^- + ^ - 6 + • • • + ^ 6 + i where the partials are evaluated at (£,£o,£i,- • • , 6 ) - The powers of Z) are defined as Dm§ = D{Dm-l§). The differential index of DAE (2.6) is then defined as follows[17]. Definition 2.2 (differential i n d e x ) The index of (2.6) is the minimal nonnegative in- teger m such that F has m continuous derivatives and the nonlinear system of "derivative equations" JD m F(t,Xo,Xi) = 0 DF(*,Xo,Xi,x 2 ) = 0 F(i,x0,Xi,x2,---,xm+1) = 0 Chapter 2. Index Reductions of DAEs and Reformulations 19 can be solved on an interval I for Xi in terms oft and Xo .' Xi = <^>(t,Xo). The differential index of a DAE defined above can be related to the sensitivity of the solutions to perturbations in the equations. In particular Hairer et.al. [42] introduce the perturbation index concept which can directly relate the perturbation index to a perturbation estimate. Definition 2.3 ( p e r t u r b a t i o n index) Equation (2.6) has perturbation index m along a solution x on [0,t], if m is the smallest integer such that, for all functions x having a defect F ( i , x , x ' ) = *(t) there exists an estimate on [0,t] ||x(t) - x(*)|| < C(||x(0) - x(0)|| + max | | ^ ) | | + • • • + max where the expression on the right-hand side is sufficiently ^ ( O H ) small. Here C denotes a con- stant which depends only on F and the length of the interval. The differential index and the perturbation index are equal when the DAE has the form: Bx.' = F1(t,x) (2.8) where B is a constant matrix. In general the following relations hold [34]: differential index < perturbation index < differential index + 1 2.1.2 I n d e x R e d u c t i o n for Higher I n d e x D A E s Higher index DAEs are ill-posed problems and index reduction should generally be made before the integration. Let us first consider several DAE examples which we will study frequently in the thesis. They are semi-explicit DAEs in Hessenberg form [42]. Chapter 2. Index Reductions of DAEs and Reformulations 20 The index-1 case: x' = f(x,y) (2.9a) 0 = g(x,y) (2.9b) where g y = <9g/<9y has a bounded inverse in a neighborhood of the solution (or in the region of interest). The initial value (x 0 ,yo), to be consistent, satisfies g(x 0 ,yo) = 0. The index-2 case: x' = f(x,y) (2.10a) 0 = g(x) (2.10b) where g x f y has a bounded inverse. A consistent initial value (x 0 ,yo) satisfies g(x 0 ) = 0 and g x ( x 0 ) f ( x 0 , y o ) = 0. The index-3 case: x' = f(x,y) (2.11a) y' = h(x,y,z) (2.11b) 0 = g(x) (2.11c) where g x fyhz has a bounded inverse. An example of (2.11) is the equations (2.2) describing the dynamics of multibody systems. To be consistent, the initial value (xo,yo,Zo) is required to satisfy g(x 0 ) = 0, g x ( x 0 ) f ( x 0 , y o ) = 0 and g x x ( f , f ) + gxfxf + gxfyh = 0 at (x0,yo,z0). Notice that all three DAEs (2.9), (2.10) and (2.11) can be written in the form BX' = F(t,X) where B is a diagonal matrix with entries 1 and 0. So for DAEs (2.9), (2.10) and (2.11) the differential index and the perturbation index are equal to each other. Chapter 2. Index Reductions of DAEs and Reformulations 21 Numerical integration of higher index DAEs usually involves two steps [12, 17]: t h e stabilization or regularization by reformulating the DAE to a lower index DAE and the discretization of the reduced index DAE. There are several methods of reformulations for higher index DAEs in the literature. For illustration we once again consider the Euler-Lagrange equations for a constrained multibody system as an example. M(p)p" = 0 = f(p,v,t)-GT(p)\ (2.12a) g(p) (2.12b) I n d e x r e d u c t i o n by differentiation If we replace (2.12) by the equation M(P)P" 0 = f(p,v,*)-GT(p)A (2.13a) = - f (2.13b) Gp(v,v) + G ( / - G T A ) then the resulting DAE (2.13) has index 1. The equation (2.13b) has a double zero eigenvalue, so g(x(t)) will have a polynomial growth as t —> oo. This is partly responsible for the growing drifts when we integrate the equations (2.13). B a u m g a r t e stabilization To reduce the drift Baumgarte[15] proposed to use f(p,v,i)-GT(p)A (2.14a) ° = £?+*! + ™ < 2 - i4b > Af(p)p" = to replace (2.12). Baumgarte's technique has been used for many practical computations because it is simple to implement and understand, and can make use of ODE solvers. Of course at each time step some linear solver has to be used in the evaluation of the corresponding right-hand side function for the state space form. If there is no other Chapter 2. Index Reductions of DAEs and Reformulations 22 stiffness the integration can be efficient as an explicit scheme can be used. One of the problems with Baumgarte stabilization lies with the choice of parameters 7;. It is not obvious in practice how to choose suitable parameters [8, 12]. Analytically for any positive 71 and 72 the ODE (2.14b) is asymptotically stable which implies g(x(£)) —• 0 as t —> 00. Mathematically, the larger 7,-, the faster g(x(t)) damps out as t —• 00. But in numerical integration large 7; can introduce stiffness into the equation. There is no general recipe for choosing the right parameters and there is no guarantee that we can always find the Baumgarte parameters so that the drift will be controlled to the desired small amount. For certain problems small Baumgarte parameters are preferred while for the others large parameters have to be used [8]. P r o j e c t e d invariants m e t h o d s Another technique is maintaining more constraints by introducing an additional multiplier^) JJL [12, 36, 37]. By using this technique DAE (2.12) can be reformulated as [37] p' = v + G T (p) / u (2.15a) Mv' = f(p,v,t)-GT(p)A (2.15b) 0 = Gv (2.15c) 0 = g(p) (2.15d) The system (2.15) is an index-2 DAE for variables (p, v, A,/u). The exact solution for JJ, of (2.15) is fj, = 0 so that (2.15) and (2.12) will share the same solutions for (p, v, A). As the numerical solution of (2.15) satisfies both the position constraint (2.12b) and velocity constraint (2.15c) the method has no drift problem. But the computation of (2.15) could be expensive as implicit schemes have to be used for the even larger dimension DAE (2.15). Chapter 2. Index Reductions 2.2 of DAEs and Reformulations 23 V e c t o r Field on a Manifold and its Stabilizations In this section we study the relationship between a DAE and its underlying ODEs. A general form of underlying ODEs for a given DAE is given. Then we consider the stability of the invariant manifold in the underlying ODEs by using the Lyapunov's direct method. This can shed light on stabilization reformulations for a given DAE. It is well known that DAEs can be regarded as ODEs on certain manifolds, so we consider ODEs with invariant manifolds. For such ODEs we give a general form of stabilized ODEs or reformulations in which the invariant manifold considered is asymptotically stable. 2.2.1 D A E s and their U n d e r l y i n g Vector Fields A DAE can be regarded (at least locally) as a vector field on a certain manifold [77]. Mathematically, this manifold is an invariant set of the vector field in both directions of time t. By invariant set we mean that an orbit of the vector field will remain in the manifold all the time if the orbit starts from it. Simple examples of invariant sets are periodic orbits in the state space and the level curves in Hamiltonian systems [41, 45]. Due to the algebraic (constraint) nature, DAEs are always associated with certain invariant manifolds. For higher index DAEs there are more invariant conditions than just defined by the visible constraints, so the invariant manifolds are more complicated. For illustration consider the constrained multibody system with holonomic constraints: q' = v (2.16a) M(q)v' = f(q,v)-GT(q)A (2.16b) 0 = g(q) (2.16c) As we noted at the beginning of this chapter, the Lagrange multiplier A can be solved by differentiating the constraints (2.16c) twice, so the state space form ODE for (q, v) Chapter 2. Index Reductions of DAEs and 24 Reformulations has the form: q' = v (2.17a) v' = M-1(q)f(q,v)-M-1(q)G,T(q)A(q,v) (2.17b) where A(q,v) = (G(q)M-1(q)G'T)-1(G(q)M-1(q)f(q,v) + Z ( q , v ) ) (2.18) In order that a solution of the state space form ODE (2.17) be a solution of DAE (2.16) the initial value for (2.17) must be consistent. For DAE (2.16), consistent initial values are those which satisfy the position and velocity constraints. We therefore define a manifold which is invariant in the state space form vector field (2.17): M = {(q, v) : g(q) = 0, G(q)v = 0} (2.19) DAE (2.16) is equivalent to the ODE (2.17) on the manifold M. in the sense that they have the same (analytical) solution set for (q, v) on M. . Therefore the existence and uniqueness of solutions of the DAE (2.16) can be seen by considering vector field (2.17) on the manifold M.. To this end suppose the initial value (qo, Vo, A0) satisfies (2.18) and (qo, v 0 ) (E M.. Through the initial value ( q 0 , v 0 ) (at t0 = 0 ), ODE (2.17) has a solution ( q ( t ) , v ( t ) ) . Let A = A(q(i),v(t)) be defined by (2.18). We claim that (q(i), v ( t ) , \(t)) is a solution of DAE (2.16). A simple calculation shows that (q(t), v ( i ) , A(t)) satisfies (2.16b). So we only need to verify that g(q) = 0. This is implied by the fact that ^g(q(0) = G(q)v' + z(q(<),v(t)) = G M - 1 ( / - G T A ) + z(q,v) = 0 and g(q 0 ) = 0, | g ( q ( < 0 ) ) = C(qo)v 0 = 0. Chapter 2. Index Reductions of DAEs and Reformulations 25 In general consider the DAE F(i,x,x') = 0 (2.20) Definition 2.4 ( u n d e r l y i n g O D E ) An ODE is called an underlying ODE of DAE (2.20), if it admits all solutions of (2.20). For a given DAE neither the underlying ODE nor the invariant manifold are unique. For instance, DAE (2.20) has invariant manifolds M. defined by (2.19) and the one defined by Mi = { ( q , v ) : G(q)v = 0} Stabilization can be made on any of the invariant manifolds. Also we can construct different underlying ODEs for a given DAE. Different index reduction formulations as well as stabilization techniques will result in different reformulations of vector fields [12, 15, 17, 25, 31]. All those formulations admit solutions of the original DAE. But the stability and other dynamical properties are not necessarily the same. Those different properties can directly affect the behavior of numerical solutions. In order to investigate the better reformulations for a given DAE, we start our study with an ODE which has an invariant manifold. This ODE can be any of the underlying ODEs of the original DAE. We present a criterion for any other vector fields which share the solution with the ODE on the manifold. According to the definition of underlying ODEs all those vector fields which satisfy the criterion will be underlying ODEs of the original DAE. Then starting from the general form of underlying ODEs we analyze the stability of the invariant manifold and therefore a better reformulation can be obtained. Consider an ODE z' = f(z) (2.21) on an invariant manifold M. = {(z) £ i ? n , h ( z ) = 0 } , where f and h are smooth vector valued functions. For any other ODE defined on M. we have the following conclusion. Chapter 2. Index Reductions of DAEs and Reformulations 26 P r o p o s i t i o n 2.1 A vector field z' = F(z) (2.22) admits all solution of (2.21) on M. if and only if ~P\M = f \M Proof If ¥\M = f \M, for any solution z(t) of (2.21) on M, z'{t) = f (z(i)) = F(z(*)). So z(t) is a solution of (2.22); Suppose (2.22) admits all solutions of (2.21) which lie on M.. Then for any (z 0 ) G A4, we have a solution of (2.21) z(t) € M., t G I such that z(£ 0 ) = z 0 . By assumption z(2) is also a solution of (2.22), so we have z'{t) = F(z(t)) f (z(<)) = z'(<) = F(z(t)), especially for t0: f (z 0 ) = F(z 0 ). So F | ^ = f \M t E I. Therefore n As an immediate conclusion of P r o p o s i t i o n 2.1, the following ODEs share the same solutions on M. with (2.21): z/ = f ( z ) - 7 j F ( t , x ) h ( z ) (2.23) where 7 is a constant and F(t, z) is a matrix and we call F the stabilization matrix. ODE (2.23) is a general form of reformulations in which a variety of choices of parameter 7 and matrix F can be made. For (2.22) we focus on the reformulations in which the manifold M. is stable. As pointed out in [12], all the index reductions and stabilizations currently used in numerical solutions of DAEs have the stability property analytically. But the numerical behavior can be different for various problems. By analyzing the stable reformulations we can find suitable reformulations which can be attractive for numerical discretization. 2.2.2 S t a b i l i t y of a D A E and its reformulations Stability (in the Lyapunov sense) is often defined for an equilibrium state solution of an ODE. It can be extended to the stability for an invariant set of an ODE. In fact, the Chapter 2. Index Reductions of DAEs and Reformulations 27 orbit of an equilibrium solution is a special case of an invariant set. Consider an ODE z' = f (z) , (2.24) and suppose M. is an invariant set of the ODE (2.24). M. is called stable if for any given e > 0, there is a 8 — 8(e,t0) > 0 such that dist(z(t),M.) Similarly we call M. asymptotically < e for t > i 0 , if dist(z0, M.) < 8. stable if in addition lim dist(z(t),M) = 0 Stability for DAEs is studied in [67] although no precise definition of stability for DAEs was given. Similar to the stability of solutions of ODEs, we call a solution of a DAE stable, if solutions of the DAE with small perturbations of initial value remain nearby the solution of the unperturbed DAE. Special care must be taken in the DAE case because small perturbations of consistent initial value may end up with inconsistent initial values. Clearly, if an ODE which admits all solutions of the DAE is (asymptotically) stable the DAE is (asymptotically) stable as well. The converse is not true, however. In fact those ODEs may have different stability properties. We therefore classify the instability in a reformulation into two sources: one is due to the DAE itself and the other arises from the reformulation. An instability due to a DAE is inherent to that DAE and we can't change it by reformulations into an underlying ODE. But an instability due to the reformulation can be improved by selecting a better underlying ODE. For the stability study of the invariant manifold we can construct a certain Lyapunov function. Suppose the invariant manifold M of the vector field (2.24) is given by M = {(z) : h(z) = 0} (2.25) We can construct the Lyapunov function L = h T ( z ) h ( z ) . Then the total derivative of L with respect to (2.24) is i|(2.24) = 2h T (*,x)h z (z)f(z) (2.26) Chapter 2. Index Reductions of DAEs and Reformulations 28 If the total derivative in (2.26) is negative definite then M. is asymptotically stable. T h e proof can be carried out in a similar way as Lyapunov stability theory for an equilibrium state in differential equations. The property of L in (2.26) provides a useful clue when we construct our stabilizations in (2.23). Combined with other considerations we will end up with new stabilizations and index reductions for higher index problems. 2.2.3 S t a b i l i z a t i o n s of U n d e r l y i n g O D E s We now come to the stabilizations of underlying vector fields. With the aid of the Lyapunov function considered above a general form of underlying vector fields with asymptotically stable invariant manifold is given. Also this result can be applied to the stabilizations of any semi-explicit higher index (> 2) DAEs of Hessenberg form. Consider a vector field z' = f(z) (2.27) 0 = h(z) (2.28) with an invariant set M. defined by We propose the following general form of stabilizations for vector field (2.27): z' = f ( z ) - 7 F ( « , z ) h ( z ) (2.29) where F = D(HD)-1 (2.30) or F = D with D(t,z) (2.31) smooth such that HD is nonsingular. If (2.31) is used we also require HD to be symmetric positive definite for all ( t , z ) . Also the smallest eigenvalue of HD is Chapter 2. Index Reductions of DAEs and Reformulations 29 bounded from below by a positive number. For choices of F in (2.30) or (2.31) we have the following stability theory. T h e o r e m 2.1 Suppose the stabilization is symmetric positive matrix F is chosen such that the matrix definite with the smallest HF(z) eigenvalue A(z) > Ao > 0 for all z. Assume that there is a positive number 70 such that ||#f(z)||<7o||h(z)|| (2.32) for all z near AA\. Then the invariant manifold AA (2.28) is asymptotically stable in the vector field (2.29) if 7 > 7o/A 0 . Proof Consider the Lyapunov function V(z) = | h r h ( z ) . Then ^||h(z(*)|| 2 = h(z(t))Tfh(z(t)) h(z(t))TH(z(t))(fo->yF(z(t))h(z(t)) = h(z(t))Tm0--fhTHFh < < < -(7A0-7o)hTh 0 Therefore the invariant manifold is asymptotically stable. R e m a r k 2.1 Recall that Hf(z) • = 0 when h(z) = 0 so the condition (2.32) in the theorem can be realized. On the other hand we cannot remove the condition (2.32). following example shows that even though all the assumptions hold, the invariant manifold can be unstable. E x a m p l e 2.1 •a-Oftsp The of the theorem except (2.32) Chapter 2. Index Reductions where z = (zi,z2)T, of DAEs and 30 Reformulations r2 = z\ + z\ and h(z) = r 2 — 1. The instability of the manifold M. = {z, h(z) = 0} is implied by the fact that h = ct3/2 where c is a constant. R e m a r k 2.2 From the proof of the theorem above the invariant following monotonic property in the stabilized vector fields manifold Ai has the (2.29): ||h(t + a ) | | < | | h ( t ) | | for any a > 0 and t. For index-2 DAEs the monotonic property of the invariant fold also hold in Baumgarte's stabilization. For index-3 DAEs, however, the property does not necessarily hold in Baumgarte's differences between Baumgarte's 2.3 and our stabilization. mani- monotonic This is one of the main stabilizations. I n d e x r e d u c t i o n a n d stabilization for i n d e x - 2 D A E s For higher index (> 2) DAEs it is desirable to reduce the index as far as we can. The index reduction can always be achieved by differentiating the constraints (or the whole DAE in general) one or more times. However, the resulting DAEs(or ODEs) may have different stability property from the original DAEs. To illustrate the possible different stability properties, consider the following index-2 DAE example. E x a m p l e 2.2 1 xx 1 x 2 = 1 2 = 1 ~2X2 0 = where xi,x2 c-* - x2X x2 ' (2.34a) XlX (2.34b) e _i + x\x2 (2.34c) and A are unknown variables, x2 > 0. If we reduce the index by differenti- ating the constraint (2.23c) with respect to i, we have: 0 = - e _ t + x2 l-x-i H x2\\ +x-i(—-X2-xr\ Chapter 2. Index Reductions of DAEs and Reformulations 31 = -\{x\ + x\) Recall that x^ > 0, so A = 0. Substituting this in (2.34a) and (2.34b) we have an O D E for (x1,x2): x1 — -x\ H 2 x'2 = (2.35a) X2 --x2 (2.35b) We can verify that the ODE (2.35) has the general solution: **> = (*<°> + ^ 5 j ) e i - 5) s 2 (i) = x 2 (0)e-2 (2 3fa) ' (2.36b) where (a;i(0), 0:2(0)) is the initial value with x2(0) > 0. We can verify that DAE (2.34) has the solution: 5i(t) = - e " 2 , x2(t) = e~2 Consider the stability of the original DAE (2.34). The solutions of DAE (2.34)(we only consider xi,x2) are those among the general solutions of the ODE (2.35), such that the initial value (xi(0),x2(0)) lies on the manifold: M = {(iCi,X2) : r |l +X-LX2 = 0,X2 > 0} Therefore the general solutions of DAE (2.34) are given by So (xi(t),X2(t)) is asymptotically stable as a solution of DAE (2.34) but unstable as a solution of ODE (2.35). We can easily see that AA. is stable but not asymptotically stable in (2.35). The reformulation can bring additional instability into the equations and there could be a serious problem especially for nonlinear DAEs. Note that if we consider the Chapter 2. Index Reductions of DAEs and Reformulations 32 variational problem of (2.34)(i.e.the linearization about the exact solution) we will have an unbounded coefficient linear system. The variational problem of (2.34) is unstable. In reformulations and index reductions for DAEs certain aspects, like stiffness of the resulting reformulation and expense of the evaluation can has to be considered. Analytically Baumgarte stabilization can be understood as adding " correction" terms into the generalized applied force so that the correction forces the trajectories back into the (invariant) manifold whenever the trajectory is away from it. However the way or the geometrical direction Baumgarte's stabilization drives the trajectories is not necessarily optimal. In this section we will study reformulations for index-2 DAEs. A stabilization technique improving Baumgarte's stabilization for semi-explicit Hessenberg form DAE is proposed. We consider the reformulations and index reduction for the semi-explicit index-2 DAEs. We can write Baumgarte's stabilization as a special case in the general form of stabilized underlying ODEs (2.29). To this end we first derive an underlying ODE by differentiation for a given DAE. For convenience we only consider the autonomous case: x' = f(x) + £ ( x ) y (2.37a) 0 = g(x) (2.37b) where x, y are vectors not necessarily of the same dimension, G(x)B(x.) and G(x) = dg/dx.. is nonsingular For DAE(2.37) the set of consistent initial values of x is given by M = {x € Rnv | g(x) = 0} (2.38) If we differentiate (2.37b) with respect to t, we have 0 = G(x)[f(x) + 5 ( x ) y ] (2.39) Substituting y in (2.37a) we have an underlying ODE: x' = F 0 (x) (2.40) Chapter 2. Index Reductions of DAEs and = f(x) - Reformulations 33 BWlGWBW-iGWW The manifold M. in (2.38) is invariant in (2.40). This can be seen from the fact that any solution x*(i) of (2.40) satisfies ^g(x*(t)) = G(x*(t))x.*'(t) = 0. Baumgarte's stabiliza- tion for (2.37) is an index-1 DAE: x' = f(x) + £ ( x ) y (2.41a) 0 = G(x)[f(x) + 5 ( x ) y ] + 7g(x) (2.41b) where 7 > 0 is a constant. The underlying vector field for x has the form: x' = F7(x) = f(x)- J B(x)[G'(x)i?(x)]- 1 [G'(x)f(x) + 7 g ( x ) ] Notice that F^\M = FQ\M f ° r ai (2.42) l 7 > 0, by P r o p o s i t i o n 2.1 the vector field (2.42) admits all solutions x(t) of DAE (2.37). It is worthwhile to consider the geometrical property of the stabilization term in (2.42). Although both G(x)B(x) and G ( x ) G T ( x ) have inverses, 5 ( x ) and GT(x) in (2.37) can be different. Taking ny = 1 for example, the existence of the inverse of G(x)B(x.) the vector B(x) means does not lie on the tangential plane of M. at x while the vector GT(x.) represents the normal direction of Jli at x. Furthermore G T (x) and -B(x) can be different in magnitude [12]. Those considerations are helpful when we choose an underlying ODE for the DAE to integrate numerically. The general form of stabilized underlying ODEs (2.29) for the index-2 DAE (2.37) for x(<) has the form: x' = F 0 ( x ) where S'(x) = SnxXny(^) 7 5-(x)g(x) is a smooth matrix and 7 > 0 is a constant. (2.43) Chapter 2. Index Reductions of DAEs and Reformulations 34 Multiplying (2.43) by G(x) and remembering that g' = G(x)x' and G(x)F0 — 0 we have ^ + 7G(x)5(x)g(x) = 0 (2.44) The purpose of the choice of matrix S(x) is clear: to make (2.44) asymptotically stable. For the Baumgarte technique S(x) = B(x)[G(x)B(x))~1 and G(x)S(x) = I. Therefore (2.44) is asymptotically stable for any positive 7. However, when -B(x) in DAE (2.37) is less pleasant [12], we may find a better reformulation instead. One of the alternatives is S(x) = GT(x)[G(x)GT(x)}-1 (2.45) S(x) = D(x)[G(x)D(x)}-1 (2.46) or in general where D(x) is a matrix such that G(x)D(x) S in (2.46) we have G(x)S(x) has a (bounded) inverse. For all choices of = I and M. is asymptotically stable. If we choose D = B then (2.46) is the Baumgarte stabilization. The vector GT(x) represents the normal direction of the manifold Ai in (2.38). When we choose S(x) in (2.45), the correction term 7»S'g drives the trajectories of (2.43) back into M. along the normal direction which is the shortest way to reach M.. So the choice of (2.45) is particularly attractive. The drifts in the numerical integration of (2.45) can be much smaller. E x a m p l e 2.3 Consider for 0 < t < 1 [12] x[ = (2 - t)vy + qi(t) x 2 = {v - l)y + 92(0 0 = (t + 2 ) n + (t2 - i)x2 + r(t) with £i(0) = ^2(0) = 1 and y(0) = —0.5. Here v > 1 is a parameter and qi(t) = (1 + z/)e* Chapter 2. Index Reductions of DAEs and 35 Reformulations v ~ 1* ft(0 = 1 + 2 T T C ' r(t) = -(i2 + t-2)e* Here we use the backward Euler method and the step size h = 0.01. We take v = 1000 so hv{= 10) is not small (the stability restriction needs h to be decreased further [12]). 'Error 1' and 'Drift 1' in Table 2.1 refer to Baumgarte. 'Error2' and 'Drift2' refer to t h e stabilization with S = GT[GGT}'1 7 0. 1. 10. 100. 1000. l.e+8 Error 1 .19e-2 .56e-3 .91e-3 .27e-4 .13e+42 .92e+74 in (2.43). Drift 1 .85e-2 .49e-2 .29e-3 .93e-8 .45e+39 .45e+58 Error2 .20e-2 .lle-2 .56e-4 .14e-4 .14e-4 .14e-4 Drift2 .85e-2 .49e-2 .31e-3 .39e-5 .40e-7 0 Error3 .19e-2 .25e-4 .14e-4 .14e-4 .14e-4 .14e-4 Drift3 .85e-2 .10e-3 .12e-5 .12e-7 .13e-9 0 Table 2.1: Comparison of Baumgarte and two improved methods We see that with Baumgarte's method small values of 7 do not lead to good results while large 7 leads to disastrous ones. In contrast, using (2.45) the error decreases (or does not increase) monotonically in 7 and, moreover, the range in 7 which yields acceptable errors is much larger. In this sense the stabilization (2.45) provides a much improved situation. We can safely choose 7 as large as needed to counterbalance F0 (there is no need to choose 7 larger). In (2.45) there is an evaluation of the inverse of G , (x)G T (x) (or at least the solution of linear systems involving GGT as well as GB). Sometimes this can be expensive. The main concern about the reformulation is the stability of the invariant manifold (2.38). For the stabilizations of (2.43) the matrix S could have different form other than (2.46). In general we have the following stability theorem for the reformulation (2.43). Chapter 2. Index Reductions of DAEs and Reformulations T h e o r e m 2.2 If the matrix S in (2-43) is chosen such that G(x)S(x) positive definite and the smallest eigenvalue of GS 36 is symmetric satisfies A^(x)>ao>0 where CTQ is a constant then the invariant manifold (2.47) (2.38) is asymptotically stable in (2.43). Proof We consider the Lyapunov function L(x) = | g T ( x ) g ( x ) for (2.43). Then ^1(2.43) = - 7 g T ( x ) G ( x ) S ' ( x ) g ( x ) Noticing that G(x)S(x) (2.48) is symmetric positive definite and satisfies (2.47) we have ^1(2.43) < - 7 0 o g T ( x ) g ( x ) = -jcr0L By a comparison theorem in [60] the invariant manifold is asymptotically stable. (2.49) • For the various choices of S in (2.46), including the Baumgarte stabilization, the condition of Theorem 2.3 is automatically satisfied with GS the identity matrix. We can consider another choice of S: x' = F 0 ( x ) - 7 G ' : r g ( x ) (2.50) \s G has full row rank GGT is indeed symmetric positive definite. For the particular case need to verify the condition (2.47). We have tested the reformulation (2.50) with the pie 2.3. 'Error3' and 'Drift3' in Table 2.1 refer to the stabilization (2.50). With the parameter setting and same scheme the numerical results show (2.50) is better mgarte stabilization. And it is even somewhat better than (2.45). Testing to notice that the stabilized underlying ODE has a close relationship ^cted invariants methods mentioned earlier. Our stabilization matrix has a Chapter 2. Index Reductions of DAEs and Reformulations 37 similar form with the projection matrix. The main advantage using stabilization matrix lies in the fact that the stabilization matrix determines the stability property of the invariant manifold and thus has a potential to end up with a better reformulation. The general reformulations (2.43) with the choices of (2.45) or (2.46) studied above can also be obtained by introducing additional Lagrange multipliers. We now reformulate DAE (2.37) as x' = f (x) + 5 ( x ) y - Q(x)fi (2.51a) 0 = G(x)(f(x) + £ ( x ) y ) (2.51b) 0 = g(x) - a P ( x ) / i (2.51c) g(x(0)) = 0 where a > 0, P(x) is a ny x ny nonsingular smooth matrix and G(x)Q(x) has an inverse for any x. Different choices of matrices P and Q correspond to different stabilizations. For example, if we let P = GGT, x' = 0 = Q = GT then we have f(x) + 5 ( x ) y - - G , T ( x ) ( G G , T ) - 1 g a G(x)(f (x) + B(x)y) (2.52a) (2.52b) g(x(0)) = 0 If we eliminate the algebraic component y then (2.52) is the reformulation (2.45) (with a = I / 7 ) . When we take P = 7, Q — GT we have the reformulation (2.50) x' = 0 = f(x) + 7 i ( x ) y - - G T ( x ) g a G(x)(f(x) + fl(x)y) g(x(0)) = 0 For a = 0, (2.51) is then a projected invariants method [12]. (2.53a) (2.53b) Chapter 2. Index Reductions of DAEs and Reformulations 38 T h e o r e m 2.3 DAE (2.51) has index 1 for a > 0 and index 2 for a = 0. The multiplier additional JJ, in (2.51) has the exact solution ji(t) = 0. So DAE (2.51) has the same solution ( x , y ) as DAE (2.37). Proof For a = 0, notice that dg/dx 0 = G(x) and differentiate (2.51c). We have = G(x)x' = G(x)(f(x) + 5 ( x ) y - Q ( x ) A i ) = -G(x)Q(x)// This implies /i = 0 as G(x)(5(x) has an inverse for any x. For a > 0, (2.51c) becomes M = -P-1(x)g(x) a (2.54) So DAE (2.54) has index 1. Now differentiating (2.54) we get fi'(t) = - ^ ( P - 1 ( x ) ) g ( x ) + ip-1(x)G'(x)(f(x) + a at a ^(p-1(x))P(x)-ip-1(x)G(x)g(x) at a fi(0) = i P " 1 ( x ) g ( x ( 0 ) ) = 0 gives fi(t) = 0. JB(x)y-g(x)) V • The stabilization can be obtained by suitably choosing the matrices P ( x ) , Q(x) and the parameter a in (2.51). Analytically smaller a corresponds to damping the drift faster but smaller a can also bring stiffness into the system. So far it is not clear how to choose the parameter a in practical computation. By choosing suitable stabilization matrices P and Q, large range of /J, can be expected as shown in Example 2.3. Usually the reformulation (2.43) with (2.46) is used and this includes Baumgarte stabilization as a special case. The difference is that in Baumgarte's stabilization only the parameter is Chapter 2. Index Reductions of DAEs and Reformulations 39 chosen while in our stabilization we may have different stabilization matrices to choose which can provide a better reformulation. And the choice of stabilization matrices can be made from the point view of geometric as well as the stability considerations. 2.4 S t a b i l i z a t i o n for I n d e x - 3 D A E s Consider the index-3 DAE x" = 0 = f(x,x') + £ ( x , x ' ) y (2.55a) g(x) (2.55b) with G(x)i?(x) nonsingular. The index-3 DAE (2.55) is more complicated than index-2 DAEs. Stabilizations for the index-3 DAE can be achieved in two ways. One is regarding the index-3 DAE as a (underlying) vector field on an invariant manifold M. like in index-2 DAE cases. Then reformulate the underlying vector field so that M. has a nice stability property. The other reformulations of index-3 DAEs can be completed in two steps: first reduce the index-3 DAE to an index-2 DAE by either projected invariant or Baumgarte stabilization, then use the reformulations for index-2 DAEs developed in the previous section. To avoid repetition we just focus on the key difference regarding the stability of invariant manifold in the reformulations for index-3 DAEs. Unlike the index-2 cases, the stability of the invariant manifold cannot be guaranteed for certain similar reformulations. Let {* ) \ gx(x)x' j We can write the underlying ODE of (2.55) which is obtained by two differentiations of (2.55b) as z' = f(z) Chapter 2. Index Reductions of DAEs and Reformulations 40 with an invariant manifold M. defined by 0 = h(z) (2.56) where h : U C K2n -> K2ny and f0 : U C Tl2n -»• ft2n are assumed to be smooth, and if(z) = h z ( z ) has a full row rank. The function fo(z) in (2.27) has the form: X f (z) = | I \ f ( x , x ' ) - P (gxi?)"1 [gxx(x',x') + gxf] / It is easy to verify that H(z)f(z) = ( ^ ) (2-57) Thus the condition (2.32) holds with 70 = 1. Therefore the invariant manifold M. (2.56) has a monotonic property in the stabilized underlying vector fields (2.29). While in Baumgarte's stabilization g satisfies the equations g" + 7 i g ' + 72g = 0 and usually ||h|| is not monotonic decreasing. For index-3 DAE (2.55) the analogue to stabilization (2.51) of index-2 DAEs is x" = f(x,x') + J B ( x , x / ) y - 2 7 G T ^ - 7 2 G ' T 7 / (2.58a) 0 = G ( / + By) + z(x,x') (2.58b) 0 = Gx'-Pi(x,x')^ (2.58c) 0 = g-P2(x,x> (2.58d) g(x(0)) = G(x(0))x'(0) = 0 where z ( x , x ' ) = G x ( x ) ( x ' , x ' ) , G = 5 g ( x ) / 5 x , P i ( x , x ' ) and P 2 ( x , x ' ) are ny x ny nonsingular smooth matrices and 7 > 0 as before. Chapter 2. Index Reductions of DAEs and Reformulations 41 Similar to the index-2 case we can prove that the solutions for fj,(t) and rj(t) in (2.58) are fi(t) = r/(t) = 0. So (2.58) and (2.55) share the same solution for ( x , y ) . For (2.58) we have dg dt dt2 = G(x)x' = Gx(x)(x',x') + C?(x)x" = z(x, x') + G(f + By- 27GT^ - 2 7 = ~21{GGr)P''dft-1\GGT)p-^ = -27(GGr)p-1^-72(G'GT)p-1g GTV) So g satisfies the following equation g" + 27(GGT)Pr1g' + 12{GGT)P^g = 0 (2.59) There are many choices of P\ and P2 in (2.58) so that (2.59) is asymptotically stable. For example if we choose Pi = P^ = GGT, then (2.59) has the differential equation which is the same as the Baumgarte reformulation: g" + 2 7 g ' + 2 7 g = 0 (2.60) (2.60) is asymptotically stable for any 7 > 0. Stabilization (2.58) with Pi = P2 = GGT is different from Baumgarte stabilization when B / GT. If we choose Pi = / , P 2 _ 1 = GGT then (2.59) becomes g" + 2 7 ( G G T ) g ' + 2 7 (GGT)2g = 0 (2.61) This is analogous to the stabilization (2.52) for index-2 DAEs. A numerical example for those two stabilizations with the choices of Pi = P 2 = GGT and Pi = / , P 2 _ 1 = GGT is given by a modification of example 2.3: x'[ = (2-t)vy + qi(t) Chapter 2. Index Reductions of DAEs and x'2' = 0 = (u-l)y Reformulations 42 + q2{t) (t + 2)xi + (t2 - 4)x 2 + r(t) where qi, q2 and r are same as in Example 2.3. We use the backward Euler method, h — 0.01 and v = 1000. Errorl and Driftl refer to Baumgarte, Error2 and Drift2 refer to P1 = P2 = GGT and Error3 and Drift3 refer to i \ = J and P2 = GGT. We can see that for larger 7 ( = 1000) Baumgarte technique produces an unreasonable result. But the stabilizations we proposed still work well. 7 0. 1. 10. 100. 1000. Errorl .82e-2 .80e-2 .16e-l .19e-l .20e+65 Driftl .72e-2 ,63e-2 .63e-3 .79e-5 .lle+61 Error2 .82e-2 .80e-2 .73e-2 .71e-2 .71e-2 Drift2 .72e-2 .63e-2 .68e-3 .16e-4 .81e-7 Error3 .82e-2 .72e-2 .71e-2 .71e-2 .71e-2 Drift3 .72e-2 .23e-3 .68e-5 .47e-6 .45e-7 Table 2.2: Baumgarte and Baumgarte improved for an index-3 DAE The stability analysis for (2.61) is not as easy as (2.44) which is for the index-2 DAE case. T h e eigenvalue criterion for stability can't be used here as (2.61) has time dependent coefficients. Although GG is positive definite, (2.61) is not always stable as the following example shows The ODE (2.62) has a solution g = it*. Definitely we can use the reformulations with the choice of P\ = P2 = GGT in (2.58). Note that in (2.60), although it is asymptotically stable for 7 > 0, the damping is not necessarily monotonic. The stabilization for index-3 DAEs can also be carried out by two steps: reducing (2.55) to an index-2 problem first, (for example we can use Baumgarte or projected Chapter 2. Index Reductions of DAEs and Reformulations 43 invariants method), then applying the stabilizations to the reduced index-2 DAE. For instance using Baumgarte stabilization we can reformulate (2.55) to an index-2 DAE: x' = v (2.63a) v' = f(x,v) + 5 ( x , v ) y (2.63b) G(x)v + 7 g (2.63c) 0 = Next we apply the stabilization proposed in the previous section to (2.63). If the reformulation (2.53) is used we have: x' = v - ( ( G x v ) T + 7 G T ) (i v' = f(x,v) + 0 = ( G x v + 7 G , )v + G ( f ( x , v ) + J B(x,v)y) (2.64c) 0 = g' + 7g (2-64d) (2.64a) fi(x,v)y-G'T/i (2.64b) g'(x(0),v(0)) + 7 g ( x ( 0 ) ) = 0 where C x = dG(x)/dx.. Notice that in (2.64) the algebraic variables // and y can be eliminated using (2.64c) and (2.64d). In general the reformulations or stabilizations for index-3 DAEs are more difficult than the index-2 case. By stabilizing the invariant manifold M 2 = { ( x , v ) : 0 = gx(x)v} (2.65) alone usually is not enough as it is equivalent to an index-2 problem. The different reformulations discussed so far guarantees the asymptotic stability for the invariant manifold analytically. It can be expected that our reformulations studied in this chapter work well for a variety of problems. Chapter 3 Stabilizations of D i s c r e t e D y n a m i c s Formulations In the previous chapter we have studied various reformulations and index reductions for a given DAE. For the semi-explicit index-2 or index-3 DAEs of Hessenberg form, by differentiation we can rewrite the DAE as a underlying ODE: z' = f(*,z) (3.1) with an invariant manifold M = {z, 0 = h(z)} (3.2) The stability of the invariant manifold M. (3.2) may directly affect the efficiency of a numerical integration when we are interested in the solutions of (3.1) on the manifold (3.2). We stabilize the invariant manifold M. by reformulating the vector field (3.1) to the following: z' = f ( * , z ) - 7 F ( z ) h ( z ) (3.3) where 7 is a positive number. The stabilization matrix F in (3.3) is chosen so that M. is asymptotically stable. The numerical integration of (3.3) depends on the choices of Baumgarte parameter 7 and the stabilization matrix F. There are many studies [46] on the optimal choice for the Baumgarte parameters 7. Numerical examples in Chapter 2 show that the numerical behavior depends on the matrix F. In [8] we point out that the optimal choice of the parameters depends on the step size and the numerical scheme as well. This is discussed and demonstrated further in [7]. Moreover, the analytical 44 Chapter 3. Stabilizations of Discrete Dynamics Formulations 45 asymptotic stability of the invariant manifold cannot guarantee that the numerical solutions have small drifts. The stability property of the invariant manifold M. may appear different after the discretization and the numerical solutions may still have large drifts. To see this consider a simple index-2 DAE example. E x a m p l e 3.1 z' = 3t2 + y (3.4a) 0 = z-t3 (3.4b) Using the Baumgarte stabilization for (3.4): = (*' - 3t 2 ) + 7(2 - t 3 ) = y+ 0 <y(z-t3) So we have the underlying ODE of (3.4) for z: z' = 3t2 - 7(2 - t3) (3.5) with the invariant manifold given by 0 = z-t 3 (3.6) In the underlying ODE (3.5) if 7 = 0 the invariant (3.4b) is stable but not asymptotically stable. For any positive 7, the asymptotic stability of M. in (3.5) is implied by the fact: ^{z-t3)2 = -1{z-t3f (3.7) Now consider the discretization of ODE (3.5) with the initial value (to, ZQ) = (0,0). Using the forward Euler method with step size h, we have zn+1 = zn + h(3t2n-7(zn-t3n)) = (1 - h~f)zn + 3h(nh)2 + h-/(nh)3 (3.8) Chapter 3. Stabilizations of Discrete Dynamics Formulations 46 When 7 = 0 the solution (3.8) has the form: zn+1=t3n+1-h3(n + l)(^n + l) (3.9) For 7 > 0 we require that |1 — hj\ < 1, otherwise the solution (3.8) will be exponentially large. Take 7 = \jh (It can be shown that in this example there are no significant differences for the drifts among those 7 with |1 — hj\ < 1 ). Then (3.8) reads zn+1=t3n+1-h3(3n + l) (3.10) Clearly for all 7 > 0, h ( t n + i ) = zn+\ — t^ + 1 = — h3(3n + 1) will grow arbitrarily for n large enough (i.e. t —> 00). So the manifold M. is not stable in the difference equation (3.8). R e m a r k 3.1 cretization, 1. The simple stable invariant manifold can become unstable after disand the drift can be quadratically 2. The asymptotically growing; stable invariant manifold can become unstable after discretiza- tion, and the drift can be linearly growing. The observations above motivate us to seek more efficient stabilizations for the integration of underlying ODEs with the invariant manifold. A discretization results in difference equations or a discrete dynamical system. Within the difference equations a new asymptotic property could appear and some stability properties of the vector field could be lost too [98]. In order to restore the stability for the invariant manifold we have to reformulate the difference equations somehow without diluting the scheme's other desirable properties such as accuracy. We call this procedure a post-stabilization. We start with a discussion of the difference equations in Section 1. Using the comparison theory we consider the stability of an invariant manifold in difference equations in Section 2. We prove that the manifold is uniformly stable in the difference equations Chapter 3. Stabilizations of Discrete Dynamics Formulations 47 with post-stabilization. The drifts of the numerical solutions are directly related to the stability of the manifold in the difference equations. In Section 3 we analyze the order of the drift in the difference equations with a post-stabilization. When the stabilization matrix is chosen in a certain way the drift is 0(/i 2 ( p + 1 ') if a scheme with local truncation error 0(hp+1) is used. This is a stronger result then the ones reported in [8] or [24]. In Section 4 we propose stabilization algorithms for holonomic constrained multibody dynamics. Among the proposed post-stabilization methods, one variant which we call PStab is our method of choice and it has been used in our numerical code for constrained multibody dynamics. In Section 5 we consider the stabilization scheme for non-holonomic constrained systems. We present other methods like coordinates projection methods in Section 6. Our stabilization for difference equations corresponds to one iteration of a certain projection method if a particular weighted norm has been used. Our stabilization methods restores in the difference equations the stability property of the invariant manifold and therefore it ensures numerical solutions with very small drifts. A variety of numerical examples in the next chapter show that PStab is efficient and robust. 3.1 D y n a m i c s of Difference E q u a t i o n s Difference equations appear in the study of discretization methods for differential equations. They also appear as natural descriptions of observed evolution phenomena because most measurements of time evolving variables are discrete and such difference equations are in their own right important mathematical models. Some results in the theory of difference equations have been obtained as more or less natural discrete analogues of corresponding results for differential equations. This is especially true in the case of the Lyapunov theory of stability [60, 1, 61, 98]. Nonetheless, the theory of difference equations is much richer than the corresponding theory of differential equations. For example, Chapter 3. Stabilizations of Discrete Dynamics Formulations 48 a simple difference equation resulting from a first order differential equation may exhibit a phenomenon often called the appearance of "ghost" solutions or existence of chaotic orbits that can only happen for higher order or higher dimensioned differential equations. Consequently, the theory of difference equations is interesting in itself and it plays an important role in numerical analysis[1, 61, 98]. N o t a t i o n : We will use the following notation in this chapter: • N = { 0 , 1 , • • •} the set of natural numbers including zero • N(a) = {a, a + 1, • • •} where a £ N • N(a, b) = {a, a + 1, • • • , b} where a < b < oo and a,b £ TV Let N be any one of the above sets. A difference equation in one independent variable k £ N and one unknown u(k) is a functional equation of the form f(k,u(k),u(k + l),---,u(k + n)) = 0 (3.11) where / is a given function. The order of (3.11) is defined to be the difference between the largest and smallest arguments explicitly involved. Equation (3.11) is said to be normal if it is of the form u(k + n) = f(k,u(k),u(k + 1),- • • ,u(k + n - I)) (3.12) We will mainly consider a system of difference equations u(fc + l) = F(Jfe,u(ife)) (3.13) where u is a 1 x n vector and F is a given vector valued function. A system of linear difference equations has the form u{k + l) = A(k)u{k) + b(k) (3.14) Chapter 3. Stabilizations of Discrete Dynamics Formulations 49 where A(k) is a given n x n matrix and h(k) is a given n x l vector. Initial value p r o b l e m s A function u(k) defined on Nn, where Nn N(a,b-l+n) ifN = N(a,b-l) N(a) if N = N(a) N ifN = N is said to be a solution of the given n-th difference equation (3.11) on A^ if the values of u(k) reduce the difference equation to an identity over JV. Similarly , for the system of difference equations (3.13) a function u(k) defined on A\ is a solution of (3.13) on N if the values of u(k) reduce the difference equations to an equality over N. The initial value problem is to find a solution of (3.13) with the initial value u(a) = u° where u° is a given value. We denote the initial value solution by u(k, a, u°) Our concern for a given difference system is the ultimate behavior of its solutions. Let us denote by B(v, S) the open ball having center at v and radius 6. For the special case v = 0 we simply use the notation Bg. Suppose (3.13) has u(k) = 0 as its solution. The stability of the zero solution is defined as follows [61]: Definition 3.1 The solution u = 0 of (3.13) is said to be 1. stable if given e > 0 , there is a 8(e,ko) such that for any u° £ B$ the u(k) u(k,k0,U0)eBe. = 2. Uniformly stable if it is stable and 6 can be chosen independent 3. Attractive solution on ko- if there is S(ko) > 0 such that for u° G B$ one has lirm;_+00 u(k) = 0. 4- Asymptotically stable if it is stable and attractive. For a nonlinear system, usually the study of one particular solution's asymptotic property is through the linearization around the solution. Consider the stability of nonlinear Chapter 3. Stabilizations of Discrete Dynamics Formulations 50 difference equations: u fc+1 =A(k)uk + f{k,uk) (3.15) where f : N x Ba —*• Ba and f (k, 0) = 0. If the function f is small in a certain sense and we regard (3.15) as a perturbation of the difference equations u fc+1 = A(k)uk (3.16) then for the stability of (3.15) we have the following theorem [1, 61]: T h e o r e m 3.1 Assume that ||/(fc,u*)||<e||u*|| for all k and all \\uk\\ sufficiently 0 of (3.16) is uniformly asymptotically small, where e > 0 is small. asymptotically (3.17) If the solution uk = stable, then the solution uk = 0 of (3.15) is stable. As a corollary of T h e o r e m 3.1, the Perron theorem holds. Corollary 3.1 (Perron). Consider the difference equations uk+1 = Auk + f(k,uk) where A is a constant matrix with all its eigenvalues inside the unit disk and li m IW-;>II = 0 u-o uniformly 3.2 (3.18) moreover (3.19) ||u|| with respect to k. Then the zero solution of (3.18) is asymptotically stable. S t a b i l i t y of Invariant Manifolds in Difference E q u a t i o n s Consider a discretization of (3.1) by a one-step (explicit) scheme y/c+i =yk + (j)(tk,yk,hk) (3.20) Chapter 3. Stabilizations of Discrete Dynamics Formulations 51 where hk is the (fcth) time step and (3.20) advances the solution of (3.1) from the approximation yk at tk to an approximation y^+i at tfc+i. Using the theory of difference equations we can study various stability and other properties of the scheme (3.20). Here we are mostly interested in the asymptotic property of manifold M. (3.2) in the difference equations (3.20). We now present the definition of stability for an invariant manifold in the difference equation. Suppose M* is an invariant manifold of difference equations (3.13). We denote the ^-neighborhood of M* in u-space by B(M;6) B(A4*,6). = {u,vMJu-v\\<6} (3.21) Definition 3.2 The invariant manifold A4* of (3.13) is said to be 1. stable if given e > 0 , there is a 8(e,ko) such that for any u° G B(M.*,S) solution u(k) 6 B(M*, the e) for all k. 2. Uniformly stable if it is stable and 8 can be chosen independent of ko • Definition 3.2 directly relates the drifts of the numerical solution (3.20) with the stability of M. in the discrete dynamical system. In order for the numerical solution to remain nearby the invariant manifold we have to require that M. be stable in the difference equations (3.20). Obviously a stable invariant manifold M. in the underlying ODE (3.1) has a better chance to be stable in the difference equations (3.20) than an unstable M.. However we cannot rely on that. We therefore propose a stabilization method which reformulates the difference equations so that M. is stable or even asymptotically stable in the resulting reformulation. Consider the difference equations (3.20) arising from a one-step discretization. We Chapter 3. Stabilizations of Discrete Dynamics Formulations 52 stabilize the manifold M. in (3.2) by reformulating (3.20) into the following set of difference equations [8]: zk+1 = zk + (j>(tk,Zk,hk) (3.22a) zk+1 = zk+1 - F(zk+1)h(zk+1) (3.22b) If the solution z(k) of (3.22) is near the manifold M,, we can regard (3.22) as a perturbation to the difference equations (3.20). Then the solutions of (3.22) will stay nearby the solution of (3.20). If the perturbation term in (3.22) (i.e. the correction term in (3.22b)) is smaller than the truncation error of the discretization (3.20), we can expect that the solutions of (3.22) have the same order of accuracy as that of (3.20). In order to give a global estimation of the difference between the solutions of the two systems, we present the discrete Gronwall inequality [1, 61] which is analogous to the well-known Gronwall inequality for differential equations. L e m m a 3.1 (Discrete Gronwall inequality) Suppose that for all k 6 N(a) the following inequality is satisfied k—1 u(k)<p(k) + q(k)J2f(l)u(l) (3.23) l=a Then k—l k—l «(*) < P(k) + q(k) £ K0/(0 II (1 + <?(T)/(T)) l=a for all k e (3.24) T=1+1 N(a). Using the discrete Gronwall inequality we can prove the following estimation theorem. T h e o r e m 3.2 For (f>(tklyk,hk) suppose there are nonnegative functions \(k) and fj,(k) such that \\<l>(tk,yk,hk) - <j>(tk,zk,hk)\\ < A ( f c ) | | y - z | | \\F{zk+1)h(zk+1)\\ < fiik) (3.25) (3.26) Chapter 3. Stabilizations of Discrete Dynamics Formulations 53 Then for the solution y^ of (3.20) and the solution z^ of (3.22) we have the following estimation: k—l k—1 /—1 z ||y* - zfc|| < (||y° - z°|| + £ MO) + E ( l l y ° - °ll + £ M) /=a ?=a j=a k—1 I I (1 + A (r)) (3.27) T=1+1 where y° and z° are the initial values. Proof We can rewrite (3.20) and (3.22) as follows: k-i Yk = y° + 5 X * ' > y > M l=a k—l xk = z° + ^ ( ^ / , z / , / i / ) - i ? ( z / ) h ( z / ) ) l=a So k—l fc—1 +£ f(z,)h(z,) From (3.25) and (3.26) we have ( A: —1 \ k—l lly o -« o ll + E M 0 ) + E A ( 0 l l y » - ^ l l (3-28) Then by the discrete Gronwall inequality (3.27) follows. • The function (f> in the discretization (3.20) is of 0{hk) and so is \{k) in (3.25). The main difference between the solutions y(k) and z(k) lies in the fact that the stability property of M. in the difference equations (3.20) and (3.22) could be different when we vary the choice of the stabilization matrix F. Different choices of F lead to different difference equations. By reformulating the difference equations so that the manifold M. is asymptotically stable we can expect the solution z(k) to stay nearby M. for any k. The stability study for an invariant manifold is more difficult than that of an equilibrium state or other simple invariant sets like periodic orbits in difference equations. Chapter 3. Stabilizations of Discrete Dynamics Formulations 54 There is not too much study reported in the literature on the stability of an invariant manifold for difference equations. Fortunately in our study the invariant manifold M. has a special form, i.e. it is defined explicitly by a certain known function h(z). monitoring the values of the functions h we are better off for the stability of M. for u(k) By But = h(zjt) where Zk is a numerical approximation for an underlying ODE at tk, the next term u(k + 1) = h(z*;+i) is determined by <j> i.e. the discretization scheme, F and u(k) altogether. There is no recursive formulation for u(k + 1) and therefore we cannot write u(&) as a discrete dynamical system. For this reason, the stability theory for difference equations like Theorem 3.1 can not be directly applied to h(z(&)). We have to turn to other qualitative theory in dynamical systems to study the stability property oiM. In the qualitative study for differential equations or difference equations, there are two powerful tools: the comparison theory and the Lyapunov direct method. For our asymptotic study of the values h(z(fc)) we can make a full use of these tools and techniques. For example, one useful comparison theorem in difference equations is the following [61] P r o p o s i t i o n 3.1 Let k £ N(k0), r > 0 and g(k,r) be a function respect to r for any fixed k. Suppose that for k > ko, the y(k + l)<g(k,y(k)) u(k + l) hold. >g(k,u(k)) Then y(k0) < u(k0) implies y(k) < u(k), for all k > ko inequalities, nondecreasing with Chapter 3. Stabilizations of Discrete Dynamics Formulations 55 We now come to the stability theorem for an invariant manifold M. in difference equations (3.22). T h e o r e m 3.3 Suppose h(z) and F(z) are sufficiently smooth near the compact M. defined in (3.2). If the stabilization manifold matrix F in (3.22) is chosen such that \\I-H(z)F(z)\\<P<l where p is a constant and H(z) = h z is the Jacobian matrix, then for sufficiently (3.29) small hk in (3.22), manifold A4 is uniformly stable in (3.22). Proof Since M. = {z : h(z) = 0} is compact, the uniform stability of M. is implied by the following statement: If for any given e > 0, there exists a 8(e) > 0 so that the solution of (3.22) satisfies ||h(zfc)|| < e for any k when hk is small enough and ||h(zo)|| < 8(e). Recall that h(z) and F(z) are sufficiently smooth near the compact manifold M. so we can find a constant M > 0 such that ||h(z - F ( z ) h ( z ) ) | | < ||h(z) - # ( z ) F ( z ) h ( z ) | | + M | | h ( z ) | | 2 (3.30) for all z near M.. Let pi — p -\- {\ — p)/2 then p\ < 1. Now for any given e > 0, we first assume the step-size in scheme (3.22a) is small such that ||h(z + ^ f c , z , ^ ) ) | | < | | h ( z ) | | + « 5 i (3.31) where • i e i l - Pi) n ^ 1 - P^ h = rDm{ '(Wi)"^-} Choosing 6(e) = m i n { e / 2 , (1 — p)/AM}, we claim that ||h(z fc )|| < 6 (3.32) Chapter 3. Stabilizations of Discrete Dynamics Formulations 56 for all k if ||h(z 0 )|| <6. In order to prove (3.32), let us first show that |h(z,)|| < ^ / (3-33) for all k. Using (3.31) we have So (3.33) holds for k = 1. Now let TV = {1, • • •, / } , so that (3.33) holds for all k e N. Then for k < I we have ||h(z f c + i)|| = \\h(zk + <f>(tk,Zk,hk))\\ < ||h(z,)||+^ = \\h(zk-F(zk)h(zk))\\+Si + M | | h ( z , ) | | 2 + *! < ||7 - H(zk)F(zk)\\\\h(zk)\\ < ( / , + M||h(z f c )||)||h(z f c )||+* 1 < = (/>+^)l|h(zfc)||+*i /?i||h(z f c )||+£i Using Proposition 3.1 ( the comparison theorem ) we have llh(zJ+1)n < ^i|h(Zl)ii + E M < HhCzOII+px-A1- pi < | | h ( z 0 ) | | + * i + />i = ||h(z 0 )|| + - ^ 1 - pi 1 - p 1- p < - -+ 4M - 4M Si I-Pi (3.34) Chapter 3. Stabilizations of Discrete Dynamics Formulations 57 So we have (/ + 1) G N and therefore (3.33) holds for any k. Repeating the above arguments we have a similar estimation of (3.34) for ||h(zjt)||: ||h(z fc )|| < ||h(z 0 )|| + - ^ 1 - px < § +§=< (3.35) Finally using (3.30) and (3.29), we have ||h(z,)|| = ||h(z f c -F(z f c )h(z f c ))|| < (p + = M\\h(zk)\\)\\h(zk)\\ Pi\\Hik)\\ Notice that p\ < 1 and we have ||h(zfc)|| < e by (3.35). R e m a r k 3.2 To determine the stability of the linear • system v{k + 1) = (I - H(zk+i )F(ik+1)) v(k) (3.36) is not an easy task as it is a time dependent system in general. the stabilization matrix F such that (3.36) is asymptotically F such that HF = / then the linear system stable. system is uniformly If we can choose (3.36) has a zero coefficient matrix. stability of a constant coefficient linear system is totally determined of its eigenvalues. Our goal is to find by the The distribution Therefore if all its eigenvalues lie inside the unit disk then the linear asymptotically stable. For time-dependent the eigenvalue criterion may not be sufficient. linear systems, however, See the following example [61]. E x a m p l e 3.2 Consider the linear system y(k + 1) = A(k)y(k) (3.37) Chapter 3. Stabilizations of Discrete Dynamics Formulations 58 where / m = o\9-{-l) 7 k 9 + (-l)*7 (3.38) 0 The eigenvalues of A{k) are ± v 2/2 for all k, and they are inside the unit disk. The zero solution of (3.37) is unstable as the general solution of (3.38) has the form: if k is even, and when k is odd. Here C\ and c% are any constants. We can see that in any case a solution will grow exponentially. • By Theorem 3.3 the drifts (away from the manifold Ai) of the difference equations (3.22) will remain fairly small globally provided the step size h is small and (3.36) is uniformly asymptotically stable. In the next section we will consider several poststabilization algorithms for the Euler-Lagrange equations and the ODEs with invariants. The solutions by the resulting algorithms will stay nearby M. and also retain the same order of accuracy as the original scheme z(k + 1) = z(k) + <f>(tk, z(&), hk). 3.3 P o s t - S t a b i l i z a t i o n s for O D E s w i t h Invariant Manifolds We have studied the stability of the manifold M. in difference equations. By suitable choices of the stabilization matrix F in (3.22b) the manifold M is uniformly stable. Therefore we can control the drift to a small amount. An estimation between solutions of the stabilized difference equations (3.22) and the original form (3.20) has been given Chapter 3. Stabilizations of Discrete Dynamics Formulations 59 by Theorem 3.2. In this section we describe the actual construction of the stabilization algorithms for an ODE with an invariant manifold. Then we propose two stabilization algorithms for Euler's equations in descriptor form for constrained multibody dynamics. Consider a stabilized ODE system z' = f(z) (3.39) with an invariant manifold M. defined by h(z) = 0 where h : 7Zn -> Km (3.40) and f : Tln -> 7ln are sufficiently smooth and H(z) = hz(z) has a full row rank. From example 3.1 we see that the analytical stability property of invariant manifold M. defined by (3.40) can play an important role in the discretization. Furthermore we can also consider the monotonicity property of the invariant manifold and this leads to the following definition. Definition 3.3 An invariant manifold M. is said to be decreasing in the vector field (3.39) if for any solution z(t) of (3.39) and a > 0 ||h(z(t + a ) ) | | < | | h ( z ( t ) ) | | (3.41) The decreasing property is an extension of the concept of integral invariant manifold. Recall that an invariant manifold M. of (3.40) is an integral invariant set in the vector field (3.39) if H{z)f(z) = 0 for all z. There are, however, invariant manifolds with decreasing property which are not necessarily integral invariant manifolds. Consider a stabilized ODE system z' = f (z) - 1HT (HHT) ~l (z)h(z) (3.42) with an invariant manifold given by (3.40), where 7 is a positive parameter. Suppose there are positive 70 and r such that ||#(z)f(z)||<7o||h(z)|| (3.43) Chapter 3. Stabilizations for all z G B(A4,r). of Discrete Dynamics Formulations 60 Then the invariant manifold M. of (3.40) is decreasing in (3.39) for all 7 > 70 as we pointed out in the previous chapter. R e m a r k 3.3 For the semi-explicit index-3 DAEs of Hessenberg form, the constant^ (3.43) is 1. So there is no need to choose 7 very large to make M. in decreasing. Theorem 3.3 ensures that the invariant manifold M. is uniformly stable in the difference equations (3.22). We now come to a practical estimate for the possible drifts when the stabilized scheme (3.22) has been used. To this end we start with an order p discretization scheme for ODE (3.39) zk+1 = zk + cf>(tk,zk,hk) (3.44) with the solution zk at tk and the step size hk = tk+i — tk. The post-stabilization scheme (3.22) can be written as z fc+ i = zk + cf)(tk, zk, hk) - F(zk + <f>h(tk, zk, hk))h(zk + cj>(tk, zk, hk)) (3.45) We have the following theorem. T h e o r e m 3.4 Suppose the manifold M. is decreasing in (3.39) and compact in lZn. we choose the stabilization If matrix F such that \\I-HF{z)\\<p<\, (3.46) then for the scheme (3.^5) we have ||h(z,)|| < M0hpk+1 where Mo is a constant of moderate size. Proof (3.47) Chapter 3. Stabilizations of Discrete Dynamics 61 Formulations Suppose h and F are smooth for z near M.. The manifold M. is compact so there exists a positive number M\ such that for any z near M. h(z - F ( z ) h ( z ) ) = ( / - HF{z)) h(z) + Ai(z) (3.48) anc ||ff(z)|| < Mx 11^(2)11 < MJIIMZJII 2 (3.49) Let pi = p + (1 — p)/2. Then 0 < p\ < 1 because ^ < 1. Since the scheme (3.44) has an order p accuracy we can find a constant M2 so that the local truncation error Tk satisfies Tk < M2hl+1 for all k if the right-hand side function in vector field (3.39) is smooth. Furthermore we choose the step size small enough so that T ^ T = 1-p 1 - p\ ( W^ . } , Applying (3.48) for the scheme (3.45) we have h(z f c + 1 ) = (/ - HF(zk+1)) h(z f c + 1 ) + A a ( z , + 1 ) (3.51) The initial condition Zo is consistent so h(zo) = 0. Therefore for i = 0 w e have MxCHMzOH + M x T ^ l a - / ? ) (3.52) Let N denote the number set that for any k £ TV (3.52) holds for all i < k. We know that 0 £ N i.e. N is not empty. For k G N let z(t) — z(t,t^,z^) be the exact solution of (3.39) with z(tk) = Zfc. Then ||h(z* + 1 )|| = \\h(z(tk+1)) + H(zi+1)(zk+1-z(tk+1))\\ < ||h(z(t J f c + 1 ))||+M 1 --n f e + 1 (3.53) Chapter 3. Stabilizations of Discrete Dynamics Formulations 62 Using the decreasing property of M. and (3.51) we have ||h(z f c + 1 )|| < | | / - f r i f ( z J H . 1 ) | | | | h ( z f c + 1 ) | | + ||A 1 (z f c + 1 )|| < /5||h(zfc+1)||+M1||h(z,+1)|| < (^ + M 1 ||h(z f c + 1 )||)||h(z J f e + 1 )|| < (p + M1(\\h(z(tk+1))\\ < (/0 + M 1 (||h(z f c )|| + M 1 r ) ) . ( | | h ( z J f e ) | | + M 1 r f c + 1 ) < (p + \(l = Pl(\\h(zk)\\ 2 + MlTk+1)) - p)) • (\\h(zk)\\ + • (||h(z(t A + 1 ))|| + Min+1) (3.54) (3.55) MlTk+1) + MlTk+1) (3.56) Therefore with consistent initial condition, ||h(zfc+1)|| < ^ +1 ||h(z 0 )|| + E ^ i ^ ^ j=o < Mi -^M2hpk+1 (3.57) Using (3.57) and (3.50) we have M 1 (||h(z Jfe+1 )||+M 1 .r fc+1 ) = 1 1- Pl ~ l Mlx 21-^r<Mf-^~. -^.l-^=12-^ 2M\ 2Pl Pl Therefore the estimation (3.52) holds for any i when we choose the local truncation error r so that (3.50) is valid. So (3.57) holds for any k. • For a semi-explicit index-3 DAE the index reduction can be easily carried out by two differentiations of the algebraic equations. In such an underlying vector field the invariant manifold defined by the constraint equations along with their time derivatives does not have the decreasing property. However the decreasing property of the invariant Chapter 3. Stabilizations of Discrete Dynamics Formulations 63 manifold can be relaxed by choosing the stabilization matrix F such that I-H(z)F(z) = 0 (3.58) If (3.58) holds we have even a better results for h(zfc): T h e o r e m 3.5 Suppose we choose the stabilization M.. If the local truncation error of scheme (3.44) matrix F so that (3.58) holds near zs 0(hp+1) then ||h(z,0|| = O ( / i 2 ( p + 1 ) ) (3.59) Proof The proof can be carried out by two steps. The first step is to prove ||h(z,)||=0(>+1) (3.60) and the second one is to prove (3.59). In fact (3.60) can be proved in similar way to that of Theorem 3.4. Recall that the key step of using the decreasing property in the proof of Theorem 3.4 is the one from (3.54) to (3.55) in order to obtain (3.56). In order to prove (3.56) in the current situation let us consider an auxiliary (stabilized) initial value problem at the integration step tk y' = f(y)-"foHT(HHTy\y)h(y) (3.61) y ( * k ) = %k Let y(t) = y(t,tk,Zk) and z{t) = z(t,tk,Zk) be the exact solutions of (3.61) and (3.39) for t 6 [tk,tk+i\. As we mentioned earlier M. is decreasing in (3.61) i.e. we have | | h ( y ( M , , z , ) ) | | < ||h(y(<*))|| = ||h(z*)|| t € [**,i*+i] Now for the exact solutions y(t) and z(t) we have an estimation in [ifc, ^fc+i] \\z(t,tk,zk) -y(Mjfc,zjfc)|| 64 Chapter 3. Stabilizations of Discrete Dynamics Formulations = || j t k (f (z) - f(y) + 7 o i f T ( M T ) - 1 h ( y ) ) <fe|| < f ||f (z) - f (y)||<fa + 70 f Jtk Jtk ftL\\Z(s)-y(s)\\ds < \\HT{HHT)-lh{y)\\ds + l0\\h(y(tk))\\ Jtk 'tk \\HT{HHTy'\\d& f Jtk where L is a Lipschitz constant. Therefore \z(t,tk,zk)-y(t,tk,zk)\\ < 7 1 rt ° J|h(z f c )|| lkL, 1 — hkL Jtk Jtk f\\HT{HHTYl\\ds for t £ [tk,tk+i]. Let yk+i be the numerical solution of (3.61) at tk+i by using the same discretization scheme as for z (3.44). For step-size h small enough we have ||h(z fc+1 )|| = ||h(y fc+1 ) + H(6k+1){zk+1 - < (||h(y(t J H .i))|H-M 1 -r y ) + yk+1)\\ \\H(ek+l)\\ (\\z(tk+1) - y(tk+1)\\ + M1(TE + ry)) < (\\h(y(tk+1))\\+M1-ry) + \\H(9k+1)\\x < M 1 ||h(z fc )||+ M(ry + Tz) (3.62) where r z and r y are local truncation errors for z and y respectively. From (3.58) we have p = 0 in (3.54). Using the estimation (3.62) we have a similar estimation as (3.56). Similar to (3.57) we have (3.60). From (3.62) and the assumption that the local truncation error is 0(hp+1) we have iih(z,+1)ii = o ( y + i ) Therefore (3.59) can be obtained as follows h(z* +1 ) = (I-HF)h(zk+1) = 0(\\h(zk+1)\n + o(\\h(zk+1)\\2) n Chapter 3. Stabilizations of Discrete Dynamics Formulations 65 Algorithm 3.1 INPUT from the one-step integration: Solve H(zk+1)D(zk+i)z Compute OUTPUT zk+i = zk + (j>(tk,zk,hk) = h(zk+i) z* = _D(zfc+i)z zk+1 = zk+1 - z* The post-stabilization with F satisfying (3.58) is our recommended method. For the vector field (3.39) with an invariant manifold (3.40) we can choose the stabilization matrix F = D ( i f D ) - 1 where D(z) is a matrix such that the inverse of the matrix HD exists. With this choice we propose A l g o r i t h m 3 . 1 . From Theorem 3.5 if the onestep integration scheme being used has 0(hp+1) Algorithm 3.1 will be 0(h2(p+^). local truncation error, the drifts with On the other hand, if I - HF ^ 0 then our theory requires the decreasing property, even though numerical results in [7, 8] indicate that such stabilizations may work anyway, although less effectively. Chapter 3. Stabilizations 3.4 of Discrete Dynamics 66 Formulations Stabilization for M u l t i b o d y D y n a m i c s w i t h H o l o n o m i c Constraints The Lagrange formulation of the equations describing a constrained multibody system can be written as q' = v M(q)v' = f(q,v)-<3(q)rA 0 = g(q) (3.63) where M ( q ) is the mass (n x n) matrix, f (q, v) is the vector of applied forces and A is the vector of dimension m of Lagrange multipliers. By differentiating the constraints twice a state space form differential equations for q and v can be obtained: q' = v M(q)v' = f(q,v)-G(q)TA (3.64) where A = ( G M ^ G ^ - ^ G M ^ f + Lv) Then the vector field (3.64) has the invariant manifold w ^ ( g(q) ) h(q,v) = The Jacobian matrix of h is given by / C(q) 0 # ( q , v ) = h(q,v)(q,v) = \L(q,v) G(q) where H has a full row rank. We first choose our stabilization matrix F to be F = HT{HHT)~X (3.65) Chapter 3. Stabilizations of Discrete Dynamics 67 Formulations Algorithm 3.2 INPUT *w= (q2>D one- step integration: (q^+i, 1 Solve T + <f>(tk,<lk,vk, M M r ( q f c + 1 , v m ) ( q T , v T f = (g T (q j t + 1 ),(G(q f c 4- i ) v * + i ) T ) T ( q T , v T = ^ ( q H l,Vfc+i )(q T ,v T ) r Compute OUTPUT ,vf + 1 ) T -(q r (q*+i> v fc+i) T = (q*+i vrf We consider Algorithm 3.2 If the local truncation error in the integration is 0(hp+1) the drift with Algorithm 3.2 is 0(h2(p+1>>). We call this algorithm S-full [8, 7]. It involves the inversion of matrix HHT which could be expensive for a large constrained multibody system. Also the evaluation of the entries of i J ( q , v ) , in particular Z(q, v ) , could be costly. In order to avoid the evaluation of L and solving with the big matrix HHT we proposed in [7] the stabilization matrix M-XGT (GM-1GT) 0 M~lGT 0 \ (3.66) (GM-lGT) The stabilization (3.66) was called S-both [8]. With the choice of (3.66) we have the amplifying matrix / — HF = G(q) L(q,v) \ I M-XGT 0 G(q) ) \ 0 0 -LM^GT (GM~lGT) 0 M-^G? -l (GM-lGT) 0 \ , [GM~1GT) -i (3.67) 0 / Chapter 3. Stabilizations of Discrete Dynamics 68 Formulations The amplifying matrix of S-both is not a constant matrix, so we cannot use its eigenvalues to guarantee the stability of the corresponding linear system (3.36). The norm condition of \\I — HF\\ < p < 1 does not hold for (3.67) either. However, a simple variation on the proof of Theorem 3.4 shows that its results still hold, essentially because / V n 0 0n \ / Ek 0 / \ £*-i 0 0 = 0 o regardless of what the m x m matrices Ek-i and Ek are. We now come to our favorite post-stabilization for constrained multibody systems. Suppose that from an order p one-step scheme we have (qk+i,vk+i). Applying (3.66) twice we obtain the following stabilization: / <ik+i \ ( qjfe+i - F{qk+1, \ V/t+i / \ Vfc+i / qk+1 \ I qfc+1 v fe+1 )h(qfc + i,vjfe + i) (3.68a) - F(qjfe + i,v fc+1 )h(qfc+i, v f c + i) (3.68b) where F is defined (3.66) and g(q) \ G(q)vJ In [7] we called the stabilization (3.68) S-both 2 . In S-both from (3.68b) we have h(q fc+ i,VA:+i) = (/-/r(qfc + i,Vjfc + i)F(qjfc + i,Vjfc+i))h(qjb + i,Vjt + i) + 0 (||h(qfc+i, v f c + i)|| 2 ) h(q f c + i,v f c + 1 ) = (I - H(qk+1,vk+1)F(qk+l,vk+1))h(qk+1,vk+1) + 0 (||h(qjfc+i, Vfc+i)||2) Using (3.68a) we have # ( q * + i , v * + i ) = /r(q J k + 1 ,v J t + 1 ) + 0 (||h(q f c + i,v f c + i)|| 2 ) Therefore h(q A+ i,Vjfc + i) = (I - H(qk+1,vk+1)F(qk+1,vk+1))h(qk+1,vk+1) = (I - H(qk+l,vk+1)F(qk+1,vk+1))2h(qk+1,vk+1) + O (||h(qjfc+i, v*+i)|| 2 ) + 0 (||h(q f c + 1 , v f c + 1 )|| 2 ) Chapter 3. Stabilizations of Discrete Dynamics Formulations 69 We can easily see that the amplifying matrix for the stabilization (3.68) is the zero matrix, i.e. ( / — HF)2 = 0. Applying Theorem 3.5 to the S-both 2 stabilization (3.68) we have the following corollary. Corollary 3.2 Suppose (q^+i, Vjt+i) is the solution of a numerical integration of (3.64) at tk+i with local truncation 0(hp+1) local truncation error 0(hp+1). error and the drift is Then the S-both2 stabilization (3.68) has 0(h2^p+1^). In [7] we have compared experimentally the post-stabilization algorithms S-full, Sboth, S-both 2 , and stabilizations of the velocity or position constraints alone [3, 12]. Rather than repeating those here we simply note that S-both 2 has clearly emerged as the best compromise between accuracy and economy. Compared to S-full of Algorithm 3.2, S-both 2 has the same order of drift term but computationally it is much cheaper than S-full. We therefore adopted S-both 2 as the method of choice, and will refer to it henceforth as PStab (for post-stabilization). Based on PStab we have developed a numerical code to integrate the Euler-Lagrange equations in constrained multibody systems. Our numerical experiments show (see next chapter) that PStab is efficient and robust for different multibody dynamics problems. E x a m p l e 3.3 A double planar pendulum with a prescribed path at its "free" end is described and tested in [4, 7]. We employ Algorithms 3.2 and 3.3 to simulate this mechanism. The onestep integration we used is an explicit Runge-Kutta scheme of order 2 with a constant step size h. In Table 3.3 we record the numerical drifts on position and velocity levels (drift-position and drift-velocity) of our runs on the time interval [0, 10]. Baum(12, 70) refers to Baumgarte's stabilization with Baumgarte parameter pair (12, 70). For this example we tried many other pairs of Baumgarte parameters and the pair (12, 70) is the "optimal" one. From Table 3.1 Algorithms 3.2 and 3.3 yield very small numerical Chapter 3. Stabilizations drifts which are 0(h2(p+1>). of Discrete Dynamics 70 Formulations We will present the equations of motion and the algebraic constraints in Section 4 of the next chapter as we will further consider the simulations of this mechanism (case I) using a numerical code. A l g o r i t h m 3.3 P S t a b 1 INPUT 2. Evaluate one-step integration M = M(q Jfc+1 ), G--= 3. Solve Af- 1GT(GM~1GT)-1q (qf+1,vf+1f = (qLvDT + G(qk+1), NSTAB 5. IF NSTAB < 1 THEN qk+1 = Qfc+1, v ^ + i = v fc+1 = NSTAB == C ( q M-l)Vfc+i + 1 GOTO 3; q*:+i and v fc+1 OUTPUT /s .01 .01 .01 .001 .001 .001 qfc, Vfc, hk) =0 = g(q* + 1 ), M-1GT(GM-1GT)-1v 4. Compute Qifc+l = <J*+1 — Cb vk+ 1 = Vfc+i — v, NSTAB ELSE ^ stabilization Baum(12, 70) Algorithm 3.2 Algorithm 3.3 Baum(12, 70) Algorithm 3.2 Algorithm 3.3 drift-velocity .72e-2 .12e-7 .67e-8 .70e-4 .15e-13 .18e-13 drift-position .14e-2 .61e-9 .15e-13 .14e-4 .40e-14 .31e-14 Table 3.1: Maximum drifts of a two-arm robot R e m a r k 3.4 There are other choices of the stabilization (3.66). When a multibody system matrix F than the one in involves heavy massed substructures and extremely Chapter 3. Stabilizations light mass components than GM~1G of Discrete Dynamics Formulations 71 we can choose the matrix GGT which can be better . In this case we take the stabilization GT(GGT)~1 conditioned matrix o N (3.69) GT(GGTy1 0 For the stabilization 3.5 matrix F (3.69) the conclusion of Theorem 3.5 holds. S t a b i l i z a t i o n for N o n h o l o n o m i c S y s t e m s In constrained multibody systems the constraints impose restrictions on the possible geometric positions of the individual bodies of the system. The so-called geometric constraints have the following form in general, g(*,q) = 0 (3.70) There can also be constraints of another type which restrict the kinematically possible motions of the system, i.e. the possible values of the velocities of the individual bodies, g(<,q,q) = 0 (3.71) The constraint (3.71) is called a kinematic constraint. Clearly every geometric constraint gives rise to a certain kinematic constraint. However a constraint on velocities in a system need not lead to a restriction on their possible positions. If the constraint (3.71) is integrable then (3.71) is essentially a geometric constraint. Nonintegrable kinematic constraints are not equivalent to geometric constraints, and they are called nonholonomic constraints. A nonholonomic multibody system is defined as a system with some nonintegrable kinematic constraints. Let us consider the following nonholonomic multibody system [68]. q' = v (3.72a) Chapter 3. Stabilizations of Discrete Dynamics Formulations 72 M(q)v' = f(q,v)-G(q)TA (3.72b) 0 = G(q)v (3.72c) where G has a full row rank. (3.72c) is an index-2 DAE. Notice that for holonomic constraint system if we differentiate the position constraints we will have equations of a similar form as (3.72c). All our conclusions here hold for the holonomic constraint system with the invariant set M.2 defined by M 2 = { ( q , v ) | , G(q)v = 0} (3.73) As before, we can eliminate the Lagrange multiplier vector A by differentiating the constraint (3.72c) with respect to t. We get the following underlying ODE for (3.72c): q' = v (3.74a) M(q)v' = f(q,v)-G(q)TA(q,v) (3.74b) The vector field (3.74) has (3.73) as its invariant set. Suppose that by an order p one-step integration of (3.74) we have a solution (c{k+i,^k+i) at tk+i- We propose the following stabilization: qfc+i\ vk+1J /qfc+i\ / 0 \ \M-iGT(GM-iGT)-1G(qk+1)vk+1) ~ \vk+1) The stabilization (3.75) can be developed as the general form stabilizations: qfc+i \ / q^+i \ = Vfc+l / -F(qjb + i,v A + i)G(qA. + i)vA.+i (3.76) V v^ + i / For (3.75) the stabilization matrix F has the form '-( ° \DiGD)-1 where D is a matrix such that GD has a bounded inverse. For example we can choose D = M~XGT or D = GT depending on the multibody system. For both of these choices of Chapter 3. Stabilizations of Discrete Dynamics Formulations 73 D, the amplifying matrix I — HF is zero for (3.75) so the drift is 0{h2^+^). A numerical test for stabilization (3.75) will be given for a rolling disc problem in the following chapter. R e m a r k 3.5 Of course any holonomic system (3.64) first and then the post-stabilization velocity constraints is not always best. particularly 3.6 can be written in the form (3.72c) (3.75) can be applied. This is the essence of stabilizing only [3j. However, results in [7] indicate that this simple The example in Section 5 of the following approach chapter is one where poor results are obtained when ignoring the position-level constraints. S t a b i l i z a t i o n s and P r o j e c t i o n M e t h o d s For the numerical integration of the an ODE z' = f(z) (3.77) 0 = h(z) (3.78) with the invariant manifold A coordinate projection method [11, 12, 16, 35] can be written as zn+i = (j)f(tk,zn,hn) (3.79a) zn+1 = z n + i - D(zn+1)\ (3.79b) = h(zn+1) (3.79c) 0 One reason for our claim [8] that the post-stabilization (3.45) is particularly attractive is that it is very close to the coordinate projection method (3.79): it can be easily verified that one Newton iteration for solving (3.79b) and (3.79c) yields (3.45). The poststabilization (3.45) guarantees the asymptotic stability of M. in the difference equations even when M. is slightly unstable in the underlying vector field. Therefore the numerical solutions will stay near ]V[ for all time integration. Chapter 4 S i m u l a t i o n of Constrained Rigid M u l t i b o d y D y n a m i c s The efficient simulation of multibody dynamics requires a concerted integration of several computational aspects [47, 81, 82]. These include selection of a data structure for the system's configuration, computerized generation of governing equations of motion, incorporation of constraint conditions and implementation of suitable solution algorithms. For the equations of motion in multibody systems, the methods used can, in general, be divided into two main approaches. In the first approach, the configuration of the system is identified using a set of Cartesian coordinates that describe the locations and orientations of the bodies in the system. This approach has the advantage that the dynamics formulation of the equations that govern the motion of the system is straightforward. Moreover, this approach, in general, allows easy additions of complex force functions and constraint equations. For each spatial rigid body in the system, six coordinates are sufficient to uniquely describe the configuration. Connectivity between different bodies can be introduced to the formulation using a set of nonlinear algebraic constraint equations. This set of constraint equations can be adjoined to the system equations of motion using Lagrange multipliers. In the second approach [47, 53, 65, 19, 52], relative or joint coordinates are used to formulate a minimum number of dynamic equations which are written in terms of the system degree of freedom. This approach leads to a complex recursive formulation based on loop closure equations. Unlike the formulation based on the Cartesian coordinates, 74 Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 75 the incorporation of general forcing functions and constraint equations or specified trajectories in the recursive formulation is difficult [19]. When there are closed-loops in a multibody system a recursive formulation has to be combined with a cut-joint technique in order to retrieve the closure conditions in the multibody systems. The equations of motion for a multibody system with closed-loops are index-3 or index-2 DAEs depending on the constraints. For a large multibody system suitable linear algebra algorithms can improve the efficiency of simulations of the multibody dynamics [65, 5]. Those linear algorithms take advantage of the sparseness structure in the equations of motion when Cartesian coordinates are used. Certain results of using the special linear algorithms on large multibody systems are reported in [65]. In this work we focus on stabilization methods and our numerical experiments are conducted on multibody systems with relatively small numbers of bodies ( < 7). For those relatively small multibody systems the special linear solvers have no immediate preference over the standard linear solvers like Gaussian elimination with partial pivoting. The efficiency of our post-stabilizations are demonstrated with simple linear solvers for those relatively small multibody systems. For our poststabilization PStab we rewrite it in a new form in this chapter and with the new form, PStab can be implemented by employing the efficient linear algebra algorithms [65, 5] designed for large systems. Combined with the efficient linear algebra algorithms, PStab can be expected to be efficient for large multibody systems as well. In this chapter we study numerical simulations for constrained multibody dynamics as an application of our post-stabilization method developed in the previous chapter. In Section 1 we consider problem formulation for multibody dynamic equations. For a given multibody system different choices of generalized coordinates result in different equations of motion. These differences between range from different dimensions, different Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 76 (algebraic) structure in mass matrix and Jacobian matrix and different index formulations. Of course for different formulations different numerical methods should be used. In Section 2 we consider an embedded Runge-Kutta scheme DOPRI(5,4) [43, 23] due to Dormand and Prince as our basic discretization. We perform PStab on the 5th order solution and use it as our next (starting) value of integration. The local error is obtained as the difference between the output (stabilized) 5th order solution and a modified (from DOPRI(5,4)) solution which is obtained by the 7th stage of function evaluation on the stabilized 5th solution. We prove that the modified solution also has the order 4 like the one in DOPRI(5,4). Furthermore this modified numerical scheme has the "first same as last" property and is a 6-stage method. Using this modified numerical scheme we discretize the underlying vector field for constrained multibody systems. For a multibody system with a symmetric positive semi-definite mass matrix, we rewrite PStab in a general form. In the new form we can implement PStab even when the mass matrix is singular. More importantly we can apply the special linear algebra algorithms [65, 5] designed for large multibody systems to our PStab and therefore make PStab efficient for large systems. In Section 3 we briefly describe our numerical simulation code, RKSTAB, for constrained multibody dynamics. RKSTAB also includes Baumgarte's stabilization as an option with a very small addition to the code. Therefore we can make a systematic comparison of PStab with Baumgarte's method on many multibody systems. In the comparison "optimal" Baumgarte parameters (found by ad hoc search) are used. A variety of numerical experiments are conducted with RKSTAB in section 4. Comparisons of RKSTAB with other codes like MEXX and DASSL on certain problems are also made. In Section 5 an ill-posted problem i.e. a singular slider-crank mechanism, is considered. Although RKSTAB is not designed for such a singular problem satisfactory simulation results are obtained with RKSTAB in this case as well. Finally in Section 6 a nonholonomic constrained mechanism is considered. Numerically nonholonomic constraint systems require Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 77 less sophistication than holonomic constraint systems as the nonholonomic systems have lower index. 4.1 4.1.1 P r o b l e m formulation Absolute coordinates The configuration of a multibody system in space can be specified by using absolute coordinates. For each body in a multibody system the absolute coordinates consist of three Cartesian coordinates Rk = (Xk, Yk, Zk)T indicating the position of center of mass of the body with respect to the global frame and the orientation of the body fixed frame with respect to the global frame. The orientation of a body is usually described by three angles (Eulerian angles) or by four Euler parameters [69, 85, 53, 65]. If Euler parameters Ok = (eg, e i i e2> e 3 ) T are use d a simple relationship exists between the components of the global angular velocity vector Qk and time derivatives of Euler parameters 0'k = (e»,e?,e»,e»)T [69], 2El^k 0'k = where Ek is a semi-transformation matrix [69] that depends linearly on Euler parameters: &k ~k ~k ~k Eu \ e 3 e 2 e e l 0 / The position variables for the multibody system are q = ( q ? V - - , q L ) T 5 qk = (RTk,oTky The velocity variables are v = (vf,---,vyT, vk=(RkT,nTky Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 78 The (body) position and velocity variables are related by ^ = , Rk\ , I Inp , for k = 1, • • • , nB. The descriptor form of equations of motion can be written as q' = MY' (4.1a) T(t,q)v = f ( i , q , v ) --GT(t,q)X 0 = (4.1b) (4.1c) g(*,q) where q has dimension up, v has dimension nv, matrix G has the form G = dg/dq-T(t, and T ( i , q ) = diag(Tk) q) is a block diagonal matrix. In (4.1) the mass matrix M is a constant block diagonal matrix [65] / M = diag(Mk), Mk \ mkI3 0 ^ 0 Jk ) where Jk is a 3 x 3 moment of inertia matrix. The mass matrix M in (4.1) is symmetric positive definite. For given consistent initial values qo and Vo one can solve the reduced index (index-1) form of (4.1) [5] q' M GfcqfWVN G(t,q) 0 )\\J = = T(t,q)v (4.2a) tt(t,^)\ \7(t,q,v)J where 7 ( t , q , v ) can be derived from twice differentiations of g ( t , q ) . The constraints g = 0 indicate the kinematic relations among the adjacent bodies in the system and the normalization condition | | ^ | | = 1 for Euler parameters. Therefore matrix G is sparse and has a block tridiagonal structure [5, 65]. Chapter 4. Simulation of Constrained Rigid Multibody Dynamics R e m a r k 4.1 Euler angles can also be used as orientation coordinates np = nv. In this case the Tk is a nonlinear square matrix function Euler parameters 79 which leads to of the Euler angles. result in a simple form of matrix Tk but have four variables and an additional algebraic condition i.e. the normalization condition is required. The equations of motion for a multibody system are obtained in a similar way by using Euler angles or Euler parameters. For notational simplicity we suppress time t in the equations of motion and assume that the matrix G is the Jacobian matrix of g(q) (i.e. we assume the matrix T in (4.2a) to be the identity matrix) as before. Our discussion on a general form can be made through minor modifications. (4.2) is a system of index-1 DAEs for variables (q, v, A). The values of A can be determined once the numerical solutions of (q, v) are computed. Numerical integrations of (q, v) in (4.2) frequently require solving linear systems with the coefficient matrix: / M G(t,q)T\ (4.3) \G(t,q) 0 ) Therefore the efficiency of numerical (explicit) integration of (4.2) depends on the performance of linear solvers, especially for large multibody systems. Direct and indirect methods for solving such a linear system have been reported in [65, 5, 16]. For large multibody systems with complex kinematic structures more work is needed for efficient and fast linear algorithms. 4.1.2 R e l a t i v e coordinates For a multibody system the degree of freedom between each pair of adjacent bodies is often low and in many cases is one, as in a chain of links of a manipulator. Therefore using absolute coordinates can result in a relatively large amount of redundancy among the generalized coordinates. Relative coordinates are used to represent the joint degree of Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 80 freedom between adjacent links [46, 47, 57, 65, 5, 81, 82]. In the case of tree-configured multibody systems, using relative coordinates leads to a formulation of equations of motion as an ODE system M r ( q ) q " = f r (q,q') (4.4) The dimension of (4.4) is equal to the multibody system's degree of freedom and there are no algebraic constraints. In control applications, it is desirable to determine joint reacting force (torques) associated with the joint relative coordinates. The relative coordinates provide the measurable variables, while the reaction forces (torques) associated with the relative joint degrees of freedom represent the actuating forces required to drive the system components [19, 20]. For a forward dynamic analysis the mass matrix M r ( q ) is dense and evaluating the force term fr can be complicated. Furthermore, formulation (4.4) can be cumbersome when dealing with a multibody system with variable kinematic constraints. When there are closed-loops in a multibody system the loop closure conditions will eventually lead the equations of motion of the multibody system in general to a DAE form rather than (4.4). Absolute and relative coordinates have their advantages and disadvantages. There are certain discussions of attempting to use the two in a mixed form of equations of motion [53, 65]. Then depending on the purpose of the dynamics analysis (inverse or forward dynamics analysis) relative or absolute coordinates are explicitly computed and the others remain intermediate values. In [53] the relationship among various recursive formulations for a multibody system with an open loop is considered in a unified way. For a multibody system with a complicated kinematic structure a lot of work may be needed to find an efficient formulation. Chapter 4. Simulation 4.2 of Constrained Rigid Multibody Dynamics 81 Implementation with Runge-Kutta Methods Consider solving an initial value ODE problem z' = f(t,z) (4.5a) z(*0) = zo (4.5b) along with an invariant manifold M. given by 0 = h(*,z) (4.6) A numerical method for (4.5) can be classified as a multistep, Runge-Kutta or extrapolation method according to the type of its underlying basic formula. A scheme for (4.5) can also be classified as a stiff or nonstiff method based on whether or not the scheme is suitable for solving stiff problems [43, 59, 27]. In our study here, however, we consider the numerical solution of (4.5) on the manifold (4.6). Because of the presence of algebraic equations (4.6) numerical integrations or behavior of existing codes on (4.5) could be different. For example, as we pointed out in previous chapters, the stability property of invariant manifold defined by (4.6) in the vector field (4.5a) can directly affect the efficiency of numerical integration of (4.5). Therefore we have to consider the vector field (4.5a) and the manifold (4.6) in a uniform way when integrating (4.5) numerically. The post-stabilization developed in the previous chapter provides a robust and efficient algorithm for such a purpose. Particularly when the ODE (4.5a) is non-stiff we can employ an explicit scheme and along with our poststabilization to integrate (4.5) on the manifold (4.6). We choose Runge-Kutta schemes as the discretization because the overhead is small with these methods. Furthermore, Runge-Kutta methods are known to be preferable when there are discontinuities in the equations [43]. The post-stabilization technique can be used with other discretization methods as well, though. Chapter 4. Simulation of Constrained Rigid Multibody 82 Dynamics Among initial value ODE solvers virtually all effective methods are adaptive, variablestep size methods [27]. In such methods a user is asked to specify accuracy parameters (both absolute and relative error tolerances). A numerical method then includes a stepchoosing strategy which is usually based on monitoring and controlling the size of the local error introduced in each step of the integration [27]. The most efficient error control strategies for explicit Runge-Kutta methods involve constructing Runge-Kutta formulas which contain two numerical approximations z n + 1 and z n + i which share the same function evaluations: s zn+1 = Zn + h^kki (4.7a) t=i s zn+i = zn + hJ2bik (4-7b) where ^1 = t(tn?zn) ki — f (t0 + Cih, zn + h ^2 i-l aijkj) i=i for i = 2, • • •, 5. z ra+ i and z n + 1 are approximations of z(tn + h) of order p and q (q > p, usually q = p + 1) [59, 43]. The difference s zn+1 - z n + i = hJ2 Eik then yields an estimate of the local error with ET = bT — bT = ( £ i , • • •, ES)T. (4.8) This is called an embedding Runge-Kutta formula [43, 59]. If the solution z ra+ i is used as the initial value for the next step, the method has order p and the method is called (p,q) or p(q) method. When z n +i is used as the next step initial value the method is called (q,p) or q(p) method. Obviously the (q,p) method has higher order than the (p,q) method but it has a less reliable local error estimation (4.8) than the (p,q) method [43, 59]. Chapter 4. Simulation of Constrained Rigid Multibody 83 Dynamics One of the most popular embedding Runge-Kutta methods is DOPRI(5,4) [43, 59, 23] due to Dormand and Prince. In DOPRI(5,4) it is the higher order solution that has the small moduli of the coefficients in the principal part of local truncation error [59]. The modified Butcher array for DOPRI(5,4) has the following form [59, 43]: 0 z 1 5 1 5 3 10 3 40 4 5 44 45 56 15 32 9 8 9 19372 6561 25360 2187 64448 6561 1 9017 3168 355 33 46732 5247 49 176 5103 18656 1 35 348 0 500 1113 125 192 2187 6784 11 84 Zra+1 5179 57600 0 7571 16695 293 640 92097 339200 187 2100 1 40 Zra+1 35 348 0 500 1113 125 192 2187 6784 11 84 0 71 57600 0 71 16695 71 1290 17253 339200 22 525 n+l _ Zn+1 9 40 212 729 (4.9) 1 40 For our purpose of implementing the post-stabilization PStab we have to modify DOPRI(5,4) (4.9). We construct our post-stabilization scheme as follows. zn+1 = zn + hj^bjkj 3=1 zn+1 = zn+1 - F(zn+i)h(zn+i) z n +i = Zn+i - F(zn+i)h(zn+i) (4.10) Unlike scheme (4.9) we take kj = f(tn + h,zn+1) (4.11) h,zn+1) (4.12) while in DOPRI(5,4) (4.9) k7 has the form &7 = f(tn + Chapter 4. Simulation of Constrained Rigid Multibody Dynamics For the error estimation we construct the following numerical solution for (4.5a) at tn+\ 6 z n + 1 = zn + h £ b3kn3 + hb7fy (4.13) We claim that z n + 1 of (4.13) is a 4th order solution of (4.5a) if the vector function f is smooth. To see this, let 7 ( 4 - 14 ) zn+1 = zn + h^W? 3=1 Then z„+i is the 4th order solution in DOPRI(5,4). Using theorem 3.5 we have ||z n+1 -z n+1 || = \\hfy(k?-k7)\\ = —\\f(tn + h,zn+1)-f(tn - ^Mo||Zn+l -Zn+l|| < ^M||h(z n + 1 )|| < A^. +1 + h,zn+1)\\ = 0 (,r) where M0, M and M0 are constants depending or f and h. Therefore z n + i (4.13) is a solution of (4.5a) with the same order as z n +i. As (4.11) is used it is easy to see that t,n _ Ixi'j Ln+l A/-I So the modified method has the "First Same As Lasf (FSAL) property and it is actually a 6-stage method. In general M ( q ) is symmetric and positive semi-definite 1 and matrix G(q) has a full row rank. We therefore rewrite PStab for constrained multibody systems as follows. 1 When a mixture of absolute and relative coordinates is used the mass matrix has the form M - ' Moo 0 0 0 where MQO is a symmetric positive definite matrix. Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 85 Starting from the 5th order solutions q n +i and v n +i, evaluate matrices M and G at z n + i and Vn+i we construct PStab in the following way M GT G 0 J I Vn+1 = Vn+i V n +1 = - J = J q«+i GT\ ' M v Qn+i I G ( M (4.15a) PG(qn+1)vn+1 (4.15b) = v n + i - /? n+ i 0 J GT G 0 M GT ^ V n + i -J G Pg(qrW-l) = q re +l - Cin+l -Pg(c[n+l) = 4i+l - «n+l (4.15c) P G ( q n + i ) v n + i = v n + 1 - J3n+i (4.15d) 0 where J and P are constant projection matrices given by J — \^npt "npXnX) R e m a r k 4.2 * = \yn\xnpi J-nX) 1. The new form of PStab (4-15) allows us to implement PStab when the mass matrix M is singular. When M is nonsingular, formulation (4-15) coincides with A l g o r i t h m 3.3 proposed in Chapter 3. 2. Through formulation (4-15) PStah has the same linear equations as the ones we have to solve at each time step. Therefore PStab can he implemented efficiently large systems once appropriate linear equation solvers are developed. By those linear solvers the cost of performing for employing PStah can he controlled to be a small amount for large systems, in the sense that it is always a small fraction of the cost of carrying out the unstabilized time step. For PStab (4.15) we have the following theorem. T h e o r e m 4.1 The manifold M. is uniformly stable in the difference equations Furthermore, PStab (4-15) has local error 0(h6) and numerical drifts 0(h12). (4-15). Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 86 In order to prove Theorem 4.1 we first prove the following lemma. L e m m a 4.1 Let G = G(q) be the Jacobian matrix o/g(q) and M ( q ) the mass matrix for a constrained multibody system. positive semi-definite Suppose G has a full row rank, and M(q) symmetric for all q. If for all a £ Rnp with G(q)a = 0 and a ^ 0 aTM(q)a > 0 Then the is (4.16) matrix M GT N (4.17) G is nonsingular 0 ) and for any vector w £ RnX we have w - M(q) G(q)J GT(q)X Pw = 0 G(q) o (4.18) y P r o o f of L e m m a 4.1 The nonsingularity of matrix (4.17) can be easily seen as follows. Suppose u = ( a T , b r ) T £ RnP+nX is a vector with Then Multiplying (4.19) by a Ma+GTb = 0 (4.19) Ga = 0 (4.20) and noticing that a G arMa = 0 = 0 we have Chapter 4. Simulation of Constrained Rigid Multibody 87 Dynamics Therefore a = 0 by (4.16) and (4.20). Replacing a with 0 in (4.19) and noticing that G has a full row rank we have b = 0. Therefore matrix (4.17) is nonsingular. To prove (4.18) let x and y be the vectors such that M GT \ / x \ /0 Pw = G 0 I w Then Mx + Gry = 0 Gx = w Therefore (4.18) follows from the fact that G J I X J = Gx = w • R e m a r k 4 . 3 Using lemma 1^.1 we have g(<Wi) - G(qn+1)an+i M GT \ g(q„+i) - GJ ^g(qn+i) G 0 and similarly we have G ( q n + i ) v n + i - G(q n + i)/? n + i = 0 P r o o f of T h e o r e m 4.1: From the PStab formulation (4.15) and Lemma 4.1 we have g(q„+i) = g(qn+i-a„+i) = g ( q n + i ) - G ( q n + i ) a n + i + 0 (||an+i| Chapter 4. Simulation = G(qn+i)vn+i of Constrained Rigid Multibody Dynamics 88 0(\\an+1\f) = G(qn+1 - an+1)(vn+1 - f3n+1) = (C(q„+l) - gqq(q„+l, Ctn+i)) (v n +i - $n+1) + O (\\&n+1 || 2 ) = G(qn+i)vn+i = G(qn+i)v n + 1 - C?(q n+1 )/? n+1 - g q q ( q „ + i , v „ + 1 ) o r i + i + O ( | | a n + i , / ? n + 1 | | 2 = -gqq(q«+l 5 V n + i ) a n + i + O (||o: n + l>/Wll 2 ) - g q q ( q r i + 1 , a n + 1 , v r i + i ) - G(qn+i)J3n+i + O (||a n + i,/3 n +i|j Let B be the matrix B = ( M(qn+1) GT(qn+l)\ V G(qn+1) Using (4.15a) we have an+\ = Bg(qn+i) 0 1 p ) and also an+1 = 0 ( g ( q n + 1 ) ) = 0(h6) Pn+1 = 0 ( G ( q n + i ) v n + 1 ) = 0(h6) Therefore / g(q»+i) \ / 0nA 0nA W g(qn+i) \ )+0 \<?(q»+i)Vn+i/ \ -gqq(qn+i,vn+1)JB Similarly we can write g(q n +i) and G(qn+1)vn+1 0nA / \ G(qn+1)vn+1 /12N , (h12) (4.21) / in terms of g(q„+i) and G(qn+1)vn+1 like (4.21). From (4.15) those an+i and /3n+1 have the same order as those of g(q n +i) and G(q n +i)v„+i, therefore the condition (3.58) in Chapter 3 is satisfied. Applying Theorem 3.5 we have that the drift is 0(h12). Finally the uniform stability can be obtained from Theorem 3.3 in Chapter 3. 4.3 • A Numerical Code: R K S T A B R K S T A B is a Fortran code for time integration of constrained (both holonomic and nonholonomic) mechanical systems. It is suitable for integration of the index-1 form of Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 89 the equations of motion in descriptor form /M(i,q) GT(t,q)\ U(i,q) 0 /q"\ _ ) \ \ ) ~ /f(t,q,vr (4.22) \-L(t,q)vt with the constraints defined by g(t,q) = 0 (4.23a) G(t,q)v = 0 (4.23b) RKSTAB uses DOPRI(5,4) as the discretization scheme and PStab as the stabilization method. The local error control strategy is adapted from the ODE solver DOPRI5 developed by Hairer et al [43]. Several subroutines from LINPACK axe used as well. RKSTAB consists of a set of subroutines and the user can interface the code by calling RKSTAB from a user written driver program. The user has to supply the problem describing functions and matrices M(t,q), f ( t , q , v ) , G ( t , q ) , g ( i , q ) (for holonomic constraints only), gt(t, q) and L(t, q, v ) . The subroutines which describe the functions and matrices above have the following form: C*** MASS MATRIX M(t,q) : C SUBROUTINE AMAT(N,NDIM,T,Q,AM) C DOUBLE PRECISION AM(NDIM.N),Q(N) C*** VALUE OF GENERALIZED FORCE f(t,q,v) : C SUBROUTINE FCN(N,T,Q,V,F) C DOUBLE PRECISION T,Q(N),V(N),F(N) C*** VALUE OF THE CONSTRAINT FUNCTION C SUBROUTINE GCN(M,N,T,q,G) C DOUBLE PRECISION T,Q(N),G(M) C*** JACOBIAN MATRIX OF g g(t,q) : WITH RESPECT TO q : Chapter 4. Simulation of Constrained Rigid Multibody Dynamics C SUBROUTINE GCNQ(M,N,T,Q,GQ) C DOUBLE PRECISION T,Q(N),GQ(M,N) C*** DERIVATIVE OF g WITH RESPECT TO C SUBROUTINE GCNT(M,N,T,Q,GT) C DOUBLE PRECISION T,Q(N),GT(M) C*** SECOND DIRECTIONAL DERIVATIVE OF 90 t . g : C (d"2/dt~2)g(t,q) = (d~2/ds~2) g(t+s,q+s*v) I s=0 C (d~2/dt~2)g(q,t) = GQ*v' + L(t,q,v) C SUBROUTINE GCNqQ(M,N,T,Q,V,L) C DOUBLE PRECISION T,Q(N),V(N),L(M) Solutions (of the position and velocity), along with the constraint forces (or the La- grange multiplier vector) at the internally selected integration points or user specified points can be printed out as output by a user supplied subroutine. The user can also monitor and let a simulation stop once certain physically meaningful violations occur by supplying additional physical constraints in the driver. The main routine RKSTAB is called by the user supplied main program as CALL RKSTAB(IHOL,N,M,MI,AMAT,FCN,GCN,GCNQ,GCNT,GCNQQ,GIN, & BAUMG,X,Q,V,XEND,H, & RT0L,AT0L,DAT0L,IT0L, & S0L0UT,I0UT, & WORK,LWORK,IWORK,LIWORK,LRCONT,LICONT,IDID) While the purpose of RKSTAB is to implement PStab, a very small addition to the code allows us to also include Baumgarte's method as an option. The task of choosing the Baumgarte parameters is left to the user. The user also needs to specify the following parameters in the main program. Chapter 4. Simulation of Constrained Rigid Multibody Dynamics IHOL the code assumes holonomic constraints unless ihol=l on input N dimension of the generalized coordinates M number of constraints BAUMG Baumgarte parameter, relevant only if Baumgarte's method is chosen T initial time Q(N) initial position value V(N) initial velocity value RTOL,ATOL relative and absolute tolerances DATOL drift absolute tolerance (Specifying DATOL=0 implies that 91 PStab is used at every step.) 4.4 N u m e r i c a l simulations of h o l o n o m i c constrained s y s t e m s In this section several multibody systems are integrated with RKSTAB. We start with Andrews squeezing mechanism which has been regarded as a benchmark for testing numerical methods for multibody systems [44, 65, 5]. A preliminary comparison of RKSTAB with another numerical code MEXX with respect to CPU time has been made for various error tolerances. Compared to Baumgarte's stabilization with "optimal" Baumgarte parameters, PStab uses fewer function evaluations and the amount of CPU time spent on the post-stabilization PStab is about 1% of the total CPU time in RKSTAB. We then consider a five-point vehicle suspension mechanism [40, 63]. This mechanism is a spatial 7-body system with 4 closed loops. Numerical experiments with this mechanism are conducted when the wheel has excitations from rough road conditions. Numerical solutions with Baumgarte's stabilization have relatively large drifts especially when the road excitation involves a discontinuity. We then turn to a two-link manipulator with prescribed path at the tip, considered in [7]. Improvement on Baumgarte's stabilization Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 92 with PStab is obtained with fewer function evulations and much smaller drifts in the numerical solutions. We also conduct simulations for a highly oscillating slider-crank mechanism considered in [100]. The rapid oscillation in the slider-crank is caused by a translational spring damping actuator. In [100] the numerical code DASSL which is designed for general DAEs and stiff ODEs, was used for several stabilized index-2 formulations. Our numerical experiments with RKSTAB show that significantly fewer function evaluations are used than those with DASSL reported in [100], suggesting that the underlying ODEs of the mechanism is not stiff even if there are rapid oscillations. The fifth example we considered is a singular slider-crank mechanism investigated in [74]. In this mechanism the Jacobian matrix G becomes rank-deficient at certain points. The optimization matrix becomes singular of those points and therefore the underlying ODEs give ill-posed problems. PStab is not designed for such a singular system. But the way that PStab enforces the numerical solutions back into the invariant manifold, does make a difference in the simulations. For a fairly small error tolerances RKSTAB can carry out the integration for a large time interval in this example. Finally we simulate a rolling disk on the horizontal plane. This is a nonholonomic constrained system with the constraints involving generalized velocity v. Unlike the holonomic case, this nonholonomic constrained system yields an index-2 DAE. Based on the post-stabilization proposed for nonholonomic constrained systems in Chapter 3 RKSTAB has an option to integrate such a system. Compared to holonomic constrained systems, nonholonomic constrained systems are easier to derive methods for its lower index nature. For parameter settings in those experiments, we take the absolute tolerance equal to the relative tolerance divided by 10 in all of our numerical experiments unless otherwise stated. Other parameters are defined as follows. N F C N is the number of function evaluations used Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 93 Figure 4.1: Andrews squeezing mechanism N S T E P is the number of steps needed to carry out the integration N A C C P is the number of steps that RKSTAB predicts and accepts N S T A B is the number of PStab steps used 4.4.1 Andrews squeezing mechanism The Andrews squeezing mechanism has been extensively used to test numerical methods or simulation codes for constrained multibody dynamics by many researchers [81, 43, 8, 16, 65, 47, 70]. Most of the simulation runs are made over a small time interval [0,0.015]. Its simulation on a larger time interval, however, becomes more difficult due to its highly oscillatory nature (see Fig. 4.2). We test our RKSTAB on the time interval [0,0.3]. A preliminary comparison of RKSTAB with another simulation code MEXX 2 We thank Professor Ch. Lubich for kindly sending us the code MEXX. 2 [65] is made Chapter 4. Simulation of Constrained Rigid Multibody -81 0 Dynamics 94 ' ' ' 0.1 0.2 0.3 Figure 4.2: Solution of Andrews squeezing mechanism as well. The Andrews squeezing mechanism consists of 7 rigid bodies denoted by K{ in Fig. 4.1 [43]. The rigid bodies are connected by joints without friction. The planar motion of this mechanism is caused by a motor located at O and a spring acting on body K3. The driving motor in the mechanism rotates at a high speed. rtol l.d-4 l.d-5 l.d-6 l.d-7 B(20,100) NFCN 11150 17132 25625 38828 PStab NFCN 10784 17024 25592 38780 B(20,100) NSTEP 1858 2855 4275 6471 PStab NSTEP 1791 2837 4265 6463 Table 4.1: Comparison of B(20,100) with PStab on squeezing mechanism Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 0.04 0.02 -0.02 -0.04 -0.06 Figure 4.3: Drifts with B(20,100) on squeezing mechanism 0 0.1 0.2 Figure 4.4: Drifts with PStab on squeezing mechanism 0.3 Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 96 Our first numerical experiment on this mechanism is to run a Baumgarte's stabilization and PStab for different tolerances. The Baumgarte parameters used are 71 = 20 and 7 2 = 100. As shown in Table 4.1, PStab always takes slightly fewer time steps and function evaluations than Baumgarte's formulation for all relative tolerances chosen. More importantly, Baumgarte's stabilization has larger drifts than PStab as shown in Fig. 4.3 and Fig. 4.4. datol l.d-5 l.d-6 l.d-7 l.d-8 l.d-9 0 NFCN 16427 17024 17066 17090 17090 17096 NSTEP 2475 2837 2844 2848 2848 2849 NACCP 2436 2464 2446 2446 2446 2446 NSTAB 5 601 2126 2821 2839 2849 CPU 14.71 15.48 16.23 16.34 16.43 16.36 STAB-CPU(%) 0.24 0.58 0.93 1.03 1.08 1.10 Table 4.2: Results of RKSTAB on Andrews squeezing mechanism The additional amount of work for implementing the stabilization PStab is small. To see this we then fix the relative tolerance = l.d-5 and try different constraint tolerances (absolute) datol for PStab. As shown in Table 4.2 the CPU time spent on the poststabilization, which is denoted by STAB-CPU(%), is about 1% of total CPU time for the whole simulation. There are several numerical codes available for simulation of constrained multibody systems [17, 16, 44, 46, 65]. One of the frequently used codes is MEXX. MEXX is based on extrapolation of a time stepping method that is explicit in the differential equations and linearly implicit in the nonlinear constraints [65]. A preliminary comparison has been made between MEXX 3 3 and RKSTAB on this 7 body system. For different relative The version of MEXX is faster than the old version for larger problems, because of different linear system solvers. But for this problem the linear algorithm used in MEXX is similar to ours and is faster than the O(N) alternative. Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 97 tolerances which range from l.d-8 to l.d-4 we compute the CPU and NFCN. As shown in Table 4.3 for tolerances > l.d-7 RKSTAB does a better job as it takes fewer function evaluations and CPU times than MEXX. For tolerance < l.d-8 MEXX is better in terms of CPUs and NFCNs. This is not surprising because for very small tolerances DOPRI(5,4) will not be as efficient as a higher order scheme [43]. For tolerances smaller than l.d-7 MEXX switches to higher order discretization and therefore it works well for those small tolerances. rtol l.d-4 l.d-5 l.d-6 l.d-7 l.d-8 MEXX CPU 18.61 25.72 33.62 40.54 46.08 RKSTAB CPU 10.08 14.41 22.53 33.99 52.18 MEXX NFCN 18724 26060 34550 41764 47935 RKSTAB NFCN 10712 16550 25688 39230 59852 Table 4.3: Comparison of MEXX and RKSTAB on squeezing mechanism 4.4.2 A vehicle s u s p e n s i o n m o d e l The simulation of mechanical systems is a field of current research in both applied mathematics and mechanical engineering [82, 47]. Particularly in robotics and vehicle dynamics, the multibody system approach provides an efficient tool for the design and analysis process [93]. The vehicle dynamics simulations can also be used for a ride quality investigation and a driving safety analysis by means of the variations of the contact forces between the tires and the road. D e s c r i p t i o n of t h e m o d e l The five-point vehicle suspension mechanism as shown in Fig. 4.5, is a spatial multibody system [40, 63]. The mechanical model consists of a wheel attached to a wheel Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 98 <t>l , V i pherical joints universal joints Figure 4.5: A five-point vehicle suspension mechanism carrier at point B by a revolute joint which is attached to a car body by five rods. The rods are connected to the car body by universal joints and to the wheel carrier by spherical joints. The suspension has four coupled closed loops as shown in Fig. 4.6. The existence of closed loops complicates the situation when we set up the equations of motion. Unlike an open-loop problem a recursive algorithm can not be used for a closed loop system directly. A well known technique to deal with closed loop problems is to open the closed loops by adding "cut-joint constraints" [47, 63]. In the suspension model, after cutting the spherical joints between the rods 2 to 5 and the wheel carrier as shown in Fig. 4.6, we have an open-loop system with fourteen degrees of freedom. The generalized coordinates used are as follows [63]. (<^i, ij>i) is the universal joint variable of rod i and car body, for i = 1, • • • 5; ( a , / ? , 7 ) is the spherical joint variable of rod 1 and the wheel carrier; Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 99 car body cut joint Figure 4.6: Suspension graph with cut-joints indicated 6 is the rotation of the wheel respect to the car body. The equations of motion for the unconstrained system can be obtained by using a recursive algorithm for the fourteen relative joint coordinates. Then D'Alembert's virtual work principle [47] is utilized to drive the equations of motion, augmented by reaction forces from the cut-joint constraints introduced via Lagrange multipliers. There are twelve constraint equations corresponding to the cuts on four spherical joints. The equations of motion for the suspension mechanism were generated by G. Leister with the help of the program NEWEUL 4 [63]. In simulations the car body is assumed to be fixed to the inertial frame and the wheel is assumed to rotate with a constant angular velocity (30m/second) [40, 63]. There is a spring damper element between rod (body) 5 and the car body which represents the shock 4 We thank Dr. B. Simeon for his kindly sending us the program which evaluates the equations of motion for the suspension model. Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 100 absorber of the suspension. Friction effects are not considered. The road is modeled as an excitation acting on the wheel. Both smooth excitation and discontinuous excitation are considered. The tire is modeled by a simple linear spring. The acting force is given [40] by F0 = -200000 -(R-U) + 4404134 [TV] (4.24) where R denotes the vertical displacement of the wheel center and U denotes the value of the excitation. Force F in (4.24) could become negative if R — U is too large, which means that the tire is not compressed in this position. The tire does not transfer a negative force which would hold the wheel on the road surface. In fact, if F is equal to zero, the wheel takes off from the road [40]. Therefore the actual force used is F = max{-200000 • (R - U) + 4404134,0}[iV] For the road condition we consider two different excitations. At first a Fourier series approximation (one mode for simplicity) of driveway measurements with an additional exponential damping factor is used as excitation. Secondly the road is modeled as a stair function with a jump at time t = 0.2s [40]. Case I: A s m o o t h e x c i t a t i o n Consider a smooth excitation given by U = 300sin[27r(* + 0.20)] • e x p - 5 0 0 ^ 0 ' 3 ) 2 We simulate the suspension model on the interval [0,1] using RKSTAB. Figs. 4.7 and 4.8 display the vertical position of the wheel and the selected solution components with relative tolerance = l.d-4. The wheels vertical position experiences the main disturbance in the interval [0.2,0.4] which corresponds to the presence of a relatively large excitation. Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 101 Figure 4.7: Vertical position of the vehicle wheel, case I 0.3 0.2 0.1 -0.1 i -0.2- -0.3 0 0.1 0.2 i 0.3 0.4 0.5 Time (s) 0.6 0.7 0.8 0.9 1 Figure 4.8: Solution (<^i, t/>i,a, ^ , 7 ) of the suspension, case I Chapter 4. Simulation of Constrained Rigid Multibody 102 Dynamics 0.04 0.02 -0.02 - -0.04 -0.06 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 4.9: Drifts with B(10,30) for the suspension model, case I The vibration is absorbed by the spring damper afterwards as shown in Figs. 4.7 and 4.8. We test RKSTAB with different relative tolerances and the results are summarized in Table 4.4. rtol l.d-3 l.d-4 l.d-5 l.d-6 l.d-7 NFCN 530 896 1298 2060 3284 NSTEP 88 149 216 343 547 NACCP 85 136 211 338 541 NREJCT 3 13 5 5 6 Table 4.4: Results of RKSTAB on the suspension mechanism, case I Notice that Baumgarte's stabilization does not work well for this model around the excitation period. As shown in Fig. 4.9 the drifts in velocities are unacceptably large right after the main excitation. In contrast, for PStab the drifts on position and velocity Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 103 x10 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 4.10: Drifts with PStab for the suspension model, case I levels are very small. The relative tolerance taken for those runs is l.d-4 and atol is l.d-7. Case II: A discontinuous e x c i t a t i o n Consider a discontinuous excitation which happens after driving over the border at time 0.2 seconds [63]. The wheel then will lose contact with the road. Due to the missing contact force the suspension drops and after some vibrations moves to a new equilibrium position. The excitation is given [63, 40] by 0 if t < 0.2 100 if t > 0.2 U The situation becomes complicated for a code to select or predict a step size when discontinuity occurs. We run RKSTAB on this discontinuous model for different tolerances and make a comparison with Baumgarte's formulation. As shown in Table 4.5 PStab always has smaller NFCN and N R E J C T than those of Baumgarte's stabilization. Chapter 4. Simulation of Constrained Rigid Multibody Dynamics -160 Figure 4.11: Vertical position of the vehicle wheel, case II 0.6 ~1 1 1 1 1 1 T" 1 0.4 Figure 4.12: Solution (<^!, ^1,01,^,7) of the suspension, case II Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 0.25 0.15- 0.05 -0.05 - 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 4.13: Drifts with B(40,400) for the suspension, case II x10'~ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 4.14: Drifts with PStab for the suspension, case II Chapter 4. Simulation rtol l.d-3 l.d-4 l.d-5 l.d-6 l.d-7 of Constrained Rigid Multibody B(40,400) NFCN 644 1040 1556 2240 3404 PStab NFCN 608 968 1442 2216 3356 B(40,400) NREJCT 19 37 41 41 47 Dynamics 106 PStab NREJCT 16 28 32 40 45 Table 4.5: Comparison of B(40,400) and PStab on suspension model, case II The drifts in Baumgarte's stabilization are complicated when the system has discontinuities. As shown in Fig. 4.13 the drifts on velocities are much larger than those of PStab shown in Fig. 4.14. Compared to case I we can see that Baumgarte's stabilization has more difficulty for discontinuous systems while PStab treats the two systems in an equal way. 4.4.3 A planar two-link m a n i p u l a t o r w i t h prescribed p a t h A two-link planar robotic system with a prescribed path at its end effector is investigated (see Fig. 4.15) in [4, 7]. The system consists of two rigid links: link 1 with its one end fixed at the origin and link 2. In link 2 one end is connected to link 1 with rotations allowed in the x-y plane and the other end (end effector) is set with a prescribed path. Two different types of prescribed path are considered. In case I, the y-coordinate of the end effector is given a time dependent motion, while the x-coordinate is determined by the forces acting on the system. Therefore the constraint can be written as y2 = g(x2) where we take g(x) = bxx2 + 62, &i and b2 are two constants. (4.25) Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 107 (*i > y) (*2>%) Figure 4.15: A two-link planar robotic system In case I I , a path for the end effector is specified by V2 = f(t) (4.26) where f(t) = sin2(u)t) and a; is a constant. To prescribe the equations of motion let 0^ be the angle that link 1 makes with the horizontal axis, and let 02 be the angle that link 2 makes with respect to link 1. The masses of the links are denoted by m,- and their lengths are denoted by U. The coordinates of the link 2 and its end effector are given by x\ = hci J/1 = /isia;2 = xx + l2cu 2/2 = 3/1 + I2S12 where C; = cos#;, s,- = sin#;, c\2 = cos(#i + #2)5 &vi = sin(#i + 92). The generalized coordinate q = (9i,02)T and / mi/2/3 + m2(/12 + q/3 + hl2c2) V m2(P2/3 + hl2c2/2) m2(P2/3 + hl2c2/2) m2P2/3 \ ) Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 108 -0.5 0 x-coordinate of the tip Figure 4.16: Constraint path of the end-effector, Case I -rmghd/2 - m2g(hc1 + / 2 c 1 2 /2) \ / m2l1l2s2/2(26162 + 0\) f = -m2gl2c12/2 J \ -m2hl2s29l/2 In the following simulations the following data are used. mi = m2 = 3Qkg, l\ = l2 — lm, g = 9.81m/.s 2 0i(0) = 70°, 02(O) = - 1 4 0 ° , ^i(O) = 62(0) = 0 Case I Consider the case where the coordinates (x2,y2) are constrained to lie on a parabola V2 = x22 - B where (3 is chosen such that the above initial conditions are consistent. For different tolerances comparisons are made for Baumgarte's stabilization and PStab. We take Tend = 10 and the Baumgarte parameters 71 = 12 and 72 = 70. Chapter 4. Simulation of Constrained Rigid Multibody Dynamics Figure 4.17: Joint constraint torques of the manipulator, Case I Figure 4.18: y-coordinate solution of the end-effector, Case I 109 Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 110 6000 4000 2000 - 1 1 / 1 1 '• > — - ^ 1, 1 / 1 ' / 1 1 1 1 1 1 1 1 •• o -2000 - 1 -4000 1 1 -6000 - -8000 1 1 2 time (s) Figure 4.19: Joint constraint torques of the manipulator for to = 0.5 Figure 4.20: y-coordinate solution of the end-effector in Case II Chapter 4. Simulation rtol l.d-4 l.d-5 l.d-6 l.d-7 l.d-8 of Constrained Rigid Multibody B(12,70) NFCN 920 1412 2180 3344 5156 PStab NFCN 716 1052 1580 2288 3578 B(12,70) NSTEP 153 253 363 557 859 PStab NSTEP 119 175 263 381 596 Dynamics 111 B(12,70) max-drift 1.3E-04 9.0d-06 5.3d-07 2.3d-08 2.3d-09 PStab max-drift 5.8d-08 9.5d-10 5.8d-12 1.5d-14 3.3d-15 Table 4.6: Comparison of B(12,70) with PStab on the manipulator in case I As shown in Table 4.6 for all specified tolerances PStab takes one third fewer function evaluations and steps than Baumgarte's stabilization. The drift on position level in Baumgarte's stabilization are much larger than those in PStab as shown in Table 4.6. The joint constraint torque are shown by Fig. 4.17. The solid line represents the joint constraint torque associated to 9\ and the dash line is associated to 02. Figs. 4.16 4.18 are the constraint path and y-coordinates of the end effector respectively. Case II A more difficult case [7] is obtained when the coordinates (x,y) are constrained to follow y = s'm2(uit) B(0,0) B(8,16) B(12,70) B(20,100) B(50,400) PStab NFCN 5416 4034 4058 4226 5258 3980 NSTEP 883 672 676 704 876 663 NREJCT 186 115 109 96 68 117 max-driftp 9.3E-03 3.9d-05 1.3d-05 9.5d-06 3.8d-06 1.3d-10 max-drift v 2.7d-03 1.0d-04 1.0d-04 8.8d-05 1.5d-04 2.4d-07 Table 4.7: Comparison of Baumgarte's parameters on the manipulator in case II where w is a parameter and the problem becomes difficult when LO is large. The Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 112 CRANK SLIDER damper ^nvvwv-': GROUND Figure 4.21: Slider-crank mechanism with a TSDA choice of Baumgarte's parameters is sensitive. For rtol = l.d-5, the results for different Baumgarte's parameters are shown in Table 4.7. In this example small Baumgarte's parameters (see B(8,16) and B(12,70)) are better than large ones. While PStab is better than all selected Baumgarte's stabilizations. The joint constraint torques in case II are shown by Fig. 4.19. Compared to case I the joint constraint torques are much larger than those in case I. Fig. 4.20 is ^/-coordinate of the end effector for UJ = 0.5. 4.4.4 Slider-Crank m e c h a n i s m w i t h a T S D A A slider-crank mechanism with a translational-spring-damping-actuator (TSDA) is investigated in [100]. As shown in Fig. 4.21 the mechanism has four bodies, GROUND, CRANK, ROD, SLIDER with three revolute joints and a translational joint [100]. Using absolute coordinates the generalized coordinates consist of 12 components (xi,yi,9i) indicating the Cartesian coordinate of center of mass and the orientation coordinate for body i (i = 1, • • • ,4). The constraint equations have the form [100] 9\ = si 92 = 3/i Chapter 4. Simulation of Constrained Rigid Multibody 93 113 Qi = 94 = xi- 9s V\ - 2/2 + sin(0 2 ) = Dynamics x2 + cos(0 2 ) 96 = x2 + 2 cos(0 2 ) — £3 + 2.5 cos(0 3 ) #7 = 2/2 + 2 sin(0 2 ) - J/3 + 2.5 sin(0 3 ) #8 = ^3 + 2.5 cos(0 3 ) - x4 99 = ?/3 + 2.5 sin(0 3 ) - y4 0io = («i - x 4 )sin(0i) - (j/! - y 4 ) c o s ( 0 i ) gii = cos(0i) sin(0 4 ) — sin(0i) cos(0 4 ) There is a constant torque (=100) in CRANK. A translational-spring-damping-actuator force element of spring coefficient 1000, damping factor 10 and zero actuator force is applied to the slider-crank mechanism. With the free-length 10, the TSDA is attached to the GROUND and the SLIDER, yielding applied forces given by fsda = ky where l0 = 10, I = (dTd)2 and / = }(dTd). ( _ /Q ) + cl The dist ance function is given by X\ + 12 cos(0i) — x4 — 2 cos(0 4 ) \ yi + 1 2 s i n ( 0 i ) - T / 4 - 2 s i n ( 0 4 ) / The other force acting on the mechanism is gravity. The equations of motion for the whole system can be written in the Euler-Lagrange equation form in which the mass matrix is a 12 x 12 constant diagonal matrix. In [100] a stabilized index-2 formulation of the constrained equations of motion [37] and its revised form were used for the crank-slider mechanism with a TSDA. A DAE solver DASSL was used for the integration in [100], and numerical simulations were carried out for 0.5 seconds. With different tolerances Chapter 4. Simulation of Constrained Rigid Multibody Dynamics Solution components of slider crank with a TSDA Time (seconds) Figure 4.22: Solutions of slider-crank with a TSDA, £=1000 Velocity components of slider crank with a TSDA Time (seconds) Figure 4.23: Velocities of slider-crank with a TSDA, fc=1000 114 Chapter 4. Simulation of Constrained Rigid Multibody Dynamics Solution components of slider crank with a TSDA 5i 1 1 1 1 Time (seconds) Figure 4.24: Solutions of slider-crank with a TSDA, £=10000 Velocity components of slider crank with a TSDA 1001 1 1 1 1 Time (seconds) Figure 4.25: Velocities of slider-crank with a TSDA, £=10000 115 Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 116 ranging from 10~ 3 to 10~ 6 numerical results reported in [100] and our similar runs with RKSTAB are contained in Table 4.8, i.e. where LM and CM represent the results of DASSL on the stabilized index-2 form and its revised form respectively. The tolerance settings for RKSTAB are the same as the corresponding runs with DASSL i.e. ATOL = RTOL = TOL. Clearly RKSTAB takes about (1/8) fewer steps than DASSL and the number of function evaluations in RKSTAB are about one third as those in DASSL. Method LM CM PStab LM CM PStab TOL 3 ioio- 3 10" 3 4 ioio- 4 IO -4 NFCN 251 242 62 328 327 86 NSTEP 84 86 10 128 119 14 TOL io- 5 10~5 io- 5 IO"6 IO"6 IO -6 NFCN 462 456 134 641 622 200 NSTEP 166 172 22 250 237 33 Table 4.8: Comparison of RKSTAB and DASSL on crank-slider with a TSDA The TSDA introduces oscillatory applied force, causing small oscillations of the numerical solutions. Figs. 4.22, 4.23 are the plots of numerical solutions for (#3, £3, j/3, #3, X4) and their velocities with spring coefficient = 1000. Figs. 4.24, 4.25 are those solutions and their velocities with spring coefficient = 10000. Clearly the larger spring coefficient causes more severe oscillation in the mechanism to appear. 4.5 R K S T A B on a singular slider-crank m e c h a n i s m A singularity occurs when the Jacobian matrix G becomes rank-deficient at some point in time. In such a case the optimization matrix (4.3) is singular and any underlying ODEs will yield an ill-posed problem. There are several methods to deal with singular descriptor form equations of motion [71, 9, 74]. Here we use RKSTAB to integrate a singular slider-crank mechanism which is considered in [74]. Although RKSTAB is not Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 117 Figure 4.26: A singular slider-crank mechanism designed for a singular system certain satisfactory results on this mechanism are obtained as well. Consider the slider-crank mechanism shown in Fig. 4.26 [74]. The equations of motion in descriptor form is given [74] by / 2.0 \cos{9d) co8(9d)\ 1.0 sm(ed)(92)2 fh' ) + 19.6 s i n ( ^ ) \ s i n ( ^ ) ( ^ ) 2 + 9.8sin(^) \92 / / sin(^) Vsin(^ 2 ; where 9& = 92 — 9\. As the lengths of the crank and slider are same the constraint has the form: 9 = - cos(0i) - cos(0 2 ) = 0 (4.27) The Jacobian matrix G = (sin(^i),sin(^ 2 )) becomes (0,0) as (91,92) takes the value (0,7r). A S in [74] we take the absolute tolerance = l.d-5 and simulate the mechanism for t G [0,3]. The numerical solutions and velocities without stabilization are plotted in Chapter 4. Simulation of Constrained Rigid Multibody Dynamics Figure 4.27: Solutions of the singular system with rtol = l.d-5 Time Figure 4.28: Solutions of the singular system with rtol = l.d-8 118 Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 119 0.35 0.25 0.15 0.05 Time Figure 4.29: Step sizes for the singular system without stabilization 0.35 / \ /' \ 0.3 / \ 0.25 0.2 0.15 / / / 0.1 / / / / x x * * * * * * * * x *X ,X><X. ''*..vx^ X xxxX XXX Xxx, X * X xxXxxv X * *. xxx . x x vy xx"^ xxx x > < * * 0.05 Figure 4.30: Step sizes for the singular system with PStab 3£ Chapter 4. Simulation of Constrained Rigid Multibody 5 i i i i Dynamics i 120 r 4 3 2 1 0 -1 -2 " 0 1 2 3 4 5 Time 6 7 8 9 10 Figure 4.31: Solution of the singular system with B(20,100), rtol=l.d-10 Fig. 4.27. The solid line and dashed line represent 9\ and 02 respectively. There are two singular points near t = 1 and t = 2.7 (where 9\ = 0) on the interval [0,3]. All the curves are smooth on the whole simulation interval with tolerance = l.d-5. If the tolerance is reduced to l.d-8 then the numerical solutions without stabilization will run into trouble as shown in Fig. 4.28. The reason that for smaller tolerances the numerical solution of the singular system could run into trouble lies in the fact that with smaller tolerances the code selects smaller step sizes and therefore increases the possibility for the mesh point and stage points to hit (close to) the singular points. We also test this problem with other tolerances. Fig. 4.29 shows the plots of step sizes used by RKSTAB for the tolerances Toll = l.d-4, Tol2 = l.d-5, Tol3 = l.d-6, Tol4 = l.d-7 and Tol5 = l.d-8. For the first singular point (i.e. near t = 1) all the runs pass it smoothly. Near the second singular moment, however, the simulation with Tol = l.d-8 has a problem as the code takes extremely small step Chapter 4. Simulation of Constrained Rigid Multibody 4 Dynamics ^ / / \ 121 r \ \ 3 \ / \ / 2 1 0 -1 -2 0 _J 1 I 2 L_ 3 _l 4 5 6 7 L_ 8 9 10 Time Figure 4.32: Solution of the singular system with PStab, rtol=l.d-10 sizes near the second singular moment. The optimization matrix will be evaluated at or near the singular point. Therefore the accuracy of the numerical solution cannot be guaranteed after the second singular moment. Stabilizations like Baumgarte's and the post-stabilization developed in the previous chapter can be useful when simulating a possibly singular mechanism although they are not developed for such a singular mechanism. Baumgarte's stabilization can be carried out for tolerances which are not smaller than l.d-10 for the time interval [0,10]. For tolerance = l.d-10, however, the simulation with Baumgarte's stabilization runs into trouble as shown in Fig. 4.31. PStab does a better job than Baumgarte's stabilization on a range of tolerances and time intervals for this singular crank-slider singular mechanism. Fig. 4.32 shows the solution with tolerance = l.d-10 and Tout = 10. Both Baumgarte's stabilization and PStab are proposed for forcing the numerical solutions to stay on or nearby the invariant manifold where the true solutions are supposed Chapter 4. Simulation of Constrained Rigid Multibody ~i r Dynamics —i i / / 0 ' » 2 4 122 \ \ 6 8 10 Time 12 14 _J L_ 16 18 20 Figure 4.33: Solution of new formulation for singular slider-crank to be on. For such a purpose Baumgarte's method is simply not as powerful as PStab. R e m a r k 4.4 The singularity in this singular slider-crank mechanism or similar ones can be cured by adapting better equations of motion in descriptor form for the mechanism. From Fig. 1^.26 we know that B1 + B2 = TT (4.28) since L\ = Li- Then using this simple, linear constraint equation the equations of motion for the mechanism can be written in the following 2.0 cos(0 d ) cos(9d) 1.0 '0,\ form /-sin(^)(^ 2 ) 2 + 19.6sin(^)\ /l\ sin(0d)(^)2 + 9.8sin(0) \l) / 01 + O2-T =O Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 123 where 9d = O2 — 61 as before. Notice that the Jacobian matrix of the new algebraic equation (4-^8) has a full row rank for any 61 and 62 therefore there is no singularity for all values of Qi and 62. Furthermore 62, our post-stabilization unlike (4-27), the algebraic equation (4-28) is linear in 6\ and PStab coincides with a projection method in this case. Fig. 4-33 shows the numerical solutions of the new formulation above using RKSTAB interval [0, 20] with absolute tolerance equal to l.d-10. on the time Of course, this is only a simple example, and in general we may not know how to reformulate the problem. Our point is that reformulations 4.6 that remove singularities may be possible in certain cases. A N o n h o l o n o m i c constrained s y s t e m Rolling D i s c on t h e Horizontal P l a n e A well-known example of a nonholonomic system is that of a circular disc set in rolling without sliding on a horizontal plane. We assume the disc has radius r and mass m. The position of the disc on the plane can be specified by generalized coordinates, the Cartesian coordinates x and y of the point of contact of the disc on the horizontal plane; the angle (f> between Ox axis and the straight line in which the plane of the disc intersects the horizontal plane; the angle 9 between the plane of the disc and Oz axis; and the angle of rotation if) of the disc measured in the plane of the disc from the radius through the point of contact to a fixed radius. The rolling without sliding assumption gives the two conditions dx = — r cos (f>dif) dy = — r sin 4>dif) where dx, dy and dip are arbitrary small displacements of the coordinates x, y and if). Chapter 4. Simulation of Constrained Rigid Multibody Dynamics Figure 4.34: A rolling disc on a horizontal plane Orbit of the rolling coin centre 0.8 0.60.40.2- 0 0 Figure 4.35: The orbit of the rolling coin center 124 Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 125 Therefore the kinematic constraints for the rolling disc are given by x = —r cos (/>?/>, y = — r sin <f>il> (4.29) The kinetic energy of the rotation motion is Tr = ]-A{92 + ft cos 2 9) + \c{^ + <> / sin Of where A and C are the principal central moments of inertia of the disc. Let (4.30) (xc,yc,zc) be the coordinate of the center of mass of the disc. The total kinetic energy of the disc is 1 T = m(xc +yc + zc ) + Tr -m(x2 + ij2 + r262 + r2<t>2 sin2 9) —mr[0 cos 9(x sin <f) — y cos <f>) -\- 4> sin 9(x cos </> + ?/ sin </>)] + ^A(92 + ft cos2 0) + l-C{^ + ^ sin 9)2 The Lagrangian function L = T — V = T — mgr cos 9. Let the generalized coordinate q = (x, y, </>,#, ?/>). Then the equation of motion for the rolling disc can be obtained by using the Euler-Lagrange equations: T d_ (dL\ dt\dq) _ dL _ / 1 0 0 0 cos, dq~[01Q0 gin( Ai (4.31) \ A2 ) The actual forces in the rolling disk problem are the gravity force and the contact force which acts on the disk and the plane through the contacting point. We set up the initial velocity and position as follows: 'o^ I 0 0.1 q(0) = 2 v(0) 0.1 0 -0.1 v0/ -0.1 ^ Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 126 drifts for post stabilization x10 0.2 0.4 0.6 0.8 ~1 T" 1 1.2 1.4 1.6 1.8 Figure 4.36: Drifts of rolling disk with post stabilization drifts for Baumgarte stabilization x10* 0 0.2 0.4 0.6 0.8 n r 1 1.2 1.4 1.6 1.8 Figure 4.37: Drifts of rolling disk with Baumgarte's stabilization Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 127 contact points in the plane 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 -0.05 0 0.05 0.1 0.15 0.2 0.25 Figure 4.38: The orbit of contact points in the plane z-coordinate of the coin center 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Figure 4.39: The height of the rolling disk center 1.8 Chapter 4. Simulation of Constrained Rigid Multibody Dynamics 128 Using RKSTAB we set up the error tolerance equal to l.d-3. We first try our poststabilization method proposed in Chapter 3 for nonholonomic constraint systems. As shown in Fig. 4.36 the drift with our post-stabilization is small. Then we switch to Baumgarte's stabilization. Fig. 4.37 is the drift with Baumgarte parameter 7 = 20. The drift for Baumgarte stabilization is especially large around t = 1.6 at which the rolling disk is going to rest on the horizontal plane. Fig. 4.35 to Fig. 4.39 are the plots of the orbit of the rolling disk's center in 3-dimension, the orbit of the contact point between the disk and the horizontal plane and the z-coordinate of the rolling disk's center. For multibody systems with nonholonomic constraints the equations of motion yield semi-explicit index-2 DAEs. Compared to holonomic constraint systems which are index3 DAEs nonholonomic constraint systems have less difficulties in numerical integrations. Furthermore, when the (nonholonomic) constraints are linear in the generalized velocity [68], our post-stabilization (3.76) proposed in Chapter 3 coincides with the following projection method at each successful advance integration step q = q v = v — Dfj, g(q)v = 0 Therefore the numerical drifts with our post-stabilization for nonholonomic constraint systems are very small. Chapter 5 P o s t - S t a b i l i z a t i o n in a D e f o r m a b l e M u l t i b o d y S y s t e m 5.1 Introduction During the past decade many analytical and numerical techniques have been developed for multibody systems that consist of interconnected rigid components [4, 5, 7, 29, 31, 46, 96, 99]. Advances in machine design frequently require higher operating speeds and greater positioning accuracy for a system with lightweight elements in it. At high speeds the rigid body assumption is often inadequate because of oscillatory elastic motion [51]. The demand for an accurate mathematical model, that accounts for the flexibility effect, has been a motivation for various analytical and numerical investigations for deformable multibody systems. The rigid body dynamics and structural analysis can be viewed as special cases of the dynamic analysis of deformable bodies that undergo large translational and rotational displacements. The study of dynamics of deformable multibody systems is complicated as the coupling between the rigid body motion and the elastic deformation has to be considered [19, 20, 51, 52]. This coupling can be significant when lightweight components, operating at high speeds, are involved. In the past, it was assumed that a small deformation has no significant effect on the gross rigid body motion. The system dynamical responses are thus identified by two different and uncoupled sets of equations, one which governs the gross rigid body motion of the system while the other describes the elastic deformation of the flexible components in the system [88, 89]. The equations 129 Chapter 5. Post-Stabilization in a Deformable Multibody System 130 of motion which define the rigid body dynamic system are first formulated and solved using a rigid body dynamic model. The calculated inertia, and sometimes the Coriolis and centrifugal forces, are then used with the externally applied forces to solve the linear structural dynamic problem for each flexible component [51, 52]. The total motion of the system is then obtained by the superpositioning of the elastic deformations of the flexible bodies on the system gross rigid motion. The accuracy of this approach is questionable since it does not provide for dynamical coupling of the rigid and elastic motions [51, 52]. Analysis procedures developed by Yoo and Haug [101] and Agarwal and Shabana [2] involve formulation of the equations of motion of each elastic body in terms of its "absolute" rigid body and elastic degrees-of-freedom. The configuration of each deformable body in a multibody system is identified by using two sets of coordinates: reference and elastic coordinates [85]. Reference coordinates define the location and orientation of a selected body reference. Elastic coordinates, on the other hand, describe the body deformation with respect to the body reference. The inertia coupling between the reference motion and the elastic deformation is identified and written in terms of the coupled sets of reference and elastic coordinates. Like that in a rigid system, the interconnection conditions of the adjacent substructures in the system are described by a set of constraint (nonlinear) equations formulated for each type of joints [46, 85]. This leads to a mixed system of differential and algebraic equations that must be solved simultaneously. In this approach all coordinates (reference and elastic) are maintained in the final form of dynamic equation. As high positioning accuracy is required in deformable multibody systems, the tolerances for a numerical drift can be very small. Therefore PStab could be an ideal algorithm for simulations of deformable multibody systems. In this chapter we apply PStab in a simulation of a deformable multibody system. In Section 2 we present the Chapter 5. Post-Stabilization in a Deformable Multibody System 131 concept of generalized coordinates for a deformable multibody system. There are basically two kinds of methods for simulating a deformable system: the assumed modes method and the finite element method [95, 85]. Using the finite element method considered in [85] we set up the equations of motion for a deformable slider-crank mechanism in Section 3. The equations of motion yield semi-explicit index-3 DAEs and the algebraic equations involve both the reference coordinates and elastic coordinates. The effects of dynamical coupling between the reference and elastic coordinates are reflected in the mass matrix and certain force terms in the equations. In Section 4 we consider simulations of the deformable slider-crank mechanism using our code RKSTAB. The motion of the deformable slider-crank mechanism is mainly caused by an initial angular velocity and an applied torque on the rigid crank in the mechanism. For relatively small initial angular velocities we also model the mechanism into a rigid system for numerical comparison. We consider the performance of Baumgarte's stabilization and PStab on both the rigid and deformable models. In Baumgarte's method the behavior of numerical drifts appears quite differently in the rigid and deformable models and large numerical drifts occur in both models. With PStab the numerical drifts are small for the whole time interval for both the rigid and deformable models. 5.2 G e n e r a l i z e d c o o r d i n a t e s in deformable m u l t i b o d y s y s t e m s In the study of dynamics of deformable multibody systems, different sets of coordinates are used to identify the configuration of a flexible multibody system. The so-called body coordinate describes the location and orientation of the deformable body with respect to the inertial frame. Usually the body coordinate is written in partitioned form as qj = (RT, 8*? (5.1) Chapter 5. Post-Stabilization in a Deformable Multibody System o 132 x Figure 5.1: A planar deformable body where Rl and 8l are the translational and rotation coordinates of the zth body axes. For simplicity of notation we consider a planar motion case. As shown in Fig. 5.1, the global location of an arbitrary point p on the deformable body i is then given by [85] r; = R + A\tfQ + u}) (5.2) where Rl is the location of the origin of the body axes X^X^ with respect to the global frame X Y and A1 is the transformation matrix from the body axes to the global axes. UQ is the position of point p in the undeformed state and iiy- represents the deformable displacement at point p measured with respect to the body frame. In general the equations of motion for a deformable body can be derived from constitutive equations (A.9) and equations of equilibrium (A.7) in Appendix A. Therefore the equations of motion are normally a set of nonautomous partial differential equations. The exact solutions of the displacement field for a deformable structure have an infinite number of degrees of Chapter 5. Post-Stabilization in a Deformable Multibody System 133 freedom. An approximation is sought through replacing the displacement field by certain shape functions Sl [85]: u) = S""(x)q)(t) (5.3) where 5*!(x) is a matrix of spatial admissible functions, x is the space variable and c?Jt) is the vector of temporal elastic generalized coordinates of body i. Using the aforementioned notations we can write the global position of an arbitrary point p on body i as r ; = R + A*(uj, + 5"q>) (5.4) Therefore generalized coordinates for body i are normally defined as those of qj, and q^. Depending on the choice of the spatial function ^ ' ( x ) , the commonly used approaches for a deformable system are: assumed modes method and finite elements method [95, 85]. In an assumed modes method, the shape function ^ ' ( x ) spans the entire body [85]. The shape function should satisfy the differentiability requirement and be consistent with certain boundary conditions. When structures have simple geometrical shapes, such as rods, beams and plates, the shape functions (5.3) in an assumed modes method can be determined without too much difficulty. In a multibody system with complex geometrical substructures, however, difficulty may be encountered in defining the assumed shape functions. In those cases the admissible functions are usually obtained by the finite element method [95]. In a finite element method, a deformable body in the system is discretized to relatively small regions called elements which are rigidly interconnected at selected nodal points [21, 85]. The shape functions S* in (5.3) are thus defined locally on each element as polynomials. Within each element the deformation can then be described by interpolating these polynomials. The coefficients of these polynomials are defined in terms of physical coordinates called the element nodal coordinates or elastic coordinates, which describe the displacements and slopes of selected nodal points on the element [21, 85]. By using Chapter 5. Post-Stabilization in a Deformable Multibody System 134 I i I 1 1f4 Figure 5.2: A two-dimensional beam element the connectivity between elements, the assumed displacement field (5.3) can be written in terms of the element shape functions and nodal coordinates of the structure. T w o - D i m e n s i o n a l B e a m E l e m e n t and its D i s p l a c e m e n t Field In the next section we will use a finite element method to discretize a deformable connecting rod in a slider-crank mechanism. The elements used in the discretization are two-dimensional beam elements. As shown in Fig. 5.2 a two-dimensional beam element has two nodal points. Each node has three nodal coordinates representing the longitudinal and transverse displacements ( q ^ q ^ or q/4?q/5 ) an d the slope (qj- 3 or qj 6 ) at the nodal point. Therefore six elastic degrees of freedom are associated with each two-dimensional beam element. Like that considered in [85, 86], the longitudinal displacement is approximated by a linear polynomial and the transverse displacement is approximated by a cubic polynomial in each beam element. The displacement field for the j t h beam element then has the form ( f i j ^ f = uj = S'qj where q} = ( q / i , q / 2 , q / 3 , q / W / 5 5 q / 6 ) T (5.5) Chapter 5. Post-Stabilization in a Deformable Multibody System 135 Figure 5.3: A slider-crank mechanism with deformable connecting rod are the nodal coordinates in the element mentioned above. The matrix S3' is the spacedependent element shape function [86] given by / 3 S = l-C 0 0 £ 0 0 where £ = x/lj and lj is the length of the j t h two-dimensional beam element. 5.3 E q u a t i o n s of m o t i o n for t h e deformable slider-crank m e c h a n i s m The slider-crank mechanism shown in Fig. 5.3 consists of a rigid crank centered at O, a deformable connecting rod AB of uniform circular section with length /, and a sliding block located at B [89]. The mass of the sliding block is ignored. The elastic connecting rod is assumed initially straight and simply supported at the ends by pins. Transverse and axial deformations of the connecting rod in the plane is allowed. There is a torque at Chapter 5. Post-Stabilization in a Deformable Multibody System 136 point 0 given by ri\ = sin(t) — <f)'. In this section we consider the equations of motion for the deformable mechanism using the finite element method [85, 86]. To avoid the tedious calculations we use MAPLE [66] to conduct the symbolic computations involved in the formulation by Lagrange's equations. We also use MAPLE to generate the necessary user specified Fortran subroutines called in our code RKSTAB. A d i s c r e t i z a t i o n of t h e c o n n e c t i n g rod in t h e m e c h a n i s m The connecting rod AB in the slider-crank mechanism can be discretized in a number of ways by using the finite element method. A simple discretization of rod AB is to divide rod AB into several two-dimensional beam elements. For a similar slider-crank mechanism where the crank rotates at a specified constant angular speeds, discretizations of the connecting rod into one, two and three equal elements are considered in [94] and good agreement in the deformations of the connecting rod is obtained using two and three finite elements. The situation here is a little bit more complicated as the rigid crank is driven by a given applied torque which depends on the derivative of the angle of the crank and therefore we have to treat the angle of the crank as an unknown variable. The gross reference motions of the connecting rod in our consideration are similar to those considered in [94]. We discretize the connecting rod AB into two equal beam elements with the length lj = 1/2 in our simulations of the deformable slider-crank mechanism. Although our discussion focuses on the two elements discretization for the connecting rod, the approach we use is general and can be used for any number of elements discretizations of the connecting rod AB. When the connecting rod is discretized into two beam elements there are 12 elastic coordinates. As shown in Fig. 5.4 the two beam elements in rod AB are rigidly connected. The connectivity conditions for the two beam elements can be written in terms of the Chapter 5. Post-Stabilization in a Deformable Multibody *5 <]f5 <*£ T6A t Of! %4 1f2 ^ 137 1 <te> tf System t - Qf| %8 <Jf5 «*! T ^T"^^ ^ qfV * Figure 5.4: Connection condition between two beam elements elastic coordinates as (5.6) q p = q} 6 To remove the redundancy among the elastic coordinates we define the total vector of nodal coordinates for the deformable connecting rod as 2 \T q = (q}i, q} 2 5 q/ 3 , q/4> q/5» q/e> q/4> q/s, q? 6 Then the nodal coordinates in each beam element can be written in terms of the total vector of nodal coordinates by introducing Boolean matrices B1 and B2 : q} = B'q where B\ = (-^6,06x3) and B2 = {06x3,h)- Similarly we write the displacement field on each beam element in terms of q u} = SjBJq Chapter 5. Post-Stabilization in a Deformable Multibody System 138 In the deformable slider-crank mechanism we assume that the connecting rod AB is simply supported at point A which implies q/i = q } 2 = q?5 = ° Thus we define the elastic generalized coordinates for the connecting rod as q/ - (q/35 q/4> q/s 5 q/6> q/4> q/e) The displacement field on each element can be written in terms of q j by introducing another Boolean matrix Bc : (u^vjf = ujf = SjBjBcqf (5.7) for j = 1,2 and i? c is given by 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 R e m a r k 5.1 The use of the Boolean matrices BJ and Bc leads to a similar of displacement field on each beam element. is written in terms of the shape function generalized coordinates. in the two-dimensional calculations the discretization Thus, the displacement field on each element S'3, (constant) Boolean matrices and the elastic When the finite elements used are identical i.e. with equal length beam elements, the shape function SJ is the same. Therefore regarding the integrations used in other elements. expression and products by MAPLE on one element can be These Boolean matrices also provide a systematic of rod AB with more elements. many approach for Chapter 5. Post-Stabilization in a Deformable Multibody System 139 M a s s m a t r i x of t h e deformable c o n n e c t i n g rod Let q r = ( i ? i , i ? 2 , ^ ) T be the set of reference coordinates that define the location of the origin of the body coordinates and the orientation of X1X2 with respect to the global frame X Y as shown in Fig. 5.3. The configuration of the connecting rod AB is determined by q r and q/. The global position of an arbitrary point p on element j of the connecting rod can be written as r3p = R + A (uj + S3B3Bcqf) where uJo = 0/ I A V-R2. cos 9 — sin 9 = sin 0 cos 9 The velocity of r3 with respect to X Y has the form = R + A89u1 + AS3B3BC q/ / R\ = (/2, Aeu3, AS3B3BC) 9 U// where u3 = UQ + S3B3Bcc\f and Ag = dA/d9. Therefore the kinetic energy of the connecting rod is given by Ti = Y.T^lX.L^W h 0TM)T. if Vi \ symm Aeu3 AS3B3BC u'Tu'' ujTIS3B3Bc Bj B3T S3T S3 B3 Bc /R\ dV3 9 Vq// Chapter 5. Post-Stabilization in a Deformable Multibody System 140 (R\ = -(RT,9,ig)M (5.8) 6 W where I = AgA is a constant matrix. Matrix M in (5.8) is called the mass matrix of the connecting rod and it has the following structure, / M \ / mRR A*u3 AS3B3BC u3Tu3 ujTIS3B3Bc mRf \ symm. dV3 Bj B^T S3T S3 B3 Bc symm. mR8 \ \ m/j / In M the coupling between the gross motion and the elastic deformation is implied through matrices mRg and m$f which depend on 9 and q/ nonlinearly. R e m a r k 5.2 The mass matrix M can be assembled by MAPLE matrices B3 and Bc. Consider the matrix m / y for mff with the use of Boolean instance. = Y, I . pBjBjTS3TS3BjBcdV3 j = YBjB3T f pSjTS3dV3B3Bc JVi 3 = where S is a 6 x 6 constant polynomials using YBIBJTSBJBC matrix which can be obtained by integrations of certain MAPLE. Elastic force and stiffness m a t r i x Using the assumptions of the elementary beam theory and neglecting rotary inertia Chapter 5. Post-Stabilization in a Deformable Multibody System 141 and shear deformation, the strain energy of the element j can be written [85, 86] as u / „„ 2 1 rh 3 „• \ 2" 2 Jo = W'l \£MY- where E3 is the modulus of elasticity, a3 is the cross-sectional area, I3 is the second moment of area. From (5.7) the strain energy for the connecting rod can be written as u = Y,uj B3 Bc<\jdx\ j l*Tf(XBlB>TK}fB3Bc)qf = 1 T (5.9) The matrix Kfj in (5.9) is called the stiffness matrix [21, 85]. The elastic force vector acting on the connecting rod is [85] 1 Q fc = - 0 0 0 ^ 0 0 0 e w 0 0 K Si / Constraints in t h e deformable slider-crank m e c h a n i s m The constraints in the deformable slider-crank mechanism are given by R\ — r cos (j) = 0 (5.10a) R2 — r sin (j) = 0 (5.10b) ^2 + (/ + q/(2) + q/(5))sin0 + q/(3)cos0 = 0 (5.10c) Chapter 5. Post-Stabilization in a Deformable Multibody System 142 where 4> is the angle of the crank. The condition (5.10c) indicates that the y-coordinate of slider block at B is zero. E q u a t i o n s of M o t i o n for t h e deformable slider-crank m e c h a n i s m The generalized coordinates for the deformable slider-crank mechanism are defined as (5.11) q = (#i,.R2,0,q/, The total kinetic energy for the mechanism is given by T = T2 + J<fi2/2 where J is the moment of inertia of the crank. The equations of motion for the deformable slider crank mechanism can be derived by using Lagrange's equation [85] |H + G T A = Qfc + Qe dt \dq) (5.12) where Q e is the vector of externally applied forces, Q^ is the vector of elastic force, G is the Jacobian matrix of the constraint (5.10), A is the vector of Lagrange multipliers. Equivalently we can write (5.12) in the following form. mRR mRe mr/ 0 (R"\ meg mef 0 6" / 0 0 Q+ GT\ (5.13) •Kff<lf \ symm. J 1 V 0 where Q is the collection of the generalized applied forces vector which includes the gravity, applied torque at O and those force terms that arise from differentiating the kinetic energy with respect to time and with respect to the generalized coordinates. 5.4 S i m u l a t i o n s of t h e deformable slider-crank m e c h a n i s m In the deformable slider-crank mechanism the connecting rod is considered to be a circular steel bar. Chapter 5. Post-Stabilization in a Deformable Multibody System 143 0.02 0.015 0.01 0.005 -0.005 - -0.01 - -0.015 -0.02 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Figure 5.5: Transverse deflection at the midpoint for </>'(0) = 180 The modulus of elasticity E = 3 0 x l 0 6 lb/in 2 ; The mass density p = 0.0007331 lb-sec 2 /in 4 ; The area moment of inertia about neutral axis is then given by / = irD4/6Ain4 where D is the diameter of the connecting rod. Other parameters of the deformable slider-crank mechanism are as follows. radius of the crank r length of connecting rod / diameter of connecting rod D moment of inertia of the crank J 6 in 12 in 0.25 in 1.0 Table 5.1: Slider-crank parameters Using RKSTAB we integrate the index-1 form of the equations of motion (5.13) along with the algebraic equations (5.10). We consider the deformations of the connecting rod Chapter 5. Post-Stabilization in a Deformable Multibody System x10 0.02 0.03 0.04 0.05 0.06 0.07 Figure 5.6: Transverse deflection at the midpoint for ^'(0) = 120 x10 0.02 0.03 0.04 0.05 0.06 0.07 Figure 5.7: Transverse deflection at the midpoint for </>'(0) = 60 Chapter 5. Post-Stabilization in a Deformable Multibody System 145 4 3 2 1 0 -1 -2 -3 -4 0 0.5 1 1.5 2 2.5 3 3.5 4 Figure 5.8: Drifts of deformable model with B(30,100) for ^'(0) = 180 in the mechanism when high velocities of the crank are involved. In the mechanism the crank must be operated at a speed such that the maximum amplitude of vibrations of the connecting rod is within the allowable design limit [94]. We therefore compute the deformations of the midpoint in the connecting rod for different initial velocities of <f>. The amplitudes of transverse vibration at the midpoint of the connecting rod for three different initial speed levels are shown in Figs. 5.5, 5.6 and 5.7. As shown in Figs. 5.5 and 5.7 the maximum viberational amplitude at initial speed </>'(0) = 180rad/s is more than ten times the maximum vibration amplitude at the initial speed <^'(0) = 60rad/s. The maximum deflection at the midpoint for <^'(0) = 120 is about three times the maximum deflection at the midpoint for <^'(0) = 60. We then consider simulations of the mechanism in the time interval [0, 4]. The absolute error and relative error tolerances in RKSTAB are l.d-4 and l.d-3. For </>'(0) = 180 the numerical drifts in Baumgarte's stabilization with the Baumgarte's parameters Chapter 5. Post-Stabilization in a Deformable Multibody System 146 .X10 T^ 1.5 2.5 3.5 Figure 5.9: Drifts of deformable model with PStab for <^'(0) = 180 (30, 100) increase and eventually become fairly large as shown in Fig. 5.8. With PStab the numerical drifts are small and appear to have no increasing tendency in the time interval as shown in Fig. 5.9. For a relatively small initial velocity <f>'(0) = 60 we also model the slider-crank mechanism into a rigid multibody system by only considering the reference variables (i?i, i?2, &, </>)• The equations of motion for the rigid model is simpler and RKSTAB takes much fewer steps for the rigid system for the same time interval and similar tolerance settings. The optimal Baumgarte's parameters are different for the deformable and rigid models and the behavior of numerical drift with Baumgarte's method is quite different. As shown in Fig. 5.11 the Baumgarte's stabilization has relatively large numerical drifts for large time t in the deformable model, while for the rigid model relatively large numerical drifts occur at early time as shown in Fig. 5.12. The different behavior of numerical drifts also indicates the sensitivity of the Baumgarte's stabilization in different problems. Chapter 5. Post-Stabilization in a Deformable Multibody System 147 -0.4 -0.8 0.5 1 2 2.5 3 3.5 4 Figure 5.10: Drifts of deformable model with B(16,90) for <^'(0) = 60 Figure 5.11: Drifts of deformable model with PStab for </>'(0) = 60 Chapter 5. Post-Stabilization in a Deformable Multibody System Figure 5.12: Drifts in the rigid model with B(10,80) for <//(0) = 60 x10 Figure 5.13: Drifts in the rigid model with PStab for ^'(0) = 60 148 Chapter 5. Post-Stabilization in a Deformable Multibody System 149 With PStab, however, the numerical drifts are small both in the rigid and deformable models. Since for the deformable model RKSTAB takes much smaller step sizes than those for the rigid model, smaller numerical drifts appear with PStab in the deformable model than those in the rigid model. The simulations of deformable multibody systems are more complicated than those for rigid systems. Also the appearance of the elastic coordinates could make the equations of motion stiff. Therefore stiff solvers could do a better job in certain cases. There are many interesting problems in the simulation of deformable multibody systems. For instance, 1. to find more efficient formulations of equations of motion for deformable multibody systems with large translations and rotations [83]; 2. to find efficient transformation to conduct velocity transformations [20] between the Cartesian coordinates and the relative coordinates; 3. to find efficient algorithms of linear algebra and other numerical schemes for deformable systems. Our investigation for deformable multibody systems has so far focussed on controlling the numerical drifts to a desirable small amount. We will continue our study in a number of directions with the goal of fast and efficient simulations for deformable multibody systems. Chapter 6 Conclusions High index DAEs can be regarded as limits of very stiff ODE systems and are usually ill-posed problems [17, 12]. In many cases index reduction with stabilization plays an important role in the numerical solution of semi-explicit higher index DAEs. In simulations of constrained multibody systems, the equations of motion are index-3 DAEs and an index reduction with stabilization yields better conditioned reformulations for numerical integrations. In particular when there is no other stiffness an explicit scheme can be used and this can lead to an efficient integration. In the reformulations the asymptotical stability of the invariant are important in the behavior of the numerical drifts [12, 8]. In the thesis we start with index reduction with stabilization for semi-explicit index-2 and index-3 DAEs. By investigating the underlying vector field of a given DAE we derive a general form of stabilization formulations with the aid of certain Lyapunov functions. Our considerations are made for any vector field with invariants and therefore can be applied to any semi-explicit high index DAEs. The stabilized forms of vector fields we proposed are different from the existing stabilization methods and improvements with our formulations are shown by numerical examples. In the stabilization formulations a choice can be made not only of certain parameters but also of the stabilization matrix. For semi-explicit index-2 DAEs the new stabilization formulations include Baumgarte's method with a special choice of the stabilization matrix. The invariant manifold defined by the algebraic equations is asymptotically stable in all of our stabilization formulations. Furthermore the invariant manifold has the monotonic property along the trajectories 150 Chapter 6. Conclusions of the stabilized vector field. 151 This monotonic property does not necessarily hold in Baumgarte's stabilization for semi-explicit index-3 DAEs. When a numerical discretization is applied the stability property of the invariant can be totally different in the difference equations. In Chapter 3 we present an example which shows that an asymptotically stable invariant manifold becomes unstable after a discretization. This is one of the reasons why large drifts can also occur with Baumgarte's stabilization, and we then propose to change this situation. In fact the difference equations resulting from a discretization of a vector field may have a much richer dynamic behavior than their continuous counterpart [98]. We regard the difference equations obtained from a discretization as a discrete dynamical system. Inspired by the stabilization for a vector field with an invariant manifold considered in Chapter 2, we propose stabilization methods for the invariant manifold in the discrete dynamical system. Since the stabilization is made after a discretization step we call it a post-stabilization. With proper choice of stabilization matrix in a post-stabilization we prove that the manifold is uniformly stable in the stabilized difference equations. Therefore a post-stabilization method restores the stability of the invariant manifold in the difference equations and thus very small numerical drifts can be expected in the numerical solutions. Combined with a discretization scheme with local error 0(hp+1) are 0(h2(p+1>) we show that the numerical drifts in certain post-stabilization schemes. For the Euler-Lagrange equations in multibody systems our choice of post-stabilization is PStab which is described in Algorithm 3.3. Numerical examples show that PStab is efficient and the amount of CPU time spent on the post-stabilization is around 1% of the total CPU time spent on the entire integration. We then develop a numerical code RKSTAB for simulations of multibody dynamics based on PStab. The discretization scheme is an embedding Runge-Kutta formula which is modified from the code DOPRI(5,4) [43]. The local error estimation in the modified Chapter 6. Conclusions 152 scheme uses the difference between the stabilized 5th order solution and a modified 4th order solution. In RKSTAB a minor modification allows us to include Baumgarte's stabilization as an option. The Baumgarte parameters are left to the user. With RKSTAB we conduct a variety of simulations of multibody systems including planar systems, spatial systems, systems involving high oscillation, systems with singularities and a system with nonholonomic constraints. These numerical experiments reported in Chapter 4 show that RKSTAB is efficient, fast and reliable. In Chapter 5 we apply PStab to a deformable multibody system. Using the finite element method we discretize a deformable connecting rod in a deformable slider-crank mechanism. The equations of motion are obtained from Lagrange equations. We use MAPLE to formulate the equations of motion and generate the necessary user-specified Fortran subroutines called in RKSTAB. Numerical integration of the equations of motion for the deformable slider-crank mechanism is carried out by RKSTAB. We compute the relative deformations at the midpoint of the connecting rod for different initial velocities of the crank in the mechanism. Larger initial velocities of the crank result in larger deformations in the connecting rod. Numerical drifts in Baumgarte's stabilization and PStab are considered for different initial velocities of the crank. For a relatively small initial velocity of the crank we compare the numerical drifts in Baumgarte's stabilizations for both the deformable and rigid models. With PStab small numerical drifts result in the simulations while large numerical drifts always occur with Baumgarte's method for both the deformable and the rigid models. This shows that the performance of Baumgarte's stabilization is very sensitive to the problem formulation. For all problems considered PStab is robust and efficient. D i s c u s s i o n of Future Work Efficient simulation of large constrained m u l t i b o d y s y s t e m s The need to simulate and control the dynamics of multibody systems arises in many Chapter 6. Conclusions 153 applications such as vehicle motion simulation, robot control, mechanism design and physically based modeling in computer graphics and vision. A particularly demanding application involves person-in-the-loop simulations, such as those involving force-reflecting haptic interfaces. In the simulation applications large systems with sparse structures have to be solved frequently. The efficiency of the simulation can be improved by incorporating suitable stabilization methods, integration schemes as well as certain sparse system algorithms. A continuation of this thesis research is needed to do efficient simulations of constrained multibody systems in order to approach the goal of real-time simulation. S i m u l a t i o n s and analysis of s y s t e m s w i t h variable structures Many industrial systems such as robots are subjected to a change in their kinematic structure during the functional operations [19, 20]. Mass captures, mass releases and additions or deletions of mechanical joints are possible sources of these mechanical structure changes. As time progresses, different sets of constraint equations are imposed, thus the dimension and rank of the system constraint Jacobian matrix depend on time, and accordingly a sequence of configuration spaces which have different dimension and basis has to be defined. The trajectories prescribed by different bodies in the multibody system depend on a set of degrees of freedom that has variable dimension. When there is a change in the kinematic structure, impacts between system components and jump discontinuities in the system velocity vector can occur [19, 20]. Simulations and dynamical analysis of multibody systems with variable kinematic structures are needed. S i m u l a t i n g c o n s t r a i n e d , deformable m u l t i b o d y s y s t e m s It is well-known that in many applications the rigid body assumption is violated when high speed and lightweight substructures are involved in the system. In the interest of geometric flexibility a finite element discretization (rather than some assumed modes methods) has to be used. The deformations at the finite element nodes become part of the generalized coordinate vector. The effect of coupling between elastic deformation and the Chapter 6. Conclusions 154 gross rigid body motion is then implied in the equations of motion. Usually the dimension of equations is large as the model is essentially a discretized system of partial differential equations. The matrices involved are all large and sparse. Drifts between the adjacent substructures can be large when using Baumgarte's stabilization. Moreover, efficient finite element formulations are needed for describing the motion of deformable bodies that undergo large angular rotations. Large systems with highly oscillatory solutions typically result. Their numerical solution is a challenge. Appendix A A Brief R e v i e w of M e c h a n i c s of D e f o r m a b l e B o d i e s The theory of elasticity is concerned with the mechanical behavior of solids on the macroscopic scale [85]. It treats material as uniformly distributed throughout regions of space. So the quantities like density, displacement, velocity can be defined as continuous functions of position as well as time. In a rigid case, the distance between any two particles does not change and the whole rigid body motion can be identified by the linear velocity of any point in the body and the angular velocity. A force acting at a point in a body can be equivalent to the same force acting at other point in the body plus a moment. In a deformable case, however, continuum mechanics is concerned with the changing shape and relative motion of the particles. A force that acts at a point is equivalent to a system of forces defined at another point [87]. The system force consists of force, a moment that depends on the relative displacement between the two points, and a set of generalized elastic forces that depends on the finite rotation of the body [87]. Strain c o m p o n e n t s Consider a fixed Cartesian coordinate system X 1 X 2 X 3 with origin 0 shown in Fig. A.l. Suppose at time t = 0 a deformable body occupies a fixed region of space BQ and a point PQ £ BQ. At time t the body occupies a new continuous region of space B and PQ moves to P. The displacement of P can be written as u = £-x 155 (A.l) Appendix A. A Brief Review of Mechanics of Deformable Bodies 156 Figure A.l: Deformed and Undeformed State of a Flexible Body Then the strain components in the deformable body are defined [85] by 1 / dui 2 \dxi + duj dx; J ^ duk duk' ,-x dxi dx. fc=i WU/l "-"J/ (A.2) where i, j = 1, 2, 3. The strain components arise naturally in the analysis of the deformation. They are nonlinear functions of the spatial derivatives of the displacement. Because of the symmetry of the strain tensor [21] it is sufficient to identify the strain vector e, e = ( e n , e22, €33, ei2, £i3, £23) (A-3) In case of small strains and rotations, the product terms in (A.2) can be neglected [21, 85]. If we switch the notation of Cartesian coordinate system to X Y Z the strain vector in (A.3) can be approximated as e = Du (A.4) Appendix A. A Brief Review of Mechanics of Deformable Bodies 157 where D is a derivative operator given by / 2 f 0 0 0 2^ 0 ox dy D 1 2 2 f 0 0 a_ dy _9_ dx 0 d dz 0 _9_ dx 0 9 3^ dy I Stress c o m p o n e n t s In continuum mechanics two forces are encountered: surface force and body force [21]. Body forces are those which are developed without physical contact and are distributed over the volume of the continuum. Typical examples of body forces are the gravitational and electromagnetic forces. Surface forces, or the other hand, include all forces acting on the boundaries of a medium through direct contact. Stresses in a medium result from surface forces which act on some portion of the medium [28]. Stresses describe the manner in which forces acting on the boundaries of the medium are transmitted through the medium. To be specific, let P be a point in the body and SS be the area of an element of the surface that consists P. The force SF which acts on SS may resolved into two components: 6Fn and SFt, the normal and the tangential components to the area. A normal stress an and a shear stress are then defined [28] as SFn an = lim ss^o SS ,. SFt at = lim —— ss^o SS (A.5) In a Cartesian coordinate system X Y Z we may consider the stresses acting on planes whose outward normals are in X, Y or Z directions. The state of stress at a point can be described completely by specifying the stresses acting on the mutually perpendicular Appendix A. A Brief Review of Mechanics of Deformable Bodies 158 Ax Figure A.2: Notation for stress planes through the point [28]. Therefore the stress at a point can be specified by a 3 x 3 matrix defined as [28] / V V XX T~xy 'xz Tyx °yy T 'zx T Vzz \ (A.6) yz zy J The equation of equilibrium of the body then can be written [85] as d<?x dx <9i yx dx drzx dx d-.xy + + dy da,yy dy dr..zy dy + + + di dz + hi <9T, yz 62 dz + fi dazz 63 dz + fl = pa 1 (A.7a) pa-2 (A.7b) pa-3 (A.7c) where a\, a2 and 03 are the components of the acceleration vector, and /&i, / ^ and fb3 are the body force components. It can be shown that the stress tensor matrix (A.6) is Appendix A. A Brief Review of Mechanics of Deformable Bodies 159 symmetric and the stress component vector is defined as [85] & — I, ^ r a i &yy j O' zz i Txy •> T~xz i Tyz ) \ A. O) Constitutive equations The stress and strain vectors are not sufficient to describe the mechanical behavior of a deformable body. Body deformations depend on the applied forces and the force-displacement relationship [21]. The deformation also depends on the material of the body. An additional equations known as constitutive equation, are required to com- plete the specification of the mechanical properties of a material. Constitutive relations serve to distinguish one material from another [85]. For an isotropic material, the stress components can be expressed in terms of the strain components as [85] °xx = n , .-U1 ^ [ ( 1 - 2v)exx + i/et] (1 + " ) ( 1 - -2v) E -2v) (1 + E [(1 - 2v)eyy + utt] ~ zz , w N [fl (1 + */)(l — 2i/) u 'xy ^^-xyi 'yz — — ^ tyz 5 2v)e,2 + vtf\ ; 'zx J — ^^zx where E, G and v are the modulus of elasticity, modulus of rigidity and Poisson's ration respectively. et = txx -f tyy + tzz. For most materials, it is found that the strains are proportional to the applied forces [85], provided the load does not exceed the elastic limit. Mathematically it can be written as a = Ee (A.9) where E is a 6 X 6 constant matrix and it is called the matrix of the elastic constant of the material. Usually (A.9) are called the c o n s t i t u t i v e equations. Bibliography [1] R. P. Agarwal, Difference Equations cations. Marcel Dekker, Inc. 1992. and Inequalities. Theory, Methods, and Appli- [2] Agrawal, 0 . P., and Shabana, A. A, Dynamic analysis of multibody systems using component modes, Computers and Structures, Vol. 21, No. 6(1985), 1303-1312. [3] T. Alishenas, Zur numrischen Behandlung, Stbilisierung durch Projektion und Modellierung mechanischer Systeme mit Nebenbedingungen und Invarianten, PhD thesis, Institut fur Numerische Analysis und Informatik, Konigliche Technische Hochschule, S-100 44 Stockholm, Schweden, 3/1992. [4] F.M.L. Amirouche, Computational Methods in Multibody Dynamics, Englewood Cliffs, New Jersey, 1992. Prentice-Hall, [5] Th. Andrzejewki, H. G. Bock, E. Eich and R. v. Schwerin, Recent Advances in the Numerical Integration of Multibody Systems, Preprint 93 - 28, Juni 1993. [6] U. Ascher, On numerical differential algebraic problems with application to semiconductor device simulation, SIAM J. Numer. Anal. 26 (1989), 517-538. [7] U. Ascher, H. Chin, L. Petzold and S. Reich, Stabilization of constrained mechanical systems with DAEs and invariant manifolds, the Journal of Mechanics of Structures and Machines, 23(2), 135-157(1995). [8] U. Ascher, H. Chin and S. Reich, Stabilization Numer. Math. 67 (1994) 131-149. of DAEs and invariant manifolds, [9] U. Ascher and P. Lin, Sequential regularization methods for higher index DAEs with constraint singularities: I. Linear index-2 case, SIAM J. Numer. Anal, to appear. [10] U. Ascher, R. Mattheij and R. Russell, Numerical Solution of Boundary lems for Ordinary Differential Equations, Prentice-Hall, 1988 Value Prob- [11] U. Ascher and L. Petzold, Projected implicit Runge-Kutta methods for algebraic equations, SIAM J. Numer. Anal. 28 (1991), 1097-1120. differential- [12] U. Ascher and L. Petzold, Stability of computational methods for constrained namics systems, SIAM J. Scient. Comp. 14 (1993), 95-120. 160 dy- Bibliography 161 [13] U. Ascher and L. Petzold, Projected collocation for higher-order higher-index differential-algebraic equations, J. Comp. Appl. Math. 43 (1992), 243-259. [14] A. K. Banerjee and M. E. Lemak, Multi-flexible body dynamics capturing induced stiffness, Transactions of the ASME. 58 (1991), 766-775. motion- [15] J. Baumgarte, Stabilization of constraints and integrals of motion in dynamical tems, Comp. Methods Appl. Mech. 1 (1972): 1-16 [16] H. G. Bock and R. v. Schwerin, An Inverse Dynamics ADAMS-Method strained Multibody Systems, Preprint 93 - 27, Juni 1993. sys- for Con- [17] K. Brenan, S. Campbell and L. Petzold, Numerical Solution of Initial-Value lems in Differential-Algebraic Equations, North-Holland, 1989. Prob- [18] S.L. Campbell and B. Leimkuhler, Differentiation of Constraints in Differential gebraic Equations, CRSC Technical Report 062498-01 January 4, 1990. Al- [19] C. W. Chang and A. A. Shabana, Spatial dynamics of deformable multibody systems with variable kinematic structure: part 1 - dynamic model, Journal of Mechanical Design, 112 (1990), 153-159. [20] C. W. Chang and A. A. Shabana, Spatial dynamics of deformable multibody systems with variable kinematic structure: part 2 - velocity transformation, Journal of Mechanical Design, 112 (1990), 160-167. [21] R. D. Cook, Concepts and applications of finite element analysis, 2nd Ed., Wiley, New York, 1981. [22] J. J. Dongarra, C. B. Moler, J. R. Bunch, G. W. Stewart, UNPACK SIAM, Philadelphia (1979). [23] J.R. Dormand and P.J. Prince, A family Comp. Appl. Math. 6 (1980), 19-26. of embedded Runge-Kutta User's Guide, formulae, J. [24] E. Eich, Convergence results for a coordinate projection method applied to mechanical systems with algebraic constraints, SIAM J. Numer. Anal., to appear. [25] E. Eich, K. Fiihrer, B. Leimkuhler and S. Reich, Stabilization and projection methods for multibody dynamics, Research report, Inst. Math., Helsinki Univ. of Technology, 1990. [26] A.Eichberger, C. Fiihrer and R. Schwertassek, The Benefits Simulation and its Application to Vehicle Dynamics of Parallel Multibody Bibliography 162 [27' W.H. Enright, Numerical methods for systems of initial value problems - The state of the art, In: Haug, E.(ed.): Computer aided analysis and optimization of mechanical system dynamics, Springer, Berlin Heidelberg New York Tokyo(1984) [28 R. W. Fox and A. T. McDonald, Introduction Sons, Inc. 1985. to Fluid Mechanics, John Wiley & [29; K. Fuhrer and B. Leimkuhler, Numerical solution of differential-algebraic for constrained mechanical motion, Numer. Math. 59 (1991), 55-69. equations [30 K. Fuhrer and B. Leimkuhler, A new class of generalized inverses for the solution of discretized Euler-Lagrange equations. In: Haug, E., Deyo, R.(ed.): NATO Advanced Research Workshop on Real-Time Integration Methods for Mechanical System Simulation, Snowbird, Utah 1989. Springer, Berlin Heidelberg New York Tokyo(1990) [31 K. Fuhrer and B. Leimkuhler, Formulation and numerical solution of the equations of constrained mechanical motion, DFVLR-FB 89-08, Munich, Germany, 1989. [32 C. Fuhrer and R. Schwertassek, Generation and solution of multibody system equations, J. Non-Linear Mechanics 25 (1990), 127-141. [33 Wei-Hsin Gau and A. A. Shabana, Effect of finite rotation on the propagation of elastic waves in constrained mechanical systems, Transactions of the ASME, 114 (1992), 384-393. [34 C.W. Gear, Differential-algebraic equations, indices and integral algebraic SIAM J. Numer. Anal. , 27 (1990). [35 C.W. Gear, Maintaining solution invariants in the numerical SIAM J. Sci. Stat. Comp., 7 (1986), 734-743. solution [36 C.W. Gear, H. H. Hsu and L. Petzold, Differential-algebraic Proc. ODE Meeting, Oberwolfach, West Germany, 1981. equations equations, of ODEs, revisited, [37 C.W. Gear, G. Gupta and B. Leimkuhler, Automatic integration of the EulerLagrange equations with constraints, J. Comput. Appl. Math. 12(1985), 77-90. [38 G.H. Golub and C.F. Van Loan, Matrix Computations, Press, Baltimore, 1983. Johns Hopkins University [39 E. Griepentrog and R. Marz, Differential-Algebraic Equations Treatment, Teubner-Texte Math. 88, Teubner, Leipzig, 1986. [40 Friedemann Grupp and Bernd Simeon, Five-Link-Suspension and their Numerical , Benchmark-Example Bibliography 163 [41] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Applied Mathematical Sciences, Vol. 42, SpringerVerlag New York Inc. 1983. [42] E. Hairer, Ch. Lubich and M. Roche, The Numerical Solution of DifferentialAlgebraic Systems by Runge-Kutta Methods, Lecture Notes in Mathematics, vol. 1409, Springer-Verlag, 1989. [43] E. Hairer S.P. N0rsett and G. Wanner, Solving Ordinary Differential Nonstiff Problems, Second Revised Edition, Springer-Verlag, 1993. Equations [44] E. Hairer and G. Wanner, Solving Ordinary Differential Differential-Algebraic Problems, Springer-Verlag, 1991. II: Stiff [45] Hale, J. K. Ordinary Differential Equations I: and Equations, Wiley: New York, 1969. [46] E.J. Haug, Computer-Aided Kinematics and Dynamics of Mechanical System Basic Methods, Allyn and Bacon, 1989. Vol.1: [47] Haug, E., Deyo, R.(ed.) NATO Advanced Research Workshop on Real-Time Integration Methods for Mechanical System Simulation, Snowbird, Utah 1989. Springer, Berlin Heidelberg New York Tokyo(1990) [48] Edwar J. Haug, Jon G. Kuhl and Fuh Feng Tsai, Virtual Prototyping for System Concurrent Engineering, 1992. [49] M. Hirsch, C. Pugh, and M. Shub, Invariant matics. No. 583, Springer-Verlag, 1976. Mechanical manifolds, Lecture Notes in Mathe- [50] W. C. Hurty, Dynamic analysis of structural system using component modes, AIAA, J. Vol.3, No 4. (1965) 678-685 [51] S. K. Ider and F. M. L. Amirouche, Nonlinear modeling of flexible multibody systems dynamics subjected to variable constraints, Transactions of the ASME, Vol. 56, (1989), 444-450. [52] S. K. Ider and F. M. L. Amirouche, Influence of geometric nonlinearities in the dynamics of flexible treelike structures , J. Guidance, Vol. 12, No. 6 (1989), 830-837. [53] A. Jain, Unified formulation of dynamics for serial rigid multibody systems, Journal of Guidance, Control, and Dynamics, 14 (1991), 531-542. [54] J. J. Jiang, C.L. Hsiao and A.A. Shabana, Calculation of non-linear variation of rotating beams by using tetrahedral and solid finite elements, J. Sound and Vibration 148(1991), 193-213. Bibliography 164 [55] J. B. Jonder and J. P. Meijaard, SPACAR-Computer program for dynamic analysis of flexible spatial mechanisms and manipulators, In Multibody System Handbook (W. Schiehlem ed.), Springer-Verlag, 1990. [56] Ben Jonker, A finite element dynamic analysis of spatial mechanisms with flexible links , Computer Methods in Applied Mechanics and Engineering 76 (1989) 17-40. [57] T. R. Kane and R. R. Ryan, Dynamics base, J. Guidance, 10 (1987), 139-151. of a cantilever beam attached to a moving [58] Y. A. Khulief and A. A. Shabana, Dynamics of multibody systems with variable kinematic structure, J. Mechanism, Transmissions, and Automation in Design, 108, (1986), 167-175. [59] J. D. Lambert, Numerical Methods for Ordinary Differential [60] J. P. Lasalle The Stability of Dynamical Applied Mathematics, SIAM, 1976. Systems, Systems, Wiley, 1991. Regional Conference Series in [61] V. Lakshmikantham and D. Trigiante Theory of Difference Equations: Numerical Methods and Applications Mathematics in science and engineering, Vol.181, Academic Press, 1986. [62] B. J. Leimkuhler, L.R. Petzold and C. W. Gear, On the consistent initialization of differential-algebraic systems of equations, SIAM J. Numer. Anal. 28 (1991), 105-226. [63] G. Leister and D. Bestle, Symbolic-numerical solution of multibody systems closed loops, Vehicle System Dynamics, 21 (1992), 129-142. with [64] P. Lotstedt and L. Petzold, Numerical solution of nonlinear differential equations with algebraic constraints I: Convergence results for backward differentiation formulas, Math. Comp., 46 (1986), 491-516. [65] Ch. Lubich, U. Nowak, U. Pohle and Ch. Engstler, MEXX - Numerical Software for the Integration of Constrained Mechanical Multibody Systems, Preprint SC 92-12 (December 1992). [66] MAPLE V, Language Reference Manual, Springer-Verlag, 1991. [67] R. Marz, Practical Liapunov Stability Criteria for Differential Preprint No. 91-28. Berlin 1991. [68] J. I. Neimark and N. A. Fufaev, Dynamics Mathematical Society, 1972. of Nonholonomic Algebraic Equations, Systems, American Bibliography 165 [69] P. E. Nikravesh and I. S. Chung, Application of Euler parameters to the dynamic analysis of three-dimensional constrained mechanical systems, J. Mechanical Design, Vol. 104, (1982), P 785-791. [70] Orn Olafsson and Taifun Alishenas, A comparative study of the numerical integration and velocity stabilization of two mechanical test problems, TRITA-NA-9210. [71] K.C. Park and J.C. Chiou, Stabilization of computational prodedures for constrained dynamical systems, J. Guidance, Control and Dynamics II (1988), 365-370. [72] L. Petzold, Differential/algebraic put., 3 (1982), 367-384. equations are not ODEs, SIAM J.Sci. Stat. Corn- [73] L. R. Petzold and P. Lotstedt, Numerical solution of nonlinear differential equations with algebraic constraints II: Practical implications, SIAM J. Sci. Stat. Comput., 7(1986), 720-733. [74] L. Petzold, Y. Ren and T. Maly, Regularization of higher-index differential-algebraic equations with rank-deficient constraints, Preprint, AHPCRC, University of Minnesota, 1993. [75] P.J. Prince and J.R. Dormand, Higher order embedded Runge-Kutta Comp. Appl. Math. 7 (1981), 67-75. formulae, J. [76] P.J. Rabier and W.C. Rheinboldt, A general existence and uniqueness theory for implicit equations, Differential and integral equations. Vol.4, No.3, May 1991, 563582. [77] W.C. Rheinboldt, Differential-algebraic systems as differential folds, Math. Comp., 43 (1984), 473-482. [78] R.E. Robersom and R. Schwertassek, Dynamics Verlag, 1988. equations on mani- of Multibody Systems, Springer- [79] Can E. Rosenthal, An order n formulation for robotic systems, The Journal of the Astronautical Sciences, 38, (1990), 511-529. [80] SD/FAST User's Manual, Version B . l . l , Symbolic Dynamics, Inc. 1991. [81] W. Schiehlen (Editor): Multibody Systems Handbook, Springer, 1990. [82] W. Schiehlen (Editor): Advanced Multibody Systems Dynamics, Simulation ware Tools, Kluwer Academic Publishers 1993. and Soft- [83] A. A. Shabana, Constrained motion of reformahle bodies, International journal for numerical methods in engineering, 32 (1991), 1813-1831. Bibliography 166 [84] A. A. Shabana, Automated analysis of constrained systems of rigid and flexible bodies, J. Vibration, Acoustics, Stress, and Reliability in Design, 107, (1985), 431-439. [85] A. A. Shabana, Dynamics 1989. of Multibody Systems, John Wiley and Sons, New York, [86] A. A. Shabana, Theory of Vibration, Springer-Verlag, New York, 1991. Vol II: Discrete and Continuous Systems, [87] A. A. Shabana and C. W. Chang, Connection forces in deformable multibody namics , Computer & Structures, 33, 1, (1989), 307-318. dy- [88] A.A. Shabana and R.A. Wehage, A coordinate reduction technique for transient analysis of spatial substructures with large angular rotation, J. Struct. Mech. 11 (1983), 401-431. [89] A.A. Shabana and R.A. Wehage, Variable degree-of-freedom component mode analysis of inertia flexible mechanical systems, J. Mech. Trans. Auto. Desi. 105. (1983). [90] A.A. Shabana and R.A. Wehage, Spatial transient analysis of inertia-variant mechanical systems, J. Mech. Tran. and Auto. Desi. 106, 1984. flexible [91] L.F. Shampine, Conservation laws and the numerical solution of ODEs, Comp. and Maths, with Appls., Part b, 12, (1986), 1287-1296. [92] L.F. Shampine and M.K. Gordon, Computer Solution of Ordinary Differential tions: The Initial Value Problem, W. H. Freeman and Company, 1975. Equa- [93] B. Simeon, F. Grupp, C. Fiihrer and P. Rentrop, A nonlinear truck model and its treatment as a multibody system, preprint. [94] J. O. Song and E. J. Haug, Dynamic analysis of planar flexible mechanisms, puter Methods in Applied Mechanics and Engineering, 24 (1980), 359-381. [95] A. Suleman, Dynamics and Control of Evolving Space Platform: Applications , Ph.D. Thesis, UBC, 1992. A Formulation Com- with [96] R.A. Wehage and E. J. Haug, Generalized coordinate partitioning for dimension reduction in analysis of constrained dynamic systems, J. of Mechanical Design 104 (1982), 247-255. [97] T.Y. Yang, Finite Element Structural Analysis, Prentice-Hall, 1986. Bibliography 167 [98] H. C. Yee, P. K. Sweby and D. F. Griffiths, Dynamical approach study of spurious steady-state numerical solutions of nonlinear differential equations. 1. The dynamics of time discretization and its implications for algorithm development in computational fluid dynamics NASA Technical Memorandum 102820. April 1990. [99] J. Yen, Constrained equations of motion in multibody dynamics as ODEs on manifolds, SIAM J. Numer. Anal. 30 (1993), 553-568. [100] J. Yen and L. R. Petzold, On the numerical solutions of constrained dynamics system, Preprint, University of Minnesota, July 15, 1994. multibody [101] W. S. Yoo and E. J. Haug, Dynamics of flexible mechanical systems using vibration and static correction modes , ASME Design Engineering Technical Conference , 1985, Cincinnati, Ohio.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Stabilization methods for simulations of constrained...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Stabilization methods for simulations of constrained multibody dynamics Chin, Hong Sheng 1995
pdf
Page Metadata
Item Metadata
Title | Stabilization methods for simulations of constrained multibody dynamics |
Creator |
Chin, Hong Sheng |
Date Issued | 1995 |
Description | The descriptor form of constrained multibody systems, and any general formulation of such systems with closed loops, yield index-3 differential-algebraic equations (DAEs). Generally, some index reduction techniques have to be used before numerical discretization can be safely applied. However, a direct differentiation of the constraints introduces mild instability. Hence one must consider index reduction with stabilization. One popular method for index reduction with stabilization is Baumgarte's technique. However the difficulty of choosing the parameters in practice makes this method's robustness unclear. Moreover our numerical experiments show that there are still large constraint drifts even with Baumgarte's stabilization. In the thesis, we employ concepts and techniques of dynamical systems in order to improve the situation. We first study the relationship between a DAE and its underlying vector fields. A general form of vector fields with stabilized invariant manifolds is given. We propose a new numerical stabilization method for semi-explicit index-2 and index-3 DAEs of Hessenberg form. Our stabilization method improves on Baumgarte's stabilization which is widely used in engineering as well as in simulations of multibody systems. We then develop a numerical code based on the new stabilization method with an adaptive step-size control for the descriptor form of the Euler-Lagrange equations in multibody systems. Numerical simulations have been conducted with our code on a variety of multibody systems including a spatial five-link-suspension model in a vehicle, Andrews squeezing mechanism, a two-arm manipulator with a prescribed motion and a mechanism with kinematic singularity. Our code is efficient, fast and therefore is more attractive for real-time simulations. When high speed and light-weight substructures are involved in the multibody system, the rigid body model is usually no longer valid. In such a case we compute the elastic deformations and oscillations of the substructures using the finite element method. Satisfactory numerical results using our method are presented for deformable multibody systems as well. |
Extent | 6204021 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-06-04 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080031 |
URI | http://hdl.handle.net/2429/8776 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1995-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1995-982691.pdf [ 5.92MB ]
- Metadata
- JSON: 831-1.0080031.json
- JSON-LD: 831-1.0080031-ld.json
- RDF/XML (Pretty): 831-1.0080031-rdf.xml
- RDF/JSON: 831-1.0080031-rdf.json
- Turtle: 831-1.0080031-turtle.txt
- N-Triples: 831-1.0080031-rdf-ntriples.txt
- Original Record: 831-1.0080031-source.json
- Full Text
- 831-1.0080031-fulltext.txt
- Citation
- 831-1.0080031.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080031/manifest