Sobre el método de Godunov para sistemas hiperbólicos no conservativos

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Sobre el m´etodo de Godunov para sistemas hiperb´olicos no conservativos M.L. Mu˜ noz Ruiz1, C. Par´ es Madro˜ nal2 1Dpto. Matem´atica Aplicada, Universidad de M´alaga, Campus de Teatinos s/n, 29071 M´alaga. E-mail: [email protected]. 2Dpto. An´alisis Matem´atico, Universidad de M´alaga, Campus de Teatinos s/n, 29071 M´alaga. E-mail: [email protected]. Palabras clave: sistemas hiperb´olicos no conservativos, m´etodo de Godunov Resumen En este trabajo se aborda la aproximaci´on num´erica del problema de Cauchy para sistemas hiperb´olicos no conservativos en dimensi´on uno. El concepto de soluci´on d´ebil de dichos sistemas se define utilizando la teor´ıa de Dal Maso, Le Floch y Murat, basada en la elecci´on de una familia de caminos en el espacio de estados. En primer lugar, establecemos hip´otesis para la elecci´on de esta familia de caminos y estudiamos sus implicaciones, en particular, la obtenci´on de una expresi´on del m´etodo de Godunov que generaliza su expresi´on cl´asica para sistemas de leyes de conservaci´on. A continuaci´on se estudian las propiedades de buen equilibrado de estos m´etodos. Finalmente, probamos la consistencia del esquema num´erico obtenido con la definici´on de soluciones d´ebiles. En concreto, probamos que, bajo la hip´otesis de variaci´on total acotada, si las aproximaciones obtenidas mediante un m´etodo de Godunov basado en una familia de caminos converge uniformemente a alguna funci´on cuando la malla se refina, entonces esta funci´on es una soluci´on d´ebil del sistema no conservativo, relativa a esa familia de caminos. Este resultado se extiende a los esquemas num´ericos basados en Resolvedores de Riemann Aproximados. 1. Introducci´on En este trabajo se estudia la aproximaci´on num´erica de sistemas hiperb´olicos no conservativos: ∂W ∂t +A(W)∂W ∂x = 0, x ∈R, t > 0,(1) donde W(x, t) toma valores en un subconjunto abierto y convexo Ω de RNyAes una funci´on regular y localmente acotada de Ω en MN(R). Suponemos que este sistema es 1 M.L. Mu˜noz, C. Par´es estrictamente hiperb´olico y que los campos caracter´ısticos asociados son o bien genuinamente no lineales o bien linealmente degenerados. En particular, los sistemas de E.D.P. de la forma ∂W ∂t +∂F ∂x (W) = B(W)∂W ∂x +S(W)dσ dx,(2) donde Wpertenece a un subconjunto abierto y convexo Ω de RN,Fes una funci´on regular de Ω en RN,Bes una funci´on regular de Ω en MN(R), Ses una funci´on de Ω en RNy σes una funci´on conocida de Ren Rse pueden reescribir bajo la forma cuasilineal (1) a˜nadiendo a (2) la ecuaci´on ∂σ ∂t = 0. El sistema (2) incluye como casos particulares: sistemas de leyes de conservaci´on (B= 0, S= 0); sistemas de leyes de conservaci´on con t´ermino fuente o leyes de equilibrio (B= 0); y sistemas acoplados de leyes de conservaci´on como los del tipo definido en [1]. 2. Soluciones d´ebiles Consideramos el problema en formulaci´on no conservativa (1), bajo las hip´otesis establecidas en la secci´on anterior. En el caso en que una soluci´on Wes discontinua, el producto no conservativo A(W)Wx no tiene sentido en el marco de las distribuciones. Sin embargo, haciendo uso de la teor´ıa desarrollada por Dal Maso, LeFloch y Murat [2], se puede dar una definici´on rigurosa de soluciones d´ebiles, definici´on asociada a la elecci´on de una familia de caminos en Ω, es decir, una aplicaci´on Φ: [0,1] ×Ω×Ω→Ω localmente lipschitziana que verifica Φ(0; WL, WR) = WLy Φ(1; WL, WR) = WR,∀WL, WR∈Ω. Bajo ciertas hip´otesis de regularidad sobre Φ, dada una funci´on W∈(L∞(R×R+)∩ BV (R×R+))N, el producto no conservativo A(W)Wxpuede ser interpretado como una medida de Borel, que denotaremos por [A(W)Wx]Φ. Si Wes regular a trozos, esta medida de Borel asociada al producto no conservativo viene dada, en un tiempo fijado t, por A(W(·, t))∂W ∂x (·, t)Φ , ϕ=ZR A(W(x, t))∂W ∂x (x, t)ϕ(x)dx +X mZ1 0 A(Φ(s;W− m, W+ m))∂Φ ∂s (s;W− m, W+ m)dsϕ(xm(t)),∀ϕ∈ C0(R).(3) La expresi´on Wxde la primera integral representa la derivada puntual de W(·, t); xm(t) representan la posici´on de las discontinuidades de Wen el tiempo t;W− myW+ mson los l´ımites de Wa izquierda y derecha, respectivamente, de la discontinuidad m-´esima en el tiempo t; y C0(R) es el conjunto de aplicaciones continuas con soporte compacto. Una funci´on W∈(L∞(R×R+)∩BV (R×R+))NC1a trozos es una soluci´on d´ebil de (1) si satisface la igualdad ∂W ∂t +A(W)∂W ∂x Φ = 0. 2 Sobre el m´etodo de Godunov para sistemas hiperb´olicos no conservativos Se dir´a que Wes una Φ-soluci´on de (1) si es una soluci´on d´ebil asociada a la familia de caminos Φ. Cuando no haya lugar a confusi´on, se omitir´a la dependencia de Φ en la notaci´on. A lo largo de una discontinuidad, una soluci´on d´ebil debe satisfacer la condici´on de Rankine-Hugoniot generalizada: Z1 0ξI − A(Φ(s;W−, W+))∂Φ ∂s (s;W−, W +)ds = 0,(4) donde ξes la velocidad de propagaci´on de la discontinuidad, Ies la matriz identidad, y W−yW+son, respectivamente, los l´ımites a izquierda y derecha de la soluci´on a lo largo de la discontinuidad. En el caso particular de un sistema de leyes de conservaci´on, es decir, cuando A(W) corresponde a la matriz jacobiana de una funci´on de flujo F(W), la definici´on de producto no conservativo como una medida de Borel es independiente de la elecci´on de caminos y la condici´on de Rankine-Hugoniot generalizada se reduce a la condici´on de Rankine-Hugoniot habitual. Al igual que ocurre en el caso conservativo, no toda discontinuidad es admisible. Por lo tanto, se hace necesario asumir un concepto de soluci´on entr´opica, como el de Lax o bien el siguiente, basado en un par de entrop´ıa (η, G) para (1), i.e. un par de funciones regulares de Ω en Rtal que ∇G(W) = ∇η(W)· A(W),∀W∈Ω, y seg´un el cual una soluci´on d´ebil se dice que es entr´opica si satisface la desigualdad ∂η(W) ∂t +∂G(W) ∂x ≤0 en el sentido de las distribuciones. 3. Elecci´on de caminos La elecci´on de la familia de caminos es importante ya que determina la velocidad de propagaci´on de las discontinuidades. La elecci´on m´as simple es la dada por la familia de segmentos, que corresponde la definici´on de productos no conservativos introducida por Volpert en [7]. En general, debe estar basada en los aspectos f´ısicos del problema; sin embargo, desde el punto de vista matem´atico es natural exigir a esta familia que satisfaga algunas hip´otesis concernientes a la relaci´on de los caminos con las curvas integrales de los campos caracter´ısticos. En este trabajo asumiremos que la familia de caminos satisface las siguientes hip´otesis: (H1) Dados dos estados WLyWRpertenecientes a la misma curva integral γde un campo linealmente degenerado, el camino Φ(·;WL, WR) es una parametrizaci´on del arco de γque une WLyWR. (H2) Dados dos estados WLyWRpertenecientes a la misma curva integral γde un campo genuinamente no lineal, Ri, y tal que λi(WL)< λi(WR), el camino Φ(·;WL, WR) es una parametrizaci´on del arco de γque une WLyWR. 3 M.L. Mu˜noz, C. Par´es (H3) Denotemos por RP ⊂ Ω×Ω el conjunto de pares (WL, WR) tales que el problema de Riemann        Wt+A(W)Wx= 0, W(x, 0) = (WLif x < 0, WRif x > 0, (5) tiene una ´unica soluci´on d´ebil autosimilar compuesta por a lo sumo Nondas simples (i.e. choques entr´opicos, discontinuidades de contacto u ondas de rarefacci´on) conectando J+ 1 estados intermedios constantes W0=WL, W1, . . . , WJ−1, WJ=WR, con J≤N. Entonces, dados (WL, WR)∈ RP, la curva descrita por el camino Φ(·;WL, WR) en Ω es igual a la uni´on de las correspondientes a los caminos Φ(·;Wj−1, Wj), j= 1, . . . , J. El establecimiento de estas hip´otesis nos permite probar las propiedades siguientes: Proposici´on Asumamos que el concepto de soluciones d´ebiles de (1) est´a definido en la base de una familia de caminos que satisface las hip´otesis (H1)-(H3). Entonces: (i) Dados dos estados WLyWRpertenecientes a la misma curva integral de un campo linealmente degenerado, la discontinuidad de contacto dada por W(x, t) = (WLif x < σt, WRif x > σt, donde σes el valor constante del autovalor correspondiente a lo largo de la curva integral, es una soluci´on d´ebil entr´opica de (1). (ii) Sea un par de estados (WL, WR)de RP y sea Wla soluci´on del problema de Riemann (5). Se tiene, para cada t > 0, la desigualdad hA(W(·, t))Wx(·, t),1i=Z1 0 A(Φ(s;WL, WR))∂Φ ∂s (s;WL, WR)ds. Por tanto, la masa total de la medida de Borel A(W(·, t))Wx(·, t)es independiente de t. (iii) Sea (WL, WR)un par de RP y sea Wjcualquiera de los estados intermedios en la soluci´on del problema de Riemann (5). Entonces Z1 0 A(Φ(s;WL, WR))∂Φ ∂s (s;WL, WR)ds =Z1 0 A(Φ(s;WL, Wj))∂Φ ∂s (s;WL, Wj)ds +Z1 0 A(Φ(s;Wj, WR))∂Φ ∂s (s;Wj, WR)ds. 4 Sobre el m´etodo de Godunov para sistemas hiperb´olicos no conservativos Es posible dar un procedimiento general para construir una familia de caminos que satisfaga las hip´otesis (H1)-(H3), al menos para la clase RP, extendiendo la teor´ıa de ondas simples de sistemas hiperb´olicos de leyes de conservaci´on y los resultados concernientes a las soluciones de los problemas de Riemann al caso de sistemas hiperb´olicos no conservativos. 4. El m´etodo de Godunov para sistemas hiperb´olicos no conservativos Consideremos el sistema (1) bajo las condiciones establecidas en la secci´on 1, con condici´on inicial W(x, 0) = W0(x), x ∈R,(6) y asumamos que la familia de caminos Φ usada en la definici´on del producto no conservativo A(W)Wxsatisface las hip´otesis (H1)-(H3). Para discretizar el sistema consideramos celdas Ii= [xi−1/2, xi+1/2] que supondremos, por comodidad, de tama˜no constante ∆xy definimos xi+1 2=i∆x. Definimos tambi´en xi= (i−1/2)∆x, el centro de cada celda Ii. Sea ∆tel paso de tiempo y definamos tn=n∆t. Denotamos por Wn ia la aproximaci´on de las medias en cada celda de la soluci´on exacta obtenida con el esquema num´erico en t=tn: Wn i∼ =1 ∆xZxi+1/2 xi−1/2 W(x, tn)dx. Supongamos conocidos Wn i. Como es usual en los m´etodos de Godunov, aproximamos la soluci´on en t=tn+1 por Wn+1 i=1 ∆xZxi xi−1/2 Wi−1/2(x, tn+1)dx +Zxi+1/2 xi Wi+1/2(x, tn+1)dx, donde Wi+1/2denota la soluci´on del problema de Riemann asociado a los estados Wn iy Wn i+1 en la intercelda x=xi+1/2. Asumiendo una condici´on CFL 1/2 obtenemos la siguiente expresi´on para el m´etodo de Godunov: Wn+1 i=Wn i−∆t ∆xZ1 0 A(Φ(s, Wn i−1/2, Wn i))∂Φ ∂s (s;Wn i−1/2, Wn i)ds +Z1 0 A(Φ(s;Wn i, Wn i+1/2))∂Φ ∂s (s;Wn i, Wn i+1/2)ds,(7) donde Wn i+1/2es el valor (constante) de Wi+1/2en la intercelda x=xi+1/2. N´otese que Wi+1/2puede ser discontinua en x=xi+1/2. En tal caso, la discontinuidad en x=xi+1/2 es estacionaria, por lo que la masa asociada al salto correspondiente es cero. Por tanto, en el esquema podemos reemplazar Wn i+1/2bien por el l´ımite de Wi+1/2a la izquierda de xi+1/2,Wn,− i+1/2, o bien por el l´ımite a la derecha, Wn,+ i+1/2. 5 M.L. Mu˜noz, C. Par´es En el caso particular de un sistema de leyes de conservaci´on, en que Aes la matriz jacobiana de una funci´on de flujo F, (7) coincide con la expresi´on habitual del m´etodo de Godunov para sistemas de leyes de conservaci´on. Este m´etodo es bien equilibrado. Recordemos que la propiedad de buen equilibrado est´a relacionada con la aproximaci´on num´erica de las soluciones estacionarias. El sistema (1) admite soluciones estacionarias no triviales s´olo si tiene campos linealmente degenerados: si W(x) es una soluci´on estacionaria regular, satisface A(W(x)) ·W0(x) = 0,∀x∈R. Si W0(x)6= 0, entonces 0 es un autovalor de A(W(x)) con W0(x) como autovector asociado. Por tanto, esta soluci´on se puede interpretar como una parametrizaci´on de una curva integral de un campo linealmente degenerado cuyo autovalor correspondiente toma el valor 0 a lo largo de la curva. Si llamamos Γ al conjunto de todas las curvas integrales γde un campo linealmente degenerado de A(W) tal que el autovalor correspondiente se anula en Γ, dada una curva γ∈Γ, un esquema num´erico se dice que es exactamente bien equilibrado (o respectivamente bien equilibrado con orden k) para γsi resuelve de modo exacto (respectivamente con orden O(∆xk)) soluciones regulares estacionarias Wtales que W(x)∈γpara cada x. Se dice que el esquema num´erico es exactamente bien equilibrado (o bien equilibrado con orden k) si estas propiedades se satisfacen para cada curva de Γ (v´ease [6] para m´as detalles). 5. Consistencia del m´etodo de Godunov En [5] se conjetur´o un teorema de convergencia de tipo Lax-Wendroff para esquemas num´ericos camino-conservativos: si las soluciones obtenidas con un esquema caminoconservativo convergen a alguna funci´on cuando se refina la malla, entonces esta funci´on debe ser una soluci´on d´ebil del sistema no conservativo si las definiciones de esquema camino-conservativo y soluci´on d´ebil hacen referencia a la misma familia de camios. El resultado principal de esta secci´on apunta en ese sentido. Asumimos que ∆t/∆xes una constante fija λcuando ∆xy ∆ttienden a 0 y una condici´on CFL 1/2. Como discretizaci´on de la condici´on inicial elegimos las medias W0 i=1 ∆xZxi+1/2 xi−1/2 W0(x)dx. Sea h= ∆xy sea Wh(x, t) la funci´on definida c.p.d. en R×[0,∞) por Wh(x, 0) = W0 i, x ∈(xi−1/2, xi+1/2),(8) y Wh(x, t) = Wi+1/2(x, t),(x, t)∈(xi, xi+1)×[tn, tn+1).(9) La condici´on CFL impuesta asegura que las soluciones de los problemas de Riemann adyacentes no interfieren una con otra antes del instante de tiempo tn+1. 6 Sobre el m´etodo de Godunov para sistemas hiperb´olicos no conservativos Teorema Dado h > 0, sea Whla aproximaci´on num´erica (9) obtenida con el m´etodo de Godunov, y sea Wh(x, 0) dado por (8). Supongamos que existe una funci´on W∈(L∞(R×[0,∞)) ∩BV (R×[0,∞)))Ntal que kWh(·, t)−W(·, t)kL∞(R)N−→ h→00en L1([0,∞)) (10) y que existe una constante Ctal que TV (Wh(·, t)) ≤C∀t∈[0,∞), h > 0.(11) Entonces Wes una soluci´on d´ebil del problema (1) con condici´on inicial (6). La prueba de este teorema es una adaptaci´on del teorema de convergencia de LaxWendroff para sistemas hiperb´olicos de leyes de conservaci´on, para la que hemos de usar un resultado de estabilidad para productos no conservativos establecido en [2]. El resultado anterior no implica que la soluci´on d´ebil obtenida al refinar la malla sea una soluci´on entr´opica. Si asumimos el concepto de soluci´on d´ebil entr´opica basado en un par de entrop´ıa debemos probar que Z∞ 0Z∞ −∞ η(W(x, t))∂ϕ ∂t (x, t)dxdt −Z∞ 0Z∞ −∞ G(W(x, t))∂ϕ ∂x (x, t)dxdt ≤ − Z+∞ −∞ η(W0(x))ϕ(x, 0) dx ∀ϕ∈ C1 0(R×R+), ϕ ≥0.(12) Para ello ser´ıa suficiente suponer que se tiene una desigualdad de entrop´ıa discreta de la forma η(Wn+1 i)≤η(Wn i)−∆t ∆xGn i+1/2− Gn i−1/2,(13) donde Gn i+1/2=GWn i−q, . . . , Wn i+p,(14) siendo Guna funci´on Lipschitz continua de Ωp+q+1 en Ω, consistente con Gen el sentido de que G(W, . . . , W) = G(W)∀W∈Ω.(15) 6. Algunas consideraciones El teorema establecido en la secci´on anterior no puede ser considerado como una extensi´on del teorema cl´asico de Lax-Wendroff para problemas no conseervativos, ya que la convergencia uniforme aqu´ı requerida es mucho m´as fuerte que la exigida en dicho resultado cl´asico. Ha de ser pues entendido ´unicamente como un resultado de consistencia del esquema num´erico con el concepto de soluci´on d´ebil elegido y, por tanto, con las condiciones de Rankine-Hugoniot correspondientes. Este teorema se puede extender al caso de una familia de esquemas num´ericos basada en la elecci´on de un resolvedor de Riemann aproximado. En efecto, en [5] se introdujo una generalizaci´on para el caso de sistemas no conservativos del concepto de resolvedor de 7 M.L. Mu˜noz, C. Par´es Riemann aproximado, basada en la elecci´on de una familia de caminos, as´ı como la definici´on de esquema num´erico camino-conservativo, que generaliza la de esquema num´erico conservativo para problemas conservativos. Este concepto de resolvedor de Riemann aproximado para sistemas no conservativos permite la construcci´on de un esquema num´erico que puede ser escrito bajo la forma de un esquema camino-conservativo. El esquema de Godunov presentado es el caso particular correspondiente a la elecci´on del resolvedor de Riemann exacto. Agradecimientos Este estudio ha sido parcialmente financiado por los proyectos BFM2003-07530-C02-02 y MTM2006-08075. Referencias [1] M.J. Castro, J. Mac´ıas, C. Par´es. AQ-scheme for a class of systems of coupled conservation laws with source term. Application to a two-layer 1-D shallow water system. Math. Mod. Num. Anal., 35 (2001), 107–127. [2] G. Dal Maso, P.G. LeFloch, F. Murat. Definition and weak stability of nonconservative products. J. Math. Pures Appl., 74 (1995), 483–548. [3] P.D. Lax, B. Wendroff. Systems of conservation laws. Comm. Pure Appl. Math., 13 (1960) 217–237. [4] M.L. Mu˜noz-Ruiz, C. Par´es. Godunov method for nonconservative hyperbolic systems. Math. Mod. Num. Anal., 41 (2007), 169–185. [5] C. Par´es. Numerical methods for nonconservative hyperbolic systems: a theoretical framework. SIAM J. Numer. Anal., 44 (2006), 300–321. [6] C. Par´es, M.J. Castro. On the well-balance property of Roe’s method for nonconservative hyperbolic systems. Applications to shallow water systems. Math. Mod. Num. Anal., 38 (2004), 821–852. [7] A.I. Volpert. The space BV and quasilinear equations. Math. USSR Sbornik, 73 (1967), 225–267. 8