Estudio micromecánico de despegues en la interfase entre fibra y matriz en materiales compuestos

320 ANALES DE l\tiECANICA DE LA FRACTURA Vol. 12 ESTUDIO MICROMECANICO DE DESPEGUES EN LA INTERFASE ENTRE FIBRA Y MATRIZ EN MATERIALES COMPUESTOS J.C. del Cañot, J. Cañas:j:, F. París:j:, V. Mantic:j:, J. Sánchez:j: t E.T.S.I. Industriales; Dpto. Resistencia de Materiales, Pº del Cauce s/n, 47011 Valladolid :j: E.T.S.I. Industriales; Dpto. Mecánica Medios Continuos, Av. Reina Mercedes s/n, 41012 Sevilla Resumen. Se presenta un estudio numérico basado en el Método de los Elementos de Contorno del problema micromecánico de despegue parcial entre fibra y matriz en materiales compuestos reforzados con fibra, bajo acción de cargas transversales. El análisis contempla la posible presencia de zonas de contacto sin fricción en la porción despegada de la interfase. Se identifica carácter singular de algunas componentes de tensión en la proximidad del fin del pegado, y se calculan numéricamente parámetros asociados a la singularidad mediante un proceso de extrapolación. La comparación de los parámetros conocidos analíticamente con los valores obtenidos resulta satisfactoria. Se sugiere la utilidad potencial de este tipo de análisis para determinar propiedades físicas de la interfase con el apoyo de datos experimentales. Abstract. A numerical study on the mícromechanícal problem of partially debonded fiber and matrix in transversely loaded fiber-reinforced composites is presented. The analysis accounts for the possibility of frictionless contact wnes along the debonded wne of the interface. The Boundary Element Method is used through. The singular nature of sorne stress components near the end of the bonded wne is ídentified, and parameters related to the singularity are numerically calculated by using an extrapolation procedure. Comparison of theoretically known parameters against its calculated values is satisfactory. The potential usefulness of this kind of numerical analysis together with experimental data in order to determine physícal properties of the interface is suggested. l. INTRODUCCION La observación mediante microscopio de grietas transversales en materiales reforzados con fibra, revela que el proceso de agrietamiento comienza por despegues en la interfase fibra-matriz, para niveles relativamente bajos de la carga transversal. Para niveles superiores de carga, estos despegues se conectan a través de la matriz formando una macrogrieta. Por tanto, las propiedades y mecanismos de fallo de la interfase juegan un papel dominante en el comportamiento· del material bajo carga transversal. El objetivo último que anima la realización del presente trabajo es el establecimiento de un procedimiento para la cuantificación de las propiedades resistentes a fractura de la interfase. Ello podría ser llevado a cabo realizando el mejor ajuste posible de los parámetros de un criterio de resistencia dado, al comportamiento observado en la experimentación. La distribución de tensiones y/o desplazamientos, que intervendrá de una u otra manera en el criterio de resistencia, debe ser pues conocida. En resumen, sería necesario contar con: -información precisa acerca del estado tensional en condiciones de despegue de la interfase. -datos experimentales acerca de las condiciones en que se inicia y propaga dicho despegue. El presente trabajo se dedica al primero de los puntos anteriores, estudiándose la viabilidad de obtener numéricamente de modo suficientemente preciso el estado tensional en el tipo de problemas considerado, utilizando discretizaciones basadas en técnicas de contorno, más ventajosas que las técnicas de dominio en el tipo de problemas en elasticidad lineal que se estudian. En la sección 2 se expone el problema a analizar. Su elección responde al criterio de facilidad de materialización y observación en laboratorio, si bien estos aspectos no son abordados en este artículo. Las secciones 3 y 4 describen brevemente las principales características de la técnica numérica que se ha implementado para obtener los resultados que se presentan en la sección 5. En la sección 6 se resumen las conclusiones del estudio. realizado. ANALES DE MECANICA DE LA FRACTURA Vol. 12 (1995) 321 2. DESCRIPCION DEL PROBLEMA El presente análisis se realiza bajo hipótesis de deformación plana en elasticidad lineal. El problema propuesto· es mostrado en la Fig. l. Consiste en una fibra única; cuyas propiedades corresponden a una fibra de vidrio, rodeada por una matriz de dimensiones suficientemente grandes como para poder ser considerada infinita. Las propiedades de la matriz corresponden a una resina epoxy. Se asumirá un ángulo de despegue 2a orientado simétricamente respecto a la dirección de la carga, y separación nula en la zona despegada cuando no actúan cargas exteriores. Es conocido, ver por ejemplo [1], que este tipo de problemas es independiente de la carga, en el sentido de que las zonas de separación y contacto (y la posible partición de esta última en zonas de adhesión y deslizamiento, en el caso de que exista fricción), no dependen del nivel de carga, suponiendo que efectivamente se aplique alguna fracción no nula de la carga. El valor de la tracción exterior no tiene por tanto una significación relevante en este estudio, puesto que todas las tensiones serán proporcionales al nivel de carga; no obstante, en la figura 1 se indica un valor tipo cómodo para la observación durante la experimentación. Solamente se analizará la mitad superior del problema, aprovechando la simetría del mismo. 3. EL METODO DE LOS ELEMENTOS DE CONTORNO EN PROBLEMAS DE CONTACTO Los fundamentos del Método de los Elementos de Contorno en elasticidad lineal están bien descritos en la literatura [3,4], y no serán completamente desarrollados aquí. Para el propósito de esta exposición considérese un . sólido bidimensional elástico sin cargas de volumen. Podernos usar el segundo teorema de Betti para generar tantas ecuaciones corno necesitemos para calcular los parámetros indeterminados de la solución aproximada que hayamos planteado previamente. Tornaremos siempre la solución fundamental de Kelvin para un medio infinito corno solución auxiliar, aunque son posibles otras elecciones. Las ecuaciones que se obtienen son de la forma : ciJ (x)u/x)+ JaDT¡)l' (x,y)uj(y)d5(y) = = lw 'Jiij(x,y)TJ(y)d5(y) (1) dónde an es el contorno del dominio D que ocupa el sólido, "y" simboliza genéricamente un punto del sólido (siempre en su contorno en esta ecuación), "x" es la posición de la carga unidad de la solución auxiliar, "i","j" son índices que varían de 1 a 2 en problemas bidimensionales, uj(y), Tj(y) son los 35x10"6 Pa zona pegada = 2n: - 2a matriz J E=2.8 x1 0"9 Pa l Poiss=0.33 zona no pegada = 2a línea de simetría Fig. l. Problema a analizar y direcciones locales de las componentes de tensión y desplazamiento. A primera vista podría pensarse que si a es pequeño no habría zona de contacto, pero en todo caso parece claro que para a grande, digamos 90º, existirá una zona de contacto entre las superficies despegadas de la fibra y la matriz en las cercanías de a. El estudio numérico debe, por tanto, ser capaz de modelar estas zonas de contacto, lo que por otra parte constituirá la aportación más importante respecto a las soluciones analític¡iS disponibles [2], tal como se comenta en la sección4. vectores de desplazamiento y tensión en , y '!f¡j(x.y), T¡jW(x,y) son los vectores desplazamiento y tensión de la solución auxiliar, actuando su carga unidad en dirección "i". Se entiende sumatorio de 1 a 2 en índices que, corno "j", se repiten en el mismo término. Si "x" se sitúa en el contorno del sólido los dos íntegrandos de (1) son singulares: la integral izquierda exi_ste en el sentido del valor principal de Cauchy, mientras la derecha contiene sólo una singularidad débil, integrable ·en el sentido convencional. Adicionalmente, si "x" está en el contorno, depende de la geometría local del contorno en punto. Si "x" está fuera o dentro del 322 ANALES DE MECANICA DE LA FRACTURA Vol. 12 sólido, los dos integrandos son regulares. Si "x" está fuera es cij=O, y si está dentro, es c¡j=Oij· La técnica de contorno más popular consiste en escribir la ecuación (1) con el punto fuente "x" situado en cada una de las posiciones nodales yP (p= l...n) de una discretización habitual en elementos de contorno. Esta es la técnica empleada principalmente en este trabajo, aunque también se usan ecuaciones obtenidas con "x" fuera del sólido para propósitos particulares que se indicarán más adelante. La solución aproximada, cuyos parámetros se pretende determinar, se plantea del modo habitual, expresando las tensiones y los desplazamientos como combinación lineal de p-funciones NP(y), no nulas en pequeñas porciones del contorno (elementos). Cada una de estas funciones vale 1 en un punto particular del contorno (el nodo p, que se aSQGia a esa función), y cero en las demás posiciones nodales. Con esta elección, los parámetros a determinar son precisamente las· componentes del desplazamiento y tensión en los nodos. La aproximación resulta: uj(y) = í:NP(y)uj(yP) p (2) p La sustitución de (2) en (l) para el punto fuente en xq conduce a: "H'l:Pu.(yP)="G'l:Pt.(yP) (3) .1:,.,¡ IJ J .1:,.,¡ IJ j p p donde: si p =t:. q (4) Se ha elegido aproximar tensiones y desplazamientos en base a tres posiciones nodales yP en cada elemento: un nodo (compartido por dos elementos) en cada extremo del elemento y otro en algún punto intermedio. Habitualmente se denominan elementos parabólicos continuos en la literatura. En cuanto a la geometría, el problema contiene sólo contornos circulares y rectos, que pueden ser descritos exactamente, no siendo necesario recurrir a interpolación. Por ello se ha implementado geometría exactamente circular para los elementos de la interfase. La ecuación (3) planteada para tantos puntos xq como nodos yP tengamos constituye un sistema de ecuaciones algebraicas que permitiría calcular las tensiones y desplazamientos desconocidos en problemas de un solo cuerpo con condiciones de contorno correctamente defmidas. En un problema de contacto entre dos sólidos A y B (Figura 2), las condiciones de contorno en la zona de contacto involucran variables de ambos sólidos, por lo que las ecuaciones (3) planteadas para los sólidos A y B, y las ecuaciones de contacto deben resolverse conjuntamente. Sólido A Fig. 2. Zona de contacto y defmición de ejes locales. Las mencionadas condiciones de contacto, para el caso sin fricción que nos ocupa, deben garantizar el equilibrio, la ausencia de interpenetraciones, y que no haya tracciones en la zona despegada (aunque puede haberlas en lazona pegada). Así, para un punto M en la zona pegada tendremos: ut (M)+ u~ CM:)= O, i=l,2 T¡A(M)=T¡B(M), i=l,2 (5) y para M en la zona despegada, pero en contacto: u~(M)+u~(M)=O T1 A (M)= T1 B (M)~ 0 Ti' (M)= T:(M) =O (6) Finalmente, en la zona de separación de la zona despegada se tiene u~ (M)+ u~ (M) ~O T¡A (M)= T¡B (M)= 0, i= 1,2 (7) En los problemas sin fricción, la solución no depende de la historia de carga, no siendo por tanto necesario plantear la aplicación incremental de la misma. Por tanto, la estrategia de solución se ha basado en la determinación iterativa de las zonas de contacto y separación. Para ello se parte de una estimación de la zona de contacto, se resuelve el conjunto de ecuaciones formado por (3), planteado para A y B, junto a las ecuaciones de contacto, y se comprueba si hay violaciones de la 2ª de las ecs. (6) o de la 1 ª de las (7). Si es así, se modifican algunas de las condiciones de ANALES DE MECANICA DE LA FRACTURA Vol. 12 (1995) 323 contacto supuestas y se reinicia el proceso. La solución se alcanza cuando no se violan las ecuaciones mencionadas. Este procedimiento es una particularización del algoritmo más general que se describe en [5]. En cuanto a la resolución conjunta de las ecuaciones (3) para los sólidos A y B y las de contacto, que debe realizarse en cada iteración del proceso anterior, se ha implementado una resolución a su vez iterativa. Se comienza resolviendo las ecuaciones (3) del sólido matriz suponiendo desplazamientos nulos en las zonas pégadas y de contacto. Las tensiones obtenidas son aplicadas a la fibra mediante las condiciones de contacto. Se obtiene la solución para la fibra bajo estas condiciones, y los desplazamientos obtenidos se aplican, nuevamente a través de las condiciones de contacto, al sólido matriz. El proceso continúa hasta que los desplazamientos de la interfase se estabilizan con un nivel de tolerancia relativa global predeterminado (0.0002 en los casos presentados). Para este tipo de problemas de contacto, en los que uno de los sólidos es mucho más rígido que el otro, el procedimiento descrito suele resultar ventajoso frente a la resolución conjunta de todas las ecuaciones del problema, porque los requerimientos de memoria y el número de operaciones suele ser menor (el proceso converge en dos o tres iteraciones en los casos estUdiados). Adicionalmente, este procedimiento no impide realizar una condensación previa (técnica descrita en [6]) de los grados de libertad que permanecen pasivos. Ha.y que notar, sin embargo que para problemas con materiales de rigidez similar, el proceso converge más lentamente, perdiendo así su atractivo. 4. CONSIDERACIONES ACERCA DEL ESTADO TENSIONAL EN EL BORDE DE GRIETA Si los sólidos fuesen del mismo material, y el ángulo a de despegue fuese muy pequeño, el problema presentaría la típica singularidad de tensión (asociada al modo I) en el lado pegado. Sin embargo, los materiales son distintos y a puede tener cualquier valor, lo que hace difícil anticipar el carácter de la solución de tensiones para la configuración general. Entre otros autores, England [7] observó que la aplicación de la teoría lineal de la elasticidad a un problema análogo al de la Fig.l, considerando bordes de grieta libres, conduce a soluciones oscilantes cerca del fin del pegado, manifestando este autor que dichas soluciones son físicamente incorrectas, puesto que implican interpenetraciones entre los sólidos. Comninou [8-10] ha planteado la resolución de problemas de grietas rectas admitiendo zonas de contacto, no presentando sus soluciones el inaceptable carácter oscilante antes descrito. En [8] se muestra que con tracción transversal a la grieta, situación asimilable a nuestro problema si a es pequeño, la zona de contacto es pequeñísima. Sin embargo, bajo carga combinada [9], situación más parecida a la que se da en nuestro problema en la vecindad del fin del pegado si a no es pequeño, la zona de contacto cubre buena parte de la longitud de la grieta. Lo anterior sugiere que teóricamente tendremos siempre zona de contacto, lo que hace esperar un papel importante del fallo por cortante en la propagación de la grieta [8, 10], con independencia de la orientación de la misma respecto de las cargas exteriores. En otro orden de cosas, la modelización numérica de la discontinuidad de t¡ en el fin del pegado usando el tipo de elemento referido en la sección 3, requiere alguna actuación especial. Se ha implementado una técnica de nodo doble, que para un sólido A consiste en colocar dos nodos en la misma posición geométrica, justo en el punto de discontinuidad, teniendo cada nodo (lA, 2A en la fig.3) defmida su función de aproximación N1A(y), N2A(y), sobre un elemento distinto. La Fig.3. muestra la disposición de nodos dobles, lA, 2A, lB, 2B, en un problema de contacto, asi como los elementos adyacentes. sólido A lA 2A ~$~ _y~ sólido B lB 2B Fig. 3. Nodos dobles en el fin del pegado. En nuestro caso, lA y lB estarán siempre pegados, siéndoles de aplicación las ecuaciones (5). Los nüQos 2A y 2B pueden estar, en la solución numérica, en contacto o bien en separación, aplicándose las ecuaciones (6) o (7). En todo caso, tenemos 16 incógnitas en la discontinuidad (cuatro por nodo), y ocho ecuaciones de contacto. Las otras ocho ecuaciones que precisamos se han obtenido mediante una de las dos opciones siguientes: 1 ª, planteando la ecuación (3) para cada sólido. con xP en lá discontinuidad (4 ecuaciones), y u¡(IA)=u¡(2A), u¡(1B)=ui(2B) (otras 4 ecuaciones), o bien 2ª, planteando dos ecuaciones del tipo (3) en cada sólido, con xP en el exterior de los mismos (8 ecuaciones). En la Fig.3. se indican con aspas, a escala orientativa, las posiciones que se han usado en este trabajo. La opción 1 ª, más simple, es apropiada sólo si 2A y 2B están en 324 ANALES DE MECANICA DE LA FRACTURA Vol. 12 (1995) separación, ya que en otro caso, la primera de las ecuaciones (6) resulta redundante. En todo caso, se obtiene, como debe ocurrir, el mismo desplazamiento en los cuatro nodos, salvo errores numéricos inherentes a la evaluación de las integrales en (4), y de la resolucióñ del sistema lineal de ecuaciones. S. RESULTADOS Las discretizaciones se realizan de modo que la primera distancia entre nodos desde el fin del pegado sea pequeña (típicamente 0.001 º), siendo cada intervalo nodal sucesivo el doble que el anterior. En todos los casos analizados se encuentran oscilaciones de la solución numérica, que afectan únicamente a los tres o cuatro nodos más próximos a a, con independencia de ló cercano que se sitúe el nodo más próximo. Estas oscilaciones se atribuyen al hecho de que en los nodos en e=a. la tensión es teóricamente infinita, y no puede representarse numéricamente. Los resultados en los nodos afectados por esta inestabilidad simplemente se han descartado. 72 o ~ M ~ 00 ~ M M M Wº ~~~~~~~=n==~e -5 -10 -15 -20 -25 -30 -35 -40 -45 Fig. 4. Tensiones calculadas para el caso a = 80º El carácter singular de cada componente de tensión se ha identificado, además de por simple inspección, mediante dos procedimientos: a) comparando la evolución al afinar la discretización, usando discretizacíones con su nodo más próximo á a a 0.01º, 0.001º, y 0.0001º y, b) comprobando la aparición de un segmento casi recto en gráficos análogos a la Fig.6, que representan log(T¡) frente a -log 1 a-e l. La Fig. 4. muestra la distribución de tensiones normales y tangenciales para despegue a = 80º, como ejemplo. Se aprecia evolución singular de la tensión tangencial T2 en la zona pegada (9>a), y de la tensión normal T 1 en la zona de contacto (Eka). La tensión normal T 1 en la zona pegada no presenta carácter singular. Para a = 60º (Fig.5) se aprecia el mismo carácter de cada componente de tensión. 59.5 10 60.5 61 61.5 o~~~~~~==~~ -10 e -20 -30 -40 -50 -60 -70 -80 T(xlOA-7) Pa Fig. 5. Tensiones calculadas para el caso a= 60º Para a mayor que unos 40º, las discretizaciones empleadas detectan alguna zona de contacto, y las soluciones presentan el aspecto de las Figs.4 y 5. Es notorio que para a menor que unos 75º coexistan una tracción (no singular) en la zona pegada con una inmediata zona de contacto en la zona despegada. Si a es menor de aproximadamente 40º, la zona de contacto no es captada por las discretizaciones empleadas, abarcando por tanto menos de 0.0001 º· Parece claro que la solución para bordes libres (numérica o analítica) podría aplicarse para propósitos prácticos en estos casos sin error apreciable. 2.5 Log( T2 x 10"-7 Pa) 2 -·- a=60º ---{j;-- a=80 º 1.5 / -Log ( (9-a) rad ) o 2 3 Fig. 6. T 2 frente a (9--a.), a escala logarítmica. La aparición de una zona recta en las proximidades de 9=a en gráficas logarítmicas de tensión frente a distancia al fin del pegado, es mostrada en la Fig.5 para dos casos representativos. La pendiente de estos segmentos rectos es el orden de la singularidad, mientras que su ordenada en el origen es fácilmente relacionable .con el factor de intensificación. El orden de la singularidad así calculado se muestra en la Tabla 1, junto con el coeficiente de correlación r2 obtenido en el ajuste por mínimos cuadrados de los segmentos "casi rectos" identificados. Se han usado seis puntos para cada ajuste. ANALES DE MECANICA DE LA FRACTURA Vol. 12 (1995) 325 Para juzgar la prec1s1on de los resultados, cabe comparar el orden calculado de las singularidades con el del problema, localmente similar, de interfase recta, con zona pegada y subsecuente zona de contacto. Este orden se calcula exactamente a partir de las expresiones dadas en [11], demostrándose que para el caso sin fricción, vale 0.5, y que no depende de las constantes elásticas. La comparación con los órdenes de singularidad obtenidos revela una diferencia menor del 1%. Adicionalmente, existe una concordancia excelente entre los resultados obtenidos cuando no se capta zona de contacto (a<40º; no se muestran estos. resultados porque su interés es anecdótico disponiéndose de solución analítica) y la correspondiente solución analítica [2]. Tabla 1. Ordenes de singularidades calculados numéricamente, y contacto detectado. Zona pegada Contacto a t¡ r2 orden 60º t, 0.9998 0.505 0.1º 80º t, 0.9996 0.504 7º 6. CONCLUSIONES Se ha reconocido, en base a los resultados numéricos y a soluciones analíticas ya disponibles, que la presencia de materiales diferentes en el problema de despegue parcial entre fibra y matriz implica la presencia de zonas de contacto en la zona despegada. Se ha mostrado que la tensión normal en la zona pegada no es singular, dándose la sorprendente posibilidad de que sea de tracción, a pesar de la presencia inmediata de lá zona de contacto. La compresión en la zona de contacto y la tensión tangencial en la zona pegada se han identificado como singulares. Este carácter de las soluciones es acorde con el presentado en [8,9] para grietas rectas. Cuando la zona de contacto es muy pequeña, la discretización puede no captarla. En estos casos parece claro que la solución analítica disponible se podría aplicar sin error apreciable, tal como se hace en [2]. De la precisión obtenida en el estudio numérico realizado, cabe concluir la validez del método seguido para obtener el estado de tensión de modo suficientemente preciso, justamente en los casos en que la aplicación de la solución analítica mencionada produciría errores importantes. Desde un punto de vista práctico, el parámetro de interés será el factor de intensificación, o el propio estado de . tensiones y desplazamientos cerca del borde de grieta. La precisión obtenida para el orden de la singularidad, permite estimar la precisión obtenida en otras variables del problema. Agradecimiento. El apoyo económico aportado por la Caja de Ahorros de Salamanca bajo su Proyecto para Jóvenes Investigadores es cordialmente agradecido. REFERENCIAS [1] Gladwell G.M.L., "Contact problems in the classical theory of elasticity", Sijthoof & Noordhoff Int. Publ., Germantown, Maryland, USA (1980). [2] Toya M. "A crack along the interface of a circular inclusion embedded in an infinite so lid", J. Mech. Phys. Solids, VoL 22,325-348 (1973). [3] Brebbia C.A., Telles J.C.F. and Browen L.C., "Boundary Element Techniques", SpringerVerlag, Berlin,(l984). [4] Hall W.S. "The boundary Element Method", Kluwer Academic Publishers, Dordretch, (Holanda) (1994) [5] París F. & Garrido J.A., "An incremental procedure for friction contact contact problems with the B.E.M.", Int. J. for Engineering Analysis with Boundary Elements, Vol 6, nº 4, 202-213 (1989). [6] París F. & Garrido J.A., "Aspectos numéricos de la aplicación del M:E:C: al problema de contacto", Rev. Int. de Métodos Numéricos para el Cálculo y Diseño en Ingeniería, Vol 2, nº 1, 43-62 (1986). [7] England A.H. "An are crack around a circular elastic inclusion", Journal of Applied Mechanics, p.637, septiembre (1966). [8] Comninou M., "The interface crack", Journal of Applíed Mechanics, p.631, diciembre (1977). [9] Comninou M. & Schmueser D.,"The interface crack in a conbined tension-compression and shear field", Journal of Applied Mechanics, p.345, junio (1979) [10] Comninou M."An overview on interface cracs", Engineering Fracture Mechanics, Vol. 37, nº 1, 197-208 (1990) [11] Dempsey J.P. & Sinclair G.B., 'On the singular behaviour at the vertex of a bi-material wedge', Journal of Elasticity, Vol.ll no. 3, 317 -327 (1981).