Detección y monitorización virtual del crecimiento de grietas en paneles de nanocompuestos reforzados con nanotubos de carbono usando XFEM

Equation Chapter 1 Section 1 Dpto. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Trabajo Fin de Máster Máster en Ingeniería Industrial Sevilla, 2020 Autor: José Vargas Cruz Tutor: Luis Rodríguez de Tembleque Solano Detección y monitorización virtual del crecimiento de grietas en paneles de nanocompuestos reforzados con nanotubos de carbono usando XFEM iii Trabajo Fin de Máster Máster en Ingeniería Industrial Detección y monitorización virtual del crecimiento de grietas en paneles de nanocompuestos reforzados con nanotubos de carbono usando XFEM Autor: José Vargas Cruz Tutor: Luis Rodríguez de Tembleque Solano Profesor Titular de Universidad Dpto. de Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 v Trabajo Fin de Máster: Detección y monitorización virtual del crecimiento de grietas en paneles de nanocompuestos reforzados con nanotubos de carbono usando XFEM Autor: José Vargas Cruz Tutor: Luis Rodríguez de Tembleque Solano El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2020 El Secretario del Tribunal vii A mis padres y a mi hermano A María A Bru ix Agradecimientos Porque solo con constancia y trabajo se consiguen las cosas. Se termina una larga etapa. Tras seis años de esfuerzo se ve la luz al final del túnel. Y nada de esto podría haber sido posible sin vuestra ayuda. Gracias mamá, porque gracias a ti he aprendido lo que es saber luchar, porque eres el mejor ejemplo de mujer luchadora que puede existir. Gracias papá, por tu apoyo y por haberme enseñado desde pequeño que con constancia y esfuerzo se consiguen las cosas mas díficiles de esta vida. Gracias hermano, por haber alcanzado tantas metas juntos y por superar cualquier obstáculo que se nos haya presentado en el camino. Por eso y mucho más, gracias. Gracías María, por ser uno de los pilares más importantes de mi vida. Por el sacrificio que has hecho por mi en muchos momentos, por esos domingos de quedarnos en casa por mis trabajos. Simplemente gracias por ser como eres conmigo. Estaré eternamente agradecido. No podía dejar atrás a mi Bruno, porque me ha ayudado en los momentos mas difíciles de mi vida. Simplemente con mirarme te das cuenta de como me siento, y me ofreces esos abracitos y esos besos que tanto me gustan. Gracias Bru. Por último, gracias al equipo de trabajo de este proyecto. Gracias por acogerme desde el primer día como un miembro más del equipo, por transmitirme vuestra sabiduría y conocimientos y por haberme dado esta magnífica oportunidad de trabajar codo con codo con vosotros. Gracias Luis. José Vargas Cruz Sevilla, 2020 4.1.2 Análisis electroestático ..................................................................................................................... 39 4.1.2.1 Influencia del tamaño de grieta para varias orientaciones. Grieta impermeable ................ 40 4.1.2.2 Influencia de la permitividad eléctrica de la grieta para varias orientaciones. Grieta permeable 43 4.1.2.3 Medición y detección de varias grietas.................................................................................... 45 4.1.3 Análisis piezorresistivo ...................................................................................................................... 50 4.2 Monitorización virtual del crecimiento de grietas .................................................................................. 50 4.2.1 Placa con grieta interna .................................................................................................................... 51 4.2.1.1 Análisis mecánico ...................................................................................................................... 51 4.2.1.2 Análisis electroestático ............................................................................................................. 56 4.2.1.3 Análisis piezorresistivo .............................................................................................................. 59 4.2.2 Viga a flexión en tres puntos con defecto de borde ....................................................................... 59 4.2.2.1 Análisis mecánico ...................................................................................................................... 59 4.2.2.2 Análisis electroestático ............................................................................................................. 63 4.2.2.3 Análisis de la variación de la resistencia eléctrica frente a la flecha en el caso de viga apoyada sometida a flexión en tres puntos ................................................................................................ 67 5 Resumen, conclusiones y trabajos futuros ........................................................................................... 69 5.1 Resumen ..................................................................................................................................................... 69 5.2 Conclusiones .............................................................................................................................................. 69 5.3 Futuros trabajos ........................................................................................................................................ 70 Referencias.................................................................................................................................................. 71 Anexo A. Códigos de ANSYS ........................................................................................................................ 73 Anexo B. Códigos de MATLAB ..................................................................................................................... 75 xvii ÍNDICE DE FIGURAS Figura 1; Clasificación general de materiales compuestos [2] 1 Figura 2: Representación de la sección longitudinal de una viga de material compuesto con capacidad para detectar delaminación y deformaciones [4] 3 Figura 3: Esquema de sistema SHM en un material multifuncional compuesto con CNTs [5] 3 Figura 4: Estructuras de nanotubos de pared simple y múltiple [8]. 6 Figura 5: Resistencia a tracción (a), resistencia a flexión (b) y módulo de Young (c) para distintos porcentajes de fibras en MWCNT/epoxy [8]. 7 Figura 6: Representación de los mecanismos de salto electrónico (a) y red conductora (b). 8 Figura 7: Variación de la conductividad eléctrica de un nanocompuesto polimérico en funcón de la fracción volumétrica de fibras [10] 8 Figura 8: Representación de modos de carga I (a), II (b) y III (c) [11]. 11 Figura 9: Coordenadas polares ligadas al vértice de la grieta [11]. 12 Figura 10: Evolución de la tenacidad en función de la temperatura en aceros [11]. 14 Figura 11: Valores típicos de 𝐾𝐼𝐶 y límite elástico para diferentes materiales [11]. 14 Figura 12: Grieta centrada en placa infinita en modo I (a), en modo II (b) y en modo III (c) [11]. 15 Figura 13: Grieta inicial (a) y grieta propagada (b) [11]. 15 Figura 14: Esquema de modo de carga mixto [11]. 18 Figura 15: Representación de una placa infinita sometida a un modo de carga mixto [11]. 20 Figura 16: Evolución de la tensión circunferencial respecto al ángulo de propagación [11]. 20 Figura 17: SIFs en el vértice de grieta antes (a) y después (b) de la propagación [11]. 21 Figura 18: Geometría del dominio 𝛺 con una discontinuidad Г𝑐 [12]. 22 Figura 19: Representación esquemática de los diferentes nodos y elementos en la discretización [12]. 23 Figura 20: Vértice de grieta dentro del elemento (izq) y en el borde del elemento (der) [13]. 24 Figura 21: Nodos en el método Phantom-Node [13]. 25 Figura 22: Cálculo y definición de los valores 𝜓 y 𝜙 [13]. 26 Figura 23: Definición del radio de enriquecimiento [13]. 27 Figura 24: Definición de snap-tolerance [13]. 27 Figura 25: Evaluación de los critgerios STTMAX y PSMAX entorno al vértice de grieta [13]. 28 Figura 26: Ley cohesiva rígida lineal [13]. 28 Figura 27: Esquema general de monitorización del crecimiento de grieta. 32 Figura 28: Diagrama de flujo de la resolución del problema mecánico con ANSYS. 33 Figura 29: Diagrama de flujo de la resolución del problema electroestático y piezorresistivo con ANSYS. 34 Figura 30: Esquema problema mecánico de placa con grieta centrada inclinada bajo esfuerzo uniforme (𝜎𝑦𝑦). 36 Figura 31: Detalle del mallado de la placa para análisis estático usando XFEM. 37 Figura 32: Estudio de convergencia de los valores analíticos 𝐾𝐼 y 𝐾𝐼𝐼. 37 Figura 33: Influencia de la orientación de la grieta (α) en los SIFs. 38 Figura 34: Influencia tamaño de grieta en el cálculo de los SIFs KI (a) y KII (b). 38 Figura 35: : Esquema problema eléctrico de placa con grieta centrada inclinada sometida a diferencia de potencial (∆ϕ). 39 Figura 36: Influencia del tamaño de la grieta (𝐿/𝑎) y orientación de la grieta (𝛼) en la variación de la resistencia eléctrica. 40 Figura 37: Influencia del tamaño de grieta sobre la distribución de potencial eléctrico para 𝛼=0 º. 41 Figura 38: Influencia del tamaño de grieta sobre la distribución de potencial eléctrico para 𝛼=45 º 42 Figura 39: Influencia de las condiciones de permitividad eléctrica en la grieta (𝜅𝑐/𝜅𝑚) y la orientación (𝛼) en la variación de la resistencia eléctrica (𝑅). 43 Figura 40: Influencia de la permitidad sobre la distribución de potencial eléctrico para 𝛼=0 º y 𝐿/𝑎=1.5. 44 Figura 41: Influencia de la permitidad sobre la distribucuión de potencial eléctrico para 𝛼=45 º y 𝐿/𝑎=1.5. 45 Figura 42: Esquema problema eléctrico para la medición y detección de dos grietas equivalentes a la original. 46 Figura 43: Influencia de la separación 𝛿𝑥 y 𝛿𝑦 en la variación de la resistencia eléctrica del problema con grieta equivalente. 46 Figura 44: Potencial eléctrico para 𝛿𝑦/𝑎2=0 (a) y 𝛿𝑦/𝑎2=0.05 (b). 47 Figura 45: Potencial eléctrico para 𝛿𝑦/𝑎2=0.1 (a) y 𝛿𝑦/𝑎2=0.5 (b). 48 Figura 46: Potencial eléctrico para 𝛿𝑦/𝑎2=1 (a) y 𝛿𝑦/𝑎2=2 (b). 49 Figura 47: Comparación entre análisis piezorresistivo y electroestático. 50 Figura 48: Esquema problema mecánico de placa con grieta inclinada bajo desplazamiento uniforme (𝑢𝑦). 52 Figura 49: Definición de áreas en problema mecánico de placa traccionada para 𝛼=45 º. 52 Figura 50: Mallado de placa sometida a tracción para 𝛼=45 º. 53 Figura 51: Isocontornos de desplazamientos 𝑢𝑦 (a) y tensión equivalente Von Mises 𝜎𝑒𝑞 (b) para 𝛼=0 º. 54 Figura 52: Isocontornos de desplazamientos 𝑢𝑦 y 𝛼=0 º para 𝑎/𝑎0=3 (a), 𝑎/𝑎0=6 (b) y 𝑎/𝑎0=9 (c). 54 Figura 53: Isocontornos de desplazamientos 𝑢𝑦 (a) y tensión equivalente Von Mises 𝜎𝑒𝑞 (b) para 𝛼=45 º. 55 Figura 54: Isocontornos de desplazamientos 𝑢𝑦 y 𝛼=45 º para 𝑎/𝑎0=3 (a), 𝑎/𝑎0=6 (b), 𝑎/𝑎0=9 (c) y 𝑎/𝑎0=12 (d). 55 Figura 55: Esquema problema eléctrico de placa con grieta inclinada genérica sometida a tracción. 56 Figura 56: Potencial eléctrico para (a) 𝛼=0 º y (b) 𝛼=45 º el problema de placa con grieta interna. 57 Figura 57: Evolución de distribución del potencial eléctrico placa a tracción para 𝛼=0 º. 57 Figura 58: Evolución de distribución del potencial eléctrico placa a tracción para 𝛼=45 º. 58 Figura 59: Análisis de la influencia del tamaño de grieta en la variación de la resistencia eléctrica en el caso de placa sometida a tracción para 𝛼=[0,15,30,45,60,75] º. 58 Figura 60: Comparación análisis piezorresistivo y electroestático para placa con grieta interna. 59 Figura 61 Esquema problema mecánico de placa con grieta inclinada sometida a flexión. 60 xix Figura 62: Definición de áreas en problema mecánico de placa sometida a flexión para 𝛼=45 º. 60 Figura 63: Mallado de placa sometida a flexión para 𝛼=45 º. 61 Figura 64: (a) Isocontornos de desplazamientos 𝑢𝑦 y (b) tensión equivalente Von Mises 𝜎𝑒𝑞 para 𝛼=0 º. 61 Figura 65: Isocontornos de desplazamientos 𝑢𝑦 y 𝛼=0 º para (a) 𝑎/𝑎0=1 , (b) 𝑎/𝑎0=2 , (c) 𝑎/𝑎0=3 y (d) 𝑎/𝑎0=4 . 62 Figura 66: (a) Isocontornos de desplazamientos 𝑢𝑦 y (b) tensión equivalente Von Mises 𝜎𝑒𝑞 para 𝛼=45 º. 62 Figura 67: Isocontornos de desplazamientos 𝑢𝑦 y 𝛼=45 º para (a) 𝑎/𝑎0=2, (b) 𝑎/𝑎0=3 y (c) 𝑎/𝑎0=4, 63 Figura 68: Esquema problema eléctrico de placa con grieta inclinada genérica sometida a flexión. 64 Figura 69: Distribución de potencial eléctrico para el problema de placa sometida a flexión. 64 Figura 70: Evolución de distribución del potencial eléctrico placa a flexión para 𝛼=0 º. 65 Figura 71: Evolución de distribución del potencial eléctrico placa a flexión para 𝛼=45 º. 66 Figura 72: Análisis de la influencia del tamaño de grieta en la variación de la resistencia eléctrica en el caso de viga sometida a flexión en tres puntos para (a) 𝛼=0 º y (b) distintos valores de 𝛼. 67 Figura 73: Variación de la resistencia eléctrica en función del desplazamiento vertical de la viga. 67 xxi Notación A Superficie de la grieta A* Semiancho de la grieta 𝑎𝑖 Vector desplazamiento asociado al dominio de la grieta 𝐶 Vector quiral dA Incremento de superficie de grieta E Módulo elástico 𝐸𝑠 Energía para crear nueva superficie de grieta ∆𝐸𝑖𝑛𝑡 Energía interna f Vector de fuerzas externas 𝐹𝑏 Fuerzas por unidad de volumen aplicadas en el sistema 𝐹𝑗 (𝑥) Funciones de singularidad G Tasa de liberación de energía 𝐺𝑐 Tasa de liberación de energía crítica H (x) Funciones Heaviside I Intensidad de corriente J Densidad de corriente K Factor de intensidad de tensiones 𝐾𝑐 Factor de intensidad de tensiones crítico L Semilargo de la placa 𝑁𝑖(𝑥) Funciones de forma asociadas a cada nodo 𝑛ℎ Número de nodos asociados a los elementos de la superficie de la grieta 𝑛𝑡 Número de nodos asociados al dominio del vértice de la grieta R Resistencia eléctrica 𝑅0 Resistencia eléctrica de la placa sin deformar r Tensor de cambio relativo en la resistividad t Espesor de la placa 𝑡𝑖1 Vector desplazamiento asociado a los elementos del vértice de grieta 𝑇1 Funciones de enriquecimiento del vértice de grieta 𝑢𝑥 Desplazamiento en el plano de la placa en dirección “X” 𝑢𝑦 Desplazamiento en el plano de la placa en dirección “Y” W Semiancho de la placa 𝛼 Ángulo de inclinación de la grieta Г Contorno del problema generalizado Г𝑢 Contorno de desplazamientos Г𝑡 Contorno de tensiones Г𝑐 Contorno de la grieta 𝛾𝑠 Energía superficial 𝜅 Permitividad eléctrica 𝜅𝑐 Permitividad eléctrica de la grieta 𝜅𝑚 Permitividad eléctrica del medio 𝜆𝑖𝑗 Componentes de la matriz piezorresistiva 𝜈 Coeficiente de Poisson Π Matriz piezorresistiva 𝜌𝑒𝑓𝑓 Matriz resistividad efectiva 𝜌𝑒𝑓𝑓 0 Matriz resistividad efectiva de la placa sin deformar ϕ Potencial eléctrico Ω Dominio del problema generalizado 𝜎 Tensor de tensiones de Cauchy 𝜎𝑥𝑥 Tensión axial dirección “X” 𝜎𝑦𝑦 Tensión axial dirección “Y” 𝜎𝑥𝑦 Tensión tangencial plano “XY” 𝜎𝜃𝜃 Tensión circunferencial 𝜎𝑟𝜃 Tensión tangencial circunferencial 1 1 INTRODUCCIÓN 1.1 Motivación De forma general, un material compuesto puede definirse como un sistema integrado por una combinación de dos o más materiales que difieren en forma y composición química, y que además son insoluble entre sí. Los materiales compuestos se han convertido en el principal foco de interés en la ingeniería ya que combinan propiedades y prestaciones de los materiales constituyentes, presentando notables mejoras. La mayor parte de los materiales compuestos están formados por una matriz continu que rodea y aglomera las fases dispersas, también conocidas como refuerzos. De esta forma, en función de la naturaleza de los refuerzos, los materiales compuestos se pueden clasificar en tres grandes grupos: compuestos con partículas como refuerzos, reforzados con fibras y compuestos estructurales. Sin embargo, los materiales compuestos también pueden clasificarse en función de la naturaleza de la matriz, existiendo los siguientes tipos: MMC (Metalic Matrix Composites, por sus siglas en inglés), CMC (Ceramic Matrix Composites) y PMC (Polymeric Matrix Composites) [1]. Figura 1; Clasificación general de materiales compuestos [2] Por otro lado, existe un concepto que se va afianzando cada vez más en la actualidad, y es el de materiales multifuncionales o materiales inteligentes. Este tipo de materiales se definen como una combinación de materiales que desarrollan sus funciones de una forma inteligente en función de los cambios que experimente el medio que los rodea. En la actualidad, existe una infinidad de aplicaciones tecnológicas que requieren el uso de distintos tipos de sensores. Estas aplicaciones requieren que se lleve a cabo un desarrollo y estudio exhaustivo de estos materiales multifuncionales que presentan capacidades sensoriales. 2 El rápido crecimiento de estas tecnologías y su interés en la industria se ha visto reflejada en el surgimiento de nuevas ramas de investigación y tecnologías que tienen un gran impacto en diversas áreas científicas. Se trata de la nanotecnología. Esta ciencia persigue crear nuevos materiales que tengan capacidad para lograr desafíos tecnológicos que han ido surgiendo en los últimos años debido al incremento de la necesidad tecnológica actual. Generalmente los sistemas inteligentes están formados por sensores, actuadores y controladores. En el caso de la fibra óptica, esta es capaz de reaccionar ante cambios de temperatura, deformación, presión e incluso campos eléctricos, lo que la convierte en un excelente sensor para ser utilizado en el sector de las telecomunicaciones. Las aleaciones con memoria de forma (AMF) son los materiales más comúnmente usados como actuadores, destacando su uso en instrumentos quirúrgicos. Sin embargo, existen casos particulares como los compuestos de matriz polimérica con nanotubos de carbono como refuerzo, que se agrupan dentro de lo que se conoce como estructuras intrínsecamente inteligentes [3]. La principal característica de estos materiales es que no requieren de dispositivos como sensores y actuadores, ya que se combinan las propiedades mecánicas, térmicas, eléctricas y piezorresistivas de manera que son capaces de detectar por sí solas daños mecánicos o térmicos, esfuerzos y deformaciones. Cuando se añaden en pequeñas concentraciones y su dispersión es adecuada, se ha podido comprobar que suponen una mejora en términos mecánicos, así como mejoran propiedades multifuncionales como la conductividad eléctrica. Estos materiales se han convertido en el epicentro de numerosos estudios científicos y aplicaciones industriales, debido a las excelentes y peculiares propiedades eléctricas, térmicas y mecánicas que presentan, entre las cuales destaca su propiedad piezorresistiva. La piezorresistividad es una propiedad que tienen ciertos materiales por la cual se genera una respuesta eléctrica ante un estímulo mecánico, tal y como se verá detalladamente en el primer capítulo del presente proyecto. Una alternativa potencial para el desarrollo de materiales piezorresistivos es la utilización de nanocompuestos formados una matriz epoxy y nanotubos como refuerzos. La cuestión es que en la actualidad existen muchos materiales metálicos y poliméricos que no cumplen con los requisitos para las aplicaciones que se demandan. Dicha propiedad piezorresistiva de este tipo de nanocompuestos hace posible plantearse su utilización para aplicaciones que abarquen industrias como la aeronática, automotriz o incluso el sector energético. Esta propiedad puede aprovecharse de forma que se utilicen como sensores para la detección de daño, ya que un componente fabricado con este tipo de material, al estar sometido a una serie de solicitaciones mecánicas experimentará unas deformaciones que se traducirán finalmente en una variación de la resistencia eléctrica. Si fuera posible relacionar dicha variación de la resistencia eléctrica con fenómenos físicos como defectos en el componente u otro tipo de daño, sería posible establecer criterios y métodos que permitieran obtener información acerca del estado de dicho componente o estructura. Actualmente existe una disciplina llamada SHM (Structural Health Monitoring, por sus siglas en inglés) que constitutye una práctica común en mútiples actividades industriales. Esta disciplina abarca la aplicación de ensayos no destructivos (NDT, por sus siglas en inglés) y detección de daño para evaluar la solidez de infraestructuras y el comportamiento del mantenimiento basado en condición que permite aumentar su esperanza de vida. Particularmente, la mayoría de las investigaciones han centrado sus esfuerzos en detectar daño en estructuras a través de variaciones en la rigidez de los componentes estructurales. Dichas variaciones a menudo se originan por procesos de degradación por envejecimiento, corrosión, fatiga, existencia de grietas o por sucesos accidentales. Sin embargo, el número de estudios en la aplicación de técnicas SHM son escasos [5]. Los nuevos materiales multifuncionales, como el hormigón “inteligente” o materiales compuestos “autosensitivos”, ofrecen una innovadora alternativa para monitorear el estado de las estructuras dentro del campo de la SHM [5]. En este contexto, el presente proyecto se desarrolla como una estrategia de monitorización de la integridad de estructuras formadas por nanocompuestos con nanotubos de carbono como refuerzos, persiguiendo el objetivo de establecer un procedimiento que permita evaluar el daño en estructuras fabricadas con este material. 3 Figura 2: Representación de la sección longitudinal de una viga de material compuesto con capacidad para detectar delaminación y deformaciones [4]. Figura 3: Esquema de sistema SHM en un material multifuncional compuesto con CNTs [5]. 1.2 Descripción del problema y objetivos En este proyecto se pretende estudiar el comportamiento de MWCNTs como refuerzos en nanocompuestos, así como interpretar los resultados obtenidos para logar llegar a conclusiones que permitan desarrollar una metodología de validación y control de paneles fabricados con nanocompuestos reforzados con nanotubos de carbono. Para ello se propone una simulación numérica que permita realizar una monitorización virtual del daño en estructuras tipo panel, poniendo en juego el efecto piezorresistivo que presenta este tipo de nanocompuestos. Se trata de un método no invasivo que podría implementarse en numerosas aplicaciones industriales, permitiendo no solo detectar defectos en estos componentes, sino también establecer una monitorización en tiempo real de dichos defectos. De esta forma se podría disponer de un regristro constante de la evolución de dichos defectos en la estructura, pudiéndose determinar criterios para proceder a su reparación o declaración de no conformidad para su uso. 11 3 FORMULACIÓN DE ELEMENTOS FINITOS PARA LA DETECCIÓN Y MONITORIZACIÓN DE GRIETAS EN PANELES DE NANOCOMPUESTOS REFORZADOS CON NANOTUBOS DE CARBONO ste tercer capítulo se centra en el desarrollo de la formulación de elementos finitos orientada a la detección y monitorización de grietas. Por un lado, se tratará la formulación del problema mecánico, trantando conceptos típicos de la Mecánica de la Fractura avanzada como la tenacidad y energía de fractura, modo de fractura, métodos analíticos para establecer criterios de crecimiento de grieta, cómo se desarrolla la formulación XFEM y cómo el software ANSYS modela este tipo de problemas. Por otro lado, se trata la formulación y resolución del problema eléctrico mediante el software de simulación ANSYS, abarcando cómo se realiza el cálculo de la respuesta piezorresistiva, cómo se resuelve el campo eléctrico y el cálculo de la resistencia eléctrica en paneles con presencia de grietas o daños. 3.1 Formulación del problema mecánico En primer lugar, se dará paso al desarrollo de la formulación del problema mecánico, en el cuál se abarcarán aspectos como la energía de fractura, tenacidad de fractura, la formulación del problema mediante la metodología XFEM y la implementación de toda la formulación en el software ANSYS. 3.1.1 Energía de fractura y tenacidad a fractura Con el objetivo de resolver el caso de una placa con una grieta inicial, Irwin estudió tres modos de carga distintos:  Modo I: también llamado modo de apertura, ya que tiende a abrir los bordes de la grieta.  Modo II: también llamado modo de deslizamiento, ya que tiende a deslizar los bordes de la grieta uno sobre otro.  Modo III: también conocido como modo de cizallamiento, desplazando la grieta fuera de su plano. (a) (b) (c) Figura 8: Representación de modos de carga I (a), II (b) y III (c) [11]. E 12 La mecánica de sólidos clásica puede ser aplicada a cada uno de los modos anteriores para calcular el campo de tensiones y deformaciones. Para ello, se utilizarán coordenadas polares vinculadas al vértice de la grieta. Los tres modos de fractura anteriores se diferencian en las condiciones de contorno impuestas. A continuación, se muestran las expresiones del campo de tensiones para cada uno de los modos descritos: MODO I 𝜎𝑥𝑥=𝐶 √𝑟𝑐𝑜𝑠𝜃2[1−𝑠𝑖𝑛3𝜃 2𝑠𝑖𝑛𝜃2]+𝑓(𝑟0) 𝜎𝑦𝑦=𝐶 √𝑟𝑐𝑜𝑠𝜃2[1+𝑠𝑖𝑛3𝜃 2𝑠𝑖𝑛𝜃2]+𝑓(𝑟0) 𝜎𝑥𝑦=𝐶 √𝑟𝑐𝑜𝑠𝜃2𝑐𝑜𝑠3𝜃 2𝑠𝑖𝑛𝜃2+𝑓(𝑟0) 𝑢𝑥=𝐶(1+𝜈) 𝐸√𝑟𝑐𝑜𝑠𝜃2[𝜅−1+2𝑠𝑖𝑛2𝜃2]+𝑓(𝑟) 𝑢𝑦=𝐶(1+𝜈) 𝐸√𝑟𝑠𝑖𝑛𝜃2[𝜅+1+2𝑐𝑜𝑠2𝜃2]+𝑓(𝑟) (3 – 1) donde 𝜅 es una expresión del coeficiente de Poisson 𝜈 y es igual a 3−𝜈 1+𝜈 en el caso de deformación plana e igual a 3−4𝜈 para deformación plana. Figura 9: Coordenadas polares ligadas al vértice de la grieta [11]. MODO II 𝜎𝑥𝑥=−𝐶 √𝑟𝑠𝑖𝑛𝜃2[2+𝑐𝑜𝑠3𝜃 2𝑐𝑜𝑠𝜃2]+𝑓(𝑟0) 𝜎𝑦𝑦=𝐶 √𝑟𝑐𝑜𝑠𝜃2𝑐𝑜𝑠3𝜃 2𝑠𝑖𝑛𝜃2+𝑓(𝑟0) 𝜎𝑥𝑦=𝐶 √𝑟𝑐𝑜𝑠𝜃2[1−𝑠𝑖𝑛3𝜃 2𝑠𝑖𝑛𝜃2]+𝑓(𝑟0) 𝑢𝑥=𝐶(1+𝜈) 𝐸√𝑟𝑠𝑖𝑛𝜃2[𝜅+1+2𝑐𝑜𝑠2𝜃2]+𝑓(𝑟) 𝑢𝑦=𝐶(1+𝜈) 𝐸√𝑟𝑐𝑜𝑠𝜃2[𝜅−1+2𝑠𝑖𝑛2𝜃2]+𝑓(𝑟) (3 – 2) 13 MODO III 𝜎𝑥𝑧=−𝐶 √𝑟𝑠𝑖𝑛𝜃2+𝑓(𝑟0) 𝜎𝑦𝑧=𝐶 √𝑟𝑐𝑜𝑠𝜃2+𝑓(𝑟0) 𝑢𝑥=4𝐶(1+𝜈) 𝐸√𝑟𝑠𝑖𝑛𝜃2+𝑓(𝑟) (3 – 3) Una vez desarrollada la formulación del campo de tensiones y deformaciones de los tres modos de fractura, nos centraremos en el modo I y vamos a analizar el término 𝜎𝑦𝑦 (término dominante), que es el campo de tensiones que caracteriza la concentración de tensiones en una placa sometida a tracción. El término dominante cerca del vértice de la grieta es 𝐶 √𝑟, que representa que existe una singularidad en el vértice de la grieta. Por esta razón el valor de la tensión no puede usarse para determinar si una determinada grieta inicial va a propargarse o no. Irwin planteó considerar “cómo de rápido la tensión tienda a infinito” cerca del vértice de la grieta, para lo cual definió el factor de intensidad de tensiones (SIF). 𝐾𝐼=𝑙𝑖𝑚𝑟→0(√2𝜋𝑟𝜎𝑦𝑦 𝑚𝑜𝑑𝑜 𝐼⎹ 𝜃=0)=𝐶√2𝜋 𝐾𝐼𝐼=𝑙𝑖𝑚𝑟→0(√2𝜋𝑟𝜎𝑥𝑦 𝑚𝑜𝑑𝑜 𝐼𝐼⎹ 𝜃=0)=𝐶√2𝜋 𝐾𝐼𝐼𝐼=𝑙𝑖𝑚𝑟→0(√2𝜋𝑟𝜎𝑦𝑧 𝑚𝑜𝑑𝑜 𝐼𝐼𝐼⎹ 𝜃=0)=𝐶√2𝜋 (3 – 4) De esta forma, se puede desarrollar la formulación genérica de la evolución de tensiones y desplazamientos cerca del vértice de grieta a partir de un SIF expresado en unidades 𝑀𝑃𝑎∙√𝑚 : 𝜎𝑚𝑜𝑑𝑜 𝑖=𝐾𝑖 √2𝜋𝑟𝑓𝑚𝑜𝑑𝑜 𝑖(𝜃) 𝑢𝑚𝑜𝑑𝑜 𝑖=𝐾𝑖√𝑟 2𝜋𝑔𝑚𝑜𝑑𝑜 𝑖(𝜃) (3 – 5) donde 𝑓 y 𝑔 son funciones definidas para cada modo independientes de la carga y geometría. De esta forma Irwin tuvo la idea de considerar el valor de los SIFs para analizar y detectar la propagación de las grietas. De hecho, experimentalmente se puede demostrar que para un determinado material que obedezca las hipótesis de LEFM (Linear Elastic Fracture Mechanics), la grieta se propaga si el valor del SIF correspondiente alcanza un cierto valor crítico 𝐾𝐶 llamado tenacidad de fractura. Con todo lo descrito anteriormente se puede establecer un criterio de propagación de grieta para el modo I: 𝐾𝐼<𝐾𝐼𝐶 : la grieta no se propaga 𝐾𝐼>𝐾𝐼𝐶 : la grieta se propaga (3 – 6) El término 𝐾𝐼𝐶 representa la tenacidad de fractura. Cabe destacar que bajo las hipótesis de LEFM, 𝐾𝐼 depende exclusivamente de la geometría y las condiciones de carga, mientras que 𝐾𝐼𝐶 depende solamente del material. Por otro lado, hay que señalar que los materiales dúctiles tienen una tenacidad de fractura mayor, aunque 14 normalmente las hipótesis de LEFM no se cumplen debido a que su comportamiento no es elástico. Algunos materiales tienen un comportamiento frágil a bajas temperaturas y un comportamiento dúctil a altas temperaturas. Para dichos materiales la tenacidad depende de la temperatura y existe una zona de transición, tal y como se muestra en la Figura 10. Figura 10: Evolución de la tenacidad en función de la temperatura en aceros [11]. Figura 11: Valores típicos de 𝐾𝐼𝐶 y límite elástico para diferentes materiales [11]. Otra cuestión interesante es cómo se puede determinar el valor de los SIFs, los cuales están asociados a la tenacidad del material. Tal y como se comentó previamente, estos factores dependen de la geometría y de las condiciones de carga. En la práctica, existen distintos métodos entre los que se encuentran el método analítico, el método de superposición, manuales específicos, método de elementos finitios (FEM) y el método de enfoque energético. El método analítico solamente puede utilizarse en cass simples en los cuales exista una grieta de tamaño 2𝑎 en una placa infinita. En estos casos los resultados analíticos para los diferentes SIFs son los que se muestran a continuación. 𝐾𝐼=𝜎∞√𝜋𝑎 𝐾𝐼𝐼=𝜏∞√𝜋𝑎 𝐾𝐼𝐼𝐼=𝜏∞√𝜋𝑎 (3 – 7) Para casos mas complejos, se utilizan otras expresiones para el cálculo de los SIFs tal y como se muestra a continuación. 𝐾𝐼=𝛽𝐼𝜎∞√𝜋𝑎 𝐾𝐼𝐼=𝛽𝐼𝐼𝜏∞√𝜋𝑎 𝐾𝐼𝐼𝐼=𝛽𝐼𝐼𝐼𝜏∞√𝜋𝑎 (3 – 8) donde 𝛽𝑖 son valores que dependen de la geometría y del tamaño de la grieta. 15 (a) (b) (c) Figura 12: Grieta centrada en placa infinita en modo I (a), en modo II (b) y en modo III (c) [11]. Tal y como se ha visto previamente, el tamaño de grieta crítico esta relacionado con la energía de fractura. Consideremos un caso de una placa infinito con una grieta centrada de longitud 2𝑎, y tomemos un sistema de coordenadas polar en uno de los vértice de la grieta en el estado inicial, es decir, antes de que la grieta se haya propagado (ver Figura 13). El campo de tensiones y desplazamientos pueden calcularse tomando los valores de 𝜃=0º y 𝑥=𝑟+𝑎: 𝜎𝑦𝑦 0(𝜃=0;𝑟=𝑥−𝑎)=𝜎∞√𝑎 √2(𝑥−𝑎) 𝑢𝑦0(𝜃=0;𝑟=𝑥−𝑎)=0 (3 – 9) Ahora consideremos que la grieta ha crecido y ha alcanzado una longitud 2(𝑎+∆𝑎). El campo de tensiones y desplazamiento arriba y debajo de los bordes de grieta se pueden calcular considerando valores 𝜃 y x tales que 𝜃=±𝜋 y 𝑥=𝑎+∆𝑎−𝑟: 𝜎𝑦𝑦 0(𝜃=𝜋 ; 𝑟=𝑎+∆𝑎−𝑥)=0 𝑢𝑦0(𝜃=±𝜋 ; 𝑟=𝑎+∆𝑎−𝑥)=±𝜎∞(1+𝜈)(𝜅+1) 𝐸√2√𝑎+∆𝑎√𝑎+∆𝑎−𝑥 (3 – 10) (a) (b) Figura 13: Grieta inicial (a) y grieta propagada (b) [11]. 16 De las expresiones anteriores se puede deducir que antes de la propagación existe un campo de desplazamiento 𝜎𝑦𝑦 0 distinto de cero pero un desplazamiento nulo, mientras que tras la propagación existe un desplazamiento 𝑢𝑦0 asociado a la apertura de la grieta pero la zona en cuestión está libre de tensiones. Por tanto, la estructura libera energía cuando el tensor de tensiones exhibe un trabajo durante la propagación de la grieta. Esta energía se calcula como se muestra a continuación: ∆𝐸𝑖𝑛𝑡=−4𝑑𝑧∫ ∫ 𝜎𝑦𝑦𝑑𝑢𝑦𝑑𝑥 𝑢𝑦 1 𝑢𝑦 0 𝑎+∆𝑎 𝑎 (3 – 11) donde el factor de 4 viene de considerar que existe una doble simetría alrededor de los ejes 𝑥 e 𝑦. Hay que destacar que existe un incremento de grieta de tamaño ∆𝑎 en cada lado de la grieta. Considerando un material elástico lineal, en el proceso de apertura de la grieta la tensión 𝜎𝑦𝑦 debe decrementarse linealmente con el desplazamiento 𝑢𝑦. La expresión una vez simplificada quedaría de la siguiente forma: ∆𝐸𝑖𝑛𝑡=−2𝑑𝑧∫ 𝜎𝑦𝑦 0𝑢𝑦1𝑑𝑥 𝑎+∆𝑎 𝑎 (3 – 12) A continuación se aplica el cambio de variable 𝑥=𝑎+∆𝑎𝑐𝑜𝑠2𝜃, con lo que la expresión de la energía liberada queda: ∆𝐸𝑖𝑛𝑡=−𝑑𝑧𝜎∞ 2√𝑎(𝑎+∆𝑎)(1+𝜈)(1+𝜅) 2𝐸 𝜋∆𝑎 (3 – 13) Esta energía se puede relacionar con el incremento de la superficie de la grieta tal que 𝑑𝐴=2∆𝑎𝑑𝑧. De esta forma se obtiene un concepto muy importante en la Mecánica de la Fractura, la tasa de liberación de energía 𝐺: 𝐺=−𝑑𝐸𝑖𝑛𝑡 𝑑𝐴 (3 – 14) Esta energía corresponde a la energía liberada por la estructura por unidad de incremento de superficie de la grieta. La anterior expresión puede modificarse tal que: 𝐺=−𝑑𝐸𝑖𝑛𝑡 𝑑𝐴 =lim ∆𝑎→0∆𝐸𝑖𝑛𝑡 2∆𝑎𝑑𝑧=𝜋𝑎𝜎∞ 2 4𝐸 (1+𝜈)(𝜅+1)=𝜋𝑎𝜎∞ 2 𝐸′=𝐾𝐼2 𝐸′ (3 – 15) 𝐸′=𝐸 para un estado de tensión plana 𝐸′=𝐸 1−𝜈2 para un estado de deformación plana Considerando un estado de tensión plana, es posible relacionar los SIFs con la tasa de liberación de energía 𝐺 , de forma que: 𝐺=𝐾𝐼2 𝐸′+𝐾𝐼𝐼2 𝐸′+𝐾𝐼𝐼𝐼 2(1+𝜈) 𝐸′ (3 – 16) 17 Por otro lado, la energía total puede asumirse como la suma de la energía elástica interna del material 𝐸𝑖𝑛𝑡 y la energía necesaria 𝐸𝑠 para crear una nueva superficie de grieta A: 𝐸=𝐸𝑖𝑛𝑡+𝐸𝑠 (3 – 20) Si se deriva la expresión anterior respecto a la superficie de grieta A y se considera que la energía total se conserva (no se considera disipación de energía): 𝜕𝐸 𝜕𝐴=𝜕𝐸𝑖𝑛𝑡 𝜕𝐴 +𝜕𝐸𝑠 𝜕𝐴=0 (3 – 21) Teniendo en cuenta la expresión de la tasa de liberación de energía 𝐺, la expresión anterior pasa a ser: 𝐺=−𝜕𝐸𝑖𝑛𝑡 𝜕𝐴 =𝜕𝐸𝑠 𝜕𝐴 (3 – 22) Si asumimos que la grieta está creciendo, el término 𝜕𝐸𝑠 𝜕𝐴 corresponde precisamente con el doble de la energía superficial 𝛾𝑠. 𝐺=−𝜕𝐸𝑖𝑛𝑡 𝜕𝐴 =𝜕𝐸𝑠 𝜕𝐴=2𝛾𝑠=𝐺𝑐 (3 – 23) La anterior expresión significa que para que una grieta crezca la estructura tiene que liberar una energía mayor o igual a la energía requerida para formar formar una grieta. Para materiales frágiles el valor de 𝐺𝑐 es dos veces la energía superficial y corresponden con la energía requerida para que se produzca clivaje. Por tanto, el criterio para que una grieta crezca se puede expresar en términos de SIF o en términos energéticos, tal y como se muestra a continuación: 𝐾𝐼≥𝐾𝐼𝐶 𝐺≥𝐺𝐶 (3 – 24) 3.1.1.1 Crecimiento de grieta. Modo de fractura mixto Una cuestión interesante que debe plantearse es: ¿en qué dirección crece la grieta? En los puntos siguientes se verán aspectos fundamentales a tener en cuenta como combinación de modos de carga, métodos de cálculo del ángulo de propagación de la grieta y relación entre la ley cohesiva y parámetros característicos de la Mecánica de fractura. 3.1.1.1.1 Combinación de modos I y II Anteriormente se ha determinado si la grieta crece o no cuando está cargada en modo I. En ese caso, por simetrñia la grieta se propaga en línea recta. Pero, ¿qué ocurriría si el modo de carga es mixto? ¿Cuál es la condición para que la grieta crezca? ¿En que dirección lo hace? Para dar respuesta a las preguntas anteriores, consideremos un caso de carga mixto tal y como se muestra en la Figura 14, en la cual la carga está orientada según una dirección que forma un ángulo 𝛽 respecto al plano de la grieta (un ángulo 𝛽=𝜋/2 correspondería al caso de modo I puro). Bajo estas condiciones de carga, la grieta se propaga según un ángulo 𝜃. 18 Figura 14: Esquema de modo de carga mixto [11]. A continuación, se va a establecer la formulación necesaria para evaluar los SIFs para las condiciones de carga previamente mencionadas. El tensor de tensiones es conocido en la dirección de la carga y puede expresarse referenciándolo a la grieta utilizando una matriz de rotación: (𝜎𝑥𝑥 𝜎𝑥𝑦 𝜎𝑥𝑦 𝜎𝑦𝑦)=(cos𝛽 sin𝛽 −sin𝛽 cos𝛽)𝑇(𝜎∞0 0 0)(cos𝛽 sin𝛽 −sin𝛽 cos𝛽) { 𝜎𝑥𝑥=𝜎∞ 𝑐𝑜𝑠2𝛽 𝜎𝑥𝑦=𝜎∞sin𝛽cos𝛽 𝜎𝑦𝑦=𝜎∞ 𝑠𝑖𝑛2𝛽 (3 – 25) Si se asume que la placa es infinita, los SIFs correspondientes a los modos I y II se obtienen a partir del campo de tensiones yy y xy. La determinación del ángulo de propagación de la grieta se puede llevar a cabo mediante dos métodos que se verán más adelante. {𝐾𝐼=𝜎∞ 𝑠𝑖𝑛2𝛽 √𝜋𝑎 𝐾𝐼𝐼=𝜎∞sin𝛽cos𝛽 √𝜋𝑎 (3 – 26) 3.1.1.1.2 Método de la tensión máxima circunferencial El método de la tensión circunferencial máxima fue propuesto por Erdogan y Sih en 1963 y está basado en las siguientes dos hipótesis: 1. La grieta crece en la dirección para la cual el SIF relacionado con el modo I es máxima. 2. La grieta crece solo si el SIF es mayor que la tenacidad medida para el modo I puro. Como el modo I en la nueva estructura está caracterizada por la nueva tensión 𝜎𝑦𝑦, por ejemplo 𝜎𝜃𝜃, el criterio de propagación de la grieta es el siguiente: (√2𝜋𝑟𝜎𝜃𝜃(𝑟,𝜃∗))≥𝐾𝐼𝐶 (3 – 27) 19 El valor de la dirección de propagación 𝜃∗ se obtiene del siguiente sistema: {𝜕𝜎𝜃𝜃 𝜕𝜃 ⎹𝜃∗=0 𝜕2𝜎𝜃𝜃 𝜕𝜃2⎹𝜃∗<0 (3 – 28) A partir de las expresiones del campo de tensiones para un problema 2D y modos I y II, aplicando el principio de superposición se obtienen las expresiones del campo de tensiones para modo mixto: 𝜎𝑥𝑥=1 √2𝜋𝑟[𝐾𝐼𝑐𝑜𝑠𝜃2−𝐾𝐼 2𝑠𝑖𝑛3𝜃 2𝑠𝑖𝑛𝜃2−2𝐾𝐼𝐼𝑠𝑖𝑛𝜃2−𝐾𝐼𝐼 2𝑐𝑜𝑠3𝜃 2𝑠𝑖𝑛𝜃2] 𝜎𝑦𝑦=1 √2𝜋𝑟[𝐾𝐼𝑐𝑜𝑠𝜃2+𝐾𝐼 2𝑠𝑖𝑛3𝜃 2𝑠𝑖𝑛𝜃2+𝐾𝐼𝐼 2𝑐𝑜𝑠3𝜃 2𝑠𝑖𝑛𝜃2] 𝜎𝑥𝑦=1 √2𝜋𝑟[𝐾𝐼 2𝑐𝑜𝑠3𝜃 2sin𝜃+𝐾𝐼𝐼𝑐𝑜𝑠𝜃2−𝐾𝐼𝐼 2𝑠𝑖𝑛3𝜃 2𝑠𝑖𝑛𝜃2] (3 – 29) Aplicando una rotación para expresar el tensor de tensiones en coordenadas polares, se obtiene la expresión de la tensión circunferencial 𝜎𝜃𝜃: 𝜎𝜃𝜃(𝜃=𝛽)=𝜎𝑦𝑦 𝑐𝑜𝑠2𝛽+𝜎𝑥𝑥 𝑠𝑖𝑛2𝛽−2𝜎𝑥𝑦sin𝛽cos𝛽 (3 – 30) Sustituyendo las expresiones del campo de tensiones en coordenadas cartesianas y realizando las convenientes agrupaciones, se puede expresar 𝜎𝜃𝜃 como se muestra a continuación: 𝜎𝜃𝜃=𝐾𝐼 √2𝜋𝑟 𝑐𝑜𝑠3𝜃2−3𝐾𝐼𝐼 2√2𝜋𝑟sin𝜃𝑐𝑜𝑠𝜃2 (3 – 31) Si ahora se introduce el ratio de carga mixta cot 𝛽∗=𝐾𝐼𝐼 𝐾𝐼, la expresión de la tensión circunferencial queda de la siguiente forma: 𝜎𝜃𝜃=𝐾𝐼 √2𝜋𝑟 [𝑐𝑜𝑠3𝜃2−3cot 𝛽∗ 2sin𝜃𝑐𝑜𝑠𝜃2] (3 – 32) Consideremos un caso de una placa infinitra sometida a un modo de carga mixto caracterizado por el ratio 𝛽∗ tal y como se muestra en la Figura 15. Se puede concluir que para el caso en que 𝛽∗=90 º la tensión circunferencial es máxima para 𝜃=0 º y la grieta se propagará en línea recta, correspondiendo a un modo I puro. Para 𝛽∗>0 la tensión circunferencial es máxima cuando 𝜃<0, lo que significa que la grieta se propaga en una dirección más cerca a la perpendicular con la carga. El ángulo de propagación de la grieta puede obtenerse de acuerdo a la siguiente expresión: 0=𝜕𝜎𝜃𝜃 𝜕𝜃 ⎹𝜃∗=𝐾𝐼 √2𝜋𝑟[−32𝑐𝑜𝑠2𝜃∗ 2sin𝜃∗ 2−3cot𝛽∗ 2𝑐𝑜𝑠3𝜃∗ 2+3cot𝛽∗𝑠𝑖𝑛2𝜃∗ 2cos𝜃∗ 2] (3 – 33) 26 Como se vió anteriormente, la función de desplazamiento que utiliza el método Singularity-Based y las funciones de singularidad 𝐹𝑗(𝑥) son las siguientes: 𝑢(𝑥)=𝑁𝐼(𝑥)∙𝑢𝐼+𝐻(𝑥)∙𝑁𝐼(𝑥)∙𝑎𝐼+𝑁𝐼(𝑥)∙∑𝐹𝑗(𝑥)∙𝑏𝐼𝑗 4 𝑗=1 𝐹𝑗(𝑥)=[√𝑟∙sin𝜃2,√𝑟∙cos𝜃2,√𝑟∙sin𝜃sin𝜃2,√𝑟∙sin𝜃cos𝜃2] (3 – 53) donde (𝑟,𝜃) representa las coordenadas de un sistema de coordenadas polares con el origen en el vértice de la grieta. Por otro lado, el conjunto de valores (𝜙,𝜓) permiten caracterizar las tensiones en el vértice de la grieta y el campo de desplazamientos en lugar de utilizar las coordenadas locales (𝑟,𝜃). La relación entre las variables anteriores se muestra a continuación: 𝑟=√𝜙2+𝜓2 𝜃=𝑡𝑎𝑛−1𝜓 𝜙 (3 – 54) Para comprender mejor que representan los valores 𝜓 y 𝜙 se muestran a continuación dos representaciones gráficas que lo clarifican. (a) (b) Figura 22: Cálculo y definición de los valores 𝜓 y 𝜙 [13]. 3.1.3.3 Análisis quasi-estático de crecimiento de grieta Con el método XFEM se asume que la simulación de crecimiento de grietas es quasi-estático. Para este análisis el método que se utiliza es el Phantom-Node, ya que ANSYS no permite realizar este tipo de análisis con el método Singularity-Based. Este método requiere llevar a cabo una serie de pasos de manera ordenada, los cuales se explican detalladamente a continuación. 27 1. Definición del modelo Este paso es común a los dos métodos. Entre las tareas involucradas destacan: a) Definición de los parámetros de enriquecimiento de la grieta, con el objetivo de definir la zona del modelo en la cual se localiza la grieta. Esta zona es mejorada utilizando nodos internos adicionales. Cabe destacar que es posible definir varias grietas en una misma zona de enriquecimiento. En este punto es necesario definir qué método queremos utilizar para el análisis de la grieta, Phantom-Node o Singularity-Based [13]. b) Definición del radio de mejora para tener en cuanto los efectos de las singularidades del vértice de grieta. Este punto evidentemente solo será utilizado en el método Singularity-Based. Este paso es necesario debido a que, por defecto, las funciones de singularidad solo se aplican al elemento donde se encuentra el vértice de la grieta, por lo que si queremos que los efectos de dicha singularidad se tengan en cuenta en los elementos de alrededor, deberemos definir este radio [13]. Figura 23: Definición del radio de enriquecimiento [13]. c) Definición de snap-tolerance. Este tercer paso solamente tiene aplicabilidad en el método Singularity-Based. Es necesario definir esta tolerancia debido a que, en un análisis típico con este método, es una buena práctica colocar el vértice de grieta en algún punto de la zona media del elemento [13]. Figura 24: Definición de snap-tolerance [13]. d) Definición de la grieta inicial. En este cuarto paso se definen los valores 𝜓 y 𝜙. En el caso de utilizar el método Phantom-Node solamente se definirá el valor 𝜙, ya que este método requiere intrínsecamente que el vértice de grieta se encuentre en el borde o cara del elemento [13]. 28 2. Definición del criterio de crecimiento de grieta ANSYS tiene implementado los dos siguientes criterios de crecimiento de grieta:  STTMAX (Criterio de máxima tensión circunferencial) Este criterio se basa en la evaluación del máximo valor de la tensión circunferencial 𝜎𝜃𝜃 cuando el vértice de grieta cambia de posición [13].  PSMAX (Criterio de tensión circunferencial basado en 𝜎𝑟𝜃=0) Este criterio en una alternativa al criterio STTMAX, que consiste en evaluar la tensión circunferencial en el punto donde 𝜎𝑟𝜃=0. Idealmente ambos criterios deben de arrojar los mismos resultados, pero debido a la discretización en elementos finitios, los resultados pueden diferir un poco [13]. Figura 25: Evaluación de los critgerios STTMAX y PSMAX entorno al vértice de grieta [13]. 3. Definición del decremento de tensiones en los segmentos de grietas creados Cuando la grieta comienza a abrirse, las tensiones cohesivas en los segmetos de grieta disminuyen gradualmente hasta cero conforme la deformación aumenta. Dicha disminución de las tensiones se modela de acuerdo a una ley cohesiva rígida lineal [13]. Figura 26: Ley cohesiva rígida lineal [13]. Cuando el desplazamiento perpendicular entre las dos superficies de los segmentos de grieta aumenta, las tensiones disminuyen linealmente desde el punto A hasta el punto B. Mientras la tensión no alcance el máximo valor especificado en el criterio de crecimiento, la grieta no crecerá. Cabe destacar también que el área encerrada bajo la curva corresponde a la tasa de liberación de energía crítica 𝐺𝐶, de manera que se puede relacionar dicho parámetro con el valor STTMAX tal que 𝐺𝐶=12𝛿𝑛∙𝑆𝑇𝑇𝑀𝐴𝑋 29 4. Especificación del comportamiento de la zona cohesiva en la grieta inicial Las grietas iniciales en el modelo pueden o no tener un comportamiento cohesivo. En caso de la grieta o grietas iniciales requieran dicho comportamiento, hay que especificarlo. Es importante tener en cuenta que XFEM solo admite un comportamiento de la zona cohesiva bilineal. Además, por defecto la interacción entre las superficies de la grieta se tiene en cuenta usando una formulación de sanción de contacto sencilla en la dirección normal. Las caras de la grieta se asumen que están en contacto libre de fricción [13]. 5. Evaluación del criterio de crecimiento de grieta Una vez definido el criterio, es necesario implementarlo para que se pueda llevar a cabo la evaluación de la tensión circunferencial mñaxima en el vértice de grieta [13]. 6. Cálculo del crecimiento de grieta Por último, se procede a calcular el criterio de crecimiento en la fase de solución, una vez que el análisis ha convergido [13]. 3.1.3.3.1 Hipótesis en la simulación XFEM de crecimiento de grieta A continuación, se exponen las hipótesis más relevantes a tener en cuenta en este tipo de análisis.  Se asume que el material presenta un comportamiento elástico lineal.  El análisis se considera que es quasi-estático, sin ser posible realizar otros tipos de análisis.  Debido a que los efectos de la singularidad del vértice de grieta no se tienen en cuenta en este tipo de análisis, el campo de tensiones y deformaciones en torno al vértice de grieta son aproximados. 3.2 Formulación y resolución del problema eléctrico con ANSYS Retomando el objetivo del proyecto, se trata de realizar un análisis numérico para ser capaces de detectar y monitorizar grietas de distintos tamaños en una placa de composites compuesto por epoxy y MWCNT (Multi Walled Carbon Nanotube). Este tipo de material presenta un comportamiento piezorresistivo, esto es, su resistividad eléctrica sufre cambios debido a la contribución de dos factores: deformación de la placa al estar sometida a unos esfuerzos mecánico y la presencia y crecimiento de grietas. Por tanto, en este punto se va a estudiar la formulación y como se resuelve el problema eléctrico con el software comercial ANSYS. 3.2.1 Cálculo de la respuesta piezorresistiva Modelar el comportamiento piezorresistivo de un material compuesto por una matriz epoxy y MWCNT requiere combinar las ecuaciones que gobiernan su comportamiento eléctrico con aquellas que describen el estado de deformaciones mecánicas. Para ello el tensor de resistividad eléctrica 𝜌𝑒𝑓𝑓 se calcla como la inversa del tensor conductividad eléctrica 𝜎𝑒𝑓𝑓. El tensor resistividad eléctrica puede escribirse en forma matricial de la siguiente manera: 𝜌𝑒𝑓𝑓=[𝜌1𝜌6𝜌5 𝜌6𝜌1𝜌4 𝜌5𝜌4𝜌3] (3 – 55) 30 Cuando el material no está sometido a deformación, MWCNTs están orientados de manera aleatoria a lo largo del material y el tensor de resistividad eléctrica toma la forma de una matriz escalar con los términos de la diagonal 𝜌0 [14]. Una vez que el material está sometido a deformación mecánica, las componentes de la matriz resistiva cambian de acuerdo a la siguiente expresión: [ 𝜌1 𝜌2 𝜌3 𝜌4 𝜌5 𝜌6 ] = [ 𝜌0 𝜌0 𝜌0 000 ] + [ ∆𝜌1 ∆𝜌2 ∆𝜌3 ∆𝜌4 ∆𝜌5 ∆𝜌6 ] o escrito más compacto, 𝝆𝑒𝑓𝑓=𝝆𝑒𝑓𝑓 0(𝑰+𝒓) (3 – 56) donde 𝜌𝑒𝑓𝑓 0 hace referencia al tensor de resistividad cuando el material no está deformado. El término 𝒓 define el tensor de cambio relativo en la resistividad y puede relacionarse con el tensor de deformación mecánica como 𝒓=𝜫∙ɛ , siendo 𝜫 la matriz piezorresistiva [14]. Si asumimos que la matriz piezorresistiva de un compuesto con CNT posee una simetría cristalina cúbica, se puede escribir: [ ∆𝜌1/𝜌0 ∆𝜌2/𝜌0 ∆𝜌3/𝜌0 ∆𝜌4/𝜌0 ∆𝜌5/𝜌0 ∆𝜌6/𝜌0 ] = [ 𝜆11 𝜆12 𝜆12 𝜆12 𝜆11 𝜆12 𝜆12 𝜆12 𝜆11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 𝜆44 0 0 0 𝜆44 0 0 0 𝜆44 ] ∙ [ ɛ1 ɛ2 ɛ3 2ɛ23 2ɛ13 2ɛ12 ] (3 – 57) Los términos 𝜆𝑖𝑗 son las componentes de la matriz piezorresistiva, por lo que 𝜆11 hace referencia al efecto piezorresistivo a lo largo de un eje cristalino principal para deformaciones aplicadas en ese eje (efecto piezorresistivo longitudinal), 𝜆12 se refiere al efecto piezorresistivo a lo largo de un eje cristalino principal para deformaciones aplicadas en un eje perpendicular (efecto piezorresistivo transversal) y 𝜆44 describe el efecto piezorresistivo en un campo eléctrico fuera del plano debido al cambio de la corriente inducida en el plano por la deformación cortante. Para calcular las componentes de 𝜫, se necesitan dos experimentos virtuales: ensayo de dilatación unilateral restringida y ensayo de distorsión [14]: [∆𝜌1/𝜌0 ∆𝜌2/𝜌0 ∆𝜌3/𝜌0]=[𝜆11 𝜆12 𝜆12 𝜆12 𝜆11 𝜆12 𝜆12 𝜆12 𝜆11]∙[00 ɛ3]=𝜫𝑑𝑖𝑙∙[00 ɛ3] [∆𝜌4/𝜌0 ∆𝜌5/𝜌0 ∆𝜌6/𝜌0]=[𝜆44 0 0 0 𝜆44 0 0 0 𝜆44]∙[00 2ɛ12]=𝜫𝑑𝑖𝑠∙[00 2ɛ12] (3 – 58) 3.2.2 Resolución del campo eléctrico Tal y como se ha comentado previamente, el objetivo es obtener los cambios que sufre la resistencia eléctrica de la placa cuando está sometida a una serie de esfuerzos mecánicos y existe una grieta que crece a lo largo del tiempo. Para realizar dicho cálculo, el procedimiento consiste en colocar estratégicamente unos electrodos en 31 los bordes de la placa, tal y como se muestra en el cuarto capítulo. Uno de esos electrodos suministrará un voltaje de valor conocido ∆𝜙 mientras que el otro electrodo estará conectado a tierra. Dichos electrodos hacen que se induzca un campo eléctrico en la placa, circulando una corriente de intensidad 𝐼 a través del material. De acuerdo a la ley de Ohm, existe una relación directa entre la intesidad y la resistencia, de forma que ∆𝜙=𝑅∙𝐼. Si el potencial ∆𝜙 se mantiene constante y teniendo en cuenta que la resistencia 𝑅 de la placa variará debido a la presencia de la grieta y a las deformaciones mecánicas del material, la intensidad I también variará. Dicho con otras palabras, existirá una densidad de corriente 𝐽 entre los bordes en los cuales han sido colocados los electrodos. El tensor de la conductividad eléctrica 𝜿 se puede obtener teniendo en cuenta que se define como la inversa de la resistividad, tal que 𝛋=𝝆−1 [14]. Por tanto, existe una relación directa entre la densidad de corriente 𝐽 y el potencial eléctrico ∆𝜙 tal que: 𝐉=−𝛋∙∇𝜙 (3 – 59) Por último, quedaría plantear la ecuación que resuelve el potencial. La distribución del potencial eléctrico 𝜙 en un dominio Ω caracterizado por el tensor de conductividad eléctrica 𝛋 puede calcularse atendiendo a la ecuación (3 - 60) como se muestra a continuación: ∇(𝛋∙∇𝜙)=0 (3 – 60) Finalmente, el problema consistirá en resolver el campo eléctrico, o la densidad de corriente 𝐽, y posteriormente proceder al cálculo de la corriente eléctrica 𝐼, ya que la intensidad de corriente no es mas que la suma de las densidad de la densidad de corriente que atraviesa una o varias secciones que tengan como vectores normales a la superficie 𝑑𝑠 [14] : 𝐼=∫𝐉∙𝐧 𝑆∙𝑑𝑠 (3 – 61) 3.2.3 Cálculo de la resistencia eléctrica en el panel dañado Teniendo en cuenta las formulaciones descritas en los dos puntos anteriores, si somos capaces de capturar la densidad de corriente que atraviesa cada uno de los elementos finitos tras la discretización, se podrá obtener la intesidad eléctrica que atraviesa el panel y por ende la resistencia eléctrica de la placa. Además, en el caso de un análisis quasi-estático de crecimiento de grieta, el problema radica en obtener el valor de la densidad de corriente en cada iteración del proceso, es decir, cada vez que la grieta aumente su tamaño. De esta forma, y teniendo en cuenta la ley de Ohm, se podrá obtener una evolución de la resistencia eléctrica de la placa en función del tamaño de la grieta, quedando resuelto el objetivo del problema de detectar y monitorizar el tamaño de la grieta en todo momento. Para realizar el cálculo anterior, es necesario disponer de las propiedades electromecánicas del material que se va a ensayar, tales como las componentes de la matriz piezorresistiva, la conductividad eléctrica del material, el módulo de elasticidad y el coeficiente de Poisson. Los valores de dichas propiedades, relativas a un material compuesto por fibra epoxy y MWCNT con un porcentaje de 0,5 𝑤𝑡% se muestran en la siguiente tabla. Componentes de la matriz piezorresistiva 𝜫 𝜎 (𝑆∙𝑚−1) 𝐸 (𝐺𝑃𝑎) 𝜈 𝜆11 𝜆12 𝜆44 6.84 7.99 1.19 1.22∙10−2 2.86 0.28 Tabla 1: Propiedades electromecánicas efectivas de MWCNT/epoxy [15] 32 3.3 Esquema de monitorización del crecimiento de grieta implementado con ANSYS En este último punto del tercer capítulo se van a incluir dos diagramas de flujo que permitan interpretar y dar a conocer al lector los pasos que se han seguido para resolver los problemas mecánicos y eléctricos que se abordan en este proyecto con el software comercial ANSYS. Figura 27: Esquema general de monitorización del crecimiento de grieta. En la Figura 27 se ha resumido de forma esquemática el proceso que se ha seguido para realizar la monitorización virtual del crecimiento de grieta con ANSYS. Tal y como puede apreciarse, el proceso comienza con un análisis mecánico, cuypo procedimiento se especifica detalladamente en la Figura 28. El análisis mecánico tiene como output un documento .txt con toda la información relativa a las coordenadas del vértice de grieta conforme esta va avanzando en el modelo creado. El documento .txt con con las coordenadas del vértice de grieta sirve como input para el análisis electroestático y piezorresistivo que se realiza posteriormente. Cabe destacar que ambos análisis están íntimamente relacionados. La diferencia entre ellos radica en que, en el análisis piezorresistivo se imponen condiciones de contorno eléctricas y mecánicas, mientras que en el análisis electroestático solo se utilizan condiciones de contorno eléctricas. Con el análisis electroestático y piezorresistivo se pretende obtener la densidad de corriente en cada uno de los pasos de carga, es decir, el análisis se realiza tantas veces como posiciones nuevas alcancen los vértices de grieta del problema mecánico, obteniendo de esta forma un valor de la densidad de corriente en cada uno de los avances de la grieta. Esta información se escribe en un documento .txt que será el output del análisis electroestático y piezorresistivo, que servirá a su vez como input para el posterior tratamiento de datos que se realizará en MATLAB para obtener finalmente los resultados que se mostrarán en el capítulo cuarto. Aunque en el esquema anterior no se ha representado, el documento con las coordenadas del vértice de grieta son también input para el postprocesado de datos. El procedimiento seguido para el análisis electroestático y piezorresistivo se detalla en la Figura 29. 33 Figura 28: Diagrama de flujo de la resolución del problema mecánico con ANSYS. 34 Figura 29: Diagrama de flujo de la resolución del problema electroestático y piezorresistivo con ANSYS. 35 4 ESTUDIOS NUMÉRICOS n este cuarto capítulo se va a abordar los estudios numéricos del proyecto, los cuales abarcan los procedimientos realizados tanto para la detección de grietas como para la monitorización del crecimiento de las mismas. En la primera parte, referente a la detección de grietas, se estudiará el caso de una placa cuadrada sometida a una carga de tracción con una grieta central ya existente. En esta parte se hará un estudio estático y se abordarán distintas configuraciones y casuísticas que se detallarán mas adelante. Por otra parte, en cuanto a la monitorización del crecimiento de grietas, se estudiarán paralelamente dos casos. El primero de ellos hace referencia al mismo problema estudiado en la primera parte de este capítulo. El segundo de ellos, abordará el caso de una placa rectangular sometida a flexión en dos puntos con una grieta inicial situada en el borde de la placa. En ambos casos, se estudiará el crecimiento de grietas predefinidas con distintas orientaciones y la relación existente entre dicho crecimiento y la resistividad de la placa. 4.1 Detección de grietas: análisis estático Tal y como se comentó previamente en la descripción de este capítulo, en este punto se va a abordar el problema de una placa cuadrada con una grieta predefinida con distintas orientaciones. A partir de esta configuración, se realizarán dos análisis:  Análisis mecánico Llamamos análisis mecánico al estudio relacionado con el cálculo de parámetros característicos de la Mecánica de Fractura Lineal. En este primer estudio se calculará analíticamente los valores de los SIF, tanto KI como KII, y se compararán los valores obtenidos con los valores teóricos.  Análisis eléctrico El análisis eléctrico hace referencia al cálculo de la resistencia eléctrica de la placa en función del ángulo de inclinación de la grieta respecto a la horizontal. En este punto se estudiará la variación de la resistencia eléctrica para distintas orientaciones con relación al caso de ausencia de grieta incial Los dos estudios anteriores se acompañarán de otros estudios secundarios:  Estudio de la influencia del tamaño de la grieta al cálculo de los SIFs.  Estudio de la influencia en la resistividad de la placa de considerar o no conductancia en la grieta.  Medición y detección de varias grietas en la placa. E 42 L/a=10 L/a=5 L/a=2.5 L/a=1.5 Figura 38: Influencia del tamaño de grieta sobre la distribución de potencial eléctrico para 𝛼=45 º 43 4.1.2.2 Influencia de la permitividad eléctrica de la grieta para varias orientaciones. Grieta permeable Los resultados mostrados anteriormente son relativos al caso de una grieta impermeable, es decir, aquella con una relación 𝜅𝑐/𝜅𝑚 nula. A continuación, se va a realizar un estudio de la influencia de las condiciones de permitividad eléctrica en la grieta (𝜅𝑐/𝜅𝑚) y la orientación de la misma (𝛼) en la variación de la resistencia eléctrica. Los resultados se muestran en la Figura 39, y son relativos al caso concreto de L/a=1.5. Figura 39: Influencia de las condiciones de permitividad eléctrica en la grieta (𝜅𝑐/𝜅𝑚) y la orientación (𝛼) en la variación de la resistencia eléctrica (𝑅). En la Figura 39 puede observarse que la permitividad eléctrica de la grieta es un parámetro que determina la resistencia eléctrica de la placa. Cuanto menor sea el valor de la relación 𝜅𝑐/𝜅𝑚, es decir, cuanto mas próximo este 𝜅𝑐 del valor nulo (caso de grieta eléctricamente impermeable), mayor es la resistencia eléctrica de la placa. A continuación, en las Figuras 40 y 41 se muestra la distribución del potencial eléctrica para el caso 𝐿/𝑎=1.5 y valores de orientación de la grieta 𝛼=0 º y 𝛼=45 º para distintas relaciones 𝜅𝑐/𝜅𝑚. Obsérvese en la Figura 41 la diferencia entre la distribuciñon de potencial eléctrico para 𝜅𝑐/𝜅𝑚=0.001 y 𝜅𝑐/𝜅𝑚=0.5, situaciones opuestass lo más similar posible al caso de grieta impermeable y permeable, respectivamente. 44 𝜅𝑐/𝜅𝑚=0.001 𝜅𝑐/𝜅𝑚=0.01 𝜅𝑐/𝜅𝑚=0.1 𝜅𝑐/𝜅𝑚=0.5 Figura 40: Influencia de la permitidad sobre la distribución de potencial eléctrico para 𝛼=0 º y 𝐿/𝑎=1.5. 45 𝜅𝑐/𝜅𝑚=0.001 𝜅𝑐/𝜅𝑚=0.01 𝜅𝑐/𝜅𝑚=0.1 𝜅𝑐/𝜅𝑚=0.5 Figura 41: Influencia de la permitidad sobre la distribucuión de potencial eléctrico para 𝛼=45 º y 𝐿/𝑎=1.5. 4.1.2.3 Medición y detección de varias grietas En este último problema se va a utilizar la misma geometría para la placa e iguales condiciones de contorno eléctrica que en los casos anteriores. Concretamente se van a modelar dos grietas horizontales impermables (𝜅𝑐=0) cuyos tamaños serán la mitad de la grieta original, es decir, si el semiancho de grieta inicial lo llamamos 𝑎1 y al semiancho de las grieta nueva como 𝑎2, se obtiene la relación 𝑎1=2𝑎2. Con este estudio se pretende analizar si los resultados obtenidos son iguales al caso de que existiera una grieta equivalente de tamaño 𝑎1. Para la localización de dichas grietas se utilizarán dos parámetros: 𝛿𝑥 y 𝛿𝑦. En el caso de 𝛿𝑥, representa la distancia entre los vértices de grieta más próximos de ambas grietas. Por otro lado 𝛿𝑦 representa la distancia entre los ejes de las dos grietas. Con la descripción anterior, parece razonable definir dos sistemas de coordenadas locales para cada una de las grietas. El sistema x’ - y’ (grieta 1) tendrá su origen en el punto (−𝛿𝑥/2 ,𝛿𝑦/2). Por otro lado, el origen del sistema x’’ - y’’ se encontrará en el punto (𝛿𝑥/2 ,−𝛿𝑦/2). Una vez descrito los parámetros a utilizar y la geometría de la placa, se procederá a utilizar diferentes relaciones 𝛿𝑦/𝑎2 para simular distintos casos en los que las grietas se encuentren a mayor o menor distancia entre sí. 46 Concretamente se simularán los casos 𝛿𝑦/𝑎2=[0,0.1,0.5,1]. Como resulta lógico, ambas grietas estáran alineadas en dirección x cuando 𝛿𝑦/𝑎2=0. Para comprender un poco mejor la filosofía de este problema, se muestra en la Figura 42 un esquema. Una vez comprendido el problema, se procede a realizar el análisis eléctrico. Los resultados del mismo se muestran a continuación en la Figura 43. Figura 42: Esquema problema eléctrico para la medición y detección de dos grietas equivalentes a la original. Figura 43: Influencia de la separación 𝛿𝑥 y 𝛿𝑦 en la variación de la resistencia eléctrica del problema con grieta equivalente. De la figura anterior puede deducirse que los menores valores de resistencia eléctrica de la placa se alcanzan cuando las grietas están separadas en dirección x un valor negativo de 0.01 m, es decir, cuando las grietas están 47 solapadas. Por otra parte, los valores más altos de resistencia eléctrica se alcanzan para valores de 𝛿𝑥=0. Por último cabe destacar el resultado obtenido para 𝛿𝑥=0 y 𝛿𝑦/𝑎2=0. Estos valores representan la situación en la que las dos grietas se encuentran alineadas según su eje mayor y en la que los vértices derecho e izquierdo de las grietas 1 y 2 coinciden, lo que es equivalente al caso de una sola grieta cuya tamaño sea 2𝑎1. Se puede comprobar la coherencia del resultado, ya que la relación 𝑅2/𝑅1 toma un valor prácticamente igual a la unidad. A continuación, las Figuras 44-46 muestran la distribución del potencial eléctrico para las distintas relaciones 𝛿𝑦/𝑎2 y para valores 𝛿𝑥/𝑎2=[−0.5 0 0.5], representando físicamente estos valores los casos de grieta “enmascarada” y grietas separadas en la dirección longitudinal de la grieta. a) b) Figura 44: Potencial eléctrico para 𝛿𝑦/𝑎2=0 (a) y 𝛿𝑦/𝑎2=0.05 (b). 48 a) b) Figura 45: Potencial eléctrico para 𝛿𝑦/𝑎2=0.1 (a) y 𝛿𝑦/𝑎2=0.5 (b). 49 a) b) Figura 46: Potencial eléctrico para 𝛿𝑦/𝑎2=1 (a) y 𝛿𝑦/𝑎2=2 (b). 50 4.1.3 Análisis piezorresistivo El modelo que se ha simulado corresponde a un material compuesto por una matriz epoxy y MWCNTs como refuerzos. Como se comentó en el capítulo tercero, este material presenta un comportamiento piezorresistivo, de forma que la resistencia eléctrica de la placa experimenta una variación al deformarse. Hasta ahora, solo se ha realizado un análisis electroestático, en el cual se ha analizado la variación que experiementa la resistencia eléctrica de la placa cuando existe una grieta inicial, independientemente de la configuración y características de dicha grieta. En este punto se va a realizar un análisis piezorresistivo para ver si existe diferencia entre realizar un análisis puramente eléctrico y un análisis introduciendo dicho comportamiento piezorresistivo del material. Para ello el análisis se ha realizado considerando el caso de una grieta impermeable, mostrándose los resultados obtenidos en la Figura 47. Figura 47: Comparación entre análisis piezorresistivo y electroestático. De la Figura 47 puede concluirse que, para un caso de grieta estática como el analizado, a penas existen diferencias entre los resultados obtenidos al considerar el efecto piezorresistivo o no, quedando justificado que los resultados obtenidos en el análisis electroestático son válidos para el caso de grieta estática. 4.2 Monitorización virtual del crecimiento de grietas En este segundo punto del capítulo cuatro, se abordará el tema de la monitorización virtual de grietas dinámicas, es decir, grietas que van creciendo progresivamente al estar sometida la placa en cuestión a unos esfuerzos determinados. Al igual que se hizo en el punto anterior, para realizar el estudio correctamente se va a llevar a cabo un análisis mecánico y un análisis eléctrico.  Análisis mecánico El análisis mecánico tiene por objetivo definir la geometría del problema, la configuración de cargas y condiciones de contorno, la definición de una grieta inicial y finalmente analizar el crecimiento de la misma. A través de este análisis seremos capaces de conocer el estado de la grieta y su tamaño en cada momento. 51  Análisis electroestático Una vez realizado el análisis mecánico, se procederá a diseñar el problema electroestático A través de este problema se pretende analizar la variación de la resistividad de la placa conforme la grieta va alcanzado una nueva posición a lo largo de su proceso de crecimiento. De esta forma seremos capaces de evaluar la resistencia eléctrica de la placa para cada tamaño de grieta conforme esta va creciendo a lo largo del plano de la placa. Para abordar este punto se han proyectado dos problemas distintos. Uno de ellos es totalmente indéntico al problema abordado en el caso de grieta estática, es decir, el problema estará compuesto por una placa cuadrada sometida a tracción. La única diferencia es que en este caso no se impondrán esfuerzos en el borde superior de la placa, sino que se impondrá un desplazamiento 𝑢𝑦. El segundo problema consiste en una placa rectangular apoyada en dos puntos sometida a flexión. Para ello el problema se ha modelado aplicando un desplazamiento 𝑢𝑦 en dirección y negativa aplicado en el punto central del borde superior de la placa. Tal y como se comentó previamente, el objetivo es monitorizar el crecimiento de la grieta y ser capaces de analizar la influencia que tiene dicho crecimiento en la variación de la resistencia eléctrica de la placa. Ambos estudios, tanto el mecánico como el eléctrico, se ejecutarán para los dos problemas mencionados. Por último, cabe destacar que los análisis se realizaran para distintas configuraciones de grieta, es decir, la grieta se definirá inicialmente con distintos ángulos de orientación y posteriormente la grieta crecerá en la dirección que tenga que hacerlo por la naturaleza del problema. 4.2.1 Placa con grieta interna 4.2.1.1 Análisis mecánico Para desarrollar como se ha llevado a cabo el análisis y la obtención de los resultados de la parte mecánica, se hará una distinción entre los dos problemas comentados previamente, ya que cada uno de ellos presenta una configuración distinta. En primer lugar, se va a describir la geometría de la placa, configuración de cargas y condiciones de contorno y la definición de la grieta inicial. Posteriormente se dará paso a la justificación del mallado y elementos utilizados. Por último, se finalizará con el análisis de los resultados obtenidos. Respecto a la geometría, se trata de una placa cuadrada (2L x 2L x t ; L=100mm y t=2 mm), la cuál estará sometida a un desplazamiento uniforme aplicado en el borde superior de la placa (𝑢𝑦=0,0025 𝑚). Además, se impondrán unas condiciones de contorno en el borde inferior de la placa, de manera que los desplazamientos en dirección y estáran impedidos en todo el borde inferior y el desplazamiento en dirección x sólo estará impedido en el punto central del borde inferior (ver Figura 48). La placa presenta una grieta recta centrada con una orientación que irá variando en cada una de las configuraciones. Los distintos ángulos adoptados para el análisis son 0 º, 15 º, 30 º, 45 º, 60 º y 75 º. En la Figura 48 se muestra un esquema que representa la configuración adoptada para el análisis mecánico. Para el mallado de la placa se han utilizado elementos tipo PLANE182. Con el obetivo de facilitar la definición de la grieta inicial y obtener unos resultados coherentes de crecimiento, la placa se ha dividido en siete áreas distintas tal y como puede verse en la Figura 49. 58 (a) (b) (c) (d) Figura 58: Evolución de distribución del potencial eléctrico placa a tracción para 𝛼=45 º. Figura 59: Análisis de la influencia del tamaño de grieta en la variación de la resistencia eléctrica en el caso de placa sometida a tracción para 𝛼=[0,15,30,45,60,75] º. 59 4.2.1.3 Análisis piezorresistivo Al igual que se hizo en el análisis estático de grietas, se va realizar un análisis piezorresistivo para ver si existen diferencias al considerar el efecto piezorresistivo del material o no. La Figura 60 muestra los resultados obtenidos de este análisis, pudiéndose establecer una comparación entre los resultados del análisis electroestático ejecutado previamente y el análisis piezorresistivo. Figura 60: Comparación análisis piezorresistivo y electroestático para placa con grieta interna. Tal y como puede observarse en la figura anterior, al introducir el efecto piezorresistivo la resistencia eléctrica de la placa aumenta para todas las inclinaciones de grieta analizadas. Sin embargo, hay que destacar que la influencia de incluir o no el efecto piezorresistivo en el análisis es mínimo. El error de considar o no el efecto resistivo se encuentra en un rango del 4 % y 6 %. Por tanto, puede llegarse a la conclusión de que los análisis previamente realizados sin considerar el efecto piezorresistivo son totalmente válidos. 4.2.2 Viga a flexión en tres puntos con defecto de borde 4.2.2.1 Análisis mecánico Al igual que se hizo con el problema anterior, primero va a realizarse una descripción de la geometría de la placa, configuración de cargas y condiciones de contorno y de la grieta inicial. Luego se justificará la malla y elementos utilizados. Finalmente se hará un análisis de los resultados obtenidos. Se trata de una placa plana rectangular (𝐿 𝑥 𝑊 𝑥 𝑡 ; 𝐿=100𝑚𝑚,𝑊=30 𝑚𝑚 𝑦 𝑡=2 𝑚𝑚), la cuál estará sometida a un desplazamiento aplicado en el punto central del borde superior de la placa (𝑢𝑦=−0,003 𝑚). Además, se impondrán unas condiciones de contorno en las esquinas inferiores derecha e izquierda de la placa, de forma que los desplazamientos en dirección y quedan impedidos (ver Figura 61). La placa presenta una grieta recta localizada en el centro del borde inferior de la placa, con una orientación que irá variando en cada una de las configuraciones. Los distintos ángulos seleccionados para ejecutar el análisis son 0 º, 15 º, 30 º, 45 º y 60 º. Dicha grieta cuenta con un sistema de coordenadas local que no coincide con el sistema de coordenadas global. El sistema de coordenadas global tiene su origen en el borde inferior izquierdo de la placa. Por otra parte, el sistema de coordenadas global está orientado según la grieta, de forma que el eje x’ tiene la dirección del eje de 60 la grieta y el eje y’ es perpendicular a la misma. Al estar ligado dicho sistema de coordenadas a la grieta, el eje x’ formará un ángulo 𝛼+ 𝜋/2 respecto al eje x global, siendo 𝛼 el ángulo de orientación de la grieta respecto a la vertical. En la Figura 61 se muestra un esquema que representa la configuración adoptada para el análisis mecánico. Para el mallado de la placa se han utilizado elementos tipo PLANE182. Al igual que en el ejemplo anterior, resultó conveniente estudiar la distribución de áreas para posteriormente mallar cada una de ellas y facilitar el cálculo y análisis de resultados En la Figura 62 se puede apreciar la distribución elegida. Figura 61 Esquema problema mecánico de placa con grieta inclinada sometida a flexión. Figura 62: Definición de áreas en problema mecánico de placa sometida a flexión para 𝛼=45 º. A continuación, se justifica la distribución y definición de las áreas mostradas en la Figura 62 tal y como se hizo en el caso anterior:  Las áreas A2 y A3 se mallarán mediante “free meshing”. En estas áreas no es vital que los elementos sean regulares, ya que la grieta no atravesará dichas zonas. Sin embargo, tal y como se verá más adelante, los elementos adjuntos a los bordes derecho e izquierdo de la placa son regulares con el objetivo de obtener con exactitud unos resultados coherentes posteriormente en el problema eléctrico.  En las áreas A5 y A6 se utilizará también “mapped meshing” ya que son regulares. Aún asi, estas áreas podrían haber sido malladas mediante “free meshing” ya que la grieta tampoco pasará por ellas.  Las áreas A7 y A8 estarán malladas utilizando “free meshing” al ser irregulares. A priori la grieta no atravesará dichas zonas, aunque en caso de que lo hicieran sería un recorrido muy corto que no afectaría a los cálculos y obtención de resultados.  En el área 4 será por donde la grieta avance en su inmensa mayoría. Sin embargo, no ha sido posible utilizar “mapped meshing” al ser un área irregular. Aún así los elementos serán lo mas regulares posible y tendrán forma cudrilátera, lo que facilitará el post-procesado de los resultados. 61  Por último, el área A1 es el más importante de todas ya que en ella se definirá la grieta inicial. Esta área se mallará mediante “mapped meshing”, con el objetivo de definir los elementos con la mayor exactitud posible y tener controlado la distribución de los mismos. Dicha área estará definida de acuerdo al sistema de coordenadas local x’ – y’ conforme a la orientación de la grieta (𝛼). A continuación, se muestra el mallado utilizado en la placa para el caso de 𝛼=45 º. El proceso es análogo para el resto de orientaciones. En la Figura 63 puede apreciarse el avance de la grieta por cada una de las zonas, quedando de esta forma justicado el mallado utilizado. Figura 63: Mallado de placa sometida a flexión para 𝛼=45 º. Una vez expuesto todo lo relativo a la geometría y mallado del problema, se dará paso al análisis de los resultados. La Figura 64 muestra los isocontornos de desplazamiento (𝑢𝑦) y la tensión equivalente de Von Misses (𝜎𝑒𝑞) para el caso de 𝛼=0 º una vez que la grieta ha alcanzado su tamaño final. Tal y como se hizo en el problema de la placa sometida a tracción, se va a mostrar una evolución de los isocontornos de desplazamiento (𝑢𝑦) para diferentes tamaño de grieta. Concretamente se mostarán las relaciones 𝑎/𝑎0=[1,2,3,4], donde 𝑎 es el valor de la grieta en un instante 𝑡 y 𝑎0 es el tamaño de la grieta inicial. Por otro lado, se van a mostrar los resultados obtenidos en cuanto a isocontornos de desplazamiento (𝑢𝑦) y tensión equivalente de Von Misses (𝜎𝑒𝑞) se refiere, para el caso de 𝛼=45 º (ver Figura 66). En la Figura 67 se puede ver la evolución de los isocontornos de desplazamiento para distintos tamaños de grieta 𝑎/𝑎0= [2,3,4]. (a) (b) Figura 64: (a) Isocontornos de desplazamientos 𝑢𝑦 y (b) tensión equivalente Von Mises 𝜎𝑒𝑞 para 𝛼=0 º. 62 (a) (b) (c) (d) Figura 65: Isocontornos de desplazamientos 𝑢𝑦 y 𝛼=0 º para (a) 𝑎/𝑎0=1 , (b) 𝑎/𝑎0=2 , (c) 𝑎/𝑎0=3 y (d) 𝑎/𝑎0=4 . (a) (b) Figura 66: (a) Isocontornos de desplazamientos 𝑢𝑦 y (b) tensión equivalente Von Mises 𝜎𝑒𝑞 para 𝛼=45 º. 63 (a) (b) (c) Figura 67: Isocontornos de desplazamientos 𝑢𝑦 y 𝛼=45 º para (a) 𝑎/𝑎0=2, (b) 𝑎/𝑎0=3 y (c) 𝑎/𝑎0=4, 4.2.2.2 Análisis electroestático En este punto abordaremos el análisis electroestático para el caso de una placa apoyada en dos puntos sometida a flexión, cuya geometría es la misma que la descrita en el apartado de análisis mecánico correspondiente. En este caso, se coloca en la placa un electrodo que suministra 10 V en el borde izquierdo de la placa y el borde derecho se conectará a tierra. Para el análisis eléctrico se ha utilizado la misma malla que en el problema mecánico. Cabe destacar que los bordes superior e inferior de la placa se encuentran aislados eléctricamente, es decir, la densidad del flujo eléctrico es nula. A continuación, en la Figura 68 se muestra una representación del problema eléctrico planteado. Una vez descrito el modelo del problema eléctrico aplicado a la placa apoyada en dos puntos sometida a flexión, se dará paso al análisis de resultados. En la Figura 69 se muestra la distribución del potencial eléctrico concretamente para una orientación de la grieta de 𝛼=0 º y 𝛼=45 º. 64 Figura 68: Esquema problema eléctrico de placa con grieta inclinada genérica sometida a flexión. (a) (b) Figura 69: Distribución de potencial eléctrico para el problema de placa sometida a flexión. Por otro lado, en las Figuras 70 y 71 se muestran los resultados relativos a la evolución de la distribución del potencial eléctrico para orientaciones de grieta 𝛼=0 º y 𝛼=45 º. Para ambas orientaciones se mostará dicha evolución para un tamaño de grieta 𝑎/𝑎0= [2,3,4]. 65 (a) (b) (c) Figura 70: Evolución de distribución del potencial eléctrico placa a flexión para 𝛼=0 º. 66 (a) (b) (c) Figura 71: Evolución de distribución del potencial eléctrico placa a flexión para 𝛼=45 º. Por último, se representará la evolución de cambios de la resistencia eléctrica de la placa (𝑅/𝑅0) en función del tamaño de longitud de grieta (𝑎/𝑎0). Dicha evolución se mostrará para orientaciones de grieta tales que 𝛼= [0,15,30,45,60] º (véase Figura 72). 67 (a) (b) Figura 72: Análisis de la influencia del tamaño de grieta en la variación de la resistencia eléctrica en el caso de viga sometida a flexión en tres puntos para (a) 𝛼=0 º y (b) distintos valores de 𝛼. 4.2.2.3 Análisis de la variación de la resistencia eléctrica frente a la flecha en el caso de viga apoyada sometida a flexión en tres puntos El último análisis llevado a cabo en el proyecto consiste en ver como varía la resistencia eléctrica de la viga sometida a flexión en tres puntos conforme la flecha aumenta. Para este análisis se ha utilizado una orientación de grieta 𝛼=0 º y considerando todos y cada uno de los pasos de carga que utiliza ANSYS para resolver. De esta forma, podemos analizar si la resistencia eléctrica varía linealmente conforme la flecha aumenta o si por el contrario en algunos momentos permanece constante ya que, a pesar de existir un desplazamiento vertical de la viga, la grieta no crece. Figura 73: Variación de la resistencia eléctrica en función del desplazamiento vertical de la viga. 74 Monitorización virtual del crecimiento de grietas Código ANSYS: Mecánico grieta interna.txt Este código modela el caso de una grieta inicial con una inclinación 𝜶 que ha de ser introducida por el usuario. El programa calcula las coordenadas de los vértices de grieta en cada paso de carga y filtra aquellos en los cuales la grieta crece. De esta forma, se almacenan en un archivo .txt que se genera automática las coordenadas X e Y de los dos vértices de grieta. Código ANSYS: Electroestático grieta interna.txt Este código utiliza como input el archivo .txt generado anteriormente con los vértices de grieta y realiza un análisis electroestático para calcular la densidad de corriente que atraviesa la placa. Para ello, el programa ha sido diseñado de forma que los elementos en los cuales se encuentre la grieta serán anulados, simulando que en esos elementos no existe material. Generará un archivo. txt con los valores de la densidad de corriente total en cada paso de carga. Código ANSYS: Piezorresistivo grieta interna.txt Este código es similar al anterior, con la diferencia de que se le atribuyen al material las propiedades eléctricas, mecánicas y piezorresistivas para simular el efecto resistivo del mismo. Código ANSYS: Mecánico viga flexión.txt Este código modela el caso de una grieta inicial con una inclinación 𝜶 que ha de ser introducida por el usuario. El programa calcula las coordenadas de los vértices de grieta en cada paso de carga y las filtra al igual que el código de análisis mecánico de grieta interna. Este código genera un archivo .txt con los valores de la densidad de corriente en cada paso de carga y otro archivo con la flecha que experimenta la viga en cada paso de carga. Código ANSYS: Electroestático viga flexión.txt El código ha sido diseñado para realizar un análisis electroestático al igual que en el caso de grieta interna, con los mismos inputs y outputs. Código ANSYS: Piezorresistivo viga flexión.txt Igual que el código anterior pero teniendo en cuenta el efecto piezorresistivo. 75 Anexo B. Códigos de MATLAB Detección de grietas Código MATLAB: SIFs.m Este código recibe como entrada el valor de los factores de intensidad de tensiones obtenidos en las simulaciones de ANSYS para diferentes ángulos de orientación de las grietas y calcula la convergencia de dichos valores a la solución analítica. Código MATLAB: Estatico_Variación_Resistencia.m Este código utiliza como input el valor del flujo eléctrico obtenido en ANSYS para diferentes relaciones L/a y calcula la variación que experimenta la resistencia eléctrica de la placa para las diferentes relaciones L/a y orientaciones de la grieta. Código MATLAB: Estatico_Variación_Resistencia_Piezorresistivo.m Este código es similar al anterior, con la diferencia de que calcula la resistencia eléctrica de la placa cuando se tiene en cuentra el efecto piezorresistivo, mostrando gráficamente la diferencia que existe entre ejecutar un análisis puramente electroestático y un análisis piezorresistivo. Código MATLAB: Estatico_Grieta_Permeable.m Este código muestra la diferencia que existe en la variación de la resistencia eléctrica de la placa para diferentes orientaciones de la grieta y diferentes relaciones 𝜿𝒄/𝜿𝒎. Código MATLAB: Estatico_Varias_Grietas.m Este código calcula la variación de la resistencia eléctrica que experimenta la placa cuando existe una grieta de tamaño a y dos griets cuyo tamaño es la mitad de la primera. Así mismo se evalúa dicha diferencia para distintas posiciones y ubicaciones de las grietas, variando la distancia en dirección X entre ambas grietas y la distancia en dirección Y entre ellas. 76 Monitorización virtual del crecimiento de grietas Código MATLAB: Grieta_Interna_Crecimiento.m Este código recibe como input el flujo de densidad de corriente para distintas orientaciones de grieta y la posición del vértice de grieta en cada paso de carga que la grieta crece. Con estos datos, calcula la variación de la resistencia eléctrica conforme la grieta va creciendo. Código MATLAB: Viga_Flexión_Crecimiento.m El código es similar al anterior pero diseñado para el problema de viga a flexión. La diferencia fundamental es que en el código anterior existían dos vértices de grieta y en este problema solo uno, por lo que el código debe ser adaptado para capturar este crecimiento de grieta. Código MATLAB: Grieta_Interna_Crecimiento_Piezorresistivo.m Este código permite realizar una comparación entre la variación que experimenta la resistencia eléctrica de la placa en el caso de análisis electroestático y análisis piezorresistivo, para orientaciones de grieta 𝜶= 𝟎 º y 𝜶=𝟒𝟓 º. Código MATLAB: Viga_Flexión_Crecimiento_Piezorresistivo.m Este código realiza la misma comparación que el código anterior pero para el caso de viga a flexión. Código MATLAB: Viga_Flexión_Crecimiento_Flecha.m Este código recibe como input el flujo de densidad de corriente, la posición del vértice de grieta durante el proceso de crecimiento y el desplazamiento en dirección Y (flecha) en cada paso de carga. De esta manera, se representa la variación que experimenta la resistencia eléctrica en función del desplazamiento experimentado.