Ajuste de campos de velocidad de aguas someras. El método de gradiente conjugado

Adjusting shallow water vector fields. The conjugate gradient method

^1*Jorge López-López, ²Lorenzo Héctor Juárez-Valencia

¹ División Académica de Ciencias Básicas. Universidad Juárez Autónoma de Tabasco, Km 1 carretera Cunduacán - Jalpa de Méndez, Cunduacán, Tabasco, México. C.P. 86690 Tel. (914)1001886 *jorge.lopez@ujat.mx

Recibido el 04 de marzo de 2014
Aceptado el 22 de agosto de 2014

RESUMEN

En este trabajo se muestra la aplicación al caso de campos de flujo en aguas someras de un modelo y un método de aproximación para ajustar campos de velocidades de viento bidimensionales. El modelo es de masa consistente y permite la reconstrucción de un campo vectorial solenoidal u a partir de un campo dado u¹. Matemáticamente el modelo se puede replantear como un problema de punto-silla, en donde se obtiene un problema de tipo Stokes degenerado para el campo de velocidad desconocido u con un multiplicador λ, asociado a la restricción de incompresibilidad. Este problema se expresa en forma operacional que se resuelve con el método de gradiente conjugado y el de elemento finito.

Palabras clave: Campo vectorial, gradiente conjugado, multiplicador de Lagrange, aguas someras, Saint-Venant, elemento finito.

ABSTRACT

Key words: Vector field, conjugate gradient, Lagrange multiplier, shallow water equations, Saint-Venant, finite element.

INTRODUCCIÓN

Para una amplia gama de aplicaciones de la dinámica de fluidos se requiere conocer, o al menos ajustar, campos de velocidades en una región dada ya sea en dos o tres dimensiones. Como ejemplos se tienen: el transporte y difusión de contaminantes en ríos, el efecto del viento en la propagación del fuego sobre estructuras, el análisis de calidad del aire (Cox et al. 2005, Wang et al. 2005, Chertock et al. 2006), y de particular interés el estudio de fenómenos asociados con el flujo en aguas someras, por su importancia para Tabasco, ya que su planicie está sujeta a inundaciones frecuentes, ocasionadas por flujos de aguas que pueden clasificarse como aguas someras (López et al. 2009).

Una forma de estudiar estos fenómenos es la modelación matemática y la simulación numérica de los mismos, en virtud de que este tipo de escenarios pueden ser modelados muy satisfactoriamente por las ecuaciones de Saint-Venant sobre un dominio espacial bidimensional Ω,

complementadas con condiciones de frontera y con una condición inicial u₀ = (u₁₀, u₂₀) y h₀ sobre Ω, donde las variables dependientes son el calado h, y los promedios verticales de las velocidades en las direcciones horizontales u₁, u₂, mientras que las variables dependientes y parámetros son la aceleración de la gravedad g, la fuerza de Coriolis C, el efecto del viento en la superficie del fluido τ, la pendiente del terreno S₀ y la fricción del fluido con el terreno S_f. En este modelo se desprecia la componente vertical, y las hipótesis principales son: i) La pendiente del fondo es pequeña, ii) El movimiento principal de las partículas ocurre en planos horizontales, iii) Las fuerzas de masa que actúan son la gravedad en la dirección vertical y la fuerza de Coriolis en el plano horizontal, iv) La curvatura de las líneas de corriente es pequeña, por lo que la distribución de la presión se considera hidrostática. Estas implican que la superficie del agua es gradualmente variada, sin cambios bruscos en su pendiente, y que las ecuaciones diferenciales parciales sean a lo sumo bidimensionales (López et al. 2009, Hodges 2013), es decir, en cada tiempo t, incluido el tiempo inicial, la velocidad u = (u₁, u₂) es un campo bidimensional.

En lo que respecta a la obtención de un campo de velocidad inicial, en la práctica es común contar solo con un número pequeño de mediciones aisladas dentro del dominio, exigiendo esto como primera necesidad, una interpolación de estos datos para extenderlos a todo el dominio, proceso que normalmente no cumple conservación de masa, siendo necesario un ajuste. En este contexto ajustar un campo significa que:

i) A partir de un campo vectorial dado se obtiene un campo vectorial (ajustado) que cumple conservación de masa (divergencia cero en todo punto del dominio);

ii) El campo ajustado no penetra las fronteras sólidas, sino solo las toca tangencialmente.

En la Figura 1 se muestra un ejemplo de dominio para una aplicación de viento y para una aplicación de aguas someras. Se ilustran ejemplos de forma del dominio, frontera artificial (el flujo la puede atravesar) y frontera natural (sólida). El dominio Ω se considera una region en R² con frontera , la cual está dada por en el caso de campos de viento, y para el caso de agua En ambos casos es la frontera natural, y el resto es la frontera artificial.

El objetivo de este trabajo es mostrar la aplicación de un modelo de masa consistente para ajuste de campos de viento (flujo no viscoso) bidimensionales, el cual ha sido aplicado en una variedad de problemas meteorológicos (Cox et al. 2005, Wang et al. 2005), al caso de ajustar campos de aguas someras (flujo viscoso y bidimensional), así como motivar la aplicación de la modelación matemática en el estudio del movimiento de aguas superficiales para la disminución de riesgo de inundaciones, y la prevención de los desastres asociados con este tipo de fenómenos.

MATERIALES Y MÉTODOS

Modelo y algoritmo de solución

i) Las características de u^I, ii) El conjunto factible de campos de velocidad donde se buscará a u, y iii) Una norma para medir cercanía entre u^I y los candidatos.

En lo que sigue se describe el modelo y el algoritmo numérico usado. Los detalles pueden ser consultados en Flores et al. (2010), definiendo el espacio de funciones factibles como:

y usando en H(div; Ω) la métrica

entonces el problema de ajuste de campos se puede formular como

se tiene que la u buscada es única y satisface

Para λ, en H¹ (Ω), λ = 0 en La integral que define la norma así como las demás que se mencionan en el texto son sobre todo Ω.

En Flores et al. (2010) se ha resuelto la ecuación (7) transformándola en un problema elíptico para y resolviendo un problema de punto silla para (u, λ). El enfoque de punto silla para resolver el problema (7), o equivalentemente (6), se basa en la metodología usual para resolver problemas de optimización sujeto a restricciones, haciendo uso del método de multiplicadores de Lagrange. El punto importante es relajar la condición de conservación de masa, introduciendo un nuevo parámetro: el multiplicador de Lagrange λ. Esto significa que se introduce un nuevo espacio de funciones vectoriales, más amplio que V, donde se busca la solución, el cual se define y denota como

Con la introducción de este nuevo espacio se considera sobre (V_N, L₂(Ω)) el Lagrangiano

Un punto estacionario de es un par (u, λ) que satisface las condiciones de primer orden, es decir,

donde λ no necesita satisfacer condiciones de frontera. En la solución (u,λ) de (10)-(11), u es el mínimo de J. Sea (u, A) solución del problema (10)-(11), y supóngase que el campo u es de la forma u = u¹ + u_λ, donde u¹ es el campo inicial y u_λ ∈ V_N representa el ajuste necesario que requiere u^I para ser la solución u. Con este supuesto, sustituyendo u = u^I + u_λ en (10) se llega a que el ajuste u_λ satisface la ecuación

Con esto, la solución del problema de punto silla (10)-(11) se ha traducido en resolver (12)-(13) para u_A. El problema (12)-(13) puede formularse como una ecuación funcional para u. Para esto se introduce el operador lineal A de L₂(Ω) en L₂(Ω) definido por

donde u_q ∈ V_N es la solución de

Con la definición de este funcional se tiene que el multiplicador A buscado satisface

El operador A es autoadjunto (simétrico) y fuertemente elíptico (definido positivo), así que para aproximar λ puede usarse gradiente conjugado, cuyo algoritmo, después de considerar la definición del operador A, es:

1. Inicio: Aλ⁰ ∈ L₂(Ω) dado, resolver para u⁰ ∈ V_N, (17)

Hacer

2. Descenso: Para k ≥ 0, conocidos λ^k, g^k, d^k y u^k calcular λ^{k + 1}, g^{k + 1} y u^{k + 1} como:

3. Paro o nueva dirección conjugada: Si tomar λ = λ^{k + 1}; u = u^{k + 1} y parar. En otro caso calcular d^{k + 1} = -g;^{k + 1} + β_kd^k : β_k = (20)

Hacer k = k + 1 y regresar al paso 2.

En el algoritmo anterior denota el producto interior usual en L₂(Ω). El algoritmo no requiere condición de frontera explícita para el multiplicador A (pero sí para las funciones tipo velocidad). Obsérvese que en el paso 2 se actualizan λ y u simultáneamente.

Las ecuaciones integrales del tipo (15), en el paso 1 y 2 del algoritmo de gradiente conjugado, se resuelven por medio del método de elemento finito con elementos lineales mixtos del tipo Bercovier-Pironneau (Bercovier y Pironneau 1979), en donde el multiplicador se aproxima por medio de elementos lineales en una malla dada y la velocidad también se aproxima por medio de elementos lineales pero en una malla doblemente fina como se muestra en la Figura 2; el objetivo de esto es satisfacer la condición inf-sup o condición LBB (Brezzi y Fortin 1991), para que la aproximación sea estable.

Se ha implementado un programa propio en Fortran con la metodología descrita y se ha aplicado para recuperar o ajustar tres campos de prueba. La descripción de esos campos y los resultados respectivos se muestran a continuación.

Campo u = (x,-y)

Para este ejemplo se tomó al campo inicial como u^I = (x,0) en Ω = (0,1) x (0,1). La frontera natural es la frontera inferior del cuadrado y la frontera artificial es el resto. Se tomó S = I. En la Tabla 1 se muestran los errores relativos y un promedio de la divergencia, obtenidas para diferentes mallas de presión (y sus respectivas mallas de velocidad), así como las iteraciones que se requirieron para que el algoritmo de gradiente conjugado alcanzara la convergencia. En las Figuras 3 y 4 se muestran el multiplicador y el contraste entre los campos exacto y ajustado, respectivamente, para la malla más gruesa usada.

Flujo invíscido con obstáculo

Se tomó aquí como campo de prueba a u = (u₁, u₂) donde: (21)

con α = 1.0, U₀ = 0.01 y como campo de inicio se tomó u¹ = (u, 0) en Ω = (—2, 2) x (0,1) menos el semicírculo con centro en (0,0) y radio uno. Se consideró S₁₁ = 1 y S₂₂ = 1. La frontera natural es la frontera inferior del dominio y la frontera artificial es el complemento. En este caso para la malla de presión más gruesa usada, h=0.42, se tuvo un error relativo de 8.91E-2, una divergencia de 4.24E-7 y se necesitaron 49 iteraciones para la convergencia de gradiente conjugado. En la Figura 5 se muestra el multiplicador y en la Figura 6 se contrastan el campo exacto con el campo aproximado, para la malla de presión más gruesa usada. Se consideraron condiciones de frontera exactas en la frontera natural y en las secciones verticales. Estos resultados están de acuerdo a los presentados en (Wang et al. 2005), quien usa una matriz S diferente a la identidad, resuelve un problema elíptico para el multiplicador.

Canal con paredes curvas

Para este ejemplo se tomó como campo base la solución numérica estacionaria u = (u₁,u₂) de un flujo de agua que entra horizontalmente en la parte izquierda (x = 0) y fluye a través de un canal con paredes curvas cuyas dimensiones son (0,20) x (0.3,10) menos los círculos con radio y centros en (10, 0.3) y (10,10), modelado con las ecuaciones de Saint-Venant en 2D. Este fenómeno involucra tanto las velocidades horizontales (u₁,u₂) como un calado h en cada punto del dominio, pero para este ejercicio solo se toma en cuenta la velocidad, olvidando el calado, ya que se está considerando un estado estacionario y en estas condiciones el calado no influye en la ecuación de continuidad. Como campo a ajustar se tomó u^I = (u^I, 0). Se consideró S₁₁ = 1, S₂₂ = 1. Aquí la frontera natural son las paredes del canal y el complemento es la artificial. Se consideraron condiciones exactas en toda la frontera. Para la malla de presión más gruesa usada, h=0.56, se tuvo un error relativo de 9.87E-2, una divergencia de -6.01E-5 con 100 iteraciones en gradiente conjugado. En las Figuras 7 y 8 se muestra el multiplicador y se contrastan el campo exacto con el campo aproximado en un acercamiento, respectivamente, para la malla de presión más gruesa usada. No se conocen resultados publicados relacionados con este ejemplo que sirvan de comparación. Sin embargo, los autores están trabajando en este mismo ejemplo, pero resolviendo un problema elíptico para A, y probando con diferentes condiciones de frontera para A. Cabe mencionar que con las ideas expuestas se está introduciendo la inquietud de recobrar campos de velocidad en el contexto de aguas someras.

En conclusión, se recuperó adecuadamente la segunda componente, a partir de la primera componente, de dos campos de velocidad que cuentan con solución analítica; esto ha servido para validación de la metodología y de los programas de cómputo, pues como muestra la tabla de convergencia, tanto el error relativo como la divergencia son excelentes y mejoran al refinar la malla. Para el ejemplo 1 Flores et al. (2010) reportan resultados similares para una malla de presión de 80x80. En Cervantes (2013) se reportan errores del mismo orden para este problema usando funciones de base radial, un método libre de mallas. Para el segundo ejemplo, los resultados son comparables cualitativamente con los presentados en (Wang et al. 2005) y de acuerdo al error y la divergencia obtenidos en éste trabajo, se ha ajustado bien a pesar de la complejidad de la frontera y la no linealidad de las componentes. En el tercer campo el ajuste se obtiene con precisión similar al ejemplo 2, lo que podría deberse a que el campo de prueba es resultado de una simulación numérica, el cual no tiene divergencia cero y además cambia abruptamente o casi discontinuamente de un nodo a otro. A pesar de ello, el campo ajustado suaviza el comportamiento del campo base y tiene una divergencia numérica muy pequeña. Además, es importante destacar que este tercer campo presenta dos áreas de recirculación que se muestran en los acercamientos en las partes derecha, superior e inferior, lo cual ha sido capturado también por el campo ajustado aun y cuando el modelo usado no considera ningún término asociado con viscosidad. Los dos últimos campos tienen características del dominio y del flujo similares a las de campos de aguas someras, y aunque el error relativo no ha sido tan espectacular como en el ejemplo 1, lo cual está asociado tanto con la complejidad del campo como del dominio, se considera que el método descrito da buenos resultados en este tipo de aplicaciones. En todos los casos el orden de la divergencia alcanzada es superior al obtenido por otras metodologías (Flores et al. 2010).

Resumiendo, se ha presentado un modelo y un algoritmo de solución que puede ser usado para recobrar campos de velocidad de flujos de aguas someras a partir de algún campo inicial. En general se pretende contar con una metodología eficiente para aplicaciones reales, como son los de aguas someras y los de los diversos campos experimentales de la dinámica de los fluidos y la meteorología. En dichos casos, el campo inicial ya no será un campo horizontal en el que se quiere recobrar la vertical, sino un campo vectorial perturbado en las distintas direcciones.

Bercovier M, Pironneau O (1979) Error estimates for the finite element method solution of the Stokes problem in the primitive variables. Numerical Mathematics 33: 211-224.         [ Links ]

Brezzi F, Fortin M (1991) Mixed and hybrid finite element methods, Springer-Verlag, New York, Inc. New York, NY, USA. 350 p.         [ Links ]

Cervantes DA, González-Casanova P, Gout C, Juárez LH, Reséndiz LR (2013) Vector field approximation using radial basis functions. Journal of Computational and Applied Mathematics 240: 163-173.         [ Links ]

Chertock A, Kurganov A, Petrova G (2006) Finite-Volume-Particule methods for models of transport of pollutants in shallow water. Journal of Scientific Computer 27: 189-199.         [ Links ]

Cox RM, Sontowski J, Dougherty CM (2005) An evaluation of three diagnostic wind models (CALMET, MCSCIPUF, and SWIFT) with wind data from the dipole pride 26 field experiments. Meteorological applications 12: 329-341.         [ Links ]

Flores C, Juárez H, Núñez MA, Sandoval ML (2010) Algorithms for vector field generation in mass consistent models. Numerical Methods for Partial Differential Equations 26: 826-842.         [ Links ]

Hodges BR (2013) Challenges in continental river dynamics. Environmental Modeling &. Software 50: 16-20.         [ Links ]

López J, Alavez J, Hernández JL (2009) Solución numérica de las ecuaciones de Saint-Venant vía volúmenes finitos, Revista de Ciencias Básicas 8: 34-53.         [ Links ]

Vázquez-Cendón ME, Cea L (2012) Analysis of a new Kolgan-type scheme motivated by the shallow water equations. Applied Numerical Mathematics 62: 489-506.         [ Links ]

Vázquez M (2008) Introducción al método de volúmenes finitos, Publicacións da USC. Colección Manuais Universitarios, vol. 10, España, 184 p.         [ Links ]

Wang Y, Williamson C, Garvey D, Chang S, Cogan J (2005) Application of a multigrid to a mass-consistent diagnostic wind model. Journal of Applied Meteorology 44: 1078-1089.         [ Links ]