Modelado numérico de la propagación de fractura por fatiga utilizando el método de los elementos de contorno

Jairo F. Useche Vivero^a

Resumen

El trabajo presenta una herramienta computacional basada en el método de los elementos de contorno, (boundary element method –BEM–), que permite estudiar el fenómeno de propagación de fractura en estructuras sometidas a esfuerzos cíclicos. El código sirve para desarrollar modelos computacionales enfocados al estudio de resistencia residual de componentes bidimensionales fracturados, basados en la mecánica de la fractura lineal elástica, (linear elastic fracture mechanics –LEFM–). Se presenta la formulación por elementos de contorno del problema, su implementación computacional y su aplicación al análisis de problemas específicos de crecimiento de fractura por fatiga. El crecimiento de la fractura es modelado a través de la teoría propuesta por Paris, aplicada al caso de elasticidad plana, y es usado el criterio de esfuerzos principales máximos para definir la dirección de propagación (Portela, 1992). La formulación BEM fue implementada en un código computacional desarrollado en Matlab 7.0 . El programa ha demostrado efectividad para predecir trayectorias de propagación y resistencia residual por fatiga, permitiendo su uso como herramienta de análisis en ingeniería.

Palabras clave: Propagación de fractura, método de elementos de contorno, mecánica de la

Abstract

The work presents a computer tool based on the boundary element method (BEM), that allows the study of the propagation of fracture phenomenon in structures which have sustained cyclical efforts. The code functions to develop computer models focused on the study of residual resistance of fractured bi-dimensional components based on linear elastic fracture mechanics (LEFM). The boundary element formulation of the problem is presented, as well as its computer implementation and its application to the specific problem analysis of fracture growth by fatigue. Fracture growth is modeled through the theory proposed by Paris, applied to the case of flat elasticity, and the criteria of maximum main efforts is used to define the direction of propagation (Portela, 1992). The BEM formulation was implemented in a computer code developed in Matlab 7.0 . The program has demonstrated effectiveness in predicting trajectories of propagation and residual resistance by fatigue, allowing its use as a tool of analysis in engineering.

Key words: Fracture growth propagation, boundary element method, fracture mechanics.

Fecha de recepción: 12 de octubre de 2007
Fecha de aceptación: 26 de noviembre de 2007

________________________

^aUniversidad Tecnológica de Bolívar. Departamento de Ingeniería Mecánica.
Correo electrónico: juseche@unitecnologica.edu.co

............................................................................................................................................................

Introducción

La falla catastrófica de estructuras es causada por la propagación de fracturas más allá de un tamaño seguro. Estas fracturas, presentes de alguna forma en todo componente estructural, son el resultado de defectos de fabricación o de daños localizados generados durante su servicio y pueden crecer lentamente debido a procesos tales como creep o fatiga, entre otros (Broek, 1986), hasta provocar una disminución de su resistencia estructural y generar en el futuro una falla catastrófica.

La determinación de la tolerancia al daño, propiedad estructural que determina la capacidad para sobrellevar de manera segura una fractura durante la vida útil de una estructura, es uno de los objetivos en el análisis y evaluación de fallas. Adicionalmente, la determinación de la velocidad de cre-cimiento y del tamaño crítico, son parámetros que determinan la estabilidad de la fractura y la vida residual de un componente estructural.

En mecánica de la fractura, la determinación de los factores de intensidad de tensiones (stress intensity factors –SIF–) son parámetros con los cuales puede establecerse la tolerancia al daño. La gran dificultad, desde el punto de vista de la mecánica de sólidos, es la ausencia de soluciones analíticas adecuadas que describan la respuesta de cuerpos con geometrías y condiciones de cargas complejas y materiales no lineales (plasticidad, creep, etc.) presentes en aplicaciones reales. Factores de intensificación de tensiones han sido encontrados utilizando las hipótesis de la teoría LEFM para problemas con geometrías y cargas simples (S. Aoki, N. Hasebe, 1987), las cuales son de utilidad limitada en ingeniería.

Si bien la mecánica de la fractura ofrece bases teóricas adecuadas para el estudio de su propagación y tolerancia al daño, su aplicación práctica en ingeniería solo es posible a través del uso de he-rramientas computacionales que permitan obtener soluciones numéricas a las ecuaciones de gobierno que rigen un problema, obteniendo así soluciones más realistas.
El trabajo presenta una formulación matemática basada en el método de los elementos de contorno, que permite modelar cuerpos bidimensionales fracturados en estado de tensión o deformación plana y estudiar el fenómeno de propagación basado en la teoría LEFM. El desarrollo teórico está basado en el trabajo de Portela (1992). La formulación ha sido implementada en un programa computacional. Se presentan tres ejemplos de aplicación en mecánica de la fractura que demuestran la efectividad del programa para predecir trayectorias de propagación, permitiendo su uso como herramienta de análisis para evaluación de daño en ingeniería.

Aspectos teóricos

A diferencia de los métodos de aproximación en el dominio, el de los elementos de contorno (BEM) trabaja con la discretización del contorno donde son definidas ecuaciones integrales de gobierno del problema (Gráfico 1). En comparación con el método de los elementos finitos (FEM), en el BEM la dimensionalidad del problema es reducida, generando de este modo sistemas de ecuaciones de menor tamaño. Igualmente, BEM es superior al FEM en el tratamiento de problemas en mecánica de la fractura, ya que no se requieren mallas altamente refinadas (y algoritmos especiales para gene-rarlo y actualizarlo si hay dificultades de propagación) alrededor de la punta de la fractura para capturar los altos gradientes de tensiones, tal como ocurre en el FEM.


Formulación integral de contorno

La formulación por elementos de contorno para problemas en elasticidad lineal parte de considerar la formulación por residuos ponderados de las ecuaciones de Navier (Brebbia y Domínguez, 1989):

donde son las componentes de desplazamiento, son las fuerzas de cuerpo, v es el módulo de Poisson y es el modulo cortante. Utilizando integración por partes se obtiene la ecuación integral:

La anterior ecuación corresponde al teorema de reciprocidad de Betti en mecánica de sólidos. En esta ecuación, son tensiones aplicadas en el contorno, son funciones de ponderación y representa el divergente del tensor de tensiones ponderado. La primera integral del lado izquierdo de esta ecuación puede ser transformada en una integral de contorno si se considera la solución fundamental de la ecuación 1, la cual es obtenida considerando una carga puntual en un dominio infinito:

De esta manera, y teniendo en cuenta la propiedad fundamental de la función Delta de Dirac: se obtiene una ecuación integral de contorno:

donde es una constante que depende de la ubicación del punto de colocación para puntos en el dominio, para puntos sobre el contorno y para puntos por fuera del contorno) son las soluciones fundamentales para tracción y desplazamiento, respectivamente, para la ecuación 1, las cuales han sido tomadas como funciones ponderadoras en la ecuación 2. La anterior ecuación es la base del método directo de los elementos de contorno. Existen técnicas especiales para el tratamiento de la integral de dominio de las fuerzas de cuerpo utilizadas para hallar una integral equivalente evaluadas en el contorno, pero no serán tratadas aquí (ver Patridge, Brebbia y Wrobel, 1992).

Formulación por elementos de contorno

Ya que en general no existe una solución analítica para la ecuación 4, esta puede hallarse a través del tratamiento numérico de las integrales con la discretización geométrica del contorno, utilizando elementos discretos y aproximando sobre ellos el campo de desplazamientos y tensiones. Varios tipos de elementos pueden ser utilizados, siendo más comunes los lineales y cuadráticos continuos o discontinuos (ver Brebbia y Domínguez, 1989).

Si se utilizan como puntos de colocación los mismos que definen los elementos en el contorno y utilizando la siguiente aproximación para los campos de desplazamientos y tensiones dentro de un elemento cualquiera se tiene:

donde φi ,i=1...3 son las funciones de forma de un elemento de borde cuadrático continuo o discontinuo. Se obtiene un sistema de ecuaciones linealmente independiente, el cual puede ser escrito como:

donde N es el número de nodos en el contorno, es el vector de desplazamientos en el nodo i, es el vector de tracciones aplicadas en el contorno en el nodo j, y H, G y B son matrices cuyos términos contienen las soluciones fundamentales de la ecuación 1. Una vez obtenidos desplazamientos y tensiones en el contorno, pueden utilizarse puntos en el dominio para encontrar soluciones internas: desplazamientos, tensiones y deformaciones. De forma compacta y considerando todas las contribuciones para un nodo dado se tiene finalmente un sistema de ecuaciones del tipo:

Modelado de fractura

El anterior procedimiento no puede ser aplicado directamente a problemas de fractura, ya que sobre el contorno de esta existirán elementos superpuestos (nodos coincidentes), lo cual generaría un sistema de ecuaciones linealmente dependientes (sistema singular). Por tanto, para analizar problemas de mecánica de la fractura mediante BEM, es utilizada una formulación alterna donde son escritas adicionalmente a las ecuaciones de desplazamientos, ecuaciones de tracción para ciertos nodos en el contorno. Dicha formulación es conocida como método dual de elementos de contorno (DBEM). En esta formulación, son escritas ecuaciones de tracción para los nodos sobre uno de los bordes de la fractura y ecuaciones de desplazamientos para los nodos en el borde opuesto. De esta manera se obtiene un sistema linealmente independiente de ecuaciones, lo cual asegura la solución del sistema de ecuaciones resultante.

Derivando la ecuación 4 con respecto a las coor-denadas espaciales y teniendo en cuenta que:

donde es el vector de tracción, es el tensor de tensiones de Cauchy, es el tensor de constantes elásticas y es el tensor de pequeñas deformaciones, se obtiene la ecuación integral de tracción (Portela, 1992):

donde es el vector normal del elemento en el punto de colocación y los kernels y son encontrados derivando las soluciones fundamentales de las ecuaciones de Navier.

Propagación de fractura

La mecánica de la fractura provee criterios que permiten determinar sus tensiones críticas (SIFs críticos), dirección de propagación y dimensión crítica. Varios criterios han sido propuestos para determinar la dirección de propagación de una fractura, entre ellos el más utilizado es el de tensión principal máxima, el cual establece que la propagación ocurrirá en la dirección perpendicular a la dirección de tensión máxima calculada en la punta de la fractura. Considerando la definición de los factores de intensidad de tensiones este criterio puede ser escrito como (Meguid, 1989):

donde es la dirección de propagación. En este trabajo, estos factores son calculados mediante extrapolación del campo de desplazamiento en los nodos próximos a la punta de la fractura:

donde son los desplaza-mientos de los nodos D, E, F y G de los elementos que definen la punta de la fractura (Aliabadi y Rooke, 1992).

Por otra parte, la extensión de la fractura ocurrirá cuando el valor principal de tensión alcance el valor que en modo I generará una propagación de la fractura:

donde representa la resistencia a la fractura en estado plano para el material. En el presente trabajo el crecimiento de fractura por fatiga es modelado a través de la ley de Paris (ver Broek, 1986), la cual está dada por:

donde los parametros y N representan el tamaño y el número de ciclos de carga respectivamente, y C y m son constantes empíricas características del material, determinadas a través de ensayos de laboratorio para fractura en modo I. En este modelo, el rango de variación del factor de intensificación de tensiones puede calcularse como:

en esta ecuación R representa la amplitud de la variación de la tensión aplicada . Otros modelos empíricos han sido propuestos los cuales pueden ser consultados en Broek (1986).

Implementación computacional

La formulación por elementos de contorno desarrollada en la sección anterior fue implementada en un código computacional utilizando Matlab 7.0 , con sistema operacional Windows XP . El programa trabaja con elementos lineales cuadráticos discontinuos para la aproximación del campo de desplazamientos y de tensiones en los elementos. Utiliza funciones de formas cuadráticas continuas para describir la geometría de los elementos. En el tratamiento numérico de las integrales de contorno se utilizó cuadratura de Gauss de seis puntos para integrales no singulares y cuadratura de Gauss logarítmica con diez puntos para integrales “débiles” (términos en la diagonal de la matriz G) (Dirgantara, 2002). Para las integrales “fuertes” que se presentan en los términos sobre la diagonal de la matriz H se utilizó una expansión asintótica por series de Taylor, e integración gaussiana de seis puntos. Se utiliza un esquema de actualización de las matrices H y G durante cada incremento de la fractura obteniéndose una disminución considerable en el costo computacional al reutilizar estas matrices (Portela, 1992).

Aplicaciones

Placa con fractura central

Se presenta el análisis de una placa rectangular con una fractura central, como verificación del programa para el cálculo de factores de intensificación de tensiones. El problema se muestra en el Gráfico 2A. Los parámetros geométricos consi-derados son: y espesor igual a 0.01 m. Es aplicada una tracción normal en los extremos de la placa tal como se muestra. Las constantes elásticas para la placa son:

En el Gráfico 2B se ve el modelo por elementos de contorno propuesto. Los contornos de la placa han sido discretizados utilizando 30 elementos cuadráticos discontinuos y 8 elementos en cada borde de la fractura. El Gráfico 2C muestra la geometría deformada de la malla de elementos de contorno. En el Cuadro 1 está la comparación de los SIFs encontrados utilizando BEM y la solución analítica (Aoke y Hasebe, 1987). Se observa una buena correlación con los resultados analíticos con errores máximos de 0.95 %.

Viga con fractura de canto

Se modela el proceso de crecimiento de una fractura de canto en una viga fabricada en acrílico utilizada en el estudio experimental de propagación de fractura plana. El Gráfico 3A muestra las características de la viga y el modelo por elementos de contornos. Se emplean 26 elementos cuadráticos discontinuos para modelar el contorno de la placa y 8 en los bordes de la fractura. El contorno de los agujeros es discretizado con ocho elementos. Se ha utilizado un incremento de extensión de fractura igual a 0.0762 m. Los resultados pueden verse en el Gráfico 3B para el quinto incremento de extensión y para el incremento final. Se observa una buena correlación con los resultados dados por Bittencourt, Wawrzynek e Ingraffea (1996) no mostrados aquí.

Fractura en placa con agujeros

Este ejemplo muestra el análisis de crecimiento de una fractura de borde de una placa con agujeros, tal como se presenta en el Gráfico 4A. La placa tiene las siguientes dimensiones: . Se ha consi-derado la longitud inicial de la fractura igual a 0.1m. Una carga uniforme de tracción es aplicada en los extremos de la placa y de manera perpendicular al eje inicial de la fractura, tal como se muestra. Se considera la placa con propiedades: E = 200000 y v = 0.25. Un proceso de propagación de fractura por fatiga ha sido modelado considerando una amplitud constante de carga con R = 2/3. Las constantes utilizadas en la ley de Paris son: C = 4.624x10-12 y m = 3.3.

El contorno de la placa ha sido discretizado inicialmente con 30 elementos cuadráticos discontinuos en el borde externo y 8 en el borde que define los agujeros. Inicialmente tres elementos modelan la fractura. Se utilizan incrementos en el tamaño de la fractura iguales a 0.06667 m. La serie de gráficos mostrados en el Gráfico 4B presentan el proceso de propagación de la fractura, donde inicialmente se observa un avance de esta hacia el agujero inferior generada por el campo de tensiones en la cercanía a este. Posteriormente la fractura se propaga casi sobre la línea media de la placa, lo cual podría esperarse si se considera que el campo de tensiones en aquella zona es casi uniforme y genera una propagación básicamente en modo I. Estos resultados presentan una alta correlación con los reportados por Portela (1992).

El diagrama del Gráfico 5A muestra la variación de los factores de intensificación de tensiones en los modos como función del número de incremento de extensión de la fractura. Hay una disminución del factor en modo I con el aumento de la longitud de la fractura, mientras que el factor en modo II permanece casi constante y con valor muy pequeño, debido básicamente a que es un problema de propagación de fractura en modo I.

Por otra parte, el Gráfico 5B, muestra la variación de la resistencia residual en función del número de ciclos de carga N calculados mediante la ecuación 13. Para esto es utilizada integración numérica basada en el método trapezoidal. Este diagrama de resistencia muestra la variación de la máxima carga que el componente estructural puede soportar sin generar una propagación inestable de la fractura a medida que la longitud de fractura aumenta. En este grafico, el valor de resistencia es representado mediante la relación:

donde representan los valores de resistencia residual y el factor equivalente de intesificación de tensiones en modo I, calculados en la configuración inicial de la fractura. Específicimente, la resistencia residual está dada por la expresión (Portela, 1992):

donde Y es un factor geométrico y es un valor de referencia para el esfuerzo aplicado.

Conclusiones

Fue presentada una formulación basada en el método de elementos de contorno para el modelado numérico de la propagación de fractura por fatiga, implementada en un código computacional desarrollado en Matlab 7.0 . Los resultados de los ejemplos analizados demuestran que los factores de intensificación de tensiones, las trayectorias de propagación y resistencia residual calculados, presentan una buena correlación con aquellos obtenidos experimentalmente y con los reportados en la literatura especializada. Lo anterior evidencia la versatilidad y confiabilidad del método de los elementos de contorno en el tratamiento de problemas en mecánica de la fractura. La herramienta computacional desarrollada permite el análisis de propagación de fractura por fatiga, estimación de trayectorias de propagación y resistencia residual de componentes fracturados, lo cual la convierte en un instrumento valioso para la toma decisiones en labores de inspección de equipos, siendo posible su implementación en sistemas de vigilancia y control de fracturas en tiempo real.

Referencias

Aliabadi,M. H.; Rooke, D. P. (1992) Numerical Fracture Mechanics, Southampton, Computa-tional Mechanics Publications, Elsevier.
Aoki, S.; Hasebe, N. (1987) Stress Intensity Factors Handbook, New York, Pergamon Press.

Brebbia, C. A.; Domínguez, J. (1989) Boundary Elements an Introductory Course, New York, McGraw-Hill.

Bittencourt, T. N., Wawrzynek, P. A., Ingraffea, A. R. (1996) “Quasi-Automatic Simulation of Crack Propagation for 2D LEFM problems”, Engineering Fracture Mechanics, vol. 55, n.o 2, pp. 321-334.

Broek, D. (1986) Elementary Engineering Fracture Mechanics, The Netherlands, Martinus Nijhoff Publishers, Dordrecht.

Dirgantara, T. (2002) Boundary Element Analysis of Cracks in Shear Deformable Plates and Shells, Southampton, WIT press.

Meguid, S. A. (1989) Engineering Fracture Mecha-nics, New York, Elsevier.

Patridge, P. W.; Brebbia, C. A.; Wrobel, L. C. (1992), The dual reciprocity boundary element method, Southampton, Computational mecha-nics publications, Elsevier.

Portela, A. (1992) Dual Boundary Element Incremental Analisis of Crack Growth [tesis doctoral], Southampton, Wessex Institute of Techonology.

Portela, A.; Aliabadi, M. A. (1992) “The Dual Boundary Element Method: Effective Implementation for Crack Problems”, International Journal for Numerical Methods in Engineering, vol. 33, pp. 1269-1297.