Análisis cinemático y dinámico de sistemas mecánicos formados por varios sólidos rígidos

4 Revista internacional de métodos numéricos para cálculo y diseño en ingeniería, Vol. 1,4,31-48 (1985) ANALISIS CINEMATICO Y DINAMICO DE SISTEIMAS MECANICOS FORMADOS POR VARIOS SOLIDOS RIGIDOS JORGE UNDA Y JAVIER GARCIA DE JfiLON Centro de Estudios e Investigaciones Ticnicas de Guipúzcoa Escuela Superior de Ingenieros Industrhdes de San Sebastián Universidad de Navaira RESUMEN En este artículo se presenta un nuevo método numérico de análisis cinemáticlo y dinámico con computador de sistemas mecánicos formados por varios só!idos rígidos. El método presentado utiliza un nuevo sistema de coordenadas no independientes fonnado por las coordenadas cartesianas de algunos puntos del mecanismo y las componentes cartesiarias de algunos vectores cinitarios solidarios al mismo, mediante los que se define la posición y el movimiento del sistema. Las ecuaciones de restricción que ligan estas coordenadas se establecen coino condiciones de sólido rígido de cada elemento y como restricciones de par. La inclusión de componentes cartesianas de vectores unitarios dentro de este sistema de coordenadas, facilita en gran manera la formulación de las restricciones de par cuando el par cinemático está asociado con una dirección particular, como sucede en los pares de revolución (R), cilíndricos (C) o prismiíticos (P). Las ecuaciones de restricción resultan lineales o cuadráticas. Las ecuaciones diferenciales del movimiclnto se obtienen fácilmente mediante el Teorema de las Potencias Virtuales. Finalmente, se presentan algunos ejemplos de análisis de mecanismos planos y tridimensionales. SUMMARY In this paper a new method for the nunierical kinematic and dynamic analysis of multi-rigid-body systems is described. The method presented uses a neur system of non independent coordinates formed by the cartesian coordinates of some points of the mechanism and by the cartesian components of unitary vectors fixed to it, which determine the position and the motion of the multi-rigidbody system. The constraint equations arise from the rigicl body condition of each, element and from the pair constraint equations. The inclusion of unitary vectors as mechanism coordinates allows an easy formulation of pair constraints when the pair is related to a particular direction, as in revolute (R) cylindrical (C) or prismatic (P) pairs. The constraint equations are always linear or quadratic. The differential equations of motion are obtained easily through the application of the Theorem of Virtual Power. Some exarnples of dynamic analysis of planar and threedirnensional mechanisms are presented. INTRODUCCION Muchos sistemas mecánicos pueden ser modelizados de forma efectiva como sistemas formados por varios sólidos rígidos unidos ent:re sí. El análisis diinámico de estos sistemas puede ser llevado a cabo rnediante mét.odos analíticos o mediante métodos numéricos. Los métodos analíticos, como los detjcritos para sistemas tridimensionales Recibido: Marzo 1985 t3 Univasitat Politkcnica de Catalunya (España) ISSN 0213-1315 3 2 JORGE UNDA Y JAVIER GARCIA DE JALON en la referencia1, han aumentado considerablemente su campo de aplicación en los últimos años, pero las dificultades teóricas y prácticas son tan grandes que la única alternativa posible para lograr la generalidad en el análisis son los métodos numéricos. Estos métodos han abierto la posibilidad de desarrollar programas de computador generales, como los bien conocidos DRAM~,~, 1MP4, ADAMS~ y DADS-~D~-~. Las características distintivas de los métodos de análisis en que se basan estos programas son, por una parte, el tipo de coordenadas elegidas para la definición de la posición y movimiento del mecanismo; por otro lado, las ecuaciones de restricción que ligan, estas coordenadas, cuya formulación está estrechamente relacionada con el tipo de i coordenadas elegido; y finalmente, el establecimiento de las ecuaciones de equilibrio1 dinámico cuya integración numérica permitirá conocer la evolución del mecanismo a lo largo del tiempo. En este artículo se presenta un método de análisis cinemático y dinámico de mecanismos basado en la utilización de un nuevo sistema de coordenadas dependientes, denominadas "coordenadas naturales", cuyas ecuaciones de restricción son de fácil formulación. Las ecuaciones diferenciales del movimiento se establecen mediante el Teorema de las Potencias Virtuales, previa definición de las matrices de inercia de los sólidos rígidos que forman el mecanismo y posterior ensamblaje de éstas en la matriz de inercia del sistema. COORDENADAS DE UN MECANISMO El punto crucial de cualquier método de análisis cinemático y dinámico de mecanismos es la definición de las cccoordenadas" del mecanismo. Dichas coordenadas vienen constituidas por un conjunto de parámetros no independientes que definen unívocamente la posición de todos y cada uno de sus elementos; no son independientes porque cualquier conjunto de parámetros cuyo número sea superior al número de grados de libertad del mecanismo, deberá satisfacer ciertas ecuaciones de compatibilidad geométrica adicionales, que se conocen con el nombre de "ecuaciones de restricción". Las ecuaciones de restricción juegan un papel de fundamental importancia en el análisis de estos sistemas, y se corresponden estrechamente con el tipo de coordenadas elegido. Por otra parte, se definen las "coordenadas generalizadas" de un mecanismo como un conjunto de coordenadas independientes, cuyo número coincide con el número de grados de libertad del mismo. Las coordenadas generalizadas no determinan la posición de todos los elementos del mecanismo sino a través de la resolución del "problema de posición", que es un problema no lineal que puede tener varias soluciones. Por esta razón, las coordenadas generalizadas no se pueden utilizar por sí solas para definir la posición, sino que se suelen usar para definir las velocidades y aceleraciones de los elementos de entrada, o para la integración numérica de las ecuaciones diferenciales del movimiento. Así pues, la definición de las coordenadas del mecanismo y de las ecuaciones de restricción constituyen el núcleo de todo método numérico de análisis cinemático y dinámico de mecanismos. Existen tres tipos principales de coordenadas de un mecanismo, que han dado lugar a tres familias diferentes de programas de análisis con computador: las coordenadas "relativas", utilizadas por los programas DRAM2*3 e IMP4, las coordenadas de "punto de referencia", utilizadas con ligeras variantes por los programas ADAMS5 y DADS3DG8; y las coordenadas "básicas" o "naturales", de las que hace uso el programa COMPAMM9"13, y a las que se hará especial referencia en este artículo. Coordenadas relativas Como se ha apuntado, las "coordenadas relativas" se usan en los programas DRAM e IMP. Con ellas, la posición de cada elemento se define con relación al elemento que le precede en la cadena cinemática, por medio de un número mínimo de parámetros que dependen de los grados de libertad de movimiento relativo permitidos por el par que une ambos elementos. Así, un par esférico introduce tres ángulos como coordenadas relativas, un par cilíndrico introduce un ángulo y una distancia, y un par de revolución introduce un único ángulo. La Figura 1 muestra un mecanismo tridimensional RSCR en el que la posición de cada uno de sus elementos queda definida por un conjunto de coordenadas relativas. 33 ANALISIS CINEMATICO Y DINAMICO Figura 1 Las ecuaciones de restricción asociadas con 1z.s coordenadas relativas se obtienen formulando vectorial o matricialmente las ecuaciones correspondientes a los lazos cerrados del mecanismo, y por tanto, para poder formularlas es necesario el reconocimiento previo de los lazos independientes. Esta tarea se puede realizar automáticamente (por medio de la teoría de grafgs, por ejemplo), o por medio del analista. Los principales inconvenientes de las coorderiadas relativas derivan del hecho de que no determinan directamente la posición, velocidad y aceleración de los elementos, con las consiguientes dificultades a la hora de plantear las ecuaciones de equilibrio ,dinámico, y la necesidad de utilizar pre y postprocesadores. Coordenadas de punto de referencia Las coordenadas de punto de referencia definen la posición de cada elemento utilizando las coordenadas cartesianas de uno de sus puntos -usualmenite el centro de gravedady la orientación angular del elemento definida respecto de unos ejes fijos de alguna determinada manera. Existen para esto diversos procedimientos. Así, el programa ADAMS~ utiliza los "ángulos de Euler"', que constituyen un mínimo sistema de coordenadas angulares independientes, pero con los que pueden surgir algunas dificultades cuando el ángulo de nu,tación se hace cero o múltiplo de T, pues en esta 34 JORGE UNDA Y JAVIER GARCIA DE JALON posición los ángulos precesión y de rotación propia no están definidos unívocamente. Este inconveniente puede ser obviado cambiando la orientación de los ejes fijos. Otra forma de determinar la orientación angular de un sólido es por medio de los "parámetros de Euler", que son cuatro parámetros no independientes que fijan la posición del eje y la magnitud de un giro. Estas coordenadas no tienen posiciones singulares, pero como contrapartida su número es mayor que con los ángulos de Euler, y también el hecho de no ser independientes presenta dificultades adicionales en" la integración numérica de las ecuaciones diferenciales del movimiento. Estas coordenadas son las utilizadas en k1 programa DADS-~D~-' y en la referencia14. Otros autores han propuesto otros tipos de variables para fijar la posición angular entre un triedro fijo y otro móvil, como por ejemplo los nueve elementos de la matriz de cosenos directores de los ejes móviles, o bien los seis elementos correspondientes a dos columnas de dicha matriz15'16. Estas coordenadas tienen los mismos inconvenientes o aún mayores que los parámetros de Euler: número elevado y relaciones de dependencia. Las ecuaciones de restricción asociadas con las coordenadas de punto de referencia surgen de las restricciones impuestas en el movimiento relativo entre dos elementos contiguos por el par cinemática que los une. El número de coordenadas de punto de referencia, para un mismo mecanismo, es siempre muy superior al de coordenadas relativas, pero sus defensores arguyen que las e~uaciones matriciales que se obtienen son muy dispersas, y que el mayor número no les resta eficacia. Sin embargo, hay que señalar que -en la prácticaesto sí que constituye un verdadero inconveniente, por la gran cantidad de ecuaciones de restricción que hay que introducir y los sistemas de ecuaciones diferenciales tan grandes que resultan. Coordenadas naturales y coordenadas básicas Las coordenadas básicas constituyen otra opción en el análisis de mecanismos planosg-ll y tridimensi~nalesl~-~~. Estas coordenadas están formadas por las coordenadas cartesianas de algunos puntos del mecanismo. La Figura 2 muestra los "puntos básicos" del mismo mecanismo espacial RSCR de la Figura 1. Cada elemento debe tener al menos dos puntos básicos en el caso plano, y tres en el tridimensional. Figura 2 ANALISIS CINEMATICO Y DINIMICO 3 5 Las coordenadas básicas pueden considerarse como una evolución de las coordenadas de punto de referencia, en las que dichos puntos se: desplazan hacia los pares del mecanismo o a otros puntos de interés. De esta manera, cada elemento tiene varios puntos de referencia, y éstos son compartidos por varios elementos, lo que tiene la importante ventaja de hacer innecesario el uso de las coordenadas angulares. Las ecuaciones de restricción asociadas con las coordenadas básicas se obtienen por un doble camino: de la condición de sólido rígido entre puntos básicos pertenecientes a un mismo elemento, y de las res.tricciones en el movimiento relativo de puntos básicos pertenecientes a elementos adyacentes irnlpuestas por el par cinemática que los une. Todas estas restricciones pueden quedar reducidas al establecimiento de tres tipos de condiciones elementales: distancia constiinte entre dos puntos básicos, área constante del triángulo formado porr tres puntos básicos, y volumeri constante del tetraedro formado por cuatro puntos l~ásicos. Las coordenadas básicas son atractivas porque dan directamente la posición, velocidad y aceleración de los puntos de ~nterés del mecanismo, y porque las ecuaciones de restricción son simples y fáciles de manejar. Como desventaja de estas coordenadas aparece el gran número de incógnitas a que dan lugar en el caso tridimensional. Por ello, estas coordenadas son especialmente adecuadas para el caso plano, en el que su número es moderado (siempre menor que con las coordenadas de punto de referencia), y la formulación extraordinariamente sencilla. Las coordenadas naturales constituyen una evolución de las coordenadas básicas dirigida a solucionar sus desventajas en el caso tiridimensional. De igual manera que éstas, las coordenadas naturales utilizan las coorde:nadas cartesianas de algunos puntos del mecanismo, pero además de éstas utilizan ts.mbién las componentes cartesianas de algunos vectores unitarios. En la Figura 3 puede verse el mismo mecanismo de las figuras anteriores, cuya configuración queda definida por medio de coordenadas naturales. Figura 3 36 JORGE UNDA Y JAVIER GARCIA DE JALON Los vectores unitarios contribuyen a definir los elementos como sólidos rígidos, y proporcionan una forma fácil y simple de considerar los pares cinemáticos asociados con una determinada dirección, como son los pares R, C y P. Así, las coordenadas naturales mantienen las ventajas de las coordenadas básicas, reduciendo el número de parámetros necesarios para definir la posición del mecanismo. Como puede verse en las Figuras 1 a 3, el mismo mecanismo necesita de seis coordenadas relativas, dieciocho coordenadas básicas, y doce coordenadas naturales. El mismoi mecanismo necesitaría dieciocho coordenadas de punto de referencia, con ángulos de Euler: y ventiuna en el caso de utilizar parámetros de Euler, para su definición. Al utilizar coordenadas naturales los elementos del mecanismo quedan definidos mediante un conjunto de puntos básicos y vectores unitarios, existiendo diferentes posibilidades para la definición de los elementos, tal como puede verse en la Figura 4. Figura 4 Las ecuaciones de restricción asociadas a las coordenadas naturales son muy sencillas de establecer, y serán explicadas con detalle en el apartado siguiente. ECUACIONES DE RESTRICCION Como se ha indicado en el apartado anterior de este artículo, existe una estrecha relación entre las coordenadas del mecanismo y el tipo de ecuaciones de restricción que las interrelaciona. Así, las coordenadas relativas dan lugar a ecuaciones de restricción de lazo cerrado, las coordenadas de punto de referencia dan lugar a ecuaciones de restricción de par, y las coordenadas básicas a ecuaciones de restricción de sólido rígido (o de elemento) y de par. A continuación se detallará brevemente el modo en que se formulan las ecuaciones de restricción en coordenadas básicas para el caso plano, y en coordenadas naturales en el caso tridimensional. En ambos casos se distinguirá entre las restricciones que derivan de la condición de sólido rígido y las que derivan de los pares cinemáticos del mecanismo. ANALISIS CINEMATICO Y' DINAMICO 3 7 Restricciones de elemento En el caso plano, la condición de que un elemento sea un sólido rígido se reduce -por lo generala la condición de que se mantengan constantes las distancias entre los puntos básicos que pertenecen a dicho elemento. En algunos casos particulares, como cuando hay tres puntos básicos que estári alineados, las tres condiciones de distancia constante no son independientes, y hay que recurrir a un nuevo tipo de condición, que es la de área constante del triángulo formado por los tres puntos básicos. Matemáticamente, estas condiciones pueden formu:iarse en la forma siendo (xi,yi), (xj ,yj), (xk ,yk) las coordenadas cartesianas de los puntos entre los que se establecen las ecuaciones de restricción, dij el valor de la distancia entre los puntos i y j, y Aijk el valor del área del triángulo formado :por los puntos i, j, y k, En el caso tridimensional las ecuaciones de restricción se formulan también de un modo muy sencillo, estableciendo condiciones .de ángulo y distancia constantes entre puntos y vectores unitarios pertenecientes al mismo elemento. Así, en un elemento tridimensional definido mediante p puntos básicos y vectores unitarios, será necesario establecer al menos (3p-6) restricciones de elemento independientes, pues cada punto o vector aporta tres coordenadas y un sólido rígido en el espacio tiene seis grados de libertad. Las condiciones de distancia y ángulo constante se formulan mediante el producto escalar entre vectores unitarios, o entre un vector formado por dos puntos básicos y un vector unitario. Así, el elemento de la figura 5a, definido mediante dos puntos Figura S 3 8 JORGE UNDA Y JAVIER GARCIA DE JALON básicos, tiene seis coordenadas naturales y cinco grados de libertad (no se considera la rotación alrededor de la recta que une los dos puntos), y por tanto hay que introducir una ecuación de restricción que es la condición de distancia constante entre ambos puntos donde C1 es una constante y rij representa el vector asociado con el segmento (i-j). Análogamente, el elemento de la Figura 5b, definido mediante dos puntos básicos y un vector unitario, da lugar a las siguientes ecuaciones de restricción: donde C1 y C2 son constantes, rij tiene el mismo significado que en la expresión anterior y u, es un vector unitario. Finalmente, en Ia Figura 5c puede verse un elemento definido mediante dos puntos básicos y dos vectores unitarios. De acuerdo con la regla establecida con anterioridad,' las ecuaciones de restricción de elemento que habrá que introducir serán las seis siguientes: donde el significado de las variables es similar al de las expresiones 4-6. Cuando los dos vectores unitarios son coplanarios, las ecuaciones de restricción anteriores no son las adecuadas. La condición de producto escalar constante entre ambos vectoñes debe ser sustituida por la condición de que uno de ellos pueda expresarse como combinación lineal del otro y del vector rij. De estas tres ecuaciones algebraicas, únicamente una es independiente de las anteriores, y ésta puede ser determinada mediante la técnica del pivotamiento utilizada en la resolución del sistema de ecuaciones lineales resultante. De modo análogo, se pueden establecer las restricciones de sólido rígido para elementos mis complicados. Restricciones de par Una vez que los elementos rígidos han sido definidos como tales, quedan por establecer las restricciones correspondientes a los pares cinemáticos. ANALISIS CINEMATICO Y DINAMICO 3 9 En el caso bidimensional se considerarán los pares de rotación R, y de traslación P. En la Figura 6a aparece un par R uniendo dos elementos contiguos. Mediante las coordenadas básicas no hace falta introducir ninguna ecuación nueva, pues la restricción de dicho par queda automáticamente introducida por compartir ambos e1ement.o~ un único punto básico. Por otra parte, en la Figura 6b aparece un par P que devenga Figura 6 las condiciones de que sea constante el área del írriángulo formado por los puntos i, j y k, y de que el vector (9-k) sea perpendicular al vector (i-j). Matemáticamente, estas ecuaciones se formulan del siguiente modo, En el caso tridimensional, las coordenadas naturales permiten formular también muy fácilmente este tipo de restricciones, pues, de hecho, una característica esencial de las coordenadas naturales es la de desplazar el énfasis de las restricciones de par a las restricciones de sólido rígido, con lo cual la formulación se simplifica considerablemente. La restricción que en el movimierito relativo entre dos elementos impone un par esférico S (Figura 7a), es tenida en cuenta de modo automático cuando ambos comparten el punto básico localizado en dicho par. De igual manera, un par de revolución R (Figura 7b) no hace tarnpoco necesaria la formulación de ninguna ecuación de restricción, pues los elementos que une comparten el punto básico situado en el par y el vector unitario asociado con el eje de giro de dicho par. Por otra parte, un par cilíndrico C (Figura 7c) se introduce compartiendo ambos elementos un vector unitario en la dirección del eje del par, e imponiendo la condición de que los puntos i y j permanezcan alineados con dicho vector unitario. Esta condición se impone por medio de la proporcionalidad de las componentes de. ambos vectores . .- donde c es una constante para cada posición del mecanismo. De un modo análogo, y siempre relativamente sencillo, puedenintroducirse las restricciones correspondientes a otros pares: prismático, junta de Cardan, etc. 46 JORGE UNDA Y JAVIER GARCIA DE JALON la barra durante cinco segundos, representándose una posición cada cuatro centésimas de segundo. Los impactos han sido analizados con un coeficiente de restitución de 0.9 y con un coeficiente de rozamiento de 0.2. Figura 9 Las Figuras loa, 10b y 10c muestran, según diferentes perspectivas, la evolución de un mecanismo tridimensional SSC bajo la acción de la gravedad. El inovimiento se representa durante 1.2 segundos, dibujándose una posición cada 0.03 segundos. Figura 10 ANALISIS CINEMATICO Y DINA!VIICO 4 7 Las Figuras 1 la y 1 lb, por otra parte, muestran bajo diferentes puntos de vista la evolución estacionaria de un giróscopo en el que la velocidad de nutación es nula. Figura 1 1 Finalmente, las Figuras 12a y 12b representan la evolución del mismo giróscopo de las Figuras 1 la y 1 lb, siendo ahora la velocidad de nutación diferente de cero. Figura 12 CONCLUSIONES En este artículo se ha descrito un nuevo método para el análisis cinemático y dinámico de mecanismos. Este método utiliza como coordenadas del mecanismo, las coordenadas cartesianas de ciertos puntos y las componentes cartesianas de algunos vectores unitarios rígidamente unidos a los elementos del mismo. Estas coordenadas y sus ,derivadas determinan la posición, velocidad y aceleración de cada elemento del mecanismo, eliminando la necesidad de considerar coordenadas angulares. Como consecuencia de esto, las ecuaciones de restricción que ligan estas coordenadas no independientes 4 8 JORGE UNDA Y JAVIER GARCIA DE JALON son siempre cuadráticas o lineales, y por tanto fáciles de formular; además, nunca contienen funciones trascendentes. Para el problema dinámico se ha desarrollado una matriz de inercia, obteniéndose las ecuaciones de equilibrio dinámico mediante el Teorema de las Potencias Virtuales. Los Multiplicadores de Lagrange, utilizados para plantear las ecuaciones del movimiento, se eliminan mediante la proyección de dichas ecuaciones sobre el subespacio nulo del jacobiano de la matriz de restricciones. El método presentado es conceptualmente simple y fácil de implementar en un computador. REFERENCIAS 1. J. Duffy. ''Analysis ofMechanisms and Robot Manipulators". Edward Arnold, (1980). 2. M. A. Chace y D. A. Smith. "DAMNDigital Computer Program for the Dynarnic Analysis of Generalized Mechanical Systems". Dans. SAE 80,969-991, (197 1). 3. D. A. Smith, M. A. Chace y A. C. Rubens. "The Automatic Generation of a Mathematical Model for Machinery Sistems". Joumal of Engineering for Industry, 629-635, (1 973). 4. P. N. Sheth y J. J. Uicker. "IMP (Integrated Mechanisms Program), A Computer-Aided Design Analysis System for Mechanisms and Linkages" . Joumal of Engineering for Industry , 454464, (1972). 5. N. Orlandea, M. A. Chace y D. A. Calahan. "A Sparsity-Oriented Approach to the Dynarnic Analysis and Design of Mechanical Systems - Part 1 and Part 2". Joumal of Engineering for Industry, 773-784, (1 977). 6. R. A. Wehage y E. J. Haug. "Generalized Coordinate Parttioning for Dimension Reduction in Analysis of Constrained Dynamic Systems". Joumal of Mechanical Design 104,247-255, (1982). 7. R. A. Wehage y E. J. Haug. "Dynamic Analysis of Mechanical Systems with Intermittent Motion". Joumal ofMechanica1 Design 104,778-784, (1982). 8. P. E. Nikravesh y 1. S. Chung. "Application of Euler Parameters to the Dynamic Analysis of Three-Dimensional Constrained Mechanical Systems". Joumal of Mechanical Design 104,785791, (1982). 9. J. García de Jalon, M. A. Serna y R. Aviles. "Computer Method for Kinematic Analysis of Lower-Pair Mechanisms-1 Velocities and Accelerations". Mechanism and Machine Theory 16, 543-556, (1981). 10. J. García de Jalón, M. A. Sema y R. Avilés. "Computer Method for Kinematic Analysis of Lower-Pair Mechanisms-11 Position Problem". Mechanism and Machine Theoty 16, 5 57-566, (1981). 1 1. M. A. Serna, R. Aviles y J. García de Jalón. "D ynamic Analysis of Plane Mechanisms with Lower Pairs in Basic Coordinates". Mechanism and Machine Theory 17,397.403, (1 982). 12. J. García de Jalón, M. A. Serna, F. Viadero y J. Flaquer. "A Simple Numencal Method for the Kinematic Analysis of Spatial Mechanisms". Joumal of Mechanical Design 104,78-82, (1982). 13. A. Aveilo, J. Unda y J. García de Jalón. "Análisis Dinámico de Mecanismos Tridimensionales en Coordenadas Naturales". 111 Congreso Nacional de Ingeniería Mecánica, Gijón, España (1984). 14. J. Wittenburg. "Dynamics of Systems of RigidBodies", (B.G. Teubner, Stuttgart, 1977). 15. R. E. Roberson. "Correction of Numerical Error when Two Direction Cosine Columns Are Used as Kinematical Variables". Computer Methods in Applied Mechanics and Engineering 46, 15 1 - 158, (1984). 16. R. E. Roberson. "Correction ofNumenca1 Error when One Direction Cosine is Known". Computer Methods in Applied Mechanics and Engineenng 46,307-3 12, (1 984). 17. C. W. Gear. 'Differential-Akebraic Equations". Computer Aided Analysis and Optimization of Mechanical System Dinarnics. Spnnger Verlag, Heilderberg, (1984). 18. C. W. Gear. "Numerical Initial Value Problems in Ordinary Differential Equations ". PrenticeHaii, Englewood Cliffs, New Jersey, (1971). 19. G . E. Forsythe , M. A. Malcolrn y C. B. Moler. "Computer Methods for Ibhthematical Computations '[ Prentice-Hall, Englewood Cliffs, New Jersey, (1 977). 20. Harwell Subroutine Library. AERE Harwell. Oxfordshire. U.K., (198 1).