Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones

Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería de Tecnologías Industriales Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones Dpto. Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Autor: Angel Hidalgo Armayones Tutor: Alfredo Navarro Robles Sevilla, 2020 iii Trabajo Fin de Grado Grado en Ingeniería de Tecnologías Industriales Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones Autor: Angel Hidalgo Armayones Tutor: Alfredo de Jesús Navarro Robles Catedrático de Universidad Dpto. de Ingeniería Mecánica y Fabricación Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2020 v Trabajo Fin de Grado: Simulación numérica de endurecimiento en ensayos biaxiales cíclicos noproporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones Autor: Ángel Hidalgo Armayones Tutor: Alfredo de Jesús Navarro Robles El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2020 El Secretario del Tribunal vii A mi familia A mis maestros ix Agradecimientos Quiero agradecer a todas las personas cercanas a mí lo mucho que me han inspirado para realizar este trabajo, y el apoyo recibido por la parte docente. Angel Hidalgo Armayones Sevilla, 2020 ÍNDICE DE FIGURAS Figura 1.1 4 Figura 1.2 5 Figura 1.3 6 Figura 1.4 7 Figura 3.1.1 15 Figura 3.1.2 16 Figura 4.1.2.1.1 21 Figura 4.1.2.1.2 21 Figura 4.1.2.1.3 22 Figura 4.1.2.1.4 22 Figura 4.1.2.1.5 23 Figura 4.1.2.1.7 23 Figura 4.1.2.1.8 24 Figura 4.1.2.2.1 24 Figura 4.1.2.2.2 25 Figura 4.1.2.2.3 25 Figura 4.1.2.2.4 26 Figura 4.1.2.2.5 26 Figura 4.1.2.2.6 27 Figura 4.1.2.2.7 27 Figura 4.1.2.3.1 28 Figura 4.1.2.3.2 28 Figura 4.1.2.3.3 29 Figura 4.1.2.3.4 29 Figura 4.1.2.3.5 30 Figura 4.1.2.3.6 30 Figura 4.1.2.3.7 31 Figura 4.1.3.1 32 Figura 4.1.3.2 32 Figura 4.2.2.1.1 36 Figura 4.2.2.1.2 36 Figura 4.2.2.1.3 37 Figura 4.2.2.1.4 37 Figura 4.2.2.1.5 38 Figura 4.2.2.1.6 38 xvii Figura 4.2.2.1.7 39 Figura 4.2.2.2.1 39 Figura 4.2.2.2.2 40 Figura 4.2.2.2.3 40 Figura 4.2.2.2.4 41 Figura 4.2.2.2.5 41 Figura 4.2.2.2.6 42 Figura 4.2.2.2.7 42 Figura 4.2.2.3.1 43 Figura 4.2.2.3.2 43 Figura 4.2.2.3.3 44 Figura 4.2.2.3.4 44 Figura 4.2.2.3.5 45 Figura 4.2.2.3.6 45 Figura 4.2.2.3.7 46 Figura 4.3.2.1.1 48 Figura 4.3.2.1.2 49 Figura 4.3.2.1.3 49 Figura 4.3.2.1.4 50 Figura 4.3.2.1.5 50 Figura 4.3.2.1.6 51 Figura 4.3.2.1.7 51 Figura 4.3.2.2.1 52 Figura 4.3.2.2.2 52 Figura 4.3.2.2.3 53 Figura 4.3.2.2.4 53 Figura 4.3.2.2.5 54 Figura 4.3.2.2.6 54 Figura 4.3.2.2.7 55 Figura 4.3.2.3.1 55 Figura 4.3.2.3.2 56 Figura 4.3.2.3.3 56 Figura 4.3.2.3.5 57 Figura 4.3.2.3.6 58 Figura 4.3.2.3.7 58 Figura 5.1.1 62 Figura 5.1.2 62 Figura 5.2.1 63 Figura 5.2.2 64 Figura 5.3.2.1 65 Figura 5.3.2.2 66 Figura 5.3.2.3 66 Figura 5.3.2.4 67 Figura 5.3.2.5 67 Figura 5.3.2.6 68 Figura 5.3.2.7 68 Figura 5.3.4 69 xix Notación 𝝈 Vector tensión según métrica establecida 𝜺 Vector deformación según métrica establecida 𝜎 Componente normal del vector tensión 𝜏 Componente tangencial del vector tensión 𝜀 Componente longitudinal del vector deformación 𝛾 Componente transversal del vector deformación |𝝈| Módulo efectivo de tensiones |𝜺| Módulo efectivo de deformaciones 𝑑( ) Incremento diferencial 𝜃 Ángulo entre dos vectores 𝑊𝑝 Trabajo plástico 𝜺𝒑 Vector deformación plástica 𝜺𝒆 Vector deformación elástica 𝐻( ) Función de endurecimiento 𝝓(|𝜎|) Módulo de endurecimiento 𝒏 Normal a la superficie de plastificación 𝑛1 Componente x de 𝒏 𝑛2 Componente y de 𝒏 𝐸 Módulo de Young del material 𝐺 Módulo de Elasticidad transversal del material 𝐾 Parámetro Ramberg-Osgood del material 𝑛 Parámetro Ramberg-Osgood del material 𝝈𝑘 Punto de inversión de carga activo 𝐶𝑘 Circunferencia de plastificación activa asociada a 𝝈𝑘 𝝈𝑐,𝑘 Centro de la circunferencia de plastificación 𝐶𝑘 𝑞 Distancia entre punto actual y punto de descarga 𝑡 Parámetro de muestreo para modelización en MATLAB 𝑡𝑐𝑎𝑟𝑔𝑎 Parámetro 𝑡 para tramo de aumento de deformaciones impuestas |𝜀𝑚𝑎𝑥| Parámetro de deformaciones efectivas impuestas en la simulación 𝑡𝑓𝑖𝑛 Valor final del parámetro 𝑡 𝑁𝐶 Número de ciclos simulados 𝑆𝑓 Tensión de fallo del material 𝜀𝑓 Deformación de fallo del material 𝑏 Parámetro curva 𝜀−𝑁 𝑎 Parámetro curva 𝜀−𝑁 1 1 TEORÍA DE PLASTICIDAD BASADA EN LA DISTANCIA as simulaciones que se realizarán en este proyecto utilizan una teoría de plasticidad basada en la distancia entre puntos representados en el espacio tensional. Esta teoría está basada en los métodos de cálculos de tensiones y deformaciones utilizados mediante la ecuación de Ramberg-Osgood para deformaciones a bajo número de ciclos en ensayos monoaxiales, permite el cálculo de estos valores para el caso multiaxial. Para ello utilizaremos el criterio de plastificación de von Mises, definiendo así la métrica del espacio tensional que se describirá a continuación. El método descrito a continuación consiste en una teoría de plasticidad cíclica, de forma que se puedan aplicar al modelo multiaxial procedimientos similares a los que se aplican en los problemas monoaxiales mediante el método de las Deformaciones locales. Para ello, se utilizará una representación espacial de tensiones y deformaciones, cuyos ejes se determinarán según el criterio de plastificación de von Mises. Esto permite utilizar el concepto de distancia entre tensiones para poder visualizar de una forma gráfica y sencilla todo el procedimiento. Esta teoría se desarrolla completamente en las referencias [2], [3] y [4], pero se ha considerado útil expresar sus ideas claves para una mayor comprensión de este proyecto por parte del lector. Para entender mejor esta teoría, empezaremos desarrollándola conforme al modelo biaxial, para más adelante aplicarla a casos multiaxiales. Para comenzar debemos definir la métrica en la cual representaremos las tensiones y deformaciones del modelo. Para ello trataremos los vectores tensión y deformación como vectores, los cuales tendrán componentes normales y tangenciales. En el caso de las tensiones, las definiremos mediante el vector 𝝈, que a su vez se compondrá de las componentes 𝜎 (tensión normal) y 𝜏 (tensión tangencial). El vector de deformaciones a su vez lo definiremos como el vector 𝜺, compuesto por las componentes 𝜀 (deformación longitudinal) y 𝛾 (deformación transversal): 𝝈=(𝜎,𝜏) (1.1) 𝜺=(𝜀,𝛾) (1.2) El criterio de plastificación de von Mises nos permite definir superficies de plastificación, lo cual significa que la plastificación comienza cuando se supera el valor crítico: |𝜎|2=23𝜎2+2𝜏2 (1.3) Como esta expresión no es la norma euclidiana, podemos decir que la base vectorial no estará formada por vectores ortonormales (perpendiculares y de módulo unidad), si no por vectores que su producto escalar sea proporcional a la expresión anterior. En el espacio vectorial de tensiones, el ángulo entre dos vectores tensión 𝝈𝟏=(𝜎1,𝜏1) y 𝝈𝟐=(𝜎2,𝜏2) se calcula mediante: cos𝜃= 𝝈𝟏 .𝝈𝟐 |𝝈𝟏 | .|𝝈𝟐 | (1.4) L Teoría de plasticidad basada en la distancia 2 Siendo el producto escalar 𝝈𝟏 .𝝈𝟐 : 𝝈𝟏 .𝝈𝟐=23𝜎1 𝜎2+2𝜏1 𝜏2 (1.5) Sin embargo, para las componentes del vector deformación es diferente, ya que debido a que la expresión para calcular el trabajo plástico es: 𝑑𝑊𝑝=𝜎𝑑𝜀𝑝+𝜏𝑑𝛾𝑝 (1.6) Para que esta expresión sea así debemos definir para el espacio vectorial de deformaciones la base dual tomando como referencia el espacio tensional. Como los vectores de las bases duales son ortonormales, podemos calcular el módulo efectivo de deformaciones con la siguiente expresión: |𝜀𝑝|=√(𝜀𝑝)2+ (𝜀𝑝 2)2+ (𝜀𝑝 2)2+ (𝛾𝑝 2)2+ (𝛾𝑝 2)2 = √32(𝜀𝑝)2+12(𝛾𝑝)2 (1.7) Avanzando en la materia, tendremos en cuenta para este método una hipótesis de endurecimiento isótropo, asumiendo que el radio de la circunferencia de plastificación según von Mises sólo es función del incremento de deformación efectiva: |𝜎|=𝐻(∫|𝑑𝜀𝑝|) (1.8) Para la mayoría de los metales que se endurecen plásticamente es una función creciente. La expresión inversa a la ecuación anterior viene dada por: ∫|𝑑𝜀𝑝|=𝐻−1(|𝜎|) (1.9) Siendo la función 𝐻−1 la función inversa a 𝐻. Derivando obtenemos: |𝑑𝜀𝑝|=𝑑𝐻−1(|𝜎|) 𝑑|𝜎|=𝝓(|𝜎|)𝑑|𝜎| (1.10) Siendo 𝝓(|𝜎|) la función que relaciona el incremento de la deformación plástica efectiva con el incremento de la tensión efectiva. La denominaremos módulo de endurecimiento. Esta ecuación implica que la deformación plástica sólo ocurre cuando aumenta el módulo efectivo de tensiones, por lo que el endurecimiento sólo depende de la distancia del punto de tensión al origen de coordenadas. El incremento de deformación plástica tendrá una dirección normal a la circunferencia |𝝈|=𝑐𝑡𝑒, según la expresión: 𝑑𝜀𝑝= 𝜙(|𝜎|)𝑑|𝜎|𝒏 (1.11) Siendo 𝒏 el vector normal a la circunferencia |𝝈|=𝑐𝑡𝑒. Esta normal se define derivando según la métrica descrita en las ecuaciones (1.3) y (1.7): 3 3 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones 𝑛1=𝜕|𝜎| 𝜕𝜎=23𝜎 |𝜎| (1.12) 𝑛2=𝜕|𝜎| 𝜕𝜏 =2𝜏 |𝜎| (1.13) Sustituyendo estas expresiones en la ecuación (1.11) y añadiendo las componentes elásticas según la Ley de Hooke, tendremos a nuestra disposición un sistema de ecuaciones diferenciales el cual describe la evaluación de las deformaciones plásticas y elásticas del material introduciendo como parámetro una evolución en las tensiones. Estas ecuaciones sólo se aplican para el caso 𝑑|𝜎|>0, es decir, cuando el módulo efectivo de tensiones aumenta. Si 𝑑|𝜎|=0 , el punto de tensión representado en el diagrama de tensión normal 𝜎 y tensión tangencial 𝜏 describirá una trayectoria circular, manteniendo la misma distancia al origen de coordenadas. Cuando esto ocurre, podemos afirmar que solo ocurren cambios en las tensiones y en las deformaciones de forma elástica. Para obtener la función de endurecimiento 𝜙(|𝜎|) partiremos de la ecuación de Ramberg-Osgood, según la referencia [1]: 𝜀=𝜎𝐸+(𝜎𝐾)1/𝑛 (1.14) En el cual la deformación plástica se define: 𝜀𝑃=(𝜎𝐾)1/𝑛 (1.15) Teniendo en cuenta que para el caso monoaxial: |𝝈|=√23𝜎 (1.16) |𝜺𝑃|=√32𝜀𝑃 (1.17) Y sustituyendo en la ecuación (1.14), finalmente tenemos: |𝜺𝑃|=(√32)(1𝑛+1)(|𝝈| 𝐾)1𝑛 (1.18) Teoría de plasticidad basada en la distancia 4 La cual podemos derivar para obtener la expresión de 𝜙(|𝝈|) según la expresión (1.10) 𝜙(|𝝈|)=(√32)(1𝑛+1)1𝑛(|𝝈| 𝐾)1𝑛1 |𝝈| (1.19) Expresión que sustituida en (1.11) se puede aplicar para el caso multiaxial, siempre que calculemos los módulos |𝝈| y |𝜺𝑃| según las métricas adoptadas previamente. Recordamos que esto es válido siempre y cuando el módulo efectivo de tensiones aumente, pero ¿qué ocurre si disminuyen? Tenemos que empezar a hablar de puntos de descarga. Los puntos de descarga se pueden entender de forma análoga a los picos de tensiones que se determinan para aplicar el método de conteo de ciclos Rainflow. Primero representemos un estado tensional cualquiera: Figura 1.1 Es fácil de entender que si el siguiente paso de 𝝈 tiene una distancia al origen mayor al radio de la circunferencia representada en la figura 1.1, decimos que en nuestro ensayo el material está en carga, por lo que pueden aplicarse las ecuaciones descritas en el apartado anterior. Pero si ocurre lo contrario, es decir, la distancia es menor, se crea un punto de descarga. El primer punto de descarga lo definiremos como 𝝈0. Este punto se define como el valor cuya distancia al origen es máxima, antes de que se produzca el primer descenso de |𝝈|. Se creará entonces una nueva circunferencia de plastificación según la siguiente imagen: 5 5 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones Figura 1.2 Entonces, debemos calcular las distancias respecto a esta nueva circunferencia 𝐶 de plastificación, cuyo origen, difiere del origen de coordenadas. Para calcular las distancias y poderlas utilizar en el modelo de ecuaciones diferenciales indicado anteriormente, debemos utilizar la siguiente expresión: 𝑞=|𝝈−𝝈0| cos𝜃 (1.20) Donde 𝜃 es el ángulo entre los segmentos que conectan 𝝈0 con el origen de coordenadas y con el punto actual de tensión 𝝈. Si el módulo efectivo de tensiones vuelve a aumentar, 𝑑|𝝈|>0, el radio de 𝐶 crecerá, siendo siempre tangente a 𝐶0 en el punto de descarga 𝝈0. Mediante el producto escalar podemos obtener el ángulo 𝜃: −𝝈𝟎∗(𝝈−𝝈𝟎)=|𝝈0| |𝝈−𝝈0|cos𝜃 (1.21) Con lo cual la ecuación (1.20) nos queda: 𝑞=|𝝈−𝝈0| cos𝜃 =|𝝈−𝝈0|2|𝝈0| −𝝈𝟎∗(𝝈−𝝈𝟎) (1.22) Simplificando: 𝑞= 23(𝜎−𝜎0)2+2(𝜏−𝜏0)2 23(−𝜎0)(𝜎−𝜎0)+2(−𝜏0)(𝜏−𝜏0)𝑞0 2 (1.23) Donde: 𝑞0 2=√23𝜎02+2𝜏02 (1.24) Representando con 𝑞0 el módulo efectivo de tensiones del punto de descarga, es decir, el tamaño de la circunferencia 𝐶0. Consideraciones teóricas del modelo de cálculo 12 Este parámetro se definirá con mayor detenimiento en el posterior apartado, donde hablaremos del modelo numérico que se creará para realizar las simulaciones. La expresión de 𝝓(|𝜎|) viene dada por la ecuación (1.19) la cual expresa el flujo plástico dado por los coeficientes de Ramberg-Osgood. 2.2 Análisis de estado límite e hipótesis principales Una vez tenemos definido claramente el sistema de ecuaciones que debemos integrar para realizar nuestras simulaciones, es interesante realizar un pequeño análisis de dichas ecuaciones. Si analizamos las expresiones (1.19) y (2.1.14), podemos tener en cuenta una serie de hipótesis que marcarán la naturaleza de las simulaciones propuestas y el transcurso del presente proyecto. Según la teoría propuesta, al mantener las deformaciones en módulo constante, el material se endurece y plastifica cada vez más, por lo que |𝜎| es cada vez mayor, y con él 𝝓(|𝜎|) se va incrementando. Al ser 𝑛 del orden de 0,2 o menor, podemos afirmar que 𝝓(|𝜎|) crecerá a un ritmo mayor que |𝜎|2, siendo el cociente 𝝓(|𝜎|) |𝜎|2 cada vez mayor a medida que aumentan los ciclos a los que está sometido el material. Es entonces cuando podemos despreciar el término elástico frente al plástico, pues: 1𝐸1𝐺≪𝝓(|𝜎|) |𝜎|2; 1𝐺≪𝝓(|𝜎|) ; 1𝐸≪𝝓(|𝜎|) (2.2.1) Esta simplificación se puede aplicar cuando 𝝓(|𝜎|) es suficientemente grande, ya que, al inicio de las simulaciones, en el cual los valores serán menores, no podemos despreciar el término elástico frente al término plástico. Sin embargo, conforme vayan avanzando las simulaciones, es bastante probable que podamos aplicar dicha simplificación. Simplificando la expresión (2.1.14) según las hipótesis de las ecuaciones (2.2.1): [𝑑𝜎 𝑑𝑡 𝑑𝜏 𝑑𝑡]=1 𝝓(|𝜎|) |𝜎|2(4𝜏2 𝐸+4𝜎2 9𝐺) [ 𝝓(|𝜎|)4𝜏2 |𝜎|2−𝝓(|𝜎|)43𝜎𝜏 |𝜎|2 −𝝓(|𝜎|)43𝜎𝜏 |𝜎|2𝝓(|𝜎|)49𝜎2 |𝜎|2 ] [𝑑𝜀 𝑑𝑡 𝑑𝛾 𝑑𝑡] (2.2.2) Simplificando: [𝑑𝜎 𝑑𝑡 𝑑𝜏 𝑑𝑡]=1 (𝜏2 𝐸+𝜎2 9𝐺)[𝜏2−13𝜎𝜏 −13𝜎𝜏 19𝜎2][𝑑𝜀 𝑑𝑡 𝑑𝛾 𝑑𝑡] (2.2.3) Una vez hemos conseguido simplificar bastante nuestro sistema de ecuaciones diferenciales, podemos ver que la matriz es singular. Esto nos lleva a un análisis muy interesante de dicha ecuación por componentes. Para ello, compararemos ambas componentes para relacionarlas con el módulo efectivo de tensiones: 13 13 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones 𝑑𝜎 𝑑𝑡=1 (𝜏2 𝐸+𝜎2 9𝐺)[𝜏2𝑑𝜀 𝑑𝑡−13𝜎𝜏𝑑𝛾 𝑑𝑡]=𝜏 (𝜏2 𝐸+𝜎2 9𝐺)[𝜏𝑑𝜀 𝑑𝑡−13𝜎𝑑𝛾 𝑑𝑡] (2.2.4) 𝑑𝜏 𝑑𝑡=1 (𝜏2 𝐸+𝜎2 9𝐺)[−13𝜎𝜏𝑑𝜀 𝑑𝑡+19𝜎2𝑑𝛾 𝑑𝑡]=−13𝜎 (𝜏2 𝐸+𝜎2 9𝐺)[𝜏𝑑𝜀 𝑑𝑡−13𝜎𝑑𝛾 𝑑𝑡] (2.2.5) Si comparamos estas dos expresiones, podemos llegar a una expresión bastante curiosa, la cual relaciona las hipótesis de simplificación que se han llevado a cabo con el valor del módulo efectivo de tensiones a lo largo del ensayo: 𝑑𝜎 𝑑𝑡 𝑑𝜏 𝑑𝑡=−3𝜏 𝜎 (2.2.6) 𝜎𝑑𝜎 𝑑𝑡+3𝜏𝑑𝜏 𝑑𝑡=0 (2.2.7) 𝜎𝑑𝜎 𝑑𝑡=−3𝜏𝑑𝜏 𝑑𝑡 (2.2.8) Si calculamos la derivada del módulo de las tensiones según la métrica del espacio tensional, expuesto en la ecuación (1.3): 𝑑|𝜎| 𝑑𝑡 =𝑑√23𝜎2+2𝜏2 𝑑𝑡 (2.2.9) 𝑑|𝜎| 𝑑𝑡 =𝑑√23𝜎2+2𝜏2 𝑑𝜎 𝑑𝜎 𝑑𝑡+𝑑√23𝜎2+2𝜏2 𝑑𝜏 𝑑𝜏 𝑑𝑡 (2.2.10) 𝑑|𝜎| 𝑑𝑡 =− 1 2√23𝜎2+2𝜏243𝜎𝑑𝜎 𝑑𝑡−1 2√23𝜎2+2𝜏24𝜏𝑑𝜏 𝑑𝑡 (2.2.11) Podemos sustituir en la expresión (2.2.11) la ecuación (2.2.8) Consideraciones teóricas del modelo de cálculo 14 𝑑|𝜎| 𝑑𝑡 =1 2√23𝜎2+2𝜏2433𝜏𝑑𝜏 𝑑𝑡−1 2√23𝜎2+2𝜏24𝜏𝑑𝜏 𝑑𝑡 (2.2.12) Por lo tanto: 𝑑|𝜎| 𝑑𝑡 =0 (2.2.13) Esto implica que llega un momento en el cual a pesar de que seguimos teniendo ciclos de carga, el módulo de las tensiones se tiene que estabilizar en un valor. En este momento, el vector tensión describiría una trayectoria circular cuyo diámetro es el valor del módulo de sigma |𝜎|. 15 3 PARAMETRIZACIÓN DEL MODELO MATEMÁTICO EN MATLAB ara determinar si las previsiones que se han determinado en el apartado anterior y las hipótesis utilizadas se cumplen, y si es así, cuándo y en qué valor se estabiliza |𝜎|,realizaremos una simulación en MATLAB integrando las ecuaciones diferenciales y representando los resultados gráficamente para analizarlos en detalle. 3.1 Consideraciones generales Para definir el tramo de aumento del módulo del vector deformación |𝜀|, podemos imponer diferentes trayectorias. Para este caso de estudio, utilizaremos las dos simplificaciones que hemos considerado más ilustradoras para el análisis de nuestro problema. 3.1.1 Aproximación lineal Por un lado, podemos afirmar que el módulo de las deformaciones |𝜀| se incrementará linealmente, y sólo debido a la componente longitudinal 𝜀. En este caso, en el material el módulo del vector deformación alcanzaría el valor de |𝜀𝑚𝑎𝑥| siguiendo el mismo procedimiento que en el caso monoaxial, produciendo a su vez un incremento lineal de la componente axial 𝜎 del vector tensión. Una vez alcanzado el punto |𝜀𝑚𝑎𝑥|, el vector deformación empezaría a dibujar una trayectoria circular, mediante la cual empezarían a aparecer deformaciones transversales 𝛾 y tensiones tangenciales 𝜏. En la siguiente ilustración podemos verlo con más detalle: Figura 3.1.1 Las expresiones del vector deformación 𝜺(𝒕) para este procedimiento de carga vienen dadas por: 𝜀(𝑡)=√23|𝜀𝑚𝑎𝑥|𝑡 𝑡𝑐𝑎𝑟𝑔𝑎 (3.1.1) P Parametrización del modelo matemático en MATLAB 16 𝜏(𝑡)=0 (3.1.2) Siendo 𝑡 el parámetro de tiempo que va desde 𝑡=0 hasta 𝑡𝑐𝑎𝑟𝑔𝑎, instante en el cual el vector deformación 𝜺 alcanza el valor |𝜀𝑚𝑎𝑥|. Para entender mejor los cálculos que haremos posteriormente, 𝑡𝑐𝑎𝑟𝑔𝑎 será múltiplo de 2𝜋. 3.1.2 Aproximación en espiral Por el otro lado, la carga del material se puede modelar hacia la estabilización en la superficie |𝜀𝑚𝑎𝑥|=𝑐𝑡𝑒 mediante una aproximación en espiral. La forma más sencilla de modelar este tipo de aproximación es mediante una circunferencia que crece linealmente con el parámetro 𝑡. Mediante este método, desde el primer instante el material está experimentando un proceso de carga multiaxial, tanto en tensiones como en deformaciones, a diferencia del caso anterior. Después de alcanzar el valor del módulo del vector de deformaciones |𝜀𝑚𝑎𝑥|, proseguiríamos como en el caso anterior, describiendo una circunferencia de radio constante igual a |𝜀𝑚𝑎𝑥|. Podemos representarlo gráficamente como: Figura 3.1.2 Las expresiones del vector deformación 𝜺 impuesto para este caso particular, serán de la forma: 𝜀(𝑡)=√23|𝜀𝑚𝑎𝑥|cos (𝑡) 𝑡 𝑡𝑐𝑎𝑟𝑔𝑎 (3.1.2.1) 𝜏(𝑡)=√2|𝜀𝑚𝑎𝑥|sin (𝑡) 𝑡 𝑡𝑐𝑎𝑟𝑔𝑎 (3.1.2.2) Los valores constantes que multiplican las variables de las expresiones anteriores vienen dados para que representen una circunferencia en la métrica de nuestro espacio de deformaciones. Siendo el parámetro 𝑡 desde 𝑡=0 hasta 𝑡=𝑡𝑐𝑎𝑟𝑔𝑎 como en el caso anterior. Mirando las expresiones (3.1.2.1) y (3.1.2.2), podemos apreciar en este momento la utilidad de definir 𝑡𝑐𝑎𝑟𝑔𝑎 como múltiplo de 2𝜋. Así podemos definir cuantos ciclos tarda el material en alcanzar un vector deformación de módulo |𝜀𝑚𝑎𝑥|. 17 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones 3.1.3 Tramo con modulo de deformaciones constante Una vez definido el vector deformación 𝜺 para el tramo de aproximación, este pasa a describir una circunferencia de radio constante y valor |𝜀𝑚𝑎𝑥|, como vimos en las figuras anteriores. Para este tramo, que consideramos ya como el endurecimiento cíclico del material, las expresiones del vector deformación 𝜺 son las siguientes: 𝜀(𝑡)=√23|𝜀𝑚𝑎𝑥|cos (𝑡) (3.1.3.1) 𝜏(𝑡)=√2|𝜀𝑚𝑎𝑥|sin (𝑡) (3.1.3.2) En este caso el parámetro 𝑡 empieza desde el valor 𝑡=𝑡𝑐𝑎𝑟𝑔𝑎 hasta el valor 𝑡=𝑡𝑓𝑖𝑛. Para determinar con mayor facilidad el número de ciclos que va a durar nuestra simulación, definimos el instante final 𝑡𝑓𝑖𝑛 como 𝑡𝑓𝑖𝑛=2𝜋𝑁𝐶, siendo 𝑁𝐶 el número de ciclos totales que se producen en el ensayo. Para comparar los mismos resultados, utilizaremos en el transcurso de las simulaciones la aproximación mediante espiral, ya que al aproximarse al segundo tramo de una forma más progresiva, podrá soportar mejor los posibles errores numéricos que un cambio de pendiente tan brusco como ocurre si realizamos el primer tramo de forma lineal.. El problema puede estar mal condicionado, puesto que el valor de entrada son deformaciones (cuyo orden de magnitud es de 0.01) y el valor de salida son tensiones (orden de magnitud 500). Es por esto por lo que debemos utilizar tolerancias de integración pequeñas, aunque esto nos aumente el tiempo de ejecución del código. Se utilizarán tolerancias de integración absolutas y relativas de 10−12. Para el muestreo de los datos y su representación gráfica, tomaremos una tolerancia ℎ=10−4 para el parámetro 𝑡. Este valor permite una resolución gráfica óptima, y a su vez no sobrecarga en exceso la memoria del ordenador que realizará los cálculos. 3.2 Parámetros y coeficientes del material En esta simulación definiremos los valores de los parámetros siguientes para introducirlos en las ecuaciones diferenciales y expresiones descritas en los apartados anteriores. Estos son los parámetros Ramberg-Osgood, Módulo de Young y coeficiente de Poisson. Utilizaremos los atributos del material Aluminio 7075 T651 según la referencia [5] Para comprobar la utilidad de los datos utilizados también disponemos de los parámetros de la curva 𝜀−𝑁 asociada, según el método de las deformaciones locales, cuya ecuación es: 𝜀=𝑆𝑓 𝐸(2𝑁𝑐)𝑏+𝜀𝑓(2𝑁𝑐)𝑐 (3.2.1) Esta ecuación relaciona el número de ciclos de rotura según una deformación efectiva 𝜀 dada, siendo 𝑆𝑓 y 𝜀𝑓 , respectivamente, la tensión y deformación de fallo del material en el ensayo de tracción. Dichos valores, para el Aluminio 7075 T651 vienen dados por [5]: Parametrización del modelo matemático en MATLAB 18 K 852 MPa n 0,074 E 71 GPa G 26,7 GPa v 0,33 1231 MPa 0.26 b -0,122 c -0,806 Parámetros Ramberg-Osgood Parámetros del ensayo de tracción Parámetros Aluminio 7075 T651 Parámetros de la curva ε-N 19 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones 4 MÉTODOS DE CÁLCULO UTILIZADOS EN MATLAB eniendo en cuenta que estamos ante un problema numérico donde la matriz que define el sistema de ecuaciones diferenciales es posible que sea singular en algunos puntos, se establecerán un total de tres modelos de cálculo. Estos modelos surgirán adaptando al primer modelo las necesidades que encontremos a lo largo de las simulaciones. 4.1 Primer modelo de cálculo. Hipótesis de estado tensional creciente. Dado que, según las previsiones teóricas, el módulo efectivo de tensiones |𝜎| va creciendo y tiende a estabilizarse en un valor, podemos decir que este valor sólo crece durante la simulación, siendo posible usar sólo las ecuaciones de carga del material, las cuales son bastante más sencillas, lo que simplifica notablemente el problema. 4.1.1 Parametrización Para evaluar los resultados, haremos varias simulaciones, definiendo diferentes valores de |𝜀𝑚𝑎𝑥| como valor de entrada. Para definir este rango, haremos uso de la ecuación de la curva 𝜀−𝑁, definiendo un rango de ciclos 𝑁𝐶 entre 50-250 ciclos, ya que el método de las deformaciones locales y la teoría establecida en la introducción de este proyecto permite predecir el comportamiento del material a bajo número de ciclos (menos de 1000 ciclos). Este rango a su vez nos permitiría estar en el rango medio que abarca desde 𝑁𝐶=1 hasta 𝑁𝐶=1000 Utilizando la ecuación 𝜀−𝑁 mediante el script curvaen.m podemos establecer un rango de valores de |𝜀𝑚𝑎𝑥| entre |𝜀𝑚𝑎𝑥| = 0.011 y |𝜀𝑚𝑎𝑥|=0.017 según la tabla: Para realizar la simulación de este primer modelo (sólo ecuaciones de carga), utilizaremos el script espiralcl. Este script se puede utilizar para calcular mediante otros modelos de cálculo más avanzados que veremos en los siguientes apartados, por lo que para este primer modelo ecuación deberemos definir los siguientes parámetros en el script: • “ecuacion=1” para utilizar el primer modelo de cálculo. • “esp=1” para aproximar el primer tramo en espiral, el cual se utilizará durante toda la simulación. • “material=1” para tomar los parámetros del Aluminio 7075 T651 • “𝑁𝐶=100” para calcular 100 ciclos en cada iteración. Esto se cambiará según los resultados para un mejor muestreo. Es decir, para |𝜀𝑚𝑎𝑥|=0.014 evaluaremos 200 ciclos y para |𝜀𝑚𝑎𝑥|=0.017 1000 ciclos. • “|𝜀𝑚𝑎𝑥|=𝑉𝐴𝐿𝑂𝑅 𝐷𝐸 𝑇𝐴𝐵𝐿𝐴 𝐴 𝐸𝑉𝐴𝐿𝑈𝐴𝑅" 0.010 0.011 0.012 0.013 0.014 0.015 0.016 0.017 T Métodos de cálculo utilizados en MATLAB 20 Dentro de dicho script, para integrar el sistema de ecuaciones diferenciales (2.1.14), debemos hacerlo por dos tramos, como explicamos anteriormente. Uno de esos tramos para el aumento de |𝜀| en trayectoria espiral y otro para cuando el vector deformación describe una trayectoria circular según la circunferencia |𝜀|=|𝜀𝑚𝑎𝑥|, parámetro previamente definido. Para ello, disponemos de dos funciones, syst.m y syst2.m, (para los tramos 1 y 2, respectivamente), las cuales toman como entrada las componentes del vector deformación según las ecuaciones (3.1.2.1) y (3.1.2.2) para el primer tramo y las ecuaciones (3.1.3.1) y (3.1.3.2.) para el segundo tramo. A su vez, devuelven las componentes 𝜎 y 𝜏 del vector tensión, a lo largo de cada instante del vector parametrizado 𝑡. Estas funciones se integrarán según la función ode45 que incluye MATLAB. Al ejecutar el script con dichos parámetros, realizará una serie de representaciones, las cuales son: • Componente normal de tensión 𝜎 frente a componente normal de deformación 𝜀. (Ciclos de histéresis en componente normal). • Componente transversal de tensión 𝜏 frente a componente transversal de deformación 𝛾. (Ciclos de histéresis en componente transversal). • Componente normal de tensión 𝜎 frente a componente transversal de tensión 𝜏. (Evolución del vector tensión 𝝈) • Componente normal de deformación𝜀 frente a componente transversal de deformación 𝛾. (Evolución del vector tensión 𝜺) • Norma de tensiones efectivas |𝜎| frente al número de ciclos 𝑁𝐶 • Vector de diferencias de |𝜎| frente al número de ciclos 𝑁𝐶. Esto nos sirve para conocer si ha habido descarga, es decir, |𝜎| ha disminuido. Si esto ocurriera, tendríamos que utilizar las ecuaciones completas (carga y descarga). • Ángulo de desfase entre los dos vectores. En primera instancia, nos interesa sobre todo si los resultados van a ser estables numéricamente o no, y entonces analizar los distintos escenarios. En este primer modelo nos centraremos en la curva 𝑁𝐶−|𝜎|. Siguiendo los valores de |𝜀𝑚𝑎𝑥| dados por la tabla 4.1, los resultados gráficos más representativos se indican a continuación. Escogeremos como valores representativos tres valores de la tabla 4.1. Estos valores serán un valor de deformaciones bajas, como |𝜀𝑚𝑎𝑥|=0.011; un valor de deformaciones intermedio, el cual será |𝜀𝑚𝑎𝑥|=0.014 y un valor de deformaciones elevado, respectivamente al intervalo de deformaciones que hemos escogido, tomado |𝜀𝑚𝑎𝑥|=0.017 Los resultados para este primer método de cálculo se muestran a continuación. 21 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones 4.1.2 Resultado de las simulaciones 4.1.2.1 Resultados para módulo de deformaciones impuesto |𝜺𝒎𝒂𝒙|=𝟎.𝟎𝟏𝟏 Figura 4.1.2.1.1 Figura 4.1.2.1.2 Métodos de cálculo utilizados en MATLAB 28 4.1.2.3 Resultado para módulo de deformaciones impuesto |𝜺𝒎𝒂𝒙|=𝟎.𝟎𝟏𝟕 Figura 4.1.2.3.1 Figura 4.1.2.3.2 29 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones Figura 4.1.2.3.3 Figura 4.1.2.3.4 Métodos de cálculo utilizados en MATLAB 30 Figura 4.1.2.3.5 Figura 4.1.2.3.6 31 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones Figura 4.1.2.3.7 4.1.3 Análisis de los resultados Procedemos a analizar las representaciones gráficas que nos da el primer modelo de cálculo, en el cual sólo se implementan las ecuaciones de carga del material, las cuales son válidas cuando el módulo del vector tensión |𝜎| aumenta, es decir, 𝑑|𝜎|>0. Es por esto por lo que las figuras que representan las diferencias divididas de |𝜎| será crucial, ya que si vemos que en algún momento |𝜎| decrece, (vector de diferencias divididas negativo) la validez de los resultados sería muy cuestionable. En principio, es interesante analizar las figuras 4.1.2.1.1, 4.1.2.1.2, 4.1.2.2.1, 4.1.2.2.2, 4.1.2.3.1 y 4.1.2.3.2, las cuales representan los lazos histéresis tensión-deformación para las componentes axiales y tangenciales, respectivamente. Parece ser que, viendo las figuras de los ciclos de histéresis, en la medida que se incrementa el número de ciclos, el material se va endureciendo cada vez más, por lo que el área que encierra el lazo de histéresis, cada vez se va estrechando más. Esto podemos apreciarlo mejor si utilizamos el comando comet en MATLAB para la representación de los resultados, en vez de utilizar el comando plot, lo que nos muestra una evolución en tiempo real de cómo se construyen los lazos de histéresis. Llega un momento en el que los lazos se estrechan tanto que el material parece que colapsa, pues se queda sin margen de plastificación. Este momento, podría ser cuando el valor de |𝜎| alcanza su máximo, según las previsiones teóricas vistas en las hipótesis del apartado 3 de esta memoria. Sin embargo, si observamos las figuras 4.1.2.1.6, 4.1.2.2.6 y 4.1.2.3.6, las cuales representan las diferencias divididas para los tres valores del módulo de deformación impuesto, respectivamente, podemos observar que el módulo efectivo de tensiones |𝜎| decrece, lo que nos hace dudar de la validez de los cálculos realizados. Esto se puede observar a su vez en las figuras 4.1.2.1.5, 4.1.2.2.5 y 4.1.2.3.5, las cuales muestran la evolución de |𝜎| frente al número de ciclos 𝑁𝑐. Se ve claramente un descenso a partir de cierto punto del valor de |𝜎|, desde el cual las ecuaciones se comportan de forma extraña. Cabe destacar que el comportamiento de las ecuaciones difiere según cambia el módulo de deformación impuesta |𝜀𝑚𝑎𝑥|. Se observa un comportamiento más inestable o “turbulento” en el caso de deformaciones más bajas, comportándose nuestro modelo de cálculo de forma más estable o “laminar” en el caso de deformaciones mayores. Podemos establecer que el caso |𝜀𝑚𝑎𝑥|=0.014 podría ser el punto de transición entre un tipo de comportamiento y el otro. Métodos de cálculo utilizados en MATLAB 32 Para comprobar qué ocurre en el momento en que las curvas de histéresis se estrechan al límite de colapsar, hemos de calcular sólo hasta el número de ciclos donde |𝜎| alcanza su valor máximo, verificando así si este es el punto en el que los lazos de histéresis colapsan. Representando para los ciclos de histéresis en coordenadas axiales (𝜀−𝜎), nos darían los siguientes resultados: Para |𝜀𝑚𝑎𝑥|=0.011: Figura 4.1.3.1 Para |𝜀𝑚𝑎𝑥|=0.014: Figura 4.1.3.2 A deformaciones impuestas superiores, como en el caso de |𝜀𝑚𝑎𝑥|=0.017, la progresión hasta cerrar el ciclo de histéresis es tan lenta que las representaciones no pueden distinguirse de las otras sin hacer zoom, o sin utilizar 33 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones el comando 𝑐𝑜𝑚𝑒𝑡. Hay que apreciar también que según la deformación impuesta |𝜀𝑚𝑎𝑥| es mayor, mayor número de ciclos son necesarios para que |𝜎| alcance el valor máximo. Se ve claramente que hay inestabilidades numéricas, que van cambiando según cambia el valor de |𝜀𝑚𝑎𝑥| impuesto. Para valores más pequeños, las fluctuaciones son rápidas y de amplitud pequeña. Sin embargo, para valores más altos, las fluctuaciones tienen una amplitud considerable, aunque tengan una frecuencia menor. Esto podría poner en duda la validez de los resultados obtenidos. Sin embargo, estas fluctuaciones ocurren después de que |𝜎| sea máximo. Al estar utilizando las ecuaciones sólo teniendo en cuenta la carga, es posible que se deba a esto, pues después del pico de |𝜎|, su valor desciende. Teniendo sólo en cuenta el tramo en el que |𝜎| va aumentando, podemos afirmar que el punto en el que |𝜎| es máximo, es cuando el ciclo de histéresis no se puede cerrar más. Según la teoría propuesta, el valor de |𝜎| se debería de estabilizar entonces. ¿Qué ocurriría si el valor de |𝜎| se estabiliza? (𝑑|𝜎|=0) Al ser la deformación plástica 𝜀𝑝 directamente proporcional a 𝑑|𝜎| según la expresión 1.11, podemos decir que, en este caso, en el que |𝜎| se estabiliza, no hay deformación plástica. Si no hay deformación plástica, la evolución del material debe ser elástica, quedando el sistema de ecuaciones diferenciales dependiendo sólo del tensor de Hooke, sin el término plástico, cuya expresión sería: [𝑑𝜎 𝑑𝑡 𝑑𝜏 𝑑𝑡]=[𝐸 0 0 𝐺][𝑑𝜀 𝑑𝑡 𝑑𝛾 𝑑𝑡] (4.1.3) Pero al usar las ecuaciones sin tener en cuenta la descarga, tenemos un problema asociado a ello, pues después de que |𝜎| alcance el máximo, este valor comienza a decrecer en todos los ensayos realizados, en unos más rápido (a menores deformaciones impuestas), y en otros de forma más progresiva (a mayores deformaciones impuestas). Esto es totalmente incompatible con este modelo, pues a partir de este momento deberíamos utilizar las ecuaciones teniendo en cuenta la descarga (efecto memoria), descritas en la teoría de plasticidad propuesta al inicio del documento, que no se han implementado en este modelo. Debemos entonces, implementar las ecuaciones teniendo en cuenta la descarga, para verificar si efectivamente después de que |𝜎| alcance el máximo, este valor se estabiliza, como hemos previsto con las hipótesis anteriormente señaladas. 4.2 Segundo modelo de cálculo. Utilización de ecuaciones de descarga En el apartado anterior hemos podido ver que nuestros resultados no podrían tener la validez necesaria, puesto que hemos utilizado un modelo de cálculo suponiendo que |𝜎| es siempre creciente, es decir, 𝑑|𝜎|>=0 y vemos que esto se incumple en las simulaciones que llevamos hasta ahora. Debemos entonces calcular las tensiones como distancias respecto a la circunferencia de plastificación 𝐶𝑘 que se genere, como hemos explicado detalladamente en el apartado introductorio de este proyecto. Pero esto tiene una problemática asociada. Para que este método funcione, en nuestro sistema de ecuaciones diferenciales en derivadas parciales, debemos detectar cuando existe un paso de integración en el que 𝑑|𝜎|<= 0 y a partir de este, almacenar este paso como primer punto de descarga 𝜎0 y calcular su circunferencia de plastificación asociada 𝐶0. Hasta aquí no hay problema, pero ¿qué ocurriría si se almacenaran sucesivos puntos de descarga, y una vez que respecto a estas circunferencias de plastificación |𝜎| vuelva a crecer y estas se borren debido al efecto memoria? Pues que los sistemas de integración optimizados que se incluyen en el paquete MATLAB no permiten volver a pasos de integración anteriores, por lo que no podríamos almacenar los puntos de descarga y seguir integrando, para volver a ellos cuando se diera el efecto memoria. Para poder implementar este complejo sistema de ecuaciones, hay una alternativa, establecer un método de integración manual paso a paso siguiendo el método Runge-Kutta de orden 4. Métodos de cálculo utilizados en MATLAB 34 Para ello se crearán las funciones syst2rk4d.m y syst2d.m, que sustituirán respectivamente a las funciones ode45.m de MATLAB y syst2.m que hemos utilizado para calcular el segundo tramo de nuestras ecuaciones en el modelo de cálculo anterior. Cabe destacar, que, en el primer tramo de cálculo, es decir, el tramo de carga en el que |𝜀| se aproxima a |𝜀𝑚𝑎𝑥| mediante una espiral, no es necesario establecer estas ecuaciones ya que en este tramo |𝜎| siempre crece y por lo tanto 𝑑|𝜎| es siempre positivo, siendo el modelo de cálculo anterior efectivo en este caso. 4.2.1 Parametrización de modelo de integración paso a paso mediante método Runge-Kutta de cuarto orden El método Runge-Kutta de cuarto orden nos permite establecer un método de integración paso a paso de una forma más exacta y optimizada que el conocido método de Euler. Explicado brevemente, consiste en establecer las siguientes expresiones para calcular cada paso de integración: 𝑌(𝑖+1)=𝑌(𝑖)+ℎ(𝐾1+2𝐾2+2𝐾3+𝐾4) (4.2.1) 𝐾1=𝑑𝑌(𝑋(𝑖),𝑌(𝑖)) (4.2.2) 𝐾2=𝑑𝑌(𝑋(𝑖)+ℎ2,𝑌(𝑖)+𝐾1ℎ2) (4.2.3) 𝐾3=𝑑𝑌(𝑋(𝑖)+ℎ2,𝑌(𝑖)+𝐾2ℎ2) (4.2.4) 𝐾3=𝑑𝑌(𝑋(𝑖)+ℎ,𝑌(𝑖)+𝐾3ℎ) (4.2.5) Expresiones en las que 𝑌 es la función incognita a obtener, 𝑖 es el paso del ultimo valor de 𝑌 obtenido, ℎ es el paso de integración establecido, y 𝑋 constituye el conjunto de variables independientes de las que depende 𝑌. Mediante este método podremos acceder a todos los puntos de integración y poder así almacenar los sucesivos puntos de descarga 𝜎𝑘 que se vayan calculando y sus circunferencias de plastificación asociadas 𝐶𝑘, calculando a su vez el centro de estas circunferencias 𝜎𝑐,𝑘 desde el cual se medirán relativamente las distancias en el espacio tensional para calcular los siguientes puntos. ¿Cómo podemos implementar entonces las ecuaciones teniendo en cuenta la descarga? Explicaremos brevemente el algoritmo utilizado en el fichero syst2rk4d.m: • Se establecen primero los valores iniciales de centros de circunferencias de plastificación (empieza en el [0,0]), índice de paso 𝑖 en 1, y el parámetro 𝑐 en 1 también, parámetro que nos ayudará a contabilizar más adelante cuál es la circunferencia de plastificación activa. • Empezamos a utilizar las ecuaciones sin tener en cuenta la descarga, integrando nuestro sistema basado en RK4 introduciendo la función syst2.m como 𝑑𝑌. Se establecerá un condicional al final de cada paso, que verificará si |𝜎| en el último paso calculado ha aumentado respecto al paso anterior. Si esto ocurre, seguiremos integrando con estas ecuaciones. Sin embargo, si ocurre lo contrario, establecemos el parámetro 𝑑 con valor 1, parámetro respecto a qué punto de descarga medimos las distancias. Al tomar valor 1, se crea el primer punto de descarga. • Si al calcular un nuevo paso de integración el parámetro 𝑑 es distinto de 0, volveremos al punto de integración anterior, calcularemos la distancia 𝑞0, respecto a la cual se calculará el siguiente punto, mediante las ecuaciones (1.23), (1.24) y (2.1.14), sustituyendo en esta última la función de endurecimiento 𝜙(|𝝈|) por 𝜙(𝑞). Esto se implementa en MATLAB introduciendo como función a integral la función syst2d.m, que toma como valores de entrada 𝝈𝑐,𝑘 , el último valor de 𝝈, y el punto de descarga activo 𝝈𝑘. 35 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones • Si este nuevo valor de |𝜎| es mayor que el anterior, se produce efecto memoria y borramos el punto de descarga, estableciendo de nuevo el parámetro 𝑑 en 0. Volvemos a utilizar las ecuaciones sin tener en cuenta la descarga, según el primer modelo de cálculo, mediante la función syst2.m. • Esto cambia radicalmente si se produce más de un punto de descarga sucesivo, pues entonces se crearía la primera circunferencia de plastificación asociada 𝐶0 y su centro 𝜎𝑐,0, respecto al que se calculara entonces la distancias 𝑞 y 𝑞0, según la expresión 1.26. Cabe destacar 𝑞𝑘 es igual al diámetro de la cirfuncerencia de plastificación 𝐶𝑘 activa. • A partir de entonces, se establece un bucle repetitivo por el cual si las distancias respecto al centro de la circunferencia de plastificación activa 𝜎𝑐,𝑘 se reduce, se crearán nuevos puntos de descarga y se incrementarán los valores de 𝑑 y de 𝑐. Cada vez que aumentan estos parámetros, se almacenan los puntos de descarga y centros correspondientes. Al contrario, si las distancias aumentan y superan el radio de la circunferencia de plastificación activa, esta y su punto de descarga asociado se borrarán, tomándose como referencia el punto de descarga y su circunferencia de plastificación anteriores. Para entender mejor esto podemos examinar detenidamente las figuras 1.1, 1.2, 1.3, y 1.4 Para mayor detalle, los cálculos están perfectamente descritos y comentados en el fiche. Tenemos entonces implementado nuestro sistema para calcular de forma más precisa nuestras simulaciones, ya que en este caso sí hemos tenido en cuenta que el material puede descargar, por lo que si vemos en las figuras que |𝝈| desciende en algún instante, los resultados no tendrían por qué perder validez. Pero no todo son ventajas, pues la exactitud de nuestros resultados dependerá de la tolerancia de integración ℎ que escojamos. Aun así, podremos establecer como valor límite una tolerancia de integración para este método paso a paso de ℎ=10−4 ya que el costo computacional de almacenar variables con pasos de integración (y por ende de muestreo) pequeños es altísimo, tanto en consumo de memoria como en tiempo de ejecución. Como simple comparativa, el fichero espiral.cl tarda en ejecutar las simulaciones del primer modelo de cálculo unos 20 segundos. En este modelo de cálculo esta cifra asciende hasta horas incluso si mantenemos el paso de integración en ℎ=10−4. Veremos más delante si hay alguna forma de corregir esto. Estableceremos las simulaciones para los mismos valores que en el apartado 4.1, representando las figuras para |𝜀𝑚𝑎𝑥|=0.011, |𝜀𝑚𝑎𝑥|=0.014 y |𝜀𝑚𝑎𝑥|=0.017. Para utilizar este modelo de cálculo, configuraremos el fichero escpiral.cl de la siguiente forma: • “ecuacion=1” para utilizar el primer modelo de cálculo. • “esp=1” para aproximar el primer tramo en espiral, el cual se utilizará durante toda la simulación. • “material=1” para tomar los parámetros del Aluminio 7075 T651 • “𝑁𝐶=100” para calcular 100 ciclos en cada iteración. Esto se cambiará según los resultados para un mejor muestreo. Es decir, para |𝜀𝑚𝑎𝑥|=0.014 evaluaremos 200 ciclos y para |𝜀𝑚𝑎𝑥|=0.017 1000 ciclos. • “|𝜀𝑚𝑎𝑥|=𝑉𝐴𝐿𝑂𝑅 𝐷𝐸 𝑇𝐴𝐵𝐿𝐴 𝐴 𝐸𝑉𝐴𝐿𝑈𝐴𝑅" El resto de la ejecución que no tenga que ver con la integración de las variables se hará análoga al procedimiento del primer modelo de cálculo descrito en el apartado 4.1 Métodos de cálculo utilizados en MATLAB 36 4.2.2 Resultados de las simulaciones 4.2.2.1 Resultado para módulo de deformaciones impuesto |𝜺𝒎𝒂𝒙|=𝟎.𝟎𝟏𝟏 Figura 4.2.2.1.1 Figura 4.2.2.1.2 37 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones Figura 4.2.2.1.3 Figura 4.2.2.1.4 Métodos de cálculo utilizados en MATLAB 44 Figura 4.2.2.3.3 Figura 4.2.2.3.4 45 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones Figura 4.2.2.3.5 Figura 4.2.2.3.6 Métodos de cálculo utilizados en MATLAB 46 Figura 4.2.2.3.7 4.2.3 Análisis de los resultados Una vez implementadas las ecuaciones de plasticidad teniendo en cuenta la existencia de puntos de inversión y del efecto memoria, hemos obtenido los resultados según las figuras (4.2.2.1.1) hasta (4.2.2.3.7), manteniendo las mismas condiciones y parámetros de las simulaciones realizadas con el primer modelo de cálculo, que no tenía estos factores. Esto nos permite comparar ambas figuras para determinar el comportamiento de ambos modelos. El factor más importante de los resultados obtenidos con este modelo es que ahora los resultados no tienen por qué invalidarse por ver un descenso en el módulo efectivo de tensiones |𝜎|, o apreciar que en las figuras de diferencias divididas de |𝜎| los valores son negativos. Aunque las previsones teóricas vistas en los apartados 1 y 2 de este documento expresan que |𝜎| no puede descender debido a que la deformación plástica se mantiene constante; esta vez sí es posible permitirnos pequeñas oscilaciones numéricas. Debemos tener en cuenta que el sistema de ecuaciones diferenciales indicado en la expresión (2.1.14) es complejo y por ello puede estar no exento de comportamientos numéricos indeseados. En los resultados del primer modelo de cálculo podemos apreciar una inestabilidad y un comportamiento errático de las ecuaciones una vez el módulo efectivo de tensiones |𝜎| llega a su máximo y se detecta el primer descenso de este valor. Este descenso hacía invalidar los resultados puesto que dichas ecuaciones no podían calcular verazmente a partir de ahí. Esto cambia con este modelo de cálculo. Ahora puede verse en las figuras (4.2.2.1.5), (4.2.2.2.5) y (4.2.2.3.5) que indican la evolución del valor |𝜎| respecto al número de ciclos totales transcurridos 𝑁𝑐, que el valor de |𝜎| llega a estabilizarse, con oscilaciones debido a un comportamiento numérico indeseado. Estas oscilaciones se producen tanto para deformaciones bajas como para deformaciones altas. Sin embargo, esto no invalida los resultados al haber tenido en cuenta las ecuaciones de descarga. Hay que tener en cuenta que hemos definido un paso de integración no muy bajo ℎ=10−4 debido al elevado consumo de recursos que requiere implementar un modelo de integración paso a paso para poder almacenar los 47 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones sucesivos puntos de descarga 𝜎𝑘 y sus circunferencias de plastificación asociadas 𝐶𝑘. Este paso de integración puede ser el responsable de las oscilaciones numéricas, siendo necesaria una simplificación de estas ecuaciones, para poder realizar un análisis con mayor profundidad. 4.3 Tercer modelo de cálculo. Simplificación de ecuaciones de descarga Volviendo al apartado 4.1.3, si verificamos que |𝜎| se estabiliza, lo que ocurre, pero con oscilaciones de naturaleza numérica; no existe deformación plástica, puesto que según la expresión 1.11 vemos que el incremento de deformación plástica 𝑑𝜀𝑝 es directamente proporcional al incremento del módulo efectivo de tensiones 𝑑|𝜎|, el cual es 0 a partir de la estabilización de |𝜎|. Entonces podríamos decir que a partir de la estabilización de |𝜎| el comportamiento es elástico al no existir plastificación. Las figuras (4.2.2.1.1), (4.2.2.1.2), (4.2.2.2.1), (4.2.2.2.2), (4.2.2.3.1), (4.2.2.3.2), nos indican que la órbita de estabilización es casi una línea recta, lo cual sería un comportamiento plástico perfecto. Esto se aprecia mejor utilizando el comando comet para las figuras en vez del comando plot, lo que nos permite apreciar el trazado de las figuras en tiempo real. Recordemos que, en el primer modelo de cálculo, el máximo de |𝜎| es el punto donde los ciclos de histéresis se estrechan al máximo y a partir de ahí |𝜎| comienza a descender. Esto sugiere la siguiente simplificación, con la cual no necesitaríamos utilizar las ecuaciones de descarga, en las cuales no podemos introducir un paso de integración más pequeño y en las cuales el tiempo de ejecución y un consumo de recursos es altísimo. Teniendo en cuenta esto último, podemos utilizar el primer modelo de cálculo, pero con un condicionante: puesto que |𝜎| comenzará a descender por comportamientos numéricos indeseados una vez alcandado el máximo, que coincide con el valor de estabilización; una vez se detecte que 𝑑|𝜎|<0 en un paso de integración, sustituiremos el sistema de ecuaciones (2.1.14) por el (4.1.3), el cual indica un comportamiento elástico, sólo para este paso de integración. La diferencia de esto con el segundo método de cálculo es que este pequeño algoritmo no necesita un sistema de integración paso a paso a mano, sino que se puede realizar utilizando las subrutinas de integración que nos ofrece el paquete MATLAB, las cuales están perfectamente optimizadas para dar unos resultados rápidos y precisos, como es el caso de la función ode45 utilizada en el primer modelo de cálculo. Esto, aparentemente, nos permitirá obtener unos resultados con mejor comportamiento numérico y con un tiempo de ejecución indiscutiblemente más rápido. 4.3.1 Parametrización La parametrización de este modelo es muy sencilla una vez definida la del primer modelo de cálculo, ya que se basa estrictamente en la de este. Sólo debemos introducir un pequeño cambio, en el script que definimos el sistema de ecuaciones diferenciales (2.1.14), syst.m, y syst2.m para el primer modelo, introduciremos la condición de que si 𝑑|𝝈| es negativo en algún paso de integración, en dicho paso de integración sustituiremos las ecuaciones (2.1.14) por las de la expresión (4.1.3), siendo este paso de integración puramente elástico. Se crean a su vez los ficheros systcl.m y syst2cl.m en el que se incluye este condicionante. Utilizaremos la subrutina de integración ode45 la cual trabajará con estas dos funciones, para el primer y segundo tramos, respectivamente. Estableceremos las simulaciones para los mismos valores que en el apartado 4.1, para |𝜀𝑚𝑎𝑥|=0.011 , |𝜀𝑚𝑎𝑥|=0.014 y |𝜀𝑚𝑎𝑥|=0.017. Para utilizar este modelo de cálculo, configuraremos el fichero espiral.cl de la siguiente forma: • “ecuacion=3” para utilizar el primer modelo de cálculo. • “esp=1” para aproximar el primer tramo en espiral, el cual se utilizará durante toda la simulación. • “material=1” para tomar los parámetros del Aluminio 7075 T651 • “𝑁𝐶=100” para calcular 100 ciclos en cada iteración. Esto se cambiará según los resultados para un mejor muestreo. Es decir, para |𝜀𝑚𝑎𝑥|=0.014 evaluaremos 200 ciclos y para |𝜀𝑚𝑎𝑥|=0.017 1000 ciclos. • “|𝜀𝑚𝑎𝑥|=𝑉𝐴𝐿𝑂𝑅 𝐷𝐸 𝑇𝐴𝐵𝐿𝐴 𝐴 𝐸𝑉𝐴𝐿𝑈𝐴𝑅" Métodos de cálculo utilizados en MATLAB 48 Los resultados que se representarán serán análogos a los vistos anteriormente en los apartados 4.1 y 4.2, para poder hacer una comparación efectiva sobre estos. 4.3.2 Resultados de las simulaciones 4.3.2.1 Resultado para módulo de deformaciones impuesto |𝜺𝒎𝒂𝒙|=𝟎.𝟎𝟏𝟏 Figura 4.3.2.1.1 49 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones Figura 4.3.2.1.2 Figura 4.3.2.1.3 Métodos de cálculo utilizados en MATLAB 50 Figura 4.3.2.1.4 Figura 4.3.2.1.5 51 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones Figura 4.3.2.1.6 Figura 4.3.2.1.7 Métodos de cálculo utilizados en MATLAB 52 4.3.2.2 Resultados para módulo de deformaciones impuesto |𝜺𝒎𝒂𝒙|=𝟎.𝟎𝟏𝟒 Figura 4.3.2.2.1 Figura 4.3.2.2.2 53 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones Figura 4.3.2.2.3 Figura 4.3.2.2.4 61 5 ANÁLISIS DE RELACIÓN ENTRE VARIABLES oncluido este análisis preliminar, se asegura que se cumplen las previsiones teóricas. Sin embargo, hay un pequeño detalle a tener en cuenta. Cuanto mayor es la deformación impuesta |𝜀𝑚𝑎𝑥|, mayor es el número de ciclos 𝑁𝑐 necesario para la estabilización de |𝜎|. Sería interesante establecer una relación de cómo varía este valor 𝑁𝑐 de establización según la deformación máxima impuesta |𝜀𝑚𝑎𝑥| en nuestro modelo. A partir de este punto utilizaremos el tercer modelo de simulaciones, visto su buen comportamiento numérico y acorde a las previsiones teóricas del apartado 2 5.1 Evolución del módulo de tensiones estabilizado frente a las deformaciones impuestas En este apartado veremos cómo varía el módulo de tensiones efectivo |𝜎| en función de las deformaciones impuestas constantes |𝜀𝑚𝑎𝑥| que definamos en la simulación. Para ello crearemos un nuevo fichero en MATLAB apoyándonos en el script principal de simulación espiral.cl. En dicho script utilizaremos los métodos de integración del tercer modelo de cálculo, el cual se describe en el apartado 4.3, ya que son los más adecuados y precisos para realizar las simulaciones, además de ser bastante más rápidos de ejecutar que los del segundo modelo de cálculo descrito en el apartado 4.2. El programa para este cometido es el fichero ordenadatos.m. Consiste en el siguiente algoritmo: • Utilizamos la estructura base del programa espiral.cl, configurando los parámetros exactamente igual que en el tercer modelo de cálculo, fijando el parámetro “ecuación” en 3. • Se definirá un parámetro ℎ𝑒 el cual nos servirá para ir iterando y almacenando distintas simulaciones con diferentes valores de |𝜀𝑚𝑎𝑥|. ℎ𝑒 será el incremento de |𝜀𝑚𝑎𝑥| entre una simulación y otra. • Se definen dos parámetros, |𝜀𝑚𝑎𝑥|𝑖 y |𝜀𝑚𝑎𝑥|𝑓, que definirán el valor inicial y final de |𝜀𝑚𝑎𝑥| de las simulaciones a almacenar; cuyos valores serán 0.075 y 0.0175, respectivamente. • En cada iteración, se almacenarán los valores 𝑄 y 𝑁𝑒𝑠𝑡, que representarán el valor en el cual se estabiliza el módulo de tensiones efectivo |𝜎| y el número de ciclos necesario para que se estabilize, respectivamente. • Se almacena esta información en un archivo de texto, tabl.txt, para su posterior tratamiento, ya que al realizar bastantes simulaciones, el tiempo de ejecución es alto. Una vez ejecutado podemos representar los datos obtenidos con las siguientes figuras: C Análisis de relación entre variables 62 62 Figura 5.1.1 Figura 5.1.2 63 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones 5.2 Ajustes gráficos Podemos observar claramente que la relación entre el valor estabilizado de |𝜎| y las deformaciones impuestas |𝜀𝑚𝑎𝑥| es lineal. Además, el número de ciclos necesario para la estabilización 𝑁𝑒𝑠𝑡 parece ser que aumenta exponencialmente a medida que aumenta |𝜀𝑚𝑎𝑥|. Para comprobar que las evoluciones de estas dos relaciones son lineal y exponencial, respectivamente, se utilizará la función ajuste.m. Esta función comparará los resultados gráficos reales con unas gráficas generadas mediante ajustes de regresión. Esto puede implementarse fácilmente gracias a la función incluida en MATLAB lsqcurvefit, la cual recibe una función deseada, y los vectores de las variables independientes y dependientes, devolviendo los coeficientes de la función introducida. Dicha función con estos coeficientes representa el ajuste más cercano a los puntos de muestreo. Esto se realizará con una función lineal para la gráfica 𝑄−|𝜀𝑚𝑎𝑥| del tipo: 𝑓(𝑥)=𝑐1+𝑐1𝑥 (5.2.1) y para la gráfica 𝑁𝑒𝑠𝑡−|𝜀𝑚𝑎𝑥| una función exponencial del tipo: 𝑓(𝑥)=𝑐1𝑐2𝑥 (5.2.2) Siendo en ambos casos 𝑐1 y 𝑐2 los coeficientes que calcurará lsqcurvefit para los datos recogidos en el fichero tabl.txt, que se deben calcular anteriormente con el ejecutable ordenadatos.m También podremos comparar ambos datos y calcular el error relativo máximo que se comete con el ajuste. Los resultados se ven a continuación: Figura 5.2.1 Análisis de relación entre variables 64 64 Figura 5.2.2 Los errores relativos máximos son en el caso del ajuste lineal representado en la figura 5.2.1 𝑒𝑚𝑎𝑥=0,53 %, y en el caso del ajuste lineal 𝑒𝑚𝑎𝑥=0,91 %. Estos valores son muy pequeños, por lo que vemos que el ajuste es excelente en ambos casos. Los coeficientes obtenidos para ambas relaciones son los siguientes, teniendo cuenta los parámetros para el Aluminio 7075 T651 según la referencia [5]: 𝑐1𝑙𝑖𝑛=47.000 ; 𝑐2𝑙𝑖𝑛=45 𝑐1𝑒𝑥𝑝=0 ; 𝑐2𝑒𝑥𝑝=1.727 (5.2.3) Siendo 𝑐1𝑙𝑖𝑛 y 𝑐2𝑙𝑖𝑛 los parametros para el ajuste líneal de la relación 𝑄−|𝜀𝑚𝑎𝑥| y 𝑐1𝑒𝑥𝑝 y 𝑐2𝑒𝑥𝑝 los parametros para el ajuste exponencial de la relación 𝑄−𝑁𝑒𝑠𝑡, 5.3 Predicciones para deformaciones plásticas superiores 5.3.1 Consideraciones generales Si nos figamos más detenidamente en la figura 5.2.2 vemos que al verificarse que el número de ciclos necesarios para que el material se estabilice, 𝑁𝑒𝑠𝑡, evoluciona exponencialmente a deformaciones altas. Por ello, a deformaciones sensiblemente más altas, la estabilización se produciría tan tarde que se saldría del campo de la fatiga a bajo número de ciclos, ya que el valor 𝑁𝑒𝑠𝑡 superaría holgadamente los 1000 ciclos. Es por ello de utilidad examinar detenidamente el comportamiento a deformaciones más altas para estudiar si podemos realizar una predicción mejor que declarar que el valor de |𝜎| se estabilizará. 65 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones Para ello simularemos de nuevo con el tercer modelo de cálculo, pero esta vez estableceremos |𝜀𝑚𝑎𝑥| en un valor bastante mayor, por ejemplo, |𝜀𝑚𝑎𝑥|=0.05, evaluando 500 ciclos de carga. Esto nos permitirá ver si |𝜎| crece muy lentamente o si podemos aproximarlo con otro tipo de ajuste. Los resultados se muestran en el siguiente apartado. 5.3.2 Resultados de las simulaciones Figura 5.3.2.1 Análisis de relación entre variables 66 66 Figura 5.3.2.2 Figura 5.3.2.3 67 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones Figura 5.3.2.4 Figura 5.3.2.5 Análisis de relación entre variables 68 68 Figura 5.3.2.6 Figura 5.3.2.7 69 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones 5.3.3 Análisis de los resultados Empezaremos a analizar los ciclos o lazos de histéresis de las figuras 5.3.2.1 y 5.3.2.2, en comparación con lo que hemos visto anteriormente para las deformaciones impuestas menores. Vemos que los lazos de histéresis son anchos y que se estrechan a un ritmo bastante lento, lo que nos indica que la supuesta estabilización, si llegase, habría llegado mucho después de la rotura del material, por lo que dicho cálculo no nos sería muy útil. El ángulo entre los vectores tensión 𝝈 y deformación 𝜺 se acerca bastante a 90º, y evoluciona muy lentamente, diciéndonos también que faltan muchos ciclos para la estabilización. Esto nos dice también que el comportamiento plástico predomina todavía mucho en material. Podemos apreciar esta información en la figura 5.3.2.7 Por último, si nos fijamos en la figura 5.3.2.5 vemos que el comportamiento se asemeja más a una gráfica logarítmica que a una asíntota, que es a lo que tendría que asemejarse dicha figura para asegurar que |𝜎| se va a estabilizar en no muchos ciclos relativamente hablando. 5.3.4 Ajuste gráfico Análogamente al procedimiento del apartado 5.2, volveremos a utilizar el comando lsqcurvefit de MATLAB para realizar un ajuste a la figura 5.3.2.5 según una función logarítmica de la forma: 𝑓(𝑥)=𝑐1𝑙𝑜𝑔(𝑥−𝑐2)+𝑐3 (5.3.4.1) Siendo 𝑐1, 𝑐2 y 𝑐3 los coeficientes a buscar para representar la función del ajuste. Análogamente al apartado 5.2, calcularemos también el error relativo máximo para este ajuste. El ajuste se realizará para el segundo tramo, ya que el tramo de carga hasta que el módulo de deformación efectiva alcanza el valor de |𝜀𝑚𝑎𝑥| es casi una línea recta, lo que invalidaría el ajuste Los resultados se muestran a continuación: Figura 5.3.4 Anexo. Código de programación necesario para las simulaciones 76 76 %REPRESENTACIÓN DE FIGURAS% figure (1) plot(epsilon,(sigma(:,1))) title('\bf Curva de histéresis axial') xlabel('\bf \epsilon [Adimensional]') ylabel('\bf \sigma [MPa]') figure (2) plot(gamma,sigma(:,2)) title('\bf Curva de histéresis transversal') xlabel('\bf \gamma [Adimensional]') ylabel('\bf \tau [MPa]') figure (3) plot(sqrt(2/3)*sigma(:,1),sqrt(2)*sigma(:,2)) title('\bf Evolución del vector tensión (métrica de von Mises)') xlabel('\bf \sigma efectivo [MPa]') ylabel('\bf \tau efectivo [MPa]') figure (4) plot(sqrt(3/2)*epsilon,sqrt(1/2)*gamma) title('\bf Evolución del vector deformación (métrica de von Mises)') xlabel('\bf \epsilon efectivo [Adimensional]') ylabel('\bf \gamma efectivo [Adimensional]') figure (5) plot(N,Q) title('\bf Evolución del módulo efectivo de tensiones según el número de ciclos') xlabel('\bf N_{c} [Ciclos]') ylabel('\bf Norma efectiva de tensiones [MPa]') figure (6) plot(N(1:end-1),dm*h/(2*pi)) title('\bf Diferencias divididas del módulo efectivo de tensiones') xlabel('\bf N_{c} [Ciclos]') ylabel('\bf Vector de diferencias de \bf \sigma [MPa]') % figure (7) % % plot(Q*1e3,flujo(K,na,sigma')) % title('EVOLUCIÓN DEL FLUJO PLÁSTICO FRENTE A LA NORMA EFECTIVA DE TENSIONES') % % REPRESENTACIÓN EN TIEMPO REAL DE LAS GRAFICAS SIGMA-TAU, EPSILONGAMMA 77 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones % figure (1) % title('Tensión axial frente a tensión tangencial') % % hold on % figure (2) % title('Deformación normal frente a deformación transversal') % hold on % % % for t=1:length(tspan) % figure (1) % % plot(sqrt(2/3)*sigma(t,1),sqrt(2)*sigma(t,2),'.r') % figure (2) % % plot(sqrt(3/2)*epsilon(t),sqrt(1/2)*gamma(t),'.b') % % end %CÁLCULO DEL DESFASE ENTRE LOS VECTORES TENSIÓN Y DEFORMACIÓN% angeps=acos(sqrt(3/2)*epsilon(l1:end)/emax); %CÁLCULO DEL ÁNGULO GIRADO POR EL VECTOR DEFORMACIONES A PARTIR DEL SEGUNDO TRAMO% angsig=acos(sqrt(2/3)*sigma(l1:end,1)./Q(l1:end)')'; %CÁLCULO DEL ÁNGULO GIRADO POR EL VECTOR TENSIONES A PARTIR DEL SEGUNDO TRAMO% angsig(1)=0; %PONEMOS A CERO LOS PRIMEROS INSTANTES, POR SI FUERAN INFINITO% angsig2=dp(angsig); %CON UNA FUNCIÓN, PASAMOS LOS ÁNGULOS EN UNIDADES DE REFERENCIA, A ÁNGULO RECORRIDO TOTAL POR EL VECTOR% angeps(1)=0; angeps2=dp(angeps); %ESTO SE HACE PARA LOS DOS ÁNGULOS (DE LAS DEFORMACIONES Y DE LAS TENSIONES% %AHORA QUE TENEMOS LOS ÁNGULOS RECORRIDOS POR LOS DOS VECTORES, PODEMOS %CALCULAMOS EL DESFASE ENTRE AMBOS VECTORES desf=angeps2-angsig2; %REPRESENTACIÓN DEL DESFASE ENTRE VECTORES DEFORMACIÓN Y TENSIÓN% figure(15) plot(N(l1:end-1),desf*360/(2*pi)) title('\bf Ángulo de desfase entre vectores deformación y tensión') ylabel('\bf \alpha [º]') xlabel('\bf N_{c} [Ciclos]') Anexo. Código de programación necesario para las simulaciones 78 78 Syst.m function dsigma=syst(t,sigma,eps,esp) global K na E G %FUNCIÓN DE INTEGRACIÓN DE ECUACIONES DIFERENCIALES SIMPLIFICADA PARA EL %TRAMO1 Q=sqrt(2/3*sigma(1)^2+2*sigma(2)^2); %SE DEFINE EL MÓDULO EFECTIVO DE TENSIONES% aux=((sqrt(3/2))^(1/na+1))/na*((1/K)^(1/na))*(Q^(1/na-3)); %SE DEFINE UNA VARIABLE AUXILIAR= PHI(Q)/Q^2 PARA REDUCIR POSIBLES PROBLEMAS NUMÉRICOS% %SE DEFINE LA MATRIZ DE LAS DERIVADAS PARCIALES% A=1/(1/E*1/G+(aux)*((4*sigma(2)^2)/E+(4*sigma(1)^2)/(9*G)))*[1/G+(4*s igma(2)^2)*aux -4/3*sigma(1)*sigma(2)*aux;-4/3*sigma(1)*sigma(2)*aux 1/E+4/9*(sigma(1)^2)*aux]; %CÁLCULO DE LA DERIVADA DEL VECTOR DEFORMACIÓN% %MÉTODO CONOCIENDO LA DERIVADA EXACTA% %CARGA LINEAL% if esp==0 deps=@(t)[eps.e;0]; %CARGA EN ESPIRAL% else deps=@(t)[eps.e*(cos(t)-t*sin(t)) ; eps.g*(sin(t)+t*cos(t))]; end % DIFERENCIACIÓN NUMÉRICA % % h=1e-10; % deps=@(t)[(eps.e(t+h)-eps.e(t-h))/(2*h);(eps.g(t+h)-eps.g(th))/(2*h)]; dsigma=A*deps(t); %CALCULAMOS LA DERIVADA DEL VECTOR SIGMA A PARTIR DE LA MATRIZ A Y LA DERIVADA DEL VECTOR DEFORMACIÒN Syst2.m function dsigma=syst2(t,sigma,eps) global K na E G %FUNCIÓN DE INTEGRACIÓN DE ECUACIONES DIFERENCIALES SIMPLIFICADA PARA EL %TRAMO2 Q=sqrt(2/3*sigma(1)^2+2*sigma(2)^2); %SE DEFINE EL MÓDULO EFECTIVO DE TENSIONES% aux=((sqrt(3/2))^(1/na+1))/na*((1/K)^(1/na))*(Q^(1/na-3)); %SE DEFINE UNA VARIABLE AUXILIAR= PHI(Q)/Q^2 PARA REDUCIR POSIBLES PROBLEMAS NUMÉRICOS% %SE DEFINE LA MATRIZ DE LAS DERIVADAS PARCIALES% A=1/(1/E*1/G+(aux)*((4*sigma(2)^2)/E+(4*sigma(1)^2)/(9*G)))*[1/G+(4*s igma(2)^2)*aux -4/3*sigma(1)*sigma(2)*aux;-4/3*sigma(1)*sigma(2)*aux 1/E+4/9*(sigma(1)^2)*aux]; 79 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones %CÁLCULO DE LA DERIVADA DEL VECTOR DEFORMACIÓN% %MÉTODO CONOCIENDO LA DERIVADA EXACTA% deps=@(t)[(-eps.e*sin(t));(eps.g*cos(t))]; % DIFERENCIACIÓN NUMÉRICA % % h=1e-10; % deps=@(t)[(eps.e(t+h)-eps.e(t-h))/(2*h);(eps.g(t+h)-eps.g(th))/(2*h)]; dsigma=A*deps(t); %CALCULAMOS LA DERIVADA DEL VECTOR SIGMA A PARTIR DE LA MATRIZ A Y LA DERIVADA DEL VECTOR DEFORMACIÒN %se pueden meter la funcion y derivarla o meter directamente la funcion %derivada % Systrk4d.m function [t,sigma]=systrk4d(tini,tfin,sigma0,eps,h,esp) d=0; %PARÁMETRO DE DESCARGA% sigma(:,1)=sigma0; %PUNTO INICIAL DE TENSIONES% t=tini:h:tfin; %INICIALIZACIÓN DEL VECTOR T A PARTIR DEL PASO DE INTEGRACIÓN H% sigmac=[0;0]; %CENTRO DE CIRCUNFERENCIA DE PLASTIFICACIÓN INICIAL% c=1; %CIRCUNFERENCIA DE PLASTIFICACIÓN ACTIVA% i=1; %INDICE DE PASO PARA CADA INSTANTE DE T% while i<=(length(t)-1) %MIENTRAS QUE NO LLEGUEMOS AL INSTANTE FINAL DE T if d>=1 %SI SE HA DETECTADO DESCARGA% sigmad(:,1)=sigma(:,i-1); %EL PRIMER PUNTO DE DESCARGA ES EL PASO ANTERIOR A LA DESCARGA% q0=2*modu(sigmad(:,1)); %DEFINIMOS EL MÓDULO EFECTIVO DE LA PRIMERA CIRCUNFERENCIA DE PLASTIFICACIÓN %RESOLVEMOS EL SISTEMA CON LAS ECUACIONES DE DESCARGA, MEDIANTE MÉTODO %RUNGE KUTTA PASO A PASO DE 4TH ORDEN K1=systd(t(i),sigma(:,i),sigmad(:,d),sigmac(:,c),q0,eps,esp); K2=systd(t(i)+h/2,(sigma(:,i)+K1/2*h),sigmad(:,d),sigmac(:,c),q0,eps,e sp); K3=systd(t(i)+h/2,(sigma(:,i)+K2/2*h),sigmad(:,d),sigmac(:,c),q0,eps,e sp); K4=systd(t(i)+h,(sigma(:,i)+K3*h),sigmad(:,d),sigmac(:,c),q0,eps,esp); sigma(:,i+1)=sigma(:,i)+h*(K1+2*K2+2*K3+K4)/6; % q0=abs((2/3*(sigma(1,i)-sigmad(1,d))^2+2*(sigma(2,i)- sigmad(2,d))^2)/(2/3*(sigmac(1,c)-sigmad(1,d))*(sigma(1,i)- Anexo. Código de programación necesario para las simulaciones 80 80 sigmad(1,d))+2*(sigmac(2,c)-sigmad(2,d))*(sigma(2,i)- sigmad(2,d)))*q0/2); i=i+1; %INCREMENTAMOS INDICE DE PASO% while modu(sigma(:,i))<modu(sigmad(:,1)) && i<=(length(t)-1) %MIENTRAS QUE LAS TENSIONES NO SEAN MAYORES AL PRIMER PUNTO DE DESCARGA% while modu(sigma(:,i)-sigmac(:,c))<modu(sigmad(:,d)-sigmac(:,c)) && i<=(length(t)-1) %MIENTRAS LA DISTANCIA DEL PUNTO DE TENSIÓN AL CENTRO ACTUAL DE PLASTIFICACIÓN SEA MENOR A LA DISTANCIA DEL PUNTO DE DESCARGA ACTIVO AL CENTRO ACTUAL DE PLSATIFICACIÓN% qaux=abs((2/3*(sigma(1,i-1)-sigmad(1,d))^2+2*(sigma(2,i-1)- sigmad(2,d))^2)/(2/3*(sigmac(1,c)-sigmad(1,d))*(sigma(1,i-1)- sigmad(1,d))+2*(sigmac(2,c)-sigmad(2,d))*(sigma(2,i-1)- sigmad(2,d)))*q0/2); %CALCULAMOS EL DIÁMETRO DE LA NUEVA CIRCUNFERENCIA PROVISIONAL DE PLASTIFICACIÓN, QUE SE USARÁ PARA CALCULAR EL NUEVO CENTRO% %CALCULAMOS EL CENTRO DE LA NUEVA CIRCUNFERENCIA PROVISIONAL DE PLASTIFICACIÓN% c=c+1; sigmac(:,c)=sigmad(:,d)-qaux/2*(sigmad(:,d)-sigmac(:,c- 1))/(modu(sigmad(:,d)-sigmac(:,c-1))); if modu(sigma(:,i)-sigmac(:,c))<modu(sigmad(:,d)-sigmac(:,c)) %SI EL NUEVO PUNTO DE TENSIONES CALCULADO ANTERIORMENTE ESTA DENTRO DE LA CIRCUNFERENCIA PROVISIONAL DE PLASTIFICACIÓN% %SE CREA UN NUEVO PUNTO DE DESCARGA% d=d+1; sigmad(:,d)=sigma(:,i-1); q0=abs(((modu(sigmad(:,d)-sigmad(:,d-1)))^2)/(2/3*(sigmac(1,c-1)- sigmad(1,d-1))*(sigmad(1,d)-sigmad(1,d-1))+2*(sigmac(2,c-1)- sigmad(2,d-1))*(sigmad(2,d)-sigmad(2,d-1)))*q0/2); %SE CALCULA EL DIÁMETRO DE LA NUEVA CIRCUNFERENCIA DE PLASTIFICACIÓN% else %SI EL NUEVO PUNTO DE TENSIONES CALCULADO ANTERIORMENTE NO ESTÁ DENTRO DE LA CIRCUNFERENCIA PROVISIONAL DE PLASTIFICACIÓN c=c-1; %SE BORRA DICHA CIRCUNFERENCIA AUXILIAR% end %CALCULAMOS UN NUEVO PUNTO DE TENSIONES USANDO EL MÉTODO RUNGE KUTTA DE %CUARTO ORDEN K1=systd(t(i),sigma(:,i),sigmad(:,d),sigmac(:,c),q0,eps,esp); K2=systd(t(i)+h/2,(sigma(:,i)+K1/2*h),sigmad(:,d),sigmac(:,c),q0,eps,e sp); K3=systd(t(i)+h/2,(sigma(:,i)+K2/2*h),sigmad(:,d),sigmac(:,c),q0,eps,e sp); 81 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones K4=systd(t(i)+h,(sigma(:,i)+K3*h),sigmad(:,d),sigmac(:,c),q0,eps,esp); sigma(:,i+1)=sigma(:,i)+h*(K1+2*K2+2*K3+K4)/6; % q0=abs(((modu(sigma(:,i)-sigmad(:,d)))^2)/(2/3*(sigmac(1,c)- sigmad(1,d))*(sigma(1,i)-sigmad(1,d))+2*(sigmac(2,c)- sigmad(2,d))*(sigma(2,i)-sigmad(2,d)))*q0/2); i=i+1; %INCREMENTAMOS INDICE DE PASO% while modu(sigma(:,i)-sigmac(:,c))>modu(sigmad(:,d)-sigmac(:,c)) && d>1 && c>1 %MIENTRAS QUE EL NUEVO PUNTO CALCULADO SALGA DE LA CIRCUNFERENCIA DE %PLASTIFICACIÓN ACTIVA% d=d-1; %BORRAMOS EL PUNTO DE DESCARGA ACTIVO Y LA CIRCUNFERENCIA DE PLASTIFICACIÓN ACTIVA, VOLVIENDO A LA ANTERIOR c=c-1; %EL PUNTO DE DESCARGA end if d>1 %SI EL PUNTO DE DESCARGA ACTIVO NO ES EL PRIMERO% %RECALCULAMOS LA CIRCUNFERENCIA DE PLASTIFICACIÓN QUE PASA A SER %LA ACTUAL q0=abs(((modu(sigmad(:,d)-sigmad(:,d- 1)))^2)/(2/3*(sigmac(1,c-1)-sigmad(1,d-1))*(sigmad(1,d)-sigmad(1,d- 1))+2*(sigmac(2,c-1)-sigmad(2,d-1))*(sigmad(2,d)-sigmad(2,d- 1)))*q0/2); else %SI EL PUNTO DE DESCARGA ACTIVO ES EL PRIMERO% %RECALCULAMOS LA CIRCUNFERENCIA DE PLASTIFICACIÓN QUE PASA A SER %LA ACTUAL q0=2*modu(sigmad(:,1)); end end end end %SI SALIMOS DEL BUCLE, ES POR QUE LAS TENSIONES HAN SUPERADO EL PRIMER %PUNTO DE DESCARGA% d=0; %PONEMOS EL PARÁMETRO DE DESCARGA A 0% if i<=(length(t)-1) %RESOLVEMOS EL SISTEMA CON LAS ECUACIONES DE DESCARGA, MEDIANTE MÉTODO %RUNGE KUTTA PASO A PASO DE 4TH ORDEN K1=syst(t(i),sigma(:,i),eps,esp); K2=syst(t(i)+h/2,(sigma(:,i)+K1/2*h),eps,esp); K3=syst(t(i)+h/2,(sigma(:,i)+K2/2*h),eps,esp); K4=syst(t(i)+h,(sigma(:,i)+K3*h),eps,esp); sigma(:,i+1)=sigma(:,i)+h*(K1+2*K2+2*K3+K4)/6; Anexo. Código de programación necesario para las simulaciones 82 82 if modu(sigma(:,i+1))<modu(sigma(:,i)) %SI LAS TENSIONES DEL NUEVO PUNTO HAN DISMINUIDO RESPECTO AL ANTERIOR% d=1; %SE CREA EL PRIMER PUNTO DE DESCARGA% end i=i+1; %INCREMENTAMOS ÍNDICE DE PASO% end end Syst2rk4d.m function [t,sigma]=syst2rk4d(tini,tfin,sigma0,eps,h) d=0; %PARÁMETRO DE DESCARGA% sigma(:,1)=sigma0; %PUNTO INICIAL DE TENSIONES% t=tini:h:tfin; %INICIALIZACIÓN DEL VECTOR T A PARTIR DEL PASO DE INTEGRACIÓN H% sigmac=[0;0]; %CENTRO DE CIRCUNFERENCIA DE PLASTIFICACIÓN INICIAL% c=1; %CIRCUNFERENCIA DE PLASTIFICACIÓN ACTIVA% i=1; %INDICE DE PASO PARA CADA INSTANTE DE T% while i<=(length(t)-1) %MIENTRAS QUE NO LLEGUEMOS AL INSTANTE FINAL DE T if d>=1 %SI SE HA DETECTADO DESCARGA% sigmad(:,1)=sigma(:,i-1); %EL PRIMER PUNTO DE DESCARGA ES EL PASO ANTERIOR A LA DESCARGA% q0=2*modu(sigmad(:,1)); %DEFINIMOS EL MÓDULO EFECTIVO DE LA PRIMERA CIRCUNFERENCIA DE PLASTIFICACIÓN %RESOLVEMOS EL SISTEMA CON LAS ECUACIONES DE DESCARGA, MEDIANTE MÉTODO %RUNGE KUTTA PASO A PASO DE 4TH ORDEN K1=syst2d(t(i),sigma(:,i),sigmad(:,d),sigmac(:,c),q0,eps); K2=syst2d(t(i)+h/2,(sigma(:,i)+K1/2*h),sigmad(:,d),sigmac(:,c),q0,eps) ; K3=syst2d(t(i)+h/2,(sigma(:,i)+K2/2*h),sigmad(:,d),sigmac(:,c),q0,eps) ; K4=syst2d(t(i)+h,(sigma(:,i)+K3*h),sigmad(:,d),sigmac(:,c),q0,eps); sigma(:,i+1)=sigma(:,i)+h*(K1+2*K2+2*K3+K4)/6; % q0=abs((2/3*(sigma(1,i)-sigmad(1,d))^2+2*(sigma(2,i)- sigmad(2,d))^2)/(2/3*(sigmac(1,c)-sigmad(1,d))*(sigma(1,i)- sigmad(1,d))+2*(sigmac(2,c)-sigmad(2,d))*(sigma(2,i)- sigmad(2,d)))*q0/2); i=i+1; %INCREMENTAMOS INDICE DE PASO% while modu(sigma(:,i))<modu(sigmad(:,1)) && i<=(length(t)-1) %MIENTRAS QUE LAS TENSIONES NO SEAN MAYORES AL PRIMER PUNTO DE DESCARGA% 83 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones while modu(sigma(:,i)-sigmac(:,c))<modu(sigmad(:,d)-sigmac(:,c)) && i<=(length(t)-1) %MIENTRAS LA DISTANCIA DEL PUNTO DE TENSIÓN AL CENTRO ACTUAL DE PLASTIFICACIÓN SEA MENOR A LA DISTANCIA DEL PUNTO DE DESCARGA ACTIVO AL CENTRO ACTUAL DE PLSATIFICACIÓN% qaux=abs((2/3*(sigma(1,i-1)-sigmad(1,d))^2+2*(sigma(2,i-1)- sigmad(2,d))^2)/(2/3*(sigmac(1,c)-sigmad(1,d))*(sigma(1,i-1)- sigmad(1,d))+2*(sigmac(2,c)-sigmad(2,d))*(sigma(2,i-1)- sigmad(2,d)))*q0/2); %CALCULAMOS EL DIÁMETRO DE LA NUEVA CIRCUNFERENCIA PROVISIONAL DE PLASTIFICACIÓN, QUE SE USARÁ PARA CALCULAR EL NUEVO CENTRO% %CALCULAMOS EL CENTRO DE LA NUEVA CIRCUNFERENCIA PROVISIONAL DE PLASTIFICACIÓN% c=c+1; sigmac(:,c)=sigmad(:,d)-qaux/2*(sigmad(:,d)-sigmac(:,c- 1))/(modu(sigmad(:,d)-sigmac(:,c-1))); if modu(sigma(:,i)-sigmac(:,c))<modu(sigmad(:,d)-sigmac(:,c)) %SI EL NUEVO PUNTO DE TENSIONES CALCULADO ANTERIORMENTE ESTA DENTRO DE LA CIRCUNFERENCIA PROVISIONAL DE PLASTIFICACIÓN% %SE CREA UN NUEVO PUNTO DE DESCARGA% d=d+1; sigmad(:,d)=sigma(:,i-1); q0=abs(((modu(sigmad(:,d)-sigmad(:,d-1)))^2)/(2/3*(sigmac(1,c-1)- sigmad(1,d-1))*(sigmad(1,d)-sigmad(1,d-1))+2*(sigmac(2,c-1)- sigmad(2,d-1))*(sigmad(2,d)-sigmad(2,d-1)))*q0/2); %SE CALCULA EL DIÁMETRO DE LA NUEVA CIRCUNFERENCIA DE PLASTIFICACIÓN% else %SI EL NUEVO PUNTO DE TENSIONES CALCULADO ANTERIORMENTE NO ESTÁ DENTRO DE LA CIRCUNFERENCIA PROVISIONAL DE PLASTIFICACIÓN c=c-1; %SE BORRA DICHA CIRCUNFERENCIA AUXILIAR% end %CALCULAMOS UN NUEVO PUNTO DE TENSIONES USANDO EL MÉTODO RUNGE KUTTA DE %CUARTO ORDEN K1=syst2d(t(i),sigma(:,i),sigmad(:,d),sigmac(:,c),q0,eps); K2=syst2d(t(i)+h/2,(sigma(:,i)+K1/2*h),sigmad(:,d),sigmac(:,c),q0,eps) ; K3=syst2d(t(i)+h/2,(sigma(:,i)+K2/2*h),sigmad(:,d),sigmac(:,c),q0,eps) ; K4=syst2d(t(i)+h,(sigma(:,i)+K3*h),sigmad(:,d),sigmac(:,c),q0,eps); sigma(:,i+1)=sigma(:,i)+h*(K1+2*K2+2*K3+K4)/6; % q0=abs(((modu(sigma(:,i)-sigmad(:,d)))^2)/(2/3*(sigmac(1,c)- sigmad(1,d))*(sigma(1,i)-sigmad(1,d))+2*(sigmac(2,c)- sigmad(2,d))*(sigma(2,i)-sigmad(2,d)))*q0/2); i=i+1; %INCREMENTAMOS INDICE DE PASO% Anexo. Código de programación necesario para las simulaciones 84 84 while modu(sigma(:,i)-sigmac(:,c))>modu(sigmad(:,d)-sigmac(:,c)) && d>1 && c>1 %MIENTRAS QUE EL NUEVO PUNTO CALCULADO SALGA DE LA CIRCUNFERENCIA DE %PLASTIFICACIÓN ACTIVA% d=d-1; %BORRAMOS EL PUNTO DE DESCARGA ACTIVO Y LA CIRCUNFERENCIA DE PLASTIFICACIÓN ACTIVA, VOLVIENDO A LA ANTERIOR c=c-1; %EL PUNTO DE DESCARGA end if d>1 %SI EL PUNTO DE DESCARGA ACTIVO NO ES EL PRIMERO% %RECALCULAMOS LA CIRCUNFERENCIA DE PLASTIFICACIÓN QUE PASA A SER %LA ACTUAL q0=abs(((modu(sigmad(:,d)-sigmad(:,d- 1)))^2)/(2/3*(sigmac(1,c-1)-sigmad(1,d-1))*(sigmad(1,d)-sigmad(1,d- 1))+2*(sigmac(2,c-1)-sigmad(2,d-1))*(sigmad(2,d)-sigmad(2,d- 1)))*q0/2); else %SI EL PUNTO DE DESCARGA ACTIVO ES EL PRIMERO% %RECALCULAMOS LA CIRCUNFERENCIA DE PLASTIFICACIÓN QUE PASA A SER %LA ACTUAL q0=2*modu(sigmad(:,1)); end end end end %SI SALIMOS DEL BUCLE, ES POR QUE LAS TENSIONES HAN SUPERADO EL PRIMER %PUNTO DE DESCARGA% d=0; %PONEMOS EL PARÁMETRO DE DESCARGA A 0% if i<=(length(t)-1) %RESOLVEMOS EL SISTEMA CON LAS ECUACIONES DE DESCARGA, MEDIANTE MÉTODO %RUNGE KUTTA PASO A PASO DE 4TH ORDEN K1=syst2(t(i),sigma(:,i),eps); K2=syst2(t(i)+h/2,(sigma(:,i)+K1/2*h),eps); K3=syst2(t(i)+h/2,(sigma(:,i)+K2/2*h),eps); K4=syst2(t(i)+h,(sigma(:,i)+K3*h),eps); sigma(:,i+1)=sigma(:,i)+h*(K1+2*K2+2*K3+K4)/6; if modu(sigma(:,i+1))<modu(sigma(:,i)) %SI LAS TENSIONES DEL NUEVO PUNTO HAN DISMINUIDO RESPECTO AL ANTERIOR% d=1; %SE CREA EL PRIMER PUNTO DE DESCARGA% end i=i+1; %INCREMENTAMOS ÍNDICE DE PASO% end end 85 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones Systd.m function dsigma=systd(t,sigma,sigmad,sigmac,q0,eps,esp) global K na E G q=abs(modu(sigma-sigmad))^2/(2/3*(sigmac(1)-sigmad(1))*(sigma(1)- sigmad(1))+2*(sigmac(2)-sigmad(2))*(sigma(2)-sigmad(2)))/2; %CÁLCULO DEL NUEVO DIÁMETRO, DIVIDIA ENTRE EL DIAMETRO ANTERIOR (PARA EVITAR PROBLEMAS NUMERICOS SI Q0 %ES CERCANO A CERO% n=[4/(3*q)*(sigma(1)-sigmad(1))-4/3*(sigmac(1)- sigmad(1));4/q*(sigma(2)-sigmad(2))-4*(sigmac(2)-sigmad(2))]; %DIRECCIÓN DE ENDURECIMIENTO MULTIPLICADO POR Q0 (PARA EVITAR PROBLEMAS %NUMÉRICOS% aux=((sqrt(3/2))^(1/na+1))/na*((1/K)^(1/na))*((q/2)^(1/na- 1))*(q0)^(1/na-3); %FLUJO PLÁSTICO DE LA NUEVA CIRCUNFERENCIA DE PLASTIFICACIÓN, DIVIDO ENTRE %Q0^2 B=[(aux*n(1)^2+1/E) (aux*n(1)*n(2));(aux*n(1)*n(2)) (aux*n(2)^2+1/G)]; %MONTAJE DE LA MATRIZ B, EL 1/Q0^2 QUE MULTIPLICA A AUX SE COMPENSA CON EL Q0 QUE %MULTIPLICA EN CADA COMPONENTE DEL VECTOR N %CÁLCULO DE LA DERIVADA DEL VECTOR DEFORMACIÓN% %DIFERENCIACIÓN SIMBÓLICA% %LINEAL% if esp==0 deps=@(t)[eps.e;0]; %ESPIRAL% else deps=@(t)[eps.e*(cos(t)-t*sin(t)) ; eps.g*(sin(t)+t*cos(t))]; end % DIFERENCIACIÓN NUMÉRICA % % h=1e-10; % deps=@(t)[(eps.e(t+h)-eps.e(t-h))/(2*h);(eps.g(t+h)-eps.g(th))/(2*h)]; dsigma=B\deps(t); Syst2d.m function dsigma=syst2d(t,sigma,sigmad,sigmac,q0,eps) global K na E G q=abs((modu(sigma-sigmad))^2/(2/3*(sigmac(1)-sigmad(1))*(sigma(1)- sigmad(1))+2*(sigmac(2)-sigmad(2))*(sigma(2)-sigmad(2))))/2; %CÁLCULO DEL NUEVO DIÁMETRO, DIVIDIA ENTRE EL DIAMETRO ANTERIOR (PARA EVITAR PROBLEMAS NUMERICOS SI Q0 %ES CERCANO A CERO% n=[4/(3*q)*(sigma(1)-sigmad(1))-4/3*(sigmac(1)- sigmad(1));4/q*(sigma(2)-sigmad(2))-4*(sigmac(2)-sigmad(2))]; Anexo. Código de programación necesario para las simulaciones 92 92 [tspan2,sigma2]=ode45(@systcl2,tspan2,[sigma(end,1);sigma(end,2)],opti ons,epsc2); %INTEGRACIÓN DE LAS ECUACIONES DE FLUJO PLÁSTICO EN DERIVADAS PARCIALES PARA EL SEGUNDO TRAMO% epsilon=[eps.e(tspan) ;eps2.e(tspan2)]; %SE CALCULAN LAS DEFORMACIONES EN COMPONENTES EPSILON Y GAMMA PARA CADA INSTANTE DEL VECTOR T% gamma=[eps.g(tspan) ;eps2.g(tspan2)]; %SE AGRUPAN A SU VEZ LOS TRAMOS 1 Y 2% sigma=[sigma; sigma2]; %SE AGRUPAN LAS TENSIONES DE LOS TRAMOS 1 Y 2 , POR COMPONENTES SIGMA Y TAU% tspan=[tspan ;tspan2]; %SE AGRUPA EL VECTOR DEL PARÁMETRO T DE LOS TRAMOS 1 Y 2% Q=modu(sigma'); %SE CALCULA EL MÓDULO EFECTIVO DE LAS TENSIONES PARA CADA PUNTO CALCULADO% dm=diff(Q); %SE CALCULA EL VECTOR DE DIFERENCIAS DEL MÓDULO DE SIGMA, PARA VER SI HA DECRECIDO EN ALGÚN PUNTO% N=tspan/(2*pi); %DESHACEMOS EL CAMBIO DE VARIABLE DE T A NÚMERO DE CICLOS TOTALES, PARA CADA INSTANTE EVALUADO% [m,i]=max(Q); %SE CALCULA EL MÁXIMO DEL MÓDULO EFECTIVO DE TENSIONES, Y EL ÍNDICE EN EL CUAL APARECE% %ESTE VALOR SERÁ MÁXIMO JUSTO ANTES DE LA ESTABILIZACIÓN, LO TOMAREMOS COMO ORDEN DE MAGNITUD% i=tspan(i)/(2*pi); %CAMBIAMOS LA VARIABLE ÍNDICE A NÚMERO DE CICLOS tabeps(pos)=emax; %GUARDAMOS LOS VALORES DE EMAX, NEST, Y QMAX PARA CADA SIMULACIÓN% Nest(pos)=i; Qmax(pos)=m; pos=pos+1; end toc figure (1) plot(tabeps,Qmax) %REPRESENTAMOS LAS TENDENCIAS E-QMAX Y E-NEST title('EVOLUCIÓN DEL NUMERO DE CICLOS EN EL CUAL SE ESTABILIZA Q') xlabel('MÓDULO DE EPSILON') ylabel('NÚMERO DE CICLOS NECESARIO') figure (2) plot(tabeps,Nest) xlabel('MÓDULO DE EPSILON') ylabel('MÓDULO DE SIGMA Q') title('EVOLUCIÓN DEL VALOR DE Q EN LA ESTABILIZACIÓN') save('tabl.txt','tabeps','Nest','Qmax','-ascii') %ALMACENAMOS LOS DATOS, PUES EL TIEMPO DE EJECUCION DE ESTE %SCRIPT ES ALTO 93 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones Ajuste.m %SCRIPT PARA AJUSTAR LAS TENDENCIAS EMAX-QMAX Y EMAX-NEST DE FORMAS LINEAL %Y EXPONENCIAL RESPECTIVAMENTE %SE DEBEN EJECUTAR DESPUES DE ORDENADATOS.M, O TENER LOS FICHEROS QUE ESTE %GENERA clear,clc close all A=load('tabl.txt'); %CARGA DE LOS DATOS ALMACENADOS POR EL SCRIPT ORDENADATOS.M tabeps=A(1,:); %ORGANIZACIÓN DE LOS DATOS PARA CADA VARIABLE% Nest=A(2,:); Qmax=A(3,:); fact=100; %FACTOR DE MULTIPLICACIÓN PARA MEJOR LA CONVERGENCIA% tabepsch=tabeps*fact; %MULTIPLICAMOS EMAX POR UN FACTOR DE 100 PARA FACILITAR %LA CONVERGENCIA DEL AJUSTE EXPONENCIAL fun=@(coef,x)coef(1)*(coef(2)).^(x); %DEFINICION DE UNA FUNCIÓN DE AJUSTE TIPO EXPONENCIAL% coef0=[10;10]; %VALORES INICIALES DE LOS COEFICIENTES% coef = lsqcurvefit(fun,coef0,tabepsch,Nest); %CÁLCULO DE LOS COEFICIENTES PARA NUESTROS DATOS (EMAX, NEST)% Nest2=fun(coef,tabepsch); %CÁLCULO DEL NUMERO DE CICLOS NEST SEGÚN EL AJUSTE% fun2=@(coef2,x)coef2(1).*x+coef2(2); %DEFINICIÓN DE UNA FUNCIÓN DE AJUSTE TIPO LINEAL% coef20=[1,1]; %VALORES INICIALES DE LOS COEFICIENTES% coef2 = lsqcurvefit(fun2,coef20,tabeps,Qmax); Qmax2=fun2(coef2,tabeps); %COMPARACIÓN GRÁFICA ENTRE EL AJUSTE REALIZADO, Y LOS VALORES NÚMERICOS %OBTENIDOS POR INTEGRACIÓN% %GRÁFICA EMAX-NEST% figure (1) plot(tabepsch/fact,Nest2,'b',tabeps,Nest,'r') title('Evolución del número de ciclos frente a deformaciones impuestas') ylabel('\bf N_{c} [Ciclos]') Anexo. Código de programación necesario para las simulaciones 94 94 xlabel('\bf Norma efectiva de deformaciones') legend('Ajuste Exponencial','Datos Numéricos','location','north') error1=max(abs(Nest2-Nest)/Nest)*100; %ERROR MÁXIMO COMETIDO CON LA APROXIMACIÓN EN TANTO POR CIENTO% er=num2str(error1); ab=[er '%']; disp(ab) %GRÁFICA EMAX-QMAX% figure (2) plot(tabeps,Qmax2,'b',tabeps,Qmax,'r') ylabel('\bf Q [MPa]') xlabel('\bf Norma efectiva de deformaciones') title('Evolución del valor de Q según deformaciones impuestas') legend('Ajuste Lineal','Datos numéricos','location','north') error2=max(abs(Qmax2-Qmax)/Qmax)*100; %ERROR MÁXIMO COMETIDO CON LA APROXIMACIÓN EN TANTO POR CIENTO% Ajuste2.m %SCRIPT QUE SE EJECUTA DESPUES DE ESPIRAL.M% %PERMITE CALCULAR UN AJUSTE LOGARITMICO EN LA EVOLUCIÓN N-Q PARA %DEFORMACIONES ALTAS EN LAS QUE EL VECTOR Q NO SE ESTABILIZA DE FORMA %TEMPRANA tc=round(tcarga/h)+1/h; N1=N(tc:end,:)'; %INICIALIZMOS LOS VECTORES N Y Q A PARTIR DE CUANDO TERMINA LA CARGA, CON UN CIERTO MARGEN Q1=Q(:,tc:end); fun=@(coef,N1)coef(1)*log(N1-coef(2))+coef(3); %DEFINICION DE UNA FUNCIÓN DE AJUSTE TIPO LOGARÍTMICA% coef0=[1,1,10]; %VALORES INICIALES DE LOS COEFICIENTES% coef = lsqcurvefit(fun,coef0,N1,Q1); %CÁLCULO DE LOS COEFICIENTES PARA NUESTROS DATOS (N1,Q1)% Q1=fun(coef,N1); Q1=real(Q1); %CÁLCULO DE LA GRAFICA N-Q SEGÚN EL AJUSTE% Q2=Q(tc:end); error=Q1-Q2; error=error./Q2; error=abs(error); [error,i]=max(error); %CÁLCULO DEL ERROR MAXIMO QUE SE COMETE CON EL AJUSTE% %COMPARACIÓN GRÁFICA ENTRE EL AJUSTE REALIZADO, Y LOS VALORES NÚMERICOS %OBTENIDOS POR INTEGRACIÓN% figure (9) 95 Simulación numérica de endurecimiento en ensayos biaxiales cíclicos no-proporcionales mediante una teoría de plasticidad basada en la idea de distancia en el espacio de tensiones plot(N,Q,'r',N1,Q1,'b') title('\bf Evolución del módulo efectivo de tensiones según el número de ciclos') xlabel('\bf N_{c} [Ciclos]') ylabel('\bf Norma efectiva de tensiones [MPa]') legend('Datos Numéricos','Ajuste Logarítmico','location','best')