Diferentes aproximaciones numéricas para el acoplamiento de no linealidad material y fractura basada en la metodología phase-field
Abstract
L os procesos de fractura son una de las principales causas de fallo que ocurren en diferentes materiales en aplicaciones en Ingeniería. Por ello, a lo largo de las últimas décadas, su estudio ha estado enfocado en tareas como el desarrollo de criterios de prevención en fase de proyección de estructuras o la modelización del comportamiento de materiales ante situaciones (estados de carga) que conlleven la iniciación y el desarrollo de grietas. En el caso concreto de la Ingeniería Civil, los materiales de naturaleza geomecánica (suelos y rocas) sufren a su vez este tipo de procesos de carácter inelástico. En el contexto de la Mecánica Computacional, se han propuesto diversas metodologías de análisis de procesos de fractura basadas en la Mecánica de la Fractura Elástica Lineal (MFEL), modelos de grieta cohesiva, entre otros muchos. Recientemente, una de las técnicas con mayor impacto en la comunidad científica ha sido la denominada Phase Field Approach of Fracture (PFAF), que es la visión energética originalmente propuesta por Griffith incluida en un enfoque variacional desde el punto de vista matemático. En este trabajo, se desarrolla la extensión de la PFAF para la modelización 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álisis de las diferentes opciones de modelización de fuerzas precursoras de fractura en el contexto de PFAF así como su correspondiente implementación numérica empleando el método de los elementos finitos.
Full text
Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Máster Máster en Ingeniería de Caminos, Canales y Puertos Diferentes aproximaciones numéricas para el acoplamiento de no linealidad material y fractura basada en la metodología phase-field Autor: Pablo Olivares Rodríguez Tutor: José Antonio Reinoso Cuevas Tutor Externo: Adrià Quintanas Corominas Dpto. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022
Trabajo Fin de Máster Máster en Ingeniería de Caminos, Canales y Puertos Diferentes aproximaciones numéricas para el acoplamiento de no linealidad material y fractura basada en la metodología phase-field Autor: Pablo Olivares Rodríguez Tutor: José Antonio Reinoso Cuevas Profesor Titular de Universidad Tutor Externo: Adrià Quintanas Corominas Doctor en Ingeniería Mecánica Dpto. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022
Trabajo Fin de Máster: Diferentes aproximaciones numéricas para el acoplamiento de no linealidad material y fractura basada en la metodología phase-field Autor: Pablo Olivares Rodríguez Tutor: José Antonio Reinoso Cuevas Tutor Externo: Adrià Quintanas Corominas El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha:
Resumen Lo s procesos de fractura son una de las principales causas de fallo que ocurren en diferentes materiales en aplicaciones en Ingeniería. Por ello, a lo largo de las últimas décadas, su estudio ha estado enfocado en tareas como el desarrollo de criterios de prevención en fase de proyección de estructuras o la modelización del comportamiento de materiales ante situaciones (estados de carga) que conlleven la iniciación y el desarrollo de grietas. En el caso concreto de la Ingeniería Civil, los materiales de naturaleza geomecánica (suelos y rocas) sufren a su vez este tipo de procesos de carácter inelástico. En el contexto de la Mecánica Computacional, se han propuesto diversas metodologías de análisis de procesos de fractura basadas en la Mecánica de la Fractura Elástica Lineal (MFEL), modelos de grieta cohesiva, entre otros muchos. Recientemente, una de las técnicas con mayor impacto en la comunidad científica ha sido la denominada Phase Field Approach of Fracture (PFAF), que es la visión energética originalmente propuesta por Griffith incluida en un enfoque variacional desde el punto de vista matemático. En este trabajo, se desarrolla la extensión de la PFAF para la modelización 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álisis de las diferentes opciones de modelización de fuerzas precursoras de fractura en el contexto de PFAF así como su correspondiente implementación numérica empleando el método de los elementos finitos. PALABRAS CLAVE: mecánica computacional, plasticidad, phase field, elementos finitos, mecánica de la fractura I
Abstract Fr acture events are one of the most common causes of failure that occur in different materials used in 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 dimensioning stages of structures as well as on 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 (LEFM), cohesive crack models, among many others, have been proposed so far. Recently, one of the techniques with great impact on the scientific community is the so-called Phase Field Approach of Fracture (PFAF), that is 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 III
4Capítulo 1. Introducción Figura 1.2 Tipos de comportamiento en modelización computacional (extraído de [12]). tratar con modelos 1-D, tal y como se ilustra en la Figura 1.3. De esta forma, con vista en la Figura 1.2, se distingue: ▼ Modelo E-D (Rotura frágil): Se alcanza antes la tensión σd , sin en ningún momento alcanzar σp . Este caso se trata de un modelo puramente de daño, con disminución progresiva de la rigidez y sin formación apreciable de deformaciones plásticas. Encontramos ejemplos de implementación en [ 13 ] y [14]. ▼ Modelo E-D-PD (Rotura cohesiva): Este caso es similar al anterior con la salvedad de que tras alcanzar σd , existe un momento en el que se llega también a σp . A partir de este momento, se desarrollan deformaciones plásticas que van acopladas al proceso de daño, observándose en (b1) de la Figura 1.2
1.2 Revisión del estado del arte 5 Figura 1.3 Definición de las tensiones σdyσpa través de la curva de comportamiento. como a medida que se desarrolla el daño las deformaciones de carácter permanente aumentan. Se pueden encontrar ejemplos y desarrollo de técnicas numéricas como el Cohesive Zone Model (CZM) [15]. ▼ Modelo E-P-D (Rotura dúctil desacoplada): En este tipo de modelo se alcanza en primer lugar la tensión σp para posteriormente alcanzar σd , simulando el recorrido de tensiones en comportamientos tipo dúctil. Al ser de tipo desacoplado, se observa en (c1) de la Figura 1.2 que primero se desarrollan las deformaciones plásticas en el tramo plástico σ=σp<σd sin ningún tipo de disminución en la rigidez del material (pendiente del tramo de carga-descarga) y, después, se da el tramo de daño σ=σd donde ya existe una bajada notoria de la rigidez del material pero sin seguir desarrollando deformaciones permanentes. Ejemplos de este tipo de modelizaciones se encuentran en [16] y [17]. ▼ Modelo E-P-PD (Rotura dúctil acoplada): En contraste con el modelo anterior, en esta tipología de carácter acoplado, tras alcanzar la tensión σd , existe una evolución conjunta de plasticidad y daño por la cual se da de forma simultánea la disminución de la rigidez y la continuación de la formación de deformaciones permanentes. Este tipo de modelos cuentan con mayor validez experimental al simular de forma más realista el comportamiento de rotura tipo dúctil en materiales. Se encuentran ejemplos de aplicación en [18], [19] y [20]. La mayoría de modelos ideados para simular procesos de fractura según la aproximación del PFAF se basan en la aproximación de tipo variacional al modelo canónico propuesto por Francfort y Marigo [ 5 ], para obtener una versión regularizada. Dicho procedimiento se conoce como la formulación tipo phase-field. Esta formulación se basa en regularizar la superficie discontinua de la grieta a través de una variable auxiliar de tipo difuso, por la que se obtiene un campo continuo de la variable daño evaluable en cada punto del dominio, véase [ 21 ]. En la Figura 1.4 se puede observar como la grieta (representada por el campo discreto Γ ) se aproxima a través de la variable phase-field, que actúa en la zona delimitada por la longitud característica lc (parámetro de control) relacionada con σd. Figura 1.4 Situación de grieta discreta y grieta regularizada según PFAF.
6Capítulo 1. Introducción La variable phase-field se define como d:Ω×[0,t]−→[0,1] . Es decir, d interpola entre el estado intacto del material d=0 y el estado totalmente dañado d=1 donde los estados intermedios suponen una pérdida progresiva de rigidez. Por otro lado, existen otras variables que entran en juego en el problema de modelización de fractura, y que se recogen en la Tabla 1.1. Las variables independientes del problema son el campo de desplazamientos u, el campo de deformaciones plásticas ε ε εpy la variable phase-field d. Tabla 1.1 Variables del problema (Extraída de [22]). Tabla de definición de variables Símbolo Variable Tipo uDesplazamiento Observable ε ε εDeformación total ε ε εpDeformación plástica Interna αVariable de endurecimento (se adopta la deformación plástica acumulada ¯ εp) ∇xαGradiente de la variable de endurecimiento dVariable de daño Interna ∇xdGradiente de la variable daño Estas variables entran dentro del montaje del modelo a través de la formulación del potencial de energía interna total, propia de formulaciones de tipo variacional. La densidad de energía presenta la siguiente forma [22]: Wint (ε ε ε,ε ε εp,d,∇xd,α,∇xα):=ψ(ε ε ε,ε ε εp,d,α,∇xα)+∆(d,∇xd,α) = =g(d)ψe(ε ε ε−ε ε εp) E +h(d)ψp(α,∇xα) H +ˆ ψf(d,∇xd) F +p(d)ˆ ψp(α) P (1.1) En un primer desarrollo la energía se divide en dos términos: la energía libre ψ y el trabajo disipativo ∆ , que puede ser local o no local para el campo de fractura en función de si se considera el término de ∇xd . Pasando a un segundo nivel contamos con cuatro términos, que siguiendo el orden de sumandos, la nomenclatura y desarrollos de [22] son: ▼E : Término correspondiente a la energía elástica degradada. Dicha energía se degrada a partir de una función g(d) cuyo valor disminuye desde 1 a 0 a medida que el valor de d aumenta. Las condiciones que debe cumplir esta función pueden consultarse en [21]. ▼H : Término correspondiente a la energía plástica degradada, suponiendo un endurecimiento isótropo del material. De nuevo existe una función h(d) que cumple el mismo propósito y condiciones que g(d) . La definición de la energía libre plástica viene establecida por: ψp(α,∇xα):=1 2Hα2+f(α)+ 1 2l2 p|∇xα|2(1.2) Donde la contribución del endurecimiento se divide en un primer término de carácter lineal, viniendo caracterizado por el módulo de endurecimiento H , y el término no lineal, que puede o no tenerse en cuenta. El último sumando corresponde al carácter no local del modelo de plasticidad, donde el ancho de la zona plástica viene delimitado por lp (parámetro de control). En [ 19 ], se considera la plasticidad no local (también conocida como Gradient-Plasticity) añadiendo una ecuación adicional al
1.2 Revisión del estado del arte 7 sistema resultante del principio de minimización. En el caso de no considerarse estaríamos tratando con modelos de plasticidad local. ▼F: Término correspondiente al trabajo disipativo de daño. En el contexto que se trata toma la forma: ˆ ψf(d,∇xd):=Gc cωω(d) lc +lc|∇xd|2(1.3) Donde Gc es la tenacidad a fractura del material, y ω(d) y cω son los coeficientes propios del funcional de Ambrosio-Tortorelli. En el caso más común (conocido como AT-2) se considera la función de densidad de superficie de grieta definida en [21]: ˆ ψf(d,∇xd):=Gc 2d2 lc +lc|∇xd|2(1.4) Puede comprobarse que utilizando este tipo de formulación se invoca el concepto de Γ -Convergence con respecto al modelo de fractura frágil de Griffith llevando a 0 la longitud característica lc , véase [ 21 , 23 ]. Sin embargo, tiene el inconveniente de que el daño 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ño [ 24 , 25 ]. También existen otros autores que emplean un funcional tipo AT-1, que desarrolla un tramo lineal previo al daño, aunque es necesario realizar ajustes para mantener la condición de irreversibilidad [ 12 ]. Dicha condición consiste en que el material no puede recuperarse del daño, es decir, las grietas creadas o desarrolladas no pueden desvanecerse, por lo que el proceso no es reversible. ▼P : Término correspondiente al trabajo disipativo plástico. Para establecer el acoplamiento entre las variables de daño y plasticidad, distintos autores utilizan la función de degradación p(d). Con estos términos, podemos ensamblar la densidad de energía e integrarla en el dominio de estudio para después obtener la formulación energética basada en: (i) la condición de estabilidad ( ˙ d>0 y ˙ α>0 ), (ii) el balance energético ( ˙ Π=0 , siendo Π el potencial total del sistema), (iii) la condición de irreversibilidad ( δΠ≥0 ). Pueden consultarse desarrollos teóricos y aplicaciones directas de estos principios en [ 26 , 27 ]. Para modelar el comportamiento tipo frágil, únicamente se consideran los términos E y F de (1.1) , ya que no existe desarrollo de deformaciones plásticas. En cambio, para el comportamiento tipo dúctil, es necesario considerar todos los términos. Mientras que en el primer caso no existe mucha variabilidad en la formulación de los términos bajo el contexto del PFAF (salvo los coeficientes propios de la formulación de Ambrosio-Tortorelli (1.3) ), la formulación de cada término en comportamiento dúctil suele variar entre autores. Esto es debido principalmente por dos razones: ▼ No existe una formulación variacional de la fractura dúctil, a diferencia que en frágil [ 5 ], por lo que es necesario idear una formulación que cumpla los propósitos de la versión regularizada, que es donde se integra la metodología phase-field [17]. ▼ No existe un método estándar para acoplar fractura y plasticidad, por lo que cada autor propone su metodología. Así, existen numerosas metodologías con el propósito de simular el comportamiento dúctil en materiales. En la Tabla 1.2 (construida a partir de [ 22 ]), se recoge un resumen del planteamiento matemático, basado en (1.1) , en el que cada autor propone su modelo en base a las justificaciones que pueden consultarse en los artículos correspondientes recogidos en la primera columna. Los comentarios acerca de las diferencias entre autores pueden consultarse en [ 22 ]. Aquí destacaremos las más importantes: ▼ En el término de energía elástica, la función g(d) varía su expresión en función del autor que se considere y las necesidades que precise, como la modelización del comportamiento post-crítico o ajuste del comportamiento de las variables. Una aplicación reseñable de la función de degradación en este término para una formulación de carácter anisótropo se encuentra en [ 29 ]. Por otra parte, es destacable el caso de Ambati et al. [ 17 ], que añade el acoplamiento entre fractura y plasticidad en esta
8Capítulo 1. Introducción Tabla 1.2 Comparativa de modelos usando la formulación (1.1). Histórico de modelos Modelo [Referencia] E g(d)ψe(ε ε ε−ε ε εp) H h(d)ψp(α,∇xα) F ˆ ψf(d,∇xd) P p(d)ˆ ψp(α) Alessi et al. [12] 1−d)2 k−(k−1)(1−d)2ψe(ε ε ε−ε ε εp) con k>0 -3Gc 8d lc+lc|∇xd|2(1−d)2nσpα con n>0 Borden et al. [20] g(d)ψe(ε ε ε−ε ε εp) con g(d) = (1−d)2(1+d(2−k)) yk>0 h(d)1 2Hα2 con h(d)≡g(d)Gc 2d2 lc+lc|∇xd|2 p(d)σpα con p(d)≡g(d) Duda et al. [16] (1−d)2ψe(ε ε ε−ε ε εp) 1 2Hα2σpα Kuhn et al. [18] (1−d)21 2Hα2 con h(d)≡g(d)(1−d)2σpα Miehe et al. [19] (1−d)21 2Hα2+f(α)+ 1 2l2 p|∇xα|2 con f(α) = (σ∞−σp)α+1 ηe−ηα yη>0,σ∞>σp1−(1−d)2ωc+ωc ξld2 lc+lc|∇xd|2 con ωc>0yξ>0 Ambati et al. [17] (1−d)2b(α)ψe(ε ε ε−ε ε εp) donde b(α) = (α/αcrit )m, con αcrit >0ym>01 2Hα2 Gc 2d2 lc+lc|∇xd|2σpα Yin et al. [28] (1−d)2ψe(ε ε ε−ε ε εp)f(ζ(t))G0 c 2d2 lc+lc|∇xd|2
1.2 Revisión del estado del arte 9 función, con el objetivo de establecer el comienzo de daño para un cierto nivel de deformación plástica acumulada. ▼ Comparando los términos de energía plástica y disipación plástica podemos conocer aquellos autores que optan por utilizar un modelo acoplado, en base al empleo de factores con variable d . En el caso de Miehe et al. [ 19 ], se acopla además con un modelo de plasticidad de carácter no local, al emplear un término proporcional a |∇xα|2y delimitar la zona de acción con el parámetro de control lp. ▼ En cuanto al término disipativo de fractura se conserva la misma estructura en la mayoría de los casos, salvo en Miehe et al. [ 19 ] que emplea una formulación que añade dos parámetros adicionales: wc , que controla el comienzo del mecanismo de fractura y sustituye a Gc ; y ξ que controla la pendiente del tramo de reblandecimiento, aludiendo a la modelización del comportamiento post-crítico. Por otro lado, Alessi et al. [ 12 ] utiliza un término de disipación de daño lineal con un AT-1. Con base en [ 22 ], la metodología de Alessi et al. conlleva la pérdida de las fronteras de la variable de phase-field d , que deja de moverse de forma natural entre 0 y 1 y necesita de algoritmos específicos para su recuperación. Por último, es interesante comentar que Yin et al. [ 28 ] opta por establecer el acoplamiento entre fractura y plasticidad en este término, a partir de la definición de una tenacidad de fractura degradable con la variable de endurecimiento. Fractura frágil Una de las motivaciones principales que persigue la modelización de procesos como el de fractura es reproducir los estudios de carácter experimental. Con este propósito, se intentan suplir limitaciones del método experimental como estados de carga complejos o el gran coste económico que pueden llevar consigo algunos ensayos. Para entender de forma adecuada la modelización de procesos de fractura frágil es preciso entender el desarrollo y consecución de las técnicas numéricas que han llevado a alcanzar el estado actual. A continuación se ofrece una pequeña síntesis de esta idea. La base de la metodología computacional en el campo que estamos tratando residió de forma primigenia en la Mecánica de la Fractura Elástica-Lineal (MFEL), ya que es una metodología bastante simple y robusta de implementar. Sin embargo, a medida que el problema se va complicando mediante la implementación de geometrías y estados de carga más complejos, los mecanismos de fractura se hacen también más complejos, pudiéndose dar varios a la vez de distintas morfologías. Este hecho es el principal problema de la metodología MFEL, ya que se le hace bastante complicado reproducir caminos de grieta fuera de los casos básicos. Normalmente, en suelos y rocas topamos con el caso más complejo, ya que recordemos que estamos ante un material ya de por sí bastante difícil de modelar, por lo que la MFEL no parece la mejor manera de tratar el problema de estudio. De esta forma, con el objetivo de solventar estas dificultades, se abren paso las metodologías asociadas al Método de los Elementos Finitos (MEF) y el Método de los Elementos de Contorno (MEC). Para nuestro caso particular nos focalizaremos en el MEF, que es la que más extendida está en el ámbito científico de modelización de procesos de fractura. Las metodologías principales que se han desarrollado en los últimos años se dividen principalmente en: ▼ Técnicas basadas en la metodología de grieta difusa: Con base en el campo de la Mecánica de Daño Continuo (MCD), se basa en la definición de variables auxiliares de daño para la modelización de los distintos mecanismos de fractura que se pueden dar en el sólido. El principal problema que presenta es que su versión local tiene una fuerte dependencia entre la solución y la malla de discretización empleada. Para ello, distintos autores han presentado alternativas como el uso de procedimientos de regularización de la formulación o metodologías no locales, como se verá más adelante [ 30 , 31 , 32 , 33 ]. ▼ Técnicas de implementación de discontinuidades fuertes: A modo de de reproducir la grieta de forma explícita, esta metodología se divide en estrategias de enriquecimiento del campo cinemático de interpolación a nivel nodal para evitar el remallado del entorno de la grieta, teniendo el eXtended FEM (X-FEM) como ejemplo; o técnicas basadas en el empleo de elementos modelados con formulaciones basadas en modos incompatibles de deformación, como el enhanced FEM (E-FEM) [ 34 , 35 , 36 , 37 , 38 ]. ▼ Técnicas de inclusión de elementos de interfaz cohesiva de carácter adaptativo o permanente en los bordes de la malla de elementos finitos [39, 40, 41, 42, 43, 44]. Por tanto, podemos ver que mientras el primer grupo realiza un enfoque de la grieta de carácter difuso, los otros grupos la aproximan de manera discreta. El enfoque de carácter discreto suele tener el inconveniente
10 Capítulo 1. Introducción de tener que buscar alternativas más complicadas para evitar la dependencia de la solución con la malla de elementos finitos. Todas estas técnicas presentan la limitación común de ser incapaces de predecir la nucleación de grieta y dar problemas en modelizaciones de procesos de ramificación y coalescencia para frentes múltiples. En este momento es cuando surge el origen de la metodología phase-field, que ofrece una nueva perspectiva frente a modelizaciones de grietas vistas hasta ahora de carácter más complejo, con carácter mixto al estar fundamentada en técnicas de grieta difusa y modelos de daño no local asociados al uso del gradiente espacial de la variable de daño. Es interesante comparar las similitudes y disparidades de los modelos asociados al gradiente con la metodología phase-phield, véase [ 45 ]. Para controlar la extensión del gradiente de daño en phase-field se establece una ecuación adicional en función de una longitud característica lc , que actúa como parámetro de control. Un ejemplo de aplicación ampliamente estudiado y tratado es la modelización de fractura frágil, objeto de este subapartado. Otra aplicación interesante aplicada específicamente en geomateriales, distinguiendo distinto comportamiento a tracción y a compresión puede consultarse en [46]. De forma particular, la base computacional de la mecánica de la fractura frágil se basa en los desarrollos de la Teoría Clásica de Griffith [ 4 ]. Dicha teoría propone un criterio por el cual la grieta se propaga si el ratio de energía disponible para el crecimiento de grieta supera cierto valor crítico. Ambos valores se consideran constantes para condiciones dadas y se denominan G y Gc , respectivamente. Al presentar los mismos problemas comentados en el párrafo anterior, se opta por utilizar métodos variacionales basados en procesos de minimización de energía para resolver el problema [ 5 ], como ya ha sido introducido anteriormente. Dicho método se basa en incluir un término asociado a Gc de carácter superficial en la expresión del potencial total de energía interna, dentro del término disipativo. Finalmente, para obtener una formulación manejable, se obtiene una versión regularizada a través de la metodología phase-field. Con esta técnica, el crecimiento y desarrollo de grietas se modela a partir de la evolución y crecimiento de una variable auxiliar de daño en una malla fija de elementos finitos. Los puntos fuertes del uso de la metodología phase-field son: (i) la Γ -Convergencia al Modelo canónico de Griffith [ 21 , 23 ] (aplicado posteriormente en (1.11) ), (ii) evitar el modelado de discontinuidades y (iii) posibilidad de un tipo de implementación de carácter robusto, flexible y efectivo mediante un solver básico. La formulación analizada en este trabajo se basa en la estructura termodinámicamente consistente de la metodología phase-field para fractura frágil que se recoge en [ 13 , 21 , 24 ]. Desde este punto, se hará uso de aquellos desarrollos que se entiendan necesarios para comprender la materia. En el caso de necesidad de profundidad de contenido pueden consultarse las referencias anteriores. La configuración de referencia de nuestro problema de estudio, con referencia en la Figura 1.4, está dado por el dominio B⊂Rn cuyo contorno viene delimitado por ∂B⊂Rn−1 . Rn hace referencia al espacio Euclídeo de orden n, siendo en este caso de aplicación el caso dimensional (n=2) y tridimensional (n=3). Los campos básicos del problema de fractura frágil son el de desplazamientos u y el de daño, que, como ya vimos en el apartado anterior, viene definido por la variable de phase field d en una formulación de tipología no local. Las condiciones de contorno vienen definidas en el contorno ∂Bu para las acciones cinemáticas y ∂Bt para las acciones estáticas. Por otro lado, como se verá más adelante, las condiciones de contorno de la variable de daño se definen en el contorno completo ∂B. En el caso del campo de desplazamientos utambién se hará uso tensor de pequeñas deformaciones: ε ε ε(u):=∇xsu=1 2∇xu+∇T xu(1.5) Uno de los aspectos más importantes de la formulación empleada [ 21 ] es la descomposición espectral del tensor de (1.5) de deformaciones que permite separar la componente positiva y negativa, que será útil para desarrollos posteriores: ε ε ε=ε ε ε++ε ε ε− ε ε ε±:= n ∑ i=1⟨εi⟩±ni⊗ni (1.6) Donde la descomposición se hace en función de los autovalores εi y los autovectores ni de ε ε ε . El operador ⟨z⟩±son los paréntesis de Macaulay que representan la función a trozos ⟨z⟩±:= (z±|z|)/2. Para obtener la formulación energética del problema es necesario establecer el funcional del sistema completo, que tiene la siguiente forma:
1.2 Revisión del estado del arte 11 Π(u,d) = Π(u,d)int −Π(u)rmext =E(u,d)+F(d,∇xd)−W(u) = =ZB ψe(ε ε ε,d)dV +ZB Gcγ(d,∇xd)dV | {z } Πint −ZB b·udV +Z∂B ¯ t·udA | {z } Πext (1.7) Donde se realiza la siguiente descomposición (los términos son análogos a los de (1.1)): ▼ Energía libre elástica (E): Caracterizada por una descomposición espectral en cuanto a su parte positiva y negativa. Esta idea se basa en el desarrollo planteado en (1.6): ψe(ε ε ε,d) = g(d)ψe+(ε ε ε)+ψe−(ε ε ε) ψe±(ε ε ε) = λ 2⟨tr[ε ε ε]⟩±2+µtr ε ε ε2 ±(1.8) Donde la función de degradación g(d) únicamente afecta a la parte positiva de la energía elástica. Esto es coherente con el hecho de que la parte positiva del tensor de deformaciones, asociada a estados tensionales de tracción, es la que cede energía para desarrollar mecanismos de fractura, a modo de disipación: g(d) = (1−d)2+K(1.9) El término cuadrático va acompañado de un parámetro residual K para evitar errores de condicionamiento numérico. ▼ Energía disipada por procesos de fractura (F): Previamente definida en (1.4) , está caracterizada por la tenacidad a fractura Gcy la función de de densidad superficial de grieta γ: γ(d,∇xd) = 1 2lc d2+lc 2|∇xd|2(1.10) Como ya se ha comentado anteriormente, a través de esta función puede obtenerse la longitud de grietra discreta Γcuando lc→0a través de la Γ-Convergencia: Γ=l´ ım lc→0ZB γ(d,∇xd)dV (1.11) ▼ Energía referida al trabajo de las fuerzas externas (W): Caracterizada por los dos últimos sumandos de (1.7), están referidos al campo de desplazamientos u. Estableciendo el PTV para definir el equilibrio entre el incremento de energía total y el trabajo realizado por las fuerzas externas a partir de (1.7) nos queda el sistema (pueden consultarse los pasos intermedios en [ 13 ]): Equilibrio: ∇x·σ+b=0 Sistema KKT: ˙ d≥0 f−Gcδdγ≤0 ˙ d(f−Gcδdγ) = 0 CDC: u=¯ uen ∂Bu σ·n=¯ ten ∂Bt ∇xd·n=0en ∂B (1.12) Que distingue entre la forma fuerte de la ecuación de equilibrio para la determinación del campo u , un sistema tipo Karush-Khun-Tucker para el campo d y las condiciones de contorno. En el caso del campo de daño, se impone la condición de contorno sobre la dirección del gradiente de daño, que debe ser perpendicular a la normal en el contorno. Por otro lado, dentro del sistema (1.12) aparece:
12 Capítulo 1. Introducción ▼Tensor de tensiones: σ σ σ:=∂ε ε εψ(ε ε ε,d) = g(d)[λ⟨tr[ε ε ε]⟩+I+2µε ε ε+]+[λ⟨tr[ε ε ε]⟩−I+2µε ε ε−] ▼Energetic force:f:=−∂dψ(ε ε ε,d) = 2(1−d)ψe+(ε ε ε) ▼Derivada funcional de la densidad superficial de grieta: δdγ=∂dγ−∇x·∂∇xdγ=1 lcd−l2 c∆d Utilizando la formulación basada una estructura del tipo Ginzburg-Landau generalizada (véase [ 24 , 25 ]) el sistema KKT puede cambiarse por: η˙ d Evolución =2(1−d)H Driving Force −Gc lc [d−l2 c∆d] R. Geométrica (1.13) Podemos observar que el mecanismo de fractura frágil por el cual se desarrollan las grietas se basa en un equilibrio, a modo de competición, entre el término de la resistencia geométrica y el término de Driving Force, gobernada por la variable histórica H. Dicha variable toma la forma: H(x,t):=m´ ax s∈[0,t]ψe+(ε ε ε(x,s)) (1.14) Por tanto, este tipo de comportamiento viene controlado por la la parte positiva de la energía libre elástica, tal y como se anticipó al principio del capítulo. A (1.14) pueden añadirse términos a modo de thresholds para delimitar el comienzo de daño, y evitar la desventaja del empleo de un AT-2. Sin embargo, el hacer uso de estos términos ocasiona la pérdida de consistencia variacional, por lo que no se garantiza que el modelo converja al canónico de Griffith. Integrando el sistema (1.12) en el dominio de aplicación es posible utilizar un algoritmo de resolución basado en desacoplar el problema en los campos u y d haciendo uso de (1.13) . A continuación se muestra un posible modo de resolver el problema basado en [13]: Rate independent Operator Split Algorithm Esquema escalonado de resolución del problema para el pseudo-tiempo [tn,tn+1]: I. Inicialización: Las variables un , dn y Hn son conocidas. Se prescriben las cargas y desplazamientos impuestos bn+1,¯ un+1y¯ tn+1. II. Cálculo de H: H= ψ+ e+(ε ε ε)si ψ+ e+(ε ε ε)>Hn Hnen otro caso III. Cálculo del campo daño (PF): d=arg´ ınf dZBGcγ(d,∇xd)−(1−d)2HdV IV. Cálculo del campo de desplazamientos: Con dcongelado u=arg´ ınf uZB [ψ(ε ε ε,d)−b·u]dV −Z∂Bt ¯ t·udA V. Repetición de los pasos II-IV hasta que se alcance la convergencia para la pareja de variables (un+1,dn+1)
1.2 Revisión del estado del arte 13 Donde rate independent hace alusión a no considerar términos viscosos en la ecuación (1.13) , por lo que el término asociado al parámetro viscoso η desaparece. La viscosidad suele ser empleada para estabilizar la solución post-crítica. La formulación en elementos finitos no resulta de interés para este problema, ya que estamos interesados en modelar el caso acoplado con la plasticidad. De cualquier forma, puede consultarse desarrollos análogos por la cual se obtiene el sistema de ecuaciones a resolver por el MEF y diversas aplicaciones de interés en [ 14 , 47 ]. Fractura dúctil La modelización de la fractura dúctil emerge de la base del conocimiento desarrollada para fractura frágil. Sin embargo, debemos partir de la idea de que el mecanismo de funcionamiento no es el mismo. Según [ 16 ], en la MFEL clásica los procesos de fractura pueden darse únicamente en las puntas de grieta, donde el crecimiento viene dominado por la integral-J en el entorno. En el caso elastoplástico, dicha integral se desvanece y en principio, no existiría dicha fuente de crecimiento. Esto no deja de ser un mal uso de la teoría, ya que en plasticidad la zona de fractura se forma delante de la punta de la grieta. Por tanto, los modelos que se desarrollen deben reflejar la existencia de dicha zona de fractura. De esta forma, la teoría desarrollada para fractura dúctil parte de utilizar la plasticidad para complementar la energía libre dentro de la integral (1.7) en el término Hde (1.1) y los desarrollos de phase-field del apartado anterior, con el objetivo de reproducir el proceso de formación de la zona de fractura. Para materiales dúctiles es necesario considerar un término disipativo adicional debido al desarrollo de deformaciones plásticas. Dicho término no puede ser aproximado desde un enfoque variacional (no existe una teoría variacional que describa el proceso de fractura de tipo dúctil), luego en la mayoría de los casos no deja de ser un artificio matemático sin sentido físico, aunque es importante dárselo para la justificación del modelo. Uno de los aspectos más importantes de este tipo de modelizaciones es la manera en la que se considera el acoplamiento entre la plasticidad y el daño, existiendo ejemplos de ello según los autores de la Tabla 1.2. Siguiendo los mismos pasos que en fractura frágil, utilizando el PTV podemos obtener la formulación compacta del problema (extraída de [ 48 ]), que puede diferenciarse en tres subproblemas según la variable principal (campo de desplazamientos u, campo de deformaciones plásticas ε ε εpy la variable phase-field d): ▼Equilibrio ∇x·σ+b=0 CDC:(σ·n=¯ ten ∂Bt u=¯ uen ∂Bu (1.15) ▼Plasticidad Sistema KKT: ˙ α≥0 ¯ fp(u,α)≤0 ˙ α¯ fp(u,α) = 0 Regla de flujo: ˙ εp=˙ αN Función de plastificación: ¯ fp(u,α):=¯ σeq(u,α)−Hα−σy (1.16) Donde la función de plastificación compara el término de ¯ σeq(u,α) con el límite elástico σy(α) . En este caso, ¯ σeq(u,α) varía según el criterio de plastificación que se utilice y σy(α) se toma lineal, con un módulo de endurecimiento H.
20 Capítulo 1. Introducción Figura 1.8 Diagrama de flujo para la obtención del vector de tensiones para un paso de tiempo dado usando el algoritmo de Return Mapping (extraído de [6]).
1.3 Organización del documento 21 Por otro lado, para derivar la matriz Cepc es preferible obtener primeramente su expresión en el espacio de tensiones principales, por simplicidad, a partir de las siguientes expresiones [6]: Cepc ppal(i,j) = ∂σi ∂εe trial j ,i,j=1,2,3 σi=ˆ σ(εe trial 1,εe trial 2,εe trial 3) (1.31) Donde básicamente se deriva la función implícita ˆ σ (obtenida a partir de (1.30) ) con respecto a sus variables, que son las deformaciones elásticas de prueba (elastic trial strains). Estas deformaciones se obtienen sumando el incremento de deformación del paso actual a las deformaciones elásticas obtenidas en el paso anterior. De forma análoga al caso anterior, en la subrutina CTMC de [ 6 ] se realiza una implementación en base a las expresiones desarrolladas para el cálculo de tensiones en el modelo de Mohr-Coulomb, para cada una de las fases de la Figura 1.8. 1.3 Organización del documento Tras el desarrollo de la base teórica necesaria para la ejecución del objetivo propuesto al principio del capítulo, el trabajo se estructura como sigue: ▼ En el Capítulo 2, se abordará la implementación de la ley constitutiva haciendo uso del modelo de plasticidad de Mohr-Coulomb a través de los algoritmos y códigos propuestos en [ 6 ]. Posteriormente, se aplicará y validará la metodología propuesta a través de un par de ejercicios. ▼ En el Capítulo 3, se propone un modelo acoplado de fractura y plasticidad, en base a una serie de metodologías recogidas en la Tabla 1.2. Mediante la implementación de dichas metodologías se persigue el estudio de las fuerzas precursoras de fractura y cómo sus parámetros influyen en la forma y comportamiento de la solución obtenida. ▼ En el Capítulo 4, se expondrán las conclusiones y líneas de investigación futuras en base al modelo desarrollado, con el objetivo de extender su aplicación y validez.
2 Implementación de la ley constitutiva Co n el fin de desarrollar el modelo de fractura considerando la no linealidad material en base a la teoría de la plasticidad, en este capítulo se abordará la implementación de la ley de comportamiento referida a los fundamentos establecidos en el apartado 1.2.2 de este documento. Al tratar con un material con posibilidad de desarrollar deformaciones plásticas, los ingredientes básicos que definen dicho comportamiento son: (i) el criterio de plastificación, (ii) la regla de flujo y (iii) la ley de endurecimiento (en el caso de que permitamos el endurecimiento de ciertas variables). La secuencia por la que pasará la implementación realizada cuenta con los siguientes bloques: ▼ Inicialización del modelo. Se declaran las subrutinas internas, variables y parámetros que se emplean en el programa. ▼Actualización del estado de tensiones. A través de un estado de tensiones de predicción y un proceso de corrección de carácter iterativo, se actualiza el estado de tensiones para cada paso, controlado por ciertos valores de tolerancia. ▼Modelización de la matriz de comportamiento elastoplástica consistente Cepc. Para obtener la matriz de comportamiento es necesario obtener la variación del estado tensional con respecto al campo de deformaciones. Dicho proceso se llevará a cabo en primer lugar en el espacio de tensiones principales, donde las expresiones a derivar son más tratables, para luego pasar al espacio global. La ley constitutiva del material se modelizará a través del Modelo de Mohr-Coulomb de carácter no asociativo, con una ley de endurecimiento aplicada a la variable de cohesión. En el espacio de tensiones principales, con σ1≥σ2≥σ3 (primer sextante), la función de plastificación toma la forma: Φ1(σ σ σ,c) = (σ1−σ3)+(σ1+σ3)sinφ−2ccosφ(2.1) El potencial plástico, basado en (2.1) para el primer sextante es: Ψp1(σ σ σ,c) = (σ1−σ3)+(σ1+σ3)sinψ(2.2) Esta formulación del potencial plástico resulta ser el caso más sencillo al únicamente sustituir el ángulo de rozamiento φ por el de dilatancia ψ . Existen otras formulaciones para el potencial plástico más complejas, como la propuesta por Menétrev y William en [ 52 ], consistente en una una hipérbola en el plano p-q y una elipse en el plano π. Por último, la ley de endurecimiento para la cohesión estará en función de la deformación plástica acumulada c(¯ εp) , pudiendo adoptar diversas formas según conveniencia (lineal, exponencial, racional). Este tipo de leyes se utilizan para modelar los comportamientos de endurecimiento (la cohesión aumenta de valor ha medida que se van acumulando deformaciones plásticas) y reblandecimiento (la cohesión disminuye ha medida que aumenta la deformación plástica acumulada). Algunos ejemplos de aplicación se muestran en la Figura 2.1. 23
24 Capítulo 2. Implementación de la ley constitutiva Ejemplos de leyes de reblandecimiento Ejemplos de leyes de endurecimiento Figura 2.1 Ejemplos de leyes de endurecimiento para la cohesión c(¯ εp). 2.1 Implementación y desarrollo del código La implementación de esta ley se llevará a cabo a través de una UMAT que se procesará en ABAQUS . La UMAT no es más que una subrutina que se encarga de establecer la ley constitutiva de un material de forma genérica. Tiene como fin obtener el tensor de tensiones y la matriz de rigidez para ensamblar el sistema de ecuaciones a resolver por el MEF. El código de la UMAT , que se utilizará como guía para explicar los conceptos desarrollados,se desarrollará en lenguaje FORTRAN. Los procedimientos y cálculos realizados toman la base y la metodología empleada en [ 6 ], de los cuales se ofrece en este trabajo un resumen conciso con las ecuaciones más significativas y útiles para nuestro objetivo. 2.1.1 Inicialización y declaración de variables y subrutinas En primer lugar, se establece la subrutina principal a partir de la cual deriva todo el código diseñado. En dicha subrutina se declaran todas las variables propias de ABAQUS y de nuestro código, que entrarán en subrutinas internas posteriores para la realización de los cálculos (Código 2.1). Código 2.1 Inicialización subrutina UMAT. Declaración de variables y parámetros del modelo. subroutine \texttt{UMAT}( & STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,DRPLDE,& DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,& CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,&
2.1 Implementación y desarrollo del código 25 PNEWDT,CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC) !------------------------------------------------------------------------- ! Plastic strain update for a elastoplastic with linear hardening ! Mohr-Coulomb soil. Applications for plane strain problems. ! Non-associated linear Mohr-Coulomb material model with input ! 1) Young’s modulus, E ! 2) Poisson’s ratio, NU ! 3) Hardening law: initial cohesion, C0 and hardening modulus, H ! 4) friction angle, PHI [degrees] ! 5) dilation angle, PSI [degrees] !----------------------------------------------------------------------- implicit none !----------------------------------------------------------------------- ! Variables defined by ABAQUS input/output !----------------------------------------------------------------------- character(len=80), intent(in) :: CMNAME integer(4), intent(in) :: NTENS integer(4), intent(in) :: NSTATV integer(4), intent(in) :: NPROPS integer(4), intent(in) :: NDI integer(4), intent(in) :: NSHR integer(4), intent(in) :: NOEL integer(4), intent(in) :: NPT integer(4), intent(in) :: LAYER integer(4), intent(in) :: KSPT integer(4), intent(in) :: KSTEP integer(4), intent(in) :: KINC real(8), intent(in) :: CELENT real(8), intent(in) :: DTIME real(8), intent(in) :: TEMP real(8), intent(in) :: DTEMP real(8), intent(in) :: STRAN(NTENS) real(8), intent(in) :: DSTRAN(NTENS) real(8), intent(in) :: PROPS(NPROPS) real(8), intent(in) :: PREDEF(1) real(8), intent(in) :: DPRED(1) real(8), intent(in) :: TIME(2) real(8), intent(in) :: COORDS(3) real(8), intent(in) :: DROT(3,3) real(8), intent(in) :: DFGRD0(3,3) real(8), intent(in) :: DFGRD1(3,3) real(8), intent(inout) :: SSE real(8), intent(inout) :: SPD real(8), intent(inout) :: SCD real(8), intent(inout) :: RPL real(8), intent(inout) :: DRPLDT real(8), intent(inout) :: STATEV(NSTATV) real(8), intent(inout) :: STRESS(NTENS) real(8), intent(inout) :: DDSDDE(NTENS,NTENS) real(8), intent(inout) :: DDSDDT(NTENS) real(8), intent(inout) :: DRPLDE(NTENS) real(8), intent(inout) :: PNEWDT !----------------------------------------------------------------------- ! Variables defined by this subroutine
26 Capítulo 2. Implementación de la ley constitutiva !----------------------------------------------------------------------- real(8) :: DGAM(2), STRES(NTENS) real(8) :: DMATX(NTENS,NTENS) real(8) :: PHI_D, PSI_D, PHI, PSI, COH0, H, E, NU logical :: LALGVA(5) real(8), parameter :: pi = 3.1415926535897932_8 Donde las variables más significativas y utilizadas en el resto del código son: ▼ABAQUS: ■Argumentos de entrada: ♦NTENS: Número de componentes del vector tensión/deformación (en notación de Voigt). ♦STRAN(NTENS): Deformación total al principio del incremento. ♦DSTRAN(NTENS): Incremento de deformación. ♦ NSTATV: Número de variables de estado. Este tipo de variables se utilizan en ABAQUS para guardar variables propias de nuestro problema que son dependientes del paso de tiempo. ♦ PROPS(NPROPS): Array con las propiedades del material a modelar, limitado por su número a partir de NPROPS. Para nuestro caso, y tal y como se detalla en el Código 2.1, para el Modelo Mohr-Coulomb es necesario definir las propiedades E y ν para la fase elástica, los ángulos φ y ψ para la función de plastificación y potencial plástico respectivamente, y parámetros para establecer la ley de endurecimiento como c0 (cohesión inicial) y H (módulo de endurecimiento) en el caso de ser de carácter lineal. ♦KSTEP: Número de paso de tiempo. ♦KINK: Número de incremento dentro de un paso de tiempo. ■Argumentos de salida: ♦DDSDDE(NTENS,NTENS): Matriz elastoplástica constitutiva. ♦STRESS(NTENS): Vector de tensiones. ♦STATEV(NSTATEV): Variables de estado. ▼Propias: ■ DGAM: Array para guardar los incrementos de los multiplicadores plásticos, incógnita dentro del algoritmo de Return Mapping. ■ STRES(NTENS), DMATX(NTENS,NTENS): Arrays para guardar e ir actualizando los argumentos STRESS y DDSDDE respectivamente durante el número de iteraciones que sea necesario para alcanzar la convergencia deseada dentro de un paso de tiempo. ■E, NU: Propiedades elásticas del material. ■COH0, H: Propiedades de endurecimiento del material. ■ PHI, PSI: Ángulos de rozamiento y dilatancia (propiedades propias de la formulación del Modelo de Mohr-Coulomb). A continuación, debemos declarar las subrutinas internas que entran dentro de la subrutina general, además de las propiedades que se mantendrán constantes durante todo el proceso de resolución (Código 2.2).
2.1 Implementación y desarrollo del código 27 Código 2.2 Inicialización subrutina UMAT . Establecimiento de subrutinas internas y definición de los parámetros del material. ! !---- Explicit interfaces ---------------------------------- interface subroutine MCSTUP(DGAM,LALGVA,RSTAVA, & STRAT,STRES,MSTRE, & KSTEP,KINC,PROPS,PHI,PSI) logical, intent(out) :: LALGVA(5) real(8), intent(out) :: DGAM(2) integer(4), intent(in) :: MSTRE real(8), intent(inout) :: RSTAVA(MSTRE+1) real(8), intent(in) :: STRAT(MSTRE) real(8), intent(out) :: STRES(MSTRE) integer(4), intent(in) :: KSTEP,KINC real(8), intent(in) :: PROPS(4), PHI, PSI end subroutine MCSTUP subroutine MCCONTANMA(DMATX,LALGVA,RSTAVA,& STRAT,STRES,MSTRE,PROPS,PHI,PSI) logical, intent(in) :: LALGVA(5) integer(4), intent(in) :: MSTRE real(8), intent(out) :: DMATX(MSTRE,MSTRE) real(8), intent(in) :: RSTAVA(MSTRE+1) real(8), intent(in) :: STRAT(MSTRE) real(8), intent(in) :: STRES(MSTRE) real(8), intent(in) :: PROPS(4),PHI,PSI end subroutine MCCONTANMA end interface ! !----------------------------------------------------------- !----------------------------------------------------------------------- ! Material Properties !----------------------------------------------------------------------- E = PROPS(1) ! E - Young’s modulus NU = PROPS(2) ! NU - Poisson’s ratio COH0 = PROPS(3) ! c - initial cohesion H = PROPS(4) ! H - hardening modulus PHI_D = PROPS(5) ! PHI - angle of friction (degrees) PSI_D = PROPS(6) ! PSI - angle of dilation (degrees) !------------------------------------------------------- ! Convertion of angles from degrees to radians !------------------------------------------------------- PHI = pi*PHI_D/180.0_8 ! [rad] Friction angle PSI = pi*PSI_D/180.0_8 ! [rad] Dilatation angle Por una parte, la subrutina MCSTUP se encargará de actualizar el estado de tensiones en los puntos de Gauss de los elementos del modelo conforme al algoritmo de Return Mapping, mientras que la subrutina MCCONTANMA realizará la construcción de la matriz de rigidez de forma consistente con la formulación de las tensiones. Por último, en la subrutina principal se llaman a las subrutinas anteriores para la realización de los cálculos. Como vector de deformaciones de entrada se utiliza las deformaciones elásticas del paso anterior más el incremento de deformación del paso actual, de acuerdo con las premisas del algoritmo. Una vez alcanzada la
28 Capítulo 2. Implementación de la ley constitutiva convergencias de cada subrutina se actualizan los argumentos principales de la UMAT : STRESS Y DDSDDE (Código 2.3). Código 2.3 Desarrollo subrutina UMAT . Llamadas a subrutinas internas y actualización de los argumentos de salida principales. !---------------------------------------------------------------------- - ! Stress update and creation of consistent constitutive matrix !----------------------------------------------------------------------- call MCSTUP(& & DGAM,LALGVA,STATEV(1:5),& & STATEV(1:4)+DSTRAN(:),STRES,NTENS,KSTEP,KINC,PROPS(1:4),PHI,PSI) call MCCONTANMA(& & DMATX,LALGVA,STATEV(1:5),STATEV(1:4)+DSTRAN(:),STRES,NTENS,PROPS(1:4), PHI,PSI) !----------------------------------------------------------------------- ! Postconditioner for ABAQUS !----------------------------------------------------------------------- STRESS = STRES ! updated stress at the Gauss point DDSDDE = DMATX ! consistent constitutive matrix !----------------------------------------------------------------------- end subroutine \texttt{UMAT} !----------------------------------------------------------------------- 2.1.2 Actualización del estado tensional Para la actualización del estado tensional seguiremos el diagrama de flujo de la Figura 1.8 conforme al esquema del Return Mapping ilustrado en la Figura 2.2. Figura 2.2 Esquema de aplicación del Return Mapping (extraído de [6]).
2.1 Implementación y desarrollo del código 29 Cuya formulación de tipo de incremental para las tensiones, de la cual deriva la actualización del estado tensional mostrada en (1.30), es [6]: σ σ σn+1=σ σ σtrial n+1−Ce:∆ε ε εp(2.3) Donde Ce hace refencia a la matriz constitutiva elástica. La aplicación al Modelo de Mohr-Coulomb mediante la regla de flujo no asociativa es [6]: σ σ σn+1=σ σ σtrial n+1−Ce: 6 ∑ i=1 ∆γiN N Ni n+1(2.4) Dependiente de un conjunto de incógnitas ∆γi , que son los incrementos de los multiplicadores plásticos para el pseudo-tiempo [tn,tn+1] . La subrutina planteada ha sido adaptada del código programado para HYPLAS de la rutina SUMC, de [ 6 ]. En primer lugar, como en el caso anterior, se declaran las variables, parámetros a utilizar y subrutinas internas (Código 2.4). Código 2.4 Inicialización subrutina MCSTUP. Declaración de parámetros y variables. subroutine MCSTUP(& DGAM,LALGVA,RSTAVA,STRAT,STRES,MSTRE, & KSTEP,KINC,PROPS,PHI,PSI) ! ------------------------------------------------------------------------ ! State update procedure for the Mohr-Coulomb type material model. ! ! This routine uses the fully implicit elastic predictor/return ! mapping algorithm as the state update procedure for the ! Mohr-Coulomb elasto-plastic material model. The essential return mapping is ! carried out here in principal stress space. ! ------------------------------------------------------------------------ implicit none ! --- Parameters for subroutines ------------- integer(4), parameter :: MCOMP=4, NDIM=2 integer(4), parameter :: MXITER=50 ! -------------------------------------------- logical, intent(out) :: LALGVA(5) real(8), intent(out) :: DGAM(2) integer(4), intent(in) :: MSTRE real(8), intent(inout) :: RSTAVA(MSTRE+1) real(8), intent(in) :: STRAT(MSTRE) real(8), intent(out) :: STRES(MSTRE) integer(4), intent(in) :: KINC, KSTEP(4) real(8), intent(in) :: PROPS(6), PHI, PSI ! -------------------------------------------- logical :: IFPLAS, SUFAIL, EDGE, APEX logical :: RIGHT, DUMMY real(8) :: E, NU, COH0, H real(8) :: EIGPRJ(MSTRE,2), PSTRS(3), STREST(3) real(8) :: DGAMA, DGAMB, EPBAR, EPBARN real(8) :: GMODU, BULK, R2G, R4G, R2BULK, R2CPHI, R1D3 real(8) :: EETV, PT, EETVD3 integer(4) :: II, JJ, MM, I, J real(8) :: PSTRS1, PSTRS2, PSTRS3 real(8) :: COHE, SMCT, PHIA, RES real(8) :: SCAPRD real(8) :: SPHSPS, CONSTA, R4C2PH real(8) :: DENOM, DDGAMA real(8) :: S1, S2, S3, DELTA
36 Capítulo 2. Implementación de la ley constitutiva La fase del vértice se considera la última antes de salir del algoritmo, según la Figura 1.8. El vértice de la pirámide de Mohr-Coulomb se encuentra en la diagonal principal del espacio de tensiones principales. De la función de plastificación, se puede deducir que además de que el estado tensional se encuentre en esta línea se debe cumplir: p=ccotφ(2.12) Donde pes la parte hidróstatica del tensor de tensiones, que en el espacio principal equivale a: p=1 3(σ1+σ2+σ3)(2.13) Con todo lo visto hasta el momento, el tensor de tensiones principales solución debe ser de la forma [ p,p,p ], al tener solo parte hidrostática. Para esta fase en particular, la incógnita del algoritmo del Return Mapping será el incremento de deformación plástica volumétrica ∆εp v , igual a la traza del tensor de deformaciones plásticas en el sistema principal, ya que la solución se encuentra sobre la diagonal principal. La ecuación residual a resolver es la siguiente (su deducción, basada en la ley de comportamiento elástica y (2.12) , y consideraciones pertinentes pueden consultarse en [6]): c(¯ εp n+∆¯ εp)cotφ−ptrial n+1+K∆εp v=0(2.14) Donde: ptrial n+1=K(εe 1n+εe 2n+εe 3n) ∆¯ εp=λ∆εp vλ=cosφ sinψ (2.15) El estado actualizado de tensiones en cada iteración se corresponde con: σik+1=ptrial n+1−K∆εp vi=1,2,3(2.16) Donde los valores actualizados en cada paso se corresponden con el subíndice n+1 y en cada iteración con k+1. Una vez alcanzada la convergencia del paso los valores con k+1pasan al estado n+1. De nuevo, para la ecuación (2.14) se utiliza el método de Newton-Raphson, obteniendo como resultado ∆εp v para su aplicación en (2.16) . En el caso de que el método falle en este punto, al no existir más fases, se considera que el incremento de deformación utilizado al comienzo del paso es demasiado elevado, tal y como se comentó previamente. En la siguiente parte del código se implementa la fase del vértice con las ecuaciones anteriores (Código 2.9). Código 2.9 Desarrollo subrutina MCSTUP. Fase de corrección. Aplicación del multi-vector Return Mapping al vértice. 50 CONTINUE ! Apply multi-vector return mapping to APEX ! --------------------------------------- ! Set initial guess for volumetric plastic strain increment DEPV DEPV=0.0_8 EPBAR=EPBARN COHE=COH0+EPBAR*H COTPHI=COS(PHI)/SIN(PHI) RES=COTPHI*COHE-PT ! Newton-Raphson iterations for DEPV DO NRITER=1,MXITER DENOM=COS(PHI)*COTPHI/SIN(PSI)*H+& BULK DDEPV=-RES/DENOM DEPV=DEPV+DDEPV EPBAR=EPBARN+COS(PHI)/SIN(PSI)*DEPV
2.1 Implementación y desarrollo del código 37 COHE=COH0+EPBAR*H P=PT-BULK*DEPV RES=COTPHI*COHE-P ! check for convergence RESNOR=ABS(RES) IF(PT.NE.0.0_8)RESNOR=RESNOR/ABS(PT) IF(RESNOR.LE.1.0e-10_8)THEN APEX=.TRUE. DGAMA=DEPV DGAMB=0.0_8 ! update principal stresses S1=P S2=P S3=P GOTO 70 ENDIF ENDDO Si se ha alcanzado la convergencia y el estado final de tensiones actualizado se ubica en el primer sextante en alguna de las fases anteriormente explicadas significa que el algoritmo de Return Mapping se ha empleado con éxito. El último paso en este apartado sería transformar el estado principal al global mediante los autovalores y autovectores obtenidos mediante la subrutina SPEC2 y actualizar todas las variables resultado (tensiones, deformaciones elásticas, plásticas y plástica acumulada) al paso n+1 . En el caso de las deformaciones plásticas, son fácilmente obtenibles a partir del incremento de deformaciones elásticas y el incremento de deformación total del paso: ε ε εp n+1=ε ε εp n+∆ε ε ε−∆ε ε εe(2.17) 2.1.3 Modelización de la matriz Cepc Tras obtener el estado actualizado de tensiones estamos en disposición de ensamblar la matriz de comportamiento elastoplástica consistente. Elastoplástica ya que tiene en cuenta el carácter plástico de la formulación empleada y consistente ya que está de acuerdo con el algoritmo de Return Mapping desarrollado en el apartado anterior. La definición de la matriz es la siguiente: Cepc =dσ σ σn+1 dε ε εn+1 =∂e σ σ σ ∂ε ε εn+1αn σ σ σn+1=e σ σ σ(αn,ε ε εn+1) (2.18) Donde la variable de endurecimiento queda fijada en la obtención de esta matriz, al igual que en la expresión de tensiones actualizadas para el algoritmo de Return Mapping (1.30) . De esta forma, podemos asemejar una ley elastoplástica (path-dependent) a una elástica no lineal con un módulo tangente (path-independent) (2.18) para el pseudo-tiempo [tn,tn+1] . Por tanto, no es necesario integrar todo el historial de tensiones, sino únicamente tener en cuenta el estado actual y su módulo tangente. Sin embargo, como ya se comentó en el capítulo anterior, es más sencillo obtener esta matriz en el espacio de tensiones principales a través de las expresiones de (1.30) y (1.31) . La relación entre (1.31) y (2.18) puede dividirse en dos partes. En primer lugar, se considera que (2.18) también está definida en el espacio de tensiones principales. En este caso, la equivalencia se produce a través de la definición de la deformación elástica de prueba ε ε εe trial n+1. A través del siguiente set de relaciones: ε ε εe trial n+1=ε ε εe n+∆ε ε ε ∆ε ε ε=ε ε εn+1−ε ε εn ε ε εn=ε ε εe n+ε ε εp n (2.19)
38 Capítulo 2. Implementación de la ley constitutiva Es fácilmente demostrable que ε ε εe trial n+1=ε ε εn+1−ε ε εp n . Por tanto, se puede establecer la siguiente equivalencia: σ σ σn+1=e σ σ σ(αn,εn+1)≡ˆ σ σ σ(αn,εe trial n+1)(2.20) Que se puede entender como un simple cambio de variable. Mediante este cambio se concluye: Cepc ppal =∂e σ σ σ ∂ε ε εn+1 =∂ˆ σ σ σ ∂ε ε εe trial n+1 (2.21) Donde se ha hecho implícito el uso de la regla de la cadena. La composición de esta matriz, para el paso n+1, queda de la siguiente forma: Cepc ppal = ∂σ1 ∂εe trial 1 ∂σ1 ∂εe trial 2 ∂σ1 ∂εe trial 3 ∂σ2 ∂εe trial 1 ∂σ2 ∂εe trial 2 ∂σ2 ∂εe trial 3 ∂σ3 ∂εe trial 1 ∂σ3 ∂εe trial 2 ∂σ3 ∂εe trial 3 (2.22) La segunda parte del proceso para llegar a (2.18) pasa por transformar la matriz (2.22) al espacio de tensiones globales. Se entiende por espacio global a la ordenación de las tensiones de la forma [σ11 σ22 σ33 σ12] , considerando problemas de deformación plana generalizada y axisimétrico. Para ello, se denota el hecho de que la matriz (2.22) calcula las derivadas de los autovalores del pseudovector σ σ σn+1 con respecto a los autovalores del pseudovector ε ε εe trial n+1 . De esta forma, conociendo los arrays de autovalores y las derivadas de un array con respecto al otro es bastante sencillo el cálculo de las derivadas de los pseudovectores de origen. El proceso de transformación se lleva a cabo mediante una subrutina que evalúa las derivadas del pseudovector de tensiones con respecto al pseudovector de deformaciones elásticas de prueba (véase la sección A.3 de [6] para más información). En resumen, tras conocer el valor de las tensiones principales, deformaciones elásticas de prueba y la fase del Return Mapping donde se ha alcanzado la convergencia, el proceso de obtención de Cepc pasa por obtener la matriz de derivadas en el espacio principal y posteriormente pasarla al espacio global. Todo este proceso, modelado para las distintas fases explicadas en el apartado anterior (fase elástica o plástica dividida en corrector al plano principal, bordes o vértice) se llevará a cabo a través de la subrutina MCCONTANMA, adaptada del código para HYPLAS de la rutina CTMC, de [ 6 ]. Cabe destacar que si estamos en fase elástica la matriz es igual a la elástica, por simple aplicación de las ecuaciones de Lamé. En cuanto al desarrollo de la subrutina, en primer lugar, se declaran las variables, parámetros y subrutinas internas a utilizar (Código 2.10). Código 2.10 Inicialización subrutina MCCONTANMA. Declaración de parámetros y variables. subroutine MCCONTANMA(& DMATX,LALGVA,RSTAVA, & STRAT,STRES,MSTRE,PROPS,PHI,PSI) !------------------------------------------------------------------------ ! Consistent tangent matrix for the Mohr-Coulomb model ! The tangent matrix computed in this routine is consistent with ! fully implicit elastic predictor/return mapping algorithm for the ! Mohr-Coulomb elasto-plastic material model. ! It contains the plane strain and axisymmetric implementations of ! the model. ! It returns either the elastic tangent or the elasto-plastic ! consistent tangent matrix depending on the input value of EPFLAG. !------------------------------------------------------------------------- implicit none !--- Parameters for subroutines -------------
2.1 Implementación y desarrollo del código 39 integer(4), parameter :: MCOMP=4, MDIM=3, NDIM=2 integer(4), parameter :: NTYPE=3 !----------------------- logical, intent(in) :: LALGVA(5) integer(4), intent(in) :: MSTRE real(8), intent(out) :: DMATX(MSTRE,MSTRE) real(8), intent(in) :: RSTAVA(MSTRE+1) real(8), intent(in) :: STRAT(MSTRE) real(8), intent(in) :: STRES(MSTRE) real(8), intent(in) :: PROPS(4),PHI,PSI !----------------------- logical :: EPFLAG, APEX, EDGE, OUTOFP,RIGHT logical :: REPEAT real(8) :: E, NU, COH0, H real(8) :: DPSTRS(MDIM,MDIM), PSTRS(MDIM), STRAC(MDIM) real(8) :: EIGPRJ(MSTRE,2), PSTRA(MDIM) real(8) :: FOID(MSTRE,MSTRE), SOID(MSTRE) real(8) :: GMODU, BULK, R2G, R4G, R2BULK, R2CPHI real(8) :: R4C2PH, R1D3, R2D3, R2GD3, R4GD3 integer(4) :: II, JJ, MM, I, J real(8) :: PSTMAX, PSTMIN real(8) :: SPHSPS, CONSTA, CONSTB real(8) :: FACTA, DRVAA, DRVAB, DRVBA, DRVBB real(8) :: AUX1, AUX2, AUX3 real(8) :: R1DDET real(8) :: COTPHI,DSIDEJ real(8) :: DENOM, B1, B2, B3 real(8) :: FACTOR integer(4) :: NSTRE !----------------------------------------------------------------------- !---- Explicit interfaces ---------------------- interface subroutine SPDEC2(EIGPRJ,EIGX,REPEAT,X, & MCOMP,NDIM) integer(4), intent(in) :: MCOMP,NDIM real(8), intent(in) :: X(MCOMP) real(8), intent(out) :: EIGPRJ(MCOMP,NDIM),EIGX(NDIM) logical, intent(out) :: REPEAT end subroutine SPDEC2 subroutine DGISO2(DEIGY,DYDX,EIGPRJ,EIGX,EIGY, & OUTOFP,REPEAT,MCOMP,MDIM,NDIM) integer(4), intent(in) :: MCOMP,MDIM,NDIM real(8), intent(in) :: DEIGY(MDIM,MDIM), EIGPRJ(MCOMP,NDIM) real(8), intent(in) :: EIGX(NDIM), EIGY(NDIM) logical, intent(in) :: OUTOFP, REPEAT real(8), intent(out) :: DYDX(MCOMP,MCOMP) end subroutine DGISO2 end interface Donde como variable nueva se presenta DMATX, que es el array que almacena la matriz elastoplástica consistente, análoga a la variable DDSDDE utilizada por ABAQUS . También se crea una variable lógica nueva, EPFLAG, análoga a IFPLAS en MCSTUP, que determina si la fase de actualización de tensiones ha sido elástica o plástica. La subrutina interna DGISO2 [ 6 ] se encargará de transformar la matriz de derivadas en principales a globales, de acuerdo con lo explicado anteriormente. La subrutina SPDEC2 se volverá a utilizar para obtener los autovalores de ε ε εe trial n+1.
40 Capítulo 2. Implementación de la ley constitutiva Al ser el proceso de obtención de la matriz bastante mecánico en todas las fases, se ilustra a continuación un ejemplo de cálculo de una de las componentes de esta en el espacio principal. Dicha componente será la (2,2) de la fase del borde izquierdo, es decir ∂σ2 ∂εe trial 2 . En primer lugar, se toman diferenciales de (1.30) para σ2 teniendo en cuenta que contamos con dos multiplicadores plásticos [∆γa∆γb] correspondientes a los vectores de flujo de los planos 1 y 2, N1yN2, quedando: dσ2=dσtrial 2+4 3G−2Ksinψd∆γa−2G1+1 3sinψ+2Ksinψd∆γb= =dσtrial 2+AUX1 d∆γa+AUX2 d∆γb (2.23) Por un lado, para obtener dσtrial 2 , o cualquier componente del vector de tensiones de prueba en principales, en función de dεe trial 2se hace uso de la ley elástica lineal, tal que: ∂σtrial 1 ∂εe trial 2 ∂σtrial 2 ∂εe trial 2 ∂σtrial 3 ∂εe trial 2 = K−2 3G K+4 3G K−2 3G (2.24) Por otro lado, para obtener [d∆γad∆γb] en función dεe trial 2 es necesario diferenciar el sistema de ecuaciones residual (2.10) del apartado anterior (las constantes y variables definidas en este desarrollo provienen de las utilizas en el sistema (2.10)): d˜ Φa(d∆γa,d∆γb) d˜ Φb(d∆γa,d∆γb) = d˜ σa+∂˜ Φa ∂∆γad∆γa+∂˜ Φa ∂∆γbd∆γb d˜ σb+∂˜ Φb ∂∆γad∆γa+∂˜ Φb ∂∆γbd∆γb ="0 0#(2.25) Que es equivalente al siguiente sistema: −∂˜ Φa ∂∆γa−∂˜ Φa ∂∆γb −∂˜ Φb ∂∆γa−∂˜ Φb ∂∆γb d∆γa d∆γb = d˜ σa d˜ σb (2.26) Donde: ∂˜ Φa ∂∆γa ∂˜ Φa ∂∆γb ∂˜ Φb ∂∆γa ∂˜ Φb ∂∆γb =−a−4Hcos2φ−b−4Hcos2φ −b−4Hcos2φ−a−4Hcos2φ ∂˜ σa ∂εe trial 2 ∂˜ σb ∂εe trial 2 = 2K−4 3G 2G+2K+2 3Gsinφ (2.27) En la que en la primera igualdad de (2.27) se ha hecho implícito el uso de la regla de la cadena para c(¯ εp n+∆¯ εp) utilizando la primera ecuación de (2.11) , mientras que en la segunda se ha hecho uso de (2.24) de nuevo. Despejando el set de incógnitas de (2.26) queda:
2.1 Implementación y desarrollo del código 41 d∆γa d∆γb =1 ∂˜ Φa ∂∆γa ∂˜ Φb ∂∆γb−∂˜ Φb ∂∆γa ∂˜ Φb ∂∆γb −∂˜ Φb ∂∆γb ∂˜ Φa ∂∆γb ∂˜ Φb ∂∆γa−∂˜ Φa ∂∆γa d˜ σa d˜ σb (2.28) Finalmente, sustituyendo la segunda componente de (2.24), (2.27) y (2.28) en la ecuación (2.23) diferenciando con respecto a dεe trial 2resulta: ∂σ2 ∂εe trial 2 =K+4 3G+1 D(AUX1"2K ∂˜ Φa ∂∆γb−∂˜ Φb ∂∆γb!+2G∂˜ Φa ∂∆γb+G 2 3 ∂˜ Φa ∂∆γb+4 3 ∂˜ Φb ∂∆γb!sinφ#+ +AUX2"2K ∂˜ Φa ∂∆γa−∂˜ Φb ∂∆γa!+2G∂˜ Φa ∂∆γa+G 2 3 ∂˜ Φa ∂∆γa+4 3 ∂˜ Φb ∂∆γa!sinφ#) (2.29) Donde D es el denominador (determinante) de (2.28) . Esta componente y el resto se encuentran implementadas en la siguiente parte del código (Código 2.11). Se muestra únicamente la fase del borde izquierdo ya que el procedimiento de resolución en el resto de fases del dominio plástico son análogas a la mostrada. Código 2.11 Desarrollo subrutina MCCONTANMA. Cálculo de la matriz elastoplástica consistente para la fase del borde izquierdo. IF(EPFLAG)THEN ! Compute elastoplastic consistent tangent ! ---------------------------------------- ! Spectral decomposition of the elastic trial strain STRAC(1)=STRAT(1) STRAC(2)=STRAT(2) STRAC(3)=STRAT(4)*0.5_8 CALL SPDEC2(EIGPRJ,PSTRA,REPEAT,STRAC,MCOMP,NDIM) PSTRA(3)=STRAT(3) ! and current total stress PSTRS(1)=STRES(1)*EIGPRJ(1,1)+STRES(2)*EIGPRJ(2,1)+2.0_8*STRES(4)* EIGPRJ(3,1) PSTRS(2)=STRES(1)*EIGPRJ(1,2)+STRES(2)*EIGPRJ(2,2)+2.0_8*STRES(4)* EIGPRJ(3,2) PSTRS(3)=STRES(3) ! Identify directions of maximum and minimum principal trial stresses ! II MAX ! JJ MIN ! MM MED II=1 JJ=1 PSTMAX=PSTRA(II) PSTMIN=PSTRA(JJ) DO I=2,3 IF(PSTRA(I).GE.PSTMAX)THEN II=I PSTMAX=PSTRA(II) ENDIF IF(PSTRA(I).LT.PSTMIN)THEN JJ=I PSTMIN=PSTRA(JJ) ENDIF ENDDO
42 Capítulo 2. Implementación de la ley constitutiva IF(II.NE.1.AND.JJ.NE.1)MM=1 IF(II.NE.2.AND.JJ.NE.2)MM=2 IF(II.NE.3.AND.JJ.NE.3)MM=3 IF(EDGE)THEN ! Tangent consistent with 2-vector return to edge SPHSPS=SIN(PHI)*SIN(PSI) CONSTA=R4G*(1.0_8+R1D3*SPHSPS)+4.0_8*BULK*SPHSPS IF(RIGHT)THEN CONSTB=R2G*(1.0_8+SIN(PHI)+SIN(PSI)-R1D3*SPHSPS)+4.0_8*BULK*SPHSPS ELSE CONSTB=R2G*(1.0_8-SIN(PHI)-SIN(PSI)-R1D3*SPHSPS)+4.0_8*BULK*SPHSPS ENDIF FACTA=R4C2PH*H DRVAA=-CONSTA-FACTA DRVAB=-CONSTB-FACTA DRVBA=-CONSTB-FACTA DRVBB=-CONSTA-FACTA AUX1=R2G*(1.0_8+R1D3*SIN(PSI))+R2BULK*SIN(PSI) AUX2=(R4GD3-R2BULK)*SIN(PSI) AUX3=R2G*(1.0_8-R1D3*SIN(PSI))-R2BULK*SIN(PSI) R1DDET=1.0_8/(DRVAA*DRVBB-DRVAB*DRVBA) ! ...returned to left edge DPSTRS(II,II)=BULK+R4GD3+(AUX1*(((R2BULK*(DRVBB-DRVAB)+(DRVAB*R4GD3+ DRVBB*R2GD3))*SIN(PHI))+DRVBB*R2G)+AUX2*(((R2BULK*(DRVBA-DRVAA) +(DRVAA*R4GD3+DRVBA*R2GD3))*SIN(PHI))+DRVBA*R2G))*R1DDET DPSTRS(II,MM)=BULK-R2GD3+(AUX1*(((R2BULK*(DRVBB-DRVAB)-(DRVAB*R2GD3+ DRVBB*R4GD3))*SIN(PHI))-DRVAB*R2G)+AUX2*(((R2BULK*(DRVBA-DRVAA) -(DRVAA*R2GD3+DRVBA*R4GD3))*SIN(PHI))-DRVAA*R2G))*R1DDET DPSTRS(II,JJ)=BULK-R2GD3+((AUX1*(DRVBB-DRVAB)+AUX2*(DRVBA-DRVAA)) *(((R2BULK+R2GD3)*SIN(PHI))-R2G))*R1DDET DPSTRS(MM,II)=BULK-R2GD3+(AUX1*(((R2BULK*(DRVAA-DRVBA)-(DRVAA*R4GD3+ DRVBA*R2GD3))*SIN(PHI))-DRVBA*R2G)+AUX2*(((R2BULK*(DRVAB-DRVBB) -(DRVAB*R4GD3+DRVBB*R2GD3))*SIN(PHI))-DRVBB*R2G))*R1DDET DPSTRS(MM,MM)=BULK+R4GD3+(AUX1*(((R2BULK*(DRVAA-DRVBA)+(DRVAA*R2GD3+ DRVBA*R4GD3))*SIN(PHI))+DRVAA*R2G)+AUX2*(((R2BULK*(DRVAB-DRVBB) +(DRVAB*R2GD3+DRVBB*R4GD3))*SIN(PHI))+DRVAB*R2G))*R1DDET DPSTRS(MM,JJ)=BULK-R2GD3+((AUX1*(DRVAA-DRVBA)+AUX2*(DRVAB-DRVBB)) *(((R2BULK+R2GD3)*SIN(PHI))-R2G))*R1DDET DPSTRS(JJ,II)=BULK-R2GD3+(AUX3*(((R2BULK*(DRVAB-DRVBB-DRVAA+DRVBA)+( DRVAA-DRVAB)*R4GD3+(DRVBA-DRVBB)*R2GD3)*SIN(PHI))+(DRVBA-DRVBB)* R2G))*R1DDET DPSTRS(JJ,MM)=BULK-R2GD3+(AUX3*(((R2BULK*(DRVAB-DRVBB-DRVAA+DRVBA)+( DRVAB-DRVAA)*R2GD3+(DRVBB-DRVBA)*R4GD3)*SIN(PHI))+(DRVAB-DRVAA)* R2G))*R1DDET DPSTRS(JJ,JJ)=BULK+R4GD3+(AUX3*(DRVAB-DRVBB-DRVAA+DRVBA)*(((R2BULK+ R2GD3)*SIN(PHI))-R2G))*R1DDET ! Take from principal to global space CALL DGISO2(DPSTRS,DMATX,EIGPRJ,PSTRA,PSTRS,& & OUTOFP,REPEAT,MCOMP,MDIM,NDIM) !-------------------------------------------------------------------------- end subroutine MCCONTANMA !------------------------------------------ !-------------------------------------------------------------------------- La componente (MM,MM) de la matriz DPSTRS corresponde al ejemplo resuelto anteriormente. Una vez calculada se llama a la subrutina interna DGISO2 para obtener la componente de la matriz elastoplástica en el espacio global.
2.2 Ejemplo de aplicación 43 2.2 Ejemplo de aplicación En este ejemplo se va a mostrar una de las capacidades más útiles por la cual se ha decido implementar el modelo de carácter elastoplástico de Mohr-Coulomb de [ 6 ] frente a otros. Dicha capacidad consiste en poder introducir una ley de endurecimiento propuesta de forma sencilla en la UMAT diseñada para este trabajo, con el fin de poder ajustarse a posteriores contrastes con resultados experimentales. Para poder ilustrar esta ventaja se van a aplicar diversos tipos de leyes de endurecimiento isótropo a un modelo simple de carácter cualitativo. El modelo propuestoconsiste en un elemento sujetoa las siguientes condiciones de contorno, esquematizado en la Figura 2.3. Figura 2.3 Esquema del modelo. La idea es obtener la curva del comportamiento en la dirección de aplicación de la condición de contorno en desplazamiento, y observar las distintas respuestas del modelo en función de la ley de endurecimiento establecida. La formulación general de la ley de endurecimiento, tal y como se describió en el Capítulo 1 para el modelo Mohr-Coulomb, es la siguiente: c=c(¯ εp)(2.30) Donde c es la cohesión del material (variable sujeta a evolución por endurecimiento) y ¯ εp la deformación plástica acumulada (variable de endurecimiento). A continuación, se detalla la lista de leyes de endurecimiento que se van a implementar, en las que la mayoría ya fueron mostradas a través de la Figura 2.1. Se realiza una clasificación en función de si la cohesión crece con la aplicación de la carga (endurecimiento), disminuye (reblandecimiento) o muestra un comportamiento mixto (mixto): ▼Leyes de endurecimiento: ■Lineal: c(¯ εp) = c0+H¯ εp ■Cuadrática: c(¯ εp) = c0+H(¯ εp)2 ■Exponencial: c(¯ εp) = c0·eH¯ εp ▼Leyes de reblandecimiento: ■Lineal: c(¯ εp) = c0−H¯ εp ■Inversa: c(¯ εp) = c0 1+H¯ εp ■Inversa al cuadrado: c(¯ εp) = c0 1+H(¯ εp)2
44 Capítulo 2. Implementación de la ley constitutiva ▼Ley mixta: ■Cuadrática: c(¯ εp) = a(¯ εp)2+b¯ εp+c En la mayoría de leyes propuestas se denota una dependencia con un par de parámetros que son c0 (cohesión inicial) y H (módulo de endurecimiento). En el caso de la ley mixta se opta por establecer una serie de coeficientes a,b,c; que podrían ser objeto de ajuste de una nube de puntos obtenida de resultados experimentales. Para este caso todos los parámetros y coeficientes tendrán un valor teórico, con el objetivo de mostrar resultados de índole cualitativo que se ajusten a un mismo rango de tensiones. Los datos de entrada para cada tipo de modelo son los siguientes: Tabla 2.1 Parámetros empleados para los modelos. Endurecimiento c0[kPa] H [kPa] Lineal 15.0 40.0 Cuadrático 15.0 100.0 Exponencial 15.0 2.0 Reblandecimiento c0[kPa] H [kPa] Lineal 15.0 -30.0 Inversa 15.0 15.0 I. al cuadrado 15.0 40.0 Mixto a [kPa] b [kPa] c [kPa] Cuadrático -2750/7 1675/14 15.0 En todos los casos se trabajará con valores de φ y ψ de carácter residual, con el objeto de que no interfieran en el cometido de este ejercicio, que es evaluar la forma de las curvas de comportamiento en función de las leyes de endurecimiento. Los valores asignados de los parámetros de la fase elástica son E=4.0·103 kPa y ν=0.3 . 2.2.1 Resultados Resolviendo el problema conforme a la UMAT planteada en la sección anterior, en la Figura 2.4 se muestra los resultados para cada tipo de ley de endurecimiento. 0 0.2 0.4 0.6 0.8 1 0 10 20 30 40 50 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 5 10 15 20 25 30 35 40 Figura 2.4 Curvas de comportamiento para cada tipo de ley de endurecimiento.
2.3 Ejemplo de verificación 45 Se observa que en todos los casos el tramo plástico de la curva de comportamiento se ajusta al carácter de la ley de endurecimiento modelada. Esto supone una amplia libertad para poder realizar ajustes experimentales, seleccionando el tipo de ley más adecuado y calibrando una serie de parámetros para acercarnos al comportamiento real del material que estemos estudiando. Como se deduce de la formulación mostrada en la sección anterior, el hecho de cambiar el tipo de ley de endurecimiento no afecta a la resolución del algoritmo de Return Mapping, ya que hacemos uso del Método de Newton-Raphson para resolver la ecuación residual, que es donde entra en juego la variable de cohesión. 2.3 Ejemplo de verificación En este ejemplo se va a plantear un problema de verificación cuyo objetivo es comparar las soluciones dadas por tres metodologías distintas, donde una de ellas es la planteada y explicada en este trabajo. A continuación se hace un pequeño resumen de cada una de ellas: ▼ Abaqus: Este software plantea un modelo propio formado por una superficie de plastificación tipo Mohr-Coulomb (1.25) y un potencial plástico distinto al descrito en este trabajo. Dicho potencial se basa en la formulación de Menétrey y William [52]. ▼ Clausen: Formulación verificada y basada en [ 3 ]. Está planteada para el caso elástico-perfectamente plástico, sin posibilidad de endurecimiento o reblandecimiento. El potencial plástico no asociativo es análogo al caso explicado en el apartado 2.1. ▼ Hyplas: Formulación adaptada al entorno de ABAQUS a través de la UMAT diseñada en este trabajo. Se permite además la modelización de endurecimiento o reblandecimiento a través de una ley de endurecimiento. El caso planteado se basa en la evaluación de la capacidad de carga de un material tipo suelo a través de la modelización del fenómeno de hundimiento de una zapata corrida infinitamente larga que descansa sobre un terreno homogéneo en condiciones ideales de de deformación plana, tal y como se muestra en la Figura 2.5. Figura 2.5 Esquema del problema. Cuyo resultado principal es la carga límite o de colapso del terreno P para las condiciones impuestas. Este problema ya ha sido resuelto por numerosos autores tal y como se detalla en [ 53 ], empleando el teorema de la cota superior. Una de las soluciones más famosas es la ideada por Prandtl & Hill [ 54 , 55 ], que predispone que la zapata es indeformable y que la carga superficial P se asemeja a un descenso uniforme del plano horizontal bajo la zapata, todo ello suponiendo un comportamiento elástico perfectamente plástico. Por otro lado, dentro de la solución de Prandtl & Hill se incluye la definición de la línea de rotura, que conforma el mecanismo de colapso. Dicho mecanismo de colapso se ilustra en la Figura 2.6.
52 Capítulo 3. Desarrollo del modelo de fractura acoplado El esquema de trabajo seguido por la UEL en cuanto al cálculo de contribuciones se detalla a través del diagrama de la Figura 3.1. Tras la llamada a la UEL para cada elemento, se obtiene el valor de las funciones de forma, matriz de deformación (derivada de las funciones de forma) y Jacobiano de la transformación (matriz de paso de coordenadas locales a naturales) en los puntos de Gauss. A través de la estimación de la variable de daño se calcula su valor y su gradiente en estos puntos. Tras calcular el campo d se llama a la UMAT para realizar el cometido explicado al principio de este capítulo. Cabe destacar que los resultados obtenidos son de carácter nominal, es decir, el estado tensional y matriz de comportamiento están afectados de la variable de daño a través de la función de degradación (1.9) . Además, dentro de la UMAT también se calculan las energías elástica y plástica necesarias para la obtención de la Driving Force. En cuanto a este último cálculo, por un lado, hay diversas formas de calcular la energía elástica (término E de (1.1) sin degradar), que para materiales isótropos y en base a las constantes de Lamé toma la forma: Ψe=λ 2(tr ε ε εe)2+µtr(ε ε εe)2(3.1) Se puede observar que (3.1) es combinación lineal de los invariantes del tensor de deformaciones elástico. A través de esta formulación se exponen las descomposiciones definidas en [59]: ▼ Descomposición de Miehe [ 24 ]: Ya introducida en (1.8) , se basa en la descomposición en la parte positiva (tracción) y negativa (compresión) haciendo uso de las deformaciones elásticas del sistema principal, que son los autovalores del tensor de deformaciones elástico (εe 1,εe 2,εe 3). Ψe±=λ 2⟨εe 1+εe 2+εe 3⟩2 ±+µ⟨εe 1⟩2 ±+⟨εe 2⟩2 ±+⟨εe 3⟩2 ±(3.2) ▼ Descomposición de Amor et al (parte esférica-desviadora): En este caso la descomposición se hace en base a la parte esférica y la parte desviadora del tensor de deformaciones elástico: ε ε εe=εe kI I I+ε ε εe des (3.3) Donde εe k es la deformación esférica media elástica, obtenida a partir de la traza del tensor. Haciendo uso de los parámetros K, módulo de compresibilidad, y G, módulo de cizalladura, la descomposición queda de la siguiente forma: Ψe=K 2(tr ε ε εe)2+Gtr(ε ε εe des)2(3.4) Que puede adoptar a su vez una descomposición en parte positiva y negativa: Ψe+=K 2⟨tr ε ε εe⟩2 ++Gtr(ε ε εe des)2 Ψe−=K 2⟨tr ε ε εe⟩2 − (3.5) La descomposición en parte positiva y negativa es útil ya que es más realista degradar únicamente la parte positiva, que es la encargada de la nucleación y desarrollo de defectos en el material. En este trabajo se hará uso del primer tipo de descomposición. Por otro lado, la energía disipativa plástica Ψpd (término H+P de (1.1) sin degradar) se calcula a través de la obtención del área bajo la curva de tensión-deformación plástica, de forma que: Ψpd =Ψpd + n ∑ i=1 σn+1 i+σn i 2(εn−1 pi−εn pi)(3.6) Siendo n el número de componentes del pseudovector de deformaciones plásticas, en el caso 2-D serían tres y en el caso 3-D seis. Tal y como se muestra en la Figura 3.2 para una componente genérica, la energía se calcula mediante la discretización de la integral o suma discreta de áreas sucesivas, añadiendo en cada paso la encerrada entre los valores de deformación plástica del paso anterior y del actual.
3.2 Análisis de metodologías 53 Figura 3.2 Obtención de la energía disipativa plástica a través de suma de áreas para una componente. Se observa que este tipo de disipación se calcula en función de las tensiones efectivas, que no están afectadas del campo daño. Se emplean este tipo de tensiones, y no las nominales (afectadas del campo de daño), ya que la energía se degrada posteriormente a través de la función de degradación. Tras la obtención de todas las variables necesarias, se añaden como contribuciones a los vectores de residuos y matrices de rigidez, a través de las expresiones definidas en el bloque de Staggered Algorithm del Capítulo 1. Dichas contribuciones se actualizan hasta alcanzar la convergencia establecida mediante las tolerancias impuestas en el modelo. 3.2 Análisis de metodologías Tras quedar claro el proceso de obtención de las componentes del sistema de ecuaciones del MEF a resolver, a continuación haremos especial énfasis en el cálculo del residuo y parte de la matriz de rigidez dependiente del campo de daño d , rd i y Kdd ij . Recordemos que estas variables dependían, además de la discretización del MEF realizada, del término de Driving Force de (1.18) a través de la variable histórica H y la función de degradación g(d). Como se comentó en su momento con referencia en [ 28 ], el término de Driving Force es objeto de cambios para el desarrollo de distintas metodologías de carácter acoplado como los mostrados en la Tabla 1.2. En este trabajo se estudiará la derivación de las metodologías de Borden et al. [ 20 ] y Miehe et al. [ 19 ], que son los que más versatilidad ofrecen a la hora de modelar este tipo de comportamiento en materiales. En ambos casos, se parte de la ecuación de evolución de la superficie de grieta (1.18) realizando una serie de modificaciones a través de parámetros y thresholds. A continuación se muestra como obtener las componentes del vector de residuos y matriz de rigidez asociados al daño para cada metodología. ▼Metodología propuesta por Borden et al. [20]: En este caso, la ecuación de evolución sufre cambios en la definición de H y g(d) , partiendo de (1.4) , quedando de la siguiente forma: Gc lcd−l2 c∆d=−g′(d)βeH+βp⟨Ψpd −W0⟩+(3.7)
54 Capítulo 3. Desarrollo del modelo de fractura acoplado Donde: H=m´ ax s∈[0,t][Ψe(ε ε εe)] g(d) =(1−d)2(1+d(2−m)),m>0 (3.8) Siendo nomenclatura asignada es la siguiente: ■Gc: Tenacidad a fractura. ■lc: Longitud característica. ■g(d): Función de degradación. ■βe , βp : Coeficientes de influencia de la parte de la energía elástica y plástica respectivamente en la Driving Force. ■H: Variable histórica del término de Driving Force. ■Ψpd: Energía disipativa plástica. ■W0:Threshold de la parte de la energía plástica. Por facilitar la comprensión en el desarrollo, se ha obviado intencionadamente la descomposición en parte positiva y negativa de la energía elástica en ambas metodologías, teniendo solo que afectar la función de degradación a la parte positiva. En base a la discusión ofrecida en [ 22 ], Borden y coautores utilizan una función de degradación cúbica con AT-2 para poder generar un tramo lineal previo a alcanzar la deformación correspondiente a la tensión pico (como comparación al comportamiento de la de carácter cuadrático consultar [ 20 ]) y así suplir el defecto del comienzo inmediato del daño con la puesta en carga propia de la combinación { g(d) cuadrática+AT-2}. En la Figura 3.3 se compara el uso de ambas funciones para la respuesta a carga uniaxial correspondiente a un caso homogéneo planteado en [ 20 ]. Se observa que antes de alcanzar la tensión pico encontramos un tramo lineal para la función cúbica, mientras que la cuadrática se aleja de un comportamiento elástico-lineal. Para conseguir mantener el mismo nivel de deformaciones εcuad d al alcanzar la tensión pico σcub d que en el caso cuadrático y mantener la semejanza con un modelo AT-2, el parámetro m se 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ón de defectos para el estado intacto del material. En [20], siguiendo diversas modelizaciones, se opta por un valor de m=10−4. Figura 3.3 Curva de comportamiento del caso comparativo para carga uniaxial extraída de [20]. Por otro lado, Borden separa la Driving Force en dos términos afectados por unos coeficientes de ponderación βe y βp . El primer término afecta a la energía elástica y el segundo a la disipación plástica controlada además por un threshold W0 . Los coeficientes de ponderación modelan la contribución de cada término al crecimiento de grieta, que es el fin último de la Driving Force.
3.2 Análisis de metodologías 55 Para obtener la componente asociada al daño del vector de residuos y matriz de rigidez es necesario utilizar la discretización propia del MEF: d=N N Nid d di ∇xd=B B Bid d di (3.9) Por la cual, sustituyendo en (3.7) y derivando se obtiene: r r rd i=ZBGc lc d+g′(d)βeH+βp⟨Ψpd −W0⟩+N N Ni+GclcB B BT i∇xddV K K Kdd i j =∂rd i ∂dj =ZBGc lc +g′′(d)βeH+βp⟨Ψpd −W0⟩+N N NiN N Nj+GclcB B BT iB B BjdV (3.10) ▼Metodología propuesta Miehe y coautores [19]: Por su parte, Miehe et al. [ 25 ] realizan la modificación del término F de (1.1) . En dicho término establece una extensión del modelo AT-2: ∆f(d,∇xd):=2ωc ξd2 lc +lc|∇xd|2+(1−g(d))ωc(3.11) Cuya nomenclatura es la siguiente: ■lc: Longitud característica. ■g(d): Función de degradación. ■ωc: Energía crítica de fractura. ■ξ: Parámetro de control de la fase postcrítica. Básicamente, se añade un término adicional controlado por la energía crítica de fractura ωc que actúa a modo de threshold. Desarrollando la energía total del sistema y derivando con respecto a la variable daño para obtener la ecuación de evolución queda: 2ωc ξd−l2 c∆d=−g′(d)H(3.12) Donde: H=m´ ax s∈[0,t]⟨Ψe+Ψpd −ωc⟩ g(d) =(1−d)2(3.13) Como puede observarse, esta formulación introduce el threshold de manera natural en el término de energía de fractura, que deriva en la ecuación de evolución, a diferencia de otros autores que optan por establecerlo en los términos de energía elástica o plástica. De esta forma, consigue asemejarse a un modelo AT-1 con un tramo previo lineal sin daño apreciable, pero con las ventajas propias de un modelo AT-2. Además, emplea un parámetro ξpara controlar la caída de la fase postcrítica de daño. De forma análoga al proceso realizado en Borden et al. [ 20 ], se obtienen las contribuciones al sistema del MEF: r r rd i=ZB2ωc ξd+g′(d)HN N Ni+2ωcl2 c ξB B BT i∇xddV K K Kdd i j =∂rd i ∂dj =ZB2ωc ξ+g′′(d)HN N NiN N Nj+2ωcl2 c ξB B BT iB B BjdV (3.14)
56 Capítulo 3. Desarrollo del modelo de fractura acoplado 3.3 Ejemplo de aplicación Para verificar las metodologías presentadas en la sección anterior con el modelo computacional desarrollado, se comprobará de forma cualitativa su efectividad a través de un ejemplo sencillo de carácter homogéneo. Dicho ejemplo supone la modelización de un elemento sometido a tracción en condiciones de deformación plana. El elemento utilizado es de tipo bilineal de 4 nodos, con 4 puntos de integración. En la Figura 3.4 se ilustra el esquema del problema planteado, análogo al planteado en la sección 2.2, detallando las condiciones de contorno impuestas. Figura 3.4 Esquema del modelo. Dentro de la modelización de la plasticidad a partir del modelo de Mohr-Coulomb se propone la siguiente ley de endurecimiento de carácter lineal para la cohesión, en función de la deformación plástica acumulada ¯ εp(variable de endurecimiento): c(¯ εp) = c0+H¯ εp(3.15) Las propiedades empleadas en la modelización son las recogidas en la Tabla 3.1. Estas propiedades modelizan de forma cualitativa un suelo tipo arcilloso en condiciones no drenadas (simulación a corto plazo), material característico de aplicación del modelo Mohr-Coloumb, con un módulo de endurecimiento H negativo ya que la tracción va asociada a una pérdida de resistencia en suelos. Tabla 3.1 Propiedades y parámetros utilizados para las metodologías de Borden et al. y Miehe et al.. 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
3.3 Ejemplo de aplicación 57 3.3.1 Resultados A continuación, se incluye un análisis pormenorizado de los resultados para ambas metodologías: ▼Metodología de Borden et at.: ■Análisis del parámetro m: En la Figura 3.5 y la Figura 3.6 se compara la influencia del parámetro m con respecto al estado no dañado. En ninguno de los casos se aplica un threshold. Mientras que en el primer caso, correspondiente al empleo de m=10−4 , se verifica la conservación de una fase no dañada antes de darse la caída de resistencia, en el caso de m=2 (correspondiente a un modelo AT-2, compruébese sustituyendo en la segunda ecuación de (3.8) ) la degradación comienza desde la fase elástica. Este hecho demuestra una de las ventajas de la metodología de Borden et al., no necesitando de un threshold para evitar el comienzo del daño desde el comienzo del estado de carga. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 25 30 Figura 3.5 Comparativa del empleo de m=10−4con el caso no dañado. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 25 30 Figura 3.6 Comparativa del empleo de m=2con el caso no dañado.
58 Capítulo 3. Desarrollo del modelo de fractura acoplado En la Figura 3.7 se compara la metodología de Borden et al. (m=10−4) con un modelo AT-2 con función de degradación cuadrática (m=2) . Uno de los aspectos más interesantes a comentar es la transición más suave a la caída del primer caso con respecto al segundo, debido al empleo de una función de degradación cúbica, en vez de una cuadrática. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 25 30 0.25 0.27 25 28 Figura 3.7 Comparativa del empleo de m=10−4(Borden et al.) con m=2(g(d)cuadrática+AT-2). ■Análisis del parámetro W0: Como refuerzo del primer análisis realizado en esta sección, en la Figura 3.8 se demuestra que la metodología de Borden et al. no necesita de threshold para desarrollar una fase sin daño previa. También se denota el crecimiento de la variable de daño d hasta un valor cercano a 1, aunque nos encontremos ante un caso homogéneo. Este hecho indica que este tipo de metodología desarrolla un fuerte crecimiento del daño una vez tiene comienzo. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 25 30 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figura 3.8 Resultados de las curvas de comportamiento (línea continua) y evolución del daño (línea discontinua) para la comparativa de Borden et al. sin threshold y el caso no dañado.
3.3 Ejemplo de aplicación 59 En la Figura 3.9 se ilustra una comparativa del valor del threshold para distintos casos. Como era de esperar, a mayor valor del threshold más se retrasa el inicio de la fase de daño en el elemento. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 25 30 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figura 3.9 Resultados de las curvas de comportamiento (línea continua) y evolución del daño (línea discontinua) para la comparativa de distintos valores de threshold. ▼Metodología de Miehe et al.: ■Análisis del parámetro ωc: En la metodología de Miehe et al. se emplea el parámetro ωc para establecer la energía crítica desencadenante del proceso de daño (similar a Gc ) y establecer un threshold, tal y como se ilustró en (3.11) . A través de la Figura 3.10 se observa como el parámetro ωc establece una fase previa al daño, ajustable a través de su valor. En contraste con la metodología de Borden, se denota la caída más suave, donde la variable de daño d alcanza aproximadamente un valor de 0.8. Más adelante se verá que la forma de la caída se puede ajustar gracias al parámetro ξ. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 25 30 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Figura 3.10 Resultados de las curvas de comportamiento (línea continua) y evolución del daño (línea discontinua) para la comparativa de Miehe et al. con ωc=1.0y el caso no dañado.
60 Capítulo 3. Desarrollo del modelo de fractura acoplado De forma análoga a la metodología de Borden et al., en la Figura 3.11 se comprueba el efecto del crecimiento de la energía crítica ωc . A mayor valor de este parámetro, mayor es la energía necesaria para desencadenar la etapa de daño. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 25 30 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figura 3.11 Resultados de las curvas de comportamiento (línea continua) y evolución del daño (línea discontinua) para la comparativa de distintos valores de ωc. ■Análisis del parámetro ξ: Por último, en la Figura 3.12 se ilustran las diferentes formas de caída propias de la fase postcrítica debido al parámetro ξ . A medida que aumenta el valor de este parámetro la pendiente va creciendo, asemejándose a una fase de daño sin mucho recorrido en deformaciones. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 25 30 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figura 3.12 Resultados de las curvas de comportamiento (línea continua) y evolución del daño (línea discontinua) para la comparativa de distintos valores de ξ.
4 Conclusión En el presente trabajo se ha realizado un estudio en profundidad de la mecánica de la fractura en el contexto del comportamiento dúctil, que ha derivado en la aplicación de leyes no lineales del material a partir de la modelización de la plasticidad con posibilidad de ser de carácter asociativa o no asociativa. 4.1 Conclusiones En esta sección se describe los desarrollos llevados a cabo en cada capítulo de este trabajo: ▼ En el primer capítulo, a modo de introducción, se ha realizado una revisión del estado del arte de la mecánica de la fractura y la plasticidad. En el apartado de la mecánica de la fractura, se ha descrito el comportamiento tanto frágil como dúctil en materiales, desde una perspectiva de la mecánica computacional a través de metodologías de tipo variacional como es el PFAF. Mediante esta metodología, se ha podido establecer la versión regularizada del modelo por el cual se obtienen los algoritmos a emplear para la resolución del problema. Centrándonos en el comportamiento dúctil, se han dispuesto las herramientas y conocimientos necesarios para desarrollar la labor de este trabajo. Por otro lado, en el apartado de plasticidad, se han descrito los fundamentos y desarrollos necesarios para el empleo de técnicas como el Return Mapping para la actualización del estado tensional y la matriz de comportamiento elastoplástico consistente. De forma particular, con miras al próximo capítulo, se ha definido el Modelo de Mohr-Coulomb, con base en [6]. ▼ En el segundo capítulo, se ha desarrollado el modelo computacional de la plasticidad a través de la creación de una subrutina UMAT utilizando el programa de elementos finitos ABAQUS . Para demostrar la eficacia y validez del modelo, se han desarrollado un par de ejemplos de aplicación. ▼ En el tercer capítulo, se ha mostrado la capacidad de varias metodologías descritas en el primer capítulo empleando el modelo computacional desarrollado. En este caso, el carácter asociativo o no asociativo de la plasticidad viene dada por una aplicación sencilla del modelo de Mohr-Coulomb. Gracias al problema planteado en este capítulo, se ha podido comprobar como los parámetros de cada metodología influyen en la solución. Con este conocimiento se podría trabajar con el ajuste de resultados experimentales. Sin embargo, cabe destacar que detrás de cada metodología y de la obtención del valor de cada parámetro debería existir un trasfondo físico que le aporte rigor científico y experimental a la modelización efectuada. Finalmente, se denota la versatilidad del modelo computacional para incluir otro tipo de leyes de comportamiento distintas a través del desarrollo de UMAT específicas, que pueden acoplarse fácilmente al entorno de la UEL implementado en ABAQUS . Además, una vez verificado los resultados para un elemento, los problemas planteados pueden asemejarse a geometrías y condiciones más realistas. 4.2 Líneas de investigación futuras En base a los conocimientos adquiridos con la realización de este trabajo, como líneas de investigación se proponen: 61
Bibliografía [1] S. Teichtmeister, D. Kienle, F. Aldakheel, y M. A. Keip, “Phase field modeling of fracture in anisotropic brittle solids,” International Journal of Non-Linear Mechanics, vol. 97, 2017. [2] “Apuntes de Mecánica de la Fractura,” Grupo de Elasticidad y Resistencia de Materiales, Departamento de Mecánica de Medios Continuos y Teoría de Estructuras, ETSI, US, Curso 2019/2020. [3] J. Clausen, Efficient Non-linear Finite Element Implementation of Elasto-Plasticity for Geotechnical Problems : Thesis. Computational Mechanics. Esbjerg Institute of Technology. Aalborg University Esbjerg, 2007. [4] A. A. Griffith, “The phenomena of rupture and flow in solids,” Masinovedenie, 1995. [5] G. A. Francfort y J. J. Marigo, “Revisiting brittle fracture as an energy minimization problem,” Journal of the Mechanics and Physics of Solids, vol. 46, 1998. [6] E. A. D. S. Neto, D. Perić, y D. R. Owen, Computational Methods for Plasticity: Theory and Applications, 2008. [7] J. Lubliner, “Plasticity Theory: Revised Edition,” Journal of Applied Mechanics, vol. 59, 2006. [8] R. I. Borja, Plasticity Modeling and Computation, 2013. [9] E. Hoek, C. Carranza, y B. Corkum, “Hoek-brown failure criterion – 2002 edition,” Narms-Tac, 2002. [10] E. Hoek y E. T. Brown, “The Hoek–Brown failure criterion and GSI – 2018 edition,” Journal of Rock Mechanics and Geotechnical Engineering, vol. 11, 2019. [11] J. Clausen y L. Damkilde, “An exact implementation of the Hoek-Brown criterion for elasto-plastic finite element calculations,” International Journal of Rock Mechanics and Mining Sciences, vol. 45, 2008. [12] R. Alessi, J. J. Marigo, C. Maurini, y S. Vidoli, “Coupling damage and plasticity for a phase-field regularisation of brittle, cohesive and ductile fracture: One-dimensional examples,” International Journal of Mechanical Sciences, vol. 149, 2018. [13] C. Miehe, M. Hofacker, y F. Welschinger, “A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits,” Computer Methods in Applied Mechanics and Engineering, vol. 199, 2010. [14] M. A. Msekh, J. M. Sargado, M. Jamshidian, P. M. Areias, y T. Rabczuk, “Abaqus implementation of phase-field model for brittle fracture,” Computational Materials Science, vol. 96, pp. 472–484, 1 2015. [15] M. Paggi y J. Reinoso, “Revisiting the problem of a crack impinging on an interface: A modeling framework for the interaction between the phase field approach for brittle fracture and the interface cohesive zone model,” Computer Methods in Applied Mechanics and Engineering, vol. 321, 2017. [16] F. P. Duda, A. Ciarbonetti, P. J. Sánchez, y A. E. Huespe, “A phase-field/gradient damage model for brittle fracture in elastic-plastic solids,” International Journal of Plasticity, vol. 65, 2014. 69
70 Bibliografía [17] M. Ambati, T. Gerasimov, y L. D. Lorenzis, “Phase-field modeling of ductile fracture,” Computational Mechanics, vol. 55, pp. 1017–1040, 5 2015. [18] C. Kuhn, T. Noll, y R. Müller, “On phase field modeling of ductile fracture,” GAMM Mitteilungen, vol. 39, 2016. [19] C. Miehe, S. Teichtmeister, y 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, vol. 374, 2016. [20] M. J. Borden, T. J. Hughes, C. M. Landis, A. Anvari, y I. J. Lee, “A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects,” Computer Methods in Applied Mechanics and Engineering, vol. 312, 2016. [21] C. Miehe, F. Welschinger, y M. Hofacker, “Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations,” International Journal for Numerical Methods in Engineering, vol. 83, 2010. [22] R. Alessi, M. Ambati, T. Gerasimov, S. Vidoli, y L. D. Lorenzis, “Comparison of phase-field models of fracture coupled with plasticity,” Computational Methods in Applied Sciences, vol. 46, 2018. [23] G. A. Francfort, B. Bourdin, y J. J. Marigo, “The variational approach to fracture,” Journal of Elasticity, vol. 91, pp. 5–148, 4 2008. [24] C. Miehe, L. M. Schänzel, y H. Ulmer, “Phase field modeling of fracture in multi-physics 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, vol. 294, 2015. [25] C. Miehe, M. Hofacker, L. M. Schänzel, y 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, vol. 294, 2015. [26] P. Rodriguez, J. Ulloa, C. Samaniego, y E. Samaniego, “A variational approach to the phase field modeling of brittle and ductile fracture,” International Journal of Mechanical Sciences, vol. 144, 2018. [27] C. Li, J. Fang, C. Wu, G. Sun, G. Steven, y Q. Li, “Phase field fracture in elasto-plastic solids: Incorporating phenomenological failure criteria for ductile materials,” Computer Methods in Applied Mechanics and Engineering, vol. 391, p. 114580, 2022. [28] B. Yin y M. Kaliske, “A ductile phase-field model based on degrading the fracture toughness: Theory and implementation at small strain,” Computer Methods in Applied Mechanics and Engineering, vol. 366, 2020. [29] J. Y. Wu, V. P. Nguyen, H. Zhou, y Y. Huang, “A variationally consistent phase-field anisotropic damage model for fracture,” Computer Methods in Applied Mechanics and Engineering, vol. 358, 2020. [30] M. Jirásek, “Nonlocal models for damage and fracture: Comparison of approaches,” International Journal of Solids and Structures, vol. 35, 1998. [31] J. Lemaitre, “A continuous damage mechanics model for ductile fracture,” Journal of Engineering Materials and Technology, Transactions of the ASME, vol. 107, 1985. [32] J. Lemaitre y R. Desmorat, Engineering damage mechanics: Ductile, creep, fatigue and brittle failures, 2005. [33] G. T. Camacho y M. Ortiz, “Computational modelling of impact damage in brittle materials,” International Journal of Solids and Structures, vol. 33, 1996. [34] M. Jirásek, “Comparative study on finite elements with embedded discontinuities,” Computer Methods in Applied Mechanics and Engineering, vol. 188, 2000. [35] N. Sukumar, N. Moës, B. Moran, y T. Belytschko, “Extended finite element method for three-dimensional crack modelling,” International Journal for Numerical Methods in Engineering, vol. 48, 2000.
Bibliografía 71 [36] M. Stolarska, D. L. Chopp, N. Mos, y T. Belytschko, “Modelling crack growth by level sets in the extended finite element method,” International Journal for Numerical Methods in Engineering, vol. 51, 2001. [37] E. Béchet, H. Minnebo, N. Moës, y B. Burgardt, “Improved implementation and robustness study of the X-FEM for stress analysis around cracks,” International Journal for Numerical Methods in Engineering, vol. 64, 2005. [38] J. Oliver, A. E. Huespe, y P. J. Sánchez, “A comparative study on finite elements for capturing strong discontinuities: E-FEM vs X-FEM,” Computer Methods in Applied Mechanics and Engineering, vol. 195, 2006. [39] T. Yang, K. M. Liechti, y R. Huang, “A multiscale cohesive zone model for rate-dependent fracture of interfaces,” Journal of the Mechanics and Physics of Solids, vol. 145, 2020. [40] N. Chandra, H. Li, C. Shet, y H. Ghonem, “Some issues in the application of cohesive zone models for metal-ceramic interfaces,” International Journal of Solids and Structures, vol. 39, 2002. [41] D. Xie y A. M. Waas, “Discrete cohesive zone model for mixed-mode fracture using finite element analysis,” Engineering Fracture Mechanics, vol. 73, 2006. [42] H. D. Espinosa, P. D. Zavattieri, y S. K. Dwivedi, “A finite deformation continuum/discrete model for the description of fragmentation and damage in brittle materials,” Journal of the Mechanics and Physics of Solids, vol. 46, 1998. [43] J. Reinoso y M. Paggi, “A consistent interface element formulation for geometrical and material nonlinearities,” Computational Mechanics, vol. 54, 2014. [44] M. Paggi y J. Reinoso, “An anisotropic large displacement cohesive zone model for fibrillar and crazing interfaces,” International Journal of Solids and Structures, vol. 69-70, 2015. [45] R. de Borst y C. V. Verhoosel, “Gradient damage vs phase-field approaches for fracture: Similarities and differences,” Computer Methods in Applied Mechanics and Engineering, vol. 312, 2016. [46] J. Ulloa, J. Wambacq, R. Alessi, E. Samaniego, G. Degrande, y S. François, “A micromechanics-based variational phase-field model for fracture in geomaterials with brittle-tensile and compressive-ductile behavior,” Journal of the Mechanics and Physics of Solids, vol. 159, 2022. [47] J. Reinoso, “Análisis numérico de procesos de fractura en rocas heterogéneas mediante el uso combinado de la metodología phase field y grietas de interfase,” Trabajo Fin de Máster, Universidad de Sevilla, 2020. [48] J. Fang, C. Wu, J. Li, Q. Liu, C. Wu, G. Sun, y Q. Li, “Phase field fracture in elasto-plastic solids: Variational formulation for multi-surface plasticity and effects of plastic yield surfaces and hardening,” International Journal of Mechanical Sciences, vol. 156, 2019. [49] R. Alessi, J. J. Marigo, y S. Vidoli, “Gradient damage models coupled with plasticity: Variational formulation and main properties,” Mechanics of Materials, vol. 80, 2015. [50] J. Fang, C. Wu, T. Rabczuk, C. Wu, C. Ma, G. Sun, y Q. Li, “Phase field fracture in elasto-plastic solids: Abaqus implementation and case studies,” Theoretical and Applied Fracture Mechanics, vol. 103, 2019. [51] R. Picón, “Teoría de la Plasticidad,” Grupo de Elasticidad y Resistencia de Materiales, Departamento de Mecánica de Medios Continuos y Teoría de Estructuras, ETSI, US, Curso 2019/2020. [52] P. Menetrey y K. J. Willam, “Triaxial failure criterion for concrete and its generalization,” ACI Structural Journal, vol. 92, 1995. [53] W. F. Chen, Limit Analysis and Soil Plasticity, 1975. [54] L. Prandtl, “Hauptaufsätze: Über die Eindringungsfestigkeit (Härte) plastischer Baustoffe und die Festigkeit von Schneiden,” ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, vol. 1, 1921.
72 Bibliografía [55] R. Hill, “The plastic yielding of notched bars under tension,” Quarterly Journal of Mechanics and Applied Mathematics, vol. 2, 1949. [56] “Guía de Cimentaciones en Obras de Carretera,” Ministerio de Fomento, Gobierno de España, 2009. [57] W. F. Chen, A. F. Saleeb, y G. J. Dvorak, “Constitutive Equations for Engineering Materials, Volume I: Elasticity and Modeling,” Journal of Applied Mechanics, vol. 50, 1983. [58] K. Terzaghi, Theoretical Soil Mechanics, 1943. [59] N. P. Dijk, J. J. Espadas-Escalante, y P. Isaksson, “Strain energy density decompositions in phase-field fracture theories for orthotropy and anisotropy,” International Journal of Solids and Structures, vol. 196-197, 2020. [60] D. Kienle, F. Aldakheel, y M. A. Keip, “A finite-strain phase-field approach to ductile failure of frictional materials,” International Journal of Solids and Structures, vol. 172-173, 2019. [61] E. E. Alonso, A. Gens, y A. Josa, “A constitutive model for partially saturated soils,” Geotechnique, vol. 40, 1990.