1 Departamento de Mecánica, Universidad Simón Bolívar, Venezuela. e-mail: cmuller@usb.ve,

² Centro de Bioingeniería, Universidad Central de Venezuela, e-mail: mcerrola@reacciun.ve

En este trabajo se diseña y desarrolla un modelo tridimensional de los huesos, en condiciones fisiológicas, para el mejor entendimiento del comportamiento óseo. Para realizar esta simulación es necesario desarrollar una metodología para la extracción de la información geométrica y propiedades de los materiales partiendo de las Tomografías Axiales Computarizadas y crear un arreglo de datos con las propiedades mecánicas del hueso que sirva de entrada para el método de elementos finitos con la versión p, que es la herramienta que se utilizará para la determinación de tensiones y deformaciones. Se creó un modelo de la tibia humana heterogéneo aplicando cargas fisiológicas para determinar los esfuerzos y deformaciones, para el cual se compararon los resultados con modelos presentes en la literatura. Los resultados demuestran que la metodología presentada es adecuada para la simulación ósea y supera algunas dificultades de los modelos presentes en la literatura.

Palabras clave: Huesos, simulación, versión p de elementos finitos, superficies suaves, heterogeneidad del hueso.

In this work a three-dimensional model of bone, under physiological conditions, is developed with the purpose of a better understanding of the bone tissue. To accomplish the simulation it is necessary to develop a methodology for the extraction of geometry and material properties from Computerize Axial Tomographies. An ASCII format properties-array is created as input data for the Finite Element Program in p-version. The stresses and displacements are calculated for a human tibia model as a heterogeneous body. The results and comparisons with the literature demonstrate that the developed methodology is adequate for bone simulation and overcomes some of the difficulties reported by other investigators.

Key words: bone, simulation, P-version of finite element method, smooth surfaces, bone heterogeneity.

Los recientes avances en la tecnología aumentan cada día la esperanza de una vida mas larga y sana, por medio de la aplicación de una medicina menos invasiva y más eficiente. Entre los avances que se pueden mencionar está la creación de equipos como tomógrafos computacionales y resonadores magnéticos, que permiten observar al cuerpo humano internamente sin tener que intervenir quirúrgicamente para realizar un diagnóstico. Por otro lado, existen operaciones con rayos láser o con aparatos sumamente pequeños como son los utilizados en las artroscopias que permiten operar y curar a un paciente con una mínima interacción con la integridad del cuerpo humano.

En este trabajo se presenta una metodología para la reconstrucción, la simulación y el análisis mecánico de los huesos utilizando Tomografías Axiales Computarizadas tanto para extraer la geometría del hueso como los materiales que lo componen. El modelo supera algunas dificultades presentes en otros modelos publicados recientemente por otros investigadores. Para lograr este objetivo, se utilizan el método de elementos de contorno, el método de Elementos Finitos en su versión clásica y en su versión p-adaptable, ésta última es una metodología bastante reciente y que no ha sido aplicada al área de biomecánica.

La conducta biomecánica del cuerpo humano es determinada por huesos, cartílagos, ligamentos, tendones, músculos y otros tejidos conjuntivos. Estos elementos se clasifican como estructuras activas o pasivas dependiendo de si ellos producen o no fuerzas. Los huesos son considerados estructuras pasivas y constituyen la mayoría del esqueleto. La matriz ósea es un material compuesto formado por una componente orgánica (aprox.65%), otra inorgánica (aprox.20%) y agua (aprox. 10%), estos valores varían dependiendo del tipo de hueso. La matriz orgánica está constituida aproximadamente por 95% de fibras de colágeno reforzadas por depósitos de calcio y sales de fosfato en forma de hidroxiapatita (Cowin, 1989). Los depósitos de calcio y fosfato son los que le proporcionan dureza, rigidez y resistencia al hueso, mientras que las fibras de colágeno son las que proporcionan flexibilidad. El hueso joven u osteoide contiene inicialmente mucha agua que será desplazada por mineral en períodos de días y/o meses. La componente orgánica le da forma al hueso y contribuye a la capacidad de resistir a la tracción, mientras que la componente inorgánica o mineral contribuye a la resistencia a la compresión. Los huesos desmineralizados (ligamentos o tendones) son flexibles y resistentes a la tracción. Con el incremento de la mineralización en el proceso de maduración ósea, el hueso se hace más rígido, mas quebradizo y menos elástico, de allí la diferencia en el comportamiento a la fractura de un hueso de un niño y un hueso de un adulto.

Los huesos en el sistema esquelético son estructuras complejas compuestas de dos tipos de hueso (ver figura 1) que tienen un comportamiento y una función bien diferenciada, ellos son:

• El hueso esponjoso o trabecular

• El hueso compacto o cortical.

Se han estudiado las propiedades mecánicas de estos dos tipos de huesos extensivamente y varios informes están disponibles en la literatura (Buckwalter et. al. 1995, Pettermann et. al 1997). La mayoría de los estudios están de acuerdo que el hueso cortical y el hueso esponjoso contienen la misma composición de la matriz y estructura, pero el hueso cortical tiene una porosidad mucho menor (1:5) que el hueso trabecular. La porosidad se define como el volumen de vacío por unidad de volumen de hueso, y representa la parte proporcional del hueso ocupado por médula ósea o material no-mineralizado que se encarga de la irrigación del hueso. El hueso compacto tiene una porosidad de aproximadamente 5 a 30% mientras que la porosidad del hueso trabecular es aproximadamente 30 a 90%. La geometría y orientación de las trabéculas contribuyen en el comportamiento anisotrópico estructural que es evidente tanto en el hueso trabecular, como en el hueso compacto. El hueso trabecular se encuentra principalmente en los huesos planos y cortos o en la epífisis de los huesos largos (ver figura 1). Los poros del hueso trabecular interconectados de forma irregular proporcionan una densidad aparente y unas propiedades mecánicas del hueso inconstantes. El módulo de elasticidad y la resistencia a la fractura del hueso cortical pueden ser diez veces mayor que aquél de un volumen similar de hueso trabecular. El hueso esponjoso tiene aproximadamente veinte veces mas superficie por unidad de volumen que el cortical, por lo cual se cree que debido a esta diferencia, el hueso esponjoso presenta habitualmente una mayor actividad metabólica y una mayor velocidad de remodelación, respondiendo mas rápidamente a las cargas mecánicas que el hueso cortical. La característica porosa del hueso esponjoso tiene una función importante en la absorción de fuerzas en las articulaciones, el reemplazo de la epífisis y del cartílago que la recubre por una prótesis, unida al hueso por medio de un cemento de polimetilmetacrilato, elimina la propiedad de absorción de impactos lo cual aumenta drásticamente las fuerzas transmitidas por la articulación.

El hueso puede ser clasificado según su geometría en tres grupos principales: huesos cortos, llanos y largos o tubulares. Los huesos cortos, como tarsianos, carpapianos y cuerpos vertebrales, miden aproximadamente lo mismo en todas las direcciones y son de forma trapezoidal, cuboidal, cuneiforme o irregular. Los huesos "planos" tienen una dimensión que es mucho menor que las otras dos, entre ellos están la escápula y las crestas del ilión. Los huesos "largos" tienen una dimensión que es mucho mayor que las otras dos, como por ejemplo el fémur, la tibia, el húmero, los metacarpianos, los metatarsianos y las falanges. Los huesos largos se describen basándose en su composición como muestra la figura 1. Estos tienen una zona tubular formada mayormente por hueso cortical denominada eje del hueso o diáfisis. En los extremos del hueso la diáfisis se transforma en una zona más amplia compuesta por hueso trabecular denominada epífisis que se articula con otros huesos y está protegida por una capa de cartílago de hialine llamado cartílago articular. La diáfisis es una estructura que rodea la cavidad medular. En el canal medular se encuentra alojada la medula ósea que tiene como función la irrigación a las células óseas, por lo que cualquier problema que presente la médula afecta las funciones del hueso. Entre la epífisis y la diáfisis hay una zona creciente durante el desarrollo que se llama el metáfisis. El grueso hueso cortical en la diáfisis proporciona resistencia a la torsión y a la flexión. En la epífisis el apoyo de la delgada capa cortical que cubre al hueso trabecular distribuye cargas mecánicas de las superficies articulares al árbol del hueso y permite mayor deformación, ayudando a absorber las carga de impacto aplicadas en las articulaciones sinoviales (figura 2). Los huesos se encuentran recubiertos por una membrana que ayuda a su crecimiento, el periostio externa y endostio internamente, este último separa el hueso de la médula.

Los huesos realizan varias funciones mecánicas (Nigg y Herzog, 1994), entre las más importantes están:

• Proporcionar apoyo al cuerpo para soportar fuerzas externas (p. ej.: gravedad).

• Actuar como un sistema de la palanca para transferir fuerzas (p. ej. las fuerzas musculares).

• Proporcionar protección para los órganos interiores vitales (p. ej. cerebro, pulmones etc.).

• Funciones fisiológicas como tomar parte en el proceso de circulación sanguínea (hematopoiesis).

• Almacenar calcio, fósforo e iones que utilizan los músculos para su funcionamiento (homeostasis mineral).

El estudio de las propiedades mecánicas del hueso permite predecir las fuerzas que el hueso es capaz de resistir, las posibles consecuencias de las enfermedades, entender el efecto de envejecimiento y otras características. Currey (1970) afirma que el hueso es casi dos veces más resistente a compresión de lo que lo es a tensión. La resistencia máxima a tensión se ha considerado como la propiedad mecánica más importante de esta estructura. Por consiguiente, la mayoría de las pruebas hechas para estudiar las propiedades mecánicas de los huesos han sido a tensión con probetas orientadas a lo largo del eje axial de los huesos.

Una gran cantidad de investigadores han hecho trabajos experimentales para determinar los valores de estas constantes, entre ellos se puede nombrar a Ashman et. al. (1984) quienes utilizaron una técnica de ondas acústicas para encontrar el promedio del módulo de Young en dirección transversal (13.4 GPa) y en dirección longitudinal (20.0 GPa); Zysset et. al. (1999) que llevaron a cabo un estudio mecánico para medir el módulo de Young del hueso trabecular y del hueso cortical, ellos encontraron que el módulo de Young en la dirección longitudinal es aproximadamente 40% mayor que el de la dirección transversal. También notaron que el módulo de Young del hueso trabecular era algo mayor que el del hueso cortical transversal y sustancialmente menor que el del hueso cortical en dirección longitudinal. Ellos estimaron el coeficiente de Poisson entre 0.2 y 0.4, obteniendo que la variación en este rango no producía mayor variación en los resultados. Algunos investigadores han publicado las propiedades del hueso en todas las direcciones, entre los mas destacados están Reilly y Burstein (1975), Yoon y Katz (1976) quienes consideraron el hueso transversalmente isotrópico; Knets et. al. (1981) y Ashman et. al. (1984) quienes consideraron el material como ortotrópico. En la tabla 1 se muestran los resultados de estas investigaciones.

Considerando que el hueso es una geometría complicada compuesta por material no homogéneo, anisotrópico, no lineal y viscoelástico, una de las mayores dificultades es medir de forma precisa la geometría y las propiedades mecánicas del hueso, por lo muchos investigadores han invertido gran cantidad de tiempo en este tema. El primero en ofrecer una ley de comportamiento óseo fue Wolff (1892), quien sostuvo que el hueso debe su forma, densidad y propiedades a un proceso evolutivo de manera que éste se forma y transforma según las cargas a las cuales está sometido, hasta poder soportar dichas cargas de una forma óptima.

Avances recientes en las mediciones de la densidad ósea han permitido determinar la relación que existe entre la densidad y las propiedades mecánicas del hueso. Se ha realizado gran cantidad de investigación en torno a relacionar estas dos propiedades por medio de las tomografías axiales computarizadas (TAC).

En esta sección se hace primero, un recuento de lo que se ha publicado en el modelaje óseo utilizando tomografías computarizadas, y luego un recuento de como se ha utilizado esta información para realizar modelos utilizando el método de los elementos finitos.

Alguno de los primeros investigadores que publicaron sobre este tema son Weaver y Chalmers (1966), quienes investigaron la influencia de la densidad ósea relacionada con la edad de diferentes pacientes y el esfuerzo de fractura a compresión del hueso trabecular de las vértebras. Ellos demostraron que el esfuerzo de fractura está relacionado con el contenido de mineral del hueso y es independiente de la edad y del sexo del paciente. Luego, Mc. Elhaney et. al. (1970) observaron una relación positiva entre la densidad aparente del hueso y el esfuerzo a compresión y módulo de elasticidad.

Carter y Hayes (1977) publicaron una relación entre la densidad aparente del hueso y las propiedades mecánicas que todavía es vigente y está referenciada en casi todas las publicaciones posteriores. Ellos hicieron pruebas con probetas de hueso trabecular humano y bovino, y sugirieron, como otros investigadores previos, que el hueso se puede clasificar como compacto o como trabecular dependiendo de su porosidad que es proporcional al volumen ocupado por el tejido no mineral. Por esta razón ellos afirmaron que las propiedades del hueso trabecular son similares al hueso compacto cuando se observa en forma microscópica. Ellos encontraron que el esfuerzo a compresión para todo tejido óseo es aproximadamente proporcional a la velocidad de deformación a la potencia de 0.06 y al cuadrado de la densidad aparente del hueso

             (1)

donde 

                         S = esfuerzo de compresión (MN/m²)
                         r = densidad aparente (gr/cm³)
                           = velocidad de deformación (1/seg)
                      S_c = esfuerzo de compresión del hueso compacto con una densidad r_c a una tasa de deformación de 1.0 1/seg
También presentaron una relación similar entre el módulo de elasticidad y la densidad aparente

       (2)

donde

              E = Módulo de Elasticidad a la compresión (MN/m²)
             Ec = Módulo de Elasticidad a la compresión del hueso compacto con una densidad r_c a una tasa de deformación  de 1.0 1/seg

Estas relaciones son válidas para cualquier rango de densidad ósea en el esqueleto, para todo hueso compacto y casi todo el poroso hueso trabecular.

Posterior a este trabajo, ha aparecido un significativo número de publicaciones (entre otros Townsend et al. 1975, Neil et al. 1983), que comprueban estas ecuaciones y agregan que se pueden ver afectadas por la estructura ósea, las propiedades microscópicas, la orientación de las trabéculas y la anisotropía del hueso. Estos investigadores también enfatizaron sobre la importancia de calibrar adecuadamente las tomografías computarizadas, ya que se puede obtener entre un 30% y un 40% de error al relacionar la densidad aparente que se lee en las imágenes y las propiedades del hueso.

En este sentido, McBroom et al. (1985) realizaron estudios cuantitativos de tomografías computarizadas in vitro de manera de investigar la precisión que esta técnica ofrece. Sus resultados demostraron que las TAC ofrecen una información precisa sobre las propiedades mecánicas, siempre y cuando éstas hayan sido correctamente calibradas.

Más tarde, Rice et. al (1988) demostraron que la constante de proporcionalidad entre el módulo de Young y la densidad aparente del hueso difiere entre el humano y el bovino. Por otra parte afirmaron que la sugerencia de Wolff (1892) en cuanto a que las propiedades del hueso compacto son las mismas que el hueso trabecular, es solo cierta si se consideran factores de composición química, geometría, y no así cuando solo se consideran las propiedades mecánicas. Ellos hicieron un estudio estadístico entre todos los datos de tomografías presentes en la literatura y expusieron que para el hueso trabecular una relación cuadrática representa mejor los datos que una representación cúbica, además de que la dirección tiene un papel importante como muestra la siguiente ecuación:

        (3)

donde
                              D=1 para la dirección transversal, D=0 para la dirección longitudinal
                              T=1 para tensión, T=0 para compresión

Rice et. al. (1988) también trabajaron con el módulo de elasticidad para las trabéculas individuales y aseguraron que éste no es igual al módulo de elasticidad del hueso cortical.

Ahsman et. al. (1989) investigaron la variación anatómica del módulo de elasticidad ortotrópico de la tibia humana utilizando técnicas de ultrasonido. Ellos encontraron que el hueso trabecular de la parte proximal de la tibia resulta muy heterogéneo, donde el módulo axial varía desde 340 a 3350 MPa, y su anisotropía es constante.

Keyak et. al. (1994) consideraron que era necesario unificar todas las relaciones entre la densidad, el módulo de elasticidad y el esfuerzo del hueso trabecular que hasta entonces se habían publicado, y llegaron a la conclusión de que todas las ecuaciones existentes en la literatura proporcionan la misma relación luego de ser convertidas a una misma base y reafirmando que las TAC son apropiadas para determinar las propiedades mecánicas siempre y cuando se calibren adecuadamente las imágenes. Existe una gran cantidad de publicaciones en la década de los 90 (Rohl et al. 1991, Hvid et. al. 1985, 1989, Keaveny 1994) que evidencia que la falla por tensión o compresión de las trabéculas es independiente a la densidad o al módulo de elasticidad, e isotrópico independiente de la orientación de la trabécula, sugiriendo que la falla de la trabécula es fundamentalmente función de la deformación a la cual el hueso está sometido y no al esfuerzo.

En los últimos años los esfuerzos se han dirigido al modelado microscópico de las trabéculas para el modelado más preciso del hueso, se ha invertido gran cantidad de tiempo en la evaluación de la fracción de volumen del hueso trabecular medida por medio de TAC microscanning. Entre los investigadores que están desarrollando esta área están Ding. et al. (1999), Kuhn et. al. (1990) Rüegsegger et. al. (1996), Müller. et. al. (1996).

Se han reportado varios modelos de huesos utilizando el método de elementos finitos, un buen resumen de los primeros trabajos está reportado en el artículo de Huiskes y Chao (1983). Otras referencias de los modelos que utilizan el método de elementos finitos combinado con el uso de tomografías axiales computarizadas son: Huang et.al (1980), Steiz y Rüegsegger, (1983), Rhodes et. al. (1985) entre otros.

Las publicaciones más recientes intentan utilizar las tomografías computarizadas para hacer modelos de elementos finitos personalizados para cada paciente. Esto todavía no se ha logrado de forma automática. Por el contrario, realizar una malla 3D donde se reproduzca las complicadas superficies del hueso requiere un enorme esfuerzo manual.

Entre los avances que se han realizado en esta dirección se pueden nombrar el trabajo de Keyak, et. al. (1990), quienes presentaron un método basado una caja de pixels o "arreglo de propiedades" (voxel based method), para generar un modelo 3D para pacientes específicos, donde tanto la geometría como las propiedades del hueso son extraídas de tomografías axiales computarizadas. El método de ellos utiliza elementos cúbicos de ocho nodos con funciones de interpolación lineales, tres grados de libertad translaciones por nodo y propiedades de material isotrópicas. Todas las celdas cuadriculadas son alineadas con la dirección de la tomografía, de manera que el material del elemento está directamente relacionado con la densidad radiológica. Todos los elementos están direccionados de la misma forma. En este trabajo, no pretendieron simular la superficie del hueso exactamente, por el contrario éste resulta bastante impreciso por causa de la forma cuadriculada de los elementos. Ellos utilizaron la formula de Carter y Hayes (1977) para relacionar las imágenes con las propiedades mecánicas del hueso.

En esta misma dirección Marom y Linden (1990) utilizaron datos de TAC para generar un modelo tridimensional del hueso. Ellos reconstruyeron las secciones transversales de imágenes pixel por pixel, de manera de crear una matriz tridimensional de información de las TAC para obtener un acceso directo a los datos estructurados. La malla de elementos finitos 3D se creó a partir de estas secciones transversales bidimensionales, de manera que en cada sección transversal se creó una malla con igual número de elementos que en la siguiente y así al conectar las dos secciones transversales se crea el modelo tridimensional. Luego se aplica una subrutina para asignar la propiedad del elemento como un promedio de las propiedades de los pixels involucrados en ambas secciones transversales. Ellos utilizaron el software de Elementos Finitos SAP IV y el ANSYS como procesador y post procesador.

Kullmer et. al. (1998) desarrollaron diferentes subrutinas para el modelado óseo a partir de TAC, proponiendo dos métodos: el primero en el cual se extrae la superficie de las imágenes y se malla el sólido con un software comercial y el segundo donde se construye una malla de hexaedros donde cada pixel representa un elemento. En este último caso, ellos generaron un algoritmo para reducir el número de elementos y suavizar la superficie. Cada uno de los métodos tiene ventajas y desventajas, y se utilizan dependiendo de las necesidades que requiera el modelo.

Viceconti et. al. (1998) propusieron unos lineamientos para un adecuado protocolo de adquisición de imágenes para el modelado de huesos utilizando el MEF. Con este procedimiento se obtienen mayor cantidad de imágenes donde el hueso es mas complicado y menor cantidad donde tanto la geometría y la composición del hueso es mas sencilla.

Zannoni et. al. (1998), desarrollaron algoritmos para la generación de malla de elementos y la adquisición de propiedades mecánicas a partir de las tomografías computarizadas. En el software desarrollado, los elementos no tienen que estar necesariamente direccionados paralelos al eje de las tomografías y esto da mas libertad al programa. Ellos también desarrollaron una técnica para asignar propiedades mecánicas a los elementos (constantes en todo el elemento) a partir de los grises de las imágenes y utilizando las relaciones propuestas por Carter y Hayes (1977). Con esta metodología lograron colocarle 624 materiales diferentes al modelo del hueso, reduciéndolo a 214 con una diferencia del 1% en la energía de densidad de deformación.

Kerner, et. al. (1999) construyeron un modelo tridimensional de elementos finitos a partir de una densitometria ósea y realizaron remodelación ósea a partir del modelo. El modelo estaba constituido por elementos tridimensionales isoparamétricos de 8 nodos y las propiedades fueron asignadas elemento a elemento como el promedio del valor de las propiedades mecánicas proveniente de los grises correspondientes.

Van Rietbergen et. al. (1999) desarrollaron un modelo de la parte proximal de un fémur canino utilizando micro tomografías. Esto fue traducido en un modelo de micro elementos finitos con 7.9 millones de elementos adecuadamente refinados para representar las trabéculas. Utilizando un software especial de elementos finitos realizaron análisis para distintos tipos de carga, de manera que superponiendo los resultados se pueden obtener los esfuerzos y las deformaciones para otras cargas diferentes. Otros análisis de los resultados permiten ver con este método el volumen de interés, en la trabécula, localizado en la parte esponjosa del hueso. El objetivo de este proyecto fue proporcionar estimados razonables en esfuerzos y deformaciones de las trabéculas y comparar los resultados con los esfuerzos que deberían aparecer en las trabéculas según la hipótesis de Wolf (1982), que dice que la arquitectura de la trabécula es tal que el esfuerzo mínimo debe aparecer donde exista el mínimo peso. Estos mismos investigadores continuaron su trabajo y luego van Rietbergen et. al (2000) publicaron modelos de hasta 97 millones de elementos, con el objetivo de comparar el comportamiento de las trabéculas de un hueso normal con un hueso osteoporoso.

El uso de TAC se ha convertido en una práctica común en todos los estudios de Biomecánica. El concepto de tomografía fue introducido a la radiología durante los años treinta. Si bien la radiología convencional produce imágenes bidimensionales de un objeto, la adquisición de imágenes tomográficas se realiza por medio de una rotación que secciona al objeto y lo organiza en imágenes paralelas sucesivas. La cantidad y la calidad de información que contiene cada valor de gris de una sección, denominada resolución espacial, incrementa con el número de mediciones de atenuaciones en diferentes ángulos. Las imágenes están formadas por cierta cantidad de elementos o pixels que al tomar en cuenta la distancia entre las imágenes se puede hablar de elementos de volumen o voxels de la imagen. A cada voxel se le asigna un valor numérico que representa el valor de atenuación, que corresponde al promedio de la irradiación absorbida por el tejido en ese pixel, por lo cual la densidad de la tomografía computarizada es directamente proporcional al coeficiente de absorción. El valor o tono de gris de cada pixel se expresa en unidades de Hounsfield (HU). Cuando el tomógrafo esta calibrado, la densidad tomográfica del agua es considerada como 0 HU y la del aire como -1000 HU. A cada uno de los diferentes tejidos del cuerpo en estudio se le asignará un valor relativo a la escala de Hounsfield que se muestra en la figura 3.

Las imágenes salen del tomógrafo en formato DICOM¹ que contiene la información en unidades Hounsfield, están son traducidas en muchos casos en formato TIFF para ser leídas por programas comerciales, pero en la transformación de las imágenes se piede mucha información por lo que es preferible trabajar con las imágenes en su formato original.

El valor de atenuación del agua (0 HU) y del aire (-1000 HU) representan puntos fijos en la escala de densidad de las TAC y no son afectadas por el voltaje del tubo irradiador del tomógrafo. Dependiendo de la efectividad del aparato irradiador, la relación de los distintos tipos de tejidos y fluidos respecto al agua variarán. Por esta razón, los valores de densidad presentados en la literatura son solo una referencia, la cual implica que una calibración del tomógrafo es indispensable para obtener datos coherentes.

Cann y Genant (1980) desarrollaron un proceso de calibración por medio de un dispositivo con cilindros de hueso de densidad conocida el cual se digitaliza con el cuerpo que se desea estudiar, ellos encontraron que las densidades óseas conocidas se relacionaban con una escala lineal a los números de Hounsfield. Las TAC utilizadas en este trabajo corresponden a un paciente masculino de 24 años de edad, fueron adquiridas en la ciudad de Caracas en la Clínica Felix Boada, Baruta, por un radiólogo especializado. El tomógrafo utilizado tiene las siguientes características: Toshiba AUKLET, 1999, Helicoidal con espesor de corte de 3 mm a 6mm. Velocidad de arranque de la mesa igual a f=1/1.8 seg. Las imágenes fueron reconstruidas con interpolación en intervalos de 2, la velocidad de mesa (pitch) fue de 6 mm/s, filtro de alta resolución FCO2, ángulo de colimación² de 30º. La resolución de la imagen es de 0.5 mm por pixel. Cada punto del fantoma de 15 mm de diámetro tiene una densidad ósea y un número de Hounsfield (HU) conocido, lo cual proporciona una curva de calibración, que para este caso es: y = 0.69767 x+10.8347287 (mg/cm³), como muestra el gráfico de la figura 4.

La conversión de TAC en modelos tridimensionales es un proceso complicado y requiere varios pasos. Cada uno de estos pasos representa una potencial fuente de error si no se realiza adecuadamente. El proceso de reconstrucción ósea depende de la utilidad que se le vaya a dar al modelo. En este proyecto se generaron modelos para el cálculo de deformaciones y tensiones utilizando el Método de Elementos de Contorno (MEC) y el Método de Elementos Finitos (MEF). Para el primero, es necesario crear un modelo superficial y para el segundo es necesario crear un modelo sólido, que serán discretizados por programas comerciales. En otros casos, se crea la malla sin antes construir unos modelos completos del hueso que se desee modelar, lo que requiere un proceso manual que es el que se utiliza para el análisis con el método de elementos finitos. En todos los casos se requiere reproducir inicialmente la geometría externa basándose en TAC tomadas a distintas distancias, creando un modelo de "alambre". El proceso de adquisición de TAC para la reconstrucción ósea debe ser debidamente realizado para obtener resultados confiables.

Para la reconstrucción de la tibia en este trabajo se digitalizó el área de la rodilla de 360 mm, por medio de 180 imágenes adquiridas a un milímetro de espesor en la parte proximal de la tibia, donde existe una fuerte presencia de hueso trabecular y 4 mm en la diáfisis, como muestran las figuras 5 y 6.

Figura 5. Proceso de segmentación de la tibia humana

La segmentación es el proceso de subdividir la imagen en regiones basándose en sus características pictóricas. La subdivisión se hace aislando los voxels correspondientes a un rango de unidades Hounsfield deseado. Ésto puede observarse en la figura 5 donde se muestra la segmentación de la tibia. La posición de cada voxel es medida desde el primer voxel de la imagen. Se pueden utilizar tres diferentes métodos en el proceso: fijación de límites de valores superiores e inferiores, detección de bordes y expansión de regiones a partir de una semilla. La segmentación se puede hacer de forma manual en programas ingenieriles CAD/CAM o se pueden utilizar programas comerciales que extraen los contornos para todas las imágenes secuenciales. Existen varios programas comerciales, el utilizado en este trabajo es el Surfdriver (Moody y Lazonoff 1999) que proporciona los contornos en polilíneas en formato IGES³ . En muchos casos, y sobre todo en la zona de las articulaciones donde existe gran cantidad de hueso esponjoso, las imágenes son difusas, los contornos deben ser cerrados y perfeccionados manualmente lo cual representa un proceso laborioso.

Los contornos son importados en programas ingenieriles de simulación CAD/CAM, donde las polilíneas son convertidas en curvas suaves (splines). Esto es lo que se denomina el modelo de "alambre" (wireframe model) como muestra la figura 6. Para que las superficies óseas queden suaves es importante reducir el número de puntos de control de los splines y comprobar que todos tengan un número similar de puntos o aumentando en forma creciente. Utilizando las herramientas de los programas CAD/CAM que permiten crear planos de trabajo (working planes) y perfiles (profiles), estas curvas son transformadas en secciones transversales planas, las superficies son creadas a partir de estos contornos cerrados (surface by boundaries) y los sólidos son creados por extrusión (lofting) (ver figura 7).

                (a) Modelo superficial                         (b) Modelo sólido
Figura 7. Generación de modelos óseos utilizando herramientas CAD/CAM

Existen varias formas de generar un modelo discreto, y la escogencia del tipo de malla dependerá del método numérico que desee utilizar.

El primer modelo que se generó en este trabajo, se hizo utilizando el MEC para el cual es necesario crear una malla superficial. Se utilizaron elementos cuadriláteros generados automáticamente por un programa comercial (Rutina FEM de Pro/Engineer) a partir de un modelo superficial de la tibia (figura 8a), los cálculos fueron realizados con un programa creado en la Universidad de Londres (Müller-Karger et. al, 1999). Este modelo tiene la desventaja que se modeló como una sola región, modelando la médula como un hueco, por lo que no incluye los diferentes materiales que constituyen el hueso y fue superado por otros análisis realizados posteriormente.

La segunda simulación de la tibia que se realizó fue mallada a partir de un modelo sólido, por medio de tetraedros con el programa comercial Pro/Engineer. El análisis de esfuerzo deformación se realizó utilizando el MEF con el programa Pro/Mechanica (ver figura 8b). Este modelo tiene la desventaja que a pesar de que se pueden caracterizar distintas regiones para asignar diferentes propiedades a zonas determinadas del hueso, no permite modelarlo como un cuerpo altamente heterogéneo.

En ambos casos mencionados, la superficie del hueso es discretizada de una forma muy aproximada donde los elementos tienen caras planas y la superficie no es suave, situación que se agrava en la epífisis de los huesos, donde las superficies son más complicadas y el hueso es muy heterogéneo.

El mejor método para simular huesos resultó ser el MEF en su versión p adaptable. Este método, combinado con funciones de transformación isoparamétricas adecuadas, permite el modelado de la superficie del hueso de forma muy precisa y suave (ver figura 8c). Para generar la malla, se parte del modelo de contorno y se crean macro elementos manualmente. En este caso se utilizó un programa que acopla el modelo realizado en el programa comercial Mechanical Desktop con un programa de cálculo de Elementos Finitos en versión p adaptable desarrollado en la Universidad Técnica de Munich (Düster et. al 2000), en el Departamento de Informática en Ingeniería Civil. Este programa, en lenguaje C++, llamado AdhoC, esta diseñado para correr solo en plataformas UNIX. Para obtener la información geométrica se utiliza la interfase desarrollada para Mechanical Desktop (2000) que funcionan en máquinas PC bajo ambiente Windows NT (Bröker, 2001). Este acoplamiento requiere que las máquinas trabajen paralelamente utilizando modeladores de UNIX en ambiente Windows. Comercialmente existen programas que crean los modelos discretizados a partir del modelo geométrico, pero normalmente estos malladores son de elementos tetraédricos que necesariamente poseen una propiedad constante para todo el elemento y no reproducen exactamente la geometría del modelo original. La principal ventaja del acoplamiento utilizado en este trabajo es la capacidad de cálculo simultaneo del sistema CAD con el programa de MEF, lo cual implica la transferencia directa de la geometría, condiciones de carga y de contorno del modelo a cada macro elemento sin necesidad de extraerlo de la fuente donde se generó la geometría. Lo dicho anteriormente implica que los datos necesarios para el cálculo con el MEF vienen sin la pérdida de información y sin la simplificación de las superficies indispensable para otros programas. Por otra parte se consigue transferir cualquier cambio realizado en el modelo en el sistema CAD directamente al programa de MEF. Este modelo tiene la ventaja que como estos elementos deben ser integrados con funciones de orden superior, cada punto de Gauss puede ser relacionado con las características mecánicas del hueso, permitiendo 1000 ó más propiedades materiales diferentes para cada elemento a partir de la matriz de propiedades proveniente de las TAC, simulando así un material totalmente heterogéneo. La desventaja de este método es que la malla debe realizarse de forma manual ya que no existe un mallador de hexaedros comercial y esto es un proceso laborioso. Para generar estos modelos se desarrollaron técnicas que facilitan el trabajo manual.

Otros autores han generado mallas manualmente creando pequeños elementos hexaédricos por cada voxel de las TAC, con la ventaja de que cada voxel representa un sólido rectangular con dimensiones y propiedades constantes y con la desventaja que para TAC de alta resolución estos modelos resultan en un gran número de elementos finitos y las superficies no resultan suaves. Algunos de estos modelos se muestran en la figura 9. Algunos autores han creado algoritmos para suavizar las superficies y disminuir el tamaño de las mallas, entre ellos se pueden nombrar Keyak et al (1990), Kullmer et. al (1998) y Zanonni et. al. (1999) entre otros.

Tendencias más novedosas (Garcia, 1999) se están investigando actualmente relacionadas con los métodos sin malla, que trabajan directamente sin elementos, utilizando únicamente una nube de puntos que reproduce la geometría del sólido en estudio.

El método de elementos finitos arroja una solución al problema que difiere de la solución exacta por un error de aproximación. Este error disminuye de tamaño si el tamaño de la subdivisión (h) se reduce o si el orden del polinomio (p) de la función que describe el elemento crece. La variación de estos parámetros hasta conseguir una solución aceptable es lo que se denomina métodos adaptables (Zienkiewicz y Taylor, 1994, Szabo y Babuska, 1991, Bathe K.J., 1990).

La versión p del método de elementos finitos ha comprobado ser una excelente herramienta para resolver ecuaciones lineales elípticas, como son la ecuación de Poisson, las ecuaciones de Lamé, etc. Muchos autores han demostrado que esta versión es superior al método de elementos finitos clásico o la versión h (Düster et. al., 2000). Este método, combinado con un diseño apropiado de la malla, puede presentar una tasa de convergencia exponencial, lo cual permite llegar a una solución más rápidamente. Por otro lado la elaboración de los macro elementos tipo p utilizando métodos de fusión de funciones (blending function method) permite describir los contornos reales del sistema de una forma más precisa sin tener que aumentar el número de elementos. Esto permite discretizar estructuras complejas utilizando una malla burda como base para el uso de todo el espectro de funciones de interpolación del cual se hará referencia mas adelante.

La implementación de la versión p del método de elementos finitos está basada en conjuntos de funciones jerárquicas de una dimensión. La principal diferencia este tipo de funciones respecto a las funciones de interpolación tradicionales es que todas las funciones de menor orden están contenidas en las funciones de orden superior. Este tipo de funciones se basa en que para generar un polinomio de orden p a lo largo de un elemento no es necesario introducir nuevos nodos, sino que se pueden utilizar parámetros sin significado físico, las cuales permiten aumentar el grado de polinomio sin necesidad de realizar todos los cálculo de nuevo. Según Szabó y Babuška (1991), las funciones jerárquicas básicas se pueden implementar para cualquier grado de polinomio deseado, la figura 10 muestra el ordenamiento de las funciones jerárquicas utilizadas.La implementación de la versión p utilizada en este trabajo permite resolver problemas estructurales tridimensionales que pueden contener cualquier tipo de superficies casi de forma arbitrarias, como son los huesos. La implementación está basada en una formulación de hexaedros con la capacidad de variar el grado del polinomio en las tres direcciones locales así como el espectro de las funciones de interpolación para las tres componentes del desplazamiento.

Figura 10. Ordenamiento de las funciones jerárquicas

Una diferencia importante entre la versión h y p del método de elementos finitos radica en los requerimientos de las transformaciones de los elementos del sistema global al sistema local. Debido a que en la versión p el tamaño de los elementos no se reduce a medida que se incrementa el grado del polinomio, la descripción de la geometría es independiente del número de elementos y esto trae como consecuencia la necesidad de construir elementos que siguen la forma exacta de la superficie de la estructura a modelar. La transformación de elementos puede ser visto como un caso especial cuando se utiliza el método de las funciones de fusión (Blending function method, Dey et. al. 1997), a diferencia de las transformaciones isoparamétricas que se utilizan comúnmente en el MEF clásico, ya que con este método se pueden crear elementos cuyas caras pueden ser superficies casi de forma arbitraria. En la figura 11 se muestra como dibujar en un espacio cartesiano x,y,z, la forma distorsionada de un elemento en coordenadas locales x ,h ,z  por medio de una transformación biunívoca Q:

            (4)

La función de fusión utilizada en este trabajo es la desarrollada por Királyfalvi y Szabo (1997) en la cual además de transformar las coordenadas, transforma los bordes y las caras de los elementos. Considerando el elemento de la figura 11, se definen:
X_i=(X_i,Y_i,Z_i) con i=1,...8 , las coordenadas globales de cada nodo
E_i=(E_ix,E_iy,E_iz) con i=1,…12 , funciones que describen las formas de los bordes
F_i=(F_ix,F_iy,F_iz) con i=1,...6 , funciones que describen la forma de las caras
La función de fusión Q^e desde las coordenadas locales a las coordenadas globales x^T=(x,y,z) es:

    (5)

El primer término de esta ecuación corresponde a la transformación estándar de los nodos, el segundo término corresponde a la transformación de las caras y el tercer término corresponde a la transformación de los bordes de un elemento. Mas detalle de estas funciones e_i y f_i se encuentra en la referencia Düster et. al (2000).

La integración de las expresiones de las matrices de los elementos debe hacerse en forma numérica, para lo cual existen muchas técnicas entre las cuales se puede nombrar: la cuadratura de Newton-Cotes y la cuadratura de Gauss.

En este trabajo, el mismo número de puntos es utilizado en todas las direcciones (ver figura 12). Las propiedades mecánicas del modelo se asignan a cada uno de estos puntos en el proceso de calculo. Esto da la posibilidad de simular modelos altamente heterogéneos utilizando hasta 1000 materiales diferentes en cada elemento, dependiendo del número de puntos de Gauss utilizados.

Varios autores han encontrado relaciones entre la densidad del hueso y el módulo de elasticidad como muestra la tabla 2. Existen distintas formas de medir la densidad del hueso, entre ellas se pueden nombrar: densidad seca (dry density, ρ_dry )⁴ , densidad húmeda (wet density, ρ_wet )⁵ , densidad de ceniza (ash density, ρ_ash )⁶ y densidad aparente de la TAC (ρ_ct ), por lo cual es difícil comparar los resultados de todas las relaciones existentes en la literatura. Keyak et. al. (1994) presentaron ecuaciones matemáticas que relacionan todos los tipos de densidades, permitiendo unificar las relaciones entre propiedades mecánicas y densidad del hueso, y comparar resultados. Ellos llegaron a la conclusión que para las mediciones donde se había calibrado el tomógrafo las ecuaciones unificadas arrojaban valores muy parecidos.Los estudios no calibrados arrojaban variaciones considerablemente grandes, por lo cual no fueron incorporados en la estandarización de los datos. La relación entre las unidades Hounsfield (HU) o valores de grises y la densidad aparente del hueso en cada píxel, se calcula por medio de una recta que el usuario suministra de acuerdo a la curva de calibración del tomógrafo (ρ=C*HU+B).En esta investigación se utiliza la ecuación cúbica desarrollada por Carter y Hayes (1977), por ser ésta la más utilizada en la literatura existe una buena posibilidad de comparar los resultados. Asumiendo una velocidad de deformación de 0.01 mm/mm.seg (Cowin, 1989), la densidad del hueso compacto igual a 1.73 g/cm³ y el módulo de Young para ese hueso de 22 GPa (Zanoni et al. 1998) se obtiene:

                                                            E = kr³     con  k=4.249 [GPa/(g/cm³)] (6)

la cual es válida para todo el esqueleto, tanto para hueso cortical como para hueso esponjoso.

Se desarrolló un software en lenguaje C++, al cual se le asignó el nombre Greyvalues y tiene por finalidad traducir un conjunto de imágenes, en formato de tomógrafo (DICOM) o en formato TIFF, en valores ASCII que representen propiedades mecánicas. El programa permite tener mayor concentración de imágenes en algunas partes del estudio y menor concentración en otras partes donde se necesite menos información. En la entrada de datos del programa, las imágenes se clasifican de acuerdo a grupos de igual espesor (concentración de imágenes). El programa interpola la información necesaria para que el archivo de salida contenga información equidistante, a la menor distancia entre los grupos de imágenes suministradas.

Para ilustrar la capacidad del programa Greyvalues, se incluye aquí un ejemplo del estudio de una rodilla, para el cual se digitalizaron imágenes tomográficas de la siguiente forma: Fémur: 120 mm cada 4mm=30 imágenes, Rodilla:120 mm, cada 1mm = 120 imágenes, Tibia: 120 mm cada 4mm = 30 imágenes.

La imagen completa tiene 512x512 pixels, que corresponden a 255x255 mm. Se desea estudiar la tibia derecha por lo que se escoge una ventana de 150x150 pixelsl ubicada en el pixel (40,150) como muestra la figura 13. Para este ejemplo solo se analizaron las últimas dos imágenes del estudio, lo cual representan las imágenes ubicadas a 356 y 360 mm desde la primera imagen. La curva de calibración para este estudio es:

                                            r = 0.69767 (HU)+10.8347287 (mg/cm³) (7)

Los nombres de las imágenes a estudiar son ct356 y ct360, el encabezado de las imágenes es de 1608 bits, la posición de la caja de propiedades tiene el mismo origen (0.0,0.0,0.0) que el modelo de Mechanical Desktop. Se desea una salida de cuatro imágenes separadas por 1 mm. Existe un solo grupo de imágenes a procesar compuesto por dos imágenes separadas por 4 mm. La salida incluye todos los archivos que contienen información sobre unidades Hounsfield, densidad y módulo de elasticidad.

Las figuras 14 y 15 muestran la salida para la imagen ct356 de la densidad y del Módulo de elasticidad respectivamente, los valores de la densidad varían en un rango de -0.908 a 1.6797 gr/cm3 y los valores de módulos de elasticidad del hueso varían entre 0 y 20.135 GPa; por razones físicas los valores negativos del módulo de elasticidad son forzados a cero. Es de notar que el rango comienza con números negativos ya que esta calibración es válida exclusivamente para componentes óseos y en la imagen existen también tejidos bandos que no están adecuadamente representados con esta escala.

El programa crea un archivo de reporte, con la información de lo que se ha calculado, que archivos se han creado y donde ha almacenado la información. El nombre de este archivo es suministrado por el usuario con la variable IDout, la tabla 4 muestra este archivo para el ejemplo de corrida.

El programa de Elemento finitos en versión p adaptable utiliza el arreglo de propiedades suministrado y recoge la propiedad del hueso en los puntos de integración, de manera que cada punto de Gauss tendrá un módulo de elasticidad diferente dependiendo de la heterogeneidad del cuerpo que se este analizando.

6.1 Homogeneización del arreglo de propiedades mecánicas

Una vez ensamblado el arreglo de propiedades, se procede al cálculo de esfuerzos y deformación con un software de Elementos Finitos en versión p. Para comprobar la veracidad de los resultados y evaluar el beneficio de la complejidad del modelo, se realizó un estudio paramétrico, en el cual se varió la cantidad de valores suministrados en el arreglo de propiedades mecánicas según el número de pixels evaluados y se homogeneizaron las imágenes a través del uso de promedio de los valores de los pixels vecinos. Los resultados se midieron con curvas de esfuerzos de von Mises y energía de deformación.

         Figura 16. Comparación entre una imagen no homogeneizada y una imagen con 3 pixels promediados  

Para los primeros modelos se obtuvieron buenos resultados en las deformaciones pero saltos en los esfuerzos por lo cual se decidió homogeneizar la matriz de manera de que no existieran valores muy distintos en el módulo de elasticidad para puntos de Gauss vecinos. La homogeneización del arreglo se realizó haciendo un promedio entre el valor de cada pixel y los pixels vecinos. En el programa, el usuario puede introducir el número de pixels vecinos que se desea promediar y el peso que este promedio debe llevar respecto al valor original, por defecto se toma el promedio de todos los pixels incluyendo el que se desea reemplazar como nuevo valor. Las figura 16 y 17 muestran la influencia de promediar pixels en la imagen tomográfica.

Se presentarán resultados del análisis de esfuerzo de los modelos de la tibia, calculados considerando en primer lugar al hueso como un cuerpo homogéneo y en segundo lugar como un cuerpo heterogéneo, donde las propiedades son extraídas de TAC. Se describen las ventajas y desventajas del método desarrollado en esta tesis respecto a los métodos utilizados por otros autores. Con la metodología desarrollada se realizaron varios modelos para los cuales se muestra los resultados del análisis de esfuerzo y las curvas de convergencia

Los primeros resultados son de un modelo de tibia de 357 mm de longitud compuesto con un solo material y con una cavidad interna para simular la medula. Este modelo se realizó con la intención de demostrar la capacidad de modelar y analizar mecánicamente los huesos con superficies suaves y pocos elementos.

Las imágenes tomográficas utilizadas para generar este modelo provienen del "Visible Human Project" donde se tienen imágenes cada un mm de un hombre completo con una resolución de 0.5mm/pixel. Para este modelo se utilizaron 97 imágenes distribuidas con una mayor concentración donde la geometría es mas complicada (epífisis) y menor concentración en el eje de la tibia.

En este ejemplo se asumió un módulo de elasticidad promedio de 17000 MPa y un coeficiente de Poisson v =0.3. Se aplico una carga constante distribuida en los cóndilos de F=2450 N. La carga se colocó en los puntos nodales que cubren áreas aproximadamente igual a las áreas del contacto siguientes (Kettlekcamp y Jacobs 1972): de 468mm2 (platillo medial) y 297 mm2 (platillo lateral). La malla para este hueso resultó en 136 elementos hexaédricos. Los esfuerzos de von Mises se calcularon con grados de polinomio, que definen el orden del los elementos, variable de p=3 hasta p=6. La tabla 5 muestra el análisis para distintos grados de polinomio, demostrando que la energía de deformación no presenta grandes cambios, por lo cual se considera que el resultado es convergente.

En la figura 18 se muestran los resultados para p=6, con un esfuerzo máximo de compresión de 104 MPa que ocurre hacia el final del eje de la tibia a una distancia de la parte proximal de 330 mm; el eje del hueso presenta una distribución suave de esfuerzos. Estos resultados concuerdan con los resultados presentados por Metha y Rajani (1995) en donde ellos obtuvieron una deflexión de 7.954 mm y un esfuerzo máximo de 147 MPa para un modelo compuesto por tres materiales. En el modelo presentado en esta investigación el esfuerzo y desplazamiento resultaron algo menores ya que se consideró la tibia sólo compuesta por hueso compacto y una cavidad para la medula. A pesar de que este modelo ya representa un avance respecto a los modelos presentes en la literatura en cuanto al modelado de las superficies óseas, ya que son mas suaves y reproducen mejor al hueso real, éste no pudo ser analizado con el arreglo de propiedades debido a su tamaño y a la capacidad computacional de la que se disponía. El arreglo de propiedades para la tibia completa tiene un tamaño de aproximadamente 200 Mb.

Otro modelo de tibia reportado en la literatura es el de Bedzinski y Scigala, (2000). Ellos presentan un modelo completo de tibia compuesto por diferentes regiones simulado con tetraédros, utilizando el MEF. Este modelo presenta su máximo valor de esfuerzo en la parte inferior a 300 mm de la articulación con un valor de 84 Mpa. Este valor no puede ser comparado cuantitativamente con los valores presentados en este trabajo ya que ellos no especifican el valor de las cargas utilizadas, pero cualitativamente se observa un comportamiento similar a los resultados obtenidos en este trabajo.

El segundo análisis que se presenta es un modelo de la diáfisis de la tibia 144 mm discretizada en 48 macro elementos, se comparan los resultados para la misma geometría utilizando, en el primer caso un modelo sólido compuesto por un solo material (E= 17200M Pa, v=0.3), en el segundo caso un modelo sólido heterogéneo utilizando tantos materiales, extraídos del arreglo de propiedades, como puntos de Gauss. En este ejemplo se utilizó diez puntos de Gauss en cada dirección, lo que representa 1000 puntos de integración y mil materiales diferentes para cada elemento. El área de la superficie superior es de 937,88 mm², y se le aplicó una tensión de 2.61N/mm².

Las figuras 19 y 20 muestran una gran diferencia entre los resultados para el modelo homogéneo y el modelo heterogéneo en cuanto a los esfuerzos de von Mises. En primer lugar el valor máximo de esfuerzo es mucho más alto para el caso heterogéneo y en segundo lugar cabe destacar la diferencia en la distribución de los esfuerzos. El caso homogéneo se comporta como una viga a compresión y a flexión donde el esfuerzo máximo ocurre en la superficie cerca del extremo vinculado. En el caso heterogéneo los esfuerzos se distribuyen de forma más homogénea en la zona periférica donde existe mayor densidad ósea, como era de esperarse.

En la figura 21 se muestra un análisis de convergencia para este modelo, para el cual se varía el grado del polinomio de 1 a 8. El gráfico de energía de deformación demuestra que los resultados para grados de polinomio bajo no son de mucha utilidad sobre todo en el caso del modelo heterogéneo.

El esfuerzo de von Mises máximo para el caso heterogéneo incrementa de forma constante a medida que se afina el polinomio, lo cual da a pensar de que existe una singularidad que se calcula con mayor precisión cuando el orden de integración aumenta. Existe posibilidad que haya un cambio brusco de propiedades en la frontera del elemento cuando se utiliza el arreglo de propiedades o bien que existan puntos de Gauss vecinos con valores de propiedades mecánicas muy distintas. Estas dos situaciones son fuentes de imprecisiones ya que pueden producir saltos en los valores obtenidos del esfuerzo. Una análisis más representativo se puede realizar comparando, no solo los esfuerzos máximos, sino también los esfuerzos en puntos particulares de los modelos como muestra la figura 21. Se escogió un corte transversal del modelo y se escogieron tres puntos de esta sección arbitrariamente. Los puntos 1 y 3 fueron elegidos en la periferia y el punto 2 se colocó en la médula del hueso. La figura 22 muestra convergencia del esfuerzo de von Mises para todos los puntos a partir del grado de polinomio p=5. La figura 23(a) muestra los esfuerzos von Mises para el caso homogéneo y la figura 23(b) los esfuerzos von Mises para el caso heterogéneo. Se puede ver claramente que en el punto 2 en el caso heterogéneo los esfuerzos son muy pequeños debido que no existe material en esa zona.

Figura 22. Puntos de comparación de resultados

7.3 Influencia de la omisión de pixels y la homogeneización en el arreglo de propiedades

Finalmente se presentan los resultados de un modelo de la tibia sin cóndilos, omitiendo dos pixels para controlar el tamaño del arreglo de propiedades, se discretizó utilizando el método de las superficies. Se comparan los resultados para la misma geometría utilizando, en el primer caso un modelo sólido compuesto por un solo material (E= 17200M Pa, v=0.3), en el segundo caso un modelo sólido heterogéneo utilizando tantos materiales, extraídos del arreglo de propiedades, como puntos de Gauss. Las figuras 24 y 25 muestran los resultados en esfuerzo de von Mises para ambos casos y la figura 26 muestra la influencia de homogeneizar la matriz en el esfuerzo máximo de von Mises y en el tiempo de cálculo. Homogeneizar la matriz en este ejemplo produce los mismo resultados en que el ejemplo anterior.

Figura 24. Esfuerzos y deformaciones para tibia sin cóndilos utilizando el arreglo de propiedades promediando 5 pixels

Figura 25. Esfuerzos y deformaciones para tibia sin cóndilos utilizando solo un material

La metodología, logros y tropiezos presentados en esta investigación muestran la complejidad del problema que representa simular un sistema de geometría y composición tan complicada como son los huesos. El trabajo realizado se puede resumir en los siguientes aspectos:

* Se creo una metodología para reconstruir y simular huesos a partir de TAC. 

*Se utilizaron diferentes métodos numéricos, como son el método de elementos de contorno y el método de elementos finitos en versión clásica y versión p, para analizar mecánicamente la tibia humana, y determinar que método es más adecuado para la simulación ósea. 

* Se utilizaron programas comerciales para generar modelos discretos de la tibia con elementos cuadrilateros superficiales de ocho nodos para calcular esfuerzos y deformaciones por MEC, y utilizando tetraedros para cálculos con el MEF. Se crearon los programas necesarios para generar los datos de entrada para los programas de cálculo. 

* Se generaron modelos discretos con macro elementos hexaédricos por medio del modelaje manual con sistema CAD, para el cálculo de esfuerzos y deformaciones utilizando la versión p del MEF. 

* Se utilizó una interfase entre sistema CAD y el programa de MEF para la transferencia por medio de funciones de fusión de la geometría, directamente como entrada de datos a programa de cálculo por MEF en versión p, lo cual permite que las caras de los macro elementos sean exactamente la superficie suave de los huesos. 

* Se desarrolló un programa que lee las imágenes tomográficas, convierte los valores de grises en propiedades mecánicas del hueso y crea un arreglo de propiedades como entrada de datos para el programa de MEF en versión p.

De los aspectos desarrollados se puede concluir:

* El método que se desarrolló en este trabajo plantea un acercamiento diferente en la simulación ósea a los métodos encontrados en la literatura. 

* El método mas adecuado para simular los huesos resultó ser el MEF en versión p adaptable, ya que permite por medio de macro elementos simular los huesos de forma sencilla con pocos elementos. 

* La simulación de los huesos por medio de la versión p del MEF conjuntamente con la transformación de fusión de los elementos, permite la discretización de los sistemas con macro elementos compuestos por casi cualquier superficie, lo que conlleva a un modelo de los huesos con superficies suaves. 

* Uno de los principales aportes de este trabajo es que la metodología utilizada permite la asignación de las propiedades mecánicas del hueso a través de los puntos de integración de Gauss. En caso de que se utilicen 10 puntos en cada dirección esto representa 1000 materiales diferentes para cada elemento. Dependiendo del número de elementos utilizado en el modelo se puede simular la heterogeneidad del hueso de forma bastante precisa. Por lo dicho anteriormente, se puede afirmar que este método difiere de los métodos macroscópicos en que permite simular internamente al hueso de forma más continua y con una gran cantidad de materiales sin necesidad de utilizar intervalos. 

* Con esta metodología se logra determinar el comportamiento óseo de forma independiente a la arquitectura de las trabéculas y por tanto el análisis es independiente de la malla que se utilice. 

* Se simula el hueso como altamente heterogeneo e isotrópico, el próximo paso es considerar la anisotropía del hueso que se pretende resolver con un proceso iteratico, de manera de obtener un modelo cuyas lineas de fuerza sigan la estructura del hueso trabecular.

La investigación presentada en este trabajo refleja solo un pequeño porcentaje del avance que representa la simulación del tejido óseo. Las perspectivas de esta investigación son poder estudiar de forma virtual el beneficio de colocar una prótesis o algún otro dispositivo traumatológico en un paciente antes de la intervención, de manera que la cirugía sea mucho más exitosa luego de un estudio teórico previo. Por otra parte podría estudiarse como se va a remodelar el tejido óseo una vez colocada la prótesis, permitiendo estudiar cuanto tiempo durará el dispositivo en el paciente sin ocasionar problemas de aflojamiento, dolor y por consecuencia una nueva cirugía.  

1. Ashman R.B., Cowin S.C., van Buskirk W.C., Rice,(1984), A continuous wave technique for measurement of elastic properties of cortical bone , J. Biomechanics, 17:349-361.        [ Links ]

2. Ashman R.B., Rho J.Y., (1988), Elastic Moduli of trabecular bone material, J. Biomechanics, 21:177-181.        [ Links ]

3. Bathe K. J.,(1982), Finite Element Procedures in Engineering Analysis. Prentice Hall, Inc.        [ Links ]

4. Bedzinski R., Scigala K., (2000), Physical and numerical models of the tibia bone, 12th Conference of the European Society of Biomechanics, Dublín 2000, p 186.        [ Links ]

5. Bröker H, (2001) Integration von geometrischer Modellierung und Berechnung nach p-Version der FEM., PhD thesis, Teschnische Universität München.        [ Links ]

6. Buckwalter J.A., Glimcher M.J. Cooper R.R., Recker R., (1995), Bone Biology. Part I: Structure, Blood Supply, Cells, Matrix, and Mineralization, J. Bone and Joint Surgery, 77A(8) :1256-1275.        [ Links ]

7. Can C.E., Genant H.K, (1980), Precise measurement of vertebral mineral contnet using computd tomography, J. Computer Assisted Tomography , 4:493:500.        [ Links ]

8. Carter, D.R. and Hayes W.C., (1977), The compressive behavior of bone as a two-phase porous structure, J Bone Joint Surg, 59A, pp. 954-62.        [ Links ]

9. Comín M., Peris J.L., Prat J.M, Dejoz J. R., Vera P. M., Hoyos J. V., (1999), Biomecánica de la fractura ósea y técnicas de reparación. Instituto de Biomecánica de Valencia, ISBN 84-923974-5-4        [ Links ]

10. Cowin S.,(1989), Bone Mechanics,CRC Press Inc. Florida USA:        [ Links ]

11. Currey (1970) The Mechanical Properties of Bone, Clinical Orthopeadics and Related Research, 73:211:231        [ Links ]

12. Ding M., Odggard A. Hvid I., (1999), Accuracy of cancellous bone volume fraction measured by micro-CT scanning, J. Biomechanics, 32:323-326.        [ Links ]

13. Düster A., H. Broker, E. Rank (2000), The p-version of the finite element method for three-dimensional curved thin walled structures, submitted to: International Journal for Numerical Methods in Engineering.        [ Links ]

14. Garcia J.M., (1999), Modelos de remodelación ósea: análisis numérico y aplicaciones al diseño de fijaciones de fracturas del fémur proximal, Tesis de Doctorado, Centro Politécnico Superior de la Universidad de Zaragoza, España.        [ Links ]

15. Hvid I, Jensen N, Bunger C., Solund K and Djurhuus J., (1985) Bone mineral assay: its relation to the mechanical strength of cancellous bone. Engng Med., 14:79-83.        [ Links ]

16. Hvid I, Bentzen S.M., Linde F., Mosekilde L., Pongsoipetch B (1989), X-ray quantitative computed tomography: the relations to physical properties of proximal tibial trabecular bone specimens, J. Biomechanics 22:837-844.        [ Links ]

17. Huang HK., Suarez F., Toridis TG, Khonozeimeh K, Ovenshire L, (1980), Utilization of computerized tomographic scans as input to finite element analysis. International Conference Procedings: Finite Element in Biomechanics, 2:797-815        [ Links ]

18. Huiskes R., Chao Eye., (1983) A survey of finite element analysis in orthopedic biomechanics: the first decade, J Biomech.,16: 385-409.        [ Links ]

19. Keaveny T, Wachtel E.F., McMahon T.A.,and Hayes W.C., (1994) Trabecular bone exhibits fully linear elastic behavior and yield at low strains, J. Biomechanics 27, 1127-1136.        [ Links ]

20. Kerner J., Huiskes R., van Lenthe G.h., Weinans H., van Rietbergen B., Engh C.A., Amis. A.A., (1999), Correlation between pro-operative periprostthetic bone density and post-operative bone loss in THA can be explained by strain-adaptative remodeling, J. of Biomechanics , 32:695-703.        [ Links ]

21. Keyak J.H., Meagher J.M., Skinner H.B., and Mote C.D. Jr., (1990), Automated three-dimensional finite element modeling of bone: a new method, J Biomed Eng, 12:389        [ Links ]

22. Keyak J.H., Lee I.Y., and Skinner H.B., (1994), Correlations between orthogonal mechanical properties and density of trabecular bone: use of different densitometry measures, J Biomed Mat Res, 28:1329-1336.        [ Links ]

23. Kettlekcamp D.B., Jacobs A.W., (1972) Tibiofemoral Contact Area-Determination and Implantations, J. of Bone and Joint Surgery, 54A:167-172.        [ Links ]

24. Kuhn J.L., Goulet R.W., Pappas M., Goldstein S.A., (1990), Morphometric and anisotropic symentries of canine distal femur. J. Ortho. Res. 8:776-780.        [ Links ]

25. Kullmer G. Weiser J., Richar A,(1998), Construction of finite element models on the basis of computed tomography data, in Computer Methods in Biomechanics & Biomedical Engineering - 2. (J. Middleton, M.L. Jones, G. N. Pande, eds). Gordon & Bread Sci. Pub, The Netherlands, 279:287.        [ Links ]

26. Lacroix D, Prendergast P.J. Li G., Marsh D. (2000), A 3D finite element model of a tibia to simulate the regenerative and resorption phases during fracture healing. 12th Conference of the European Society of Biomechanics, Dublin, 2000, p.60        [ Links ]

27. Marom S.A., and Linden M.J., (1990), Computer aided stress analysis of long bones utilizing computed tomography, J. Biomechanics, 23(5): 399-404        [ Links ]

28. McBroom R.J., Hayes W.C., Edwards W.T., Goldberg R.P.,(1985), Prediction of vertebral body compressive fracture using quantitative computed tomography, J. Bone Joint Surg, Vol. 67-A, pp. 1206-14        [ Links ]

29.  Mc. Elhaney J.H., Alem N.M., Roberts VL.,(1970), A porous block model for cancellous bone, New York, AM. Soc. Of Mech. Engineers publication 70 WA/BHF-2; 1-9.        [ Links ]

30. Metha B:V, Rajani S., (1995), Finite element analysis of the human tibia, Computer Simulations in Biomedicine, Editors H. Power., R.T. Hart., Computational Mechanics Publications, 309-316.        [ Links ]

31. Moody D., Lazonoff S, (1999), SURFdriver, User Manual        [ Links ]

32. Müller-Karger, CM., Gonzalez C. Aliabadi, M, Cerrolaza M., (1999) Boundary Element Análisis of Tibial Plateau, International Conference on Boundary Elements, Londres 6 al 10 Julio 1999.        [ Links ]

33. Müller-Karger C.M., Bröker H., Rank E, Cerrolaza M, (2000), 3D Geometric modeling and analysis of bone using the p- version finite element method. In proceedings of European Congress on Computational Applied Sciences and Engineering, Barcelona 2000.        [ Links ]

34. Müller-Karger, C.M., Rodríguez R., Marques C., Cerrolaza M. (2000) Análisis tridimensional de la tibia para el estudio mec'anico de la rodilla": V Congreso Internacional de Métodos Numéricos en Ingeniería y Ciencias Aplicadas (CIMENICS´2000), Puerto la Cruz, Venezuela entre el 20-24 de abril de 2000.        [ Links ]

35. Neil J.L., Demos T.C., Stone J.L, Hayes W.C., (1983), Tensile and Compressive Properties of Vertebral Trabecular Bone, Trans Orthop Res. Soc. 8: 344.        [ Links ]

36. Nigg B., Herzog W.,(1994), Biomechanics of the Musculo-Skeletal System, John Wiley & Sons Ltd.        [ Links ]

37. Pettermann H.E., Reiter T.J., Rammerstorfer F.G., (1997), Computational Simulation of Internal Bone Remodeling, Archives of Computational Methods in Engineering, 4(4):295-323.        [ Links ]

38. Polikeit A, Orr T.E., Nolte L-P.(2000), Factors influencing the stresses in the lumbar vertebra after insertion of interbody cages: finite element analyses. 12th Conference of the European Society of Biomechanics, Dublin, 2000, p.168.        [ Links ]

39. Reilly D.T., Burstein A. H., (1975), The elastic and ultimate properties of compact bone tissue . J. Biomechanics, 8: 393-405.        [ Links ]

40. Rice J.C., Cowin S.C., Bowman J. A., (1988), On the dependence of the elasticity and strength of cancellous bone on apparent density , J. Biomech., 21 (2): 155-168.        [ Links ]

41. Rohl L. Larsen E., Linde F., Odgarrd A and Jorgensen J (1991) Tensile and compressive properties of cancellous bone , J. Biomechanics, 24, 1143-1149.        [ Links ]

42. Rüegsegger P., Koller B., Müller R., (1996), A microtomographic system for the nondestructive evaluation of bone architecture. Calcification Tissue International 58, 24-29.        [ Links ]

43. Szabó B., Babuška I, (1991), Finite Element Method. John Wiley and Sons, Inc. New York.        [ Links ]

44. Steiz P, Ruegsegger P, (1983), Fast Contour detection algorithm for high precision quantitative CT, IEEE Trans Md Img; MI2(3): 136-141.        [ Links ]

45. Townsend P.R., Raux P., Rose R. M., Miegel R. E., Radin E. L., (1975), The distribution and anisotropy of the Stiffness of Cancellous Bone in the Human Patella, J. Biomechanics, 8 :363-367.        [ Links ]

46. van Rietbergen B., Müller R., Ulrich D., Rüegsegger P., Huiskes R., (1999), Tissue stresses and strain in trabeculae of a canine proximal femur can be quantified from computer reconstructions, J. of Biomechanics, 32:165-173.        [ Links ]

47. van Rietbergen, F. Eckstein, B. Koller, R. Huiskes, F.P.T. Baaijens, P. Rüegsegger (2000), Trabecular bone tissue strains in the healthy and osteoporotic human femur, 12th Conference of the European Society of Biomechanics, Dublin, 2000, p.19.        [ Links ]

48. Viceconti M, Zannoni C, Baruffaldi F., Pierotti L, Toni A., Cappello A, (1998) CT-scan data acquisition to generate biomechanical model of bone structure , in Computer Methods in Biomechanics & Biomedical Engineering - 2. (J. Middleton, M.L. Jones, G. N. Pande, eds). Gordon & Bread Sci. Pub , The Netherlands, 279:287, 1998        [ Links ]

49. Weaver and Chalmers (1966), Cancellous Bone: Its Strength and changes with aging ad an evaluation of some methods for measurements its mineral content. I. Age changes in cancellous bone,. J. Bone and Joint Surg. 48-A: 289-299,        [ Links ]

50. Wolff (1872) Beiträge zür Lehre von der Heilung der Frakturen. Arch. Klin. Chri., 14:389-453.        [ Links ]

51. Yoon, H.S., y Katz J.L. (1976) Ultrasonic wave propagation in human cortical bone: Measurement of elastic properties ans micro hardness. J. Biomechanics, 9:459-465        [ Links ]

52. Zannoni C, Mantovani R, Viceconti M. (1998) ;Material properties assignment to finite element models of bone structures: a new method. Med Eng Phys 20(10):735-40.        [ Links ]

53. Zienkiewicz O. C., Morgan K, (1983) Finite Elements and Approximation, John Wiley and Sons, INC.        [ Links ]

54. Zienkiewicz O. C., Taylor R.L. (1994) El método de los elementos finitos, 4ta Edición, Mc Graw Gill/Interamenricana de España Vol 1.        [ Links ]

55. Zysset P.K., Guo X.E., Hoffer C.E., Moore K.E., Goldstein S.A., (1999), Elastic modulus and hardness of cortical and trabecular bone lamellae measured by nanoindentation in the human femur, J. of Biomechanics, 32:1005-1012.        [ Links ]

Notas pie de pagina.
1 Digital Imaging and Communications in Medicine3
2 Ángulo de aberturadel haz de rayos 

3 Initial Graphics Exchange Specification. Formato universal que puede guardar información sobre contornos y superficies.
4 Hueso seco tratado en baño de alcohol y éter.
5 Hueso sumergido en baño de agua al vacío por 24 h.
6 Hueso tratado en un horno a 600ºC por 24 horas