Full text
DIFERENTES ENFOQUES PARA LA FRACTURA D ´ UCTIL BASADOS EN LA METODOLOG´ IA PHASE-FIELD: APLICACI ´ ON A MODELOS DE PLASTICIDAD ASOCIADA Y NO ASOCIADA P. Olivares-Rodr´ ıguezab*, A. Quintanas-Corominascd, E. Mart´ ınez-Pa˜ nedad, J. Reinosoa, P. Durandb aGrupo de Elasticidad y Resistencia de Materiales, E.T.S. de Ingenier´ ıa, Universidad de Sevilla Camino de los Descubrimientos s/n, 41092 Sevilla, Espa˜ na bDepartamento de Estructuras de Edificaci´ on e Ingenier´ ıa del Terreno, E.T.S. de Arquitectura, Universidad de Sevilla Av. Reina Mercedes 2, 41012 Sevilla, Espa˜ na cAMADE, Escuela Polit´ ecnica, Universidad de Girona, Campus Montilivi s/n, 17073 Girona, Espa˜ na dImperial College, Exhibition Rd, South Kensington, SW7 2BX Londres, Reino Unido *Persona de contacto: poliv[email protected] RESUMEN Los procesos de fractura son una de las principales causas de fallo en diferentes materiales y estructuras en diferentes aplicaciones en Ingenier´ ıa. Por ello, a lo largo de las ´ ultimas d´ ecadas, su estudio ha estado enfocado en tareas como el desarrollo de criterios de prevenci´ on en fase de proyecci´ on de estructuras o la modelizaci´ on del comportamiento de materiales ante situaciones (estados de carga) que conlleven la iniciaci´ on y el desarrollo de grietas. En el caso concreto de la Ingenier´ ıa Civil, los materiales de naturaleza geomec´ anica (suelos y rocas) sufren a su vez este tipo de procesos de car´ acter inel´ astico. En el contexto de la Mec´ anica Computacional, se han propuesto diversas metodolog´ ıas de an´ alisis de procesos de fractura basadas en la Mec´ anica de la Fractura El´ astica Lineal, modelos de grieta cohesiva, entre otros muchos. Recientemente, una de las t´ ecnicas con mayor impacto en la comunidad cient´ ıfica ha sido la denominada Phase Field Approach of Fracture (PFAF), basada en la visi´ on energ´ etica originalmente propuesta por Griffith incluida en un enfoque variacional desde el punto de vista matem´ atico. En este trabajo, se desarrolla la extensi´ on de la PFAF para la modelizaci´ on de la fractura conjuntamente con leyes materiales no lineales siguiendo diversas formulaciones propuestas en la literatura y aplicadas a modelos de fractura con flujo asociativo y no asociativo. El trabajo concierne un profundo an´ alisis de las diferentes opciones de modelizaci´ on de fuerzas precursoras de fractura en el contexto de PFAF as´ ı como su correspondiente implementaci´ on num´ erica empleando el m´ etodo de los elementos finitos. PALABRAS CLAVE: mec´ anica computacional, plasticidad, phase field, elementos finitos, mec´ anica de la fractura ABSTRACT Fracture events are one of the most common causes of failure that occur in different materials used in different engineering applications. Therefore, over the last decades, the comprehensive understanding of these phenomena has been received a great deal of attention with focus on the prevention criteria development in the designing and dimesnioning stages of structures or the modelling of the material behaviour in situations (related to load states) that lead the development of cracks. In Civil Engineering applications, geomaterials (soils and rock) are prone to develop this type of inelastic process. Within the context of Computational Mechanics, several fracture methodologies based on Linear Elastic Fracture Mechanics, cohesive crack models, among many others, have been proposed so far. Recently, one of the techniques with the greatest impact on the scientific community has been the so-called Phase Field Approach of Fracture (PFAF), based on the energetic vision originally proposed by Griffith included in a variational approach from a mathematical point of view. In this work, the extension of the PFAF for material nonlineal effects prior fracture occurence is developed following several formulations proposed in the literature and applied to fracture models with associative and non-associative flow. The work concerns with an in-depth analysis of the different options for modelling fracture driving forces in the context of PFAF as well as their corresponding numerical implementation using the finite element method. KEYWORDS: computational mechanics, plasticity, phase field, finite element method, fracture mechanics La motivaci´ on que persigue este trabajo desde el punto de vista computacional se basa en desarrollar un modelo para la aplicaci´ on de diversas formulaciones de fractura siguiendo la metodolog´ ıa PFAF acoplada con una ley constitutiva de car´ acter elastopl´ astico. La mayor´ ıa de modelos ideados para modelar procesos de fractura seg´ un la aproximaci´ on de tipo variacional PFAF siguen el modelo can´ onico propuesto por Francfort y Marigo [1], para obtener una versi´ on regularizada. Dicho procedimiento se conoce como la formulaci´ on tipo phase- Revista de Mecánica de la Fractura Vol.3 2022 101
Figura 1: Situaci´on de grieta discreta y grieta regularizada seg´un PFAF. field. Esta formulaci´ on se basa en regularizar la superficie discontinua de la grieta a trav´ es de una variable auxiliar de tipo difuso, por la que se obtiene un campo continuo de la variable da˜ no evaluable en cada punto del dominio, ve´ ase [2]. En la Figura 1 se puede observar como la grieta (representada por el campo discreto Γ) se aproxima a trav´ es de la variable phase-field, que act´ ua en la zona delimitada por la longitud caracter´ ıstica lc(par´ ametro de control). La variable de phase-field se define como d:Ω×[0,t]→ [0,1]. Es decir, dinterpola entre el estado intacto del material d=0y el estado totalmente da˜ nado d=1, donde los estados intermedios suponen una p´ erdida progresiva de rigidez. Por otro lado, los fen´ omenos no lineales de tipo material se tienen en cuenta en este trabajo a partir de un modelo elasto-pl´ astico [3] definido por: (i) criterio de plastificaci´ on, (ii) regla de flujo y (iii) la correspondiente ley de endurecimiento. En este trabajo se ha implementado el modelo de Mohr-Coulomb con una ley de endurecimien- to afectando a la variable de cohesi´ on, siguiendo las directrices de [4]. El car´ acter asociativo o no asociativo del modelo pl´ astico reside en la diferencia entre los valores del ´ angulo de rozamiento ϕy el ´ angulo de dilatancia ψ, par´ ametros de la funci´ on de plastificaci´ on y el potencial pl´ astico respectivamente. Por tanto, las variables independientes del problema son el campo de desplazamientos u, el campo de deformaciones pl´ asticas εpy la variable phase-field d.Estas variables entran dentro del montaje del modelo de car´ acter acoplado a trav´ es de la formulaci´ on de la densidad de energ´ ıa interna total, propia de formulaciones de tipo variacional. La densidad de energ´ ıa presenta la siguiente forma [5]: Wint(ε,εp,d,∇d, α, ∇α) :=ψ(ε, εp,d, α, ∇α)+ ∆(d,∇d, α)= =g(d)ψe(ε−εp) E +h(d)ψp(α, ∇α) H + ∆f(d,∇d) F + +p(d)∆p(α) P (1) En la que αes la variable de endurecimiento. En un primer desarrollo la energ´ ıa se divide en dos t´ erminos: la energ´ ıa libre ψy el trabajo disipativo ∆, que puede ser local o no local para el campo de fractura en funci´ on de si se considera el t´ ermino de ∇d. Pasando a un segundo nivel contamos con cuatro t´ erminos, cuyo desarrollo puede consultarse en [5], enfoc´ andonos en este documento en el t´ ermino de fractura F: E: T´ ermino correspondiente a la energ´ ıa el´ astica degradada. Dicha energ´ ıa se degrada a partir de una funci´ on g(d) cuyo valor disminuye desde 1 a 0 a medida que el valor de daumenta. H: T´ ermino correspondiente a la energ´ ıa pl´ astica degradada. De nuevo existe una funci´ on h(d) que cumple el mismo prop´ osito y condiciones que g(d). F: T´ ermino correspondiente al trabajo disipativo de da˜ no. En el contexto que se trata toma la forma: ∆f(d,∇d) :=Gc cω ω(d) lc +lc|∇d|2!(2) Donde Gces la tenacidad a fractura del material, y ω(d) y cωson los coeficientes propios del funcional de Ambrosio-Tortorelli. En el caso m´ as com´ un (conocido como AT-2) se considera la funci´ on de densidad de superficie de grieta definida en [2]: ∆f(d,∇d) :=Gc 2 d2 lc +lc|∇d|2!(3) Puede comprobarse que utilizando este tipo de formulaci´ on se da la Γ-Convergencia al modelo de fractura fr´ agil de Griffith llevando a 0 la longitud caracter´ ıstica lc, v´ ease [2, 6]. Sin embargo, tiene el inconveniente de que el da˜ no empieza a entrar en juego desde que se inicia el estado de carga. Esta desventaja se suele suplir mediante el uso de thresholds, por los que se impone un umbral para el comienzo del desarrollo del da˜ no [7, 8]. Tambi´ en existen otros autores que emplean un funcional tipo AT-1, que desarrolla un tramo lineal previo al da˜ no, aunque es necesario realizar ajustes para mantener la condici´ on de irreversibilidad [9]. Revista de Mecánica de la Fractura Vol.3 2022 102
P: T´ ermino correspondiente al trabajo disipativo pl´ astico, degradado a partir de la funci´ on p(d). Con estos t´ erminos, se ensambla la densidad de energ´ ıa y se integra en el dominio de estudio para despu´ es obtener la formulaci´ on energ´ etica propia del m´ etodo de los elementos finitos. La organizaci´ on del art´ ıculo se estructura como sigue: En primer lugar, se describe y se detalla el modelo computacional desarrollado (Secs. 1 y 2), para posteriormente extraer una serie de resultados de car´ acter cualitativo a trav´ es de un ejemplo sencillo para una serie de metodolog´ ıas reconocidas (Sec. 3). Finalmente se resumen las principales conclusiones en la Sec. 4. 1. MODELO COMPUTACIONAL El modelo computacional desarrollado se basa en el m´ etodo de los elementos finitos y fue desarrollado en ABAQUS. La forma de resolver el problema acoplado pasa por implementar un procedimiento tipo Staggered Algorithm [10]. En ABAQUS, se trabajar´ a con los entornos personalizables tipo UEL y UMAT. Una UEL (User Element) es un tipo de interfaz que permite implementar elementos finitos seg´ un una metodolog´ ıa propuesta por el usuario. Las aportaciones de este entorno son las contribuciones al vector de residuos y a la matriz de rigidez del sistema de ecuaciones del MEF. Por otro lado, la UMAT se encarga de definir el modelo de comportamiento del material a trav´ es de la actualizaci´ on del estado de tensiones y de la matriz de rigidez en base al modelo pl´ astico implementado. Como variables generales de entrada se detalla el valor estimado del campo de desplazamientos y de da˜ no, que deben converger en cada paso de tiempo por el car´ acter desacoplado del algoritmo empleado. Dentro del procedimiento, uno de los aspectos m´ as destacables es el c´ alculo de las energ´ ıas. Para la energ´ ıa el´ astica (t´ ermino Esin degradar) la formulaci´ on en materiales is´ otropos toma la forma [11]: Ψe=λ 2(tr εe)2+µtr h(εe)2i(4) Se puede observar que la Ec.(4) es combinaci´ on lineal de los invariantes del tensor de deformaciones el´ astico. En este trabajo se utilizar´ a la descomposici´ on de Miehe [7], que se basa en la diferenciaci´ on de la parte positiva (tracci´ on) y negativa (compresi´ on) haciendo uso de las deformaciones el´ asticas del sistema principal, que son los autovalores del tensor de deformaciones el´ astico (εe 1, εe 2, εe 3): Ψe±=λ 2⟨εe 1+εe 2+εe 3⟩2 ±+µ⟨εe 1⟩2 ±+⟨εe 2⟩2 ±+⟨εe 3⟩2 ±(5) El operador ⟨z⟩±son los par´ entesis de Macaulay, que representan la funci´ on a trozos ⟨z⟩±:=(z± |z|)/2. Por otro lado, la energ´ ıa disipativa pl´ astica Ψpd (t´ ermino H+Psin degradar) se calcula a trav´ es de la obtenci´ on del ´ area bajo la curva de tensi´ on-deformaci´ on pl´ astica, de forma que: Ψpd = Ψpd + n X i=1 σn+1 i+σn i 2(εn−1 pi−εn pi) (6) Siendo nel n´ umero de componentes del pseudovector de deformaciones pl´ asticas. Tal y como se muestra en la figura 2 para una componente gen´ erica, la energ´ ıa se calcula mediante la discretizaci´ on de la integral o suma discreta de ´ areas sucesivas, a˜ nadiendo en cada paso la encerrada entre los valores de deformaci´ on pl´ astica del paso anterior y del actual. Figura 2: Obtenci´on de la energ´ıa disipativa pl´astica. Se observa que este tipo de disipaci´ on se calcula en funci´ on de las tensiones efectivas, que no est´ an afectadas del campo da˜ no. Se emplean este tipo de tensiones, y no las nominales (afectadas del campo de da˜ no), ya que la energ´ ıa se degrada posteriormente a trav´ es de la funci´ on de degradaci´ on. Tras la obtenci´ on de todas las variables necesarias, se a˜ naden como contribuciones a los vectores de residuos y matrices de rigidez, a trav´ es de las expresiones definidas del procedimiento tipo Staggered Algorithm en [10]. 2. AN ´ ALISIS DE METODOLOG´ IAS En esta secci´ on se har´ a especial enf´ asis en el c´ alculo del residuo y parte de la matriz de rigidez dependiente del campo de da˜ no d,rd iyKdd ij . Estas variables dependen, adem´ as de la discretizaci´ on del MEF realizada, del t´ ermino de Driving Force. Dicho t´ ermino tiene su origen en la ecuaci´ on de evoluci´ on local del da˜ no [8]: Gc lc [d−l2 c∆d] R. Geom´ etrica =2(1 −d)H Driving Force (7) Donde el crecimiento de grieta se basa en un equilibrio, a modo de competici´ on, entre el t´ ermino de la resistencia geom´ etrica y la Driving Force, gobernada por la variable Revista de Mecánica de la Fractura Vol.3 2022 103
hist´ orica H, responsable de la condici´ on de irreversibilidad del da˜ no en modelos AT-2. Diferentes autores, en aras de desarrollar un modelo de fractura d´ uctil que acople los campos de plasticidad y da˜ no, han propuesto y desarrollado diversas formulaciones para el t´ ermino de Driving Force. En este trabajo se estudiar´ a la derivaci´ on de las metodolog´ ıas de Borden [12] y Miehe [13], que son los que m´ as versatilidad ofrecen a la hora de modelar este tipo de comportamiento en materiales. En ambos casos, se parte de la ecuaci´ on de evoluci´ on de la superficie de grieta (7) realizando una serie de modificaciones a trav´ es de par´ ametros y thresholds. A continuaci´ on, se muestra como obtener las componentes del vector de residuos y matriz de rigidez asociados al da˜ no para cada metodolog´ ıa. Metodolog´ ıa propuesta por Borden et al.: En este caso, la ecuaci´ on de evoluci´ on sufre cambios en la definici´ on de Hyg(d), partiendo de la Ec. (3), quedando de la siguiente forma: Gc lcd−l2 c∆d=−g′(d)hβeH+βp⟨Ψpd −W0⟩+i (8) Donde: H=m´ ax s∈[0,t]Ψe(εe) g(d)=(1 −d)2(1 +d(2 −m)),m>0 (9) Por facilitar la comprensi´ on en el desarrollo, se ha obviado intencionadamente la descomposici´ on en parte positiva y negativa de la energ´ ıa el´ astica en ambas metodolog´ ıas, teniendo solo que afectar la funci´ on de degradaci´ on a la parte positiva. En base a la discusi´ on ofrecida en [5], Borden y coautores utilizan una funci´ on de degradaci´ on c´ ubica con AT-2 para poder generar un tramo lineal previo a alcanzar la deformaci´ on correspondiente a la tensi´ on pico (como comparaci´ on al comportamien- to de la de car´ acter cuadr´ atico consultar [12]) y as´ ı suplir el defecto del comienzo inmediato del da˜ no con la puesta en carga propia de la combinaci´ on {g(d) cuadr´ atica+AT-2}. Para conseguir mantener el mismo nivel de deformaciones al alcanzar la ten- si´ on pico que en el caso cuadr´ atico y mantener la semejanza con un modelo AT-2, el par´ ametro mse debe ajustar. Dicho ajuste se corresponde a un valor mucho menor que la unidad, pero sin alcanzar su nulidad ya que no se podr´ ıa producir la nucleaci´ on de defectos para el estado intacto del material. En [12], siguiendo diversas modelizaciones, se opta por un valor de m=10−4. Por otro lado, Borden separa la Driving Force en dos t´ erminos afectados por unos coeficientes de ponderaci´ on βeyβp. El primer t´ ermino afecta a la energ´ ıa el´ astica y el segundo a la disipaci´ on pl´ astica controlada adem´ as por un threshold W0. Los coeficientes de ponderaci´ on modelan la contribuci´ on de cada t´ ermino al crecimiento de grieta, que es el fin ´ ultimo de la Driving Force. Para obtener la componente asociada al da˜ no del vector de residuos y matriz de rigidez es necesario utilizar la discretizaci´ on propia del MEF: d=Nid d di ∇d=Bid d di (10) Por la cual, sustituyendo en (8) y derivando se obtiene: rd i=ZB("Gc lc d+g′(d)βeH+βp⟨Ψpd− −W0⟩+)]Ni+GclcBT i∇dodV Kdd i j =∂rd i ∂dj =ZB("Gc lc +g′′(d)βeH+βp⟨Ψpd− −W0⟩+)]NiNj+GclcBT iBjodV (11) Metodolog´ ıa propuesta Miehe y coautores: Por su parte, Miehe et al. [8] realizan la modificaci´ on del t´ ermino Fde (1). En dicho t´ ermino establece una extensi´ on del modelo AT-2: ∆f(d,∇d) :=2ωc ξ d2 lc +lc|∇d|2!+(1 −g(d))ωc (12) En el que a˜ nade un t´ ermino adicional controlado por la energ´ ıa cr´ ıtica de fractura ωcque act´ ua a modo de threshold. Desarrollando la energ´ ıa total del sistema y derivando con respecto a la variable da˜ no para obtener la ecuaci´ on de evoluci´ on queda: 2ωc ξd−l2 c∆d=−g′(d)H(13) Donde: H=m´ ax s∈[0,t]h⟨Ψe+ Ψpd −ωc⟩i g(d)=(1−d)2 (14) Como puede observarse, esta formulaci´ on introduce el threshold de manera natural en el t´ ermino de energ´ ıa de fractura, que deriva en la ecuaci´ on de evoluci´ on, a diferencia de otros autores que optan por establecerlo en los t´ erminos de energ´ ıa el´ astica o pl´ astica. De esta forma, consigue asemejarse a un modelo AT-1 con un tramo previo lineal sin da˜ no apreciable, pero con las ventajas propias de un modelo AT-2. Adem´ as, emplea un par´ ametro ξpara controlar la ca´ ıda de la fase postcr´ ıtica de da˜ no. De forma an´ aloga al proceso realizado en Borden et al. [12], se obtienen las contribuciones al sistema Revista de Mecánica de la Fractura Vol.3 2022 104
0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 0.25 0.27 25 28 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 0 0.2 0.4 0.6 0.8 1 0 0.002 0.004 0.006 0.008 0.01 0 5 10 15 20 25 30 0 0.2 0.4 0.6 0.8 1 0 0.002 0.004 0.006 0.008 0.01 0 5 10 15 20 25 30 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 0.25 0.27 25 28 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 30 0.25 0.27 25 28 Figura 3: Resultados de las curvas de comportamiento (l´ınea continua) y evoluci´on del da˜no (l´ınea discontinua). del MEF: rd i=ZB("2ωc ξd+g′(d)H#Ni+ +2ωcl2 c ξBT i∇d)dV Kdd i j =∂rd i ∂dj =ZB("2ωc ξ+g′′(d)H#NiNj+ +2ωcl2 c ξBT iBj)dV (15) 3. RESULTADOS Para verificar las metodolog´ ıas presentadas en la secci´ on anterior con el modelo computacional desarrollado, se comprobar´ a de forma cualitativa su efectividad a trav´ es de un ejemplo sencillo. Dicho ejemplo supone la modelizaci´ on de un elemento sometido a tracci´ on en condiciones de deformaci´ on plana. El elemento utilizado es de tipo bilineal de 4 nodos, con 4 puntos de integraci´ on. Dentro de la modelizaci´ on de la plasticidad a partir del modelo de Mohr-Coulomb se propone la siguiente ley de endurecimiento de car´ acter lineal para la cohesi´ on, en funci´ on de la deformaci´ on pl´ astica acumulada ¯εp(variable de endurecimiento): c(¯εp)=c0+H¯εp(16) Las propiedades empleadas en la modelizaci´ on son las recogidas en la Tabla 1. Estas propiedades modelizan de forma cualitativa un suelo tipo arcilloso en condiciones no drenadas (simulaci´ on a corto plazo), material caracter´ ıstico de aplicaci´ on del modelo Mohr-Coloumb, con un m´ odulo de endurecimiento H negativo ya que la tracci´ on va asociada a una p´ erdida de resistencia en suelos. Material E[kPa] νGc[kPa·m] lc[m] Arcilla 8.0·1030.2 0.5 0.2 Material c0[kPa] H [kPa] ϕ[º] ψ[º] Arcilla 15.0 -10.0 1.0 1.0 Tabla 1: Propiedades y par´ametros utilizados para las metodolog´ıas de Borden y Miehe. A continuaci´ on, a la vista de los resultados ilustrados en la Figura 3 puede concluirse: En la comparativa I de Borden se denota la diferencia entre la aproximaci´ on de Borden et al. (m = 10−4) y un AT-2 (m =2) sin threshold, con referencia al estado no da˜ nado. Se observa que mientras que la aproximaci´ on de Borden et al. conserva una fase no da˜ nada antes de darse la ca´ ıda s´ ubita de resistencia, en AT-2 la degradaci´ on comienza desde la fase el´ astica. Adem´ as, gracias al zoom en la gr´ afica se puede apreciar una transici´ on m´ as suave a la ca´ ıda en Borden et al. que en AT-2. En la comparativa II de Borden se observa el efec- to del threshold en la parte pl´ astica de la Driving Force. Como era de esperar, a mayor valor del threshold m´ as se retrasa el inicio de la fase de da˜ no en el elemento. De forma an´ aloga al caso anterior, en la comparativa I de Miehe se comprueba el efecto de la energ´ ıa Revista de Mecánica de la Fractura Vol.3 2022 105
cr´ ıtica ωc. A mayor valor de este par´ ametro mayor es la energ´ ıa necesaria para desencadenar el proceso de da˜ no. Por ´ ultimo, en la comparativa II de Miehe observamos las diferentes formas de ca´ ıda propias de la fase post-cr´ ıtica debido al par´ ametro ξ. A medida que aumenta el valor de este par´ ametro la pendiente va creciendo, asemej´ andose a una fase de da˜ no sin mucho recorrido en deformaciones propia de la fractura de tipo fr´ agil. 4. CONCLUSIONES El presente trabajo muestra la capacidad de las distintas metodolog´ ıas descritas empleando el modelo computacional desarrollado. En este caso, el car´ acter asociado o no asociado de la plasticidad viene dada por una aplicaci´ on sencilla del modelo de Mohr-Coulomb. Gracias al problema planteado, se ha podido comprobar como cada par´ ametro influye en la soluci´ on. Con este conocimiento se podr´ ıa trabajar con el ajuste de resultados experimentales. Sin embargo, cabe destacar que detr´ as de cada metodolog´ ıa y de la obtenci´ on del valor de cada par´ ametro deber´ ıa existir un trasfondo f´ ısico que le aporte rigor cient´ ıfico a la modelizaci´ on efectuada. Finalmente, se denota la versatilidad del modelo computacional para incluir otro tipo de leyes de comportamiento distintas a trav´ es del desarrollo de UMAT espec´ ıficas, que pueden acoplarse f´ acilmente al entorno de la UEL implementado en ABAQUS.Adem´ as, una vez verificado los resultados para un elemento, los problemas planteados pueden asemejarse a geometr´ ıas y condiciones m´ as realistas. AGRADECIMIENTOS JR, PD y POR agradecen su apoyo al Ministerio de Ciencia e Innovaci´ on y a la Agencia Estatal de Investigaci´ on a trav´ es de la subvenci´ on del proyecto PID2019-109723GB- I00. Adem´ as, POR agradece el apoyo de la Universidad de Sevilla a trav´ es del contrato predoctoral otorgado en la convocatoria VIPPIT-2020-II.2. Por su parte, AQC agradece su apoyo al Ministerio de Universidades del Gobierno Espa˜ nol y a los fondos del programa Next Generation EU de la Uni´ on Europea. REFERENCIAS [1] G. A. Francfort and J. J. Marigo. Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids, 46, 1998. [2] C. Miehe, F. Welschinger, and M. Hofacker. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations. International Journal for Numerical Methods in Engineering, 83, 2010. [3] Jacob Lubliner. Plasticity theory: Revised edition. Journal of Applied Mechanics, 59, 2006. [4] E. A. De Souza Neto, D. Peri´ c, and D. R.J. Owen. Computational Methods for Plasticity: Theory and Applications. 2008. [5] R. Alessi, M. Ambati, T. Gerasimov, S. Vidoli, and L. De Lorenzis. Comparison of phase-field models of fracture coupled with plasticity. Computational Methods in Applied Sciences, 46, 2018. [6] Gilles A. Francfort, Blaise Bourdin, and Jean Jacques Marigo. The variational approach to fracture. Journal of Elasticity, 91:5–148, 4 2008. [7] Christian Miehe, Lisa Marie Sch¨ anzel, and Heike Ulmer. Phase field modeling of fracture in multiphysics problems. part i. balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids. Computer Methods in Applied Mechanics and Engineering, 294, 2015. [8] C. Miehe, M. Hofacker, L. M. Sch¨ anzel, and F. Aldakheel. Phase field modeling of fracture in multiphysics problems. part ii. coupled brittle-to-ductile failure criteria and crack propagation in thermoelastic-plastic solids. Computer Methods in Applied Mechanics and Engineering, 294, 2015. [9] Roberto Alessi, Jean Jacques Marigo, Corrado Maurini, and Stefano Vidoli. Coupling damage and plasticity for a phase-field regularisation of brittle, cohesive and ductile fracture: One-dimensional examples. International Journal of Mechanical Sciences, 149, 2018. [10] Jianguang Fang, Chengqing Wu, Timon Rabczuk, Chi Wu, Conggan Ma, Guangyong Sun, and Qing Li. Phase field fracture in elasto-plastic solids: Abaqus implementation and case studies. Theoretical and Applied Fracture Mechanics, 103, 2019. [11] N. P.van Dijk, J. J. Espadas-Escalante, and P. Isaksson. Strain energy density decompositions in phasefield fracture theories for orthotropy and anisotropy. International Journal of Solids and Structures, 196- 197, 2020. [12] Michael J. Borden, Thomas J.R. Hughes, Chad M. Landis, Amin Anvari, and Isaac J. Lee. A phasefield formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects. Computer Methods in Applied Mechanics and Engineering, 312, 2016. [13] C. Miehe, S. Teichtmeister, and F. Aldakheel. Phase-field modelling of ductile fracture: A variational gradient-extended plasticity-damage theory and its micromorphic regularization. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 374, 2016. Revista de Mecánica de la Fractura Vol.3 2022 106