Simulación numérica del conformado metálico. Eficiencia numérica en procesos de laminación. Parte II

Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería. Vol. 13,3, 297-320(1997) PROCESOS DE LAMINACIÓN. PARTE 11 M.A. GUTIÉRREZ* M. OJANGUREN* Y J.J. ANZA** * Departamento de Mecánica y CAD-CAM, LABEIN Cuesta de Olabeaga 16, 48013 Bilbao, España Tel.: + 34-4-489 2400 Fax: + 34-4-411 1749 ** Departamento de Matemática Aplicada Universidad del PazS Vasco Alameda de Urquijo s/n, 48013 Bilbao, España Tel.: + 34-4-427 8055 Fax: + 34-4-441 4041 E-mail: [email protected] RESUMEN La resolución aproximada mediante técnicas de elementos finitos de la simulación del conformado requiere un gran esfuerzo computacional contradictorio con la interactividad requerida en la aplicación industrial. En la segunda parte de este trabajo se analiza la eficiencia computacional de los algoritmos combinados de la elastoplasticidad, deformación finita y contacto, incidiendo especialmente en la resolución óptima del sistema de ecuaciones lineales que es preciso afrontar numerosas veces dentro del esquema incremental iterativo propio del análisis no lineal. Se presenta en esta línea un método basado en la descomposición de dominios cuya aplicación a problemas de contacto y en especial a la simulación de proceso de laminación resulta muy ventajosa. NUMERICAL SIMULATION OF METAL FORMING. NUMERICAL EFFICIENCY IN ROLLING PROCESSES. PART 11 SUMMARY The numerical simulation of metal forming processes follows a highly non linear analysis where general aspects as elastoplasticity, finite deformation and contact mechanics are combined. Approximated solutions obtained by finite element techniques require strong computational effort, that contradicts the need of interactivity for industrial applications. In this work, computational efficiency is analysed, especially the optimum solution of the linear systems to be solved through the incremental-iterative schemes associated with non linear analysis. A domain decomposition method, especially adapted to contact problems is presented, as well as the improved performance obtained in the application to hot rolling simulation. Recibido: Diciembre 1996 QUniversitat Polithcnica de Catalunya (España) ISSN 0213-1315 En el marco del análisis no lineal del sólido deformable discretizado mediante elementos finitos, el estudio de la eficiencia del proceso incremental-iterativo de resolución implícita, admite dos planteamientos o perspectivas. En la primera se consideran los aspectos relacionados con la formulación teórica y su aproximación, que afectan al coste del cálculo elemental y pueden condicionar el tamaño de cada incremento y el número de iteraciones por incremento1. En la segunda se considera el proceso iterativo de resolución de un incremento; dada una formulación y su aproximación, la eficiencia global del proceso depende de la estrategia de actualización y factorización de la rigidez tangente en cada iteración4. El primer planteamiento afecta al cálculo de los vectores de fuerzas internas y de las matrices de rigidez tangente. El reciente desarrollo de los procedimientos implícitos ha sido importante en este sentido y en la actualidad existen métodos robustos que permiten alargar los incrementos, manteniendo su resolución en pocas iteraciones. Ello ha sido posible a partir de la mejora de la integración de las ecuaciones constitutivas de la elastoplasticidad en deformación finita mediante algoritmos de retorno y formulaciones corrotacionales o hiperelástica~~~~, que solventan la integración objetiva de las ecuaciones planteadas como relación de tasas temporales y permiten establecer rigideces tangentes consistentes que se conservan la convergencia cuadrática del método de Newton5-?. El segundo planteamiento se refiere a la sucesiva actualización y resolucián de los grandes sistemas lineales de ecuaciones asociados con la aproximación mediante elementos finitos. En el campo de la resolución de sistemas lineales de ecuaciones. el esfuerzo investigador trasciende el ámbito de aproximación numérica del sólido e incluso del continu~~~~~ tanto en relación con la computación secuencia1 o paralela como con respecto a los métodos directos o iterativos. En la línea del segundo planteamiento y como continuación de la parte primera1, en esta segunda parte del trabajo se profundiza en la eficiencia computacional de procesos de conformado con grandes deformaciones elastoplásticas y contacto. En el contexto de las estructuras numéricas asociadas con el tratamiento de contacto mecánico mediante multiplicadores de Lagrange, los desarrollos que se proponen están orientados al tratamiento diferenciado de todas aquellas partes físicas y de interaccióri que lo requieran, procediéndose a la actualización y la factorización del sistema de ecuaciones en bloques, mediante los métodos de resolución directos o iterativss más adecuados para cada una de ellos. Los resultados obtenidos en una estación Indigo2 R3000 se discuten en los diversos ejemplos presentados, en los que el comportamiento elastoplástico está calculado en formato hipoelástico corrotacional con integración en la configuración intermediaL. CONTACTO Resolución del sistema lineal de ecuaciones El sistema de ecuaciones lineal que surge en un problema de contacto sin deslizamiento entre dos cuerpos deformables mediante multiplicadores de Lagrange tiene la forma' donde K1 y K2 son las matrices de rigidez de cada uno de los sólidos y C1 y C2 son las matrices de acomplamiento que introducen la restricción cinemática del contacto. El sitema es no singular, y por tanto tiene solución única, pero la matriz de coeficientes no es definida positiva y no es posible aplicar directamente el método del gradiente conjugado. Al estar los multiplicadores concentrados en las últimas ecuaciones, el complemento de Schur es definido negativo y el sistema puede resolverse mediante una descomposición LDL~. Si las matrices de rigidez de cada uno de los cuerpos en contacto son invertibles, eliminando los desplazamientos u1 y u2, el sistma anterior puede expresarse en función de los multiplicadores de Lagrange, obteniéndose el siguiente condensado en la interfase de contacto LA = b donde L = C~KT~CY + ~2~2~~2 y b = ClKl1f1+ CZK~~~Z - d (2) La matrix L es simétrica, definida positiva y llena. Su dimensión coincide con el número de grados de libertad en contacto, y por tanto es en general menor que la del sistema original. Existen varias formas de resolver este problema, dependiendo del método elegido para obtener KT~ y Kyl y para resolver el sistema condensado. Resuelto éste, los desplazamientos pueden recuperarse, ya que El sistema condensado (2) puede resolverse con un método directo, pero para ello es necesario calcular de forma explícita la matrix de coeficientes L. El método del gradiante conjugado16 (GC), resuelve el sistema realizando solamente productos de la inatriz L por las direcciones p de búsqueda y esta operación puede realizarse sin calcular de forma explícita la matriz L. Una vez factorizadas cada una de las submatrices de rigidez, el producto Lp se realiza de forma eficiente mediante multiplicación de matrices dispersas por un vector y sustituciones progresivas y regresivas. Para cada uno de los cuerpos en contacto se procede de la forma siguiente: 1) Multiplicación del vector p por c:; 2) Sustitución progresiva y regresiva en la factorización de Kj; 3) Multiplicación del vector resultante por Cj. Del mismo modo, los desplazamientos (3) pueden obtenerse una vez calculados los multiplicadores de Lagrange mediante sustituciones progresivas y regresivas en la factorización de las submatrices de rigidez. Para la factorización de las submatrices de rigidez puede elegirse cualquier método directo. En este trabajo se ha utilizado el método clásico de eliminación gaussiana basado en columnas activas (ACTCOL1'). Denominaremos mixto al método que divide la matriz de rigidez en los bloques inducidos por el contacto con multiplicadores de Lagrange y combina un método directo (ACTCOL) para la resolución de las submatrices de los cuerpos en contacto y un método iterativo (GC) para resolver la interacción entre ellos. Las ventajas proporcionadas por este método frente a la resolución directa del sistema de ecuaciones están relacionadas con la reducción del almacenamiento y del tiempo de cálculo. Los escasos elementos no nulos de las matrices de contacto Cj afectan negativamente al ancho de banda de matriz global y al formato standard skyline que implica el almacenamiento de un gran número de ceros. En el método mixto la matrix de rigidez se trata por bloques que se almacenan de forma independiente pudiéndose así emplear un modelo skyline para las submatrices de rigidez de los sólidos y un almacenamiento disperso para las matrices de contacto. Los sistemas que van a triangularizarse ahora con métodos directos para obtener K;' tendrán tamaño y ancho de banda menores que el problema original y el esfuerzo necesario para la factorización de todas las matrices locales será sustancialmente menor que el empleado para la matriz completa. Para que el método mixto proporcione una solución más eficiente en su conjunto, es necesario que el orden de convergencia del método del gradiante conjugado en la interfase sea suficientemente bueno. El número de iteraciones empleadas por el método del gradiente conjugado depende en gran medida del número de condición de la matriz del sistma a resolver y es como mucho igual al número de ecuaciones. La matriz del sistema condensado es menor y está mejor condicionada que la matriz original, por ello GC resulta un método efectivo de resolución incluso sin precondi~ionador~~. Otra ventaja adicional se tiene cuando el tratamiento del contacto con fricción conduce a una matriz no simétrica, cuya resolución directa presenta los inconvenientes de aumento de almacenamiento y tiempo de cálculo que ello conlleva. Con el método mixto es posible mantener la simetría de las matrices locales y diferenciar tan sólo la parte de las restricciones, siendo el aumento de almacenamiento mucho menor, ya que estas matrices se almacenan de forma dispersa. Será necesario utilizar un método de resolución para matrices no simétricas tan sólo en el sistema condensado en la interfase. Tratamiento de los sistemas singulares. Gradiante conjugado proyectado Aunque el problema global sea estable, es decir, no tenga movimientos de sólido rígido, puede ocurrir que algunos de los cuerpos componentes sólo lo sea como consecuencia de las restricciones de contacto. En ese caso, las matrices locales de estos cuerpos son singulares y su factorización requiere un tratamiento especial. Sean n cuerpos en contacto y supongamos que uno de ellos denominado genéricamente (j) tenga impedidos los movimientos de sólido rígido sólo a través de las restricciones de contacto. Entonces, en la aplicación del método mixto descrito en el punto anterior es necesario resolver un sistema singular de la forma SIMULACIÓN NUMÉRICA DEL CONFORMADO METÁLICO 301 Como la matriz Kj es simétrica, el sistema anterior tendrá al menos una solución, si y sólo si el término independiente no tiene componentes en el espacio nulo de Kj, es decir, fj - C;X es ortogonal a una base Rj del espacio nulo de Kj La condición (5) de ortogonalidad expresa que los elementos de la base Rj represetan desplazamientos del cuerpo j que no aportan energía de deformación sobre las fuerzas externas que actúan sobre el mismo, condición que es satisfecha por los modos de movimiento de sólido rígido del cuerpo j, que constituyen una base del espacio nulo de Kj . Existe por tanto una matriz K: pseudoinversa, verificando K~K~K~ = Kj de forma que la solución de (4) viene dada por donde a es un vector real de dimensión igual al número de movimientos de sólido rígido del cuerpo j. Si el cuerpo j tiene n grados de libertad y r movimientos de sólido rígido, la matriz Kj tiene rango p = n - r y puede reordenarse de forma que las singularidades aparezcan al final La matrix K: es una pseudoinversa de Kj y la matriz Rj es una base del espacio nulo de Kj. De esta forma si se tiene un sistema de contacto con nf cuerpos flotantes, la sustitución de los desplazamientos (6) en la ecuación de los multiplicadores (2) introduce como incógnita adicional al vector a. Añadiendo la condición de ortogonalidad (5), el problema quedará definido, obteniéndose por el método mixto un sistema condensado en la interfase, simétrico y no singular, que admite una única solución (A, a) y que determina unívocamente los desplazamientos de cada uno de los n cuerpos del sistema. El sistema condensado en la interfase es F -G X [-GT 0 ] [a] = [be] donde ."J F=CCjK:C) G=(CIR i...CnfRnf) K:=K-' 3 si no es flotante j=l n. r (8) Pero la matriz del sistema (7) ahora es indefinida por lo que el método de gradiante 302 M.A. GUTIÉRREZ, M. OJANGUREN Y J.J. ANZA conjugado no puede aplicarse de forma directa. El sistema anterior es equivalente al problema de minimización con restricciones algebraicas 1 min #(A) = -AT~A + ATb sujeto a GTA = e 2 Como F es una matriz simétrica y defininido positiva, el método de gradiante conjugado puede emplearse para obtener una solución del problema sin restricciones. Este algoritmo convergerá a la solución restringida si y sólo si se modifica de forma que la restricción se verifique en cada iteración18~lg. Esto se consigue proyectando todas las direcciones de búsqueda sk en el espacio nulo de GT. De esta forma, si Ak-l verifica la restricción, entonces Ak = Ak-i + PPsk también la verifica. La única diferencia que se introduce en el método de gradiente conjugado es que los residuos deben ser proyectados antes de utilizarlos en el cálculo de la dirección de búsqueda. La elección de A. está limitada únicamente por la restricción G~X~ = e siendo una opción habitual: A. = G(GTG)-'e. Es importante destacar que sólo se realiza una proyección en cada iteración. Esta proyección es relativamente barata cuando se utiliza un proyector como P = 1 - G(GTG)-lGT, debido a la pequeña dimensión de la matriz G, que coincide con el número de movimientos de sólido rígido. Antes de empezar a iterar se construye y factoriza GTG; cada proyección, así como el cálculo Ao, suponen una sustitución progresiva y regresiva en dicha factorización. Una vez encontrado el vector A, el vector a que contiene los coeficientes de los modos de sólido rígido se calcula según a = (GT G)-l GT (FA - d) . Ejemplo 1 Se trata de la indentación de un bloque elástico en una base también elástica con fricción infinita. La geometría y condiciones de contorno (desplazamiento prescrito d) del problema junto con los valores elásticos (relativos) empleados se muestran en la Figura 1. Los multiplicadores de Lagrange se definen sobre la superficie de contacto de la base. d=OL , E=l.c7 r=ü.3J SI A ! E=l e6 i =O 35 l 23 v LIDAAAAAA 2.3 I'igura 1. Indentación con bloque El objetivo de este ejemplo es mostrar el comportamiento del método de resolucióc mixto frente al global. El problema es lineal por tratarse de materiales elásticos con deformación infinitesimal y superficie de contaco conocida. Los desplazamientos máximos obtenidos con el método directo y el método mixto para distintas tolerancias del gradiente conjugado se resumen en la Tabla 1. Tabla 1. Indentación con bloque 2D: comparación de resultados Método de resolución ACTCOL Mixto CG tol = 1.e-5 Mixto CG tol = 1.e-4 Mixto CG tol = 1.e-3 Como se puede observar, los resultados coinciden hasta el sexto dígito significativo para tolerancias de 0,00001 y 0,0001, lo que supone 8 y 7 iteraciones de CG: respectivamente. Para 0,001 el CG sólo realiza 5 iteraciones y el método pierde precisión. En adelante utilizaremos un valor de la tolerancia de 0,0001, Los resultados de tiempo CPU correspondientes a una malla uniforme de 20 x 8 elementos en el bloque y distintas discretizaciones uniformes de la base se recogen en la Tabla 11, donde se indica el número de incógnitas adicionales (A) de contacto. Desplazamiento max. en Y 0,0162096798 0,0162096445 0,0162096987 0,0162102129 Tabla 11. Indentación con bloque 2D 1 Desplazamiento max. en 2 -0,0826635224 -0,0826634461 1 1 -0,0826633574 -0,0826506147 Malla 50 x 30 70 x 50 90 x 70 Se observa la reducción del tiempo empleado en la resolución del sistema, obteniéndose una mejora global que aumenta con el número de grados de libertad (4 %, 9 %y 12 %). En tres dimensiones las ventajas computacionales obtenidas en el mismo ejemplo son sensiblemente mayores, apreciándose una mejora aún más importante con el aumento de los grados de libertad del problema. Se recogen en la Tabla 111 los resultados obtenidos con distintas discretizaciones, manteniendo la relación entre el número de elementos de la base y del bloque. No X 22 30 3 8 Tiempo total ACTCOL 3,82 12,57 30,66 Tiempo sol. sistema hlixto I 3,66 11,34 27,04 ACTCOL . 0,82 5,63 18,17 Mixto I 0,50 4,91 15,74 Tabla 111. Indentación con bloque 3D: comparación tiempos cálculo Malla Algoritmo Iterativo. Actualización del sistema de ecuaciones En los métodos iterativos tipo Newton se tratan las no linealidades de contacto conjuntamente con las derivadas de la geometría y del comportamiento del material. En ese contexto, denominaremos método mixto I al proceso iterativo de Newton-Raphson en el que los sucesivos sistemas lineales planteados se resuelven mediante el método mixto descrito en el apartado anterior. Denominaremos método mixto 11 al proceso iterativo que distingue en la matriz de rigidez la parte del material de la del contacto, actualizando tan sólo aquella parte que lo precise. Para problemas de contacto entre cuerpos elásticos lineales, el método de NewtonRaphson aplicado a la matriz global supone la refactorización completa de la matriz en cada iteración, cuando sólo la parte de contacto ha sufrido modificaciones. El esquema mixto permite calcular en cada iteración sólo las matrices de contacto, resolver el problema condensado y mediante sustituciones progresivas y regresivas sobre la primera y única factorización de las matrices de rigidez Ki, obtenerlos nuevos desplazamiento^'^. De igual manera si se tiene un cuerpo elástico lineal en contacto con uno plástico es posible en cada iteración, recalcular y factorizar sólo la parte plástica y las matrices de contacto. Esta idea es especialmente interesante para el problema de la laminación donde el comportamiento del rodillo es elástico y el comportamiento del desbaste es elastoplástico. Un esquema similar a éste fue planteado en15, aunque sin distinguir los cuerpos en contacto, es decir, considerando como un único bloque las partes de la matriz que representan la deformación mecánica. Se introducen dos niveles de iteraciones: las iteraciones menores o de contacto que se realizan dentro del bucle de iteraciones mecánicas (iteraciones mayores). Para el problema de laminación esta idea no parece tan interesante, puesto que las no linealidades de contacto y las de la deformación plástica tienen lugar en la misma zona del desbaste, la que justamente se encuentra debajo del rodillo. El método propuesto en este trabajo no separa las iteraciones de contacto de las mecánicas, sino que ambas se realizan en un único bucle. Cada uno de los cuerpos en contacto constituye un grupo de elementos con su matriz de rigidez y su residuo. Una variable de actualización controla cuando dicha matriz va a ser recalculada y factorizada No X . Tiempo total ACTCOL Mixto I Reduc. Tiempo sol. sistema ACTCOL Mixto 11 dentro del proceso iterativo. Calculados los residuos, se asignará un estado a la variable de actualización de cada cuerpo según el criterio establecido. El criterio considerado en este trabajo consiste en analizar los residuos de cada cuerpo y en base a ellos decidir las matrices que se van a recalcular. Cuando la norma del residuo de un cuerpo es menor que la norma del residuo total por una determinada tolerancia, se considera que para dicho cuerpo la matriz de rigidez calculada en la iteración anterior puede emplearse en la iteración siguiente (actualización de un método Newton modificado). El valor asignado a la tolerancia es importante para la eficacia computacional del método, pero no resulta crítico para la convergencia. Con valores muy bajos de la tolerancia se recalculan y refactorizan todas las matrices, mientras que con valores demasiado altos aumenta el número de iteraciones necesarias en cada incremento para alcanzar la convergencia, no obteniéndose en ambos casos las ventajas deseadas. La experiencia indica que valores de tolerancia entre 0,01 y 0,001 resultan aceptables. Esta estrategia permite combinar de la forma más ventajosa posible un método de Newton con un Newton modificado para un problema general de contacto con plasticidad y deformación finita. El método mixto propuesto permite en potencia la utilización de esquemas de actualización de las factorizaciones de la matriz tangente tipo Broyden o BFGS4-15, que pueden reducir el coste total del análisis. Estos métodos no son aplicables directamente sobre la formulación convencional del problema de contacto mediante el método de los multiplicadores de Lagrange, debido a que la actualización de las restricciones altera la estructura del sistema que debe ser calculado de nuevo. Ejemplo 2 Consiste en la indentación de deformación finita de un cilindro elástico (E = le 7, y = 0,35) de radio 0,4 en una base elastoplástica (E = le 6, y = 0,35, ay = 3333, H = 1010) de dimensiones 0,8 x 0,3, cuya discretización se observa en la Figura 2. Para un desplazamiento prescrito de 0,05, el análisis se realiza en 10 incrementos aplicándose en primer lugar el esquema de Newton puro bajo las dos variantes descritas anteriormente para la resolución del sistema lineal; mediante el método directo (ACTCOL) y mediante el método mixto (mixto I); en segundo lugar se aplica el sistema de Newton parcialmente modificado conjuntamente con el método mixto (mixto II). Los tiempos de cáclulo empleados para cada uno de ellos se resumen en la Tabla IV. Figura 2 muestran los resultados obtenidos para una discretización de 29279 grados de libertad, que tiene el mismo número de elementos (240 x 60) que el caso anterior. Tiempo CPU-red. Iter. CG No ec. adic. No mov. sólido 1 dominio 105,8 2 dominios 64,6 - 39 % 18 242 O 3 dominios 61,l - 42 % 36 484 O 4 dom. (2 x 2) 81,2 - 23 % 34 486 5 6 dom. (3 x 2) 72,4 - 31 % 47 730 8 9 dom. (3 x 3) 67,3 - 36 % 45 976 17 16 dom. (4 x 4) 52,9 - 50 % 47 1470 35 25 dom. (5 x 5) 47,O - 55 % 46 1968 59 36 dom. (6 x 6) 43,l - 59 % 46 2470 89 64 dom. (8 x 8) 43,7 - 58 % 46 3486 167 100 dom. (10 x 10) 45,5 - 57 % 45 4518 269 Tabla VIII. Viga en voladizo 120 x 120: comparación de tiempos de cálculo Tiempo CPU Iter. CG No ec. adic. No mov. sólido 1 dominio 46,42 2 dom. en 2 43,53 23 482 O 3 dom. en 2 57,02 7 1 964 O 2X x 22 59,24 52 606 5 2X x 32 70,69 108 1090 8 3X x 22 59,32 51 730 11 3X x 32 60,09 82 1216 17 4X x 42 57,67 95 1830 35 5X x 52 52,04 102 2448 59 6X x 62 52,35 105 3070 8 9 4X x 22 57,17 45 845 17 8X x 22 52,92 48 1350 4 1 12X x 32 45,40 50 2350 98 Tabla IX. Viga en voladizo 240 x 60: comparación de tiempos de cálculo Se observa cómo los dominios estrechos producen sistemas condensados en la interfase peor condicionados; con la misma dimensión y con menos movimientos de sólido rígido se necesitan más iteraciones de CG para resolver las particiones en 2 que para las particiones en X. Las particiones en dominios cuadrados (8X x 22 y 12X x 32) introducen tantos grados de libertad adicionales que aunque el sistema está bien condicionado no se obtienen mejoras apreciables. SIMULACIÓN NUMÉRICA DEL CONFORMADO METÁLICO Para el mismo ejemplo en tres dimensiones se consideran discretizaciones de 16 x 16 x 16 y 32 x 16 x 8 elementos (13855 y 14671 grados de libertad, respectivamente), obteniéndose los resultados que se recogen en la Tabla X y la Tabla XI, respectivamente. Tiempo CPU-red. Iter. CG No ec. adic. No mov. sól. ríg. 1 dominio 835,7 - - 2 dom. (2 x 1 x 1) -- 25 7 O 4 dom. (2 x 2 x 1) 162,9 - 80 % 53 1836 O 8 dom. (2 x 2 x 2) 157,l - 81 % 57 2916 20 32 dom. (4 x4 ~4) 111,O - 87 % 76 10800 280 Tabla X. Viga en voladizo 16 x 16 x 16: comparación de tiempos de cálculo Tiempo CPU-red. Iter. CG No ec. adic. No mov. sól. ríg. 1 dominio 232,9 2 dom. (2 x 1 x 1) 239,9 18 459 4 2 dom. (1 x 2 x 1) 150,6 - 35 % 52 891 O 2 dom. (1 x 1 x 2) 130,O - 44 % 28 1683 O 4 dom. (2 x 2 x 1) 177,9 - 24 % 72 1404 8 4 dom. (2 x 1 x 2) 188,5 - 19 % 72 2244 1 O 4 dom. (1 x 2 x 2) 120,l - 48 % 79 2772 O 8 dom. (2 x 2 x 2) 154,6 - 34 % 125 1404 20 16 dom. (2 x 4 x 2) 170,3 - 27 % 252 5700 40 16 dom. (4 x 2 x 2) 117,O - 50 % 69 4644 68 32 dom. (4 x 4 x 2) 114,5 - 51 % 112 7200 136 64 dom. (4 x 4 x 4) 200,4 - 14 % 293 12528 280 64 dom. (8 x 4 x 2) 109,7 - 53 % 65 10200 328 Tabla XI. Viga en voladizo 32 x 16 x 8: comparación de tiempos de cálculo Como era de esperar, las reducciones de tiempo de cálculo son mucho mayores en el caso tridimensional aún siendo el número de grados de libertad mucho menor, ya que las reducciones en el ancho de banda son sustancialmente mayores. En la geometría cuadrada se consiguen mejoras de hasta el 60 % en 2 dimensiones y 87 % en 3 dimensiones. Para la geometría rectangular en 2 dimensiones no se consiguen mejoras apreciables, mientras que en 3 dimensiones se reduce el tiempo de cálculo hasta un 53 9%. Ejemplo 5 l Se trata de indentación de un bloque elástico (E = 20000, y = 0,35) en una base elastoplástica (E = 10000, y = 0,35, oy = 500, H = 1111) mediante un gran desplazamiento prescrito de valor igual a la altura del bloque. El análisis se realiza en 10 incrementos de carga con hipótesis de grandes deformaciones y con la discretización que se muestra en la Figura 7. Figura 7. Indentación con bloque. Malla inicial Se muestran en la Tabla XII los tiempos de cálculo obtenidos para la resolución de problema con distintas particiones (X x Y x 2) en dominios de la base. La dirección X corresponde al largo (24 elementos), la Y al ancho (12 elementos) y la Z a la altura (8 elementos). El bloque se considera siempre un dominio excepto en el primer caso, directo, en el que no se aplica el método. En la primera parte de este trabajo se han presentado la formulación, los resultados y la metodología de validación utilizada' para la simulación numérica de un caso de laminación plana de un desbaste de 300 mm de longitud (X) y 350 mm de anchura (Y) cuyo espesor (2) de 100 mm se reduce en un 13,2 % con rodillos de 300 mm de radio. Se completan a continuación aquelllos resultados con los relativos a la eficiencia computacional, obtenidos con los distintos métodos desarrollados para la resolución del proceso iterativo tipo Newton de cada incremento. Rendimiento del método mixto de resolución del contacto En la Tabla XIII se muestran los tiempos de CPU correspondientes a la laminación 2D para distintas discretizaciones del desbaste manteniendo constante la discretización del rodillo. Se observa como el método mixto I resulta ventajoso sólo para sistemas de ecuaciones de tamaño elevado y sólo mediante la estrategia iterativa (mixto 11) se obtienen resultados eficientes también con discretizaciones pequeñas, que tienden a desaparecer con el aumento de los grados de libertad del desbaste, por los motivos aducidos en los ejemplos de contacto anteriores. Tabla XIII. Laminación 2D: rendimiento de los métodos mixtos Tabla XIV. Laminación 3D: rendimiento de los métodos mixtos En la Tabla XIV se muestran los tiempos de CPU correspondientes a la laminación 3D para distintas discretizaciones, donde al aumentar el número de grados de libertad del desbaste, también mejora la eficiencia computacional. De nuevo, para problemas 316 M.A. GUTIÉRREZ, M. OJANGUREN Y J.J. ANZA pequeños la ventaja proviene principalmente de la estrategia iterativa, mientras que para problemas grandes proviene de la resolución del sistema. Rendimiento del método de descomposición de dominios En el caso 2D, el dominio a descomponer se reduce a un rectángulo. Sólo particiones a lo alto del desbaste reducen el ancho de banda de las submatrices, pero a costa de un elevado número de multiplicadores y produciendo dominios con mala relación de aspecto. Particiones a lo largo (Figura 8) dan lugar a sistemas condensados bien condicionados, pero no reducen el ancho de banda de las matrices a factorizar, por lo que no son eficientes. En la Tabla XV se muestran los resultados obtenidos para distintas discretizaciones del desbaste y distintos procedimientos de resolución. Figura 8. Laminación 2D: partición en dominios El coste de las ecuaciones constitutivas elastoplásticas se refleja en la comparación del tiempo de cálculo de la rigidez tangente con el tiempo de resolución del sistema. Estos tiempos sólo se igualan en el último caso con el aumento del número de grados de libertad del problema. El coste de la resolución del sistema crece con el número de dominios, debido al incremento de los grados de libertad adicionales, sin embargo disminuye el coste del cálculo de la rigidez tangente, al no recalcularse las matrices de los dominios que están alejados del rodillo. Puede apreciarse que en 2D ambas tendencias se compensan y el balance final no es especialmente beneficioso. En simulaciones 3D, particiones a lo ancho del desbaste permiten reducir el ancho de banda de las submatrices sin producir dominios mal condicionados y sin un número de nodos interfase muy elevado. De esta forma es posible combinar mediante particiones a lo ancho y a lo largo las ventajas propias del método de descomposición de dominios puros (DD-1) con las introducidas por la estrategia iterativa (DD-11). En la Tabla XVI se muestran los resultados obtenidos para una malla de 18 x 10 x 6 elementos en el desbaste y 100 elementos en el rodillo. En el primer caso, desbaste y rodillo se resuelven conjuntamente mediante el método directo. En los demás casos, las diversas particiones consideran al rodillo como uno de los dominios. En la última columna se indica porcentualmente la mejora conseguida frente al método directo. De manera análoga, en la Tabla XVII se muestran los resultados obtenidos con la estrategia iterativa (DD-11) sobre una discretización más fina: 24 x 12 x 8 elementos del desbaste y 100 elementos en el rodillo. Pueden apreciarse mejoras importantes de hasta el 50 % frente al método directo standard. Tabla XV. Laminación 2D: rendimiento de la descomposición de dominios Tabla XVI. Laminación 3D. D. Dominios (1) 7 dominios 3 x 2 x 1 (11) 97615 25,O % 2 x 3 x 1 (11) 63820 51,O % 9 dominios 2 x 4 x 1 (11) 61663 52,6 % 10 dominios 3 x 3 x 1 (11) 68531 47,4 % 13 dominios 3 x 4 x 1 (11) 66408 49,O % Tabla XVII. Laminación 3D. D. Dominios (11) CONCLUSIONES En este trabajo se ha abordado la resolución eficiente en ordenadores secuenciales del sistema de ecuaciones de un problema de contacto con multiplicadores de Lagrange y la combinación eficiente con un método de descomposición de dominios semejante, también basado en multiplicadores de Lagrange (FETI). Se han obtenido resultados satisfactorios que justifican la metodología desarrollada para resolver roblem mas sólidos no lineales. El tratamiento diferenciado de las distintas zonas del dominio puede reducir sensiblemente los anchos de banda de las matrices locales y conducir a un sistema condensado en la interfase que se resuelve eficazmente mediante el método del gradiente conjugado. En la simulación de procesos de laminación se han obtenido mejoras computacionales de un 50 % en problemas en 3 dimensiones y gran número de grados de libertad. En laminaciones de productos con forma donde la modelización de la geometría del rodillo precisa discretizaciones más finas, caben esperar mejoras más importantei. Los métodos desarrollados permiten potencialmente la utilización de métodos BFGS y precondicionadores del gradiente conjugado, que pueden suponer sendas ventajas a explorar. Además, la descomposición en dominios se amolda especialmente a las nuevas posibilidades que ofrece la computación paralela, lo que puede suponer mejoras adicionales significativas. AGRADECIMIENTOS El segundo autor agradece al Departamento de Educación Universidades e Investigación del Gobierno Vasco por la beca del Programa de Formación de Investigadores concedida para la realización de este trabajo. REFERENCIAS 1. M.A. Gutiérrez, J.J. Anza, "Simulación numérica del conformado metálico. Eficiencia numérica en procesos de laminación. Parte 1", Reu. Int. Mét. Num. Cálc. y Dis. en Ing., Vol. 13, 2, pp. 211-239, (1997). 2. J.C. Simo, "Topics on the Numerical Analysis and Simulation of Plasticity", draft (1994) for Handbook of Numerical Analysis, Vol. 111, to be edited by P.G. Ciarlet y J.L. Lions, Elsevier Science Publishers B.V., North Holland, 3. J.P. Ponthot, "Traitement unifie des milieux continus solides en grandes transformations par la methode des elements finis", Tesis doctoral, Universidad de Lieja, (1994-1995). 4. H. Matties y G. Strang, "The Solution of Nonlinear Finite Element Equations". Int. J. Num. Meth. Engng., Vol. 14, pp. 1613-1626, (1979). 5. J.C. Simo y R.L. Taylor, "Consistent Tangent Operator for Rate-Independent Elastoplasticity", Comp. Meth. Appl. Mech. Engng., Vol. 48, pp. 101-118, (1985). 6. A.F. Arif y N. Zabaras, "On the Performance of Two Tangent Operators for Finite Element Analysis of Large Deformation Inelastic Problems", Int. J. Num. Meth. Engng., Vol. 35, pp. 369-389, (1992). 7. H. Parisch, "A Consistent Tangent Stiffness Operator for Three-dimensional Nonlinear Contact Analysis", Int. J. Num. Meth. Engng., Vol. 28, pp. 1803-1812, (1989). 8. J.J. Dongarra, I.S. Duff y H. Van der Vorst, "Solving Linear Systems on Vector and Shared Memory Computers" , SIAM, Philadelphia, (1991). 9. R. Barret, M. Berry, T. Chan, J. Demmel, J. Donato, J.J. Dongarra, V. Eijkour, R. Pozo, C. Romine y H. Van der Vorst, "Templates for the Solutions of Linear Systems: Building Blocks for Iterative Methods" , SIAM, (1993). 10. G.H. Golub y C.F. Van Loan, "Matriz Computations", 2nd edition, The Johns Hopkins University Press, Baltimore, (1989). 11. R. Glowinski, G.H. Golub, G.A. Meurant y J. Periaux, (Eds), Proceeding of the First International Symposium on Domain Decomposition Method for Partial Differential Equations, SIAM, (1988). 12. T.F. Chan, R. Glowinski, J. Periaux y O. Widlung (Eds), Proceeding of the 2nd International Symposium on Domain Decomposition Method for Partial Differential Equations, SIAM, (1989). 13. T.F. Chan, R. Glowinski, J. Periaux y O. Widlung (Eds), Proceeding of the 2nd International Symposium on Domain Decomposition Method for Partial Differential Equations, SIAM, (1990). 14. F.J. Gallego y J.J. Anza, "A Mixed Finite Element Model for the Elastic Contact Problem", Int. J. Num. Meth. Engng., Vol. 28 , pp. 1249-1264, (1989). 15. Nour-Omid y P. Wriggers, "A Two Levels Iteration Method for the Solution of Contact Problems", Comp. Meth. Appl. Mech. Engng., 54, pp. 131-144, (1986). 16. G.H. Golub y D. O'Leary, "Some History of the Conjugate Gradient and Lanczos Methods", SIAMS, Vol. 31, (1989). 17. O.C. Zienkiewicz y R.L. Taylor, "The Finite Element Method", 4th edition, Vol. 1, Mc.Graw Hill, (1989). 18. R. Flecher, "Practica1 Methods of Optimization", John Willey and Sons, New York, (1981). 19. P.F. Gil1 y W. Murray, "Numerical Methods for Constrained Optimization", Academic Press, London, (1974). 20. P.E. Bjordstad y O. Wildlung, "Iterative Methods for Solving Elliptic Problems on Regions Particionated into Substructures", SIAM J. Num. Anal., Vol. 32, No 6, (1986). 21. F.X. Roux, "Domain Decomposition Methods for Static Problems", Rech. Aérosp., No, (1990-1). 22. Y.H. De Roeck, "Résolution sur ordinateurs multi-processeurs de problkmes d'klasticité par décomposition de domaines", Tesis Doctoral, Université Paris IX Dauphine VER Mathématiques de la Décision, (1991). 23. Y.H. De Roeck, "Nonlinear Elasticity Solved by a Domain Decomposition Method in a Hypercube", "Appl. Num. Math., Vol. 12, (1993). 24. P. Le Tallec, "Domain Decomposition Methods in Computational Mechanics", Comp. Mech. Advances, Vol. 1, K0 2, (1994). 25. C. Farhat y F.X. Roux, "Implicit Parallel Processing in Structural Mechanics", Comp. Mech. Advances, Vol. 2, pp. 1-124., (1994). 26. P.L. Lions, "On the Schwarz Alternating Method I", Proceedings of the First International Symposium on Domain Decomposition Methods for Partial Diflerential Equations, R. Glowsinski et al., SIAM, pp. 1-42, Philadelphia, (1994). 27. C. Farhat y F.X. Roux, "A Method of Finite Element Tearing and Interconnecting and its Parallel Solution Algorithm", Int. J. Num. Meth. Engng., Vol. 32, pp. 1205-1227, (1991). 28. C. Farhat, J. Mande1 y F.X. Roux, "Optimal Convergence Properties of the FETI Domain Decomposition Method", Comp. Meth. Appl. Mech. Engng., Vol. 115, pp. 365-385 (1994). 29. TOP/DOMDEC, Totally Object-Oriented Program for Automatic DOMain DECompositionn .