Un esquema de alto orden de tipo MUSTA para la resolución numérica de 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) Un esquema de alto orden de tipo MUSTA para la resoluci´on num´erica de sistemas hiperb´olicos no conservativos. M.J. Castro1, A. Pardo1, C. Par´ es1 1Depto. de An´alisis m´atematico, U. de M´alaga. E-mail: [email protected] Palabras clave: Sistemas hiperb´olicos no conservativos, bien equilibrado, MUSTA, GFORCE, WENO, shallow-water Resumen En este trabajo se presenta la extensi´on de los esquemas MUSTA (Multi-Stage) de alto orden a problemas no conservativos. En [7] se presentaron los esquema de tipo MUSTA para leyes de conservaci´on, usando como resolvedor de Riemann aproximado tanto en las etapas de predicci´on como de correcci´on el esquema GFORCE. Estos esquemas destacan por su simplicidad y su bajo coste computacional. En este trabajo formulamos el esquema GFORCE en el marco de los esquemas num´ericos Ψ-conservativos (“camino-conservativo”) introducidos por Par´es en [5]. El esquema MUSTA se reinterpreta como un esquema de reconstrucci´on de estados, donde el operador de reconstrucci´on usado est´a ligado a la resoluci´on de los problemas de Riemann asociados a cada intercelda. Finalmente se propone el uso de un operador de reconstrucci´on de estados de alto orden, que combinado con la estrategia MUSTA, resulta un esquema de alto orden de tipo MUSTA. Se presentan adem´as algunos ensayos num´ericos para el sistema de ecuaciones de aguas someras. 1. Introducci´on Consideremos el sistema hiperb´olico no conservativo: ∂W ∂t +A(W)∂W ∂x = 0, x ∈lR, t > 0,(1) donde la inc´ognita W(x, t) toma valores en un conjunto abierto convexo Ω ⊂lRN. Supondremos que el sistema es estrictamente hiperb´olico, esto es, ∀W∈Ω la matriz A(W) 1 M.J. Castro, A. Pardo, C. Par´es tiene Nautovalores reales distintos λ1(W)< . . . < λN(W). N´otese que (1) tiene soluciones estacionarias no triviales si posee alg´un campo linealmente degenerado asociado al autovalor 0. Los sistemas de leyes de conservaci´on con productos no conservativos y t´erminos fuente ∂w ∂t +∂F ∂x (w) = −B(w)∂w ∂x +S(w)dσ dx,(2) siendo σ(x) funci´on conocida, son un caso particular de (1). Efectivamente basta definir A(W) = J(w)−S(w) 0 0 ,(3) donde J(w) = ∂F ∂w (w) + B(w), W =w σ. Supondremos adem´as que J(w) tiene N−1 autovalores reales y distintos λ1(w)< . . . < λN−1(w),(4) con autovectores asociados rj(w), j= 1, . . . , N −1. Si ninguno de estos autovalores se anula, (1) es un sistema estrictamente hiperb´olico: A(W) tiene Nautovalores reales y distintos λ1(w), . . . , λN−1(w),0,(5) con autovalores asociados R1(W), . . . , RN(W),(6) dados por Ri(W) = ri(w) 0, i = 1, . . . , N −1; RN(W) = J(w)−1·S(w) 1.(7) Obs´ervese que el campo Nes linealmente degenerado y est´a asociado al autovalor 0. Las curvas integrales del mismo vienen dadas por el sistema de ecuaciones: dW ds =RN(W).(8) 2. Esquemas num´ericos Ψ-conservativos En el contexto de la aproximaci´on num´erica para los sistemas hiperb´olicos no conservativos (1), Par´es introdujo en [5] la siguiente definici´on: 2 Un esquema MUSTA de alto orden para problemas hiperb´olicos no conservativos. Definici´on 1 Dada una familia de caminos Ψ, un esquema num´erico es Ψ-conservativo si puede escribirse de la siguiente manera: Wn+1 i=Wn i−∆t ∆xD+ i−1/2+D− i+1/2,(9) donde D± i+1/2=D±Wn i−q, . . . , Wn i+p, siendo D−yD+dos funciones continuas de Ωp+q+1 aΩque cumplen D±(W, . . . , W) = 0,∀W∈Ω,(10) y D−(W−q, . . . , Wp) + D+(W−q, . . . , Wp) = Z1 0AΨ(s;W0, W1)∂Ψ ∂s (s;W0, W1)ds, para cada Wi∈Ω, i =−q, . . . , p. Esta definici´on generaliza a la de esquema num´erico conservativo para sistemas de leyes de conservaci´on, de hecho se demuestra que si (1) es un sistema de leyes de conservaci´on (esto es, si en (1) Aes la matriz jacobiana de una funci´on flujo F), entonces todo esquema Ψ-conservativo para alguna familia de caminos Ψ es consistente y conservativo en el sentido habitual. Rec´ıprocamente, un esquema conservativo en el sentido usual es Ψ-conservativo para cualquier familia de caminos Ψ. 2.1. Esquemas num´ericos bien equilibrados La aproximaci´on num´erica del equilibrio, esto es de las soluciones estacionarias, est´a´ıntimamente ligada a la propiedad de bien equilibrado. Definimos el conjunto Γ de todas las curvas integrales γde un campo linealmente degenerado de A(W) tal que el autovalor correspondiente es nulo en Γ. En [5] se demuestr´o que para obtener un esquema num´erico Ψ-conservativo bien equilibrado para una curva γ∈Γ, debe ocurrir: Z1 0A(Ψ(s;W0, W1))∂Ψ ∂s ds = 0 ∀W0, W1∈γ. (11) 3. Esquema GFORCE para sistemas hiperb´olicos no conservativos Presentamos en esta secci´on la formulaci´on de un esquema de tipo GFORCE para sistemas hiperb´olicos no conservativos. La deducci´on del mismo se detalla en [3]. En primer lugar, es necesario introducir el concepto de linealizaci´on de Roe propuesta por Toumi en el marco de sistemas hiperb´olicos no conservativos: Definici´on 2 Dada una familia de caminos Ψ, una funci´on AΨ: Ω ×Ω7→ MN(lR) es una linealizaci´on de Roe si cumple las siguientes propiedades: 3 M.J. Castro, A. Pardo, C. Par´es Para cada WL, WR,AΨ(WL, WR)tiene Nautovalores reales y distintos. ∀W, AΨ(W, W ) = A(W). ∀WL, WR,AΨ(WL, WR)(WL−WR) = Z1 0A[Ψ(s;WL, WR)]∂Ψ ∂s (s;WL, WR)ds. En [3] se propone la siguiente generalizaci´on del esquema GFORCE para sistemas hiperb´olicos no conservativos: Wn+1 i=Wn i−∆t ∆xD+ i−1/2+D− i+1/2, D+ i+1/2=1−ω 2 ∆x ∆tb Ii+1/2(Wi+1 −Wi) + ω 2 ∆t ∆xA2 i+1/2(Wi+1 −Wi)+ +1 2Ai+1/2(Wi+1 −Wi), D− i+1/2=−1−ω 2 ∆x ∆tb Ii+1/2(Wi+1 −Wi)−ω 2 ∆t ∆xA2 i+1/2(Wi+1 −Wi)+ +1 2Ai+1/2(Wi+1 −Wi), (12) donde ω=1 1 + CFL es un peso que garantiza propiedades de monoton´ıa del esquema num´erico. b Ii+1/2es una matritz, que se construye para garantizar el car´acter bien equilibrado del mismo. En [3] se muestra que si la matriz de Roe en el sentido de Toumi para Ψ es tal que el campo Nes linealmente degenerado y est´a asociado al autovalor 0, una elecci´on posible es: b Ii+1/2=Ki+1/2Id0 0 0 K−1 i+1/2(13) donde Ki+1/2tiene por columnas los autovectores de la matriz de Roe Ai+1/2=AΨ(Wi, Wi+1). Obs´ervese que en este caso 0 es autovector de b Ii+1/2y el autovector coincide con el de la matriz de Roe Ai+1/2. De esta forma es posible probar que el esquema GFORCE resultante es bien equilibrado. En particular, para leyes de equilibrio con productos no conservativos podemos usar la siguiente reformulaci´on: 4 Un esquema MUSTA de alto orden para problemas hiperb´olicos no conservativos. wn+1 i=wn i−∆t ∆xnD+ i−1/2+D− i+1/2o; D− i+1/2=1 2FGF wn i, wn i+1+Bi+1/2(wn i+1 −wn i)−Si+1/2(σn i+1 −σn i)+ −ω 2 ∆t ∆xJi+1/2+1−ω 2 ∆x ∆tJ−1 i+1/2Si+1/2(σn i+1 −σn i); D+ i+1/2=1 2−FGF wn i, wn i+1+Bi+1/2(wn i+1 −wn i)−Si+1/2(σn i+1 −σn i)+ +ω 2 ∆t ∆xJi+1/2+1−ω 2 ∆x ∆tJ−1 i+1/2Si+1/2(σn i+1 −σn i). (14) donde FGF =ωF LW + (1 −ω)FLF , FLW wn i, wn i+1=1 2F(wn i) + F(wn i+1)−1 2 ∆t ∆xJ2 i+1/2wn i+1 −wn i, FLF wn i, wn i+1=1 2F(wn i) + F(wn i+1)−∆x 2∆twn i+1 −wn i, (15) siendo Ji+1/2,Bi+1/2ySi+1/2los bloques correspondientes de la matriz de Roe Ai+1/2 dada en (3). Aunque formalmente necesitemos la inversa de la matriz Ji+1/2, en la pr´actica lo ´unico que usamos es el producto J−1 i+1/2Si+1/2. Cuando un autovalor de Ji+1/2se anula (problemas resonantes), se usar´a una aproximaci´on del producto J−1 i+1/2Si+1/2. 4. Esquemas MUSTA En [7] se describen los esquemas de tipo MUSTA (MUlti STAge). Se trata de un esquema de tipo predictor-corrector para aproximar la soluci´on del problema de Riemann en cada intercelda y definir as´ı el flujo num´erico en cada intercelda. Para que el esquema resultante sea eficiente desde el punto de vista computacional, en cada etapa del algoritmo se suele usar resolvedores de Riemann aproximados de bajo coste computacional como es el caso del esquema de tipo GFORCE. A fin de poder extender los esquemas MUSTA a problemas no conservativos, en [3] se interpretan como operadores de reconstrucci´on, esto es, cada etapa del algoritmo MUSTA se interpreta como una “reconstrucci´on”de la soluci´on a ambos lados de la intercelda. Finalmente, para obtener un esquema de alto orden podemos combinar la estrategia MUSTA con un operador de reconstrucci´on de alto orden, resultando un esquema MUSTA de alto orden. En este trabajo hemos usado el operador de reconstrucci´on PHM (piecewise hyperbolic method) (ver [4]), combinado con un esquema MUSTA de dos etapas, una de predicci´on y otra de correcci´on. 5. Ensayos num´ericos Consideramos en esta secci´on el sistema de ecuaciones de aguas someras unidimensional que modelan la evoluci´on de fluido homog´eneo no viscoso en un canal recto poco profundo 5 M.J. Castro, A. Pardo, C. Par´es de secci´on rectangular constante:          ∂h ∂t +∂q ∂x = 0, ∂q ∂t +∂ ∂x q2 h+g 2h2=ghdH dx . (16) La variable xhace referencia al eje del canal, tal tiempo, q(x, t) y h(x, t) representan al flujo m´asico y la altura de la columna de agua, H(x) la profundidad medida desde un nivel de referencia y ges la gravedad. Para aproximar num´ericamente las soluciones del problema anterior hemos usado un esquema MUSTA de dos etapas usando en cada etapa el esquema GFORCE dado en (14), donde Ji+1/2=0 1 ghi+1/2−u2 i+1/22ui+1/2, Bi+1/2= 0, Si+1/2=0 ghi+1/2, siendo hi+1/2=hi+hi+1 2yui+1/2=√hiui+phi+1ui+1 √hi+phi+1 . Estos estados intermedios corresponden a una linealizaci´on de Roe asociada a una familia de segmentos. En este caso particular, el producto J−1 i+1/2Si+1/2= 1 1−Fi+1/2 0!, donde Fi+1/2=u2 i+1/2 ghi+1/2. Cuando alguno de los autovalores de Ji+1/2se anula Fi+1/2= 1, con lo cual el producto J−1 i+1/2Si+1/2no est´a definido. Para evitar estos problemas, proponemos sustituir el vector J−1 i+1/2Si+1/2por signo(1 −Fi+1/2) 0. El esquema GFORCE as´ı construido verifica que es exactamente bien equilibrado para el agua en reposo, aproximando el resto de soluciones estacionarias con orden 1. Otra posibilidad ser´ıa usar una t´ecnica de reconstrucci´on hidrost´atica generalizada (ver [2], [3]). 5.1. Test num´ericos Test 1: El objetivo de este test es comprobar el orden de los esquemas ROE (v´ease [6]), GFORCE, MUSTA y PHM-MUSTA. Para ello consideramos un canal de 10 m de largo y una topografia dada por H(x) = 1 −0,5e−3(x−5)2sen2πx 10 . Como datos iniciales imponemos q= 0 y h(x) = H(x)+2−0,3(1−e−0,5(x−5)2) y dejamos evolucionar la soluci´on hasta t= 0,03. Se ha tomado CFL = 0,7. Hemos tomado como soluci´on de referencia una 6 Un esquema MUSTA de alto orden para problemas hiperb´olicos no conservativos. ROE GFORCE Celdas error horden herror qorder qerror horden herror qorder q 320 1,16E−24,24E−21,49E−25,31E−2640 5,84E−30,997 2,16E−20,973 6,84E−31,129 2,46E−21,109 1280 2,92E−31,000 1,09E−20,989 3,26E−31,070 1,18E−21,060 MUSTA PHM-MUSTA Celdas error horden herror qorder qerror horden herror qorder q 320 1,28E−24,69E−27,59E−44,39E−3640 6,04E−31,085 2,24E−21,065 1,21E−42,651 7,37E−42,575 1280 2,93E−31,045 1,09E−21,036 1,82E−52,732 1,11E−42,737 Tabla 1: Test 1, error y orden en t=0.03 soluci´on calculada con un mallado de 20480 puntos. Los errores en norma L1y el orden se muestran en la tabla 1. Test 2: En este test se comprueba como se comportan los esquemas ROE, GFORCE, MUSTA y PHM-MUSTA para tiempos grandes y la convergencia hacia una soluci´on estacionaria que posea una transici´on y un choque. Se toma en el intervalo [0,25] la funci´on de fondo H(x) =    0,05(x−10)2,if 8 < x < 12; 0,2,en caso contrario. Y los datos iniciales h= 0,33, q= 0,18, con condiciones de frontera q(0, t) = 0,18, h(25, t) = 0,33. Se ha tomado CFL= 0,9. Y el tiempo final es t= 200. Se ha calculado la soluci´on exacta con un mallado de 3200 puntos. La tabla 2 muestra el error en norma L1. ROE GFORCE Cells error horden herror qorder qerror horden herror qorder q 200 4,94E−35,09E−38,19E−23,41E−2400 1,55E−31,669 2,87E−30,825 4,64E−20,818 1,75E−20,964 800 6,98E−41,156 1,39E−31,049 2,77E−20,743 9,33E−30,905 MUSTA PHM-MUSTA Cells error horden herror qorder qerror horden herror qorder q 200 4,75E−21,93E−25,64E−36,41E−3400 2,75E−20,790 1,09E−20,821 2,05E−31,458 4,03E−30,669 800 1,46E−20,914 5,39E−31,020 7,36E−41,479 1,86E−31,113 Tabla 2: Test 2, error en t=200. Test 3: El siguiente test muestra el comportamiento de los esquemas ROE, GFORCE, MUSTA y PHM-MUSTA en presencia de frentes seco/mojado, usando para ello la t´ecnica 7 M.J. Castro, A. Pardo, C. Par´es descrita en [1]. Para ello usamos un test propuesto por Gallou¨et en que aparece una zona seca sobre un fondo escalonado. En el intervalo [0,25] definimos el fondo y los datos iniciales H(x) =    13,si 25/3< x < 25/2; 14,en otro caso. q(x) =    −300,si 50 3≤x; 300,si 50 3> x. yh(x) = H(x)−4. Se ha tomado CFL = 0,9 y ∆x= 0,125. La figura 1 muestran los resultados de los diferentes esquemas ante una soluci´on de referencia calculada sobre una mallado fino compuesto por 5000 vol´umenes en t= 0,25. 0 5 10 15 20 25 −14 −12 −10 −8 −6 −4 −2 ROE 0 5 10 15 20 25 −14 −12 −10 −8 −6 −4 −2 GFORCE 0 5 10 15 20 25 −14 −12 −10 −8 −6 −4 −2 MUSTA 0 5 10 15 20 25 −14 −12 −10 −8 −6 −4 −2 PHM−MUSTA Figura 1: Test 3, altura, comparaci´on con una soluci´on de refencia en t=0.25 Agradecimientos Este trabajo ha sido financiado parcialmente mediante el proyecto de investigaci´on MTM2006-08075. Referencias [1] M. Castro, A. Ferreiro, J.A. G´arc´ıa-Rodriguez, J. M. Gonz´alez, J. Mac´ıas, C. Par´es y M.E. V´azquezCendon. On the numerical treatment of wet/dry fronts in shallow flows: application to one-layer and two-layer 1-D shallow water system. Math. Comp. Model. 42 (3-4): 419-439, 2005. [2] M.J. Castro, A. Pardo, C. Par´es, Well-Balanced numerical schemes based on a generalized hydrostatic reconstruction technique. Mathematical Models and Methods in Applied Sciences. Aceptado en M3AS, 2007. [3] M.J. Castro, A. Pardo, C. Par´es, E.F. Toro, Non conservative well-balanced MUSTA schemes and high order extension. Submitted. [4] A. Marquina. Local piecewise hyperbolic reconstructions for nonlinear scalar conservation laws.. SIAM J. Sci. Comp 15:892-915, 1994. [5] C. Par´es. Numerical methods for nonconservative hyperbolic systems: a theoretical framework. SIAM J. Num. Anal. 44(1): 300-321, 2006. [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. M2AN, Vol. 38, N◦5, pp. 821-852, 2004. [7] E.F. Toro, V.A. Titarev, MUSTA Schemes for Systems of Conservation Laws. J. Comput. Phys., 216(2): 403 - 429, 2006. 8