Modelo de predicción del fallo en materiales compuestos bajo cargas transversales a las fibras. Interacción entre los despegues de interfase fibra-matriz y daño en la matriz

28 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 Diseño Avanzado en Ingeniería Mecánica Modelo de predicción del fallo en materiales compuestos bajo cargas transversales a las fibras. Interacción entre los despegues de interfase fibra-matriz y daño en la matriz. Autor: José Miguel Pimentel Medina Tutor: Luis Aristides Távara Mendoza Dpto. de Mecánica de Medios Continuos y Teoría de Estructuras Grupo de Elasticidad y Resistencia de Materiales Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Proyecto Fin de Máster Máster en Diseño Avanzado en Ingeniería Mecánica Modelo de predicción del fallo en materiales compuestos bajo cargas transversales a las fibras. Interacción entre los despegues de interfase fibra-matriz y daño en la matriz. Autor: José Miguel Pimentel Medina Tutor: Luis Aristides Távara Mendoza Dpto. de Mecánica de Medios Continuos y Teoría de Estructuras Grupo de Elasticidad y Resistencia de Materiales Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Proyecto Fin de Máster: Modelo de predicción del fallo en materiales compuestos bajo cargas transversales a las fibras. Interacción entre los despegues de interfase fibra-matriz y daño en la matriz. Autor: José Miguel Pimentel Medina Tutor: Luis Aristides Távara Mendoza 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 Pr imero que todo, agradezco el apoyo incondicional de mi familia por apoyarme durante todo este tiempo en mi decisión de continuar estudios superiores en España; su dedicación y el esfuerzo familiar serán honrados. Al profesor Luis Távara, por dejarme formar parte de este proyecto y por su orientación y consejos brindados durante todo el periodo académico. Extiendo mi agradecimiento para Alejandro Estefani, por su orientación técnica durante el seminario de elementos finitos del departamento. A todo el Grupo de Elasticidad y Resistencia de Materiales en la escuela, en particular al Profesor Vladislav Mantič: su rigor ha sido ejemplo para impulsar mi interés en los aspectos teóricos y numéricos de la mecánica de medios continuos. A Natalia García y Antonio, por los buenos momentos compartidos, su hospitalidad y amabilidad conmigo siendo compañeros de piso en Sevilla. José Miguel Pimentel Sevilla, 2022 I Resumen El uso extendido de los materiales compuestos en estructuras aeronáuticas primarias ha llevado a la necesidad de entender a mayor profundidad los mecanismos de fallo en dichos materiales heterogéneos. En esta tesis, se considera tanto el fallo inter-fibra como el de la matriz en la escala micromecánica para las láminas de 90 ◦ bajo carga transversal. Para alcanzar este objetivo, el comportamiento de la interfase fibra-matriz y el de la matriz se han implementado en una subrutina de material (UMAT) del programa comercial de elementos finitos Abaqus; las inestabilidades post-pico se capturan usando el método secuencial lineal programado en Python. Se han realizado diversos estudios en un modelo unifibra considerando la influencia de la topología del mallado de la matriz, el criterio de fallo de la matriz, niveles de densidad de energía críticos y estados de tensión biaxiales sobre el camino de propagación de la grieta; los resultados indican una mayor energía disipada para valores en incremento de la densidad de energía crítica mientras que la tensión crítica σ∞ que inicia el mecanismo de fallo se mantiene invariante. Los criterios de fallo para la matriz basados en tensiones muestran, en la mayoría de casos que el mecanismo de fallo empieza en la matriz alrededor de la circunferencia de la fibra, con la aparición de bifurcaciones espurias de grieta. Los estudios realizados en modelos multifibra para distintas configuraciones de laminados [ 0,90n ] s , demuestran la capacidad de los criterios de fallo en producir una grieta transversal que atraviesa todo el espesor de la lámina antes de empezar el proceso de delaminación. El modelo de predicción es capaz de reproducir la bifurcación de la grieta transversal principal en los modelos integrados multiescala, cuando el frente de grieta en la matriz alcanza el borde inferior de la interfase de una fibra. Palabras clave: Propagación de grieta; Materiales compuestos; Análisis por elementos finitos; Despegues de interfase; Fallo de matriz; Análisis secuencial lineal. III 1 Introducción El uso creciente de los materiales compuestos en diferentes sectores industriales como el aeroespacial, automotriz, de las energías renovables, etc., ha llevado a un estudio má profundo del daño y los mecanismos de fallo que afectan a estos. Este proyecto se centrará en el estudio del fallo de la interfase fibra-matriz, y en la implementación de un modelo de interfase en el código de elementos finitos ABAQUS que permita el estudio del inicio y la propagación de las roturas en las interfases fibra-matriz, junto con el modelado del daño de la matriz en modelos unifibra y multifibra en la escala micromecánica. El uso de los materiales compuestos reforzados con fibras de vidrio o carbono trae actualmente problemas que se derivan del grado de conocimiento que se tiene del fallo de estos materiales, en particular, en relación al que se tiene de los materiales metálicos. Por un lado está el menor tiempo en vigencia desde su uso comercial y por otro lado hay que considerar el número y complejidad de los mecanismos de daño en los materiales compuestos, así como la interacción de los mismos, todo esto derivado de la naturaleza heterogénea de estos materiales. Desde el punto de vista mecánico, son varios los mecanismos de daño que pueden afectar el comportamiento mecánico (en servicio) de dichos materiales. Los escenarios característicos de daño son generalmente categorizados como: 1. Eventos de daño intralaminar como e.g., fractura de la fibra, fallo de la matriz, decohesión fibra-matriz [14, 34, 48], entre otros. 2. Fallo interlaminar, i.e., delaminación o despegue [ 68 ] entre juntas estructurales unidas con adhesivo en aplicaciones ingenieriles [50, 10]. Durante más de una década, el Grupo de Elasticidad y Resistencia de Materiales (GERM) de la Universidad de Sevilla ha desarrollado un modelo de interfase denominado "Linear Elastic Brittle Interface Model" (LEBIM) [ 35 , 45 ], que ha permitido modelar p.ej., la nucleación y el crecimiento de grietas de interfase en inclusiones cilíndricas embebidas en una matriz finita, empleando un modelo numérico de elementos de contorno (BEM) implementado por L. Távara y colaboradores [ 63 ], además de poder modelizar la nucleación y crecimiento de los despegues de interfase en compuestos unidireccionales reforzados con fibra bajo cargas transversales bi-axiales [ 62 ], este último estudio también usando modelos numéricos de elementos de contorno. D. Castillo [ 11 ] implementó el modelo de interfase LEBIM en el programa de elementos finitos Abaqus, logrando modelizar los despegues en los estudios numéricos, aunque sin poder capturar el comportamiento de "snap-through" (i.e., no se capturan las inestabilidades post-pico después del primer despegue) en los modelos micromecánicos unifibra y multifibra. Luego, E. Paloma [ 39 ] implementa el método de solución secuencial lineal "Sequential Linear Analysis" [ 64 , 38 ] en el código de elementos finitos de Abaqus, permitiendo entonces capturar las inestabilidades post-pico de la interfase frágil en los sistemas fibra-matriz propios de una lámina de material compuesto. En la actual literatura se han desarrollado ciertos enfoques para abordar el problema de la propagación de grietas en materiales frágiles (ver p.ej., Wingender y Balzani [ 71 ], Pandolfi y Ortiz [ 40 ]), sin embargo en dichos enfoques se requiere la existencia de una fisura inicial cuya dimensión tiende a ser del orden de las dimensiones estructurales del componente, dejándose fuera la escala micromecánica del daño y procediendo entonces a degradar las propiedades resistentes en ciertos elementos para promover el inicio de grieta [71]. 1 2Capítulo 1. Introducción El trabajo propuesto aquí pretende avanzar en una línea de investigación de la predicción del fallo en la escala micromecánica de láminas de compuesto bajo carga transversal, implementando un modelo de interfase débil (LEBIM), en conjunto con diversos criterios de fallo en la matriz y el método de solución secuencial lineal (SLA), para elucidar cuestiones relativas al comienzo y propagación de grietas, junto al valor de las tensiones y deformaciones (σ2yε2en ejes de la lámina) que inician el proceso de daño. Este nuevo aporte permitirá: •Modelizar el despegue de la interfase fibra-matriz junto con el daño de la matriz. • Emplear distintos criterios de fallo para la matriz y analizar cuáles son capaces de reproducir con mayor fidelidad la topología de fractura experimental en las láminas de 90 grados bajo cargas transversales. • Conocer el valor de tensiones y deformaciones que inician el mecanismo de fallo en modelos multifibras y en modelos integrados multiescala de laminados [ 0,90n ] s en presencia de una grieta transversal inicial. 1.1 Motivación El conocimiento actual acerca del daño y los mecanismos de fallo en los materiales compuestos no está lo suficientemente avanzado como para desarrollar criterios de base física para varios tipos de fallo. Actualmente, existe una demanda en incremento para desarrollar herramientas de modelado computacional que permitan una exacta predicción de los distintos mecanismos de fallo que ocurren en los compuestos reforzados por fibras. El uso de modelos alternativos basados en supuestos distintos a los que ofrece la Mecánica de la Fractura Lineal Elástica (LEFM) [ 65 , 62 ], parece ser prometedor, ya que primero, permite complementar el entendimiento de los distintos mecanismos de fallo en los materiales compuestos, y segundo, permitiría ayudar con el cálculo de las cargas de fallo en el diseño de estructuras de material compuesto. 1.2 Daño en los Materiales Compuestos Los materiales compuestos reforzados con fibras son habitualmente diseñados para que sus fibras trabajen en la dirección longitudinal, y esto los vuelve más vulnerables al fallo en la dirección transversal. En particular, el fallo transversal se observa en el caso de laminados multidireccionales que contienen láminas de 90 grados y también en el caso de laminados sometidos a impactos. A nivel macromecánico, el principal mecanismo de fallo es la propagación de grietas de despegue entre las diferentes láminas de un laminado. A nivel micromecánico, se distinguen el fallo de la interfase fibra-matriz y el fallo de la matriz [ 14 , 59 ], como se muestra en la figura 1.1, mientras que en la figura 1.2 se muestra un diagrama del plano de fallo asociado al daño fibra-matriz. El denominado fallo inter-fibras/matriz ya ha sido tratado en estudios anteriores [14, 45, 65] para el caso de modelos de una y varias fibras, donde se ha puesto en relieve la importancia de la interfase fibra-matriz en el desarrollo del mecanismo de daño. Esos estudios se sostienen en la hipótesis de que el fallo transversal empieza con la aparición de pequeñas grietas en la interfase y analizan la propagación a través de ésta, pudiéndose identificar las siguientes etapas (mostradas en la figura 1.3): •Etapa 1: nucleación de la grieta de interfase. •Etapa 2: crecimiento simétrico e inestable de la grieta de interfase, hasta un semiángulo de despegue entre 60-70 grados. • Etapa 3: la grieta a lo largo de la interfase cambia abruptamente su dirección, penetrando en la matriz en dirección perpendicular a la carga aplicada. 1.2 Daño en los Materiales Compuestos 3 Figura 1.1 (a) Despegues de fibra-matriz y (b) propagación hacia una grieta transversal. Fuente: [18]. Figura 1.2 Diagrama del plano de fallo asociado al despegue inter-fibras y matriz en una lámina unidireccional. Figura 1.3 Desarrollo esquemático del daño; (a) grietas pequeñas de interfase originadas por la tensión radial; (b) crecimiento inestable de grieta hasta alcanzar un semi-ángulo de despegue de 60-70 grados; (c) penetración de la grieta de interfase en la matriz; (d) Crecimiento inestable de la grieta en la matriz; (e) Daño real en un compuesto de fibra bajo carga transversal. Fuente: [45]. 4Capítulo 1. Introducción 1.3 Métodos de la Mecánica de la Fractura no Singular Los métodos basados en la mecánica de la fractura lineal elástica (LEFM) [ 36 ], permiten estudiar el crecimiento y propagación de una grieta ya existente, pero no el comienzo de la misma. Además, cabe notar que en la escala micromecánica de una lámina de compuesto, la mecánica de la fractura no tiene respuestas tan satisfactorias como la que pueden ofrecer los métodos energéticos en componentes sin fisuras [20]. 1.3.1 Modelos de zona cohesiva Entre los métodos nolineales de la fractura, enraizados en los conceptos establecidos por Barenblatt y Dugdale durante la década del 1960 [ 6 , 17 ] (ver también a Hillerborg [ 24 ]), el conocido modelo de zona cohesiva (CZM, también denominado modelo de grieta cohesiva) ha probado ser uno de los enfoques más exitosos para la caracterización del daño en sólidos. La idea clave de este modelo yace en el hecho de que el campo de tensiones infinitas alrededor del vértice de la grieta no es realista, ya que en esa zona las tensiones elásticas no superan la tensión de fluencia del material. Consecuentemente, para evitar esas singularidades, los autores introdujeron el concepto de tensiones cohesivas a lo largo de la denominada zona del proceso de fractura. La región justo después del borde de la grieta, bajo la cual micro-fisuraciones toman lugar, es considerada como una extensión de la grieta física en la cual una presión de cierre σY es aplicada, concepto que se ilustra en la figura 1.4. Figura 1.4 Placa infinita sometida a una tensión remota σ , con una grieta de ancho inicial 2a . Las zonas de fluencia Dugdale corresponden a las longitudes R0 bajo el campo de tensiones σY , haciendo la longitud efectiva de grieta 2(a+R0). Fuente: [58]. El elemento principal que describe la respuesta cohesiva es la denominada ley de tracción-separación. Esta ley gobierna la respuesta de decohesión del modelo de interfase en términos de máxima tracción, tenacidad a la fractura y criterio de fallo en modo mixto. La diferencia entre las distintas leyes de tracción-separación consiste en la forma de la evolución de las tensiones con el historial de desplazamiento, tal como se ilustra en la figura 1.5. En la literatura, el distinto número de formulaciones cohesivas que obedecen a distintas leyes de tracciónseparación es extenso. Las leyes de tracción-separación disponibles se extienden desde modelos bilineales [ 68 ] hasta los exponenciales [ 46 ]. Sin embargo, una estrategia de modelado razonable es el caso frágil de los modelos de zona cohesiva, usualmente referidos como modelo de interfase lineal-elástica frágil, ó en su abreviatura anglosajona, LEBIM. 1.4 Modelo de Interfase Lineal Elástica Frágil (LEBIM) 5 Figura 1.5 Modelos de zona cohesiva. (a) Dugdale, (b) modelo de ablandamiento lineal, (c) modelo trapezoidal, (d) modelo exponencial. 1.4 Modelo de Interfase Lineal Elástica Frágil (LEBIM) Como en este proyecto se ha de modelar el fallo del despegue inter-fibras en la escala micromecánica, junto a la posibilidad del daño de la matriz adyacente a las fibras, es conveniente primero darle cabida a una descripción del modelo de interfase débil. La idea central del modelo de interfase se basa en el empleo de la mecánica de medios continuos, para modelar una (finita) capa de adhesivo de pequeño espesor, h , que une a dos sólidos y que se ilustra en la figura 1.6. La definición de una interfase de espesor finito evita la implementación de estrategias numéricas requeridas para insertar elementos de interfase de espesor nulo. Figura 1.6 El modelo de interfase LEBIM, como una distribución continua de muelles a lo largo de una interfase; Fuente: [65]. 1.4.1 Ley constitutiva de la capa de adhesivo Con el LEBIM, la interfase de adhesivo no dañada se modela como una distribución continua de muelles, donde cada punto en ella está gobernado por la siguiente ley constitutiva, Interfase Lineal Elástica (puntos sin daño): σ(x) = knδn(x) τ(x) = ktδt(x),(1.1) donde σ(x) y τ(x) denotan respectivamente, a las componentes normal y tangencial del tensor de tensiones en un punto x de la interfase; kn y kt representan las rigideces normal y tangencial de la distribución de resortes, y finalmente, δn(x) y δt(x) son los desplazamientos relativos normales y tangenciales de puntos opuestos de la interfase. 6Capítulo 1. Introducción Mientras que para un punto x de la interfase que esté roto (dañado), la ley constitutiva se define de la siguiente manera: σ(x) = kn⟨δn(x)⟩− τ(x) = 0,(1.2) donde el operador ⟨⟩− denota la parte negativa de un número real. El comportamiento descrito por las ecuaciones 1.1 y 1.2, y que se ilustra en la figura 1.7 se puede analizar como sigue: cuando el índice de liberación de energía (ERR) en modo mixto, G(x) , supera el valor critico Gc(x) (i.e., la tenacidad a la fractura del material), ocurre la rotura del punto. Cabe notar que el valor Gc(x) es dependiente de los valores críticos de las componentes normales y tangenciales del tensor de tensiones, σc(x) y τc(x) . En un punto roto, las tensiones tangenciales se vuelven nulas, mientras que las tensiones normales pueden existir si el punto está trabajando en compresión, manteniendo su rigidez inicial, kn (i.e., se considera la existencia de condición de contacto tipo penalty). Figura 1.7 Ley tensión-deformación del LEBIM. Izquierda: comportamiento en modo-1; Derecha: comportamiento en modo 2. Fuente: [65]. 1.4.2 Criterio de fallo de la interfase El criterio de fallo de la interfase tiene una base energética, donde la fractura de la interfase en un punto material x , se alcanza cuando el correspondiente índice de liberación de energía en modo mixto G(x) alcanza su correspondiente valor critico, i.e., G(x) = Gc(ψ(x)) . bajo condiciones generales de carga, el ángulo mixto de fractura energético ψse define como: tanψ=pξ−1tanψσ,(1.3) donde ξ=kt/kn y tanψσ=τ/σ , ψσ denota el ángulo mixto de fractura basado en tensiones en el punto x . Cabe notar que en la presente formulación, se asume que tan2ψ=tan2ψG=GII/GI para valores de σ>0 . Asumiendo una contribución aditiva de cada modo individual de fractura (ver p.ej., [ 35 , 29 , 32 ]), el índice de 1.4 Modelo de Interfase Lineal Elástica Frágil (LEBIM) 7 liberación de energía, ERR, de la interfase lineal elástica se define como G(x) = GI(x)+GII (x), con: GI(x) = h⟨σ(x)⟩+⟨εn(x)⟩+ 2,GII(x) = hτ(x)εn(x) 2,(1.4) verificándose GI(x) = 0para εn(x)≤0. La dependencia de la tenacidad a la fractura Gc con el ángulo mixto del modo de fractura se define como el reportado por Hutchinson y Suo [25]: Gc=¯ GIc[1+tan2((1−λ)ψ)],(1.5) donde: ¯ GIc =¯ σc¯ εnch 2=¯ σ2 c 2kn =kn¯ ε2 nch 2(1.6) representa la energía de fractura en modo puro de apertura I . En la ecuación 1.5, λ es el parámetro de sensibilidad del modo de fractura (normalmente obtenido del ajuste de los resultados experimentales). ¯ σc=σc(ψ=0)>0 , y ¯ εnc =εnc(ψ=0)>0 son las tracciones y deformaciones criticas asociadas al modo I, respectivamente, mientras que τc y εtc corresponden a los valores críticos de tensión y deformación tangencial, respectivamente. En relación al parámetro λ , un rango típico de 0.2≤λ≤0.3 caracteriza interfases con una dependencia moderadamente intensa del modo de fractura [ 25 ]. De la figura 1.8 se puede observar el desarrollo de Gc GIc en función del ángulo ψ para varios valores de λ ; cabe notar la mayor intensidad de Gc con λ=0 . En la actual formulación del LEBIM, tanto el criterio energético como el de tensiones son equivalentes. Entonces, las siguientes relaciones entre ¯ GIc,¯ σc,E′′ , y el espesor, h, del adhesivo se obtienen: ¯ GIc =¯ σ2 c 2E′′ ,(1.7) donde es interesante notar que E′′ =kn en la presente formulación. Una discusión más detallada de estos tópicos es dada por V. Mantič y colaboradores (sección 2.3) en [ 35 ]. Más detalles de la formulación del criterio de fallo previamente establecido se brinda en [61, 63, 64]. Figura 1.8 Energía de fractura Gcen función de ψG, para distintos valores de λ. Examinando la anterior derivación, cabe notar que uno de los aspectos más atractivos del presente modelo de interfase yace en el hecho de que el LEBIM se caracteriza por sólo cuatro parámetros independientes, expresamente: ¯ σc (el valor de la tensión normal crítica), ¯ GIc (la tenacidad a fractura en modo I), ξ (el cociente de las rigideces tangencial y normal de la interfase, kt kn ), y λ (el parámetro de sensibilidad del modo de 8Capítulo 1. Introducción fractura). Se debe notar que en el presente modelo LEBIM, el espesor de la interface, h, no tiene influencia en los resultados ( G y Gc son proporcionales a h , revisar las ecuaciones 1.4-1.6), para capas de interfase flexibles y muy delgadas. Sin embargo, hay que dejar expreso que en contraste a las formulaciones de zona cohesiva que usan definiciones alternas de las leyes tracción-separación, en el presente modelo LEBIM sólo puede producirse un daño parcial de la interfase, ya que la variable interna del daño sólo caracteriza estados de daño binarios (i.e., ninguno ó daño completo). Finalmente, el LEBIM puede ser considerado como un límite frágil de un modelo de zona cohesiva bilineal, y este es limitado a configuraciones donde el tamaño de la zona de proceso/ablandamiento es pequeña en comparación con el tamaño de la grieta. Ha sido demostrado por Wang [ 70 ], que el efecto de la forma de la ley cohesiva puede ser descartado si los valores de la tenacidad a la fractura y la tensión crítica son iguales, y las zonas de proceso/cohesivas son pequeñas. Escribiendo el criterio de propagación de la grieta a lo largo de la interfase lineal elástica frágil, la siguiente expresión para la tensión normal critica en función del ángulo ψGes obtenida: σc(ψG) = ¯ σcq1+tan2((1−λ)ψG)cosψG δnc(ψG) = σc(ψG) kn (1.8) De manera similar, la tensión tangencial crítica en el modo mixto, se puede expresar en términos de ¯ σc y ψG , como: τc(ψG) = skt kn ¯ σcq1+tan2((1−λ)ψG)sinψG δtc(ψG) = τc(ψG) kt (1.9) Entonces, la condición que regirá el daño de la interfase se expresa como: f actor =sG(x) Gc ,(1.10) donde qG(x) Gc mantiene una interacción lineal del estado tensional de la interfase (σ,τ) con la tensión critica ¯ σc (ver de nuevo las ecuaciones 1.8 y 1.9, con el término p1+tan2((1−λ)ψG) despejado); de la definición en la ecuación 1.10, el daño material se produce si f actor ≥1. El gráfico de las curvas de fallo de la interfase, parametrizado por las ecuaciones 1.8 y 1.9, se muestra en la figura 1.9 para 0≤ψG≤π 2 , donde una relación kt kn=0.25 ha sido considerada. Las tensiones críticas normales y tangenciales en modo mixto han sido normalizadas por la tensión crítica de rotura en el modo I, ¯ σc. 1.5 Algoritmo de Solución Empleado 9 Figura 1.9 Envolventes de fallo de la interfase, en el plano (σ,τ) para 0≤ψG≤π 2 y distintos valores de λ . 1.5 Algoritmo de Solución Empleado Usualmente, los modelos de grieta de interfase presentan inestabilidades numéricas post-fractura que tienden a ser difíciles de capturar usando los algoritmos de solución de Newton-Raphson y el método Arc-Length modificado de Riks; esto es debido a que después de los primeros despegues en adhesivos ó microfisuraciones en p.ej., el concreto, se puede producir una abrupta pérdida de rigidez estructural, generando comportamientos del tipo "snap-back" con un perfil de ablandamiento de dientes de sierra en la curva de equilibrio del sistema (ver p.ej., [ 2 ]). Una manera de evitar este problema es el uso de un algoritmo de solución secuencial lineal, conocido como el SLA ("Sequential Linear Analysis") en abreviatura anglosajona. Como este método será implementado en este trabajo para estudiar el fallo en modelos unifibra y multifibras, se procederá a su descripción bajo el siguiente apartado. 1.5.1 Algoritmo de solución secuencial lineal (SLA) R. de Borst [ 15 ] ya había realizado estudios numéricos del comportamiento post-fallo en estructuras con ley de ablandamiento lineal, usando modificaciones del método arc-length. Sin embargo, este enfoque implica el uso de los controles de iteración y convergencia necesarios en un programa de elementos finitos no lineal. Luego, J. Rots hacia el 2001 [ 53 ], había implementado un algoritmo de solución secuencial lineal para modelizar la fractura en fachadas de concreto. Este modelo aproxima la curva de ablandamiento tensión-deformación en tracción como un diagrama en dientes de sierra. Después de un análisis lineal, el elemento crítico cuyo estado tensional está mas cercano al daño, es removido. Seguido, la rigidez de dicho elemento se reduce y el proceso se repite de nuevo. Cabe notar que este enfoque había tomado inspiración de ciertos autores para modelar la heterogeneidad de compuestos con base de cemento como una estructura de celosías. [55]. El algoritmo actual analiza tanto los puntos de integración de la interfase así como aquellos de la matriz para cada iteración del paso de carga. En este proyecto, se usa Abaqus para construir los modelos y resolver el sistema lineal de ecuaciones. El LEBIM es programado en lenguaje FORTRAN, siendo implementado en Abaqus por medio de una subrutina del material UMAT [ 11 ]. Una vez que el comportamiento de la interfase y la matriz se han definido, el algoritmo analiza cada punto del sólido identificando si está dañado o intacto. Procedimiento global La estructura se discretiza usando los elementos continuos lineales estándar. El módulo elástico E , la relación de Poisson ν , y los parámetros del criterio de fallo correspondiente son asignados a los elementos. Seguido, los siguientes pasos se ejecutan en orden secuencial: 16 Capítulo 2. Criterios de Fallo de la Matriz + Implementación Abaqus UMAT linealidad del problema elástico para realizar dicha superposición de energías. Acorde con lo anterior: Ut=Uv+Ud(2.4) donde la componente volumétrica se define como: Uv=1 2I1(σi j)·εv=1−2ν 6E(σxx +σyy +σzz)2(2.5) Por otra parte, la energía critica Ucr, de la matriz se define como: Ucr =1 2 S2 t E(2.6) cuya unidad en el sistema internacional se da en J/m3 . Sin embargo, y de manera similar a como se establece en la ecuación 1.10 del capítulo 1 para el factor que define el fallo de la interfase, el criterio energético que define el fallo de la matriz tiene que estar elevado al exponente 1/2 para que se pueda mantener una relación lineal con las tensiones, manteniendo entonces una equivalencia entre el criterio energético y el de tensiones (se recomienda ver p.ej., [31, 56]), de manera que el fallo se producirá cuando: f actor =sUv Ucr ≥1(2.7) El factor de la matriz en la anterior ecuación se expresa de dicha forma debido a que tanto el criterio energético como el de tensiones forman condiciones necesarias para la fractura (dígase, la nucleación de grieta), pero ninguno de ellos es suficiente por sí solo [31]. 2.1.5 La densidad de energía total de deformación Expandiendo del apartado anterior, se toma en cuenta la componente de distorsión (aquella que provoca los cambios de forma) de la energía de deformación en la ecuación 2.4, donde Udse define como: Ud=1+ν 3E[σ2 I+σ2 II +σ2 III −(σIσII +σII σIII +σIσIII)] (2.8) Similar a la relación 2.7, este criterio establece que el fallo de la matriz se producirá cuando la energía total almacenada alcance un valor crítico: f actor =sUv+Ud Ucr ≥1(2.9) Con esto termina la exposición relativa a las consideraciones de los criterios de fallo para la matriz. Ahora, se pasará a establecer un preámbulo describiendo las convenciones de las componentes de tensión que Abaqus usa dentro de la subrutina UMAT y su equivalencia con las tensiones actuantes en un plano micromecánico del fallo de una lámina unidireccional de compuesto. 2.2 Convención de las Componentes de Tensión y Deformación 17 2.2 Convención de las Componentes de Tensión y Deformación Cuando se definen las propiedades del material, la convención usada para las componentes de tensión y deformación en Abaqus son las siguientes: •σ11: Tensión normal en la dirección-1 •σ22: Tensión normal en la dirección-2 •σ33: Tensión normal en la dirección-3 •τ12: Tensión de cortante en el plano 1-2 •τ13: Tensión de cortante en el plano 1-3 •τ23: Tensión de cortante en el plano 2-3 Las direcciones 1-, 2-, y 3- dependen del tipo de elemento escogido. En nuestro caso, para los elementos sólidos continuos como los elementos de deformación plana usados para los análisis en este trabajo, dichas direcciones por defecto corresponden las direcciones espaciales globales. Como se estudian modelos unifibra tanto en una matriz finita e infinita, así como en modelos multifibras de una lámina de 90 grados usualmente comprendida entre láminas de cero grados, la nomenclatura de las tensiones en el plano transversal a las fibras se reproduce en los modelos de la figura 2.3, i.e, se representan las tensiones que actúan en el plano 2−3 en ejes de una lámina; cabe notar que la dirección de estas tensiones no tienen nada de relación con las direcciones espaciales globales, hasta el momento. Cuando se genera un modelo de elementos finitos en deformación plana, donde las direcciones 2-3 de una lámina coinciden con las direcciones espaciales globales 1-2 (el plano X−Y ), entonces ya se tiene referencia de las componentes de tensión en un plano transversal a las fibras, para relacionarlas con la dirección de las tensiones en el plano global X−Y de Abaqus (cabe notar que Abaqus ordena el definir los elementos de deformación plana en el plano X−Y[1]), concepto que se ilustra en la figura 2.4. 18 Capítulo 2. Criterios de Fallo de la Matriz + Implementación Abaqus UMAT Figura 2.3 Plano de fallo (transversal a las fibras), y dirección de las tensiones asociadas a dicho plano [ 44 ]. Figura 2.4 Modelo de elemento finito en deformación plana. En Abaqus, el número de componentes del tensor de tensiones asociado a un elemento de deformación plana son cuatro, definidas como sigue: •S11: XX, componente normal de tensión. •S22: YY , componente normal de tensión. •S33: ZZ, componente normal de tensión. 2.3 Variables de la UMAT 19 •S12: XY, componente a cizalladura de la tensión. Por tanto, la convención a seguir para asociar las variables del vector de tensiones de la subrutina UMAT con las tensiones ubicadas en un plano transversal a las fibras es como sigue: •ST RESS(1) = S11:→equivale a σ22 de la lámina •ST RESS(2) = S22:→equivale a σ33 de la lámina •ST RESS(3) = S33:→equivale a σ11 de la lámina •ST RESS(4) = S12:→equivale a σ23 de la lámina La convención a usar para las variables de deformacion sigue el mismo orden que la usada para definir las componentes del tensor de tensiones en un elemento de deformacion plana. 2.3 Variables de la UMAT En este apartado se describirán las variables principales cuya definición dentro de la subrutina es esencial para el funcionamiento de la misma y que a partir de las cuales se modela el comportamiento de la interfase y de la matriz de manera independiente; estas variables podrían tener un valor inicial asignado que luego la subrutina actualiza para ser devuelto a Abaqus. Variables a ser definidas en todos los casos: •DDSDDE(NTENS,NTENS) : define la matriz jacobiana del modelo constitutivo, ∂∆σ/∂∆ε , donde ∆σ es el incremento de las tensiones y ∆ε es el incremento de las deformaciones. DDSDDE(I,J) define el cambio de la componente- I de tensión al final del incremento de tiempo (i.e., carga) causada por una perturbación infinitesimal de la componente- J del vector de incremento de deformaciones. A menos que se invoque la componente no simétrica del modelo de material definido por el usuario, Abaqus sólo usará la parte simétrica de DDSDDE. •STRESS(NTENS) : este vector se pasa como el tensor de tensiones al inicio del incremento de carga y se tiene que actualizar en esta subrutina para pasar a ser el tensor de tensiones al final del incremento. El tamaño (ó longitud) de este vector depende del valor NTENS , que para los modelos de deformación plana, pasa a ser NTENS=4 . La medida de tensión usada por Abaqus es la "real" o la llamada tensión de Cauchy. •STATEV(NSTATV) : es un vector que almacena las variables de estado dependientes de la solución (como p.ej., las tensiones y deformaciones). En cualquier caso STATEV debe ser pasado como los valores de las variables al final del incremento. -Variables a ser pasadas como información: •STRAN(NTENS): un vector que contiene las deformaciones totales al inicio del incremento de carga. •DSTRAN(NTENS): es el vector incremento de deformaciones. •NTENS : es el número de componentes del vector de tensiones y deformaciones ( NDI+NSHR ), siendo NDI el número de componentes normales y NSHR el número de componentes de cizalladura. En deformación plana, NTENS=4; en tensión plana, NTENS=3; para elementos sólidos tridimensionales, NTENS=6. •NSTATV : es el número de variables dependientes de la solución que están asociadas con un material definido por el usuario. En la interfaz gráfica de Abaqus, esta variable es el primer cuadro de diálogo dentro del comportamiento de material "Depvar" bajo la pestaña "General" y que se muestra en la figura 2.5. El "Depvar" tambien se puede añadir y modificar en el fichero de entrada (.inp) que genera Abaqus para ser resuelto por su programa de elementos finitos, de la siguiente manera (dentro del encabezado de los materiales definidos): ** MATERIALS ** *Material, name=Material-1 *Depvar 3, 20 Capítulo 2. Criterios de Fallo de la Matriz + Implementación Abaqus UMAT Figura 2.5 Localización por medio de la interfaz gráfica de Abaqus, de la variable NSTATV. •NOEL: número del elemento. •NPT: número del punto de integración •KINC: número del incremento. •KSTEP: número del paso. 2.4 Funcionamiento General de la Subrutina En este apartado, se dará una explicación general del funcionamiento de cómo son procesadas la variables en la subrutina. Sin embargo y previo a esto, se debe entender la manera en que Abaqus resuelve los modelos de elementos finitos mediante un análisis lineal estático, que es el caso en este proyecto ya que los efectos de inercia no intervienen en el modelo físico propuesto. Entre los pasos básicos para analizar un modelo de elementos finitos, se encuentra la escritura del fichero de entrada (en este caso un archivo del tipo .inp). Cuando se envía un "job" asociado con un modelo para su análisis, Abaqus/CAE genera primero un fichero de entrada representando el modelo creado, y entonces Abaqus/Standard realiza el análisis con los contenidos de este fichero. De manera aparte, se le puede pedir a Abaqus que genere exclusivamente el fichero de entrada haciendo clic en la pestaña superior derecha del Job Manager, como se muestra en la figura 2.6; Abaqus/CAE escribe entonces un fichero en formato ASCII, que se puede ver y editar en el directorio de trabajo. El concepto del fichero de entrada se introduce en este apartado, a razón de que el nombre de este archivo será requerido en el fichero de la subrutina UMAT y como dato de entrada para el algoritmo SLA, como se verá en el próximo capítulo. Para que Abaqus pueda resolver adecuadamente el sistema lineal de ecuaciones que ensambla, siguiendo las instrucciones del algoritmo secuencial lineal (ver sección 1.5.1), el usuario tiene que poner en el módulo del Step (para el procedimiento: Static, General) dentro de las opciones del periodo de tiempo y el máximo 2.4 Funcionamiento General de la Subrutina 21 Figura 2.6 Cómo generar el fichero de entrada en el Job Manager de Abaqus. número de incrementos del paso de carga el valor unidad ( 1 ), para que se resuelva el sistema de ecuaciones [K]{u}={F} con la magnitud total de la carga externa {F} aplicada al modelo de elementos finitos, para cada paso del SLA. Dichas opciones se ilustran en la figura 2.7, y se escriben en el fichero de entrada .inp del modelo de elementos finitos con el siguiente formato: ** STEP: Step-1 ** *Step, name=Step-1, nlgeom=NO, inc=1 *Static 1., 1., 1., 1. (a) periodo de tiempo del "step" (b) opciones del incremento Figura 2.7 Opciones del paso lineal estático para el modelo en Abaqus; (a) periodo de tiempo y (b) las opciones del incremento que le serán asignadas valores unitarios. 22 Capítulo 2. Criterios de Fallo de la Matriz + Implementación Abaqus UMAT En cada paso de carga que lanza el SLA, Abaqus llama a la subrutina UMAT, donde se forma la matriz de constantes elásticas del material para cada elemento, Cel , que posee las propiedades mecánicas asignadas en el modelo de comportamiento definido por el usuario. En cada paso de carga, la subrutina inicializa en cero las siguientes variables que se pasan como información de cada elemento: •DDS: matriz de elasticidad isótropa del elemento. •STRESS: vector de tensiones del elemento. •ψG: ángulo mixto del modo de fractura. •GI: índice de liberación de energía en modo I de fractura. •GII: índice de liberación de energía en modo II de fractura. y determina si el elemento que se evalúa en dicho paso pertenece a la interfase o la matriz, pasando a leer las propiedades del modelo de material definido por el usuario ( PROPS en la UMAT) para este elemento, evaluando entonces la matriz de elasticidad y el criterio de fallo establecidos para el material, seguido del cálculo de las variables dependientes del estado ( SDV ) en cada punto de integración del elemento. Las variables de estado se pasarán a Abaqus para ser guardadas en el fichero de resultados (un archivo del tipo .odb). 2.5 Estructura de la Subrutina La subrutina ha sido programada en el lenguaje de programación FORTRAN. A pesar de que las sentencias en las líneas de código usadas para la programación de la subrutina de este proyecto no es muy compleja, se recuerda al usuario que en FORTRAN (en contraste con otros lenguajes de programación como MATLAB), antes de usar una variable, se la tiene que declarar con el tipo adecuado, ya sea p.ej., un número entero ( INTEGER ), o de doble precisión ( REAL*8 ). Los enlaces introductorios "Fortran Variables" 1 y "Fortran Best Practices" 2 se dejan aquí para que el usuario se pueda familiarizar con los tipos de variables y las mejores prácticas de programación en el lenguaje FORTRAN. La sección de encabezado de la subrutina UMAT implementada se muestra en las siguientes líneas de código: Códigos 2.1 Encabezado de la UMAT. SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT, 2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, 3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT, 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC) implicit none ! maximum number of integration points and elements logical lee integer nmax_pt, nmax_el parameter (nmax_pt=4, nmax_el=600000) integer dama(nmax_pt,nmax_el) common lee, dama data lee /.true./ ! Como no cargamos el archivo 'ABA_PARAM.INC'debemos de inicializar las variables CHARACTER*80 CMNAME REAL*8 STRESS(NTENS), 1 DDSDDE(NTENS,NTENS), 2 DDSDDT(NTENS),DRPLDE(NTENS),STATEV(nstatv), 1https://fortran-lang.org/en/learn/quickstart/variables/ 2https://fortran-lang.org/en/learn/best_practices/ 2.5 Estructura de la Subrutina 23 3 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), 4 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3) INTEGER PT,ELEM INTEGER error_ap,error_lec INTEGER k1,k2,k,j,i,ndi,nshr,ntens,nprops,noel,npt INTEGER layer,kspt,kstep,kinc,nstatv,mtype,damage REAL*8 DSTRESS(4),DDS(4,4),psig,lambda,GIcb REAL*8 sigmacb,Pi,Gi,Gii,Gtot,h,ktkn,Gc,psiGcrit REAL*8 E,xnue,ebulk3,eg2,elam,trval,Gs,Gn,Knn,Ktt,Kss,K33 REAL*8 sse,spd,scd,rpl,drpldt,dtime,temp,dtemp,celent,pnewdt REAL*8 factor,sigmaprin1,sigmaprin2,VM,sigmayt,Suc,Ut,Uv,Ud,Uc En el código 2.1, se comenta que como no se incluye el archivo ABA_PARAM.INC , entonces se establecen valores iniciales a las variables STRESS , DDS y PSIG . En las últimas líneas, se definen las variables de tipo numérico REAL*8 que reciben las propiedades mecánicas de la ley constitutiva del material tanto de los elementos de la interfase así como los de la matriz. Los comandos INTEGER y REAL*8 en FORTRAN declaran variables numéricas de tipo entero y del tipo doble precisión de 8 bytes con 15 dígitos de exactitud, respectivamente. El bloque de inicialización de las variables se detalla en las siguientes líneas como: Códigos 2.2 Inicialización de las variables de la UMAT. ! Inicializacion Mutex call MutexInit(1) ! Inicializacion variables DDS=0.d0 psig=0.d0 STRESS=0.d0 Pi=ACOS(-1.d0) Gi=0.d0 Gii=0.d0 Donde se inicializan en cero ( 0.d0 ) los valores de la matriz de elasticidad DDS , el vector de tensiones STRESS , el ángulo del modo mixto de fractura de base energética psig , y los índices de liberación de energía del modo I, Gi, y del modo II, Gii. A diferencia de anteriores trabajos que sólo emplearon las propiedades y los criterios de fallo de la interfase para el modelo del material [ 11 , 39 , 65 ], en este proyecto, para modelar el daño de la matriz, se añade el comportamiento material definido por el usuario ("User Material") para los elementos de la matriz, y por tanto ya no se puede emplear el modelo de comportamiento lineal elástico isotrópico definido por defecto en Abaqus, siendo entonces la sintaxis que define el modelo material de la matriz dentro del fichero .inp como sigue: ** MATERIALS ** *Material, name=Epoxy *Depvar NSTATV, *User Material, constants=n PROPS(1), PROPS(2), PROPS(3), PROPS(4), ...,PROPS(n) donde NSTATV es el número de variables dependientes de la solución para cada punto de integración en cada elemento de la matriz, y n es el número de constantes del material usadas para definir la ley constitutiva y que será dependiente del criterio de fallo empleado, como se presentará en la sección 2.7. 24 Capítulo 2. Criterios de Fallo de la Matriz + Implementación Abaqus UMAT 2.6 Implementación del LEBIM en UMAT Bajo esta sección se darán a conocer los detalles de implementación del modelo de interfase LEBIM dentro de la subrutina. Sin embargo, primero tenemos que presentar la sintaxis de las variables de comportamiento pertinentes al LEBIM dentro del fichero .inp que genera Abaqus, para tener un claro entendimiento de su significado. La sintaxis de la tarjeta del material de la interfase en el fichero de entrada se estructura como sigue: *Material, name=InterLEBIM *Depvar 3, *User Material, constants=6 PROPS(1), PROPS(2), PROPS(3), PROPS(4), PROPS(5), PROPS(6) donde: •PROPS(1) : corresponde a la variable mtype de la UMAT; es un número entero que siempre es igual a uno para la interfase; mayor e igual a dos para la matriz. •PROPS(2): es el valor numérico de la tensión crítica de fallo en el modo I, ¯ σc. •PROPS(3): corresponde al índice de liberación de energía crítico de la interfase en modo I, ¯ GIc. •PROPS(4): es el parámetro de sensibilidad del modo de fractura, λ. •PROPS(5): es la relación de las rigideces de la interfase, kt kn(usualmente es 0.25). •PROPS(6): es el valor numérico del espesor de la interfase LEBIM, h. El código 2.3 muestra cómo la subrutina procesa el tipo de material y las propiedades que se le asignan dependiendo del valor de mtype. Códigos 2.3 Lectura y asignación de las propiedades mecánicas de la interfase (mtype=1) y la matriz (mtype≥2) mediante la subrutina UMAT. ! !READ THE MATERIAL PROPERTIES FROM THE INPUT FILE mtype=PROPS(1) if (mtype.eq.1) then sigmacb=PROPS(2) GIcb=PROPS(3) lambda=PROPS(4) ktkn=PROPS(5) h=PROPS(6) else if (mtype.gt.1) then E=PROPS(2) XNUE=PROPS(3) EBULK3=E/(1-2*XNUE) EG2=E/(1+XNUE) ELAM=(EBULK3-EG2)/3 El pseudocódigo que implementa el modelo de material de la interfase LEBIM en la subrutina UMAT se detalla en 2.8. 2.7 Implementación del Daño de la Matriz en UMAT 25 Input: parámetros ( ¯ σc,¯ GIc,λ,ξ,h) de la interfase, junto a los datos (∆εn,εn,σn,(SDV1)n) Result: Variables dependientes del estado SDV1, SDV2, SDV3 foreach punto de integración n=1:4del elemento NOEL do if (SDV1)n<1then Empieza el LEBIM sin daño: Actualizar el tensor de elasticidad Cinterf Actualizar el incremento de tensiones, ∆σn=DDS(k)·DSTRAN(k) Actualizar el tensor de tensiones σn+1=σn+∆σn Evaluar la variable de daño (SDV 1)n Calcular los ERRs de la interfase GI,GI,GcyψG Evaluar la tenacidad a la fractura Gc(x)[Ec. 1.4] Calcular el f actor =pGc(x/Gc)[Ec. 1.10] if f actor ≥1then damage = 1 else damage = 0 else if (SDV1)n=1then Empieza el LEBIM con daño: f actor =0 Actualizar la matriz de elasticidad, Cinterf, y el tensor de tensiones if εn(1)<0then σ(k)=0 for k=2, 3, 4 else if εn(1)>0then σ(k)=0 for k=1, 2, 3, 4 Actualizar la matriz de rigidez tangente DDSDDE(i,j)=DDS(i,j) Actualizar los valores de las variables dependientes del estado: SDV1 = damage SDV2 = factor SDV3 = NOEL Figura 2.8 Pseudocódigo de la implementación LEBIM en la UMAT. 2.7 Implementación del Daño de la Matriz en UMAT El daño de los elementos de la matriz se implementa en la UMAT con la variable mtype≥2 . La sintaxis de los datos en la tarjeta del material de la matriz en el fichero de entrada (.inp) se escribe como sigue: *Material, name=Epoxy *Depvar 3, *User Material, constants=4 PROPS(1), PROPS(2), PROPS(3), PROPS(4),..., PROPS(n) donde: •PROPS(1) : es la variable mtype ; es un número entero igual y mayor a dos para los elementos de la matriz. •PROPS(2): valor numérico que corresponde al módulo elástico de la matriz, E. •PROPS(3): valor numérico que corresponde al coeficiente de Poisson de la matriz, ν. •PROPS(n) , para n≥4 : corresponden a las propiedades mecánicas (resistencias o energías) relacionadas a la hipótesis de fallo de la matriz definido por el valor de mtype. 32 Capítulo 3. Implementación SLA + Python sus actuales propiedades resistentes. La carga global es entonces escalada de acuerdo al salto de daño en dicho punto de integración. El modelo es sencillo y siempre brinda una solución: la rigidez tangente de los dientes de sierra siempre es positiva, de manera que la divergencia en la solución no ocurre en un análisis secuencial lineal. Una explicación física de este modelo es que la fractura consiste en un proceso gradual de separación donde la sección neta que conecta el material, y por tanto la rigidez, es gradualmente reducida. En los siguientes apartados, se introduce la herramienta SLA, desarrollada en Python, y que ya ha sido implementada para el programa Abaqus en anteriores trabajos (ver capítulo 3 de [39], y en [65]). 3.1 Componentes del SLA y Funcionamiento General El algoritmo secuencial lineal (SLA) de este proyecto ha sido elaborado para funcionar en los servidores de cálculo del Grupo de Elasticidad y Resistencia de Materiales (GERM) de la Universidad de Sevilla, incorporando ciertos elementos para su adecuación. A pesar de que se use Abaqus/CAE en un ordenador personal para elaborar los modelos de elementos finitos, la llamada del código SLA y la solución del sistema lineal de ecuaciones serán realizadas desde el servidor de cálculo. Los componentes del algoritmo secuencial lineal en el servidor de cálculo son los siguientes: •El fichero .inp que genera Abaqus. •La subrutina de material UMAT: es un archivo del tipo Fortran Source con extensión .f1. • El fichero inputSLA.txt : es un archivo de texto que contiene información que será leída por un archivo llamado main_abaqus.py. • El fichero main_abaqus.py : es el archivo Python que lee el fichero .inp y lanza el job de Abaqus con la subrutina en el servidor de cálculo. • El fichero postprocess.py : un archivo Python llamado por main_abaqus.py en cada iteración del SLA y que post-procesa las variables dependientes de estado (SDV) de la subrutina UMAT. • Un fichero run-test : un archivo de texto que se carga en la carpeta del usuario del servidor de cálculo, y que contiene la información del ordenador y el número de procesadores para resolver el job .inp de Abaqus. 1https://fortranwiki.org/fortran/show/File+extensions 3.2 Fichero Principal: main_abaqus.py 33 En la figura 3.3 se muestra la implementación en pseudocódigo del método secuencial lineal, y cuyos componentes serán explicados en los siguientes apartados. Input: nombre del fichero .inp, subrutina UMAT, número de iteraciones del SLA, frecuencia de grabado de resultados .ODB, puntos de integración a dañar por iteración. Result: ficheros Abaqus .ODB, fichero de texto con los datos de la variable factor (SDV2) por cada iteración del SLA. Lanzar el script de Pyhon/Abaqus /* run-test lanza el fichero main_abaqus.py en servidor de cálculo */ forall step=1 : número iteraciones SLA, do comienza iteración del SLA Crear el fichero .inp_numiteracionesSLA_step Lanzar el programa Abaqus para resolver el problema lineal elástico [K]{u}={F}de .inp_numiteracionesSLA_step /* main_abaqus realiza el llamado de postprocess.py */ foreach punto de integración xde la interfase LEBIM y de la matriz, do Leer el archivo datos_procesar.txt que incluye los puntos de integración x previamente dañados if xno está en el archivo datos_procesar.txt, then Evaluar la variable f actor en x (punto de integración sin daño), [Ecs. 2.1,2.2,2.3,2.7,2.9] else Definir factor = 0 (punto de integración con daño) Crear una tabla ordenada con la variable f actor (de mayor a menor) en cada iteración del SLA Imprimir el valor máximo del f actor en cada iteración del SLA en un archivo de texto Figura 3.2 Pseudocódigo del método secuencial lineal (SLA) aplicado a un modelo con ley de comportamiento de la interfase y matriz implementados en una subrutina UMAT; la variable factor corresponde con el múltiplo de carga que resulta del criterio de fallo empleado. 3.2 Fichero Principal: main_abaqus.py El main_abaqus.py es la columna vertebral del código SLA, ya que tomando los datos que suple el fichero inputSLA.txt , lanza el job de Abaqus y postprocesa los resultados requeridos para cada iteración. Para que este código funcione adecuadamente y de acuerdo con los datos suplidos por inputSLA.txt , se procede a explicar sus componentes. inputSLA.txt es el archivo de texto que contiene información del fichero de entrada de Abaqus, la subrutina UMAT y las opciones del método secuencial lineal. Su sintaxis se presenta a continuación: Códigos 3.1 Contenidos del archivo inputSLA.txt. name_INP. Nombre del fichero INP singlefibre name_UMAT. Nombre fichero UMAT LEBIMdefplana num_iteraciones. N iteraciones que resuelve SLA 400 paso. Frecuencia de guardado ficheros ODB (0 solo guarda el ltimo) 10 tol_o_ptosint. Si < 1 Modo Tolerancia. Si >=1 Modo n ptos integraci 4 cpu_num. Numero de procesadores (Debe coincidir con ppn en pythonpbs) 2 donde: •name_INP: es el nombre (sin escribir la extensión) del archivo .inp que genera Abaqus/CAE. 34 Capítulo 3. Implementación SLA + Python •name_UMAT : es el nombre de la subrutina (sin incluir la extensión) UMAT que será usada por Abaqus en cada iteración del SLA. •num_iteraciones : es el número (entero) de pasos lineales de carga con los cuales el usuario controla la cantidad de veces que el problema elástico es resuelto de forma secuencial. •paso : es un número entero positivo que asigna la frecuencia de grabado para los archivos .ODB que Abaqus genera. •tol_o_ptosint : un valor numérico, que si es mayor o igual que 1, representa el número de puntos de integración a dañar por iteración del SLA, y si es menor de 1, se impone un rango de tolerancia dado por el usuario y que daña los puntos de integración dentro del rango (1-tolerancia) ×max(f actor) . En este trabajo se ha empleado mayormente un tol_o_ptosint≥ 1, omitiendo los casos para el modo tolerancia. •cpu_num: un número entero positivo, por defecto se puede dejar en 4. 3.2.1 Lanzamiento del job Abaqus En el fichero main_abaqus.py, hay cuatro comandos que en cada iteración del SLA, se encargan de la lectura, creación, lanzamiento del job y retornar el control del programa tras la finalizacion del job Abaqus: mdb.ModelFromInputFile() mdb.Job() mdb.jobs[].submit() mdb.jobs[].waitForCompletion() La definición del job de Abaqus se realiza con el siguiente comando: mdb.Job(atTime=None, contactPrint=OFF, description='', echoPrint=OFF, explicitPrecision=SINGLE, getMemoryFromAnalysis=True, historyPrint=OFF, memory=90, memoryUnits=PERCENTAGE, model=name_INP, modelPrint=OFF, multiprocessingMode=THREADS, name=name_files, nodalOutputPrecision=SINGLE, numCpus=int(cpu_num), numDomains=int(cpu_num), numGPUs=0, queue=None, resultsFormat=ODB, scratch='',type=ANALYSIS, userSubroutine='LEBIMkinkpstrain_uni.f', waitHours=0, waitMinutes=0) Se ha verificado que el programa produce los resultados correctos si el número de CPUS se ajusta en 1, 2, o 4, aunque el tiempo de ejecución es esencialmente el mismo. Se hace notar que el número de dominios (numDomains=int()) debe ser un múltiplo del numero de CPUS. Como principal observación, el único elemento que se recomienda modificar dentro del fichero main_- abaqus.py es la variable userSubroutine dentro del objeto mdb.Job() , donde a userSubroutine se le pasa escrito entre comillas (") el nombre de la subrutina UMAT con su extensión. Un ejemplo de ello es: LEBIMkinkpstrain_uni.f , que se ilustra en el código 3.2. Es decir, si se carga una subrutina UMAT al servidor de cálculo para ser usada en el job de Abaqus, su nombre tiene que ser cambiado en el fichero main_abaqus.py. Después de entregar el job de Abaqus al servidor de cálculo y habiendo finalizado la solución de la iteración, main_abaqus.py llama a la función postprocess.py. Códigos 3.2 Ejecución del job por main_abaqus.py. #--------------------------------------------------------------- # BUCLE PRINCIPAL #--------------------------------------------------------------- ofile = open("output.txt",'w') pp_time = 0.0 # Post-processing time flagsize=0 for step in range(num_iteraciones): step_time = time.time() 3.3 Postprocesamiento de Resultados: el Fichero postprocess.py 35 # Dar nombre a los archivos de salida en la iteracion if (float(tol_o_ptosint) < 1): name_files = name_INP + '_'+str( num_iteraciones) + 'steps'+'_'+str(step + 1) + '_'+ tol_o_ptosint + ' tol' else: name_files = name_INP + '_'+str(num_iteraciones) + 'steps'+'_'+ str(step + 1) + '_'+ tol_o_ptosint + 'ptos_danar' #--------------------------------------------- # Ejecutar inp + UMAT #--------------------------------------------- # Read the model mdb.ModelFromInputFile(inputFileName=name_INP+'.inp', name= name_INP) # Create the job mdb.Job(atTime=None, contactPrint=OFF, description='', echoPrint=OFF, explicitPrecision=SINGLE, getMemoryFromAnalysis=True, historyPrint=OFF, memory=90, memoryUnits=PERCENTAGE, model=name_INP, modelPrint=OFF, multiprocessingMode=THREADS, name=name_files, nodalOutputPrecision=SINGLE, numCpus=int(cpu_num), numDomains=int(cpu_num), numGPUs=0, queue=None, resultsFormat=ODB, scratch='',type=ANALYSIS, userSubroutine='LEBIMkinkpstrain_uni.f', waitHours=0, waitMinutes=0) # Run the job mdb.jobs[name_files].submit(consistencyChecking=OFF) # Do not return control till job is finished running mdb.jobs[name_files].waitForCompletion() 3.2.2 Llamada de la función postprocesamiento Después de haber finalizado la solución del problema elástico para la iteración step del SLA, main_abaqus llama a la función postprocess.py, pasándole cinco argumentos mostrados en el código 3.3. Códigos 3.3 Llamada del postprocesamiento en main_abaqus.py. #--------------------------------------------- # POSTPROCESADO #--------------------------------------------- # Postprocesado ODB (modulo postprocess) time1 = time.time() postprocess.postprocess(name_files, tol_o_ptosint, name_txt_evolucion, str(step), working_directory) 3.3 Postprocesamiento de Resultados: el Fichero postprocess.py El fichero postprocess.py se encarga de extraer información de los resultados en el archivo .ODB de Abaqus para cada iteración del SLA. Recibe 5 argumentos en su llamada dentro de main_abaqus.py. Las primeras líneas de comando se escriben en el código 3.4, con los argumentos de llamado explicados como sigue: • name_files: es el nombre del archivo de entrada (sin extensión alguna) correspondiente a la iteración actual del SLA. 36 Capítulo 3. Implementación SLA + Python •tol_o_ptosint: número de puntos de integración a dañar por paso del SLA. • name_txt_evolucion: una cadena de caracteres con información del número de iteraciones, la tolerancia o puntos de integración a dañar. •str(step): número que corresponde con la iteración del SLA. •working_directory: representa el directorio de trabajo actual del usuario en el servidor de cálculo. El archivo de resultados ODB que genera Abaqus finalizada cada iteración del SLAse abre con el siguiente comando: odb = openOdb(arg1 + '.odb') Previo a introducir las operaciones que guardan e imprimen los resultados del fichero .odb, se introduce en el siguiente apartado los distintos repositorios del objeto de resultados en Abaqus, con el fin de entender el uso que postprocess.py hace de dichos objetos. Códigos 3.4 Importación de módulos y comandos iniciales de postprocess.py. #----------------------------------------------------------- # DEFINICION DE FUNCIONES PARA POSTPROCESADO #----------------------------------------------------------- from odbAccess import * from step import * from operator import itemgetter from os import chdir, path, remove, mkdir, system, listdir, getcwd def postprocess(arg1, arg2, arg3, arg4, arg5): odb = openOdb(arg1 + '.odb') 3.3.1 Descripción de los objetos Odb en Abaqus Dentro del fichero postprocess.py se encuentran variables que hacen uso de los objetos fieldOutputs y historyRegions de Abaqus, y que se describen dentro de este apartado para aclarar su uso y finalidad. Para extraer resultados de las variables de campo y de ciertas entidades del modelo (ya sea de un set ó conjunto) definidas en un step de Abaqus, se hace uso de los objetos fieldOutputs y historyRegions, respectivamente; ambos se encuentran dentro del contenedor steps (ver la figura 3.4). El objeto fieldOutputs, ilustrado en la figura 3.5 contiene los datos de campo para una determinada variable de salida del modelo de elementos finitos, y cuyo formato es el siguiente: session.odbs[name].steps[name].frames[i].fieldOutputs[name] donde el argumento name en fieldOutputs[name] hace referencia al nombre de la variable de campo que puede ser p.ej., tensiones ( 'S11' ), deformaciones ( 'E11' ), componente de desplazamiento ( 'U1' ), etc. Para acceder a las entidades y valores asociados con dicha variable, se hace uso del objeto FieldValue que representa los datos de la variable de campo en un punto del modelo, y que sigue el siguiente formato: session.odbs[name].steps[name].frames[i].fieldOutputs[name].values[i] El objeto FieldValue contiene los siguientes miembros (por brevedad se introducen aquellos usados en el fichero postprocess.py): •integrationPoint: un número entero que especifica el punto de integración dentro del elemento. • elementLabel: un número entero que especifica la etiqueta (o numeración) del elemento asociado con la variable de campo. • data: es la entidad del tipo Float que indica el valor de la variable de campo asociada con el objeto FieldOutputs. 3.3 Postprocesamiento de Resultados: el Fichero postprocess.py 37 Figura 3.3 Los distintos repositorios del objeto de resultados (odb) en Abaqus. Figura 3.4 El objeto fieldOutput. Por otro lado, el objeto historyRegions permite acceder a las entidades de resultados del historial de salida (History Output) en un step de Abaqus. Dependiendo del tipo de resultado esperado, el objeto historyRegions se puede definir para uno de los siguientes: •un nodo •un punto de integración •una región •el modelo completo Para acceder a los elementos del historyRegions, se usa un método de repositorios llamado keys(), cuya descripción se compila en la tabla 3.1 junto a otros métodos. Un ejemplo del método keys() se ilustra con los siguientes comandos: from abaqus import * from abaqusConstants import * from odbAccess import * 38 Capítulo 3. Implementación SLA + Python Tabla 3.1 Algunos métodos de indagación de las entidades en los objetos Odb. Método Descripción keys() Devuelve una lista con las claves (entidades) en el repositorio values() Devuelve una lista de objetos en el repositorio items() Devuelve una lista del conjunto [claves, lista] del repositorio myodb=openOdb('singlefibre_biax3_500steps_10_4ptos_danar.odb') myodb.steps['Step-1'].historyRegions.keys() y que devuelve: ['Element MATRIZ-1.152 Int Point 1','Element MATRIZ-1.152 Int Point 2', 'Element MATRIZ-1.152 Int Point 3','Element MATRIZ-1.152 Int Point 4', 'Element MATRIZ-1.9859 Int Point 1','Element MATRIZ-1.9859 Int Point 2', 'Element MATRIZ-1.9859 Int Point 3','Element MATRIZ-1.9859 Int Point 4', 'Node INTERFACE-1.12','Node INTERFACE-1.14'] Es decir, keys() retorna un array con los elementos del objeto historyRegions en orden alfanumérico y separados por comas. Para acceder a dichos elementos se usan índices numéricos, donde [0] y [-1] representan el primer y último elemento del array, respectivamente. Finalmente, con el objeto historyOutputs (figura 3.6) se permite extraer información del historial de salida en un punto del modelo para la variable de campo especificada. Se accede a los valores numéricos de una variable en un step lineal estático mediante el miembro data[1][1] , y que significan (frameValue, value) donde frameValue es tanto tiempo, frecuencia, o modo, y value es el valor de la variable especificada en frameValue. Figura 3.5 El objeto historyOutputs. Un ejemplo que explica cómo pedir la componente de tensión S11 en un punto de integración de un elemento del modelo es como sigue: RF1 = path_HisReg.historyOutputs['S11'].data[1][1] y que se lee: del objeto historyRegions path_HisReg , captura los datos de la componente de tensión S11 y asignárselo a la variable RF1. 3.3.2 Declaración de las variables de campo e historial de regiones en postprocess.py Después de que postprocess.py abre el archivo .odb de la iteración −n del SLA, se crean las variables que guardan los datos del objeto factor (SDV2) y número del elemento (SDV3); key_HisReg guarda el array con los elementos del historial de regiones del modelo Abaqus. Dichas variables se muestran en el siguiente 3.3 Postprocesamiento de Resultados: el Fichero postprocess.py 39 segmento de código 3.5: Código 3.5 Variables de campo e historial de regiones. key_step = odb.steps.keys() path_SDV2 = odb.steps[key_step[0]].frames[-1].fieldOutputs['SDV2']. values path_SDV3 = odb.steps[key_step[0]].frames[-1].fieldOutputs['SDV3']. values key_HisReg = odb.steps[key_step[0]].historyRegions.keys() path_HisReg = odb.steps[key_step[0]].historyRegions[key_HisReg[0]] path_HisRegEl2 = odb.steps[key_step[0]].historyRegions[key_HisReg[4]] path_HisRegp = odb.steps[key_step[0]].historyRegions[key_HisReg[8]] path_HisRegn = odb.steps[key_step[0]].historyRegions[key_HisReg[-1]] Donde: •key_step : variable que almacena el objeto Step de Abaqus; se corresponde al nombre del paso lineal estático que aplica las cargas externas al modelo y que por defecto es llamado 'Step-1' , y el índice que lo selecciona en Python es [0]. •path_SDV2: almacena los objetos de la variable dependiente del estado, SDV2 = factor. •path_SDV3: almacena los objetos de la variable dependiente del estado, SDV3 = NOEL. 3.3.3 Creación de tablas y datos a imprimir Luego de crear las variables de campo y del historial de regiones, la información de la variable SDV2 es colocada en un array llamado tabla_all , que luego se ordena de mayor a menor tomando el valor numérico de SDV2 como referente (código 3.6). En todo caso, los componentes de tabla_ord se pasan al array tabla_print. Código 3.6 Creacion de tablas y ordenamiento de datos. tabla_all=[] tabla_print=[] for kin range(len(path_SDV2)): tabla_all.append([path_SDV2[k].integrationPoint, path_SDV2[k]. elementLabel, path_SDV3[k].data, path_SDV2[k].data]) # Ordenar de mayor a menor la tabla anterior en funcion de la cuarta columna (SDV2 == factor) tabla_ord=sorted(tabla_all, key=itemgetter(3), reverse=True) # En el caso de aplicar "tolerancia" if (float(arg2) < 1): for kin range(len(tabla_ord)): if (((tabla_ord[0][3] - tabla_ord[k][3])/tabla_ord[0][3]) > float( arg2)): break tabla_print.append([tabla_ord[k][0], tabla_ord[k][1], tabla_ord[k ][2], tabla_ord[k][3]]) # En el caso de aplicar "numero de puntos de integracion a danar" else: tabla_print=tabla_ord 40 Capítulo 3. Implementación SLA + Python Finalmente se crean las variables que almacenan el valor numérico del objeto historyOutputs con el campo requerido (S11, S22, U1, U2, ..., etc.), y que ya han sido previamente definidas en el History Output de Abaqus : #--------------------------------------------------- #||||||| OBTENER RF1, RF2, U1, U2 ||||||| #--------------------------------------------------- RF1 = path_HisReg.historyOutputs['S11'].data[1][1] SEl2 = path_HisRegEl2.historyOutputs['S11'].data[1][1] U1 = path_HisRegp.historyOutputs['U1'].data[1][1] U1p = path_HisRegn.historyOutputs['U1'].data[1][1] #--------------------------------------------------- En postprocess.py se crean tres archivos de datos que recogen la información del daño y de los valores del historyOutputs en cada iteración del SLA, y son los siguientes: •datos_procesar.txt : el fichero de texto donde se graba el punto de integración dañado junto a la etiqueta del elemento. Este se guarda en el fichero archivos_procesar creado por postprocess.py. •evolucion_rotura_.txt : un fichero de texto que recoge información de la evolución de la rotura por cada iteración del SLA. •tabla_grafica.txt : almacena en forma de tabla los valores numéricos del objeto historyOutputs y el valor máximo de la variable SDV2 por cada iteración del SLA. Con los datos en este fichero se puede construir p.ej., la curva carga-desplazamiento que captura las inestabilidades post-pico propias del modelo. El grueso del código postprocess.py se publica bajo el apéndice C. En la siguiente sección, se explicarán algunos detalles de la subrutina UMAT relacionados al código SLA. 3.4 Observaciones en la Subrutina UMAT Relacionadas al SLA El primer detalle a tomar en cuenta consiste en los valores que se asignan a los parámetros nmax_pt y nmax_el , que deben ser mayor o igual al máximo número de puntos de integración por elemento (usualmente cuatro)y que el número de elementos del modelo, respectivamente. ! maximum number of integration points and elements logical lee integer nmax_pt, nmax_el parameter (nmax_pt=4, nmax_el=600000) integer dama(nmax_pt,nmax_el) common lee, dama data lee /.true./ Segundo, cada vez que se cambie el nombre del fichero de elementos finitos (el archivo .inp de Abaqus), el usuario debe escribir el nombre de éste, sin extensión, en el directorio raíz del usuario, para que la subrutina pueda leer el archivo datos_procesar.txt, como sigue: OPEN (7,FILE='/home/jpimentel/singlefibre/ +archivos_procesar/datos_procesar.txt', +STATUS='OLD',ACTION='READ',IOSTAT=error_ap) 3.5 El Fichero run-test Como se había mencionado en la sección 3.1, en este archivo se escribe el nombre del fichero de entrada Abaqus (sin extensión), ubicado dentro de la carpeta PREPARASLA del usuario en el servidor de cálculo. Este fichero, que se muestra bajo el código 3.7, es el que hace la llamada de main_abaqus.py en el servidor, y las variables a modificar por el usuario son las siguientes: 3.5 El Fichero run-test 41 •--nodes : indica el número de bandejas de cálculo en el cluster para resolver la simulación, por defecto se deja en 1. •--nodelist: indica el número de la bandeja de cálculo que manejará la simulación. •--ntasks-per-node: corresponde al número de procesadores por bandeja de cálculo. •--time : indica el tiempo máximo de la simulación en el servidor de cálculo, su formato es hh:mm:ss. •--job-name : aquí se escribe el nombre del fichero de entrada de Abaqus ubicado en el fichero PREPARASLA. Código 3.7 El fichero run-test. #!/bin/bash #SBATCH --nodes=1 #SBATCH --nodelist=nodo7 #SBATCH --ntasks-per-node=4 #SBATCH --time=60:00:00 #SBATCH --job-name=singlefibreNloff JOB_NAME=$SLURM_JOB_NAME DIRECTORIO=$JOB_NAME WORK_DIR=$HOME/$DIRECTORIO unset SLURM_GTIDS mkdir -p $WORK_DIR cd $WORK_DIR echo "El directorio de trabajo actual es `pwd`" echo "-------------------" MAQUINA=`tail $PBS_NODEFILE -n 1` echo "Su trabajo se est ejecutando en el $MAQUINA" echo "-------------------" NUM_PROCS=$SLURM_NPROCS echo "Funcionando en $NUM_PROCS procesador/es." echo "La ejecuci ha comenzado el: `date`" echo "-------------------" ssh master "cd /home/$USER/PREPARASLA; tar cfz rts.tgz *" scp master:/home/$USER/PREPARASLA/rts.tgz . tar xfz rts.tgz rm -f rts.tgz ssh master rm /home/$USER/PREPARASLA/rts.tgz source /opt/intel/oneapi/setvars.sh export LD_PRELOAD=/opt/intel/oneapi/compiler/2021.1.1/linux/compiler/lib/ intel64_lin/libirc.so #python maincluster.py >salida.out /opt/abaqus/Commands/abq6141 cae noGui=main_abaqus.py > salida.out echo "-------------------" echo "FIN de la ejecuci el: `date`" 48 Capítulo 4. Modelos Unifibra Lo que de manera tácita demuestran los resultados en dichas figuras, es que la topología de discretización de la matriz sí influye en el camino de grieta, p.ej., la malla estructurada del tipo "mapped" ((b) en las figuras 4.8, 4.9) produce un camino preferencial de grieta inclinado en 45 grados con la dirección de la carga σ∞ , y dicho camino no es congruente con la evidencia de las fractografías para este tipo de carga (ver a [ 43 ] y el capítulo 1 de este documento). En contraste, la evolución de la variable de daño, SDV1, para las mallas (a) y (c) produce un camino de grieta que tiende a formar un ángulo recto con la dirección de la carga, reproduciendo entonces el daño en las fractografías del cap. 1. Recomendaciones •Se deben evitar las direcciones preferenciales de los elementos en el mallado de la matriz. • Se recomienda usar un algoritmo de mallado que genere elementos cuyos bordes no formen un patrón, y cuya organización sea lo más libre posible, p.ej., (a) en la figura 4.7. 4.1 Modelo de Matriz Finita 49 (a) Mallado libre (b) Estructurada 1 (c) Estructurada 2 Figura 4.8 Estudios de objetividad de la malla de la matriz; usando el criterio de fallo en tensión principal máxima σ1, con St=80 MPa. 50 Capítulo 4. Modelos Unifibra Paso (a) Mallado libre (b) Estructurada 1 (c) Estructurada 2 (1) (2) (3) (4) (5) Figura 4.9 Estudios de objetividad de la malla de la matriz; usando el criterio de la componente dilatacional de la energía de deformación, Uv, con Ucrit =1J/m3. 4.1 Modelo de Matriz Finita 51 4.1.5 Efectos de los distintos criterios de fallo empleados sobre el camino de grieta en la matriz finita Bajo este apartado, se documentarán los resultados numéricos para los distintos criterios de fallo de la matriz en los casos en que fue posible, y dicho esto debido a que para las propiedades mecánicas de referencia en la interfase LEBIM ( ¯ σc =90MPa, ¯ GIc =2 J/m2 , λ =0.25, kt/kn=0.25 y h=0.01 µm ) y la resistencia a tracción Sut =80 MPa del material de la matriz, no se produce el fallo de la interfase cuando se usa la tensión principal máxima y el criterio de von Mises (ambos bajo control de carga en tensiones), debido a que se emplea una resistencia de la matriz con menor valor que la tensión critica de la interfase (i.e., Sut < ¯ σc ); el daño de la matriz inicia en la periferia externa de la interfase, extendiéndose hacia los extremos y en dirección normal a la aplicación de la carga (ver figura 4.10). (a) Inicio del daño en la matriz (b) Daño final Figura 4.10 Fallo exclusivo de la matriz con el criterio de la tensión principal máxima, Sut =80 MPa y control de carga en tensiones. La figura 4.11 muestra la evolución de la variable de daño usando tres criterios de fallo en la matriz y control de carga en tensión. El criterio de Mohr-Coulomb así como el de la energía de deformación en dilatación empiezan el fallo por la interfase que luego progresa hacia la matriz, con la salvedad de que Mohr-Coulomb presenta ramificaciones dispersas del daño en contraste con el criterio energético que presenta un camino de grieta menos abrupto. Usando un control de carga en desplazamiento, la figura 4.12 presenta la variable de daño para los pasos finales de la progresión de grieta y para los distintos criterios de fallo empleados en la matriz. Tanto la tensión principal máxima, von Mises y la energía de deformación total producen un fallo dominante de la matriz; tanto von Mises como la energía total de deformación producen un camino de grieta que no toca la interfase y que son similares entre si, indicando que cuando se usa la energía total de deformación la componente de distorsión de esta última domina el daño de la matriz. Por otra parte no es menos importante notar que por el tipo de condiciones de contorno aplicadas, la morfología del camino de grieta es distinta de los modelos con el control de carga en tensiones (ya que se aplica tracción sin desplazamiento relativo entre los nodos de una cara). Finalmente, en la figura 4.13 se grafican las curvas de equilibrio para los modelos (b) y (c) de la fig. 4.11; no se produce una notable diferencia entre las tensiones críticas del despegue fibra-matriz para aquellos criterios de fallo de la matriz que garanticen la iniciación de grieta en la interfase LEBIM. 52 Capítulo 4. Modelos Unifibra Paso (a) mytpe=2, criterio σ1(b) mtype=3, Mohr-Coulomb (c) mtype=5, energético Uvol (1) (2) (3) (4) (5) Figura 4.11 Resultados de la evolución de la variable de estado SDV1, para distintos criterios de fallo de la matriz; el control de carga se realiza bajo tensión σ∞ x , con las propiedades mecánicas Sut =80 MPa, Suc=120 MPa, y Ucr=1J/m3. 4.1 Modelo de Matriz Finita 53 Criterio de fallo SDV1 Observaciones Tensión principal máxima, Sut =80 MPa Fallo dominante de la matriz Mohr-Coulomb, Sut =80 MPa ySuc=120 MPa Fallo de interfase y cambio hacia la matriz von Mises, Sut =80 MPa Fallo dominante de matriz, sin tocar la interfase Componente dilatacional, Uv;Ucr=1J/m3Fallo de interfase seguido de la matriz Energía total, U ; Ucr =1 J/m3Daño similar al de von Mises Figura 4.12 Variable de daño SDV1, del modelo unifibra para los distintos criterios de fallo evaluados; control de carga realizado en desplazamiento, u=0.01µm. 54 Capítulo 4. Modelos Unifibra Figura 4.13 Curvas de equilibrio para distintos criterios de fallo, en el modelo de matriz finita. 4.1.6 Efectos de las propiedades mecánicas interfase LEBIM/matriz Usando un control de carga en tensiones σ∞ x y el criterio de la componente dilatacional de la energía de deformación, la figura 4.14 muestra resultados de la evolución del daño para los siguientes valores de energía critica: 1, 2, y 5 J/m3 . A mayor capacidad de absorción de energía de deformación en la matriz, se alcanza un mayor ángulo de despegue de la interfase θd . El mismo fenómeno ocurre bajo un control de carga en desplazamientos, y cuando se usa el criterio de fallo de la energía total de deformación Ut ((b) y (c) en la figura 4.15), con la salvedad de que para el valor crítico de 1 J/m3 , el daño es dominado por la matriz y sin tocar la interfase LEBIM. Dejando constante en 1 J/m3 la energía crítica de la matriz, y usando el criterio de la componente de dilatación Uv , se modifican las propiedades de tenacidad y resistencia de la interfase, tal como se muestra en la figura 4.16. La morfología de la grieta para ¯ GIc = 2J/m2 , ¯ σc =60MPa es similar al caso (a) en 4.14. En contraste, para valores más elevados del índice de liberación de energía crítico ¯ GIc =10 J/m2 se produce un camino de grieta asimétrico y discontinuo entre la interfase y la matriz ((b) y (c) de 4.16). Dichos resultados de la evolución del daño (SDV1) demuestran que la capacidad de absorción de energía de deformación de la matriz Ucr tiene un mayor efecto en el semi-ángulo de despegue de la interfase θd , que la variación de las propiedades mecánicas en la interfase LEBIM. En las curvas de equilibrio σ∞ x−CMOD ilustradas en la figura 4.17, se demuestra que (tanto para el control en tensión así como en desplazamiento) los valores críticos de energía de deformación de la matriz no influyen en la tensión crítica remota σ∞ x necesaria para el despegue de grieta en la interfase; aunque dichos valores sí influyen en el inicio del comportamiento "snap-back" post-fisuración, ocurriendo esto para aperturas de grieta de interfase (CMOD) menores con valores más bajos de Ucr. Finalmente, y usando el criterio de la componente volumétrica de la energía de deformación (debido a su capacidad de iniciar el daño en la interfase, ver fig. 4.14), se ilustra en la figura 4.18 distintas curvas de equilibrio obtenidas para cuatro pares distintos de propiedades en la interfase LEBIM. A mayor índice de liberación de energía crítico en modo I ( ¯ GIc =10 J/m2 ), mayor es la tensión crítica necesaria para el despegue de grieta y mayor es la apertura radial de grieta previo al fenómeno "snap-back". De dichas curvas, cabe notar la presencia de inestabilidades en forma de dientes de sierra que se producen post-fisuración de la interfase, y que el algoritmo SLA ha sido capaz de capturar. 4.1 Modelo de Matriz Finita 55 Paso (a) Ucr=1J/m3(b) Ucr=2J/m3(c) Ucr=5J/m3 (1) (2) (3) (4) (5) Figura 4.14 Resultados de la evolución de la variable de estado SDV1, bajo el criterio de la componente dilatacional de energía de deformación, y para distintos valores de energía crítica de la matriz, (Ucr=1, 2 y 5J/m3); el control de carga se realiza bajo tensión σ∞ x. 56 Capítulo 4. Modelos Unifibra Paso (a) Ucr=1J/m3(b) Ucr=2J/m3(c) Ucr=5J/m3 (1) (2) (3) (4) (5) Figura 4.15 Resultados de la evolución de la variable de daño SDV1, bajo el criterio de la energía total de deformación, y para distintos valores de energía crítica de la matriz, ( Ucr =1, 2 y 5 J/m3 ); el control de carga se realiza en desplazamientos. 4.1 Modelo de Matriz Finita 57 Paso (a): ¯ GIc=2J/m2,¯ σc=60MPa (b): ¯ GIc=10J/m2,¯ σc=90MPa (c): ¯ GIc=10J/m2,¯ σc=60MPa (1) (2) (3) (4) (5) Figura 4.16 Evolución de la variable de estado SDV1, para distintos valores de propiedades de la interfase LEBIM; el control de carga se realiza bajo tensión σ∞ x. 64 Capítulo 4. Modelos Unifibra Tabla 4.5 Resultados de las tensiones remotas críticas, σ∞ c, para el despegue de la interfase. Número de elementos de interfase σ∞∥errorabs(σ∞)∥,% 1 el./grado 59.947 MPa — 2 el./grado 59.945 MPa 0.0033 Figura 4.24 Curvas de equilibrio σ∞- CMOD, para distinto número de elementos en la interfase LEBIM. 4.2.3 Efectos del tipo de control de carga Para estudiar los efectos del tipo de control de carga realizado sobre la matriz, se decide representar las curvas de equilibrio ( σ∞ - CMOD) obtenidas por el algoritmo SLA, debido a que con estas curvas se visualizan las diferencias entre las tensiones críticas de despegue de la interfase que una tensión o desplazamiento remoto aplicado en la matriz pudieran provocar. En la figura 4.25 se grafican dichas curvas, demostrando que la tensión crítica remota σ∞ x , necesaria para el crecimiento inestable de la grieta de interfase permanece invariante ante los dos distintos controles de carga: la tension y desplazamiento remoto ( σ∞ y u∞ , respectivamente). Esto confirma la hipótesis del principio de Saint-Venant realizada en la sección 4.2. El gráfico (b) de 4.25 demuestra el comportamiento post-fisuración de la interfase, donde el "snap-back" se mantiene invariante bajo un determinado criterio de fallo de la matriz. 4.2 Modelo de Matriz Infinita 65 (a) (b) Inestabilidades post-fractura Figura 4.25 Curvas σ∞- CMOD para distintos controles de carga en la matriz infinita. 4.2.4 Efectos de los distintos criterios de fallo de la matriz La figura 4.26 ilustra los estados avanzados de la variable de daño para los distintos criterios de fallo de la matriz. De manera similar a lo que ocurre en el modelo unifibra de matriz finita, para dos de las hipótesis de fallo en tensiones, la tensión principal máxima y de von-Mises producen un daño iniciado y dominado por la matriz (usando una resistencia Sut =80MPa [ 57 ]). En contraste, el criterio de Mohr-Coulomb y los criterios energéticos Uv y Ut permiten el inicio del fallo en la interfase, que luego progresa y se extiende en la matriz siguiendo una dirección normal con la carga remota aplicada. Las curvas σ∞ - CMOD graficadas en la figura 4.27 demuestran que la tensión crítica del despegue de la interfase no es alterada por aquellos criterios de fallo que garanticen la iniciación del daño en la interfase. En cambio, el criterio de fallo de von Mises produce un mecanismo de fallo dominado por el daño en la matriz, sin nucleación de la grieta de interfase (esto se nota por la línea roja en la figura 4.27). 66 Capítulo 4. Modelos Unifibra Criterio de fallo SDV1 Observaciones Tensión principal máxima, Sut =80 MPa Fallo dominante de la matriz Mohr-Coulomb, Sut =80 MPa ySuc=120 MPa Fallo de interfase y matriz von Mises, Sut =80 MPa Fallo dominante de matriz Componente dilatacional, Uv;Ucr=1J/m3 Energía total, Ut; Ucr=2J/m3 Figura 4.26 Variable de daño SDV1, del modelo unifibra para los distintos criterios de fallo evaluados; propiedades de la interfase LEBIM: ¯ GIc=2J/m2,¯ σc=90 MPa, λ=0.25, kt/kn=0.25, h=0.01µm. 4.2 Modelo de Matriz Infinita 67 Figura 4.27 Curvas de equilibrio σ∞- CMOD, para distintos criterios de fallo de la matriz. 4.2.5 Efecto de las propiedades mecánicas de la interfase LEBIM/matriz Manteniendo constante la energía crítica de la matriz en 1 J/m3 y usando el criterio de la componente volumétrica de la energía de deformación, se procede a variar las propiedades de resistencia ( ¯ σc ) y el índice de liberación energético ( ¯ GIc ) del LEBIM; dichos efectos se cuantifican en las curvas σ∞ x - CMOD de la figura 4.28. Se concluye que para un valor fijo de ¯ GIc , incrementar la resistencia de la interfase ¯ σc resulta en una mayor rigidez global del modelo. Usando los criterios de fallo energéticos y propiedades estándar de la interfase, en la figura 4.29 se ilustra la evolución de la variable de daño para dos distintos valores de la energía crítica de la matriz: 1 J/m3 y 5 J/m3 . El ángulo de despegue de la interfase θd incrementa cuando la matriz posee una mayor capacidad de absorción de energía, similar a lo que sucede en los modelos de matriz finita (ver figura 4.14). Un aspecto del modelo que queda mejor expuesto con la variación del valor de la energía crítica de la matriz, es que las tensiones remotas aplicadas que promueven el daño de la matriz incrementan su valor cuando se eleva Ucr , como se muestra en la figura 4.30. Figura 4.28 Curvas de equilibrio σ∞ - CMOD, para distintos valores de las propiedades mecánicas de la interfase LEBIM. 68 Capítulo 4. Modelos Unifibra Paso SDV1, Uvol:Ucr=1J/m3SDV1, Uvol:Ucr=5J/m3 1 2 3 SDV1, criterio Utotal:Ucr=2J/m3↓SDV1, criterio Utotal :Ucr=5J/m3↓ (Final) Figura 4.29 Variable de daño SDV1, del modelo unifibra para distintos valores de la energía crítica de la matriz; propiedades de la interfase LEBIM: ¯ GIc =2 J/m2 , ¯ σc =90 MPa, λ =0.25, kt/kn =0.25, h=0.01µm. 4.2 Modelo de Matriz Infinita 69 (a) Dilatación, Uvol (b) Energía total, Ut Figura 4.30 Curvas de equilibrio σ∞ x - CMOD, para distintos valores de la energía crítica de la matriz Ucr , usando (a) la componente de dilatación y (b) la energía total de deformación. 5 Modelos Multifibra: un Análisis Numérico del Camino de Grieta 5.1 Introducción El llamado fallo inter-fibra/matriz que aparece en los laminados bajo cargas transversales y que se ha introducido en la sección 1.2, ha sido objeto de estudios previos [ 13 , 14 , 45 , 65 ]. De acuerdo con Talreja y colaboradores [ 4 , 5 ], el inicio de este fallo puede estar asociado con la cavitación de la matriz o la aparición de un despegue en la interfase fibra-matriz. En estos estudios previos [ 45 , 42 ], los autores hicieron la hipótesis de que el fallo empieza con la aparición de pequeños despegues en la interfase, y este ha sido el supuesto de base para realizar los análisis de la morfología de grieta en el anterior capítulo de los modelos unifibra. El uso de los modelos multifibra ha sido empleado por distintos autores para estudiar el comportamiento de laminados en material compuesto [ 3 , 21 , 49 , 69 ]. En el presente capítulo se estudian varios modelos multifibra a nivel micromecánico, para concentrarse en la evolución del mecanismo de daño desde la nucleación y el crecimiento de la grieta de interfase fibra-matriz, hasta su bifuración y crecimiento a través de la matriz. Particularmente, el objetivo aquí en este capítulo es presentar la evolución del daño y las inestabilidades capturadas por el código SLA en la lámina de 90 ◦ de laminados [ 0,90n ] s , modelando toda una sección del laminado. A pesar de que en laminados con alto volumen de fibras Vf , implica una gran cercanía entre fibras (ver figura 5.1) y que por ende la bifurcación de la grieta de interfase hacia la matriz no ocurre en todas las fibras, este estudio dejará de lado dicha observación, y empleará mayormente criterios de fallo energéticos de la matriz para evaluar la nucleación y propagacion del daño a través de esta última. Figura 5.1 Grieta transversal en un espécimen de fibra de carbono-epxy bajo carga transversal [37]. 71 72 Capítulo 5. Modelos Multifibra: un Análisis Numérico del Camino de Grieta 5.2 Descripción de los Modelos Los modelos analizados han sido desarrollados en Abaqus/CAE, representando laminados de compuesto en tres capas del tipo [ 0,90n ] s , donde las capas principales están claramente distinguidas como se muestra en las figuras 5.2 y 5.3. La lámina de 90 ◦ se compone de una matriz con una distribución aleatoria de fibras, mientras que las capas superior e inferior de las láminas de 0 ◦ son modeladas como un material lineal elástico ortótropo. (a) modelo de 40µm, izquierda (I) (b) modelo de 40µm, derecha (D) Figura 5.2 Modelos multifibra, lámina de espesor 40µm. 5.2 Descripción de los Modelos 73 Figura 5.3 Modelo multifibra, lámina de 90◦con 60µm de espesor. Tabla 5.1 Dimensiones de los modelos multifibra. Modelo Ancho, w Alto, h Espesor lámina 90◦ 40I 100µm 190µm 40µm 40D 100µm 190µm 40µm 60D 100µm 210µm 60µm Las cotas dimensionales de los modelos se introducen en la tabla 5.1. Se deja notar que los modelos (a) y (b) en la figura 5.2 se distinguen no tanto por la dirección de las condiciones de contorno en desplazamiento impuestas en el laminado, pero más claramente por la distinta distribución de las fibras en la lámina de 90◦ que se muestra en la figura 5.4; esto consigue estudiar la influencia de distinta distribución de fibras sobre los posibles caminos de grieta. Para los tres modelos en 5.1, se han aplicado condiciones de contorno de simetría, y se asume la hipótesis de deformación plana en los modelos numéricos. (a) (b) Figura 5.4 Distribución de fibras para los modelos 40I (a) y 40D (b). En cada modelo la matriz y la interfase han sido discretizadas con elementos cuadriláteros de orden lineal con integración de Gauss con 4 puntos (se corresponde al elemento CPE4 de Abaqus). Las fibras se asumen circulares e iguales con un diámetro de 5.2µm, mientras que la interfase se modela como un anillo circular de muy fino espesor h=0.01µm. El mallado se ha elaborado para tener el mismo número de elementos en la periferia de la fibra y en la interfase, como muestra el detalle en la figura 5.5. Respecto de la matriz epoxy en la lámina de 90 ◦ , esta se modela con una ley constitutiva lineal elástica isótropa, con un módulo de Young Ede 2.79 GPa y la relación de Poisson νcon 0.33. La discretización de 80 Capítulo 5. Modelos Multifibra: un Análisis Numérico del Camino de Grieta Evol. ↓damage,Ucr =1J/m3damage,Ucr =2J/m3damage,Ucr =5J/m3 (1) (2) (3) (4) (5) Figura 5.11 Desarrollo de la variable de daño, modelo 40D. Influencia del valor crítico de la componente dilatacional de la energía de deformación. 5.3 Resultados 81 (a) εxvs. paso SLA (b) σx-εx,total (c) σx-εx(daño) Figura 5.12 Curvas tensión-deformación en la lámina de 90 grados, modelo 40D, y para distintos valores de la energía critica de la matriz. 82 Capítulo 5. Modelos Multifibra: un Análisis Numérico del Camino de Grieta 5.3.3 Influencia del número de puntos de integración a dañar por paso del SLA: modelo 40D Manteniendo fija en 1 J/m3 la energía crítica de la matriz, los contornos de la figura 5.13 muestran el desarrollo de la variable de daño en el modelo 40D, para dos, cuatro, y ocho puntos de integración por paso del SLA. Se nota que la morfología del camino de grieta no es afectada por este parámetro, sin embargo se concluye que el número de puntos de integración a dañar por paso actúa como un acelerador del daño y por ende, se recomienda discreción en que su valor sea acorde a la escala física y el número de elementos del modelo. La curvas en (a) de la figura 5.14 validan los comentarios anteriormente hechos, ya que el nivel de tensiones aplicadas decae más rápido por iteración del SLA para el caso de 8 puntos de integración. 5.3 Resultados 83 Paso. ↓SDV1, 2 ptos. SDV1, 4 ptos. SDV1, 8 ptos. (40) (80) (100) (120) (200) Figura 5.13 Evolución del daño, modelo multifibras 40D; efectos del número de puntos de integración a dañar por paso SLA. 84 Capítulo 5. Modelos Multifibra: un Análisis Numérico del Camino de Grieta (a) σx-paso SLA (b) σx−εx,total Figura 5.14 Efectos del número de puntos de integración a dañar por paso del SLA; modelo 40D. 5.3.4 Modelo 60D Influencia de los valores críticos de la energía de dilatación de la matriz. Usando dos valores distintos de la energía critica de dilatación para la matriz, en 1- y 5- J/m3 , en la figura 5.15 se ilustra la progresión de la grieta transversal en la lamina de 90 ◦ del modelo 60D. Ambos conjuntos de resultados se presentan usando el mismo factor de escala en desplazamientos. En ambos casos, el mecanismo inicial de daño se debe al despegue de la interfase fibra-matriz. Para la energía crítica de 5 J/m3 , el alineamiento de los despgues de interfase se realiza formando un ángulo que ya deja de ser 90 ◦ con la carga aplicada. Lo más notorio de este último caso es que se retarda el proceso de delaminación. Con 1 J/m3 , la fractura en la lámina sigue un camino lo más rectilíneo posible, en dirección ortogonal con la carga aplicada hasta alcanzar la frontera entre láminas para entonces iniciar el proceso de delaminación. De las curvas de deformación aplicada εx por paso del SLA (el gráfico (a) de 5.16), ya se produce una ligera diferencia entre las deformaciones críticas que inician el daño, siendo mayor para el caso con 5 J/m3 (2.229%, en contraste con εcrit x =1.86% para 1 J/m3 ). Sin embargo, la tensión promedio aplicada a la lámina 5.3 Resultados 85 Evol. ↓SDV1,Ucr=1J/m3SDV1,Ucr=5J/m3 (1) (2) (3) (4) (5) Figura 5.15 Progresión del daño, modelo 60D; criterio de la componente dilatacional de la energía de deformación con Ucr=1- y 5-J/m3. de 90 ◦ que inicia el despegue de la interfase se mantiene igual ( ≈ 68 MPa) respecto a los modelos 40I y 40D (ver b y c de 5.12). 86 Capítulo 5. Modelos Multifibra: un Análisis Numérico del Camino de Grieta (a) εx-Paso SLA (b) σxvs. Deformación total (c) σxvs. Deformación por daño Figura 5.16 Curvas tensión-deformacion en la lámina de 90 grados, modelo 60D, y para distintos valores de la energía crítica de la matriz. 5.3 Resultados 87 Influencia de los valores críticos de la energía total deformación de la matriz. Usando el criterio de fallo de la energía total de deformación para la matriz Ut , los contornos en la figura 5.17 muestran la progresión del daño para dos distintos valores de densidad de energía critica: 1- y 5 J/m3 . En ambos casos, el mecanismo de daño se inicia en la matriz, conjunto a la aparición de modos espurios de deformacion en algunos elementos (varios han sido removidos para clarificar los gráficos). Para la densidad de energía critica de 1 J/m3 , el camino de grieta en la matriz tiende a seguir el contorno externo de las fibras, para luego llegar a la frontera de la lámina e iniciar la delaminación (sin atravesar todo el espesor de 60 µ m). Se establece que una de las principales razones para obtener un fallo dominado por la matriz con el valor de 1 J/m3 reside en que la variable f actor (ec. 2.9) que calcula la subrutina UMAT es bastante superior a la que se obtiene para la interfase, resultando esto en un daño casi exclusivo de la matriz. Cuando se tiene una matriz con mayor capacidad de absorción de densidad de energía, cambian los mecanismos de daño (i.e, ya se producen los despegues en la interfase LEBIM) y la morfología de la grieta a través del espesor de la lámina. Sin embargo, basado en los hallazgos de Asp y Talreja [ 4 ], de que en regiones con fuerte tri-axialidad en las tensiones (i.e., cercano a las fibras) se produce el fallo frágil por cavitación de la matriz, se razona que el criterio de la densidad de energía total ( Ut ) puede arrojar resultados espurios con relación a la variable de daño en regiones ricas en fibras. Finalmente, se observa de las curvas en la figura 5.18, que este criterio tiende a subestimar el valor de la deformación εx aplicada (en relación a εx que resulta de usar la componente volumétrica de la densidad de energía), necesaria para iniciar el crecimiento inestable de grieta, por lo cual el uso de este criterio se puede limitar a regiones ricas en resina de la matriz. 88 Capítulo 5. Modelos Multifibra: un Análisis Numérico del Camino de Grieta Evol. ↓SDV1,Ucr=1J/m3SDV1,Ucr=5J/m3 (1) (2) (3) (4) (5) Figura 5.17 Evolución del daño, multifibras modelo 60D, para el criterio de la componente total de la energía de deformación, Ut. 5.3 Resultados 89 (a) εx-Paso SLA (b) Detalle εx-Paso SLA Figura 5.18 Curvas εx -Paso SLA, modelo 60D, usando el criterio de la energía total de deformación para la matriz. 96 Capítulo 6. Modelos Multiescala: Análisis Numérico del Daño en Laminados 0/90 Evol. ↓SDV1,Ucr =1J/m3SDV1,Ucr =2J/m3SDV1,Ucr =5J/m3 (1) (2) (3) (4) (5) Figura 6.7 Progresión del daño, modelo multiescala con 100 fibras. Para facilitar su visualización, los resultados de la variable SDV1 en cada fila se reportan para distintos pasos del SLA. 6.3 Resultados 97 (1) ε∞ y-SLA (2) Detalle, ε∞ y-SLA Figura 6.8 Deformación aplicada por paso del SLA, modelo de 100 fibras para distintos valores de la densidad de energía crítica. 98 Capítulo 6. Modelos Multiescala: Análisis Numérico del Daño en Laminados 0/90 Figura 6.9 Curvas σ∞ ydeformación por daño, en el modelo de 100 fibras. 6.3.2 Carga impuesta en las láminas de 0 grados, modelo de 100 fibras Aplicando un desplazamiento impuesto |U2| = 0.1 µm en las caras del laminado de 0 ◦ , los contornos en la figura 6.10 representan la evolución de la variable de daño en la celda de 90 ◦ , usando el criterio de la componente volumétrica de la densidad de energía Uvol , con tres valores distintos energía crítica: 1, 2, y 5 J/m3 . Los resultados obtenidos son muy parecidos a los presentados en la pasada figura 6.7, con lo cual se razona que el efecto de aplicar la carga en el laminado de 0 ◦ no afecta en manera cualitativa a la morfología de la grieta en la matriz, en contraste con la carga en toda la cara de un laminado [0,90n]s. Las curvas (a) y (b) de la figura 6.11 confirman los anteriores comentarios acerca de la evolución del daño en el modelo con la celda de 100 fibras. Para un determinado valor de la energía crítica, la deformación que inicia el mecanismo de fallo no cambia respecto del modelo multiescala que aplica la carga en todo el laminado [0,90n]s. Finalmente, la figura 6.12 presenta curvas de equilibrio σ∞ y−ε (daño) (introducidas en la sección 5.3.2). Como es sólo el laminado de 0 grados quien soporta toda la carga aplicada y en virtud de la orientación de las fibras en dicha capa, las tensiones que inician el mecanismo de fallo para los tres niveles de energía crítica es bastante superior (por un factor mayor de diez) respecto de aquellas reportadas por el modelo que aplica la carga en toda la cara del laminado (volver a ver la figura 6.9). Esto es beneficioso hasta cierto nivel, ya que implica la capacidad de las láminas de 0 ◦ de arrestar la progresión del daño, desde el punto de vista de las tensiones aplicadas. 6.3 Resultados 99 Evol. ↓SDV1,Ucr =1J/m3SDV1,Ucr =2J/m3SDV1,Ucr =5J/m3 (1) (2) (3) (4) (5) Figura 6.10 Progresión del daño, modelo multiescala bajo carga en el laminado de 0◦. 100 Capítulo 6. Modelos Multiescala: Análisis Numérico del Daño en Laminados 0/90 (a) ε∞ y-paso SLA (b) detalle ε∞ y-paso SLA Figura 6.11 Deformación aplicada por paso del SLA, modelo de 100 fibras bajo carga en laminado de 0 ◦ y ante distintos valores de la densidad de energía de dilatación. crítica. 6.3 Resultados 101 (a) σ∞ y−ε(daño) (b) detalle, σ∞ y−ε(daño) Figura 6.12 Curvas σ∞ y−ε(daño), modelo de 100 fibras bajo carga en laminado de 0◦. 6.3.3 Modelo de 160 fibras: carga en toda la cara del laminado 0/90 Usando el criterio de la componente volumétrica de la densidad de energía Uvol , los contornos de la figura 6.13 representan la evolución de la variable de daño en la celda del modelo con 160 fibras (FIB160). Para Ucr =1J/m3 , la grieta se propaga en la matriz produciendo el daño por los despegues de interfase que se encuentran en el frente de la grieta. Se puede apreciar la bifurcación que toma el camino del daño en la matriz antes de alcanzar la frontera con la lámina de 0 ◦ para iniciar el proceso de delaminación. Con Ucr =5J/m3 , el mecanismo de fallo se produce en su mayoría por los despegues de la interfase fibra-matriz, cuyo alineamiento deja de tomar un camino rectilíneo a lo largo de la celda. El proceso de delaminación se retarda, notándose que se producen despegues de interfase en dirección paralela con la carga aplicada y cercano a la frontera con la lámina de 0◦. Las curvas (1) y (2) en la figura 6.14 representan la evolución de la deformación aplicada al laminado por paso del SLA, y las curvas tensión-deformación por daño, respectivamente. Se nota de las curvas εy [%] - paso SLA (1), que la magnitud de la deformación remota aplicada al laminado incrementa con el avance del daño (en particular para 1 J/m3 ). Sin embargo, respecto a las cotas en el eje vertical de ambos gráficos 102 Capítulo 6. Modelos Multiescala: Análisis Numérico del Daño en Laminados 0/90 es recomendable que se limiten a valores de la deformacion (y tensión) transversal al fallo ε2T , para ser coherente con las resistencias de una lámina (i.e., tabla 1 en [ 57 ]). Las deformaciones remotas que inician el primer fallo son (ligeramente) menores que las reportadas para el modelo de 100 fibras (ver figura 6.8). Evol. Ucr =1J/m3Ucr =5J/m3 (1) (2) (3) (4) (5) Figura 6.13 Progresión de la variable de daño en la celda de 160 fibras, y para dos distintos valores de la energía crítica de la matriz. 6.3 Resultados 103 (1) εy[%] - paso SLA (2) σ∞ y-ε(daño) Figura 6.14 Curvas tensión-deformaciones en el modelo multiescala de 160 fibras. 6.3.4 Modelo de 160 fibras: carga en el laminado de 0 grados La figura 6.15 muestra la evolución del daño en la celda del modelo FIB160, usando el criterio de la componente volumétrica de la densidad de energía de deformación. En ambos casos (para Ucr =1 J/m3 y Ucr =5 J/m3 ), la morfología de grieta tiene gran similaridad con la representada en la pasada figura 6.13, lo cual implica que el camino prospectivo de la grieta no se altera al modificar las condiciones de contorno en carga de un laminado [ 0,90n ] s . Sólo el cambio en los valores de la energía crítica de la matriz produce cambios en los mecanismos de fallo en la escala micromecánica. Por otro lado, las deformaciones remotas ε∞ y aplicadas a las láminas de 0 ◦ presentadas en (1) de la figura 6.16 y que inician el primer fallo (0.06% para Ucr =1 J/m3 y 0.126% para Ucr =1 J/m3 ), se corresponden con las presentadas en (1) de la figura 6.14, indicando con esto que el proceso de daño es guiado por la deformación aplicada al laminado. Finalmente, como son las láminas de 0 ◦ que en este caso reciben la carga aplicada en dirección de las fibras, no es sorpresa observar que las tensiones iniciadoras del primer fallo (2 y 3 en 6.16) son mayores ( ≈ 100 MPa para Ucr =1 J/m3 y 220 MPa para Ucr =5 J/m3 ) en comparación con las que captura un elemento en la lámina de 90 ◦ adyacente con la aplicación de las cargas en el laminado (ver 104 Capítulo 6. Modelos Multiescala: Análisis Numérico del Daño en Laminados 0/90 Evol. SDV1,Ucr =1J/m3SDV1,Ucr =5J/m3 (1) (2) (3) (4) (5) Figura 6.15 Progresión del daño en el modelo multiescala con 160 fibras, bajo carga en el laminado de 0 ◦ , y para dos valores distintos de la densidad de energía crítica Ucr. p.ej., (2) en la figura 6.14). 6.3 Resultados 105 (1) εy[%] - paso SLA (2) σ∞ y-ε(daño) (3) Detalle, σ∞ y-ε(daño) Figura 6.16 Curvas tensión-deformaciones en el modelo multiescala de 160 fibras, bajo carga en el laminado de 0 grados. 112 Apéndice A. Subrutina UMAT ! Inicializacion variables DDS=0.d0 psig=0.d0 STRESS=0.d0 Pi=ACOS(-1.d0) Gi=0.d0 Gii=0.d0 ! Explicación de time(1): ! TIME(1): Current value of step time. ! TIME(2): Current value of total time. ! !READ THE MATERIAL PROPERTIES FROM THE INPUT FILE mtype=PROPS(1) if (mtype.eq.1) then sigmacb=PROPS(2) GIcb=PROPS(3) lambda=PROPS(4) ktkn=PROPS(5) h=PROPS(6) else if (mtype.gt.1) then E=PROPS(2) XNUE=PROPS(3) EBULK3=E/(1-2*XNUE) EG2=E/(1+XNUE) ELAM=(EBULK3-EG2)/3 !c JMP (Sept-2022) ************************************************** !c[ The next lines define the failure criteria selected for the ] !c[ composite laminae matrix material as a function of the input ] !c[ number (an integer) given to the epoxy resin material card ] !c[ (the first row of the user defined properties in the abaqus ] !c[ user material type). ] !c ****************************************************************** !c -------------------mtype choice:------------------------- !c 2=Maximum principal stress failure criteria ] !c 3=Mohr-Coulomb brittle failure criteria ] !c 4=von-Mises yield failure criteria ] !c 5=Critical volumetric strain energy density criteria ] !c 6=Critical total strain energy density criteria ] !c --------------------------------------------------------- if (mtype.eq.2) then ! Maximum principal stress failure criteria: sigmayt=PROPS(4) else if (mtype.eq.3) then ! Mohr-Coulomb failure criteria: sigmayt=PROPS(4) Suc=PROPS(5) else if (mtype.eq.4) then ! von-Mises yield failure criteria: sigmayt=PROPS(4) else if (mtype.eq.5) then ! Critical volumetric strain energy density criteria: Uc=PROPS(4) else if (mtype.eq.6) then ! Critical total strain energy density criteria: Uc=PROPS(4) 113 endif endif ! ********************************************************* ! Leer fichero archivos_procesar.txt la primera vez para ! construir el array dama() call MutexLock(1) ! Abrimos Mutex if (lee) then dama(:,:) = 0 OPEN (7,FILE='/home/jpimentel/singlefibre/ +archivos_procesar/datos_procesar.txt', +STATUS='OLD',ACTION='READ',IOSTAT=error_ap) error_lec=0 DO WHILE ((error_ap.eq.0).and.(error_lec.eq.0)) READ (7,*,IOSTAT=error_lec) PT, ELEM IF(error_lec.eq.0) THEN dama(pt,elem)=1 ENDIF ENDDO CLOSE(7) lee=.false. end if damage = float(dama(npt,noel)) call MutexUnLock(1) ! Cerramos Mutex ! ***************************************************** ! LEBIM Behaviour mtype=1 ! ***************************************************** ! LTM (20/04/21) if (mtype.eq.1) then ! LTM (20/04/21) !cccccccccccccccccccccccccccccccccccccccccccccccc ! LEBIM interface option ! ***************************************************** ! Se entra en LEBIM sin daño ! Se entra en el bucle principal IF (damage.lt.1) THEN ! Rigideces del resorte Knn=h*(sigmacb**2)/(2*GIcb) Kss=Knn/1d18 Ktt=Knn*ktkn K33=Knn/1d18 ! Actualización de la matriz de elasticidad DDS(1,1)=Knn DDS(2,2)=Kss DDS(3,3)=K33 DDS(4,4)=Ktt ! DSS(1,1) es la radial y DSS(2,2) es en la dirreción teta ! Determination del incremento del tensor de tensión DO k=1,4 DSTRESS(k)=DDS(k,k)*DSTRAN(k) ENDDO ! Implementación del tensor de tensión DO k=1,4 STRESS(k)=STRESS(k)+DSTRESS(k) 114 Apéndice A. Subrutina UMAT ENDDO ! Cálculo de Gi, Gii, psig Gi=h*(STRESS(1))**2.d0/(2.d0*Knn) Gii=h*(STRESS(4))**2.d0/(2.d0*Ktt) psig=datan2(STRESS(4)*dsqrt(Knn/Ktt),STRESS(1)) !cccccccccccccccccccccccccccccccccccccccccccccccccccc ! Modificaciones LTM (16/03/15) IF (STRESS(1).GT.0.D0)THEN Gtot=Gi+Gii ELSE Gtot=Gii ENDIF Gc=GIcb*(1.d0+(dtan(psig*(1.d0-lambda)))**2.d0) psiGcrit=pi/(2.d0*(1.d0-lambda)) IF(abs(psig).ge.psiGcrit) THEN Gc=GIcb*1.d8 ENDIF factor = dsqrt(Gtot/Gc) IF (factor.ge.1.d0) THEN damage=1 else damage=0 ENDIF ! Modificaciones LTM (16/03/15) !cccccccccccccccccccccccccccccccccccccccccccccccccccc ! Fin del bucle principal ELSE IF (damage.eq.1) THEN ! ************************************************ ! Actualización de la matriz de elasticidad, antes no tenía encuenta que si estaba en compresión todavía existia ! rigidez a compresión en la dirreción radial factor=0.d0 IF (DSTRAN(1).lt.0.d0) THEN Knn=h*(sigmacb**2)/(2*GIcb) Ktt=E/1e18 !reduced stiffness terms! Kss=E/1e18 K33=E/1e18 DDS(1,1)=Knn DDS(2,2)=Kss DDS(3,3)=K33 DDS(4,4)=Ktt STRESS(1)=STRESS(1)+DDS(1,1)*DSTRAN(1) DO k=2,4 STRESS(k)=0.d0 ENDDO ELSE Knn=E/1e18 Ktt=E/1e18 Kss=E/1e18 K33=E/1e18 DDS(1,1)=Knn DDS(2,2)=Kss DDS(3,3)=K33 DDS(4,4)=Ktt DO k=1,4 STRESS(k)=0.d0 115 ENDDO ENDIF ENDIF ! Determinación de la matriz TANGENTE DO i=1,4 DO j=1,4 DDSDDE(i,j)=DDS(i,j) ENDDO ENDDO statev(1)=damage statev(2)=factor statev(3)=noel ! LTM (20/04/21) !cccccccccccccccccccccccccccccccccccccccccccccccc ! LEBIM interface option END !cccccccccccccccccccccccccccccccccccccccccccccccc ELSE IF (mtype.gt.1) THEN !cccccccccccccccccccccccccccccccccccccccccccccccc ! Discrete damage at continuum media starts !cccccccccccccccccccccccccccccccccccccccccccccccc ! No damage starts ! Isotropic Elasticity matrix update IF (damage.eq.0) THEN DO k1=1,3 DO k2=1,3 DDS(k2,k1)=ELAM ENDDO DDS(k1,k1)=EG2+ELAM ENDDO DDS(4,4)=EG2 ! Stress increment tensor update ! Trace of the strain tensor: TRVAL=DSTRAN(1)+DSTRAN(2)+DSTRAN(3) ! Lame equations for a linear elastic material: DO k=1,3 DSTRESS(k)=EG2*DSTRAN(k)+ELAM*TRVAL ENDDO ! Stress increment: DSTRESS(4)=EG2*DSTRAN(4) ! Stress tensor update DO k=1,4 STRESS(k)=STRESS(k)+DSTRESS(k) ENDDO !cc JMP (Sept-2022) ******************************************* if (mtype.eq.2) then ! Failure criteria: maximum principal stress, sigma1: sigmaprin1=((STRESS(1)+STRESS(2))/2)+ @(dsqrt(STRESS(1)**2 + STRESS(2)**2 + 4*STRESS(4)**2 - @(2*STRESS(1)*STRESS(2))))/2 factor=sigmaprin1/sigmayt IF (factor.ge.1.d0) THEN damage=1 else damage=0 ENDIF else if (mtype.eq.3) then 116 Apéndice A. Subrutina UMAT ! Failure criteria: Mohr-Coulomb 3D sigmaprin1=((STRESS(1)+STRESS(2))/2)+ @(dsqrt(STRESS(1)**2 + STRESS(2)**2 + 4*STRESS(4)**2 - @(2*STRESS(1)*STRESS(2))))/2 factor=(sigmaprin1/sigmayt)-(STRESS(3)/Suc) IF (factor.ge.1.d0) THEN damage=1 else damage=0 ENDIF else if (mtype.eq.4) then ! Failure criteria: von-Mises biaxial sigmaprin1=((STRESS(1)+STRESS(2))/2)+ @(dsqrt(STRESS(1)**2 + STRESS(2)**2 + 4*STRESS(4)**2 - @(2*STRESS(1)*STRESS(2))))/2 sigmaprin2=((STRESS(1)+STRESS(2))/2)- @(dsqrt(STRESS(1)**2 + STRESS(2)**2 + 4*STRESS(4)**2 - @(2*STRESS(1)*STRESS(2))))/2 VM=dsqrt((sigmaprin1)**2 + (sigmaprin2)**2 - @(sigmaprin1*sigmaprin2)) factor=(VM/sigmayt) IF (factor.ge.1.d0) THEN damage=1 else damage=0 ENDIF else if (mtype.eq.5) then ! Failure criteria: matrix volumetric strain energy density component Uv=((1-2*xnue)/(6*E))*(STRESS(1)+STRESS(2)+STRESS(3))**2 factor=dsqrt(Uv/Uc) IF (factor.ge.1.d0) THEN damage=1 else damage=0 ENDIF else if (mtype.eq.6) then ! Failure criteria: matrix total strain energy density sigmaprin1=((STRESS(1)+STRESS(2))/2)+ @(dsqrt(STRESS(1)**2 + STRESS(2)**2 + 4*STRESS(4)**2 - @(2*STRESS(1)*STRESS(2))))/2 sigmaprin2=((STRESS(1)+STRESS(2))/2)- @(dsqrt(STRESS(1)**2 + STRESS(2)**2 + 4*STRESS(4)**2 - @(2*STRESS(1)*STRESS(2))))/2 Uv=((1-2*xnue)/(6*E))*(STRESS(1)+STRESS(2)+STRESS(3))**2 Ud=(((1+xnue)/(3*E))*(sigmaprin1**2+sigmaprin2**2+STRESS(3)**2 - @((sigmaprin1*sigmaprin2)+(sigmaprin2*STRESS(3))+ @(sigmaprin1*STRESS(3))))) Ut=Uv+Ud factor=dsqrt(Ut/Uc) IF (factor.ge.1.d0) THEN damage=1 else damage=0 ENDIF endif ! ********************************************************************** ENDIF 117 ! No damage ends !cccccccccccccccccccccccccccccccccccccccccccccccc ! Continuum media (isotropic) damage starts IF (damage.eq.1) THEN DDS(1,1)=E/1.d18 !1.D18 is a double-precision exponent DDS(2,2)=E/1.d18 !thus, E/1.D18 is a reduced stiffness DDS(3,3)=E/1.d18 !term DDS(4,4)=E/1.d18 DO k=1,4 STRESS(k)=0.d0 ENDDO !the expression above can be enriched !by including stiffness terms under !compressive strains ENDIF ! LTM (20/04/21) ! Determinación de la matriz TANGENTE DO i=1,4 DO j=1,4 DDSDDE(i,j)=DDS(i,j) ENDDO ENDDO statev(1)=damage statev(2)=factor statev(3)=noel ENDIF !cccccccccccccccccccccccccccccccccccccccccccccccc ! End of discrete damage at continuum media !cccccccccccccccccccccccccccccccccccccccccccccccc RETURN END Apéndice B Fichero main_abaqus.py Código B.1 Fichero principal del SLA. #************************************************************* # SCRIPT: main_abaqus.py # Autor: Alfonso Gijon, a partir de una version anterior de # Luis Tavara, Enrique Paloma y Laura Moreno # Fecha: Febrero 2022 # ************************************************************ #Importacion modulos from os import chdir, path, remove, mkdir, system, listdir, getcwd from shutil import move, rmtree from glob import glob from operator import itemgetter from sys import argv import time import postprocess from abaqus import * from abaqusConstants import * from odbAccess import * from part import * from material import * from section import * from assembly import * from step import * from interaction import * from load import * from mesh import * from optimization import * from job import * from sketch import * from visualization import * from connectorBehavior import * # Initialization of the execution time and start_time = time.time() #--------------------------------------------------- 119 120 Apéndice B. Fichero main_abaqus.py # NOMBRES DE NUESTROS ARCHIVOS DE TRABAJO (LEER de inputSLA.txt) #--------------------------------------------------- fileSLA = open('inputSLA.txt','r') flag = fileSLA.readline() name_INP = fileSLA.readline().rstrip() flag = fileSLA.readline() name_UMAT = fileSLA.readline().rstrip() #------------------------------------------------------------- #------------------------------------------------------------- # Nuevas entradas por fichero #------------------------------------------------------------- # numero de iteraciones flag = fileSLA.readline() num_iteraciones = int(fileSLA.readline().rstrip()) #Frecuencia de guardado de ficheros ODB (si paso = 0 solo se guarda el ultimo) flag = fileSLA.readline() paso = int(fileSLA.readline().rstrip()) if (paso==0): respuesta_ODB = 'n' else: respuesta_ODB = 'y' #Modo: tolerancia o numero de ptos int a danar flag = fileSLA.readline() tol_o_ptosint = fileSLA.readline().rstrip() #Numero de cpus a usar flag = fileSLA.readline() cpu_num = fileSLA.readline().rstrip() fileSLA.close() #------------------------------------------------------------- # Orden para indicar cual es el directorio actual de trabajo (para poder usar en cluster) working_directory = getcwd() # Indicar el directorio de trabajo y crear directorio para el archivo ' datos_procesar.txt' chdir(working_directory) process_directory = working_directory+'/archivos_procesar' # Remove previous results if they exist if path.exists(process_directory): rmtree(process_directory) for file in listdir(working_directory): if file.endswith("_danar"): rmtree(file) continue if file.startswith(name_INP+"_") or file.endswith(".txt"): remove(file) mkdir(process_directory) #mkdir('C:\\archivos_procesar') #Dar nombre al archivo txt con el numero de iteraciones y la tolerancia if (float(tol_o_ptosint) < 1): name_txt_evolucion = 'evolucion_rotura'+'_'+ str(num_iteraciones) + 'steps'+'_'+ tol_o_ptosint + 'tolerancia' 121 else: name_txt_evolucion = 'evolucion_rotura'+'_'+str(num_iteraciones) + ' steps'+'_'+ tol_o_ptosint + 'ptos_danar' #--------------------------------------------------------------- # BUCLE PRINCIPAL #--------------------------------------------------------------- ofile = open("output.txt",'w') pp_time = 0.0 # Post-processing time flagsize=0 for step in range(num_iteraciones): step_time = time.time() # Dar nombre a los archivos de salida en la iteracion if (float(tol_o_ptosint) < 1): name_files = name_INP + '_'+str( num_iteraciones) + 'steps'+'_'+str(step + 1) + '_'+ tol_o_ptosint + ' tol' else: name_files = name_INP + '_'+str(num_iteraciones) + 'steps'+'_'+ str(step + 1) + '_'+ tol_o_ptosint + 'ptos_danar' #--------------------------------------------- # Ejecutar inp + UMAT #--------------------------------------------- # Read the model mdb.ModelFromInputFile(inputFileName=name_INP+'.inp', name= name_INP) # Create the job mdb.Job(atTime=None, contactPrint=OFF, description='', echoPrint=OFF, explicitPrecision=SINGLE, getMemoryFromAnalysis=True, historyPrint=OFF, memory=90, memoryUnits=PERCENTAGE, model=name_INP, modelPrint=OFF, multiprocessingMode=THREADS, name=name_files, nodalOutputPrecision=SINGLE, numCpus=int(cpu_num), numDomains=int(cpu_num), numGPUs=0, queue=None, resultsFormat=ODB, scratch='',type=ANALYSIS, userSubroutine='LEBIMkinkpstrain_uni.f', waitHours=0, waitMinutes=0) # Run the job mdb.jobs[name_files].submit(consistencyChecking=OFF) # Do not return control till job is finished running mdb.jobs[name_files].waitForCompletion() #--------------------------------------------- # POSTPROCESADO #--------------------------------------------- # Postprocesado ODB (modulo postprocess) time1 = time.time() postprocess.postprocess(name_files, tol_o_ptosint, name_txt_evolucion, str(step), working_directory) time2 = time.time() pp_time = pp_time + time2 - time1 # Obtener size de fichero datos_procesar if (step==0): flagsize=0 128 Índice de Figuras 4.1 (a) Problema de la inclusión circular (1) inmersa en una matriz (2); (b) grieta de interfase alrededor de la inclusión circular. Fuente: [33] 43 4.2 Geometría del modelo de la inclusión circular inmersa en una matriz finita. Fuente: [65] 44 4.3 Regiones de interfase para obtener el carácter asimétrico del despegue de grieta. Fuente: [65] 45 4.4 Mallas genéricas usadas para el estudio de convergencia de tensiones. (a) 1 elemento/grado de interfase; (b) 2 elementos/grado de interfase; (c) 3 elementos/grado de interfase 45 4.5 Curvas de equilibrio σ∞ x -CMOD de modelos unifibra con matriz finita. El control de carga se realiza en tensiones 46 4.6 Concepto de la apertura radial de grieta, CMOD 47 4.7 Distintos mallados de la matriz: (a) algoritmo de malla libre, (b) algoritmo tipo "mapped", (c) algoritmo libre cercano a la interfase 47 4.8 Estudios de objetividad de la malla de la matriz; usando el criterio de fallo en tensión principal máxima σ1, con St=80 MPa 49 4.9 Estudios de objetividad de la malla de la matriz; usando el criterio de la componente dilatacional de la energía de deformación, Uv, con Ucrit =1J/m350 4.10 Fallo exclusivo de la matriz con el criterio de la tensión principal máxima, Sut =80 MPa y control de carga en tensiones 51 4.11 Resultados de la evolución de la variable de estado SDV1, para distintos criterios de fallo de la matriz; el control de carga se realiza bajo tensión σ∞ x , con las propiedades mecánicas Sut =80 MPa, Suc=120 MPa, y Ucr=1J/m352 4.12 Variable de daño SDV1, del modelo unifibra para los distintos criterios de fallo evaluados; control de carga realizado en desplazamiento, u=0.01µm 53 4.13 Curvas de equilibrio para distintos criterios de fallo, en el modelo de matriz finita 54 4.14 Resultados de la evolución de la variable de estado SDV1, bajo el criterio de la componente dilatacional de energía de deformación, y para distintos valores de energía crítica de la matriz, (Ucr=1, 2 y 5J/m3); el control de carga se realiza bajo tensión σ∞ x55 4.15 Resultados de la evolución de la variable de daño SDV1, bajo el criterio de la energía total de deformación, y para distintos valores de energía crítica de la matriz, (Ucr=1, 2 y 5J/m3); el control de carga se realiza en desplazamientos 56 4.16 Evolución de la variable de estado SDV1, para distintos valores de propiedades de la interfase LEBIM; el control de carga se realiza bajo tensión σ∞ x57 4.17 Curvas de equilibrio σ∞ x -CMOD para el sistema fibra-matriz. Efecto de los valores críticos de la energía (a) de dilatación Uvol y (b) total Ut58 4.18 Efecto de variación en las propiedades de la interfase LEBIM, modelo de matriz finita 59 4.19 Configuración del problema de la inclusión (fibra) bajo tensión biaxial remota sin despegue de interfase (a) y con despegue de la interfase (b). Fuente: [35] 59 4.20 Evolución de la variable de estado SDV1, para distintos valores del parámetro de carga biaxial η ; el control de carga se realiza bajo tensión, con σ∞ x=1 MPa 61 4.21 Curvas (adimensionales) de equilibrio para distintos valores del parámetro de carga biaxial, η ; σ∞ x =1 MPa, δn= apertura radial de grieta, 2a = diámetro de la fibra 62 4.22 Esquema del modelo unifibra de matriz infinita 63 4.23 Discretización del modelo de elementos finitos unifibra de matriz infinita 63 4.24 Curvas de equilibrio σ∞- CMOD, para distinto número de elementos en la interfase LEBIM 64 4.25 Curvas σ∞- CMOD para distintos controles de carga en la matriz infinita 65 4.26 Variable de daño SDV1, del modelo unifibra para los distintos criterios de fallo evaluados; propiedades de la interfase LEBIM: ¯ GIc=2J/m2,¯ σc=90 MPa, λ=0.25, kt/kn=0.25, h=0.01µm 66 4.27 Curvas de equilibrio σ∞- CMOD, para distintos criterios de fallo de la matriz 67 4.28 Curvas de equilibrio σ∞ - CMOD, para distintos valores de las propiedades mecánicas de la interfase LEBIM 67 4.29 Variable de daño SDV1, del modelo unifibra para distintos valores de la energía crítica de la matriz; propiedades de la interfase LEBIM: ¯ GIc=2J/m2,¯ σc=90 MPa, λ=0.25, kt/kn=0.25, h=0.01µm 68 4.30 Curvas de equilibrio σ∞ x - CMOD, para distintos valores de la energía crítica de la matriz Ucr , usando (a) la componente de dilatación y (b) la energía total de deformación 69 5.1 Grieta transversal en un espécimen de fibra de carbono-epxy bajo carga transversal [37] 71 5.2 Modelos multifibra, lámina de espesor 40µm 72 5.3 Modelo multifibra, lámina de 90◦con 60µm de espesor 73 Índice de Figuras 129 5.4 Distribución de fibras para los modelos 40I (a) y 40D (b) 73 5.5 Detalle del mallado para el modelo 40D 74 5.6 Mallado de la matriz 74 5.7 Morfología del inicio del daño en el modelo 40D, para el paso 100 del SLA, criterio de la componente volumétrica de la energía de deformación para el fallo de la matriz 76 5.8 Evolución del daño en el modelo 40D; criterio de la componente dilatacional de la energía de deformación con Ucr=1J/m377 5.9 Evolución del daño en el modelo 40I; criterio de la componente dilatacional de la energía de deformación con Ucr=1J/m378 5.10 Deformaciones y tensiones críticas de los modelos de 40µm 79 5.11 Desarrollo de la variable de daño, modelo 40D. Influencia del valor crítico de la componente dilatacional de la energía de deformación 80 5.12 Curvas tensión-deformación en la lámina de 90 grados, modelo 40D, y para distintos valores de la energía critica de la matriz 81 5.13 Evolución del daño, modelo multifibras 40D; efectos del número de puntos de integración a dañar por paso SLA 83 5.14 Efectos del número de puntos de integración a dañar por paso del SLA; modelo 40D 84 5.15 Progresión del daño, modelo 60D; criterio de la componente dilatacional de la energía de deformación con Ucr=1- y 5-J/m385 5.16 Curvas tensión-deformacion en la lámina de 90 grados, modelo 60D, y para distintos valores de la energía crítica de la matriz 86 5.17 Evolución del daño, multifibras modelo 60D, para el criterio de la componente total de la energía de deformación, Ut88 5.18 Curvas εx-Paso SLA, modelo 60D, usando el criterio de la energía total de deformación para la matriz 89 6.1 Modelo integrado multiescala de un laminado [0,90n]s. Fuente: [43] 91 6.2 Escenarios de grietas en una lámina de [ 0,90n ] s : (a) acercándose a la interfase, (b) terminando en la interfase, (c) extendiéndose a lo largo de la interfase, (d) acercándose a una interfase fisurada. Fuente: [9] 92 6.3 Modelos multiescala usados para el análisis; la unidad de las dimensiones del dibujo es el µm92 6.4 Condiciones de contorno y ubicación de la grieta transversal en el modelo multiescala 93 6.5 Distribución de fibras en la celda de los modelos multiescala 93 6.6 Detalles del mallado de las celdas en los modelos multiescala 94 6.7 Progresión del daño, modelo multiescala con 100 fibras. Para facilitar su visualización, los resultados de la variable SDV1 en cada fila se reportan para distintos pasos del SLA 96 6.8 Deformación aplicada por paso del SLA, modelo de 100 fibras para distintos valores de la densidad de energía crítica 97 6.9 Curvas σ∞ ydeformación por daño, en el modelo de 100 fibras 98 6.10 Progresión del daño, modelo multiescala bajo carga en el laminado de 0◦99 6.11 Deformación aplicada por paso del SLA, modelo de 100 fibras bajo carga en laminado de 0 ◦ y ante distintos valores de la densidad de energía de dilatación. crítica 100 6.12 Curvas σ∞ y−ε(daño), modelo de 100 fibras bajo carga en laminado de 0◦101 6.13 Progresión de la variable de daño en la celda de 160 fibras, y para dos distintos valores de la energía crítica de la matriz 102 6.14 Curvas tensión-deformaciones en el modelo multiescala de 160 fibras 103 6.15 Progresión del daño en el modelo multiescala con 160 fibras, bajo carga en el laminado de 0 ◦ , y para dos valores distintos de la densidad de energía crítica Ucr 104 6.16 Curvas tensión-deformaciones en el modelo multiescala de 160 fibras, bajo carga en el laminado de 0 grados 105 7.1 Laminado [0/90]scon daño producido por una carga transversal. Fuente: [9] 108 7.2 Modos de fractura. Fuente: [72] 108 7.3 Sección representativa de una lámina bajo las tensiones σ1,σ2,τ12 109 Índice de Tablas 3.1 Algunos métodos de indagación de las entidades en los objetos Odb 38 4.1 Propiedades mecánicas del conjunto fibra-matriz 44 4.2 Propiedades mecánicas de la interfase LEBIM 44 4.3 Convergencia de las tensiones remotas críticas para el despegue de la interfase 46 4.4 Valores pico de σ∞ x/¯ σcen función de η, con σ∞ x=1 MPa 62 4.5 Resultados de las tensiones remotas críticas, σ∞ c, para el despegue de la interfase 64 5.1 Dimensiones de los modelos multifibra 73 5.2 Número de elementos en la matriz, fracción de volumen (Vf) y número de fibras por modelo 75 5.3 Propiedades mecánicas de la lámina ortótropa 75 6.1 Dimensiones generales, de la fibra y la interfase, en la celda de los modelos multiescala 93 6.2 Propiedades mecánicas de la lámina ortótropa de 0◦94 131 Bibliografía [1] Abaqus, Abaqus 6.11, Dassault Systemes Simulia Corporation, Providence, RI, USA (2011). [2] GaMAC Alfano and MA1011 Crisfield, Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues, International journal for numerical methods in engineering 50 (2001), no. 7, 1701–1736. [3] A Arteiro, G Catalanotti, AR Melro, P Linde, and Pedro P Camanho, Micro-mechanical analysis of the in situ effect in polymer composite laminates, Composite Structures 116 (2014), 827–840. [4] L Asp, Lars A Berglund, and R Talreja, Prediction of matrix-initiated transverse failure in polymer composites, Composites Science and Technology 56 (1996), no. 9, 1089–1097. [5] LE Asp, Lars A Berglund, and Ramesh Talreja, Effects of fiber and interphase on matrix-initiated transverse failure in polymer composites, Composites science and technology 56 (1996), no. 6, 657–665. [6] Grigory Isaakovich Barenblatt, The mathematical theory of equilibrium cracks in brittle fracture, Advances in applied mechanics 7(1962), 55–129. [7] Zdeněk P Bažant and Byung H Oh, Crack band theory for fracture of concrete, Matériaux et construction 16 (1983), no. 3, 155–177. [8] Zhuming Bi, Chapter 8 - applications—solid mechanics problems, Finite Element Analysis Applications (Zhuming Bi, ed.), Academic Press, 2018, pp. 281–339. [9] A. Blázquez, V. Mantič, F. París, and N.L. McCartney, Bem analysis of damage progress in 0/90 laminates, Engineering Analysis with Boundary Elements 33 (2009), no. 6, 762–769. [10] A Carpinteri, P Cornetti, and N Pugno, Edge debonding in frp strengthened beams: Stress versus energy failure criteria, Engineering Structures 31 (2009), no. 10, 2436–2447. [11] D Castillo, Implementación de un modelo de interfase en el programa de elementos finitos abaqus, Aplicación a materiales compuestos, Proyecto fin de carrera, Universidad de Sevilla (2014). [12] Wai-Fah Chen and Da-Jian Han, Plasticity for structural engineers, J. Ross Publishing, 2007. [13] E Correa, E Kristofer Gamstedt, F Paris, and V Mantič, Effects of the presence of compression in transverse cyclic loading on fibre–matrix debonding in unidirectional composite plies, Composites Part A: Applied Science and Manufacturing 38 (2007), no. 11, 2260–2269. [14] E Correa, F París, and V Mantič, Fiber–matrix debonding in composite materials: transverse loading, Modeling Damage, Fatigue and Failure of Composite Materials, Elsevier, 2016, pp. 97–116. [15] R. de Borst, Computation of post-bifurcation and post-failure behavior of strain-softening solids, Computers and Structures 25 (1987), no. 2, 211–224. [16] Matthew J DeJong, Max AN Hendriks, and Jan G Rots, Sequentially linear analysis of fracture under non-proportional loading, Engineering Fracture Mechanics 75 (2008), no. 18, 5042–5056. 133 134 Bibliografía [17] Donald S Dugdale, Yielding of steel sheets containing slits, Journal of the Mechanics and Physics of Solids 8(1960), no. 2, 100–104. [18] EK Gamstedt and BA Sjögren, Micromechanisms in tension-compression fatigue of composite laminates containing transverse plies, Composites Science and Technology 59 (1999), no. 2, 167–178. [19] IG García, V Mantič, and E Graciani, Debonding at the fibre–matrix interface under remote transverse tension. one debond or two symmetric debonds?, European Journal of Mechanics-A/Solids 53 (2015), 75–88. [20] Emmanuel E Gdoutos, Fracture mechanics criteria and applications, vol. 10, Kluwer Academic Publishers, 1990. [21] T Guillén-Hernández, A Quintana-Corominas, IG García, J Reinoso, M Paggi, and A Turon, In-situ strength effects in long fibre reinforced composites: A micro-mechanical analysis using the phase field approach of fracture, Theoretical and Applied Fracture Mechanics 108 (2020), 102621. [22] Zvi Hashin, Failure criteria for unidirectional fiber composites, (1980). [23] Max AN Hendriks and Jan G Rots, Sequentially linear versus nonlinear analysis of rc structures, Engineering Computations (2013). [24] Arne Hillerborg, Fracture mechanics concepts applied to moment capacity and rotational capacity of reinforced concrete beams, Engineering Fracture Mechanics 35 (1990), no. 1, 233–240, Special Issue Fracture and Damage of Concrete and Rock. [25] John W Hutchinson and Zhigang Suo, Mixed mode cracking in layered materials, Advances in applied mechanics 29 (1991), 63–191. [26] S Invernizzi, D Trovato, MAN Hendriks, and AV Van de Graaf, Sequentially linear modelling of local snap-back in extremely brittle structures, Engineering Structures 33 (2011), no. 5, 1617–1625. [27] Christos Kassapoglou, Design and analysis of composite structures, 2010. [28] Jochen Kopp and Walter Michaeli, Dimensioning of thick laminates using new iff strength criteria and some experiments for their verification, Spacecraft Structures, Materials and Mechanical Engineering, vol. 386, 1996, p. 305. [29] Steen Krenk, Energy release rate of symmetric adhesive joints, Engineering Fracture Mechanics 43 (1992), no. 4, 549–559. [30] Lothar Kroll and Werner Hufenbach, Physically based failure criterion for dimensioning of thick-walled laminates, Applied composite materials 4(1997), no. 5, 321–332. [31] Dominique Leguillon, Strength or toughness? a criterion for crack onset at a notch, European Journal of Mechanics-A/Solids 21 (2002), no. 1, 61–72. [32] Stefano Lenci, Analysis of a crack at a weak interface, International Journal of Fracture 108 (2001), no. 3, 275–290. [33] V Mantič, Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. application of a coupled stress and energy criterion, International journal of Solids and Structures 46 (2009), no. 6, 1287–1304. [34] V Mantič and IG García, Crack onset and growth at the fibre–matrix interface under a remote biaxial transverse load. application of a coupled stress and energy criterion, International Journal of Solids and Structures 49 (2012), no. 17, 2273–2290. [35] Vladislav Mantič, Luis Távara, Antonio Blázquez, Enrique Graciani, and Federico París, A linear elastic-brittle interface model: application for the onset and propagation of a fibre-matrix interface crack under biaxial transverse loads, International Journal of Fracture 195 (2015), no. 1, 15–38. [36] Shaker A Meguid, Engineering fracture mechanics, Springer, 1989. Bibliografía 135 [37] Laura Moreno, Luis Távara, and Elena Correa, Numerical analysis of the crack paths produced by fibre–matrix interface failure in cross-ply lfrp laminates, Composite Structures 284 (2022), 115222. [38] Laura Moreno, Luis Távara, Elena Correa, and Federico París, Sequential linear analysis for the prediction of the symmetrical or non-symmetrical character of the debond onset and propagation along a fiber-matrix interface, Journal of Multiscale Modelling 10 (2019), no. 03, 1842004. [39] Enrique Paloma Castro, Implementación de un método de solución lineal secuencial en el programa de elementos finitos abaqus para el modelo de interfase elástica lineal frágil, (2016). [40] Anna Pandolfi and Michael Ortiz, An eigenerosion approach to brittle fracture, International Journal for Numerical Methods in Engineering 92 (2012), no. 8, 694–714. [41] F París, Teoría de la elasticidad, Grupo de Elasticidad y Resistencia de Materiales, Sevilla, España, 2000. [42] F Parıs, E Correa, and J Cañas, Micromechanical view of failure of the matrix in fibrous composite materials, Composites Science and Technology 63 (2003), no. 7, 1041–1052. [43] F París, ML Velasco, and E Correa, The scale effect in composites: An explanation physically based on the different mechanisms of damage involved in failure, Composite Structures 257 (2021), 113089. [44] Federico París, A study of failure criteria of fibrous composite materials, (2001). [45] Federico París, Elena Correa, and Vladislav Mantič, Kinking of transversal interface cracks between fiber and matrix, (2007). [46] Kyoungsoo Park, Glaucio H Paulino, and Jeffery R Roesler, A unified potential-based cohesive model of mixed-mode fracture, Journal of the Mechanics and Physics of Solids 57 (2009), no. 6, 891–908. [47] A. Puck, J. Kopp, and M. Knops, Guidelines for the determination of the parameters in puck’s action plane strength criterion, Composites Science and Technology 62 (2002), no. 3, 371–378. [48] A Puck and H Schürmann, Failure analysis of frp laminates by means of physically based phenomenological models, Failure criteria in fibre-reinforced-polymer composites, Elsevier, 2004, pp. 832–876. [49] Martin Rädel, Christian Willberg, and Daniel Krause, Peridynamic analysis of fibre-matrix debond and matrix failure mechanisms in composites under transverse tensile load by an energy-based damage criterion, Composites Part B: Engineering 158 (2019), 18–27. [50] J Reinoso, A Blázquez, A Estefani, F París, J Cañas, Edward Arévalo, and F Cruz, Experimental and three-dimensional global-local finite element analysis of a composite component including degradation process at the interfaces, Composites Part B: Engineering 43 (2012), no. 4, 1929–1942. [51] Jan G Rots and Stefano Invernizzi, Regularized sequentially linear saw-tooth softening model, International journal for numerical and analytical methods in geomechanics 28 (2004), no. 7-8, 821–856. [52] Jan Gerrit Rots, Computational modeling of concrete fracture, (1988). [53] JG Rots, Sequentially linear continuum model for concrete fracture, Fracture mechanics of concrete structures 13 (2001), 831–839. [54] Hiroshi Saito, Hiroki Takeuchi, and Isao Kimpara, Experimental evaluation of the damage growth restraining in 90 layer of thin-ply cfrp cross-ply laminates, Advanced Composite Materials 21 (2012), no. 1, 57–66. [55] E Schlangen and JGMi Van Mier, Experimental and numerical analysis of micromechanisms of fracture of cement-based composites, Cement and concrete composites 14 (1992), no. 2, 105–118. [56] G C Sih, A three-dimensional strain energy density factor theory of crack propagation, Mechanics of fracture initiation and propagation, Springer, 1991, pp. 23–56. [57] P D Soden, M J Hinton, and A S Kaddour, Lamina properties, lay-up configurations and loading conditions for a range of fibre reinforced composite laminates, Failure criteria in fibre-reinforced- polymer composites, Elsevier, 2004, pp. 30–51. 136 Bibliografía [58] C.T. Sun and Z.-H. Jin, Chapter 6 - crack tip plasticity, Fracture Mechanics (C.T. Sun and Z.-H. Jin, eds.), Academic Press, Boston, 2012, pp. 123–169. [59] Ramesh Talreja, Matrix and fiber–matrix interface cracking in composite materials, Modeling Damage, Fatigue and Failure of Composite Materials, Elsevier, 2016, pp. 87–96. [60] L Távara and V Mantič, Crack paths formed by multiple debonds in lfrp composites, Mechanics Research Communications 84 (2017), 148–154. [61] L Távara, V Mantic, E Graciani, J Cañas, and F París, Analysis of a crack in a thin adhesive layer between orthotropic materials: an application to composite interlaminar fracture toughness test, Computer Modeling in Engineering & Sciences(CMES) 58 (2010), no. 3, 247–270. [62] L Távara, V Mantič, E Graciani, and F París, Modelling interfacial debonds in unidirectional fibrereinforced composites under biaxial transverse loads, Composite Structures 136 (2016), 305–312. [63] L Távara, Vladislav Mantič, E Graciani, and F París, Bem analysis of crack onset and propagation along fiber–matrix interface under transverse tension using a linear elastic–brittle interface model, Engineering Analysis with Boundary Elements 35 (2011), no. 2, 207–222. [64] Luis Távara, Damage initiation and propagation in composite materials. boundary element analysis using weak interface and cohesive zone models, Ph.D. thesis, PhD Thesis. Universidad de Sevilla: Sevilla, 2010. [65] Luis Távara, José Reinoso, David Castillo, and Vladislav Mantič, Mixed-mode failure of interfaces studied by the 2d linear elastic–brittle interface model: Macro-and micro-mechanical finite-element applications in composites, The Journal of Adhesion 94 (2018), no. 8, 627–656. [66] M Toya, A crack along the interface of a circular inclusion embedded in an infinite solid, Journal of the Mechanics and Physics of Solids 22 (1974), no. 5, 325–348. [67] Stephen W Tsai and Edward M Wu, A general theory of strength for anisotropic materials, Journal of composite materials 5(1971), no. 1, 58–80. [68] A Turon, Pedro Ponces Camanho, J Costa, and CG Dávila, A damage model for the simulation of delamination in advanced composites under variable-mode loading, Mechanics of materials 38 (2006), no. 11, 1072–1089. [69] TJ Vaughan and CT McCarthy, Micromechanical modelling of the transverse damage behaviour in fibre reinforced composites, Composites Science and Technology 71 (2011), no. 3, 388–396. [70] John T Wang, Investigating some technical issues on cohesive zone modeling of fracture, Journal of engineering materials and technology 135 (2013), no. 1, 011003. [71] Dennis Wingender and Daniel Balzani, Simulation of crack propagation based on eigenerosion in brittle and ductile materials subject to finite strains, Archive of Applied Mechanics 92 (2022), no. 4, 1199–1221. [72] Alan T. Zehnder, Modes of fracture. book: Encyclopedia of tribology, pp. 2292–2295, Springer US, Boston, MA, 2013.