Esquemas numéricos bidimensionales de alto orden para resolución numérica del sistema de aguas someras con transporte inerte de un contaminante

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Esquemas num´ericos bidimensionales de alto orden para resoluci´on num´erica del sistema de aguas someras con transporte inerte de un contaminante M.J. Castro D´ ıaz1, E.D. Fern´ andez Nieto2, A.M. Ferreiro Ferreiro3, J.A. Garc´ ıa Rodr´ ıguez3, C. Par´ es Madro˜ nal1 1Dpto. de An´alisis Matem´atico, Univ. de M´alaga. E-mail: [email protected]. 2Dpto. de Matem´atica Aplicada I, Univ. de Sevilla. Sevilla. E-mail: [email protected]. 3Dpto. de Matem´aticas, E.P.S., Univ. de A Coru˜na. E-mail: [email protected]. Palabras clave: Reconstrucciones de estado, esquema bien-equilibrado, m´etodo de vol´umenes finitos, transporte de contaminantes, campos linealmente degenerados Resumen En este trabajo se aplican esquemas num´ericos bidimensionales de alto orden (ver [5]) empleando el m´etodo de vol´umenes finitos, para la modelizaci´on del transporte inerte de una sustancia en un fluido, como por ejemplo, el vertido de contaminantes. El modelo matem´atico consiste en un sistema acoplado de aguas someras y una ecuaci´on que modela de transporte inerte de una sustancia. Dicho acoplamiento introduce un nuevo campo linealmente degenerado en el sistema. Adem´as, si el vertido ocupa s´olo una porci´on de fluido, la frontera del mismo se propaga como una discontinuidad de contacto (ver [11]). En consecuencia, para aproximar con precisi´on la evoluci´on de un vertido es necesario desarrollar m´etodos num´ericos que capturen adecuadamente las discontinuidades de contacto. En este trabajo presentamos un m´etodo de vol´umenes finitos de alto orden bidimensional basado en t´ecnicas de reconstrucciones de estados. Concretamente hemos empleado una t´ecnica de tipo MUSCL discontinua que garantiza orden 2. Finalmente, se presentan algunos tests num´ericos. 1. Ecuaciones de aguas poco profundas con transporte de contaminantes El modelo matem´atico que modela el transporte de contaminantes resulta de acoplar las ecuaciones de aguas someras y una ecuaci´on de transporte: 1 M.J. Castro D´ıaz, E.D. Fern´andez Nieto, A.M. Ferreiro Ferreiro                              ∂h ∂t +∂qx ∂x +∂qy ∂x = 0, ∂qx ∂t +∂ ∂x µq2 x h+1 2gh2¶+∂ ∂y ³qxqy h´=gh∂H ∂x , ∂qy ∂t +∂ ∂x ³qxqy h´+∂ ∂y Ãq2 y h+1 2gh2!=gh∂H ∂y , ∂hC ∂t +∂qxC ∂x +∂qyC ∂y =σCσ, (1) donde las inc´ognicas del problema son la altura de la columna del agua h(x, t), el caudal q(x, t) = (qx(x, t), qy(x, t)) y la concentraci´on de contaminante C(x, t), donde Hes la batimetr´ıa del fondo medida desde un nivel de referencia, σ(x, t) son las fuentes emisoras (medida en m2/seg) y Cσla concentraci´on de sustancia en dichas fuentes. El sistema (1) puede escribirse como un sistema 2D no conservativo, ∂W ∂t +A1(W)∂W ∂x +A2(W)∂W ∂y = 0,(2) donde W(x, t) : O×(0, T)→Ω⊂RN,Oes un dominio acotado de R2, Ω es un conjunto convexo de RN,Ai: Ω →MN×Nfunciones regulares y localmente acotadas. Dado un vector unitario η= (ηx, ηy)∈R2se define: A(W, η) = Ax(W)ηx+A2(W)ηy. Suponemos que el sistema (2) es hiperb´olico tal que, para todo W∈Ω⊂RNy∀η∈R2, la matriz A(W, η) tiene Nautovalores reales: λ1(W, η)≤. . . ≤λN(W, η),siendo Rj(W, η), j= 1, . . . , N los autovectores asociados. En consecuencia A(W, η) es diagonalizable: A(W, η) = K(W, η)L(W, η)K−1(W, η),donde L(W, η) es la matriz diagonal cuyos coeficientes son los autovalores de A(W, η) y K(W, η) es la matriz cuyas columnas son los vectores Rj(W, η), j= 1, . . . , N. En este caso los autovalores del sistema son λ1=u·η−√ghkηk,λ2=u·η,λ3=u·η, λ4=u·η+√ghkηk. N´otese que el sistema acoplado (1) posee dos campos linealmente degenerados. 2. M´etodos de vol´umenes finitos de alto orden basados en reconstrucciones de estado Para discretizar el sistema (2), descomponemos el dominio computacional en celdas o vol´umenes de control, Vi⊂R2, que supondremos pol´ıgonos cerrados. Se usa la siguiente notaci´on: dado un vol´umen finito Vi,Nies el conjunto de ´ındices jtales que Vjes el vecino de Vi,Eij es la arista com´un de dos celdas vecinas ViyVj, y |Eij|su longitud, ηij = (ηij,x, ηij,y) es el vector unitario normal a la arista Eij y que apunta hacia la celda Vj(ver Figura 1). La discretizaci´on del sistema (2) se lleva a cabo mediante un esquema de vol´umenes finitos (ver [6]). Si W(x, t) es la soluci´on exacta denotaremos por Wn iel promedio de la soluci´on en tn, Wn i=1 |Vi|ZVi W(x, tn)dx, 2 Esquemas 2D de alto orden para acoplamiento de ecuaciones de transporte y ecuaciones de aguas someras Figura 1: Volumen finito general. Wn irepresentar´a la aproximaci´on de Wn ien tn, es decir, Wn i≃Wn i. Dado un volumen Videnotaremos por Piel operador de reconstrucci´on sobre el volumen. En concreto, Pidepender´a de una familia de valores {Wj}j∈Bi, donde Bies un conjunto de ´ındices de vol´umenes de control vecinos o pr´oximos a Vi(ver [5]). Cuando los valores de la sucesi´on dependan del tiempo, denotaremos al operador de reconstrucci´on por Pt i. Dado el vector ηij que apunta al volumen Vj, denotaremos por W− ij (t, s) y W+ ij (t, s) ∀s∈Eij, al l´ımite de Pt i(Pt jrespectivamente) cuando xtiende a sa trav´es de Vi(Vj respectivamente): l´ım x→s x·ηij < kij Pt i(x) = W− ij (t, s),l´ım x→s x·ηij > kij Pt j(x) = W+ ij (t, x).(3) 2.1. M´etodo 2D de vol´umenes finitos de alto orden Consideramos el siguiente esquema num´erico de alto orden, independiente del operador de reconstrucci´on empleado (ver detalles en [5]): W0 i(t) = −1 |Vi| X j∈Ni |Eij| n(¯r) X l=1 wlA− ij,l(W+ ij,l(t), W− ij,l(t), ηij)(W+ ij,l(t)−W− ij,l(t)) +ZViµA1(Pt i(x))∂Pt i ∂x (x) + A2(Pt i(x))∂Pt i ∂y (x)¶dx¸, (4) donde wl,l= 1, . . . , n(¯r), denota los pesos de una f´ormula de cuadratura asociada a la integral 1D sobre la arista Eij. Si mediante xlse denotan los puntos sobre la arista Eij de la f´ormula de cuadratura, entonces W± ij,l(t) = W± ij (t, xl). En la pr´actica esta f´ormula se elige en funci´on del orden del operador de reconstrucci´on de estado: si por ¯rdenotamos el orden de la f´ormula de cuadratura, entonces ¯r > p, siendo pel orden del operador de reconstrucci´on en las aristas de los vol´umenes finitos. Mediante Aij se denota la matriz de Roe asociada al problema 1D no conservativo proyectado sobre la arista Eij. El producto no conservativo A(W)Wxdel problema 1D proyectado hace dif´ıcil la definici´on de soluci´on d´ebil para esta clase de sistemas. En la teor´ıa desarrollada por Dal Maso, LeFloch y Murat (ver [4]) se introduce una definici´on de 3 M.J. Castro D´ıaz, E.D. Fern´andez Nieto, A.M. Ferreiro Ferreiro productos no conservativos como medidas de Borel, basada en la selecci´on de una familia de caminos en el espacio de fases. Tambi´en debe elegirse en este caso una familia de caminos para definir la matriz de Roe (ver [8]). En este trabajo hemos considerado un operador de reconstrucci´on de estado de orden dos para mallas no estructuradas, propuesto en [5]. 2.2. Reconstrucciones de orden dos sobre mallas no estructuradas bidimensionales Por simplicidad supongamos que tenemos vol´umenes finitos de tipo arista (ver [6]). Un volumen de tipo arista lo podemos escribir como Vi=Ti,1∪Ti,2∪Ti,3∪Ti,4, donde Ti,k,k= 1,2,3,4, son tri´angulos (ver Figura 2) definidos por Ci(punto medio de la arista sobre la que se construye el volumen finito de tipo arista) y las cuatro aristas del volumen de control. Sean bi,k,k= 1,2,3,4, los baricentros de estos tri´angulos, respectivamente. Figura 2: Tri´angulos que forman el volumen finito de tipo arista Vi. Se considera el siguiente operador de reconstrucci´on Pi(x) de orden dos (ver [5]), definido como Pi(x) = Wi+p(x),con p(x) = ∇Wi(x−Ni).(5) donde ∇Wies una aproximaci´on al menos de primer orden del gradiente de la soluci´on W(x), y Nial punto definido por: Ni= 4 X k=1 |Ti,k| |Vi|bi,k. Si por ∇W|Tjdenotamos la estimaci´on del gradiente sobre Tj, se considera la siguiente aproximaci´on, de segundo orden, del gradiente de la soluci´on en Ni(ver [5]), ∇Wi≈ 4 X j=1 |Tj|∇W|Tj 4 X j=1 |Tj| .(6) Es frecuente que en la aproximaci´on num´erica de sistemas hiperb´olicos como los considerados en este trabajo, la soluci´on presente discontinuidades. A fin de obtener un operador de reconstrucci´on que aproxime con orden dos la soluci´on en zonas regulares y que al 4 Esquemas 2D de alto orden para acoplamiento de ecuaciones de transporte y ecuaciones de aguas someras mismo tiempo capture las regiones donde W(x) sea discontinua, es necesario modificar el operador de reconstrucci´on (5), usando para ello una funci´on limitadora de pendiente. 0 10 20 30 40 50 60 70 0 0.1 0.2 0.3 0.4 0.5 Sol. exacta Roe Rec. orden dos (a) 0 seg. 0 10 20 30 40 50 60 70 0 0.1 0.2 0.3 0.4 0.5 Sol. exacta Roe Rec. orden dos (b) 4 seg. 0 10 20 30 40 50 60 70 0 0.1 0.2 0.3 0.4 0.5 Sol. exacta Roe Rec. orden dos (c) 5 seg. 0 10 20 30 40 50 60 70 0 0.1 0.2 0.3 0.4 0.5 Sol. exacta Roe Rec. orden dos (d) 8 seg. Figura 3: Evoluci´on de la concentraci´on de contaminante en diferentes instantes. Secci´on longitudinal central. 3. Experimentos num´ericos 3.1. Transporte de contaminante en un canal rectangular con fondo plano Este test bidimensional estudia la evoluci´on de una mancha de contaminantes de forma circular que se arrastra con velocidad constante en un canal rectangular de dimensiones 75 m×30 m. Las condiciones iniciales son: h(x, y, 0) = 2, qx(x, y, 0) = 10, qy(x, y, 0) = 0; y la concentraci´on inicial de contaminante viene dada por, C(x, y, 0) = (0,5 si (x−15)2+ (y−15)2≤36, 0 en otro caso. Suponemos adem´as que no existen fuentes emisoras. La condici´on CFL es igual a 0,8. Empleamos una malla estructurada de 9000 vol´umenes finitos. Se impone el caudal q= (10,0). Como condiciones de contorno se imponen condiciones libres en las fronteras correspondientes a la recta x= 0 y x= 75. En las paredes 5 M.J. Castro D´ıaz, E.D. Fern´andez Nieto, A.M. Ferreiro Ferreiro laterales se impone la condici´on de deslizamiento q·η= 0. Con estos datos se deja evolucionar el experimento hasta t= 10 segundos. La soluci´on exacta viene dada por un c´ırculo que se desplaza a velocidad constante en la direcci´on del eje X,vx=qx h= 5 m/s. En la Figura 3 se presenta una comparativa en la secci´on longitudinal central entre la soluci´on exacta, la soluci´on aproximada empleando un m´etodo de orden 1 (Roe) y la soluci´on aproximada mediante el esquema de la secci´on 2.2. Se aprecia como con el esquema de orden 2 se consigue mantener la concentraci´on inicial 0,5 y se preserva mejor la frontera que delimita el contaminante. En la Figura 4 se presenta una comparativa mediante curvas de nivel entre la difusi´on del contaminante empleando Roe y el m´etodo 2.2. 0 10 20 30 40 50 60 70 0 5 10 15 20 25 30 Esquema de orden 2 0 1 0 2 0 30 4 0 50 60 7 0 0 5 10 15 20 25 30 Roe (a) Tiempo 5 seg. 0 10 20 30 40 50 60 70 0 5 10 15 20 25 30 0 1 0 2 0 30 4 0 50 60 7 0 0 5 10 15 20 25 30 Esquema de orden 2 Roe (b) Tiempo 10seg. Figura 4: Evoluci´on del contaminante (de arriba abajo Roe vs. reconstrucci´on de segundo orden): curvas de nivel. 3.2. Transporte de contaminante en un canal rectangular con fondo variable Este test se lleva a cabo en un canal rectangular de dimensiones 75m×30m. El fondo del canal presenta un bump que viene dado por la siguiente funci´on: B(x, y) = e−0,075·(x−37,5)2.(7) Se supone inicialmente un contaminante ocupando un c´ırculo de centro (15,15) y radio 6 m. Suponemos inicialmente un flujo constante de q= (1,0) y una superficie libre constante de 3 m. de altura en su zona m´as profunda. Como condiciones de contorno se impone una condici´on de deslizamiento q·n= 0 en las fronteras y= 0 e y= 30, y condiciones libres en x= 0 y x= 75. Se considera CFL = 0,8 y se deja evolucionar la simulaci´on hasta el tiempo t= 120 seg. En la Figura 5 se presentan los resultados obtenidos empleando Roe 6 Esquemas 2D de alto orden para acoplamiento de ecuaciones de transporte y ecuaciones de aguas someras (a) Roe 65 seg. (b) Roe 120 seg. (c) Esq. de orden2, t=65 seg. (d) Esq. de orden2, t=120 seg. Figura 5: Evoluci´on de la concentraci´on de contaminante en diferentes instantes. 0 10 20 30 40 50 60 70 0 0.05 0.1 0.15 0.2 0.25 0.3 Roe Rec. orden dos (a) Tiempo 65 seg. 0 10 20 30 40 50 60 70 0 0.05 0.1 0.15 0.2 0.25 0.3 Roe Rec. orden dos (b) Tiempo 120seg. Figura 6: Evoluci´on de la concentraci´on de contaminante: secci´on longitudinal central (Roe vs. esquema de orden 2). y el esquema 2.2. Las figuras 5-(a) y 5-(b) corresponden al resultado obtenido empleando el m´etodo de Roe. En las figuras 5-(c) y 5-(d) se muestra la evoluci´on del contaminante empleando el esquema de reconstrucciones de orden 2. En la Figura 6 se muestra la comparativa entre el esquema de Roe y el esquema 2.2 a lo largo de una secci´on longitudinal central. En l´ınea punteada se muestra el resultado obtenido con el esquema de Roe y en l´ınea continua se muestra el resultado empleando el 7 M.J. Castro D´ıaz, E.D. Fern´andez Nieto, A.M. Ferreiro Ferreiro esquema de reconstrucciones de estado de orden 2. Agradecimientos Este trabajo ha sido parcialmente financiado por los proyectos de investigaci´on MTM200608075 y MTM2006-01275, financiados por el Gobierno Espa˜nol. Referencias [1] M.O. Bristeau, B.Perthame. Transport of pollutant in shallow water using kinetic schemes. CEMRACS 1999, ESAIM Proceedings 10, Soc. Math. Appl. Indust., Paris 1999, 9-21. [2] M.J. Castro, J.M. Gallardo, Carlos Par´es. Finite volume schemes based on weno reconstruction of states for solving nonconservaative hyperbolic systems. Applications to shallow water systems. Mathematics of Computation, 75(255), 1103–1134, 2006. [3] M. J. Castro, J. A. Garc´ıa , J. M. Gonz´alez, C. Par´es. A parallel 2D finite volume scheme for solving systems of balance laws with nonconservative products: application to shallow flows. Computer Methods in Applied Mechanics and Engineering. 195(19–22), 2006. [4] G. Dal Maso, P.G. LeFloch, F. Murat. Definition and weak stability of nonconservative products. J. Math. Pures Appl. 74: 483–548, 1995. [5] A.M. Ferreiro Ferreiro. Desarrollo de t´ecnicas de post-proceso de flujos hidrodin´amicos, modelizaci´on de problemas de transporte de sedimentos y simulaci´on num´erica mediante t´ecnicas de vol´umenes finitos. Tesis Doctoral. Universidad de Sevilla. 2006. [6] J.A. Garc´ıa Rodr´ıguez. Paralelizaci´on de esquemas de vol´umenes finitos: aplicaci´on a la resoluci´on de sistemas de tipo aguas someras. Tesis Doctoral. Universidad de M´alaga. 2005. [7] A. Marquina. Local piecewise hyperbolic reconstruction of numerical fluxes for nonlinear scalar conservation laws. SIAM, Journal Sci. Comput. 15(4), 892-915, 1994. [8] C. Par´es, M.J. Castro. On the well-balance property of Roe’s method for nonconservative hyperbolic systems. Applications to shallow-water systems. ESAIM: M2AN, 38(5):821–852, 2004. [9] H. J. Schroll and F. Svensson. A Bihyperbolic Finite Volume Method for Quadrilateral Meshes. SIAM: J. Sci. Comput., 26(2):237–260, 2006. [10] Susana Serna. A class of extended limiters applied to piecewise hyperbolic methods. SIAM: J. Sci. Comput., 28(1):123–140, 2006. [11] E. F. Toro. Shock-capturing methods for free-surface shallow flows. Wiley, 2001. [12] Zhengfu Xu, Chi-Wang Shu. Anti-diffusive finite difference weno methods for shallow water with transport of pollutant. Journal of Computational Mathematics, 24(3):239–251, 2006. 8