El método de los elementos de contorno axisimétrico aplicado al análisis de grietas de interfase en materiales compuestos reforzados con fibras

EL MÉTODO DE LOS ELEMENTOS DE CONTORNO AXISIMÉTRICO APLICADO AL ANÁLISIS DE GRIETAS DE INTERFASE EN MATERIALES COMPUESTOS REFORZADOS CON FIBRAS E. Graciani 1 , V. Manti 1 , F. París 1 , J. Varna 2 1 Grupo de Elasticidad y Resistencia de Materiales, Escuela Superior de Ingenieros, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092-Sevilla, España. E-mail: [email protected], [email protected], [email protected] Tfno: 954 48 73 00. Fax: 954 46 16 37 2 Division of Polymer Engineering, Luleå University of Technology, SE-97187-Luleå, Sweden. [email protected] RESUMEN En el presente artículo se muestra un análisis de los diferentes estados elásticos que aparecen durante la realización del ensayo de fragmentación de fibra única, empleando el Método de los Elementos de Contorno. El objetivo fundamental del trabajo es el análisis de la aparición y propagación de una grieta de despegue entre la fibra y la matriz, y de su influencia en las tensiones axiales a lo largo de la fibra. Se han empleado dos enfoques diferentes en el análisis: el enfoque de grieta abierta (el cual, aunque no es realista en este caso, ofrece algunos resultados aceptables) y el enfoque de grieta cerrada (considerando o no el efecto de la fricción entre las caras de la grieta). Se analiza, además, la solución asintótica de las tensiones y los desplazamiento en el entorno de dos grietas diferentes: primero, una grieta que parte la fibra y, segundo, una grieta de despegue. ABSTRACT An analysis of the micromechanical elastic fields arising in the single fibre fragmentation test, carried out employing the Boundary Element Method, is presented in this paper. The main objective of the work is to study the initiation and growth of a debond crack along the fibre-matrix interface, and its influence on the axial stresses along the fibre. Two different models have been considered: the open crack model (which, although is clearly unrealistic in this particular case, give some acceptable results) and the closed crack model (excluding and including the effect of friction between crack faces). The asymptotic behaviour of stresses and displacements in the vicinity of the crack tips are studied, first, for a crack splitting the fibre and, second, for a debond crack. ÁREAS TEMÁTICAS PROPUESTAS: Fractura de Materiales Compuestos. PALABRAS CLAVE: Método de los Elementos de Contorno, Materiales Compuestos, Grietas de Interfase. 1. INTRODUCCIÓN El ensayo de fragmentación de fibra única [1], es una técnica experimental ampliamente usada en la caracterización de la interfase fibra-matriz en materiales compuestos. Para ello se emplea una probeta con una única fibra rodeada de matriz polimérica sometida a tracción. Como resultado de esta tracción, la fibra se fragmenta sucesivamente y aparecen grietas de despegue que nacen del final de las grietas de la fibra y crecen, de forma estable, a lo largo de la interfase. Dado que, por un lado, la diferencia entre los coeficientes de expansión térmica entre la fibra y la matriz induce una compresión radial de la fibra durante el enfriamiento desde la temperatura de solidificación hasta la temperatura ambiente, y, por otro, la diferencia entre sus coeficientes de Poisson hace que la compresión radial de la fibra aumente al aplicar la carga, la grieta de despegue se encuentra cerrada en toda su longitud, creciendo, por tanto en modo II puro. Debido a la gran diferencia entre los tamaños relativos del diámetro de la probeta, el diámetro de la fibra y el tamaño de las grietas de despegue, se ha considerado que el Método de los Elementos de Contorno (MEC) es la herramienta más adecuada para llevar a cabo el análisis numérico del problema. Dado que existe simetría de revolución respecto al eje de la fibra, el estudio se ha llevado a cabo mediante un algoritmo que implementa el MEC para el problema de contacto elástico con fricción entre sólidos con simetría axial. Con el objeto de poder emplear discretizaciones no conformes en las interfases y zonas de contacto, las ecuaciones de equilibrio y compatibilidad se imponen mediante una formulación débil [2]. Para llevar a cabo el análisis se han tenido en cuenta los dos enfoques más extendidos en el análisis de grietas de interfase [3]: el enfoque de Williams o de grieta abierta [4,5], que supone que las caras de la grieta están libres Anales de Mecánica de la Fractura Vol. II (2006) 559 de tensiones, y el enfoque de Comninou o de grieta cerrada [6,7], que asume que en el entorno del vértice existe una zona de contacto. Aunque, a priori, en el presente caso, dado que la grieta está cerrada en toda su longitud, el enfoque de grieta abierta no es claramente el más adecuado, se han incluido en el presente trabajo alguno de los resultados obtenidos con el mismo con el objeto de comprobar su validez. 2. DESCRIPCIÓN DEL MODELO Considerando que los fragmentos en que se divide la fibra son suficientemente largos como para que lejos del plano de rotura de la fibra la solución esté regularizada, que en cada fragmento existe simetría de la solución respecto a dicho plano, y que el estado tensional en la probeta tiene simetría axial respecto al eje de la fibra, es posible reducir el análisis numérico del ensayo a la sección radial del tramo de probeta correspondiente a la mitad de un fragmento de la fibra. Con lo cual la geometría y las condiciones de contorno serán las que se muestran en la figura 1. Para el caso considerado, en que la probeta consiste en una fibra de vidrio de radio r f .=.5.m rodeada por una matriz cilíndrica de radio r m .=.55.r f , se ha comprobado que para longitudes L f . .40.r f la solución en el entorno de las fisuras es idéntica en todos los casos, lo cual valida la hipótesis de solución repetitiva en el extremo de todos los fragmentos. El tamaño de la grieta de despegue se ha variado en el rango 0. .a/r f . .30. Las propiedades termoelásticas empleadas para la fibra han sido: E f .=.70.GPa, f .=.0.2 y f .=.7.K -1 ; y las de la matriz: E m .=.3.5.GPa, m .=.0.3 y m .=.50.K -1 . La temperatura de solidificación durante el se ha supuesto a 80.K por encima de la temperatura ambiente. fo Lu r z f L f r m r pegado o T o T fo Lu r z f L f r m r pegado r z f L f r m r pegado o T o T Figura. 1. Condiciones de contorno. Las condiciones de contorno empleadas son las que se esquematizan en la figura 1, es decir, simetría respecto al plano de la rotura de la fibra, alargamiento uniforme del plano medio del fragmento, fibra libre de tensiones en su plano de rotura y contacto entre las caras de la grieta de despegue (salvo en los casos correspondientes al enfoque de Williams en que se ha considerado que están libres de tensiones). Dado que se pretende calcular la evolución de las tensiones axiales a lo largo del eje de la fibra, se ha empleado en el contorno de la misma una longitud máxima de los elementos igual a r f. /3.=.1.7.m. Dado que en la matriz no es preciso refinar tanto, se han empleado discretizaciones no conformes en la interfase, tal y como se observa en la figura 2, donde se muestra una discretización típica, en la cual, en la zona de interfase no dañada hay 216 elementos del lado de la fibra y 91 del lado de la matriz. El empleo del MEC hace que sea muy sencillo obtener una discretización extremadamente fina en el entorno del vértice de la fisura (siendo la longitud del menor elemento empleado de 10 -10 .m). 0 50 100 150 200 250 050100150200250 Matriz Fibra 0 10 20 0102030 Figura. 2. Discretización. 3. FRAGMENTACIÓN DE LA FIBRA Durante la primera etapa del ensayo de fragmentación, el único efecto que puede observarse es la rotura de la fibra en fragmentos cada vez más pequeños. Esta rotura se produce por los puntos más débiles de la misma y progresa instantáneamente hasta alcanzar la interfase fibra-matriz. De acuerdo a [8], el análisis de la solución asintótica en una esquina multimaterial en deformación plana desarrollado en [9] puede aplicarse en el caso análogo de simetría axial. Con lo cual, se van a usar esos resultados semi-analíticos para comprobar la exactitud de la solución obtenida mediante el MEC en el entorno del vértice de la grieta existente en el plano z.=.0 del problema de la figura 1 (cuando a.=.0). De esta forma, el comportamiento asintótico esperado cuando la distancia al vértice, , disminuye ( . .0), es u i .=.O( ) para los desplazamientos y ij .=.O( -1 ) para las tensiones, siendo =.0.1869. En la figura 3 se muestran, conjuntamente, los resultados obtenidos con el MEC en el entorno del vértice de la grieta en la fibra para los desplazamientos (relativos al vértice) a lo largo de la cara de la grieta y las compo- Anales de Mecánica de la Fractura Vol. II (2006) 560 nentes singulares del vector tensión en el fondo de grieta y la interfase. Como puede observarse, el comportamiento asintótico cuando . .0 se ajusta con excelente precisión a una recta en coordenadas log-log. Los resultados mostrados corresponden al caso de L f .=.100.r f y una deformación aplicada a la probeta de valor .=.2%. 1 E-2 1 E-1 1 E-5 1 E-4 1 E-3 1 E-2 Distancia al vértice [ m] 1 E+2 1 E+3 1 E+4 1 E+5 1 E+6 1 E+7 rr rz zz u z u r Figura. 3. Comportamiento asintótico de la solución en el entorno del vértice de la grieta en la fibra. En la tabla 1 se recoge el valor de obtenido a partir de cada una de las variables mostradas en la figura 3, empleando un ajuste de mínimos cuadrados de los resultados numéricos en el intervalo 10 -5 .m.<. .<.10 -4 .m. Como puede observarse, la precisión obtenida es excelente en todos los casos. [9] u r u z rr rz zz .1869 .1889 .1868 .1928 .1870 .1869 Err (%) 1.07 0.05 3.16 0.05 0.00 Tabla 1. Orden de singularidad de la solución numérica. 4. DESPEGUE ENTRE FIBRA Y MATRIZ Una vez alcanzado un cierto nivel de carga en la probeta, pueden surgir grietas de despegue que nacen del vértice de las grietas de la fibra y crecen a lo largo de la interfase. A continuación se va a estudiar la solución elástica del ensayo durante el proceso de crecimiento de dichas grietas, despreciando el efecto de la fricción (el cual se incluirá en el apartado siguiente). El análisis se va a centrar en el estudio de la solución asintótica en el entorno del vértice de la grieta de despegue, en el efecto que tienen dichas grietas sobre la tensión axial en la fibra y, por último, en su propagación. 4.1. Solución asintótica En la figura 4 se muestra (para el caso particular en que L f .=.100.r f , a.=.5.r f y .=.2%) la solución obtenida para las componentes del vector tensión a lo largo de la interfase, en el entorno del vértice de la grieta, empleando los dos modelos descritos. Como puede observarse, en el modelo de contacto la grieta está cerrada en toda su longitud, siendo singulares la presión de contacto, rr,C , y la tensión tangencial, rz,C+ , en el fondo de grieta, por contra, la tensión normal, rr,C+ , es acotada, tal como predice la teoría [6]. -400 -300 -200 -100 0 20 25 30 z [ m] vértice rz ,C+ rz ,A+ rr ,C+ rr ,A+ rr ,C Figura. 4. Solución en el entorno del vértice. Aunque en la figura 4 parece que las diferencias entre el modelo abierto (A) y el modelo de contacto (C) se producen fundamentalmente en las tensiones normales, en la figura 5, puede comprobarse que el carácter asintótico de la solución obtenida en el fondo de grieta es diferente en ambos modelos, presentando las tensiones normales, rr,A , y tangenciales, rz,A+ , del modelo abier- to una singularidad oscilante, tal y como predice la teoría [4,5]. -200 -150 -100 -50 0 50 1 E-7 1 E-6 1 E-5 1 E-4 1 E-3 1 E-2 1 E-1 Distancia al vértice de la grieta, [ m] rz ,C+ rz ,A+ rr ,C+ rr ,A+ Figura. 5. Comportamiento asintótico de la solución en el entorno del vértice de la grieta de despegue. Puede observarse que las tensiones rr,A+ (del modelo abierto) son de compresión en el entorno del vértice de la grieta, incluso para distancias tan próximas que dejan de tener sentido físico, con lo cual, sea cual sea la escala que se elija, definirán un factor de intensificación de tensiones K I .<.0 y, en consecuencia, no se cumplen las condiciones de validez de dicho modelo [5]. Por ello, en lo que sigue, salvo que se indique lo contrario, se mostrarán únicamente resultados del modelo de contacto. 4.2. Tensiones axiales en la fibra Para ver la influencia de la aparición de las grietas de despegue en la tensión axial en la fibra, se muestran en la figura 6 los resultados obtenidos para zz a lo largo del eje de la fibra en una serie de problemas en los cuales L f .=.100.r f , 0.<.a.<.25.r f y .=.2%). Anales de Mecánica de la Fractura Vol. II (2006) 561 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 100 200 300 400 z [ m] a= 25 rf a= 20 rf a= 15 rf a= 10 rf a= 5 rf a= r f a= 0 Figura. 6. Tensiones axiales en la fibra. Como puede observarse, cuanto mayor es la longitud de la grieta de despegue, menos eficiente es la transferencia de carga de la matriz a la fibra y, por tanto, mayor es la zona de fibra que se encuentra descargada. En consecuencia, al incrementar la carga aplicada a la probeta, cuanto mayor sea el tamaño de las grietas de despegue, menor será el volumen de fibra sometido a altas tensiones y, por tanto, susceptible a una nueva rotura, con lo cual, a interfases débiles les corresponderán, una vez terminado el ensayo, longitudes de fragmento mayores que si la interfase es más resistente. Este efecto se pone de manifiesto en la figura 7, en la que se muestra la cota z 95% a partir de la cual se recupera el 95% de la tensión previa a la rotura de la fibra. Como puede observarse, para tamaños pequeños de la grieta de despegue la cota z 95% aumenta rápidamente, tendiendo para tamaños mayores a una relación z 95% (a) lineal, bastante más suave. 0 20 40 60 0 10 20 30 Tamaño de la grieta de despegue, a/rf Figura. 7. Efecto de la grieta de despegue en el tamaño de la zona de exclusión. 4.2. Propagación de la grieta de despegue Para estudiar la propagación de la grieta de interfase se ha calculado el índice de liberación de energía, G, empleando una generalización de la técnica del cierre virtual de grieta propuesto por Irwin al caso de grietas de interfase. Para ello se ha efectuado una integración numérica, usando una cuadratura de Gauss-Jacobi con una densidad apropiada y ocho puntos de integración, en un intervalo de tamaño a.=.10 -3 .m del entorno del vértice de la fisura. En la figura 8 se muestran los resultados obtenidos con el modelo de contacto, G T , C .=.G II , C , correspondiente a un crecimiento en modo II puro, junto con los obtenidos a partir de la solución del modelo abierto, con el cual se obtienen unos resultados globales, G T , A , muy similares a los del modelo de contacto, a pesar de que la división entre modo I, G I , A , y modo II, G II , C , es en realidad ficticia y dependiente de la longitud de integración a. Los resultados mostrados corresponden al caso particular en que L f .=.100.r f , no obstante, para longitudes de fragmento L f . .40.r f se obtienen valores de G prácticamente idénticos a los mostrados, lo cual valida la suposición de que no existía interacción entre las grietas de despegue de ambos extremos del fragmento. 1 E+0 1 E+1 1 E+2 1 E+3 1 E-2 1 E-1 1 E+0 1 E+1 1 E+2 Tamaño de la grieta de despegue, a [ m] GI,A GII ,A GT,A= GI,A+ GII ,A GT,C = G II ,C Figura. 8. Variación del índice de liberación de energía con el tamaño de la grieta de despegue. Como puede verse en la figura, G no está acotada cuando a. .0, lo cual está de acuerdo con el comportamien- to asintótico G.=.O( -1 ).=.O( 0.626 ) obtenido generalizando el análisis de [10]. Ajustando los resultados numéricos obtenidos con el modelo de contacto para las tres longitudes de grieta más pequeñas se obtiene G T , C .= O( 0.609 ) lo cual da una idea de la excelente precisión de los resultados obtenidos. El hecho de que G no esté acotada cuando a. .0 implica que la grieta de despegue aparecerá de forma inmediata cuando la grieta de la fibra alcance la interfase. No obstante, dado que, para valores pequeños de a, G decrece rápidamente cuando progresa la fisura, el crecimiento inicial será estable. Si se emplea el modelo de contacto, con un criterio del tipo G II .=.G C , para estudiar la propagación de la grieta de despegue se obtienen los resultados mostrados en la figura 9 para dos valores característicos de la tenacidad a fractura de la interfase. Como puede observarse, la propagación es muy lenta al principio, acelerándose progresivamente hasta llegar a una propagación inestable, que se corresponde con la asíntota horizontal de la curva G(a), cuando a. . , que se observa en la figura 8. Anales de Mecánica de la Fractura Vol. II (2006) 562 0 40 80 120 160 01234 [%] 2 2 J/m 100 J/m 50 C C G G Figura. 9. Propagación de la grieta de interfase. Ajustando los resultados experimentales de a( ) con una curva como la mostrada en la figura 9 permitiría la obtención de valores cuantitativos de la tenacidad a fractura de la interfase. No obstante, para ello, sería preciso incluir el efecto de la fricción en la propagación de la grieta. 5. EFECTO DE LA FRICCIÓN Para incluir el efecto de la fricción en el análisis es obvio que debe emplearse el modelo de contacto, ya que es imposible incluirlo en el modelo abierto que supone que las caras de la grieta están libres de tensiones. A continuación se presentan los resultados obtenidos para diferentes valores del coeficiente de fricción, , entre las caras de la grieta de despegue. La solución ha sido obtenida con las condiciones de contorno indicadas en la figura 1 y en todos los casos se ha obtenido como resultado que la grieta está cerrada y deslizando en toda su longitud, con lo cual la solución del problema (para un tamaño fijo de la grieta de despegue) es lineal con la deformación aplicada. 5.1. Solución asintótica Comparando los resultados del modelo de contacto con y sin fricción se obtienen los resultados que se muestran en la figura 10 (para el caso particular de L f .=.100.r f , a.=.5.r f , .=.2% y .=.0.5). -300 -200 -100 0 20 25 30 z [ m] vértice rr ,C ( =0.0) rr ,C ( =0.5) rz ,C ( =0.5) rr ,C+ ( =0.0) rr ,C+ ( =0.5) rz ,C+ ( =0.0) rz ,C+ ( =0.5) Figura. 10. Solución en el entorno del vértice. Como puede observarse, debido a la fricción aparece una tensión tangencial singular, rz,C .=. rr,C , en la zona despegada y, por lo demás, el resto de tensiones en el entorno del vértice disminuyen. Además, el orden de singularidad de la solución se desvía del clásico .=.0.5, obteniéndose un valor que depende de según la ecuación tg( ).=.( ) 1 , donde es el segundo parámetro de Dundurs [7]. Como muestra de la precisión con la que los resultados obtenidos cumplen el comportamiento asintótico predicho, en la figura 11 se recogen, junto con la solución analítica, la solución numérica obtenida para el orden de singularidad y el error porcentual cometido su cálculo (efectuado a partir de u z y rz de forma análoga a como se explicó en el caso de la grieta en la fibra), en una serie de problemas en que se ha variado el coeficiente de fricción en el intervalo 0. . . .1. Se observa que el error cometido es prácticamente nulo en el cálculo de a partir de u z y del 1% como máximo en el cálculo a partir de rz . 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 Coeficiente de fricción -1 0 1 ( rz ) ( uz) Err (rz ) Err (uz) Figura. 11. Influencia de la fricción en el orden de singularidad. Solución numérica, , y error cometido, Err . 5.2. Tensiones axiales en la fibra El efecto que tiene la fricción en las tensiones axiales en la fibra puede observarse en las figuras 12 y 13, donde se muestra la solución obtenida para zz en el eje de la fibra en dos series de problemas. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 100 200 300 400 z [ m] = a=25 rf, = 1.0 a=25 rf, = 0.8 a=25 rf, = 0.6 a=25 rf, = 0.4 a=25 rf, = 0.2 a=25 rf, = 0 Figura. 12. Efecto de en zz (0,z): a.=.25r f , 0....1. En primer lugar, figura 12, se ha empleado un tamaño constante de la grieta, a.=.25r f , y se ha variado el coefi- Anales de Mecánica de la Fractura Vol. II (2006) 563 ciente de fricción, observándose cómo la presencia de la fricción hace que el aumento de la tensión axial de la fibra sea más paulatino, presentando una zona inicial en que las tensiones aumentan linealmente con una pendiente proporcional al coeficiente de fricción. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 100 200 300 400 z [ m] a= 0 a= 30 rf, =0.5 a= 20 rf, =0.5 a= 10 rf, =0.5 a= 30 rf, =0 a= 20 rf, =0 a= 10 rf, =0 Figura. 13. Efecto de en zz (0,z): .=..5, 0. .a. .30r f . En el segundo caso, figura 13, se ha mantenido constante el coeficiente de fricción, .=.0.5, y se ha variado el tamaño de la grieta, observándose, que (una vez superada la zona lineal) la diferencia entre la solución con y sin fricción es pequeña si el tamaño de la grieta de despegue es pequeño y aumenta progresivamente con dicho tamaño. Aunque a la vista de los resultados anteriores es claro que la fricción eleva el estado tensional de la fibra en la zona de transferencia, no debe tener, sin embargo, una gran influencia en el proceso de fragmentación de la misma, dado que la cota a la cual se recupera un nivel de carga suficiente como para que se pueda volver a producir la rotura de la fibra en el interior del fragmento no depende significativamente del coeficiente de fricción, como se observa en la figura 10, en la que se muestra claramente que la influencia de en la cota z 95% a partir de la cual se recupera el 95% de la tensión previa a la rotura de la fibra es pequeña durante la etapa en que comienzan a progresar las grietas de despegue. 0 20 40 60 0 10 20 30 Tamaño de la grieta de despegue, a/rf =0 =0.5 Figura. 10. Efecto de la fricción en el tamaño de la zona de exclusión. 5.3. Propagación Dado que, debido a la presencia de la fricción, el deslizamiento relativo de las caras de la grieta introduce una disipación de energía que es imposible de recuperar, si se quiere incluir este efecto, el análisis de la propagación de la fisura debe de hacerse siguiendo paso a paso el crecimiento de la grieta de despegue, lo cual conlleva un complejo proceso de remallado, es por ello que no se ha incluido dicho análisis en el presente estudio. CONCLUSIONES Se ha desarrollado una herramienta numérica basada en el MEC para el análisis tensional de problemas con grietas de interfase en materiales compuestos que permite calcular con excelente precisión la solución asintótica en entornos extremadamente próximos al vértice, empleando discretizaciones sencillas y teniendo en cuenta el posible contacto con fricción entre sus caras. Del análisis del ensayo de fragmentación se han obtenido, por un lado, una serie de curvas que permiten relacionar cualitativamente la tenacidad a fractura de la interfase con la longitud media de los fragmentos resultante tras el ensayo y, por otro, una serie de curvas que se espera sirvan para la cuantificación de dicha tenacidad mediante ajuste con los resultados experimentales. REFERENCIAS [1] Kelly, A, y Tyson W.R., “Tensile properties of fiber-reinforced metals: copper/tungsten and copper / molybdenum”, Journal of the Mechanics and Physics of Solids, 13: 329-350, 1965. [2] Graciani, E., Manti , V., París, F. y Blázquez, A., “Weak formulation of axi-symmetric frictionless contact problems with boundary elements. Application to interface cracks”, Computers & Structures, 83: 836-855, 2005. [3] Manti , V., Blázquez, A., Correa, E. y París, F., “Analysis of interface cracks with contact in composites by 2D BEM”, Fracture and Damage of Composites, M. Guagliano y M.H. Aliabadi eds., WIT Press, 189-248, 2006. [4] Williams M.L., “The stress around a fault or crack in dissimilar media”, Bulletin of the Seismological Society of America, 49: 199-204, 1959. [5] Rice J.R., “Elastic fracture mechanics concepts for interfacial cracks”, Journal of Applied Mechanics, 55: 98-103, 1988. [6] Comninou M., “The interface crack”, Journal of Applied Mechanics, 44: 631-636, 1977. [7] Comninou M., “Interface crack with friction in the contact zone”, Journal of Applied Mechanics, 44: 780-781, 1977. [8] Zak A.R., “Stresses in the vicinity of boundary discontinuities in bodies of revolution”, Journal of Applied Mechanics, 31: 150-152, 1964. [9] Barroso A., Mantic V. y París F., “Singularity analysis of anisotropic multimaterial corners”, International Journal of Fracture, 119: 1-23, 2003. [10] Leguillon D. y Sánchez-Palencia E., “Fracture in heterogeneous materials - weak and strong singularities”, New Advances in Computational Structural Mechanics, P. Ladeveze y O.C. Zienkiewicz eds., 423-434, 1992. Anales de Mecánica de la Fractura Vol. II (2006) 564