Localización de defectos en hueso cortical empleando el método de los elementos de contorno y algoritmos genéticos.

David Ojeda¹, Eduardo Divo², Alain Kassab³, Miguel Cerrolaza⁴

¹ Departamento de Diseño Mecánico y Automatización, Facultad de Ingeniería, Universidad de Carabobo, Naguanagua, Venezuela, dojeda@uc.edu.ve

² Department of Engineering Technology, University of Central Florida, Orlando, FL, USA, edivo@mail.ucf.edu

³ Dep. of Mechanical, Materials, and Aerospace Engineering, University of Central Florida, Orlando, FL, USA, kassab@mail.ucf.edu

⁴ Instituto Nacional de Bioingeniería, Universidad Central de Venezuela, Caracas, Venezuela, mcerrola@gmail.com

Resumen

Se presenta una solución eficiente para la detección de cavidades en hueso cortical a partir de una técnica de superposición de cargas puntuales empleando el método de elementos de contorno. Se muestran dos ejemplos donde se logra la detección de una cavidad a partir de una imagen tomográfica de un hueso cortical, validándose la solución a partir de resultados numéricos experimentales usando el Método de los Elementos de Contorno.

Palabras claves: Método de los elementos de contorno (MEC); detección de defectos; algoritmos genéticos.

Locating defects in cortical bone using the boundary element method and genetic algorithms.

ABSTRACT

An efficient solution for cavity detection in cortical bone, using a point load superposition technique employing the boundary element method is presented in this paper. Two examples where the cavity is located using a TAC is presented. Results of cavity presence problems simulated using numerical experiments validate the approach in cortical bone with a single cavity.

Keywords: Cortical bone; boundary element method (BEM); defect detection; genetic algorithms.

Recibido: 13/06/07  Revisado: 08/10/07 Aceptado: 26/11/07

1. Introducción

Se presenta la resolución de un problema de ensayo no destructivo de materiales que permite la localización de una cavidad en una sección transversal de hueso cortical perteneciente a un fémur humano, empleando el método de los elementos de contorno y algoritmos genéticos. Se utiliza un método de superposición singular de cargas puntuales para resolver el problema geométrico inverso permitiendo simular, de manera eficiente, la presencia de cavidades internas (Ojeda et al 2007a; Gámez et al 2007).

Para determinar el campo de variables que se obtienen de las entradas, en un problema directo se especifica: a. La ecuación que rige el campo variable, b. Las propiedades físicas, c. Las condiciones de contorno, d. Las condiciones iniciales y e. La geometría del sistema. Por el contrario en un problema inverso, que es el caso de esta publicación, se especifica: (1) parte de las condiciones anteriormente

numeradas (a, b, c, d, e) y (2) una condición especificada (conocida como la condición expuesta del contorno). La razón de resolver un problema inverso es encontrar alguna condición desconocida numerada de la (a) a la (e) en un problema directo. Normalmente en un problema inverso, las condiciones especificadas se obtienen de alguna medición interna del campo de variables a través de algún sensor (Divo et al 2000, 2004; Ulrich et al 1996; Kassab et al 1994; Kassab et al 1995) o, como el caso de esta investigación, por medio del resultado de un análisis numérico.

El método de los elementos de contorno (MEC), por su propia naturaleza, ofrece la alternativa perfecta (Ojeda et al 2007c; Gámez et al 2007; Divo et al 2000, 2004; Ulrich et al 1996; Kassab et al 1994; Kassab et al 1995; Cerrolaza et al 2000; Müller et al 2001; Annicchiarico et al 2002; Annicchiarico et al 2004; Annicchiarico et al 2005; Martínez et al 2006) para la resolución de este tipo de problemas debido a que no requiere discretización del dominio ni regeneración completa del mallado así como la de su geometría (Brebbia y Domínguez et al 1989; Kane 1994).

En este artículo se determinaron las características de las cavidades encontradas en una sección transversal de la diáfisis de fémur, mediante el desarrollo de un método de superposición de singularidades (Ojeda et al 2007a, 2007b, 2007c; Gámez et al 2007).

Este método propone una solución para la localización de cavidades o tumores en zonas cercanas a la presencia de un dispositivo metálico, como un clavo, fijador externo o placa, que produzca brillo cuando se realice una tomografía axial computarizada o rayos X.

2. PROBLEMA DIRECTO Y MEC EN ELASTICIDAD

La solución de problemas directos en elastostática provee el campo de los desplazamientos dentro del sistema de ecuaciones de gobierno, las condiciones de contorno y el sistema geométrico. El campo de los desplazamientos para un medio isótropo, homogéneo y linealmente elástico sujeto a una distribución de fuerzas internas o de campo b_i es regida por la ecuación de Navier:

   (1)

donde u_i y u_j representan los campos de desplazamientos, n es la relación de Poisson y m es el módulo de rigidez. En un problema directo, la solución de esta ecuación, normalmente, es posible cuando se prescriben correctamente dos posibles tipos de condiciones de borde:

   (2)

donde, la sobrebarra denota una cantidad especificada, t_i es el campo de tracciones y n_j es el vector normal unitario hacia afuera del dominio W rodeado por el contorno G = G_u U G_t. El campo de esfuerzos s_ij está definido como:

   (3)

Con el campo de deformaciones relativas e_ij como:

  (4)

La técnica del método de los elementos de contorno consiste en transformar las ecuaciones que describen el comportamiento de la variable en el interior y en el contorno del dominio, en ecuaciones integrales referidas solamente al contorno y encontrar la solución numérica de estas ecuaciones. De manera tal, que todas las aproximaciones numéricas se efectúan para valores en el contorno, trayendo como consecuencia la reducción de la dimensionalidad del problema (Inglessis y Cerrolaza 2003). En una formulación del método de elementos de contorno basada en la ecuación integral de Somigliana, queda planteada una relación integral entre la deformación en un punto de colocación “p” y las deformaciones u_i y tracciones t_i en todo el contorno G. Adicionalmente, las fuerzas de campo b_i, permanecen relacionadas a través de una integral de dominio W, como sigue:

  (5)

donde, G_ij y H_ij son las soluciones fundamentales de desplazamientos y tracciones, respectivamente, y  es una constante geométrica, tal que:  si  p Î W ; si  p Ï W y  si p Î G en el caso de un contorno suave (Brebbia y Domínguez, 1989; Kane, J., 1994).

Estableciendo que el campo de fuerzas internas b_i está constituido únicamente por cargas puntuales (Ojeda et al 2007b), tal que:

  (6)

donde NL es el número de cargas puntuales, es la intensidad de cada carga l y  es la función de Dirac ubicada en el punto de impacto de cada carga . Utilizando las propiedades de la función de Dirac, el último término de la ecuación integral (5) se reduce a:

(7)

donde  es la solución fundamental de desplazamiento evaluada en el punto de impacto de cada carga x^l. Esta substitución permite reducir, de manera simple, la ecuación integral (7) a integrales de contorno como:

(8)

De igual forma se puede definir una ecuación integral de contorno que relaciona el campo de esfuerzos  en un punto de colocación p como:

   (9)

donde los términos  y  son las soluciones fundamentales de esfuerzos. Esta ecuación integral de contorno se puede emplear para calcular el esfuerzo  en cualquier punto p del dominio W e incluso se puede llevar al contorno pero se deben utilizar métodos especiales de integración ya que el término  se vuelve hipersingular cuando se integra sobre el punto de colocación p.

Discretizando la ecuación integral de contorno de desplazamientos (9) utilizando NE elementos de contorno con NN nodos independientes en cada elemento, se obtiene el siguiente arreglo matricial:

[H]{u}=[G]{t}+{q}   (10)

donde las matrices [H] y [G], de dimensiones N´N, contienen los coeficientes de influencia que relacionan los desplazamientos {u} y tracciones {t} en el contorno (Brebbia y Domínguez, 1989; Kane, J., 1994). Donde N = d´NE´NN, d es el número de dimensiones espaciales (2 o 3), NE es el número de elementos y NN es el número de nodos. Lo más resaltante de este desarrollo es el hecho que todos los efectos generados por las cargas puntuales quedan remitidos al vector {q}, lo que esencialmente permite alegar que si de alguna forma se logra simular la presencia de una cavidad utilizando un conjunto de cargas puntuales, entonces no se necesitaría re-discretizar el contorno en cada paso de un eventual proceso de optimización para la detección de esta cavidad ya que  su efecto  simplemente estaría afectando al sistema algebraico con la suma de un vector {q} conocido.

A la relación de coeficientes (10) se introducen las condiciones de borde  y  para arreglar el sistema algebraico en la forma:

[A]{x} ={b}+{q}   (11)

donde el vector {x} contiene los valores nodales desconocidos de {u} y {t} en donde no fueron prescritas condiciones de borde (Brebbia y Domínguez, 1989; Kane, J., 1994). La solución de este sistema sigue métodos estándar. Cabe destacar que para esta investigación se emplearon solo elementos de contorno cuadráticos discontinuos iso-paramétricos, donde tanto la geometría como las cantidades de campo {u} y {t} son aproximadas utilizando funciones de forma cuadráticas pero los nodos de las cantidades de campo se ubican dentro del elemento un 12.5%, (Divo et al 2000, 2004; Ulrich et al 1996; Kassab et al 1994; Kassab et al 1995) y no en los puntos extremos del elemento como es el caso de los nodos geométricos. Una de las ventajas para el uso de este tipo de nodo es para evitar las singularidades en las esquinas (Brebbia y Domínguez, 1989).

3. MÉTODO DE OPTIMIZACIÓN: ALGORITMOS GENÉTICOS (AGs)

Tres características resaltantes que definen a los algoritmos genéticos (AGs) de otros algoritmos evolutivos: (1) representación binaria de la solución, (2) la proporcionalidad del método de selección, y (3) mutación y cruce como método primario para producir variaciones. En general, las propiedades de un organismo se describen como una serie de genes (variables de diseño) en un cromosoma. Por lo tanto, si se trata de simular la naturaleza usando un computador, se debe codificar el diseño de las variables de una manera conveniente. Se adoptará un modelo haploide usando un modelo de vector binario que simule un cromosoma sencillo (Ojeda et al 2007a, 2007b, 2007c; Gámez et al 2007; Divo et al 2000, 2004). Se determina la longitud del vector por el número y precisión requerida de cada variable de diseño. Cada variable de diseño es limitada entre un valor máximo y mínimo. Adicionalmente, se establece la precisión de la variable de diseño asignando un número de bits por variable. La información contenida en la cadena de vectores que compone el cromosoma caracteriza a cada individuo en su población. Además, cada individuo, dependiendo de su diseño, posee una aptitud o capacidad de supervivencia. El proceso de optimización comienza con la suposición aleatoria de un grupo de posibles soluciones denominadas población que contienen un tamaño inicial o número de individuos. Cada individuo se define por la combinación de parámetros  y se representa como una

cadena binaria o cromosoma, ver figura 1. Ya que los AGs están planteados para maximizar y no minimizar, se puede formular una función de aptitud, Z (Ojeda et al 2007a), como el inverso de la función objetivo S como:

(12)

Esta función de aptitud Z se evalúa para cada individuo en la población actual definiendo que tan apto es este individuo para la supervivencia. En cada iteración de los AGs se realiza una operación de selección, cruce y mutación para actualizar el arreglo de la población. La selección se realiza para determinar pares de individuos que cruzarán su material genético para producir un descendiente para la nueva generación. Este operador selecciona aleatoriamente a la pareja de individuos, sin embargo, esta aleatoriedad es proporcionada con el nivel de aptitud de cada individuo para garantizar que los individuos más aptos tendrán más probabilidad de transmitir su carga genética a la nueva generación. Adicionalmente, un operador de mutación aleatorio afecta la información obtenida al cruzar el material genético de los individuos. Este paso es importante para un mejoramiento continuo de las generaciones. Una serie de parámetros son fijados inicialmente en los AGs y determinan y afectan la eficiencia del proceso de optimización.

Los parámetros que controlan el proceso de optimización son el número de parámetros por individuos o variables de optimización, el número de bits por parámetro o longitud de los cromosomas que define cada individuo, el número de individuos o tamaño de la población por generación, el número de hijos por cada cruce, la probabilidad de cruce y la probabilidad de mutación. Los operadores de los AGs se aplican generación tras generación hasta satisfacer un criterio de convergencia o hasta alcanzar un número máximo de generaciones. Una vez que el tamaño de la población sea fijada, el algoritmo inicializa todos los cromosomas. Esta operación se lleva a cabo aleatoriamente. Después de inicializar la población, se realiza una evaluación de la aptitud de cada individuo computando la función de aptitud. Teniendo los valores de la función de aptitud para cada individuo, el proceso de selección puede comenzarse. La probabilidad de selección de cada individuo es calculada utilizando la relación mostrada en la ecuación 13.

  (13)

donde TP es el tamaño de la población, ui es el miembro i de la población y Zi es la aptitud del miembro i de la población. Se genera una ruleta de peso en la cual a cada miembro se le asigna una porción de la rueda en proporción a su probabilidad de selección. La ruleta es girada tantas veces como hayan individuos en la población para producir el número adecuado de descendientes. Acá algunos cromosomas serán seleccionados más de una vez y los mejores cromosomas tendrán más oportunidades de transmitir su material genético a la nueva generación. El cruce y mutación ocurren una vez que se haya aplicado la selección para garantizar la diversidad genética en la población de la nueva generación. Una vez estimadas las cargas puntuales mediante el proceso de optimización de los AGs, un nuevo proceso de optimización se inicializa y lleva a cabo para ubicar la línea o contorno de cero tracciones al rededor de las cargas. Este segundo proceso de optimización se realiza utilizando los AGs con las mismas propiedades de la etapa anterior pero con una nueva función (Ojeda et al 2007a) de aptitud tal como:

  (14)

La optimización de esta función de aptitud utilizando los AGs da como resultado la ubicación y forma final aproximada de la superficie de la cavidad en cuestión.

4. METODOLOGÍA

Para la localización de la cavidad, primero se somete el hueso a cargas externas de tracción o compresión. Se miden los desplazamientos en el contorno y se introducen en el programa de algoritmos genéticos de acuerdo al esquema mostrado en la figura 2. Puede observarse que el procedimiento lee los desplazamientos a partir de un archivo, posteriormente distribuye las cargas puntuales que simulan la presencia de cavidad en el dominio. Carga la función evaluación e inicia la evaluación de la primera función objetivo (ecuación 12) hasta obtener un valor menor a un error fijado previamente (e). Una vez concluido esto, se continúa con el siguiente paso que corresponde a la determinación de la forma de la cavidad. Para ello se carga la segunda función objetivo (ecuación 14) para determinar las tracciones cero, concluyendo con la ubicación y forma del defecto.

5. VALIDACIÓN NUMÉRICA

Los resultados numéricos de problemas directos en 2D con cavidades usando el MEC ofrecen una alternativa para generar las condiciones de contorno en el exterior de la superficie que servirá para validar el programa. Se usan elementos de contorno discontinuos y cuadráticos en todos los casos.

La figura 3.a muestra un primer ejemplo. Corresponde a una sección transversal de un hueso cortical perteneciente a un fémur de una persona adulta. Las dimensiones en promedio son 0,04 m de diámetro externo y diámetro interno de 0,0254 m. Posee una cavidad elíptica rotada 45º, cuyos ejes mayor y menor son 0,0045 m y 0,0026 m, respectivamente. En la figura 3.b se observa la discretización de la sección con 128 elementos (80 en los lados y 48 para la sección interna). Las condiciones de borde establecidas en el modelo se observa en la figura 3.c. El hueso está empotrado en uno de los lados y tiene una carga de compresión de 100 MPa distribuida en el resto de los lados.

Un segundo ejemplo se muestra en la figura 4.a. Las dimensiones del hueso similares al ejemplo 1, con la diferencia de que la cavidad es circular con un diámetro aproximado de 0,002 m. Se observa en la mitad de la sección un empotramiento y la otra mitad una carga distribuida perpendicular al hueso de 10 MPa. En la figura 4.b se muestra la discretización del hueso similar al ejemplo 1. Las figuras 5.a y 5.b muestran los desplazamientos para ambos ejemplos, usando problemas directos y MEC. Estos valores serán usados para validar el proceso de simulación (Ojeda, D., et al 2007a, 2007b, 2007c; Gámez, B., et al 2007).

Después de resolver el problema inverso para el ejemplo 1, se observa en las figuras 6.a y 6.b, los resultados de la evolución de las cargas puntuales desde la primera generación hasta lograr una convergencia, detectando finalmente la ubicación de la cavidad. Esta convergencia es alcanzada por la primera función objetivo (ecuación 12) después de 1000 generaciones del AGs. Posteriormente se obtiene la convergencia de la segunda función objetivo (ecuación 14) después de 50 generaciones, permitiendo obtener la forma de la cavidad. La evolución de las dos funciones objetivos para el ejemplo 1 se muestran en las figuras 7.a y 7.b.

Similar ocurre con el ejemplo 2. En las figuras 8.a y 8.b se observa el movimiento de las cargas puntuales, localizando la cavidad, hasta alcanzar la convergencia en 500 generaciones del AGs. La forma de la cavidad se obtiene después de 20 generaciones de AGs. La evolución de las dos funciones objetivos se observa en las figuras 9.a y 9.b.

6. CONCLUSIONES

Se ha propuesto un método superposición singular de cargas puntuales para la solución de problemas geométricos inversos que localiza de manera eficiente las cavidades en problemas que involucra una configuración geométrica irregular como lo es la sección transversal de la diáfisis de hueso de fémur. Se ha integrado el método de los AGs como herramienta de optimización con el MEC, validando la precisión de la solución con una cavidad simple de forma circular y una elíptica, sometido a dos diferentes arreglos de cargas de compresión. Los resultados numéricos han provenido de los problemas directos haciendo uso del MEC, aprovechando la ventaja de que con el mismo no se requiere remallado en el interior del dominio; lo que ha permitido la reducción del tiempo computacional.

7. AGRADECIMIENTOS

Los autores quieren expresar su más sincero agradecimiento a la Universidad de Carabobo (Venezuela), University of Central Florida (USA) e Instituto Nacional de Bioingeniería (Venezuela) por el financiamiento y apoyo institucional para el desarrollo de esta investigación.

8. REFERENCIAS BIBLIOGRÁFICAS

1. Annicchiarico W., Cerrolaza M., (2004). “A 3D boundary element optimization approach based on genetic algorithms and surface modelling”, Eng. Anal. With Bound. Elem, Vol. 28(11), pp. 1351-1361.        [ Links ]

2. Annicchiarico W., Cerrolaza M., (2002). “An Evolutionary Approach for the Shape Optimization of General Boundary Elements Models”. Electronic Journal of Boundary Elements, Vol.2.        [ Links ]

3. Annicchiarico W., Martínez G., Cerrolaza M., (2005) “Boundary elements and b-spline modelling for medical applications”, J. of App. Math. Mod. Vol. 31(2), p.p. 194 – 208.        [ Links ]

4. Brebbia C.A., Dominguez J., (1989). “Boundary element, An introductory course”. Computational mechanics Publ., Boston, co-published with McGraw-Hill, New York. pp. 134-250        [ Links ]

5. Cerrolaza M., Annicchiarico W. and Martinez M., (2000). “Optimization of 2D boundary element models using b-splines and genetic algorithms”, Engineering Anal. with Bound. Elem, 24(5): 427-440.        [ Links ]

6. Divo E.A., Kassab A.J., Rodríguez F., (2004). “An efficient singular superposition technique for cavity detection and shape optimization”. Numerical Heat Transfer, Part B, 46: 1 - 30, Copyright Taylor & Francis Inc.        [ Links ]

7. Divo E.A., Kassab A.J., Rodríguez F., (2000).“Characterization of space dependent thermal conductivity with a BEM-Based genetic algorithm”. Institute for Computational Engineering, University of Central Florida, Orlando, Florida, 32816-2450, Numerical Heat Transfer, Part A: Applications, Vol. 37, No. 8, pp 845 - 877.        [ Links ]

8. Gámez, B., Ojeda, D., Divo, E., Kassab, A., Cerrolaza, M. (2007). “Cavity detection and crack analysis in biomechanics using BEM”, Proceeding of IX International Conference on Computational Plasticity COMPLAS IX. E. Oñate and D. R. J. Owen (Eds) Ó CIMNE, Barcelona.        [ Links ]

9. Inglessis P., Cerrolaza M. (2003). Una introducción al método de los elementos de contorno en ingeniería y ciencias aplicadas. CDCH-UCV, Caracas, Venezuela        [ Links ]

10. Kassab A.J., Moslehy F.A., Daryapurkar A.B., (1994). “Nondestructive detection of cavities by an inverse elastostatics boundary element method”. J. Engineering Analysis with Boundary Elements 13. 45 - 55.        [ Links ]

11. Kassab A.J., Moslehy F.A., Ulrich T.W., (1995). “Inverse boundary element solution for locating subsurface cavities in thermal and elastostatic problems”. In Proc. IABEM-95, Computational Mechanics '95, Hawaii, July 30-August 3 (ed. Atluri, Yagawa and Cruse), pp. 3024-3029, Springer, Berlin.        [ Links ]

12. Kane J.H. (1994). “Boundary element analysis in engineering continuum mechanics”. Prentice Hall, Englewood Cliffs, New Jersey 07632.        [ Links ]

13. Martínez G. and Cerrolaza M., (2006). “A bone adaptation integrated approach using BEM”, J. Eng. Anal. With Bound. Elem. 30, 107-115.        [ Links ]

14. Müller-Karger C., González C., Aliabadi M.H. and Cerrolaza M., (2001). “Three dimensional BEM and FEM stress analysis of the human tibia under pathological conditions”, J. of Comp. Mod. In Eng. and Sciences, 2(1): 1-13.        [ Links ]

15. Ojeda, D., Divo, E., Kassab, A. Cerrolaza, M., (2007a). “Cavity detection in biomechanics by an inverse evolutionary point load bem technique”. Inverse Problems, Design and Optimization Symposium (IPDO-2007), Miami, Florida. USA.        [ Links ]

16. Ojeda, D., Divo, E., Kassab, A., Cerrolaza, M., (2007b). “Superposición de singularidades para similar la presencia de cavidades en problemas de elastostática usando el método de los elementos de contorno”. In Press, Acta Científica Venezolana. IVIC, Caracas, Venezuela.        [ Links ]

17. Ojeda, D., Gámez, B., Divo, E., Kassab, A. and Cerrolaza, M., (2007c). “Singular superposition elastostatic BEM/GA algorithm for cavity detection”, Brebbia, C. and Popov, V. (eds), Boundary elements and other mesh reduction methods XXIX, pp. 313 – 322, WIT Press, UK.        [ Links ]

18. Ulrich T.W., Moslehy F.A. Kassab A.J., (1996). “A BEM based pattern search solution for a class of inverse elastostatic problems”. Int. J. Solids Structures, Vol. 33, No. 15, pp. 2123-2131.        [ Links ]