Estudio a Nivel Mesomecánico De La Fractura Del Hueso Trabecular

ESTUDIO A NIVEL MESOMECÁNICO DE LA FRACTURA DEL HUESO TRABECULAR. C.M. López*, O. Roa**, I. Carol* *

ETSECCPB, Universitat Politècnica de Catalunya, 08034 Barcelona ** FONACIT-UNET, Venezuela

Resumen. En este artículo se presenta un estudio numérico mesomecánico del comportamiento del hueso trabecular en el cual se emplean mallas 2D generadas aleatoriamente mediante la teoría de Voronoï–Delaunay. El modelo, originalmente desarrollado y verificado para muestras de hormigón, incorpora el comportamiento de fractura utilizando "elementos junta" de espesor cero. La estructura interna del material se representa explícitamente y las interfases se insertan a priori entre los elementos estándar del medio continuo, a lo largo de las principales potenciales trayectorias de fisuración. Para la obtención de las mallas se ha desarrollado un generador automático que permite representar mesoestructuras óseas anisótropas (también isótropas, como caso particular) con determinados parámetros geométricos como la porosidad y el tensor de fábrica. En primer término se analiza el comportamiento mecánico de geometrías isótropas bajo acciones de tracción y compresión uniaxial. Los resultados numéricos obtenidos muestran modos de rotura y curvas promedio tensión-deformación realistas cuando se considera no linealidad geométrica que permite capturar los efectos de inestabilidad y pandeo de las trabéculas individuales. Posteriormente, se consideran geometrías anisótropas. Se han realizado ensayos numéricos de tracción uniaxial en muestras cortadas con distinta orientación de los ejes principales del tensor de fábrica, para diferentes valores de la porosidad. Se presentan resultados de la variación del módulo de elasticidad y de la tensión máxima en función del ángulo de giro, así como las tendencias de fisuración y deformación obtenidas. Abstract. In this article we present a mesomechanic study of trabecular bone behavior employing 2D random meshes generated by Voronoï–Delaunay theory. The model, firstly developed and veryfied for concrete specimens, incorporates crack behavior using zero thickness “joint elements”. The internal material structure is explicitly represented, and interfaces are inserted between the predetermined standard elements of the continuum media, along the main potential crack paths. An automatic generator has been developed to obtain the meshes, which allows to represent anisotropic bone mesostructures (also isotropic, as a particular case) with some geometric parameters like porosity and fabric tensor. Firstly, the mechanical behavior of isotropic geometries under uniaxial tension and compression is analysed. Numerical results obtained show crack modes and realistic mean stress-strain curves when geometrical non-linearity that allows to capture instability effects and bending of individual trabeculae. Afterwards, anisotropic geometries are considered. Numerical tests were carried out considering uniaxial tension in specimens cut with different fabric tensor direction, for different values of porosity. Results obtained on variation of elasticity module and maximum tension with regard to rotation angle, and crack and deformation paths are presented. 1. INTRODUCCIÓN Existen dos tipos básicos de tejido óseo, el hueso cortical (compacto) y el hueso poroso o trabecular. Independientemente de la anatomía microscópica, todos los segmentos de esqueleto se componen de una zona externa de hueso compacto (cortical) y una zona interna de hueso poroso. El hueso cortical es homogéneo, altamente organizado y orientado regularmente y muy resistente. El hueso poroso, en cambio, presenta una estructura poco organizada y orientada aleatoriamente, denominada trabecular (del latín “trabeculae” = vigas pequeñas). En 2D al entramado óseo se lo denomina estructura de colmena cuando presenta una orientación preferencial acentuada. En 3D estructura en forma de celdas tridimensionales, que al conectarse presentan un aspecto esponjoso. La diferencia principal entre ambos

tejidos óseos viene dada por la porosidad que en el hueso cortical varía entre el 5% y el 30%, mientras que en el hueso trabecular varía entre el 30 y el 90% [1]. Los ensayos experimentales en huesos son técnicamente complicados. A diferencia de otros materiales de ingeniería, el hueso trabecular no puede ser confeccionado fácilmente en configuraciones arbitrarias, las dimensiones de las muestras son pequeñas y, adicionalmente se presentan dificultades para la preservación de muestras. Por estos motivos, en los últimos años se ha extendido el empleo de métodos numéricos sobre todo basados en los elementos finitos. En este trabajo el estudio se centra en el análisis a nivel mesoescala del hueso trabecular. En primer término se describen sintéticamente los aspectos más importantes de la modelización desarrollada. Luego, se presentan resultados numéricos correspondientes a ensayos de

tracción y compresión uniaxial en mallas isótropas. Finalmente, se muestran y analizan los resultados numéricos obtenidos sobre muestras anisótropas sometidas a cargas de tracción con distinta orientación respecto de los ejes principales del tensor de fábrica y para diferentes valores de la porosidad. 2. MODELIZACIÓN MESOESTRUCTURAL El comportamiento de materiales heterogéneos como el hormigón o el hueso poroso está claramente influenciado por la geometría y las propiedades mecánicas de los componentes de la estructura interna del material. Por ello, una manera más adecuada de profundizar el estudio de estos materiales consiste en representar explícitamente la microestructura de primer nivel (a veces también llamada mesoestructura), introducir leyes de comportamiento de sus componentes básicos, e intentar reproducir el comportamiento “macro” del material como resultado del análisis mecánico a nivel “meso”. Las figuras 1a y 1b muestran las imágenes de una probeta de hormigón de 15 cm. de diámetro y una muestra de hueso trabecular de 1x1x1 cm.3, tomada de [2]. En la probeta de hormigón se puede identificar fácilmente una mesoestructura constituida por partículas de áridos de un tamaño mayor, embebidas en una matriz formada por el mortero y los áridos de menor tamaño. El hueso poroso presenta una estructura entramada de barras y placas de tejido calcificado, orientadas aleatoriamente (arquitectura trabecular). En la muestra de la figura 1b se puede distinguir una estructura porosa no homogénea. Se puede apreciar en estas figuras que salvo por la diferencia de escala de longitud (cm. en un caso y mm en el otro), a este nivel de observación las geometrías de ambos materiales no es muy diferente si se reemplazan los áridos por huecos. Teniendo en cuenta la analogía entre ambas estructuras internas, en este trabajo se modeliza una muestra de hueso trabecular tomando como punto de partida un modelo mesomecánico desarrollado previamente para hormigón [3-4]. Para ello, una muestra 2D se representa mediante una matriz correspondiente a la

fase sólida de material óseo que rodea a una fase de vacíos (figura 1c). La matriz de material óseo se discretiza mediante elementos finitos triangulares a los cuales se les asigna un comportamiento constitutivo elástico lineal. A fin de posibilitar la no linealidad de comportamiento, mecanismos de rotura y fisuración, etc., la malla de elementos finitos incluye elementos junta sin espesor, distribuidos entre algunos contactos de la matriz de modo de representar las principales direcciones potenciales de fisuración (figuras 1d y 1e). Para el comportamiento de las interfases se utiliza una ley constitutiva no lineal con reblandecimiento por trabajo de fractura, formulada en términos de las tensiones en su plano medio y los desplazamientos relativos entre ambos lados de la junta [3, 5, 6]. Para la generación de la geometría se utiliza un procedimiento aleatorio basado en la teoría de Voronoï/Delaunay. Se ha desarrollado un generador automático de mallas que permite representar mesoestructuras óseas isótropas y anisótropas con determinados parámetros geométricos, como la porosidad (establece la proporción entre huecos y trabécula) y el tensor de fábrica (tensor de segundo orden que caracteriza la forma y orientación de los huecos). Para representar la anisotropía del hueso trabecular se emplean diversos parámetros tensoriales, como la longitud media de intercepción (mean intercept length, MIL). La MIL representa un promedio de la distancia intertrabecular medida a lo largo de una línea y es función del ángulo de inclinación de la línea a lo largo de la cual se hace la medición (figura 2a). Whitehouse [7] ha observado que al representar un diagrama polar con los valores de la MIL, el diagrama en 2D se aproxima a una elipse (lo que es equivalente a un tensor positivo de segundo orden). La relación entre las longitudes principales de la elipse es representativa de la anisotropía de la muestra. En las figuras 2b y 2c se muestran la representación gráfica de la longitud de intercepción media y del tensor de fábrica resultante correspondientes a muestras anisótropa e isótropa, respectivamente.

Figura 1. a) Muestra de hormigón de 15 cm de diámetro, b) muestra de hueso trabecular de 1 cm de lado, c) ejemplo de la discretización de la matriz, d) disposición de las junta, e) detalle de los elementos junta.

Figura 2. Ejemplo del cálculo del tensor de fábrica: a) obtención de la longitud de intercepción media, b) diagrama del tensor de fábrica para malla anisótropa, y c) diagrama del tensor de fábrica para malla isótropa. 3. ANÁLISIS ISÓTROPAS

DE

MESOESTRUCTURAS

Los estudios iniciales en mallas isótropas se enfocaron en sus propiedades elásticas y el error de discretización [8]. En este apartado se analiza el comportamiento no lineal y en fractura bajo acciones de tracción y compresión uniaxial. Los resultados numéricos se comparan con curvas experimentales disponibles. En cada ensayo se imponen los desplazamientos verticales uniformes en todos los nodos de la cara superior de la malla (los desplazamientos verticales de la cara inferior son fijos), mientras que los desplazamientos transversales se dejan libres excepto en un punto fijo. La deformación a lo largo de la dirección de aplicación de la carga se toma como el desplazamiento impuesto dividido por la longitud de la muestra. La tensión nominal se obtiene como la suma de las coacciones en la cara de la probeta, dividida por su ancho. Las magnitudes de los parámetros constitutivos se han obtenido a partir de diferentes ensayos numéricos realizados con el fin de ajustar los valores que permitan encontrar la mejor aproximación global a la respuesta experimental. Los resultados de los ensayos experimentales a tracción y compresión fueron tomados de [9] y [2], respectivamente (figuras 3 y 4). La figura 3 resume los resultados obtenidos en tracción. En el gráfico de la izquierda se incluyen las curvas promedio tensión-deformación experimentales y las obtenidas del cálculo considerando y sin considerar los efectos de la no linealidad geométrica ("grandes deformaciones"). En el gráfico del centro se muestra la malla deformada y en el de la derecha, la malla con la representación de la energía disipada en cada junta en el proceso de fractura (el espesor de las líneas es representativo del trabajo disipado durante el proceso de fractura en cada punto de la junta). Se observa que en general los cálculos reproducen bien la curva experimental. En el caso de considerar el efecto geométrico no lineal en el análisis, la curva decrece

hasta el valor de tensión cero de manera similar a lo que ocurre en la curva experimental, mientras que sin considerar estos efectos la evolución de la curva tensión-deformación tiende a una meseta horizontal con una fricción residual ficticia debida a que el cálculo ignora la creciente inclinación de los puentes de material entre los dos lados que se separan. Se puede apreciar que el modo de falla es aproximadamente horizontal, tal como era de esperar en correspondencia con la tendencia de fisuración comúnmente observada en tracción. Los resultados en compresión uniaxial se recogen en la figura 4. En el gráfico de la izquierda se aprecia que si no se consideran los efectos de la no linealidad geométrica la curva tensión deformación crece ilimitadamente sin llegar a un pico en el rango de deformaciones analizado. Esto sucede a pesar de haber incluido las mismas juntas susceptibles de fisuración que en el análisis en tracción. La curva obtenida con efectos geométricos no lineales parece mucho más realista, con una clara tendencia hacia una meseta horizontal. En los gráficos del centro y derecha se presenta un estado avanzado del proceso de carga (punto B de la figura. 4a) en términos de la energía de fractura gastada en las juntas y la malla deformada, respectivamente. Una observación detallada de estas figuras permite detectar una limitación que aparece en los cálculos realizados. En las zonas señaladas con círculos en la figura 4c, se presentan solapamientos de la malla en correspondencia con las juntas que desarrollan los mayores desplazamientos relativos de corte. En estas juntas la energía disipada en modo mixto es muy superior a la energía gastada en las demás juntas, fundamentalmente en modo I (figura 4b). Estos solapes podrían ser físicamente inaceptables por lo que para seguir avanzando en la historia de carga, sería necesario tener en cuenta en el análisis el problema de los nuevos contactos que se forman entre elementos inicialmente no relacionados, alrededor del perímetro de cada vacío.

Figura 3. a) Curva tensión-deformación nominal en tracción uniaxial, b) malla deformada y, c) estado de fisuración en término del trabajo disipado, para el estado final del ensayo con efectos geométricos no lineales.

Figura 4. a) Curva tensión-deformación nominal en compresión uniaxial, b) estado de fisuración en términos de la energía disipada en las juntas, y c) malla deformada, en el punto B. 4. ANÁLISIS ANISÓTROPAS

DE

MESOESTRUCTURAS

Se presentan y analizan los resultados obtenidos en muestras anisótropas frente a cargas de tracción con distinta orientación respecto de los ejes de ortotropía. Tanto el módulo de elasticidad como la resistencia dependen de la dirección de la carga aplicada para la misma morfología trabecular [10], y varían considerablemente en función de los valores de los principales parámetros geométricos como la porosidad y el tensor de fábrica [11]. Se ha realizado una serie de ensayos numéricos de tracción uniaxial en ambas direcciones de muestras cuadradas cortadas con distinta orientación de los ejes principales del tensor de fábrica (figura 5). Se han analizado seis muestras cortadas cada 15° de giro con respecto al eje horizontal (dirección principal con menor valor de tensor de fábrica) para cinco diferentes valores de la porosidad (la relación entre los valores principales del tensor de fábrica H1/H2 se mantiene constante en una relación de aproximadamente 0.7 ± 0.01 al variar la porosidad). En las figuras 6 y 7 se presentan los diagramas polares del módulo de elasticidad y de la tensión máxima, respectivamente, en función del ángulo de giro. En ambos gráficos los resultados del primer cuadrante corresponden a los ensayos en la dirección transversal

(horizontal) de la muestra para los distintos ángulos de giro analizados, mientras que los del segundo cuadrante corresponden a la carga en dirección longitudinal (vertical). En la figura 6 se puede apreciar que el módulo de elasticidad aumenta para igual porosidad desde la dirección horizontal hasta la vertical, en acuerdo con los valores principales del tensor de fábrica. Asimismo, se observa que a medida que disminuye la porosidad, se incrementa el módulo en cada orientación, si bien la relación entre los módulos a 90° y 0° decrece en términos relativos (de más del doble para porosidad de 0,74 hasta aproximadamente 1.4 para 0.34).

Figura 5. Generación de muestras giradas con respecto a las direcciones principales.

obtenido se asemeja al que cabe esperar de un medio continuo, con desplazamientos y efecto Poisson en las direcciones principales del tensor de fábrica, y con desplazamientos de corte perpendiculares a la dirección de carga para otras orientaciones, como se aprecia en la malla girada 45°. Para la malla de mayor porosidad puede apreciarse una localización más acentuada. En las direcciones principales la deformación se produce a lo largo de la muestra, mientras que para otros ángulos de giro (véase la representación de 45°) se observa una zona localizada en la que las trabéculas giran y tienden a orientarse según la dirección de carga, mientras que la otra parte de la muestra prácticamente permanece sin deformarse. 5. CONSIDERACIONES FINALES Figura 6. Diagrama polar del módulo elástico en función del ángulo de giro, para distintos valores de la porosidad. En el diagrama de tensiones máximas (figura 7) puede apreciarse una característica diferente: la resistencia en función del ángulo presenta una discontinuidad; mientras que se mantiene constante en un cierto rango de ángulos cercanos a la dirección principal de valor más bajo (trayectoria aproximadamente circular), luego, en el otro rango cercano al mayor valor sigue una tendencia distinta tipo “elipse”. Este tipo de comportamiento ya ha sido señalado en la literatura para otros materiales, como suelos o algunos metales.

Un modelo mesomecánico inicialmente desarrollado y verificado ampliamente para hormigón, se ha adaptado y aplicado en este trabajo para hueso trabecular. Los trabajos presentados forman parte de la tesis doctoral de Olga Roa, actualmente en su fase final [12]. En el análisis no lineal de mallas isótropas se han obtenido buenos resultados con modos de rotura y curvas promedio tensión-deformación realistas. En estos estudios se ha visto la necesidad de considerar la no linealidad geométrica para poder capturar los efectos de inestabilidad y pandeo de las trabéculas individuales. En el ámbito anisótropo, se han realizado ensayos numéricos de tracción uniaxial sobre distintas mallas, obtenidas a partir de una misma muestra cortada con distintas orientaciones de los ejes principales del tensor de fábrica y para diferentes valores de la porosidad. Se han presentado interesantes resultados de las relaciones del módulo de elasticidad y de la tensión máxima en función del ángulo de giro, así como de las tendencias de fisuración y deformación. En hueso poroso los ensayos experimentales son técnicamente complicados debido a las pequeñas dimensiones de las trabéculas y a la necesidad de preservación de las muestras. La muestra no puede ser confeccionada fácilmente, y adicionalmente, debido a la heterogeneidad, resulta imposible obtener distintas muestras (para ensayar en las distintas direcciones) de la misma porción de hueso. Por lo tanto, la simulación desarrollada puede resultar de una gran importancia en diversas aplicaciones biomédicas. 6. AGRADECIMIENTOS

Figura 7. Diagrama polar de la tensión máxima en función del ángulo de giro, para distintos valores de la porosidad. En las figuras 8 y 9 se presentan las mallas deformadas con porosidades de 0.44 y 0.74, respectivamente, para tres distintas orientaciones de la carga de tracción. Para la malla de menor porosidad el comportamiento

Los autores agradecen al MCYT (Madrid) la ayuda recibida a través de los proyectos de investigación MAT2000-1007 y MAT2003-02481. El segundo autor agradece la beca doctoral financiada por FONACIT y UNET (Venezuela). Finalmente, se agradece a la ETSECCPB-UPC la bolsa de viaje concedida para asistir al XXI Encuentro del Grupo Español de Fractura.

Figura 8. Mallas deformadas correspondientes a una porosidad de 0.44 para tres ángulos de giro de 0° (izq.), 45° (medio) y 90° (der.).

Figura 9. Mallas deformadas correspondientes a una porosidad de 0.74 para tres ángulos de giro de 0° (izq.), 45° (medio) y 90° (der.). 7. REFERENCIAS [1] Newton, C. and D. Nunamaker, D., “Biomechanics pertinent to fracture etiology, reduction and fixation”, Textbook of small animal orthopaedics, cap. 12, J.B. Lippincott Company (1985). [2] Pini, M., “Modellazione ed analisi di un tessuto osseo spongioso”. Tesi di Laurea, Politecnico di Milano, Italia. , (1998). [3] López, C. M., “Análisis microestructural de la fractura del hormigón utilizando elementos finitos tipo junta. Aplicación a diferentes hormigones”. Tesis Doctoral. ETSECCPB, Universitat Politècnica de Catalunya, Barcelona, España (1999). [4] López, C., Carol, I., Aguado, A., “Modelización microestructural de la fractura del hormigón. I: tracción y compresión uniaxial y ensayo brasileño", XIX Encuentro del Grupo Español de Fractura, Gerona, 13-15 de marzo (2002). [5] Carol, I., Prat, P. C. y López, C. M., “A normal/shear cracking model. Application to discrete crack analysis”, J. of Engineering Mechanics, 123, No 8 (1997). [6] Carol, I., López C. M. y Roa, O., “Micromechanical analysis of quasi-brittle materials using fracture-based interface elements”,

Int. J. Numer. Meth. in Engrg. Vol 52, 193-215 (2001). [7] Whitehouse, J., “The quantitative morphology of anisotropic trabecular bone”, J. Microscopy. Vol 101, 153-168 (1974). [8] Roa, O., Lopez, C.M., Maranzana, F., Carol, I. and Contro, R., “On the use of generated geometries for the numerical analysis of trabecular bone specimens”. Middleton, J., 2002,Computer Methods in Biomechanics and Biomedical Engineering (Proceedings of the Int. Conf. in Rome, Oct. 2001),WSAtkins Science and Technology. [9] Keaveny, T., Borchers, R., and Gibson, L, “Theoretical analysis of the experimental artefact in trabecular bone compressive modulus”, J. Biomechanics. Vol 26 (4/5), 599-607 (1993). [10] Yuehuei, H.A., R. A. Draughn, R.A., “Mechanical testing of bone and the bone-implant interface”. CRC Press (2000). [11] Cowin, S.C. “Bone Mechanics.” CRC Press, (2001). [12] Roa, O., “Análisis microestructural del comportamiento mecánico del hueso trabecular”. Tesis Doctoral en curso. ETSECCPB, Universitat Politècnica de Catalunya, Barcelona, España. (2004).