Source: https://www.scribd.com/document/81053219/7-CAPITULO-6
Timestamp: 2016-10-23 07:17:13+00:00

Document:
BrowseBrowseInterestsBiography & MemoirBusiness & LeadershipFiction & LiteraturePolitics & EconomyHealth & WellnessSociety & CultureHappiness & Self-HelpMystery, Thriller & CrimeHistoryYoung AdultBrowse byBooksAudiobooksComicsSheet MusicBrowse allUploadSign inJoinBooksAudiobooksComicsSheet MusicDISEÑO DE UN SILO CONFORME AL EUROCÓDIGOCAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 72 CAPÍTULO 6 6. MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA Ansys puede usar diferentes métodos matemáticos para la resolución de las ecuaciones numéricas de equilibrio necesarias para resolver los problemas de elementos finitos planteados. En nuestro caso estamos interesados en la resolución de las ecuaciones de pandeo que se han planteado en el capítulo 4. Ansys dispone de métodos numéricos distintos según se esté resolviendo el problema lineal (Pandeo de Autovalores, apropiado para obtener una primera carga crítica de manera más rápida) o el problema no lineal (la carga crítica bajará al introducir no-
linealidades en el problema). En este capítulo se presentan dichos métodos, introduciendo el cálculo de pandeo de autovalores primero y luego entrando más en detalle en el cálculo no lineal por ser este último más complejo y tener más problemas de convergencia. Los métodos numéricos de DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 73 cálculo no lineal se basan sobre todo en el método de Newton-Raphson y en variaciones del mismo. Se verá cómo el programa los lleva a cabo, con las diferentes órdenes asociadas y los parámetros de que dispone para controlar estos métodos. Por último se hará especial hincapié en el método del “arc-length” que es el más apropiado para el caso que se estudia en el presente proyecto. 6.1. Pandeo de Autovalores El análisis de pandeo de autovalores predice la carga crítica de pandeo de una estructura ideal lineal y elástica (es decir el punto de bifurcación en la curva carga-
desplazamiento de la estructura). Este método se corresponde con el método clásico del análisis de pandeo elástico que se puede encontrar en los libros; por ejemplo, al hacer el análisis de pandeo de autovalores de un cilindro hueco de espesor de pared constante sometido a compresión pura los resultados se corresponden con la solución clásica de Euler (4.36). Sin embargo, las imperfecciones y otras no linealidades hacen que la mayoría de las estructuras que podemos encontrar en la vida real no alcancen esta carga crítica clásica de pandeo, sino que la carga real de pandeo será menor de la prevista por el análisis elástico. Por esto este análisis sólo debe hacerse como paso previo de un análisis no lineal o como una primera aproximación a la carga crítica de una estructura ya que proporciona valores de la carga crítica mayores a los reales y que son, por tanto, no seguros. Ansys realiza el análisis de pandeo de autovalores usando un modelo linealizado de estabilidad elástica. El punto de bifurcación se corresponde, en este caso, con el crecimiento no acotado de un nuevo patrón de deformación. Una estructura con una curva de comportamiento carga-desplazamiento como la mostrada en la figura 6.1.a puede ser analizada bien por este método, mientras que si muestra el comportamiento mostrado en la figura 6.1.b no se analiza bien. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 74 Figura 6.1: “Validez de los resultados de un análisis de pandeo de autovalores” En cualquier cálculo de pandeo de una estructura se debe tener en cuenta el fenómeno de “Stress stiffening”. Este es el fenómeno por el cual las fuerzas de membrana (axiles) influyen en la rigidez a deflexión lateral; por ejemplo, una cuerda de guitarra se hace más rígida al tensarla (su frecuencia natural de vibración es más alta), pero una carga axial de compresión reduce la rigidez a flexión haciendo que la estructura se haga inestable al alcanzar cierta carga: la carga crítica de pandeo. Este fenómeno de 2º orden conocido como “Stress stiffening” se tiene en cuenta al realizar el análisis de pandeo de autovalores de Ansys incluyendo un cálculo tensional previo para incluir los efectos de las tensiones de membrana en la matriz de rigidez geométrica (“stress stiffness matrix”) que se suma a la obtenida con la teoría de primer orden. ss
= + (6.1) La matriz de rigidez geométrica depende de: - Geometría del elemento. - Desplazamientos. - Valores de los esfuerzos de membrana. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 75 6.1.1. Proceso del cálculo de Pandeo de Autovalores Cuando se analiza el pandeo lineal en Ansys, el problema de pandeo es formulado como un problema de autovalores siguiendo el siguiente proceso: 1. Se carga la estructura hasta un nivel arbitrario de referencia de cargas exteriores,
F . Se hace un análisis estático lineal para calcular los esfuerzos de membrana en los elementos ref ref
K D F ⋅ = (6.2) De la ecuación (6.2) calculamos los desplazamientos ref
D y después las tensiones y esfuerzos correspondientes al nivel de carga ref
F . Se calcula la matriz de rigidez geométrica para ese nivel de carga ,ref
, para ello se hace el análisis estático previo antes mencionado, para el cual se introduce en el programa el código expuesto en la figura 6.2: Figura 6.2: “Código de Ansys para el cálculo previo de la matriz de rigidez geométrica” 2. Para otro nivel de cargas exteriores se tiene: , ref ref
F F K K
λ λ = ⋅ ⇒ = ⋅ (6.3) /SOLU ANTYPE,STATIC PSTRES,ON ! Aplicar cargas
SOLVE FINISH DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 76 Al ser un análisis lineal, se supone que al multiplicar las cargas por λ también se multiplican las tensiones por el mismo factor, pero sin cambiar la distribución de tensiones. Los desplazamientos para el nuevo valor de carga son ( )
,ref ref
K K D F
λ λ + ⋅ ⋅ = ⋅ (6.4) 3. Si se aplica una perturbación sobre la carga δF, se produce un movimiento δD. El pandeo (bifurcación del equilibrio) sucede cuando es posible un δD 0 con δF=0: ( ) ( )
λ δ λ + ⋅ ⋅ + = ⋅ (6.5) Restando (6.5) y (6.4) se obtiene el problema de autovalores al que se alude al realizar el análisis: ( )
λ δ + ⋅ ⋅ = (6.6) El menor de los autovalores nos proporciona la carga crítica de pandeo cr
λ . Su autovector es el modo de pandeo. El problema de autovalores de la ecuación (6.6) se resuelve mediante alguno de los procedimientos explicados en 6.1.2. Los autovectores están normalizados de manera que el componente más grande es 1,0. Por tanto, las tensiones (cuando se calculan como resultado del problema) deben ser interpretadas como una distribución relativa de tensiones y no directamente como el estado tensional. Una muestra del código que se introduce en Ansys para resolver el problema de pandeo por autovalores se muestra en la figura 6.3: DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 77 Figura 6.3: “Código de Ansys para realizar el cálculo del problema de Pandeo por Autovalores” Si el primer autovalor es negativo (lo cual indica que el pandeo se produce para cargas aplicadas en el sentido opuesto), Ansys da un warning y el programa termina. En el método del subespacio, se muestra un error “Number of stress-stiffness DOF is less than requested modes”. Para evitar esta dificultad, aplicar un desplazamiento inicial de la búsqueda (con el comando BUCOPT) hasta un valor próximo al autovalor esperado. 6.1.2. Problema de extracción de autovalores Ansys dispone de varios métodos para la extracción de autovalores que se enumeran en la tabla 6.1. Los que se usan en el análisis del pandeo de autovalores están recuadrados en rojo en dicha tabla. /SOLU !Debe haberse resuelto previamente un problema estático con PSTRES,ON ANTYPE,BUCKLE ! Análisis de pandeo !por autovalores. BUCOPT,LANB,10 ! Define el método de !cálculo de autovalores y ! el nº de modos a calcular. MXPAND,1 ! Expande la forma !de los 10 modos calculados. SOLVE FINISH DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 78 Tabla 6.1: “Métodos de extracción de Autovalores en Ansys 8.0” ANSYS ofrece varios métodos para la obtención de autovalores y autovectores en un sistema de N g.d.l.: • REDUCED – Emplea matrices reducidas, asociadas a una serie de g.d.l. maestros que el usuario (o el propio ANSYS) debe seleccionar previamente. El proceso de solución es más rápido que con las opciones SUBSPACE pero menos exacto, ya que las matrices reducidas sólo permiten obtener una solución aproximada. • SUBSPACE – Permiten obtener un determinado número M (reducido: M<<N) de autovalores y autovectores. No resulta necesario definir g.d.l. maestros. • BLOCK LANCZOS – Válido en los mismos casos que Subspace, pero tiene una convergencia más rápida. Aplicable en problemas de autovalores con matrices simétricas grandes. Es el más eficiente en el caso de que haya autovalores próximos DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 79 (es el caso del pandeo por autovalores, por tanto es el método que usaremos). Como se ha mencionado anteriormente, los métodos de extracción de ecuaciones más frecuentemente usados por Ansys en la extracción de autovalores del problema de Pandeo de autovalores son el método del Subespacio y el de Block Lanczos y por ello los veremos con más detenimiento. Método del Subespacio El método iterativo del subespacio (que se activa con el comando BUCOPT,SUBSP) se describe con detalle en Bathe [5]. Las mejoras de este método sugeridas por Wilson e Itoh [16] también se han incluido en Ansys. El algoritmo básico consiste en los siguientes pasos: 1. Se define el pivote inicial s : • En un análisis de pandeo (ANTYPE,BUCKLE), s SHIFT = en el comando BUCOPT (su valor por defecto es 0,0). 2. Inicializar los vectores iniciales [
X ]. 3. Triangularizar la matriz pivotada [ ] [ ]
K K s M ( = +
(6.7) donde: - [ ] K es la matriz de rigidez ensamblada. - [ ] M es la matriz de masas (o de rigidez geométrica) ensamblada. Se hace una comprobación de “secuencia de Sturm” si este es un punto de pivote distinto del inicial y se ha requerido (mediante los comandos Strmck=ALL, que se usa por defecto, o PART en el comando SUBOPT). 4. Para cada iteración n del subespacio (desde 1 hasta M
N ), se repiten los pasos del 5 al 14: DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 80 donde M
N es el máximo número de iteraciones del subespacio (se pasa como entrada en NUMSSI del comando SUBOPT). 5. Se forma [ ] [ ][ ]
= y se escala [ ] F por { }
son los autovalores previamente estimados. 6. Se resuelve para n
K X F ( ( =
(6.8) Estas ecuaciones se resuelven usando el frontal direct equation solver de Ansys (EQSLV,FRONT) o el iterative PCG solver (EQSLV,PCG). 7. Se escalan los vectores n
por ( ) { }
s λ λ
− 8. Se ortogonalizan los vectores a los vectores previamente convergidos (Ortogonalización de Gram-Schmidt). 9. Se definen las matrices de subespacio K (
y M (
K X K X ( ( ( =
(6.9) [ ]
M X M X ( ( ( =
(6.10) 10. Se ajusta para el pivote *
K K s M ( ( ( = +
11. Computar los autovalores y vectores del subespacio usando la iteración generalizada de Jacobi: [ ] [ ]{ }
K Q M Q λ ( ( =
(6.11) donde: - [ ] Q son los autovectores del subespacio. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 81 - { }
λ son los autovalores actualizados. 12. Actualizar la aproximación a los autovectores: [ ] [ ]
X X Q ( =
(6.12) 13. Si se encuentra algún modo redundante o negativo, se eliminan y se crea un nuevo vector aleatorio. 14. Se comprueba la convergencia (con la tolerancia definida por el analista en Ansys): • Si todos los modos requeridos convergen, ir al paso 15. • Si se requiere un nuevo pivote (como se indica más abajo), ir al paso 3. • Ir a la nueva iteración, paso 4. 15. Se realiza una comprobación de la secuencia de Sturm si se pide (Strmck=ALL, por defecto en el comando SUBOPT). Los pasos del 5 al 12 sólo se realizan en los vectores no convergidos: Una vez que un autovalor ha convergido, el autovector asociado no se vuelve a iterar. El procedimiento de ortogonalización de Gram-Schmidt (paso 8) asegura que los autovectores que no han convergido permanezcan ortogonales a los vectores convergidos sobre los que no se itera. Método de Block Lanczos El método de extracción de autovalores de Block Lanczos (es el algoritmo de Lanczos por bloques al que se accede con el comando BUCOPT,LANB) está disponible para problemas de gran tamaño de autovalores con simetría. Típicamente, este solucionador es aplicable al tipo de problemas resolubles con el método de autovalores del subespacio, pero consigue una tasa de convergencia más rápida. El algoritmo de pivote por bloques de Lanczos , como se detalla en Grimes [17], es la base teórica de este extractor de autovalores. El método usado por el análisis modal emplea una estrategia de cambio de pivote automática, combinada con comprobaciones DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 82 de las secuencias de Sturm, para extraer el número de autovalores requerido. La comprobación de las secuencias de Sturm, también asegura que se encuentra el número de autovalores requeridos más allá de la frecuencia proporcionada por el analista como pivote inicial (FREQE en el comando MODOPT) se encuentran sin obviar ningún modo. El algoritmo de Block Lanczos es una variación del algoritmo de Lanczos clásico, donde las recursiones de Lanczos se efectúan usando un bloque de vectores, en lugar de un solo vector. Más detalles sobre el método clásico de Lanczos se pueden encontrar en Rajakumar y Rogers [18]. El uso del método de Lanczos por bloques (o método Block Lanczos) para resolver grandes modelos (100.000 gdl, por ejemplo) con muchas restricciones puede requerir una cantidad de memoria del ordenador bastante significativa. Por esta razón, Ansys utiliza los Multiplicadores de Lagrange para tratar las ecuaciones de restricción en el extractor de Lanczos por bloques en lugar de eliminar dichas ecuaciones explícitamente antes de escribir las matrices completas en el archivo de cálculo. Para más detalles sobre la formulación de la teoría de los Multiplicadores de Lagrange ver Cook [19]. Comparación: Método de Block Lanczos vs Método del Subespacio El Método de extracción de autovalores de Block Lanczos es el que usa Ansys por defecto en el cálculo de pandeo por autovalores. Este método es tan preciso como el Método del Subespacio y además es más rápido. El método de Block Lanczos es especialmente potente cuando busca frecuencias de autovalores en partes determinadas del espectro de autovalores de un sistema (en el problema de pandeo de autovalores normalmente los autovalores están próximos entre sí). La tasa de convergencia de las frecuencias de autovalores en las partes media y alta del espectro es casi la misma que cuando se extraen los modos más bajos, por tanto, cuando se usa una frecuencia dada como pivote (especificada en el comando FREQB) para extraer n modos más allá de los hallados con el pivote inicial (valor inicial de FREQB), el algoritmo extrae estos n modos adicionales prácticamente a la misma velocidad que extrajo los n modos más bajos. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 83 Por su parte, el método del Subespacio es muy preciso debido a que usa las matrices [ ] K y [ ] M completas, aunque precisamente esa sea la causa de que sea más lento que el método Reducido o el de Block Lanczos. Por ello, el método del Subespacio sólo se usa cuando sea necesaria una gran precisión o cuando no sea práctico seleccionar g.d.l maestros (como se indicó anteriormente el método Reducido usa estos g.d.l maestros). A modo de resumen se muestra la tabla 6.2 en la que se comparan esquemáticamente los distintos métodos de extracción de autovalores usados por Ansys y que se han comentado en esta sección. Tabla 6.2: “Comparación esquemática de varios extractores de autovalores de Ansys” 6.2. Pandeo No-lineal 6.2.1. Método de Newton-Raphson: Visión general El proceso de discretización de los elementos finitos produce una serie de ecuaciones simultáneas: [ ]{ } { }
K u F = (6.13) DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 84 donde: - [ ] K = Matriz de coeficientes - { } u = Vector de incógnitas GDL (grados de libertad) - { }
F = Vector de cargas aplicadas Si la matriz de coeficientes [ ] K es ella misma función de los grados de libertad que son incógnita (o de sus derivadas) entonces la ecuación (6.13) es una ecuación no lineal. El método de Newton-Raphson es un proceso iterativo para resolver las ecuaciones no lineales que puede ser escrito como (Bathe[5]): { } { } { }
T a nr
K u F F ( ∆ = −
(6.14) { } { } { }
= + ∆ (6.15) donde: - T
= Jacobiano de la matriz (matriz tangente) - i = Subíndice que representa la iteración de equilibrio actual - { }
F = Vector de las fuerzas restauradoras correspondiente a las cargas elementales internas Ambos T
F se evalúan basándose en los valores dados por { }
u . La parte derecha de la ecuación (6.14) es el residuo o vector de cargas desequilibradas; i.e., la cantidad en que el sistema esta fuera del equilibrio. En la figura 6.4 se describe gráficamente la iteración i-ésima del algoritmo de Newton-Raphson aplicado a un problema de un solo grado de libertad. En un análisis estructural, T
es la matriz de rigidez tangente, { }
u es el vector de desplazamientos y { }
F es el vector de fuerzas restauradoras calculado a partir de las tensiones elementales. En un análisis térmico, T
es la matriz de conductividad, { }
u es el vector de temperaturas y { }
F es el vector de cargas resistentes calculado a partir de los flujos de calor elementales. En un DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 85 análisis electromagnético, T
es la matriz de Dirichlet, { }
u es el vector de potenciales magnéticos y { }
F es el vector de cargas resistentes calculado a partir de los flujos magnéticos elementales. En un análisis transitorio, T
es la matriz de coeficientes efectiva y { }
F es el vector de cargas efectivas aplicadas que incluye los efectos de inercia y de amortiguamiento. Como se puede observar en las siguientes figuras, se requiere más de una iteración del método de Newton-Raphson para obtener una solución que converja. El algoritmo general es como sigue: 1. Se supone{ }
u . { }
u es normalmente la solución convergida del “time step” (paso de carga) previo. En el primer “time step”, { }
u = { } 0 . 2. Calcular la matriz tangente modificada T
y la carga restauradora { }
F a partir de la configuración de los { }
u . 3. Calcular { }
u ∆ mediante la ecuación (6.14). 4. Añadir { }
u ∆ a { }
u para obtener la nueva aproximación { }
(ecuación (6.15)). 5. Repetir los pasos desde el 2 al 4 hasta obtener la convergencia. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 86 Figura 6.4: "Solución de Newton-Raphson - Una iteración" En la figura 6.5, se muestra la solución de la siguiente iteración ( ) 1 i + a la del ejemplo de la figura 6.4. Las iteraciones subsecuentes se realizarían de una manera análoga. La solución obtenida al final del proceso iterativo correspondería al nivel de carga { }
F . La solución final convergida estaría en equilibrio, de forma que el vector de cargas restauradoras { }
F (computado a partir del estado actual de tensiones, flujos de calor…) igualaría al vector de cargas aplicadas { }
F (o al menos estaría dentro de cierta tolerancia). Ninguna de las soluciones intermedias estaría en equilibrio. Cada uno de estos pasos intermedios se denomina “iteración de equilibrio” (equilibrium iteration) en Ansys. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 87 Figura 6.5: "Solución de Newton-Raphson - Siguiente iteración" Si el análisis incluye no linealidades dependientes de la trayectoria (como por ejemplo plasticidad), entonces el proceso de solución requiere que ciertos pasos intermedios estén en equilibrio para poder seguir correctamente la trayectoria de carga. Esto se consigue, efectivamente, especificando un análisis incremental paso a paso; i.e., el vector final de carga { }
F se alcanza aplicando la carga en incrementos y aplicando las iteraciones de Newton-Raphson en cada paso: { } { } { }
n i n i n n i
(6.16) donde: - ,
= Matriz tangente para el paso n, y la iteración i. - { }
F = Vector de fuerzas totales aplicadas en el paso n. - { }
F = Vector de las fuerzas restauradoras para el paso n, y la iteración i. Cada paso se denomina “paso de carga” (load step) en Ansys. Además cada paso de carga puede subdividirse a su vez en “sub-pasos de carga” (substeps) en Ansys, al final de los cuales también se alcanza el equilibrio. Este proceso es el “Procedimiento Incremental de Newton-Raphson” y se muestra en la figura 6.6. El procedimiento de Newton-Raphson garantiza la convergencia, si y solo si la solución en cualquier iteración { }
u esta “cerca” de la solución exacta. Por tanto, incluso de no haber presente una no DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 88 linealidad dependiente de la trayectoria, la aproximación incremental (i.e., aplicando las cargas en incrementos mas pequeños) es a veces necesaria para obtener una solución que se corresponda con el nivel final de carga. Figura 6.6: “Procedimiento Incremental de Newton-Raphson” Cuando la matriz de rigidez se recalcula en cada iteración (como se indica en la ecuación (6.14) y en la ecuación (6.16)) el proceso se denomina “procedimiento de solución completa de Newton-Raphson” (Full Newton-Raphson solution procedure), en Ansys: ( NROPT,FULL o NROPT,UNSYM). Alternativamente, la matriz de rigidez podría ser recalculada con menor frecuencia usando el “Procedimiento Modificado de Newton-
Raphson”, en Ansys: (NROPT,MODI). Especialmente, para análisis estáticos o transitorios, se recalcularía sólo durante la primera o la segunda iteración de cada subpaso, respectivamente. El uso del “Procedimiento de Rigidez Inicial”, en Ansys: (NROPT,INIT), previene cualquier actualización o recálculo de la matriz de rigidez, como se muestra en la figura 6.7. Si existiese un elemento “multistatus” en el modelo, sin embargo, la matriz de rigidez sería recalculada en la iteración en la cual se cambiara de estado independientemente de la opción de Newton-Raphson que se haya usado. Los procedimientos modificado y de rigidez inicial de Newton-Raphson, convergen más lentamente que el procedimiento de solución completa de Newton-Raphson, pero requieren menos reformulaciones e inversiones matriciales. Unos pocos elementos forman una matriz de rigidez tangente aproximada por lo que las características de convergencia son en cierto modo diferentes. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 89 Figura 6.7: “Rigidez-Inicial Newton-Raphson” 6.2.2. Convergencia El proceso iterativo descrito en la sección anterior continúa hasta que se consigue la convergencia. Se asume que hay convergencia cuando: { }
R R ε < (Convergencia residual) (6.17) y/o: { }
i u ref
u u ε ∆ < (Convergencia incremental de los GDL) (6.18) donde { } R es el vector residual: { } { } { }
R F F = − (6.19) DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 90 que es el término de la derecha de la ecuación de Newton-Raphson(6.14). { }
u ∆ es el vector de incrementos de los GDL, R
ε y u
ε son tolerancias (TOLER en el comando CNVTOL de Ansys) y ref
u son valores de referencia (VALUE en el comando CNVTOL de Ansys). • es la norma de un vector; esto es, una medida escalar de la magnitud del vector (más adelante se definen los tres tipos de norma que usa Ansys). La convergencia, por tanto, se obtiene cuando la norma del residuo (desequilibrio) es menor que una tolerancia multiplicada por un valor de referencia y/o cuando el tamaño de los incrementos de los GDL es menor que una tolerancia multiplicada por un valor de referencia. Por defecto, el programa Ansys sólo comprueba la convergencia residual, la de la ecuación (6.17). El valor por defecto de las tolerancias es 0,001 (tanto para
ε como para
ε ). Hay tres normas disponibles entre las que elegir (NORM en el comando CNVTOL de Ansys): 1. Norma infinita: { } max
= 2. Norma L1: { }
R R = ∑ 3. Norma L2: { } ( )
R R = ∑ Para la convergencia incremental de los GDL, basta con sustituir u ∆ por R en las anteriores ecuaciones. La norma infinita es simplemente el máximo valor en el vector (máximo residuo o máximo incremento de GDL), la norma L1 es la suma de los valores absolutos de los términos, y la norma L2 es la raíz cuadrada de la suma de los cuadrados (SRSS) de los valores de los términos, también llamada norma Euclídea. Por defecto se usa la norma L2. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 91 El valor por defecto de referencia de residuo ref
R es { }
F . Para los GDL con condiciones de contorno de desplazamientos impuestos, los { }
F en esos GDL se usan para computar ref
R . Para GDL estructurales, si { }
F cae por debajo de 1,0, entonces ref
R toma 1,0 como valor. Esto ocurre a menudo en análisis del movimiento de cuerpos rígidos (por ej., rotación libre). Para GDL térmicos, si { }
F cae por debajo de 1,0E-6, entonces ref
R toma 1,0E-6 como valor. Para el resto de GDL, ref
R toma el valor 0,0. El valor por defecto de referencia de ref
u es { } u . El número máximo de ecuaciones de equilibrio permitidas, N en cada subpaso, se fija con el comando NEQIT. Si después de realizar NEQIT iteraciones de equilibrio sin alcanzar la convergencia, Ansys da un mensaje de error y termina el cálculo. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 92 6.2.3. Variaciones del método de Newton-Raphson y opciones que usa el programa Predictor La solución usada para el comienzo de cada paso de carga n { }
u es generalmente igual a la solución actual de GDL { }
. La matriz de rigidez tangente ,0 n
F se basan en esta configuración. La opción “predictor” de Ansys (comando PRED), extrapola la solución de GDL usando la historia previa para una mejor estimación de la siguiente solución. En análisis estáticos, la predicción se basa en usar los incrementos de los desplazamientos acumulados sobre el paso de carga previo, multiplicados como factor por el tamaño del paso de carga para hallar el { }
u : { } { } { }
,0 1 n n n
u u u β
= + ∆ (6.20) donde: - { }
u ∆ = Incremento del desplazamiento acumulado sobre el paso de carga previo. - n = Paso de carga actual. { } { }
(6.21) β se define como: 1
(6.22) donde: - n
t ∆ = Tamaño del paso de carga actual. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 93 - 1 n
∆ = Tamaño del paso de carga precedente. No se permite que β sea mayor que 5. En análisis transitorios, la predicción se basa en las velocidades y aceleraciones actuales usando las fórmulas de Newmark para GDL estructurales: { } { } { } { }
,0 1 1 1
u u u t u t α
= + ∆ + − ∆
! !! (6.23) donde: - { } { } { }
! !! = Desplazamientos, velocidades y aceleraciones actuales. - n
t ∆ = Tamaño de paso de carga actual. - α = Parámetro de Newmark (se introduce con el comando TINTP). Para sistemas térmicos, magnéticos y otros sistemas de primer orden, la predicción se basa en la fórmula trapezoidal: { } { } ( ){ }
u u u t θ
= + − ∆ ! (6.24) donde: - { }
= Temperaturas actuales (o potenciales magnéticos). - { }
! = Tasas actuales de estas cantidades. - θ = Parámetro de tiempo de integración trapezoidal (se introduce con el comando TINTP). Las subsecuentes iteraciones de equilibrio proporcionan los incrementos de los GDL { } u ∆ con respecto a los valores predichos de los mismos { }
u , de ahí que este sea un algoritmo predictor-corrector. Descenso Adaptativo El descenso adaptativo (Adptky en el comando NROPT) es una técnica que cambia a una matriz “más rígida” si se encuentran dificultades de convergencia, y vuelve DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 94 a cambiar a la matriz de rigidez tangente completa al converger la solución, dando como resultado una tasa de convergencia más rápida, lo que es una ventaja. (Eggert [6]). La matriz usada en la ecuación de Newton-Raphson, (ecuación(6.14)) se define como la suma de dos matrices: ( ) 1
K K K ξ ξ ( ( ( = + −
(6.25) donde: - S
= Matriz secante (o más estable). - T
= Matriz tangente. - ξ = Parámetro de descenso. El programa ajusta el parámetro de descenso( ) ξ adaptándolo durante las iteraciones de equilibrio de la siguiente manera: 1. Empieza cada substep (subpaso) usando la matriz tangente( ) 0 ξ = . 2. Observa el cambio en el residuo { }
R en las iteraciones de equilibrio: Si aumenta (indicando una posible divergencia): • Elimina la solución actual si 1 ξ < , establece 1 ξ = y rehace la iteración usando la matriz secante. • Si ξ ya es igual a 1, continua iterando. Si decrece (indicando una solución convergente): • Si 1 ξ = (matriz secante) y el residuo ha decrecido en tres iteraciones seguidas (o 2 si ξ se incrementó a 1 durante el proceso de iteración de equilibrio por el caso arriba indicado), entonces se reduce ξ por un factor de 1/4 (estableciéndolo en 0,25) y se continúa iterando. • Si 1 ξ < , se vuelve a disminuir por un factor de 1/4 y se continúa iterando. Una vez que ξ sea menor que 0,0156, se toma 0, 0 ξ = (se usa la matriz tangente). DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 95 3. Si se encuentra un mensaje de “pivote negativo” (lo cual indica un matriz mal condicionada): • Si 1 ξ < , se elimina la solución actual, se toma 1 ξ = y se rehace la iteración usando la matriz secante. • Si 1 ξ = , se bisecta el paso de carga si el comando “automatic time stepping” está activado, en caso de que no lo esté se termina la ejecución. Las no linealidades que hacen uso del descenso adaptativo (esto es, que forman una matriz secante si 0 ξ > ) incluyen: plasticidad, contacto, rigidez de tensiones con grandes deformaciones, campos magnéticos no lineales que usan la formulación potencial escalar, el elemento concreto SOLID65 con el KEYOPT(7) = 1, y el elemento de membrana laminar SHELL41 con el KEYOPT(1) = 2. El descenso adaptativo se usa por defecto en estos casos a no ser que las opciones de “búsqueda lineal” o “arc-length” estén activadas. Sólo está disponible con el “Procedimiento de solución completa de Newton-Raphson”, en el que la matriz se recalcula en cada iteración. El “Procedimiento de solución completa de Newton-Raphson” es también el que se usa por defecto para plasticidad, contacto y no linealidades debidas a grandes deformaciones. Búsqueda Lineal La opción de búsqueda lineal (a la que se accede con el comando LNSRCH) trata de mejorar una solución de Newton-Raphson { }
u ∆ escalando el vector solución por un valor escalar denominado el “Parámetro de búsqueda lineal”. Considerando la ecuación (6.15) de nuevo: { } { } { }
= + ∆ (6.26) DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 96 En algunas situaciones de la solución, el uso del { }
u ∆ completo conduce a inestabilidades de la solución. De ahí que, si se usa la opción de la búsqueda lineal, la ecuación (6.26) sea modificada para escribirse: { } { } { }
= + ∆ (6.27) donde: - s = Parámetro de búsqueda lineal; 0, 05 1, 0 s < < . s se determina automáticamente minimizando la energía del sistema, lo que se reduce a encontrar el cero de la ecuación no lineal: { } { } { } ( ) { } ( )
g u F F s u = ∆ − ∆ (6.28) donde: - s
g = Gradiente de la energía potencial con respecto a s . Un esquema de resolución iterativa basado en el “Método de Regula Falsi” (el cual, nos permite encontrar la raíz de la ecuación f(x)=0, donde f(x) es una función continua definida en el intervalo [a,b], con f(a) y f(b) de signos diferentes)se usa para resolver la ecuación (6.28) (Schweizerhof y Wriggers [7]). Las iteraciones continúan hasta que: 1. s
g sea menor que 0
0, 5g , donde 0
g es el valor de la ecuación (6.28) en 0, 0 s = (esto es, usando { }
en lugar de { } ( ) { }
F s u ∆ . 2. s
g no cambie significativamente entre iteraciones. 3. Se hayan ejecutado 6 iteraciones. Si 0, 0
g > , no se realiza ninguna iteración y s se toma como 1,0. No se permite un s menor de 0,05. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 97 La solución escalada { }
u ∆ se usa para actualizar los valores actuales de los GDL { }
en la ecuación (6.15), que se ejecutan en la siguiente ecuación de equilibrio. 6.2.4. Método del Arc-Length Los análisis con detalle de estructuras geométricamente no lineales requieren la creación de modelos matemáticos que incluyan con precisión las condiciones específicas de carga y soporte; y, aun más importante, que modelen la rigidez y la respuesta dada por la estructura. La principal característica del pandeo no lineal frente al pandeo de autovalores que calcula Ansys, es que el fenómeno de pandeo no lineal es capaz de resolver una región de inestabilidad durante la trayectoria de post-pandeo, mientras que el pandeo de autovalores sólo incluye el comportamiento lineal de la región de pre-pandeo hasta la bifurcación (punto de carga crítica). Esto se ilustra en la figura 6.8. En dicha figura también se observa que cuando se realiza el análisis no lineal, más ajustado a la realidad, el pandeo comienza a una carga crítica menor a la obtenida con un análisis lineal o, como lo llama Ansys, con un cálculo de pandeo de autovalores. Figura 6.8: ”Comportamiento de Pandeo No lineal vs Pandeo de autovalores” DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 98 Las estructuras geométricamente no lineales contienen a menudo puntos límite en los cuales la trayectoria de equilibrio tiene tangente horizontal, como se observa en la figura 6.8. La región comprendida entre dos puntos críticos es inestable, porque la recta tangente a la trayectoria de equilibrio es negativa, lo que indica un aumento de los desplazamientos al disminuir las cargas. Si el proceso de carga se hace sin control en desplazamientos, como es habitual (salvo en ensayos de laboratorio), cuando la carga crece superando el primer punto límite, la estructura sufre un “salto” brusco de desplazamientos, que la llevan a una nueva situación estable, como se muestra en la figura 6.9. A este fenómeno se le denomina “snap-through”. Figura 6.9: “Fenómenos de snap-through y snap-back. Respuesta en carga y descarga” Al descargar esta estructura, cuando la carga cae por debajo del segundo punto límite, se produce otro salto brusco de desplazamientos que se denomina “snap-back”. El fenómeno de “snap-through” se ilustra muy bien en el caso del arco rebajado cargado puntualmente de la figura 6.10. Figura 6.10: ”Snap Through en el Pandeo” DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 99 La presencia de puntos críticos de estabilidad y trayectorias de equilibrio inestables son las principales dificultades que las soluciones numéricas deben superar para capturar completamente la respuesta no lineal. Aun existen ciertos obstáculos para encontrar respuestas no lineales de estructuras sometidas a grandes cambios geométricos. Predecir “snap-through” y “snap-back” de ciertas estructuras es difícil y resulta computacionalmente costoso. También es difícil hallar cuánta carga adicional puede soportar con seguridad una estructura bajo estas circunstancias. El método del “arc-length” se presenta como el adecuado para predecir la respuesta correcta de estructuras con comportamientos complejos del tipo “snap-through” y seguir la trayectoria de equilibrio no lineal a través de los puntos límite. Introducción Las estructuras pueden experimentar condiciones de carga que causen grandes desplazamientos que cambien de forma significativa su geometría y requieran, por tanto, que las ecuaciones de equilibrio se planteen en la geometría deformada. Las grandes deflexiones son descritas por ecuaciones diferenciales no lineales que pueden ser resueltas usando técnicas incrementales, como el Método de Newton-
Raphson anteriormente presentado. En análisis no lineales la matriz de rigidez tangente sustituye a la matriz de rigidez que se usa en los análisis lineales. Se usan pasos de carga iterativos para aplicar la carga a la estructura en forma de pequeñas cargas incrementales y hallar cada vez los correspondientes desplazamientos incrementales. La representación de estos resultados define la curva de la trayectoria de equilibrio de la estructura bajo las cargas aplicadas. Una estructura que sufra grandes cambios en su geometría a menudo presenta puntos límite con una respuesta inestable de “saltos” (“snap-through” y “snap-back”) durante un colapso estático. La solución a estas inestabilidades estructurales es difícil de hallar con los métodos comunes de resolución de ecuaciones, tales como el método de Newton-Raphson. Estos métodos fallan cuando tengan lugar comportamientos de “snap-
back” a lo largo de la trayectoria de carga. Esto sucede porque debido a la naturaleza de DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 100 la ecuación empleada en el método iterativo de Newton-Raphson, es decir, la ecuación (6.14), no puede haber convergencia si la matriz de rigidez tangente (la pendiente de la curva fuerza-desplazamiento en cualquier punto) se hace cero. Ver figura 6.11. Figura 6.11:”Divergencia en el Método de Newton-Raphson” Los investigadores han estudiado continuamente estos problemas y se han obtenido mejoras en el proceso que han sido gradualmente introducidas en los programas comerciales de elementos finitos. El análisis geométricamente no lineal llevado a cabo mediante elementos finitos debe ser capaz de hallar todas las posibles respuestas durante la aplicación de grandes cargas. Pasar puntos críticos durante la respuesta geométricamente no lineal es retador. Dos tipos de puntos críticos hallados en este tipo de comportamiento son: - Puntos límite de carga que se alcanzan cuando la trayectoria de la respuesta tiene un “snap-through” local; - Y puntos límite de control que definen un “snap-back” local. En uno de estos puntos límite de control, la carga puede cambiar de dirección al cambiar las deflexiones cuando se atraviesa un máximo local. Existe una familia importante de métodos de resolución de ecuaciones no lineales llamados genéricamente métodos del “arc-length” (longitud del arco) que pueden superar estos puntos críticos. El método general del “arc-length” usado para análisis estructurales, DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 101 surge como una variación del método general de Newton-Raphson ideada para superar las dificultades de éste para pasar por puntos críticos. Fué desarrollado originalmente por Riks (1972; 1979) y Wempner (1971) [10]. La técnica se asemeja mucho al método de Newton-Raphson descrito en Riks [8] y [9] excepto en que en este caso el incremento de carga aplicado pasa a ser una incógnita adicional en el problema. No son pocas las variaciones que se desarrollaron sobre el trabajo original de Riks y Wempner. De hecho la forma del método del arc-length que usa el programa Ansys 8.0. para resolver ecuaciones no es la de Riks y Wempner sino la variación del método introducida por Forde y Stiemer [11] en 1987. Por este motivo y, dado que no es el objetivo de este proyecto hacer un análisis demasiado exhaustivo de la gran cantidad de métodos del “arc-length” de resolución de ecuaciones no lineales, se presentarán sólo dos: el planteamiento de Riks y Wempner, por ser la primera históricamente y base de las demás, y la forma de Forde y Stiemer, usada por Ansys. Análisis estático geométricamente no lineal mediante el método del “Arc-
Length” de Riks-Wempner En el método de Riks-Wempner la longitud de un vector tangente a la trayectoria de equilibrio se usa para hallar un nuevo punto que es intersección del plano normal a la tangente. Una carga dada por el analista estimará las magnitudes del incremento inicial de carga para cada paso. Se llega a la terminación del método cuando se cumplen las condiciones de máximo factor proporcional de carga o máximo desplazamiento nodal dados. El proceso también se acaba cuando se alcanza el máximo número de incrementos de carga en un paso determinado. El método del arc-length de Riks-Wempner traza la trayectoria no lineal de equilibrio usando un proceso iterativo que empieza computando los desplazamientos iniciales provocados por un incremento de carga dado. El método pasa a hallar el nuevo punto de equilibrio desde el punto inicial; i, como se puede apreciar en la figura 6.12. (Para mas detalle consultar Crisfield [12], Owen [13] y Riks [8] y [9]). En la figura 6.12 se muestra la curva carga-desplazamiento para un sistema de un solo grado de libertad. El vector tangente a la curva en i puede escribirse como: DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 102 i
(6.29) donde i
λ ∆ es la carga incremental aplicada en i, i
K es la matriz de rigidez tangente evaluada en i y i
es el vector incremental de desplazamiento hallado resolviendo i
K u F λ ∆ = ∆
(6.30) El vector normal, i
también se muestra en la figura 6.12 y se puede expresar como: i
(6.31) Figura 6.12: “Método del arc-length de Riks-Wempner” DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 103 Figura 6.13: “Inicio del método del arc-length de Riks-Wempner” Al comenzar el método del “arc-length” de Riks-Wempner un incremento inicial de carga,
λ ∆ , se usa para computar el primer vector desplazamiento,
, y la longitud del primer vector tangente
. Estas variables se muestran en la representación de carga – desplazamiento de la figura 6.13. Para hallar los desplazamientos iniciales,
, se usan triángulos semejantes. Durante este incremento inicial la rigidez tangente es la misma que la rigidez lineal. El incremento de carga es un parámetro dado por el analista que divide la carga total aplicada en incrementos iguales. Un incremento de carga dado empieza el proceso y a partir de él se hallan los desplazamientos 0
u ∆ usando la matriz de rigidez tangente
K . Los desplazamientos iniciales 0
u ∆ se hallan resolviendo 0
λ λ ∆
(6.32) donde 1 λ = y tot
se deduce de la expresión 0
K u F λ ∆ =
(6.33) La longitud del vector tangente
, que será la que determine el radio del arc-
length, a lo largo de la trayectoria de equilibrio se puede calcular como 2
s t t u u λ ∆ = ⋅ = ∆ + ∆ ∆
(6.34) DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 104 A través del resto de iteraciones la longitud de arco (arc-length) es constante o puede ser escalada por la introducción de un parámetro por parte del analista en la siguiente expresión 1
(6.35) El analista decide sobre el número de iteraciones requeridas,
, y sobre el número de iteraciones deseadas,
I . Método del arc-length usado por Ansys El método del arc-length usado por Ansys (que se activa con el comando ARCLEN,ON) es apropiado para soluciones no lineales estáticas de problemas inestables. Este método usa iteraciones esféricas explícitas para mantener la ortogonalidad entre el radio del arc-length y las direcciones ortogonales, como se describe en Forde y Stiemer (1987). Se asume que todas las magnitudes de carga están controladas por un único parámetro escalar (i.e., el factor de carga total). Al ser este un método “quasi-Newton” todas las opciones en el programa del método de Newton-Raphson son aún la base de la solución del arc-length. Ya que ahora, los vectores de desplazamiento y el factor escalar de carga se tratan como incógnitas, el método del arc-length es él mismo un método que automatiza los tamaños de los pasos de carga (No necesitará el comando AUTOTS,ON). Para problemas con giros bruscos en la curva carga-desplazamiento o materiales dependientes de la trayectoria de carga, es necesario limitar el radio del arc-length (el tamaño del incremento de carga del arc-length) usando el radio inicial o de referencia (con el comando NSUBST) y el rango de variación permitido. Durante la solución, el método del arc-length variará el radio del mismo en cada substep (subpaso) de acuerdo con el grado de no linealidades que estén incluidas. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 105 El rango de variación del radio del arc-length se limita por los multiplicadores máximo y mínimo del radio de referencia del arc-length (MAXARC y MINARC en el comando ARCLEN). En el procedimiento de resolución del método del arc-length usado por Ansys, la ecuación no lineal (6.14), se replantea asociada al factor de carga total λ : { } { } { }
K u F F λ ( ∆ = −
(6.36) donde λ está comprendido normalmente en el rango 1, 0 1, 0 λ − ≤ ≤ . Es interesante recalcar que el factor de carga λ sólo multiplica a las fuerzas aplicadas, ya que el otro término de fuerzas son las fuerzas internas que se calculan a partir de las aplicadas. Si escribimos el factor proporcional de carga λ de forma incremental, tenemos en el substep n y la iteración i (ver figura 6.14): { } { } ( ){ } { } { }
T a a nr
K u F F F R λ λ λ ( ∆ − ∆ = + − = −
(6.37) donde: - λ ∆ = Factor de carga incremental (tal y como se muestra en la figura 6.14). Figura 6.14: ”Aproximación mediante el arc-length de Forde-Stiemer con el método completo de Newton-Raphson” DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 106 El desplazamiento incremental { }
u ∆ puede escribirse en dos partes siguiendo la ecuación(6.37): { } { } { }
u u u λ ∆ = ∆ ∆ + ∆ (6.38) donde: - { }
u ∆ = Desplazamiento debido a un factor de carga unidad. - { }
u ∆ = Desplazamiento incremental proveniente del método convencional de Newton-Raphson. Estos se definen: { } { }
( ∆ =
(6.39) { } { }
( ∆ = −
(6.40) En cada iteración del arc-length, es necesario usar las ecuaciones (6.39) y (6.40) para resolver { }
u ∆ y { }
u ∆ . Además, al haber introducido una nueva incógnita en el problema, el factor de carga incremental λ ∆ , el método introduce una nueva restricción para cada substep que se conoce como la “ecuación del arc-length” que puede ser escrita, por ejemplo, en la iteración i (ver figuras 6.14 y 6.15)(es equivalente a la ecuación (6.34) del método de Riks-Wempner) { } { }
u u λ β = + ∆ ∆ # (6.41) donde: - β = Factor de escala (con unidades de desplazamiento
) usado para asegurar el correcto escalamiento en las ecuaciones (ya que λ es adimensional). - n
u ∆ = Suma de todos los desplazamientos incrementales i
u ∆ de esta iteración. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 107 Así el método de Forde-Stiemer usado por Ansys puede verse más claramente en la figura 6.15. Figura 6.15: ”Método del arc-length de Forde-Stiemer. Un substep” El radio del arc-length i
# se fuerza, durante las iteraciones de un mismo substep, a ser idéntico al radio 1
# de la primera iteración de ese substep, es decir, se impone que se mantenga constante a lo largo de un substep dado. 1 1
i i−
= = = # # # (6.42) Por lo tanto el método del arc-length de Forde-Stiemer que usa Ansys permite a la carga y al desplazamiento que varíen durante un substep o subpaso de carga tal y como se muestra en la figura 6.16. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 108 Figura 6.16: ”Comportamiento de la convergencia en Ansys usando el arc-length” El radio del arc-length 0
# en la primera iteración del primer substep de un paso de carga o load step se denomina “radio inicial o de referencia del arc-length” (definido por el comando NSUBST ya que este radio de referencia se calcula como el cociente entre la carga total aplicada en el paso de carga correspondiente y el número de substeps especificados en dicho comando para ese mismo paso de carga). Los radios del arc-
length en los substeps subsiguientes del mismo paso de carga o load step son calculados por el programa a partir del radio del arc-length del substep previo y el comportamiento de la solución y teniendo en cuenta que estos radios deben estar comprendidos dentro del rango límite que determinan el radio de referencia para ese paso de carga y los multiplicadores del mismo (dados en las etiquetas MAXARC y MINARC del comando ARCLEN). Para el método del arc-length, por tanto, no se usa el “Automatic time stepping”. 0 0
MINARC MAXARC ⋅ ≤ ≤ ⋅ # # # (6.43) Las ecuaciones (6.38) y (6.41) conjuntamente, determinan, ellas solas, el vector solución ( ) ,
u λ ∆ ∆ . Sin embargo, hay muchas formas de resolver λ ∆ aproximadamente. El método de la iteración esférica explícita se usa para asegurar la ortogonalidad (Forde y Stiemer (pág.174 del libro indicado en la bibliografía)). En este método, el residuo DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 109 requerido i
r (un escalar) para la iteración explícita en una esfera, se calcula de antemano como producto escalar de los vectores normal i
y tangencial i
. Los autores sugieren la selección de un vector normal i
con una dirección arbitraria con respecto al vector i
que es tangente a la configuración de carga-desplazamiento en un punto dado. Los vectores normal y tangencial constan de m componentes provenientes del vector desplazamiento y una proveniente del factor de carga. Estas componentes se combinan usando el factor escalar β formando vectores de m+1 componentes que pueden escribirse como: i
(6.44) i
∆ ¦ ¹
(6.45) Por consiguiente, el factor de carga incremental del arc-length se determina mediante la fórmula: { } { }
(6.46) Esta expresión puede ser simplificada para casos particulares de ortogonalidad. El método funciona bien incluso cuando se encuentra en las proximidades de puntos críticos en los que hay cambios bruscos en la solución. Finalmente, reseñar que los vectores solución se actualizan de acuerdo con (ver figuras 6.14 y 6.15): { } { } { } { }
1 i n n i
= + ∆ + ∆ (6.47) y 1 i n i
= + + ∆ (6.48) donde: - n = Número del substep actual. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 110 Los valores de n
λ y λ ∆ están disponibles en el módulo POST26 (comando SOLU) y se corresponden con las etiquetas ALLF y ALDLF, respectivamente. El radio normalizado del arc-length se encuentra en la etiqueta ARCL (comando SOLU) y se corresponde con el valor 0
# es el radio inicial o de referencia del arc-length en el substep i definido en el comando NSUBST, y que se relaciona con el factor de carga y los desplazamientos a través de la ecuación (6.41). En el caso en el que las cargas aplicadas sean mayores o menores que las máximas o mínimas cargas críticas, el arc-length continuará con las iteraciones en ciclos ya que λ no se acerca a la unidad. Es recomendable, por tanto, acabar con las iteraciones del arc-length imponiendo un criterio de terminación (puede hacerse usando los comandos ARCTRM o NCNV). 6.3. Resumen Una vez presentada la base teórica de los métodos de resolución de ecuaciones que usa Ansys y visto con más detalle el Método del Arc-length, convendría enumerar los puntos más importantes que se han presentado en este capítulo, haciendo hincapié sobre todo en aquellos que resulten más prácticos desde el punto de vista operativo en la resolución del problema no lineal planteado. En este pequeño resumen, se incluyen una serie de consejos prácticos, provenientes de bibliografía diversa, analistas y de la propia ayuda del programa, que orientan sobre cómo realizar un análisis no lineal en Ansys. • En análisis no lineales es muy recomendable aplicar la carga en forma de incrementos (pasos de carga o load steps) para poder seguir bien la trayectoria de la curva carga-desplazamiento. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 111 • El método de resolución numérica de ecuaciones más usado por su sencillez y antigüedad es el “Método iterativo de Newton-Raphson” base de casi todos los demás y que se describe en la ecuación (6.16): { } { } { }
n i i n n i
. • En análisis no lineales en los que se esperen comportamientos de “snap through” y “snap back” o en los que se quiera trazar la trayectoria de la solución en la región de post-pandeo, el método de Newton-Raphson antes mencionado presenta problemas debido a la existencia de puntos críticos de pendiente horizontal en los que se anula la matriz de rigidez. En estos casos es aconsejable usar el Método del Arc-length que puede pasar bien esos puntos críticos. Una comparativa ilustrativa de ambos métodos puede observarse en la figura 6.17. Figura 6.17:”Newton-Raphson vs Arc-length”. • El método del arc-length que usa Ansys es la forma de Forde y Stiemer del mismo (puede observarse en las figuras 6.14 y 6.16). • El procedimiento básico de resolución del problema usando el método del arc-
length en Ansys puede resumirse en los siguientes pasos: - El método se activa con el comando ARCLEN,ON. - La carga se aplica a la estructura incrementalmente a través de una serie de pasos de carga o load steps siendo la carga aplicada en cada uno de ellos un valor definido por el analista. A su vez cada paso de carga se resuelve a través de un número de subpasos de carga o substeps que se determinan a través DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 112 del comando NSUBST. Los substeps se resuelven mediante una serie de iteraciones que siguen hasta que: la solución converja en el substep determinado (ya que se impone que la solución esté en equilibrio en cada substep) o bien se sobrepase sin converger el número máximo de iteraciones permitidas en un substep dado por el comando NEQIT. - El “radio del arc-length de referencia” es el radio del arc-length de la primera iteración del primer substep de un paso de carga (en realidad es el radio de todo el primer substep porque Ansys impone que el radio sea constante en las iteraciones de un mismo substep). Este radio lo define el analista mediante el comando NSUBST y la carga que aplica en el “loadstep” ya que se calcula con la expresión: 0
" " " " i
Fracción de la carga total aplicada en el LS
Número de substeps del LS definidos en NSUBST
= # (6.49) - Los radios de los siguientes substeps, ya no los define directamente el analista, si no que Ansys los calcula basándose en el radio del arc-length del substep previo y el comportamiento no lineal de la solución. Ahora bien, el analista determina el rango en el que se moverán estos valores imponiendo límites de variación al radio del arc-length de referencia (que por tanto sirve como referencia para el tamaño del resto de los radios del mismo paso de carga) a través de las etiquetas MAXARC y MINARC del comando ARCLEN. Los límites de variación del radio del arc-length en los substeps de un loadstep determinado son: 0
* Límite inferior MINARC = # (6.50) 0
* Límite superior MAXARC = # (6.51) donde 0
# es el radio de referencia del arc-length (el del primer substep) para ese paso de carga. - Cada paso de carga o loadstep, se resuelve de la manera indicada hasta que: La solución converja. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 113 Se alcance alguna de las condiciones de salida del análisis del programa que pueden especificarse en muchos comandos pero sobre todo en: ARCTRM, NEQIT o NCNV. (Si la solución no consigue converger en el número de iteraciones de equilibrio indicadas en NEQIT, el programa bisectará automáticamente el radio del arc-length y continuará con el análisis. La bisección continuará hasta que se obtenga una solución convergida o se alcance el límite inferior del arc-length en ese LS). • Según se definen los radios de las diferentes iteraciones del método, la carga aplicada en un substep determinado diferente del primero en la resolución de un paso de carga no es fácil de conocer de forma exacta, ya que el valor de la carga aplicada varía a lo largo del arco esférico. Sin embargo con otros métodos como el de Newton-Raphson la carga sí es perfectamente conocida. Esto puede apreciarse por ejemplo en la figura 6.17. De ahí que, aunque también pueden obtenerse cargas críticas con el método del Arc-
length ajustando bien sus parámetros, siempre que sea posible se recomienda que se obtengan mediante otros métodos como el de Newton-Raphson cuya precisión será mayor. Por todo esto, cuando el objetivo es determinar la carga crítica de pandeo de una estructura, se aconseja realizar un análisis lineal previo de pandeo de autovalores antes de realizar el análisis de pandeo no lineal. • El arc-length es un método muy potente para hacer frente a los análisis no lineales mediante elementos finitos en el cálculo de estructuras, pero su correcta aplicación requiere cierta pericia, ya que hay que elegir los parámetros que lo gobiernan con mucho cuidado, se debe hacer una buena estimación del radio de referencia que se usará para resolver un paso de carga llegando a un compromiso entre tiempo de ejecución y precisión. • Cuando se usa el arc-length, es mejor basar el criterio de convergencia en fuerzas [CNVTOL,F] y no en desplazamientos [CNVTOL,U]. • El factor de carga total del arc-length, λ , (que puede encontrarse en la etiqueta ALLF del comando SOLU) puede ser tanto positivo como negativo. De manera análoga, el valor de TIME en la solución, que en un análisis del arc-length se relaciona con el factor DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA 114 de carga total, λ , también puede ser positivo o negativo. La interpretación de valores negativos en ALLF o TIME indica, simplemente que el método del arc-length está aplicando la carga en la dirección contraria a la definida en el problema para mantener la estabilidad en la estructura. Cuando esto ocurre, se corre el peligro de que el arc-length vuelva atrás siguiendo la misma trayectoria de carga, pero en sentido inverso, produciéndose el fenómeno denominado “drifting back”. El fenómeno de “drifting back” suele deberse a una mala elección del tamaño del radio del arc-length (radios demasiado grandes o demasiado pequeños). Los valores negativos de ALLF o TIME se encuentran de forma habitual cuando se dan comportamientos de “snap-through” en la estructura. DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA cálculo no lineal se basan sobre todo en el método de Newton-Raphson y en variaciones del mismo. Se verá cómo el programa los lleva a cabo, con las diferentes órdenes asociadas y los parámetros de que dispone para controlar estos métodos. Por último se hará especial hincapié en el método del “arc-length” que es el más apropiado para el caso que se estudia en el presente proyecto.
Pandeo de Autovalores El análisis de pandeo de autovalores predice la carga crítica de pandeo de una
estructura ideal lineal y elástica (es decir el punto de bifurcación en la curva cargadesplazamiento de la estructura). Este método se corresponde con el método clásico del análisis de pandeo elástico que se puede encontrar en los libros; por ejemplo, al hacer el análisis de pandeo de autovalores de un cilindro hueco de espesor de pared constante sometido a compresión pura los resultados se corresponden con la solución clásica de Euler (4.36). Sin embargo, las imperfecciones y otras no linealidades hacen que la mayoría de las estructuras que podemos encontrar en la vida real no alcancen esta carga crítica clásica de pandeo, sino que la carga real de pandeo será menor de la prevista por el análisis elástico. Por esto este análisis sólo debe hacerse como paso previo de un análisis no lineal o como una primera aproximación a la carga crítica de una estructura ya que proporciona valores de la carga crítica mayores a los reales y que son, por tanto, no seguros. Ansys realiza el análisis de pandeo de autovalores usando un modelo linealizado de estabilidad elástica. El punto de bifurcación se corresponde, en este caso, con el crecimiento no acotado de un nuevo patrón de deformación. Una estructura con una curva de comportamiento carga-desplazamiento como la mostrada en la figura 6.1.a puede ser analizada bien por este método, mientras que si muestra el comportamiento mostrado en la figura 6.1.b no se analiza bien.
E.T.S. INGENIEROS, UNIVERSIDAD DE SEVILLA
DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA
Figura 6.1: “Validez de los resultados de un análisis de pandeo de autovalores”
En cualquier cálculo de pandeo de una estructura se debe tener en cuenta el fenómeno de “Stress stiffening”. Este es el fenómeno por el cual las fuerzas de membrana (axiles) influyen en la rigidez a deflexión lateral; por ejemplo, una cuerda de guitarra se hace más rígida al tensarla (su frecuencia natural de vibración es más alta), pero una carga axial de compresión reduce la rigidez a flexión haciendo que la estructura se haga inestable al alcanzar cierta carga: la carga crítica de pandeo.
Este fenómeno de 2º orden conocido como “Stress stiffening” se tiene en cuenta al realizar el análisis de pandeo de autovalores de Ansys incluyendo un cálculo tensional previo para incluir los efectos de las tensiones de membrana en la matriz de rigidez geométrica (“stress stiffness matrix”) que se suma a la obtenida con la teoría de primer orden.
K ss = K + Kσ
La matriz de rigidez geométrica depende de: - Geometría del elemento. - Desplazamientos. - Valores de los esfuerzos de membrana.
DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. ESTUDIO DE LA CONVERGENCIA 6.1.1. Proceso del cálculo de Pandeo de Autovalores Cuando se analiza el pandeo lineal en Ansys, el problema de pandeo es formulado como un problema de autovalores siguiendo el siguiente proceso: 1. Se carga la estructura hasta un nivel arbitrario de referencia de cargas exteriores, Fref .
Se hace un análisis estático lineal para calcular los esfuerzos de membrana en los elementos
K ⋅ Dref = Fref
De la ecuación (6.2) calculamos los desplazamientos Dref y después las tensiones y esfuerzos correspondientes al nivel de carga Fref .
Se calcula la matriz de rigidez geométrica para ese nivel de carga Kσ ,ref , para ello se hace el análisis estático previo antes mencionado, para el cual se introduce en el programa el código expuesto en la figura 6.2:
/SOLU ANTYPE,STATIC PSTRES,ON ! Aplicar cargas SOLVE FINISH
Figura 6.2: “Código de Ansys para el cálculo previo de la matriz de rigidez geométrica”
2. Para otro nivel de cargas exteriores se tiene:
F = λ ⋅ Fref ⇒ Kσ = λ ⋅ Kσ ,ref
Los desplazamientos para el nuevo valor de carga son
(K + λ ⋅ K
σ . ESTUDIO DE LA CONVERGENCIA Al ser un análisis lineal.4)
3.2. pero sin cambiar la distribución de tensiones.5) y (6. se supone que al multiplicar las cargas por λ también se multiplican las tensiones por el mismo factor.1.0. El pandeo (bifurcación del equilibrio) sucede cuando es posible un δD 0 con δF=0:
realizar el análisis:
σ . UNIVERSIDAD DE SEVILLA
.4) se obtiene el problema de autovalores al que se alude al
autovector es el modo de pandeo. ref
) ⋅δ D = 0
(6. Si se aplica una perturbación sobre la carga δF. Por tanto.S. Los autovectores están normalizados de manera que el componente más grande es 1.3:
E. se produce un movimiento δD.ref
)⋅ D
= λ ⋅ Fref
(6.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. INGENIEROS. Su
El problema de autovalores de la ecuación (6. Una muestra del código que se introduce en Ansys para resolver el problema de pandeo por autovalores se muestra en la figura 6.0.6)
El menor de los autovalores nos proporciona la carga crítica de pandeo λcr .5)
Restando (6.T.6) se resuelve mediante alguno de los procedimientos explicados en 6.ref
)⋅( D
+ δ D ) = λ ⋅ Fref
(6. las tensiones (cuando se calculan como resultado del problema) deben ser interpretadas como una distribución relativa de tensiones y no directamente como el estado tensional.
0. Para evitar esta dificultad.1. ESTUDIO DE LA CONVERGENCIA
/SOLU !Debe haberse resuelto previamente un problema estático con PSTRES.T.1.S. aplicar un desplazamiento inicial de la búsqueda (con el comando BUCOPT) hasta un valor próximo al autovalor esperado. BUCOPT.BUCKLE ! Análisis de pandeo !por autovalores.1 ! Expande la forma !de los 10 modos calculados. MXPAND.
Problema de extracción de autovalores Ansys dispone de varios métodos para la extracción de autovalores que se
enumeran en la tabla 6.
E. UNIVERSIDAD DE SEVILLA
.10 ! Define el método de !cálculo de autovalores y ! el nº de modos a calcular. Los que se usan en el análisis del pandeo de autovalores están recuadrados en rojo en dicha tabla. SOLVE FINISH
Figura 6.3: “Código de Ansys para realizar el cálculo del problema de Pandeo por Autovalores”
Si el primer autovalor es negativo (lo cual indica que el pandeo se produce para cargas aplicadas en el sentido opuesto).DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. En el método del subespacio.LANB. INGENIEROS. se muestra un error “Number of stress-stiffness DOF is less than requested modes”.ON ANTYPE. Ansys da un warning y el programa termina.
BLOCK LANCZOS – Válido en los mismos casos que Subspace. asociadas a una serie de g.d. ya que las matrices reducidas sólo permiten obtener una solución aproximada. INGENIEROS.d. maestros que el usuario (o el propio ANSYS) debe seleccionar previamente.S. No resulta necesario definir g.T.l.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. El proceso de solución es más rápido que con las opciones SUBSPACE pero menos exacto.
SUBSPACE – Permiten obtener un determinado número M (reducido: M<<N) de autovalores y autovectores.: •
REDUCED – Emplea matrices reducidas.l. Es el más eficiente en el caso de que haya autovalores próximos
E. ESTUDIO DE LA CONVERGENCIA
Tabla 6. UNIVERSIDAD DE SEVILLA
.0. maestros. pero tiene una convergencia más rápida.1: “Métodos de extracción de Autovalores en Ansys 8.d. Aplicable en problemas de autovalores con matrices simétricas grandes.0”
ANSYS ofrece varios métodos para la obtención de autovalores y autovectores en un sistema de N g.
BUCKLE). Las mejoras de este método sugeridas por Wilson e Itoh [16] también se han incluido en Ansys. que se usa por defecto. Como se ha mencionado anteriormente. [M ]
es la matriz de masas (o de rigidez geométrica) ensamblada.0.
Se hace una comprobación de “secuencia de Sturm” si este es un punto de pivote distinto del inicial y se ha requerido (mediante los comandos Strmck=ALL. los métodos de extracción de ecuaciones más frecuentemente usados por Ansys en la extracción de autovalores del problema de Pandeo de autovalores son el método del Subespacio y el de Block Lanczos y por ello los veremos con más detenimiento. Método del Subespacio El método iterativo del subespacio (que se activa con el comando BUCOPT. Inicializar los vectores iniciales [ X 0 ]. Se define el pivote inicial s : • En un análisis de pandeo (ANTYPE. se repiten los pasos del 5 al 14:
E. Para cada iteración n del subespacio (desde 1 hasta N M ). ESTUDIO DE LA CONVERGENCIA (es el caso del pandeo por autovalores. 4. 3.SUBSP) se describe con detalle en Bathe [5]. UNIVERSIDAD DE SEVILLA
.T.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. s = SHIFT en el comando BUCOPT (su valor por defecto es 0. INGENIEROS.7)
[ K ] es la matriz de rigidez ensamblada.0). Triangularizar la matriz pivotada
K *  = [K ] + s [M ]  
(6. 2. El algoritmo básico consiste en los siguientes pasos: 1. o PART en el comando SUBOPT). por tanto es el método que usaremos).S.
M  =  X n  [M ]  X n       
10. Se resuelve para  X n  :  
K*   X n  = [F ]   
E.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. 7. Se ajusta para el pivote  K *  =  K  + s  M        11.9) (6. Computar los autovalores y vectores del subespacio usando la iteración generalizada de Jacobi:
 K *  [Q ] =  M  [Q ]{λn }    
(6. 6. Se escalan los vectores  X n  por   8. ESTUDIO DE LA CONVERGENCIA donde N M es el máximo número de iteraciones del subespacio (se pasa como entrada en NUMSSI del comando SUBOPT).FRONT) o el iterative PCG solver (EQSLV. Se forma [ F ] = [ M ][ X n −1 ] y se escala [ F ] por {λn −1} donde {λn −1} son los autovalores previamente estimados. Se definen las matrices de subespacio  K  y  M  :    
K  =  X n     
[K ]  X n   
(6.PCG). 5. Se ortogonalizan los vectores
{( λ
− s ) λn −1}
vectores previamente convergidos
(Ortogonalización de Gram-Schmidt).11)
[Q ] son los autovectores del subespacio. 9. INGENIEROS. UNIVERSIDAD DE SEVILLA
.T.S.8)
Estas ecuaciones se resuelven usando el frontal direct equation solver de Ansys (EQSLV.
12. Se realiza una comprobación de la secuencia de Sturm si se pide (Strmck=ALL. Ir a la nueva iteración. como se detalla en Grimes [17]. INGENIEROS. es la base teórica de este extractor de autovalores. ESTUDIO DE LA CONVERGENCIA -
{λn } son los autovalores actualizados. ir al paso 3. El procedimiento de ortogonalización de Gram-Schmidt (paso 8) asegura que los autovectores que no han convergido permanezcan ortogonales a los vectores convergidos sobre los que no se itera. Si se encuentra algún modo redundante o negativo. Si se requiere un nuevo pivote (como se indica más abajo). ir al paso 15.S. pero consigue una tasa de convergencia más rápida. Método de Block Lanczos El método de extracción de autovalores de Block Lanczos (es el algoritmo de Lanczos por bloques al que se accede con el comando BUCOPT. UNIVERSIDAD DE SEVILLA
. combinada con comprobaciones
E.LANB) está disponible para problemas de gran tamaño de autovalores con simetría. el autovector asociado no se vuelve a iterar.
15.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. Se comprueba la convergencia (con la tolerancia definida por el analista en Ansys): • • • Si todos los modos requeridos convergen. Actualizar la aproximación a los autovectores:
[ X n ] =  X n  [Q ]  
13. Típicamente. El método usado por el análisis modal emplea una estrategia de cambio de pivote automática. se eliminan y se crea un nuevo vector aleatorio.T. 14. paso 4. este solucionador es aplicable al tipo de problemas resolubles con el método de autovalores del subespacio. El algoritmo de pivote por bloques de Lanczos . Los pasos del 5 al 12 sólo se realizan en los vectores no convergidos: Una vez que un autovalor ha convergido. por defecto en el comando SUBOPT).
0. UNIVERSIDAD DE SEVILLA
. Este método es tan preciso como el Método del Subespacio y además es más rápido. cuando se usa una frecuencia dada como pivote (especificada en el comando FREQB) para extraer n modos más allá de los hallados con el pivote inicial (valor inicial de FREQB). La comprobación de las secuencias de Sturm. también asegura que se encuentra el número de autovalores requeridos más allá de la frecuencia proporcionada por el analista como pivote inicial (FREQE en el comando MODOPT) se encuentran sin obviar ningún modo.S. La tasa de convergencia de las frecuencias de autovalores en las partes media y alta del espectro es casi la misma que cuando se extraen los modos más bajos.
E. Ansys utiliza los Multiplicadores de Lagrange para tratar las ecuaciones de restricción en el extractor de Lanczos por bloques en lugar de eliminar dichas ecuaciones explícitamente antes de escribir las matrices completas en el archivo de cálculo. Por esta razón. El método de Block Lanczos es especialmente potente cuando busca frecuencias de autovalores en partes determinadas del espectro de autovalores de un sistema (en el problema de pandeo de autovalores normalmente los autovalores están próximos entre sí). para extraer el número de autovalores requerido. El algoritmo de Block Lanczos es una variación del algoritmo de Lanczos clásico.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. por tanto. Comparación: Método de Block Lanczos vs Método del Subespacio El Método de extracción de autovalores de Block Lanczos es el que usa Ansys por defecto en el cálculo de pandeo por autovalores. Más detalles sobre el método clásico de Lanczos se pueden encontrar en Rajakumar y Rogers [18]. donde las recursiones de Lanczos se efectúan usando un bloque de vectores. Para más detalles sobre la formulación de la teoría de los Multiplicadores de Lagrange ver Cook [19]. en lugar de un solo vector.000 gdl.T. El uso del método de Lanczos por bloques (o método Block Lanczos) para resolver grandes modelos (100. el algoritmo extrae estos n modos adicionales prácticamente a la misma velocidad que extrajo los n modos más bajos. ESTUDIO DE LA CONVERGENCIA de las secuencias de Sturm. por ejemplo) con muchas restricciones puede requerir una cantidad de memoria del ordenador bastante significativa. INGENIEROS.
2. A modo de resumen se muestra la tabla 6. 6.0. Por ello.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.13)
.2: “Comparación esquemática de varios extractores de autovalores de Ansys”
6.l maestros). INGENIEROS.1.d.2.d. el método del Subespacio sólo se usa cuando sea necesaria una gran precisión o cuando no sea práctico seleccionar g.S.
Tabla 6.2 en la que se comparan esquemáticamente los distintos métodos de extracción de autovalores usados por Ansys y que se han comentado en esta sección. el método del Subespacio es muy preciso debido a que usa las matrices [ K ] y [ M ] completas. aunque precisamente esa sea la causa de que sea más lento que el método Reducido o el de Block Lanczos.
Pandeo No-lineal Método de Newton-Raphson: Visión general El proceso de discretización de los elementos finitos produce una serie de
ecuaciones simultáneas:
[ K ]{u} = { F a }
E.l maestros (como se indicó anteriormente el método Reducido usa estos g. ESTUDIO DE LA CONVERGENCIA Por su parte. UNIVERSIDAD DE SEVILLA
(6.T.
{ui } es el vector de desplazamientos y
es el vector de fuerzas
restauradoras calculado a partir de las tensiones elementales. En un análisis térmico. {ui } es el vector de temperaturas y { Fi nr } es el  
vector de cargas resistentes calculado a partir de los flujos de calor elementales. La
parte derecha de la ecuación (6. INGENIEROS. la cantidad en que el sistema esta fuera del equilibrio.T.0.  K iT  es la matriz de   rigidez tangente. El método de Newton-Raphson es un proceso iterativo para resolver las ecuaciones no lineales que puede ser escrito como (Bathe[5]):
 K iT  {∆ui } = { F a } − { Fi nr }   {ui +1} = {ui } + {∆ui }
(6.14) (6. En la figura 6.14) es el residuo o vector de cargas desequilibradas. i.15)
 K iT  = Jacobiano de la matriz (matriz tangente)  
i = Subíndice que representa la iteración de equilibrio actual
Vector de las fuerzas restauradoras correspondiente a las cargas
elementales internas Ambos  K iT  y  
se evalúan basándose en los valores dados por {ui } . En un
E. En un análisis estructural..S. UNIVERSIDAD DE SEVILLA
.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.13) es una ecuación no lineal.e.4 se describe gráficamente la iteración i-ésima del algoritmo de Newton-Raphson aplicado a un problema de un solo grado de libertad.
 K iT  es la matriz de conductividad. ESTUDIO DE LA CONVERGENCIA donde: -
[ K ] = Matriz de coeficientes
{u} = Vector de incógnitas GDL (grados de libertad)
{F } = Vector de cargas aplicadas
Si la matriz de coeficientes [ K ] es ella misma función de los grados de libertad que son incógnita (o de sus derivadas) entonces la ecuación (6.
15)).0.S. En el primer “time step”. {u0 } = {0} . El algoritmo general es como sigue: 1. ESTUDIO DE LA CONVERGENCIA análisis electromagnético. 5. 4. Como se puede observar en las siguientes figuras. En un análisis transitorio. {ui } es el vector de potenciales   magnéticos y
nr i nr
es el vector de cargas resistentes calculado a partir de los flujos
magnéticos elementales. UNIVERSIDAD DE SEVILLA
. Se supone {u0 } .
E.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. 2. se requiere más de una iteración del método de Newton-Raphson para obtener una solución que converja.T. Añadir {∆ui } a {ui } para obtener la nueva aproximación {ui +1} (ecuación (6.  K iT  es la matriz de coeficientes   efectiva y
es el vector de cargas efectivas aplicadas que incluye los efectos de
inercia y de amortiguamiento. Repetir los pasos desde el 2 al 4 hasta obtener la convergencia. INGENIEROS. Calcular {∆ui } mediante la ecuación (6. Calcular la matriz tangente modificada  K iT  y la carga restauradora Fi nr   a partir de la configuración de los {ui } .  K iT  es la matriz de Dirichlet.
{u0 }
es normalmente la solución convergida del “time
step” (paso de carga) previo.14). 3.
INGENIEROS. de forma que el vector de cargas
restauradoras Fi nr
{ } (computado a partir del estado actual de tensiones.
La solución obtenida al final del proceso iterativo correspondería al nivel de carga
{F } .S.
E.4. UNIVERSIDAD DE SEVILLA
.0.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.T. Ninguna de las soluciones intermedias estaría en equilibrio.4: "Solución de Newton-Raphson . La solución final convergida estaría en equilibrio. Cada uno de estos pasos intermedios se denomina “iteración de equilibrio” (equilibrium iteration) en Ansys.5. ESTUDIO DE LA CONVERGENCIA
Figura 6. se muestra la solución de la siguiente iteración ( i + 1) a la del ejemplo de la figura 6.Una iteración"
En la figura 6. flujos de calor…) {F }
igualaría al vector de cargas aplicadas
(o al menos estaría dentro de cierta
tolerancia). Las iteraciones subsecuentes se realizarían de una manera análoga.
Esto se consigue.T. incluso de no haber presente una no
E.6.Siguiente iteración"
Si el análisis incluye no linealidades dependientes de la trayectoria (como por ejemplo plasticidad). y la iteración i. efectivamente. si y solo si la solución en cualquier iteración
{ui } esta “cerca” de la solución exacta. entonces el proceso de solución requiere que ciertos pasos intermedios estén en equilibrio para poder seguir correctamente la trayectoria de carga.i
Cada paso se denomina “paso de carga” (load step) en Ansys. UNIVERSIDAD DE SEVILLA
donde: T  K n .e.
{F } = Vector de las fuerzas restauradoras para el paso n. INGENIEROS. Además cada paso de carga puede subdividirse a su vez en “sub-pasos de carga” (substeps) en Ansys.i  {∆un . el vector final de carga F a se alcanza aplicando la carga en incrementos y aplicando las iteraciones de Newton-Raphson en cada paso:
T  K n.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.S. i. especificando un análisis incremental paso a paso. ESTUDIO DE LA CONVERGENCIA
Figura 6.i  
(6.i } = { Fna } − { Fnnr } .0. al final de los cuales también se alcanza el equilibrio. y la iteración i.
nr n ..  
{F } = Vector de fuerzas totales aplicadas en el paso n.i  = Matriz tangente para el paso n. Este proceso es el “Procedimiento Incremental de Newton-Raphson” y se muestra en la figura 6. El procedimiento de Newton-Raphson garantiza la convergencia.5: "Solución de Newton-Raphson . Por tanto.
la aproximación incremental (i. Si existiese un elemento “multistatus” en el modelo. El uso del “Procedimiento de Rigidez Inicial”. la matriz de rigidez podría ser recalculada con menor frecuencia usando el “Procedimiento Modificado de NewtonRaphson”. previene cualquier actualización o recálculo de la matriz de rigidez.INIT). Especialmente. sin embargo.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.14) y en la ecuación (6. en Ansys: (NROPT.UNSYM). la matriz de rigidez sería recalculada en la iteración en la cual se cambiara de estado independientemente de la opción de Newton-Raphson que se haya usado. pero requieren menos reformulaciones e inversiones matriciales. para análisis estáticos o transitorios. ESTUDIO DE LA CONVERGENCIA linealidad dependiente de la trayectoria.MODI). respectivamente. Los procedimientos modificado y de rigidez inicial de Newton-Raphson.
Figura 6.16)) el proceso se denomina “procedimiento de solución completa de Newton-Raphson” (Full Newton-Raphson solution procedure). Unos pocos elementos forman una matriz de rigidez tangente aproximada por lo que las características de convergencia son en cierto modo diferentes. aplicando las cargas en incrementos mas pequeños) es a veces necesaria para obtener una solución que se corresponda con el nivel final de carga.0. Alternativamente.e. se recalcularía sólo durante la primera o la segunda iteración de cada subpaso.
E.6: “Procedimiento Incremental de Newton-Raphson”
Cuando la matriz de rigidez se recalcula en cada iteración (como se indica en la ecuación (6. convergen más lentamente que el procedimiento de solución completa de Newton-Raphson. en Ansys: ( NROPT. como se muestra en la figura 6. en Ansys: (NROPT. UNIVERSIDAD DE SEVILLA
.FULL o NROPT..T.S.7. INGENIEROS.
{∆ui } < ε u uref
(Convergencia incremental de los GDL)
(6.T.S.
Se asume que hay convergencia cuando:
{R} < ε R Rref
(Convergencia residual)
El proceso iterativo descrito en la sección anterior continúa hasta que se consigue la convergencia.2.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0. UNIVERSIDAD DE SEVILLA
. INGENIEROS.19)
E.7: “Rigidez-Inicial Newton-Raphson”
6. ESTUDIO DE LA CONVERGENCIA
Figura 6.2.18)
donde { R} es el vector residual:
{ R} = {F a } − { F nr }
la norma L1 es la suma de los valores absolutos de los términos.
{∆ui } es
vector de incrementos de los GDL. La convergencia. INGENIEROS. • es la norma de un vector. Por defecto.S. una medida escalar de la magnitud del vector (más adelante se definen los tres tipos de norma que usa Ansys). Norma infinita: 2.
Hay tres normas disponibles entre las que elegir (NORM en el comando CNVTOL de Ansys): 1. esto es. ε R y ε u son tolerancias (TOLER en el comando CNVTOL de Ansys) y Rref y uref son valores de referencia (VALUE en el comando CNVTOL de Ansys). ESTUDIO DE LA CONVERGENCIA que es el término de la derecha de la ecuación de Newton-Raphson(6. también llamada norma Euclídea. basta con sustituir ∆u por R en las anteriores ecuaciones.17).0. Norma L1: 3. y la norma L2 es la raíz cuadrada de la suma de los cuadrados (SRSS) de los valores de los términos. UNIVERSIDAD DE SEVILLA
.14).DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. la de la ecuación (6.001 (tanto para ε R como para ε u ). Por defecto se usa la norma L2.
La norma infinita es simplemente el máximo valor en el vector (máximo residuo o máximo incremento de GDL). por tanto. El valor por defecto de las tolerancias es 0.
E. el programa Ansys sólo comprueba la convergencia residual.T. se obtiene cuando la norma del residuo (desequilibrio) es menor que una tolerancia multiplicada por un valor de referencia y/o cuando el tamaño de los incrementos de los GDL es menor que una tolerancia multiplicada por un valor de referencia. Norma L2:
{ R} ∞ = max Ri
{R} 1 = ∑ Ri
{R} 2 = ( ∑ Ri2 ) 2
Para la convergencia incremental de los GDL.
{ } en esos GDL se usan
{F } cae por debajo de 1.0 como valor.0E-6. si
entonces Rref toma 1. Si después de realizar NEQIT iteraciones de equilibrio sin alcanzar la convergencia.0. los F nr para computar Rref . UNIVERSIDAD DE SEVILLA
.. Ansys da un mensaje de error y termina el cálculo.0E-6 como valor. entonces
{F } cae por debajo de 1. El valor por defecto de referencia de uref es
.T.S. Para GDL estructurales. Para el resto de GDL. ESTUDIO DE LA CONVERGENCIA El valor por defecto de referencia de residuo Rref es
{F } . Para los GDL con
condiciones de contorno de desplazamientos impuestos. Esto ocurre a menudo en análisis del movimiento de cuerpos
rígidos (por ej.0.0. rotación libre).
Rref toma 1. Rref toma el valor 0.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. Para GDL térmicos. INGENIEROS. se fija con el comando NEQIT.
El número máximo de ecuaciones de equilibrio permitidas. N en cada subpaso.
extrapola la solución de GDL usando la historia previa para una mejor estimación de la siguiente solución.2.0 n
(6. La opción
“predictor” de Ansys (comando PRED).0.3.
E.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. UNIVERSIDAD DE SEVILLA
La matriz de rigidez tangente
 K n .20)
{∆ u n } =
previo. En análisis estáticos.21)
β se define como:
∆tn ∆tn −1
Incremento del desplazamiento acumulado sobre el paso de carga
n = Paso de carga actual.0
basan en esta configuración. ESTUDIO DE LA CONVERGENCIA
{un−1} . INGENIEROS. programa
Variaciones del método de Newton-Raphson y opciones que usa el
La solución usada para el comienzo de cada paso de carga n generalmente igual a la solución actual de GDL
{u } es
n .T.S. la predicción se basa en usar los incrementos de los desplazamientos acumulados sobre el paso de carga previo.
{∆un } = ∑ {∆ui }
∆tn = Tamaño del paso de carga actual. multiplicados como factor por el tamaño del paso de carga para hallar el un .0  y la carga restauradora  
{F } se
n .0 :
{u } = {u } + β {∆u }
Las subsecuentes iteraciones de equilibrio proporcionan los incrementos de los GDL {∆u} con respecto a los valores predichos de los mismos un . magnéticos y otros sistemas de primer orden. la predicción se basa en las velocidades y aceleraciones actuales usando las fórmulas de Newmark para GDL estructurales:
! {u } = {u } + {u } ∆t
n .∆tn −1 = Tamaño del paso de carga precedente. UNIVERSIDAD DE SEVILLA
∆tn = Tamaño de paso de carga actual. y vuelve
E. de ahí que este sea un algoritmo predictor-corrector. Descenso Adaptativo
El descenso adaptativo (Adptky
en el comando NROPT) es una técnica que
cambia a una matriz “más rígida” si se encuentran dificultades de convergencia.0 .
Para sistemas térmicos.
! {un −1} = Tasas actuales de estas cantidades. {un−1} .T.24)
{un −1} = Temperaturas actuales (o potenciales magnéticos). velocidades y aceleraciones actuales.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.S. En análisis transitorios.
θ = Parámetro de tiempo de integración trapezoidal (se introduce con el
comando TINTP). ESTUDIO DE LA CONVERGENCIA .0 n −1 n −1
1  2 !! +  − α  {un −1} ∆tn 2  
(6. No se permite que β sea mayor que 5. la predicción se basa en la fórmula trapezoidal:
! {u } = {u } + (1 − θ ){u } ∆t
n . {un−1} = Desplazamientos.0.0 n −1 n −1
! !! {un −1} .
α = Parámetro de Newmark (se introduce con el comando TINTP). INGENIEROS.
se vuelve a disminuir por un factor de 1/4 y se continúa iterando.
El programa ajusta el parámetro de descenso (ξ ) adaptándolo durante las iteraciones de equilibrio de la siguiente manera: 1. (ecuación(6. Si ξ ya es igual a 1. entonces se reduce ξ por un factor de 1/4 (estableciéndolo en 0. Si ξ = 1 (matriz secante) y el residuo ha decrecido en tres iteraciones seguidas (o 2 si ξ se incrementó a 1 durante el proceso de iteración de equilibrio por el caso arriba indicado). establece ξ = 1 y rehace la iteración usando la matriz secante. La matriz usada en la ecuación de Newton-Raphson. • Si ξ < 1 . 0 (se usa la matriz tangente). Una vez que ξ sea menor que 0. continua iterando.0156. Empieza cada substep (subpaso) usando la matriz tangente (ξ = 0 ) . UNIVERSIDAD DE SEVILLA
Si decrece (indicando una solución convergente):
.  
ξ = Parámetro de descenso. lo que es una ventaja. ESTUDIO DE LA CONVERGENCIA a cambiar a la matriz de rigidez tangente completa al converger la solución. E.0. INGENIEROS.    K T  = Matriz tangente.25)
 K S  = Matriz secante (o más estable).25) y se continúa iterando. Observa el cambio en el residuo
{R} 2 en las iteraciones de equilibrio:
Si aumenta (indicando una posible divergencia): • • • Elimina la solución actual si ξ < 1 . (Eggert [6]).14)) se define como la suma de dos matrices:
 K iT  = ξ  K S  + (1 − ξ )  K T       
(6. se toma ξ = 0.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.T. dando como resultado una tasa de convergencia más rápida.S. 2.
Búsqueda Lineal La opción de búsqueda lineal (a la que se accede con el comando LNSRCH) trata de mejorar una solución de Newton-Raphson {∆ui } escalando el vector solución por un valor escalar denominado el “Parámetro de búsqueda lineal”. Si se encuentra un mensaje de “pivote negativo” (lo cual indica un matriz mal condicionada): • • Si ξ < 1 .26)
E. que forman una matriz secante si ξ > 0 ) incluyen: plasticidad. contacto.S. en el que la matriz se recalcula en cada iteración.15) de nuevo:
{ui +1} = {ui } + {∆ui }
Las no linealidades que hacen uso del descenso adaptativo (esto es. se toma ξ = 1 y se rehace la iteración usando la matriz secante. Si ξ = 1 . INGENIEROS. ESTUDIO DE LA CONVERGENCIA 3.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. UNIVERSIDAD DE SEVILLA
. en caso de que no lo esté se termina la ejecución. El “Procedimiento de solución completa de Newton-Raphson” es también el que se usa por defecto para plasticidad. se elimina la solución actual. El descenso adaptativo se usa por defecto en estos casos a no ser que las opciones de “búsqueda lineal” o “arc-length” estén activadas.0. campos magnéticos no lineales que usan la formulación potencial escalar. se bisecta el paso de carga si el comando “automatic time stepping” está activado. el elemento concreto SOLID65 con el KEYOPT(7) = 1. y el elemento de membrana laminar SHELL41 con el KEYOPT(1) = 2.
Considerando la ecuación (6. rigidez de tensiones con grandes deformaciones. Sólo está disponible con el “Procedimiento de solución completa de Newton-Raphson”.T. contacto y no linealidades debidas a grandes deformaciones.
Un esquema de resolución iterativa basado en el “Método de Regula Falsi” (el cual.27)
s se determina automáticamente minimizando la energía del sistema.
3.28) (Schweizerhof y Wriggers [7]). si se usa la opción de la búsqueda lineal.5 g 0 . −
g s no cambie significativamente entre iteraciones. 0 .DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. ESTUDIO DE LA CONVERGENCIA En algunas situaciones de la solución.0.26) sea modificada para escribirse:
{ui +1} = {ui } + s {∆ui }
donde: s = Parámetro de búsqueda lineal.0.0 (esto es. No se permite un s menor de 0. Las iteraciones continúan hasta que:
1. Se hayan ejecutado 6 iteraciones.b]. con f(a) y f(b) de signos diferentes)se usa para resolver la ecuación (6.
(6. donde g 0 es el valor de la ecuación (6.
g s sea menor que 0. 0 . 05 < s < 1. la ecuación (6. De ahí que. donde f(x) es una función continua definida en el intervalo [a. el uso del
{∆ui } completo
inestabilidades de la solución.05.T. INGENIEROS. nos permite encontrar la raíz de la ecuación f(x)=0. usando Fnnr1 en lugar de F nr s {∆u} .
en s = 0. 0.S. lo que se
reduce a encontrar el cero de la ecuación no lineal:
g s = {∆ui }
({F } − {F ( s {∆u })})
a nr i
(6. Si g s > 0. no se realiza ninguna iteración y s se toma como 1.28)
g s = Gradiente de la energía potencial con respecto a s .
aun más importante.
Figura 6.T. Esto se ilustra en la figura 6. es que el fenómeno de pandeo no lineal es capaz de resolver una región de inestabilidad durante la trayectoria de post-pandeo. más ajustado a la realidad.
La principal característica del pandeo no lineal frente al pandeo de autovalores que calcula Ansys.2. que modelen la rigidez y la respuesta dada por la estructura. Método del Arc-Length Los análisis con detalle de estructuras geométricamente no lineales requieren la creación de modelos matemáticos que incluyan con precisión las condiciones específicas de carga y soporte.0. UNIVERSIDAD DE SEVILLA
.15).8. mientras que el pandeo de autovalores sólo incluye el comportamiento lineal de la región de pre-pandeo hasta la bifurcación (punto de carga crítica).8: ”Comportamiento de Pandeo No lineal vs Pandeo de autovalores”
6. INGENIEROS. En dicha figura también se observa que cuando se realiza el análisis no lineal. como lo llama Ansys.S. y. ESTUDIO DE LA CONVERGENCIA
La solución escalada {∆ui } se usa para actualizar los valores actuales de los GDL
{ui +1} en la ecuación (6. que se ejecutan en la siguiente ecuación de equilibrio. con un cálculo de pandeo de autovalores. el pandeo comienza a una carga crítica menor a la obtenida con un análisis lineal o.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.4.
Figura 6.S. El fenómeno de “snap-through” se ilustra muy bien en el caso del arco rebajado cargado puntualmente de la figura 6. que la llevan a una nueva situación estable. lo que indica un aumento de los desplazamientos al disminuir las cargas.0.T. como se observa en la figura 6. porque la recta tangente a la trayectoria de equilibrio es negativa. como se muestra en la figura 6. INGENIEROS. La región comprendida entre dos puntos críticos es inestable. A este fenómeno se le denomina “snap-through”.10. como es habitual (salvo en ensayos de laboratorio).10: ”Snap Through en el Pandeo”
. cuando la carga crece superando el primer punto límite.9.8. cuando la carga cae por debajo del segundo punto límite.9: “Fenómenos de snap-through y snap-back. Si el proceso de carga se hace sin control en desplazamientos. Respuesta en carga y descarga”
Al descargar esta estructura.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. la estructura sufre un “salto” brusco de desplazamientos. ESTUDIO DE LA CONVERGENCIA Las estructuras geométricamente no lineales contienen a menudo puntos límite en los cuales la trayectoria de equilibrio tiene tangente horizontal. se produce otro salto brusco de desplazamientos que se denomina “snap-back”.
El método del “arc-length” se presenta como el adecuado para predecir la respuesta correcta de estructuras con comportamientos complejos del tipo “snap-through” y seguir la trayectoria de equilibrio no lineal a través de los puntos límite. También es difícil hallar cuánta carga adicional puede soportar con seguridad una estructura bajo estas circunstancias.S. Predecir “snap-through” y “snap-back” de ciertas estructuras es difícil y resulta computacionalmente costoso. UNIVERSIDAD DE SEVILLA
. por tanto. Introducción
Las estructuras pueden experimentar condiciones de carga que causen grandes desplazamientos que cambien de forma significativa su geometría y requieran. Se usan pasos de carga iterativos para aplicar la carga a la estructura en forma de pequeñas cargas incrementales y hallar cada vez los correspondientes desplazamientos incrementales. como el Método de NewtonRaphson anteriormente presentado. Estos métodos fallan cuando tengan lugar comportamientos de “snapback” a lo largo de la trayectoria de carga. INGENIEROS. Aun existen ciertos obstáculos para encontrar respuestas no lineales de estructuras sometidas a grandes cambios geométricos.T. Esto sucede porque debido a la naturaleza de
E.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. ESTUDIO DE LA CONVERGENCIA La presencia de puntos críticos de estabilidad y trayectorias de equilibrio inestables son las principales dificultades que las soluciones numéricas deben superar para capturar completamente la respuesta no lineal.0. La representación de estos resultados define la curva de la trayectoria de equilibrio de la estructura bajo las cargas aplicadas. Las grandes deflexiones son descritas por ecuaciones diferenciales no lineales que pueden ser resueltas usando técnicas incrementales. Una estructura que sufra grandes cambios en su geometría a menudo presenta puntos límite con una respuesta inestable de “saltos” (“snap-through” y “snap-back”)
durante un colapso estático. que las ecuaciones de equilibrio se planteen en la geometría deformada. La solución a estas inestabilidades estructurales es difícil de hallar con los métodos comunes de resolución de ecuaciones. tales como el método de Newton-Raphson. En análisis no lineales la matriz de rigidez tangente sustituye a la matriz de rigidez que se usa en los análisis lineales.
11:”Divergencia en el Método de Newton-Raphson”
Los investigadores han estudiado continuamente estos problemas y se han obtenido mejoras en el proceso que han sido gradualmente introducidas en los programas comerciales de elementos finitos. El método general del “arc-length” usado para análisis estructurales. la carga puede cambiar de dirección al cambiar las deflexiones cuando se atraviesa un máximo local. . Ver figura 6.Y puntos límite de control que definen un “snap-back” local. El análisis geométricamente no lineal llevado a cabo mediante elementos finitos debe ser capaz de hallar todas las posibles respuestas durante la aplicación de grandes cargas. INGENIEROS.
Pasar puntos críticos durante la respuesta geométricamente no lineal es retador.11.Puntos límite de carga que se alcanzan cuando la trayectoria de la respuesta tiene un “snap-through” local. no puede haber convergencia si la matriz de rigidez tangente (la pendiente de la curva fuerza-desplazamiento en cualquier punto) se hace cero.T.S. la ecuación (6. ESTUDIO DE LA CONVERGENCIA la ecuación empleada en el método iterativo de Newton-Raphson.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. En uno de estos puntos límite de control.
Figura 6.0. es decir. Dos tipos de puntos críticos hallados en este tipo de comportamiento son:
. Existe una familia importante de métodos de resolución de ecuaciones no lineales llamados genéricamente métodos del “arc-length” (longitud del arco) que pueden superar estos puntos críticos.
usada por Ansys. para resolver ecuaciones no es la de Riks y Wempner sino la variación del método introducida por Forde y Stiemer [11] en 1987.12. y la forma de Forde y Stiemer. 1979) y Wempner (1971) [10]. dado que no es el objetivo de este proyecto hacer un análisis demasiado exhaustivo de la gran cantidad de métodos del “arc-length” de resolución de ecuaciones no lineales. El método pasa a hallar el nuevo punto de equilibrio desde el punto inicial.0. Se llega a la terminación del método cuando se cumplen las condiciones de máximo factor proporcional de carga o máximo desplazamiento nodal dados. De hecho la forma del método del arc-length que usa el programa Ansys 8. En la figura 6. ESTUDIO DE LA CONVERGENCIA surge como una variación del método general de Newton-Raphson ideada para superar las dificultades de éste para pasar por puntos críticos.0.S. La técnica se asemeja mucho al método de Newton-Raphson descrito en Riks [8] y [9] excepto en que en este caso el incremento de carga aplicado pasa a ser una incógnita adicional en el problema. (Para mas detalle consultar Crisfield [12]. Owen [13] y Riks [8] y [9]).DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. El proceso también se acaba cuando se alcanza el máximo número de incrementos de carga en un paso determinado.T. por ser la primera históricamente y base de las demás.12 se muestra la curva carga-desplazamiento para un sistema de un solo grado de libertad. i. El vector tangente a la curva en i puede escribirse como:
E. Por este motivo y. El método del arc-length de Riks-Wempner traza la trayectoria no lineal de equilibrio usando un proceso iterativo que empieza computando los desplazamientos iniciales provocados por un incremento de carga dado. No son pocas las variaciones que se desarrollaron sobre el trabajo original de Riks y Wempner. se presentarán sólo dos: el planteamiento de Riks y Wempner. como se puede apreciar en la figura 6. Fué desarrollado originalmente por Riks (1972. Análisis estático geométricamente no lineal mediante el método del “ArcLength” de Riks-Wempner
En el método de Riks-Wempner la longitud de un vector tangente a la trayectoria de equilibrio se usa para hallar un nuevo punto que es intersección del plano normal a la tangente. Una carga dada por el analista estimará las magnitudes del incremento inicial de carga para cada paso. INGENIEROS. UNIVERSIDAD DE SEVILLA
T. ESTUDIO DE LA CONVERGENCIA
" "  ∆ui  ti =    ∆λi 
(6.S.29)
donde ∆λi es la carga incremental aplicada en i.12 y se puede expresar como:
" "  ∆ui  ni =    −∆λi 
Figura 6. ni también se muestra en la figura 6.12: “Método del arc-length de Riks-Wempner”
. INGENIEROS. KTi es la matriz de rigidez tangente evaluada en i y ∆ui es el vector incremental de desplazamiento hallado resolviendo
" " KTi ∆ui = ∆λi F "
(6.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.0.30)
El vector normal.
.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. Para hallar los desplazamientos iniciales. ∆u0 .T.0. ∆λ0 . se usa para computar el primer vector desplazamiento.13.13: “Inicio del método del arc-length de Riks-Wempner”
Al comenzar el método del “arc-length” de Riks-Wempner un incremento inicial de carga. se usan triángulos semejantes. INGENIEROS. El incremento de carga es un parámetro dado por el analista que divide la carga total aplicada en incrementos iguales.34)
E. a lo largo de la trayectoria de equilibrio se puede calcular como
" " "T " ∆s0 = t0 ⋅ t0 = ∆λ 2 + ∆u0 ∆u
(6. ∆u0 .S.32)
" " KT0 ∆utot = λ F "
(6. Estas variables se muestran en la representación de carga – desplazamiento de la figura 6. Un incremento de carga dado empieza el proceso y a partir de él se hallan los desplazamientos ∆u0 usando la matriz de rigidez tangente KT0 . Durante este incremento inicial la rigidez tangente es la misma que la rigidez lineal. y la longitud del primer vector tangente t0 . Los desplazamientos iniciales ∆u0 se hallan resolviendo
∆λ0 λ " = " ∆u0 ∆utot
donde λ = 1 y ∆utot se deduce de la expresión
(6. ESTUDIO DE LA CONVERGENCIA
Figura 6. que será la que determine el radio del arclength.33)
La longitud del vector tangente t0 .
Durante la solución. Este método usa iteraciones esféricas explícitas para mantener la ortogonalidad entre el radio del arc-length y las direcciones ortogonales.0. es necesario limitar el radio del arc-length (el tamaño del incremento de carga del arc-length) usando el radio inicial o de referencia (con el comando NSUBST) y el rango de variación permitido. I des .e.
Al ser este un método “quasi-Newton” todas las opciones en el programa del método de Newton-Raphson son aún la base de la solución del arc-length.
Método del arc-length usado por Ansys El método del arc-length usado por Ansys (que se activa con el comando ARCLEN.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. el método del arc-length es él mismo un método que automatiza los tamaños de los pasos de carga (No necesitará el comando AUTOTS.T. Se asume que todas las magnitudes de carga están controladas por un único parámetro escalar (i. el método del arc-length variará el radio del mismo en cada substep (subpaso) de acuerdo con el grado de no linealidades que estén incluidas.ON).S. el factor de carga total).. Para problemas con giros bruscos en la curva carga-desplazamiento o materiales dependientes de la trayectoria de carga. UNIVERSIDAD DE SEVILLA
Ya que ahora.35)
El analista decide sobre el número de iteraciones requeridas.
E.ON) es apropiado para soluciones no lineales estáticas de problemas inestables. ESTUDIO DE LA CONVERGENCIA A través del resto de iteraciones la longitud de arco (arc-length) es constante o puede ser escalada por la introducción de un parámetro por parte del analista en la siguiente expresión
 I 2 ∆si = ∆si −1  des   I i −1 
(6. los vectores de desplazamiento y el factor escalar de carga se tratan como incógnitas. como se describe en Forde y Stiemer (1987). I i −1 . y sobre el número de iteraciones deseadas. INGENIEROS.
S. UNIVERSIDAD DE SEVILLA
. se replantea asociada al factor de carga total λ :
 K iT  {∆ui } = λ { F a } − { Fi nr }  
(6. Es interesante recalcar que el factor de carga λ sólo multiplica a las fuerzas aplicadas. INGENIEROS. Si escribimos el factor proporcional de carga λ de forma incremental.36)
donde λ está comprendido normalmente en el rango −1.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.14).
En el procedimiento de resolución del método del arc-length usado por Ansys.14):
 K iT  {∆ui } − ∆λ { F a } = ( λn + λi ) { F a } − { Fi nr } = − { Ri }  
∆λ = Factor de carga incremental (tal y como se muestra en la figura 6. ya que el otro término de fuerzas son las fuerzas internas que se calculan a partir de las aplicadas. tenemos en el substep n y la iteración i (ver figura 6. la ecuación no lineal (6. 0 .T.14: ”Aproximación mediante el arc-length de Forde-Stiemer con el método completo de Newton-Raphson”
E.0. ESTUDIO DE LA CONVERGENCIA El rango de variación del radio del arc-length se limita por los multiplicadores máximo y mínimo del radio de referencia del arc-length (MAXARC y MINARC en el comando ARCLEN).
Figura 6. 0 ≤ λ ≤ 1.14).
{∆u } = Desplazamiento incremental proveniente del método convencional de
∆un = Suma de todos los desplazamientos incrementales ∆ui de esta iteración.
Estos se definen:
{∆u } =  K  {F }   {∆u } = −  K  {R }  
(6. al haber introducido una nueva incógnita en el
' i '' i
problema.39) y (6.38)
{∆u } = Desplazamiento debido a un factor de carga unidad.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.39) (6. en la iteración i (ver figuras 6. el factor de carga incremental ∆λ . es necesario usar las ecuaciones (6.37):
El desplazamiento incremental {∆ui } puede escribirse en dos partes siguiendo la
{∆ui } = ∆λ {∆ui' } + {∆ui'' }
(6.40) para resolver
{∆u } y {∆u } .0. por ejemplo.15)(es equivalente a la ecuación (6.S. INGENIEROS.40)
En cada iteración del arc-length. Además. UNIVERSIDAD DE SEVILLA
. el método introduce una nueva restricción para cada substep que se conoce como la “ecuación del arc-length” que puede ser escrita.T.
E. ESTUDIO DE LA CONVERGENCIA ecuación(6.14 y 6.34) del método de Riks-Wempner)
# 2 = λi2 + β 2 {∆un } {∆un } i
β = Factor de escala (con unidades de desplazamiento-1) usado para asegurar
el correcto escalamiento en las ecuaciones (ya que λ es adimensional).
Por lo tanto el método del arc-length de Forde-Stiemer que usa Ansys permite a la carga y al desplazamiento que varíen durante un substep o subpaso de carga tal y como se muestra en la figura 6. durante las iteraciones de un mismo substep.15.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.16.15: ”Método del arc-length de Forde-Stiemer.T.
Figura 6. Un substep”
El radio del arc-length # i se fuerza. se impone que se mantenga constante a lo largo de un substep dado..S.
E.0. ESTUDIO DE LA CONVERGENCIA Así el método de Forde-Stiemer usado por Ansys puede verse más claramente en la figura 6.
# i = # i −1 = . es decir.. UNIVERSIDAD DE SEVILLA
. a ser idéntico al radio #1 de la primera iteración de ese substep. INGENIEROS.
0 0 MINARC ⋅ #1 ≤ # i ≤ MAXARC ⋅ #1
(6. Sin embargo. ellas solas. hay muchas formas de resolver ∆λ aproximadamente. Para el método del arc-length. determinan.0.T. no se usa el “Automatic time stepping”.S. INGENIEROS. el vector solución ( ∆ui .43)
Las ecuaciones (6. ∆λ ) . el residuo
Figura 6.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. Los radios del arclength en los substeps subsiguientes del mismo paso de carga o load step son calculados por el programa a partir del radio del arc-length del substep previo y el comportamiento de la solución y teniendo en cuenta que estos radios deben estar comprendidos dentro del rango límite que determinan el radio de referencia para ese paso de carga y los multiplicadores del mismo (dados en las etiquetas MAXARC y MINARC del comando ARCLEN).16: ”Comportamiento de la convergencia en Ansys usando el arc-length”
0 El radio del arc-length #1 en la primera iteración del primer substep de un paso de
carga o load step se denomina “radio inicial o de referencia del arc-length” (definido por el comando NSUBST ya que este radio de referencia se calcula como el cociente entre la carga total aplicada en el paso de carga correspondiente y el número de substeps especificados en dicho comando para ese mismo paso de carga).
El método de la iteración esférica explícita se usa para asegurar la ortogonalidad (Forde y Stiemer (pág.38) y (6.41) conjuntamente. por tanto. UNIVERSIDAD DE SEVILLA
.174 del libro indicado en la bibliografía)). En este método.
44) (6. INGENIEROS.46)
Esta expresión puede ser simplificada para casos particulares de ortogonalidad. se calcula de antemano como producto escalar de los vectores normal ni y tangencial ti .45)
Por consiguiente.15):
{ui +1} = {un } + {∆un } + {∆ui }
(6. el factor de carga incremental del arc-length se determina mediante la fórmula:
ri − {∆un }
{∆u } λ + {∆u } {∆u }
T '' i T i n ' i
(6. El método funciona bien incluso cuando se encuentra en las proximidades de puntos críticos en los que hay cambios bruscos en la solución.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. UNIVERSIDAD DE SEVILLA
λi +1 = λn + λi + ∆λ
(6. ESTUDIO DE LA CONVERGENCIA requerido ri (un escalar) para la iteración explícita en una esfera.14 y 6.S. Los autores sugieren la selección de un vector normal ni con una dirección arbitraria con respecto al vector ti que es tangente a la configuración de carga-desplazamiento en un punto dado. Los vectores normal y tangencial constan de m componentes provenientes del vector desplazamiento y una proveniente del factor de carga.
E. reseñar que los vectores solución se actualizan de acuerdo con (ver figuras 6.0. Estas componentes se combinan usando el factor escalar β formando vectores de m+1 componentes que pueden escribirse como:
" " ui  ti =    βλi  " "  ∆u  ni =    β ∆λ 
n = Número del substep actual.T. Finalmente.
DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. que orientan sobre cómo realizar un análisis no lineal en Ansys. donde # 0 es el radio inicial o de referencia del arc-length en i 0 #i
el substep i definido en el comando NSUBST. En el caso en el que las cargas aplicadas sean mayores o menores que las máximas o mínimas cargas críticas.T.41). convendría enumerar los puntos más importantes que se han presentado en este capítulo.3. UNIVERSIDAD DE SEVILLA
E. y que se relaciona con el factor de carga y los desplazamientos a través de la ecuación (6.S. provenientes de bibliografía diversa. Resumen
Una vez presentada la base teórica de los métodos de resolución de ecuaciones que usa Ansys y visto con más detalle el Método del Arc-length. analistas y de la propia ayuda del programa. respectivamente.0. Es recomendable. haciendo hincapié sobre todo en aquellos que resulten más prácticos desde el punto de vista operativo en la resolución del problema no lineal planteado. ESTUDIO DE LA CONVERGENCIA Los valores de λn y ∆λ están disponibles en el módulo POST26 (comando SOLU) y se corresponden con las etiquetas ALLF y ALDLF. el arc-length continuará con las iteraciones en ciclos ya que λ no se acerca a la unidad. por tanto. INGENIEROS. • En análisis no lineales es muy recomendable aplicar la carga en forma de incrementos (pasos de carga o load steps) para poder seguir bien la trayectoria de la curva carga-desplazamiento.
En este pequeño resumen. El radio normalizado del arc-length se encuentra en la etiqueta ARCL (comando SOLU) y se corresponde con el valor
#i . se incluyen una serie de consejos prácticos. 6. acabar con las iteraciones del arc-length imponiendo un criterio de terminación (puede hacerse usando los comandos ARCTRM o NCNV).
14 y 6.T.16). INGENIEROS. UNIVERSIDAD DE SEVILLA
.ON.17. . ESTUDIO DE LA CONVERGENCIA • El método de resolución numérica de ecuaciones más usado por su sencillez y antigüedad es el “Método iterativo de Newton-Raphson” base de casi todos los demás y
T que se describe en la ecuación (6.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. . A su vez cada paso de carga se resuelve a través de un número de subpasos de carga o substeps que se determinan a través E.17:”Newton-Raphson vs Arc-length”.i  
• En análisis no lineales en los que se esperen comportamientos de “snap through” y “snap back” o en los que se quiera trazar la trayectoria de la solución en la región de post-pandeo. el método de Newton-Raphson antes mencionado presenta problemas debido a la existencia de puntos críticos de pendiente horizontal en los que se anula la matriz de rigidez.S.i  {∆ui } = Fna − Fnnr .
• El método del arc-length que usa Ansys es la forma de Forde y Stiemer del mismo (puede observarse en las figuras 6. • El procedimiento básico de resolución del problema usando el método del arclength en Ansys puede resumirse en los siguientes pasos:
Figura 6. En estos casos es aconsejable usar el Método del Arc-length que puede pasar bien esos puntos críticos.16):  K n . Una comparativa ilustrativa de ambos métodos puede observarse en la figura 6.La carga se aplica a la estructura incrementalmente a través de una serie de pasos de carga o load steps siendo la carga aplicada en cada uno de ellos un valor definido por el analista.El método se activa con el comando ARCLEN.0.
se resuelve de la manera indicada hasta que:
La solución converja.50) (6.
.Los radios de los siguientes substeps. Los substeps se resuelven mediante una serie de iteraciones que siguen hasta que: la solución converja en el substep determinado (ya que se impone que la solución esté en equilibrio en cada substep) o bien se sobrepase sin converger el número máximo de iteraciones permitidas en un substep dado por el comando NEQIT.S. Ahora bien. si no que Ansys los calcula basándose en el radio del arc-length del substep previo y el comportamiento no lineal de la solución.49)
E.51)
0 donde # 1 es el radio de referencia del arc-length (el del primer substep) para
ese paso de carga. el analista determina el rango en el que se moverán estos valores imponiendo límites de variación al radio del arc-length de referencia (que por tanto sirve como referencia para el tamaño del resto de los radios del mismo paso de carga) a través de las etiquetas MAXARC y MINARC del comando ARCLEN.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. UNIVERSIDAD DE SEVILLA
.El “radio del arc-length de referencia” es el radio del arc-length de la primera iteración del primer substep de un paso de carga (en realidad es el radio de todo el primer substep porque Ansys impone que el radio sea constante en las iteraciones de un mismo substep). ya no los define directamente el analista.T. Este radio lo define el analista mediante el comando NSUBST y la carga que aplica en el “loadstep” ya que se calcula con la expresión:
#0 = i
Fracción de la carga total aplicada en el " LS " Número de substeps del " LS " definidos en NSUBST
(6. Los límites de variación del radio del arc-length en los substeps de un loadstep determinado son:
0 Límite inferior = MINARC * #1 0 Límite superior = MAXARC * #1
(6.Cada paso de carga o loadstep. ESTUDIO DE LA CONVERGENCIA del comando NSUBST.0. INGENIEROS.
. INGENIEROS. La bisección continuará hasta que se obtenga una solución convergida o se alcance el límite inferior del arc-length en ese LS). cuando el objetivo es determinar la carga crítica de pandeo de una estructura. ya que hay que elegir los parámetros que lo gobiernan con mucho cuidado. (Si la solución no consigue converger en el número de iteraciones de equilibrio indicadas en NEQIT.T. • El arc-length es un método muy potente para hacer frente a los análisis no lineales mediante elementos finitos en el cálculo de estructuras. se debe hacer una buena estimación del radio de referencia que se usará para resolver un paso de carga llegando a un compromiso entre tiempo de ejecución y precisión.17. se aconseja realizar un análisis lineal previo de pandeo de autovalores antes de realizar el análisis de pandeo no lineal. aunque también pueden obtenerse cargas críticas con el método del Arclength ajustando bien sus parámetros. Sin embargo con otros métodos como el de Newton-Raphson la carga sí es perfectamente conocida. (que puede encontrarse en la etiqueta ALLF del comando SOLU) puede ser tanto positivo como negativo. la carga aplicada en un substep determinado diferente del primero en la resolución de un paso de carga no es fácil de conocer de forma exacta.DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8. es mejor basar el criterio de convergencia en fuerzas [CNVTOL. • Cuando se usa el arc-length. que en un análisis del arc-length se relaciona con el factor
E.0.S. ya que el valor de la carga aplicada varía a lo largo del arco esférico. siempre que sea posible se recomienda que se obtengan mediante otros métodos como el de Newton-Raphson cuya precisión será mayor. Esto puede apreciarse por ejemplo en la figura 6. De manera análoga. pero su correcta aplicación requiere cierta pericia. De ahí que. el valor de TIME en la solución. ESTUDIO DE LA CONVERGENCIA Se alcance alguna de las condiciones de salida del análisis del programa que pueden especificarse en muchos comandos pero sobre todo en: ARCTRM.F] y no en desplazamientos [CNVTOL. NEQIT o NCNV. • Según se definen los radios de las diferentes iteraciones del método. el programa bisectará automáticamente el radio del arc-length y continuará con el análisis. λ .U]. Por todo esto. • El factor de carga total del arc-length.
Los valores negativos de ALLF o TIME se encuentran de forma habitual cuando se dan comportamientos de “snap-through” en la estructura.0. simplemente que el método del arc-length está aplicando la carga en la dirección contraria a la definida en el problema para mantener la estabilidad en la estructura. se corre el peligro de que el arc-length vuelva atrás siguiendo la misma trayectoria de carga. también puede ser positivo o negativo.T.
E. ESTUDIO DE LA CONVERGENCIA de carga total. λ . produciéndose el fenómeno denominado “drifting back”. UNIVERSIDAD DE SEVILLA
. INGENIEROS. pero en sentido inverso.S. La interpretación de valores negativos en ALLF o TIME indica. Cuando esto ocurre. El fenómeno de “drifting back” suele deberse a una mala elección del tamaño del radio del arc-length (radios demasiado grandes o demasiado pequeños).DISEÑO DE UN SILO CONFORME AL EUROCÓDIGO CAPÍTULO 6: MÉTODOS NUMÉRICOS USADOS POR ANSYS 8.
AboutBrowse booksSite directoryAbout ScribdMeet the teamOur blogJoin our team!Contact UsPartnersPublishersDevelopers / APILegalTermsPrivacyCopyrightSupportHelpFAQAccessibilityPressPurchase helpAdChoicesMembershipsJoin todayInvite FriendsGiftsCopyright © 2016 Scribd Inc. .Terms of service.Accessibility.Privacy.Mobile Site.Site Language: English中文EspañolالعربيةPortuguês日本語DeutschFrançaisTurkceРусский языкTiếng việtJęzyk polskiBahasa indonesia7.CAPÍTULO 6 by camello375 viewsEmbedDownloadRead on Scribd mobile: iPhone, iPad and Android.Copyright: Attribution Non-Commercial (BY-NC)List price: $0.00Download as PDF, TXT or read online from ScribdFlag for inappropriate contentMore informationShow less
More from camelloarches elementsby camelloLa Estructura de Las Revoluciones Cientificasby camelloFuerza de Un Cilindroby camelloTop Nonfiction on ScribdThe Innovators: How a Group of Hackers, Geniuses, and Geeks Created the Digital Revolutionby Walter IsaacsonThe Emperor of All Maladies: A Biography of Cancerby Siddhartha MukherjeeHard Choicesby Hillary Rodham ClintonTop Fiction on ScribdThe Light Between Oceans: A Novelby M. L. StedmanBrooklyn: A Novelby Colm TóibínThe Rosie Project: A Novelby Graeme Simsion

References: resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución