Un esquema de alto orden de tipo MUSTA para la resolución numérica de sistemas hiperbólicos no conservativos

XX Cong eso de Ecuaciones Di e enciales y Aplicaciones X Cong eso de Ma em´ a ica Aplicada Se illa, 24-28 sep iemb e 2007 (pp. 1–8) Un esquema de al o o den de ipo MUSTA pa a la esoluci´on num´e ica de sis emas hipe b´olicos no conse a i os. M.J. Cas o1, A. Pa do1, C. Pa ´ es1 1Dep o. de An´alisis m´a ema ico, U. de M´alaga. E-mail: [email p o ec ed] Palab as cla e: Sis emas hipe b´olicos no conse a i os, bien equilib ado, MUSTA, GFORCE, WENO, shallow-wa e Resumen En es e abajo se p esen a la ex ensi´on de los esquemas MUSTA (Mul i-S age) de al o o den a p oblemas no conse a i os. En [7] se p esen a on los esquema de ipo MUSTA pa a leyes de conse aci´on, usando como esol edo de Riemann ap oximado an o en las e apas de p edicci´on como de co ecci´on el esquema GFORCE. Es os esquemas des acan po su simplicidad y su bajo cos e compu acional. En es e abajo o mulamos el esquema GFORCE en el ma co de los esquemas num´e icos Ψ-conse a i os (“camino-conse a i o”) in oducidos po Pa ´es en [5]. El esquema MUSTA se ein e p e a como un esquema de econs ucci´on de es ados, donde el ope ado de econs ucci´on usado es ´a ligado a la esoluci´on de los p oblemas de Riemann asociados a cada in e celda. Finalmen e se p opone el uso de un ope ado de econs ucci´on de es ados de al o o den, que combinado con la es a egia MUSTA, esul a un esquema de al o o den de ipo MUSTA. Se p esen an adem´as algunos ensayos num´e icos pa a el sis ema de ecuaciones de aguas some as. 1. In oducci´on Conside emos el sis ema hipe b´olico no conse a i o: ∂W ∂ +A(W)∂W ∂x = 0, x ∈lR, > 0,(1) donde la inc´ogni a W(x, ) oma alo es en un conjun o abie o con exo Ω ⊂lRN. Su- pond emos que el sis ema es es ic amen e hipe b´olico, es o es, ∀W∈Ω la ma iz A(W) 1 M.J. Cas o, A. Pa do, C. Pa ´es iene Nau o alo es eales dis in os λ1(W)< . . . < λN(W). N´o ese que (1) iene soluciones es aciona ias no i iales si posee alg´un campo lineal- men e degene ado asociado al au o alo 0. Los sis emas de leyes de conse aci´on con p oduc os no conse a i os y ´e minos uen e ∂w ∂ +∂F ∂x (w) = −B(w)∂w ∂x +S(w)dσ dx,(2) siendo σ(x) unci´on conocida, son un caso pa icula de (1). E ec i amen e bas a de ini A(W) = J(w)−S(w) 0 0 ,(3) donde J(w) = ∂F ∂w (w) + B(w), W =w σ. Supond emos adem´as que J(w) iene N−1 au o alo es eales y dis in os λ1(w)< . . . < λN−1(w),(4) con au o ec o es asociados j(w), j= 1, . . . , N −1. Si ninguno de es os au o alo es se anula, (1) es un sis ema es ic amen e hipe b´olico: A(W) iene Nau o alo es eales y dis in os λ1(w), . . . , λN−1(w),0,(5) con au o alo es asociados R1(W), . . . , RN(W),(6) dados po Ri(W) =  i(w) 0, i = 1, . . . , N −1; RN(W) = J(w)−1·S(w) 1.(7) Obs´e ese que el campo Nes linealmen e degene ado y es ´a asociado al au o alo 0. Las cu as in eg ales del mismo ienen dadas po el sis ema de ecuaciones: dW ds =RN(W).(8) 2. Esquemas num´e icos Ψ-conse a i os En el con ex o de la ap oximaci´on num´e ica pa a los sis emas hipe b´olicos no conse - a i os (1), Pa ´es in odujo en [5] la siguien e de inici´on: 2 Un esquema MUSTA de al o o den pa a p oblemas hipe b´olicos no conse a i os. De inici´on 1 Dada una amilia de caminos Ψ, un esquema num´e ico es Ψ-conse a i o si puede esc ibi se de la siguien e mane a: Wn+1 i=Wn i−∆ ∆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 unciones con inuas 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, pa a cada Wi∈Ω, i =−q, . . . , p. Es a de inici´on gene aliza a la de esquema num´e ico conse a i o pa a sis emas de leyes de conse aci´on, de hecho se demues a que si (1) es un sis ema de leyes de conse aci´on (es o es, si en (1) Aes la ma iz jacobiana de una unci´on lujo F), en onces odo esquema Ψ-conse a i o pa a alguna amilia de caminos Ψ es consis en e y conse a i o en el sen ido habi ual. Rec´ıp ocamen e, un esquema conse a i o en el sen ido usual es Ψ-conse a i o pa a cualquie amilia de caminos Ψ. 2.1. Esquemas num´e icos bien equilib ados La ap oximaci´on num´e ica del equilib io, es o es de las soluciones es aciona ias, es ´a´ın i- mamen e ligada a la p opiedad de bien equilib ado. De inimos el conjun o Γ de odas las cu as in eg ales γde un campo linealmen e degen- e ado de A(W) al que el au o alo co espondien e es nulo en Γ. En [5] se demues ´o que pa a ob ene un esquema num´e ico Ψ-conse a i o bien equilib ado pa a una cu a γ∈Γ, debe ocu i : Z1 0A(Ψ(s;W0, W1))∂Ψ ∂s ds = 0 ∀W0, W1∈γ. (11) 3. Esquema GFORCE pa a sis emas hipe b´olicos no con- se a i os P esen amos en es a secci´on la o mulaci´on de un esquema de ipo GFORCE pa a sis emas hipe b´olicos no conse a i os. La deducci´on del mismo se de alla en [3]. En p ime luga , es necesa io in oduci el concep o de linealizaci´on de Roe p opues a po Toumi en el ma co de sis emas hipe b´olicos no conse a i os: De inici´on 2 Dada una amilia de caminos Ψ, una unci´on AΨ: Ω ×Ω7→ MN(lR) es una linealizaci´on de Roe si cumple las siguien es p opiedades: 3 M.J. Cas o, A. Pa do, C. Pa ´es Pa a cada WL, WR,AΨ(WL, WR) iene Nau o alo es eales y dis in os. ∀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 p opone la siguien e gene alizaci´on del esquema GFORCE pa a sis emas hipe b´olicos no conse a i os: Wn+1 i=Wn i−∆ ∆xD+ i−1/2+D− i+1/2, D+ i+1/2=1−ω 2 ∆x ∆ b Ii+1/2(Wi+1 −Wi) + ω 2 ∆ ∆xA2 i+1/2(Wi+1 −Wi)+ +1 2Ai+1/2(Wi+1 −Wi), D− i+1/2=−1−ω 2 ∆x ∆ b Ii+1/2(Wi+1 −Wi)−ω 2 ∆ ∆xA2 i+1/2(Wi+1 −Wi)+ +1 2Ai+1/2(Wi+1 −Wi), (12) donde ω=1 1 + CFL es un peso que ga an iza p opiedades de mono on´ıa del esquema num´e ico. b Ii+1/2es una ma i z, que se cons uye pa a ga an iza el ca ´ac e bien equili- b ado del mismo. En [3] se mues a que si la ma iz de Roe en el sen ido de Toumi pa a Ψ es al que el campo Nes linealmen e degene ado y es ´a asociado al au o alo 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/2 iene po columnas los au o ec o es de la ma iz de Roe Ai+1/2=AΨ(Wi, Wi+1). Obs´e ese que en es e caso 0 es au o ec o de b Ii+1/2y el au o ec o coincide con el de la ma iz de Roe Ai+1/2. De es a o ma es posible p oba que el esquema GFORCE esul an e es bien equilib ado. En pa icula , pa a leyes de equilib io con p oduc os no conse a i os podemos usa la siguien e e o mulaci´on: 4 Un esquema MUSTA de al o o den pa a p oblemas hipe b´olicos no conse a i os. wn+1 i=wn i−∆ ∆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 ∆ ∆xJi+1/2+1−ω 2 ∆x ∆ J−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 ∆ ∆xJi+1/2+1−ω 2 ∆x ∆ J−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 ∆ ∆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∆ wn i+1 −wn i, (15) siendo Ji+1/2,Bi+1/2ySi+1/2los bloques co espondien es de la ma iz de Roe Ai+1/2 dada en (3). Aunque o malmen e necesi emos la in e sa de la ma iz Ji+1/2, en la p ´ac ica lo ´unico que usamos es el p oduc o J−1 i+1/2Si+1/2. Cuando un au o alo de Ji+1/2se anula (p oble- mas esonan es), se usa ´a una ap oximaci´on del p oduc o J−1 i+1/2Si+1/2. 4. Esquemas MUSTA En [7] se desc iben los esquemas de ipo MUSTA (MUl i STAge). Se a a de un esquema de ipo p edic o -co ec o pa a ap oxima la soluci´on del p oblema de Riemann en cada in e celda y de ini as´ı el lujo num´e ico en cada in e celda. Pa a que el esquema esul an e sea e icien e desde el pun o de is a compu acional, en cada e apa del algo i mo se suele usa esol edo es de Riemann ap oximados de bajo cos e compu acional como es el caso del esquema de ipo GFORCE. A in de pode ex ende los esquemas MUSTA a p oblemas no conse a i os, en [3] se in e p e an como ope ado es de econs ucci´on, es o es, cada e apa del algo i mo MUSTA se in e p e a como una “ econs ucci´on”de la soluci´on a ambos lados de la in e celda. Finalmen e, pa a ob ene un esquema de al o o den podemos combina la es a egia MUSTA con un ope ado de econs ucci´on de al o o den, esul ando un esquema MUSTA de al o o den. En es e abajo hemos usado el ope ado de econs ucci´on PHM (piecewise hype bolic me hod) ( e [4]), combinado con un esquema MUSTA de dos e apas, una de p edicci´on y o a de co ecci´on. 5. Ensayos num´e icos Conside amos en es a secci´on el sis ema de ecuaciones de aguas some as unidimensional que modelan la e oluci´on de luido homog´eneo no iscoso en un canal ec o poco p o undo 5 M.J. Cas o, A. Pa do, C. Pa ´es de secci´on ec angula cons an e:          ∂h ∂ +∂q ∂x = 0, ∂q ∂ +∂ ∂x q2 h+g 2h2=ghdH dx . (16) La a iable xhace e e encia al eje del canal, al iempo, q(x, ) y h(x, ) ep esen an al lujo m´asico y la al u a de la columna de agua, H(x) la p o undidad medida desde un ni el de e e encia y ges la g a edad. Pa a ap oxima num´e icamen e las soluciones del p oblema an e io hemos usado un esquema MUSTA de dos e apas usando en cada e apa 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 . Es os es ados in e medios co esponden a una linealizaci´on de Roe asociada a una amilia de segmen os. En es e caso pa icula , el p oduc o 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 au o alo es de Ji+1/2se anula Fi+1/2= 1, con lo cual el p oduc o J−1 i+1/2Si+1/2no es ´a de inido. Pa a e i a es os p oblemas, p oponemos sus i ui el ec o J−1 i+1/2Si+1/2po signo(1 −Fi+1/2) 0. El esquema GFORCE as´ı cons uido e i ica que es exac amen e bien equilib ado pa a el agua en eposo, ap oximando el es o de soluciones es aciona ias con o den 1. O a posibilidad se ´ıa usa una ´ecnica de econs ucci´on hid os ´a ica gene alizada ( e [2], [3]). 5.1. Tes num´e icos Tes 1: El obje i o de es e es es comp oba el o den de los esquemas ROE ( ´ease [6]), GFORCE, MUSTA y PHM-MUSTA. Pa a ello conside amos un canal de 10 m de la go y una opog a ia dada po H(x) = 1 −0,5e−3(x−5)2sen2πx 10 . Como da os iniciales imponemos q= 0 y h(x) = H(x)+2−0,3(1−e−0,5(x−5)2) y dejamos e oluciona la soluci´on has a = 0,03. Se ha omado CFL = 0,7. Hemos omado como soluci´on de e e encia una 6 Un esquema MUSTA de al o o den pa a p oblemas hipe b´olicos no conse a i os. ROE GFORCE Celdas e o ho den he o qo de qe o ho den he o qo de q 320 1,16E−2- 4,24E−2- 1,49E−2- 5,31E−2- 640 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 e o ho den he o qo de qe o ho den he o qo de q 320 1,28E−2- 4,69E−2- 7,59E−4- 4,39E−3- 640 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: Tes 1, e o y o den en =0.03 soluci´on calculada con un mallado de 20480 pun os. Los e o es en no ma L1y el o den se mues an en la abla 1. Tes 2: En es e es se comp ueba como se compo an los esquemas ROE, GFORCE, MUSTA y PHM-MUSTA pa a iempos g andes y la con e gencia hacia una soluci´on es a- ciona ia que posea una ansici´on y un choque. Se oma en el in e alo [0,25] la unci´on de ondo H(x) =    0,05(x−10)2,i 8 < x < 12; 0,2,en caso con a io. Y los da os iniciales h= 0,33, q= 0,18, con condiciones de on e a q(0, ) = 0,18, h(25, ) = 0,33. Se ha omado CFL= 0,9. Y el iempo inal es = 200. Se ha calculado la soluci´on exac a con un mallado de 3200 pun os. La abla 2 mues a el e o en no ma L1. ROE GFORCE Cells e o ho den he o qo de qe o ho den he o qo de q 200 4,94E−3- 5,09E−3- 8,19E−2- 3,41E−2- 400 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 e o ho den he o qo de qe o ho den he o qo de q 200 4,75E−2- 1,93E−2- 5,64E−3- 6,41E−3- 400 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: Tes 2, e o en =200. Tes 3: El siguien e es mues a el compo amien o de los esquemas ROE, GFORCE, MUSTA y PHM-MUSTA en p esencia de en es seco/mojado, usando pa a ello la ´ecnica 7 M.J. Cas o, A. Pa do, C. Pa ´es desc i a en [1]. Pa a ello usamos un es p opues o po Gallou¨e en que apa ece una zona seca sob e un ondo escalonado. En el in e alo [0,25] de inimos el ondo y los da os iniciales H(x) =    13,si 25/3< x < 25/2; 14,en o o caso. q(x) =    −300,si 50 3≤x; 300,si 50 3> x. yh(x) = H(x)−4. Se ha omado CFL = 0,9 y ∆x= 0,125. La igu a 1 mues an los esul ados de los di e en es esquemas an e una soluci´on de e e encia calculada sob e una mallado ino compues o po 5000 ol´umenes en = 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 Figu a 1: Tes 3, al u a, compa aci´on con una soluci´on de e encia en =0.25 Ag adecimien os Es e abajo ha sido inanciado pa cialmen e median e el p oyec o de in es igaci´on MTM2006-08075. Re e encias [1] M. Cas o, A. Fe ei o, J.A. G´a c´ıa-Rod iguez, J. M. Gonz´alez, J. Mac´ıas, C. Pa ´es y M.E. V´azquez- Cendon. On he nume ical ea men o we /d y on s in shallow lows: applica ion o one-laye and wo-laye 1-D shallow wa e sys em. Ma h. Comp. Model. 42 (3-4): 419-439, 2005. [2] M.J. Cas o, A. Pa do, C. Pa ´es, Well-Balanced nume ical schemes based on a gene alized hyd o- s a ic econs uc ion echnique. Ma hema ical Models and Me hods in Applied Sciences. Acep ado en M3AS, 2007. [3] M.J. Cas o, A. Pa do, C. Pa ´es, E.F. To o, Non conse a i e well-balanced MUSTA schemes and high o de ex ension. Submi ed. [4] A. Ma quina. Local piecewise hype bolic econs uc ions o nonlinea scala conse a ion laws.. SIAM J. Sci. Comp 15:892-915, 1994. [5] C. Pa ´es. Nume ical me hods o nonconse a i e hype bolic sys ems: a heo e ical amewo k. SIAM J. Num. Anal. 44(1): 300-321, 2006. [6] C. Pa ´es, M.J. Cas o. On he well-balance p ope y o Roe´s me hod o nonconse a i e hype bolic sys ems. Applica ions o Shallow-Wa e Sys ems. M2AN, Vol. 38, N◦5, pp. 821-852, 2004. [7] E.F. To o, V.A. Ti a e , MUSTA Schemes o Sys ems o Conse a ion Laws. J. Compu . Phys., 216(2): 403 - 429, 2006. 8