Proyecto Eliptica 2

Proyecto Fin de Carrera
RENDIMIENTO ENERGÉTICO
DEL SISTEMA
CUERPO HUMANO/MECANISMO DE 4 BARRAS
Autor: Rubén Rosado Ávila
Tutor: José Ángel Acosta Rodríguez
INGENIERÍA INDUSTRIAL
DEPARTAMENTO DE INGENIERÍA DE SISTEMAS Y AUTOMÁTICA
ESCUELA TÉCNICA SUPERIOR DE INGENIERÍA
UNIVERSIDAD DE SEVILLA
Agradecimientos
A mi padre, mi madre y mi hermano,
por darme todo su apoyo y confianza.
A mi familia en general.
A los que no se encuentran entre nosotros,
siempre estaréis presentes.
A mis amigos y compañeros, con los que he compartido estos magníficos años en Sevilla.
A aquellos que me ayudaron a seguir hacia adelante.
A aquellos profesores que siempre me apoyaron.
A mi tutor José Ángel A. R.
por su orientación, seguimiento y supervisión del mismo.
1
2
ÍNDICE:
CAPÍTULO 1: INTRODUCCIÓN……………………...…………………………………………………………………PAG. 5
1.1 MOTIVACIÓN………………………………………………………………………………………………..PAG. 5
1.2 OBJETIVOS…………………………………………………………………………………………….……..PAG. 7
CAPÍTULO 2: TEORÍA DE MECANISMOS ……………………………………………………….……………….PAG. 9
2.1 CINEMÁTICA NUMÉRICA PLANA………………………………………………………..…………PAG. 9
2.1.1 INTRODUCCIÓN…………………..…………………………..………………………..…….PAG. 9
2.1.2 COORDENADAS CARTESIANAS O DE REFERENCIA….………………………..PAG. 9
2.1.3 RESTRICCIONES IMPUESTAS POR LOS PARES CINEMÁTICOS ………....PAG. 11
2.1.3.1 Par de revolución…………………………………………………………..PAG. 11
2.1.3.2 Par prismático………………………………………..……………………PAG. 12
2.1.4 PROBLEMA DE POSICIÓN…………………………………………………………….PAG. 13
2.1.5 SIMULACIÓN CINEMÁTICA…………………………………….…………………….PAG. 13
2.1.5.1 Problema de posición………………………………………………….PAG. 14
2.1.5.2 Problema de velocidad…………………………..…………………..PAG. 14
2.1.5.3 Problema de aceleración………………………..…………………..PAG. 15
2.1.5.4 Matriz jacobiana y matriz Bq…………………..…………………..PAG. 15
2.2 DINÁMICA NUMÉRICA PLANA……………………………………………………………………PAG. 17
2.2.1 INTRODUCCIÓN………………………………………………………………………….PAG. 17
2.2.2 DINÁMICA CON RESTRICCIONES…………………………………………………PAG. 17
2.2.2.1 Principio de los Trabajos Virtuales……………………………..PAG. 18
2.2.2.2 Método de los multiplicadores de Lagrange..…………….PAG. 18
2.3 PROBLEMA INVERSO………………………………………………………………………………...PAG. 23
CAPÍTULO 3: MODELADO Y DESCRIPCIÓN DEL SISTEMA BICICLETA ELÍPTICA……….………..PAG. 25
3.1 DEFINICIÓN DEL MECANISMO……………………………………………………………………PAG. 25
3.2 DATOS BICICLETA ELÍPTICA………..………………………………………………………………PAG. 27
3
3.3 CÁLCULO DE MASAS Y MOMENTOS DE INERCIAS…………………………………….PAG. 28
3.4 ESTUDIO CINEMÁTICO……………………………………………………………………………..PAG. 32
3.4.1 PROBLEMA DE POSICIÓN………………………………………………………….PAG. 35
3.4.2 PROBLEMA DE VELOCIDAD……………………………………………………….PAG. 41
3.4.3 PROBLEMA DE ACELERACIÓN……………………………………………………PAG. 41
3.5 ESTUDIO DINÁMICO…………………………………………………………………………………PAG. 44
CAPÍTULO 4: DESCRIPCIÓN DEL SOFTWARE Y USO……………..…………………………………………PAG. 47
4.1 FUNCIONES PROBLEMA POSICIÓN………………………………..…………………………PAG. 47
4.2 FUNCIONES PROBLEMA VELOCIDAD……………………………..…………………………PAG. 60
4.3 FUNCIONES PROBLEMA ACELERACIÓN…………………………………………………….PAG. 60
4.4 FUNCIONES CÁLCULO DINÁMICO…………………………………..………………………..PAG. 66
4.5 REPRESENTACIÓN Y SIMULACIÓN DEL MECANISMO……..………………………...PAG. 68
CAPÍTULO 5: ESTUDIO ECONÓMICO…………..…………………………………………………………………PAG. 75
5.1 INTRODUCCIÓN………….……………………………………………………………………………..PAG. 75
5.2 ELECCIÓN DEL GENERADOR……………………………………………………………………….PAG. 75
5.3 ELECCIÓN DEL CABLE..……………………………………………………………………………….PAG. 80
5.4 ELECCIÓN DE LA BATERÍA…….……………………………………………………………………PAG. 84
CAPÍTULO 6: RESULTADOS……………..…………………………………………………………………………….PAG. 87
CAPÍTULO 7: CONCLUSIÓN……………..…………………………………………………………………………….PAG. 91
CAPÍTULO 8: BIBLIOGRAFÍA………………..………………………………………………………………………..PAG. 93
CAPÍTULO 9: ANEXO………………..……………………………………………………………………………………PAG. 95
9.1 ANEXO A………….……………..………………………………………………………………………..PAG. 95
9.2 ANEXO B…………………………..…………………………………………………………………….PAG. 109
4
CAPÍTULO 1: INTRODUCCIÓN
Este capítulo es una breve introducción a la motivación para la realización de este proyecto fin
de carrera. También se comenta los objetivos finales que se pretenden alcanzar.
1.1 MOTIVACIÓN
Las “bicicletas elípticas“, o simplemente “elípticas”, son unas máquinas de uso indoor o de
gimnasio exclusivamente. Se trata de dos plataformas donde se apoyan los pies y dos barras
verticales donde agarrarse con las manos. Ambos movimientos, el de las manos y los pies van
acompasados de forma armoniosa. El esfuerzo es erguido, y a pesar de que se le llama
“bicicletas”, poco o nada tiene que ver con ellas, ya que es más afín con el esquí de fondo que
con el primero. En realidad es una especie de mezcla, sin llegar a ser ninguno de ellos, del esquí,
pedaleo, marcha y correr. El movimiento de los brazos recuerda mucho al esquí nórdico de
fondo, lo que lo diferencia especialmente del resto de máquinas aeróbicas, en las que es el tren
inferior el involucrado en el esfuerzo.
Las bicicletas elípticas son tras las bicicletas estáticas (spinning) y las cintas de corres, son las
máquinas más populares para el entrenamiento aeróbico en los gimnasios.
Al igual que ocurre con el ciclismo, el ejercicio en máquinas elípticas no produce impacto,
minimizando el riesgo de lesión. Son especialmente valoradas dentro del campo de la
rehabilitación de ambas extremidades. Tanto para lesiones deportivas, o personas mayores que
necesitan recuperar movilidad en alguna articulación.
El consumo calórico es superior en las máquinas elípticas que en las bicis estáticas debido a que
el movimiento implica mayor número de grupos musculares.
Los brazos y las piernas en su movimiento simultáneo hacen que la cantidad de calorías
quemadas en la unidad de tiempo sea sensiblemente mayor, siendo por ello una opción muy
adecuada para aquellos que pretenden perder peso.
Comparativa Bicicletas elípticas vs Bicis estáticas o spinning.
Diferencias entre bicicletas elípticas y las estáticas o de spinning:
 No disponen de sillín.
 El usuario está de pie y no sentado.
 El movimiento de los pedales es elíptico en vez de circular.
 Se ejercitan también los brazos y la parte superior del cuerpo, no sólo las piernas.
 No provoca molestias en la espalda ni en la zona de apoyo del sillín.
 Prácticamente no tiene semejanzas directas con ningún otro deporte, salvo el esquí de
fondo, que es muy relativa.
Similitudes entre las bicicletas elípticas y las estáticas:
 Movimiento de muy bajo impacto sobre las articulaciones.

 Alto consumo calórico.
5
 Sistema de resistencia variable para ajustar la intensidad del ejercicio al estado de forma
del usuario.
Así pues, las máquinas elípticas reúnen en un todo en uno, las ventajas de la cinta de correr
y las de la bici estática (spinning). Se realizan dos movimientos, uno circular y otro de
traslación, en el que se utilizan los músculos del tren superior e inferior.
Ventajas de las bicicletas elípticas:
Una bicicleta elíptica es en realidad una máquina de acondicionamiento físico con una
particularidad, el movimiento que realizas es la combinación de varios otros movimientos
deportivos a modo de síntesis.
Por tanto los movimientos como correr, marchar, pedalear, esquiar, escalar se reducen a
uno solo. Esta particularidad hace que sea un aparato de los llamados aeróbicos, que
entrena a la vez, todos estos gestos deportivos.
Como incorpora varios movimientos a la vez, trabaja más cantidad de grandes grupos
musculares simultáneamente, dando como resultado mayor cantidad de calorías quemadas,
por lo que disminuye tu grasa corporal y por consiguiente tu abdomen.
Desventajas de las bicicletas elípticas:
La desventaja de la bicicleta elíptica es que si eres deportista, como por ejemplo corredor,
entrenando en el elíptico no puedes transferir el gesto específico del deporte que practicas,
y pierdes eficacia en tu puesta en forma.
Sin embargo, suelen darse períodos de sobreentrenamiento entre los deportistas que los
obliga a suspender por un breve tiempo o reemplazar el plan programado, siendo la
bicicleta elíptica un excelente recurso para seguir trabajando el sistema aeróbico utilizando
un gesto deportivo distinto.
Cualquier rutina de ejercicios físicos debe incluir un plan aeróbico, si te gusta trabajar más
en la parte de fuerza muscular, la bicicleta elíptica es el perfecto complemento para tu
programa, debido a las diferentes posiciones que puedes adoptar para entrenar diversos
grupos musculares.
Por tanto, la bicicleta elíptica termina siendo una máquina versátil, completa y accesible
para que puedas optar, a la hora de iniciar ese postergado programa de acondicionamiento
físico en tu vida, que te permita perder peso, tonificando tu cuerpo.
A todo lo comentado anteriormente, se le añade además la motivación extra de aprovechar

la energía mecánica producida por la acción humana en la bicicleta elíptica, y convertirla en
energía eléctrica a través de un generador eléctrico.
Por ello, este proyecto se dedicará a modelar y simular la dinámica de la bicicleta elíptica a
distintas velocidades y para diferentes pesos de personas, para así calcular la energía
eléctrica que es capaz de obtener. Con ello se analizarán los resultados y así poder ver el
6
ahorro económico que todo ello supondría en un gimnasio considerando un número
variable de bicicletas elípticas funcionando unas determinadas horas al día.
1.2 OBJETIVOS
El objetivo principal que se pretende conseguir con este Proyecto Fin de Carrera es el
estudio del rendimiento energético del cuerpo humano en una bicicleta elíptica.
Para ello, se van a utilizar un modelo de bicicleta elíptica, exactamente la KETTLER VISO XS y
simular el movimiento de la misma, para así poder realizar un análisis tanto cinemático
como dinámico ante acciones externas. Esto es, ante las distintas variables que puede haber
en una bicicleta (peso de la persona, frecuencia de pedaleo y tiempo).
Los objetivos que desea cumplir este trabajo son:
1. Modelado de la bicicleta elíptica.
2. Simulación de la cinemática y dinámica de la bicicleta elíptica para las diferentes

combinaciones de variables, y posterior evaluación de los resultados.
3. Análisis del ahorro económico que supondría la colocación de un generador eléctrico a

nuestro mecanismo.
7
8
CAPÍTULO 2: TEORÍA DE MECANISMOS
En este capítulo se describen las herramientas que utilizaremos para poder realizar los cálculos
cinemáticos y dinámicos de nuestro mecanismos y que implementaremos en nuestro programa
de cálculo Matlab.
2.1 CINEMÁTICA NUMÉRICA PLANA.
En la siguiente sección se describirá las herramientas que se usan en el análisis computacional

de sistemas multicuerpo.
2.1.1 INTRODUCCIÓN
Se denominan sistemas multicuerpo a los conjuntos de sólidos rígidos y/o flexibles que se
encuentran unidos por pares cinemáticos, muelles, amortiguadores, actuadores, etc. Los
ejemplos de sistemas multicuerpo son numerosos: máquinas herramientas, robots, vehículos,
tanto terrestres como aeroespaciales, etc. En algunas aplicaciones de Biomecánica se estudia el
propio cuerpo humano o algunas partes de él como un sistema multicuerpo.
Aquí se describirán algunas herramientas utilizadas para resolver el problema cinemático en

mecanismo planos con métodos computacionales, lo que exige a la metodología que sea
completamente sistemática.
El problema cinemático consiste en conocer las velocidades y aceleraciones de todas las barras
de un mecanismo conocidas las velocidades y aceleraciones de algunas de ellas. En este sistema
se elegirán las coordenadas cartesianas o de referencia para definir el movimiento de un cuerpo
rígido.
A continuación se estudiarán las restricciones que imponen los pares más comunes que se
encuentran en mecanismos planos. Finalmente se mostrará el planteamiento general del
problema de simulación cinemática, analizando los problemas de posición, velocidad y
aceleración.
2.1.2 COORDENADAS CARTESIANAS O DE REFERENCIA
En la siguiente figura 2.1 se muestra un sólido rígido, i, sobre el que se han definido unos ejes
locales
FIGURA 2.1. COORDENADAS CARTESIANAS O DE REFERENCIA.
9
solidarios a él, (xi,yi). La posición del origen de este sistema local queda definida por el vector R i,
medido desde un sistema de referencia inercial (X,Y). La posición de un punto cualquiera del
espacio, P, será Rp vista desde el sistema inercial, al que a partir de este momento nos
referiremos como sistema global, y desde el sistema local. Ambos vectores de posición vienen
ligados por la siguiente ecuación:
Rp = Ri + Ai·
Donde Ai es la matriz de giro que expresa el vector local en coordenadas globales. Véase la
figura 2.2.
Y yi
θi
Rp
i
xi
Ri
FIGURA 2.2. SISTEMAS DE COORDENADAS DE REFERENCIA.
La expresión anterior se puede desarrollar en componentes quedando de la siguiente forma:
( )=( )+( )·( ) (2.1.1)
Se observa en la ecuación (2.1.1), que las coordenadas globales de un punto cualquiera, y en

particular de un punto del sólido, se pueden expresar en función de las coordenadas de
referencia del sólido, (Xi,Yi,θi) y de las coordenadas locales del punto, ( ).
La velocidad de un punto del sólido vendrá dada al derivar respecto al tiempo la ecuación
(2.1.1). Al realizar esta derivada hay que tener en cuenta que ahora P es un punto perteneciente
al sólido i, y, por tanto, es un vector constante. Así:
̇ p = ̇ i + ̇ i· = ̇ I + ̇ i· · (2.1.2)
10
Donde designa a la matriz que se obtiene al derivar la matriz Ai componente a componente
con respecto a θi:
=( ) (2.1.3)
La aceleración del punto P viene dada la derivar de nuevo con respecto al tiempo la ecuación
(2.1.2):
̈ p = ̈ i + ̈ i· · + ̇ i· ̇ · = = ̈ i + ̈ i· · - ̇ ·Ai· (2.1.4)
El segundo término del último miembro representa la aceleración tangencial relativa, mientras
que el último término es la aceleración normal relativa del punto P respecto al origen de
coordenadas local.
2.1.3 RESTRICCIONES IMPUESTAS POR LOS PARES CINEMÁTICOS MÁS COMUNES.
En esta sección se verá cómo se expresa matemáticamente, utilizando las coordenadas de

referencia, las restricciones impuestas por los pares más comunes. En nuestro mecanismo sólo
intervendrán pares de revolución.
2.1.3.1 Par de revolución
Un par de revolución es una conexión entre dos o más barras que permite un movimiento
rotacional entre ellas. Obsérvese la figura 2.3.
FIGURA 2.3.PAR DE REVOLUCIÓN
11
La condición que tiene que cumplir un par de revolución es la siguiente:
= (2.1.5)
Que si la desarrollamos nos queda:
( )+( )·( ) -( ) + ( )·( ) =( ) (2.1.6)
2.1.3.2 Par prismático
Un par de revolución es una conexión entre dos o más barras que permite un movimiento de
traslación o prismático entre ellas. Obsérvese la figura 2.4.
Las condiciones cinemáticas que tienen que cumplir son las siguientes:
(2.1.7)
( ) (2.1.8)
Pi
Ri
hi
Ri
Qi
Ri
Mj
i Ri
Ri
j
Ri
FIGURA 2.4.PAR PRISMÁTICO
Todo ello teniendo en cuenta, que nuestro vector de coordenadas de referencia es:
=( ) (2.1.9)
12
2.1.4 PROBLEMA DE POSICIÓN
Sea un mecanismo en el que se definen n coordenadas de referencia, entre las que se

establecen m ecuaciones de restricción impuestas por los pares cinemáticos. El número de
grados de libertad del mecanismo es n-m. El problema de posición consiste en dar valores a n-m
de esas coordenadas y resolver el sistema de m ecuaciones de restricción con m incógnitas que
resulta y que se suele escribir:
C(qd,qi)=[C1(qd,qi) C2(qd,qi )…Cm( qd,qi)]T=0 (2.1.10)
Donde qd es el vector de coordenadas dependientes, vector de incógnitas y qi es el vector de

coordenadas independientes, la posición de las barras de entrada, que están prescritas.
El anterior es un sistema de ecuaciones generalmente no lineal, que se puede resolver

numéricamente mediante el algoritmo de Newton- Raphson. Partiendo de un estimado inicial
qd0, se puede obtener una aproximación de la solución a partir de un estimado anterior
mediante la expresión:
qdj+1 =qdj-Cqd-1(qdj,qi)·C(qdj,qi) (2.1.11)
donde Cqd es la matriz jacobiana de las restricciones, dada por:
Cqd = = (2.1.12)
( )
Las iteraciones se detendrán cuando se considere que se ha alcanzado la convergencia. Para ello
se pueden considerar dos condiciones:
‖ ‖ (2.1.13)
‖ ‖
Las ecuaciones (2.1.13) imponen que la variación entre dos iteraciones sea pequeña,
estableciendo para ello una tolerancia Є1 y/o que se cumplan las ecuaciones de restricción con
una tolerancia Є2.
2.1.5 SIMULACIÓN CINEMÁTICA
En una simulación cinemática se pretende conocer las sucesivas posiciones que adoptan las
barras de un mecanismo cuando se mueve una serie de ellas, según una ley horaria definida. Las
ecuaciones de movilidad son las leyes horarias que establecen la variación temporal de las
coordenadas de referencia independientes, presentes en un número igual al de grados de
13
libertad. Conocida la evolución de las posiciones de todas las barras también se pueden analizar
las velocidades y aceleraciones.
2.1.5.1 Problema de posición
Las ecuaciones de movilidad se escriben:
qi – f(t) =0 (2.1.14)
siendo las coordenadas independientes funciones conocidas del tiempo. Al sistema de

ecuaciones de restricción (2.1.10) se le añaden las ecuaciones de movilidad resultando el
siguiente sistema de n ecuaciones con n incógnitas.
C(q,t)=[C1(q,t) C2(q,t )…Cn( q,t)]T=0 (2.1.15)
Donde ya no se distingue entre coordenadas dependientes e independientes. Todas se

consideran dependientes y el tiempo es la única variable independiente. Para un determinado
instante tk, el sistema de ecuaciones que resulta se puede resolver también por el método de
Newton-Raphson.
qj+1 =qj-Cq-1(qj,tk)·C(qj,tk) (2.1.16)
debiéndose plantear un nuevo problema de posición para el instante siguiente tk+1. Así,
dividiendo el intervalo de tiempo total [t1,t2] en subintervalos y resolviendo un problema de
posición en cada extremo de los subintervalos se puede conocer la posición de las barras del
mecanismo en el intervalo total.
No hay que olvidar que pueden darse casos de posiciones singulares en los mecanismos.
Las posiciones singulares del mecanismo son aquellas en las que no existe la inversa de la matriz
jacobiana, porque su determinante (también llamado jacobiano) es nulo. Los puntos muertos y
los puntos de bloqueo son algunas configuraciones singulares. Como se sabe, en los puntos
muertos el movimiento está indeterminado y esto es una consecuencia de que el jacobiano sea
nulo, como se verá enseguida. En lo que se refiere a la posición, que el jacobiano sea nulo sólo
impide resolver el problema de posición mediante el algoritmo de Newton- Raphson, aunque
puede resolverse por otros procedimientos.
2.1.5.2 Problema de velocidad
El sistema de ecuaciones (2.1.15) se deriva con respecto al tiempo y se obtienen las ecuaciones
de velocidad:
̇ (q,t) = 0 → · + =0 (2.1.17)
Ésta última ecuación se suele escribir:
Cq· ̇ + Ct = 0 (2.1.18)
14
Una vez resuelto el problema de posición se constituye un sistema de n ecuaciones con n
incógnitas lineal en las incógnitas ̇ . La solución a este problema es única si y sólo si el jacobiano
es no nulo, es decir, si el mecanismo no está en una configuración singular, como un punto de
bloqueo o un punto muerto. Por esa razón se dice que en un punto muerto el movimiento está
indeterminado, porque la solución al problema de velocidad no es única.
2.1.5.3 Problema de aceleración
Si se deriva con respecto al tiempo el sistema de ecuaciones de velocidad (2.1.18) se obtiene el

sistema de ecuaciones de aceleración.
(Cq· ̇ + Ct)=0; (2.1.19)
(Cq· ̇ + Ct)q ̇ + (Cq· ̇ + Ct)=0 (2.1.20)
(Cq· ̇ )q· ̇ + Ctq· ̇ + Cqt· ̇ + Cq· ̈ + Ctt=0 (2.1.21)
Cq· ̈ + (Cq· ̇ )q· ̇ + 2Cqt· ̇ + Ctt=0 (2.1.22)
Éste es un sistema de ecuaciones lineal que se puede reescribir:
Cq· ̈ = Qd con Qd =- (Cq· ̇ )q· ̇ - 2Cqt· ̇ - Ctt (2.1.23)
Si se ha resuelto el problema de posición y de velocidad todos los términos de Q d serán

conocidos y el problema de aceleración tendrá solución si existe la inversa de la matriz
jacobiana. El primer término de Qd se denominará
(Cq· ̇ )q = Bq (2.1.24)
2.1.5.4 Matriz jacobiana y matriz Bq
Ahora se muestra como es la estructura que tiene tanto la matriz jacobiana como la matriz B q
aplicada al sistema de ecuaciones de restricción de los pares cinemáticos.
Par de revolución
Matriz jacobiana.
La matriz jacobiana correspondiente a las ecuaciones de restricción impuesta por los pares de
rotación (2.1.6) es la siguiente:
(– )
Cq = [Cqi Cqj] = ( ) (2.1.25)
( ) –
15
Que reescribiendo queda:
Cq = [ · · ] (2.1.26)
Matriz Bq.
La matriz Bq (2.1.24) correspondiente a las ecuaciones de restricción impuesta por los pares de
rotación (2.1.6) sería:
Cq· ̇ = [Cqi Cqj]·( ̇̇ ) = ̇ I + ̇ i· · - ̇ j - ̇ j· · (2.1.27)
Bq = (Cq· ̇ ) = [0 -Ai· ·θi 0 Aj· ·θj] (2.1.28)
Par prismático
Matriz jacobiana.
La matriz jacobiana correspondiente a las ecuaciones de restricción impuesta por los prismáticas
(2.1.7) y (2.1.8) es la siguiente:
( ) (2.1.29)
( ) ( ) ( ) ( )
donde los términos que aparecen son los siguientes:
( ) ( ) ( ) ( )
Matriz Bq
La matriz Bq (2.1.24) correspondiente a las ecuaciones de restricción impuesta por los pares
prismáticos (2.1.7) y (2.1.8) sería:
̇ ̇
Cq· ̇ = [Cqi Cqj]·( ̇̇ )=(
( ) ̇ ̇ ( ) ̇ ( ) ( ) ̇ )
(2.1.30)
Bq = (Cq· ̇ )=( ̅ ̅ ̅ ̅ ̅ ̅ ) (2.1.31)
donde los términos que aparecen son los siguientes:
̅ ̅ ( ) ̇
̅ ̅ ( ) ̇
̅ ( ) ( ̇ ̇ ) ̇ ( ) ( ) ̇
16
y en esta ecuación se obtiene derivando
( ) ( ) ( ) · ( )
· ( )
dado que ( ) · ( ) , finalmente queda:
( ) ( )
El último término es:
̅ ( ) ( ) ̇ ̇ ( ) ( ) ̇ ( ) ( )
̇
2.2 DINÁMICA NUMÉRICA PLANA.
En esta sección se va a explicar la forma de plantear y resolver numéricamente problemas

dinámicos en mecanismos planos. Se utilizará también para su planteamiento la formulación en
coordenadas cartesianas de referencia.
2.2.1 INTRODUCCIÓN
El principio de D’Alembert establece que cualquier sólido está en equilibrio de fuerzas y

momentos si entre las fuerzas y momentos que actúan sobre él se incluyen las fuerzas y
momentos de inercia. El planteamiento general para sistemas multicuerpo con n sólidos (o
mecanismos con n barras) será:
̈ (2.2.1)
Donde M es la matriz de masa; el vector de fuerzas externas en el sólido i; el vector de las

reacciones en el sólido i; el vector de las fuerzas centrífugas del sólido i
2.2.2 DINÁMICA CON RESTRICCIONES
A continuación se mostrará cómo se pueden incluir en las ecuaciones de la dinámica el hecho de

que las coordenadas que describen el sistema estén relacionadas mediante ecuaciones de
restricción. Hay varios métodos, pero nosotros nos centraremos en el método de los
multiplicadores de Lagrange. Antes de pasar a describir el método se introducirá el Principio de
Trabajos Virtuales, que es una forma de establecer equilibrio de un sistema multicuerpo.
17
2.2.2.1 Principio de los Trabajos Virtuales.
El Principio de Trabajos Virtuales (PTV) establece que: “En cualquier sistema en equilibrio, el
trabajo desarrollado por las fuerzas aplicadas al sistema es nulo en cualquier desplazamiento
virtual, compatible con las ligaduras”. En un problema dinámico hay que incluir entre las fuerzas
aplicadas al sistema las fuerzas de inercia. En mecanismos, las ligaduras vienen impuestas por
los pares cinemáticos, así que un desplazamiento virtual puede ser cualquiera compatible con
las restricciones impuestas por los pares. EL PTV se expresa matemáticamente de la siguiente
forma:
∑ (2.2.2)
donde: es el trabajo virtual de las fuerzas de inercias
es el trabajo virtual de las fuerzas externas
es el trabajo virtual de las fuerzas de reacción
∑ ̈
∑ (2.2.3)
∑ =0
Donde: Qv es el vector que contiene a las fuerzas centrífugas
Qe es el vector de las fuerzas externas
Qc es el vector de reacciones
Siendo nulo el trabajo virtual de las reacciones en virtud del principio de acción y reacción.
Y por tanto la ecuación (2.2.2) queda:
̈ (2.2.4)
habiendo desaparecido las reacciones de las ecuaciones. Hay que observar que el paréntesis no
es nulo, dado que no todas las componentes de δq son independientes, ya que están
relacionadas por las ecuaciones de restricción impuestas por los pares.
2.2.2.2 Método de los multiplicadores de Lagrange.
La forma para resolver el problema dinámico con restricciones que vamos a utilizar es mediante
los multiplicadores de Lagrange. Este método permite además calcular las reacciones entre los
pares, que serán unas incógnitas más del problema y que están relacionadas con los
multiplicadores de Lagrange, como se verá a continuación. Así, al problema que permite obtener
tanto las aceleraciones como las reacciones se le denomina problema aumentado.
18
Fuerzas generalizadas
Una forma de plantear las ecuaciones de equilibrio es mediante el uso de las ecuaciones de
Lagrange:
( ̇
) (2.2.5)
En la Qi representa la fuerza generalizada asociada al desplazamiento virtual δqi y T la energía

cinética.
(2.2.6)
Considérese en primer lugar una fuerza F aplicada en un punto cualquiera P de un sólido (véase
figura 2.5). El trabajo virtual asociado a esta fuerza y a un desplazamiento virtual δq del sólido,
que provoque un desplazamiento virtual del punto de aplicación de la fuerza δRp, será:
(2.2.7)
FIGURA 2.5 FUERZA APLICADA EN UN PUNTO CUALQUIERA P DE UN SÓLIDO Y MOMENTO

APLICADO AL MISMO.
El desplazamiento virtual del punto de aplicación de la fuerza se obtiene diferenciando la

expresión que da su posición:
( ) ( ) (2.2.8)
De esta forma,
(2.2.9)
19
de la que se puede despejar QF
( ) ( )
( ) (2.2.10)
( ) ( )
Se puede observar que la última componente de QF representa el momento que da la fuerza

aplicada, F, respecto al origen de coordenadas local, con lo cuál, si este origen se hace coincidir
con el centro de gravedad del sólido, el vector QF contiene todas las acciones externas que
intervienen en el teorema de conservación de la cantidad de movimiento y del momento
cinético.
Si se aplica un momento M = M ⃗ , como también se muestra en la figura 2.5, el trabajo virtual

debido a dicho momento será
(2.2.11)
Conexión rígida
Si dos cuerpos se conectan rígidamente, se transmiten entre ellos una fuerza y un momento,
como se representa en la figura 2.6. En este caso, las ecuaciones de restricción expresadas en
nuestro sistema de referencia, serían:
(2.2.12)
y la matriz jacobiana Cq correspondiente al sistema de ecuaciones anterior (2.2.12) quedaría:
Cq=( ) (2.2.13)
Los multiplicadores de Lagrange están relacionados con las reacciones que se transmiten los
cuerpos a través de la conexión rígida, de la siguiente forma:
( ) ( ) (2.2.14)
Se puede definir una fuerza generalizada asociada a las reacciones, que se pueden considerar
unas fuerzas externas más, y esa definición se puede hacer tanto en el sólido i como en el j.
Según las ecuaciones (2.2.10) y (2.2.11) dicha fuerza generalizada será:
20
( ) ( ) (2.2.15)
FIGURA 2.6. DOS CUERPOS CONECTADOS RÍGIDAMENTE ENTRE SÍ.
El vector de fuerzas generalizadas debido a las reacciones se puede expresar en función de la

matriz jacobiana de las restricciones y los multiplicadores de Lagrange:
( ) ( ) (2.2.16)
( )
Par de revolución
El par de revolución permite el giro e impide los dos desplazamientos relativos entre los sólidos.
Éstos se transmiten entre sí una fuerza, de manera que el vector de multiplicadores de Lagrange
es:
λ=( ) (2.2.17)
La fuerza generalizada asociada a esta reacción es:
( ) ( ) ( ) (2.2.18)
21
Par prismático
En el par prismático se impide el giro relativo y el desplazamiento perpendicular al eje de

deslizamiento de un sólido respecto a otro. Al impedir estos grados de libertad aparecen un
momento y una reacción normal al contacto, que proporcionan los multiplicadores de Lagrange:
λ=( ) (2.2.19)
Las restricciones impuestas por los pares son
( ) ( ) (2.2.20)
y la matriz jacobiana
( ) (2.2.21)
( ) ( ) ( ) ( )
Se puede demostrar que entre el vector de reacciones y los multiplicadores de Lagrange existe
una relación idéntica a la de otros pares:
( ) (2.2.22)
Aplicación del método
A continuación se demostrará que los vectores Qc definidos en función de los multiplicadores de

Lagrange (reacciones) y la matriz jacobiana de las restricciones, son los mismos que aparecen en
la ecuación 2.2.1, es decir, que definiendo estos , cumplen:
̈ (ecuación 2.2.1)
En efecto, el PTV establece que
̈ (ecuación 2.2.4)
Por otro lado, al cumplirse las ecuaciones de restricción
(2.2.23)
Introduciendo esta ecuación en (2.2.4) resulta
̈ (2.2.24)
Esta ecuación se puede separar en componentes dependientes e independientes,
( ( ) ( ̈̈ ) ( ) ( ) ( ) ) (2.2.25)
22
que al desarrollarse queda:
( ̈ ̈ )+
( ̈ ̈ ) (2.2.26)
La matriz es cuadrada, de dimensiones mxm. Si además es no singular, como ocurre si el

mecanismo no adopta una configuración singular, tal como un punto de bloqueo o un punto
muerto, se puede invertir y es posible encontrar un conjunto de multiplicadores de Lagrange, un
vector λ, que haga que el segundo paréntesis se anule:
λ= ̈ ̈ (2.2.27)
Si se anula el segundo paréntesis de la ecuación (2.2.26), el primero también lo hará, ya que la

ecuación (2.2.25) debe de cumplirse para cualquier desplazamiento virtual. Esto implica que si
los multiplicadores de Lagrange se eligen como indica la ecuación (2.2.27), el paréntesis de la
ecuación (2.2.24) se anula y por tanto,
̈ (2.2.28)
Se ha demostrado entonces que si - la ecuación anterior reproduce la ecuación

(2.2.1). Es decir, que los multiplicadores de Lagrange que son solución de la ecuación (2.2.27)
deben de estar relacionados con las reacciones como ha quedado descrito, ya que así
reproducen la misma ecuación de equilibrio. En la formulación aumentada con multiplicadores
de Lagrange estas ecuaciones se completan con las relaciones cinemáticas de aceleración,
resultando finalmente el sistema de ecuaciones
( ) ( ̈) ( ) (2.2.29)
que constituye un sistema de ecuaciones diferenciales y algebraicas que permite calcular tanto
la respuesta del sistema como las reacciones, que están relacionadas con los multiplicadores de
Lagrange de la forma que se ha analizado anteriormente. La matriz de coeficientes del sistema
anterior es cuadrada de dimensiones n+m y simétrica, lo que simplifica su tratamiento desde el
punto de vista numérico.
Este sistema de ecuaciones que sirve tanto para resolver los problemas dinámicos directos, en
los que las incógnitas son las aceleraciones, como los dinámicos inversos, en los que las
incógnitas son las fuerzas. En nuestro caso se trata de un problema inverso ya que lo que no
conocemos son las fuerzas.
2.3 PROBLEMA INVERSO
En el problema inverso es conocido el estado de movimiento del mecanismo y las fuerzas

resistentes y se pretende calcular la fuerza o fuerzas motrices necesarias para conseguir el
movimiento deseado venciendo las fuerzas resistentes.
23
El movimiento queda definido por medio de la ley horaria qi, que define las coordenadas
independientes. Para conocer la evolución de las demás coordenadas es necesario resolver el
problema de posición y el problema de velocidad, definidos por
C(q,t)=0 → q(t)
̇ → ̇ (2.3.1)
La última ecuación de (2.3.1) se sustituye por la ecuación de la formulación aumentada. Cuando

en ésta se sustituye la posición y la velocidad de todas las barras calculadas anteriormente y las
fuerzas resistentes, se pueden despejar las aceleraciones y los multiplicadores de Lagrange.
( ̈) ( ) ( ) (2.3.2)
Los multiplicadores de Lagrange asociados a las ecuaciones de movilidad son las fuerzas
motrices que producen el movimiento prescrito. Considérese por ejemplo, el mecanismo de
cuatro barras de la figura 2.7 en el que se desea calcular el par motriz, M2 necesario para que la
barra 2 gire a velocidad constante. La ecuación de movilidad y la fila de la matriz Cq asociada
serán
| | | T
(2.3.3)
Si dicha ecuación de movilidad corresponde a la última de las restricciones, el vector de

multiplicadores de Lagrange será
( | | | | (2.3.4)
y al hacer aparece el par motriz, en la ecuación de equilibrio de momentos de la barra

2.
FIGURA 2.7: PAR MOTRIZ EN UN CUADRILÁTERO ARTICULADO
24
CAPÍTULO 3: MODELADO Y DESCRIPCIÓN DEL SISTEMA BICICLETA
ELÍPTICA
En este capítulo se describe el mecanismo de la bicicleta elíptica y el modelo matemático
cinemático y dinámico descrito en el capítulo anterior particularizado para nuestro caso.
3.1 DEFINICIÓN DEL MECANISMO
Aquí en esta sección se va a definir el mecanismo de la bicicleta elíptica y representar las

ecuaciones utilizadas para resolver su cinemática y dinámica.
Nuestro mecanismo se puede dividir en dos partes, la parte derecha y la parte izquierda de la
bicicleta. Ambas partes forman un mecanismo de 4 barras.
El mecanismo de la parte derecha sería el disco de inercia (barra 2), la barra 3 que es la que
llevaría la plataforma para colocar el pie derecho, y la barra 4 que es la barra donde
colocaríamos la mano derecha. Y el mecanismo de la parte izquierda sería el disco de inercia
(barra 2), la barra 5 que es la que llevaría la plataforma para colocar el pie izquierdo, ya la barra
6 que es la barra donde colocaríamos la mano izquierda. Para que haya una armonía en el
movimiento ambos lados están unidos al disco de inercia (barra 2) en puntos opuestos, para así
asegurar que cuando la plataforma de un pie esté abajo, la del otro pie quede arriba. Todo lo
descrito se observa en la figura en la figura 3.1.
Como estamos modelando en un plano, entonces nuestro mecanismo estará formado en total
por 6 barras, el cual tiene 7 pares de rotación, que son los siguientes:
O2 : unión entre la barra 2 y la barra fija. Se trata de un par de rotación fijo.
A: unión entre la barra 2 y la barra 3.
B: unión entre la barra 3 y la barra 4.
C: unión entre la barra 2 y la barra 5.
D: unión entre la barra 5 y la barra 6.
O4: unión entre la barra 4 y la barra fija. Se trata de un par de rotación fijo.
O4: unión entre la barra 6 y la barra fija. Se trata de un par de rotación fijo.
Como se observa en la siguiente figura 3.1, en el par de rotación fijo O4 se unen tres barras.
25
FIGURA 3.1: DESCRIPCIÓN Y NUMERACIÓN DE BICICLETA ELÍPTICA
Para poder realizar el análisis de nuestro mecanismo, también hay que definir un sistema de
coordenadas de referencia.
En la siguiente figura 3.2 se representa las coordenadas de referencia utilizadas:
FIGURA 3.2: COORDENADAS DE REFERENCIA DE BICICLETA ELÍPTICA
26
3.2 DATOS BICICLETA ELÍPTICA
A continuación se detallan los datos de la bicicleta elíptica utilizados para el estudio. Estos datos
lo hemos obtenidos del catálogo de una marca de bicicletas elípticas llamada KETTLER.
En la siguiente figura 3.3 se representa los vectores definidos para la realización del modelo
matemático
FIGURA 3.3: VECTORES DE BICICLETA ELÍPTICA
El modelo elegido es KETTLER VISO XS, cuyos datos son los siguientes:
Numeración de las barras
Barra 2: La que gira vueltas completas. Es un disco de inercia.
Barra 3: La siguiente a la barra 2. Es la flotante, donde va puesta la plataforma para colocar el

pie derecho.
Barra 5: La siguiente a la barra 2i. Es la flotante, donde va puesta la plataforma para colocar el
pie izquierdo.
27
Barra 4: La que va desde la barra 3 hasta el punto fijo O4, y de O4 hasta P
Barra 6: La que va desde la barra 5 hasta el punto fijo O4, y de O4 hasta Q.
Longitud de las barras
Barra 2. Radio 0,24 m
Barra 3=Barra 5: longitud 0,78 m
Barra 4=Barra 6: longitud 1,37 m
Sección de las barras
Barra 2: espesor 1 cm
Barra 3=Barra 5: sección rectangular (5x2,5 cm)
Barra 4=Barra 6: cilíndrica hueca (diámetro 3 cm)
Espesor del acero: 2mm
Otros datos de interés
Distancia entre puntos fijos (0,97 m) y ángulo de la recta de unión entre puntos fijos y la
horizontal (0,578 rad).
Ambos pedales están opuestos en la barra 2, es decir, cuando la barra 3 está abajo, la barra 5
está arriba.
Material de las barras: Acero (ρ=7850 kg/m3)
Con estos datos podemos calcular las masas y las inercias de todas las barras.
3.3 CÁLCULO DE MASAS Y MOMENTOS DE INERCIAS.
En esta sección vamos a proceder a calcular las masas y los momentos de inercias de todas las
barras del mecanismo.
Cálculo de Masas
m=ρ·V siendo ρ=densidad; V=volumen
28
Barra 2: Disco de inercia. Véase figura 3.4
FIGURA 3.4: DIMENSIONES BARRA 2
m=
Barra 3=Barra 5: Sección rectangular. Véase figura 3.5.
FIGURA 3.5: DIMENSIONES SECCIÓN BARRAS 3 Y 5
29
m=[ ] =
[ ]
Barra 4=Barra 6: Sección circular. Véase figura 3.6.
FIGURA 3.6: DIMENSIONES SECCIÓN BARRAS 4 Y 6
m=[( )]
[( )]
Cálculo del Momento de inercia de una distribución de masa continúa.
∫ siendo dm un elemento de masa situado a una distancia x del eje de rotación
 Momento de inercia de una barra
Vamos a calcular el momento de inercia de una varilla de masa y longitud L respecto de un

eje perpendicular a la varilla que pasa por el centro de masas, como muestra la figura 3.7.
FIGURA 3.7
30
La masa dm del elemento de longitud de la varilla comprendido entre x y x+dx es
El momento de inercia de una barra es:
Aplicando el teorema de Steiner, podemos calcular el momento de inercia de la varilla respecto

de un eje perpendicular a la misma que pasa por uno de sus extremos. Véase figura 3.8.
FIGURA 3.8
( ) siendo la inercia en el centro de la barra.
 Momento de inercia de un disco
Vamos a calcular el momento de inercia de un disco de masa M y radio R respecto de un

eje perpendicular al plano del disco y que pasa por su centro.
Tomamos un elemento de masa que dista x del eje de rotación. El elemento es un anillo
de radio x y de anchura dx. Si recortamos el anillo y lo extendemos, se convierte en un
rectángulo de longitud 2x y anchura dx, véase figura 3.9.
FIGURA 3.9
La masa es:
31
El momento de inercia del disco es:
Entonces particularizando en las barras del mecanismo, obtenemos:
Barra 2:
Barra 3=Barra 5:
donde es la masa de la persona modelada de forma que en cada lado irá la mitad de su
masa, y d es la distancia de donde está aplicada la masa al eje de centro de gravedad.
Dicha masa será una variable de nuestro estudio.
Barra 4=Barra 6:
3.4 ESTUDIO CINEMÁTICO
La resolución del problema de optimización comienza por el estudio cinemático. Es decir, el

cálculo de la posición, velocidad y aceleración de una serie de puntos característicos del
mecanismo.
Sea "q" el vector de coordenadas naturales de los puntos característicos del mecanismo y "C" el
conjunto de restricciones geométricas que debe cumplir el mecanismo durante su
funcionamiento.
Las restricciones geométricas del mecanismo se pueden escribir, de forma compacta, como:
C(q , t) = 0
La resolución del problema de posición consiste en determinar el vector "q" de coordenadas

naturales que cumpla con las condiciones de restricción, para una determinada posición del
eslabón de entrada.
Como las condiciones de restricción normalmente son no lineales, se utiliza en su resolución el

método de linealización iterativo de Newton-Raphson. Con este método, se obtiene el vector de
32
coordenadas naturales para una posición del mecanismo que cumple las restricciones
geométricas, para una determinada posición del eslabón de entrada.
Para iniciar el método de Newton-Raphson se debe partir de un vector de coordenadas

naturales aproximadas. Según sea ese vector inicial puede que el método no converja a una
solución aceptable; en cuyo caso, se debe probar con otro vector inicial de coordenadas
naturales, y así sucesivamente hasta conseguir converger a una solución que represente una
posición real del mecanismo. Un buen vector inicial suele ser el correspondiente a una posición
real del mecanismo y fácil de determinar, que sea próxima a la posición que se desea calcular.
Para obtener este vector inicial lo que hemos hecho es resolver la ecuación de lazo de cada
parte del mecanismo, es decir, primero para la parte derecha y luego para la parte izquierda.
La siguiente figura 3.10 define el mecanismo a estudiar. Los datos de partida son las longitudes
de las cuatro barras, las cuales son indeformables. El eslabón O2A es el eslabón motor, y en
función a cada posición (conocida) de este se va a calcular la posición de los eslabones AB y O4B.
FIGURA 3.10: MECANISMO 4 BARRAS SIMPLE
Para estudiar la posición del mecanismo, es necesario definir sus ecuaciones de cierre. Estas
ecuaciones definen la geometría de la máquina mediante un número de ecuaciones no
dependientes igual al número de incógnitas. Para ello, se hará la descomposición cinemática del
mecanismo.
A continuación, se representa el mecanismo de 4 barras que constituye cada lado (izquierdo y

derecho) de la bicicleta elíptica y se definen los ángulos de cada barra y los vectores que
utilizaremos para plantear las ecuaciones. Véase la figura 3.11.
33
FIGURA 3.11: MECANISMO 4 BARRAS BICICLETA ELÍPTICA
A continuación, hay que plantear las ecuaciones de cierre del mecanismo, que son las siguientes:
donde equivale a la mitad del vector de la figura 3.2, es decir =
Es un sistema de dos ecuaciones con dos incógnitas no lineal, en el que tenemos como datos las
longitudes de las barras, es decir, z1, z2, z3, z4; y el ángulo ϕ en cada instante. Las incógnitas son
Ω y ψ en función del ángulo del eslabón motor ϕ .
Por lo que si le damos un valor cualquiera a ϕ, en este caso le damos 0 rad, que es una posición
fácil de calcular, y resolvemos el sistema, obtenemos:
Ω=6,1706 rad
ψ=5,0295 rad
Resolviendo el mismo sistema de ecuaciones para el mecanismo de 4 barras del lado izquierdo,
dándole en este caso un valor a ϕ = π rad, obtenemos:
Ω’=6,2073 rad
ψ’=4,2762 rad
34
Una vez obtenido estos ángulos para una posición determinada ya podemos calcular el vector
de coordenadas naturales para dicha posición inicial “ ”. Para ello hay que tener en cuenta el
sistema de coordenadas de referencias elegido anteriormente. Véase figura 3.2.
Realizando las operaciones nos queda:
Coordenadas barra 1 (fija):
=0; =0; =0
Coordenadas barra 2:
=0; =0; =0
= ; = ; =
= ; = ; =
= ; = ; =
= ; = ; =
Por lo que queda:
(0 0 0 0 0 0 (z2+(z3/2))T θ3 (z2+z3+(z4/2))T θ4 (-z2+(z5/2))T θ5 (-z2+z5+(z6/2))T θ6)T
Con todo ello ya tenemos un valor de q0 bastante real para poner como vector inicial y así poder
iniciar el método de Newton-Raphson.
3.4.1 PROBLEMA DE POSICIÓN
En nuestro mecanismo las ecuaciones de restricción, que son sólo las correspondientes a los
pares de revolución son las siguientes:
También aclarar que las tres primeras ecuaciones de restricciones son referidas a la colocación
del sistema de referencia global, que en este caso lo vamos a colocar en el par O2, entonces:
(3.4.1)
; (3.4.2)
35
(3.4.3)
Para los pares de revolución de nuestro mecanismo se tiene que cumplir como se ha descrito en
la teoría la ecuación 2.1.6:
Sustituyendo la ecuación 2.1.6 en cada par de revolución, obtenemos:
Par de revolución O2
= siendo RO2 = Ri + Ai·
( )+( )·( ) -( ) - ( )·( ) =( )
como ( ) ( ) ; ( ) ( )
Por lo que queda:
(3.4.4)
(3.4.5)
Hacemos lo mismo con los siguientes pares de revolución:
Par de revolución A
= siendo RA = Ri + Ai·
( )+( )·( ) -( ) - ( )·( ) =( )
como ( ) ( ) ; ( ) ( )
Por lo que queda:
( ) (3.4.6)
( ) (3.4.7)
36
Par de revolución B
= siendo RB = Ri + Ai·
( )+( )·( ) -( ) - ( )·( ) =( )
como ( ) ( ) ; ( ) ( )
Por lo que queda:
( ) (3.4.8)
( ) (3.4.9)
( )+( )·( ) -( ) - ( )·( ) =( )
como ( ) ( ) ; ( ) ( )
Por lo que queda:
(3.4.10)
(3.4.11)
Par de revolución C
= siendo RC = Ri + Ai·
( )+( )·( ) -( ) - ( )·( ) =( )
como ( ) ( ) ; ( ) ( )
37
Por lo que queda:
( ) (3.4.12)
( ) (3.4.13)
Par de revolución D
= siendo RD = Ri + Ai·
( )+( )·( ) -( ) - ( )·( ) =( )
como ( ) ( ) ; ( ) ( )
Por lo que queda:
( ) (3.4.14)
( ) (3.4.15)
( )+( )·( ) -( ) - ( )·( ) =( )
como ( ) ( ) ; ( ) ( )
Por lo que queda:
(3.4.16)
(3.4.17)
Nos faltaría la ecuación de movilidad, para poder cerrar nuestro sistema de ecuaciones, que en
nuestro caso se la imponemos a la barra 2 la cual va a ser nuestra barra motriz:
38
(3.4.18)
Por lo que nuestro sistema de ecuaciones de restricción queda:
( )
( )
( )
( )
C(q,t)=
( )
( )
( )
( )
La resolución del problema de posición consiste en determinar el vector "q" de coordenadas

naturales que cumpla con las condiciones de restricción, para una determinada posición del
eslabón de entrada.
39
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
40
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0
𝑙
0 0 0 1 0 -𝑙 𝑠𝑒𝑛𝜃 -1 0 - 𝑠𝑒𝑛𝜃 0 0 0 0 0 0 0 0 0
𝑙
0 0 0 0 1 𝑙 𝑐𝑜𝑠𝜃 0 -1 𝑐𝑜𝑠𝜃 0 0 0 0 0 0 0 0 0
𝑙 𝑙
0 0 0 0 0 0 1 0 - 𝑠𝑒𝑛𝜃 -1 0 - 𝑠𝑒𝑛𝜃 0 0 0 0 0 0
𝑙 𝑙
0 0 0 0 0 0 0 1 𝑐𝑜𝑠𝜃 0 -1 𝑐𝑜𝑠𝜃 0 0 0 0 0 0
𝐶𝑞
-1 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 1 0 0 0 0 0 0 0
𝑙
0 0 0 1 0 𝑙 𝑠𝑒𝑛𝜃 0 0 0 0 0 0 -1 0 - 𝑠𝑒𝑛𝜃 0 0 0
𝑙
0 0 0 0 1 𝑙 𝑐𝑜𝑠𝜃 0 0 0 0 0 0 0 -1 𝑐𝑜𝑠𝜃 0 0 0
𝑙 𝑙
0 0 0 0 0 0 0 0 0 0 0 0 1 0 - 𝑠𝑒𝑛𝜃 -1 0 - 𝑠𝑒𝑛𝜃
𝑙 𝑙
0 0 0 0 0 0 0 0 0 0 0 0 0 1 𝑐𝑜𝑠𝜃 0 -1 𝑐𝑜𝑠𝜃
tiempo hay que aplicar el método de Newton-Raphson (véase ecuación 2.1.11):
-1 0 𝑑𝑥 𝑠𝑒𝑛𝜃 𝑑𝑦 𝑐𝑜𝑠𝜃 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0
0 -1 𝑑𝑥 𝑐𝑜𝑠𝜃 𝑑𝑦 𝑠𝑒𝑛𝜃 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
Para poder determinar el vector “q” de nuestro sistema de ecuaciones en cada instante de
Para ello, tenemos que calcular el jacobiano “Cq” de nuestro sistema C. Véase ecuación 2.1.12.
Ahora ya sí, se puede aplicar el método y obtener “q” para cualquier instante de tiempo. Todo ello se
ha realizado con el programa de cálculo numérico Matlab.
3.4.2 PROBLEMA DE VELOCIDAD
Una vez resuelto el problema de posición, el sistema de ecuaciones (2.1.15) se deriva con respecto al
tiempo y se obtienen las ecuaciones de velocidad (2.1.17) y que se suele escribir de forma más
abreviada como ecuación (2.1.18):
Entonces de esta ecuación (2.1.18) sólo nos queda por conocer la matriz Ct . Resultando:
Por lo que, el cálculo de la velocidad es:
3.4.3 PROBLEMA DE ACELERACIÓN
Una vez obtenidas las expresiones para el cálculo de la velocidad, volviendo a derivar las ecuaciones de
restricción respecto del tiempo se obtiene la ecuación (2.1.22):
A continuación calculamos las matrices que nos quedan por conocer. Resultando:
2Cqt· ̇ =0
41
̇
̇ ̇
̇ ̇ ̇ ̇
̇ ̇ ̇ ̇
̇ ̇ ̇ ̇
̇ ̇ ̇ ̇
̇= ̇
̇ ̇
̇ ̇ ̇
̇ ̇ ̇ ̇
̇ ̇ ̇ ̇
̇ ̇ ̇ ̇
̇ ̇ ̇ ̇
̇ ̇ ̇
̇ ̇ ̇
Este último sistema se deriva parcialmente respecto a las coordenadas de referencia y obtenemos la
matriz que nos falta para poder calcular la aceleración.
Que puede expresarse en forma compacta como:
42
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
43
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 -𝜃̇ 𝑙 𝑐𝑜𝑠𝜃 0 0 -𝜃̇ 𝑙 𝑐𝑜𝑠𝜃 0 0 0 0 0 0 0 0 0
0 0 0 0 0 -𝜃̇ 𝑙 𝑠𝑒𝑛𝜃 0 0 -𝜃̇ 𝑙 𝑠𝑒𝑛𝜃 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 -𝜃̇ 𝑙 𝑐𝑜𝑠𝜃 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 -𝜃̇ 𝑙 𝑐𝑜𝑠𝜃 0 0 0 0 0 0 0 0 0
𝑙
0 0 𝜃̇ 𝑑𝑥 𝑐𝑜𝑠𝜃 𝑑𝑦 𝑠𝑒𝑛𝜃 0 0 0 0 0 0 0 0 -𝜃̇ 𝑐𝑜𝑠𝜃 0 0 0 0 0 0
𝐵𝑞
𝑙
0 0 𝜃̇ 𝑑𝑥 𝑠𝑒𝑛𝜃 𝑑𝑦 𝑐𝑜𝑠𝜃 0 0 0 0 0 0 0 0 -𝜃̇ 𝑠𝑒𝑛𝜃 0 0 0 0 0 0
𝑙
0 0 0 0 0 𝜃̇ 𝑙 𝑐𝑜𝑠𝜃 0 0 0 0 0 0 0 0 -𝜃̇ 𝑐𝑜𝑠𝜃 0 0 0
𝑙
0 0 0 0 0 𝜃̇ 𝑙 𝑠𝑒𝑛𝜃 0 0 0 0 0 0 0 0 -𝜃̇ 𝑠𝑒𝑛𝜃 0 0 0
𝑙 𝑙
0 0 0 0 0 0 0 0 0 0 0 0 0 0 -𝜃̇ 𝑐𝑜𝑠𝜃 0 0 -𝜃̇ 𝑐𝑜𝑠𝜃
𝑙 𝑙
0 0 0 0 0 0 0 0 0 0 0 0 0 0 -𝜃̇ 𝑠𝑒𝑛𝜃 0 0 -𝜃̇ 𝑠𝑒𝑛𝜃
0 0 𝜃̇ 𝑑𝑥 𝑐𝑜𝑠𝜃 𝑑𝑦 𝑠𝑒𝑛𝜃 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 𝜃̇ 𝑑𝑥 𝑠𝑒𝑛𝜃 𝑑𝑦 𝑐𝑜𝑠𝜃 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Entonces una vez conocido las matrices, el cálculo de la aceleración sería: ̈ ̇
3.5 ESTUDIO DINÁMICO
Una vez resuelto el problema cinemático, se estudia el problema dinámico, es decir, el estudio de las
ecuaciones que relacionan las masas con la cinemática del mecanismo y con las fuerzas.
Debido a que el conjunto de coordenadas naturales no son independientes, se introducen los

multiplicadores de Lagrange en las ecuaciones que relacionan las masas con las fuerzas y las
aceleraciones. El sistema de ecuaciones (2.3.2) es el que tendremos que resolver.
A continuación se procede a representar las matrices necesarias para el cálculo dinámico.
Nuestra matriz M, al estar colocados los ejes locales en los centros de gravedad de cada barra, resulta
una matriz triangular, siendo:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 m2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 m2 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 I2 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 m3 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 m3 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 I3 0 0 0 0 0 0 0 0 0
M=
0 0 0 0 0 0 0 0 0 m4 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 m4 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 I4 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 m5 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 m5 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 I5 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 m6 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 m6 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 I6
44
El vector Q de fuerzas externas de nuestro sistema estará formado por la suma de la fuerza
gravitatoria, fuerza centrífuga y la fuerza de fricción.
Q= Qgrav + Qcentrigufa + Qfriccion
Qgrav=(0 mi·g mi·g· · )T ; Qcentrigufa= (mi· · · 0)T;

g=(0 -9.81)T;
=(0 0)T; ya que los ejes locales están colocados en los centros de gravedad de cada barra
0 0 0
0 0 0
0 0 0
0 0 0
- m2·9.81 0 0
0 0 c·w
0 0 0
-m3·9.81 0 0
0 0 0
Qgrav= Qcentrigufa= Qfriccion=
0 0 0
-m4·9.81 0 0
0 0 0
0 0 0
-m5·9.81 0 0
0 0 0
0 0 0
-m6·9.81 0 0
0 0 0
Como se observa la fuerza externa Qcentrifuga= 0 ya que hemos colocado los ejes locales de referencia en
su centro de gravedad.
45
También comentar que el término c·w que aparece en el vector Qfriccion representa la resistencia que va
a provocar el generador, el cual irá colocado en la barra 2.
Los multiplicadores de Lagrange asociados a las ecuaciones de movilidad son las fuerzas motrices que
producen el movimiento prescrito. Considerando nuestro mecanismo en el cual se desea calcular el
par motriz, M2 necesario para que la barra 2 gire a velocidad constante. La ecuación de movilidad y la
fila de la matriz Cq asociada serán
| | | | | T
Si dicha ecuación de movilidad corresponde a la última de las restricciones, el vector de

multiplicadores de Lagrange será:
λ= (0 0 0
)T;
y al hacer aparece el par motriz, en la ecuación de equilibrio de momentos de la barra 2.
46
CAPÍTULO 4: DESCRIPCIÓN DEL SOFTWARE Y USO.
En este capítulo se detalla el proceso de cálculo que nos llevará a la solución. Se hará mediante la
ayuda del software de cálculo numérico MATLAB.
La programación se ha hecho muy esquemática, para así facilitar su lectura. Para ello se ha
programado varias funciones que se explican a continuación. Todas las funciones serán llamadas por
un script en el cual se le introduce todos los datos de nuestro mecanismo y las condiciones iniciales, y
el script te calcula tanto la cinemática como la dinámica y finalmente representa y simula el
mecanismo.
El programa lo dividiremos en cuatro apartados, que serán:
1. Funciones problema posición

2. Funciones problema velocidad
3. Funciones problema aceleración
4. Funciones cálculo dinámico
Antes de explicar cada apartado, a modo de introducción, hay que decir que cada función estará
gobernada por un archivo llamado simulacionfinal.m.
Este archivo lanza la orden de llamar a todos las funciones.
4.1 FUNCIONES PROBLEMA POSICIÓN
Esta sección engloba las funciones necesarias para el cálculo del problema de posición de nuestro
mecanismo.
Función NewtonRaphsonJacob.m
Finalidad del archivo
El objetivo de esta función es calcular nuestro problema de posición mediante algoritmo de Newton-
Raphson explicado anteriormente para la cinemática numérica plana. Para ello se introducirá la
ecuación a resolver, los datos de entrada necesarios y las restricciones necesarias.
La ecuación a resolver es la vista anteriormente, la cual es la siguiente:
qj+1 =qj-Cq-1(qj,tk)·C(qj,tk) (ecuación 2.1.9)
Argumentos de entrada y de salida
Algunos de los datos necesarios para este archivo serán proporcionados de forma automática por
simulacionfinal.m. El resto de datos, (condiciones de contorno y ecuaciones a resolver) deberán
introducirse manualmente abriendo el fichero y escribiendo en él.
Para ello hay que darle unos datos de entrada, que son los siguientes:
47
Las restricciones de nuestro problema. Que llamamos fun, y que como veremos en simulacionfinal.m
será una función llamada Restricciones.m
El jacobiano de nuestro problema. Que llamamos jac, y que como veremos en simulacionfinal.m será
una función llamada Jacobiano.m
Y unas condiciones iniciales aproximada x0.
El archivo te devuelve nuestro sistema de coordenadas x exacto.
Algoritmo y codificación con Matlab
% Se leen los datos de entradas

function x = NewtonRaphsonJacob(fun,jac,x0)
%Se introduce los parámetros necesarios del NEWTONRAPHSON
AbsErr=10^(-3);
MaxIter=100;
Niter=0;
ok=0;
while ((ok==0)&&(Niter<MaxIter))
Niter=Niter+1;
CC = feval(fun,x0);
Cq = feval(jac,x0);
%Se evalúa la convergencia de la solución
%Se comprueba que el max(abs(CC)) sea inferior al error impuesto
if (max(abs(CC))<AbsErr)
ok=1;
end
%Se realiza la ITERACIÓN DE NEWTON-RAPHSON
x=x0-Cq\CC;
x0=x;
end
Función Restricciones.m
El objetivo de esta función es generar las ecuaciones de restricción de nuestro problema (C(q,t)=0),
para ello hay que darle como dato de entrada nuestro sistema de coordenadas q.
48
Los datos necesarios para el archivo es nuestro sistema coordenadas q.
Este archivo genera nuestro sistema de ecuaciones de nuestro mecanismo C.

function C = Restricciones(q)
%Se establecen las variables globales

global n nr np R P DR DP;
global t w;
%Se genera la matriz C correspondiente a las filas de las restricciones de

los pares cinemáticos
C = RestrPares(q);
%Se genera la última fila de la matriz correspondiente a la RESTRICCIÓN DE

MOVILIDAD
C(3*n,1) = q(6)+w*t;
Función RestrPares.m
El objetivo de esta función es calcular las ecuaciones de restricciones impuestas por los pares
cinemáticos.
Los datos necesarios para el archivo es nuestro sistema coordenadas q.
Este archivo genera las filas correspondientes a las restricciones de los pares cinemáticos de nuestro
sistema de ecuaciones de nuestro mecanismo C.

function C = RestrPares(qdep)

global IND DEP qind;
%Se genera las coordenadas independientes del mecanismo

for i=1:length(IND)
q(IND(i))=qind(i);
end
%Se genera las coordenadas dependientes del mecanismo
49
for i=1:length(DEP)
q(DEP(i))=qdep(i);
end
%Se genera las tres primeras filas de nuestro sistema de ecuaciones

correspondiente a la Barra fija
C(1,1)=q(1);
C(2,1)=q(2);
C(3,1)=q(3);
%Se genera las filas correspondientes a cada par cinemático de revolución de

nuestro mecanismo
for k=1:nr;
i=R(1,k);
j=R(2,k);
Ri=[q((i-1)*3+1) q((i-1)*3+2)]';
teti=q((i-1)*3+3);
Rj=[q((j-1)*3+1) q((j-1)*3+2)]';
tetj=q((j-1)*3+3);
ri=[DR(1,k) DR(2,k)]';
rj=[DR(3,k) DR(4,k)]';
Rr = Rotula(Ri,teti,ri,Rj,tetj,rj);
C(3+2*(k-1)+1,1)=Rr(1);
C(3+2*(k-1)+2,1)=Rr(2);
end
%%Se genera las filas correspondientes a cada par cinemático prismático de

nuestro mecanismo
for k=1:np;
i=P(1,k);
j=P(2,k);
Ri=[q((i-1)*3+1) q((i-1)*3+2)]';
teti=q((i-1)*3+3);
Rj=[q((j-1)*3+1) q((j-1)*3+2)]';
tetj=q((j-1)*3+3);
ri=[DP(1,k) DP(2,k)]';
hi=[DP(3,k) DP(4,k)]';
rj=[DP(5,k) DP(6,k)]';
Rp = Prism(Ri,teti,ri,hi,Rj,tetj,rj);
C(3+2*nr+2*(k-1)+1,1)=Rp(1);
C(3+2*nr+2*(k-1)+2,1)=Rp(2);
end
50
Función Rotula.m
El objetivo de esta función es devolver las ecuaciones de restricciones de los pares de rotación (2.1.5),
es decir, introduciendo Ri,θi,ri,Rj,θj,rj lo que hace internamente es lo siguiente:
Argumentos de entrada y salida
Esta función toma como argumentos de entrada la posición de la rótula respecto a cada barra
implicada, es decir, Ri,θi,ri,Rj,θj,rj .
Generaría como salida una matriz C correspondientes a las ecuaciones de los pares de revolución.

function C = Rotula(Ri,teti,ri,Rj,tetj,rj);
% Se genera la matriz A de cada barra implicada llamando a la función Rotmat

que explicaremos a continuación
Ai = Rotmat(teti);
Aj = Rotmat(tetj);
% Se genera la matriz C correspondiente a los pares de revolución

correspondiente a la ecuación (1.6)
C = Ri+Ai*ri-Rj-Aj*rj;
Función RotMat.m
El objetivo de esta función es devolverte la matriz Ai introduciendo los ángulos respectivos, donde Ai es
la matriz de giro que expresa el vector local en coordenadas globales.
Ai=( )
Esta función toma como argumentos de entrada la posición del ángulo de la barra implicada.
Generaría como salida la matriz Ai correspondiente a la barra i.

function A=Rotmat(x)
% Se genera la matriz A
51
A(1,1)=cos(x);
A(1,2)=-sin(x);
A(2,1)=sin(x);
A(2,2)=cos(x);
Función Prism.m
El objetivo de esta función es devolver las ecuaciones de restricciones de los pares prismáticos (2.1.8) y
(2.1.9), es decir, introduciendo Ri,θi,ri,hi,Rj,θj,rj lo que hace internamente es lo siguiente:
Esta función toma como argumentos de entrada la posición del par prismático respecto a cada barra
implicada, es decir, Ri,θi,ri,hi,Rj,θj,rj.
Generaría como salida una matriz C correspondientes a las ecuaciones de los pares prismáticos.

function C = Prism(Ri,teti,ri,hi,Rj,tetj,rj)

Ai = Rotmat(teti);
Aj = Rotmat(tetj);
% Se genera la fila de la matriz C correspondiente a los pares prismáticos

implicados. Ecuación (1.8)
C(1)=teti-tetj;
% Se genera la fila de la matriz C correspondiente a los pares prismáticos

implicados. Ecuación (1.9)
C(2)=(Ai*hi)'*(Ri+Ai*ri-Rj-Aj*rj);
En nuestro mecanismo no existen pares prismáticos. Por lo que no devolverá nada.
Función Jacobiano.m
El objetivo de esta función es calcular el jacobiano (Cq(q,t)) del sistema del ecuaciones de restricción
de nuestro problema, para ello hay que darle como dato nuestro sistema de coordenadas q.
52
Cq = =
( )
A esta función hay que darle como dato nuestro sistema de coordenadas q.
Se genera la matriz jacobiana Cq de nuestro mecanismo.
% Se leen los datos de entrada

function Cq = Jacobiano(q)

global n;
global t w;
%Se genera la matriz Cq a través de la llamada de otra función llamada

JacobPares.m
Cq = JacobPares(q);
%Se añade la última fila de la matriz correspondiente a la RESTRICCIÓN DE

MOVILIDAD de nuestro problema
Cq(3*n,6) = 1;
Función JacobPares.m
Esta función te calcula el jacobiano de las ecuaciones de restricciones impuestas por los pares
cinemáticos, tanto los de revolución como los prismáticos.
Hay que darle como dato nuestro sistema de coordenadas q.
Se genera la submatriz Cq correspondiente a las ecuaciones impuestas por los pares cinemáticos
existentes.

function Cqdep = JacobPares(qdep)

53
%Se genera las coordenadas independientes del mecanismo

for i=1:length(IND)
q(IND(i))=qind(i);
end
%Se genera las coordenadas dependientes del mecanismo

for i=1:length(DEP)
q(DEP(i))=qdep(i);
end
%JACOBIANO RESPECTO A TODAS LAS COORDENADAS
%Se calcula el Número de restricciones dependiendo del número de pares

cinemáticos que contenga nuestro mecanismo
m = 2*(nr+np)+3;
%Se genera la matriz Cq completa de ceros

Cq = zeros(m,3*n);
%Se rellena las casillas de la matriz Cq correspondientes a la Barra fija

Cq(1,1)=1;
Cq(2,2)=1;
Cq(3,3)=1;
%Se rellena las casillas de la matriz Cq correspondientes a los pares de

revolución. Rótulas
for k=1:nr;
%Localizamos las barras que actúan en el par nr correspondiente
i=R(1,k);
j=R(2,k);
%Localizamos la posición del centro de gravedad de cada barra según los

datos leídos
Ri=[q((i-1)*3+1) q((i-1)*3+2)]';
teti=q((i-1)*3+3);
Rj=[q((j-1)*3+1) q((j-1)*3+2)]';
tetj=q((j-1)*3+3);
%Localizamos la posición del par respecto a los ejes locales de cada

barra
ri=[DR(1,k) DR(2,k)]';
rj=[DR(3,k) DR(4,k)]';
% Generamos la matriz Cq referidas a las restricciones impuesta por los

pares de revolución
Cqr=Jac_rotula(Ri,teti,ri,Rj,tetj,rj);
% Rellenamos las casillas de la matriz Cq correspondientes a los pares de

revolución
Cq(3+2*(k-1)+1:3+2*(k-1)+2,3*(i-1)+1:3*(i-1)+3)=Cqr(1:2,1:3);
Cq(3+2*(k-1)+1:3+2*(k-1)+2,3*(j-1)+1:3*(j-1)+3)=Cqr(1:2,4:6);
end
%Se rellena las casillas de la matriz Cq correspondientes a los pares

Prismáticos
54
for k=1:np;
%Localizamos las barras que actúan en el par np correspondiente
i=P(1,k);
j=P(2,k);
%Localizamos la posición del centro de gravedad de cada barra según los

datos leídos
Ri=[q((i-1)*3+1) q((i-1)*3+2)]';
teti=q((i-1)*3+3);
Rj=[q((j-1)*3+1) q((j-1)*3+2)]';
tetj=q((j-1)*3+3);
%Localizamos la posición del par respecto a los ejes locales de cada barra
ri=[DP(1,k) DP(2,k)]';
hi=[DP(3,k) DP(4,k)]';
rj=[DP(5,k) DP(6,k)]';
% Generamos la matriz Cq referidas a las restricciones impuesta por los

pares prismáticos
Cqp=Jac_prism(Ri,teti,ri,hi,Rj,tetj,rj);
% Rellenamos las casillas de la matriz Cq correspondientes a los pares

prismáticos
Cq(3+2*nr+2*(k-1)+1:3+2*nr+2*(k-1)+2,3*(i-1)+1:3*(i-1)+3)=Cqp(1:2,1:3);
Cq(3+2*nr+2*(k-1)+1:3+2*nr+2*(k-1)+2,3*(j-1)+1:3*(j-1)+3)=Cqp(1:2,4:6);
end
%EXTRAE EL JACOBIANO RESPECTO A LAS COORDENADAS DEPENDIENTES
Cqdep = zeros(m,m);
for i=1:length(DEP)
Cqdep(:,i) = Cq(:,DEP(i));
end
Función Jac_rotula.m
Esta función te devuelve el jacobiano correspondiente a las ecuaciones de restricciones de los pares de
revolución.
Ejemplo:
Par de revolución A
= siendo RA = Ri + Ai·
55
( )+( )·( ) -( ) - ( )·( ) =( )
como ( ) ( ) ; ( ) ( )
Por lo que queda:
( )
( )
Entonces
( )
Cq= =( )=( )
( )
Esta función toma como argumentos de entrada la posición de la rótula respecto a cada barra
implicada, es decir, Ri,θi,ri,Rj,θj,rj .
Generaría como salida la matriz jacobiana Cq correspondiente a las ecuaciones de restricción de los
pares de revolución.

function Cq = Jac_rotula(Ri,teti,ri,Rj,tetj,rj)
% Se genera la matriz de ceros Cq correspondientes a las restricciones de

los pares de revolución
Cq=zeros(2,6);
% Se genera la matriz Ateta de cada barra implicada llamando a la función

RotTet que explicaremos a continuación
Ateti = RotTet(teti);
Atetj = RotTet(tetj);
% Vamos montando las casillas de la matriz que sabemos que va a resultar 1 o

-1
Cq(1:2,1:2)=eye(2);
Cq(1:2,4:5)=-eye(2);
% Montamos las restantes casillas que son las correspondientes a las

derivadas respecto de los ángulos.
Cq(1:2,3)=Ateti*ri;
Cq(1:2,6)=-Atetj*rj;
56
Función RotTet.m
El objetivo de esta función es devolverte la matriz introduciendo los ángulos respectivos, donde Ai
es la matriz de giro que expresa el vector local en coordenadas globales.
Ai=( )
y donde designa a la matriz que se obtiene al derivar la matriz Ai componente a componente con
respecto a θi:
=( )
Esta función toma como argumentos de entrada la posición del ángulo de la barra implicada.
Generaría como salida la matriz correspondiente a la barra i.

function Atet = RotTet(x)
% Se genera la matriz Atet
Atet(1,1)=-sin(x);
Atet(1,2)=-cos(x);
Atet(2,1)=cos(x);
Atet(2,2)=-sin(x);
Función Jac_prism.m
Esta función te devuelve el jacobiano correspondiente a las ecuaciones de restricciones de los pares
prismáticos. Véase ecuación 2.1.29.
Realiza prácticamente los mismos cálculos que la función Jac_rotula.m pero en este caso con las
restricciones de los pares prismáticos.
Esta función toma como argumentos de entrada la posición del par prismático respecto a cada barra
implicada, es decir, Ri,θi,ri,hi,Rj,θj,rj
57
Generaría como salida, la matriz jacobiana Cq correspondiente a las ecuaciones de restricción de los
pares de prismáticos existentes.

function Cq=Jac_prism(Ri,teti,ri,hi,Rj,tetj,rj);
% Se genera la matriz de ceros Cq correspondientes a las restricciones de

los pares prismáticos
Cq=zeros(2,6);
Ai = Rotmat(teti);
Aj = Rotmat(tetj);
% Se genera la matriz Ateta de cada barra implicada llamando a la función

RotTet
% Vamos montando las casillas de la matriz que sabemos que va a resultar 1 o

-1
Cq(1,3)=1;
Cq(1,6)=-1;
% Vamos montando las casillas de la matriz que sabemos que va a resultar

(Ai*hi)’ o –(Ai*hi)’
Cq(2,1:2)=(Ai*hi)';
Cq(2,4:5)=-(Ai*hi)';
% Montamos las restantes casillas que son las correspondientes a las

derivadas respecto de los ángulos.
Cq(2,3)=(Ateti*hi)'*(Ri+Ai*ri-Rj-Aj*rj)+(Ai*hi)'*(Ateti*ri);
Cq(2,6)=(Ai*hi)'*(-Atetj*rj);
Aclarar que en nuestro mecanismo no hay pares prismáticos, así que estas funciones no van a influir.
A continuación se representa un esquema para ver con mayor claridad la interacción de las distintas
funciones anteriormente comentadas para realizar el cálculo del problema de posición de nuestro
mecanismo mediante el algoritmo de Newton-Raphson. Véase figura 4.1
58
Entrada
qinicial NewtonRaphsonJacob.m
Devuelve (q)
Entrada
Entrada
Restricciones.m Jacobiano.m
Devuelve(C(q,t)) Devuelve Cq(q,t)
Llama a Llama a
RestrPares.m JacobPares.m
Devuelve C(q) Devuelve Cq(q)
Llama a
Llama a
Rotula.m Prism.m Jac_rotula.m Jac_prism.m
Devuelve C(q) Devuelve C(q) Devuelve Cq(q) Devuelve Cq(q)
(pares (pares (pares (pares

revolución) prismáticos) revolución) prismáticos)
Llama a
Llama a
RotTet.m Llama a
RotMat.m
Devuelve Devuelve
FIGURA 4.1. ESQUEMA PROGRAMA NEWTON-RAPHSON
59
4.2 FUNCIONES PROBLEMA VELOCIDAD
Este apartado engloba las funciones necesarias para el cálculo del problema de velocidad de nuestro
mecanismo.
Una vez resuelto el problema de posición, resolvemos el problema de velocidad, que para ello
tenemos que resolver la ecuación.
Como se observa solo tenemos que programar una función para que nos devuelva la matriz Ct.
La función que realiza el cálculo es la siguiente:
Función DtRestr.m
El objetivo de esta función es generar la matriz Ct = de nuestro mecanismo.
Esta función no toma nada como argumentos de entrada.
Generaría como salida la matriz Ct de nuestro mecanismo.

function Ct = DtRestr
% Se establece las variables generales

global t w;
% Se genera las filas de la matriz correspondientes a las restricciones de

los pares, las cuales al no depender del tiempo serán 0.
Ct = zeros(3*n,1);
%Se genera la última fila de la matriz, correspondiente a la RESTRICCIÓN DE

MOVILIDAD
Ct(3*n,1) = -w;
4.3 FUNCIONES PROBLEMA ACELERACIÓN
A continuación, este apartado engloba las funciones necesarias para el cálculo del problema de
aceleración de nuestro mecanismo.
60
Después del problema de velocidad hay que calcular el problema de aceleración, para ello tenemos
que resolver la ecuación:
̈ ̇
En este apartado nos queda por conocer la matriz Bq y la matriz Ctt.
Dichas matrices se calculan con las siguientes funciones.
Función MatrizBq.m
Te devuelve la matriz Bq de nuestro sistema, para ello hay que introducirle la posición (q) y la velocidad
( ̇ ).
El objetivo de esta función es calcular la matriz Bq del sistema del ecuaciones de restricción de nuestro
problema, siendo:
Cq = = Cq· ̇
( )
Bq = (Cq· ̇ )
Como argumentos, a esta función hay introducirle la posición q y la velocidad ( ̇ ).
Se genera la matriz Bq de nuestro mecanismo.

function Bq = MatrizBq(q,v)
% Se establece las variables globales

global n;
global t w;
% Se genera la matriz Bq, para ello se llama a la función MatBqPares que

explicaremos a continuación
Bq = MatBqPares(q,v);
%Por último se genera la última fila de la matriz correspondiente a la

RESTRICCIÓN DE MOVILIDAD
Bq(3*n,:) = zeros(1,3*n);
61
Función MatBqPares.m
Esta función te calcula la matriz Bq de las ecuaciones de restricciones impuestas por los pares
cinemáticos, tanto los de revolución como los prismáticos.
Hay que darle como dato la posición q y la velocidad ( ̇ ).
Se genera la submatriz Bq’ correspondiente a las ecuaciones impuestas por los pares cinemáticos
existentes.

function Bq=MatrizBq(q,v)
% Se establecen las variables globales

global n nr np R P DR DP w;
global t;
%Se calcula el Número de restricciones dependiendo del número de pares

cinemáticos que contenga nuestro mecanismo
m = 2*(nr+np)+3;
%Se genera la matriz Bq completa de ceros

Bq=zeros(m,3*n);
%Se rellena las casillas de la matriz Bq correspondientes a la Barra fija

Bq(1,1)=0;
Bq(2,2)=0;
Bq(3,3)=0;
%Se rellena las casillas de la matriz Bq correspondientes a los pares de

revolución. Rótulas
for k=1:nr;
i=R(1,k);
j=R(2,k);
%Localizamos la posición y velocidad del centro de gravedad de cada barra

según los datos leídos
Ri=[q((i-1)*3+1) q((i-1)*3+2)]';
Vi=[v((i-1)*3+1) v((i-1)*3+2)]';
teti=q((i-1)*3+3);
dteti=v((i-1)*3+3);
Rj=[q((j-1)*3+1) q((j-1)*3+2)]';
Vj=[v((j-1)*3+1) v((j-1)*3+2)]';
tetj=q((j-1)*3+3);
dtetj=v((j-1)*3+3);

barra
ri=[DR(1,k) DR(2,k)]';
rj=[DR(3,k) DR(4,k)]';
62
% Generamos la matriz Bq referidas a las restricciones impuesta por los
pares de revolución
Bqr=Bq_rotula(Ri,teti,Vi,dteti,ri,Rj,tetj,Vj,dtetj,rj);
% Rellenamos las casillas de la matriz Bq correspondientes a los pares de

revolución
Bq(3+2*(k-1)+1:3+2*(k-1)+2,3*(i-1)+1:3*(i-1)+3)=Bqr(1:2,1:3);
Bq(3+2*(k-1)+1:3+2*(k-1)+2,3*(j-1)+1:3*(j-1)+3)=Bqr(1:2,4:6);
end
%Se genera las casillas de la matriz Bq correspondientes a los pares

Prismáticos
for k=1:np;
i=P(1,k);
j=P(2,k);
%Localizamos la posición y velocidad del centro de gravedad de cada barra

según los datos leídos
Ri=[q((i-1)*3+1) q((i-1)*3+2)]';
Vi=[v((i-1)*3+1) v((i-1)*3+2)]';
teti=q((i-1)*3+3);
dteti=v((i-1)*3+3);
Rj=[q((j-1)*3+1) q((j-1)*3+2)]';
Vj=[v((j-1)*3+1) v((j-1)*3+2)]';
tetj=q((j-1)*3+3);
dtetj=v((j-1)*3+3);

barra
ri=[DP(1,k) DP(2,k)]';
hi=[DP(3,k) DP(4,k)]';
rj=[DP(5,k) DP(6,k)]';
% Generamos la matriz Bq referidas a las restricciones impuesta por los

pares prismáticos
Bqp=Bq_prism(Ri,teti,Vi,dteti,ri,hi,Rj,tetj,Vj,dtetj,rj);
% Rellenamos las casillas de la matriz Bq correspondientes a los pares

prismáticos
Bq(3+2*nr+2*(k-1)+1:3+2*nr+2*(k-1)+2,3*(i-1)+1:3*(i-1)+3)=Bqp(1:2,1:3);
Bq(3+2*nr+2*(k-1)+1:3+2*nr+2*(k-1)+2,3*(j-1)+1:3*(j-1)+3)=Bqp(1:2,4:6);
end
Función Bq_rotula.m
Esta función te devuelve los valores de la matriz Bq correspondiente a las ecuaciones de restricciones
de los pares de revolución. Véase ecuación (2.1.29)
63
Esta función toma como argumentos de entrada la posición y velocidad de la rótula respecto a cada
barra implicada, es decir, Ri,θi, ̇ , ̇ ,ri,Rj,θj, ̇ , ̇ ,rj .
Generaría como salida los valores de la matriz Bq correspondiente a las ecuaciones de restricción de
los pares de revolución.

function Bq=Bq_rotula(Ri,teti,Vi,dteti,ri,Rj,tetj,Vj,dtetj,rj)
% Se genera la matriz de ceros Bq correspondientes a las restricciones de

Bq=zeros(2,6);

Ai = Rotmat(teti);
Aj = Rotmat(tetj);
% Observando la ecuación (1.29) rellenamos las casillas

Bq(1:2,1:2)=zeros(2);
Bq(1:2,4:5)=zeros(2);
Bq(1:2,3)=-dteti*Ai*ri;
Bq(1:2,6)=dtetj*Aj*rj;
Función Bq_prism.m
Esta función te devuelve los valores de la matriz Bq correspondiente a las ecuaciones de restricciones
de los pares prismáticos. Véase ecuación (2.1.31).
Esta función toma como argumentos de entrada la posición y velocidad del par prismático respecto a
cada barra implicada, es decir, Ri,θi, ̇ , ̇ ,hi,ri,Rj,θj, ̇ , ̇ ,rj ,hj.
Generaría como salida los valores de la matriz Bq correspondiente a las ecuaciones de restricción de
los pares prismáticos

function Bq=Bq_prism(Ri,teti,Vi,dteti,ri,hi,Rj,tetj,Vj,dtetj,rj)
% Se genera la matriz de ceros Bq correspondientes a las restricciones de

Bq=zeros(2,6);
64
Ai = Rotmat(teti);
Aj = Rotmat(tetj);

% Observando la ecuación (1.35) rellenamos las casillas

Bq(2,3)=(-Ateti*hi)'*Vi+((-Ai*hi)'*(Ri+Ai*ri-Rj-
Aj*rj)+2*(Ateti*hi)'*(Ateti*ri)+(Ai*hi)'*(Ai*ri))*dteti+(-
Ateti*hi)'*Vj+((Ateti*hi)'*(-Atetj*rj)+(Ai*hi)'*(Aj*rj))*dtetj;
Bq(2,6)=(Ateti*hi)'*(-Atetj*rj)*dteti+(Ai*hi)'*(Aj*rj)*dtetj;
Función DttRestr.m
El objetivo de esta función es generar la matriz Ctt = de nuestro mecanismo.
Esta función no toma nada como argumentos de entrada.
Generaría como salida la matriz Ctt de nuestro mecanismo.

function Ctt = DttRestr
% Se establece las variables generales

global t w;
% Se genera las filas de la matriz correspondientes a las restricciones de

los pares, las cuales al no depender del tiempo serán 0.
Ctt = zeros(3*n,1);
%Se genera la última fila de la matriz, correspondiente a la RESTRICCIÓN DE

MOVILIDAD
Ctt(3*n,1) = 0;
4.4 FUNCIONES CÁLCULO DINÁMICO.
En este apartado se explica las funciones utilizadas para poder realizar el cálculo dinámico de nuestro
mecanismo.
65
( )·( ̈ )=( ̇ ̇ ̇
)
En el anterior sistema de ecuación tenemos que programar las funciones para que nos calcule la matriz
de masa M, y los vectores de fuerzas externas necesarios. Las funciones programadas son las
siguientes:
Función MatrizMasa.m
Esta función tiene como objetivo generar la matriz de Masa M del mecanismo.
siendo
M =( )
Como argumento de entrada hay que introducirle la posición (q).
Y genera como argumento de salida la matriz M del mecanismo.

function M = MatrizMasa(q)
% Se establecen los variables globales

global inercias;
% Se genera la matriz cero M, siendo n el número de barras

M = zeros(3*n,3*n);
for i = 1:n,
% Se localiza el ángulo de cada barra y generamos la variable teti
teti = q(3*(i-1)+3,1);
% Se genera la matriz Atet de cada barra
Atet = RotTet(teti);
% Se localiza los valores de las masas y generamos la variable mi
mi = inercias(1,i);
% Se localiza la posición del eje local de cada barra elegido, con
respecto al centro de gravedad de las barras, y generamos el vector rG
rG = inercias(3:4,i);
% Se calcula las inercias de cada barra respecto al eje elegido

I = inercias(2,i)+mi*(rG(1)^2+rG(2)^2); %Teorema de Steiner
% Vamos calculando los valores correspondientes de la matriz Qv y

colocándolos en sus respectivas posiciones
M(3*(i-1)+1:3*(i-1)+2,3*(i-1)+1:3*(i-1)+2) = mi*eye(2);
M(3*(i-1)+3,3*(i-1)+3) = I;
66
M(3*(i-1)+1:3*(i-1)+2,3*(i-1)+3) = mi*Atet*rG;
M(3*(i-1)+3,3*(i-1)+1:3*(i-1)+2) = mi*rG'*Atet';
end
Función FCentrifuga.m
Esta función te devuelve la fuerzas centrífugas Qcentrifugas, para ello hay que introducirle la posición (q).
Siendo Qcentrigufa= (mi· · · 0)T;
Y genera como argumento de salida la matriz Qcentrifugas del mecanismo.

function Qv = FCentrifuga(q)

global inercias;
% Se genera la matriz cero Qv, siendo n el número de barras

Qv = zeros(3*n,1);
for i = 1:n,
teti = q(3*(i-1)+3,1);
% Se genera la matriz A de cada barra llamando a la función Rotmat
A = Rotmat(teti);
mi = inercias(1,i);
respecto al centro de gravedad de las barras, y generamos la variable rG
% Vamos calculando los valores correspondientes de la matriz Qv y
Qv(3*(i-1)+1:3*(i-1)+2,1) = mi*teti^2*A*rG;
end
Función FGravedad.m
Esta función te devuelve la fuerzas gravitatorias Qgrav, para ello hay que introducirle la posición (q).
67
Siendo Qgrav=(0 mi·g mi·g· · )T y g=(0 -9.81)T;
Y genera como argumento de salida la matriz Qgrav del mecanismo.

function Qgrav = FGravedad(q)

global inercias g;
% Se genera la matriz cero Qgrav, siendo n el número de barras

Qgrav = zeros(3*n,1);
for i = 1:n,
teti = q(3*(i-1)+3,1);
% Se genera la matriz Atet de cada barra llamando a la función RotTet
mi = inercias(1,i);
respecto al centro de gravedad de las barras, y generamos la variable rG
% Vamos calculando los valores correspondientes de la matriz Qgrav y
Qgrav(3*(i-1)+1:3*(i-1)+2,1) = mi*g;
Qgrav(3*(i-1)+3,1) = mi*rG'*Atet'*g;
end
4.5 REPRESENTACIÓN Y SIMULACIÓN DEL MECANISMO
Una vez resuelto la programación tanto de la cinemática como de la dinámica de nuestro mecanismo,
en este apartado se mostrará las funciones y el script utilizado para poder representar los resultados
obtenidos y simular el mecanismo.
Función PosicionPunto.m
Esta función te devuelve la posición de cualquier punto (P) respecto de nuestro sistema referencia
global. El cálculo que realiza es el siguiente:
68
Como argumento de entrada hay que introducirle la posición (q) y el vector .
Y genera como argumento de salida el vector Rp, que sería la posición en ejes globales de cada punto
del mecanismo.

function Rp = PosicionPunto(qi,rp)
% Se localiza las coordenadas de cada barra y se genera el vector de

posición Ri y la matriz A
Ri = qi(1:2,1);
Ai = Rotmat(qi(3,1));
% Se genera el vector Rp que indica la posición de cada punto en ejes

globales
Rp = Ri + Ai*rp;
Simulacionfinal.m
Script del programa ejecutable.
El script es un archivo-m que contiene una serie de comandos que se ejecutarán al ejecutar dicho
archivo en MatLab. En dicho script se introducirán los datos de entradas necesarios. Una vez
introducidos el script irá llamando a las funciones anteriormente explicadas, e irá calculando todas
nuestras ecuaciones. Tanto el problema cinemático (posición, velocidad, aceleración) como el
problema dinámico. También hará una representación y simulación del mecanismo. Por último sacará
por pantalla los gráficos del par motor, potencia instantánea y trabajo motor respecto al tiempo.
global t w;
global inercias g;
%A continuación se introducen los datos de la posición las barras en una

%condición inicial calcula en la apartado 3.4 ESTUDIO CINEMÁTICO
z2=[0.24*cos(0) 0.24*sin(0)]';%vector barra 2
z3=[0.78*cos(353.549*pi/180) 0.78*sin(353.549*pi/180)]';%vector barra 3
z4=[1.3*cos((108.169)*pi/180) 1.3*sin((108.169)*pi/180)]';%vector barra 4
z1=[0.97*cos(0.578) 0.97*sin(0.578)]';%vector barra fija 1
n=6;%indica el número de barras del mecanismo
nr=7;%indica el número de pares de rotación del mecanismo
np=0;%indica el número de pares prismáticos del mecanismo
%Son vectores que nos indician las barras que intervienen en cada rótula
69
R(1:2,1)=[1 2]';%por ejemplo en la rótula 1 interviene la barra 1 y la 2
R(1:2,2)=[2 3]';
R(1:2,3)=[3 4]';
R(1:2,4)=[4 1]';
R(1:2,5)=[2 5]';
R(1:2,6)=[5 6]';
R(1:2,7)=[6 1]';
%Estos vectores indican la distancia de la rótula al eje local definido en

%cada barra
DR(1:2,1)=[0 0]';
DR(3:4,1)=[0 0]';
DR(1:2,2)=[norm(z2) 0]';
DR(3:4,2)=[-norm(z3)/2 0]';
DR(1:2,3)=[norm(z3)/2 0]';
DR(3:4,3)=[-norm(z4)/2 0]';
DR(1:2,4)=[0 0]';
DR(3:4,4)=[0.97*cos(0.578) 0.97*sin(0.578)]';
DR(1:2,5)=[-norm(z2) 0]';
DR(3:4,5)=[-norm(z5)/2 0]';
DR(1:2,6)=[norm(z5)/2 0]';
DR(3:4,6)=[-norm(z6)/2 0]';
%
DR(1:2,7)=[0 0]';
DR(3:4,7)=[0.97*cos(0.578) 0.97*sin(0.578)]';
%PROPIEDADES INERCIALES DE LOS SÓLIDOS
%(masa, momento de inercia respecto a G y posición de G en globales)
inercias(:,1) = [0 0 [0 0]]';
inercias(:,2) = [14.205 ((norm(z2)^2)*14.205/2+(0.046)) [0 0]]';
inercias(:,3) = [1.738 ((norm(z3)^2)*1.738/12 + 50*(0.2^2)) [0 0]]';
inercias(:,4) = [1.892 ((norm(z4)^2)*1.892/12) [0 0]]';
inercias(:,5) = [1.738 ((norm(z5)^2)*1.738/12 + 50*(0.2^2)) [0 0]]';
inercias(:,6) = [1.892 ((norm(z6)^2)*1.892/12) [0 0]]';
%Dirección y sentido de la gravedad
g = [0 -9.81]';
%%%%%%%%%%%%%%%%%%%%%%%%%
% SIMULACIÓN CINEMÁTICA %
%%%%%%%%%%%%%%%%%%%%%%%%%
IND = [];
DEP = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18]';
%Representamos los ángulos de cada barra en la condicional inicial elegida
fi20 = atan2(z2(2),z2(1))
psi20 = atan2(z4(2),z4(1))
70
teta20 = atan2(z3(2),z3(1))
psid20 = atan2(z6(2),z6(1))
tetad20 = atan2(z5(2),z5(1))
w = 95*(2*pi/60); %VELOCIDAD ANGULAR (95 rpm) (en este caso)

tspan = 0:0.05:(9*pi/w); %TIEMPO DE INTEGRACIÓN
%Estimacion inicial
q0=[0 0 0 0 0 fi20 (z2+(z3./2))' teta20 (z2+z3+(z4./2))' psi20 (-

z2+(z5./2))' tetad20 (-z2+z5+(z6./2))' psid20]';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BUCLE DE SOLUCIÓN DEL PROBLEMA DE POSICIÓN %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
q = zeros(3*n, length(tspan));
v = zeros(3*n, length(tspan));
a = zeros(3*n, length(tspan));
lam = zeros(3*n, length(tspan));
for i=1:length(tspan),
t = tspan(i);
%PROBLEMA DE POSICION
%Resolvemos el problema de posición mediante el algoritmo de
%Newton-Raphson empleando la siguiente función
q(:,i) = NewtonRaphsonJacob(@Restricciones,@Jacobiano,q0);
Cq = Jacobiano(q(:,i));%Calcula el jacobiano
Ct = DtRestr;
%CALCULO DE VELOCIDADES
v(:,i) = -Cq\Ct;
Bq = MatrizBq(q(:,i),v(:,i));
Ctt = DttRestr;
%CALCULO DE ACELERACIONES Y MULTIPLICADORES DE LAGRANGE

gam = -(Bq*v(:,i)+Ctt);
M = MatrizMasa(q(:,i));
Qv = FCentrifuga(q(:,i));
Qgrav = FGravedad(q(:,i));
Qfriccion = [0 0 0 0 0 0.95*w*1 0 0 0 0 0 0 0 0 0 0 0 0]';
AA = zeros(6*n,6*n);
bb = zeros(6*n,1);
AA(1:3*n,1:3*n) = M;
AA(1:3*n,3*n+1:6*n) = Cq';
AA(3*n+1:6*n,1:3*n) = Cq;
bb(1:3*n,1) = Qv+Qgrav+Qfriccion;
bb(3*n+1:6*n,1) = gam;
xx = AA\bb;
a(:,i) = xx(1:3*n,1);
71
lam(:,i) = xx(3*n+1:6*n,1);%multimplicadores de Lagrange
q0 = q(:,i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% POSTPROCESO DE RESULTADOS %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% REPRESENTACION GRAFICA DEL MECANISMO EN TODAS LAS POSICIONES ANALIZADAS
figure(1);
hold on;
axis equal;
RO2 = [0 0]';
RA = PosicionPunto(q(3*(2-1)+1:3*(2-1)+3,i),[norm(z2) 0]');%Esta función
%te va calculando el punto que se desee, en este caso A con respecto al
%eje de coordenadas local de la barra en cada instante de tiempo
RB = PosicionPunto(q(3*(3-1)+1:3*(3-1)+3,i),[norm(z3)/2 0]');
RO4 = PosicionPunto(q(3*(4-1)+1:3*(4-1)+3,i),[0 0]');
RP=PosicionPunto(q(3*(4-1)+1:3*(4-1)+3,i),[norm(z4)/2 0]');
RO5=z1;
RC = PosicionPunto(q(3*(5-1)+1:3*(5-1)+3,i),[-norm(z5)/2 0]');
RD = PosicionPunto(q(3*(6-1)+1:3*(6-1)+3,i),[-norm(z6)/2 0]');
RQ=PosicionPunto(q(3*(6-1)+1:3*(6-1)+3,i),[norm(z6)/2 0]');
X=[RO2(1) RA(1) RB(1) RO4(1) RP(1)];
Xd=[RO2(1) RC(1) RD(1) RO4(1) RQ(1)];
Y=[RO2(2) RA(2) RB(2) RO4(2) RP(2)];
Yd=[RO2(2) RC(2) RD(2) RO4(2) RQ(2)];
plot(X,Y,'b',Xd,Yd,'g');
end
% ANIMACIÓN
RO2 = [0 0]';
RA = PosicionPunto(q(3*(2-1)+1:3*(2-1)+3,i),[norm(z2) 0]');
RO5 = z1;
X=[RO2(1) RA(1) RB(1) RO4(1) RP(1)];
Y=[RO2(2) RA(2) RB(2) RO4(2) RP(2)];
Xd=[RO2(1) RC(1) RD(1) RQ(1)];
Yd=[RO2(2) RC(2) RD(2) RQ(2)];
72
if i ==1
figure(2);
hold on;
axis equal;
axis([-2 4 -2 2]);
h = plot(X,Y,'XDataSource','X','YDataSource','Y');
h2 = plot(Xd,Yd,'g','XDataSource','Xd','YDataSource','Yd');
else
refreshdata(h,'caller')
refreshdata(h2,'caller')
pause(0.1);
drawnow;
end
end
%Representación Par Motor
figure(3);
plot(tspan,lam(18,:));
title('Par Motor ');
ylabel('Par (Nm)');
xlabel('Tiempo (s)');
%Representación Potencia Instantánea

figure(4);
plot(tspan,lam(18,:)*w);
title('Potencia Instantánea');
ylabel('Potencia (W)');
%Representación Trabajo motor

Wmot(1) = 0;
for i = 2:length(tspan)
Wmot(i) = Wmot(i-1)+(tspan(i)-tspan(i-1))*(lam(18,i)*w);
end
figure(5);
plot(tspan,Wmot);
title('Trabajo Motor');
ylabel('Trabajo (J)');
media_par=mean(lam(18,:))
media_potencia=mean(lam(18,:)*w)
rms_par=sqrt(mean((lam(18,:)).^2))
%Cálculo RMS de la potencia. Este es el valor que vamos a utilizar para
%el estudio económico
rms_potencia=sqrt(mean((lam(18,:)*w).^2))
73
74
CAPÍTULO 5: ESTUDIO ECONÓMICO
En este capítulo se realizará un estudio económico sobre lo que supondría la colocación de un
generador eléctrico a nuestro mecanismo.
5.1. INTRODUCCIÓN.
En este capítulo del proyecto se trata de analizar la rentabilidad del mismo proyecto. La energía
mecánica generada por la fuerza humana se va a aprovechar para convertirla a energía eléctrica. Para
ello vamos a tener que colocar un generador eléctrico en cada bicicleta elíptica de manera que todos
ellos vayas conectados mediante los cables a una batería.
Para realizar el estudio vamos a tener en cuenta una serie de variables de nuestro mecanismo que va a
ser clave para la obtención de los resultados. Dichas variables son la frecuencia angular (rpm) de la
barra motriz, es decir, la barra 2 y la masa de la persona que se monte en la bicicleta (kg). Dándole
valores a las variables y combinándolas, obtendremos la potencia mecánica generada en nuestro
mecanismo (W).
Una vez obtenida la potencia mecánica generada buscaremos un generador que se adapte a nuestras
condiciones lo máximo posible, y con ello obtendremos la potencia eléctrica generada, la cual
utilizaremos para realizar nuestro estudio.
5.2 ELECCIÓN DEL GENERADOR
Para poder convertir la energía mecánica generada por la acción humana a través de nuestro
mecanismo a energía eléctrica necesitamos un generador eléctrico que se adapte lo máximo posible a
nuestras condiciones, especialmente de par (N·m) y frecuencia angular (rpm).
Para ello hemos estado barajando tres opciones:
Dinamo: La dinamo es una máquina eléctrica que, absorbiendo energía mecánica, genera una
corriente eléctrica pulsante, que en la práctica puede considerarse como continua, cuya tensión
depende de la velocidad de rotación: aumentando el número de revoluciones aumenta también la
tensión. Dado que la instalación eléctrica de los automóviles trabaja a tensión fija (6, 12 ó 24 V), la
tensión de la dinamo debe mantenerse constante mediante un sistema de regulación.
Las dinamos o generadores eléctricos para ruedas siempre han estado presentes en la historia del
ciclismo, principalmente para proveer de la energía necesaria a infinidad de modelos de luces para
bicicletas.
Alternador: Aparato o dispositivo electromagnético que genera corriente eléctrica alterna mediante
la transformación de energía mecánica en energía eléctrica gracias a la inducción producida por un
imán que se mueve en el interior de una bobina.
Motor de corriente continua: El motor de corriente continua es una máquina que convierte la energía
eléctrica continua en mecánica, provocando un movimiento rotatorio.
75
Esta máquina de corriente continua es una de las más versátiles en la industria. Su fácil control de
posición, par y velocidad la han convertido en una de las mejores opciones en aplicaciones de control y
automatización de procesos. Los motores de corriente continua se siguen utilizando en muchas
aplicaciones de potencia (trenes y tranvías) o de precisión (máquinas, micro motores, etc.)
La principal característica del motor de corriente continua es la posibilidad de regular la velocidad

desde vacío a plena carga.
Su principal inconveniente, el mantenimiento, muy caro y laborioso.
Una máquina de corriente continua (generador o motor) se compone principalmente de dos partes, un
estator que da soporte mecánico al aparato y tiene un hueco en el centro generalmente de forma
cilíndrica. En el estator además se encuentran los polos, que pueden ser de imanes permanentes o
devanados con hilo de cobre sobre núcleo de hierro. El rotor es generalmente de forma cilíndrica,
también devanado y con núcleo, al que llega la corriente mediante dos escobillas.
Los motores y los generadores de corriente continua están constituidos esencialmente por los mismos
elementos, diferenciándose únicamente en la forma de utilización. Por reversibilidad entre el motor y
el generador se entiende que si se hace girar al rotor, se produce en el devanado inducido una fuerza
electromotriz capaz de transformarse en energía en el circuito de carga. En cambio, si se aplica una
tensión continua al devanado inducido del generador a través del colector de delgas, el
comportamiento de la máquina ahora es de motor, capaz de transformar la fuerza contraelectromotriz
en energía mecánica. En ambos casos el inducido está sometido a la acción del campo inductor
principal.
De entre estos tres aparatos, decidimos utilizar un motor de corriente continua ya que después de
buscar información, fue el que más se adaptaba a nuestras condiciones de potencia y frecuencia
angular, con un rendimiento admisible de nuestro mecanismo. Además resultaba algo más barato.
El motor elegido es un motor de eje utilizado en bicicletas eléctricas. Véanse figuras 5.1 y 5.2.
Marca GoldenMotor
Model:PW-12H
Voltage:24V
Power:180W-300W
Weight: 7.2 Kgs
FIGURA 5.1: MOTOR
76
FIGURA 5.2: CARACTERÍSTICAS TÉCNICAS DEL MOTOR
Con estos datos ya podemos saber cuál es la energía eléctrica que genera la acción humana en este
tipo de bicicletas.
Para calcular la energía vamos a tener como variables la frecuencia angular (rpm) de la barra motriz
(barra 2) y la masa de la persona.
Vamos a tomar los siguientes valores
w(rpm)= 70, 75,80,85,90 y 95 rpm
masa(kg)= 55,60,65,70,80,90 y 100 kg
Hay que aclarar que la resistencia que opone el generador, que denominaremos “c”, va a variar
dependiendo de la frecuencia angular. Véase tabla 5.1.
77
Generador:
M=c·w → c=M/w
TABLA 5.1.
GENERADOR ELÉCTRICO
Frecuencia Angular (rpm) 70 75 80 85 90 95
Par Motor (N·m) 24,5 22,4 19,34 16,13 12,76 9,5
Constante "c" 3,34 2,85 2,31 1,81 1,35 0,95
En el anexo B, se muestra las gráficas tipos que obtenemos al realizar las simulaciones de cada
combinación. La potencia generada que utilizaremos será la RMS de la potencia representada.
También añadir que, como se observa en la tabla 5.1, el rendimiento del generador va a variar
dependiendo de la frecuencia angular empleada.
Las fórmulas que se emplearán para el estudio son:
Potencia eléctrica obtenida (kW)) =potencia generada (kW) x rendimiento generador
Energía eléctrica diaria unitaria (kWh)= horas trabajo al día x potencia eléctrica obtenida(kW).
Ahorro diario unitario = coste del kWh x energía eléctrica diaria unitaria (kWh)
Para los cálculos hemos considerado que las bicicletas trabajarán 10 horas al día y el coste del kWh de
0,16 €/kWh.
Realizando una combinación entre las dos variables de nuestro problema (frecuencia angular y masa)
obtenemos los resultados.
En las siguientes tablas, se muestran los resultados obtenidos para las diferentes masas de personas
elegidas.
TABLA 5.2. MASA 55 Kg.
Potencia Energía
Potencia Ahorro diario
Velocidad Horas Precio Rendimiento eléctrica eléctrica
generada unitario sin
(rpm) trabajo/día Kw·h generador obtenida diaria
(kW) perdidas
(kW) unitaria(kWh)
70 10 0,160 € 0,180 57,00% 0,102 1,024 0,164 €
75 10 0,160 € 0,176 60,00% 0,106 1,057 0,169 €
80 10 0,160 € 0,163 63,75% 0,104 1,038 0,166 €
85 10 0,160 € 0,145 66,97% 0,097 0,969 0,155 €
90 10 0,160 € 0,122 70,37% 0,086 0,858 0,137 €
95 10 0,160 € 0,098 72,89% 0,072 0,716 0,115 €
78
TABLA 5.3. MASA 60 Kg
Potencia Consumo
Velocidad Horas Precio Rendimiento eléctrica eléctrico
(rpm) trabajo/día Kw·h generador obtenida diario
(kW) perdidas
(kW) unitario(kWh)
70 10 0,160 € 0,180 57,00% 0,103 1,025 0,164 €
75 10 0,160 € 0,176 60,00% 0,106 1,058 0,169 €
80 10 0,160 € 0,163 63,75% 0,104 1,040 0,166 €
85 10 0,160 € 0,145 66,97% 0,097 0,972 0,155 €
90 10 0,160 € 0,123 70,37% 0,086 0,863 0,138 €
95 10 0,160 € 0,099 72,89% 0,072 0,724 0,116 €
Potencia Consumo
(kW) perdidas
(kW) unitario(kWh)
70 10 0,160 € 0,180 57,00% 0,103 1,026 0,164 €
75 10 0,160 € 0,177 60,00% 0,106 1,059 0,170 €
80 10 0,160 € 0,164 63,75% 0,104 1,042 0,167 €
85 10 0,160 € 0,146 66,97% 0,098 0,976 0,156 €
90 10 0,160 € 0,124 70,37% 0,087 0,869 0,139 €
95 10 0,160 € 0,101 72,89% 0,074 0,735 0,118 €
Potencia Consumo
(kW) perdidas
(kW) unitario(kWh)
70 10 0,160 € 0,180 57,00% 0,103 1,027 0,164 €
75 10 0,160 € 0,177 60,00% 0,106 1,061 0,170 €
80 10 0,160 € 0,164 63,75% 0,105 1,045 0,167 €
85 10 0,160 € 0,146 66,97% 0,098 0,980 0,157 €
90 10 0,160 € 0,125 70,37% 0,088 0,877 0,140 €
95 10 0,160 € 0,103 72,89% 0,075 0,750 0,120 €
Potencia Consumo
(kW) perdidas
(kW) unitario(kWh)
70 10 0,160 € 0,181 57,00% 0,103 1,030 0,165 €
75 10 0,160 € 0,178 60,00% 0,107 1,066 0,171 €
80 10 0,160 € 0,165 63,75% 0,105 1,053 0,169 €
85 10 0,160 € 0,148 66,97% 0,099 0,993 0,159 €
90 10 0,160 € 0,128 70,37% 0,090 0,899 0,144 €
95 10 0,160 € 0,108 72,89% 0,079 0,787 0,126 €
79
Potencia Consumo
(kW) perdidas
(kW) unitario(kWh)
70 10 0,160 € 0,181 57,00% 0,103 1,034 0,165 €
75 10 0,160 € 0,179 60,00% 0,107 1,072 0,171 €
80 10 0,160 € 0,167 63,75% 0,106 1,064 0,170 €
85 10 0,160 € 0,151 66,97% 0,101 1,011 0,162 €
90 10 0,160 € 0,132 70,37% 0,093 0,927 0,148 €
95 10 0,160 € 0,115 72,89% 0,084 0,835 0,134 €
Potencia Consumo
(kW) perdidas
(kW) unitario(kWh)
70 10 0,160 € 0,182 57,00% 0,104 1,039 0,166 €
75 10 0,160 € 0,180 60,00% 0,108 1,079 0,173 €
80 10 0,160 € 0,169 63,75% 0,108 1,076 0,172 €
85 10 0,160 € 0,154 66,97% 0,103 1,032 0,165 €
90 10 0,160 € 0,137 70,37% 0,096 0,962 0,154 €
95 10 0,160 € 0,122 72,89% 0,089 0,892 0,143 €
Como se puede observar con una bicicleta elíptica se puede obtener como máximo un ahorro diario de
0,173 €. Por lo que vamos a centrar el estudio, en cuantas bicicletas se necesitará para ahorrarnos 1 €
al día. Para ello el cálculo que realizamos es hacer la media de todos los ahorros diarios unitarios
obtenidos, que resulta 0,154 €.
Entonces para saber el número de bicicletas necesario para ahorrarnos 1 € al día hacemos:
1/0,154 = 6,49 → el número de bicicletas es 7.
Sabiendo el número de bicicletas vamos a proceder a la elección del cable, que irá desde los
generadores de cada bicicleta a la batería, y después procederemos a la elección de la batería.
5.3 ELECCIÓN DEL CABLE
En este apartado se va a elegir el cable que mejor se ajuste a nuestra instalación, ya que en el tramo
que va desde los generadores a la batería se producirá una caída de tensión, y por tanto de potencia,
la cual tiene que ser lo más pequeña posible. Para ello se hará un pequeño cálculo, consistente en
hallar un cable con la sección suficiente para que la caída de tensión sea lo más pequeña posible.
80
El cálculo de secciones de líneas eléctricas es un método de cálculo para obtener la sección idónea del
conductor a emplear, siendo este capaz de:
 transportar la potencia requerida con total seguridad;
 que dicho transporte se efectúe con un mínimo de pérdidas de energía;
 mantener los costes de instalación en unos valores aceptables.
La caída de tensión ( ) se produce como consecuencia de la resistencia de los conductores. Como

regla general, en España, se permite una ( ) máxima de:
 3% para cualquier circuito interior de viviendas.

 3 % en instalaciones de alumbrado.
 5 % en el resto de instalaciones.
La ecuación a emplear es la siguiente:
siendo: S:sección del cable ; I:intensidad en el tramo
ρ:resistividad del cable ; :caída de tensión tramo
L:Longitud del cable
Para realizar los cálculos se colocará las bicicletas en la posición más desfavorable, es decir, una al lado
de otra separadas una determinada distancia, como se representa en la siguiente figura 5.3.
FIGURA 5.3.
siendo L’: separación entre bicicletas;L’=1,5 m
Lb: longitud desde la última bicicleta a la batería; Lb= 5 m.
81
Nuestra instalación es una línea con diferentes cargas distribuidas uniformemente, entonces hay que
calcular antes el momento eléctrico.
El momento eléctrico de una línea es el producto de la carga eléctrica por la distancia hasta el
origen. Puede considerarse como el equivalente de la línea constituido por un único tramo de línea
con una única carga en su extremo.
 En corriente continua:
=L·I
Y para cargas uniformemente repartidas sería:
Como intensidad generada para nuestros cálculos, vamos a tomar la más alta, siendo:
I (A)= = 7,59 A, entonces:
Entonces para el número de bicicletas n= 7, el momento eléctrico resulta:
M= L·I+L·2·I+L·3·I+…+(n-1)·L·I+n·Lb·I =504,73 Am
Tomando como caída de tensión máxima el 3 %, queda:
= 0,03·24=0,72 V
ρ= 1,71·10-8 Ω·m
Obtenidos todos los parámetros, la sección del cable es:
=23,97 mm2≈ 24 mm2
Por lo que necesitamos un cable de 24 mm2 para nuestra instalación para que la pérdida sea lo más
pequeña posible.
El cable elegido es:
El cable TopFlex MS TRI-RATED
Conductor: Cobre electrolítico, clase 5 (flexible) según EN60228 y BS 6360
Aislamiento: PVC de alta temperatura de servicio tipo TI3 según norma UNE 21031/HD 21 y
Clase 43 según UL 1581. El material especial utilizado para el aislamiento proporciona buenas
propiedades de deslizamiento al cable.
82
Embalaje: Las secciones pequeñas (de 0,75 mm2 hasta 6 mm2) se suministran en cajas de alta
resistencia (ver tabla inferior). Las secciones medias (de 10 mm2 hasta 35 mm2) se suministran
en rollos con film retractilado. Las secciones mayores (> 35 mm2) se suministran en bobinas.
STYLE SECCIÓN EMBALAJE
1015 0,75-6 mm2 Cajas de alta resistencia
1028 10 mm2 Rollos retractilados
1283 16-35 mm2 Rollos retractilados
1283 >50 mm2 Bobinas
Norma Nacional/ Europea: UNE-EN 60332-1
Norma Internacional: IEC 60332-1 / UL 2556
ITC-BT: 30
Una foto del cable se muestra en la figura 5.4
FIGURA 5.4: CABLE
83
5.3 ELECCIÓN DE LA BATERÍA.
En este apartado se va a proceder a la elección de la batería. La batería considerada tiene que ser de
24 V y con la capacidad suficiente para cumplir nuestras condiciones. La intensidad que llega a la
batería consideraremos la más desfavorable, es decir, la que generaría todas las bicicletas a máxima
potencia, por lo que:
Ibateria>n·I=7·7,59= 53,13 A
Para ello colocaremos 5 baterías de 100 Ah en paralelo, ya que para esa intensidad cada batería se
cargará en 2 horas, y al estar las bicicletas funcionando 10 horas al día cada una, entonces
necesitaremos 5 baterías para garantizar que toda esa energía se almacenará en 10 horas.
Otra opción, hubiera sido colocar una batería por cada batería. En este caso la batería tendría que ser
de 70 Ah para que con la intensidad de 7,59 A se cargara en 10 horas las baterías. También nos
ahorraríamos de colocar el cable. Pero realizando un estudio económico resulta que la instalación de 5
baterías de 100 Ah en paralelo resulta algo más rentable.
Las características de la batería elegida se muestra en la siguiente tabla 5.9:
TABLA 5.9
Marca TopBand
Model: TB-24100F
Type: LiFePO4 battery
Nominal voltage: 24 V
Typical capacity: 100Ah
Max continuous discharge current 60 A
Max pulse discharge current : 150 A
Work voltage range: 18.4~29.2V
Charging temperatura 0°C~45°C
Discharging temperatura -20°C~60°C
Storage temperatura -20°C~45°C
Size: 360*300*196mm
Weight: approx.34kg
Uso: Coches,sillas de ruedas, ups, solar…
Ciclo de vida: Más de 2000 ciclos
84
Una foto de la batería se muestra en la figura 5.5.
FIGURA 5.5 BATERÍA
85
86
CAPÍTULO 6: RESULTADOS
En este capítulo vamos a mostrar los resultados obtenidos definitivamente, teniendo en cuenta las
pérdidas en el cable y la batería. Las pérdidas serán un 3% en el cable como se ha explicado
anteriormente y un 2 % en la batería, por lo que vamos a tener unas pérdidas totales del 5%.
Además de mostrar el ahorro anual eléctrico que supondrá la colocación de los generadores, también
mostraremos la inversión que habría que hacer para conseguir los generadores, cables y batería. Una
vez que sabemos el ahorro anual eléctrico y la inversión necesaria, mostraremos cuánto sería el
tiempo de amortización de dicha inversión.
Cálculos realizados para la obtención de los resultados:
Los cálculos realizados para obtener estos resultados son los siguientes:
- Ahorro diario unitario = Ahorro diario sin pérdidas x pérdidas
siendo las pérdidas del 5%, por lo tanto habrá que multiplicarlo por 0,95.
- Ahorro diario = ahorro diario uniatrio x número de bicicletas

- Coste inversión = coste generador + coste cable + coste batería
siendo sus respectivos costes:
Coste generador= 90 €/ud.; Coste batería = 68,3 €/ud.; Coste cable = 2,2 €/metro.
Por tanto como tenemos 7 bicicletas y vamos a necesitar 30 metros de cable:
Coste inversión = 90·7 + 68,3·5 + 2,2·30= 1037,50 €
- Ahorro anual = ahorro diario x número de días laborales al año
siendo el número de días laborales considerado 252.
Amortización= Coste inversión/ahorro anual
En las siguientes tablas, se muestran los resultados obtenidos para las diferentes masas de personas
elegidas.
Ahorro diario Ahorro

Velocidad Número de Ahorro Coste Ahorro Amortiza Amortización
unitario sin diario
(rpm) bicicletas diario total inversión anual ción en años
perdidas unitario
70 0,164 € 0,156 € 7,00 1,09 € 1.037,50 € 274,66 € 3,78 3años y9meses

75 0,169 € 0,161 € 7,00 1,12 € 1.037,50 € 283,43 € 3,66 3años y8meses
80 0,166 € 0,158 € 7,00 1,10 € 1.037,50 € 278,42 € 3,73 3años y9meses
85 0,155 € 0,147 € 7,00 1,03 € 1.037,50 € 259,82 € 3,99 3años y12meses
90 0,137 € 0,130 € 7,00 0,91 € 1.037,50 € 230,14 € 4,51 4años y6meses
95 0,115 € 0,109 € 7,00 0,76 € 1.037,50 € 192,04 € 5,40 5años y5meses
87

unitario sin diario
(rpm) bicicletas diario inversión anual ción en años
perdidas unitario
70 0,164 € 0,156 € 7,00 1,09 € 1.037,50 € 274,84 € 3,77 3años y9meses

75 0,169 € 0,161 € 7,00 1,13 € 1.037,50 € 283,69 € 3,66 3años y8meses
80 0,166 € 0,158 € 7,00 1,11 € 1.037,50 € 278,87 € 3,72 3años y9meses
85 0,155 € 0,148 € 7,00 1,03 € 1.037,50 € 260,55 € 3,98 3años y12meses
90 0,138 € 0,131 € 7,00 0,92 € 1.037,50 € 231,33 € 4,48 4años y6meses
95 0,116 € 0,110 € 7,00 0,77 € 1.037,50 € 194,17 € 5,34 5años y4meses

unitario sin diario
perdidas unitario
70 0,164 € 0,156 € 7,00 1,09 € 1.037,50 € 275,08 € 3,77 3años y9meses

75 0,170 € 0,161 € 7,00 1,13 € 1.037,50 € 284,05 € 3,65 3años y8meses
80 0,167 € 0,158 € 7,00 1,11 € 1.037,50 € 279,49 € 3,71 3años y9meses
85 0,156 € 0,148 € 7,00 1,04 € 1.037,50 € 261,58 € 3,97 3años y12meses
90 0,139 € 0,132 € 7,00 0,92 € 1.037,50 € 233,03 € 4,45 4años y5meses
95 0,118 € 0,112 € 7,00 0,78 € 1.037,50 € 197,19 € 5,26 5años y3meses

unitario sin diario
perdidas unitario
70 0,164 € 0,156 € 7,00 1,09 € 1.037,50 € 275,38 € 3,77 3años y9meses

75 0,170 € 0,161 € 7,00 1,13 € 1.037,50 € 284,52 € 3,65 3años y8meses
80 0,167 € 0,159 € 7,00 1,11 € 1.037,50 € 280,29 € 3,70 3años y8meses
85 0,157 € 0,149 € 7,00 1,04 € 1.037,50 € 262,89 € 3,95 3años y11meses
90 0,140 € 0,133 € 7,00 0,93 € 1.037,50 € 235,22 € 4,41 4años y5meses
95 0,120 € 0,114 € 7,00 0,80 € 1.037,50 € 201,06 € 5,16 5años y2meses

unitario sin diario
perdidas unitario
70 0,165 € 0,157 € 7,00 1,10 € 1.037,50 € 276,17 € 3,76 3años y9meses

75 0,171 € 0,162 € 7,00 1,13 € 1.037,50 € 285,76 € 3,63 3años y8meses
80 0,169 € 0,160 € 7,00 1,12 € 1.037,50 € 282,39 € 3,67 3años y8meses
85 0,159 € 0,151 € 7,00 1,06 € 1.037,50 € 266,38 € 3,89 3años y11meses
90 0,144 € 0,137 € 7,00 0,96 € 1.037,50 € 241,03 € 4,30 4años y4meses
95 0,126 € 0,120 € 7,00 0,84 € 1.037,50 € 211,14 € 4,91 4años y11meses
88

unitario sin diario
perdidas unitario
70 0,165 € 0,157 € 7,00 1,10 € 1.037,50 € 277,20 € 3,74 3años y9meses

75 0,171 € 0,163 € 7,00 1,14 € 1.037,50 € 287,39 € 3,61 3años y7meses
80 0,170 € 0,162 € 7,00 1,13 € 1.037,50 € 285,16 € 3,64 3años y8meses
85 0,162 € 0,154 € 7,00 1,08 € 1.037,50 € 270,96 € 3,83 3años y10meses
90 0,148 € 0,141 € 7,00 0,99 € 1.037,50 € 248,64 € 4,17 4años y2meses
95 0,134 € 0,127 € 7,00 0,89 € 1.037,50 € 223,99 € 4,63 4años y8meses

unitario sin diario
perdidas unitario
70 0,166 € 0,158 € 7,00 1,11 € 1.037,50 € 278,49 € 3,73 3años y9meses

75 0,173 € 0,164 € 7,00 1,15 € 1.037,50 € 289,42 € 3,58 3años y7meses
80 0,172 € 0,164 € 7,00 1,15 € 1.037,50 € 288,58 € 3,60 3años y7meses
85 0,165 € 0,157 € 7,00 1,10 € 1.037,50 € 276,58 € 3,75 3años y9meses
90 0,154 € 0,146 € 7,00 1,02 € 1.037,50 € 257,88 € 4,02 4años y0meses
95 0,143 € 0,136 € 7,00 0,95 € 1.037,50 € 239,16 € 4,34 4años y4meses
En la siguiente figura 6.1 se representa una gráfica que muestra el ahorro eléctrico anual frente a la
frecuencia angular de cada masa de la persona elegida.
AHORRO ELÉCTRICO ANUAL

350,00 €
300,00 €
ahorro eléctrico anual (€)
250,00 € MASA 55 Kg
MASA 60 Kg
200,00 €
MASA 65 Kg
150,00 €
MASA 70 Kg
100,00 € MASA 80 Kg
50,00 € MASA 90 Kg
MASA 100 Kg
0,00 €
70 75 80 85 90 95
frecuencia angular (rpm)
FIGURA 6.1
89
También se muestra una gráfica de la amortización frente a la frecuencia angular de cada masa de
persona elegida. Véase figura 6.2
AMORTIZACIÓN
6,00
5,00
MASA 55 Kg
4,00
Amortización
MASA 60 Kg
3,00 MASA 65 Kg
MASA 70 Kg
2,00
MASA 80 Kg
1,00 MASA 90 Kg
MASA 100 Kg
0,00
70 75 80 85 90 95
frecuencia angular (rpm)
FIGURA 6.2
Como se puede observar en las gráficas el ahorro eléctrico es importante al cabo de un año para
cualquier frecuencia angular y masa, siendo poca la diferencia de ahorro entre masas para
revoluciones bajas. Según va aumentando la frecuencia angular la masa va tomando importancia, y por
tanto la inercia también aumentará, entonces transmitirá un par mayor a la barra motriz (barra 2), y
por consecuencia mayor la potencia mecánica generada.
Al generar más potencia mecánica, se genera más potencia eléctrica, y por tanto un mayor ahorro
anual. Como el coste de inversión es constante, entonces la amortización será menor.
90
CAPÍTULO 7: CONCLUSIÓN
En este capítulo se comentan las conclusiones finales del proyecto.
Este proyecto como se observa ha dado unos resultados bastante positivos. Las bicicletas elípticas, son
unas máquinas de uso indoor o de gimnasio exclusivamente. Aprovechando la acción humana que las
personas ejercen voluntariamente para tonificar su cuerpo, conseguimos aproximadamente de media,
que en tres años y pocos meses amortizar unos gastos de inversión que después supondrá reducir
nuestros gastos eléctricos una cantidad bastante considerable. Por tanto, colocar los generadores a las
bicicletas elípticas es una buena solución para gimnasios que contengan un número elevado de
bicicletas.
En la siguiente tabla 7.1 se muestra el ahorro y la amortización que se obtiene al aumentar el número
de bicicletas existentes.
Tabla 7.1
Ahorro
diario Ahorro
Número de Coste Ahorro Amortiza Amortización
diario total
unitario bicicletas inversión anual total ción en años
medio
medio
0,146 € 7 1,02 € 1.037,50 € 257,75 € 4,03 4años y0meses
0,146 € 8 1,17 € 1.197,70 € 294,58 € 4,07 4años y1meses
0,146 € 9 1,32 € 1.357,90 € 331,40 € 4,10 4años y1meses
0,146 € 10 1,46 € 1.518,10 € 368,22 € 4,12 4años y1meses
0,146 € 11 1,61 € 1.678,30 € 405,04 € 4,14 4años y2meses
0,146 € 12 1,75 € 1.838,50 € 441,87 € 4,16 4años y2meses
0,146 € 14 2,05 € 2.075,00 € 515,51 € 4,03 4años y0meses
0,146 € 21 3,07 € 3.112,50 € 773,26 € 4,03 4años y0meses
Entonces en la gráfica de la figura 7.1 se puede evaluar el ahorro eléctrico y el coste de inversión
frente al número de bicicletas.
3.500,00 €
3.000,00 €
2.500,00 €
2.000,00 €
Ahorro eléctrico anual
1.500,00 € medio
1.000,00 € Coste inversión
500,00 €
0,00 €
7 8 9 10 11 12 14 21
Número de bicicletas
FIGURA 7.1
91
También se representa gráficamente la amortización frente al número de bicicletas. Véase figura 7.2.
Amortización
4,18
4,16
4,14
Amortización (años)
4,12
4,10
4,08
Amortización
4,06
4,04
4,02
4,00
0 5 10 15 20 25
Número de bicicletas
FIGURA 7.2.
Como se observa en la tabla y gráficas, hemos ido aumentando el número de bicicletas para hacer una
comparativa de cómo influye dicho número. La conclusión final es que mientras mayor sea el número
de bicicletas instaladas en el gimnasio, mayor será nuestro coste de inversión (ya que hay que comprar
más generadores y metros de cable), pero por consiguiente se obtiene un ahorro económico superior,
siendo finalmente el periodo de amortización similar para múltiplos de 7. Este hecho se produce ya
que cada 7 bicicletas sólo colocaríamos 5 baterías de 100 Ah en paralelo con sus respectivos metros de
cable. Mientras que si, por ejemplo, colocamos 8 bicicletas, entonces a la bicicleta número 8
tendríamos que colocarle una batería de 70 Ah, y dejar a las otras 7 con la instalación calculada
anteriormente. Esto hace que el coste aumente un poco más de lo normal, produciendo que la
amortización aumente un poco.
Aun así, aprovechar esa fuerza humana voluntaria a la larga nos proporciona beneficio.
92
CAPÍTULO 8: BIBLIOGRAFÍA
En este capítulo se muestran las referencias utilizadas para el estudio y elaboración del proyecto fin de
carrera.
(1) F. Reuleaux.” The kinematics of Machinery”. Dover, Nueva York (1963)
(2)J.N. Kozhevnikov. “Mecanismos”. Gili, B arcelona (1970)
(3) J.A Hrones y G.L Nelson. “Analysis of the Four-bar Linkage”.John Wiley, Nueva York (1951)
(4) Apuntes cinemática y dinámica de máquinas. (2011)
(5) J. Nieto “Síntesis de mecanismos”
Referencias web
www.goldenmotor.es
spanish.alibaba.com
www.topcable.com
www.wikipedia.es
nuestrolifestyle.com/tag/bicicletas-elipticas
www.kettler.net
93
94
CAPÍTULO 9: ANEXO
En este capítulo se muestra el código empleado para poder realizar la simulación cinemática y
dinámica, y así obtener la potencia mecánica que es capaz de generar la acción humana.
9.1. ANEXO A
En este anexo se muestra el código explicado en el capítulo 4 para que sea fácil su referencia y de
forma ordenada.
Script del programa ejecutable
global t w;
global inercias g;
%ESTUDIO CINEMÁTICO
z2=[0.24*cos(0) 0.24*sin(0)]';%vector barra 2
z1=[0.97*cos(0.578) 0.97*sin(0.578)]';%vector barra fija 1
n=6; nr=7; np=0;
%Son vectores que nos indician las barras que intervienen en cada rótula
R(1:2,1)=[1 2]';%por ejemplo en la rótula 1 interviene la barra 1 y la 2
R(1:2,2)=[2 3]';
R(1:2,3)=[3 4]';
R(1:2,4)=[4 1]';
R(1:2,5)=[2 5]';
R(1:2,6)=[5 6]';
R(1:2,7)=[6 1]';
%Estos vectores indican la distancia de la rótula al eje local definido en

%cada barra
DR(1:2,1)=[0 0]';
DR(3:4,1)=[0 0]';
DR(1:2,2)=[norm(z2) 0]';
DR(3:4,2)=[-norm(z3)/2 0]';
DR(1:2,3)=[norm(z3)/2 0]';
DR(3:4,3)=[-norm(z4)/2 0]';
DR(1:2,4)=[0 0]';
DR(3:4,4)=[0.97*cos(0.578) 0.97*sin(0.578)]';
95
DR(1:2,5)=[-norm(z2) 0]';
DR(3:4,5)=[-norm(z5)/2 0]';
DR(1:2,6)=[norm(z5)/2 0]';
DR(3:4,6)=[-norm(z6)/2 0]';
%
DR(1:2,7)=[0 0]';
DR(3:4,7)=[0.97*cos(0.578) 0.97*sin(0.578)]';
%PROPIEDADES INERCIALES DE LOS SÓLIDOS
%(masa, momento de inercia respecto a G y posición de G en globales)
inercias(:,1) = [0 0 [0 0]]';
inercias(:,2) = [14.205 ((norm(z2)^2)*14.205/2+(0.046)) [0 0]]';
inercias(:,3) = [1.738 ((norm(z3)^2)*1.738/12 + 50*(0.2^2)) [0 0]]';
inercias(:,4) = [1.892 ((norm(z4)^2)*1.892/12) [0 0]]';
inercias(:,5) = [1.738 ((norm(z5)^2)*1.738/12 + 50*(0.2^2)) [0 0]]';
inercias(:,6) = [1.892 ((norm(z6)^2)*1.892/12) [0 0]]';
%Dirección y sentido de la gravedad
g = [0 -9.81]';
%%%%%%%%%%%%%%%%%%%%%%%%%
% SIMULACIÓN CINEMÁTICA %
%%%%%%%%%%%%%%%%%%%%%%%%%
IND = [];
DEP = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18]';
%Representamos los ángulos de cada barra en la condicional inicial elegida
fi20 = atan2(z2(2),z2(1))
psi20 = atan2(z4(2),z4(1))
teta20 = atan2(z3(2),z3(1))
psid20 = atan2(z6(2),z6(1))
tetad20 = atan2(z5(2),z5(1))
w = 95*(2*pi/60); %VELOCIDAD ANGULAR (95 rpm)

tspan = 0:0.05:(9*pi/w); %TIEMPO DE INTEGRACIÓN
%Estimacion inicial
q0=[0 0 0 0 0 fi20 (z2+(z3./2))' teta20 (z2+z3+(z4./2))' psi20 (-

z2+(z5./2))' tetad20 (-z2+z5+(z6./2))' psid20]';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BUCLE DE SOLUCIÓN DEL PROBLEMA DE POSICIÓN %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
q = zeros(3*n, length(tspan));
v = zeros(3*n, length(tspan));
a = zeros(3*n, length(tspan));
lam = zeros(3*n, length(tspan));
96
t = tspan(i);
%CÁLCULO DE POSICIONES
q(:,i) = NewtonRaphsonJacob(@Restricciones,@Jacobiano,q0);
Cq = Jacobiano(q(:,i));%Calcula el jacobiano
Ct = DtRestr;
%CALCULO DE VELOCIDADES
v(:,i) = -Cq\Ct;
Bq = MatrizBq(q(:,i),v(:,i));
Ctt = DttRestr;
%CALCULO DE ACELERACIONES Y MULTIPLICADORES DE LAGRANGE

gam = -(Bq*v(:,i)+Ctt);
M = MatrizMasa(q(:,i));
Qv = FCentrifuga(q(:,i));
Qgrav = FGravedad(q(:,i));
Qfriccion = [0 0 0 0 0 0.95*w*1 0 0 0 0 0 0 0 0 0 0 0 0]';
AA = zeros(6*n,6*n);
bb = zeros(6*n,1);
AA(1:3*n,1:3*n) = M;
AA(1:3*n,3*n+1:6*n) = Cq';
AA(3*n+1:6*n,1:3*n) = Cq;
bb(1:3*n,1) = Qv+Qgrav+Qfriccion;
bb(3*n+1:6*n,1) = gam;
xx = AA\bb;
a(:,i) = xx(1:3*n,1);
lam(:,i) = xx(3*n+1:6*n,1);%multimplicadores de Lagrange
q0 = q(:,i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% POSTPROCESO DE RESULTADOS %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% REPRESENTACION GRAFICA DEL MECANISMO EN TODAS LAS POSICIONES ANALIZADAS
figure(1);
hold on;
axis equal;
RO2 = [0 0]';
RA = PosicionPunto(q(3*(2-1)+1:3*(2-1)+3,i),[norm(z2) 0]');%Esta función
%te va calculando el punto que se desee, en este caso A con respecto al
97
%eje de coordenadas local de la barra en cada instante de tiempo
RO5=z1;
X=[RO2(1) RA(1) RB(1) RO4(1) RP(1)];
Xd=[RO2(1) RC(1) RD(1) RO4(1) RQ(1)];
Y=[RO2(2) RA(2) RB(2) RO4(2) RP(2)];
Yd=[RO2(2) RC(2) RD(2) RO4(2) RQ(2)];
plot(X,Y,'b',Xd,Yd,'g');
end
% ANIMACIÓN
RO2 = [0 0]';
RA = PosicionPunto(q(3*(2-1)+1:3*(2-1)+3,i),[norm(z2) 0]');
RO5 = z1;
X=[RO2(1) RA(1) RB(1) RO4(1) RP(1)];
Y=[RO2(2) RA(2) RB(2) RO4(2) RP(2)];
Xd=[RO2(1) RC(1) RD(1) RQ(1)];
Yd=[RO2(2) RC(2) RD(2) RQ(2)];
if i ==1
figure(2);
hold on;
axis equal;
axis([-2 4 -2 2]);
h = plot(X,Y,'XDataSource','X','YDataSource','Y');
h2 = plot(Xd,Yd,'g','XDataSource','Xd','YDataSource','Yd');
else
refreshdata(h,'caller')
refreshdata(h2,'caller')
pause(0.1);
drawnow;
end
end
%Representación Par Motor
figure(3);
plot(tspan,lam(18,:));
title('Par Motor ');
ylabel('Par (Nm)');
%Representación Potencia Instantánea

figure(4);
plot(tspan,lam(18,:)*w);
title('Potencia Instantánea');
98
ylabel('Potencia (W)');
%Representación Trabajo motor

Wmot(1) = 0;
for i = 2:length(tspan)
Wmot(i) = Wmot(i-1)+(tspan(i)-tspan(i-1))*(lam(18,i)*w);
end
figure(5);
plot(tspan,Wmot);
title('Trabajo Motor');
ylabel('Trabajo (J)');
media_par=mean(lam(18,:))
media_potencia=mean(lam(18,:)*w)
rms_par=sqrt(mean((lam(18,:)).^2))
%Cálculo RMS de la potencia. Este es el valor que vamos a utilizar para
%el estudio económico
rms_potencia=sqrt(mean((lam(18,:)*w).^2))
FUNCIONES PROBLEMA POSICIÓN
Función NewtonRaphsonJacob.m
function x = NewtonRaphsonJacob(fun,jac,x0)
%PARÁMETROS DEL NEWTONRAPHSON
AbsErr=10^(-3);
MaxIter=100;
Niter=0;
ok=0;
while ((ok==0)&&(Niter<MaxIter))
Niter=Niter+1;
CC = feval(fun,x0);
Cq = feval(jac,x0);
% CONVERGENCIA
% max(abs(CC))
if (max(abs(CC))<AbsErr)
ok=1;
end
% ITERACIÓN DE NEWTON-RAPHSON
x=x0-Cq\CC;
99
x0=x;
end
Función Restricciones.m
function C = Restricciones(q)
global t w;
C = RestrPares(q);
%RESTRICCIÓN DE MOVILIDAD
C(3*n,1) = q(6)+w*t;
Función RestrPares.m
function C = RestrPares(qdep)
for i=1:length(IND)
q(IND(i))=qind(i);
end
for i=1:length(DEP)
q(DEP(i))=qdep(i);
end
%Barra fija
C(1,1)=q(1);
C(2,1)=q(2);
C(3,1)=q(3);
%Rótulas
for k=1:nr;
i=R(1,k);
j=R(2,k);
Ri=[q((i-1)*3+1) q((i-1)*3+2)]';
teti=q((i-1)*3+3);
Rj=[q((j-1)*3+1) q((j-1)*3+2)]';
tetj=q((j-1)*3+3);
ri=[DR(1,k) DR(2,k)]';
rj=[DR(3,k) DR(4,k)]';
Rr = Rotula(Ri,teti,ri,Rj,tetj,rj);
100
C(3+2*(k-1)+1,1)=Rr(1);
C(3+2*(k-1)+2,1)=Rr(2);
end
%Prismáticos
for k=1:np;
i=P(1,k);
j=P(2,k);
Ri=[q((i-1)*3+1) q((i-1)*3+2)]';
teti=q((i-1)*3+3);
Rj=[q((j-1)*3+1) q((j-1)*3+2)]';
tetj=q((j-1)*3+3);
ri=[DP(1,k) DP(2,k)]';
hi=[DP(3,k) DP(4,k)]';
rj=[DP(5,k) DP(6,k)]';
Rp = Prism(Ri,teti,ri,hi,Rj,tetj,rj);
C(3+2*nr+2*(k-1)+1,1)=Rp(1);
C(3+2*nr+2*(k-1)+2,1)=Rp(2);
end
Función Rotula.m
function C = Rotula(Ri,teti,ri,Rj,tetj,rj);
Ai = Rotmat(teti);
Aj = Rotmat(tetj);
C = Ri+Ai*ri-Rj-Aj*rj;
Función RotMat.m
function A=Rotmat(x)
A(1,1)=cos(x);
A(1,2)=-sin(x);
A(2,1)=sin(x);
A(2,2)=cos(x);
Función Prism.m
function C = Prism(Ri,teti,ri,hi,Rj,tetj,rj)
Ai = Rotmat(teti);
Aj = Rotmat(tetj);
101
C(1)=teti-tetj;
C(2)=(Ai*hi)'*(Ri+Ai*ri-Rj-Aj*rj);
Función Jacobiano.m
function Cq = Jacobiano(q)
global n;
global t w;
Cq = JacobPares(q);
Cq(3*n,6) = 1;
Función JacobPares.m
function Cqdep = JacobPares(qdep)
for i=1:length(IND)
q(IND(i))=qind(i);
end
for i=1:length(DEP)
q(DEP(i))=qdep(i);
end
%JACOBIANO RESPECTO A TODAS LAS COORDENADAS
m = 2*(nr+np)+3; %Número de restricciones
Cq = zeros(m,3*n);
%Barra fija
Cq(1,1)=1;
Cq(2,2)=1;
Cq(3,3)=1;
%Rótulas
for k=1:nr;
i=R(1,k);
j=R(2,k);
Ri=[q((i-1)*3+1) q((i-1)*3+2)]';
teti=q((i-1)*3+3);
102
Rj=[q((j-1)*3+1) q((j-1)*3+2)]';
tetj=q((j-1)*3+3);
ri=[DR(1,k) DR(2,k)]';
rj=[DR(3,k) DR(4,k)]';
Cqr=Jac_rotula(Ri,teti,ri,Rj,tetj,rj);
Cq(3+2*(k-1)+1:3+2*(k-1)+2,3*(i-1)+1:3*(i-1)+3)=Cqr(1:2,1:3);
Cq(3+2*(k-1)+1:3+2*(k-1)+2,3*(j-1)+1:3*(j-1)+3)=Cqr(1:2,4:6);
end
%Prismáticos
for k=1:np;
i=P(1,k);
j=P(2,k);
Ri=[q((i-1)*3+1) q((i-1)*3+2)]';
teti=q((i-1)*3+3);
Rj=[q((j-1)*3+1) q((j-1)*3+2)]';
tetj=q((j-1)*3+3);
ri=[DP(1,k) DP(2,k)]';
hi=[DP(3,k) DP(4,k)]';
rj=[DP(5,k) DP(6,k)]';
Cqp=Jac_prism(Ri,teti,ri,hi,Rj,tetj,rj);
Cq(3+2*nr+2*(k-1)+1:3+2*nr+2*(k-1)+2,3*(i-1)+1:3*(i-1)+3)=Cqp(1:2,1:3);
Cq(3+2*nr+2*(k-1)+1:3+2*nr+2*(k-1)+2,3*(j-1)+1:3*(j-1)+3)=Cqp(1:2,4:6);
end
%EXTRAE EL JACOBIANO RESPECTO A LAS COORDENADAS DEPENDIENTES
Cqdep = zeros(m,m);
for i=1:length(DEP)
Cqdep(:,i) = Cq(:,DEP(i));
end
Función Jac_rotula.m
function Cq = Jac_rotula(Ri,teti,ri,Rj,tetj,rj)
Cq=zeros(2,6);
Cq(1:2,1:2)=eye(2);
Cq(1:2,4:5)=-eye(2);
103
Cq(1:2,3)=Ateti*ri;
Cq(1:2,6)=-Atetj*rj;
Función Jac_prism.m
function Cq=Jac_prism(Ri,teti,ri,hi,Rj,tetj,rj);
Cq=zeros(2,6);
Ai = Rotmat(teti);
Aj = Rotmat(tetj);
Cq(1,3)=1;
Cq(1,6)=-1;
Cq(2,1:2)=(Ai*hi)';
Cq(2,4:5)=-(Ai*hi)';
Cq(2,3)=(Ateti*hi)'*(Ri+Ai*ri-Rj-Aj*rj)+(Ai*hi)'*(Ateti*ri);
Cq(2,6)=(Ai*hi)'*(-Atetj*rj);
Función RotTet.m
function Atet = RotTet(x)
Atet(1,1)=-sin(x);
Atet(1,2)=-cos(x);
Atet(2,1)=cos(x);
Atet(2,2)=-sin(x);
FUNCIONES PROBLEMA VELOCIDAD
Función DtRestr.m
function Ct = DtRestr
global t w;
Ct = zeros(3*n,1);
Ct(3*n,1) = -w;
104
FUNCIONES PROBLEMA ACELERACIÓN
Función MatrizBq.m
function Bq = MatrizBq(q,v)
global n;
global t w;
Bq = MatBqPares(q,v);
Bq(3*n,:) = zeros(1,3*n);
Función MatBqPares.m
function Bq=MatrizBq(q,v)
global n nr np R P DR DP w;
global t;
m = 2*(nr+np)+3; %Número de restricciones
Bq=zeros(m,3*n);
%Barra fija
Bq(1,1)=0;
Bq(2,2)=0;
Bq(3,3)=0;
%Rótulas
for k=1:nr;
i=R(1,k);
j=R(2,k);
Ri=[q((i-1)*3+1) q((i-1)*3+2)]';
Vi=[v((i-1)*3+1) v((i-1)*3+2)]';
teti=q((i-1)*3+3);
dteti=v((i-1)*3+3);
Rj=[q((j-1)*3+1) q((j-1)*3+2)]';
Vj=[v((j-1)*3+1) v((j-1)*3+2)]';
tetj=q((j-1)*3+3);
dtetj=v((j-1)*3+3);
ri=[DR(1,k) DR(2,k)]';
rj=[DR(3,k) DR(4,k)]';
Bqr=Bq_rotula(Ri,teti,Vi,dteti,ri,Rj,tetj,Vj,dtetj,rj);
Bq(3+2*(k-1)+1:3+2*(k-1)+2,3*(i-1)+1:3*(i-1)+3)=Bqr(1:2,1:3);
Bq(3+2*(k-1)+1:3+2*(k-1)+2,3*(j-1)+1:3*(j-1)+3)=Bqr(1:2,4:6);
105
end
%Prismáticos
for k=1:np;
i=P(1,k);
j=P(2,k);
Ri=[q((i-1)*3+1) q((i-1)*3+2)]';
Vi=[v((i-1)*3+1) v((i-1)*3+2)]';
teti=q((i-1)*3+3);
dteti=v((i-1)*3+3);
Rj=[q((j-1)*3+1) q((j-1)*3+2)]';
Vj=[v((j-1)*3+1) v((j-1)*3+2)]';
tetj=q((j-1)*3+3);
dtetj=v((j-1)*3+3);
ri=[DP(1,k) DP(2,k)]';
hi=[DP(3,k) DP(4,k)]';
rj=[DP(5,k) DP(6,k)]';
Bqp=Bq_prism(Ri,teti,Vi,dteti,ri,hi,Rj,tetj,Vj,dtetj,rj);
Bq(3+2*nr+2*(k-1)+1:3+2*nr+2*(k-1)+2,3*(i-1)+1:3*(i-1)+3)=Bqp(1:2,1:3);
Bq(3+2*nr+2*(k-1)+1:3+2*nr+2*(k-1)+2,3*(j-1)+1:3*(j-1)+3)=Bqp(1:2,4:6);
end
Función Bq_rotula.m
function Bq=Bq_rotula(Ri,teti,Vi,dteti,ri,Rj,tetj,Vj,dtetj,rj)
Bq=zeros(2,6);
Ai = Rotmat(teti);
Aj = Rotmat(tetj);
Bq(1:2,1:2)=zeros(2);
Bq(1:2,4:5)=zeros(2);
Bq(1:2,3)=-dteti*Ai*ri;
Bq(1:2,6)=dtetj*Aj*rj;
Función Bq_prism.m
function Bq=Bq_prism(Ri,teti,Vi,dteti,ri,hi,Rj,tetj,Vj,dtetj,rj)
Bq=zeros(2,6);
Ai = Rotmat(teti);
Aj = Rotmat(tetj);
106
Bq(2,3)=(-Ateti*hi)'*Vi+((-Ai*hi)'*(Ri+Ai*ri-Rj-
Aj*rj)+2*(Ateti*hi)'*(Ateti*ri)+(Ai*hi)'*(Ai*ri))*dteti+(-
Ateti*hi)'*Vj+((Ateti*hi)'*(-Atetj*rj)+(Ai*hi)'*(Aj*rj))*dtetj;
Bq(2,6)=(Ateti*hi)'*(-Atetj*rj)*dteti+(Ai*hi)'*(Aj*rj)*dtetj;
Función DttRestr.m
function Ctt = DttRestr
global t w;
Ctt = zeros(3*n,1);
Ctt(3*n,1) = 0;
FUNCIONES CÁLCULO DINÁMICO.
Función MatrizMasa.m
function M = MatrizMasa(q)
global inercias;
M = zeros(3*n,3*n);
for i = 1:n,
teti = q(3*(i-1)+3,1);
mi = inercias(1,i);
I = inercias(2,i)+mi*(rG(1)^2+rG(2)^2); %Teorema de Steiner
M(3*(i-1)+1:3*(i-1)+2,3*(i-1)+1:3*(i-1)+2) = mi*eye(2);
M(3*(i-1)+3,3*(i-1)+3) = I;
M(3*(i-1)+1:3*(i-1)+2,3*(i-1)+3) = mi*Atet*rG;
M(3*(i-1)+3,3*(i-1)+1:3*(i-1)+2) = mi*rG'*Atet';
end
Función FCentrifuga.m
function Qv = FCentrifuga(q)
global inercias;
107
Qv = zeros(3*n,1);
for i = 1:n,
teti = q(3*(i-1)+3,1);
A = Rotmat(teti);
mi = inercias(1,i);
Qv(3*(i-1)+1:3*(i-1)+2,1) = mi*teti^2*A*rG;
end
Función FGravedad.m
function Qgrav = FGravedad(q)
global inercias g;
Qgrav = zeros(3*n,1);
for i = 1:n,
teti = q(3*(i-1)+3,1);
mi = inercias(1,i);
Qgrav(3*(i-1)+1:3*(i-1)+2,1) = mi*g;
Qgrav(3*(i-1)+3,1) = mi*rG'*Atet'*g;
end
REPRESENTACIÓN Y SIMULACIÓN DEL MECANISMO
Función PosicionPunto.m
function Rp = PosicionPunto(qi,rp)
Ri = qi(1:2,1);
Ai = Rotmat(qi(3,1));
Rp = Ri + Ai*rp;
108
9.2. ANEXO B
En este anexo se muestra las gráficas tipos que obtenemos con la simulación del programa. En este
caso la única combinación mostrada es la que hemos utilizado para realizar los cálculos, es decir, la
que obtiene mayor ahorro económico.
0.8
0.6
0.4
0.2
-0.2
-0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2
FIGURA 9.1. POSICIÓN MECANISMO EN DIFERENTES INSTANTES DE TIEMPO
1.5
0.5
-0.5
-1
-1.5
-2
-2 -1 0 1 2 3 4
FIGURA 9.2. SIMULACIÓN MECANISMO
109
Par Motor
32
30
28
26
Par (Nm)
24
22
20
18
16
0 0.5 1 1.5 2 2.5 3 3.5 4
Tiempo (s)
FIGURA 9.3. PAR MOTOR RESPECTO AL TIEMPO
Potencia Instantánea
260
240
220
Potencia (W)
200
180
160
140
120
0 0.5 1 1.5 2 2.5 3 3.5 4
Tiempo (s)
FIGURA 9.4. POTENCIA GENERADA RESPECTO AL TIEMPO
110
Trabajo Motor
700
600
500
Trabajo (J)
400
300
200
100
0
0 0.5 1 1.5 2 2.5 3 3.5 4
Tiempo (s)
FIGURA 9.5. TRABAJO MOTOR RESPECTO AL TIEMPO
111
112

Proyecto Eliptica 2

Cargado por

Copyright:

Formatos disponibles

Proyecto Eliptica 2

Cargado por

Información del documento

Título original

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Proyecto Eliptica 2

Cargado por

Copyright:

Formatos disponibles

Proyecto Fin de Carrera

Autor: Rubén Rosado Ávila

Tutor: José Ángel Acosta Rodríguez

DEPARTAMENTO DE INGENIERÍA DE SISTEMAS Y AUTOMÁTICA

ESCUELA TÉCNICA SUPERIOR DE INGENIERÍA

A mi padre, mi madre y mi hermano,

por darme todo su apoyo y confianza.

A los que no se encuentran entre nosotros,

siempre estaréis presentes.

A aquellos que me ayudaron a seguir hacia adelante.

A aquellos profesores que siempre me apoyaron.

A mi tutor José Ángel A. R.

por su orientación, seguimiento y supervisión del mismo.

CAPÍTULO 2: TEORÍA DE MECANISMOS ……………………………………………………….……………….PAG. 9

2.1 CINEMÁTICA NUMÉRICA PLANA………………………………………………………..…………PAG. 9

2.1.2 COORDENADAS CARTESIANAS O DE REFERENCIA….………………………..PAG. 9

2.1.3 RESTRICCIONES IMPUESTAS POR LOS PARES CINEMÁTICOS ………....PAG. 11

2.1.3.1 Par de revolución…………………………………………………………..PAG. 11

2.1.3.2 Par prismático………………………………………..……………………PAG. 12

2.1.4 PROBLEMA DE POSICIÓN…………………………………………………………….PAG. 13

2.1.5 SIMULACIÓN CINEMÁTICA…………………………………….…………………….PAG. 13

2.1.5.1 Problema de posición………………………………………………….PAG. 14

2.1.5.2 Problema de velocidad…………………………..…………………..PAG. 14

2.1.5.3 Problema de aceleración………………………..…………………..PAG. 15

2.1.5.4 Matriz jacobiana y matriz Bq…………………..…………………..PAG. 15

2.2 DINÁMICA NUMÉRICA PLANA……………………………………………………………………PAG. 17

2.2.2 DINÁMICA CON RESTRICCIONES…………………………………………………PAG. 17

2.2.2.1 Principio de los Trabajos Virtuales……………………………..PAG. 18

2.2.2.2 Método de los multiplicadores de Lagrange..…………….PAG. 18

2.3 PROBLEMA INVERSO………………………………………………………………………………...PAG. 23

CAPÍTULO 3: MODELADO Y DESCRIPCIÓN DEL SISTEMA BICICLETA ELÍPTICA……….………..PAG. 25

3.1 DEFINICIÓN DEL MECANISMO……………………………………………………………………PAG. 25

3.2 DATOS BICICLETA ELÍPTICA………..………………………………………………………………PAG. 27

3.4 ESTUDIO CINEMÁTICO……………………………………………………………………………..PAG. 32

3.4.1 PROBLEMA DE POSICIÓN………………………………………………………….PAG. 35

3.4.2 PROBLEMA DE VELOCIDAD……………………………………………………….PAG. 41

3.4.3 PROBLEMA DE ACELERACIÓN……………………………………………………PAG. 41

3.5 ESTUDIO DINÁMICO…………………………………………………………………………………PAG. 44

CAPÍTULO 4: DESCRIPCIÓN DEL SOFTWARE Y USO……………..…………………………………………PAG. 47

4.1 FUNCIONES PROBLEMA POSICIÓN………………………………..…………………………PAG. 47

4.2 FUNCIONES PROBLEMA VELOCIDAD……………………………..…………………………PAG. 60

4.3 FUNCIONES PROBLEMA ACELERACIÓN…………………………………………………….PAG. 60

4.4 FUNCIONES CÁLCULO DINÁMICO…………………………………..………………………..PAG. 66

4.5 REPRESENTACIÓN Y SIMULACIÓN DEL MECANISMO……..………………………...PAG. 68

CAPÍTULO 5: ESTUDIO ECONÓMICO…………..…………………………………………………………………PAG. 75

5.2 ELECCIÓN DEL GENERADOR……………………………………………………………………….PAG. 75

5.3 ELECCIÓN DEL CABLE..……………………………………………………………………………….PAG. 80

5.4 ELECCIÓN DE LA BATERÍA…….……………………………………………………………………PAG. 84

9.1 ANEXO A………….……………..………………………………………………………………………..PAG. 95

9.2 ANEXO B…………………………..…………………………………………………………………….PAG. 109

Comparativa Bicicletas elípticas vs Bicis estáticas o spinning.

Diferencias entre bicicletas elípticas y las estáticas o de spinning:

Similitudes entre las bicicletas elípticas y las estáticas:

 Movimiento de muy bajo impacto sobre las articulaciones.

Ventajas de las bicicletas elípticas:

Desventajas de las bicicletas elípticas:

A todo lo comentado anteriormente, se le añade además la motivación extra de aprovechar

Los objetivos que desea cumplir este trabajo son:

1. Modelado de la bicicleta elíptica.