MODEL OF FINITE ELEMENTS FOR LINEAL ANALYSIS OF A POROELASTIC MATERIAL TRANSVERSE BY ISOTROPIC
DANIEL H. CORTÉS C.
Escuela de Ingeniería Mecánica. Universidad del Valle AA 25360 Cali
JOSÉ JAIME GARCÍA A.
Escuela de Ingeniería Civil y Geomática. Universidad del Valle. AA 25360 Cali
GILBERTO GÓMEZ
Escuela de Ingeniería Civil y Geomática. Universidad del Valle. AA 25360 Cali
Recibido para revisar 10 de Febrero de 2004, aceptado 23 de Junio de 2004, versión final 20 de Agosto de 2004
RESUMEN Se presentan las ecuaciones de la poroelasticidad lineal y su solución mediante el método de los elementos finitos considerando simetría axial e isotropía transversal para la fase sólida. Este modelo puede ser usado para describir el comportamiento mecánico de materiales poroelásticos como el cartílago articular y algunos suelos. Los resultados del programa, montado en la plataforma de dominio público Scilab, mostraron una buena correlación con las soluciones analíticas para el caso de compresión no confinada.
]]> PALABRAS CLAVES: Poroelasticidad, materiales bifásicos, tejidos biológicos, anisotropía.ABSTRACT: The equations of linear poroelasticity and their solution by finite elments are presented considering axial symmetry and a transversely isotropic model for the solid phase. This model can be used to describe the mechanical behavior of poroelastic materials like articular cartilage and some type of soils. The results of the computer program, implemented in the public domain platform Scilab, showed good correlation with the analytical solutions for the case of non confined compression.
KEY WORDS: Poroelasticity, biphasic materials, biologic weaves, anisotropy
1. INTRODUCCIÓNLa poroelasticidad estudia el comportamiento de los materiales elásticos cuyos poros están ocupados por un fluido. A raíz del estudio de problemas de mecánica de suelos la poroelasticidad tuvo un desarrollo importante desde inicios del siglo pasado, especialmente con los trabajos emprendidos en la Universidad de Viena por Terzaghi y Fillunger, descritos detalladamente en un artículo de De Boer [1]. Cabe anotar como hecho anecdótico que Terzaghi y Fillunger fueron antagonistas académicos y personales, disputa que terminó con el infortunado suicidio de Fillunger y su señora.
Se puede considerar que la poroelasticidad se incorporó a la biomecánica a raíz del artículo de Mow et al [2] quienes tomaron elementos teóricos ampliamente utilizados en la mecánica de suelos [3] para el análisis de las pruebas de creep y relajación en muestras del cartílago articular, y denominaron al material poroelástico como material “bifásico”. Desde el trabajo de Mow et al [2] se han realizado muchas investigaciones para aplicar las teorías bifásicas o poroelásticas a la descripción del comportamiento mecánico de los tejidos biológicos como el cartílago articular, los tendones y aún los huesos (Armstrong et al [4], Lai et al [5] y Cowín et al [6]).
Una de las motivaciones de estos estudios es conocer los efectos mecánicos en la etiología de las enfermedades. En particular, existen indicios (Donohue et al [7], Newberry et al [8]) de que los impactos en las articulaciones pueden ocasionar microfisuras en el cartílago articular, las cuales alteran los esfuerzos y el cartílago se degrada paulatinamente hasta que desaparece, lo que se denomina osteoartrosis, enfermedad que ocasiona en el paciente la pérdida de movilidad en la articulación y dolores intensos. Los costos asociados con el tratamiento de esta enfermedad fueron estimados en $54 billones de dólares anuales en los Estados Unidos (Atkinson et al [9]). Por tanto, en los países desarrollados se destinan sumas importantes para estas investigaciones. El estudio de esta enfermedad se hace de una manera controlada mediante experimentos con animales; en algunos de estos experimentos se efectúan pruebas de impacto sobre la articulación patelo-femoral de conejos y se estudian luego sus efectos, contando, por ejemplo, el número de condrocitos muertos y haciendo seguimiento a la condición del cartílago (Newberry et al [10]). Para entender las observaciones y tratar de establecer correlaciones con las variables mecánicas es necesario estudiar la distribución de esfuerzos en el tejido.
Las soluciones analíticas con modelos poroelásticos se pueden obtener solamente para geometrías y casos de carga sencillos. Por ello se han propuesto formulaciones mediante elementos finitos, como las de Donzelli y Spilker [11], Spilker et al [12], Almeida y Spilker [13], para resolver problemas más reales. Son pocos los programas comerciales de computador que permiten utilizar materiales poroelásticos. Algunos con esta opción, como Abaqus (Hibbit, Karlsson & Sorensen, Inc, Pawtucket, RI, USA) y Marc (MARC Análisis Research Corporation, Palo Alto, CA, USA) son muy costosos para nuestro medio.
Es por ello que uno de los objetivos de nuestro grupo es el de generar herramientas computacionales propias para estudiar el comportamiento mecánico de los materiales poroelásticos y establecer correlaciones con los experimentos realizados en conejos en el Laboratory of Orthopaedic Research de Michigan State University. En este documento se describen las ecuaciones de la poroelasticidad lineal, luego se explica su formulación en elementos finitos y se presentan los resultados del problema de compresión no confinada para materiales isótropos y transversalmente isótropos.
2. ECUACIONES DE LA POROELASTICIDAD LINEALA continuación se presentan las ecuaciones constitutivas generales que han sido desarrolladas por diferentes investigadores como Mow et al [2], Simon et al [14], Oomens [15]. Un artículo de De Boer [16] presenta una recopilación y revisión general del estado del arte sobre materiales porosos. El material bifásico se representa por la fase sólida y el fluido, los cuales se supone que son intrínsecamente incompresibles. La compresibilidad del tejido completo se explica entonces por la cantidad de fluido que ha salido de él. Los esfuerzos totales en el tejido, sij, se descomponen en los que toma el fluido, sij(f) y los que toma la fase sólida sij(s) así
donde dij es el delta de Kronecker y f (f) es la fracción de volumen del fluido, igual a la relación entre el volumen del fluido y el volumen total. Si se supone que el tejido está saturado, la suma de las fracciones de volumen del fluido y del sólido, f (s) , es igual a 1. Los esfuerzos en la fase sólida sij(s) se descomponen en una fracción de la presión de poros proporcional a f (s) más los denominados esfuerzos efectivos, s ij(e), frecuentemente utilizados para establecer los criterios de falla en la mecánica de suelos
En la descomposición presentada en (1) se considera como si las fases sólida y fluida llenaran simultáneamente todo el volumen, es por ello que en (2) y (3) los esfuerzos en el fluido y sólido toman una fracción de la presión real del fluido o presión de poros p. Si se reemplazan (2) y (3) en (1) se obtiene otra descomposición de los esfuerzos totales en el tejido
.
La ecuación de continuidad o de conservación de masa se puede escribir así
donde la coma indica derivada respecto a las coordenada xi, y vi(s)y vi(f)son las velocidades del sólido y fluido, respectivamente.
También deben cumplirse las ecuaciones de equilibrio para cada fase, las cuales se pueden formular así, para la fase sólida y fluida, respectivamente:
]]>
donde se desprecian los efectos inerciales y se considera que la interacción entre las fases se da por la diferencia entre la velocidad del fluido vi(f) y del sólido vi(s) , donde K es una constante de arrastre.
También se debe establecer una ecuación constitutiva para la fase sólida. Por ejemplo, si se considera un comportamiento lineal elástico tenemos:
donde Cijkm es un tensor de cuarto orden.
Si las seis componentes de esfuerzo y deformación se organizan en vectores, la ecuación (8) puede relacionarse mediante el uso de una matriz [D] 6 x 6, la cual se presenta a continuación en términos de las cinco constantes elásticas de un material transversalmente isótropo:
donde
además E, n y G son el módulo de elasticidad, el coeficiente de Poisson y el módulo a cortante, respectivamente; los subíndices 1 y 2 son direcciones dentro del plano de isotropía y la dirección 3 es perpendicular al plano de isotropía.
]]> Para su implementación en elementos finitos, las ecuaciones anteriores se pueden expresar en términos de los desplazamientos de la fase sólida y la presión de poros de la siguiente forma.Si se suman (6) y (7) y se substituyen (1) y (4) se obtiene:
.
Por otra parte, si se reemplazan en (7) los esfuerzos en el fluido dados en (2), se obtiene la velocidad del fluido en términos de la velocidad del sólido y la presión de poros, así
que es la ley de Darcy según la cual la diferencia entre las velocidades del fluido y sólido es proporcional al gradiente de presión.
El reemplazo de (10) en (5) da como resultado
donde k es la permeabilidad del tejido, definida en términos del coeficiente de arrastre K y de la fracción de fluido así .
y la relación entre la divergencia de la velocidad del sólido y la tasa de deformación infinitesimal
donde el punto arriba del tensor de deformaciones infinitesimales expresa derivada respecto al tiempo, el sistema (9) y (11), de cuatro ecuaciones, se puede expresar y resolver para los tres desplazamientos de la fase sólida y la presión de poros. Un planteamiento igual al anterior se puede obtener a partir de las ecuaciones de poroelasticidad de Biot, presentadas por Smith [17].
3. SOLUCIÓN MEDIANTE ELEMENTOS FINITOS PARA SIMETRÍA AXIAL
3.1. FORMULACIÓN DEL PROBLEMA MEDIANTE RESIDUOS PONDERADOS
A continuación se describe el planteamiento de (9) y (11) mediante el método de desplazamientos virtuales o de residuos ponderados.
Si se multiplica (9) por los desplazamientos virtuales ui (identificados como virtuales por la barra abajo) y se integra sobre el volumen V del cuerpo, se obtiene
.
.
Mediante la aplicación del teorema de divergencia se obtiene
donde son las tracciones definidas sobre el contorno St del cuerpo. Por tanto, el reemplazo de (17) en (15) da como resultad
donde el gradiente de desplazamientos virtual ui,,j se reemplazó por el tensor de deformaciones infinitesimales virtuales, ya que el tensor esfuerzo es simétrico.
Después de la substitución de (14) en (11), la multiplicación por la presión virtual p y su integración sobre el volumen V del cuerpo se obtiene
.
Luego se aplica el teorema de divergencia al segundo término de (19) así:
]]>
donde g es el flujo definido sobre una parte del contorno, el cual, para este tipo de problemas, suele ser cero; lo cual significa que, según esta formulación, se admite flujo sobre los contornos, pero ese flujo no está definido a priori como una condición de borde. Con lo anterior (19) se convierte en
Considerando una primera aplicación para la prueba de compresión confinada, las ecuaciones (18) y (21) se discretizaron espacialmente mediante el método de los elementos finitos como se explica en el apéndice A y se obtuvieron las ecuaciones matriciales
donde [Ks], [C] y [Kp] son matrices definidas en las ecuaciones (A13), donde (u) es el vector de desplazamientos, (p) son las presiones de poros y (u) es el vector de velocidades. La ecuación (23) se multiplicó por (-1) para garantizar la simetría de la matriz de coeficientes.
3.3. DICRETIZACIÓN EN EL TIEMPOLas ecuaciones (22) y (23) contienen el vector de desplazamientos (u) y su derivada respecto al tiempo. Por tanto, para resolver el sistema es necesario aplicar un algoritmo de integración paso a paso en el tiempo, tal como se explica en el apéndice B. Mediante este algoritmo se obtiene el sistema de ecuaciones:
.
donde (Du), (Dp) son los vectores que contienen los cambios de desplazamiento y presión en el intervalo dt. Los resultados de un problema específico se obtienen entonces mediante la solución sucesiva de (24) y (25), para hallar los desplazamientos y presión de poros. Este algoritmo se montó en un programa desarrollado mediante Scilab (http://scilabsoft.inria.fr/), el cual se pone a disposición de la comunidad científica colombiana.
4.VERIFICACIÓN DEL PROGRAMA Y ANÁLISIS PARA COMPRESIÓN NO CONFINADA
]]> 4.1. MÉTODOSCon el objeto de verificar el programa se compararon sus resultados con los obtenidos mediante la solución analítica del problema de compresión no confinada, desarrollada por Armstrong et al [3] para un material isótropo. Se utilizó un modelo con cuatro elementos rectangulares de 1 mm de altura en la dirección axial y 0.5 mm en la dirección radial, de tal forma que el radio total fue de 2 mm y la altura 1 mm (Figura 1). Las propiedades elásticas consideradas se presentan en la Tabla 1.
Figura 1. Malla de elementos finitos.
Figure 1. Finite elements mesh
Tabla 1. Propiedades elásticas de la fase sólida utilizadas en la verificación [18].
Table 1. Elastic properties of solid phase used to verify
La presión de poros en el extremo libre se hizo igual a cero, de tal manera que el flujo se produce en la dirección radial. A este modelo se le aplicó una presión uniforme de 1 KPa en la dirección axial, en la cara superior, y los nodos inferiores fueron restringidos en la dirección axial. La presión se aplicó en una rampa de 1 s y se mantuvo hasta el equilibrio. Este problema se analizó también con un modelo transversalmente isótropo (TI), con las propiedades reportadas en la Tabla 1.
4.2. RESULTADOS Y DISCUSIÓN
]]> Se observó una buena correlación entre los resultados del programa y los analíticos (Figuras 2 y 3). El desplazamiento vertical del nodo 1 (ubicado en la superficie y en el centro) tiende al valor esperado de 0.14 mm cuando se alcanza el equilibrio. Este desplazamiento es igual al cociente de la presión externa (1KPa) sobre el módulo de elasticidad de la fase sólida (0.69 MPa), ya que en equilibrio la presión de poros es cero y toda la carga externa la toma la fase sólida. En la primera parte de las curvas de presión (Figura 3) se observan algunas diferencias que podrían ser atribuidas a que en la solución analítica se aplica la carga instantáneamente y en el programa se aplica en una rampa de un segundo de duración.
Figura 2. Desplazamiento vertical del nodo 1.
Figure 2. Vertical displacement of node 1
Figura 3. Presión de poros en el nodo 1.
Figure 3. Porous presion in node 1
Durante la aplicación rápida de la carga la presión de poros aumenta rápidamente (Figura 3) y ayuda de manera significativa a soportar la carga externa. Este modelo corrobora así la función protectora del fluido cuando se aplican cargas rápidas sobre el cartílago articular, como las producidas por impactos en las articulaciones ocasionadas por saltos o caídas.
El desplazamiento inicial vertical del modelo TI fue aproximadamente el 40% del desplazamiento del modelo isótropo (Figura 4); y el pico de la presión de poros en el modelo TI fue aproximadamente tres veces mayor que el observado en el modelo isótropo (Figura 5). Estas diferencias, documentadas en la literatura (García et al. [19]), se deben a que el modelo TI es significativamente más rígido en la dirección radial, dado el mayor valor del módulo de elasticidad en el plano de isotropía (5.8 Mpa vs 0.68Mpa), lo cual implica que la deformación lateral inicial del modelo TI es menor que la del modelo isótropo, por lo cual existe un mayor confinamiento para el fluido y se genera una mayor presión de poros. Este modelo TI para el cartílago articular representa una mejoría con respecto al modelo isótropo ya que, como se ha documentado recientemente (Wang et al [20]), las propiedades a tensión en la dirección superficial son mucho mayores que en la dirección perpendicular al elemento.
]]>
Figura 4. Desplazamiento vertical nodo 1
Figure 4. Vertical displacement of node 1
Figura 5. Presión de poros en nodo 1
Figure 5. Porous presion in node 1
En estos materiales poroelásticos, el creep, o aumento de la deformación con el tiempo del tejido sometido a carga, se puede explicar por el movimiento de fluido dentro del tejido. Inicialmente, cuando se aplica una carga rápida y dada la baja permeabilidad del tejido, el fluido no puede escapar del esqueleto sólido y adquiere entonces una presión importante. A medida que transcurre el tiempo, el flujo hace que la presión disminuya paulatinamente hasta que llega a cero. Esta es la situación de equilibrio en la cual la fase sólida toma toda la carga externa.
5. CONCLUSIONESSe describieron las ecuaciones generales para materiales bifásicos presentadas por Mow et al [2] y, mediante la discretización espacial con elementos finitos y la solución en el tiempo mediante un procedimiento de integración paso a paso, se generó un programa válido para deformaciones infinitesimales con el cual se obtuvieron resultados muy similares a los dados por la solución analítica de Armstrong et al [4]. La rigidez inicial del tejido sometido a carga rápida es mayor en el modelo TI que en el isótropo, lo cual se debe a la generación de una presión de poros mayor por el mayor grado de confinamiento en la dirección transversal.
Este trabajo es un primer paso hacia el desarrollo de un programa más general que considere inicialmente deformaciones finitas, como las que puede experimentar un tejido biológico durante su funcionamiento, y modelos no lineales para el comportamiento de la fase sólida, como el hiperelástico o hipoelástico. En una etapa posterior también deben incorporarse los efectos viscoelásticos de la fase sólida, que pueden ser importantes en ciertas condiciones de los tejidos biológicos, como lo han demostrado DiSilvestro et al [21] para el cartílago articular.
APÉNDICE A
En un elemento cuadrilátero de cuatro nodos, las funciones desplazamiento u y v, en las direcciones radial y axial, respectivamente, se aproximan como
]]>
donde N1, N2, N3, N4 son las funciones de forma o interpolación y u1, u2, u3, u4, v1, v2, v3 y v4 son los desplazamientos nodales en las direcciones radial y axial; la letra u aplica para la dirección radial y la letra v para la axial. Además, los corchetes encierran matrices y los paréntesis encierran vectores. La presión de poros se aproxima en (24) mediante las mismas funciones de forma y los valores nodales de presión contenidos en el vector (p).
Para el caso de simetría axial las componentes no triviales del tensor de de deformaciones infinitesimales son las deformaciones ez, er y eq y la angular grz, las cuales se organizan en el vector (e). Para coordenadas cilíndricas este vector se calcula así en términos de los desplazamientos:
El tensor esfuerzo en su forma vectorial se relaciona con las deformaciones infinitesimales como aparece en (A4)
donde [D] se expresa en la ecuación (8-a).
La traza del tensor deformación se calcula así para coordenadas cilíndricas
La derivada respecto al tiempo de la traza del tensor infinitesimal se obtiene así
donde la matriz [B]2 es la definida en (A7). El reemplazo de las ecuaciones (A1) - (A6) en (18) y (21) da como resultado el sistema de ecuaciones (A9) y(A10).
donde el vector contiene las tracciones especificadas sobre una parte del dominio. El volumen Ve se refiere al del elemento. Las ecuaciones (A9) y (A10) deben cumplirse para cualquier conjunto de vectores virtuales (u) y (p). Por tanto:
donde:
Se divide el tiempo en intervalos dt y se denominan (u)(i) y (p)(i) los vectores de desplazamiento y presión, respectivamente, en el tiempo ti. La idea consiste en determinar los vectores de desplazamiento y presión en el tiempo (i+1) a partir de los valores conocidos en el tiempo (i). Por tanto en el tiempo (i+1) se deben satisfacer las ecuaciones
Se supone ahora que la velocidad en cada intervalo es constante y obedece a la ecuación:
donde t es el tiempo que varía desde cero hasta dt y q es un parámetro entre cero y uno. Si es cero, la velocidad es igual a la inicial, y si es uno, la velocidad es igual a la final.
Mediante integración en el tiempo de la ecuación (B3) se obtienen los desplazamientos
Si se evalúa (B4) en t = dt se obtiene
Y este resultado se substituye en (B2)
Si se suma a ambos lados de la ecuación anterior el término [Kp](p)(i) se simplifica y reorganiza se obtiene:
Donde (Du), (Dp) son los vectores que contienen los cambios de desplazamiento y presión en el intervalo dt.
En forma incremental, (A10) queda así
.
REFERENCIAS
]]>[1] R. de Boer, ”Highlights in the historical development of the porous media theory: Toward a consistent macroscopic theory”, Appl Mech Rev 32, pp. 201-262, 1996. [ Links ]
[2] V.C. Mow, S.C. Kuei, W.M. Lai y C.G. Armstrong CG, “Biphasic creep and stress relaxation of articular cartilage in compression: theory and experiments”, J Biomech Eng (ASME) vol. 102, pp. 73-84, 1980. [ Links ]
[3] M.A. Biot “Theory of elasticity and consolidation for a porous anisotropic solid”, J Appl Phys, vol. 26, pp.182-185, 1955. [ Links ]
[4] C.G. Armtrong, W.M. Lai y V.C. Mow, “An analysis of the unconfined compression of articular cartilage”, J Biomech Eng, vol. 106, pp. 165-173, 1984. [ Links ]
[5] W.M. Lai, J.S. How, y V.C. Mow VC “A triphasic theory for the swelling and deformation behaviors of articular cartilage”, J Biomech Eng, vol. 113, pp. 245-258, 1981. [ Links ]
[6] S.C. Cowin, S. Weinbaum, y Y. Zeng Y “A case for bone canaliculi as the anatomical site of strain generated potentials”, J Biomech, vol. 28, pp.1281-1297, 1995. [ Links ]
[7] J.M.Donahue, D. Buss, T.R. Oegema y R.C. Thompson “The effects of indirect blunt trauma on adult canine articular cartilage”, The Journal of Bone and Joint Surgery, vol. 65(A), pp. 948-957, 1983. [ Links ]
[8] W.N.Newberry, D. Zukosky, R.C.Haut,” Subfracture insult to a knee joint causes alterations in bone and in the functional stiffnes of overlying cartilage”, Journal of Orthopaedic Research, vol. 15, 450-455, 1997. [ Links ]
[9] P.J. Atkinson, J.J.García, N.J. Altiero y R.C. Haut “The influence of impact interface on human knee injury: implications for instrument panel design and the lower extremity injury criterion”, 41st Stapp Car Crash Conference Proccedings, 1997. [ Links ]
[10] W.N. Newberry, C.D. Mackenzie y R.C. Haut, “Blunt impact causes changes in bone and cartilage in regularly exercised animal model”, Journal of Orthopaedic Research, vol. 16, pp.348-354, 1998. [ Links ]
[11] P.S. Donzelli PS y R.L. Spilker, “A finite element formulation for contact of biphasic materilals: evaluation for plane problems”, 1993 Advances in Bioengineering, J. Tarbell, ed, New York, ASME, pp. 47-50, 1993. [ Links ]
[12] R.L. Spilker, E.S. Almeida y P.S. Donzelli, “Finite element methods for the biomechanics of soft hydrated tissues: non linear analysis and adaptive control of meshes”, High Performance Computing in Biomedical Research, T.C. Pñikington et al., ed, Boca Raton, Florida, USA, CRC Press, pp. 227-261, 1993. [ Links ]
[13] E.S. Almeida y R.L. Spilker, “Mixed and penalty finite element models for the nonlinear behavior of biphasic soft tissues in finite deformation: Part I-Alternative formulations”, CMMBBE, vol. 1, pp. 25-46, 1997. [ Links ]
[14] B.R. Simon, M.V Kufmann, M.A. McAfee y A.L. Baldwin, “Porohyperelastic theory and finite element models for soft tissues with application to arterial mechanics”, P.S. Selvadurai (Ed), Mechanics of Poroelastic Media, 245-261, 1996. [ Links ]
[15] C.W. Oomens, D.H. Van Campen y H.J.Grootenboer HJ, ”A mixture approach to the mechanics of skin, Journal of Biomechanics, vol. 20, pp. 877-885, 1987. [ Links ]
[16] R. de Boer, ”Contemporary progress in porous media theory”, Appl Mech Rev, vol. 53, pp. 323-370, 2000. [ Links ]
[17] I.M. Smith, D.V.Griffiths DV, Programming the finite element method, Third Edition, John Wiley & Sons, Chichester, UK, 1998. [ Links ]
[18] B. Cohen, W.M. Lai, V.C. Mow, “A transversely isotropic biphasic model for unconfined compression growth plate and chondroepiphysis”, ASME Journal of Biomedical Engineering, Vol. 120, pp 491-496. 1998. [ Links ]
[19] J.J. García, N.J. Altiero y R.C.Haut RC,”Estimation of in situ elastic properties of biphasic cartilage based on a transversely isotropic hypo-elastic model”, Journal of Biomechanical Engineering, vol. 122, pp.1-8, 2000. [ Links ]
[20] C. C-B. Wang, N.O Chahine, C. T. Hung, G.A. Ateshian, “Optical determination of anisotropic material properties of bovine articular cartilage in compression”, Journal of Biomechanics, vol. 36, pp. 339-353. 2003. [ Links ]
[21] M.R. DiSilvestro, Q.Zhu, M. Wong, J.S. Jurvelin, J.F. Suh, “Biphasic poroviscoelastic simulaton of the unconfined compression of articular cartilage: I-Simultaneous prediction of reaction force and lateral displacement”, ASME Journal of Biomechanical Engineering, vol. 123, pp. 191-197, Abril, 2001. [ Links ] ]]>