El método de Glimm

EL METODO DE GLIMM GUILLERMO MARSHALL* + Y ANGEL MENENDEZ+ Centro de Cálculo Cientfico Comisión Nacional de Energía Atómica 1429 Buenos Aires Laboratorio de Hidráulica Aplicada - INCYTH Casilla de Correo 21 1802 Aeropuerto Ezeiza RESUMEN En contraste con otros métodos numéricos tales como diferencias finitas o elementos finitos, el método de Glimm resuelve las ondas de choque y otras discontinuidades con relativa facilidad y sin necesidad de adicionar términos de viscosidad artificial. Desafortunadamente la literatura existente sobre dicho método es, en general, difícil de entender por aquellos que encaran su estudio por primera vez. Trataremos aquí de presentar una introducción sencilla del método de Glimm y sus principales aplicacio~ies. SUMMARY Contrarily to other numerical methods such as finite differences or finite elements, the Glimm's Method resolves impact waves or other discontinuities with relative facility and without being necessary to add artificial viscosity terms. Unfortunately, the existing litterature on this method is, en general, difficult to understand for those who study it for the first time. Here, we shall try to present a simple introduction on the Glimm's Method and its principal applications. INTRODUCCION Y CONSIDERACIONES GENERALES Glimm2 introdujo un método aproximado para la construcción de la solución' de sistemas hiperbólicos de leyes de conservación. Esta construcción es la base para su elegante teorema de existencia. A partir de entonces hubo numerosas tentativas, sin éxito aparente, de utilizar dicha construcción como una herramienta computacional, ya que se intuía su utilidad para cierto tipo de problemas. Una década después Chorinl presenta los primeros resultados computacionales exitosos basados en la construcción de Glimm para problemas de dinámica de gases. Desde entonces se realizaron gran cantidad de trabajos y lo que se dió en llamar el método de Glimm (o el método de la Recibido: Octubre 1985 * Consejo Nacional de Investigaciones. + Universidad de Buenos Aires. Universitat Politecnica de Catalunya (España) ISSN 0213-1315 Revista internacional de métodos numéricos para cálculo y diseño en ingeniería, Vol. 2,3,231-252 (1 986) 232 G. MARSHALL Y A. MENENDEZ DISWNTINUIüAD ONDA DE CHOQUE t Figura 1.2. Problema de la explosión en atmósfera libre. elección aleatoria o el método del muestreo uniforme) resultó ser una técnica de gran precisión en el tratamiento numérico de problemas con discontinuidades. El método de Glimm consiste, en esencia, en realizar un muestreo estadístico de soluciones teóricas locales de problemas de Riemann. El poder de esta técnica radica en que las soluciones teóricas locales contienen información extensiva sobre el fenómeno de interacción elemental de ondas, y el muestreo estadístico evita la aparición de la viscosidad numérica típica de los métodos de diferencias finitas y elementos finitos. En particular, el método de Glimm permite el tratamiento automático de la formación espontánea y la evolución de discontinuidades (típicas de sistemas hiperbólicos no lineales). Figura 1 .l. Flujo de un gas en una tobera convergente-divergente. A continuación, se discuten brevemente algunos problemas físicos para los cuales el niétodo de Glimm es especialmente útil. La Figura 1.1 muestra el flujo estacionario de gas en una tobera. Debido a las variaciones del área transversal existen cambios en el régimen de movimiento. En la zona divergente de la tobera se produce (por efecto de las condiciones a la salida) una transición entre los regímenes supersónico y subsónico a través de una onda de choque. El problema de una explosión se ilustra en la Figura 1.2. Inicialmente se tiene una zona de alta presión en el interior de una esfera de radio pequeño. En el proceso de igualación de presiones, se produce una onda de choque principal que se aleja del origen, una discontinuidad de contacto que avanza en la misma dirección y una onda de choque secundaria que se mueve en sentido contrario. Esta última es un "efecto de la curvatura". Otro problema de interés en dinámica de gases es el de la difracción de ondas de choque. (Ver Figura 5.1 ). e I f-- S' L Figura 1.3. Onda de frente abrupto generada por el cierre parcial de una compuerta. EL METODO DE GLIMM 233 l También en hidráulica existen problemas aptos para ser tratados por esta técnica. Tal es el caso del resalto hidráulico o de las ondas de frente abrupto generadas por la operación de compuertas (ver. Figura 1.3) o la rotura de una presa de embalse. Finalmente el método de Glimm se utiliza exitosamente en la simulación del desplazamiento de un fluído bifásico en un medio poroso. Este problema se presenta, por ejemplo, en la recuperación secundaria y terciaria en reservorios d,e petróleo. EL METODO DE GLIMM. INTRODUCCION La forma general de las ecuaciones a resolver por el método de Glimm es la siguiente: W, + F(W), + G(W,x) = O (2.1) donde F(W) es una "densidad de flujo" y G(W) un "término fuente". Se introducirá el método aplicándolo a la ecuación de onda escalar (homogénea) caracterizada por: W=u ; F(W)=au ; G(W,x)=O (2.2) donde a = constante. A partir de dadas condiciones iniciales, el método de Glimm construye la solución por medio de una sucesion de dos pasos fraccionarios de igual extensión. En el instante de tiempo t, = nk (k es el intervalo temporal de discretización), la solución u(x,tn) se aproxima por una series de estados constantes que se extienden sobre intervalos de longitud h, tal cual se ilustra en la Figura 2.1. Utilizando como valor inicial esa función continua a trozos, se construye la solución teórica de la ecuación (2.1), entre tn y tn+i,,. Esta solución teárica surge de superponer las soluciones locales a los sucesivos "problemas de Riemann" (uno para cada punto xi+t/, = (i+ ) h). El proble- Figura 2.1. AproximaciOn de u(x,tn)  Figura 2.2. SoluciOn del problema por una funci6n constante a trozos.  de Riemann para la ecuaciOn (2.1)- (2.2). u'•1  u l X. G O U, /1,20  ih  (1 • 1/20 ii.11h  x• ih 234  G. MARSHALL Y A. MENENDEZ ma de Riemann para la ecuaci6n (2.1) queda definido por las siguientes condiciones iniciales: u1 = const. para x < 0 u(x,0) =   (2.3) u r = const. para x . > 0 Si u r * 14 1 , las condiciones (2.3) corresponden a una discontinuidad inicial. La soluciOn de este problema de Riemann es trivial: la discontinuidad se traslada, sin deformarse, a lo largo de la curva caracteristica con pendiente dx/dt = a que pass por x = 0. Esto determina las soluciones locales, si se asocia u 1 y u r - con u i n y 1,, respectivamente, para cada valor 'dej. La construcciOn de la soluciOn numerica para 4 0 _ 1 / 2 , es decir, la determinaciOn de los valores de u a ser asignados a cada nodo x i + 1 /2 , , se realiza por medio de un muestreo estadistico de la soluciOn teOrica  + 1 /2 ), dentro de cada i.ntervalo espacial [ih,  Esto significa tomar n + 1/ 2 u i  u( X i , tn+ 1 /2) =  (2.4) +i./2 donde X i = (i+ 1 1 2 + O i /2) h, siendo O i un ninnero generado al azar en el intervalo [-1,1]. Entonces, u n + l /2 "  /2  si  0 1 h < ak ,n 1 + - 1 /2 1 i n_t_1/2 =  1+1 Un  si  O i h> ak (2.5) tal cual se observa en la Figura 2.2 Mediante un procedimiento id6ntico se avanza la soluciOn desde t n +1 /2 hasta tn, centrändose nuevamente la soluciOn numerica en nodos con Indice entero i. LOgicamente, la superposiciOn de las soluciones de los problemas de Riemann correspondientes a cada nodo, puede llevarse a cabo sienrpre y cuando las discontinuidades no interactnen entre si. Esto implica que debe verificarse la condiciOn de CourantFriedrichs-Lewy k > 1 a  (2.6) Una idea de la bondad del método la brinda la siguiente comparación. De acuerdo a la solución teórica del problema, cualquier discontinuidad'inicial se desplazará, luego de un tiempo T, una distancia X = aT respecto de su posición inicial, sin deformarse. La correspondiente distancia Xn calculada por el método de Glimm puede expresarse como donde n = T/k y qj son variables aleatorias independientes, pero idénticamente distribuidas, con una densidad de probabilidad que, de acuerdo a la Figura 2.2, está dada por h-ak Prob [qj = -h/2] = - 2h h+ak Prob [qj = h/2] = - 2h De las ecuaciones (2.7) y (2.8) surge que el valor esperado y la varianza de X, son, respectivamente, k Var [X,] =-(q2 -a2)~ a donde q = hlk. Se observa que el valor esperado coincide exactamente con la solución teórica. Es más, la ecuación (2.5) muestra que la propagación se realiza sin deformación. Es decir que este método no produce atenuación ni dispersión numéricas, a diferencia de los métodos numéricos dependientes de la malla de discretización, tal como diferencias finitas. Finalmente se ve que, para q y T fijos, Var [X] tiende a cero para k tendiendo a cero, lo cual prueba la convergencia del método. Chorinl desarrolló las estratégias para la elección de los valores de Bi, que resultaron cruciales para el éxito del método como herramienta de predicción. En primer lugar, Chorin introdujo la idea de elegir di = 8 = const. para cada paso de tiempo, lo cual evita la aparición de estados constantes espúrios. En segundo lugar, propuso dividir el intervalo [-1,1] en m2 subintervalos (m2 < n) y elegir el valor de 8 para cada paso de tiempo en cada uno de los subintervalos ni (i = 1,2,. ., m2) de acuerdo a la fórmula ni+ 1 = (ml +,ni) mod m2, donde ml y m2 (ml < m2) son enteros primos y no (no < m2) es arbitrario.'Este procedimiento reduce considerablemente la varianza de la solución (respecto, por ejemplo, del valor dado por la ecuación (2.9)). Una leve modificación de esta estratégia es especialmente efectiva para tener en cuenta las condiciones de borde: los subintervalos se separan en dos grupos, de acuerdo a si pertenecen al semi-intervalo [-1,0] o al [0,1], y los sucesivos valores de 8 se eligen alternadamente en cada grupo. La Figura 2.3 muestra resultados numéricos para la evolución de una discontinuidad, obtenidos con distintos valores de ml y m2. Se observa el carácter estadístico de la EL METODO DE GLIMM 235 Figura 2.3. Evolución de una discontinuidad, calculada por el mdtodo de Glimm. t;olución numérica. En efecto, la posición calculada de la discontinuidad fluctúa alrededor del valor exacto. Los valores óptimos de m, y m2 parecerían ser, de entre los probados, ml = 5 y m2 = 1 1. Estos son los utilizados por los autores en los cálculos posteriores. La extensión natural del método de Glimm para resolver ecuaciones inhomogéneas consiste en hallar la solución teórica del correspondiente problema de Riemann. Como esto puede ser, en general, complicado o, incluso, imposible Marshall y Menéndez7 introdujeron la idea de utilizar la idea de utilizar la solución teórica aproximada a primer orden en k. Este procedimiento se ilustrará con la ecuación escalar caracterizada por donde b = constante > O. La ecuación (2.1 ) - (2.10) puede ser reescrita en forma característica como La solución al problema de Riemann para la ecuacion (2.1 1) se muestra cualitativamente en la Figura 2.4. Como resultado de la no lingalidad de la ecuación, la discontinuidad inicial se transforma en una onda simple centrada. Si u, < u,, la onda simple es una onda de rarefacción; si u, > u,, es una onda de choque que se propaga con la velocidad La integración de la ecuación (2.1 1) puede aproximarse por 1 236 G. MARSHALL Y A. MENENDEZ EL METODO DE GLIMM Figura 2.4. Solución del problema de Riemann para la ecuación (2.1) - (2.10) a%!&-l teórica Figura 2.5. Evolución de una onda de choque, calculada por el método de Glimm: 238 G. MARSHALL Y A. MENENDEZ donde u. = u(x(O),o). La ecuación (2.13) muestra que tomar k = O es equivalente a tomar b = O. Esto significa que u. coincide con la solución de la ecuación homogénea, que es conocida. Más aún, de la ecuación (2.14) surge, entonces, que las curvas caracteristicas coinciden, a primer orden en k, con las correspondientes a la ecuación homogénea, que también son conocidas. Con esta información, la solución u(x,k) puede calcularse para cualquier valor dado de x. Resultados numericos obtenidos con esta metodología para el caso de propagación de una onda de choque se muestran en la Figura 2.5, para distintos valores de b. Se observa un buen acuerdo con la solución teórica. Un procedimiento alternativo para tratar ecuaciones inhomogéneas fue propuesto por Sod:consistente en una técnica de desdoblamiento de dos pasos. En el primer paso, el término inhomogéneo es removido y la ecuación homogénea resultante es resuelta por el método de Glimm. En el segundo paso, se resuelve por diferencias finitas la ecuación diferencial ordinaria resultante de considerar solamente el término inhomogéneo, es decir W, + G(W,x) = O (2.1 5) utilizando como condición inicial la solución obtenida en el paso anterior. Si bien para el presente ejemplo (ecuación (2.1) - (2.10)) ambos métodos son equivalentes, el método de Sod, aunque menos riguroso, resulta, en general, más simple que el propuesto por Marshall y Menéndez. Más aún, con la metodología de Sod es posible tratar sistemas de ecuaciones no estrictamente hiperbólicas. Tal es el caso, por ejemplo, del modelo para el desplazamiento unidimensional de un fluido bifásico inmiscible en un medio poroso, caracterizado por la ecuación (2.1 O), con donde u es la saturación de la fase mojante, x y t las variables espacial y temporal, respectivamente, a una constante y k(u) un coeficiente genérico de difusión que se supone "pequeño" (ver Laggiard y Marchal14). Las posibles soluciones al problema del Riemann correspondiente, tomando G(W,x) = O, se muestran cualitativamente en la Figura 2.6. Debido a que F(W) no es convexa, la discontinuidad inicial evoluciona como una perturbación compuesta por una onda de choque y una de rarefacción. Resultados numéricos para la evolución de una discontinuidad inicial se presentan en la Figura 2.7. Se observa como la presencia del término de difusión G(W,x) con k(u) r 0.01 7 suaviza el frente de onda. Figura 2.6. Solución de problema de Riemann para la ecuación (2.1) - (2.16) con G(W,x) = O, a) ul >u, y b) u1 <u,. EL METODO DE GLIMM 239 Figura 2.7. Evolución de una onda de choque con difusión calculada con el metodo de Glimm. A u(x,~) 1.0,~ 0.9-e 0.8 0.7 0.6 a5 04 FLUJO EN AGUAS POCO PROFUNDAS Dx* 1140 O Te0.X) DT0.01 X T*O7O x 6 9 04 • T* 1.00 1 . A* l - - - *e. O*.... - 0 o 0 o XxxxX O 0 XXIxxx - O x I X O x Las ecuaciones unidimensionales para aguas poco profundas, o ecuaciones de Saint Venant, pueden expresarse de la forma (2.1), con donde u es la velocidad del agua, h la profundidad, g la aceleración de la gravedad y R(W) un término que representa los efectos de la fricción y la pendiente del fondo. Las coordenadas x,t. adquieren ahora el significado habitual de espacio y tiempo. La aplicación del método de Glimm al sistema homogéneo (es decir, considerando G(W,x) 0) fue realizada por Marshall y Menénde9, utilizando la analogia con dinámica de gases. Más tarde, Marshall y Menéndez6 desarrollaron la aplicación directa, que se discute brevemente a continuación. La aplicación del método de Glimm a cualquier sistema hiperbólico requiere resolver el correspondiente problema de Riemann. Para el sistema (2.1) - (3.1), el problema de Riemann está caracterizado por las siguientes condiciones iniciales: W, = const. para x < O W(x,O) = W, = const. para x > O La solución de las ecuaciones (2.1) - (3.1) - (3.2) con G(W,x) = O, se ilustra en la Figura 3.1. La discontinuidad inicial~~tr~forma en "ondas de depresión, y/u "ondas G. MARSHALL Y A. MENENDEZ 1 Figura 4.4. b) mdtodo de Glimni. Figuia 4.5. Lineas isobáricas en el plano x-t. a) mitodo de Sod. EL METODO DE GLIMM 24 7 10 20 Figura 4.5. b) método de Glimm. FLUJOS SUPERSONICOS BIDIMENSIONALES ESTACIONARIOS Los problemas de difracción de ondas de choque pueden ser descriptos en forma aproximada por un sistema hiperbólico de leyes de conservaciones de la forma donde r, z y t son las variables independientes espacio-temporales, y Aquí p es la densidad, m = pu y n = pv, u y v son los componentes de la velocidad en las direcciones r y z, respectivamente, E es la energía total, P es la presión termodinámica y d es la dimensión del espacio. Estas ecuaciones están suplementadas con la ecuación de estado Asumimos que el gas es politrópico, de modo que donde e es la energía específica. El sistema (5.1) a (5.4) va acompañado de condiciones iniciales y de contornos apropiadas. Muchos problemas interesantes pueden ser descriptos por el sistema no evolucionario asociado, es decir El sistema (5.5) puede ser: globalmente hiperbdlico en el caso de flujo supersónico, mixto hiperbólico y elíptico en el caso de flujo supersónico con zona subsónica, y globalmente elíptico en el caso de flujo subsónico. Esta distinción es muy importante pues el comportamiento de dichos flujos difiere entre sí notablemente. Una variada gama de problemas estacionarios son puramente supersónicos y pueden ser descriptos, en consecuencia, por un sistema hiperbólico de leyes de conservación. Existe una analogía entre flujo estacionario supersónico y flujo evolucionario unidimensional. Este último puede ser descripto por un sistema de ecuaciones como el (2.1) - (4.1). Esta analogía se ilustra en la Figura 5.1. El flujo evolucionario unidimensional mostrado en la Figura 5.1 (a) corresponde al problema de un pistón que se acelera y luego se detiene bruscamente. El flujo estacionario supersónico equivalente, mostrado en la Figura 5.1 (b), está producido por una cuña plana bidimensional. En ambos casos aparecen un onda de choque y una de rarefacción. La segunda produce la difracción de la primera. 248 G. MARSHALL Y A. MENENDEZ Figura 5.1. Analogía entre flujo evolucionario unidimensional y flujo supersónico estacionario bidimensional: a) flujo producido por la aceleración y desaceleración bruscas de un pistón, b) flujo supersónico alrededor de una cuña plana. El ejemplo precedente pone en evidencia que cualquier procedimiento de cálculo de "marcha en el tiempo" puede también ser utilizado para flujo estacionario supersónico en dos dimensiones espaciales. Sólo se requiere considerar a una de las variables espaciales fr ó z) como una variable temporal ficticia. En particular, el método de Glimm resulta especialmente apto para este tipo de problemas. Su utilización requiere formular y resolver el correspondiente problema de Riemann. Los detalles de este procedimiento, un tanto engorrosos, se omitirán del presente trabajo (ver Marshall y Plohrs ). Resultados numéricos para el problema de la cuña se presentan en la Figura 5.2. Los obtenidos por el método de Glimm se comparan con los calculados por diferencias finitas sobre las características (de trabajosa implementación), observándose un buen acuerdo. Figura 5.2. Flujo supersónico sobre una cuña plaria, a) método de las características. Figura 5.2. b) método de Giimm. REFERENCIAS 1. A.J. Chorin. J. Comput. Phys. 23,5 17, (1976). 2. J. Glimm. Comm. Pure Appl. Math. 18,697, (1 965). 3. J. Glimm, G. Marshall and B. Plohr, Advances in Appl. Math. 5, 1. (1984). 4. E. Laggiard y G. Marshall. CNEA-NT 9/82, (1982). 5. G. Marshall and R. Menéndez. J. Compt. Phys. 39 1, (1981). 6. G. Marshall and A.N. Menéndez. Adv. Water Resources 4, 125, (1981a). 7. G. Marshall and A.N. Menéndez.J. Comput. Phys. 44,167, (1981b). 8. G. Marshall and B. Plohr. J. Comput. Phys., en prensa, (1984). 9. G. Sod. J. FluidMech. 83, pt. 4,787, (1977). Como ilustración, se presenta un listado en lenguaje FORTRAN para la resolución de la ecuación u, + u u, = O con condiciones iniciales: u = ul para x< O y u = u, para x > O, utilizando el método de Glimm. Las variables de entrada son: N = número de nodos. N1 = número de nodos con u = u,. 11, Iz = valores iniciales para generar la secuencia de números aleatorios. M = número de pasos de tiempo. MED = control de impresión. . EL METODO DE GLIMM 25 1 DIMENSION U(101) READ(6,100)N, NI, I1,12, M, MED, NI, Ml, M2 FORMAT(I5) READ(6,l lO)UO, UF, Q FORMAT(F 1 O. O) N2=2*N N12=N1*2 NIU=N2+ 1 NIA=N21 NI=NI-M1 DO 10 I=2, N2,2 U(I)= u0 IF(1. GT. N1 2)U(I)=UF CONTINUE DO50 J=l,M L= o NI=NI+ M1 IF(N1. GE. M2)NIzNI-M2 HN=NI HM=M2 R= RAN(I1,72) R=(2. * (HN+ R)/HM - 1 .)*Q W RITE(6,lO 1 )R FORMAT(12F6. 2) IF(L. EQ. 1) GO TO 42 DO 20 1=3 NIA, 2 CALL RiIEMAN(U(1-1), J(I+ 1 ), R, U(!)? CONTINUE U(I)-U0 U(?iIU)= UF L= 1 IF(MED. GE. 0)GO TO 40 WRITE(6,lOl) (U(I), I= 1, NIU, 2) GO TO 40 DO 30 I=2, N2,2 CALL RIEMAN(U(1l), U(I+ 1), R, U(1)) CONTINUE IF(MED) 43,43,44 IE= J/IABS(MED) IF(MED*IE-J) 50,43,50 WRITE(6, 101)(U(I), I=2, N2, 2) CONTINUE END SLJBROUTINE RIEMAN (UL, UR, R, UM) IF (UL. GT. UR) GO TO 10 IF (R. GT. UR) GO TO 30 IF (R. LT. UL) GO TO 20 UM= R GO TO 40 UP= 5 *(UR+ UL) IF (R. GT.UP) GO TO 30 UM=UL GO TO 40 UM=UR CONTINUE RETURN END