Dinámica de Mecanismos

Jesús Meneses Alonso

 ÍNDICE
 INTRODUCCIÓN. Dinámica inversa Vs. Dinámica directa
 COORDENADAS. Conjunto libre Vs. Conjunto ligado
 ECUACIONES DE NEWTON-EULER
► Sólido Rígido
► Mecanismo
 Ejemplo de aplicación: Mecanismo Biela-Manivela-Corredera
 Coordenadas Locales-Globales. Matriz de cambio

 ECUACIONES DE LAGRANGE
► Conjunto ligado de coordenadas. Multiplicadores de Lagrange
► Matriz de masas
► Vector de fuerzas generalizadas
 Simulación dinámica inversa
 Simulación dinámica directa.
2
 DINÁMICA INVERSA Vs. DIRECTA

Dinámica: estudio del movimiento considerando las fuerzas

 Dinámica INVERSA (DI):
► conociendo el movimiento y acaso algunas fuerzas aplicadas (pesos, etc)
► determinar el resto de fuerzas aplicadas
T(t)? (t) F(t)

 Dinámica DIRECTA (DD):

► conociendo las fuerzas aplicadas
► determinar el movimiento T(t) (t)? F(t)

 En un mecanismo también puede ser objetivo conocer las fuerzas

de enlace entre eslabones (son incógnitas tanto en DI como en DD)
3
 CONJUNTO DE COORDENADAS
Libre Vs. ligado

 Conjunto LIBRE: las coordenadas son independientes (n = g):

☺ Reducido número de ecuaciones
 No aparecen las fuerzas de enlace (entre eslabones)
 Poca generalidad
 Conjunto LIGADO: las coordenadas son dependientes (n > g) y
existen r ecs. de restricción (r = n – g)
 Mayor número de ecuaciones
 Aparecen las fuerzas de enlace
☺ Gran generalidad

4
 ECUACIONES DE NEWTON-EULER
Sólido Rígido

 Sólido rígido, caso general (6 ecs.): ෍ 𝐹Ԧ = 𝑚𝑎𝐺

► Ԧ resultante de las fuerzas externas al sólido

σ 𝐹:
෍ 𝑀𝐺 = 𝐼𝐺 𝛼Ԧ + 𝜔 × 𝐼𝐺 𝜔
► 𝑎𝐺 : aceleración del centro de masas del sólido
► σ 𝑀𝐺 : resultante de los momentos externos y momentos de las fuerzas
externas respecto al centro de masas del sólido
► 𝐼𝐺 : Tensor de inercia del sólido en el sistema de coordenadas local
Ԧ : Velocidad y aceleración angulares del sólido
► 𝜔, 𝛼

 Sólido rígido en movimiento plano (3 ecs.): ෍ 𝐹Ԧ = 𝑚𝑎𝐺

► Los ptos. realizan trayectorias planas (planos // entre sí)
෍ 𝑀𝐺 = 𝐼𝐺𝑧 𝛼
►Eje z ⊥ al plano de movimiento: 𝜔 = 𝜔𝑘; 𝛼Ԧ = 𝛼𝑘
Ԧ = 𝐹𝑥 𝑖Ԧ + 𝐹𝑦 𝑗Ԧ; 𝑀 = 𝑀𝑘
►𝐹

5
 ECUACIONES DE NEWTON-EULER
Mecanismo

 Mecanismo:
► Aplicar a cada eslabón i por separado (diagrama de sólido libre)
► Considerar las fuerzas, 𝐹𝑗𝑖 , y pares, 𝑇𝑗𝑖 , o esfuerzos de ligadura, ejercidos
por eslabones j distintos al i; además de los externos al mecanismo, 𝐹𝑖 , 𝑇𝑖 , o
esfuerzos aplicados.
► Considerar la ley de acción-reacción: 𝐹𝑘𝑖 = −𝐹𝑖𝑘 ; 𝑇𝑘𝑖 = −𝑇𝑖𝑘 .Por ejemplo,
eliminando los esfuerzos de subíndices ki, con k>i:

𝐹𝑖 + ෍ 𝐹𝑗𝑖 − ෍ 𝐹𝑖𝑘 = 𝑚𝑖 𝑎𝐺 𝑖
𝑗<𝑖 𝑘>𝑖

𝑇𝑖 + ෍ 𝑇𝑗𝑖 − ෍ 𝑇𝑖𝑘 + ෍ 𝑟𝑖𝑗 × 𝐹𝑗𝑖 − ෍ 𝑟𝑖𝑘 × 𝐹𝑖𝑘 = 𝐼𝐺 𝑖 𝛼𝑖 + 𝜔𝑖 × 𝐼𝐺 𝑖 𝜔𝑖

𝑗<𝑖 𝑘>𝑖 𝑗<𝑖 𝑗<𝑖

6
 Ejemplo 2D DI
MANIVELA-BIELA-CORREDERA

 Ecuaciones de Newton-Euler para los eslabones 1, 2 y 3:

𝐹01𝑥 − 𝐹12𝑥 = 𝑚1 𝑥ሷ 𝐺1 1(x1,y1)

G1(xG1,yG1) G2(xG2,yG2) g
𝐹01𝑦 − 𝐹12𝑦 = 𝑚1 𝑔 + 𝑚1 𝑦ሷ 𝐺1 1
𝐹12𝑥 − 𝐹23𝑥 = 𝑚2 𝑥ሷ 𝐺2  2

T(t)? (t) 3
𝐹12𝑦 − 𝐹23𝑦 = 𝑚2 𝑔 + 𝑚2 𝑦ሷ 𝐺2 2(x2,0)
A(0,0) F(t)
𝐹23𝑥 = −𝐹 + 𝑚3 𝑥ሷ 2 0
𝐹03𝑦 + 𝐹23𝑦 = 𝑚3 𝑔
𝑇 𝑡 − 𝑥𝐺1 𝐹01𝑦 + 𝑦𝐺1 𝐹01𝑥 − 𝑥1 − 𝑥𝐺1 𝐹12𝑦 + 𝑦1 − 𝑦𝐺1 𝐹12𝑥 = 𝐼1 𝛽ሷ
𝑥1 − 𝑥𝐺2 𝐹12𝑦 − 𝑦1 − 𝑦𝐺2 𝐹12𝑥 − 𝑥2 − 𝑥𝐺2 𝐹23𝑦 −𝑦𝐺2 𝐹23𝑥 = 𝐼2 𝛾ሷ

 DI→ Sistema lineal de 8 ecuaciones y 8 incógnitas: las fuerzas entre eslabones y el par:
F01x, F01y, F03y, F12x, F12y, F23x, F23y, T(t)

7
 Ejemplo 2D DI
MANIVELA-BIELA-CORREDERA
 En forma matricial:
1 0 0 −1 0 0 0 0 𝐹01𝑥 𝑚1 𝑥ሷ 𝐺1
0 1 0 0 −1 0 0 0 𝐹01𝑦 𝑚1 𝑔 + 𝑚1 𝑦ሷ 𝐺1
0 0 0 1 0 −1 0 0 𝐹03𝑦 𝑚2 𝑥ሷ 𝐺2
0 0 0 0 1 0 −1 0 𝐹12𝑥 𝑚2 𝑔 + 𝑚2 𝑦ሷ 𝐺2
0 0 0 0 0 1 0 0 = −𝐹 + 𝑚3 𝑥ሷ 2
𝐹12𝑦
0 0 1 0 0 0 1 0 𝐹23𝑥 𝑚3 𝑔
𝑦𝐺1 −𝑥𝐺1 0 𝑦1 − 𝑦𝐺1 𝑥1 − 𝑥𝐺1 0 0 1 𝐹23𝑦 𝐼2 𝛾ሷ
0 0 0 − 𝑦1 − 𝑦𝐺2 𝑥1 − 𝑥𝐺2 −𝑦𝐺2 − 𝑥2 − 𝑥𝐺2 0 𝑇 𝐼1 𝛽ሷ
𝐴(𝑡) 𝑓(𝑡) 𝑏(𝑡)

 Solución (para cada paso de tiempo):

𝑓 𝑡 = 𝐴(𝑡)\b(𝑡)

 Se supone la cinemática resuelta, con q=(x1,y1,x2,,) y sus derivadas, pero también

necesitamos xG1, yG1; xG2, yG2 y sus derivadas.
► Opción 1: introducirlas como variables en q (o mejor, utilizar variables de punto de referencia en los CM)
► Opción 2: mediante la matriz de cambio de coordenadas locales-globales
8
 COORDENADAS LOCALES-GLOBALES
Matriz de cambio, caso 2D

 Objetivo: obtener las coordenadas globales (x,y) de cualquier P a partir de las locales (X,Y),
en función de las coordenadas globales de los puntos básicos 1(x1,y1) y 2(x2,y2).
𝑥 𝑥1 𝑋 𝑥1 cos 𝜗 − sin 𝜗 𝑋
𝑟Ԧ = 𝑟1 + 𝑅 ⇒ 𝑦 = 𝑦 + 𝑀𝜗 = 𝑦 + =
1 𝑌 1 sin 𝜗 cos 𝜗 𝑌
y 𝑋 𝑌
𝑥1 1 𝑥2 − 𝑥1 𝑦1 − 𝑦2 𝑥1 + 𝑥2 − 𝑥1 + 𝑦1 − 𝑦2
𝑋 𝐿 𝐿
= 𝑦 + 𝑦 −𝑦 𝑥2 − 𝑥1 = =
P 1 𝐿 2 1 𝑌 𝑋 𝑌
(x2,y2) 𝑦1 + 𝑦2 − 𝑦1 + 𝑥2 − 𝑥1
𝑟Ԧ 2 𝐿 𝐿
𝑅
 𝑋 𝑌 𝑋 𝑌 𝑋 𝑌 𝑋 𝑌 𝑥1
𝑥1 1 − + 𝑦1 + 𝑥2 − 𝑦2 1− − 𝑦1
1(x1,y1) 𝐿 𝐿 𝐿 𝐿 𝐿 𝐿 𝐿 𝐿
= = 𝑥2
𝑟1 𝑌 𝑋 𝑌 𝑋 𝑌 𝑋 𝑌 𝑋
−𝑥1 + 𝑦1 1 − + 𝑥2 + 𝑦2 − 1− 𝑦2
x 𝐿 𝐿 𝐿 𝐿 𝐿 𝐿 𝐿 𝐿
𝐶Ӗ 𝑞ത

𝒓 = 𝑪𝒒
 C es independiente del tiempo (sí cambia de un pto a otro del sólido), por tanto:
𝒓ሶ = 𝑪𝒒ሶ 𝒓ሷ = 𝑪𝒒ሷ
 Así podemos calcular las coordenadas de los CM y sus derivadas en todo momento
9
 Ejemplo 2D DI
MANIVELA-BIELA-CORREDERA
 Coordenadas globales de G1 y G2 y sus derivadas: “extender” las matrices C(1) y C(2)

1(x1,y1)
G1(xG1,yG1) 1 2 G2(xG2,yG2) 𝑋𝐺1 𝑌𝐺1 𝑋𝐺1 𝑌𝐺1 0
1− −
𝑥𝐺1 𝐿1 𝐿1 𝐿1 𝐿1 0
𝑦𝐺1 = 𝑌𝐺1 𝑋𝐺1 𝑌𝐺1 𝑋𝐺1 𝑥1
− 1− 𝑦1
𝐿1 𝐿1 𝐿1 𝐿1
2(x2,0) ⇒
𝑋𝐺2 𝑌𝐺2 𝑋𝐺2 𝑌𝐺2 𝑥1
A(0,0) 1− −
3
𝑥𝐺2 𝐿2 𝐿2 𝐿2 𝐿2 𝑦1
𝑦𝐺2 = 𝑌𝐺2 𝑋𝐺2 𝑌𝐺2 𝑋𝐺2 𝑥2
0 − 1− 0
0 𝐿2 𝐿2 𝐿2 𝐿2

𝑋𝐺1 𝑌𝐺1
− 0
𝐿1 𝐿1
𝑥𝐺1 𝑌𝐺1 𝑋𝐺1 𝒒𝐺 = 𝑪𝒒
0 𝑥1
𝑦𝐺1 𝐿1 𝐿1
⇒ 𝑥𝐺2 =
𝑋𝐺2 𝑌𝐺2 𝑋𝐺2
𝑦1 C independiente del tiempo  𝒒ሶ 𝐺 = 𝑪𝒒ሶ
1− 𝑥2
𝑦𝐺1 𝐿2 𝐿2 𝐿2
𝑌𝐺2 𝑋𝐺2 𝑌𝐺2 𝒒ሷ 𝐺 = 𝑪𝒒ሷ
− 1−
𝐿2 𝐿2 𝐿2
𝐶Ӗ

10
 Ejemplo 2D DI
MANIVELA-BIELA-CORREDERA
 si G1 y G2 en los ptos. medios de las barras: XG1/L1 = XG2/L2 = 1/2; YG1 = YG2 = 0 
1 0 0
1 0 1 0
𝐶=
2 1 0 1
0 1 0

 Como la cinemática está resuelta para x1, y1, x2 y  en función de  y sus derivadas,
se obtienen finalmente las aceleraciones de G1 y G2:

𝑥𝐺1 1 0 0 𝑥ሶ 𝐺1 1 0 0 𝑥ሷ 𝐺1 1 0 0
𝑥1 𝑥ሶ 1 𝑥ሷ 1
𝑦𝐺1 1 0 1 0 𝑦ሶ 1 0 1 0 𝑦ሷ 𝐺1 1 0 1 0
= 𝑦1 ; 𝐺1 = 𝑦ሶ 1 ; = 𝑦ሷ 1
𝑥𝐺2 2 1 0 1 𝑥ሶ 𝐺2 2 1 0 1 𝑥ሷ 𝐺2 2 1 0 1
𝑥2 𝑥ሶ 2 𝑥2ሷ
𝑦𝐺2 0 1 0 𝑦ሶ 𝐺2 0 1 0 𝑦ሷ 𝐺2 0 1 0
𝐶Ӗ 𝐶Ӗ 𝐶Ӗ

11
 Códiogo MatLab para resolver DI-Newton
MANIVELA-BIELA-CORREDERA
%Cálculo de velocidades %Dinámica inversa Newtoniana
A=Phiq(:,1:4); C=0.5*[1 0 0;
bv=-Phiq(:,5)*betap; 0 1 0;
vel=A\bv; 1 0 1;
0 1 0];
x1p=vel(1); qG=C*q(1:3);qppG=C*ace(1:3);
y1p=vel(2);
x2p=vel(3); Nw=[1 0 0 -1 0 0 0 0;
gammap=vel(4); 0 1 0 0 -1 0 0 0;
0 0 0 1 0 -1 0 0;
%Cálculo de aceleraciones 0 0 0 0 1 0 -1 0;
if abs(tan(beta))<1 0 0 0 0 0 1 0 0;
Phiqp=[2*x1p^2+2*y1p^2; 0 0 1 0 0 0 1 0;
2*(x2p-x1p)^2+2*y1p^2; qG(2) -qG(1) 0 y1-qG(2) x1-qG(1) 0 0 1;
-L2*gammap^2*sin(gamma); 0 0 0 -(y1-qG(4)) x1-qG(3) -qG(4) -(x2-qG(3)) 0];
L1*betap^2*sin(beta)];
else bf=[m1*qppG(1);
Phiqp=[2*x1p^2+2*y1p^2; m1*(g+qppG(2));
2*(x2p-x1p)^2+2*y1p^2; m2*qppG(3);
-L2*gammap^2*sin(gamma); m2*(g+qppG(4));
L1*betap^2*cos(beta)]; -F+m3*x2pp;
end m3*g;
I2*gammapp;
ba=-Phiq(:,5)*betapp-Phiqp; I1*betapp];
ace=A\ba;
f=Nw\bf;
x1pp=ace(1); %Construcción de los vectores posición, velocidad y
y1pp=ace(2); aceleración de la
x2pp=ace(3); %corredera y del par sobre la manivela de entrada
gammapp=ace(4); X2(i)=x2;X2p(i)=x2p;X2pp(i)=x2pp;Par(i)=f(8);T(i)=t;12
 ECUACIONES DE LAGRANGE
Conjunto LIBRE de coordenadas generalizadas
𝑑 𝜕𝑇 𝜕𝑇
Ppio. D’Alembert + Ppio. Trabajos Virtuales  ෍
𝑑𝑡 𝜕𝑞ሶ 𝑗
−
𝜕𝑞𝑗
− 𝑄𝑗 𝛿𝑞𝑗 = 0 (1)
𝑗

 {qj}: Conjunto de coord. generalizadas Ec. fund. de la DINÁMICA GENERALIZADA

 qj: desplazamientos virtuales: infinitesimales, atemporales y compatibles con ligaduras

 T: Energía cinética, 𝑇 = 𝑇(𝑞𝑗 , 𝑞ሶ 𝑗 ; 𝑡)
𝑎𝑝𝑙 𝜕𝑟Ԧ
 Qj: Fuerza generalizada asociada a la coord. qj;; 𝑄𝑗 = σ𝑖 𝐹Ԧ𝑖 · 𝑖
𝜕𝑞𝑗

 Si {qj} es LIBRE (j = 1…g): los qj son arbitrarios  en (1) todos los [ ] = 0:
𝑑 𝜕𝑇 𝜕𝑇
− − 𝑄𝑗 = 0; 𝑗 = 1, … , 𝑔 ECUACIONES DE LAGRANGE
𝑑𝑡 𝜕𝑞ሶ 𝑗 𝜕𝑞𝑗

𝑑 𝜕𝐿 𝜕𝐿
− − 𝑄𝑗𝑁𝐶 = 0; 𝑗 = 1, … , 𝑔 𝐿 =𝑇−𝑉
𝑑𝑡 𝜕𝑞ሶ 𝑗 𝜕𝑞𝑗

 QjNC: Fuerza generalizada no conservativa

 V: Energía potencial, V = 𝑉(𝑞𝑗 ; 𝑡)
13
 ECUACIONES DE LAGRANGE
Conjunto LIGADO de coordenadas. Multiplicadores de Lagrange

𝑑 𝜕𝑇 𝜕𝑇
෍ − − 𝑄𝑗 𝛿𝑞𝑗 = 0
𝑗
𝑑𝑡 𝜕𝑞ሶ 𝑗 𝜕𝑞𝑗 Ec. fundamental de la DINÁMICA GENERALIZADA

 Si {q1…qr,qr+1…qn} es LIGADO (q1…qr dependientes de los qr+1…qn):

► los qj no son independientes  NO todos los [ ] = 0
► Existen r = n – g ecuaciones de ligadura,  𝜙𝑘 (𝑞𝑗 ; 𝑡) = 0, 𝑘 = 1…𝑟

 la variación virtual de todas ellas es nula: 𝑛

𝜕𝜙𝑘
𝛿𝜙𝑘 = ෍ 𝛿𝑞 = 0, ∀𝑘 = 1 … 𝑟
𝜕𝑞𝑗 𝑗
𝑗
Para cualquier conjunto de valores 1…r también se cumplirá:
𝑟 𝑛
► 𝜕𝜙𝑘
෍ ෍ 𝜆𝑘 𝛿𝑞 = 0
Y, por tanto: 𝜕𝑞𝑗 𝑗
𝑘 𝑗
𝑟
𝑑 𝜕𝑇 𝜕𝑇 𝜕𝜙𝑘
෍ − − ෍ 𝜆𝑘 − 𝑄𝑗 𝛿𝑞𝑗 = 0;
𝑑𝑡 𝜕𝑞ሶ 𝑗 𝜕𝑞𝑗 𝜕𝑞𝑗
𝑗 𝑘=1

14
 ECS. DE LAGRANGE PARA CONJUNTOS LIGADOS
“Ajuste” de los Multiplicadores de Lagrange

𝑟
𝑑 𝜕𝑇 𝜕𝑇 𝜕𝜙𝑘
෍ − − ෍ 𝜆𝑘 − 𝑄𝑗 𝛿𝑞𝑗 = 0;
𝑑𝑡 𝜕𝑞ሶ 𝑗 𝜕𝑞𝑗 𝜕𝑞𝑗
𝑗 𝑘=1

Separando términos:
𝑟 𝑛

෍ 𝛿𝑞𝑗 + ෍ 𝛿𝑞𝑗 = 0
𝑗=1 𝑗=𝑟+1

𝑟 𝑛
 Los 1os r desplazamientos 𝛿𝑞𝑗 son dependientes; los restantes g, 𝛿𝑞𝑗 , independientes.
1 𝑟+1
 𝜆𝑘 arbitrarios   𝜆𝑘 𝑟
1 tales que los primeros r términos se anulan: = 0, 𝑗 = 1 … 𝑟
𝑛
 Los 𝑛 − 𝑟 = 𝑔 restantes multiplican a los 𝛿𝑞𝑗 INDEPENDIENTES: = 0, 𝑗 = 𝑟 + 1 … 𝑛
𝑟+1
 Por tanto, las ecs. de LAGRANGE para CONJUNTOS LIGADOS son: = 0, ∀𝑗 = 1 … 𝑛

𝑟
𝑑 𝜕𝑇 𝜕𝑇 𝜕𝜙𝑘
− − ෍ 𝜆𝑘 − 𝑄𝑗 = 0; 𝑗 = 1 … 𝑛
𝑑𝑡 𝜕 𝑞ሶ 𝑗 𝜕𝑞𝑗 𝜕𝑞𝑗
𝑘=1

15
 ECUACIONES DE LAGRANGE
Conjunto LIGADO de coord. Multiplicadores de Lagrange

 n ecuaciones diferenciales con n + r incógnitas (los r multiplicadores de Lagrange

“ajustados” son desconocidos)
 El sistema se resuelve incluyendo las propias r ecuaciones de restricción:
𝑟
𝑑 𝜕𝑇 𝜕𝑇 𝜕𝜙𝑘
− − ෍ 𝜆𝑘 − 𝑄𝑗 = 0; 𝑗 = 1 … 𝑛
𝑑𝑡 𝜕𝑞ሶ 𝑗 𝜕𝑞𝑗 𝜕𝑞𝑗
𝑘=1
𝜙𝑘 𝑞𝑗 = 0; 𝑘 = 1 … 𝑟

 Sistema de n ecs. DIFERENCIALES y r ecs. ALGEBRÁICAS (DAE). En notación compacta:

𝑑 𝜕𝑇 𝜕𝑇
− − 𝝓𝑇𝒒 𝝀 − 𝑸 = 0
൞𝑑𝑡 𝜕𝒒ሶ 𝜕𝒒
𝝓=0

 Objetivo Dinámica Máquinas/mecanismos: expresar T y Q en coordenadas naturales

16
 Multiplicadores - Fuerzas de ligadura

 Si se sustituyeran las r ligaduras, 𝜙𝑘 = 0, por las correspondientes

fuerzas, 𝑅𝑗 , el conjunto 𝑞𝑗 sería LIBRE y las ecs. de movimieto serían:

𝑑 𝜕𝑇 𝜕𝑇
− = 𝑄𝑗 + 𝑅𝑗
𝑑𝑡 𝜕𝑞ሶ 𝑗 𝜕𝑞𝑗

 Por lo tanto, se pueden obtener las fuerzas generalizadas de ligadura a

partir de los multiplicadores de Lagrange:
𝑟
𝜕𝜙𝑘
𝑅𝑗 = ෍ 𝜆𝑘
𝜕𝑞𝑗
𝑘

17
 Energía cinética en coordenadas naturales
MATRIZ DE CAMBIO COORDENADAS LOCALES - GLOBALES

 Las coordenadas globales (x,y) de cualquier P a partir de las locales (X,Y), dadas las
coordenadas globales de los PB 1(x1,y1) y 2(x2,y2) (Coord. Naturales)
𝑥 𝑥1 𝑋 𝑥1 cos 𝜗 − sin 𝜗 𝑋
𝑟Ԧ = 𝑟1 + 𝑅 ⇒ 𝑦 = 𝑦 + 𝑴𝝑 = 𝑦 + =
1 𝑌 1 sin 𝜗 cos 𝜗 𝑌
y 𝑋 𝑌
𝑥1 1 𝑥2 − 𝑥1 𝑦1 − 𝑦2 𝑥1 + 𝑥2 − 𝑥1 + 𝑦1 − 𝑦2
𝑋 𝐿 𝐿
= 𝑦 + 𝑦 −𝑦 𝑥2 − 𝑥1 = =
P 1 𝐿 2 1 𝑌 𝑋 𝑌
(x2,y2) 𝑦1 + 𝑦2 − 𝑦1 + 𝑥2 − 𝑥1
𝑟Ԧ 2 𝐿 𝐿
𝑅
 𝑋 𝑌 𝑋 𝑌 𝑋 𝑌 𝑋 𝑌 𝑥1
𝑥1 1 − + 𝑦1 + 𝑥2 − 𝑦2 1− − 𝑦1
1(x1,y1) 𝐿 𝐿 𝐿 𝐿 𝐿 𝐿 𝐿 𝐿
= = 𝑥2
𝑟1 𝑌 𝑋 𝑌 𝑋 𝑌 𝑋 𝑌 𝑋
−𝑥1 + 𝑦1 1 − + 𝑥2 + 𝑦2 − 1− 𝑦2
x 𝐿 𝐿 𝐿 𝐿 𝐿 𝐿 𝐿 𝐿
𝑪 𝒒

𝑟Ԧ = 𝑪𝒒
 La matriz C es independiente del tiempo (cambia de un pto a otro del sólido), por tanto:

𝑟Ԧሶ = 𝑪𝒒ሶ 𝑟Ԧሷ = 𝑪𝒒ሷ

18
 MATRIZ DE MASAS DEL SÓLIDO
Energía cinética en coordenadas naturales, caso plano.

1 1 1
 La energía cinética de todo el sólido es: 𝑇 = න 𝑟Ԧሶ 𝑇 𝑟Ԧሶ 𝑑𝑚 = න 𝒒ሶ 𝑇 𝑪𝑇 𝑪 𝒒𝑑𝑚
2 𝑆 2 𝑆
ሶ = 𝒒ሶ 𝑇 න 𝑪𝑇 𝑪 𝑑𝑚 𝒒ሶ
2 𝑆
 Con: 2𝑋 𝑋 2 + 𝑌 2 𝑋 𝑋2 + 𝑌2 𝑌
1− + 0 − −
𝐿 𝐿2 𝐿 𝐿2 𝐿
2𝑋 𝑋 2 + 𝑌 2 𝑌 𝑋 𝑋 + 𝑌2
2
0 1− + −
𝑇
𝑪 𝑪= 𝐿 𝐿2 𝐿 𝐿 𝐿2
𝑋 𝑋2 + 𝑌2 𝑌 𝑋 + 𝑌2
2
− 0
𝐿 𝐿2 𝐿 𝐿2
𝑌 𝑋 𝑋2 + 𝑌2 𝑋2 + 𝑌2
− − 0
𝐿 𝐿 𝐿2 𝐿2

𝑋 𝑚𝑋𝐺 𝑌 𝑚𝑌𝐺 𝑋2 + 𝑌2 𝐼1
 Empleando: න 𝑑𝑚 = 𝑚; න
𝐿
𝑑𝑚 =
𝐿
;න
𝐿
𝑑𝑚 =
𝐿
; න
𝐿2
𝑑𝑚 =
𝐿2
, se tiene:
𝑆 𝑆 𝑆 𝑆

2𝑚𝑋𝐺 𝐼1 𝑚𝑋𝐺 𝐼1 𝑚𝑌𝐺

𝑚− + 2 0 − 2 −
𝐿 𝐿 𝐿 𝐿 𝐿
2𝑚𝑋𝐺 𝐼1 𝑚𝑌𝐺 𝑚𝑋𝐺 𝐼1
1 0 𝑚− + 2 − 2 1 𝑇
𝑇 = 𝒒ሶ 𝑇
𝑚𝑋𝐺 𝐼1
𝐿
𝑚𝑌𝐺
𝐿 𝐿
𝐼1
𝐿 𝐿 𝒒ሶ
𝑇 = 𝒒ሶ 𝑀𝒒ሶ
2
𝐿
− 2
𝐿 𝐿 𝐿2
0 2
−
𝑚𝑌𝐺 𝑚𝑋𝐺 𝐼1
− 2 0
𝐼1 M: Matriz de masas del sólido
𝐿 𝐿 𝐿 𝐿2
𝑴 (independiente del tiempo) 19
 MATRIZ DE MASAS DEL MECANISMO

 La energía cinética de todo el mecanismo se obtiene sumando las energías cinéticas de

sus eslabones: 1 𝑇 1 𝑇
𝑇tot = ෍ 𝑇(𝑖) = ෍ 𝒒ሶ (𝑖) 𝑀(𝑖) 𝒒ሶ (𝑖) = 𝒒ሶ 𝑀𝑡𝑜𝑡 𝒒ሶ
2 2
𝑖 𝑖
donde q(i) es el vector de las 4 coordenadas del sólido i.
𝑖 𝑖 𝑖 𝑖
𝑀11 𝑀12 𝑀13 𝑀14 𝑥𝑗ሶ
𝑖 𝑖 𝑖 𝑖
𝑀21 𝑀22 𝑀23 𝑀24 𝑦𝑗ሶ
y 𝑥𝑗ሶ 𝑦ሶ𝑗 𝑥ሶ 𝑘 𝑦ሶ 𝑘 𝑖 𝑖 𝑖 𝑖
=
𝑀31 𝑀32 𝑀33 𝑀34 𝑥ሶ 𝑘
𝑖
𝑀41 𝑖
𝑀42 𝑖
𝑀43 𝑖
𝑀44 𝑦ሶ 𝑘

𝑥𝑗 𝑦𝑗 𝑥𝑘 𝑦𝑘
K (xk,yk) 0 ⋮ ⋮ 0 ⋮ ⋮ 0
𝑥𝑗 𝑖
⋯ 𝑀11 𝑖
𝑀12 𝑖
⋯ 𝑀13 𝑖
𝑀14 ⋯ 𝑥ሶ 1
J (xj,yj) 𝑦𝑗 𝑖 𝑖 𝑖 𝑖 𝑦ሶ 1
⋯ 𝑀21 𝑀22 ⋯ 𝑀23 𝑀24 ⋯
= 𝑥ሶ 1 𝑦ሶ 1 … 𝑥ሶ 𝑛 𝑦ሶ𝑛 0 ⋮ ⋮ 0 ⋮ ⋮ 0 ⋮
x 𝑥𝑘 𝑖 𝑖 𝑖 𝑖 𝑥ሶ 𝑛
⋯ 𝑀31 𝑀32 ⋯ 𝑀33 𝑀34 …
𝑦𝑘 𝑖
⋯ 𝑀41 𝑖
𝑀42 𝑖
… 𝑀43 𝑖
𝑀44 ⋯ 𝑦ሶ𝑛
0 ⋮ ⋮ 0 ⋮ ⋮ 0

 Matriz de masas del mecanismo, Mtot (nxn): suma de las matrices extendidas de masas,
M(i)ext (nxn), de sus eslabones.
20
 MATRIZ DE MASAS DEL MECANISMO
Ejemplo
x1 y1 x2 y2
x1 y1
5𝑚2 𝑚2 𝑚2

x1
𝑚1 𝑚1 0 −
0 0 𝑚1 12 12 4
3 6 0 0 0
𝑚1 𝑚1 5𝑚2 𝑚2 𝑚2

y1
0 0 3 0
3 6 𝑒𝑥𝑡 𝑚1 12 4 12 𝑒𝑥𝑡
𝑀 = → 𝑀(1) = 0 0 0 𝑀 2 = = 𝑀(2)
1 𝑚1 𝑚1 3 𝑚2 𝑚2 5𝑚2
y1 x1

0 0 0

y2 x2
6 3 0 0 0 0 12 4 12
𝑚1 𝑚1 0 0 0 0 𝑚2 𝑚2 5𝑚2
0 0 − 0
6 3 4 12 12

x2 y2
2 𝑚3 𝑚3

y2 x2
0 0
3 6 0 0 0 0
𝑚3 𝑚3 0 0 0 0
L/2 0 0 𝑚3
1 (x1,y1) 𝑀 3 = 𝑚3
3
𝑚3
6 𝑒𝑥𝑡
→ 𝑀(3) = 0 0 0
L 0 0 3
2 (x2,y2) 6 3 𝑚3
𝑚3 𝑚3 0 0 0
1 0 0 3
6 3
3
A B 5𝑚2 𝑚1 𝑚2 𝑚2
+ 0 −
12 3 12 4
5𝑚2 𝑚1 𝑚2 𝑚2
0 +
𝑀𝑡𝑜𝑡 = 12 3 4 12
𝑚2 𝑚2 5𝑚2 𝑚3
+ 0
12 4 12 3
𝑚2 𝑚2 5𝑚2 𝑚3
− 0 +
4 12 12 3
21
 FUNCIÓN “MsMx” genera la matriz de masas
extendida de un eslabón
%Extended Mass Matrix for a link having basic points P1,P2.
%P1,P2 natural numbers for basic points, enter 0 when it is a fixed point
%n=tot number of coordinates in the mechanism; m=mass of the link; L=distance P1-P2;
%X,Y: Centre of mass local coordinates; I=Moment of Inertia of the link about P1;

function M=MsMx(P1,P2,n,m,L,X,Y,I)
M=zeros(n);
if P1==P2
M(2*P1-1,2*P1-1)=m; M(2*P1,2*P1)=m;
else
if P1>0
M(2*P1-1,2*P1-1)=m-2*m*X/L+I/L^2; M(2*P1,2*P1)=M(2*P1-1,2*P1-1);
end
if P1>0 && P2>0 2𝑚𝑋𝐺 𝐼1 𝑚𝑋𝐺 𝐼1 𝑚𝑌𝐺
𝑚− + 2 0 − 2 −
M(2*P1-1,2*P2-1)=m*X/L-I/L^2; 𝐿 𝐿 𝐿 𝐿 𝐿
2𝑚𝑋𝐺 𝐼1 𝑚𝑌𝐺 𝑚𝑋𝐺 𝐼1
M(2*P1-1,2*P2)=-m*Y/L; 0 𝑚− + 2 − 2
M(2*P1,2*P2-1)=m*Y/L; 𝐿 𝐿 𝐿 𝐿 𝐿
𝑚𝑋𝐺 𝐼1 𝑚𝑌𝐺 𝐼1
M(2*P1,2*P2)=M(2*P1-1,2*P2-1); − 2 0
𝐿 𝐿 𝐿 𝐿2
M(2*P2-1,2*P1-1)=M(2*P1-1,2*P2-1); 𝑚𝑌𝐺 𝑚𝑋𝐺 𝐼1 𝐼1
M(2*P2-1,2*P1)=M(2*P1,2*P2-1); − − 2 0
𝐿 𝐿 𝐿 𝐿2
M(2*P2,2*P1-1)=M(2*P1-1,2*P2);
M(2*P2,2*P1)=M(2*P1,2*P2);
end
if P2>0
M(2*P2-1,2*P2-1)=I/L^2; M(2*P2,2*P2)=I/L^2;
end
end

22
 ECUACIONES DE LAGRANGE
Para conjuntos ligados. Forma compacta

 Utilizando la matriz de masas, los términos de las ecuaciones de Lagrange son:

𝑛
1 1
𝑇 = 𝒒ሶ 𝑇 𝑴𝒒ሶ = ෍ 𝑀𝑖𝑘 𝑞ሶ 𝑖 𝑞ሶ 𝑘
2 2
𝑖,𝑘
𝑛 𝑛 𝑛 𝑛
𝜕𝑇 1 1 1
= ෍ 𝑀𝑖𝑘 𝛿𝑖𝑗 𝑞ሶ 𝑘 + 𝑞ሶ 𝑖 𝛿𝑘𝑗 = ෍ 𝑀𝑗𝑘 𝑞ሶ 𝑘 + ෍ 𝑀𝑖𝑗 𝑞ሶ 𝑖 = ෍ 𝑀𝑖𝑗 𝑞ሶ 𝑖
𝜕𝑞ሶ 𝑗 2 2 2
𝑖,𝑘 𝑘 𝑖 𝑖
𝑛
𝑑 𝜕𝑇 𝜕𝑇
= ෍ 𝑀𝑖𝑗 𝑞ሷ 𝑖 ; =0
𝑑𝑡 𝜕𝑞ሶ 𝑗 𝜕𝑞𝑗
𝑖
 Y las ecuaciones de Lagrange para un sistema ligado resultan:
𝑛 𝑟
𝜕𝜙𝑘
෍ 𝑀𝑖𝑗 𝑞ሷ 𝑖 = ෍ 𝜆 + 𝑄𝑗 ; 𝑗 = 1, … , 𝑛
𝜕𝑞𝑗 𝑘
𝑖 𝑘
𝜙𝑘 = 0; 𝑘 = 1, … , 𝑟
 O en forma compacta:
𝑴𝒒ሷ = 𝝓𝑇𝒒 𝝀 + 𝑸 n ecs. Diferenciales
൝ DAE
𝝓=0 r ecs. algebráicas
23
 VECTOR FZAS. GENERALIZADAS
Un sólido. Caso plano. Fza. puntual

y F
 Fuerza puntual aplicada en el punto P
 La Qj asociada a qj es el coeficiente del P (x2,y2)
desplazamiento virtual qj en la expresión del 𝑟Ԧ 2
trabajo virtual, empleando vectores: 
1(x1,y1)
𝛿𝑊 = 𝑸𝑇 𝛿𝒒 𝑟1
x
 Si hemos elegido coord. cartesianas de 1 y 2:
𝑟Ԧ = 𝑪𝒒 ⇒ 𝛿 𝑟Ԧ = 𝑪𝛿𝒒 Entonces: 𝛿𝑊 = 𝑭𝑇 𝛿 𝑟Ԧ = 𝑭𝑇 𝑪𝛿𝒒
Y por comparación: 𝑸𝑇 = 𝑭𝑇 𝑪 ⇒ 𝑄 = 𝑪𝑇 𝑭
𝑋 𝑌
1− 𝐹𝑥 − 𝐹𝑦
𝐿 𝐿
𝑄1 1 − 𝑋 Τ𝐿 − 𝑌 Τ𝐿 𝑌 𝑋
𝐹𝑥 𝐹𝑥 + 1 − 𝐹
𝑄=
𝑄2
=
𝑌 Τ𝐿 1 − 𝑋 Τ𝐿 𝐿 𝐿 𝑦
𝑄3 𝑋 Τ𝐿 𝑌 Τ𝐿 𝐹𝑦 = 𝑋 𝑌
𝑄4 𝐹𝑥 + 𝐹𝑦
− 𝑌 Τ𝐿 𝑋 Τ𝐿 𝐿 𝐿
𝑌 𝑋
− 𝐹𝑥 + 𝐹𝑦
𝐿 𝐿
24
 VECTOR FZAS. GENERALIZADAS
Mecanismo. Caso plano. Fza. en PB. Momento

 Vector de fuerzas generalizadas del mecanismo: suma de los correspondientes

vectores Q extendidos
 Fuerza aplicada en un punto básico (PB): las fuerzas generalizadas asociadas a las
coordenadas del punto son las componentes correspondientes de la fuerza
 Momento puntual: aplicar un par de fuerzas puntuales equivalente
 Ejemplo: 2 𝐹Ԧ 𝑋 𝑌
 1− 𝐹𝑥 −𝐹
𝐿 𝐿 𝑦
L2/2 𝑄1 𝑌 𝑋
1 (x1,y1) 𝐹𝑥 + 1 − 𝐹
L2 𝑄𝑖 =
𝑄2
= 𝐿 𝐿 𝑦
2 (x2,y2) 𝑄3 𝑋 𝑌
L1 𝐹𝑥 + 𝐹𝑦
m2g L3 𝑄4 𝐿 𝐿
1 3 𝑌 𝑋
M m1g − 𝐹𝑥 + 𝐹𝑦
m3g 𝐿 𝐿
A B

𝑚2 𝑔 𝐹 𝑚2 𝑔 𝑀 𝐹
− sin 𝜗 − sin 𝛽 − sin 𝜗
𝑀 2 4 𝐿1 2
4 0 − sin 𝛽
0 𝑚2 𝑔 𝐿1 𝐹 𝑚1 + 𝑚2 𝑔 𝑀 𝐹
𝑚1 𝑔 − 0 cos 𝜗 − + cos 𝛽 + cos 𝜗
− 2 𝑀 2 2 𝐿1 2
𝑸= 2 + 𝑚2 𝑔 + 0 + cos 𝛽 + =
𝑚3 𝑔 𝐿1 𝐹 𝑚2 𝑔 𝐹
0 − 𝐹 cos 𝜗 + sin 𝜗 − + 𝐹 cos 𝜗 + sin 𝜗
4 − 0 2
0 𝑚2 𝑔 2 4 2
− 0 𝐹 𝑚2 + 𝑚3 𝑔 𝐹
2 − cos 𝜗 + 𝐹 sin 𝜗 − − cos 𝜗 + 𝐹 sin 𝜗
2 2 2 25
 FUNCIÓN “QF” genera el vector extendido de
Fzas. generalizadas de una F sobre un eslabón

%Extended Generalized force vector for a point force on a link having P1,P2
%P1,P2: natural numbers for basic points, enter 0 when it is a fixed point;
%n=Total number of coordinates in the mechanism; L=Distance between P1 and P2;
X,Y: Force actuation point in local coordinates; Fx,Fy: Global force components

function Q=QF(P1,P2,n,Fx,Fy,L,X,Y)
Q=zeros(n,1);
if P1==P2
Q(2*P1-1)=Fx; 𝑋 𝑌
1− 𝐹𝑥 − 𝐹𝑦
Q(2*P1)=Fy; 𝐿 𝐿
else 𝑌 𝑋
𝐹 + 1− 𝐹
if P1>0 𝐿 𝑥 𝐿 𝑦
Q(2*P1-1)=(1-X/L)*Fx-Y/L*Fy; 𝑋 𝑌
𝐹𝑥 + 𝐹𝑦
Q(2*P1)=Y/L*Fx+(1-X/L)*Fy; 𝐿 𝐿
end 𝑌 𝑋
− 𝐹𝑥 + 𝐹𝑦
if P2>0 𝐿 𝐿
Q(2*P2-1)=X/L*Fx+Y/L*Fy;
Q(2*P2)=-Y/L*Fx+X/L*Fy;
end
end

26
 FUNCIÓN “QW” genera el vector extendido de
Fzas. generalizadas del peso de un eslabón

%Extended Generalized force vector for the weight of a link having P1,P2
%P1,P2 natural numbers for basic points, enter 0 when it is a fixed point;
%n=tot number of coordinates in the mechanism;
%m=mass of the link;
%L=distance between P1 and P2;
%X,Y: Centre of mass local coordinates;

function W=QW(P1,P2,n,m,L,X,Y)
g=9.81;
if P1==P2
W=QF(P1,P2,n,0,-m*g);
else
W=QF(P1,P2,n,0,-m*g,L,X,Y);
end

27
 DISPOSITIVOS MECÁNICOS
Resorte (k,L0), Amortiguador (c), Actuador (a); entre A y B

𝑟Ԧ𝐵 − 𝑟Ԧ𝐴
 Resorte: fza. prop. a la elongación 𝑓Ԧ𝐴 = 𝑘 𝑑𝐴𝐵 − 𝐿0 𝑑𝐴𝐵
= −𝑓Ԧ𝐵

𝑟Ԧ𝐵 − 𝑟Ԧ𝐴 𝑟Ԧ𝐵 − 𝑟Ԧ𝐴

 Amortiguador: fza. prop. a la vel. relativa 𝑓Ԧ𝐴 = 𝑐 𝑣Ԧ𝐵 − 𝑣Ԧ𝐴 ·
𝑑𝐴𝐵 𝑑𝐴𝐵
= −𝑓Ԧ𝐵

𝑟Ԧ𝐵 − 𝑟Ԧ𝐴
 Actuador: fza. prescrita 𝑓Ԧ𝐴 = 𝑎 = −𝑓Ԧ𝐵
𝑑𝐴𝐵

 Es ventajoso añadir la coord. de la distancia s entre los puntos que une el dispositivo
(o el ángulo , en caso de dispositivo angular). Entonces se incluye directamente en
la posición del vector de fzas. generalizadas correspondiente a s (ó ):

Resorte: 𝑓𝑠 = −𝑘 𝑠 − 𝐿0 𝑓𝛽 = −𝑘 𝛽 − 𝛽0
Amortiguador: 𝑓𝑠 = −𝑐 𝑠ሶ 𝑓𝛽 = −𝑐 𝛽ሶ
Actuador: 𝑓𝑠 = 𝑎 𝑓𝛽 = 𝑎

28
 SIMULACIÓN CINEMÁTICA-DINÁMICA INVERSA
DIAGRAMA DE FLUJO

Ctes.:
Coord. aprox. Inicial: gdl: xr+1= xr+10;…; xn= xn0
Geom. q=[x1;…;xn] dt, tfin
x1=x10;…; xr=xr0 y sus derivadas
Masa

definir  y t=0,dt,…,tfin FIN

err=norm() gdl en func. de t
(for)
redefinir  y
Iteración t
err=norm()
(while)
[Link].
Iteración NO Cálculo Cálculo Almac.
err>E? M,Q
N-R velocs. acels. datos
Ec Mov
SI
definir A=q(:,1:r); b=-
resolver sol=A\b; OBTENCIÓN DE
q=q+[sol;zeros(g,1)] GRÁFICAS

29
 ECUACIÓN DINÁMICA INVERSA
Obtención del Par T sobre la manivela de entrada (máq 1gdl)

 El vector de coordenadas contiene el ángulo de la manivela: q=(x1,y1,…,yN,)T

 La ec. Dinámica es: 𝑄1
𝜆1
𝑇 𝑇 ⋮
𝑴𝒒ሷ = 𝑸 + 𝝓𝑞 𝝀⟶ 𝝓
ถ 𝑞 ⋮ = 𝑴
ณ 𝒒
ณሷ −
𝑄𝑛−1
𝑛×𝑟 𝜆 𝑟 𝑛×𝑛 𝑛×1
𝑛×1
𝑇
𝑟×1
𝑛×1
𝑛×1
 Las incógnitas son i y T. Las últimas fila y columna de M están compuestas de ceros,
ya que la coordenada angular no interviene. Por tanto, la última ecuación es:
𝜙1𝑛 𝜆1 + ⋯ + 𝜙𝑟𝑛 𝜆𝑟 = −𝑇 ⇒ 𝜙1𝑛 𝜆1 + ⋯ + 𝜙𝑟𝑛 𝜆𝑟 + 𝑇 = 0
 De modo que las incógnitas se pueden agrupar en un solo vector:
0 𝜆1 𝑄1
⋮ ⋮ ⋮
𝝓𝑞 𝑇 = 𝑴
ณ ณ 𝒒ሷ − ⟶ 𝑨𝒙 = 𝒃 ⇒ 𝒙=𝑨\𝒃
0 𝜆𝑟 𝑛×𝑛 𝑛×1
𝑄𝑛−1
1 𝑇 𝑛×1
0
𝑛×(𝑟+1) (𝑟+1)×1 𝑛×1
𝑛×1
30
 Códiogo MatLab para resolver DI-Lagrange
MANIVELA-BIELA-CORREDERA
%Cálculo de velocidades %Dinámica inversa Lagrangiana
A=Phiq(:,1:4); %Matriz de masas
bv=-Phiq(:,5)*betap; M1=MsMx(0,1,5,m1,L1,XG1,YG1,I1);
vel=A\bv; M2=MsMx(1,2,5,m2,L2,XG2,YG2,I2);
x1p=vel(1);y1p=vel(2); M3=MsMx(2,2,5,m3);
x2p=vel(3);y2p=vel(4); M=M1+M2+M3;
%Cálculo de aceleraciones %Vector de Fuerzas generalizadas
if abs(tan(beta))<1 QW1=QW(0,1,5,m1,L1,XG1,YG1);
Phiqp=[2*x1p^2+2*y1p^2; QW2=QW(1,2,5,m2,L2,XG2,YG2);
2*(x2p-x1p)^2+2*(y1p-y2p)^2; QW3=QW(2,2,5,m3);
0; QF3=QF(2,2,5,Fx,Fy);
L1*betap^2*sin(beta)]; Q=QW1+QW2+QW3+QF3;
else %Resolución de la ecuación
Phiqp=[2*x1p^2+2*y1p^2; L=[Phiq;0 0 0 0 1]';
2*(x2p-x1p)^2+2*(y1p-y2p)^2; bL=M*[ace;betapp]-Q;
0; LPar=L\bL;
L1*betap^2*cos(beta)]; %Construcción de los vectores
end posición, velocidad y aceleración de la
ba=-Phiq(:,5)*betapp-Phiqp; %corredera y del par sobre la
ace=A\ba; manivela de entrada
x1pp=ace(1);y1pp=ace(2); X2(i)=x2;X2p(i)=x2p;X2pp(i)=x2pp;
x2pp=ace(3);y2pp=ace(4); Y2(i)=y2;Y2p(i)=y2p;Y2pp(i)=y2pp;
Par(i)=LPar(5);T(i)=t;

31
 SIMULACIÓN DINÁMICA DIRECTA
DAE → ODE 1ª aproximación

𝑴𝒒ሷ = 𝝓𝑇𝒒 𝝀 + 𝑸 n ecs. Diferenciales

൝ DAE
𝝓=0 r ecs. algebráicas

 1ª aproximación: derivar dos veces las ecuaciones de restricción

𝝓 = 0 →ሶ 𝝓ሶ = 𝝓𝒒 𝒒ሶ = 𝟎 →ሶ 𝝓ሷ = 𝝓ሶ 𝒒 𝒒ሶ + 𝝓𝒒 𝒒ሷ = 𝟎, con lo que:

𝑴𝒒ሷ − 𝝓𝑇𝒒 𝝀 = 𝑸 𝑴𝒒ሷ − 𝝓𝑇𝒒 𝝀 = 𝑸

൝ →ቐ
𝝓=0 𝝓𝒒 𝒒ሷ = −𝝓ሶ 𝒒 𝒒ሶ

𝑴 −𝝓𝑇𝒒 𝒒ሷ 𝑸
=
𝝓𝒒 𝟎 𝝀 −𝝓ሶ 𝒒 𝒒ሶ
La solución es INESTABLE, ya que 𝝓ሷ = 𝟎 no implica 𝝓 = 𝟎
32
SIMULACIÓN DINÁMICA DIRECTA
DAE → ODE 1ª aproximación

EJEMPLO: PÉNDULO

33
 SIMULACIÓN DINÁMICA DIRECTA
DAE → ODE 1ª aproximación

EJEMPLO: PÉNDULO
%Dinámica Directa Péndulo 1ª aprox %1ªAproximación: derivar 2 veces
clear all A=[M,-phiq';phiq,zeros(r)];
close all b=[Q;-phiqpqp];
%CONSTANTES %Obtención vector de aceleraciones
L=1;%Longitud péndulo qppLambda=A\b;
m=1;%Masa extremo %n ecs dif 2ºO. en 2n ecs dif 1er O.
n=2;%número de coordenadas %Cambio variables q a variables X. definición
r=1;%número ecs. de restricción de las derivadas de X
%MATRIZ DE MASAS Xp(1:n)=qp;
M=zeros(2); Xp(n+1:2*n)=qppLambda(1:n);
M(1,1)=m; M(2,2)=m; %función handle de derivadas de X. Resolución
%VARIABLES SIMBÓLICAS deriv=odeFunction(Xp,[x(t) y(t) xp(t) yp(t)]);
syms t x(t) y(t) xp(t) yp(t) [T,X]=ode45(deriv,[0 5],[L;0;0;0]);
assume(t,'real') %SIMULACIÓN
q=[x(t);y(t)];assume(q,'real') for i=1:length(T)
qp=[xp(t);yp(t)];assume(qp,'real') pause(0.01)
%VECTOR DE FUERZAS GENERALIZADAS circunferencia(X(i,1),X(i,2),0.05*L,'b');
Q=[0;-9.81*m]; line([0,X(i,1)],[0,X(i,2)])
%MATRICES axis([-1.2*L 1.2*L -1.2*L 1.2*L])
phi=x(t)^2+y(t)^2-L^2; axis square
phiq=[2*x(t) 2*y(t)]; hold off
phiqpqp=2*(xp(t)^2+yp(t)^2); end
34
 SIMULACIÓN DINÁMICA DIRECTA
DAE → ODE Método de BAUMGARTE

 Sustituir las ecs. de restricción 𝝓 = 0 por una expresión en 𝝓,ሷ 𝝓ሶ y 𝝓

que hace que 𝝓 tienda “rápidamente” a cero: 𝝓ሷ + 2𝜉𝜔𝝓ሶ + 𝜔2 𝝓 = 0
Con =1 (amortiguamiento crítico) y =10

𝑴𝒒ሷ − 𝝓𝑇𝒒 𝝀 = 𝑸 𝑴𝒒ሷ − 𝝓𝑇𝒒 𝝀 = 𝑸

൝ ⟹ቐ
𝝓ሷ + 2𝜉𝜔𝝓ሶ + 𝜔 𝝓 = 0
2 𝝓𝒒 𝒒ሷ = −𝝓ሶ 𝒒 𝒒ሶ − 2𝜉𝜔𝝓𝒒 𝒒ሶ − 𝜔2 𝝓

𝑴 −𝝓𝑇𝒒 𝒒ሷ 𝑸
=
𝝓𝒒 𝟎 𝝀 −𝝓ሶ 𝒒 𝒒ሶ − 2𝜉𝜔𝝓𝒒 𝒒ሶ − 𝜔2 𝝓

35
 SIMULACIÓN DINÁMICA DIRECTA
DAE → ODE Método de BAUMGARTE

EJEMPLO: PÉNDULO
%Dinámica Directa Péndulo 1ª aprox %aproximación Baumgarte
clear all psi=1;omega=10;
close all A=[M,-phiq';phiq,zeros(r)];
%CONSTANTES b=[Q;-phipqp-2*psi*omega*phiq*qp-omega^2*phi];
L=1;%Longitud péndulo %Obtención vector de aceleraciones
m=1;%Masa extremo qppLambda=A\b;
n=2;%número de coordenadas %n ecs dif 2ºO. en 2n ecs dif 1er O.
r=1;%número ecs. de restricción %Cambio variables q a variables X. definición
%MATRIZ DE MASAS de las derivadas de X
M=zeros(2); Xp(1:n)=qp;
M(1,1)=m; M(2,2)=m; Xp(n+1:2*n)=qppLambda(1:n);
%VARIABLES SIMBÓLICAS %función handle de derivadas de X. Resolución
syms t x(t) y(t) xp(t) yp(t) deriv=odeFunction(Xp,[x(t) y(t) xp(t) yp(t)]);
assume(t,'real') [T,X]=ode45(deriv,[0 5],[L;0;0;0]);
q=[x(t);y(t)];assume(q,'real') %SIMULACIÓN
qp=[xp(t);yp(t)];assume(qp,'real') for i=1:length(T)
%VECTOR DE FUERZAS GENERALIZADAS pause(0.01)
Q=[0;-9.81*m]; circunferencia(X(i,1),X(i,2),0.05*L,'b');
%MATRICES line([0,X(i,1)],[0,X(i,2)])
phi=x(t)^2+y(t)^2-L^2; axis([-1.2*L 1.2*L -1.2*L 1.2*L])
phiq=[2*x(t) 2*y(t)]; axis square
phipqp=2*(xp(t)^2+yp(t)^2); hold off
end
36
 SIMULACIÓN DINÁMICA DIRECTA
Método de PENALIZADOR

 Sustituir los multiplicadores de Lagrange por fuerzas proporcionales al

incumplimiento de las ecs. de restricción, 𝝓 (en su forma “críticamente amortiguada”:
𝝓ሷ + 2𝜉𝜔𝝓ሶ + 𝜔2 𝝓; con =1 , =10)
 Se toma un alto valor para la cte. de proporcionalidad o PENALIZADOR, 106-107
𝑴𝒒ሷ = 𝑸 + 𝝓𝑇𝒒 𝝀 → 𝑴𝒒ሷ = 𝑸 + 𝝓𝑇𝒒 𝛼 𝝓ሷ + 2𝜉𝜔𝝓ሶ + 𝜔2 𝝓 ⇒
⇒ 𝑴𝒒ሷ = 𝑸 + 𝝓𝑇𝒒 𝛼 𝝓𝒒 𝒒ሷ + 𝝓ሶ 𝒒 𝒒ሶ + 2𝜉𝜔𝝓𝒒 𝒒ሶ + 𝜔2 𝝓 ⇒

𝑴 − 𝛼𝝓𝑇𝒒 𝝓𝒒 𝒒ሷ = 𝑸 + 𝛼𝝓𝑇𝒒 𝝓ሶ 𝒒 𝒒ሶ + 2𝜉𝜔𝝓𝒒 𝒒ሶ + 𝜔2 𝝓

 Equivale a sustituir un eslabón por un sistema MASA-RESORTE-AMORTIGUADOR

con valores muy altos de M,k,c
k

c
M
 SIMULACIÓN DINÁMICA DIRECTA
Método de PENALIZADOR

EJEMPLO: PÉNDULO
%Dinámica Directa Péndulo 1ª aprox %aproximación Penalizador
clear all psi=1;w=10;p=5e6;
close all A=M-p*(phiq'*phiq);
%CONSTANTES b=Q+p*phiq'*(phipqp+2*psi*w*phiq*qp+w^2*phi);
L=1;%Longitud péndulo %Obtención vector de aceleraciones
m=1;%Masa extremo qppLambda=A\b;
n=2;%número de coordenadas %Conversión: n ec dif 2º O en 2n ec dif 1er O
r=1;%número ecs. de restricción %Cambio variables q a variables X. definición
%MATRIZ DE MASAS de las derivadas de X
M=zeros(2); Xp(1:n)=qp;
M(1,1)=m; M(2,2)=m; Xp(n+1:2*n)=qppLambda(1:n);
%VARIABLES SIMBÓLICAS %función handle de derivadas de X. Resolución
syms t x(t) y(t) xp(t) yp(t) deriv=odeFunction(Xp,[x(t) y(t) xp(t) yp(t)]);
assume(t,'real') [T,X]=ode45(deriv,[0 5],[L;0;0;0]);
q=[x(t);y(t)];assume(q,'real') %SIMULACIÓN
qp=[xp(t);yp(t)];assume(qp,'real') for i=1:length(T)
%VECTOR DE FUERZAS GENERALIZADAS pause(0.01)
Q=[0;-9.81*m]; circunferencia(X(i,1),X(i,2),0.05*L,'b');
%MATRICES line([0,X(i,1)],[0,X(i,2)])
phi=x(t)^2+y(t)^2-L^2; axis([-1.2*L 1.2*L -1.2*L 1.2*L])
phiq=[2*x(t) 2*y(t)]; axis square
phipqp=2*(xp(t)^2+yp(t)^2); hold off
end
38
 SIMULACIÓN DINÁMICA DIRECTA
MANIVELA-BIELA-CORREDERA I
%Dinámica Directa Manivela-biela-corredera
close all
clear all

%Constantes
L1=0.25;%Longitud de la manivela
L2=0.75;%Longitud de la biela
Y2=0;%Distancia trayectoria corredera-pto. fijo manivela (descentrado)
m1=2;%Masa de la manivela
m2=5;%Masa de la biela
m3=1;%Masa de la corredera
xG2=0.3;%coord del CM de la biela según eje 1-2 (1: articulación manivela, 2: articulación corredera)
yG2=0;%Coord del CM de la biela según el eje a 90º, repecto al anterior
I1=0.5;%Momento de inercia de la manivela respecto a su extremo
I2=2.5;%Momento de inercia de la biela respecto a su extremo

n=3;%número de coordenadas
r=2;%número de ecuaciones de restricción (hay redundancia)

%Matriz de masas
M=zeros(n);
M(1,1)=I1/L1^2+m2+I2/L2^2-2*m2*xG2/L2; M(2,2)=M(1,1);M(3,3)=I2/L2^2+m3;
M(1,3)=m2*xG2/L2-I2/L2^2;M(2,3)=m2*yG2/L2;
M(3,1)=M(1,3);M(3,2)=M(2,3);

%Variables simbólicas
syms t x1(t) y1(t) x2(t) x1p(t) y1p(t) x2p(t)
q=[x1(t);y1(t);x2(t)];
qp=[x1p(t);y1p(t);x2p(t)];
assume(t,'real')
assume(q,'real')
assume(qp,'real')

Fuerza=50;%fuerza horizontal sobre la corredera(vale Fuerza(t))

Torque=10;%Par sobre la manivela(vale Torque(t))
Q=[-Torque*y1(t)/L1^2;Torque*x1(t)/L1^2;Fuerza]; 39
 SIMULACIÓN DINÁMICA DIRECTA
MANIVELA-BIELA-CORREDERA II
%Matrices
phi=[x1(t)^2+y1(t)^2-L1^2;
(x2(t)-x1(t))^2+(Y2-y1(t))^2-L2^2];
phiq=[2*x1(t) 2*y1(t) 0;
-2*(x2(t)-x1(t)) -2*(Y2-y1(t)) 2*(x2(t)-x1(t))];
phiqpqp=[2*(x1p(t)^2+y1p(t)^2);
2*((x1p(t)-x2p(t))^2+y1p(t)^2)];

%1ªAproximación: doble derivada de phi=0

A1=[M -phiq';phiq zeros(r)];
b1=[Q;-phiqpqp];
qppLambda=A1\b1;
X1p(1:n)=qp;
X1p(n+1:2*n)=qppLambda(1:n);
deriv1=odeFunction(X1p,[x1(t) y1(t) x2(t) x1p(t) y1p(t) x2p(t)]);
[T1,X1]=ode45(deriv1,[0 5],[L1;0;L1+L2;0;0;0]);

%Método Baumgarte
psi=1;omega=10;
bB=[Q;-phiqpqp-2*psi*omega*phiq*qp-omega^2*phi];
qppLambdaB=A1\bB;
XBp(1:n)=qp;
XBp(n+1:2*n)=qppLambdaB(1:n);
derivB=odeFunction(XBp,[x1(t) y1(t) x2(t) x1p(t) y1p(t) x2p(t)]);
[TB,XB]=ode45(derivB,[0 5],[L1;0;L1+L2;0;0;0]);

%Método penalizadores
alfa=5*10^6;%penalizador
AP=M-alfa*(phiq)'*phiq;
bP=Q+alfa*phiq'*(phiqpqp+2*psi*omega*phiq*qp+omega^2*phi);
qppP=AP\bP;
XPp(1:n)=qp;
XPp(n+1:2*n)=qppP;
derivP=odeFunction(XPp,[x1(t) y1(t) x2(t) x1p(t) y1p(t) x2p(t)]);
[TP,XP]=ode45(derivP,[0 5],[L1;0;L1+L2;0;0;0]);
40
 SIMULACIÓN DINÁMICA DIRECTA
MANIVELA-BIELA-CORREDERA III

%simulación 1ª Aproximación
for i=1:length(T1)
figure(1)
subplot(1,3,1)
if i==length(T1)
pause(0.01)
else
pause(T1(i+1)-T1(i))
end
circunferencia(0,0,0.05*L1,'r');hold on
circunferencia(0,0,0.1*L1,'r')
circunferencia(X1(i,1),X1(i,2),0.05*L1,'b');
circunferencia(X1(i,3),Y2,0.05*L1,'b')
line([0,X1(i,1),X1(i,3)],[0,X1(i,2),Y2])
axis([-1.2*L1 1.2*(L1+L2) -1.2*(2*L1+L2)/2
1.2*(2*L1+L2)/2])
title('1ªAPROXIMACIÓN')
axis square
hold off
end

41
 SIMULACIÓN DINÁMICA DIRECTA
MANIVELA-BIELA-CORREDERA IV

%simulación método Baumgarte

for i=1:length(TB)
subplot(1,3,2)
if i==length(TB)
pause(0.01)
else
pause(TB(i+1)-TB(i))
end
circunferencia(0,0,0.05*L1,'r');hold on
circunferencia(0,0,0.1*L1,'r')

circunferencia(XB(i,1),XB(i,2),0.05*L1,'b');
circunferencia(XB(i,3),Y2,0.05*L1,'b')
line([0,XB(i,1),XB(i,3)],[0,XB(i,2),Y2])
axis([-1.2*L1 1.2*(L1+L2) -1.2*(2*L1+L2)/2
1.2*(2*L1+L2)/2])
title('BAUMGARTE')
axis square
hold off
end

42
 SIMULACIÓN DINÁMICA DIRECTA
MANIVELA-BIELA-CORREDERA V

%simulación método Penalizador

for i=1:length(TP)
subplot(1,3,3)
if i==length(TP)
pause(0.01)
else
pause(TP(i+1)-TP(i))
end
circunferencia(0,0,0.05*L1,'r');hold on
circunferencia(0,0,0.1*L1,'r')
circunferencia(XP(i,1),XP(i,2),0.05*L1,'b');
circunferencia(XP(i,3),Y2,0.05*L1,'b')
line([0,XP(i,1),XP(i,3)],[0,XP(i,2),Y2])
axis([-1.2*L1 1.2*(L1+L2) -1.2*(2*L1+L2)/2
1.2*(2*L1+L2)/2])
title('PENALIZADOR')
axis square
hold off
end

43