Nuevo elemento finito para singularidades logarítmicas de tensión en grietas de interfase tipo Winkler en modo III

NUEVO ELEMENTO FINITO PARA SINGULARIDADES LOGAR´ ITMICAS DE TENSI ´ ON EN GRIETAS DE INTERFASE TIPO WINKLER EN MODO III A. V´ azquez, V. Mantiˇ c, M. Mu˜ noz-Reja, L. T´ avara∗ Grupo de Elasticidad y Resistencia de Materiales, Escuela T´ ecnica Superior de Ingenier´ ıa Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092 Sevilla, Espa˜ na Persona de contacto: ltav[email protected] RESUMEN El objetivo de este trabajo es aumentar la precisi´ on de las soluciones obtenidas usando el m´ etodo de los elementos finitos (MEF) para grietas de interfases d´ ebiles (tipo muelle) entre adherentes el´ asticos lineales. Este modelo de fractura de interfases se denomina Modelo de Interfase El´ astico Lineal Fr´ agil (LEBIM por sus siglas en ingl´ es), y puede utilizarse, para analizar la fractura de uniones adhesivas que incluyen una capa fina de adhesivo. Para ello, se propone un elemento finito singular especial triangular de 5 nodos. Las funciones de forma singular desarrolladas reproducen el comportamiento de una soluci´ on asint´ otica anal´ ıtica desarrollada previamente en direcci´ on radial. Estos elementos especiales se implementan en un c´ odigo MEF escrito en Matlab. El nuevo elemento permite modelar grietas de interfase sin necesidad de utilizar mallas excesivamente refinadas, que es una de las desventajas actuales del LEBIM cuando se consideran interfases muy r´ ıgidas. Los resultados num´ ericos con mallas uniformes muestran que el nuevo elemento proporciona resultados significativamente m´ as precisos que cuando se usan solamente elementos est´ andar, especialmente para mallas relativamente gruesas. PALABRAS CLAVE: Modelo de Interfase El´ astica-Lineal Fr´ agil (LEBIM), muelle, interfase, elemento singular, MEF ABSTRACT The aim of this research is to increase the accuracy of Finite Element Method (FEM) solutions for cracks propagating in spring-type interfaces between linear elastic adherents. This computational model is sometimes called Linear Elastic Brittle Interface Model (LEBIM) and it can be used to analyse the fracture behaviour of adhesive joints. Thus, a special 5-node triangular crack-tip finite element is proposed. The developed special singular shape functions reproduce the radial behaviour of a previously developed analytical asymptotic solution. These special elements are implemented in a FEM code written in Matlab. The new element allows to model interface cracks without the need of using excessively refined FEM meshes, which is one of the current disadvantages in the use of LEBIM when very stiffinterfaces are considered. Numerical results using uniform meshes show that the use of the new singular element provides significantly more accurate results than those obtained using the standard elements, especially for relatively coarse meshes. KEYWORDS: Linear Elastic Brittle Interface Model (LEBIM), spring, interface, singular element, FEM 1. INTRODUCCI ´ ON La soluci´ on de la ecuaci´ on de Laplace se utiliza para diferentes tipos de problemas f´ ısicos tales como: transferencia de calor, electrost´ atica, conducci´ on el´ ectrica, magnetost´ atica, etc. [1]. En el presente trabajo, la soluci´ on de dicha ecuaci´ on se utiliza para describir el comportamien- to mec´ anico de una grieta en modo III puro (conocido en mec´ anica de la fractura como problema antiplano) entre materiales el´ asticos lineales. La figura 1 muestra las condiciones de contorno para una grieta de interfase bajo las condiciones antes mencionadas. Una forma de modelar las grietas que se propagan a lo largo de una interfase d´ ebil (que puede representar, por ejemplo, una capa adhesiva) es mediante el Modelo de Interfase El´ astico Lineal Fr´ agil (LEBIM por sus siglas en ingl´ es) [2, 3, 4, 5]. En este modelo la interfase se modela como una distribuci´ on continua de muelles con una ley el´ astica lineal fr´ agil que relaciona el desplazamien- to relativo (separaci´ on de materiales en el modo I) y los componentes del vector tensi´ on [6]. En cualquier caso, el crecimiento de la grieta en LEBIM necesita un c´ alculo preciso de la Tasa de Liberaci´ on de Energ´ ıa (ERR por sus siglas en ingl´ es), que depende de los valores del vector tracci´ on y de los desplazamientos relativos en el v´ ertice de la grieta [6]. Por ello, se requieren mallas muy finas cuando se utiliza una aproximaci´ on num´ erica, especialmente para interfases muy r´ ıgidas [7, 8, 9]. El LEBIM puede considerarse como un caso l´ ımite del conocido Modelo de Zona Cohesiva (CZM por sus siglas en ingl´ es) [5, 10, 11, 12, 13, 14]. En el pasado, hubo numerosos intentos para encontrar la serie asint´ otica de soluciones singulares en las proximidades del v´ ertice de la grieta en modelos de mec´ anica de fractura no cl´ asicos como el LEBIM [4, 15, 16]. Recientemente se introdujo una solu- Revista de Mecánica de la Fractura Vol.5 (2023) 213 2a zx y k (a) 2L x y k (b) Figura 1: Condiciones de contorno asociadas a un problema de grieta de interfase modelizado en modo de fractura III puro, es decir, problema antiplano, (a) vista 3D y (b) vista 2D. La parte sombreada representa el subdominio que puede modelizarse una vez aplicadas las condiciones de simetr´ıa. ci´ on general para este tipo de interfases en [17], y luego se particulariz´ o para una grieta trabajando en Modo III en [18]. El objetivo principal de este trabajo es maximizar la precisi´ on de las soluciones num´ ericas, especialmente en la predicci´ on de la propagaci´ on de grietas que aparecen en diversos materiales y sus uniones, y minimizar los recursos computacionales utilizados, evitando la necesidad de un refinamiento excesivo de la discretizaci´ on de s´ olidos, que a menudo es necesario para lograr la convergencia de la soluci´ on num´ erica en la actualidad. Los elementos especiales desarrollados (basados en la soluci´ on asint´ otica de [18]) implementados en un c´ odigo MEF permitir´ an la discretizaci´ on de s´ olidos con grietas sin necesidad de mallas refinadas, que es uno de los problemas actuales en el uso de LEBIM y CZMs y que suelen conllevar problemas de convergencia num´ erica de la soluci´ on. 2. ELEMENTO ESPECIAL PARA SINGULARIDAD LOGAR´ ITMICA DE TENSI ´ ON Tomando como base la forma de la soluci´ on local, se desarrolla un elemento finito capaz de reproducir ese comportamiento. As´ ı, en primer lugar se describen las funciones de forma 1-D y, a continuaci´ on, se presentan las funciones de forma 2-D. 2.1. Funciones de forma 1D Para deducir el conjunto de funciones de forma de un elemento 1D, con nodos en r1=0, r2=0.5 y r3=1, se consideran las siguientes funciones base definidas en el intervalo [0,1]: e N1(r)=1; e N2(r)=r;e N3(r)=rln (r) (1) A continuaci´ on, las nuevas funciones de forma del elemento 1D se obtienen mediante una combinaci´ on lineal de las funciones de base: Ni(r)=Ci j e Nj(r),(2) imponiendo la conocida propiedad de interpolaci´ on de Lagrange Ni(rj)=δi j,(3) donde δi j es el delta de Kronecker. Resolviendo el sistema de ecuaciones, se obtiene el siguiente conjunto de funciones: N1(r)=1−r−rln (r) ln 1 2,(4) N2(r)=2rln (r) ln 1 2,(5) N3(r)=r−rln (r) ln 1 2.(6) 2.2. Funciones de forma 2D Una vez obtenidas las funciones de forma en el espacio unidimensional, es necesario adaptarlas a un espacio de referencia bidimensional ξ-η. Para ello se genera un elemento de referencia con 6 nodos que incluye un comportamiento lineal en ηy con funciones de forma espec´ ıficas adaptadas a ξ, sustituyendo rpor ξen las Ecs. (4-6). Una vez obtenida la evoluci´ on de cada funci´ on de forma en el eje ξ, se pueden adaptar a un espacio bidimensional incluyendo una evoluci´ on lineal de dos nodos en el eje η, como se observa en las siguientes funciones de forma: N1(ξ, η)=(1 −η)     1−ξ−ξln (ξ) ln (1 2)      (7) N2(ξ, η)=(1 −η)     ξ−ξln (ξ) ln (1 2)      (8) N3(ξ, η)=η     ξ−ξln (ξ) ln (1 2)      (9) N4(ξ, η)=η     1−ξ−ξln (ξ) ln (1 2)      (10) N5(ξ, η)=(1 −η)     2ξln (ξ) ln (1 2)      (11) N6(ξ, η)=η     2ξln (ξ) ln (1 2)      (12) tomando ξyηvalores entre [0,1]. Funciones de forma singulares: Finalmente partiendo de las funciones de forma definidas para el espacio ξ−ηy bas´ andose en [19], el elemento singular cuadril´ atero se colapsa en un elemento triangular donde todos los nodos del lado con ξ=0 se colapsan en un nodo solamente. 3. RESULTADOS NUM ´ ERICOS Y ESTUDIO DE CONVERGENCIA 3.1. Definici´on de la malla y el problema analizado La Fig. 2 muestra el problema de referencia que se utilizar´ a para analizar el comportamiento del elemento singular presentado en este trabajo. En esta configuraci´ on se aplican condiciones de contorno tipo Neumann con σyz =1 en el borde superior, condiciones de contorno tipo Robin (muelle) en la mitad derecha del borde inferior y los bordes restantes son bordes libres de tensi´ on. Adem´ as, para simplificar el problema se ha utilizado condiciones de simetr´ ıa en el plano que corta el s´ olido por el centro de la grieta. x y 2k NR N (libre) N (carga uniforme) N (Plano de simetría) Figura 2: Condiciones de contorno del problema analizado. Se emplea un refinamiento tipo hde mallas uniformes para analizar el comportamiento de convergencia del nuevo elemento triangular singular para modelar grietas que se propagan a lo largo de interfases de muelles. Se utilizan elementos triangulares lineales est´ andar en todo el dominio, excepto en la proximidad del v´ ertice de la grieta, donde se utilizan los elementos singulares. A efectos comparativos, tambi´ en se emplean las mismas mallas utilizando ´ unicamente elementos triangulares lineales est´ andar en todo el dominio, incluida la vecindad de la punta de la grieta. Para estos estudios, se utilizar´ an mallas cuadradas compuestas de elementos triangulares. Los resultados se obtienen mediante una implementaci´ on del MEF en el software Matlab. Se evaluar´ an los resultados obtenidos con diferentes n´ umeros de elementos de la malla que se muestran en la Fig. 3. Figura 3: Malla estudiada para N=4. El n´ umero total de elementos se asocia a la variable Nque corresponde al n´ umero de elementos a lo largo del borde inferior con condiciones de contorno Robin (la mitad derecha). Entonces, Ntambi´ en se puede definir como: N=a h(13) donde aes la zona libre de tensiones (normalmente asociada a la semi-longitud de la grieta) y hla longitud caracter´ ıstica de un elemento. 3.2. Resultados num´ericos Se probaron diferentes valores de rigidez de la interfase variando el valor de δ, manteniendo constante la geometr´ ıa modelada. δ=2kLch µ=4ka µ(14) donde kes la rigidez del muelle (en el problema comple- to), Lch es una longitud caracter´ ıstica del problema (aqu´ ı se utiliz´ o la longitud de grieta 2a) y µes el m´ odulo de cizalladura. Los valores δanalizados son 0.1, 1, 10, 100 y 1000. La Fig. 4 muestra el campo de desplazamiento obtenido a lo largo del borde inferior para la malla analizada y para dos valores de δ. Es interesante ver que la soluci´ on tiende a la soluci´ on de interfase perfecta cuando la rigidez de la interfase crece. Tambi´ en se analiza la mejora obtenida cuando se utilizan elementos singulares. Dado que no se tiene una soluci´ on anal´ ıtica de este problema, se ha calculado el error tomando como referencia los resultados obtenidos con la malla Revista de Mecánica de la Fractura Vol.5 (2023) 215 m´ as fina (es decir el mayor n´ umero de elementos en el borde Robin) que en este caso es N=512. Los resultados muestran que la mejora obtenida cuando se utilizan elementos singulares es relevante. El error relativo en el v´ ertice de la grieta se define como: erel =τtip(N)−τtip(Nmax) τtip(Nmax)(15) -0.5 -0.25 0 0.25 0.5 x 2a 1.8 2 2.2 2.4 2.6 2.8 u7 2a= N=2 sing. elem. N=512 sing. elem. N=2 std. elem. N=512 std. elem. (a) -0.5 -0.25 0 0.25 0.5 x 2a -0.1 0 0.1 0.2 0.3 0.4 0.5 u7 2a= N=2 sing. elem. N=512 sing. elem. N=2 std. elem. N=512 std. elem. (b) Figura 4: Soluciones en desplazamiento obtenidas para (a) δ=1y (b) δ=1000 En la Fig. 5, se presenta la variaci´ on del error respecto al aumento del n´ umero de elementos, N, en el v´ ertice de la grieta para δ=0.1, 10 y 1000. Los resultados muestran que cuando la rigidez de la interfaz es peque˜ na (δ=0.1), el error relativo es similar para ambos tipos de elementos analizados, est´ andar y singular. Las principales diferencias se producen para los valores m´ as peque˜ nos de N. Por otra parte, cuando la rigidez de la interfaz es mayor (δ≥10), el error relativo entre las distintas configuraciones se hace m´ as evidente. Adem´ as, las principales diferencias se producen para los valores m´ as altos de N e interfases r´ ıgidas (δ=1000). Tambi´ en es interesante observar que los errores obtenidos al utilizar elementos singulares son siempre menores que los obtenidos utilizando ´ unicamente elementos est´ andar. 0 0.5 1 1.5 2 2.5 log(N) -6 -5 -4 -3 -2 -1 0 log(erel) elem. singular elem. estándar (a) 0 0.5 1 1.5 2 2.5 log(N) -5 -4 -3 -2 -1 log(erel) elem. singular elem. estándar (b) 0 0.5 1 1.5 2 2.5 log(N) -3.5 -3 -2.5 -2 -1.5 log(erel) elem. singular elem. estándar (c) Figura 5: Variaci´on del error relativo con respecto al aumento del n´umero de elementos en el v´ertice de la grieta para (a) δ=0.1, (b) δ=10 y (c) δ=1000. Revista de Mecánica de la Fractura Vol.5 (2023) 216 4. CONCLUSIONES En el presente trabajo se ha desarrollado un nuevo elemento finito. Este elemento es capaz de mejorar la preci- si´ on de las soluciones MEF para grietas que crecen a lo largo de interfases tipo Winkler (resorte) entre adherentes el´ asticos lineales. Las funciones de forma propuestas se basan en la soluci´ on el´ astica asint´ otica de tensiones con singularidad logar´ ıtmica en las proximidades del v´ ertice de la grieta de interfase tipo muelle trabajando en modo III de fractura. Adem´ as, el nuevo elemento puede aplicarse a cualquier otro problema f´ ısico regido por la ecuaci´ on de Laplace. El elemento reproduce el comportamiento radial de la soluci´ on asint´ otica, es triangular y est´ a formado por 5 nodos. El elemento se implementa en un c´ odigo MEF escrito en Matlab. Los resultados num´ ericos obtenidos muestran que el nuevo elemento permite modelar grietas de interfase sin necesidad de utilizar mallas de elementos finitos excesivamente refinadas, incluso para interfaces bastante r´ ıgidas. Un an´ alisis de refinamiento tipo h, utilizando mallas uniformes, mostr´ o que el nuevo elemento singular proporciona resultados m´ as precisos que los elementos finitos est´ andar, especialmente cuando se utilizan mallas gruesas. El uso del elemento propuesto permitir´ a minimizar los recursos computacionales al modelar grietas de interfase relativamente r´ ıgidas. AGRADECIMIENTOS Este estudio ha contado con el apoyo del Ministerio de Ciencia e Innovaci´ on de Espa˜ na y del Fondo Europeo de Desarrollo Regional (Proyectos PGC2018-099197-B-I00 y PID2021-123325OB-I00) y de la Junta de Andaluc´ ıa y el Fondo Social Europeo (Proyecto P18-FR-1928). REFERENCIAS [1] F. Par´ ıs and J. Ca˜ nas. Boundary Element Method, Fundamentals and Applications. Oxford University Press: Oxford, 1997. [2] L. Prandtl. A thought model for the fracture of brittle solids. Zeitschrift f¨ur Physikr Angenwandte Mathematik und Mechanik, 13(2):129–133, 1933. [3] V.M. Entov and RL. Salganik. On the Prandtl brittle fracture model. Mechanics of Solids, 3:79– 89, 1968. (translated from Russian). [4] S. Lenci. Analysis of a crack at a weak interface. International Journal of Fracture, 108(3):275–290, 2001. [5] P. Cornetti, V. Mantiˇ c, and A. Carpinteri. Finite Fracture Mechanics at elastic interfaces. International Journal of Solids and Structures, 49:1022–1032, 2012. [6] V. Mantiˇ c, L. T´ avara, A. Bl´ azquez, E. Graciani, and F. Par´ ıs. A linear elastic - brittle interface model: Application for the onset and propagation of a fibre-matrix interface crack under biaxial transverse loads. International Journal of Fracture, 195:15– 38, 2015. [7] L T´ avara, V Mantic, E Graciani, J Ca˜ nas, and F Par´ ıs. Analysis of a crack in a thin adhesi- ve layer between orthotropic materials: an application to composite interlaminar fracture toughness test. Computer Modeling in Engineering &Sciences(CMES), 58(3):247–270, 2010. [8] M. Mu˜ noz Reja, L. T´ avara, V. Mantiˇ c, and P. Cornetti. Crack onset and propagation at fibre-matrix elastic interfaces under biaxial loading using finite fracture mechanics. Composites Part A, 82:267– 278, 2016. [9] M Mu˜ noz-Reja, L T´ avara, and Vladislav Mantiˇ c. Convergence of the BEM solution applied to the CCFFM for LEBIM. In Key Engineering Materials, volume 774, pages 355–360. Trans Tech Publ, 2018. [10] M. Bialas and Z. Mr´ oz. Modelling of progressive interface failure under combined normal compression and shear stress. International Journal of Solids and Structures, 42(15):4436–4467, 2005. [11] N. Valoroso and L. Champaney. A damagemechanics-based approach for modelling decohesion in adhesively bonded assemblies. Engineering Fracture Mechanics, 73(18):2774–2801, 2006. [12] M.E. Jim´ enez, J. Ca˜ nas, V. Mantiˇ c, and J.E. Ortiz. Numerical and experimental study of the interlaminar fracture test of composite-composite adhesively bonded joints. (in Spanish). Materiales Compuestos 07, Asociaci´on Espa˜nola de Materiales Compuestos, Universidad de Valladolid, pages 499–506, 2007. [13] P. Cornetti, A. Sapora, and A. Carpinteri. Short cracks and V-notches: Finite Fracture Mechanics vs. Cohesive Crack Model. Engineering Fracture Mechanics, 168:2–12, 2016. [14] R. Dimitri, P. Cornetti, V. Mantiˇ c, M. Trullo, and L. De Lorenzis. Mode-I debonding of a double cantilever beam: A comparison between cohesive crack modeling and finite fracture mechanics. International Journal of Solids and Structures, 124:57 – 72, 2017. [15] F. Erdogan. Fracture mechanics of interfaces, In: Damage and Failure of Interfaces. Balkema Publishers: Rotterdam, 1997. [16] G. Mishuris. Interface crack and nonideal interface concept (mode III). International Journal of Fracture, 107:279–296, 2001. [17] S. Jim´ enez-Alfaro, V.r Villalba, and V. Mantiˇ c. Singular elastic solutions in corners with spring boun- Revista de Mecánica de la Fractura Vol.5 (2023) 217 dary conditions under anti-plane shear. International Journal of Fracture, 223(1):197–220, 2020. [18] S. Jim´ enez-Alfaro and V. Mantiˇ c. FEM benchmark problems for cracks with spring boundary conditions under antiplane shear loadings. Aerotecnica Missili &Spazio, 99:309–319, 2020. [19] R.S. Barsoum. On the use of isoparametric finite elements in linear fracture mechanics. International journal for numerical methods in engineering, 10(1):25–37, 1976. Revista de Mecánica de la Fractura Vol.5 (2023) 218