A 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, los esfuerzos en el fluido, que son esfuerzos promedio sobre todo el volumen, son iguales a una fracción de la presión real del fluido, o presión de poros p

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 v_i^(s)y v_i^(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 v_i^(f) y del sólido v_i^(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 C_ijkm 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í .

Dadas las relaciones entre las deformaciones infinitesimales y los desplazamientos

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

.

Esta integral se puede reorganizar así

.

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.

Las 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ÉTODOS

Con 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.

Se 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 e_z, e_r y e_q y la angular g_rz, 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

con lo que su aproximación es

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:

APÉNDICE B

Se divide el tiempo en intervalos dt y se denominan (u)⁽ⁱ⁾ y (p)⁽ⁱ⁾ los vectores de desplazamiento y presión, respectivamente, en el tiempo t_i. 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

De (B5) se despeja la velocidad en el tiempo (i+1)

Y este resultado se substituye en (B2)

Si se suma a ambos lados de la ecuación anterior el término [Kp](p)⁽ⁱ⁾ 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 ]