Implementación en Abaqus de la minimización de la energía total bajo una condición de tensión para predecir el inicio y la propagación de grietas en materiales homogéneos frágiles

IMPLEMENTACI ´ ON EN ABAQUS DE LA MINIMIZACI ´ ON DE LA ENERG´ IA TOTAL BAJO UNA CONDICI ´ ON DE TENSI ´ ON PARA PREDECIR EL INICIO Y LA PROPAGACI ´ ON DE GRIETAS EN MATERIALES HOMOG ´ ENEOS FR ´ AGILES K. Ambikakumari Sanalkumar, 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: [email protected] RESUMEN Para la predicci´ on del inicio y del crecimiento de grietas en materiales fr´ agiles se utiliza el criterio acoplado de la Mec´ anica de Fractura Finita, que considera que las grietas pueden avanzar por pasos finitos, y requiere que adem´ as del criterio de energ´ ıa se cumpla tambi´ en un criterio tensional. En el procedimiento num´ erico implementado mediante una rutina de usuario UINTER en ABAQUS, los avances finitos de la grieta, que cumplen el criterio tensional, se puentean mediante una distribuci´ on continua de muelles el´ astico-lineales, que permite definir una variable de da˜ no a lo largo de estos potenciales avances finitos de grietas. As´ ı, la energ´ ıa total (suma de la energ´ ıa potencial y la energ´ ıa disipada) es convexa por separado respecto al campo de los desplazamientos y a la variable de da˜ no, que permite aplicar el algoritmo de minimizaci´ on alternante para la minimizaci´ on de la energ´ ıa total. Este enfoque, que considera los pasos de carga, es m´ as vers´ atil que la formulaci´ on cl´ asica del criterio acoplado, proporcionando una herramienta que en el futuro pueda permitir resolver problemas complejos de relevancia industrial. Se presenta un ejemplo num´ erico, que muestra que el procedimiento num´ erico implementado proporciona predicciones precisas para el crecimiento de grietas. PALABRAS CLAVE: CCFFM, PMTE-SC, FEM, UINTER en ABAQUS ABSTRACT The coupled criterion of Finite Fracture Mechanics is used for the prediction of crack initiation and propagation in brittle materials. It is considered that cracks propagate by finite steps, requiring that both an energy and a stress criterion must be fulfilled. In the numerical procedure, implemented by means of a UINTER user subroutine in ABAQUS, the finite crack advances that fulfil the stress criterion are bridged by a continuous distribution of linear-elastic springs, which allows defining a damage variable along these potential finite crack advances. Thus, the total energy (sum of potential and dissipated energy) is separately convex with respect to the field of displacements and the damage variable, which allows the alternating minimisation algorithm to be applied for the minimisation of the total energy. This approach, which considers loading steps, is more versatile than the classical formulation of the coupled criterion, providing a tool that in the future may allow solving complex problems of industrial relevance. A numerical example is presented, showing that the implemented numerical procedure provides accurate predictions for the crack propagation. KEYWORDS: CCFFM, PMTE-SC, FEM, UINTER in ABAQUS 1. INTRODUCCI ´ ON La presencia de grietas en materiales fr´ agiles puede reducir dr´ asticamente la vida ´ util de una estructura o provocar fallos catastr´ oficos. A pesar que en el pasado se han realizado muchas investigaciones para estudiar este tipo de fractura, la iniciaci´ on de grietas bajo cargas cuasi-est´ aticas es a´ un hoy en d´ ıa un tema controvertido. La Mec´ anica de Fractura El´ astica Lineal (LEFM por sus siglas en ingl´ es) no puede predecir el inicio de la grieta en un concentrador de tensiones o en una singularidad d´ ebil, ni la desviaci´ on de la trayectoria de la grieta en el interior de un material cuando alcanza una interfase. En las ´ ultimas dos d´ ecadas, los investigadores se han centrado en desarrollar modelos de la mec´ anica de la fractura eficientes y precisos para la predicci´ on num´ erica de iniciaci´ on de grietas utilizando nuevos enfoques. El modelo de zona cohesiva (CZM por sus siglas en ingl´ es) [1, 2] es el m´ as utilizado entre los modelos no cl´ asicos de la mec´ anica de la fractura. Aplica hip´ otesis diferentes a las utilizadas en el LEFM para evitar la singularidad de tensiones en el v´ ertice de la grieta y no requiere una grieta inicial para predecir el origen de la fractura. Adem´ as, puede predecir varias trayectorias del da˜ no a la vez cuando se introducen varios elementos cohesivos [3, 4]. Por otro lado, en el denominado modelado de la fractura por “Phase Field”, desarrollado a partir del trabajo de Francfort y Marigo [5], se considera un modelo de gradiente de da˜ no que predice el inicio y la propagaci´ on de la grieta cuando un par´ ametro de longitud se define Revista de Mecánica de la Fractura Vol.5 (2023) 225 adecuadamente como una longitud caracter´ ıstica del material. Este modelo ha ganado una gran popularidad recientemente, pero desde la comunidad cient´ ıfica se han expresado preocupaciones acerca de las posibilidades de las funciones de degradaci´ on y de la energ´ ıa de fractura y del par´ ametro de longitud para representar adecuadamente situaciones complejas que involucran, por ejemplo, a compresiones, contacto, anisotrop´ ıa, iniciaci´ on de grietas, etc. Por otro lado, en los ´ ultimos a˜ nos, se ha hecho popular la Mec´ anica de la Fractura Finita (FFM por sus siglas en ingl´ es), la cual permite considerar el inicio o el avance de la fractura como un incremento finito de la longitud de la grieta. Por tanto, no mantiene la hip´ otesis de Griffith, que establece que el crecimiento de grietas es infinitesimal. En el marco de la FFM, Leguillon [6] introdujo el criterio acoplado de la FFM (CCFFM) que es la base de este trabajo. En el criterio acoplado, los criterios de tensi´ on y de energ´ ıa son dos condiciones necesarias pero no suficientes separadamente para permitir el inicio o el avance instant´ aneo de una grieta. El criterio acoplado se ha utilizado de manera efectiva para predecir el da˜ no en diferentes materiales, incluidos los materiales compuestos [7, 8, 9, 10, 11, 12, 13]. En este marco, Mantiˇ c [14] desarroll´ o una nueva formulaci´ on del criterio acoplado, donde se investigaron varios aspectos y aplicaciones de la minimizaci´ on de la energ´ ıa total para el inicio y la propagaci´ on de grietas en materiales fr´ agiles, asumiendo la evoluci´ on de la grieta de forma cuasi-est´ atica. El Principio de Energ´ ıa Total M´ ınima bajo una Condici´ on de Ten- si´ on (PMTE-SC por sus siglas en ingl´ es), introducido en [14], es m´ as vers´ atil debido a su capacidad para ser utilizado en una implementaci´ on computacional gen´ erica con una t´ ecnica de paso de tiempo. Adem´ as, el funcional de la energ´ ıa total es convexo por separado respecto a las variables de da˜ no y de desplazamientos, lo que permite que las t´ ecnicas de optimizaci´ on sean eficientes y robustas para minimizar la energ´ ıa total. 2. PMTE-SC - PRINCIPIO DE ENERG´ IA TOTAL M´ INIMA BAJO UNA CONDICI ´ ON DE TENSI ´ ON La hip´ otesis de Leguillon [6] postula que la combinaci´ on de los dos criterios de tensi´ on y de energ´ ıa puede predecir el inicio de una nueva grieta que tiene una longitud (cuando estamos estudiando un problema 2D) o un ´ area (cuando estamos estudiando un problema 3D) finita. El criterio de tensi´ on usado en este trabajo requiere que la tensi´ on normal al plano de la posible grieta finita ∆Sc supere un valor cr´ ıtico σc. Este criterio debe cumplirse en cada punto xde la nueva superficie del material ∆Sc. Siguiendo la formulaci´ on en [14, 15, 16], se pude escribir una expresi´ on bastante general del criterio tensional mediante un funcional homog´ eneo de grado 1 respecto al campo tensional, f(σi j,∆Sc): f(σi j,∆Sc)≥σc.(1) Se debe tener en cuenta que este criterio es evaluado antes del inicio de la fractura en todo el dominio, Ω, y permite obtener todas las posibles superficies de grieta, ∆Sc⊂Ω. Posteriormente el criterio energ´ etico es evaluado bas´ andose en el balance energ´ etico incremental, ya que se eval´ ua considerando el estado previo y posterior a la fractura finita definida por ∆Sc: ∆Π(∆Sc)+ ∆Ek+ ∆R(∆Sc)=0,(2) donde ∆Π(∆Sc) es la variaci´ on de la energ´ ıa potencial, ∆R(∆Sc) es la energ´ ıa disipada en la formaci´ on de la nueva superficie de grieta, ∆Sc, y ∆Ekes el incremento de energ´ ıa cin´ etica involucrada en el proceso, que siempre ser´ a mayor que cero. Por tanto, al igual que [6], que desprecia el intercambio de calor en el balance energ´ etico, se puede expresar el balance energ´ etico como: ∆Π(∆Sc)+ ∆R(∆Sc)≤0,(3) que indica que la suma del incremento de energ´ ıa potencial y disipada debe permanecer constante o disminuir para que se produzca la grieta finita. Alternativamente y de forma b´ asicamente equivalente a la formulaci´ on tradicional del CCFFM de Leguillon, el PMTE-SC, introducido en [14], tambi´ en predice el inicio y/o el crecimiento de grietas por incrementos finitos. En este trabajo, al igual que en el CCFFM cl´ asico, el PMTE-SC tambi´ en asume la evoluci´ on del problema de forma cuasiest´ atica sin tener en cuenta las fuerzas inerciales. Por tanto, el PMTE-SC se puede escribir como la siguiente minimizaci´ on bajo una condici´ on de tensi´ on: m´ ın ∆Sc [∆Π(∆Sc)+ ∆R(∆Sc)],(4) donde la condici´ on de tensi´ on significa que la minimizaci´ on de energ´ ıa total s´ olo es considerada en aquellas nuevas superficies de grietas para las que se verifica el criterio tensional para unas condiciones de contorno dadas. El concepto clave que subyace en esta formulaci´ on es que una nueva superficie de grieta puede aparecer solo en aquellas regiones donde aparecen tensiones suficientemente altas antes de la fractura, lo que representa que la formaci´ on de grietas se inhibe si las tensiones son demasiado peque˜ nas, aunque haya una cantidad suficiente de energ´ ıa disponible para liberarse. La Eq. (4) es una formulaci´ on (incremental) del principio de m´ axima disminuci´ on de la energ´ ıa total del sistema, el cual puede dar lugar a un problema de minimizaci´ on no convexo que es muy dif´ ıcil de resolver usando un algoritmo de optimizaci´ on de la energ´ ıa total. En consecuencia, se ha utilizado un algoritmo de minimizaci´ on alterna (AMA) que separa el problema de minimizaci´ on del funcional de energ´ ıa en dos t´ erminos independientes: la Revista de Mecánica de la Fractura Vol.5 (2023) 226 variable da˜ no y el campo de desplazamientos, bas´ andose en que este funcional sea convexo respecto a cada una de estas variables por separado. Esta t´ ecnica de optimizaci´ on ha mostrado ser eficiente y robusta para minimizar la energ´ ıa total. Por otro lado, en este trabajo se ha utilizado un enfoque simplificado para modelar las discontinuidades de tensiones que aparecen al producirse la fractura del material. Para ello una vez que se ha definido la posible superficie de la grieta que satisface el criterio tensional, ∆Sc, se define una distribuci´ on continua de muelles lineales con espesor cero, que se conoce como Interfaz de Muelles El´ astica Lineal (LESI por sus siglas en ingl´ es). A diferencia de los elementos con espesor finito, en los elementos de espesor cero no existe una distancia entre las superficies y, por lo tanto, la rigidez de contacto de estos elementos no depende de su espesor. Es por ello que se ha determinado que esta t´ ecnica es la m´ as adecuada para simular grietas en un medio homog´ eneo cuando se usa el M´ etodo de Elementos Finitos (FEM por sus siglas en ingl´ es). 3. RESULTADOS Para analizar la capacidad de predicci´ on del PMTE-SC se ha utilizado un problema cl´ asico de Griffith, donde se modela una placa con altura 2L=80 mm, anchura 2W=60 mm y una grieta inicial de longitud 2acentra- da en el eje x. Como se muestra en la Fig. 1, debido a la naturaleza sim´ etrica del problema y en aras de la eficiencia computacional, solo se modela la mitad derecha del esp´ ecimen. Adem´ as de las condiciones de contorno ilustradas en la figura, se asumen las condiciones de deformaci´ on plana y se usa un sistema de coordenadas cartesiano (x,y). La placa analizada esta hecha de Zirconia con las siguientes propiedades mec´ anicas: el m´ odulo de Young E=210 GPa, el coeficiente de Poisson ν=0.22, la resistencia a tracci´ on del material σc=564.19 MPa y la energ´ ıa de fractura Gc=476 Jm−2. Para realizar el mallado de la placa se ha utilizado la longitud caracter´ ıstica del material que se puede determinar a partir de la siguiente relaci´ on, lc=1 πKc σu2que da como resultado 0.1 mm para el material estudiado. Cabe destacar que el avance finito de la grieta obtenido usando el criterio acoplado suele ser en torno a esta longitud caracter´ ıstica o su fracci´ on, por lo que se debe definir una malla (suficientemente fina) que pueda capturar con precisi´ on esta longitud de avance de la grieta. En este trabajo se ha usado una malla uniforme de 300 elementos a lo largo del avance potencial de la grieta ∆as, determinado por el criterio tensional. Cabe mencionar que se ha comprobado que con una malla mucho mas basta se obtienen resultados similares. Se ha utilizado el elemento cuadr´ atico CPE8 de la librer´ ıa de ABAQUS. Para simular el inicio de la fractura con una longitud finita, ABAQUS realiza un an´ alisis FEM en el primer m´ odulo del algoritmo, es decir, el m´ odulo del criterio tensional, lo que da como resultado la m´ axima longitud inicial de la grieta ∆aspropagando desde la grieta inicial a, tal y como muestra la Fig. 1. Posteriormente, se incluye una Interfaz de Muelles El´ astica Lineal (LESI) a lo largo de las superficies producidas por la presencia de ∆as. El comportamiento de estos muelles se incluye a trav´ es del m´ odulo UINTER de ABAQUS. Posteriormente, el m´ odulo que representa el criterio energ´ etico predice la longitud de grieta resultante, ∆ac, utilizando el algoritmo AMA. Como consecuencia del modelado de la posible grieta como una distribuci´ on de muelles entre las dos superficies nuevas, las mallas usadas en los m´ odulos de los criterios tensional y energ´ etico son diferentes. Sin embargo, dentro del modulo del criterio energ´ etico no existen diferencias en las mallas de las diferentes configuraciones de avance de grieta, ya que las longitudes de las grietas en las iteraciones de AMA se pueden alterar ((activando)) y((desactivando)) los muelles presentes a lo largo de las superficies de la grieta. En consecuencia, mantener la misma topolog´ ıa de malla permite la reducci´ on de errores num´ ericos en los c´ alculos de incremento de energ´ ıa potencial. Δas a W L y x σ σ Figura 1: Ilustraci´on del un modelo de grieta en el centro con las dimensiones y las condiciones de contorno utilizadas para la simulaci´on num´erica. En la Fig. 2 se analizan los resultados obtenidos con el algoritmo descrito. En el eje de abscisas se representa la longitud de la grieta a medida que esta propaga desde su longitud inicial, y en el eje de ordenadas se representa la tensi´ on remota que produce un incremento de grieta. En las iteraciones del algoritmo que implementa el PMTESC se analizan varios avances potenciales de grietas que satisfacen el criterio tensional para una carga remota determinada σ, pero que no satisfacen el criterio energ´ etico. Algunos ejemplos de estos casos est´ an representados Revista de Mecánica de la Fractura Vol.5 (2023) 227 por las cruces rojas en la Fig. 2, el algoritmo no admite los correspondientes avances potenciales de grietas como aceptables, si no que aumenta la carga en el siguiente paso para que el sistema pueda almacenar m´ as energ´ ıa potencial y as´ ı llegar a satisfacer ambos criterios, tal y como se se˜ nala en la Fig. 2 como puntos verdes. Los resultados num´ ericos obtenidos con el PMTE-SC son comparados con una curva anal´ ıtica obtenida a partir del CCFFM cl´ asico (basado en la formulaci´ on presentada en [6]), se observa una muy buena coincidencia. El algoritmo puede predecir que la grieta crecer´ a inestablemente bajo un control de carga remota fija, tal y como se observa en los experimentos. Figura 2: Comparaci´on entre las soluciones anal´ıticas a trav´es del m´etodo cl´asico del CCFFM y los resultados num´ericos del nuevo m´etodo del PMTE-SC. 4. CONCLUSIONES En este trabajo, se ha implementado un nuevo m´ etodo, propuesto en [14], para caracterizar el inicio y la propagaci´ on de grietas basado en el Criterio Acoplado de la Mec´ anica de Fractura Finita (CCFFM), para el crecimiento de grietas en un medio homog´ eneo. Este nuevo m´ etodo permite modelar la aparici´ on simult´ anea de varias grietas en un mismo problema. El algoritmo computacional implementado permite simular un proceso de fractura en problemas planos basado en la nueva formulaci´ on del PMTE-SC, utilizando la subrutina UINTER de ABAQUS y controlada a partir de unos scripts en Python. Adem´ as, se introdujo un modelo de grieta a trav´ es de una distribuci´ on continua de muelles lineales entre las superficies de las posibles grietas, el cual permite que la minimizaci´ on de la energ´ ıa total se simplifique sustancialmente, ya que el funcional de energ´ ıa total (a priori no convexo) se puede considerar como un funcional convexo por separado en t´ erminos de las variables de desplazamiento y de da˜ no. Para evaluar la fiabilidad del algoritmo se ha analizado el problema cl´ asico de una grieta centrada en una placa infinita sometida a tracci´ on, desde un punto de vista num´ erico (con el nuevo enfoque del PMTE-SC) y anal´ ıtico (con el enfoque cl´ asico del CCFFM), arrojando resultados id´ enticos considerando tambi´ en tama˜ nos muy peque˜ nos de la grieta inicial. Estos resultados abren una nueva v´ ıa para estudiar problemas m´ as complejos de iniciaci´ on y crecimiento de m´ ultiples grietas. Sin embargo, se necesita una investigaci´ on m´ as profunda para optimizar el algoritmo, ya que es necesaria la mejora del m´ odulo del criterio tensional para poder ser utilizado en casos de m´ ultiples grietas potenciales. AGRADECIMIENTOS Este proyecto ha recibido financiaci´ on del Programa de Investigaci´ on e Innovaci´ on Horizonte 2020 de la Uni´ on Europea bajo el Acuerdo de Subvenci´ on Marie Skłodowska Curie No. 861061. – Proyecto NewFrac ITN. Adem´ as, tambi´ en ha sido financiado por el Ministerio de Ciencia e Innovaci´ on (PID2021-123325OB-I00), Consejer´ ıa de Econom´ ıa y Conocimiento de la Junta de Andaluc´ ıa (Proyecto P18-FR-1928 y Contrato US-1266016- Programa Operativo FEDER Andaluc´ ıa 2014-2020), Junta de Andaluc´ ıa y Fondo Social Europeo (Acciones de transferencia del conocimiento AT17-5908-USE). REFERENCIAS [1] G.I. Barenblatt. The mathematical theory of equilibrium cracks in brittle fracture. Advances in applied mechanics, 7:55–129, 1962. [2] D.S. Dugdale. Yielding of steel sheets containing slits. Journal of the Mechanics and Physics of Solids, 8(2):100–104, 1960. [3] A. Needleman and V. Tvergaard. Finite element analysis of localization in plasticity, In: Finite Elements: Special Problems in Solid Mechanics (Eds. J.T. Oden and G.F. Carey), pages 94–267. Prentice- Hall, 1984. [4] A. Needleman. Some issues in cohesive surface modeling. Procedia IUTAM, 10:221–246, 2014. [5] G.A Francfort and J.-J. Marigo. Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids, 46(8):1319–1342, 1998. [6] D. Leguillon. Strength or toughness? a criterion for crack onset at a notch. European Journal of Mechanics-A/Solids, 21(1):61–72, 2002. [7] D. Leguillon, E. Martin, O. ˇ Seveˇ cek, and R. Bermejo. Application of the coupled stress-energy criterion to predict the fracture behaviour of layered ceramics designed with internal compressive stresses. European Journal of Mechanics-A/Solids, 54:94– 104, 2015. Revista de Mecánica de la Fractura Vol.5 (2023) 228 [8] I.G. Garc´ ıa, B.J. Carter, A.R. Ingraffea, and V. Mantiˇ c. A numerical study of transverse cracking in cross-ply laminates by 3d finite fracture mechanics. Composites Part B: Engineering, 95:475–487, 2016. [9] A. Doitrand, C. Fagiano, N. Carrere, V. Chiarutti- ni, and M. Hirsekorn. Damage onset modeling in woven composites based on a coupled stress and energy criterion. Engineering Fracture Mechanics, 169:189–200, 2017. [10] J. Reinoso, A. Arteiro, M. Paggi, and P.P. Camanho. Strength prediction of notched thin ply laminates using finite fracture mechanics and the phase field approach. Composites Science and Technology, 150:205–216, 2017. [11] J. Li and D. Leguillon. Finite element implementation of the coupled criterion for numerical simulations of crack initiation and propagation in brittle materials. Theoretical and Applied Fracture Mechanics, 93:105–115, 2018. [12] D. Leguillon, E. Martin, O. ˇ Seveˇ cek, and R. Bermejo. What is the tensile strength of a ceramic to be used in numerical models for predicting crack initiation? International Journal of Fracture, 212:89–103, 2018. [13] E. Martin, D. Leguillon, O. ˇ Seveˇ cek, and R. Bermejo. Understanding the tensile strength of ceramics in the presence of small critical flaws. Engineering Fracture Mechanics, 201:167–175, 2018. [14] V. Mantiˇ c. Prediction of initiation and growth of cracks in composites. coupled stress and energy criterion of the finite fracture mechanics (Keynote lecture). In ECCM-16th Europ. Conf. on Composite Mater, pages 1–16, 2014. [15] V. Mantiˇ c. Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. application of a coupled stress and energy criterion. International journal of Solids and Structures, 46(6):1287–1304, 2009. [16] I.G. Garc´ ıa, M. Paggi, and V. Mantiˇ c. Fiber-size effects on the onset of fiber–matrix debonding under transverse tension: A comparison between cohesive zone and finite fracture mechanics models. Engineering Fracture Mechanics, 115:96–110, 2014. Revista de Mecánica de la Fractura Vol.5 (2023) 229 Revista de Mecánica de la Fractura Vol.5 (2023) 230