Influencia de las variaciones en los puntos de inserción de los músculos de la cadera en la distribución de densidades del fémur

i Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería de las Tecnologías Industriales Influencia de las variaciones en los puntos de inserción de los músculos de la cadera en la distribución de densidades del fémur Autor: Jesús Serrano Marcos Tutor: Joaquín Ojeda Granja Dep. Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, septiembre de 2018 i Trabajo Fin de Grado Grado en Ingeniería de las Tecnologías Industriales Influencia de las variaciones en los puntos de inserción de los músculos de la cadera en la distribución de densidades del fémur Autor: Jesús Serrano Marcos Tutor: Joaquín Ojeda Granja Profesor Contratado Doctor Interino Dep. Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, septiembre de 2018 iii Trabajo Fin de Grado: Influencia de las variaciones en los puntos de inserción de los músculos de la cadera en la distribución de densidades del fémur Autor: Jesús Serrano Marcos Tutor: Joaquín Ojeda Granja El tribunal nombrado para juzgar el Trabajo arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, septiembre de 2018 El Secretario del Tribunal v A mi familia A mis maestros vii Agradecimientos Para comenzar, las primeras palabras van dedicadas a mis padres, Jesús y Amparo. A ellos les debo todas las oportunidades que he tenido en la vida y la educación y formación que he tenido la fortuna de recibir desde pequeño y hasta el día de hoy. También estoy muy agradecido por todo a mi hermana Amparo, que siempre me ha allanado el camino y servido como guía y ejemplo a seguir. Me gustaría mostrar mi gratitud a mi tutor, Joaquín Ojeda, sin cuyas directrices y conocimientos me habría sido imposible realizar este trabajo. Mencionar también a Ezequiel Martín y agradecer su ayuda y paciencia en los primeros pasos del proyecto, puesto que sin él no habría podido ni comenzar. Por último, gracias a todos los amigos y compañeros que han estado a mi lado durante la carrera y especialmente a Marta por compartir conmigo esta travesía. Jesús Serrano Marcos Sevilla, septiembre de 2018 Figura 3.2. Parámetros de entrada para generar los archivos con las coordenadas del modelo. 41 Figura 3.3. Obtención de la distancia entre dos nodos. 42 Figura 3.4. Dirección perpendicular a la superficie, que será descartada en el análisis. 42 Figura 3.5. Imprecisión introducida al proyectar los puntos de ensayo sobre nodos de la superficie del modelo. 43 Figura 3.6. Caso nominal y ensayos considerados a priori para un punto de inserción. 43 Figura 3.7. Planos por los que se han realizado los cortes para ilustrar los resultados. 45 Figura 3.8. Comparación de los canales medulares obtenidos en los casos mostrados del psoas. 46 Figura 3.9. Resultados del caso nominal del psoas. Distribución de tensiones en el instante inicial, cuando la densidad es homogénea (a); y distribución final de densidades (b). 47 Figura 3.10. Resultados del ensayo en la dirección Y en sentido negativo del psoas. Distribución de tensiones en el instante inicial, cuando la densidad es homogénea (a); y distribución final de densidades (b). 47 Figura 3.11. Comparación de los canales medulares obtenidos en los casos mostrados del glúteo medio. 48 Figura 3.12. . Resultados del caso nominal del glúteo medio. Distribución de tensiones en el instante inicial, cuando la densidad es homogénea (a); y distribución final de densidades (b). 49 Figura 3.13. Resultados del ensayo en la dirección Y en sentido negativo del glúteo medio. Distribución de tensiones en el instante inicial, cuando la densidad es homogénea (a); y distribución final de densidades (b). 49 Figura 3.14. Comparación de los canales medulares obtenidos en los casos mostrados del aductor largo. 50 Figura 3.15. . Resultados del caso nominal del aductor largo. Distribución de tensiones en el instante inicial, cuando la densidad es homogénea (a); y distribución final de densidades (b). 51 Figura 3.16. Resultados del ensayo en la dirección X en sentido negativo del aductor largo. Distribución de tensiones en el instante inicial, cuando la densidad es homogénea (a); y distribución final de densidades (b). 51 Figura 3.17. Resultados del ensayo con las inserciones de la resonancia magnética. Distribución de tensiones en el instante inicial, cuando la densidad es homogénea (a); y distribución final de densidades (b). 52 Tabla 1. Músculos del modelo y fuerza muscular máxima en el ciclo de marcha (caso nominal). 13 Tabla 2. Valores de las constantes del modelo de remodelación empleado en las simulaciones. 30 Tabla 3. Fuerza máxima del psoas. 46 Tabla 4. Densidad y masa finales del fémur (ensayos del psoas). 46 Tabla 5. Fuerzas del glúteo medio cuando la del actuador 1 es máxima. 48 Tabla 6. Densidad y masa finales del fémur (ensayos del glúteo medio). 48 Tabla 7. Fuerza máxima del aductor largo. 50 Tabla 8. Densidad y masa finales del fémur (ensayos del aductor largo). 50 Tabla 9. Densidad y masa finales del fémur (ensayos nominal y de resonancia). 52 xv 1 1 INTRODUCCIÓN La bioingeniería se puede definir como la disciplina que estudia las aplicaciones de la ingeniería a la medicina o a la biología. Ésta se basa en aplicar conceptos y métodos físico-matemáticos y las metodologías analíticas y sintéticas de la ingeniería a organismos vivos. En una era en la que los progresos tecnológicos e ingenieriles se suceden a un ritmo vertiginoso, cada vez se integran en mayor medida los trabajos de distintos campos de la ciencia a fin de profundizar más en los conocimientos de los mismos desde distintas perspectivas. La bioingeniería no es una excepción, aunando los avances tanto en ingeniería, como en biología y medicina. Desde que comenzara su desarrollo hacia la mitad del siglo pasado, hasta el día de hoy, son innumerables los estudios que se han desarrollado y que se siguen llevando a cabo en esta disciplina. Esta ciencia multidisciplinar abarca un amplio abanico de ramas, entre las que se encuentra la que atañe a este proyecto: la biomecánica. La biomecánica se erige como uno de los campos de la mecánica con mayor potencial de desarrollo en la actualidad. Una definición que se podría dar de la biomecánica es: “el estudio de los modelos, fenómenos y leyes que sean relevantes en el movimiento de un sistema biológico a través de los métodos de la mecánica” [1]. Esta rama de la bioingeniería tiene múltiples aplicaciones de diversa índole, desde la evaluación y tratamiento de distintas patologías que afectan al aparato locomotor hasta material y equipamiento deportivo, pasando la ergonomía de distintos elementos cotidianos para mejorar el confort en determinadas actividades. En el año 1974, H. Hatze definió la biomecánica como la “ciencia que estudia el movimiento de los seres vivos usando conceptos de la mecánica clásica” [2], [3]. A día de hoy la biomecánica no limita su ámbito de investigación únicamente al movimiento, por lo que en esa definición sería adecuado sustituir la palabra movimiento por comportamiento. De hecho, en relación al párrafo anterior, el presente trabajo, si bien se ha nutrido del análisis del movimiento de un sujeto durante un ciclo de marcha; tiene por objetivo final el análisis del comportamiento del fémur (de su distribución de densidades) ante las solicitaciones en él generadas durante la marcha y el consecuente estado tensional. Para obtener la relación entre ambos fenómenos – movimiento y comportamiento estructural del hueso –, se hace uso de la teoría de la remodelación ósea, según la cual la estructura del hueso se adapta en respuesta a las alteraciones de carga que sufre. Introducción 2 1.1 Estado del arte El primer autor en proponer una relación entre la estructura del hueso y las tensiones que éste soporta fue Wolff, en 1892 [4]. Muchos son los autores que posteriormente han investigado el fenómeno de la remodelación ósea, desarrollando diferentes modelos y obteniendo resultados que permiten conocer el comportamiento del fémur a fin de desarrollar aplicaciones biomédicas. Aunque algunos ejemplos simplificados han hecho uso de métodos analíticos, los métodos de aproximación numérica son los que se emplean usualmente a la hora de cuantificar las propiedades mecánicas del fémur. Esto es debido a la complejidad de la geometría del hueso, la descripción de las propiedades del material, la distribución y variación temporal de las condiciones de contorno y la solicitación mecánica. Las aplicaciones más comunes para el desarrollo de la teoría de remodelación ósea se basan en el método de elementos finitos, el cual emplea modelos computacionales descritos mediante la subdivisión de formas complejas en una colección de elementos pequeños, pero cuantificables, que constituyen una “malla” [5]. Beaupré et al. desarrollaron un estudio [6] en el que se llevó a cabo la simulación computacional del modelo de remodelación que proponían, con elementos finitos bidimensionales. Se concluyó que una reducción de la solicitación en la historia de carga tenía como consecuencia la atrofia del hueso en determinadas zonas (reabsorción) y que el restaurar posteriormente la actividad normal en la historia de carga resultaba en una distribución diferente a la que se obtuvo calculando para un estado estático de cargas. Con un enfoque similar, Doblaré y García, hicieron también uso métodos computacionales. Desarrollaron un modelo de remodelación ósea [7] basado en la mecánica del daño, obteniendo resultados muy similares a los experimentales al aplicarlo en la evolución del fémur en el extremo proximal. Behrens et al. [8] llevaron a cabo un estudio numérico con el objetivo de determinar la utilidad de considerar únicamente un estado estático de cargas en los cálculos con elementos finitos. La conclusión que sacaron fue que un único estado estático de cargas es insuficiente para representar todo el ciclo de marcha de un sujeto y conocer la distribución de densidades del fémur. Esto puso de manifiesto la importancia de las condiciones de contorno, así como el tener en cuenta las cargas que se producen a lo largo de las distintas fases del ciclo de marcha. En la tesis doctoral de J. Ojeda [9] se propuso el modelo de remodelación ósea que se ha empleado en este trabajo, como se verá más adelante. En el estudio se llevó a cabo un análisis de sensibilidad de los principales parámetros que afectan al proceso remodelatorio y se obtuvieron resultados similares en la parte proximal a los obtenidos en la resonancia realizada al sujeto de estudio. También se tuvo en cuenta en este estudio la importancia de la dinámica de la marcha, frente a los estudios que obtienen resultados preliminares aplicando cargas estáticas. Las aplicaciones biomédicas como puede ser la implantación de prótesis [10] también suponen un recurrente objeto de estudio. En [8] se estudia el fenómeno de remodelación ósea tras una artroplastia total de cadera; y en [11] se aborda de forma numérica el diseño de fijaciones de fracturas del fémur proximal. Para poder aplicar los modelos de remodelación y llevar a cabo simulaciones mediante software de elementos finitos, se plantea el problema de resolver tanto los análisis cinemático y dinámico inversos de un ciclo de marcha humana. La solución a dicho problema la proporciona la 3 dinámica de los sistemas multicuerpo [12]. Un sistema multicuerpo es la modelización de un conjunto de sólidos rígidos o flexibles y las uniones entre los mismos. Dicho conjunto constituye un sistema físico cuya dinámica se puede estudiar en base a sus condiciones de contorno. Es necesario precisar que los modelos de sistemas multicuerpo que se utilizan para estimar las condiciones de contorno y solicitaciones mecánicas del fémur no son idénticos a los del sujeto estudiado o el modelo de elementos finitos. Por ejemplo, el programa de uso libre OpenSim [13], realiza un escalado para que el modelo se asemeje en la mayor medida posible al sujeto de estudio. No obstante, hay que tener en cuenta que el escalado que se lleva a cabo está basado en aproximaciones matemáticas y en estudios realizados a cadáveres, por lo que el resultado que se obtiene aún sigue siendo muy estandarizado. A este respecto, se han llevado a cabo estudios [14] en los que se realizan modificaciones en los puntos de origen e inserción de los músculos, cuyas posiciones se determinan mediante resonancias magnéticas, a fin de asemejarse al sistema real. Sin embargo, si bien se conoce que el modelo no reproduce de forma exacta el comportamiento del hueso del sujeto particular, no se ha hallado en la literatura un estudio que ponga de manifiesto en qué medida las diferencias o imprecisiones de la ubicación de los puntos de inserción de los músculos afectan a los resultados de la distribución de densidades. Es decir, resultaría de interés establecer la importancia de la precisión del modelo en el fenómeno de remodelación ósea. 1.2 Objetivos Según lo establecido en el párrafo anterior, la motivación del presente Trabajo de Fin de Grado será estudiar las consecuencias que tienen las variaciones del modelo en los resultados de distribución de densidad obtenidos con la aplicación de un estado de cargas. De esta forma, el objetivo será llevar a cabo un análisis de sensibilidad en el que se varíen los puntos de inserción en el fémur de los músculos de la cadera que tengan mayor influencia en el estado tensional; y analizar, cualitativa y cuantitativamente, las consecuencias de esta variación sobre la estructura del hueso, a nivel de densidad. Asimismo, se busca realizar una discusión sobre qué músculos de los analizados tienen mayor influencia en el proceso de remodelación ósea. De esta forma, se podrá determinar la influencia que pequeñas variaciones en las características del modelo empleado respecto al sujeto real tienen sobre los resultados obtenidos para dicho fenómeno. Introducción 4 5 2 METODOLOGÍA En el apartado 1.2 se ha expuesto que el fin último del estudio es comprobar cómo y en qué medida afectan los cambios de posición de los puntos de inserción en el fémur a la distribución de densidades de éste. Asimismo, se remarcó que todo el análisis se llevará a cabo mediante simulaciones con un programa de elementos finitos: Abaqus®. Con ello en mente, es procedente definir una metodología de trabajo que defina el camino a seguir para llegar a obtener los datos de entrada que habrá que proporcionarle al software: geometría del hueso a analizar y estado de carga al que se encuentra sometido durante la simulación. Hay varios archivos de los que se habrá de partir para llegar a los datos de entrada que permitirán llevar a cabo las simulaciones y el análisis de sensibilidad. En primer lugar, un modelo de elementos finitos del fémur, que proporcionará tanto la geometría como el mallado de éste [9]. Por otra parte, tablas de datos en forma de hojas de cálculo que contienen información del modelo del sistema músculo-esquelético (SME), el módulo de las fuerzas musculares (FM) y las componentes de las fuerzas de reacción articular (FRA), todo ello obtenido de OpenSim y parte del trabajo desarrollado en [14]. La Figura 2.1 presenta, pues, un diagrama que ilustra la metodología que se seguirá para obtener el estado de carga que permitirá implementar las simulaciones. En ella aparecen en verde los pasos que requieren del empleo de OpenSim y en rojo los procedimientos que se abordan con Abaqus®. En ambos casos, se emplean rutinas de programación desarrolladas en Matlab® y hojas de cálculo de Microsoft Excel® para trabajar con los datos de forma ordenada. Las sucesivas secciones de este capítulo están dedicadas al desarrollo en detalle de cada uno los bloques de esta metodología. Metodología 6 Figura 2.1. Diagrama de bloques de los pasos a seguir en la metodología. 2.1 Trabajo previo En este apartado se va a abordar todo lo relacionado con los archivos que fueron facilitados como punto de partida para llevar a cabo este estudio. Como se destacó en el texto introductorio del capítulo, estos archivos se pueden agrupar en: • Modelo de elementos finitos • Hojas de Excel® con datos procesados del registro de la marcha en OpenSim A continuación, se procede a ahondar en dicho trabajo previo a partir del cual se desarrollarán el resto de pasos descritos en la metodología. 2.1.1 Elementos finitos El modelo de elementos finitos del fémur proporciona la geometría y el mallado del hueso que se va a analizar. Este modelo se llevó a cabo en el desarrollo de la tesis de J. Ojeda [9]. En ella se estudió el fémur de un hombre adulto sin patologías en la marcha. A dicho sujeto se le realizó una resonancia del fémur mediante tomografía computacional, de modo que se pudiera reconstruir la geometría del hueso con precisión. 7 El modelo se empleará para obtener las tensiones y finalmente la distribución de densidades del fémur, aplicándole los estados de cargas a los que está sometido durante el ciclo de marcha. De esta forma, en [9] se desarrolló el modelo de elementos finitos a partir de los datos de la tomografía computacional del fémur derecho del sujeto. Dichos datos se emplearon para reconstruir la geometría en 3D y posteriormente, se realizó el mallado pertinente. Figura 2.2. Modelo de elementos finitos del fémur de estudio [9]. Tras ello, se llegó al modelo que se muestra en la Figura 2.2, que consta de 339168 elementos tetraédricos (4 nodos por elemento, C3D4 de la librería de Abaqus®) y 64757 nodos. Este modelo constituirá la base del estudio, pues sobre él se realizarán las simulaciones al aplicarle los estados de carga que se verán más adelante. 2.1.2 OpenSim y registro de la marcha Teniendo ya la geometría del fémur, para lanzar las simulaciones con el modelo de remodelación ósea, resta establecer los distintos estados de carga a los que éste se ve sometido durante el tiempo que dura el ciclo de marcha. A ello se llega siguiendo la metodología del esquema de la Figura 2.1; partiendo de un registro de la marcha de un sujeto y del tratamiento de dichos datos en OpenSim, tal y como se explica en las líneas siguientes de esta sección. Uno de los pilares fundamentales para llevar a cabo el registro de la marcha, así como para obtener los estados de carga del fémur, es el modelo del sistema músculo-esquelético que se va a emplear. En este caso, el trabajo previo [14] se basó en uno de los modelos proporcionados por OpenSim: “gait2392_simbody.osim” (Figura 2.3). La toma de datos se llevó a cabo empleando un pasillo de marcha consistente en dos plataformas de fuerza para cuantificar las fuerzas y momentos de reacción con el suelo. Asimismo, se contó con seis cámaras de visión infrarroja para la captación de los marcadores ubicados sobre el sujeto de análisis. Se empleó un protocolo de 39 marcadores [15]. Si bien, como se ha explicado en un párrafo anterior, se usó un modelo proporcionado por OpenSim, para que los resultados se ajusten fielmente a la realidad, fue preciso llevar a cabo un escalado de dicho modelo para que se ajustase a las medidas del sujeto, en base al protocolo de marcadores utilizado. Metodología 14 2.2.2.1 Componentes de las FM Del trabajo previo desarrollado en el apartado 2.1 se tiene la evolución temporal del módulo de las fuerzas musculares (Figura 2.7). Sin embargo, para llevar a cabo los ensayos en el modelo tridimensional de elementos finitos, es preciso conocer las tres componentes de las fuerzas musculares, puesto que la dirección en la que se producen estas solicitaciones tiene un papel fundamental en el estado tensional. Figura 2.7. Evolución temporal del módulo de las fuerzas musculares (Glúteo medio 1). El objetivo es, por tanto, obtener tres archivos con la evolución temporal de cada una de las componentes de las fuerzas musculares a partir del archivo que contiene los módulos y de la geometría del modelo de sistemas multicuerpo. Para ello se ha hecho uso de tres rutinas de programación en Matlab®. La primera de las rutinas simplemente lee los datos de la hoja de Excel® en la que se encuentran los módulos de las fuerzas musculares y los almacena como una variable de Matlab®, para así poder trabajar con ellos. A continuación, la segunda rutina carga los datos de la geometría del modelo de OpenSim y a partir de ésta resuelve el problema cinemático de cada uno de los músculos, de forma que se obtiene la dirección en la que inciden todas las fuerzas musculares en el fémur. A partir de los vectores unitarios que definen las direcciones; y de los módulos, se obtienen las tres componentes de las fuerzas musculares de cada músculo para el caso que se estudia (Figura 2.8). 15 Figura 2.8. Evolución temporal de las componentes de la fuerza muscular (Glúteo medio 1). La tercera y última rutina tiene una función similar a la primera: escribir los valores obtenidos como tablas en hojas de cálculo de Excel®, para facilitar el posterior uso de éstos. Así, los valores de las componentes de las fuerzas musculares se encuentran en disposición de ser introducidos más adelante como datos de entrada en Abaqus® para llevar a cabo las simulaciones que permitirán conocer la distribución de densidades. Sin embargo, antes, hace falta conocer los nodos del modelo sobre los que se aplicarán los valores obtenidos. Metodología 16 2.2.2.2 Nodos de aplicación de las FM Del modelo escalado de OpenSim se obtiene una hoja de cálculo que contiene la información de las coordenadas de los puntos de origen e inserción de todos los músculos que intervienen en el caso estudiado. Este archivo de Microsoft Excel® muestra las tres componentes del vector posición de los puntos de interés de cada músculo, especificando asimismo el sistema de referencia en el que se proporcionan. De este modo, se tienen las coordenadas de los puntos del modelo de OpenSim en los que actúan los músculos, en el sistema de referencia local del fémur de OpenSim. Como se recalcó al comienzo del capítulo, al pasar de trabajar con OpenSim a hacerlo con Abaqus, es preciso llevar a cabo un cambio de coordenadas. Tras ello, se tendrán las coordenadas de los puntos del fémur en los que habrá aplicadas fuerzas musculares, en el sistema de referencia local del fémur de Abaqus. Resta, finalmente, obtener los nodos del modelo de elementos finitos sobre los que se aplicarán las fuerzas musculares. Para ello, conocidas las posiciones de todos los nodos que conforman el modelo (se obtiene del archivo .inp) y las de los puntos de inserción en el fémur, se desarrolló una rutina en Matlab® que devuelve los nodos más próximos a los puntos introducidos, como ilustra la Figura 2.9. En ella se muestra cómo, para el punto de inserción cuya posición viene dada por 𝒓𝑖𝑛𝑠𝑒𝑟𝑐𝑖ó𝑛 , el nodo seleccionado sería el más cercano de los del modelo. En la figura, el “nodo 1”. Figura 2.9. Croquis de la selección del nodo de inserción de un músculo. 17 Cabe destacar que, si bien el modelo de OpenSim está escalado para asemejar sus dimensiones al caso que se estudia (debe tener las dimensiones escaladas a las del modelo de elementos finitos), hay detalles de la geometría, como es la curvatura del hueso en la diáfisis, que no coinciden en ambos modelos. Esto provoca que algunos de los nodos obtenidos con el método recién descrito no se encuentren sobre la superficie del fémur e incluso, en los casos más críticos, se encuentren relativamente alejados de la posición que cualitativamente deberían ocupar. Debido a lo anterior, una vez obtenida la lista de nodos mediante la rutina de Matlab®, es necesario llevar a cabo un ajuste de las posiciones: proyectar los nodos sobre la superficie y reubicar aquellos cuya posición no coincida con el punto de inserción real. Todo este ajuste se lleva a cabo, eso sí, manteniendo la mayor fidelidad al vector que se obtuvo inicialmente de OpenSim. Figura 2.10. Diagrama de flujos para la selección del nodo de inserción de un músculo. 2.2.3 Fuerzas de reacción articular En un sistema como el abordado en este trabajo, sobre el que los músculos representan una fuerza externa, el concepto de fuerzas de reacción articular es definido como el conjunto de las fuerzas y momentos netos generados por el contacto hueso-hueso entre segmentos corporales adyacentes. Estas cargas son fruto tanto de la fuerza de reacción con el suelo como de la acción de los músculos y son soportadas como fuerzas internas por la estructura de la articulación en cuestión. El estudio está más orientado a analizar los efectos de los cambios de los puntos de inserción en la parte proximal del fémur, por lo que, para lo que se pretende abarcar, las fuerzas de reacción Metodología 18 articular que interesan son las de la articulación de la cadera. Por su parte, las fuerzas de reacción articular en la rodilla no se modelarán, sino que se incluirán en nodos de la epífisis distal del hueso condiciones de contorno en desplazamientos. 2.2.3.1 Articulación de la cadera La cadera es una enartrosis compuesta por el cótilo y la cabeza femoral. La cabeza femoral constituye el componente convexo de la articulación. Forma dos tercios de esfera y está aplanada proximalmente. El cótilo o cavidad cotiloidea es el componente cóncavo de la articulación. La cabeza femoral se encuentra encajada en el cótilo, estando retenida por una banda de tejido conectivo, llamado ligamento redondo y una cápsula articular fibrosa a modo de manguito. Figura 2.11. Fuerzas de reacción articular de la cadera sobre la cabeza del fémur. Adaptada de [16]. Cabe destacar que la superficie articular anatómica es más extensa que la superficie de contacto entre la cabeza femoral y la cavidad cotiloidea. Asimismo, dentro de la superficie de contacto existe una zona de cartílago donde tiene lugar la transmisión de la fuerza resultante de la cadera sobre la cabeza femoral, conocida como superficie articular portante (Figura 2.11). Si bien el análisis de sensibilidad se va a llevar a cabo modificando únicamente las fuerzas musculares, las fuerzas de reacción articular de la cadera no dejan de ser importantes en el estado de cargas. Su magnitud es mayor que la de las fuerzas musculares, por lo que su aportación al estado tensional del fémur habrá de ser tenida en cuenta. Es por ello que se lleva a cabo el modelado de éstas. Como recoge la literatura [16], la distribución de cargas sobre la cabeza del fémur tiene la forma que se muestra en la Figura 2.11, con la resultante aplicada en su centro. Sin embargo, es 19 razonable modelar esta solicitación con una distribución en forma de paraboloide de revolución (Figura 2.12), siempre que la resultante sea la misma. 2.2.3.2 Distribución: paraboloide de revolución Como se expuso en el apartado referente al trabajo previo, se cuenta con el vector resultante de las fuerzas de reacción articular. De la misma forma, se conoce también el nodo de aplicación dicho vector, ya que al hallarse el punto de aplicación en el centro de la cabeza femoral, no será otro que el nodo que se estableció como origen de coordenadas del sistema de referencia local del fémur en el apartado de los sistemas de referencia. Figura 2.12. Distribución de la fuerza de reacción articular como un paraboloide de revolución. Adaptada de [16]. A continuación, dicho nodo de aplicación se proyecta sobre la superficie del modelo según la dirección de la resultante de las fuerzas articulares, obteniéndose un nodo que, junto con el origen del sistema de referencia local, definirán el eje de revolución del paraboloide. Para modelar la distribución parabólica de cargas, el siguiente paso consistirá en definir un nuevo sistema de referencia del paraboloide, en el cual la dirección Z coincida con el eje de revolución que se acaba de definir, de vector unitario 𝒖𝑟𝑒𝑣 (𝒖𝑟𝑒𝑣𝑓 en coordenadas del fémur en Abaqus® y 𝒖𝑟𝑒𝑣𝑝=[ 0 0 1 ]𝑇 en coordenadas del paraboloide). Dicho sistema de referencia se obtiene mediante sendos giros en torno a los ejes X e Y. Se define el sistema de ecuaciones que permitirá obtener los ángulos y matrices de giro que efectuarán el cambio de un sistema de coordenadas a otro. Denominando 𝜃 al ángulo girado en torno al eje X; y 𝜑 al ángulo girado en torno al eje Y, las matrices de giro vienen definidas por: Metodología 20 𝑹𝑥(𝜃)=[ 1 0 0 0 cos𝜃 −sin𝜃 0sin𝜃 cos𝜃] (2.3) 𝑹𝑦(𝜑)=[cos𝜑 0 sin𝜑 0 1 0 −sin𝜑 0 cos𝜑] (2.4) 𝑹𝑥𝑦(𝜃,𝜑)=𝑹𝑥(𝜃)𝑹𝑦(𝜑)=[cos𝜃 0 sin𝜑 sin𝜃sin𝜑 cos𝜃 −sin𝜃cos𝜑 cos𝜃sin𝜑sin𝜃 cos𝜃cos𝜑 ] (2.5) Quedando el sistema de ecuaciones finalmente: 𝑹𝑥𝑦[00 1 ]=𝒖𝑟𝑒𝑣𝑓 (2.6) De donde se obtienen los ángulos girados 𝜃 y 𝜑, así como las matrices de giro. Conocidos estos datos, se cambia de sistema de coordenadas al del paraboloide para proceder a definir las fuerzas. Para poder dar un valor de fuerza aplicada a cada nodo, en primer lugar se establece un diámetro máximo del paraboloide, que se producirá en el corte con la superficie de la cabeza del fémur y delimitará el área de la zona portante de carga. El valor escogido para dicho diámetro fue de 3 cm. Con el diámetro máximo elegido, el siguiente paso es realizar en Abaqus® una selección de los nodos que pertenecen a la cabeza del fémur y que, partiendo del nodo del eje de revolución sobre la superficie, incluya a todos los nodos que se encuentren en un diámetro de al menos 3 cm (se escogen algunos más alejados, para asegurar que ningún nodo queda fuera de la selección realizada). Para definir las cargas, se asumirá la superficie del fémur como un plano horizontal en el sistema de referencia del paraboloide, de modo que solo se tendrán en consideración para la obtención del valor de las cargas las componentes X e Y de las posiciones de los nodos. A partir de dichas posiciones, se define un paraboloide cuyo valor máximo en la componente Z sea unitario. La ecuación empleada es, por consiguiente, la (2.7). 𝑧𝑝=1−( 𝑥𝑝2 𝑟2+𝑦𝑝2 𝑟2 ) (2.7) Donde, como ya se ha establecido en esta sección, 2𝑟=3 𝑐𝑚. Este paraboloide unitario, se ha definido como base para obtener la distribución de cargas final que se empleará en el modelo. Para ello solo resta escalarlo para que la fuerza resultante de esta distribución se corresponda con el vector que se obtuvo de OpenSim. Se sabe que el vector resultante de una distribución parabólica continua como la que se ha 21 definido, se hallaría aplicado en la dirección del eje del paraboloide y sobre el punto de corte de éste con el plano horizontal. Su módulo sería la integral del volumen del mismo. Si bien la distribución no es continua, puesto que se trabaja con nodos discretos, es razonable asumir que el punto de aplicación de la resultante se hallará sobre el nodo que definimos como origen del sistema de referencia del paraboloide para tal efecto. Por otro lado, al tener todas las cargas la misma dirección y sentido, el módulo de la resultante será igual a la suma de los módulos de todos los vectores fuerza que componen la distribución. Esto es: 𝐹𝑢𝑛𝑖𝑡𝑎𝑟𝑖𝑜=∑𝑓𝑖 𝑛𝑝 𝑖=1 (2.8) Donde 𝐹𝑢𝑛𝑖𝑡𝑎𝑟𝑖𝑜 es el módulo de la resultante de las fuerzas del paraboloide unitario, 𝑓𝑖 son los módulos de cada una de las fuerzas que componen la distribución unitaria y 𝑛𝑝 es el número de nodos que la componen. Por tanto, siendo 𝐹𝑐𝑎𝑑𝑒𝑟𝑎 el módulo de la fuerza de reacción articular en la cadera, obtenido de OpenSim, el factor por el que habrá de ser multiplicado el paraboloide unitario para representar la distribución de cargas modelada no es otro que 𝑘, que se define como: 𝑘=𝐹𝑐𝑎𝑑𝑒𝑟𝑎 𝐹𝑢𝑛𝑖𝑡𝑎𝑟𝑖𝑜 (2.9) Tras aplicar el factor de semejanza a la distribución de fuerzas unitaria obtenida, se tiene la distribución de fuerzas que se aplicará, en las coordenadas del paraboloide, que se ilustra en la Figura 2.13. Así, restaría únicamente realizar el cambio de coordenadas. Metodología 22 Figura 2.13. Distribución de las fuerzas de reacción articular (en las coordenadas del paraboloide). Finalmente, se deshace el giro para obtener la distribución las fuerzas de reacción articular de la cadera en el sistema de coordenadas de Abaqus®. Teniéndose en este sistema de referencia, ya se pueden incluir las cargas en el archivo “.inp” de entrada para llevar a cabo las simulaciones. 2.3 Remodelación ósea 2.3.1 Fundamento teórico En la literatura se pueden encontrar distintos modelos de remodelación ósea empleados para estimar la distribución de densidad de un hueso a partir de su distribución de tensiones. A lo largo de los últimos años, varios estudios se han llevado a cabo de esa forma [8]. Uno de dichos modelos, desarrollado por Doblaré et al. [7], es el que se utilizó en la tesis de J. Ojeda [9] y se empleará en este trabajo para estimar la distribución de densidad en el fémur. Se trata de un modelo fenomenológico que considera el comportamiento anisótropo del hueso y se basa en el modelo de remodelación de Stanford [6]. 2.3.1.1 Modelo de remodelación de Stanford En este modelo se asumen ciertas ideas básicas con relevancia para el modelo aquí propuesto. En primer lugar, se establece una condición homeostática para el comportamiento del hueso a nivel 23 local. Esto quiere decir que el tejido óseo puede regular su densidad y, en consecuencia, sus propiedades mecánicas, basándose en una cierta magnitud fisiológica denominada estímulo mecánico. El estímulo mecánico local, ψ𝑡, es una medida de cuán intensamente está siendo cargado el hueso y considera los diferentes casos de carga en la actividad diaria normal. ψ𝑡=(∑𝑛𝑖𝜎𝑡𝑖 𝑚 𝑁 𝑖=1 )1/𝑚 (2.10) Donde 𝑁 es el número de diferentes casos de carga, 𝑛𝑖 es el número medio de ciclos de carga por salto de tiempo para cada caso de carga 𝑖, 𝑚 es un parámetro empírico que varía entre 3 y 8 y 𝜎𝑡𝑖 es la tensión efectiva a nivel tisular. Beaupré et al. [6] relacionaron la tensión a nivel de tejido, 𝜎𝑡, con la tensión a nivel continuo, 𝜎, siguiendo un procedimiento de homogeneización estándar y apoyado por datos experimentales [17]. Se comprobó que la tensión de rotura del hueso trabecular a nivel continuo, 𝜎𝑢𝑙𝑡, se relaciona a través de la densidad aparente con la tensión de rotura a nivel tisular, 𝜎𝑢𝑙𝑡𝑡: 𝜎𝑢𝑙𝑡(𝜌)=(𝜌𝜌)2𝜎𝑢𝑙𝑡𝑡 (2.11) Extrapolando dicha relación al estrés en el hueso a nivel continuo y tisular, se define la ecuación (2.12). 𝜎𝑢𝑙𝑡(𝜌)=(𝜌𝜌)2𝜎𝑢𝑙𝑡𝑡 (2.12) Entonces, el estímulo mecánico a nivel tisular se puede expresar como función de la tensión a nivel continuo: ψ𝑡=(∑𝑛𝑖1/𝑚 𝑁 𝑖=1 (𝜌𝜌)2𝜎𝑖) (2.13) Con 𝜎𝑖 definido por 𝜎𝑖=√2𝐸𝑈 (2.14) Siendo 𝐸 el módulo de Young y 𝑈 la densidad de energía de deformación, ambos a nivel continuo. También es posible relacionar el estímulo mecánico a nivel de tejido con el estímulo a nivel continuo mediante la ecuación (2.15). ψ𝑡=(𝜌𝜌)2ψ (2.15) La variación de tensiones debido a cargas cíclicas puede ocurrir en cuestión de segundos, mientras que la respuesta remodelatoria puede tomar días. Esta idea sugiere que el hueso está, de alguna manera, “recordando” la información del estado tensional y respondiendo pasado cierto tiempo, en base a la información acumulada durante los días previos. Carter propuso que se Metodología 30 𝑟󰇗= { −𝑐1𝑔𝑟 𝜌2−𝛽 2 𝑠𝑖 𝑔𝑟≥0 𝑦 𝑔𝑓<0 0 𝑠𝑖 𝑔𝑟<0 𝑦 𝑔𝑓<0 𝑐2𝑔𝑓 𝜌2−𝛽/2 𝑠𝑖 𝑔𝑟<0 𝑦 𝑔𝑓≥0 (2.39) Esta respuesta remodelatoria es equivalente a la respuesta del modelo de Stanford y presenta una limitación importante. De acuerdo a esta ecuación, un estímulo mecánico por encima del estímulo de referencia produce formación de hueso y a una elevada velocidad si dicho estímulo es muy alto. Como consecuencia, la rigidez aumentará muy rápidamente. Sin embargo, el modelo no considera que el hueso podría estar dañado como resultado de cargas muy elevadas y este daño podría llevar a reabsorción ósea con el fin de repararlo. El modelo anisótropo basado en la mecánica del daño no considera las diferencias en el comportamiento del hueso cuando el tejido se encuentra a tracción o a compresión. Esta limitación se debe a la ecuación que define el estímulo (2.31), cuadrática en deformaciones y, por tanto, independiente del signo de éstas. Parámetro Valor ψ𝑡∗ 50 𝑀𝑃𝑎/𝑑í𝑎 𝑐1 0.02 𝜇𝑚/𝑑í𝑎 𝑐2 0.02 𝜇𝑚/𝑑í𝑎 𝑤 12.5 𝑀𝑃𝑎/𝑑í𝑎 𝜔 0.1 𝑛 6000 𝑐𝑖𝑐𝑙𝑜𝑠/𝑑í𝑎 𝑘 1 𝑚 4 Tabla 2. Valores de las constantes del modelo de remodelación empleado en las simulaciones. 31 2.3.2 Proceso en Abaqus® A lo largo de todo el texto de la memoria se ha hecho referencia en múltiples ocasiones al software de elementos finitos con el que se ha trabajado para llevar a cabo el trabajo: Abaqus®. Esta sección está dedicada, en primer lugar, a introducir de forma breve este programa y, a continuación, la forma de llevar a cabo los procesos abordados con él y una descripción de los archivos y herramientas empleados para llegar a los resultados de la remodelación ósea. Abaqus® es un programa CAE (Computed Aided Engineering) de cálculo por elementos finitos de propósito general, parte de la plataforma SIMULIA de Dassault Systemes. Abaqus®/CAE es un entorno de trabajo cuya interfaz integra modelado, análisis, administración y visualización de las simulaciones [24]. Existe la posibilidad de trabajar directamente desde la interfaz Abaqus/Viewer, en la que se puede modelar desde la geometría hasta el estado de carga y llevar a cabo las simulaciones. Sin embargo, en este proyecto se ha trabajado introduciendo, mediante la ventana de comandos, Abaqus Command, los archivos de extensión “.inp” (input). Estos contienen la información de los nodos y elementos que conforman la malla de elementos finitos, del material del modelo y de las cargas y restricciones aplicadas sobre el mismo. A continuación se describe el contenido de los archivos “.inp” que se emplean como input en cada una de las simulaciones llevadas a cabo. En primer lugar hay un “.inp” al que se le podría denominar principal. Es el que contiene la información de los nodos y elementos, así como del material. También es el que pone en marcha la simulación, llamando a otro “.inp” de fase de carga, el cual contiene los nodos de aplicación y las componentes de todas las fuerzas que actúan sobre el modelo. De la implementación por la ventana de comandos del input principal, tras la simulación, se obtiene, entre otros, un archivo de extensión “.odb” (output database) que permitirá el post procesado de los resultados de la simulación en Abaqus/Viewer. Figura 2.18. Inputs y outputs de la primera iteración por ventana de comandos. La forma de proceder en las simulaciones fue partir de un estado inicial isotrópico y homogéneo, con densidad inicial 𝜌0=0.5 𝑔/𝑐𝑚3, como se ha realizado en trabajos de otros autores (). Se procede a aplicar 50 veces la misma fase (se traduce en simular 50 días de solicitación en el mismo estado de carga) sobre el estado inicial, desde el “.inp” principal. Posteriormente se emplearía un último “.inp” de reinicio (restart), cuya función es retomar la simulación desde el Metodología 32 punto en el que se había dejado, aplicándose otras 50 veces la fase de carga. Cabe destacar que para poder retomar la simulación mediante el archivo de reinicio mencionado, es preciso que en el directorio de trabajo se encuentren todos los archivos generados tras la última simulación realizada. Sin embargo, para el análisis de los resultados, será suficiente con conservar los “.odb” de cada iteración del reinicio. Por ello, debido al espacio que ocupa en memoria todo el volumen de archivos que se genera en el proceso, a medida que se van realizando iteraciones, se van conservando todos los output únicamente de la última simulación realizada, mientras que de las anteriores solo se guarda el “.odb”. Figura 2.19. Inputs y outputs de cada reinicio por ventana de comandos. 2.3.2.1 Archivos “.inp” En este apartado se detallará el contenido de cada uno de los inputs que se implementarán desde la ventana de comandos. En pos de la explicación, se empleará para los distintos archivos los nombres empleados durante esta sección (principal, fase de carga, reinicio). Sin embargo, como es evidente, en la práctica los nombres de estos archivos contendrán el nombre del músculo y la dirección en la que se está modificando el punto de inserción en cada grupo de archivos del análisis de sensibilidad. • Principal Es el input que implementa la primera iteración del proceso de remodelación ósea. En primer lugar contiene la definición de todos los nodos, definiéndose, para los 64757 nodos, en primer lugar la etiqueta con el número que corresponde a cada nodo y a continuación las coordenadas de éste en el sistema de referencia global de Abaqus® (Figura 2.20). 33 Figura 2.20. Definición de los nodos. A continuación, se definen los elementos que conforman el modelo a partir de los nodos generados al comienzo del código. Esto se lleva a cabo seleccionando en primer lugar el tipo de elemento (en este caso, C3D4, elementos tetraédricos con 4 nodos por elemento). Tras especificar el tipo, se indica, para los 339168 elementos, la etiqueta con el número de cada elemento y los cuatro nodos que lo conforman, como se muestra en la Figura 2.21. Figura 2.21. Definición de los elementos. Metodología 34 Una vez definidos los nodos y elementos, ya se tendría la geometría creada en Abaqus®. A continuación se definen una serie de constantes del material que determinarán su comportamiento en el proceso de remodelación, así como el número de ciclos que se aplican por cada step en cada fase. Por último, se incluye el estado de carga al que se someterá al hueso. Dicho estado se introducirá, como se ha explicado con anterioridad, mediante otro “.inp” que incluirá los valores de carga que se aplicarán en cada nodo solicitado, tal como se verá a continuación. La llamada a la fase de carga, aparece resaltada en rojo en la Figura 2.22. En este trabajo se ha empleado para cada análisis una única fase de carga, que se ha repetido 50 veces. Por tanto, el output de la implementación de la primera iteración representará el resultado de la aplicación de un mismo estado de carga durante un tiempo que equivale aproximadamente a 50 días (al considerarse un día como la repetición de los 6000 ciclos de cada escalón o step). Figura 2.22. Constantes del material y fase de carga. • Fase de carga En este archivo se incluyen las cargas tanto de las fuerzas musculares como las de reacción articular. Éste será al que se acceda cuando el “.inp” principal solicite la inclusión de un input. La Figura 2.23 muestra parte del contenido de uno de los archivos de fase de carga empleados. El primer número que se incluye es el nodo sobre el que se aplica la carga a definir. A continuación, aparece un 1, 2 ó 3, que hace referencia a la componente de la fuerza a la que corresponderá el valor introducido. De esta forma, para las fuerzas musculares se incluyen tres 35 líneas por cada nodo, introduciéndose en cada una el valor correspondiente a cada una de las tres componentes obtenidas previamente en la sección 2.2.2. Con las fuerzas de reacción articular se procede exactamente de la misma forma, introduciéndose los valores sacados de la distribución como paraboloide de revolución, de la sección 0. Figura 2.23. Definición de la fase de carga. • Reinicio Por último, el “.inp” de reinicio contiene la llamada a la lectura del último archivo generado en el proceso para continuar con la simulación detenida (“RESTART, READ”) y las fases de carga introducidas de la misma manera que en el principal. De nuevo, en el restart se introduce la misma fase 50 veces. Figura 2.24. Input de reinicio. Metodología 36 2.3.2.2 Post-procesado y convergencia de la remodelación Para saber si la remodelación ha convergido es necesario llevar a cabo un post procesado de los resultados obtenidos en cada iteración del proceso. Abriendo los archivos “.odb” que se obtienen en la interfaz Abaqus/Viewer, se puede obtener un listado de los resultados de interés. En este caso, teniendo en cuenta el criterio adoptado para saber si el proceso ha convergido, interesa sacar los resultados de las densidades aparentes, así como los volúmenes de cada elemento. Las densidades y volúmenes de cada elemento se obtendrán, como se ha dicho, desde el “.odb” en la interfaz de visualización. En ella, se habrá de seleccionar, de la pestaña “Report”, la opción “Field output” (rodeada en rojo en la Figura 2.25). Figura 2.25. Interfaz Abaqus/Viewer. Esta acción hará que aparezca una ventana emergente como las que se muestran en la Figura 2.26. En ellas se ilustra cómo se ha de seleccionar la variable cuyos datos queremos almacenar en el documento que se generará (en este caso, volumen o densidad), así como las opciones que se ofrecen para nombrar el archivo y especificar los datos a incluir en el archivo. Una vez se han seleccionado las opciones pertinentes, al hacer click sobre “OK” se generará un archivo de extensión “.rpt” (report) cuyo contenido incluirá en una columna los números con los que están etiquetados los elementos y en otra los respectivos valores de la variable deseada. 37 Figura 2.26. Opciones para generar el “.rpt”. A partir de dichos valores se podrá conocer la densidad media, 𝜌, en cada iteración, mediante la fórmula de la ecuación (2.40). 𝜌=∫𝜌𝑑𝑉 𝑉𝑜𝑙 ∫𝑑𝑉 𝑉𝑜𝑙 (2.40) El reinicio explicado anteriormente se realiza hasta que la densidad media del modelo se haya estabilizado. Se ha considerado que el proceso de remodelación ha convergido cuando entre una simulación y la anterior la variación de la densidad media es menor del 0.5% (se logra en todos los casos para 500 días de solicitación, o 9 reinicios). La Figura 2.27 muestra la evolución temporal de la densidad media del modelo de elementos finitos del fémur para el caso nominal, es decir, para el caso en el que no se ha llevado a cabo ningún cambio en los puntos de inserción. Se puede apreciar cómo el valor de la densidad media converge en 0.8482 𝑔/𝑐𝑚3. Teniendo ya la distribución de densidades final del fémur para el estado de carga aplicado, en Abaqus/Viewer se representan también los diagramas que se deseen ilustrar, para los que se presenta una amplia variedad de opciones. Metodología 38 Figura 2.27. Evolución temporal de la densidad media para el caso nominal. 39 3 ANÁLISIS DE SENSIBILIDAD A lo largo del capítulo anterior se ha desarrollado la metodología para determinar el estado de carga sobre el fémur, en un instante particular del ciclo de marcha y para unas ubicaciones concretas de los puntos de inserción de cada uno de los músculos en el fémur. Asimismo, se han detallado los pasos para, a partir de dichas solicitaciones y de la geometría del hueso, alcanzar la distribución de densidades en el fémur, solicitado según las cargas generadas en ese instante y para esas inserciones. El contenido del presente capítulo y el objetivo principal del estudio residen en realizar un análisis de sensibilidad para comprobar cómo y en qué medida afecta a la distribución final de densidades el hecho de modificar levemente la ubicación de alguna de las inserciones de los músculos. Por tanto, se habrán de repetir los pasos descritos en el capítulo 2 tantas veces como puntos de inserción se modifiquen. Es preciso remarcar que al cambiar el punto de inserción de alguno de los músculos, evidentemente, no solo cambia el nodo sobre el que se aplica la fuerza de ese músculo, sino que se están alterando las condiciones iniciales del problema de tal forma que todas las fuerzas, tanto musculares como de reacción articular, se ven modificadas y deben ser obtenidas de nuevo. 3.1 Modificación de las fuerzas musculares En [14] se realizó un análisis de sensibilidad en el que se estudiaba cómo varían los módulos de las fuerzas musculares a lo largo del ciclo de marcha al desplazar tanto los puntos de origen como los de inserción de los músculos de la cadera. Para cada uno de los músculos se realizaron 6 ensayos en cada una de las 3 direcciones, tanto para los puntos de origen como para los de inserción. De esta forma, se tienen 36 archivos de Microsoft Excel® conteniendo las evoluciones temporales de las fuerzas musculares para cada uno de los músculos estudiados. Es preciso conocer la forma en la que se llevó a cabo dicho análisis de sensibilidad, pues los desplazamientos de los puntos de inserción son los mismos que se han tenido que llevar a cabo para sacar las evoluciones temporales a lo largo del ciclo de marcha de cada una de las tres componentes de las fuerzas musculares. Para cada dirección de cada músculo, en origen e inserción, se realizaron, como ya se ha dicho, 6 ensayos. En cada uno de los grupos de 6 ensayos, se seleccionó una variación de la posición ∆𝑑. De esta forma se realizó desde el punto del caso Análisis de Sensibilidad 46 3.4.1 Psoas La distancia seleccionada para los ensayos del psoas es: ∆𝑑=1.333 𝑐𝑚 De los 140 frames del ciclo de marcha, la fuerza máxima producida por el psoas en el caso nominal se da en el frame número 62. Fuerza máxima psoas (𝑁) Componente X (𝑁) Componente Y (𝑁) Componente Z (𝑁) 782.70 27.52 -472.30 623.54 Tabla 3. Fuerza máxima del psoas. Caso Densidad media final (𝑔/𝑐𝑚3) Masa final (𝑘𝑔) Nominal 0.8315 0.4916 Dirección Y (+) 0.7832 0.4630 Dirección Y (-) 0.7351 0.4346 Dirección Z (+) 0.8543 0.5051 Dirección Z (-) 0.7624 0.4508 Tabla 4. Densidad y masa finales del fémur (ensayos del psoas). Figura 3.8. Comparación de los canales medulares obtenidos en los casos mostrados del psoas. 47 Figura 3.9. Resultados del caso nominal del psoas. Distribución de tensiones en el instante inicial, cuando la densidad es homogénea (a); y distribución final de densidades (b). Figura 3.10. Resultados del ensayo en la dirección Y en sentido negativo del psoas. Distribución de tensiones en el instante inicial, cuando la densidad es homogénea (a); y distribución final de densidades (b). Análisis de Sensibilidad 48 3.4.2 Glúteo medio La distancia seleccionada para los ensayos del glúteo medio, para los tres actuadores, es: ∆𝑑=0.667 𝑐𝑚 De los 140 frames del ciclo de marcha, la fuerza máxima producida por el actuador 1 del glúteo medio en el caso nominal se da en el frame número 51. Actuador Fuerza máxima gmed (𝑁) Comp. X (𝑁) Comp. Y (𝑁) Comp. Z (𝑁) 1 591.52 128.84 -150.75 557.29 2 210.53 113.19 12.18 177.10 3 139.53 103.39 37.44 85.90 Tabla 5. Fuerzas del glúteo medio cuando la del actuador 1 es máxima. Caso Densidad media final (𝑔/𝑐𝑚3) Masa final (𝑘𝑔) Nominal 0.8711 0.5150 Dirección Y (+) 0.8874 0.5247 Dirección Y (-) 1.0580 0.6256 Dirección Z (+) 0.8609 0.5090 Dirección Z (-) 0.8717 0.5154 Tabla 6. Densidad y masa finales del fémur (ensayos del glúteo medio). Figura 3.11. Comparación de los canales medulares obtenidos en los casos mostrados del glúteo medio. 49 Figura 3.12. . Resultados del caso nominal del glúteo medio. Distribución de tensiones en el instante inicial, cuando la densidad es homogénea (a); y distribución final de densidades (b). Figura 3.13. Resultados del ensayo en la dirección Y en sentido negativo del glúteo medio. Distribución de tensiones en el instante inicial, cuando la densidad es homogénea (a); y distribución final de densidades (b). Análisis de Sensibilidad 50 3.4.3 Aductor largo La distancia seleccionada para los ensayos del aductor largo es: ∆𝑑=0.667 𝑐𝑚 De los 140 frames del ciclo de marcha, la fuerza máxima producida por el aductor largo en el caso nominal se da en el frame número 69. Fuerza máxima addlong (𝑁) Componente X (𝑁) Componente Y (𝑁) Componente Z (𝑁) 165.61 57.78 -29.84 152.30 Tabla 7. Fuerza máxima del aductor largo. Caso Densidad media final (𝑔/𝑐𝑚3) Masa final (𝑘𝑔) Nominal 0.6944 0.4105 Dirección X (+) 0.6950 0.4109 Dirección X (-) 0.6919 0.4100 Dirección Z (+) 0.6945 0.4106 Dirección Z (-) 0.6936 0.4101 Tabla 8. Densidad y masa finales del fémur (ensayos del aductor largo). Figura 3.14. Comparación de los canales medulares obtenidos en los casos mostrados del aductor largo. 51 Figura 3.15. . Resultados del caso nominal del aductor largo. Distribución de tensiones en el instante inicial, cuando la densidad es homogénea (a); y distribución final de densidades (b). Figura 3.16. Resultados del ensayo en la dirección X en sentido negativo del aductor largo. Distribución de tensiones en el instante inicial, cuando la densidad es homogénea (a); y distribución final de densidades (b). Análisis de Sensibilidad 52 3.4.4 Resonancia magnética Este último ensayo surge de comparar el caso nominal con el de un modelo en el que se han modificado las inserciones del tensor de la fascia lata y los glúteos (mayor, medio y menor). Estos cambios se han realizado de acuerdo a las posiciones reales de las inserciones del sujeto, determinadas mediante el estudio de éste con resonancias magnéticas. El instante escogido para la carga estática del modelo del fémur ha sido aquél en el que las fuerzas de reacción articular sobre la cabeza del fémur son máximas (frame número 61). Caso Densidad media final (𝑔/𝑐𝑚3) Masa final (𝑘𝑔) Nominal 0.8482 0.5015 Resonancia magnética 0.8268 0.4889 Tabla 9. Densidad y masa finales del fémur (ensayos nominal y de resonancia). Figura 3.17. Resultados del ensayo con las inserciones de la resonancia magnética. Distribución de tensiones en el instante inicial, cuando la densidad es homogénea (a); y distribución final de densidades (b). 53 3.4.5 Discusión de los resultados Con los resultados mostrados en las páginas previas, se tiene suficiente información para llevar a cabo de forma razonada una discusión de los resultados obtenidos. En todos los casos que han sido objeto de estudio y, en particular, en los casos que se han ilustrado en la presente memoria, se aprecia cualitativamente la relación entre el estado tensional al que se somete al hueso y la distribución de densidades del mismo. En aquellos ensayos en los que aparecen zonas cuya solicitación es menor que en el caso nominal, se observa una disminución también en la densidad final en dichas áreas. Por tanto, como cabía esperar de acuerdo con la teoría, las distribuciones terminan siendo cualitativamente semejantes, puesto que en las zonas con mayor solicitación la densidad aumenta y en las zonas en desuso, se produce reabsorción. Como consecuencia de realizar las simulaciones empleando cargas estáticas y seleccionando para cada músculo estudiado el instante en el que su carga es máxima, los valores de densidad en los que converge el caso nominal para cada uno de los músculos es distinto, al tratarse de distintas fases del ciclo de marcha. Esto es, para las mismas ubicaciones de todos los puntos de inserción, solo que en instantes diferentes, se han obtenido distribuciones de densidad diferentes. Entre el instante de máxima fuerza de reacción articular y los de fuerza muscular máxima del psoas y del glúteo medio, dicha diferencia no es significativa. Sin embargo, en el caso del aductor largo, la variación de densidad aparente final respecto a la que se obtiene en el instante de máxima fuerza de reacción articular es del 18%. Un cambio tan notable debe ser considerado y a priori podría ser indicativo de la importancia de tener en cuenta todo el ciclo de marcha aplicando las cargas de una forma pseudo-dinámica (ciclos de distintas fases de cargas), y no solo cargas estáticas en determinados instantes, para conocer la distribución de densidades final del hueso. En el estudio del psoas, en el ensayo para el que se ha obtenido una mayor variación en el valor final de la densidad aparente del fémur respecto al caso nominal, ésta ha sido de un 11.6% para un desplazamiento del punto de inserción de 1.333 cm. En cuanto al glúteo medio, la mayor variación que se ha obtenido ha sido del 21.5% para un desplazamiento de 0.667 cm. Por otra parte, para el aductor largo la variación máxima obtenida es de apenas un 0.3% para un desplazamiento de también 0.667 cm. Observando las distribuciones de densidades en las secciones transversales mostradas se pueden sacar algunas conclusiones en relación a la formación más o menos clara de un canal medular en el modelo. Este canal medular aparece en el modelo debido a cargas torsionales y se manifiesta como una zona interior de hueso trabecular o esponjoso en la diáfisis, en contraste con el hueso cortical, más denso, que se presenta en la superficie. En el caso del psoas, para el caso nominal aparece de una forma más o menos difusa un canal medular en el centro, alrededor del cual se va haciendo el hueso cada vez más denso hasta la superficie, de donde se deduce que dominan las fuerzas de torsión. Sin embargo en el ensayo ilustrado, al haber modificado el punto de inserción del psoas, cuya fuerza tiene gran importancia en el instante de estudio, la distribución de densidades se hace más propia de cargas de flexión. Para el glúteo medio parece ocurrir precisamente el efecto contrario, además del aumento de densidad que se produce en toda la sección. Las secciones de ambos ensayos del aductor corto muestran la misma distribución, presentando un canal medular fruto de cargas a torsión sobre el fémur, con escasa diferencia entre ambas, como se deducía de la casi nula diferencia entre las densidades aparentes finales. Análisis de Sensibilidad 54 En cuanto a la escasa diferencia entre los casos nominal y con los datos de la resonancia magnética del sujeto en el instante de máxima fuerza de reacción articular, esto se debe probablemente a que, si bien se están modificando cuatro músculos distintos, las mayores diferencias en los nodos de inserción se dan para actuadores cuya fuerza no tiene una influencia dominante. Asimismo, el nodo del actuador 1 del glúteo medio (el que tiene una mayor importancia en ese instante en cuanto a carga se refiere) apenas se ve modificado, al igual que el tensor de la fascia lata. Por ello, pese a que se modifiquen varios músculos, no parece desencaminado el hecho de que la variación en las densidades sea más bien leve. 55 4 CONCLUSIONES Y TRABAJO FUTURO Llegado este punto, el presente apartado se destinará a recopilar y analizar el trabajo realizado, expuesto a lo largo de toda la memoria, y a concluir el alcance que se ha podido conseguir con el estudio. Tal y como se estableció al comienzo de la memoria, el objetivo del trabajo ha sido comprobar cómo y en qué medida influyen los desplazamientos de los puntos de inserción en el modelo de elementos finitos sobre la distribución final de densidades tras el proceso de remodelación ósea. 4.1 Conclusiones En este Trabajo Fin de Grado se ha propuesto una metodología para modelar las cargas a aplicar –tanto fuerzas musculares como de reacción articular– sobre un modelo de fémur en Abaqus®, a partir de la toma de datos de la marcha y su procesado mediante un modelo de sistemas multicuerpo de OpenSim. Empleando dicha metodología, se ha podido llevar a cabo un análisis de sensibilidad en el que se han variado los nodos del modelo sobre los que se aplican dichas fuerzas para observar las diferencias en las distribuciones de densidad de los diferentes ensayos. En vista de los resultados obtenidos, se puede concluir que desplazamientos del orden de un centímetro en el punto de inserción de una fuerza dominante en el estado de carga, dan lugar a variaciones importantes en la distribución final de densidad y, por consiguiente, en el valor de la densidad media del hueso, si la inserción se encuentra en una zona del hueso con una geometría crítica, como lo son los trocánteres. Sin embargo, en la diáfisis, al ser la geometría más suave y la fuerza del músculo estudiado de menor magnitud, las variaciones del orden del centímetro no conllevan un cambio apreciable en la remodelación del hueso. Todo lo anterior pone de manifiesto un hecho que se viene presentando ya desde los primeros apartados: la introducción de errores en el modelado del problema al emplear los modelos de sistemas multicuerpo; los cuales, si bien pueden ser objeto de un escalado de sus dimensiones, no dejan de presentar una geometría fruto del promediado del estudio de cadáveres de sujetos distintos al de estudio. Los resultados evidencian que el error que ello introduce supone potencialmente una imprecisión importante en la distribución de densidades obtenida. De ahí la importancia de modificar los puntos de inserción de acuerdo al problema en concreto y no solo realizar el escalado mencionado más arriba.