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

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 Proyecto Fin de Máster Máster en Ingeniería de Caminos, Canales y Puertos 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 Autor: José Reinoso Cuevas Tutor: Percy Durand Neyra Dpto. Estructuras de Edificación e Ingeniería del Terreno Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 Proyecto Fin de Máster Máster en Ingeniería de Caminos, Canales y Puertos 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 Autor: José Reinoso Cuevas Tutor: Percy Durand Neyra Profesor Titular de Universidad Dpto. Estructuras de Edificación e Ingeniería del Terreno Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 Proyecto Fin de Máster: 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 Autor: José Reinoso Cuevas Tutor: Percy Durand Neyra 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: Agradecimientos Este trabajo es la culminación de un objetivo personal que se remonta bastante tiempo atrás, pero que sin el apoyo de muchas personas hubiera sido realmente imposible. Por este motivo, no puedo dejar de recordar a las personas que con su apoyo y aliento incansable y sincero me han ayudado a llegar hasta aquí. En el ámbito académico me gustaría agradecer a los profesores que han contribuido a mi formación intelectual en diferentes ámbitos. Desde los que me recibieron en los primeros años de mis estudios de Ingeniería Industrial hasta los que he tenido la suerte de tener en una etapa posterior en mi formación como Ingeniero de Caminos, Canales y Puertos. Sin embargo, me gustaría aprovechar estas líneas para de veras agradecer a las personas más importantes de mi vida y que como no puede ser de otra forma son mi familia y amigos. A mis padres y a mis hermanas, les diría que "lo siento", tengo un examen para tal día. Es casi lo que ha sido mi excusa o frase recurrente en los últimos dos años. Os he robado demasiado tiempo y pido disculpas por ello, pero también expresaros mi profundo agradecimiento por creer siempre en mí. Esta forma de pedir perdón y agradecer al mismo tiempo es extensible a mis amigos, quienes me han sufrido bastante y que siempre me han apoyado también. Pero para ser honestos, el poder terminar esta nueva andadura se lo debo a Rocío. Su paciencia a prueba de bombas, su compresión y todo el cariño mostrado durante todos estos años ha sido una fuente inagotable para no desistir en el camino. Estos valores me han impulsado hasta límites que ni siquiera yo era consciente que pudieran ser alcanzados por alguien con mis limitaciones. Por todo ello, mis gracias se convierten en infinitas. Finalmente, como siempre en mi recuerdo, agradecer a la persona más íntegra que he tenido la suerte de conocer en mi vida: Caridad García Ávila, siempre presente en mi corazón. I Resumen El análisis de los fenómenos de fractura en rocas tiene implicaciones relevantes en diferentes aplicaciones en Ingeniería Civil y Geológica. El comportamiento complejo de tales materiales generalmente abarca la evolución de diferentes escenarios que implican ramificación de grietas, coalescencia, entre muchos otros. En este trabajo, la aplicación del enfoque variacional de fractura phase field se ha aplicado con éxito para la simulación de procesos de grietas en rocas, incluyendo la presencia de posibles defectos. Además, esta técnica se combina eficientemente con los métodos cohesivos de grietas de interfase con el fin de simular procesos de fractura en rocas. en primer lugar, la metodología actual se ha validado con los resultados del denominado Brazilian test en rocas homogéneas, con excelente precisión, y se examinan sus capacidades predictivas en situaciones incluyendo defectos únicos y múltiples. En segundo lugar, el estudio numérico comprende el análisis de rocas heterogéneas con la presencia de una interfases que separa diferentes litologías. Los resultados actuales exhiben excelentes capacidades potenciales de la presente técnica para simular fenómenos de grietas complejas en geo-materiales, lo que permite reproducir patrones de grietas tortuosas con mínima intervención del usuario. .. III 2Capítulo 1. Introducción y objetivos grietas en macizos rocosos ha sido ampliamente tratado en los últimos años mediante herramientas numéricas así como procedimientos experimentales. Esta relevancia se debe al hecho de que la presencia de estos defectos iniciales constituye el origen nuevas grietas. Figura 1.1 Macizo rocoso con presencia de grietas o defectos naturales [1] . En este sentido, con el objetivo de lograr una comprensión profunda de estos fenómenos, se ha llevado a cabo una ardua labor investigadora con el fin de poder caracterizar la propagación de grietas en las rocas. A continuación se realizará un breve repaso a algunos de los trabajos de investigación más relevantes en este contexto, desde los puntos de vista experimental y numérico en las Secciones 1.1 y 1.2, respectivamente. 1.1 Trabajos experimentales En lo concerniente a los procedimientos experimentales, el estudio en el laboratorio de los procesos de fractura en rocas ha sido un tema extensamente tratado en las últimas décadas. En este contexto en un trabajo de especial relevancia en la literatura correspondiente, Lajtai [ 3 ] llevó a cabo numerosos ensayos sobre probetas de yeso blanco (Plaster of Paris) con únicos defectos cuyas orientaciones con respecto a la carga nominal de compresión correspondía entre 0 º y 90 º con intervalos de 5 º . En estos ensayos se identificaron las diferentes etapas que puede sufrir un macizo rocoso al fracturarse bajo carga de compresión uniaxial, las cuales comprenden de forma general (Fig. 1.2): (a) fractura bajo tracción: las primeras en aparecer, las cuales tienden a tener un camino curvo, (b) fractura debido a acciones tangenciales, las cuales se propagan en dirección perpendicular a la carga principal, (c) extensión de las zonas de fractura de tracción y tangencial, y (d) aparición de zonas de fractura inclinadas. Otros ensayos de compresión uniaxial sobre piedra caliza y especímenes de grano-diorita con defectos inclinados fueron llevados a cabo por Ingraffea y Heuze [ 4 ] bajo cargas de tracción. En los resultados derivados de este estudio, se observó la secuencia de grietas primarias y secundarias durante la aplicación de la carga hasta la fractura final. Así, en la literatura relacionada, la gran mayoría de investigaciones consultadas ponen de manifiesto la dependencia del comportamiento en términos de fractura de la muestra de roca considerada con relación al número de defectos, la inclinación respectiva de cada uno de ellos así como la carga 1.2 Trabajos numéricos 3 Figura 1.2 Etapas en el proceso de fractura de un macizo rocoso con un defecto inicial identificados por Lajtai (1974) [3]. aplicada. En este sentido, y en aras de contextualizar el presente trabajo, es preciso destacar que muchos de los estudios existentes fueron llevados a cabo sobre macizos rocosos con comportamiento frágil. Sin embargo, en aplicaciones prácticas y dependiendo del tipo de roca estudiada, el alto nivel de las tensiones in-situ generalmente provocan la generación de fenómenos de plastificación de ciertas zonas alrededor de la excavación. La consideración de este aspecto está aún en desarrollo debido a la complejidad de provocar la aparición de tales fenómenos inelásticos en las condiciones de ensayo en un laboratorio. Recientemente, varios investigadores han estudiado la respuesta de la macizos rocosos por su capacidad de resistir grietas, considerando diversos aspectos tales como: ( i ) la cantidad de defectos iniciales, ( ii ) la inclinación relativa de los defectos con respecto a la carga principal y ( iii ) las diferentes cargas externas aplicadas. Estudios sobre métodos experimentales han analizado de forma pormenorizada el comportamiento de diversos macizos rocosos bajo cargas de compresión uniaxial y biaxial en rocas, a través de las cuales es posible identificar que estos materiales pueden experimentar una combinación de inicio de deslizamiento de defectos y propagación de grietas tipo wing, seguido del desarrollo de grietas secundarias (dominadas shear) que conducen a un potencial coalescencia posterior [ 5 , 6 , 7 , 8 , 9 , 10 ]. En concreto, Sagong y Bobet [ 11 ] distinguieron 9 tipos de grietas de fusión en las masas rocosas, evidenciando así la enorme variabilidad de comportamiento. Asimismo, vale la pena mencionar los estudios experimentales sobre este tema que han sido ampliamente desarrollados especialmente en las últimas dos décadas. Debido al hecho de que los métodos experimentales denominados como "estándar" presentan a menudo limitaciones relevantes en términos de introducir estados de carga complejos y caracterizar la respuesta mecánica de las muestras ante tales estados, se han venido empleando nuevas técnicas experimentales. Así, en tales circunstancias, los procedimientos experimentales tradicionales se han complementado con el uso de técnicas de correlación de imágenes digitales (DIC) y métodos de impresión 3D [ 12 ], fomentando una nueva perspectiva para la generación de muestras con los atributos deseados. 1.2 Trabajos numéricos Como se ha mencionado anteriormente, muchos de los estudios experimentales anteriores se centraron en el estudios de muestras rocosas sometidas a cargas de compresión uniaxial o biaxial, ya que las condiciones generales de carga son difíciles de reproducir en el laboratorio. Este hecho motivó el desarrollo de métodos computacionales que potencialmente permiten capturar y analizar 4Capítulo 1. Introducción y objetivos mecanismos de fractura complejos en rocas. Estos procedimientos numéricos además permiten minimizar costes económicos asociados a los dispositivos y equipos de ensayo, preparación de los testigos de rocas, entre otros, así como el proporcionar un mayor nivel de entendimiento de la problemática asociada a la propagación de grietas en macizos rocosos. La modelización de procesos de fractura en medios rocosos ha sido tradicionalmente llevado a cabo siguiendo metodologías basadas en la MFEL debido a su relativa simplicidad y robustez, y derivado del hecho de que la presión de confinamiento y la velocidad de carga en aplicaciones prácticas puede considerarse como relativamente lenta (es decir, se adoptan condiciones cuasiestática de carga). Desde un punto de vista computacional, muchos de los trabajos experimentales citados anteriormente han sido tomados como referencia para el desarrollo de técnicas numéricas para simular tales procesos de fractura. Estos procedimientos numéricos también buscan generalizar el estudio sobre el análisis de Mecánica de Rocas para diferentes solicitaciones externas además del caso de compresión simple, como son el caso de tracción, bi-tracción, bi-compresión, carga tangencial, entre otros muchos, ya que son más difíciles de llevar a cabo experimentalmente que los casos de compresión simple o el denominado ensayo triaxial. Sin embargo, los métodos basados en la MFEL presentan notables dificultades para capturar la amplia variedad de modos de fractura de dichos materiales bajo diferentes casos de carga, dando lugar a numerosos procesos de fractura que pueden potencialmente darse al mismo tiempo, así como desarrollarse grietas en localizaciones desconocidas a priori. Además, el seguimiento de tales caminos de grieta se complicará aún más por el desarrollo de ramificaciones y evoluciones de coalescencia, que son escenarios típicos que se encuentran en geo-materiales. Con el fin de solventar tales limitaciones y dificultades, en la últimas décadas se han desarrollado múltiples metodologías de cálculo computacional, especialmente en el contexto del método de los elementos finitos (MEF) y el método de los elementos de contorno (MEC). Entre las técnicas basadas en el MEC, son destacables los estudios llevados a cabo por Bobet and Einstein [ 13 ], Gonçalves da Silva and Einstein [ 14 ] para la simulación de propagación de fractura en rocas en Modos I y II en especímenes con uno o varios defectos iniciales o flaws. Posteriormente, surgieron técnicas alternativas con estos objetivos como son las denominadas: cracking particles method (CPM) [ 15 , 16 , 17 ], Peridynamics [ 18 ] y dual-horizon peridynamics [ 19 , 20 ], las cuales no requieren una representación explícita de las discontinuidades en el campo de desplazamientos debido a la fractura ni algoritmos que sean capaces de determinar el camino de la grieta durante su propagación, denominados crack tracking y que son implementados como un algoritmo adicional en el proceso de resolución numérica. Sin embargo, debido a su versatilidad y robustez, el MEF puede identificarse como el método por excelencia para la simulación de procesos de fractura en sólidos y estructuras. Dentro de las aproximaciones para simular fractura basadas en el MEF se pueden distinguir las siguientes categorías: • Técnicas basadas en grietas difusas fundamentadas en la Mecánica de Daño Continuo (MDC) a partir de la definición de una o varias variables de daño asociadas a los posibles mecanismos de fractura [ 21 ], las cuales en su versión local sufren de patologías dependientes de la malla de discretización. Es decir, la solución depende de la densidad de la malla de elementos finitos de discretización espacial. Esta limitación se debe al hecho de que una vez comienza el deterioro de las propiedades de rigidez del sólido, las ecuaciones diferenciales de gobierno del MEF pierden su carácter elíptico. Para solventar tales problemas, se han propuesto diferentes ténicas numéricas como es el uso de procedimientos de regularización integrales o metodologías no locales [ 22 , 23 , 24 , 25 , 26 ], o procedimientos alternativos de regularización según lo propuesto en [27] combinados con técnicas locales remallado local. 1.2 Trabajos numéricos 5 • Metodologías de grietas explícitas o discontinuidad fuerte. Estas técnicas se basan en un enriquecimiento bien del campo cinemático de interpolación a nivel nodal (eXtended FEM, X-FEM), fundamentados en el concepto del Partion of Unity Methods (PUM) [ 28 , 29 , 30 ], o bien formulaciones de enriquecimiento de elementos mediante modos incompatibles de deformación (enhanced FEM, E-FEM), es decir no vienen derivados del campo de desplazamientos [31, 32, 33, 34] (Fig. 1.3). • Inserción adaptativa de elementos de interfaz cohesivos durante el cálculo o implementación previa entre los bordes que delimitan la malla de elementos finitos [ 36 , 37 , 38 , 39 , 40 , 41 , 42 ]. Figura 1.3 Procedimientos numéricos de discontinuidad fuerte [35]. Aunque las anteriores estrategias se han aplicado con éxito a muchos problemas diferentes en el contexto de la Mecánica de la Fractura, todas presentan limitaciones con respecto a la predicción de inicio de grietas y en procesos de ramificación y coalescencia para frentes múltiples. Recientemente, el notable aumento de las capacidades de cálculo ha motivado el desarrollo de metodologías de simulación basadas en el MEF que retoman el concepto de daño no local y grieta difusa, considerando en su formulación una ecuación diferencial adicional referida a la denominada longitud característica del material. Especialmente enfocadas para su uso en materiales frágiles son las denominadas técnicas de phase field (PF) [ 43 , 44 ]. Este nuevo enfoque comparte algún argumento con técnicas basadas en la Mecánica de Daño Continuo (MDC) pero incorpora una formulación no local y está especialmente orientado para la modelización de fractura frágil. Los fundamentos de esta nueva metodología para fractura frágil se remontan a la formulación básica de Griffith basado en conceptos energéticos [ 45 ]. En concreto, este enfoque se desarrolla a través de la introducción de un funcional de energía total que es la suma de las contribuciones de fractura y energía elástica. La minimización de este funcional permite capturar la nucleación, propagación y coalescencia de grietas en el medio continuo de forma elegante y operativamente robusta. En este sentido, son notables las contribuciones debidas a Marigo y coatures [ 43 , 46 ], mientras que el tratamiento integral del llamado concepto Γ -convergencia ha sido abordado de forma exhaustiva en [47, 48, 49] y en las referencias citadas en ese trabajo. Las formulaciones de PF para fractura frágil han sido específicamente propuestas y rigurosamente desarrolladas por Bourdin et al. [ 46 , 50 ] y dentro de un marco termodinámicamente consistente Miehe [ 51 , 52 ], Kuhn y Müller [ 53 ], y Borden et al. [ 54 ]. Estudios recientes han ampliado esta estrategia de modelización para estructuras laminares de pequeño espesor [ 55 , 56 , 57 ], fractura dúctil [ 58 , 59 ], grietas cohesivas [ 60 ], fractura dinámica [ 61 , 62 ], y aplicaciones multi-físicas [63, 64]. Muy recientemente esta técnica ha sido empleada para la modelización de procesos de fractura en rocas con múltiples defectos iniciales [ 66 , 67 , 68 ], siendo a su vez aplicable la formulación 6Capítulo 1. Introducción y objetivos unificada propuesta en [69, 70, 71]. 1.3 Objetivos y organización 1.3.1 Objetivos Aunque durante los últimos años se han extendido de forma profusa en la literatura asociada los estudios de fractura en medios rocosos, muy pocos trabajos se han centrado en dos aspectos, los cuales serán abordados en el presente trabajo: • El uso de la técnica PF para la simulación de fractura en macizos rocosos con múltiples defectos sometidos a diferentes estados de carga. • La combinación de PF y modelos de grieta cohesiva (CZM) para la modelización de fenómenos de fractura en macizos rocosos heterogéneos compuestos de varios estratos. En este sentido se extenderá la técnica desarrollada por el autor [ 72 ], la cual introduce la interacción entre el modelo PF en el medio continuo y el modelo interfase cohesiva (que puede ser equipada con aspectos geo-mecánicos [ 73 , 74 , 75 ]), la cual puede identificarse en la separación entre las posibles litologías de macizo rocoso [76]. El uso de la metodología de simulación de procesos de fractura denominada phase field es escogida en este trabajo debido a que presenta numerosos aspectos ventajosos y de hecho su aplicación a macizos rocosos con diferentes litologías separadas por juntas/interfases que pueden inducir discontinuidades es completamente novedosa. De esta forma, la motivación del presente trabajo se fundamenta en el análisis y la comprensión del papel de la presencia de interfases en rocas, siendo ésta una temática que ha sido objeto de un gran impacto en la comunidad investigadora debido a su importancia en la práctica. Esto se debe al hecho de que las juntas en macizos rocosos poseen una cierta cantidad de cohesión debido a fenómenos de interlocking y rugosidad, conjuntamente con procesos de fricción en el caso de grandes deslizamientos, los cuales modifican la respuesta de la masa de roca [75]. Por lo tanto, el objetivo de esta investigación es abordar las capacidades predictivas de la metodología PF-CZM para su aplicación en la fractura de medios rocosos. En base a este interés, la técnica propuesta se valida primero mediante la reproducción de resultados experimentales relacionados con la fractura de roca homogénea con y sin defectos, mientras que posteriormente se extiende al análisis de medios heterogéneos. Finalmente, es preciso comentar que las simulaciones recogidas en el presente trabajo pretenden constituir un punto de partida para el desarrollo de una nueva línea de investigación en el ámbito de la simulación de procesos de fractura en geo-materiales. Este aspecto se debe al hecho de que la técnica empleada es en todo caso pionera en este campo de aplicación y que además permite la modelización de una forma versátil y sin necesidad de la intervención del usuario de fenómenos de fractura complejos, los cuales pueden tener lugar en múltiples aplicaciones prácticas dentro de la Ingeniería Civil. 1.3.2 Organización El documento se organiza de la siguiente forma. La Sección 2 resume los conceptos básicos de la técnica PF para la modelización de fractura frágil en medios continuos rocosos, mientras que las simulaciones correspondientes se presentan en la Sección 3. Posteriormente, los detalles generales de acoplamiento de grietas de interfase y de continuo siguiendo el enfoque general PF-CZM para medios heterogéneos son expuestos en la Sección 4. La Sección 5 muestra algunos resultados representativos de simulación de fenómenos de fractura en medios rocosos heterogéneos. Finalmente, las conclusiones y desarrollos futuros se discuten en la Sección 6. 2 Aspectos fundamentales de la metodología de fractura variacional de phase field En esta Capítulo se presentan de forma sucinta los fundamentos del método denominado como Phase Field (PF) para la simulación de procesos de fractura frágil, el cual es de aplicación en diversas tipologías de macizos rocosos. La formulación desarrollada en este Capítulo se fundamenta en los trabajo debidos a Miehe, Bourdin y coautores, véase [ 46 , 51 ]. Por lo tanto, diversos aspectos son omitidos en el presente documento por cuestiones de brevedad. En cualquier caso, es preciso comentar, que los fundamentos de la metodología de fractura variacional o Phase Field tiene su origen en la formulación clásica de la Mecánica de la Fractura de Griffith en términos energéticos. 2.1 Conceptos básicos En el contexto de un análisis multi-dimensional (Fig. 2.1), considérese un cuerpo arbitrario Ω∈Rnd en el espacio Euclideo de dimensión nd , cuyos puntos materiales se denotan por el vector de posición x en coordenadas Cartesianas. Las fuerzas prescritas por unidad de volumen se identifican por el vector fv:Ω−→ Rnd . La superficie que caracteriza el contorno del sólido se denota como ∂Ω∈Rnd−1 , donde ∂Ωu y ∂Ωt identifican las porciones del contorno donde están prescritas las acciones cinemáticas y estáticas, respectivamente, verificando: ∂Ωt∪∂Ωu=∂Ω y ∂Ωt∩∂Ωu=/0 . El campo de desplazamientos se denota por u , mientras que el tensor de tensiones de Cauchy se identifica como σ σ σ . Las condiciones de contorno sobre el sólido bajo estudio pueden expresarse como: u=uon ∂Ωu,and t=σ σ σ·non ∂Ωt,(2.1) donde nel vector normal exterior al sólido. En el modelo PF para grieta difusa, la grieta, la cual es habitualmente representada por una discontinuidad discreta, se regulariza mediante una variable de campo que hace las funciones de variable de daño (que cuantifica la degradación de la rigidez del sólido) d , donde d:Ω× [0,t]−→ [0,1] [ 50 ]. Esta variable de phase field cuantifica la degradación (o perdida de rigidez) del material debido a fenómenos de fractura siendo definida a nivel de punto material. De esta forma y siguiendo una evolución continua, dicha variable toma valores entre 0 y 1, donde el estado intacto se corresponde con el valor d=0yd=1identifica el estado completamente fracturado. La aproximación 1D de la función de PF se corresponde con: d=e−|x|/l , donde l indica la longitud característica del material, la cual está relacionada con la resistencia aparente del material 7 8Capítulo 2. Aspectos fundamentales de la metodología de fractura variacional de phase field 1 3 2 1 3 2  1 3  2 Figura 2.1 Comparativa entre la representación de grieta explícita (izquierda) y grieta difusa (derecha). y gobierna el ancho de la zona donde el daño es difuso. Siguiendo los desarrollos llevados a cabo en [51], las ecuaciones asociadas al problema de PF tienen la forma: d−l2∆d=0in Ωand ∇Xd·n=0in ∂Ω(2.2) donde ∆des el operador Laplaciano de d, y ∇xdes es gradiente espacial de dicha variable. En este punto, es preciso mencionar las destacable contribuciones de Wu y coautores [ 69 , 70 , 71 ], las cuales han extendido el campo de aplicación a fractura cuasi-frágil y con predicciones independientes de la longitud de característica. Con el fin de llevar a cabo la regularización difusa del campo de grieta, se define el funcional de densidad de grieta γ(d,∇xd): Γc(d):=ZΩ γ(d,∇xd)dΩ.(2.3) Como se analizó en [ 48 , 44 , 43 ], la Ec. (2.3) converge a la representación superficial de la grieta Γcen el limite l→0. Además, con el fin de preservar la consistencia termodinámica del modelo, es decir, evitar la autoreparación, se imponen la condición de evolución creciente monótona del campo de PF: ˙ Γc(d)≥0 . La forma particular del funcional de segundo orden adoptada en este trabajo tiene la expresión [ 51 ]: γ(d,∇xd) = 1 2ld2+l 2|∇xd|2.(2.4) La definición anterior de γ(d,∇xd) permite aproximar la integral de superficie original en la formulación de Griffith por la integral de volumen sobre el sólido: ZΓc Gcd∂Ω≈ZΩ Gcγ(d,∇xd)dΩ.(2.5) De esta forma, el funcional del sistema completo referido al sólido arbitrario bajo consideración puede ser dividido en las correspondientes contribuciones externas e internas, Πint(u,d) y Πext(u) , respectivamente: Π(u,d) = Πint(u,d)+Πext(u) = ZΩ ψ(ε ε ε,d)dΩ+ZΩ Gcγ(d,∇xd)dΩ+Πext(u), (2.6) donde el primer término representa la energía elástica de deformación mientras que el segundo denota la energía disipada en la grieta. Además, ψ(ε ε ε,d) es la función de energía elástica, que depende de la variable de daño d y del estado de deformaciones ε ε ε . Por otro lado, Gc es la tenacidad 2.2 Implementación numérica 9 a fractura en el sentido de Griffith y γ(d,∇xd) es el denominado funcional de densidad de grieta, que depende de la variable de daño dy su gradiente espacial ∇xd. Con respecto a la función de energía elástica almacenada en el sólido, ésta admite una descomposición espectral que permite diferenciar su parte positiva-negativa. Así, la parte positiva de la función de energía elástica depende del campo tensional asociado a las tracciones, mientras que la parte negativa a las compresiones. Según [51], tal descomposición puede expresarse: ψ(ε ε ε,d) = g(d)ψe +(ε ε ε)+ψe −(ε ε ε),(2.7a) ψe +(ε ε ε) = λ 2(htr[ε ε ε]i+)2+µtr[ε ε ε2 +],(2.7b) ψe −(ε ε ε) = λ 2(htr[ε ε ε]i−)2+µtr[ε ε ε2 −],(2.7c) donde λ y µ son las constantes de Lamé; ε ε ε+ y ε ε ε− son respectivamente las partes positivas y negativas del tensor de deformaciones. En las Ecs. (2.7) , el símbolo tr[•] representa el operador traza de un tensor, mientras h•i±es el operador de Macaulay: h•i±= (•±|•|)/2. Finalmente, la función de degradación g(d) , que introduce el deterioro o pérdida de rigidez debido al la evolución del proceso de fractura puede definirse como: g(d) = (1−d)2+K,(2.8) donde K es una valor residual de la rigidez con el fin de evitar problemas numéricos de convergencia en estados cercanos a d≈1 . Nótese además que teniendo en cuenta la Ec. (2.7a) , la función de degradación únicamente afecta a la parte positiva de la energía de deformación (Ec. (2.8) . Por lo tanto, el tensor de tensiones de Cauchy para la formulación de PF puede expresarse como σ σ σ:=∂ˆ ψ ∂ε ε ε=g(d)σ σ σ++σ σ σ−;donde σ σ σ±=λ(htr[ε ε ε]i±)1+2µε ε ε±,(2.9) siendo 1la matriz identidad de segundo orden, y σ± identifican las partes positivas y negativas del tensor de tensiones. 2.2 Implementación numérica En esta Sección se detalla la formulación del modelo de PF para la modelización de fractura frágil de rocas homogéneas en el contexto del MEF, empleando metodologías de simulación no lineal. en lo sucesivo, los operadores correspondientes al medio continuo serán identificadas con el super-índice b. Adoptando una aproximación isoparamétrica, las funciones de interpolación Lagrangianas a nivel de elemento, NI(ξ ξ ξ) , son definidas en el espacio paramétrico ξ ξ ξ={ξ1,ξ2} . Dicho espacio se emplea para la la interpolación de la geometría ( x ), el campo de desplazamientos, ( u ), su variación (δu) y su linealización correspondiente (∆u): x∼ = n ∑ I=1 NIexI=Nex;u∼ = n ∑ I=1 NIdI=Nd;δu∼ = n ∑ I=1 NIδdI=Nδd;∆u∼ = n ∑ I=1 NI∆dI=Nδd, (2.10) donde n indica el número de nodos a nivel elemental; xI y dI identifican respectivamente las coordenadas y desplazamientos nodales (es decir sus vectores correspondientes), los cuales son agrupados de forma matricial en los operadores ex y d . De forma similar, las funciones de forma de interpolación son agrupadas en la matriz N. El proceso de interpolación del campo de deformaciones ( ε ε ε ), su variación ( δε ε ε ) y su linealización ( ∆ε ε ε ) se lleva a cabo a partir de un operador de compatibilidad en su forma discreta, Bd , mediante el 10 Capítulo 2. Aspectos fundamentales de la metodología de fractura variacional de phase field campo de desplazamientos: ε ε ε∼ =Bdd;δε ε ε∼ =Bdδd;∆ε ε ε∼ =Bd∆d.(2.11) En relación a la interpolación de la variable de phase field a nivel de elemento, se emplean las mismas funciones de aproximación en el espacio paramétrico ξ ξ ξ={ξ1,ξ2} . Así, puede aproximarse a nivel de punto de cuadratura tanto la variable de phase field como su correspondiente variación y linealización: d∼ = n ∑ I=1 NIdI=Nd;δd∼ = n ∑ I=1 NIδdI=Nδd;∆d∼ = n ∑ I=1 NI∆dI=Nδd,(2.12) donde dI son los valores nodales de la variable de phase field, que se agrupa en el vector d . El gradiente espacial de dicha variable ( ∇xd ), su variación ( ∇xδd ) y su linealización ( ∇x∆d ) se interpola mediante el operador Bd: ∇xd∼ =Bdd;∇x(δd)∼ =Bdδd;∇x(∆d)∼ =Bd∆d.(2.13) Teniendo en cuenta los desarrollos anteriores, la discretización de la Ec. (2.6) a nivel elemental (siendo identificado con el super-índice el) puede expresarse: δ˜ Πel b(d,δd,d,δd) =δdTZΩel h(1−d)2+KBT dσ σ σ++BT dσ σ σ−idΩ−Z∂Ωel NTtd∂Ω−ZΩel NTfvdΩ +δdTZΩel −2(1−d)NTψe +(ε ε ε)dΩ+ZΩel Gb clBT d∇xd+1 l2NTddΩ =δdTfb d+δdTfb d(2.14) donde el vector residual fb d=fb d,ext −fb d,int , en la cual , fb d,int y fb d,ext representan las contribuciones internas y externas, respectivamente, vienen dadas por: fb d,int =ZΩel h(1−d)2+KBT dσ σ σ++BT dσ σ σ−idΩ,(2.15) fb d,ext =Z∂Ωel NTtd∂Ω+ZΩ NTfvdΩ.(2.16) Un proceso equivalente de discretización puede realizarse para la ecuación residual correspondiente al campo de phase field. De esta forma, el vector residual correspondiente, fb d , tiene la siguiente expresión: fb d=ZΩel −2(1−d)NTψe +(ε ε ε)dΩ+ZΩel Gb clBT d∇x+d1 l2NTddΩ.(2.17) En este trabajo, se ha empleado un proceso de resolución numérica incremental-iterativa del conjunto acoplado de ecuaciones no lineales siguiendo un esquema Newton-Raphson según el esquema presentado en [ 77 ]. De esta forma, una vez efectuada la correspondiente linealización, se obtiene siguiente sistema de ecuaciones acopladas [59]: Kb dd Kb dd Kb ddKb dd∆d ∆d=fb d,ext 0−fb d,int fb d.(2.18) La forma particular de los operadores Kb dd , Kb dd , Kb dd y Kb dd son omitidas por brevedad, véase [ 78 ]. La implementación numérica se ha llevado a cabo en los softwares de elementos finitos FEAP y ABAQUS, mediante subroutinas de usuario. 3 Aplicaciones representativas de simulaciones numéricas de fractura frágil en rocas homogéneas En este Capitulo se presentan una serie de ejemplos de aplicación en la que se modeliza procesos de fractura frágil en rocas homogéneas. 3.1 Validación de la metodología de modelización de fractura frágil phase field: aplicación al Brazilian test La metodología presentada en el Capitulo anterior, denominada como "Phase Field" (PF), es validada mediante su aplicación al denominado Brazilian test. Este ensayo consiste consiste en la aplicación de una carga de compresión uniaxial a una muestra de roca. En concreto, el proceso de validación presentado en esta Sección se lleva a cabo a través de la directa correlación numérica experimental con los datos del ensayo reportados en [66] para especímenes homogéneos. La definición geométrica de la muestra de roca ensayada consiste en una probeta con sección transversal circular cuyo diámetro es igual a 52 mm, ver Fig.3.1.a. Las propiedades materiales son: módulo de Young E = 31.5 GPa, coeficiente de Poisson ν = 0.25 y tenacidad a fractura Gb c = 100 J/m 2 . Además se estima el valor del parámetro de escala de phase field, el cual puede ser relacionado con la tensión de rotura σcdel material a partir de la siguiente expresión: l=27EGb c 256σ2 c .(3.1) De esta forma, asignamos l = 1 mm. Con el fin de simplificar el análisis en términos computacionales, se hace uso de la simetría del ensayo en términos materiales, geométricos y carga. Por lo tanto, unicamente se modeliza la mitad de espécimen, siendo discretizado con 55452 elementos cuadriláteros de primer orden. Las simulaciones fueron llevadas a cabo imponiendo un desplazamiento vertical que implica la compresión del componente. La Fig.3.1.b muestra el patrón de fractura obtenido obtenido mediante la implementación de la metodología PF en el software ABAQUS [ 66 ]. Durante la simulación se observó que la grieta se originó inicialmente en el centro de la probeta, siendo un resultado en consonancia con los resultados experimentales. En etapas posteriores de carga, esta grieta inicial se propagó hacia el eje superior del dominio bajo estudio. Desde un punto de vista cualitativo, la Fig. 3.2 muestra la correlación numérica-experimental del diagrama fuerza-desplazamiento donde se observa un excelente acuerdo con los resultados obtenidos por [ 66 ]. La respuesta de la probeta se caracteriza por una primera etapa lineal hasta que 11 18 Capítulo 4. Acoplamiento entre la metodología de fractura frágil en rocas phase field y grietas de interfase según el modelo de fractura cohesiva: aplicación a medios heterogéneos rocosos 1 3 2 1 3 2 Figura 4.1 (a) Coexistencia entre grieta en el medio continuo y de interfase en el contexto de la metodología PF-CZM. Uno de los aspectos fundamentales de la metodología de simulación PF-CZM consiste en asumir una descomposición aditiva de la parte dispositiva del funcional anterior. De esta forma, la energía de fractura del sistema Gc puede descomponerse en dos contribuciones asociadas con: ( i ) la fractura disipada en el medio continuo, en este caso en cada estrato del macizo rocoso,estando caracterizada por Gb c y estando tales procesos simulados con la metodología PF, y ( ii ) la energía de fractura disipada en las potenciales interfases, caracterizada por la energía de fractura Gi , y que se modeliza con la técnica de grieta cohesiva (CZM). Por lo tanto, la Ec. (4.1) puede reescribirse de forma: Π(u,Γb,Γi) = ΠΩ+ΠΓb+ΠΓi=ZΩ\Γ ψe(ε ε ε)dΩ+ZΓb Gb c(u,d)dΓ+ZΓi Gi(g,h,d)dΓ,(4.2) donde g identifica la diferencia de desplazamientos entre dos puntos materiales entre los flancos de la interfase, h es una variable de estado para evitar que la interfase recupere su rigidez una vez fracturada [60], y des la variable de phase field en el continuo. En este trabajo, por simplicidad, asumimos un comportamiento de interfase tipo tensión cut-off [ 72 ], cuya rigidez aparente puede acoplarse con la variable de phase field en el continuo vía d . Es más, nótese que la energía disipada en la interfase, Ec. (4.2) , puede descomponerse a su vez en la suma de las contribuciones asociadas a las contribuciones de los Modos ce Fractura I y II en un contexto de modelización 2D, estando representados por las energías correspondientes Gi IC yGi IIC. En la formulación derivada en [ 72 ], se asume que la apertura crítica de la interfase ( gc ) sigue una ley de dependencia lineal con respecto al estado de daño en el continuo adyacente a la interfase. dicha dependencia puede postularse: gc(d) = (1−d)gc,0+dgc,1 , donde gc,0=gc(d=0) y gc,1=gc(d=1) . Así, la ley tensión-desplazamientos en la interfase (véase Fig. 4.2 para el Modo I de Fractura), se expresa por: σ=     kn gn gnc ,if 0<gn gnc <1; 0,if gn gnc ⩾1,(4.3) τ=     kt gt gtc ,if 0<gt gtc <1; 0,if gt gtc ⩾1,(4.4) donde σ y τ son las tensiones asociadas a los Modos I y II, respectivamente, g son los desplazamientos relativos en la interfases entre los flancos de la interfase, haciendo referencia los subíndices nythacen referencia a los Modos I y II, respectivamente. 4.2 Formulación variacional y discretización mediante elementos finitos 19 Como resultado, la rigidez aparente de la interfase k de ambos Modos I y II depende del daño en el continuo den el continuo adyacente: kn=kn,0gnc,0 gnc 2 ;kt=kt,0gtc,0 gtc 2 ,(4.5) donde k0 y g0 identifican la rigidez y el campo de desplazamiento relativo en la interfase para el estado intacto, i.e. d=0. Adicionalmente, se adopta un criterio cuadrático para identificar el fallo de la interfase una vez que se disipa la energía de fractura correspondiente: Gi I Gi IC 2 +Gi II Gi IIC 2 =1,(4.6) donde Gi I y Gi II con las tasas de liberación de energía en la interfase, las cuales adquieren la forma: Gi I(d) = 1 2nt,0g2 n g2 nc,0 (1−d)gnc,0+dgnc,12;Gi II(d) = 1 2kt,0g2 t g2 tc,0 (1−d)gtc,0+dgtc,12.(4.7) Tales tasas de liberación de energía son comparadas con sus valores críticos Ec. (4.6) , donde Gi IC yGi IIC los cuales son parámetros físicos y que son independientes del daño en el continuo: Gi IC =1 2g2 nc,0kn,0;Gi IIC =1 2g2 tc,0kt,0.(4.8) Finalmente, en el presente trabajo se ha extendido el comportamiento de la interfase para tener en consideración aspectos geo-mecánicos. En concreto, se ha seguido la aproximación llevada a cabo en [ 79 ], en el la rigidez de la interfase bajo compresión está relacionada con los parámetros geo-mecánicos de la misma: gn≤0→σ=kn,0gn 1−gn/gm ,(4.9) donde gn es el cierre de la interfase, kn,0 es la rigidez inicial, gm el es cierre máximo de la interfase. Los valores kn,0 y gm , pueden relacionarse con parámetros (JCS), que puede encontrarse en la literatura asociada para diversas naturalezas rocosas (Fig. 4.2): kn,0=-7.15+1.75 JRC+0.02JCS a0-0.2510 (4.10) gm=-7.15+1.75 JRC+0.02JCS a0-0.2510 (4.11) donde JCS es la rigidez a compresión de la junta, véase Barton and Choubey [ 73 ]y[ 74 ], donde a0=JRC0/50 . Estos coeficientes pueden obtenerse de los ensayos realizados en el laboratorio para una muestra de roca. 4.2 Formulación variacional y discretización mediante elementos finitos En esta sección se deriva la formulación variacional de la contribución de la interfase en el potencial total del sistema, Ec. (4.2) . Siguiendo un procedimiento de Galerkin en el contexto del MEF, la 20 Capítulo 4. Acoplamiento entre la metodología de fractura frágil en rocas phase field y grietas de interfase según el modelo de fractura cohesiva: aplicación a medios heterogéneos rocosos Figura 4.2 Modificación de la rigidez de la interfase a compresión según parámetros geo-mecánicos gn. forma débil de la contribución de la interfase al sistema puede expresarse como: δΠΓi(u,δu,d,δd) = ZΓi∂Gi(u,d) ∂uδu+∂Gi(u,d) ∂dδddΓ,(4.12) δΠΓi(u,δu,d,δd) = ZΓi∂Gi(u,d) ∂uδu+∂Gi(u,d) ∂dδddΓ,(4.13) donde δu es el vector de funciones de aproximación del campo de desplazamientos cinemáticamente admisible ( Vu=δu|u=uon ∂Ωu,u∈H1 ), y δd corresponde a las funciones de aproximación de phase field (Vd=δd|δd=0on Γb,d∈H0). La discretización del presente modelo de interfase se realiza en el contexto del MEF empleando elementos isoparamétricos de primer orden con interpolación lineal del campo de desplazamientos, aunque su extensión a elementos de órdenes superiores puede llevarse a cabo de forma directa. Análogamente al modelo de PF en el continuo, (Sección 2.2), d denota el vector de desplazamientos nodales y ¯ d el vector de phase field nodales a nivel de elemento. Por lo tanto, la forma discreta de la Ec.(4.12) para cada elemento de interfase Γel i(Γi∼SΓel i) adquiere la forma: δ˜ Πel Γi(d,δd,¯ d,δ¯ d) = ZΓel i∂Gi(d,¯ d) ∂dδd+∂Gi(d,¯ d) ∂¯ dδ¯ ddΓ,(4.14) donde las expresiones de Gi=Gi I+Gi II, se han detallado anteriormente. El vector de desplazamientos relativos entre los flancos de la interfase g para cualquier punto de Γel i es el resultado de la diferencia entre los desplazamientos de puntos opuestos en la interfase, los cuales pueden ser obtenidos a partir de los desplazamientos nodales d multiplicados por el operador L cuya definición puede encontrase en [ 40 ] (y que se omite en el presente documento por cuestiones 4.2 Formulación variacional y discretización mediante elementos finitos 21 de brevedad): g=NLd =ˆ Bdd,(4.15) donde N es una matriz que incluye las funciones Lagrangianas de interpolación de los desplazamientos y ˆ Bd=NL identifica el operador de compatibilidad en la interfase. Para aplicar la relación de constitutiva de la interfase invocando el concepto de modelo cohesivo (CZM), siendo tal ley expresada en el sistema de referencia local definido por los vectores unitarios normal y tangencial a la interfase [ 39 , 42 ], el vector de desplazamientos relativos en el sistema global en Ec. (4.15) se multiplica por el tensor de rotación Rde forma gloc: gloc ∼ =Rg =Rˆ Bbd.(4.16) De forma similar, las siguientes expresiones y operadores son introducidos para el cálculo del valor de la variable de phase field promedio en la interfase d,Γel i, a nivel elemental: d∼ =NdMd¯ d=ˆ Bd¯ d,(4.17) donde Md es un operador para calcular la media de la variable de phase field a través de la interfase entre los dos flancos y ˆ Bd=NdMd identifica el operador de compatibilidad. La forma particular de estos operadores para formulaciones 2D y 3D es derivada en [40, 42]. Por lo tanto, la forma débil discreta correspondiente a la interfase adquiere la forma: δ˜ Πel Γi(d,δd,¯ d,δ¯ d) = δdTZΓel i∂Gi(d,¯ d) ∂dT dΓ+δ¯ dTZΓel i∂Gi(d,¯ d) ∂¯ dT dΓ =δdTZΓel i ˆ BT dRT∂Gi(d,¯ d) ∂gloc T dΓ+δ¯ dTZΓel i ˆ BT d∂Gi(d,¯ d) ∂¯ dT dΓ (4.18) dando lugar a los vectores residuales asociados tanto al campo de desplazamientos como de phase field: fi d=ZΓel i ˆ BT dRT∂Gi(d,¯ d) ∂gloc T dΓ,(4.19a) fi d=ZΓel i ˆ BT d∂Gi(d,¯ d) ∂dT dΓ.(4.19b) Finalmente, se calcula la linealización consistente de los vectores residuales a partir del concepto de derivada direccional, obteniéndose los siguientes operadores: Ki dd =∂fd ∂d=ZΓel i ˆ BT dRTCi ddRˆ BddΓ,(4.20a) Ki dd=∂fd ∂d=ZΓel i ˆ BT dRTCi ddˆ BddΓ,(4.20b) Ki dd=∂fd ∂d=ZΓel i ˆ BT dCi ddRˆ BddΓ,(4.20c) Ki dd =∂fd ∂d=ZΓel i ˆ BT dCi dd ˆ BddΓ,(4.20d) 22 Capítulo 4. Acoplamiento entre la metodología de fractura frágil en rocas phase field y grietas de interfase según el modelo de fractura cohesiva: aplicación a medios heterogéneos rocosos donde las matrices tangentes del la interfase toman la forma: Ci dd =ˆ αkn0 0ˆ βkt,(4.21a) Ci dd="gnkn ∂ˆ α ∂d,gtkt ∂ˆ β ∂d#,(4.21b) Ci dd=   gnkn ∂ˆ α ∂d gtkt ∂ˆ β ∂d   ,(4.21c) Ci dd =1 2g2 nkn ∂2ˆ α ∂d2+1 2g2 tkt ∂2ˆ β ∂d2.(4.21d) En las expresiones anteriores ˆ αyˆ βson particularizados: ˆ α=g2 nc,0 (1−d)gnc,0+dgnc,12,(4.22a) ˆ β=g2 tc,0 (1−d)gtc,0+dgtc,12.(4.22b) De forma análoga a Ec. (2.18) , el sistema de ecuaciones acoplados que se resuelve en el contexto del MEF no lineal mediante un procedimiento de solución incremental-iterativa siguiendo un esquema de Newton-Raphson multi-dimensional. El sistema resultante puede por tanto expresarse como: Ki dd Ki dd Ki ddKi dd∆d ∆d=fi d fi d.(4.23) 5 Aplicaciones representativas de simulaciones numéricas de fractura frágil en rocas heterogéneas empleando la metodología PF-CZM 5.1 Aspectos Generales La presencia de interfases en propagación de grietas en macizos rocosos puede considerarse especialmente significativa en el caso de rocas sedimentarias compuestas de capas de lutita y lutitas. El análisis de eventos de fractura en materiales y estructuras con múltiples capas ha sido objeto de investigación recurrente en los últimos años, como se describe en [72]. De forma similar a lo que ocurre en el caso de los materiales compuestos, hay tres aspectos fundamentales de crucial importancia en la fracturación a través de interfaces, llamadas también juntas en la literatura de Mecánica de Rocas: ( i ) la resistencia de la junta, ( ii ) las propiedades mecánicas de las capas y (iii) las condiciones de carga. Como se discutió en [ 80 ], cuando una grieta que se aproxima a contra una interfaz resistente, esta grieta generalmente tiende a propagarse en la sólido adyacente sin experimentar desviaciones significativas de su trayectoria inicial. Por el contrario, en el caso de interfases débiles, éstas son mucho más propensas a fallar por lo tanto conllevando fenómenos de desviación de la grieta inicial (delaminación). Este escenario básico puede reproducirse de manera directa mediante el enfoque PF-CZM descrito en la sección anterior. En este sentido, la Fig. 5.1 muestra de forma cualitativa un sistema de longitud L = 50 mm con las propiedades mecánicas descritas en los ejemplos tratado anteriormente. Dicho sistema está sometido a la acción de una carga uniformemente distribuida, la cual se introduce en el sistema mediante la imposición de los desplazamientos laterales ¯ux . En primer lugar, asignamos el ratio entre las tenacidades de fractura del medio continuo y de la interfase como Gb c/Gi c = 50 como caso de interfase débil, mientras que Gb c/Gi c = 0.05 para el caso de interfase tenaz. En este gráfico, puede observarse que la metodología desarrollada PF-CZM reproduce de forma fidedigna el comportamiento básico descrito en [ 80 ]. En concreto, se observa como en el caso de interfase débil se obtiene claramente una discontinuidad del campo de desplazamientos verticales, evidenciando un proceso de delaminación. Por el contrario, en el caso de interfase tenaz se predice un patrón de fractura perpendicular a la carga impuesta sin observarse ningún proceso de deflexión de la grieta inicial. Sin embargo, para sistemas bimateriales con interfases débiles, se pueden encontrar situaciones mucho más complejas según los estudios basados en la MFEL, véase el estudio llevado a cabo en 23 24 Capítulo 5. Aplicaciones representativas de simulaciones numéricas de fractura frágil en rocas heterogéneas empleando la metodología PF-CZM Figura 5.1 Sistema sometido a desplazamiento uniaxial. Predicciones utilizando PF-CZM para dos tipos de interfaces entre materiales similares: interfases/juntas débiles con eventos de delaminación a lo largo de la interfaz e interfases/juntas tenaces con trayectoria de grieta continua. [ 72 ]. En particular, se pueden encontrar tres escenarios de propagación diferentes basados en los parámetros de Dundurs del sistema con una grieta incidente en una interfase: (i) desviación simple a lo largo de la interfase, (ii) desviación doble a lo largo de la interfase y (iii) penetración. Figura 5.2 Geometría y condiciones de contorno para el problema bimaterial. Escenarios de patrones de grietas para una interfase frágil. Las gráficas muestran el mapa del campo de desplazamiento vertical correspondiente a tres casos diferentes: A (doble deflexión), B (deflexión simple), C (penetración). 5.2 Aplicación a rocas salinas con interfases inclinadas 25 Basándonos en las capacidades potenciales que muestra la metodología PF-CZM, las siguientes Secciones abordan la aplicabilidad del dicho método para simular eventos de fractura en masas de rocas compuestas de varios estratos. 5.2 Aplicación a rocas salinas con interfases inclinadas Una de las aplicaciones potenciales del marco numérico desarrollado en este trabajo es la investigación de la respuesta de las rocas salinas con interfases. Estas rocas muestran un comportamiento mecánico diferente de la litología correspondiente, que puede obedecer a configuraciones con una o varias capas. En consonancia con el estudio llevado a cabo en [ 81 ], los litotipos mixtos presentan lechos horizontales o inclinados de anhidrita que se alternan con halita, dando lugar a rocas compuestas de anhidrita-halita. Analizando de las evidencias experimentales de [ 81 ], se puso de manifiesto que en tales configuraciones la presecia de interlayers tienden a promover en inicio de grietas en tales localizaciones. Estas interlayers son identificadas con láminas con propiedades mecánicas menos rígidas, resistentes y tenaces que los estratos adyacentes. Este hecho además viene condicionado por la presencia de interfases con geometrías inclinadas. Estos autores investigaron el comportamiento bajo carga uniaxial y triaxial de compresión del proceso de fractura en rocas salinas con presencia de interlayers,ver Fig. 5.3, donde se aprecian la preparación de los testigos de rocas analizados así como la presencia de la interlayer (donde se incluyen diferentes espesores). Figura 5.3 Muestra rocosa salina con interlayer: muestras y definición de interlayer. La respuesta bajo carga uniaxial de compresión de estos especímenes se recoge en la Fig. 5.4 con sus correspondientes patrones de fractura. Más concretamente, estos autores discutieron que los patrones de fractura típicos en estas rocas siguen diferentes fases: ( i ) cerca/a lo largo de la región de la interfase, los micro-defectos iniciales derivados de las diferencias entre sedimentación y consolidación evolucionan generando microgrietas; ( ii ) aumentando la carga aplicada, tales micro-grietas se propagan a lo largo de la interfase 26 Capítulo 5. Aplicaciones representativas de simulaciones numéricas de fractura frágil en rocas heterogéneas empleando la metodología PF-CZM Figura 5.4 Respuesta rocosa salina con interlayer. y dentro de la región de roca salina; ( iii ) finalmente, se produce la coalescencia de diferentes grietas, rompiendo la muestra. Esta secuencia se aprecia en la Fig. 5.5. Figura 5.5 Representación de procesos de fractura de muestras de rocosa salina con interlayer. La metodología numérica desarrollada en este trabajo se aplica para estudiar el comportamiento de la deformación y las características de fractura (iniciación y coalescencia) en rocas salinas con interfaces inclinadas bajo compresión uniaxial. Esta configuración fue analizada en [ 81 ]. Para este propósito, se analizan numéricamente escenarios de compresión uniaxiales de rocas compuestas 5.2 Aplicación a rocas salinas con interfases inclinadas 27 de sal con capas intermedias de 20 mm de espesor, que mostraron un comportamiento mecánico cercano al carácter frágil. Debido a la falta de datos con respecto a la geometría de los especímenes analizados en [ 81 ], el presente estudio se restringe meramente a un análisis cualitativo del papel del interlayer central en el proceso de fractura, ver Fig. 5.6. De esta forma, se define un dominio de 100 × 200 mm con una capa intermedia central de 20 mm de espesor siguiendo una representación macroscópica con dos pequeños defectos cerca de la interfase. Este dominio se discretiza empleando 84388 elementos que combinan que combinan la metodología PF y CZM. Las propiedades mecánicas y de fractura de ambos materiales se enumeran en la Tabla 5.1. Es a su vez preciso mencionar que para la interfase dispuesta entre roca salina y la anhidrita, elegimos un valor intermedio entre las propiedades de grietas del límite transgranular reportadas en [82], de modo que se establece Gi c=1.2715 J/m2. Figura 5.6 Muestra rocosa compuesta de roca salina bajo carga de compresión uniaxial. Geometría y condiciones de contorno, patrón de fractura experimental y predicciones numéricas. Tabla 5.1 Propiedades mecánicas y de fractura: halita y anhidrita. Material E (GPa) νGb c(Jm−2) Halita 36.87 0.254 1.155 Anhidrita 74.36 0.269 1.805 La Fig. 5.6 representa la correlación cualitativa numérica-experimental del patrón de grieta predicha por la técnica numérica de predicción PF-CZM y la observada en los experimentos. como puede desprenderse de este gráfico, se observa un ajuste cualitativo muy satisfactorio. En concreto, durante las simulaciones computacionales, se observa que las grietas verticales principales se inician en los defectos cerca de la interfase y se propagan a través de la capa intermedia; posteriormente, ambas grietas se propagan a lo largo interfase originándose un fenómeno de coalescencia. El patrón de fractura final en esta localización coincide desde un punto de vista cualitativo con el reportado en [ 81 ]. Sin embargo, tenga en cuenta que el desarrollo de los caminos de grietas tipo wing yshear se desvían ligeramente de las evidencias experimentales, por lo que sería conveniente acudir a formulaciones alternativas de phase field tales como las desarrolladas en [67, 68]. 34 Capítulo 5. Aplicaciones representativas de simulaciones numéricas de fractura frágil en rocas heterogéneas empleando la metodología PF-CZM Figura 5.14 Espécimen con tres defectos inclinados e interfase vertical. (a) Condiciones de carga y geometría. (b) Mapa de daño para ¯uy=−0.1 ¯ux. (c) Mapa de daño para ¯uy=−¯ux. . Figura 5.15 Espécimen con tres defectos inclinados e interfase vertical. (a) Curva cargadesplazamiento para ¯uy=−0.1 ¯ux . (b) Curva carga-desplazamiento para ¯uy=−¯ux . . 6 Conclusiones y desarrollos futuros 6.1 Conclusiones En el presente trabajo se ha desarrollado una metodología de simulación de procesos de fractura en macizos rocosos basadas en el método de PF en combinación con el método de interfase cohesiva (CZM). La presenta técnica se ha implementado en el contexto del Método de los Elementos Finitos (MEF) no lineal bajo deformaciones infinitesimales. Para tal efecto, se han programado dos elementos finitos de desarrollo propio en los softwares FEAP y ABAQUS. La generación de la malla se ha realizado empleando ABAQUS, por lo que se han generado una serie de scripts en Matlab para poder tener la correspondencia de datos entre ambos softwares. Adicionalmente, se ha complementado el trabajo mediante la modificación del código previamente existente para incluir en el comportamiento a compresión normal del macizo rocoso parámetros geo-mecánicos. La presente metodología se ha validado en primera instancia para la modelización de procesos de fractura en medios rocosos homogéneos mediante el denominado Brazilian test, obteniéndose resultados muy satisfactorios. Posteriormente, el método desarrollado se ha aplicado a la simulación de eventos de fractura complejos en masas de roca heterogéneas con juntas/interfases débiles y múltiples defectos. Así pues, a tenor de los resultados obtenidos, se ha puesto de manifiesto que a partir de la presente técnica de modelización es posible caracterizar fenómenos de coalescencia, y bifurcaciones de grietas. A su vez, en combinación con el modelo de interfase, se ha comprobado como es posible modelizar la concomitancia de fenómenos de grietas de diversa índole, suponiendo un avance significativo en el estado actual de la técnica. 6.2 Desarrollos futuros Como potenciales desarrollos futuros, el estudio bibliográfico realizado en este trabajo ha puesto de manifiesto la importancia de desarrollar herramientas numéricas capaces de predecir fenómenos de fractura en macizos rocosos con defectos iniciales. Como planes de futuro inmediato, se realizará la correlación numérico-experimental de los resultados de simulaciones empleando la presente técnica con resultados existentes en la literatura. Seguidamente se extenderá el análisis para diferentes tipologías de interfase. Como plan adicional, se tiene en mente la producción de probetas de PMMA u otros materiales con el fin de poder realizar ensayos experimentales en los Lab de la Universidad de Sevilla. Tales probetas incluirán defectos e interfases y se procederá a su monitorización para capturar su respuesta mecánica, así como el uso de técnicas de video-correlación para caracterizar el camino de las grietas 35 36 Capítulo 6. Conclusiones y desarrollos futuros para diferentes estados de carga. De esta forma dichas probetas serán a su vez modeladas con la presente técnica para evaluar su bondad y capacidad predictiva. Finalmente desde un punto de vista más asociado tanto a los desarrollos teórico-prácticos, direcciones futuras considerarán la extensión del método propuesto para rocas anisotrópicas, fractura hidráulica, medios porosos y formulaciones modificadas rotura bajo cargas de cortadura. En esta dirección, se ha elaborado una propuesta de proyecto presentada al Plan Nacional de Investigación para extender el alcance de esta técnica a suelos parcialmente saturados. El objetivo fundamental es la caracterización de los procesos de fractura frágil en margas y en arcillas sobreconsolidadas, en condiciones parcialmente saturadas, con gran interés en la aplicación de resultados para la correcta previsión de fallos frágiles de estos materiales en estructuras y construcciones civiles. Estas acciones conllevarán el desarrollo de un marco de análisis computacional de procesos de fractura hidro-mecánica en el contexto de PF. Así, se llevará a cabo la formulación de un nuevo marco teórico-numérico aplicable a los suelos bajo estudio en este proyecto así como su correspondiente formulación numérica empleando el MEF y su posterior validación con ensayos experimentales. Lista de Figuras 1.1 Macizo rocoso con presencia de grietas o defectos naturales [1] 2 1.2 Etapas en el proceso de fractura de un macizo rocoso con un defecto inicial identificados por Lajtai (1974) [3] 3 1.3 Procedimientos numéricos de discontinuidad fuerte [35] 5 2.1 Comparativa entre la representación de grieta explícita (izquierda) y grieta difusa (derecha) 8 3.1 (a) Condiciones de carga: Brazilian splitting test. (b) Patrón de fractura. 12 3.2 Brazilian splitting test: diagrama fuerza-desplazamiento. 12 3.3 Testigo de roca con defectos inclinados: geometría y cargas 13 3.4 Caso de compresión uniaxial: mapa de daño inicial (a) y en una etapa posterior (b) 14 3.5 Diferentes estados de carga para testigo de roca con defectos inclinados: patrones de fractura. (a) Caso de compresión uniaxial. (b) Compresión vertical y tracción horizontal con ¯uy/¯ux= -1. (c) Compresión vertical y tracción horizontal con ¯uy/¯ux= -2 15 3.6 Diferentes estados de carga para testigo de roca con defectos inclinados: patrones de fractura. (a) Tracción vertical y tracción horizontal con ¯uy/¯ux= 1. (b) Tracción vertical y tracción horizontal con ¯uy/¯ux= 2. (c) Tracción vertical y tracción horizontal con ¯uy/¯ux= 10 15 4.1 (a) Coexistencia entre grieta en el medio continuo y de interfase en el contexto de la metodología PF-CZM 18 4.2 Modificación de la rigidez de la interfase a compresión según parámetros geo-mecánicos gn20 5.1 Sistema sometido a desplazamiento uniaxial. Predicciones utilizando PF-CZM para dos tipos de interfaces entre materiales similares: interfases/juntas débiles con eventos de delaminación a lo largo de la interfaz e interfases/juntas tenaces con trayectoria de grieta continua 24 5.2 Geometría y condiciones de contorno para el problema bimaterial. Escenarios de patrones de grietas para una interfase frágil. Las gráficas muestran el mapa del campo de desplazamiento vertical correspondiente a tres casos diferentes: A (doble deflexión), B (deflexión simple), C (penetración) 24 5.3 Muestra rocosa salina con interlayer: muestras y definición de interlayer 25 5.4 Respuesta rocosa salina con interlayer 26 5.5 Representación de procesos de fractura de muestras de rocosa salina con interlayer 26 5.6 Muestra rocosa compuesta de roca salina bajo carga de compresión uniaxial. Geometría y condiciones de contorno, patrón de fractura experimental y predicciones numéricas 27 37 38 Lista de Figuras 5.7 Propiedades de fractura intra- e intergranular en materiales rocosos [82]. Tabla 5 en tal publicación 28 5.8 Propiedades de fractura intra- e intergranular en materiales rocosos [82]. Tabla 6 en tal publicación 29 5.9 (a) Espécimen con dos defectos inclinados e interfase horizontal. Caso ¯uy=−0.2 ¯ux: (b) campo de displacementos horizontales; (c) variable de daño. 30 5.10 Espécimen con dos defectos inclinados e interfase horizontal. Caso ¯uy=0.5 ¯uxtraccióntracción: (a) desplazamientos verticales al final de la simulación; (b) mapa de fractura al final de la simulación; (c) Evolución fuerza-desplazamiento con incluyendo patrones de fractura en diferentes fases 31 5.11 Espécimen con dos defectos inclinados e interfase horizontal. Caso ¯uy=¯uxtraccióntracción: (a) desplazamientos verticales al final de la simulación; (b) mapa de fractura al final de la simulación; (c) Evolución fuerza-desplazamiento incluyendo patrones de fractura en diferentes fases 32 5.12 (a) Espécimen con tres defectos inclinados e interfase horizontal. (b) Patrón de fractura (mapa de la variable phase field) al final de la simulación 32 5.13 Espécimen con tres defectos inclinados e interfase horizontal. Curva carga-desplazamiento incluyendo patrones de fractura en diferentes fases y mapa de desplazamientos horizontales 33 5.14 Espécimen con tres defectos inclinados e interfase vertical. (a) Condiciones de carga y geometría. (b) Mapa de daño para ¯uy=−0.1 ¯ux. (c) Mapa de daño para ¯uy=−¯ux. 34 5.15 Espécimen con tres defectos inclinados e interfase vertical. (a) Curva carga-desplazamiento para ¯uy=−0.1 ¯ux. (b) Curva carga-desplazamiento para ¯uy=−¯ux. 34 Lista de Tablas 5.1 Propiedades mecánicas y de fractura: halita y anhidrita 27 5.2 Propiedades de fractura de la interfase 29 5.3 Propiedades geo-mecánicas de la interfase 29 39 Bibliografía [1] Hoek, E. and Bieniawski. (1984) Brittle Fracture Propagation in Rock under Compression. International Journal of Fracture 26. Pp. 276-294. [2] Hoek, E. Brittle Failure of Rock. Imperial College Rock Mechanics Research Report No 2. 1967. [3] Lajtai, E.Z. (!974) Brittle Fracture in Compression. International Journal of Fracture, Vol 10 N4, pp. 525-536. [4] Ingraffea, A.R., Heuze, F.E. (1977) An Analysis of Discrete Fracture Propagation in Rock Loaded in Compression. Proceedings of the 1 8 th U.S. Symposium on Rock Mechanics, Keystone, Colorado, pp. 2A4-1 to 2A4-7, 1977. [5] L.D. Suits, T.C. Sheahan, L.N.Y. Wong, H.H. Einstein (2009) Using high speed video imaging in the study of cracking processes in rock, Geotech. Test. J. 32 (2) 164–180. [6] A. Bobet, Numerical simulation of initiation of tensile and shear cracks, in: DC Rocks (2001) the 38th US Symposium on Rock Mechanics (USRMS), American Rock Mechanics Association. [7] S.Q. Yang, D.S. Yang, H.W. Jing, Y.H. Li, S.Y. Wang, An experimental study of the fracture coalescence behaviour of brittle sandstone specimens containing three fissures, Rock Mech. Rock Eng. 45 (4) (2012) 563–582. [8] O. Reyes, H.H. Einstein, Failure mechanisms of fractured rock - a fracture coalescence model (1991) in: 7th International Conference on Rock Mechanics, International Society for Rock Mechanics, 333–340. [9] L. Wong, H. Einstein, Crack coalescence in molded gypsum and carrara marble: Part 1. Macroscopic observations and interpretation. Rock Mech. Rock Eng. 42 (3) (2009) 475–511. [10] R. Wong, K. Chau, C. Tang, P. Lin, Analysis of crack coalescence in rock-like materials containing three flaws–part i: experimental approach, Int. J. Rock Mech. Min. Sci. 38 (7) (2001) 909–924. [11] Sagong, M., Bobet, A. (2002) Coalescence of multiple flaws in a rock-model material in uniaxial compression, Int. J. Rock Mech. Min. Sci. 39 (2) (2002) 229–241. [12] Ma et al. (2018) An investigation of non-straight fissures cracking under uniaxial compression, Eng. Fracture Mechanics 191, 15: 300-310 41 42 Bibliografía [13] Bobet, A., Einstein, H.H., (1998b). Fracture coalescence in rock-type materials under uniaxial and biaxial compression. Int. J. Rock Mech. Min. Sci. 35 (7), 863–888. [14] Gonçalves da Silva, B., Einstein, H.H., (2013) Modeling of crack initiation, propagation and coalescence in rocks. Int. J. Fract. 182 (2), 167–186. [15] Rabczuk, T. Belytschko T. (2004) Cracking particles: A simplified meshfree method for arbitrary evolving cracks, Int. J. Numer. Meth. Eng. 61 (13) (2004) 2316–2343. [16] Rabczuk, T. Belytschko T. (2007) A three-dimensional large deformation meshfree method for arbitrary evolving cracks, Comput. Methods Appl. Mech. Eng. 196 (29-30) 2777–2799. [17] Rabczuk, T., Zi, G., Bordas, S., Nguyen-Xuan, H. (2010) A simple and robust three-dimensional cracking-particle method without enrichment, Comput. Methods Appl. Mech. Eng. 199 (37) 2437–2455. [18] Rabczuk, T. Ren, H. (2017) A peridynamics formulation for quasi-static fracture and contact in rock, Eng. Geol. [19] Ren, H. Zhuang, X., Cai, Y. Rabczuk, T. (2016) Dual-horizon peridynamics, Int. J. Numer. Meth. Eng. 108 (12) (2016) 1451–1476. [20] Ren, H. Zhuang, X., Rabczuk, T. (2017) Dual-horizon peridynamics: A stable solution to varying horizons, Comput. Methods Appl. Mech. Eng. 318 762–782. [21] J. Lemaitre, R. Desmorat (2005) Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures, Springer-Verlag, Berlin. [22] Z.P. Bažant, T.G.P. Pijaudier-Cabot (1988) Nonlocal continuum damage, localization instability and convergence, J Appl. Mech. 55:287–293. [23] M. Jirásek (1998) Nonlocal models for damage and fracture: Comparison of approaches, Int. J. Solids Struct. 35:4133–4145. [24] S. Forest (2009) Micromorphic approach for gradient elasticity, viscoplasticity, and damage, J. Engnr. Mech., 135:117–131. [25] B.J.Dimitrijevic, K. Hackl (2011) A regularization framework for damage-plasticity models via gradient enhancement of the free energy, Int. J. Numer. Methods Biom. Engrg. 27:1199–1210. [26] R. Peerlings, M. Geers, R. de Borst, W. Brekelmans (2001) A critical comparison of non local and gradient-enhanced softening continua, Int. J. Solids Struct. 38:7723–7746. [27] P. Areias, M.A. Msekh, T. Rabczuk (2016) Damage and fracture algorithm using the screened Poisson equation and local remeshing, Eng. Fracture Mech. 158:116–143. [28] N. Moës, J. Dolbow, T. Belytschko (1999) A finite element method for crack growth without remeshing, Int. J. Numer. Methods Engrg. 46(1):131–150. [29] J. Dolbow, N. Moeës, T. Belytschko (2001) An extended finite element method for modeling crack growth with contact, Comput. Meth. Appl. Mech. Engng. 190:6825–6846. [30] T.P. Fries, T. Belytschko (2010) The extended/generalized finite element method: an overview of the method and its applications, Int. J. Numer. Methods Engrg. 84(3):253–304. [31] J.C. Simo, J. Oliver, F. Armero (1993) An analysis of strong discontinuities induced by strainsoftening in rate-independent inelastic solids, Comp. Mech. 12:277–296. Bibliografía 43 [32] C. Linder, F. Armero (2007) Finite elements with embedded strong discontinuities for the modeling of failure in solids, Int. J. Numer. Methods Engrg. 72:1391–1433. [33] J. Oliver, A. Huespe, S. Blanco, D. Linero (2006) Stability and robustness issues in numerical modeling of material failure with the strong discontinuity approach, Computer Methods in Applied Mechanics and Engineering, 195:7093–114. [34] F. Armero, C. Linder (2008) New finite elements with embedded strong discontinuities for finite deformations, Comput. Methods Appl. Mech. Engrg. 197:3138–3170. [35] V. Carollo (2018) A variational approach to brittle fracture and cohesive delamination: modelling and technological applications. PhD Thesis, IMT Lucca- Universidad de Sevilla. [36] X.P.Xu, A. Needleman (1994) Numerical simulation of fast crack growth in brittle solids, J. Mech. Phys. Solids 42:1397–1434. [37] M. Ortiz, A. Pandolfi (1999) Finite deformation irreversible cohesive elements for threedimensional crack-propagation analysis, Int. J. Num. Meth. Engng. 44:1267–1282. [38] K. Park, G.H. Paulino, and J.R. Roesler (2009) A unified potential-based cohesive model of mixed-mode fracture, J. Mech. Phys. Solids 57(6):891–908. [39] M. Paggi, P. Wriggers (2012) Stiffness and strength of hierarchical polycrystalline materials with imperfect interfaces, J. Mech. Phys. Solids 60:557–572. [40] J. Reinoso, M. Paggi (2014) A consistent interface element formulation for geometrical and material nonlinearities, Comp. Mech. 54:1569–1581. [41] A. Infuso, M. Corrado, M. Paggi (2015) Image analysis of polycrystalline solar cells and modelling of intergranular and transgranular cracking, J. Eur. Ceramic Society 34:2713–2722. [42] M. Paggi, J. Reinoso (2015) An anisotropic large displacement cohesive zone model for fibrillar and crazing of interfaces, Int. J. Solids Struct. 69-70:106–120. [43] G.A. Francfort, J.J. Marigo (1998) Revisiting brittle fracture as an energy minimization problem, J. Mech. Phys. Solids 46(8):1319–1342. [44] H. Amor, J.J. Marigo, C. Maurini (2009) Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments, J. Mech. Phys. Solids 57:1209–1229. [45] A.A Griffith (1921) The phenomena of rupture and flow in solids, Philos. Trans. Royal Soc. Lond. A 221:163–198. [46] B. Bourdin, G.A. Francfort, J-J. Marigo (2000) Numerical experiments in revisited brittle fracture, J. Mech. Phys. Solids 48(4):797–826. [47] L. Ambrosio, V.M. Tortorelli (1990) Approximation of functional depending on jumps by elliptic functional via C-convergence, Commun. Pure Appl. Math. 43:999–1036. [48] L. Ambrosio, V.M. Tortorelli (1992) On the approximation of free discontinuity problems, Boll. Un. Mat. Ital. B(7)6(1):105–123. [49] G. Dal Maso (1993) Introduction to Γ -convergence. Progress in nonlinear differential equations and their applications, Birkhäuser. [50] B. Bourdin, G.A. Francfort, J.J. Marigo (2008) The variational approach to fracture, J. Elast. 91(1-3):5–148.