Algoritmo numérico para problemas de elastoplasticidad con endurecimiento mixto

Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería. Vol. 6 ,3, 437-448( 1990) ALGORITMO NUMERICO PARA PROBLEMAS DE ELASTOPLASTICIDAD CON ENDURECIMIENTO MIXTO JOSE ANGEL HUIDOBRO Y ARTURO ROBLES Departamento de Matemáticas, Escuela T.S.I.I. de Gijón - Universidad de Oviedo, Carretera de Castiello s/n, Gijón 33204, Asturias. RESUMEN En este trabajo se presenta un método numérico para la resolución de un problema de elastoplasticidad con endurecimiento mixto no lineal. El método consiste en un algoritmo basado en técnicas del lagrangiano aumentado combinado con un método de elementos finitos mixto. La tensión y los parámetros de endurecimiento se aproximan por funciones constantes a trozos y los desplazamientos mediante funciones lineales a trozos. SUMMARY In this paper we present a method for the numerical resolution of a elastoplasticity problem with nonlinear mixed hardening. The method is based on Augmented Langragian thechniques. We use a mixed finit element method: stress and hardening parameters are approxiinated by piecewise constant functions and displacements by piecewise linear functions. INTRODUCCION La idea de endurecimiento mixto fue sugerida por Prager en 1935 y consiste en combinar los conceptos de endurecimiento isotrópico y cinemático. Está especialmente indicado en la modelización de procesos de endurecimiento en los que no puede prescindirse del efecto Bauschinger, como sucede por ejemplo en casos de carga cíclica. El modelo que se considera en este trabajo se enmarca en la teoría inñnitesimal y viene definido por: superficie de fluencia de Von Mises, regla de flujo de Prantdl-Reuss, endurecimiento isotrópico por deformación y regla de endurecimiento cinemático de Prager. Este modelo ha sido analizado por Axelsson y Samuelsson', donde se formula el problema mediante la construcción de la matriz elastoplástica. Para su resolución Recibido: Enero 1989 OUniversitat Politecnica de Catalunya (España) ISSN 0213-1315 437 J. HUIDOBRO Y A. ROBLES se tiene dh . dF --i) = Adrl drl y para la ley de Prager 1. dF --a = Ac da obteniéndose entonces que las ecuaciones (4)-(5)-(6) pueden expresarse como X>O, FSO, FX=O lo que indica que el vector (8, -$&, -hl(q)i)) es normal saliente a la superficie de fluencia en el punto (a, a, q). A partir de esta relación y para el caso de endurecimiento isotrópico lineal, Johnsonla formula variacionalrnente el problema y obtiene resultados de existencia y unicidad de solución. Para el caso no lineal el tratamiento del problema resulta más complejo tanto en el aspecto teórico como práctico. En este trabajo se introduce un cambio de variable que simplifica la condición de normalidad anteriormente mencionada. Sea dh + v = g(q) siendo g(7) = 1/' (d;) ds Mediante estas relaciones la ley de endurecimiento isotrópico puede expresarse como . dF -i = Acon F(a, a, v) = ~(a, a, g-l(v)) ay De esta forma la evolución de las variables plásticas queda definida por izo, F~O, FX = O o también, .(T-o,p-a,p-v)LO V(T,~,~)EB B = {(r,P,p) E M, x M, x R : F(7,P,p) < 0) (11) Esta inecuación permite obtener fácilmente una formulación variacional del problema. ALGORITMO NUMERICO PARA PROBLEMAS DE ELASSOPLASTICIDAD 441 FORMULACION VARIACIONAL Se consideran los espacios de funciones V = {WE(H'(~))~ : wlrl = O) , 1 (V Por la desigualdad de Korn, para w E V, Ic(w)l (norma en L2(Q, M,) induce sobre V una norma equivalente7 a la natural de V. Se supondrá que f E ~~(0, T; (~~(fl))~) y g E ~~(0, T; (H-~(I'))~) y se define sobre V el operador lineal continuolg Lt(w) = (f(t), w) + (g(t), w)r, wrv (e, a) producto escalar en (L2(Q))3 (-, -)rl dualidad en (H +(F))~ Se introduce el conjunto de equilibrio E(t) = {(T, ,f3, p) E H : (T, E(w)) = Lt(w) VW E V) y el subconjunto convexo y cerrado de H fJ = ((7, P,p) E H : F(7, (x),P(z), ~(2)) 1 O C-P.~. Q) Finalmente se considera la forma bilineal a(a,~) = Aijkhflij(z)~kh(z)d;i: , para U,T E L~(Q, M,), J, el operador Aa = a(u, .) y la norma asociada 1 - 1,. Obsérvese que esta norma es equivalente a la natural de L2(Q, M,) debido a la hipótesis (2b). Denotando por v la velocidad de desplazamiento, el problema puede formularse variacionalmente como: Hallar (a, a, v, u) : [O, TI -+ H x V verificando c.p.d. en [O, T] J. HUIDOBRO Y A. ROBLES a(O) = O, a(O) = O, u(O) = O (14) La existencia y unicidad de solución ha sido demostrada en [lo] bajo las hipótesis: H1. P rl E(t) # q5 c.p.d. en [O,T] H2. 3 X: R x [O,T] -+ M, tal que (~(t),r(w)) = Lt(w) Vw E V Y 121 < C c.p.d. (O, T) x a, j = 1,2 H3. f, fi L2(0, T; (L2(fl))3) PROBLEMA DISCRETO Se supondrá en lo que sigue que Q es un abierto poliédrico. Se considera una familia regular de elementos finitos (triángulos o tetraedros): ñ= " T TETh y se introducen los espacios de dimensión finita donde P,,(T) denota los polinomios de grado n sobre T. Sean además, Eh(t) = ((7, P, p) E Hh : (7, E@)) = Lt(w) VW E Vs) Ph = PnHh Se formula el problema aproximado (discreto en la variable espacial): Hallar (u~,~~,u~,v~) : [O,T] + Hh x Vh verificando c.p.d. en [O,T] La existencia y unicidad de solución se obtienen de forma análoga al problema en continuo si se dispone de hipótesis del tipo H1, H2, H3 para el problema aproximado. Con este objeto se prueba siguiendo a Johnson13, el lema: ALGORITMO NUMERICO PARA PROBLEMAS DE ELASTOPLASTICIDAD 443 LEMA: Sea lTh la proyección ortogonal de H en Elh. Entonces nh(P n E(t)) c (Ph n (Eh(t)) Demostración Sea (T,/?, p) E P n E(t). Su proyección sobre Hh viene dada por: (IT 1 medida del elemento T) Puesto que E(W) es constante en cada elemento T, (7, ~(w)) = Lt(w) * (~h, ~(w)) = Lt(w) W E Vh Y por tanto (~h, Ph vh ) E Eh ( t ) Por otra parte, puesto que T = g-loh es cóncava, la desigualdad de Jensen conduce a lo que prueba que (T~, Ph, vh) E Ph. La convergencia del método de elementos finitos queda establecida en el siguiente teorema: TEOREMA. Existe una constante C independiente de h tal que \(gh, ah, vh) - (g,a,v)12,~ < C inf {Iv - ~12,~ , w E L2(o, T; vh)) (1 . J2,w norma en L~(o,T; w)) Demostración Introduciendo la solución del problema aproximado en (13) y la proyección de la solución del problema en continuo en (16), y sumando ambas expresiones se obtiene 1. a(+) - uh(t), ~(t) - gh(t)) + -(M - irh(t), - ah(t))+ C +(~(t) - Uh(t)> v(t) - ~h(t)) 5 (E(v(~)) - ~('~h(t)), ~(t) - ~h(t)) de donde J. HUIDOBRO Y A. ROBLES 5 Ic(v(t)) - ~(w) l2 + lu(t) - uh(t) 1 5 y por la desigualdad de Korn <cI~(t)-~l:,v+U(t)-Uh(t)1~ VwEVh (C cte positiva) obteniéndose la estimación del teorema aplicando el lema de Gronwall y teniendo en cuenta la equivalencia entre la norma de L2(R, M,) y ( la. Finalmente, para obtener el problema discreto en la variable temporal, se introduce un método en diferencias finitas: Sea N un número natural y k = T/N; para n = 0,1,. . . , N se denota y se plantea el problema discreto: Hallar (u;, a:, u;) E Hh y V; E V para n = 0,1,. . .,N, verificando Este problema equivale a la minirnización sobre el conjunto convexo y cerrado Ph n E: del funcional convexo y coercivo Para asegurar la existencia del mínimo es suficiente probar que el conjunto convexo no es vacío y esto se obtiene a partir de las hipótesis H1, H2 '913. La convergencia del problema discreto al problema aproximado se realiza mediante la técnica'standard de "estimación a pri~ri"~?"~~. ALGORITMO NUMERICO PARA PROBLEMAS DE ELASTOPLASTICIDAD 445 METODO ITERATIVO Antes de introducir el esquema iterativo se expresará la inecuación (19) mediante el concepto de subdiferencial de la función indicadora del convexo de elasticidad. Con este objeto se introducen algunos resultados del análisis convexo (para más detalles ver [6l). Sea W espacio de Hilbert y K un subconjunto cerrado y convexo de W. Se define la función indicadora de K por: Identificando W con su dual, se dice que z está en la subdiferencial de IK en x si: (z,y - x) 5 O tly E K y se denota t E ¿3IK(z) (21) La aproximante Yosida de la subdiferencial de la función indicadora viene expresada por: donde I denota aplicación identidad y TIK la proyección de W en K. Teniendo en cuenta (21) el sistema (18)-(19)-(20) puede expresarse como: Hallar (u;, a;, u;) E Hh y u; E V para n = 0,1,. . .,N, verificando 1 (€(VE) - AW, --Sa:, c -Su:) E ¿31~(o:, a:, u:) (23) La ecuación multivoca (23) marca el carácter altamente no lineal de los problemas de plasticidad y su tratamiento constituye un punto crucial en la resolución numérica del problema. Diversos métodos basados en técnicas de penalización-dualidad han sido desarrollados. Citemos por ejemplo a Mercierl' y Bermudez Viaño5 para problemas de plasticidad perfecta y a Johnsonl', Samuelsson-Froier18 y Bermúdez-Viaño5 en plasticidad con endurecimiento lineal. El algoritmo que aquí se desarrolla aplica este tipo de técnicas a problemas con endurecimiento mixto no lineal, siendo las velocidades de la deformación plástica y de los parámetros de endurecimiento las que juegan el papel de multiplicador de Lagrange. Designando con (ph, qh, S;) el elemento de la subdiferencial que realiza la igualdad en (23), esta ecuación puede expresarse como J. HUIDOBRO Y A. ROBLES Estas ecuaciones permiten introducir el siguiente esquema numérico (se suprime el subíndice h en la notación): Sea (PO, qo, so) arbitrario. Conocidos (pm, q;, S", se determinan u:, a;, v:, v; resolviendo El nuevo multiplicador se calcula mediante Se observa que el sistema (28)-(29) es básicamente un problema de elasticidad lineal, cuya matriz de rigidez es constante en todas las iteraciones, y donde el multiplicador plástico juega el papel de un segundo miembro variable en cada iteración. En la resolución de la ecuación (30) hay que realizar el cálculo de la proyección; puesto que la tensión y los parámetros de endurecimiento son constantes por elemento, ésta se reduce a una proyección en R6 x R6 x R. El sistema (28)-(29)-(30) queda incluido en el marco funcional abstracto desarrollado en [3], en el que ha sido analizada su convergencia. EJEMPLO NUMERICO El método numérico anteriormente descrito ha sido implementado utilizando la biblioteca de elementos finitos MODULEF-INRIA, siendo necesaria la construcción de un nuevo módulo. Como ejemplo numérico se ha analizado la propagación de la zona plástica en una placa de aluminio entallada (Figura 1)) que tiene un comportamiento elastoplástico con endurecimiento mixto no lineal. Para un test de tensión uniaxial, la relación tensión-deformación en el rango plástico es del tipo Ramberg-Osgoodl: u e=-+- - 3u ( u 1"-' E 7E u0.7 Los valores .de la tensión de fluencia, u0.7 y n han sido analizadas en [2], considerándose en este caso 320 MPa, 366 MPa y 37 respectivamente. Las constantes elásticas son E = 70000 MPa y v = 0.2. En la Figura 2 se representa la malla utilizada que consta de 71 nodos y 105 triángulos. En la Figura 3 se representa la evolución de la zona plástica según los valores de la carga. ALGORITMO NUMERICO PARA PROBLEMAS DE ELASTOPLASTICIDAD 447 Figura 3. Propagación de la frontera plástica. AGRADECIMIENTOS Agradecemos al profesor Carlos Moreno sus valiosas sugerencias durante la elaboración de este trabajo. REFERENCIAS 1. K. Axelsson y A. Sainuelsson, "Finite Element Analysis of elastic-plastic niaterials displaying mixed hardening", Int. J. for Num. Meth. in Eng., Vol. 14, pp. 211-225, (1979). 2. H. Armen, A. Pifko y H.S. Levine, "Finite Element analysis of structures in the plastic range", NASA CR 1649, (1971).