Repositorio Institucional de Documentos

1 Trabajo Fin de Máster Cálculo no lineal de estructuras planas de barras. Obtención de las curvas de equilibrio en el análisis de postpandeo por el método de longitud del arco de Riks. Autor José Manuel Lou Pina Director Miguel Ángel Martínez Barca Escuela de Ingeniería y Arquitectura 2.014 2 RESUMEN. Este trabajo pretende ser una herramienta de apoyo en el ejercicio de mi actividad profesional. Es de todos conocido el hecho de que las estructuras deben diseñarse cumpliendo ciertos requisitos de bienestar los cuales quedan cubiertos manteniendo las deformaciones de éstas por debajo de ciertos valores normativos. Así, en el diseño de estructuras, raramente son las limitaciones tensionales las que marcan el dimensionamiento sino más bien son las flechas o deformaciones las que lo controlan. Con todo, en este trabajo se pretende estudiar el postpandeo de una tipología particular de estructuras, obteniendo las trayectorias de equilibrio de los grados de libertad de interés, de forma que se puedan visualizar los puntos límite de dichas trayectorias y obtener así la respuesta completa de las mismas. Desde un punto de vista práctico, este trabajo se pretende aplicar fundamentalmente para la obtención de las cargas criticas de pandeo de pórticos rígidos planos, a partir de las cuales se puedan obtener las longitudes de pandeo correctas de los pilares de la estructura además de en algunos casos particulares, poder seguir su respuesta de las mismas más allá del pandeo. Esto último, aunque en ocasiones puede no tener mucho sentido práctico puede servir para, en algún caso, discutir de forma razonada el comportamiento real de esta tipología estructural. Por otro lado, el hecho de abordar este problema me servirá de base técnica para, en un futuro, intentar atacar otros problemas no menos interesantes, como puede ser el tema de la respuesta de estructuras con no linealidad del material, de forma y manera que pueda obtener soluciones razonadas que pueda defender y que, por supuesto, como en este caso, garanticen la seguridad estructural. Las estructuras analizadas se resuelven mediante métodos matriciales bajo una formulación corrotacional, empleando elementos barra de Euler Bernoulli y el algoritmo empleado para la solución del problema de postpandeo es el método de longitud de arco de Riks. 3 INDICE RESUMEN. 2 1. MOTIVACIÓN. 4 2. INTRODUCCIÓN. 5 3. EXPOSICIÓN DEL TRABAJO. 6 4. ALCANCE DEL CODIGO PROGRAMADO. 6 5. TEORIA BÁSICA. 8 5.1. Descripción del movimiento. 8 5.2. Variantes de la formulación Lagrangiana. 9 5.2.1. Formulación Lagrangiana Total. 9 5.2.2. Formulación Lagrangiana Actualizada. 9 5.2.3. Formulación Lagrangiana Corrotacional. 9 5.3. Métodos de solución. 10 5.3.1. Método de Newton-Raphson. 10 5.3.2. Puntos críticos. 15 5.3.3. Métodos de longitud del arco. 16 5.4. Determinación del tamaño del salto de carga. 20 5.5. Ajuste automático de la longitud del salto de carga. 21 5.6. Cálculo de puntos de inestabilidad. 21 5.7. Criterios de convergencia. 22 5.7.1. Criterio de convergencia de desplazamientos 22 5.7.2. Criterio de convergencia de fuerzas 22 5.7.3. Criterio de convergencia de la energía 23 5.8. Matriz de rigidez tangente de la barra. 23 5.8.1. Formulación de la matriz de rigidez tangente de la barra. 23 6. PROBLEMAS RESUELTOS CON EL CODIGO. 30 6.1. Barra inclinada (caso 1). 30 6.2. Barra inclinada (caso 2). 33 6.3. Pórtico simple. 36 6.4. Arco rebajado con carga excéntrica. 38 6.5. Pilar simple. 43 6.6. Pórtico plano. 47 7. FUNCIONAMIENTO DEL PROGRAMA. 51 8. BIBLIOGRAFÍA. 62 4 1. MOTIVACIÓN. Como calculista de estructuras en activo y dada la responsabilidad que implican mis decisiones y consejos, me siento en la obligación de mantener una formación continuada, cuestión esta que, aun siendo siempre gratificante, no resulta en ocasiones fácil. Comparando desde un punto de vista únicamente técnico el análisis lineal y no lineal de estructuras, así como en el análisis lineal de estructuras, tanto si este se lleva a cabo con el empleo de programas propios como si se hace con programas comerciales, "no se requiere" de un control exhaustivo por parte del ingeniero puesto que las técnicas empleadas para la obtención de la solución son técnicas matriciales por todos bien conocidas y fácilmente chequeables con las que se obtiene un único resultado final, esto es sustancialmente diferente cuando lo que se persigue es la obtención de la solución de un problema no lineal. En este último caso, conceptos como matriz de rigidez tangente, tolerancia, control en carga o en desplazamientos, convergencia de la solución, residuo, etc... hacen que el técnico deba gobernar el análisis fijando valores de ciertos parámetros de los cuales, para hacer un uso cabal de los programas empleados, debe de conocer su correcto significado y el modo en el que los valores fijados para estos parámetros, afectan al resultado final obtenido. La necesidad del análisis no lineal de estructuras, exigido en ciertos casos incluso a nivel normativo, me obliga por tanto a conocer la teoría en la que se fundamentan las técnicas de solución de este tipo de problemas. Son muchos los programas comerciales existentes en el mercado que resuelven el problema no lineal y de los que habitualmente hacemos uso los ingenieros. Estos programas emplean algoritmos de solución los cuales resultan ser en algunos casos auténticas cajas negras no siendo posible incluso llegar a manipular. Por otra parte, estos programas comerciales lógicamente no resuelven todo tipo de problemas y se hace necesario conocer cuáles son las limitaciones de los mismos, no siendo esto último posible sin conocer previamente las técnicas que emplean y a través de las cuales, nos proporcionan la solución. Por otro lado, el conocimiento de la teoría en la que se fundamentan estos algoritmos nos puede servir de base para la ejecución de programas propios con los cuales obtener soluciones aceptables y completamente rigurosas, sin necesidad de adquirir programas comerciales los cuales siempre llevan asociados costes importantes. Con el presente trabajo se pretende encontrar la solución para el problema del postpandeo en estructuras planas de barras de la tipología expuesta en el punto 3 y que los programas comerciales de los que dispongo no me dan la solución. Aunque el código que aparece en el trabajo es bastante elemental, y sin el ánimo de inventar nada, las soluciones que se obtienen permiten ver claramente las cargas de pandeo, así como las trayectorias de equilibrio de los puntos de interés incluso después del pandeo de la estructura. El hecho de obtener este tipo de soluciones a través de un código propio, además de gratificante, resulta del mayor interés didáctico al tener la necesidad de conocer más o menos en profundidad la formulación del problema. 5 2. INTRODUCCIÓN. En el análisis lineal de estructuras, la formulación de las ecuaciones de equilibrio y compatibilidad se lleva a cabo en la geometría indeformada de la estructura admitiendo en todo momento que los desplazamientos son pequeños y que el material tiene un comportamiento lineal. Estas hipótesis dan lugar a un sistema de ecuaciones lineales en el que es válido el principio de superposición. Así, dadas una acciones exteriores (𝒇= ∑𝒇𝑖) , mayoradas por unos coeficientes (𝛾𝑖), las solicitaciones 𝒔= 𝒔(∑𝒇𝑖) y los movimientos 𝒖= 𝒖(∑𝒇𝑖) las podemos poner como 𝒔= 𝒔(∑𝛾𝑖∗𝒇𝑖)=∑𝛾𝑖∗𝒔(𝒇𝑖) (Ec1.) 𝒖= 𝒖(∑𝛾𝑖∗𝒇𝑖)=∑𝛾𝑖∗𝒖(𝒇𝑖) (Ec2.) Como es sabido, este tipo de análisis no captura aquellos desplazamientos que son consecuencia de la propia deformación de la estructura (efecto P- ∆ ) y que en estructuras esbeltas con valores altos de axiles pueden resultar significativos, lo que conlleva a un cálculo de esfuerzos en contra de la seguridad. En el análisis de estructuras con no linealidad geométrica, se debe considerar la modificación de la geometría inicial de la estructura (debida a su deformación), al formular las ecuaciones de equilibrio y compatibilidad. Sin embargo, para las estructuras comunes de edificación, caben simplificaciones en ese proceso. Aunque se abandona la hipótesis de pequeños desplazamientos, se pueden suponer giros pequeños (𝜃≪1) lo cual redunda en una simplificación muy notable del cálculo matemático. Se habla entonces de teoría de segundo orden. En este caso la expresión que se toma para la curvatura es la expresión aproximada dada por 𝝌= 1/𝑟=𝑦′′=𝑑𝜃 𝑑𝑥=𝑑2𝑦 𝑑𝑥2 (Ec3.) Cuando en el análisis de la estructura se toma para la curvatura su valor exacto (Ec4.), se habla de teoría de tercer orden. 𝝌= 1/𝑟=𝑦′′ (1+𝑦′2)3/2 (Ec4.) El análisis en teoría de segundo orden, requerido incluso a nivel normativo, resulta suficiente en la mayoría de los casos prácticos. No obstante, dado que con este análisis no se tienen en cuenta los efectos por grandes deformaciones, acota la aplicación del método a aquellas estructuras en las que los movimientos y rotaciones no sean elevados. Como el análisis P-∆ es no lineal, no podemos aplicar el principio de superposición y obtener el resultado para una combinación de carga como la suma de resultados de los estados de carga elementales. Para estimar los esfuerzos y desplazamientos se ha de realizar un estudio sobre la estructura en posición deformada (teoría de segundo orden), lo que requiere de procedimientos iterativos de cálculo basados en la teoría de primer orden, o métodos matriciales que incorporan la influencia que tienen en los desplazamientos del sistema los esfuerzos axiles de las barras y los cambios de posición de los nudos. La aplicación de la teoría de segundo orden al análisis de pórticos con imperfecciones evita las incertidumbres que se presentan en el diseño basado en el análisis lineal de esfuerzos. Por ello se recomienda, con algunas excepciones, aplicar el análisis de segundo orden a los pórticos translacionales. ****************************************************************************************************************************** 6 3. EXPOSICIÓN DEL TRABAJO. Se trata de programar en Mathematica las rutinas que nos van a servir para analizar estructuras con no linealidad geométrica, dentro de un tipo particular de tipología estructural como es la del pórtico plano. En el código presentado solo se formula la solución del problema no lineal debido a la geometría y para esta tipología de estructura, suponiendo en todo momento un comportamiento elástico lineal del material (pequeñas deformaciones) y sin entrar en los casos en los que la no linealidad del problema se debe a las condiciones de contorno (como es el caso de los problemas de contacto o de fuerzas seguidoras). El problema general en el que intervienen los tres tipos de no linealidades (por geometría, por material y por condiciones de contorno), presenta una gran dificultad y carga de programación, no siendo el objetivo de este trabajo. El trabajo se plantea como una puesta en práctica de lo visto en el Máster de Mecánica Aplicada, en particular, en la asignatura de Teoría de barras, placas y láminas. Solo para los casos en los que la no linealidad del problema se deba a condiciones geométricas, se pretende obtener las trayectorias o curvas de equilibrio de los nudos de interés de las estructuras que se analicen, y obtener los puntos límite de estas trayectorias (puntos de bifurcación del equilibrio, puntos críticos y puntos de retroceso), para lo cual se aplica el método de Riks de longitud del arco, el cual se comenta en los siguientes apartados. ***************************************************************************************************************************** 4. ALCANCE DEL CODIGO PROGRAMADO. Este trabajo trata únicamente de estructuras planas de barras del tipo pórtico de nudos rígidos. Los apoyos de las estructuras analizadas son los únicos nudos que admiten articulación. Para la resolución del sistema de ecuaciones que surge al plantear el equilibrio de la estructura se utiliza el cálculo matricial. Las diferentes barras de las estructuras se discretizan en un número adecuado de elementos barra formulados como elemento viga de Euler-Bernoulli, elemento que admite esfuerzos de tracción-compresión y flexión. Este elemento presenta dos nodos definidos por las coordenadas de sus extremos las cuales se deben de introducir en la etapa de preproceso. Además de las coordenadas de los nodos de los elementos se deberán introducir las secciones de éstos y sus inercias. Con el elemento programado, no es posible obtener la respuesta de la estructura fuera del plano de la misma lo que limita el uso de este código a estructuras planas con la tipología de pórtico rígido plano. Por el tipo de formulación empleada, la programación del método permite el estudio de estructuras planas y esbeltas en las que pueden tener lugar grandes desplazamientos, resolviendo por lo tanto el problema de la no linealidad geométrica. No obstante, para limitar la complejidad de la formulación, se ha admitido en todo momento que se está en pequeñas deformaciones, no siendo posible con este código resolver problemas en los que esta última hipótesis no sea válida. Por otro lado, dado que se admite en todo momento que el material trabaja en régimen elástico lineal, aquellos problemas en los que se presentan no linealidades debidas al material, no pueden ser abordados con este código. Se admite en todo momento que las condiciones de apoyo no se modifican a lo largo del proceso de carga, por lo que aquellos problemas en los que se presentan no linealidades debidas a las condiciones de contorno, no pueden ser abordados con este código. Dada la velocidad de los ordenadores actuales y la particular tipología estructural que se estudia en este trabajo, para simplificar el código, las barras solo admiten cargas externas (fuerzas o fuerzas y momentos) 7 aplicadas en los nudos de la malla. Para aquellas estructuras, en las existan barras que presenten cargas distribuidas en parte o toda su longitud, se podrán subdividir estas en un número suficiente de elementos, al objeto de poder capturar el efecto real de la carga distribuida sobre la barra y por tanto, sobre la estructura. Asimismo, dado que no tenemos limitación en el número de elementos, podríamos con este código y tras una división adecuada de las barras, analizar estructuras en las que se presentasen barras con inercia variable. ***************************************************************************************************************************** 8 5. TEORIA BÁSICA. 5.1. Descripción del movimiento. En el análisis de estructuras geométricamente no lineales, es fundamental la elección de un sistema de descripción matemática del movimiento y del estado de tensión y deformación de los puntos de la estructura en las distintas configuraciones deformadas por las que pasa en el proceso de carga. En teoría lineal o de deformaciones infinitesimales, no se establece distinción entre las distintas configuraciones de la estructura, ya que se supone que la geometría y características mecánicas de la estructura son invariables. Lo que caracteriza al análisis geométricamente no lineal es, precisamente, la distinción entre las distintas configuraciones a la hora de establecer las ecuaciones que describen el comportamiento mecánico de la estructura. Desde el punto de vista de la mecánica de sólidos, un cuerpo es un conjunto infinito de partículas, cada una de las cuales ocupa una posición en el espacio. Estas posiciones son variables en el tiempo y al conjunto de ellas en un instante dado se denomina configuración. Una configuración puede entenderse de manera intuitiva como una foto fija del solido a lo largo del proceso de carga. Considérese un cuerpo deformable en el espacio. En la configuración inicial, cada partícula del solido ocupa una posición en un sistema rectangular de coordenadas. Nos referimos a esta configuración como configuración original o inicial. Cuando el cuerpo se deforma, cada partícula ocupa una nueva posición que viene descrita por unas nuevas coordenadas. Llamaremos a esta configuración deformada o actual. De esta forma, por ejemplo, una partícula ºP de posición inicial (ºx1, ºx2, ºx3) se mueve a una posición deformada tP de coordenadas (tx1, tx2, tx3). El vector desplazamiento ºPtP de la partícula tendrá como componentes 𝑢1= 𝑥1 𝑡− 𝑥1 0 𝑢2= 𝑥2 𝑡− 𝑥2 0 𝑢3= 𝑥3 𝑡− 𝑥3 0 (Ec5.) Si se conoce el desplazamiento de cada partícula del sólido, puede construirse su configuración deformada. Esto sería posible si ( 𝑥1 𝑡,𝑥2 𝑡,𝑥3 𝑡) fueran funciones conocidas de las coordenadas iniciales (𝑥1 0,𝑥2 0,𝑥3 0). 𝑥1 𝑡= 𝑥1 𝑡(𝑥1 0,𝑥2 0,𝑥3 0) 𝑥2 𝑡= 𝑥2 𝑡(𝑥1 0,𝑥2 0,𝑥3 0) 𝑥3 𝑡= 𝑥3 𝑡(𝑥1 0,𝑥2 0,𝑥3 0) (Ec6.) En la mecánica de sólidos se supone que estas funciones son continuas y biunívocas, con lo que se están excluyendo de esta forma, interpenetraciones de unas partes del sólido en otras o discontinuidades producidas por fisuras. Esto implica también que las funciones ( 𝑥1 𝑡,𝑥2 𝑡,𝑥3 𝑡) tienen inversas. 𝑥1 0= 𝑥1 0(𝑥1 𝑡,𝑥2 𝑡,𝑥3 𝑡) 𝑥2 0= 𝑥2 0(𝑥1 𝑡,𝑥2 𝑡,𝑥3 𝑡) 𝑥3 0= 𝑥3 0(𝑥1 𝑡,𝑥2 𝑡,𝑥3 𝑡) (Ec7.) 9 El vector de desplazamiento puede asociarse a la partícula en la posición inicial: 𝑢1=𝑢1 (𝑥1 0,𝑥2 0,𝑥3 0)= 𝑥1 𝑡(𝑥1 0,𝑥2 0,𝑥3 0)− 𝑥1 0 𝑢2= 𝑢2 (𝑥1 0,𝑥2 0,𝑥3 0)= 𝑥2 𝑡(𝑥1 0,𝑥2 0,𝑥3 0)− 𝑥2 0 𝑢3= 𝑢3 (𝑥1 0,𝑥2 0,𝑥3 0)= 𝑥3 𝑡(𝑥1 0,𝑥2 0,𝑥3 0)− 𝑥3 0 (Ec8.) o a la partícula en la posición deformada: 𝑢1= 𝑢1 (𝑥1 𝑡,𝑥2 𝑡,𝑥3 𝑡)=𝑥1 𝑡− 𝑥1 0(𝑥1 𝑡,𝑥2 𝑡,𝑥3 𝑡) 𝑢2= 𝑢2 (𝑥1 𝑡,𝑥2 𝑡,𝑥3 𝑡)=𝑥2 𝑡− 𝑥2 0(𝑥1 𝑡,𝑥2 𝑡,𝑥3 𝑡) 𝑢3= 𝑢3 (𝑥1 𝑡,𝑥2 𝑡,𝑥3 𝑡)=𝑥3 𝑡− 𝑥3 0(𝑥1 𝑡,𝑥2 𝑡,𝑥3 𝑡) (Ec9.) A la primera forma de describir el movimiento se le conoce como descripción o formulación Lagrangiana o material y es la usual en Mecánica de Sólidos, mientras que la segunda se utiliza con más frecuencia en Mecánica de Fluidos y se conoce como formulación Euleriana o espacial. Se puede decir que la formulación Lagrangiana se ocupa de lo que le sucede a una partícula material mientras que la formulación Euleriana expresa lo que le sucede a una cierta posición del espacio. Los desarrollos que se siguen se basan en la formulación Lagrangiana. 5.2. Variantes de la formulación Lagrangiana. Dentro de la descripción Lagrangiana del movimiento de un sólido caben varias posibilidades a la hora de definir la configuración que se toma como referencia. Es importante no confundir configuración de referencia con configuración inicial, ya que esta distinción es lo que permite precisamente caracterizar las distintas variantes de la formulación Lagrangiana. Estas variantes son las siguientes: 5.2.1. Formulación Lagrangiana Total. La configuración de referencia se mantiene cte. durante el proceso completo de deformación. La configuración de referencia es la configuración de la estructura original o indeformada. Tanto las tensiones como las deformaciones han de ser medidas con respecto a esta configuración. 5.2.2. Formulación Lagrangiana Actualizada. Si el proceso de análisis se lleva a cabo, como es frecuente, mediante incrementos sucesivos de carga, resulta interesante a veces utilizar como configuración de referencia la configuración alcanzada en el incremento anterior. Cada configuración pasa a ser, por tanto, referencia de la siguiente. 5.2.3. Formulación Lagrangiana Corrotacional. En el caso de estructuras formadas por elementos esbeltos, como barras, placas o laminas, se pueden producir desplazamientos y rotaciones considerables aunque las deformaciones se mantengan moderadas. En estos casos, la componente de sólido rígido del movimiento es muy importante y en este caso resulta conveniente utilizar una configuración de referencia que elimine dicha componente. Por ejemplo, en el caso de barras se toma como referencia la barra indeformada pero trasladada y girada en el espacio mediante un movimiento de sólido rígido, tal que la posición de sus extremos coincida con la posición de los extremos de la barra deformada en el instante 16 Los métodos iterativos de solución que hemos visto en el punto anterior, fallan por tanto en la proximidad de los puntos límite. Numéricamente se observa que, al acercarse el valor de las cargas al valor máximo, la convergencia se reduce notablemente. Si se ha previsto un valor de carga superior al máximo, el proceso deja de ser convergente y no se puede obtener una solución, tal y como se representa en la figura 6. Figura 6: Con λ = cte. El método iterativo no converge. Otro tipo de puntos que aparecen en la figura 5 son los llamados puntos de retroceso. En estos puntos, uno o varios grados de libertad disminuyen de valor a medida que avanza el proceso de carga. En la figura 5, se representa también un punto de bifurcación del camino de equilibrio. En este tipo de puntos aparecen varias ramas posibles para la curva de equilibrio. El que la estructura siga un camino determinado depende de ciertos parámetros de imperfección inicial. Estos parámetros pueden ser bien fuerzas o bien imperfecciones geométricas o mecánicas. 5.3.3. Métodos de longitud del arco. Con objeto de poder seguir el comportamiento de la estructura en todo su ámbito de respuesta, interesa desarrollar métodos que no fallen en la proximidad de puntos críticos. Ya se ha visto que el proceso iterativo de corrección puede no ser convergente si se mantiene constante el factor multiplicador de las cargas λ. Una manera de evitar el inconveniente citado consiste en tratar el factor λ como una variable más del problema, de forma que el sistema de ecuaciones de equilibrio tenga siempre solución. Hará falta para ello imponer una condición adicional a las n ecuaciones de equilibrio, ya que existe una nueva incógnita λ a añadir a los n grados de libertad del problema. Si se considera un espacio (𝑼,λ) de dimensión n+1, esta nueva condición puede ser entendida como una superficie que constituye el lugar geométrico de los puntos (𝑈𝑖,λ𝑖) obtenidos en el proceso iterativo de corrección. A su vez, la curva de carga es el lugar geométrico de los puntos que cumplen las ecuaciones de equilibrio. 17 La condición usual de los métodos incrementales de control de carga es hacer λ=cte. Esta condición es la ecuación de un plano ortogonal al eje 𝑂−λ que puede cortar o no a la curva de carga. Si el nivel de carga corresponde a cargas menores que las máximas que puede soportar la estructura, el plano corta a la curva de carga y las iteraciones convergen. Si las cargas son mayores, el plano no corta a la curva de carga y el proceso es divergente. A pesar de que existen otras ecuaciones de restricción, en nuestro código hemos empleado la que se conoce como vínculo esférico. Una manera de conseguir la intersección de la curva de carga con la superficie que representa la condición adicional y, por tanto, una manera de garantizar la convergencia del proceso iterativo de corrección consiste en imponer como condición adicional superficies cerradas alrededor del punto de equilibrio determinado anteriormente. La más sencilla de estas superficies es una (híper) esfera de radio Δ𝑙 (Figura 7) cuya ecuación viene dada por: �Δ𝑼ik�TΔ𝑼ik + �α Δλik�2= Δ𝑙2 (Ec33.) En esta ecuación α es un coeficiente numérico que se introduce al objeto de homogeneizar las dimensiones y magnitud numérica de los términos de la ecuación pero en la mayoría de las ocasiones se toma como 1. Esta condición fue propuesta por Riks (ref.14) y ha sido el origen de diversos métodos de solución de problemas no lineales. Figura 7: Δ𝑙 = cte. Vinculo esférico. Del mismo modo a lo hecho en las ecuaciones Ec25 a Ec32, en este caso el residuo será función no solo de 𝑼 sino también de 𝜆 y podremos volver a escribir la ecuación Ec25, que al final de las k iteraciones del salto de carga i 𝑹𝑖𝑘= 𝜆𝑖𝑘 𝑷−𝑸𝑖𝑘(𝑼𝑖𝑘) = 0 (Ec25.) 18 Desarrollando el residuo en serie de potencias alrededor del punto de iteración anterior (k-1) se obtiene 𝑹𝑖𝑘≈ 𝑹𝑖𝑘−1+� 𝜕𝑹𝑖𝑘−1 𝜕𝜆 �δ𝜆𝑘+� 𝜕𝑹𝑖𝑘−1 𝜕𝑼 �δ𝑼𝑘= 0 (Ec34.) En esta expresión δ𝑼𝑘 es el incremento de deformación producido en la iteración k del salto de carga i, y del mismo modo δ𝜆𝑘 es el incremento de deformación producido en la iteración k de ese mismo salto de carga. Las derivadas necesarias son: 𝜕𝑹𝑖𝑘−1 𝜕𝜆 =𝑷 (Ec35.) y 𝜕𝑹𝑖𝑘−1 𝜕𝑼 =−𝜕𝑸𝑖𝑘−1 𝜕𝑼 =−𝑲𝑖𝑘−1 (Ec36.) Luego la ecuación de iteración a resolver, que es la que resolvemos en nuestro código es: 𝑹𝑖𝑘≈ 𝑹𝑖𝑘−1+� 𝜕𝑹𝑖𝑘−1 𝜕𝜆 �δ𝜆𝑘+� 𝜕𝑹𝑖𝑘−1 𝜕𝑼 �δ𝑼𝑘=𝑹𝑖𝑘−1+𝑷 δ𝜆𝑘−𝑲𝑖𝑘−1 δ𝑼𝑘= 0 𝑲𝑖𝑘−1 δ𝑼𝑘=𝑹𝑖𝑘−1+𝑷 δ𝜆𝑘 (Ec37.) Despejando el incremento de deformación se obtiene δ𝑼𝑘=�𝑲𝑖𝑘−1�−1 𝑹𝑖𝑘−1+�𝑲𝑖𝑘−1�−1𝑷 δ𝜆𝑘 (Ec38.) Que podemos escribir como: δ𝑼𝑘=�𝑲𝑖𝑘−1�−1 𝑹𝑖𝑘−1+�𝑲𝑖𝑘−1�−1𝑷 δ𝜆𝑘 (Ec39.) Por lo tanto la ecuación a resolver en la iteración k es δ𝑼𝑘= 𝑼 �𝑖𝑘+𝑼 �𝑖𝑘 δ𝜆𝑘= 0 (Ec40.) El primer sumando del incremento del desplazamiento se puede calcular con facilidad y representa el desplazamiento producido por la parte del residuo no equilibrado en la iteración anterior 𝑼 �𝑖𝑘= �𝑲𝑖𝑘−1�−1 𝑹𝑖𝑘−1 (Ec41.) El segundo sumando no puede evaluarse hasta no conocer el valor de 𝜆 pero su coeficiente puede evaluarse con sencillez y representa la deformación producida por las fuerzas básicas 𝑷. 19 𝑼 �𝑖𝑘= �𝑲𝑖𝑘−1�−1𝑷 (Ec42.) Obsérvese que si se emplea el método de Newton-Raphson modificado, no es necesario calcular este término a cada paso de la iteración, sino que puede mantenerse el del primero. Suponiendo por el momento conocido el valor de δ𝜆𝑘, y por tanto el valor del incremento de desplazamientos en esta iteración δ𝑼𝑘, se procede a actualizar los valores de las incognitas. Para los desplazamientos la actualización es: Δ𝑼𝑖𝑘 = Δ𝑼𝑖𝑘−1+δ𝑼𝑘 (Ec43.) Donde Δ𝑼𝑖𝑘−1 es el incremento de deformación acumulado a lo largo de las (k-1) iteraciones anteriores. De forma similar Δ𝑼𝑖𝑘 es el incremento de deformación acumulado tras la iteración k. De manera análoga se actualiza el parámetro de carga 𝜆: λ𝑖𝑘 = λ𝑖𝑘−1+δλ𝑘 (Ec44.) El cálculo de δλ𝑘 se efectua introduciendo una ecuación de restricción que imponga la condición de distancia máxima recorrida en este salto de carga, es decir que se limita el incremento de desplazamiento acumulado en todas las iteraciones efectuadas en este salto de carga. Si se denomina 𝛥𝑠 a la distancia máxima a recorrer, la condición es: (𝛥𝑠)2 = (Δ𝑼𝑖𝑘)𝑇 Δ𝑼𝑖𝑘 (Ec45.) Sustituyendo los incrementos por sus valores y operando se obtiene (𝛥𝑠)2 = (Δ𝑼𝑖𝑘−1+δ𝑼𝑘)𝑇 (Δ𝑼𝑖𝑘−1+δ𝑼𝑘) (𝛥𝑠)2 = (Δ𝑼𝑖𝑘−1+𝑼 �𝑖𝑘+𝑼 �𝑖𝑘 δ𝜆𝑘)𝑇 (Δ𝑼𝑖𝑘−1+𝑼 �𝑖𝑘+𝑼 �𝑖𝑘 δ𝜆𝑘) (𝛥𝑠)2 = (Δ𝑼𝑖𝑘−1+𝑼 �𝑖𝑘)𝑇 �Δ𝑼𝑖𝑘−1+𝑼 �𝑖𝑘�+ 2 (Δ𝑼𝑖𝑘−1+𝑼 �𝑖𝑘)𝑇 𝑼 �𝑖𝑘 δ𝜆𝑘+ (𝑼 �𝑖𝑘)𝑇 𝑼 �𝑖𝑘 (δ𝜆𝑘)2 La ecuación anterior es una ecuación de segundo grado en λ : a1(δ𝜆𝑘)2+ a2 δ𝜆𝑘+ a3= 0 donde a1= (𝑼 �𝑖𝑘)𝑇 𝑼 �𝑖𝑘 a2= 2 (Δ𝑼𝑖𝑘−1+𝑼 �𝑖𝑘)𝑇 𝑼 �𝑖𝑘 a3= (Δ𝑼𝑖𝑘−1+𝑼 �𝑖𝑘)𝑇 �Δ𝑼𝑖𝑘−1+𝑼 �𝑖𝑘�−(𝛥𝑠)2 Resolviendo esta ecuación se obtienen dos raíces δ𝜆1 𝑘 y δ𝜆2 𝑘. De entre ellas se elige aquella que producirá un incremento de desplazamiento acumulado más próximo al incremento de desplazamiento acumulado en la iteración anterior. Para ello, en primer lugar se determina cuál sería el incremento de desplazamiento producido por cada una de las soluciones Δ𝑼𝑖(1) 𝑘 = Δ𝑼𝑖𝑘−1+𝑼 �𝑖𝑘+𝑼 �𝑖𝑘 δ𝜆1 𝑘 Δ𝑼𝑖(2) 𝑘 = Δ𝑼𝑖𝑘−1+𝑼 �𝑖𝑘+𝑼 �𝑖𝑘 δ𝜆2 𝑘 20 A continuación se calcula la proyección de dichos incrementos de desplazamiento sobre el incremento de desplazamiento de la iteración anterior, que será un escalar, que de alguna manera estima el ángulo entre ambos vectores: (𝛥𝑠)2cos𝜑1 = (Δ𝑼𝑖𝑘−1)𝑇 Δ𝑼𝑖(1) 𝑘= (Δ𝑼𝑖𝑘−1)𝑇 (Δ𝑼𝑖𝑘−1+𝑼 �𝑖𝑘+𝑼 �𝑖𝑘 δ𝜆1 𝑘) (𝛥𝑠)2cos𝜑2 = (Δ𝑼𝑖𝑘−1)𝑇 Δ𝑼𝑖(2) 𝑘= (Δ𝑼𝑖𝑘−1)𝑇 (Δ𝑼𝑖𝑘−1+𝑼 �𝑖𝑘+𝑼 �𝑖𝑘 δ𝜆2 𝑘) (Ec46.) Eligiendo aquella solución que produzca menor ángulo, es decir el mayor valor del cos𝜑. 5.4. Determinación del tamaño del salto de carga. La idea de tratar el factor multiplicador de las cargas como una incógnita adicional en el proceso de corrección de la solución, puede además ser utilizada para decidir la estrategia de avance en el proceso de seguimiento de la curva de carga. De esta forma, en lugar de utilizar saltos de carga Δλi𝑷 determinados de antemano, puede llevarse a cabo un control del proceso de incrementación mediante otro criterio, tal como fijar la longitud del vector tangente (Δλi0,Δ𝑼i0) correspondiente a la primera aproximación dentro de cada salto mediante la ecuación (Ec25.) particularizada. (Δ𝑼i0)TΔ𝑼i0 + (α Δλi0)2= Δ𝑙2 = 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡𝑒 (Ec47.) Ya se comprende que en la práctica resulta difícil estimar un valor adecuado para el valor de la longitud del arco Δ𝑙, por lo que en la práctica se suele fijar un valor inicial de Δλi0 para el primer incremento. A partir de este valor se puede calcular Δ𝑙 despejando: Δ𝑙 = �(Δ𝑼i0)TΔ𝑼i0 + (α Δλi0)2 (Ec48.) Siendo Δ𝑼i0=Δλi0 𝑼𝑖∗ y calculándose 𝑼𝑖∗ resolviendo el sistema de ecuaciones: 𝑲0 𝑼𝑖∗=𝑷 (Ec49.) Una vez se ha determinado en el primer salto de carga el valor de la longitud del arco Δ𝑙, en los saltos posteriores, se puede deducir de manera inversa el valor del salto de carga que da lugar a un vector tangente cuya longitud valga Δ𝑙 a partir de la siguiente expresión: Δλi0 = ±�Δ𝑙2 (𝑼i∗0)T𝑼i∗0 + α2 (Ec50.) En la expresión anterior se toma el signo positivo o negativo de acuerdo con el signo del determinante de la matriz de rigidez 𝑲𝑖 utilizada al principio de cada salto para obtener 𝑼𝑖∗. De esta forma es posible seguir las ramas descendentes de las trayectorias de equilibrio, en las que es preciso utilizar factores de carga negativos. Aunque nosotros en nuestro código calculamos el determinante de la matriz de rigidez tangente de la estructura con la orden Det [K] incorporada en Mathematica, el signo del determinante de la matriz de rigidez se puede calcular sin necesidad de calcular el propio determinante. Si se utilizan algoritmos de resolución del 21 sistema de ecuaciones basados en la factorización de la matriz de rigidez en el producto de una matriz triangular inferior, con términos unidad en la diagonal principal, por otra triangular superior (factorización de Crout o también llamada factorización de Cholesky modificada), el signo del determinante se puede calcular de manera muy eficiente. La matriz de rigidez se descompone según: 𝑲=𝑳𝑼 (Ec51.) Siendo 𝑳 la matriz triangular inferior y 𝑼 la triangular superior. El signo del determinante coincide con el signo del producto de los términos de la diagonal de la matriz triangular superior: 𝑠𝑖𝑔𝑛𝑜�𝑑𝑒𝑡(𝑲)�=𝑠𝑖𝑔𝑛𝑜�𝑑𝑒𝑡(𝑳)∗𝑑𝑒𝑡(𝑼)�=𝑠𝑖𝑔𝑛𝑜�∏𝑈𝑖𝑖 𝑛 1 |∏𝑈𝑖𝑖 𝑛 1|�=�𝑈𝑖𝑖 |𝑈𝑖𝑖| 𝑛 1 (Ec52.) 5.5. Ajuste automático de la longitud del salto de carga. En el apartado anterior se ha supuesto que la longitud del vector tangente Δ𝑙 se mantenía constante en el proceso de carga. Se comprende fácilmente que manteniendo constante Δ𝑙, el proceso de corrección iterativo requerirá normalmente mayor número de iteraciones para conseguir una cierta aproximación en zonas de la trayectoria de equilibrio con gran curvatura que en otras zonas de curvatura menor. Nótese que la curvatura de la curva de equilibrio está relacionada con la variación de la rigidez de la estructura, de manera que a mayores variaciones de rigidez corresponderán mayor número de iteraciones. Si se quiere que el número de estas sea más o menos constante se puede variar Δ𝑙. Se trata de modificar esta longitud en función del número de iteraciones 𝜃𝑖−1 que han hecho falta para conseguir una cierta aproximación en el salto de carga anterior y el número de iteraciones que se desea mantener constante 𝜃𝑑𝑒𝑠. Una fórmula debida a Ramm (ref. 15) es: Δ𝑙𝑖=Δ𝑙𝑖−1 ( 𝜃𝑑𝑒𝑠 𝜃𝑖−1 )𝜸 (Ec53.) con γ igual a 0,5. Este sencillo procedimiento produce automáticamente pequeños saltos de carga en zonas de comportamiento altamente no lineal y saltos de carga mayores en zonas de comportamiento “casi” lineal. 5.6. Cálculo de puntos de inestabilidad. Los métodos que se acaban de exponer permiten el seguimiento completo del comportamiento no lineal de una estructura. El objetivo práctico de tal seguimiento es, muchas veces, estimar el valor de las cargas máximas que puede soportar la estructura. En una parte importante de los problemas de no linealidad geométrica esta estimación puede llevarse a cabo sin necesidad de trazar completamente la curva de equilibrio, evitándose así un proceso de elevado coste en términos de cálculo. Tanto los puntos límite como los puntos de bifurcación del equilibrio significan desde el punto de vista numérico puntos críticos de la curva de equilibrio en los que el sistema de ecuaciones de equilibrio esta indeterminado. En estos puntos, la matriz de rigidez tangente deja de ser definida positiva. La detección de las situaciones en las que la matriz de rigidez se hace singular constituye la base de los métodos de cálculo de puntos críticos. 22 5.7. Criterios de convergencia. Con objeto de decidir cuándo se termina el proceso iterativo se deben establecer ciertos criterios de convergencia los cuales pueden clasificarse dependiendo de las cantidades que intervienen en la comparación. Se tienen así los dos criterios más usados: 5.7.1. Criterio de convergencia de desplazamientos Con este criterio, en cada salto de carga, el proceso iterativo finaliza cuando la variación de desplazamientos δuk que se obtiene es menor que un cierto valor preestablecido. Ya que δuk es una magnitud vectorial de tantas componentes como grados de libertad tiene la estructura, para poder establecer comparaciones se utilizan normas de tres tipos. • 1. Norma infinita ‖δ𝑼‖∞= 𝑚𝑎𝑥 |δ𝑢𝑖| (Ec54.) • 2. Norma L1 ‖δ𝑼‖1= �|δ𝑢𝑖| (Ec55.) • 3. Norma L2 ‖δ𝑼‖2= ��δ𝑢𝑖2�1/2 (Ec56.) La primera de estas normas, ‖δU‖∞ es simplemente el valor máximo de las componentes del vector variación de desplazamiento. ‖δ𝑼‖1 es la suma de los valores absolutos de dichas componentes. ‖δ𝑼‖2 es la raíz cuadrada de la suma de los cuadrados de las componentes y equivale al módulo del vector. A esta última se la llama también norma Euclídea y es por lo general la más usada. Estas normas se comparan con una fracción de la norma correspondiente de los desplazamientos totales o del incremento de los desplazamientos. Por ejemplo si se utiliza la norma L2 se tendría: ‖δ𝑼‖2<𝜺𝒖‖𝑼‖2 (Ec57.) o si se compara con el incremento de desplazamiento: ‖δ𝑼‖2<𝜺𝒖‖Δ𝑼‖2 (Ec58.) siendo εu un valor que frecuentemente se fija en 0.01 ó 0.001. 5.7.2. Criterio de convergencia de fuerzas Se utiliza como cantidad a comparar una norma de las fuerzas residuales 𝑹. Como cantidad de referencia se puede utilizar la norma correspondiente de las fuerzas exteriores totales o la del incremento de carga. Utilizando por ejemplo la norma L2 del vector de incrementos de carga se tiene el siguiente criterio: 23 ‖𝑹‖2< 𝜺𝒈 ‖Δλ 𝑷‖2 (Ec59.) Un criterio también utilizado con frecuencia toma como cantidad de referencia el valor máximo de los residuos en el incremento: ‖𝑹‖2< 𝜺𝒈 𝑚𝑎𝑥 (‖𝑹‖2) (Ec60.) 5.7.3. Criterio de convergencia de la energía Otro criterio es el que considera en forma conjunta desplazamientos y fuerzas. Se basa en utilizar como cantidad a comparar al valor absoluto del producto escalar de la variación de desplazamientos por las fuerzas residuales. Como referencia se toma el valor absoluto del producto escalar de los desplazamientos totales por las cargas exteriores: ‖𝛿𝑼𝑇 𝑹‖2< 𝜺𝒆 |𝑼𝑇 𝑷| (Ec61.) 5.8. Matriz de rigidez tangente de la barra. Una vez visto lo anterior, para poder determinar la respuesta de la estructura se hace necesario ver como determinamos la matriz de rigidez tangente que según sabemos, en los problemas con algún tipo de no linealidad, y en particular de la geométrica depende del nivel de desplazamientos al que se ve sometida la estructura. 5.8.1. Formulación de la matriz de rigidez tangente de la barra. Sea la barra 1-2 de la figura 1, que inicialmente forma un ángulo β0 con el eje x global. A esta barra, de longitud inicial L0, tras aplicarle cargas en sus extremos (fuerzas y momentos), se deforma pasado dichos extremos a ocupar la posición 1´- 2´. Fig 1. En principio las incógnitas del problema son precisamente los desplazamientos nodales (movimientos U, V y 𝜃 de sus nudos extremos), medidos estos en ejes globales. 24 Los movimientos de los nudos extremos de la barra los podemos agrupar en un vector columna que denominamos vector desplazamiento y que definimos del siguiente modo 𝑼= {U1, V1,θ1, U2, V2,θ2}T (Ec62.) Por otro lado, definimos los ángulos ϑ1 y ϑ2 como ϑ1=α+𝜗1 ϑ2=α+𝜗2 (Ec63.) donde α es el giro de la barra como sólido rígido. Para ese movimiento, el alargamiento de la barra valdrá u�= L −L0 (Ec64.) y la deformación ingenieril la podremos expresar como ε=L−L0 L0=u� L0 (Ec65.) Así, el axil en la barra será por tanto N = σ∗A0=ε∗E∗A0=u�∗E∗A0 L0 (Ec66.) La variación virtual del alargamiento u� de la barra, se puede expresar en función de los movimientos de los nudos extremos de la barra como Fig 2. δu�=δU2cos β+δV2sen β−δU1cos β−δV1 sen β (Ec67.) 25 o bien δu�={−cos β,−senβ , 0, cos β,sen β, 0} δ𝑼 (Ec68.) δu�=𝒓𝑇 δ𝑼 (Ec69.) donde β es el ángulo que forma la barra con el eje x global. De la misma manera, en incremento virtual del ángulo (δα) lo podemos expresar como δα=δβ=δv�2 L0=1 L0∗{senβ,−cosβ,0,−senβ,cosβ, 0} δ𝑼 (Ec70.) δα=1 L0∗𝒛𝑇δ𝑼 (Ec71.) Por otro lado δϑ1=δα+δ𝜗1 δϑ2=δα+δ𝜗2 (Ec72.) con lo que �δ𝜗1 δ𝜗2�=�δϑ1 δϑ2�−�δα δα� (Ec73.) �δ𝜗1 δ𝜗2�=�001 000 0 0 0 0 0 1�δ𝑼−1 L0∗�𝒛𝑇 𝒛𝑇� δ𝑼=𝑨 δ𝑼 (Ec74.) donde 𝑨=�0 0 1 0 0 0 0 0 0 0 0 1�−1 L0∗�𝒛𝑇 𝒛𝑇� (Ec75.) Llamando vector de movimientos virtuales al vector δp� δ𝒑 �=�δu� δ𝜗1 δ𝜗2�=�𝒓𝑇 𝑨�δ𝑼=𝑩 δ𝑼 (Ec76.) definimos la matriz B como 𝑩=�𝒓𝑇 𝑨� (Ec77.) 32 • A partir del punto (-5, 0), para poder seguir aumentando el descenso del punto de aplicación de la carga debemos disminuir esta (empujar hacia abajo). En este punto, la estructura va rigidizándose a medida que va tomando carga como consecuencia de la tracción que le provoca el descenso del punto de aplicación de la carga. Tenemos que señalar que con la geometría particular de este primer caso, la trayectoria de equilibrio obtenida describe una curva suave. En este caso, el código programado no ha dado ningún problema a la hora de obtener los sucesivos puntos de equilibrio y no se han producido puntos de oscilación, raíces imaginarias ni cualquier otro tipo de inestabilidad. El hecho de que la curva sea suave es indicativo de que el pandeo se ha producido a nivel global de la estructura y no a nivel local de barra. Notemos que con esta geometría, el acortamiento máximo de la barra se presenta cuando esta toma la posición horizontal. Inicialmente, la longitud de la barra vale 𝐿0=�𝑑𝑥2+𝑑𝑦2=�1002+ 52=100,1249 𝑐𝑚 Mientras que la longitud de esta cuando el descenso es de 5 cm vale y admitiendo que permanece recta es 𝐿𝑓=𝑑𝑥=100 𝑐𝑚 Con esto, es acortamiento máximo de la barra vale u = 𝐿𝑓−𝐿0=100 −100,1249 =−0,1249 𝑐𝑚 por lo que el alargamiento de sus fibras es de ε=u 𝐿0=−0,1249 100,1249 =−1,247 ∗10−3 lo que da lugar a una carga de compresión de valor N = σ∗A = ε∗E∗A = −1,247 ∗10−3∗2.100.000 ∗20 = −52.392,6 𝑘𝑝 Por otro lado, como se trata de una barra apoyada-apoyada, la carga crítica de pandeo de la barra (estructura) es: N𝑐𝑟=−𝜋2∗𝐸∗𝐼 𝐿𝑘2=−3,1412∗2.100.000 ∗49,9125 100,1249 =−103.191,6 𝑘𝑝 es decir, el axil producido como consecuencia del acortamiento de la barra al descender el punto de aplicación de la carga 5 cm (N = -52.393 kp), es menor que la carga critica de pandeo de Euler de la barra. Por lo tanto, en este caso, mucho antes de que se produzca el pandeo de Euler, la estructura ha perdido su rigidez global siendo en este caso ésta la que gobierna la física del problema. No ocurre lo mismo si la coordenada y del extremo derecho de la barra vale 10 cm en lugar de los 5 cm de este primer caso. 33 6.2. Barra inclinada (caso 2). En este segundo problema, contamos con el mismo tipo de perfil que en el caso anterior, con el mismo material y sección pero ahora, la geometría inicial es diferente. Ahora dy son 10 cm mientras que dx se mantiene constante e igual a 100 cm como en el caso anterior. Las condiciones de contorno son las mismas que antes por lo que la única diferencia es la longitud de la barra (estructura). Por tanto, en este segundo caso, la barra al igual que en el caso anterior, tiene un área de 19,8 cm2 y una inercia de 49,9125 cm4. El extremo izquierdo de la barra se encuentra apoyado (sin posibilidad de movimiento ni en dirección x ni en dirección y) mientras que el extremo derecho solo tiene impedido el movimiento en dirección x, pudiendo dicho punto moverse en dirección vertical. Ambos extremos pueden rotar. Para la obtención de la trayectoria de equilibrio, en este caso se han llevado a cabo 265 incrementos de carga con los siguientes parámetros (𝛾 = 0,5, 𝜆0 = 10, 𝜃𝑑𝑒𝑠 = 5, tol = 0.001, máx. nº iter = 40). Figura 13: Geometría del modelo. Figura 14: Discretización del modelo en Mathematica. La barra se divide también en 10 elementos de igual longitud y con esta geometría y condiciones de contorno, se aplica al extremo derecho una carga vertical variable –P, la cual se grafica frente a la posición vertical del punto de aplicación de la carga. En este caso se observa lo siguiente. La curva de equilibrio que se obtiene adopta la forma de la figura 15. 20 40 60 80 100 2 4 6 8 10 34 Figura 15: Gráfica y - 𝜆. Posición del punto de aplicación de la carga - Carga aplicada. En la figura 15 puede apreciarse lo siguiente: • Los puntos (y, 𝜆) obtenidos con el código, son puntos de equilibrio y la línea que los une, es la “trayectoria de equilibrio de la estructura”. En realidad la trayectoria de equilibrio es una curva alabeada en un espacio n+1 dimensional, siendo n el nº de grados de libertad del sistema. • Inicialmente nos encontramos en el punto (10, 0). La estructura está descargada y la barra se encuentra en un estado libre de esfuerzos. En esta situación, la barra se encuentra en un estado de equilibrio estable. • Conforme disminuimos el valor de la carga P, el punto de aplicación de esta va descendiendo y el problema, que inicialmente es prácticamente lineal, se hace cada vez más no lineal observándose claramente en la figura la no linealidad entre cargas aplicadas y descensos. • A medida que la carga va aumentando, la estructura va perdiendo rigidez como consecuencia del axil que acumula debido a su acortamiento y cuando la carga es aproximadamente de unos 7.280 kp (𝜆 = -7.280), aparece un punto de bifurcación del equilibrio, en el que la estructura bien puede pandear localmente a nivel de barra o bien si la barra no pierde la rectitud, puede seguir incrementándose la carga para describir un comportamiento similar al del caso 1 con un pandeo global. Notas: En este caso, para cazar la trayectoria a partir del punto de bifurcación, se ha hecho necesario introducir una pequeña perturbación capaz de modificar la rectitud de la barra y de esa manera excitar el modo de pandeo local de la barra. Dicha perturbación se ha introducido como una carga puntual de pequeña magnitud, aplicada en el nodo central de la partición y en dirección vertical. Si no introducimos la perturbación, al llegar al punto de bifurcación el código oscila no siendo capaz de determinar las nuevas posiciones de equilibrio y este hecho hace que necesitemos de cierto conocimiento previo del modo de pandeo que puede tener lugar. En este caso es sencillo al tratarse de una sola barra y la primera forma de pandeo de ésta en sencillamente la de la figura 16. En problemas más complejos (con más barras), la solución puede no ser tan obvia. 35 Figura 16: Primer modo de pandeo de la barra articulada-articulada. Notemos que en este caso, la carga crítica de pandeo de la barra vale: N𝑐𝑟=−𝜋2∗𝐸∗𝐼 𝐿𝑘2=−3,1412∗2.100.000 ∗49,9125 100,4987 =−102.425,0 𝑘𝑝 Para que se produzca ese axil en la barra, la deformación de la barra ha de ser la que sigue N = −102.425,0 𝑘𝑝=σ∗A = ε∗E∗A ⇒ε= N E∗A=−102.425,0 2.100.000 ∗20 =−2,4387 ∗10−3 de donde se deduce que el acortamiento de la barra tiene que valer u = ε∗𝐿0=−2,4387 ∗10−3∗100,4987 =−0,2451 𝑐𝑚 Para que se produzca ese acortamiento u en la barra, el punto de aplicación de la carga se tiene que mover verticalmente una cantidad (y-y0) y la longitud final de la barra valdrá 𝐿= 𝐿0+𝑢=100,4987 −0,2451 =100,2536 𝑐𝑚 y de aquí, como 𝑑𝑦=�𝐿2−𝑑𝑥2=�100,25362+1002= 7,1263 𝑐𝑚 el descenso Δ𝑦=𝑑𝑦−𝑑𝑦0= 7,1263 −10 =−2,8737 𝑐𝑚 Es decir, la carga critica de pandeo local de la barra (Ncr = -102.425 kp), aparece cuando el descenso del punto de aplicación de la carga es de 2,8737 cm lo cual se corresponde con el desplazamiento vertical en el punto de bifurcación que aparece en la figura 15. Para esa posición de la barra, el ángulo 𝛽 que forma con la horizontal es de tan 𝛽=dy 𝑑𝑥=7,1263 100 = 0,07126 𝛽=arctg(0,07126)= 4,0762 𝑜 Por lo que la carga vertical aplicada vale 𝑃= N ∗sen 𝛽=−102.425,0 ∗𝑠𝑒𝑛 4,0762 =−7.280,0 𝑘𝑝 La cual coincide con lo que se observa en la figura 15. Además, esta carga es menor que la carga de pandeo global de la estructura por lo que cuando dy = 7,1263 cm, lo que dará lugar a que la carga en la barra sea de 102.425 kp, si la barra pierde la rectitud ideal se producirá un pandeo local de esta que marcará la trayectoria de equilibrio indicada en la figura 15. Si la barra no pierde la rectitud, la trayectoria de equilibrio será similar a la de la figura 12 del caso 1. 36 • A partir del punto de bifurcación (7,1263, -7.280), para encontrar sucesivas posiciones de equilibrio el valor absoluto de la carga P debe de disminuir. De hecho, las sucesivas posiciones de equilibrio están contenidas en una recta que pasando por el (0, 0) llega hasta la posición simétrica (-7,1263, 7.280). Durante todo este tramo recto, la rigidez tangente de la estructura es negativa. La carga de pandeo de la barra permanece constante (esta no cambia) pero la carga vertical aplicada P aumenta hasta anularse en el (0, 0) y continua creciendo hasta llegar al punto de simetría (-7,1263, 7.280). • En el punto (-7,1263, 7.280) la rigidez de la estructura pasa a ser nuevamente positiva. Para poder seguir aumentando el descenso del punto de aplicación de la carga debemos de disminuir la carga vertical aplicada P hasta el momento en el que la barra adopte la posición simétrica respecto de su posición original, momento este en el que la carga se hace nula y que corresponde al punto (-10, 0). A partir del punto (-10, 0), para poder seguir aumentando el descenso del punto de aplicación de la carga debemos de aumentar esta. En este punto, la estructura a medida que va tomando carga va rigidizándose como consecuencia de la tracción que le provoca el descenso del punto de aplicación de la carga. 6.3. Pórtico simple. En este tercer problema, se compara el resultado obtenido en nuestro código con el que obtienen P. Kotronis y F. Collin en su artículo Implementation of path following techniques into the finite element code Lagamine. En este ejemplo se trata un pórtico rígido plano formado por dos barras ortogonales de igual longitud, las cuales cuentan con las dimensiones y condiciones de contorno indicadas en la figura 17. Las barras son de un material cuyo módulo de elasticidad vale 70.000 kp/cm2 (correspondiente a un tipo de aluminio). Tienen un área de 6,0 cm2 y una inercia de 2,0 cm4. Cada barra tiene una longitud de 120 cm y se dividen en 5 elementos iguales cada uno de los cuales presenta una longitud de 24 cm. Para la obtención de la trayectoria de equilibrio, en este caso se han llevado a cabo 900 incrementos de carga con los siguientes parámetros (𝛾 = 0,5, 𝜆0 = 2, 𝜃𝑑𝑒𝑠 = 5, tol = 0.001, máx. nº iter = 20). Figura 17: Geometría del modelo. Figura 18: Discretización del modelo en Mathematica. En este caso se ha considerado el sentido positivo del eje y, hacia abajo. 20 40 60 80 100 120 20 40 60 80 100 120 37 Se aplica en el punto indicado en la figura una carga vertical variable P, la cual se grafica frente a la posición vertical del punto de aplicación de esta, obteniendo la siguiente curva de equilibrio. Figura 19: Gráfica y - 𝜆 según Kotronis y Collin. Posición del punto de aplicación de la carga - Carga aplicada. Figura 20: Gráfica y – 𝜆 según nuestro código. Posición del punto de aplicación de la carga - Carga aplicada. En la figura 20 se muestra la trayectoria de equilibrio obtenida con nuestro código: • Los puntos (y, 𝜆) obtenidos en el código, son puntos de equilibrio y la línea que los une, es la “trayectoria de equilibrio de la estructura”. • Inicialmente nos encontramos en el punto (120, 0). La estructura está descargada y las barras se encuentran en un estado libre de esfuerzos. • Conforme aumentamos el valor de la carga vertical P, el punto de aplicación de esta va descendiendo y el problema, que inicialmente tiene una rama lineal, se hace cada vez más no lineal observándose claramente en la figura la no linealidad entre cargas aplicadas y descensos. • La trayectoria presenta dos puntos límites y dos puntos de retroceso o puntos límite de desplazamiento. 60 80 100 120 -1 -0.5 0.5 1 1.5 38 • A medida que la carga va aumentando, la estructura va perdiendo rigidez y cuando la carga es aproximadamente de unos 1,89 KN (𝜆 =1,89), la estructura pierde toda su rigidez produciéndose el pandeo global de la misma. Notas: En este caso no hemos tenido que perturbar para obtener las sucesivas posiciones de equilibrio. El código no ha dado problemas de oscilaciones, raíces imaginarias... ni otros problemas numéricos. El pandeo global del sistema nos proporciona una trayectoria suave, sin presencia de puntos de bifurcación. Alcanzado el punto de pandeo global (71, 1,89), para obtener los sucesivos estados de equilibrio tenemos inicialmente que reducir la carga llegando al punto (58,5, 1,21). Este punto es un punto de retroceso dado que al atravesarlo se produce una disminución del desplazamiento acompañado de una disminución de carga. El otro punto de retroceso es el (68, -0,50). El otro punto límite es el (62, -1,00). A partir de este punto y en adelante, la rigidez tangente pasa a ser positiva y para producir aumentos en el descenso del punto de aplicación de la carga, necesitamos aumentar progresivamente la carga. • Dado que la propia deformación de la estructura excita modos de pandeo globales, no aparecen puntos de bifurcación que pudieran conducir a modos de pandeo de carácter local. 6.4. Arco rebajado con carga excéntrica. El cuarto problema que se estudia se trata de un arco rebajado con la geometría indicada en la figura 21. Los resultados obtenidos con nuestro código se comparan con los obtenidos para el mismo en el artículo de P. Kotronis y F. Collin Implementation of path following techniques into the finite element code Lagamine. Figura 21: geometría del modelo. El arco es de un material cuyo módulo de elasticidad vale 20.000 kp/cm2. Tiene un área transversal de 100,0 cm2 y una inercia de 10.000,0 cm4. El arco cubre una luz de 10 m y se encuentra apoyado en ambos extremos. Presenta una altura de 0,5 m en su punto medio y se carga de forma excéntrica tal y como se indica en la figura 21. Para el trazado del arco, se dibuja en CAD un arco de círculo que pasa por los tres puntos (los extremos y el centro) y se trazan líneas verticales paralelas y equidistantes 500 mm, de modo que se puedan obtener las 39 cotas de los diferentes puntos del arco. Dado que la posición de la carga al centro del arco vale 200 mm, se tiene en cuenta para la determinación de las coordenadas de estos puntos. El arco, se divide en 20 elementos barra y su geometría en Mathematica es la que se muestra en la siguiente figura. Figura 22: Discretización del modelo en Mathematica. Para la obtención de la trayectoria de equilibrio, en este caso se han llevado a cabo 230 incrementos con los siguientes parámetros (𝛾 = 0,5, 𝜆0 = 10, 𝜃𝑑𝑒𝑠 = 5, tol = 0.001, máx. nº iter = 15). Se aplica en el punto indicado de la figura una carga vertical variable P, la cual se grafica frente al descenso vertical del punto de aplicación de esta, obteniendo la siguiente curva de equilibrio. Figura 23: Gráfica 𝛿𝑦 - 𝜆 según Kotronis y Collin. Descenso del punto de aplicación de la carga - Carga aplicada. Según nuestro código: 200 400 600 800 1000 10 20 30 40 50 40 Figura 24: Gráfica 𝛿𝑦 - 𝜆 según nuestro código. Descenso del punto de aplicación de la carga - Carga aplicada. Figura 25: Gráfica 𝛿𝑦 - 𝜆 según nuestro código. Descenso del punto de aplicación de la carga - Carga aplicada. Según vemos en este ejemplo, se presentan cuatro puntos límites y dos puntos de retroceso o puntos límites de desplazamiento lo que da lugar a una trayectoria de equilibrio bastante enrevesada. Los resultados obtenidos con nuestro código, son muy parecidos a los aportados por P. Kotronis y por F. Collin lo que demuestra que nuestro código está correcto. En este último ejemplo, el código tampoco ha dado problemas a la hora de encontrar las sucesivas posiciones de equilibrio, lo que parece indicativo de que el problema planteado no presenta puntos de bifurcación del equilibrio que es lo que hace que nuestro código oscile. Si en lugar de aplicar la carga en el punto dado, la aplicamos en el punto central del arco, en este caso, si no aplicamos carga perturbadora, el código programado presenta problemas numéricos en las proximidades del punto de pandeo (donde la pendiente de la gráfica se hace horizontal). Aparecen raíces imaginarias, mensajes de matriz mal condicionada... La siguiente gráfica muestra el resultado obtenido para este caso. 200 400 600 800 1000 -500 500 1000 1500 200 400 600 800 1000 -500 500 1000 41 Figura 26: Gráfica 𝛿𝑦 - 𝜆 según nuestro código para el caso de carga centrada. Descenso del punto de aplicación de la carga - Carga aplicada. Como puede verse en la figura 26, en el punto de pandeo parece como si la trayectoria de equilibrio retornase hacia atrás y buscase el camino de bifurcación más o menos lineal. Al final, por más que aumentemos el nº máximo de iteraciones, el código oscila y no es capaz de seguir más allá del último punto encontrado. Si además de la carga anterior introducimos un par de fuerzas perturbadoras, el resultado es el que se muestra en la figura 28 (la perturbación se introduce según se indica en la figura 27 como un par de fuerzas de igual módulo y sentido contrario, las cuales actúan en puntos simétricos respecto del punto central del arco y de forma que el resultado de estas equivale a un par aplicado en el nodo central el cual obliga al arco a adoptar una geometría deformada no simétrica). Figura 27: Introducción de las fuerzas perturbadoras. 200 400 600 800 -250 250 500 750 1000 1250 48 Figura 35: Geometría del modelo en SAP-2000. El pórtico lo cargamos según la figura 34, con un valor para 𝑷 de 𝑷=−20.000 kg Tras un análisis plano en SAP, obtenemos que la carga de pandeo del pórtico vale: λcrit= 8,89 Es decir, tendríamos que multiplicar las cargas por 8,89 para que el pórtico colapsase por pandeo. La diferencia con el resultado obtenido por Argüelles puede deberse entre otras causas a que SAP emplea para su formulación elementos finitos. Figura 36: Carga crítica de pandeo con SAP-2000. 49 Con nuestro código, una vez introducidos el tipo de material, geometría, condiciones de contorno, propiedades de las secciones, cargas…, el resultado es el siguiente. Figura 37: Carga crítica con nuestro código. Desplazamiento vertical de la esquina superior izquierda (punto de control) - Carga aplicada. La figura 37 muestra el valor de 𝜆 frente al desplazamiento vertical del punto superior izquierdo del pórtico. Tal y como podemos ver en la figura, el valor obtenido de 𝜆 = 9,03 es algo superior al obtenido con SAP y algo inferior al obtenido para este caso por Argüelles. Esto puede deberse al hecho de que SAP divide internamente las diferentes barras a la vez que la formulación de los elementos es distinta al ser, como ya hemos comentado previamente, con elementos finitos. En la figura 37 podemos ver que la trayectoria de equilibrio de la estructura ha de hacerse horizontal con tangente en torno al valor de 𝜆𝑐𝑟𝑖𝑡 = 9,1 que es el valor dado por Argüelles. Lógicamente en el cálculo de estructuras, una de las limitaciones existentes es la deformación de estas. Las estructuras tienen que proveer cierta sensación de bienestar y esto se consigue limitando sus deformaciones a valores normativos. En este último caso vemos que para cuando la carga se aproxima a la carga de pandeo, la deformación vertical del punto de control (esquina superior izquierda) es de 35 cm lo cual, desde el punto de vista estructural es del todo inasumible, advirtiendo además que para esta deformación, la falla de la estructura puede llegar con anterioridad como consecuencia del estado tensional sobre las diferentes piezas (nuestro análisis es elástico lineal). Así, en este caso y con un nº máximo de incrementos de carga de 300, nuestro código no ha encontrado todavía el punto en el que la estructura pierde su rigidez y por lo tanto el valor aportado es un valor aproximado que podría mejorarse aumentando el nº máximo de incrementos pero que en este caso y por la deformación observada, carece de sentido. Con el valor de 𝜆 obtenido, podemos determinar las longitudes de pandeo de los pilares a través de la expresión 𝑷𝑐𝑟𝑖𝑡=𝜆𝑐𝑟𝑖𝑡 𝑵 que junto con 𝑷𝑐𝑟𝑖𝑡=𝜋2∗E∗I 𝑙𝑘2 podemos deducir el valor de la longitud de pandeo de los pilares 𝑙𝑘. En la ecuación anterior 𝑵 es el axil que obra en cada uno de los pilares. -35 -30 -25 -20 -15 -10 -5 5.5 6.5 7 7.5 8 8.5 9 50 𝑙𝑘=�𝜋2∗E∗I 𝑷𝑐𝑟𝑖𝑡 =�𝜋2∗E∗I 𝜆𝑐𝑟𝑖𝑡 𝑵 ⇒ 𝑙𝑘=β∗𝑙 y por tanto β = 𝑙𝑘 𝑙 En este caso, el β del pilar de la izquierda vale βizq = 1,444 y para el pilar de la derecha βdch = 1,338 NOTA: estos resultados están obtenidos en base a 𝜆𝑐𝑟𝑖𝑡= 9,1 , 𝑵𝑖𝑧𝑞=24.875 𝑘𝑝 , 𝑵𝑑𝑐ℎ=20.126 𝑘𝑝 , 𝑙𝑖𝑧𝑞=500 𝑐𝑚 y 𝑙𝑑𝑐ℎ=600 𝑐𝑚. *********************************************************************************************************** 51 7. FUNCIONAMIENTO DEL PROGRAMA. El código funciona de la siguiente manera. En primer lugar se define el material del que están formadas las barras de la estructura. Para ello se emplea la variable me (módulo de elasticidad), la cual es una constante del problema. En este particular al tratarse de acero el módulo de elasticidad que empleamos es 𝑚𝑒= 2.141.404 kg/cm2 En segundo lugar definimos la geometría de la estructura a analizar. Para ello introducimos las coordenadas de los puntos de la estructura (en el siguiente ejemplo las jácenas de la estructura se modelan con un único elemento barra, los pilares inferiores con tres y los superiores con dos). Así, la estructura tiene un total de 13 barras (b=13) y 13 nudos (n=13). Figura 38: Geometría del pórtico a estudiar. Figura 39: Geometría pórtico en Mathematica. Figura 40: Coordenadas (en cm) de los puntos de la estructura. Los puntos de introducen mediante la definición de las siguientes variables 𝑥0𝑖 e 𝑦0𝑖 con 𝑖= 1, 2, 3 … . 𝑛 200 400 600 800 100 200 300 400 500 600 700 52 Como ejemplo, en este caso tenemos 𝑥01= 0 ; 𝑦01= 0 𝑥03= 0 ; 𝑦03=266 Figura 41: Numeración de los puntos de la estructura. Para definir los extremos de las b barras, se definen las variables 𝒃𝒂𝒓𝒓𝒂𝑖 como un vector de dos componentes cuya primera componente es la numeración del punto inicial de la barra y la segunda componente es la numeración del punto final. Figura 42: Numeración de las barras de la estructura. 𝒃𝒂𝒓𝒓𝒂𝑖={a, b} con 𝑖= 1, 2, 3 … . 𝑏 En este ejemplo 𝒃𝒂𝒓𝒓𝒂1={1, 2} 𝒃𝒂𝒓𝒓𝒂7={8, 9} 53 Se definen las condiciones de contorno mediante la variable 𝑵𝒄. En este caso los puntos unidos a cimentación son el 1 y el 10. El punto 1 está empotrado y el 10 está articulado. Al tratarse de un estudio plano, cada punto (nodo) tiene tres posibles componentes de desplazamiento que son 𝛿𝑥, 𝛿𝑦 y el giro 𝜃 . Para este ejemplo en el que tenemos 13 puntos, esto hace un total de 39 componentes de desplazamiento. El grado de indeterminación cinemática de la estructura son precisamente esas 39 componentes menos las componentes de desplazamiento impedidas. En este ejemplo, las componentes de desplazamiento impedidas son los movimientos 𝛿𝑥, 𝛿𝑦 y el giro 𝜃 del punto 1 y los movimientos 𝛿𝑥, 𝛿𝑦 del punto 10. La variable 𝑵𝒄 se define como un vector de dimensión igual al nº de grados de libertad coartados en la estructura y que nos permite conocer las restricciones de la misma. En este ejemplo en particular esta variable es un vector de 5 componentes y resulta ser 𝑵𝒄={1, 2, 3, 28, 29} El vector 𝑵𝒄 así definido nos indica sin posibilidad de error las ligaduras de la estructura. Por otro lado, la variable 𝑫𝑵𝒄 se define como un vector de dimensión igual al grado de indeterminación cinemática de la estructura y nos indica los grados de libertad en la misma. En este ejemplo, este vector tiene 16 componentes y resulta ser 𝑫𝑵𝒄={4, 5, 6, 7, 8, 9, 10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,30,31,32,33,34,35,36,37,38,39} A continuación se definen las propiedades de sección e inercia de cada barra, asignando valores numéricos (****) a las variables siguientes 𝛺𝑖 y 𝑚𝑖 respectivamente. En nuestro ejemplo: 𝛺𝑖=149 cm2 i = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,13 𝛺11=𝛺12=198 cm2 𝑚i=25670 cm4 i = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,13 𝑚11=𝑚12=57680 cm4 Para la introducción de cargas se define un vector de cargas básicas 𝑷𝒂 de tantas componentes como grados de libertad tiene la estructura. Por cada punto no coaccionado tenemos tres posibles componentes de fuerza que son Fx, Fy y M. Así, el vector de cargas tiene la dimensión del vector 𝑫𝑵𝒄. 54 Figura 43: Estado de cargas básicas (Pa) sobre la estructura. Con las cargas apuntadas en la figura 43, 𝑷𝒂 toma la forma siguiente: 𝑷𝒂={0, 0, 0, 0, 0, 0, 0.01,−1, 0, 0, 0, 0, 0.01,−1, 0, 0, −1, 0, 0, 0, 0, 0, −1, 0, 0, 0, 0, 0, 0, 0, 0, 0, −1, 0 } Notemos que las cargas básicas están referidas a los ejes globales. El nº en rojo del vector 𝑷𝒂, es el momento externo aplicado en el punto 10, que es nulo. Con los datos introducidos, se resuelve linealmente la estructura obteniendo un vector de desplazamientos 𝑼 que en nuestro caso toma el valor que se indica la tabla de resultados siguiente. Este vector de desplazamientos o movimientos de los puntos de la estructura, también está referido a los ejes globales. Para ello se ha tenido que montar la matriz de rigidez global de la estructura 𝑲 y resolver el sistema de ecuaciones 𝑲 𝑼=𝑷𝒂 55 Tabla 1. Según vemos en el resultado, para este estado básico de carga, el movimiento horizontal del punto 6 (punto de control) vale 0,0000555412 cm. Para resolver el problema no lineal, se definen una serie de parámetros los cuales van a gobernar el análisis. Estos parámetros son los siguientes: • Exponente (𝛾) para la actualización de la longitud del arco (Ec53) 𝛾= 0,5 • Nº de iteraciones deseadas (𝜃𝑑𝑒𝑠) en cada salto de carga, para la actualización de la longitud del arco (Ec53) 𝜃𝑑𝑒𝑠= 5 • Valor inicial del parámetro de carga (𝜆), para el cálculo del vector de carga inicial 𝜆=100 • Valor de la variable (𝑡𝑜𝑙), la cual indica cuando se alcanza la convergencia 56 𝑡𝑜𝑙= 0,001 • Valor de la variable (𝑖𝑡𝑒𝑟𝑎1), la cual indica el nº máximo de iteraciones permitido dentro de un salto de carga 𝑖𝑡𝑒𝑟𝑎1=40 • Valor de la variable (𝜙𝑚𝑎𝑥), la cual indica el nº máximo de saltos de carga que van a realizarse en el análisis. 𝜙𝑚𝑎𝑥=120 • Se define la variable (𝐶𝑜𝑛𝑡𝑟𝑜𝑙). Con esta variable se indica el punto de la estructura sobre el que se quiere hacer el seguimiento de uno de sus grados de libertad. En este caso la trayectoria de equilibrio que queremos seguir es la del desplazamiento horizontal del punto 6 de la estructura. Así: 𝐶𝑜𝑛𝑡𝑟𝑜𝑙= 6 • Se define la variable (Δ𝑼𝜙). Esta variable es un vector nulo de tantas componentes como grados de libertad tiene la estructura (34) y representa el desplazamiento acumulado al inicio del primer salto de carga. Este vector se va actualizando durante el proceso iterativo dentro de cada salto de carga. Δ𝑼𝜙= {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0} • Se define la variable (cont2). Es simplemente un contador que nos marca el nº de iteraciones que llevamos dentro de un salto de carga. Inicialmente, al principio de cada salto de carga 𝑐𝑜𝑛𝑡2= 0 Asignados los valores numéricos a los parámetros o variables anteriores, se procede a calcular en primer lugar la carga inicial que actúa sobre la estructura. Esto se hace definiendo el vector 𝑷𝑷𝒂𝒂 que se calcula como 𝑷𝑷𝒂𝒂=λ∗𝑷𝒂 El parámetro de carga λ, va variando a medida que avanza el análisis. En segundo lugar, se calcula de forma lineal el vector de desplazamientos de la estructura para ese estado inicial de carga. Se calcula la primera predicción para el vector de desplazamientos dentro del primer salto de carga Δ𝑼10=𝜆∗𝑼 A continuación, el programa evalúa la longitud inicial de arco haciendo uso de la ecuación Ec48 Δ𝑙 = 𝜆∗�(Δ𝑼10)T Δ𝑼10 Se define la variable Trayectoria. Esta variable nos va a permitir acumular los puntos (𝛿𝑥, 𝜆) en el equilibrio al final de cada salto de carga y así poder luego trazar la gráfica del camino de equilibrio. Inicialmente y en este caso particular definimos esta variable como 𝑇𝑟𝑎𝑦𝑒𝑐𝑡𝑜𝑟𝑖𝑎 = {{0, 0}} 57 que representa el punto inicial de la trayectoria que deseamos obtener. Inicialmente para 𝜆= 0 (con la estructura descargada), el desplazamiento horizontal del punto 6 es nulo (𝛿𝑥= 0). Conocido el vector de desplazamientos dado por la predicción Δ𝑼10, para iniciar el análisis no lineal se procede a recalcular la geometría de la estructura a través de la definición de las variables siguientes 𝑥𝑖 = 𝑥0𝑖+ 𝛿𝑥𝑖 𝑖= 1, 2, … , 𝑛 e 𝑦𝑖 = 𝑦0𝑖+ 𝛿𝑦𝑖 𝑖= 1, 2, … , 𝑛 Con esta nueva geometría, se evalúan las nuevas longitudes (𝐿𝑗) de las barras y los nuevos ángulos (𝛽𝑗) que cada una de las barras forma con el eje x global. 𝐿𝑗 = �(𝑦𝑏−𝑦𝑎)𝑗2+ (𝑥𝑏−𝑥𝑎)𝑗2 𝑗= 1, 2, … , 𝑏 y 𝛽𝑗 = 𝐴𝑟𝑐𝑡𝑎𝑔 (𝑦𝑏−𝑦𝑎)𝑗 (𝑥𝑏−𝑥𝑎)𝑗 𝑗= 1, 2, … , 𝑏 Calculadas las nuevas longitudes y ángulos, pasamos a calcular con la Ec78 las matrices 𝑩𝑗 de cada una de las barras. De este modo, a partir de la ecuación Ec76, se puede evaluar para cada una de las barras el vector de movimientos δ𝒑 �𝑗=�δu� δ𝜗1 δ𝜗2�𝑗 =𝑩𝑗 δ𝑼𝑗 donde δ𝑼𝑗= {𝛿𝑥𝑗𝑎 , 𝛿𝑦𝑗𝑎 , 𝜃𝑗𝑎 , 𝛿𝑥𝑗𝑏 , 𝛿𝑦𝑗𝑏 , 𝜃𝑗𝑏} con a y b la numeración de los nudos extremos de la barra. Con el vector de movimientos de cada barra y a través de la ecuación Ec84, calculamos para cada barra el vector de fuerzas internas 𝑸 �𝑗. 𝑸 �𝑗=�N M1 M2�𝑗=𝑪 𝒑 �𝑗 Estas fuerzas internas están referidas a los ejes locales de la barra. Haciendo uso de la ecuación Ec82, obtenemos los vectores de fuerzas internas (𝑸) de cada una de las barras pero, referidas estas a las coordenadas globales. 𝑸𝑗=�N∗ M1∗ M2∗�𝑗=𝑩𝑗𝑇 𝑸 �𝑗 Conocidos estos vectores, sumando adecuadamente las componentes de todas las barras que concurren en un nudo, podemos montar el vector de fuerzas internas de toda la estructura 𝑽𝑭𝒏𝒐𝒅, vector que tendrá la