Esquemas 2D de alto orden basados en reconstrucciones de estado, para sistemas hiperbólicos no conservativos. Aplicación a problemas de transporte de sedimentos

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Esquemas 2D de alto orden basados en reconstrucciones de estado, para sistemas hiperb´olicos no conservativos. Aplicaci´on a problemas de transporte de sedimentos. M.J. Castro D´ ıaz1, E.D. Fern´ andez Nieto2, A.M. Ferreiro Ferreiro1 1Dpto. de An´alisis Matem´atico, Univ. de M´alaga. Campus de Teatinos. 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. U. Arquitectura T´ecnica, Univ. de A Coru˜na. E-mail: [email protected]. Palabras clave: Reconstrucciones de estado, bien-equilibrado, m´etodo de vol´umenes finitos, transporte de sedimentos Resumen En este trabajo se aborda la aproximaci´on num´erica del problema de arrastre de sedimentos causada por la evoluci´on del agua. Para la componente hidrodin´amica se consideran las ecuaciones de aguas poco profundas. La componente morfodin´amica se define mediante una ecuaci´on de continuidad, dada en funci´on del caudal s´olido. Ambas componentes constituyen un sistema acoplado que puede reescribirse como un sistema hiperb´olico no conservativo (ver [4]). Se propone un esquema 2D generalizado de Roe con reconstrucciones de estado para sistemas hiperb´olicos no conservativos (ver [4]), mediante esquemas de vol´umenes finitos y el m´etodo de l´ıneas (ver [5]), extendiendo los esquemas de alto orden para el caso 1D propuesto en [2]. La reconstrucci´on de estado empleada es de tipo MUSCL (ver [1]), que proporciona orden dos para mallas no estructuradas de vol´umenes finitos de tipo arista. Finalmente se presenta un test num´erico en el que se estudia la evoluci´on del ´angulo de expansi´on de una monta˜na de arena (ver [4]). 1. Introducci´on En este trabajo estudiamos m´etodos de alto orden de vol´umenes finitos basados en reconstrucciones de estado. Dichos m´etodos se aplican a problemas de arrastre de sedimentos. La componente hidrodin´amica se modela mediante las ecuaciones de aguas poco 1 M.J. Castro D´ıaz, E.D. Fern´andez Nieto, A.M. Ferreiro Ferreiro profundas, y la componente morfodin´amica mediante una ecuaci´on de continuidad. El sistema final puede escribirse como un sistema 2D no conservativo ∂W ∂t +A1(W)∂W ∂x1 +A2(W)∂W ∂x2 = 0,(1) 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 η= (η1, η2)∈R2se define: A(W, η) = A1(W)η1+A2(W)η2. Suponemos que el sistema (1) es estrictamente hiperb´olico, es decir, para todo W∈Ω⊂ RNand ∀η∈R2, la matriz A(W, η) tiene Nautovalores reales y distintos: λ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. 2. Ecuaciones de aguas poco profundas con arrastre de sedimentos Consideramos el sistema,                            ∂h ∂t +∂q1 ∂x +∂q2 ∂x = 0, ∂q1 ∂t +∂ ∂x1q2 1 h+1 2gh2+∂ ∂x2q1q2 h=gh ∂S ∂x1−ghSf,x1, ∂q2 ∂t +∂ ∂x1q1q2 h+∂ ∂x2q2 2 h+1 2gh2=gh ∂S ∂x2−ghSf,x2, ∂S ∂t −ξ∂qb,x1 ∂x1−ξ∂qb,x2 ∂x2 = 0, (2) donde las inc´ognicas del problema son la altura de la columna del agua h(x, t), el caudal q(x, t)=(q1(x, t), q2(x, t)) y S(x, t) = H(x)−zb(x, t), donde Hes la batimetr´ıa y zbes el espesor de la capa de sedimento. Por qb,x1yqb,x2denotamos los caudales del transporte s´olido para las componentes x1yx2, respectivamente. La definici´on de qbdepende del modelo considerado. En este trabajo consideramos el modelo de Grass y el modelo de Meyer-Peter&M¨uller. Sf,x1ySf,x2son las leyes de fricci´on de Manning. Y ξ= 1/(1 −ρ0) donde ρ0es la porosidad de la capa de sedimento. El sistema de ecuaciones (2) puede escribirse como un sistema hiperb´olico no conservativo (1). La definici´on de qb,x1yqb,x2para los casos particulares de Grass y Meyer-Peter&M¨uller son: 1. Modelo de Grass: qb,x1=Agu1(u2 1+u2 2), qb,x2=Agu2(u2 1+u2 2); (3) 2 Esquemas 2D de alto orden para sistemas hiperb´olicos no conservativos donde u1=q1 hyu2=q2 h. Siendo Agla constante de interacci´on entre el fluido y el sedimento (0 ≤Ag≤1). 2. Modelo de Meyer-Peter&M¨uller: qb,x1= 8q(G−1)gd3 isgn(u1) (τ∗,x1−τ∗,c)3/2, qb,x2= 8q(G−1)gd3 isgn(u2) (τ∗,x2−τ∗,c)3/2; (4) donde τ∗yτ∗cson la tensi´on tangencial adimensional y la tensi´on tangencial cr´ıtica adimensional, respectivamente. La definici´on de τ∗es τ∗,x1=γη2|u1|pu2 1+u2 2 (γs−γ)diR1/3 h , τ∗,x2=γη2|u2|pu2 1+u2 2 (γs−γ)diR1/3 h , donde γdenota el peso espec´ıfico del fluido, γsdenota el peso espec´ıfico del sedimen- to, G=γs/γ,ηes el coeficiente de Manning correspondiente a la capa de sedimento, dies el di´ametro medio de grano. τ∗cse suele tomar igual a 0,047. Por Rhdenotamos el radio hidr´aulico. 3. M´etodos de vol´umenes finitos de alto orden basados en reconstrucciones de estado Consideremos un sistema hiperb´olico no conservativo ∂W ∂t +A1(W)∂W ∂x1 +A2(W)∂W ∂x2 = 0.(5) Nuestro objetivo es usar un m´etodo de vol´umenes finitos de alto orden basado en reconstrucciones de estado para (5). 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,1, ηij,2) es el vector unitario normal a la arista Eij y que apunta hacia la celda Vj(ver Figura 1). La discretizaci´on del sistema (5) se lleva a cabo mediante un esquema de vol´umenes finitos (ver [5]). Wn irepresentar´a la aproximaci´on Wn i∼ =1 |Vi|ZVi W(x, tn)dx. 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 [4]). Cuando los valores de la sucesi´on dependan del tiempo, denotaremos al operador de reconstrucci´on por Pt i. 3 M.J. Castro D´ıaz, E.D. Fern´andez Nieto, A.M. Ferreiro Ferreiro Figura 1: Volumen finito general. 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, s),(6) donde la arista Eij =Vi∩Vjest´a contenida en la recta, x·ηij =x1ηij,1+x2ηij,2=kij. Proponemos el siguiente esquema num´erico (ver [4]): W0 i(t) = −1 |Vi| X j∈Ni |Eij| n(¯r) X l=1 wlA− ij,l(W+ ij,l, W − ij,l, ηij)(W+ ij,l −W− ij,l) +ZViA1(Pi(x))∂Pi ∂x1 (x) + A2(Pi(x))∂Pi ∂x2 (x)dx, (7) donde wl,l= 1, . . . , n(¯r), denota los pesos de una f´ormula de cuadratura 1D, la cual viene asociada a una integral sobre la arista Eij. Si mediante xlse denotan los puntos sobre la arista Eij de la f´ormula de cuadratura, entonces W± ij,l =W± ij (xl). En la pr´actica esta f´ormula se elige en funci´on del orden del operador de reconstrucci´on de estado: si por ¯r denotamos el orden de la f´ormula de cuadratura, entonces ¯r > p, siendo pel orden de aproximaci´on del operador de reconstrucci´on. Mediante Aij se denota la matriz de Roe asociada al problema 1D no conservativo proyectado sobre la arista Eij. Los productos no conservativos A(W)Wxdel problema 1D proyectado no tienen sentido en el marco de la teor´ıa de distribuciones. En este trabajo usamos la teor´ıa desarrollada por Dal Maso, LeFloch y Murat donde los productos no conservativos se interpretan como medidas de Borel, siendo necesario para su definici´on la elecci´on de una familia de caminos en el espacio de fases. Tambi´en ser´a necesario la elecci´on de una familia de caminos para definir la matriz de Roe (ver [7]). 4 Esquemas 2D de alto orden para sistemas hiperb´olicos no conservativos 4. Reconstrucciones de orden dos sobre mallas no estructuradas bidimensionales Por simplicidad supongamos que tenemos vol´umenes finitos de tipo arista (ver [5]). El 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 [4]), definido como Pi(x) = Wi+p(x),con p(x) = ∇Wi(x−Ni).(8) 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 [4]), ∇Wi≈ 4 X j=1 |Tj|∇W|Tj 4 X j=1 |Tj| .(9) 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 mismo tiempo capture las regiones donde W(x) sea discontinua sin introducir oscilaciones, es necesario modificar el operador de reconstrucci´on (8), usando para ello una funci´on limitadora de pendiente (ver detalles en [4]). 5 M.J. Castro D´ıaz, E.D. Fern´andez Nieto, A.M. Ferreiro Ferreiro 5. Simulaci´on bidimensional de la evoluci´on de una duna de arena c´onica El siguiente test (ver [6]) estudia el comportamiento de la evoluci´on de una duna de arena c´onica en un canal de dimensiones 1000 m×1000 m. Las condiciones iniciales son (ver Figura 3), h(x, y, 0) = 10,1−zb(x, y, 0), qx(x, y, 0) = 10, qy(x, y, 0) = 0; y la capa inicial de sedimento es una duna de arena con forma c´onica, zb(x, y, 0) =      0,1 + sin2π(x−300) 200 sin2π(y−400) 200 si 300 ≤x≤500, 400 ≤y≤600, 0,1 en otro caso. En primer lugar, estudiamos la evoluci´on de la capa de sedimentos al modelizar el arrastre mediante la f´ormula de Grass (3), para ello empleamos el esquema num´erico (7) con reconstrucciones de estado tipo MUSCL. Consideramos una interacci´on d´ebil entre el fluido y la capa de sedimentos, en concreto, imponemos Ag= 0,001. La porosidad del sedimento se ha tomado igual a 0,4. La condici´on CFL es igual a 0,8. Empleamos una malla de 7600 vol´umenes finitos de tipo arista. Se imponen el caudal q= (10,0) y el espesor de la capa de sedimentos zb= 0,1 en la frontera correspondiente a la recta x= 0 y condici´on de salida libre en la contenida en la recta x= 1000. En las paredes laterales se impone la condici´on de deslizamiento q·η= 0 (ver Figura 3). Con estos datos se deja evolucionar el experimento hasta 100 horas. (a) Vista desde arriba. (b) Vista lateral. Figura 3: Condici´on inicial. En la Figura 5 se muestran distintos instantes de la evoluci´on de la capa de sedimentos. Se observa como la capa de sedimentos se deforma gradualmente hacia un patr´on con forma de estrella, expandi´endose a lo largo del tiempo con un cierto ´angulo (ver [6]). De Vriend obtiene en [3] una aproximaci´on anal´ıtica del ´angulo, asumiendo que la interacci´on de la 6 Esquemas 2D de alto orden para sistemas hiperb´olicos no conservativos Figura 4: Estimaci´on del ´angulo de expansi´on (empleando 20 curvas de nivel). capa de sedimentos y el fluido es lenta, esto es, Ag<0,01. La aproximaci´on anal´ıtica de dicho ´angulo de expansi´on viene dada por, tan α=3√3(mg−1) 9mg−1,por lo que, para mg= 3 el valor de dicho ´angulo de expansi´on es: α= tan−1 3√3 13 != 21,786789o.En la Figura 4 se presenta el ´angulo de expansi´on estimado al comparar diferentes instantes en la evoluci´on de la duna de arena. El ´angulo estimado es aproximadamente 23o. Agradecimientos Este trabajo ha sido parcialmente financiado por los proyectos de investigaci´on MTM2006- 08075 y MTM2006-01275, financiados por el Gobierno Espa˜nol. Referencias [1] B. Van Leer. MUSCL. A new approach to numerical gas dynamics. Computing in plasma physics and astrophysics, Max-Planck-Institut fur plama physik.. Carchung, Germany, April 1976. [2] M.J. Castro, J.M. Gallardo and 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. Comp., 75: 1103-1134. 2006. [3] H.J. De Vriend. 2DH Mathematical Modelling of Morphological Evolutions in Shallow Water Coastal Engineering 11, 1–27, 1987. [4] 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. [5] 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. [6] J. Hudson. Numerical technics for morphodynamic modelling. Tesis doctoral. University of Whiteknights, 2001. [7] C. Par´es, M.J. Castro. On the well-balanced property of Roe’s method for nonconservative hyperbolic systems. Applications to shallow water systems. ESAIM: M2AN 38(5): 821-852, 2004. 7 M.J. Castro D´ıaz, E.D. Fern´andez Nieto, A.M. Ferreiro Ferreiro (a) 0 seg. (b) 75000 seg. (c) 95000 seg (d) 185000 seg (e) 260000 seg (f) 360000 seg Figura 5: Evoluci´on de la capa de sedimentos en diferentes instantes. Vista de las curvas de nivel (20 curvas de nivel). 8