Estimación de propiedades elásticas y resistentes de tejido óseo cortical a partir de ensayos experimentales y modelos de elementos finitos

i UNIVERSITAT POLITÈCNICA DE VALÈNCIA Departamento de Ingeniería Mecánica y de Materiales Trabajo Fin de Máster en Ingeniería Mecánica y de Materiales ESTIMACIÓN DE PROPIEDADES ELÁSTICAS Y RESISTENTES DE TEJIDO ÓSEO CORTICAL A PARTIR DE ENSAYOS EXPERIMENTALES Y MODELOS DE ELEMENTOS FINITOS Presentada por: D. Carlos Tudela Desantes Dirigida por: Dr. D. Eugenio Giner Maravilla Valencia, septiembre de 2016 V A L È N C I A V N I V E R S I T A T P O L I T È C N I C A E X T E C H N I C A P R O G R E S S I O ix x AGRADECIMIENTOS Gracias a mi tutor Eugenio por su infinita paciencia, y por confiar en mi. A Ricardo por todo lo que me ha enseñado y todo el tiempo que me ha dedicado. A Miriam porque me invita a ser mejor cada día con su ejemplo. A mi familia por todo el apoyo y ayuda que he recibido siempre en el ámbito académico y personal. xi xii RESUMEN En la presente tesina final de máster, se presentan los resultados de ensayos a flexión en 3 puntos en muestras de tejido óseo cortical de pequeñas dimensiones, junto con su correlación modelada con elementos finitos. Algunas de las propiedades elásticas y resistentes del material han sido inferidas mediante ensayos experimentales y simulaciones. Es más, se ha simulado la iniciación y crecimiento de microgrietas en las probetas ensayadas usando un modelo de daño basado en la deformación máxima principal. Las simulaciones muestran una buena correlación con los resultados experimentales. Los ensayos experimentales se han llevado a cabo en hueso cortical de tibia de oveja, usando microprobetas preparadas de alrededor de 20 mm de largo y 1 mm de grosor con la finalidad de ensayarlas en un test a flexión en 3 puntos. Las microprobetas se han teñido para su observación con microscopio óptico para distinguir su morfología microestructural. Experimentalmente, se ha observado que las microgrietas evolucionan a lo largo de las líneas cementantes y através del material intersticial pero sin cruzar las osteonas, tal como el modelo de elementos finitos ha predicho. La geometría de las muestras se ha considerado usando elementos finitos. Se han calibrado las propiedas elásticas y resistentes del modelo numérico de manera que su respuesta se corresponda con los resultados experimentales de fuerza-desplazamiento. De la misma manera, mediante análisis inverso, es posible estimar propiedades pobremente documentadas en la literatura. Palabras clave: Hueso cortical, microgrietas, método de los elementos finitos. ix ABSTRACT In this final master thesis, the results from 3-point bending tests on cortical tissue samples of small dimensions are presented, together with its correlation modeled by finite elements. Some of the most relevant elastic and toughness material properties have been inferred from experimental tests and simulations. Further, the initiation and growth of microcracks in the tested samples is simulated by using a damage model based on the maximum principal strain. Simulations show a good correlation with the experimental results. Experimental tests have been carried out on sheep tibia cortical bone, using prepared microsamples of about 20 mm long and 1 mm thick in order to test them on 3-point bending test conditions. Microsamples have been dyed for its observation under optic microscope in order to distinguish the micro scale morphology. Experimentally, it is observed that microcracks evolve along the cement lines, through the interstitial material but without crossing osteons, as the finite elements models predicted. Sample geometries have been considered using finite elements to calibrate the elastic and toughness properties of the FE model so that its response corresponds to the force-displacement experimental results. This way, by inverse analysis, it is possible to estimate propierties poorly documented in the literature. Keywords: Cortical bone, microcracks, finite element method. x RESUM Aquesta tesina final de màster, els resultats d’un assaig a flexió en 3 punts en mostres òssies corticals de xicotetes dimensions son presentats, junt amb la seua correlació modelada amb elements finits. Algunes de les propietats elàstiques i resistents més importants han sigut inferides a partir dels assajos experimentals i les simulacions. A més, se simulen la iniciació i creixement de les microesquerdes en les mostres emprant un model de dany basat en la deformació principal máxima. Les simulacions evidencien una bona correlaciò amb els resultats experimentals. Els assajos experimentals s’han dut a terme en os cortical de tíbia d’ovella, usant micromostres de, aproximadament, 20 mm de largària i 1 mm de grosor amb la finalitat d’ensajar-los en condicions de flexiò a 3 punts. Les micromostres s’han tenyit per a la seua observació amb microscopi òptic per a distingir la morfología microestructural. Experimentalment, s’observa que que les microesquerdes evolucionen al llarg de les línies cementants i a través del material intersticial però sense creuar les osteones, tal como el model d’elements finits ha predit. La geometria de les probetes s’ha considerat usant elements finits per a calibrar les propietats elásticas i resistents del model numèric de manera que la seua resposta corresponga amb els resultats experimentals de forçadesplaçament. De la mateixa manera, per anàlisi invers, es posible estimar propietats poc documentades en la literatura. Keywords: os cortical, microesquerdes, mètode dels elements finits. xi ix Contenido 1 INTRODUCCIÓN .......................................................................................... 1 1.1 MOTIVACIÓN ............................................................................................. 1 1.2 MICROESTRUCTURA ÓSEA .......................................................................... 1 1.3 MECÁNICA DE LA FRACTURA Y MÉTODO DE LOS ELEMENTOS FINITOS ........ 5 1.4 OBJETIVOS ................................................................................................. 7 1.5 ORGANIZACIÓN DE LA TESIS DE MÁSTER .................................................. 8 2 DESCRIPCION DE LOS ENSAYOS REALIZADOS ................................ 9 2.1 OBTENCIÓN DE LAS MUESTRAS ................................................................... 9 2.2 ENSAYO MECÁNICO Y RESULTADOS ......................................................... 11 2.2.1 Observación de las microgrietas ................................................... 13 3 MODELO NUMÉRICO ............................................................................... 16 3.1 DESCRIPCIÓN DEL MODELO ...................................................................... 16 3.1.1 Geometría ...................................................................................... 16 3.1.2 Mallado y condiciones de contorno. .............................................. 21 3.2 PROPIEDADES ELÁSTICAS ......................................................................... 23 4 RESULTADOS NUMÉRICOS Y CALIBRACIÓN DE PROPIEDADES 25 4.1 RESPUESTA FUERZA DESPLAZAMIENTO .................................................... 25 4.2 CALIBRACIÓN PROPIEDADES ELÁSTICAS .................................................. 30 4.3 CONCLUSIÓN ....................................................................................... 30 5 MODELO DE DAÑO ................................................................................... 31 5.1 MÉTODO DE LA MECÁNICA DE DAÑO CONTINUO. ..................................... 31 5.2 MODELO DE DAÑO EMPLEADO ....................................................... 31 5.3 IMPLEMENTACIÓN EN ABAQUS ................................................................. 33 5.4 RESULTADOS ........................................................................................... 33 6 CONCLUSIONES Y DESARROLLOS FUTUROS .................................. 39 6.1 CONCLUSIONES ........................................................................................ 39 6.2 DESARROLLOS FUTUROS .......................................................................... 40 7 BIBLIOGRAFÍA .......................................................................................... 41 8 ANEXOS ........................................................................................................ 43 8.1 SUBRUTINA FORTRAN ABAQUS ................................................................ 43 8.2 ARCHIVO INPUT ABAQUS .......................................................................... 44 1 INTRODUCCIÓN 7 1.4 OBJETIVOS El objetivo de este trabajo es la estimación de propiedades elásticas y resistentes del tejido óseo cortical a nivel microestructural y, más concretamente, de la línea cementante que separa las osteonas (unidad básica del tejido cortical) de la matriz intersticial gracias apartir de ensayos experimentales y su posterior modelado numérico empleando el método de elementos finitos (FEM). Son muy escasas las referencias en la literatura a estas propiedades resistentes y, sin embargo, la evidencia experimental demuestra que juegan un papel en el comportamiento a fractura de las microgrietas que aparecen en el tejido óseo cortical [4]. Asimismo se pretende simular la iniciación y el crecimiento de microgrietas en algunas de las probetas ensayadas, empleando para ello un modelo de daño cuadrático basado en deformaciones principales máximas. 1 INTRODUCCIÓN 8 1.5 ORGANIZACIÓN DE LA TESIS DE MÁSTER En primer lugar, en el capítulo 2 se detalla la obtención de las muestras usadas para obtener las curvas experimentales. A continuación, se describen los ensayos experimentales y sus condiciones de contorno, siguiendo con las gráficas fuerza-desplazamiento obtenidas y su análisis. Al final del capítulo se procede a la observación de las microgrietas aparecidas en las microprobetas durante los ensayos experimentales En el capítulo 3 se explica la creación del modelo numérico, empezando por el modelado geométrico y la generación del mallado. Después se enumeran los materiales y sus propiedades elásticas a partir de la bibliografía. En el capítulo 4, se analizan los resultados fuerza-desplazamiento obtenidos en la simulación numérica y se calibran las propiedades en función de los resultados experimentales; después se comparan los resultados con la bibliografía y se arrojan unas conclusones. En el capitulo 5, se define el método de la mecánica de daño y se establece el modelo de daño basado en deformaciones. A continuación se detalla la implementación de una subrutina en el programa Abaqus con el fin de adaptar el programa a los requerimientos del método. Posteriormente se fija un valor incial de deformación crítica y se obtienen por iteración los resultados finales. Por último se analizan los resultados de deformación crítica así como la evolución de la grieta. Las conclusiones extraídas de todo el trabajo realizado se recogen en el capítulo 6. La memoria recoge una lista de referencias bibliográficas en el capitulo 7. Y termina el documento en el capitulo 8, el anexo que recoge las partes más importante de los archivos de Abaqus y Fortran utilizados a lo largo de esta Tesina. 2 DESCRIPCIÓN DE LOS ENSAYOS REALIZADOS 9 2 DESCRIPCION DE LOS ENSAYOS REALIZADOS En el presente capítulo se hará una breve descripción sobre los ensayos realizados por el Instituto de Biomecánica de Valencia (IBV), empezando por la obtención y el tratamiento de las muestras, continuando con las condiciones y parámetros del ensayo mecánico y finalizando con los resultdos de los ensayos y la observación de las microgriestas preoducidas. 2.1 OBTENCIÓN DE LAS MUESTRAS En primer lugar se realizaron cortes transversales en la diáfisis de tibia de oveja con el fin de obtener secciones de tejido cortical de 1 mm de espesor como la mostrada en la figura 2.1 izda. Seguidamente, tras eliminar los tejidos ascociados al periostio y endostio, se obtuvieron cuatro muestras en las zonas de menor curvatura, conformando las cuatro microprobetas a ensayar. En cada una de estas probetas se realizó una microentalla de aproximadamente 0.5 mm de anchura (véase figura 2.2). Las muestras se conservaron en frío y en las condiciones de humedad apropiadas. Cabe aclarar que la morfología y las características de un hueso de oveja son similares a las del hueso de un ser humano. 2 DESCRIPCIÓN DE LOS ENSAYOS REALIZADOS 10 Figura 2.1 Izda.: Sección de tibia de oveja (sin tratar y con tinción de contraste). Dcha.: Ensayo de flexión en tres puntos de la microprobeta extraída de la sección de tibia de oveja. Figura 2.2: Microprobeta del ensayo 2 teñida con tinción de contraste después del ensayo. 2 DESCRIPCIÓN DE LOS ENSAYOS REALIZADOS 11 2.2 ENSAYO MECÁNICO Y RESULTADOS Las microprobetas fueron ensayadas a flexión en tres puntos en una máquina universal electromecánica. La figura 2.1 muestra el utillaje utilizado para el apoyo y aplicación de la carga. La aplicación de la carga se llevó a cabo a una velocidad de 1 µm/seg, de manera que los ensayos se pueden considerar cuasiestáticos. Para cada muestra, se registró una fuerza aplicada en función del desplazamiento del punto de aplicación. Los resultados se recogen en la figura 2.3 y muestran una repetibilidad aceptable. Resulta evidente la gran linealidad observada hasta la aplicación de una fuerza del orden de 4-5 N. Hasta ese rango de fuerzas aplicadas, el comportamiento se puede considerar elástico lineal. A partir de ese valor de fuerza, comienza a aparecer daño generalizado en las zonas más solicitadas (en el entorno de la entalla) y tiene lugar el crecimiento de microgrietas, lo que ocasiona la pérdida de rigidez apreciable en las distintas curvas de la figura 2.3. 2 DESCRIPCIÓN DE LOS ENSAYOS REALIZADOS 12 Figura 2.3: Respuesta fuerza desplazamiento registrada para los 4 ensayos realizados. La realización de los ensayos se llevó a cabo sin alcanzar la rotura en dos mitades de las microprobetas, con el fin de observar a posteriori el camino seguido por las microgrietas. El siguiente paso consistió en la preparación de las muestras con tinción de contraste (véase figura 2.2) y la obervación de los detalles morfológicos (osteonas, canales de havers, etc.) con el fin de delimitar las osteonas a través de sus líneas cementantes así como las microgrietas existentes (Fig. 2.4) Para la tinción de contraste se siguió el procedimiento de Bain, Impeduglia y Rubin [6]. Este procedimiento implica deshidratar las muestras, por lo que sólo se aplicó después de la realización de los ensayos. El colorante de línea cementante está compuesto por azul de toluidina (1gr.), en solución de ácido fórmico al 0.1% (100 mL), con pH 2.6. La tinción con el colorante de línea cementante se realiza durante 20 segundos. 2 DESCRIPCIÓN DE LOS ENSAYOS REALIZADOS 13 Posteriormente, se deshidrata en alcohol tert-butílico (2-metil-2-propanol) durante 30 seg. En la medida de lo posible, es importante reducir la aparición de artefactos de tinción (manchas y restos) aumentando el tiempo de deshidratación en alcohol tert-butílico. Finalmente, se realiza un aclarado rápido en xileno. La figura 2.4 muestra la distribución de osteonas a través de la microscopía óptica para la microprobeta del ensayo 2. Se observa la concentración de colorante en la matriz intersticial, dejando definidas en tono más claro las osteonas y el canal de Havers. El diámetro medio de las osteonas es de alrededor de 100 µm, aunque en el corte algunas presentan formas claramente elípticas cuyo eje mayor es del orden de 150-200 µm. 2.2.1 Observación de las microgrietas Se observa la existencia de microgrietas que evolucionan principalmente a lo largo de los contornos de las osteonas (es decir, a lo largo de las líneas cementantes). El origen de esta microgrieta se encuentra en la zona de máxima concentración de tensiones (fondo de entalla). La evolución de las microgrietas se observa claramente en el corte 2 (figura 2.4). La grieta avanza rodeando las osteonas y no atravesándolas. Figura 2.4: Microprobeta del ensayo 2. Observación de osteonas al microscopio óptico (izda. 25x, dcha. 50x) Dicha tendencia observada en el corte 2 se repite en el corte 3 (figura 2.5). Aunque en este caso se trata una grieta de menor tamaño, se observa 2 DESCRIPCIÓN DE LOS ENSAYOS REALIZADOS 14 con bastante claridad como la grieta nace en el fondo de la entalla y se propaga rodeando la primera osteona de gran tamaño por su izquierda. Figura 2.5: Microprobeta del ensayo 3. Observación de osteonas al microscopio óptico (izda. 50x, dcha. 100x) En el corte 4 se repite la tendencia de forma muy clara (figura 2.6). En este caso la grieta es de mayor tamaño y en todo momento rodea las osteonas por la línea cementante. De hecho se observa como el origen de la microgrieta parece estar en la línea cementante de una osteona que intersecta con el fondo de la entalla en su parte central, lo que confirma la hipótesis de que las grietas tienden a seguir la línea cementante más que a atravesar las osteonas. 2 DESCRIPCIÓN DE LOS ENSAYOS REALIZADOS 15 Figura 2.6: Microprobeta del ensayo 4. Observación de osteonas al microscopio óptico (50x) En el corte 1 las osteonas se visualizan de manera muy deficiente de forma que no es muy útil desde el punto de vista de la observación de la grieta (figura 2.7). Sin embargo parece compartir la misma tendencia que el resto de probetas. Figura 2.6: Microprobeta del ensayo 1. Observación de osteonas al microscopio óptico (100x) 3 MODELO NUMÉRICO 16 3 MODELO NUMÉRICO En el presente capítulo se hará una descripción sobre el proceso de generación de los modelos geométricos y su posterior mallado, así como las condiciones de la simulación. A continuación se indicarán y explicarán las propiedas elásticas empleadas en los distintos materiales del modelo 3.1 DESCRIPCIÓN DEL MODELO 3.1.1 Geometría Las geometrías a modelar son las cuatro microprobetas obtenidas tal como está explicado en el capítulo 2.1. Todos los pasos que se describen a continuación se han realizado con las cuatro microprobetas obteniendo cuatro modelos geométricos. En el proceso de modelado geométrico se distinguen dos etapas: modelado de la macroestructura y modelado de la microestructura. Para modelar la macroestructura se realiza una fotografía de alta resolución de la microprobeta sobre papel milimetrado como escala (véase figura 3.1). A continuación la imagen se procesa con el software Plot Digitizer, que permite digitalizar imágenes o gráficos escaneados. En el presente trabajo, se utiliza para obtener los puntos que determinen el contorno de la microprobeta, habiendo escalado previamente la imagen con tres puntos que definen el eje x y el eje y gracias a la referencia del papel milimetrado (véase figura 3.2). Como resultado se obtiene una tabla con las coordenadas x e y de cada punto (véase tabla 3.1). 3 MODELO NUMÉRICO 23 3.2 PROPIEDADES ELÁSTICAS Se consideran tres materiales diferentes en el modelo, según se trate de la osteona, matriz intersticial o línea cementante. El modelo de elementos finitos incluye distintas propiedades elásticas para las osteonas, matriz intersticial y línea cementante [5]. Es conocido que el tejido óseo presenta un grado de anisotropía importante. Sin embargo en el plano transversal bajo estudio en este trabajo no hay grandes diferencias en cuanto a propiedades, por lo que los tres materiales se pueden suponer isótropos. Un resumen de las propiedades élasticas y resistentes se puede encontrar en [7]. Osteona. En la literatura existen diferentes propuestas de distribución de los grupos de fibrilos de colágeno dentro de la lamela, siendo una de las más aceptadas actualmente la de 5 subcapas con espesores que varían entre 0.1 y 1.8 µm (véase [8] y las referencias allí contenidas). En [8] se calcula la matriz constitutiva del tejido lamelar, agrupando las subcapas en lamela delgada de 0.8 µm y lamela gruesa de 2.4 µm. A partir de los resultados de [8], se toma como valor inicial homogeneizado e isótropo una una rigidez de la osteona 𝐸𝑜𝑠𝑡 = 11.18 GPa y como coeficiente de Poisson 𝜐𝑜𝑠𝑡 = 0.3. Matriz intersticial. Está asociada con segmentos de antiguas osteonas que han sido sustituidas en el continuo proceso de la remodelación ósea [1]. Sus propiedades elásticas corresponden a las de un tejido más mineralizado que el de las osteonas secundarias. Existe consenso en la bibliografía en que el tejido intersticial tiene una rigidez aproximadamente entre 10 y 15% mayor que la de las osteonas [3,5]. Se ha supuesto que la rigidez es un 10% superior, por lo que 𝐸𝑖𝑛𝑡𝑒𝑟𝑠 = 12.298 GPa y 𝜐𝑖𝑛𝑡𝑒𝑟𝑠 = 0.3. Línea cementante. La línea cementante se puede considerar la interfase entre la osteona y la matriz intersticial. Es muy delgada (aproximadamente 1 µm) comparada con áreas lamelares y su composición es aún un tema discutido en la bibliografía. Existen autores que sugieren que las fibras de colágeno no atraviesan las líneas cementantes y por ello representan el eslabón más débil del tejido óseo cortical. Según este planteamiento, se explica el fenómeno observado experimentalmente por el que las microgrietas tienden a seguir las líneas cementantes más que a atravesar osteonas o tejido intersticial [9], como sucede en este trabajo. En [4] se realiza un análisis de los constituyentes de la línea cementante y se sugiere que es una región de mineralización reducida. Esta composición es 3 MODELO NUMÉRICO 24 consistente con la hipótesis de que la línea cementante constituye una interfaz relativamente débil con el tejido intersticial. Una interfaz débil entre la matriz (tejido intersticial) y las fibras (osteonas) conlleva poca capacidad de transmisión de esfuerzos, lo que promueve la iniciación de grietas en la línea cementante aunque ralentiza su crecimiento en el tejido intersticial. En la misma dirección, Thurner et al. Analizan en [10] el comportamiento de la línea cementante en el tejido cortical mediante un modelo tridimensional y se afirma que para el caso de flexión las deformaciones aumentan en las interfaces pero las tensiones disminuyen. Estos autores consideran en [10] que 𝐸𝑐𝑙 = 0.088 GPa y 𝜐𝑐𝑙 = 0.3. Sin embargo, de acuerdo con M.G Ascenzi et al. [11] se ha supuesto que 𝐸𝑐𝑙 = 3.3 GPa y 𝜐𝑐𝑙 = 0.3 que son los valores inicialmente considerados en este trabajo. Otros autores suponen valores de rigidez aún mayores [7]. 4 RESULTADOS NUMÉRICOS Y CALIBRACIÓN DE PROPIEDADES 25 4 RESULTADOS NUMÉRICOS Y CALIBRACIÓN DE PROPIEDADES 4.1 RESPUESTA FUERZA DESPLAZAMIENTO El primer objetivo del análisis es verificar la respuesta fuerzadesplazamiento del punto de aplicación de la carga en el rengo lineal. Un análisis preliminar de elementos finitos con los valores de rigidez propuestos en el apartado anterior muestra que el modelo numérico se comporta de forma mucho más rígida que los resultados obtenidos experimentalmente (el modelo numérico es hasta 32 veces más rígido). Teniendo en cuenta la simplicidad de las condiciones de contorno y que las dimensiones geométricas se han modelado con precisión, se debe variar la rigidez de las propiedades elásticas. 4 RESULTADOS NUMÉRICOS Y CALIBRACIÓN DE PROPIEDADES 26 Figura 4.1: Respuesta fuerza-desplazamiento registrada en el ensayo 1 y respuesta del modelo numérico tras calibración y ajuste de propiedades. 0 1 2 3 4 5 6 7 00,1 0,2 0,3 0,4 Fuerza (N) Desplazamiento (mm) ENSAYO 1 ABAQUS 4 RESULTADOS NUMÉRICOS Y CALIBRACIÓN DE PROPIEDADES 27 Figura 4.2: Respuesta fuerza-desplazamiento registrada en el ensayo 2 y respuesta del modelo numérico tras calibración y ajuste de propiedades. 0 1 2 3 4 5 6 7 00,1 0,2 0,3 0,4 Fuerza (N) Desplazamiento (mm) ENSAYO 2 ABAQUS 4 RESULTADOS NUMÉRICOS Y CALIBRACIÓN DE PROPIEDADES 28 Figura 4.3: Respuesta fuerza-desplazamiento registrada en el ensayo 3 y respuesta del modelo numérico tras calibración y ajuste de propiedades. 0 1 2 3 4 5 6 7 00,1 0,2 0,3 0,4 Fuerza (N) Desplazamiento (mm) ENSAYO 3 ABAQUS 4 RESULTADOS NUMÉRICOS Y CALIBRACIÓN DE PROPIEDADES 29 Figura 4.4: Respuesta fuerza-desplazamiento registrada en el ensayo 4 y respuesta del modelo numérico tras calibración y ajuste de propiedades 0 1 2 3 4 5 6 7 00,1 0,2 0,3 0,4 Fuerza (N) Desplazamiento (mm) ENSAYO 4 ABAQUS 4 RESULTADOS NUMÉRICOS Y CALIBRACIÓN DE PROPIEDADES 30 4.2 CALIBRACIÓN PROPIEDADES ELÁSTICAS Tras el ajuste de propiedades, los nuevos módulos de elasticidad calibrados para cada una de las regiones de las cuatro probetas son los siguientes (tabla 4.1): 𝑬𝒐𝒔𝒕 (GPa) 𝑬𝒄𝒍 (GPa) 𝑬𝒊𝒏𝒕𝒆𝒓𝒔 (GPa Ensayo 1 0.7900 0.1041 0.8700 Ensayo 2 0.8244 0.1086 0.9068 Ensayo 3 0.8420 0.1110 0.9300 Ensayo 4 0.8384 0.1105 0.9222 Tabla 4.1: Módulos de elasticidad calibrados El ajuste de la respuesta es ahora prácticamente el esperado, como se puede observar en las figuras 4.1, 4.2, 4.3 y 4.4. Se concluye que el ensayo realizado en condiciones húmedas presenta una rigidez del tejido óseo notablemente menor que la documentada en literatura. Hay que tener en cuenta también que las osteonas presentan mayor rigidez en la dirección axial, dirección que no interviene en el ensayo realizado, ya que aquí únicamente es relevante la rigidez transversal de las osteonas. 4.3 CONCLUSIÓN En cualquier caso, los valores obtenidos son relativamente cercanos a los valores proporcionados por Nobakhti, Limbert y Thurner [10], excepto para el caso de la rigidez de la matriz intersticial. En este trabajo se habla de un rigidez asociada a las osteonas de 𝐸𝑜𝑠𝑡,𝑡𝑟𝑎𝑛𝑠𝑣𝑒 𝑟𝑠𝑎𝑙 = 0.15 GPa, 𝐸𝑜𝑠𝑡,𝑙𝑜𝑛𝑔𝑖𝑡𝑢𝑑𝑖𝑛𝑎𝑙 = 5 GPa, 𝐸𝑖𝑛𝑡𝑒𝑟𝑠𝑡𝑖𝑐𝑖𝑎𝑙 = 13.7 GPa y 𝐸𝑐𝑙 = 0.088 GPa, que para el caso que de 𝐸𝑜𝑠𝑡,𝑡𝑟𝑎𝑛𝑠𝑣𝑒𝑟𝑠𝑎𝑙 y 𝐸𝑐𝑙 son incluso menores que los aquí estimados tras la calibración, lo que demuestra que los valores obtenidos en estas dos regiones se pueden considerar admisibles. 5 MODELO DE DAÑO 31 5 MODELO DE DAÑO 5.1 MÉTODO DE LA MECÁNICA DE DAÑO CONTINUO. La mecánica de daño continuo es una teoría constitutiva que describe la pérdida progresiva de la integridad del material debido a la propagación y unión de microgrietas, microporos y efectos similares. Estos cambios en la microestructura llevan a una degradación de la rigidez del material observada en la macroescala [14]. 5.2 MODELO DE DAÑO EMPLEADO El daño de hueso cortical es modelado con la degradación de sus propiedades mecánicas cuando se alcanzan valores críticos usando una subrutina de Abaqus. Se ha propuesto que el proceso de fallo óseo está controlado por la deformación [12] y que la estructura ósea es fuertemente jerárquica [1,2], por lo que diferentes deformaciones críticas de daño han sido propuestas para cada componente. En régimen cuasi-estático, la relación isotrópica de elasticidad basada en la mecánica de daño es expresada con [13]: 𝜎𝑖𝑗 = 1− 𝐷 𝐶𝑖𝑗𝑘𝑙 𝜀𝑘𝑙 (1) Donde D es la variable de daño, 𝜎𝑖𝑗 son los componentes de la tensión, 𝜀𝑘𝑙 las deformaciones y 𝐶𝑖𝑗𝑘𝑙 son los componentes tensores de elasticidad. La subrutina de Abaqus implementada tiene una variable de estado (f), que es el criterio de daño. Ha sido elegida para depender de la deformación máxima principal (ver ecuación 2). Cuando los valores críticos de deformación son alcanzados f ≥ 1 y las variables de daño cambian de no dañado a dañado, resultando en la degradación de la rigidez del material hasta el 1% (gráfica 4.1), para representar la “smeared crack”, modelando microgrietas [14]. 𝑓=𝜀𝑝𝑚𝑎𝑥 ,𝑖 𝜀𝑖 𝑐 (2) 5 MODELO DE DAÑO 32 Figura 5.1: Degradación de la rigidez del material al 1% al alcanzar la deformación crítica Donde 𝜀𝑝𝑚𝑎𝑥 ,𝑖 es la máxima deformación principal actual y 𝜀𝑖 𝑐 la deformación crítica de daño de cada componente. El método de modelado de daño implementado es dependiente de la malla y del número de incrementos, por lo que se ha llevado a cabo un análisis de sensibilidad atendiendo a estas dos variables. En cuanto a la dependencia de los incrementos de tiempo, los incrementos deben ser suficientemente pequeños para que la solución no varíe, mientras que la dependencia de la malla implica suficiente refinamiento de malla. σ ε 𝜺𝒄 6 CONCLUSIONES Y DESARROLLOS FUTUROS 39 6 CONCLUSIONES Y DESARROLLOS FUTUROS Resumen del capítulo: En este capítulo se recogen las principales conclusiones extraídas de los capítulos anteriores y se enumeran posibles desarrollos futuros que continúan la línea de investigación en la que se enmarca la Tesina aquí presentada. 6.1 CONCLUSIONES  Se ha llevado a cabo un ensayo a flexión de 3 puntos en muestras de tejido óseo cortical, junto con su modelado en elementos finitos.  La geometría de las probetas y la morfología de las osteonas en la región de interés se han tenido en cuenta en el modelo numérico, distinguiendo entre tres componentes a nivel microscópico: osteonas, material intersticial y líneas cementantes.  La respuesta fuerza-desplazamiento registrada experimentalmente ha permitido calibrar las propiedades elásticas de los microcomponentes considerados.  Se ha demostrado que las propiedades elásticas calibradas son similares a las que arroja la documentación consultada, lo cual es útil por la dificultad para la estimación de estas propiedades experimentalmente y la vasta diversidad de valores encontrados en la literatura, para los cuales actualmente no hay un consenso.  También se ha calibrado la deformación crítica de daño para cada microcomponente (osteona, material intersticial y línea cementante) resultando alrededor de 17% para las osteonas y líneas cementantes y sobre un 7% para el material intersticial. 6 CONCLUSIONES Y DESARROLLOS FUTUROS 40  Se ha correlacionado con el método de la mecánica de daño continuo (para material isótropo, lineal y régimen cuasiestático) la iniciación de daño y su evolución con lo observado en los test experimentales.  Tanto los resultados experimentales como las simulaciones numéricas han permitido verificar la hipótesis de que la aparición de microgrietas en tejido óseo cortical ensayado a flexión en 3 puntos está condicionada por las deformaciones, debido a los diferentes valores de rígidez de cada componente de la microestructura.  Se ha evidenciado que la grieta evoluciona a lo largo de las líneas cementantes y cruzando el material intersticial, pero no las osteonas, reforzando la idea de que las líneas cementantes constituyen un elemento protector de las osteonas y la vascularización ósea. 6.2 DESARROLLOS FUTUROS  Repetir ensayos de flexión en 3 puntos en un mayor número de muestras para tener un mayor número de resultados experimentales.  Realizar más modelos numéricos y sus correspondientes calibraciones de propiedades elásticas así como deformaciones críticas de daño. De esta manera se reforzarían las conclusiones aquí extraídas.  Asociar el mapa de daño obtenido con una propiedad intrínseca del material como es el ratio de liberación de energía.  Implementar otros modelos de daño de manera que el proceso de daño no sea binario sino gradual, con diversos escalones de degradación de la rigidez.  Realizar el mismo estudio partiendo de otras tipologías de ensayos experimentales que permitan obtener las propiedades a tracción o torsión. 7 BIBLIOGRAFÍA 41 7 BIBLIOGRAFÍA [1] S. C. Cowin, Bone mechanics handbook, CRC Press, Boca Ratón, Florida, 2001. [2] D. Taylor, J.G. Hazenberg, T.C. Lee, Living with cracks: Damage and repair in human bone. Nature Mater 6, 263-268, 2007. [3] J. Rho, J. Currey. The anisotropic Young's modulus of equine secondary osteones and interstitial bone determined by nanoindentation. J Exper Biology 204, 1775-1781, 2001. [4] D.B. Burr, M.B. Schaffler, R.G. Frederickson. Composition of the cement line and its possible mechanical role as a local interface in human compact bone. J Biomechanics 21, 939-945, 1988. [5] S. Li, A. Abdel-Wahab, E. Demirci, V.V. Silberschmidt. Fracture process in cortical bone: X-FEM analysis of microstructured models. Int J Fracture 184 (2013) 43-55. [6] S.D. Bain, T.M. Impeduglia, C.T. Rubin. Cement line staining in undecalcified thin sections of cortical bone. Stain technology 65 (1990) 1-5. [7] C. Arango, E. Giner, A. Vercher, F.J. Fuenmayor. Modelado tridimensional del fallo progresivo de la microestructura del hueso cortical mediante elementos finitos. Actas del XXXI Encuentro del Grupo Español de Fractura, San Lorenzo del Escorial (2014). [8] A. Vercher, E. Giner, C. Arango, J.E. Tarancón, F.J. Fuenmayor. Homogenized stiffness matrices for mineralized collagen fibrils and lamellar bone using unit cell finite element models with periodic boundary conditions. Biomech Model Mechanobiol 13 (2014) 437-449. [9] J.G. Skedros, J.L. Holme, E.G. Vajda, R.D. Bloebaum. Cement lines of secondary osteons in human bone are not mineral-deficient: new data in a historical perspective. Anat Rec A Discov Mol Cell Evol Biol 286 (2005) 781-803. 7 BIBLIOGRAFÍA 42 [10] S. Nobakhti, G. Limbert, P.J. Thurner. Cement lines and interlamellar areas in compact bone as strain amplifiers contributors to elasticity, fracture toughness and mechanotransduction. J Mech Behav Biomed Mat 29 (2014) 235-251. [11] M.G. Ascenzi, N.P. Kawas, A. Lutz, D. Kardas, U. Nackenhorst, J.H. Keyak. Individual specific multiscale finite element simulation of cortical bone of human proximal femur. J Comput Phys 244 (2013) 298-311. [12] Nalla, R.K., Kinney, J.H., Ritchie, R.O., 2003. Mechanistic fracture criteria for the failure of human cortical bone. Nat Mater 2, 164-168. [13] Lemaitre, J., 1985. A continuous damage mechanics model for ductile fracture. J Eng Mater Technol 107, 83-89. [14] Jirsek, M., 2011. Damage and smeared crack models. In Hofstetter, G. and Meschke, G., Numerical modeling of concrete cracking, CSIM courses and lectures, SpringerWienNewYork 532, 1-49. 8 ANEXOS 43 8 ANEXOS 8.1 SUBRUTINA FORTRAN ABAQUS Subrutina USDFLD de la probeta 2, utilizada para implementar el método de daño basado en deformaciones en Abaqus. Para cambiar el caso calculado se modifican los valores de deformación crítica marcados en negrita. Programada en FORTRAN. SUBROUTINE USDFLD(FIELD,STATEV,PNEWDT,DIRECT,T,CELENT,TIME,DTIME, 1 CMNAME,ORNAME,NFIELD,NSTATV,NOEL,NPT,LAYER,KSPT,KSTEP,KINC, 2 NDI,nshr,coord,jmac,jmtyp,matlayo,laccflg) C INCLUDE 'ABA_PARAM.INC' C CHARACTER*80 CMNAME,ORNAME CHARACTER*3 FLGRAY(15) DIMENSION FIELD(NFIELD),STATEV(NSTATV),DIRECT(3,3),T(3,3),TIME(2), * coord(*),jmac(*),jmtyp(*) DIMENSION ARRAY(15),JARRAY(15) C C Empieza aquí C MATERIAL AND STRENGTH PARAMETERS PARAMETER(ECC=17.4e-2,ECI=8.45e-2,ECO=17e-2) C C C INITIALIZE FAILURE FLAGS FROM STATEV. F = STATEV(1) !Variable de estado C C GET STRESSES AND STRAINS FROM PREVIOUS INCREMENT CALL GETVRM('EP',ARRAY,JARRAY,FLGRAY,jrcd, $ jmac, jmtyp, matlayo, laccflg) EP3 = ARRAY(3) C C IF (F .LT. 1) THEN IF (EP3 .GT. 0.) THEN IF (CMNAME .EQ. 'CEMENTLINE') THEN F=(EP3/ECC) ELSE IF (CMNAME .EQ. 'INTERSTICIAL') THEN 8 ANEXOS 44 F=(EP3/ECI) ELSE IF (CMNAME .EQ. 'OSTEONA') THEN F=(EP3/ECO) ENDIF STATEV(1)=F ENDIF ENDIF C C C IF (F .LT. 1) THEN FIELD(1)=0 ! si no falla, field=FV1=0 ELSE FIELD(1)=1 ! si falla, field=FV1=1 ENDIF C RETURN END 8.2 ARCHIVO INPUT ABAQUS Archivo input de abaqus para el modelo de la probeta 2. Para cambiar el caso calculado se modifican los valores de Material, destacados en negrita. Los comandos escritos llaman a la subrutina y la subrutina devuelve los valores correspondientes. Se han incluido únicamente aquellas partes que se consideran importantes para la explicación del modelo de daño y las condiciones de la simulación. *Heading ** Job name: Sensibilidad1 Model name: Model-1 ** Generated by: Abaqus/CAE 6.13-3 *Preprint, echo=NO, model=NO, history=NO, contact=NO ** ** MATERIALS ** *Material, name=CEMENTLINE 8 ANEXOS 45 *Depvar 1, *Elastic, dependencies=1 0.000108689, 0.3, , 0. 0.000001086, 0.3, , 1. *User Defined Field *Material, name=INTERSTICIAL *Depvar 1, *Elastic, dependencies=1 0.000906793, 0.3, , 0. 0.00000906793, 0.3, , 1. *User Defined Field *Material, name=OSTEONA *Depvar 1, *Elastic, dependencies=1 0.0008243574, 0.3, , 0. 0.000008243574, 0.3, , 1. *User Defined Field ** ---------------------------------------------------------------- ** ** STEP: Step-1 ** *Step, name=Step-1, nlgeom=NO, inc=10000 *Static, direct 0.001, 1 ** ** BOUNDARY CONDITIONS ** 8 ANEXOS 46 ** Name: BC-1 Type: Displacement/Rotation *Boundary Apoyo_1, 1, 1 Apoyo_1, 2, 2 ** Name: BC-2 Type: Displacement/Rotation *Boundary Apoyo_2, 2, 2 ** ** LOADS ** ** Name: Load-1 Type: Concentrated force *Cload Load_Point, 2, -5. ** *EL PRINT, POSITION=INTEGRATION POINTS, FREQUENCY=10 EP3,SDV,FV ** OUTPUT REQUESTS ** *Restart, write, frequency=0 ** ** FIELD OUTPUT: F-Output-1 ** *Output, field *Node Output RF, U *Element Output SDV,FV S,E ** ** HISTORY OUTPUT: H-Output-1 8 ANEXOS 47 ** *Output, history *Node Output, nset=Load_Point CF2, U2 *End Step