Un esquema de volúmenes finitos de alto orden para las ecuaciones de aguas someras con topografía y áreas secas

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matemtica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Un esquema de vol´umenes finitos de alto orden para las ecuaciones de aguas someras con topograf´ıa y ´areas secas J.M. Gallardo1, M. Castro1, C. Par´ es1 1Dpto. An´alisis Matem´atico, Universidad de M´alaga, Campus de Teatinos s/n, 29080-M´alaga. E-mails: [email protected], [email protected], [email protected]. Palabras clave: Sistemas hiperb´olicos, productos no conservativos, esquemas bien equilibrados, m´etodo de Roe, esquemas de alto orden, sistemas de aguas someras Resumen Presentamos un esquema de vol´umenes finitos para la resoluci´on de las ecuaciones de aguas someras con t´ermino fuente debido a la topograf´ıa del fondo. Se trata de un esquema de alto orden, bien equilibrado y capaz de afrontar situaciones en las que aparecen zonas secas. El esquema se ha desarrollado en un marco no conservativo general, y se basa en reconstrucciones hiperb´olicas de estados. El tratamiento de situaciones seco/mojado se lleva a cabo mediante la resoluci´on de problemas de Riemann no lineales en las interceldas donde se detecta una transici´on. 1. Introducci´on Las ecuaciones de aguas someras son ampliamente utilizadas para modelizar flujos en r´ıos, pantanos, ´areas costeras, etc. En la forma considerada en este trabajo, dichas ecuaciones conforman un sistema hiperb´olico de leyes de conservaci´on con un t´ermino fuente debido a la topograf´ıa del fondo. En los ´ultimos a˜nos ha habido un ingente desarrollo de esquemas num´ericos de alto orden para resolver las ecuaciones de aguas someras. Dichos esquemas calculan soluciones con un alto grado de precisi´on (tanto espacial como temporal) en las zonas donde dichas soluciones son regulares, mientras que, al mismo tiempo, las discontinuidades son adecuadamente capturadas. Sin embargo, cuando las ecuaciones contienen un t´ermino fuente, debe satisfacerse un cierto balance entre el flujo y el t´ermino fuente para que las soluciones estacionarias o casi estacionarias sean aproximadas correctamente. Los esquemas que satisfacen dicho balance se conocen con el nombre de esquemas bien equilibrados. Recientemente, en [1] se present´o un esquema de alto orden y bien equilibrado, desarrollado en un marco abstracto no conservativo. El esquema utiliza reconstrucciones de 1 J.M. Gallardo, M. Castro, C. Par´es estados y se basa en los esquemas de Roe generalizados introducidos en [8], y cuyas propiedades de bien equilibrado fueron analizadas en [6]. En particular, se aplic´o con ´exito a la resoluci´on de las ecuaciones de aguas someras con topograf´ıa, usando reconstrucciones WENO de quinto orden en espacio, y un esquema Runge-Kutta TVD de tercer orden para avanzar en tiempo. Una dificultad que aparece en la simulaci´on de flujos con superficie libre es la aparici´on de zonas secas, bien debidas a las condiciones iniciales o como resultado de la evoluci´on del fluido. Los ejemplos son numerosos: inundaciones, roturas de presas, tsunamis, etc. Si no se realizan modificaciones, los esquemas habituales pueden fallar ante la presencia de una situaci´on de seco/mojado, produciendo resultados sin significado f´ısico. En la literatura pueden encontrarse diversos m´etodos que solucionan el problema, aunque son a lo sumo de segundo orden. Cuando se aplican a las ecuaciones de aguas someras, los esquemas introducidos en [6] pierden sus propiedades de bien equilibrado ante la presencia de transiciones seco/mojado. De hecho, pueden producir valores negativos del espesor de la capa de agua en las proximidades del frente seco/mojado. Recientemente, en [2] se ha introducido una nueva t´ecnica para el tratamiento de frentes seco/mojado en el contexto de los esquemas de Roe generalizados. Dicha t´ecnica se basa en reemplazar, en las interceldas donde se ha detectado una transici´on seco/mojado, el correspondiente problema de Riemann lineal por otro no lineal, que se construye de acuerdo con el tipo de transici´on que aparezca. La idea central del presente trabajo consiste en combinar adecuadamente el esquema de alto orden presentado en [1] con la t´ecnica introducida en [2] para el tratamiento de frentes seco/mojado. Dicha combinaci´on no es en absoluto trivial, ya que aparecen numerosas dificultades. En particular, los flujos num´ericos deben ser modificados consistentemente en funci´on del tipo de transici´on seco/mojado que tengamos. Asimismo, las variables a reconstruir deben ser elegidas cuidadosamente, tanto para preservar el car´acter de bien equilibrado del esquema, como para mantener la positividad del espesor de la capa de agua. Por ´ultimo, hemos construido una extensi´on del esquema para resolver las ecuaciones de aguas someras bidimensionales sobre mallas estructuradas. 2. Ecuaciones de aguas someras, sistemas hiperb´olicos no conservativos y esquemas de Roe generalizados Consideremos el sistema unidimensional (∂th+∂xq= 0, ∂tq+∂xq2/h +gh2/2=ghH′(x),(1) que est´a formado por las ecuaciones que gobiernan un flujo de aguas someras a lo largo de un canal recto con secci´on rectangular constante. La variable xhace referencia al eje del canal y tes el tiempo; h(x, t) y q(x, t) son, respectivamente, el espesor y el caudal; g es la constante gravitatoria; por ´ultimo, H(x) es la funci´on de profundidad medida desde un nivel de referencia fijado. Se supone que el fluido es homog´eneo y no viscoso. Por simplicidad, los t´erminos que modelizan la fricci´on con el fondo y los efectos del viento no han sido considerados. 2 Un esquema de vol´umenes finitos para las ecuaciones de aguas someras Si a˜nadimos a (1) la ecuaci´on trivial ∂tH= 0, obtenemos el sistema no conservativo ∂tW+A(W)∂xW= 0, x ∈R, t > 0.(2) La variable de estado es W= [h, q, H]T, mientras que la matriz A(W) viene dada por A(W) =   1 0 0 −u2+c22u−c2 0 0 0  , donde u=q/h yc=√gh. Mientras no aparezcan zonas secas, la variable Wtoma sus valores en el conjunto abierto y convexo Ω = {[h, q, H]T:h > 0, q ∈R, H ∈R}. Aunque m´as adelante estaremos interesados en el caso en que Wpertenezca al conjunto Ω∗= Ω ∪{[0,0, H]T:H∈R}, supondremos en esta secci´on que se verifica h > 0. En general, el producto no conservativo A(W)Wxno est´a bien definido en el marco de la teor´ıa de distribuciones. De manera m´as precisa, el t´ermino ghHxno define una distribuci´on si tanto hcomo Hposeen discontinuidades en los mismos puntos. Sin embargo, la teor´ıa desarrollada en [3] permite interpretar dicho producto no conservativo como una medida de Borel. Muchos de los esquemas num´ericos dise˜nados para resolver sistemas de leyes de conservaci´on pueden ser adaptados para la discretizaci´on del sistema no conservativo (2). Este es el caso de los esquemas de Roe, cuya definici´on se basa en el concepto de linealizaci´on de Roe introducido en [8] (v´ease tambi´en [6]). Aqu´ı, nos restringiremos al caso particular en el que la matriz de Roe viene dada, para estados arbitrarios W0yW1en Ω, por A(W0, W1) =   0 1 0 −˜u2+ ˜c22˜u−˜c2 0 0 0  , donde ˜uy ˜cson los estados de Roe cl´asicos: ˜u=√h0u0+√h1u1 √h0+√h1 ,˜c=rgh0+h1 2.(3) Consideremos una descomposici´on en celdas Ii= [xi−1/2, xi+1/2], donde xi+1/2=i∆x; xi= (i−1/2)∆xes el centro de la celda Ii. Sea ∆tel paso temporal y definamos tn=n∆t. Denotemos por Wn ila aproximaci´on al promedio de la soluci´on exacta en la celda Ii, en el tiempo tn. Dados los estados Wn iyWn i+1, la matriz de Roe asociada se define como Ai+1/2=A(Wn i, Wn i+1). Entonces, el esquema num´erico puede escribirse como sigue: Wn+1 i=Wn i−∆t ∆xA+ i−1/2(Wn i−Wn i−1) + A− i+1/2(Wn i+1 −Wn i).(4) En lo que ata˜ne a la estabilidad, hay que considerar una condici´on de tipo CFL. Por ´ultimo, tambi´en se aplica el m´etodo de correcci´on entr´opica de Harten-Hyman. 3 J.M. Gallardo, M. Castro, C. Par´es Como es bien sabido, la presencia de t´erminos fuente o de acoplamiento en sistemas de leyes de conservaci´on puede afectar a la calidad de la soluci´on num´erica cuando se aproximan estados estacionarios o casi estacionarios. Los esquemas bien equilibrados, considerados por numerosos autores, solucionan dicho problema. Para una definici´on de la propiedad de bien equilibrado en un marco no conservativo, v´ease [6]. En el contexto de las ecuaciones de aguas someras (1), diremos que un esquema es bien equilibrado con orden γsi aproxima las soluciones estacionarias con orden γ, y es exactamente bien equilibrado si las calcula de forma exacta. En [6] se prueba que el esquema (4) es bien equilibrado con orden dos para soluciones estacionarias generales, y es exactamente bien equilibrado para soluciones correspondientes a agua en reposo. 3. Esquemas de alto orden Los esquemas de Roe generalizados introducidos en la secci´on anterior son tan s´olo de primer orden. Para capturar adecuadamente las caracter´ısticas de flujos complejos, es necesario considerar esquemas con mayor orden de precisi´on. En [1] se ha construido una extensi´on de alto orden del esquema (4), basada en reconstrucciones de los estados W. Veamos brevemente c´omo se construye dicho esquema. Fijado un tiempo t, se consideran funciones de reconstrucci´on de orden p,R−,t i+1/2(x) yR+,t i+1/2(x), definidas en [xi, xi+1/2] y [xi+1/2, xi+1] respectivamente. Sean W± i+1/2(t) las reconstrucciones laterales en la intercelda xi+1/2. Entonces el esquema puede escribirse en forma semidiscreta como sigue: W′ i(t) = −1 ∆xA+ i−1/2(W+ i−1/2(t)−W− i−1/2(t)) + A− i+1/2(W+ i+1/2(t)−W− i+1/2(t)) +Zxi xi−1/2A(R+,t i−1/2(x)) d dxR+,t i−1/2(x)dx +Zxi+1/2 xiA(R−,t i+1/2(x)) d dxR−,t i+1/2(x)dx,(5) donde Ai+1/2es la matriz de Roe asociada a los estados W− i+1/2(t) y W+ i+1/2(t). El esquema (5) resulta ser de alto orden (teorema 3.2 en [1]), exactamente bien equilibrado para soluciones de agua en reposo, y bien equilibrado con alto orden para soluciones estacionarias generales. Por ´ultimo, se considera un integrador en tiempo est´andar de orden adecuado. 4. El esquema MRoe para las ecuaciones de aguas someras En [2] se modific´o el esquema de Roe (4) para realizar un tratamiento adecuado de situaciones seco/mojado cuando se consideran las ecuaciones de aguas someras (1). En cada intercelda donde se detecta una transici´on, se sustituye el problema de Riemann aproximado habitual por un problema no lineal adecuado, que se resuelve de forma exacta. Los flujos num´ericos deben ser tambi´en modificados: v´ease [2] para los detalles. El esquema resultante se denomina esquema MRoe, y resulta ser bien equilibrado, incluso para soluciones estacionarias que contienen ´areas secas. Adem´as, se asegura la positividad del espesor hcalculado en las interceldas donde se ha detectado una transici´on seco/mojado, salvo en casos muy especiales, en los que ser´a necesario restringir la condici´on CFL. 4 Un esquema de vol´umenes finitos para las ecuaciones de aguas someras 5. Extensi´on de alto orden del esquema MRoe En esta secci´on presentamos un esquema bien equilibrado de alto orden para resolver (1), capaz de tratar adecuadamente situaciones de seco/mojado. Los ingredientes de dicho esquema han sido introducidos en las secciones precedentes. La idea es combinar el esquema de alto orden (5) con la t´ecnica del esquema MRoe para el tratamiento de zonas secas. El esquema resultante se denomina esquema HMRoe. Es necesario elegir adecuadamente las variables que ser´an reconstruidas. Para que el esquema resultante sea exactamente bien equilibrado para las soluciones de agua en reposo, consideraremos la elevaci´on η=h−H. Por otra parte, para evitar la aparici´on de valores negativos del espesor h, reconstruiremos dicha variable usando un operador que preserve la positividad. En resumen, las variables reconstruidas ser´an (h, q, η). Demos una breve descripci´on del esquema HMRoe. Supongamos conocidos los promedios Wi= [hi, qi, Hi]Ten un instante tdado. Entonces: Se definen los promedios ηi=hi−Hi. Si hi< hε, donde hεes una tolerancia fijada, hacemos hi= 0 y qi= 0 (esta puesta a cero se realiza tan s´olo en el proceso de reconstrucci´on: los valores originales de los promedios no son modificados). A partir de los promedios en una plantilla adecuada, reconstruimos las variables (h, q, η) y definimos el fondo reconstruido usando la igualdad H=h−η: esto nos proporciona los operadores de reconstrucci´on R± i+1/2(x) y los estados reconstruidos en las interceldas W± i+1/2. Si una plantilla contiene una celda seca, tomamos W− i+1/2= Wn ioW+ i+1/2=Wn i+1, seg´un el caso. En el caso de fondo mojado, los flujos num´ericos son los mismos que en (5). Si aparece una transici´on seco/mojado, se aplica la t´ecnica introducida en [2] para modificar tanto los estados intermedios como los flujos num´ericos. Una vez que el esquema semidiscreto ha sido definido, lo integramos en tiempo usando un esquema est´andar de orden apropiado. Seg´un se demuestra en [1] y [2], el esquema HMRoe es exactamente bien equilibrado para soluciones de agua en reposo (con o sin ´areas secas), y bien equilibrado con el mismo orden que el operador de reconstrucci´on utilizado, para soluciones estacionarias generales. En la pr´actica hemos usado reconstrucciones de tipo hiperb´olico (v´ease [5]). Las razones para esta elecci´on son m´ultiples: monotonicidad (evitamos la aparici´on de valores negativos del espesor h); tercer orden de precisi´on en toda la celda (el esquema resultante es de orden tres); plantilla compacta (robustez); y, por ´ultimo, el hecho de que las hip´erbolas poseen un mecanismo natural para reducir las oscilaciones cerca de un choque. Respecto a la discretizaci´on en tiempo, hemos utilizado un m´etodo Runge-Kutta TVD de tercer orden. 6. Extensi´on al caso bidimensional En [7] se ha desarrollado una extensi´on natural del m´etodo de reconstrucci´on hiperb´olica sobre mallas no estructuradas en dos dimensiones. Dicha extensi´on, denominada reconstrucci´on bi-hiperb´olica, considera en cada celda una combinaci´on de dos hip´erbolas, cada 5 J.M. Gallardo, M. Castro, C. Par´es una de ellas asociada a una direcci´on espacial. Cuando se aplica sobre mallas rectangulares (no necesariamente uniformes), se obtiene una reconstrucci´on de tercer orden en cada celda computacional. El m´etodo de reconstrucciones bi-hiperb´olicas permite extender el esquema HMRoe a un contexto bidimensional sobre mallas rectangulares, ya que la t´ecnica para el tratamiento de frentes seco/mojado puede aplicarse componente a componente. 7. Experimentos num´ericos En los experimentos num´ericos que siguen se ha tomado CFL = 0,8, hε= 10−6y g= 9,81. Los detalles de cada experimento pueden consultarse en [4]. 7.1. Generaci´on de fondo seco: la vara de Mois´es En este test consideramos un fondo plano con un obst´aculo rectangular. Inicialmente, la altura del agua es constante y el caudal es −350 si x < 50/3 y 350 si no. Dicha discontinuidad provoca la aparici´on, a la derecha del escal´on, de una zona seca entre dos ondas de rarefacci´on que viajan en direcciones opuestas. La onda que viaja hacia la izquierda interact´ua con el obst´aculo, produciendo la aparici´on de ondas en la superficie del agua. Los resultados obtenidos para diversos tiempos se representan en la figura 1. (a) Superficie libre. (b) Caudal. Figura 1: Generaci´on de un ´area seca sobre fondo no plano. 7.2. Drenado sobre fondo no plano En este experimento el fondo contiene un mont´ıculo e, inicialmente, el agua tiene altura constante. A la derecha se impone una condici´on de drenado. El flujo alcanza un estado estacionario, con agua a la izquierda del mont´ıculo y una zona seca a su derecha (v´ease la figura 2). 6 Un esquema de vol´umenes finitos para las ecuaciones de aguas someras (a) Superficie libre. (b) Caudal. Figura 2: Drenado sobre fondo no plano. 7.3. Inundaci´on en un canal abierto Consideramos un canal rectangular con tres mont´ıculos. Se supone que las fronteras superior e inferior son paredes s´olidas, y se impone un flujo entrante por la izquierda de la siguiente forma: hasta t= 300, h= 0,5 y u= 1,0; desde t= 300, h= 1,0 y u= 1,0. Los resultados obtenidos se representan en la figura 3, siendo bastante consistentes f´ısicamente. Agradecimientos Este trabajo ha sido parcialmente financiado por el proyecto BFM2003-07530-C02-02. Referencias [1] M. Castro, J.M. Gallardo, C. Par´es, High order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products. Applications to shallow water systems. Math. Comput. 75 (2006), 1103–1134. [2] M. Castro, J.M. Gonz´alez, C. Par´es, Numerical treatment of wet/dry fronts in shallow flows with a modified Roe scheme. Math. Mod. Meth. App. Sci. 16 (2006), 897–931. [3] G. Dal Maso, P. LeFloch, F. Murat, Definition and weak stability of nonconservative products. J. Math. Pures Appl. 74 (1995), 483–548. [4] J.M. Gallardo, C. Par´es, M. Castro, On a well-balanced high-order finite volume scheme for shallow water equations with topography and dry areas. Preprint (2007). [5] A. Marquina, Local piecewise hyperbolic reconstructions for nonlinear scalar conservations laws. SIAM J. Sci Comput. 15 (1995), 892–915. [6] C. Par´es, M. Castro, On the well-balanced property of Roe’s method for nonconservative hyperbolic systems. Applications to shallow water flows. ESAIM: M2AN 38 (2004), 821–852. [7] H.J. Schroll, F. Svensson, A bi-hyperbolic finite volume method on quadrilateral meshes. J. Sci. Comp. 26 (2006), 237–260. [8] I. Toumi, A weak formulation of Roe’s approximate Riemann solver. J. Comput. Phys. 102 (1992), 360–373. 7 J.M. Gallardo, M. Castro, C. Par´es Figura 3: Inundaci´on en un canal abierto. Resultados para t= 8, 30, 300 y 900. 8