Approximation of the transport of pollutants with reaction terms in shallow waters using finite elements

Latest revision as of 15:46, 27 October 2020

Resumen

En este artículo presentamos la aproximación del modelo acoplado de las ecuaciones del movimiento de un fluido en aguas poco profundas con la ecuación convección-difusión-reacción (CDR) del transporte de contaminantes. Dicha aproximación se realiza mediante elementos finitos de alto orden y usando métodos variacionales estabilizados de subescalas. El sistema acoplado de ecuaciones, previamente discretizado en el tiempo y linealizado, lo escribimos como una ecuación vectorial transitoria de CDR. Los métodos estabilizados de elementos finitos utilizados son los conocidos métodos de subescalas ASGS y OSS, los mismos que nos permiten usar igual interpolación para todas las incógnitas, así como tratar con flujos de convección y reacción dominantes. Consideramos la posibilidad de no linealidad tanto en el término convectivo como en el de reacción. No consideraremos el posible desarrollo de choques en la solución. Con el fin de examinar la precisión y robustez de los métodos ASGS y OSS, presentamos cuatro casos de prueba: convergencia en malla, transporte de un contaminante en una cavidad cuadrada, transporte de un contaminante en el golfo de Roses y en la desembocadura del río Guadalquivir, y el modelo depredador-presa, que puede escribirse como una ecuación vectorial de CDR transitoria con no linealidad en el término de reacción.

Palabras clave: Aguas poco profundas, ablandamiento y estabilización, transporte de contaminantes, modelo depredador-presa

Abstract

In this article we present the approximation of the coupled model of the equations of motion of a fluid in shallow waters with the convection-diffusion-reaction (CDR) equation of pollutant transport. This approximation is carried out using high order finite elements and using stabilised variational sub-scale methods. We write the coupled system of equations, previously discretised in time and linearised, as a transient vector equation of CDR. The stabilised finite element methods used are the known ASGS and OSS sub-scale methods, the same ones that allow us to use the same interpolation for all unknowns, as well as to deal with dominant convection and reaction flows. We consider the possibility of non-linearity in both the convective and reaction terms. We will not consider the possible development of shocks in the solution. In order to examine the accuracy and robustness of the ASGS and OSS methods, we present four test cases: mesh convergence, transport of a pollutant in a square cavity, transport of a pollutant in the Gulf of Roses and at the river Guadalquivir mouth, and the predator-prey model, which can be written as a transient CDR vector equation with non-linearity in the reaction term.

Keywords: Shallow waters, softening and stabilization, transport of pollutants, predator-prey model

1. Introducción

Un contaminante es aquel componente que está presente en el agua a niveles perjudiciales para la vida de los seres humanos, plantas y animales ([1], véase también los informes de la U.S. Environmental Protecction Agency, https://www.epa.gov/). La simulación del transporte de contaminantes para la predicción de la concentración de contaminantes en ríos, lagos, lagunas y regiones costeras tiene importancia estratégica en el análisis y diseño de soluciones de los problemas de contaminación ambiental del planeta.

En el fenómeno físico del movimiento de contaminantes se deben distinguir tres procesos: la difusión, la convección y la reacción. La difusión es el proceso físico debido al cual el contaminante se mueve como resultado del movimiento intermolecular de las partículas de ambas sustancias, el fluido que la transporta (agua) y el contaminante. La convección es el movimiento del soluto (contaminante) debido al movimiento del agua, por lo cual si el agua permanece en reposo no hay convección. Finalmente, la reacción tiene en cuenta el posible efecto de crecimiento o decrecimiento de un contaminante por factores externos (calor, especies con las que puede combinarse); en el caso de haber distintas especies de contaminantes, la reacción puede modelar la interacción entre ellas.

La formulación del problema del transporte de contaminantes se fundamenta, como la de todos los fenómenos físicos, en las ecuaciones de equilibrio y las ecuaciones constitutivas. Las ecuaciones del transporte de sustancias disueltas o en suspensión en el flujo se basan en el principio de conservación de la masa de dichas sustancias en caso de que no haya reacción o en el modelo de crecimiento o decrecimiento en caso de que la haya. En el problema que queremos abordar, la dinámica del transporte de contaminantes involucra los modelos de flujos del campo de velocidades en aguas someras y del transporte, difusión y reacción de contaminantes, los cuales tienen características de la ecuación transitoria de convección, difusión y reacción (CDR).

En este trabajo vamos a utilizar un modelo matemático que acopla las ecuaciones del movimiento de un fluido en aguas poco profundas con la ecuación de CDR del transporte de contaminantes. Las ecuaciones del movimiento de un fluido en aguas poco profundas que vamos a utilizar son las obtenidas a partir de las ecuaciones de Navier-Stokes para fluidos viscosos, y que se encuentran descritas por ejemplo en [2], acopladas a la ecuación de CDR del transporte de contaminantes como en [3,4,5], entre otras referencias. El sistema acoplado de ecuaciones se puede escribir en forma conservativa en términos de la profundidad del fluido, de las velocidades y de las concentraciones de contaminantes promediadas en el plano donde el flujo medio tiene lugar. La presencia del término fuente en la ecuación de momento, la consideración de la batimetría variable, la fricción en el fondo del lecho, la tensión en la superficie libre del agua debida al viento y el efecto de Coriolis son características que se pueden tomar en cuenta para obtener modelos físicos más realistas del flujo de aguas someras. El resultado es un sistema no lineal acoplado de ecuaciones diferenciales parciales. Para simplificar el modelo, supondremos que los contaminantes son pasivos, en el sentido de que no afectan la densidad del fluido y por lo tanto el flujo; en caso de no ser así, se podría introducir fácilmente algún modelo de flotación que tuviera en cuenta el flujo generado por la variabilidad en la concentración de contaminantes.

En la literatura existe mucha investigación numérica en el área del transporte de contaminantes en aguas poco profundas, aplicando métodos de diferencias finitas [6], métodos de volúmenes finitos [7,8,11,12,13], métodos de elementos finitos [4,9,10], métodos tipo-Godunov [18] y métodos de Boltzmann [14,15,16].

En este trabajo vamos a realizar la aproximación de las ecuaciones del movimiento de un fluido en aguas poco profundas acopladas a la ecuación de CDR del transporte de contaminantes mediante los métodos VMS (Variational Muli-Scale) llamados ASGS (Algebraic SubGrid Scale) y OSS (Orthogonal Subgrid Scale) (véase [17]). El propósito de este trabajo es tanto aplicar estas técnicas, ya conocidas, a problemas no lineales como el transporte de contaminantes en aguas poco profundas y el modelo depredador-presa, como experimentar con elementos finitos de orden elevado (hasta el cuarto orden).

Veremos en primer lugar que mediante un cambio de variables apropiado y linealizando el sistema acoplado de ecuaciones en aguas poco profundas con las ecuaciones del transporte de contaminantes, cada iteración del proceso no lineal puede escribirse como una ecuación vectorial de CDR transitoria lineal. Plantearemos la aproximación en el tiempo y el espacio de esta ecuación vectorial de CDR transitoria, para luego particularizarla al sistema acoplado de las ecuaciones del flujo en aguas poco profundas con el transporte de contaminantes.

En lo concerniente a la discretización, primero procederemos a discretizar en el tiempo mediante un esquema de integración temporal en diferencias finitas y luego realizaremos la aproximación en el espacio mediante elementos finitos, aunque el enfoque más común es proceder a la inversa, es decir, primero discretizar en el espacio y luego aproximar el sistema resultante de ecuaciones diferenciales en el tiempo. Sin embargo, el enfoque de discretizar primero en el tiempo nos permite desacoplar los errores provenientes de la discretización temporal de los correspondientes a la discretización espacial. Aunque ambos caminos para el método estandar de Galerkin dan el mismo problema discreto, esto no es así para todos los métodos basados en VMS. En lo referente a la discretización temporal, usaremos la regla trapezoidal generalizada, que es el método en diferencias finitas más sencillo; sin embargo, cualquier otro método en diferencias finitas o incluso en elementos finitos podría aplicarse. En particular, en los ejemplos numéricos usamos esquemas de integración temporal de mayor orden que la regla trapezoidal.

En cuanto a la discretización espacial, es bien conocido que el método estándar de Galerkin presenta inestabilidades numéricas cuando los términos convectivo o de reacción son dominantes con respecto al término difusivo, además de por la necesidad dada por la condición inf-sup de usar diferentes interpolaciones de elementos finitos para la profundidad, las velocidades promediadas y la concentración de contaminante. Estos dos inconvenientes pueden superarse utilizando una formulación estabilizada. Aquí vamos a usar las formulaciones basadas en subescalas introducidas por Hughes y otros en [20,21] para la ecuación CDR escalar estacionaria y generalizado a sistemas lineales estacionarios en [22,23], y a sistemas lineales transitorios en [24,25,26], dando lugar a las formulaciones ASGS y OSS para resolver el problema de CDR vectorial transitorio, como en [27]. La idea básica es aproximar el efecto de la componente de la solución continua que no puede ser resuelta por la malla de elementos finitos en términos de la solución en el espacio de elementos finitos (método ASGS), o en términos de una subescala ortogonal respecto del producto escalar de ${\textstyle L^{2}}$ al espacio de elementos finitos [28,29] (método OSS). Hay que hacer notar que no consideraremos el posible desarrollo de choques en la solución, es decir, discontinuidades que pueden llegar a producirse en el límite invíscido. Para tomar en cuenta esto es necesario desarrollar técnicas de captura de discontinuidades que no consideraremos en este trabajo, así como tampoco el tratamiento de subescalas dinámicas.

El artículo está organizado como sigue. En la siguiente sección se describe la física del problema a ser resuelto y la obtención de la ecuación vectorial de CDR transitoria que representa el problema acoplado de las ecuaciones del movimiento de un fluido en aguas poco profundas con las ecuaciones del transporte de contaminantes. Las aproximaciones y formulaciones de los métodos variacionales estabilizados ASGS y OSS para la ecuación vectorial de CDR transitoria están descritos en la sección 3. En la sección 4 se presentan pruebas numéricas de convergencia en malla con soluciones analíticas conocidas, un ejemplo clásico de la literatura del transporte de un contaminante en una cavidad cuadrada, dos ejemplos prácticos del transporte de un contaminante, uno en el golfo de Roses y otro en la desembocadura del río Guadalquivir, y un ejemplo del modelo depredador-presa. Finalmente, en base a los resultados obtenidos, en la Sección 5 se detallan algunas conclusiones.

2. Planteamiento del problema del transporte de contaminantes

2.1 Ecuaciones de gobierno

Las ecuaciones del movimiento de un fluido en aguas poco profundas se obtienen por integración vertical de las ecuaciones del movimiento, bien sea de fluidos viscosos (ecuaciones de Navier-Stokes) o invíscidos (ecuaciones de Euler); véase por ejemplo [2,19]. En este trabajo vamos a utilizar las ecuaciones obtenidas a partir de las ecuaciones de Navier-Stokes para un fluido Newtoniano incompresible e isotérmico, incluyendo la tensión superficial provocada por los efectos del viento así como también los efectos de Coriolis, al igual que batimetría variable, fricción en el fondo del lecho y distribución hidrostática de presiones. Vamos a considerar el problema de mecánica de fluidos acoplado con el de transporte de contaminantes. Esta será en esencia la diferencia del presente trabajo con respecto a [27], por lo que muchos de los desarrollos que presentamos aquí se encuentran también en esta referencia. Omitiremos aquellos desarrollos que no afectan la compresión de este artículo, aunque para permitir la lectura auto-contenida se repetirán muchos otros.

Sea ${\textstyle \Omega }$ un dominio acotado ${\textstyle \mathbb {R} ^{2}}$ con coordenadas cartesianas ${\textstyle \mathbf {x} =(x_{1},x_{2}).}$ Consideremos un fluido que se mueve durante el intervalo de tiempo ${\textstyle [0,T]}$ en el dominio de ${\textstyle \mathbb {R} ^{3}}$ dado por ${\textstyle \Omega \times (-H(\mathbf {x} ),\eta (\mathbf {x} ,t)),}$ donde la coordenada ${\textstyle x_{3}=-H(\mathbf {x} )}$ representa la batimetría del terreno y ${\textstyle x_{3}=\eta (\mathbf {x} ,t)}$ , que será incógnita del problema, es la superficie libre del fluido con ${\textstyle \mathbf {x} \in \Omega }$ , ${\textstyle t\in [0,T]}$ . Consideraremos que ${\textstyle h(\mathbf {x} ,t):=H(\mathbf {x} )+\eta (\mathbf {x} ,t)}$ es mucho menor que el diámetro de ${\textstyle \Omega }$ , por lo que supondremos válida la aproximación de las ecuaciones de movimiento en aguas poco profundas. Dichas ecuaciones acopladas a la ecuación del transporte de contaminantes pueden escribirse en forma indicial como

$\partial _{t}\left(hU_{i}\right)+\partial _{j}\left(hU_{i}U_{j}\right)+\partial _{i}{\frac {1}{2}}g\left(h^{2}-H^{2}\right)-\partial _{j}\left({\frac {h}{\rho }}{\bar {\tau }}_{ji}\right)+Q_{i}=0,$	(1)
$\partial _{t}h+\partial _{i}\left(hU_{i}\right)=0,$	(2)
$\partial _{t}\left(h\varphi _{\alpha }\right)+\partial _{i}\left(hU_{i}\varphi _{\alpha }\right)-\partial _{i}\left(hk_{ij}^{\alpha }\partial _{j}\varphi _{\alpha }\right)+S_{\alpha \beta }\varphi _{\beta }=0,$	(3)

con

${\overline {\tau }}_{ij}=\mu _{H}\left(\partial _{j}U_{i}+\partial _{i}U_{j}-{\frac {2}{3}}\delta _{ij}\partial _{k}U_{k}\right),$	(4)
$Q_{i}=h{\hat {f}}\left(-\delta _{i1}U_{2}+\delta _{i2}U_{1}\right)-g\left(h-H\right)\partial _{i}H+{\frac {h}{\rho }}\partial _{i}p_{a}-{\frac {1}{\rho }}\tau _{3i}^{s}+gU_{i}U_{0}\mathrm {Ch^{-1}} h^{-2},\quad i=1,2,$	(5)

y

(6)

donde ${\textstyle U_{i}(\mathbf {x} ,t)}$ ${\textstyle \left(i=1,2\right)}$ son las componentes de la velocidad promediadas en la profundidad, ${\textstyle \delta _{ij}}$ es la delta de Kronecker, ${\textstyle \partial _{t}}$ denota la derivada temporal y ${\textstyle \partial _{i}}$ la derivada con respecto a la ${\textstyle i}$ -ésima coordenada cartesiana ${\textstyle (i=1,2)}$ en ${\textstyle \Omega }$ . En las ecuaciones anteriores y en lo que sigue, índices repetidos indican la notación de sumación de Einstein, y los subíndices ${\textstyle i,}$ ${\textstyle j,}$ ${\textstyle k}$ recorren desde 1 hasta 2. En las expresiones previas, ${\textstyle g}$ es la aceleración de la gravedad, ${\textstyle \rho }$ la densidad del fluido, ${\textstyle {\bar {\tau }}_{ij}}$ son las tensiones viscosas promediadas, ${\textstyle \mu _{H}}$ es el coeficiente de viscosidad dinámico, ${\textstyle {\hat {f}}}$ es el parámetro de Coriolis, ${\textstyle p_{a}}$ es la presión atmosférica, ${\textstyle \tau _{3i}^{s}}$ ${\textstyle \left(i=1,2\right)}$ son las tensiones en la superficie libre del agua debidas al viento, Ch es el coeficiente de Chézy para el rozamiento en el fondo y ${\textstyle U_{0}}$ es la norma euclídea de ${\textstyle \left(U_{1},U_{2}\right)}$ . Usaremos letras griegas como subíndices para caracterizar los contaminantes, cuyo número supondremos que es ${\textstyle N}$ . Así, la concentración del contaminante ${\textstyle \alpha }$ promediada en la profundidad se representa por ${\textstyle \varphi _{\alpha }}$ . ${\textstyle S_{\alpha \beta }}$ son las componentes de la matriz de términos de reacción de la ecuación de conservación de contaminantes, con ${\textstyle \alpha }$ , ${\textstyle \beta =1,\dots ,N}$ ; en los casos de reacción no lineal, ${\textstyle S_{\alpha \beta }}$ puede depender de las concentraciones de contaminantes, es decir, ${\textstyle S_{\alpha \beta }=S_{\alpha \beta }(\varphi _{1},\dots ,\varphi _{N})}$ . ${\textstyle k_{ij}^{\alpha }}$ es el coeficiente empírico de difusión (dispersión) para el contaminante ${\textstyle \alpha }$ ; nótese que por simplificar hemos supuesto que la dispersión del contaminante ${\textstyle \alpha }$ no se ve afectada por el resto de contaminantes, aunque esta situación podríamos incluirla fácilmente.

Las ecuaciones (1), (2) y (3) representan un sistema de ${\textstyle (3+N)}$ ecuaciones, cuyo vector de incógnitas es

Dichas ecuaciones constituyen el modelo acoplado del movimiento de un fluido en aguas poco profundas y el transporte de contaminantes. Las mismas se deben completar con las condiciones iniciales y de contorno. Los valores iniciales para la elevación del agua, las velocidades promediadas y la concentración del contaminante deben ser conocidos como condiciones iniciales. Con respecto a las condiciones de contorno, es bien conocido que el problema está bien puesto si la velocidad promediada se especifica en el contorno de entrada del dominio computacional ${\textstyle \Omega }$ , mientras que en el resto del contorno se prescribe la componente normal de las tensiones viscosas y la elevación del fluido. Para las concentraciones de contaminantes, además de la condición inicial es necesario conocer como condiciones de contorno sus valores en el contorno de entrada del dominio de cálculo y los flujos en el contorno de salida; esta situación permitiría en particular considerar nulas las difusiones de contaminantes.

Debido a la complejidad del modelo descrito anteriormente, en la literatura existen muchas variantes de estas ecuaciones con distintas simplificaciones, como por ejemplo eliminar el término de difusión considerando convección pura; este es el caso de [3,4,5]. En nuestro trabajo utilizaremos las ecuaciones (1), (2) y (3) sin eliminar o hacer simplificaciones de ninguno de sus términos. Mediante un cambio de variables transformaremos las ecuaciones para poder escribirlas como una ecuación vectorial de CDR transitoria de la forma

(7)

para la que más adelante, en la Sección 3, presentaremos su aproximación numérica por elementos finitos mediante los métodos variacionales de subescalas ASGS y OSS. En las ecuaciones anteriores, ${\textstyle \mathbf {u} }$ es el vector de incógnitas, ${\textstyle \mathbf {F} }$ el vector de fuerzas conocido y ${\textstyle \mathbf {K} _{ij}}$ , ${\textstyle \mathbf {A} _{i}}$ , ${\textstyle \mathbf {S} }$ matrices dadas, que pueden depender de ${\textstyle \mathbf {u} }$ en problemas no lineales.

2.2 Cambio de variables

Realizando un cambio de variables, ${\textstyle \varphi _{a,\alpha }:=h\varphi _{\alpha }}$ , ${\textstyle u_{i}:=hU_{i}}$ y ${\textstyle P:={\frac {1}{2}}g\left(h^{2}-H^{2}\right)}$ , con ${\textstyle \partial _{t}P=hg\partial _{t}h,}$ las ecuaciones (1), (2) y (3) se pueden escribir de la siguiente manera:

$\partial _{t}u_{i}+\partial _{j}\left(U_{j}u_{i}\right)+\partial _{i}P-\partial _{j}{\frac {\mu _{H}}{\rho }}\left[\left(\partial _{i}u_{j}+\partial _{j}u_{i}-{\frac {2}{3}}\delta _{ij}\partial _{k}u_{k}\right)-\left(U_{j}\partial _{i}h+U_{i}\partial _{j}h-{\frac {2}{3}}\delta _{ij}U_{k}\partial _{k}h\right)\right]$
$\qquad +{\hat {f}}\left(-\delta _{i1}u_{2}+\delta _{i2}u_{1}\right)+gU_{0}\mathrm {Ch} ^{-1}h^{-3}u_{i}-g\left(h-H\right)\partial _{i}H+{\frac {h}{\rho }}\partial _{i}p_{a}-{\frac {1}{\rho }}\tau _{3i}^{s}=0,$	(8)
${\frac {1}{hg}}\partial _{t}P+\partial _{i}u_{i}=0,$	(9)
$\partial _{t}\varphi _{a,\alpha }+\partial _{i}\left(U_{i}\varphi _{a,\alpha }\right)-\partial _{i}\left(k_{ij}^{\alpha }\partial _{j}\varphi _{a,\alpha }-k_{ij}^{\alpha }\varphi _{\alpha }\partial _{j}h\right)+S_{\alpha \beta }h^{-1}\varphi _{a,\beta }=0.$	(10)

Las ecuaciones (8), (9) y (10) se pueden escribir fácilmente en el formato de la ecuación (7). Sin embargo esto lo haremos en el caso del problema linealizado y discretizado en el tiempo.

2.3 Discretización en el tiempo y linealización

Para la discretización temporal del modelo acoplado de las ecuaciones del movimiento de un fluido en aguas poco profundas y del transporte de contaminantes, consideremos una partición uniforme del intervalo de tiempo ${\textstyle 0=t^{0}<t^{1}<\dots {<}t^{M}=T}$ , con ${\textstyle \delta t=t^{n+1}-t^{n}}$ , donde ${\textstyle \delta t}$ es el tamaño del paso de tiempo, considerado constante. Aquí y en adelante, usaremos el superíndice ${\textstyle n}$ para indicar el nivel de tiempo ${\textstyle n}$ .

Aunque la discretización en el tiempo puede realizarse mediante cualquier aproximación de la derivada temporal en diferencias finitas o elementos finitos, aquí utilizaremos la regla trapezoidal generalizada. Para una función genérica ${\textstyle f}$ , sea

(11)

Aproximaremos

Esta aproximación es de segundo orden en ${\textstyle \delta t}$ en el caso en que ${\textstyle \theta =1/2}$ (método de Crank-Nicolson) y de primer orden si ${\textstyle \theta \not =1/2}$ . Para ${\textstyle \theta =1}$ corresponde al método de Euler implícito o el método de diferencias hacia atrás de primer orden (BDF1).

En cuanto a la linealización, utilizaremos el método de punto fijo de Picard. En cada iteración de cada paso de tiempo, el problema a resolver es de la forma:

$\delta _{t}^{n}u_{i}+\partial _{j}\left(a_{j}u_{i}\right)+\partial _{i}p-\partial _{j}{\hat {\tau }}_{ji}+S_{u,ij}u_{j}+b_{i}=0,$	(12)
${\frac {1}{h_{0}g}}\delta _{t}^{n}p+\partial _{i}u_{i}=0,$	(13)
$\delta _{t}^{n}\varphi _{a,\alpha }+\partial _{j}\left(a_{j}\varphi _{a,\alpha }\right)-\partial _{i}\left(k_{ij}^{\alpha }\partial _{j}\varphi _{a,\alpha }-k_{ij}^{\alpha }\varphi _{0}^{\alpha }\partial _{j}h_{0}\right)+S_{\alpha \beta }h^{-1}\varphi _{a,\beta }=0.$	(14)

Usando un superíndice adicional para el contador de iteraciones, las variables introducidas en (12), (13) y (14) para la iteración ${\textstyle i}$ -ésima en el paso de tiempo ${\textstyle n}$ son

{\hat {\tau }}_{ij}={\frac {\mu _{H}}{\rho }}\left(\partial _{i}u_{j}+\partial _{j}u_{i}-{\frac {2}{3}}\delta _{ij}\partial _{k}u_{k}\right),\quad \tau _{ij}^{*}={\frac {\mu _{H}}{\rho }}\left(a_{j}\partial _{i}h_{0}+a_{i}\partial _{j}h_{0}-{\frac {2}{3}}\delta _{ij}a_{k}\partial _{k}h_{0}\right),

El super-índice ${\textstyle i}$ es el contador de las iteraciones debidas a la no linealidad. Obsérvese que la forma linealizada de las tensiones viscosas ${\textstyle {\overline {\tau }}_{ij}}$ dada en la ecuación (4) corresponde a ${\textstyle {\hat {\tau }}_{ij}-\tau _{ij}^{*},}$ donde ${\textstyle {\hat {\tau }}_{ij}}$ está escrito en términos de las variables de la iteración actual, mientras que ${\textstyle \tau _{ij}^{*}}$ se calcula con los valores de las variables correspondientes a la iteración anterior. Nótese que hemos usado el símbolo ${\textstyle S_{u}}$ para indicar las componentes de la matriz de términos reactivos provenientes de la ecuación de conservación de la cantidad de movimiento, para distinguirlas de las provenientes de la ecuación de conservación de contaminantes, aunque ambas contribuirán a la matriz de términos reactivos en la forma vectorial que presentamos a continuación.

2.4 Forma vectorial

Para escribir las ecuaciones (12), (13) y (14) en la forma vectorial (7), aunque discretizada en el tiempo, escribimos las mismas como un sistema de ${\textstyle (3+N)}$ ecuaciones con ${\textstyle (3+N)}$ incógnitas. Las ecuaciones correspondientes a las velocidades y altura de la lámina de agua son las mismas que en [27], mientras que la ecuación correspondiente a la concentración de los contaminantes es:

\delta _{t}^{n}\varphi _{a,\alpha }-k_{11}^{\alpha }\partial _{1}^{2}\varphi _{a,\alpha }-k_{12}^{\alpha }\partial _{1}\partial _{2}\varphi _{a,\alpha }-k_{21}^{\alpha }\partial _{2}\partial _{1}\varphi _{a,\alpha }-k_{22}^{\alpha }\partial ^{2}\varphi _{a,\alpha }+a_{1}\partial _{1}\varphi _{a,\alpha }+a_{2}\partial _{2}\varphi _{a,\alpha }+\left(\partial _{1}a_{1}+\partial _{2}a_{2}\right)\varphi _{a,\alpha }

En esta ecuación hay suma para el índice ${\textstyle \beta }$ , pero no para ${\textstyle \alpha }$ , que es el contador de contaminante considerado.

Para simplificar la notación, consideremos un único contaminante, caracterizado por el índice genérico ${\textstyle \alpha }$ . La forma vectorial (7) del sistema de ecuaciones final, previamente discretizadas en el tiempo y linealizadas con el método de punto fijo de Picard, está dada por:

\left[{\begin{array}{c}\delta _{t}^{n}u_{1}\\\delta _{t}^{n}u_{2}\\{\begin{array}{c}{\frac {1}{h_{0}g}}\delta _{t}^{n}p\\\delta _{t}^{n}\varphi _{a,\alpha }\end{array}}\end{array}}\right]-{\frac {\mu _{H}}{\rho }}\left[{\begin{array}{cccc}{\frac {4}{3}}&0&0&0\\0&1&0&0\\0&0&0&0\\0&0&0&k_{11}^{*}\end{array}}\right]\partial _{1}^{2}\left[{\begin{array}{c}u_{1}\\u_{2}\\{\begin{array}{c}p\\\varphi _{a,\alpha }\end{array}}\end{array}}\right]-{\frac {\mu _{H}}{\rho }}\left[{\begin{array}{cccc}0&{\frac {1}{6}}&0&0\\{\frac {1}{6}}&0&0&0\\0&0&0&0\\0&0&0&k_{12}^{*}\end{array}}\right]\partial _{1}\partial _{2}\left[{\begin{array}{c}u_{1}\\u_{2}\\{\begin{array}{c}p\\\varphi _{a,\alpha }\end{array}}\end{array}}\right]

\qquad -{\frac {\mu _{H}}{\rho }}\left[{\begin{array}{cccc}1&0&0&0\\0&{\frac {4}{3}}&0&0\\0&0&0&0\\0&0&0&k_{22}^{*}\end{array}}\right]\partial _{2}^{2}\left[{\begin{array}{c}u_{1}\\u_{2}\\{\begin{array}{c}p\\\varphi _{a,\alpha }\end{array}}\end{array}}\right]-{\frac {\mu _{H}}{\rho }}\left[{\begin{array}{cccc}0&{\frac {1}{6}}&0&0\\{\frac {1}{6}}&0&0&0\\0&0&0&0\\0&0&0&k_{21}^{*}\end{array}}\right]\partial _{2}\partial _{1}\left[{\begin{array}{c}u_{1}\\u_{2}\\{\begin{array}{c}p\\\varphi _{a,\alpha }\end{array}}\end{array}}\right]

\qquad +\left[{\begin{array}{cccc}a_{1}&0&1&0\\0&a_{1}&0&0\\1&0&0&0\\0&0&0&a_{1}\end{array}}\right]\partial _{1}\left[{\begin{array}{c}u_{1}\\u_{2}\\{\begin{array}{c}p\\\varphi _{a,\alpha }\end{array}}\end{array}}\,\right]+\left[{\begin{array}{cccc}a_{2}&0&0&0\\0&a_{2}&1&0\\0&1&0&0\\0&0&0&a_{2}\end{array}}\right]\partial _{2}\left[{\begin{array}{c}u_{1}\\u_{2}\\{\begin{array}{c}p\\\varphi _{a,\alpha }\end{array}}\end{array}}\right]

\qquad +\left[{\begin{array}{cccc}\partial _{1}a_{1}+\partial _{2}a_{2}+S_{u,11}&S_{u,12}&0&0\\S_{u,21}&\partial _{1}a_{1}+\partial _{2}a_{2}+S_{u,22}&0&0\\0&0&0&0\\0&0&0&\partial _{1}a_{1}+\partial _{2}a_{2}+S_{\alpha \alpha }h_{0}^{-1}\end{array}}\right]\left[{\begin{array}{c}u_{1}\\u_{2}\\{\begin{array}{c}p\\\varphi _{a,\alpha }\end{array}}\end{array}}\right]=\left[{\begin{array}{c}-Q_{1}\\-Q_{2}\\{\begin{array}{c}0\\-Q_{4}\end{array}}\end{array}}\right],

(15)

con

En caso de considerar ${\textstyle N}$ contaminantes, el coeficiente 4-4 de la matriz de términos reactivos sería de hecho una matriz ${\textstyle N\times N}$ , con el coeficiente ${\textstyle \alpha \beta }$ dado por ${\textstyle (\partial _{1}a_{1}+\partial _{2}a_{2})\delta _{\alpha \beta }+S_{\alpha \beta }h_{0}^{-1}}$ .

La ecuación (15) corresponde a la ecuación vectorial de CDR transitoria (7) de las ecuaciones acopladas del movimiento de un fluido en aguas poco profundas y del transporte de contaminantes, discretizadas en el tiempo y linealizadas. Las matrices ${\textstyle \mathbf {K} _{ij}}$ , ${\textstyle \mathbf {A} _{i},}$ ${\textstyle \mathbf {S} }$ y ${\textstyle \mathbf {F} }$ , de nuevo para el caso de un contaminante, son

\mathbf {K} _{11}={\frac {\mu _{H}}{\rho }}\left[{\begin{array}{cccc}{\frac {4}{3}}&0&0&0\\0&1&0&0\\0&0&0&0\\0&0&0&k_{11}^{*}\end{array}}\right],\quad \mathbf {K} _{12}={\frac {\mu _{H}}{\rho }}\left[{\begin{array}{cccc}0&{\frac {1}{6}}&0&0\\{\frac {1}{6}}&0&0&0\\0&0&0&0\\0&0&0&k_{12}^{*}\end{array}}\right],\quad \mathbf {K} _{12}={\frac {\mu _{H}}{\rho }}\left[{\begin{array}{cccc}0&{\frac {1}{6}}&0&0\\{\frac {1}{6}}&0&0&0\\0&0&0&0\\0&0&0&k_{21}^{*}\end{array}}\right],

\mathbf {K} _{22}={\frac {\mu _{H}}{\rho }}\left[{\begin{array}{cccc}1&0&0&0\\0&{\frac {4}{3}}&0&0\\0&0&0&0\\0&0&0&k_{22}^{*}\end{array}}\right],\quad \mathbf {A} _{1}=\left[{\begin{array}{cccc}a_{1}&0&1&0\\0&a_{1}&0&0\\1&0&0&0\\0&0&0&a_{1}\end{array}}\right],\quad \mathbf {A} _{2}=\left[{\begin{array}{cccc}a_{2}&0&0&0\\0&a_{2}&1&0\\0&1&0&0\\0&0&0&a_{2}\end{array}}\right],

\mathbf {S} =\left[{\begin{array}{cccc}\partial _{1}a_{1}+\partial _{2}a_{2}+S_{u,11}&S_{u,12}&0&0\\S_{u,21}&\partial _{1}a_{1}+\partial _{2}a_{2}+S_{u,22}&0&0\\0&0&0&0\\0&0&0&\partial _{1}a_{1}+\partial _{2}a_{2}+S_{\alpha \alpha }h_{0}^{-1}\end{array}}\right],\quad \mathbf {F=} \left[{\begin{array}{c}-Q_{1}\\-Q_{2}\\{\begin{array}{c}0\\-Q_{4}\end{array}}\end{array}}\right].

Obsérvese que las matrices ${\textstyle \mathbf {K} _{ij}}$ son simétricas. Los coeficientes de las matrices ${\textstyle \mathbf {A} _{i},}$ ${\textstyle \mathbf {S} }$ y ${\textstyle \mathbf {F} }$ no todos son constantes, pues algunos dependen de las velocidades del flujo, elevación del nivel del agua, su batimetría y sus gradientes, cuyos valores corresponden a la iteración anterior resultado de la linealización. Asimismo, en el caso de reacciones no lineales entre contaminantes los coeficientes ${\textstyle S_{\alpha \beta }}$ pueden depender de las concentraciones de estos evaluadas en la iteración anterior del actual paso de tiempo.

En el caso de ${\textstyle N}$ contaminantes, el vector de incógnitas es

3. Aproximación espacial y estabilización numérica

Una vez hemos conseguido escribir las ecuaciones del problema en la forma (7), el proceso es exactamente el mismo que el seguido en [27] para llevar a cabo la aproximación espacial, por lo que aquí nos centraremos solamente en el problema discreto final y la forma de tratar la ecuación de transporte de contaminantes.

La aproximación de Galerkin del problema es estándar. Dicha aproximación parte de la forma variacional del problema. Una vez discretizado el dominio de cálculo en una partición de elementos finitos ${\textstyle \{\Omega ^{e}\}}$ de ${\textstyle \Omega }$ , con ${\textstyle e=1,\dots ,n_{el}}$ , siendo ${\textstyle n_{el}}$ el número de elementos de la partición, podemos construir el espacio de elementos finitos que aproxima el espacio donde se encuentra la solución del problema continuo. El método de Galerkin consiste en tomar las funciones de test de la forma variacional en el mismo espacio.

En el problema que nos ocupa, el método de Gakerkin puede ser inestable por tres motivos. En primer lugar, la interpolación de la altura de la lámina de agua no puede ser independiente de la de la velocidad promediada: ambas interpolaciones tienen que cumplir la condición inf-sup que garantiza que el problema está bien puesto desde el punto de vista matemático. En segundo lugar, si los términos viscosos (o difusivos) son pequeños comparados con los convectivos, pueden aparecer oscilaciones en todo el dominio de cálculo, como es bien conocido. Este problema afecta tanto a las ecuaciones mecánicas, para altura de agua y velocidad promedio, como a las ecuaciones de transporte de contaminantes. Finalmente, si los términos reactivos son dominantes respecto de los difusivos, pueden aparecer también oscilaciones numéricas, esta vez localizadas cerca de los contornos y de las capas límite. Aunque este problema pudiera parecer poco importante, estas oscilaciones localizadas pueden ser muy perjudiciales en problemas no lineales como el que nos ocupa, pues pueden impedir que los esquemas iterativos converjan.

Todas las inestabilidades descritas se pueden evitar usando métodos de elementos finitos estabilizados. Como se ha dicho anteriormente, en este trabajo usaremos métodos basados en el concepto de VMS. En [17] se presenta una revisión de los mismos aplicados a problemas de mecánica de fluidos, mientras que en [27] se describe su aplicación a un problema formalmente análogo al que nos ocupa, aunque sin las ecuaciones de transporte de contaminantes.

En el problema que estamos considerando, la forma final discreta del problema a resolver es:

\left(\delta _{t}^{n}\mathbf {u} _{h},\mathbf {v} _{h}\right)+B\left(\mathbf {u} _{h},\mathbf {v} _{h}\right)+{\underset {e=1}{\overset {n_{el}}{\sum }}}\int _{\Omega ^{e}}\left(-{\mathcal {L}}^{*}\left(\mathbf {v} _{h}\right)\right)^{T}{\boldsymbol {\tau }}^{e}\mathbf {P} \left[\partial _{t}\mathbf {u} _{h}+{\mathcal {L}}\left(\mathbf {u} _{h}\right)-\mathbf {F} \right]{\,{d}}\Omega =\left\langle \mathbf {F} ,\mathbf {v} _{h}\right\rangle {.}

(16)

Para explicar los distintos términos de este problema discreto, empecemos diciendo que se entiende que esta ecuación corresponde a un cierto paso de tiempo de la discretización temporal, aunque hemos omitido el superíndice del paso de tiempo para no cargar la notación. Asimismo, hemos visto que estas ecuaciones corresponden también a una cierta iteración del proceso iterativo resultado de la linealización del problema.

El subíndice ${\textstyle h}$ hace referencia a funciones en el espacio de elementos finitos (que aquí nada tiene que ver con la altura de la lámina de agua)). Así, ${\textstyle \mathbf {u} _{h}}$ es la aproximación a la incógnita ${\textstyle \mathbf {u} }$ , mientras que ${\textstyle \mathbf {v} _{h}}$ es una función de test arbitraria, cumpliendo las adecuadas condiciones de contorno. Hemos usado el símbolo ${\textstyle \langle \cdot ,\cdot \rangle }$ para indicar la integral sobre el dominio de cálculo del producto de dos funciones, reemplazándolo por ${\textstyle (\cdot ,\cdot )}$ cuando estas funciones son de cuadrado integrable, como sucede con el término temporal.

El operador ${\textstyle {\mathcal {L}}}$ es el asociado al problema (7), es decir

mientras que ${\textstyle {\mathcal {L}}^{*}}$ denota su adjunto (sin condiciones de contorno), dado por

La forma ${\textstyle B}$ en la ecuación (16) es la asociada al problema escrito en forma débil, es decir, la del método de Galerkin. Suponiendo que las matrices ${\textstyle \mathbf {A} _{i}}$ se descomponen en una parte que no se integra por partes, ${\textstyle \mathbf {A} _{i}^{c}}$ , y otra que sí, ${\textstyle \mathbf {A} _{i}^{f}}$ (y que por consiguiente contribuye a los flujos de contorno), y que las condiciones de contornos son homogéneas, para simplificar, en nuestro caso tenemos que

Solo nos queda por describir el término llamado de estabilización en la ecuación (16). En primer lugar, el operador ${\textstyle \mathbf {P} }$ define el método que estamos considerando. En el caso en que ${\textstyle \mathbf {P} =\mathbf {I} }$ , la identidad, nos referiremos a este método como ASGS, mientras que si ${\textstyle \mathbf {P} }$ es la proyección ortogonal al espacio de elementos finitos, el método resultante es el llamado OSS. Los detalles de cómo motivar dichos métodos y cómo se particularizan al problema que nos ocupa (aunque sin las concentraciones de contaminantes), pueden encontrarse en [17,27]. Solo falta describir el cálculo de la llamada matriz de parámetros de estabilización ${\textstyle {\boldsymbol {\tau }}^{e}}$ , calculada en cada elemento ${\textstyle e}$ de la partición de elementos finitos. Tomaremos dicha matriz como diagonal, de la forma

(17)

con

$\tau _{1}^{e}=\left[c_{1}{\frac {\nu _{H}}{\left({\frac {h^{e}}{d^{2}}}\right)^{2}}}+c_{2}{\frac {\|{\boldsymbol {a}}\|}{\frac {h^{e}}{d}}}+c_{3}\|S_{u,11}\|+c_{4}\|S_{u,12}\|\right]^{-1},{\mbox{ }}\tau _{2}^{e}={\frac {\left({\frac {h^{e}}{d}}\right)^{2}}{c_{1}\tau _{1}^{e}}},$	(18)
$\tau _{2+\alpha }^{e}=\left[c_{1}{\frac {\min\{k_{11}^{\alpha },k_{22}^{\alpha }\}}{\left({\frac {h^{e}}{d^{2}}}\right)^{2}}}+c_{2}{\frac {\|{\boldsymbol {a}}\|}{\frac {h^{e}}{d}}}+c_{3}\|S_{\alpha \alpha }\|h_{0}^{-1}\right]^{-1},$	(19)

donde ${\textstyle c_{1},}$ ${\textstyle c_{2},}$ ${\textstyle c_{3}}$ y ${\textstyle c_{4}}$ son parámetros numéricos, y que para los ejemplos numéricos hemos adoptado ${\textstyle c_{1}=12,}$ ${\textstyle c_{2}=2}$ y ${\textstyle c_{3}=c_{4}=1}$ . En estas expresiones, ${\textstyle \nu _{H}={\frac {\mu _{H}}{\rho }}}$ es el coeficiente de viscosidad cinemática, ${\textstyle |{\boldsymbol {a}}|}$ es la norma de la velocidad y ${\textstyle d}$ es el orden polinomial de interpolación.

La justificación de los parámetros de estabilización usados, así como algunas técnicas numéricas que hemos empleado en elementos de alto orden, puede consultarse en [17,30] (véase también [31,32]).

4. Resultados de los experimentos numéricos

Con el propósito de evaluar las formulaciones de los métodos variacionales de subescalas ASGS y OSS con elementos finitos de alto orden en la ecuación vectorial de CDR transitoria, con no linealidad en los términos convectivo y de reacción, hemos realizado diferentes experimentos numéricos, tales como pruebas de convergencia en malla para el problema del transporte de un contaminante, la aproximación de un ejemplo numérico del transporte de un contaminante en una cavidad cuadrada que se encuentra en la literatura, un ejemplo práctico de distribución del transporte de un contaminante en el golfo de Roses y en la desembocadura del río Guadalquivir, y un ejemplo del modelo depredador-presa.

El método de linealización utilizado es el de Picard. El criterio de convergencia para la no linealidad se ha fijado en la norma ${\textstyle L^{2},}$ es decir, el proceso iterativo se considera convergido cuando

donde ${\textstyle i}$ es el contador de iteraciones en un paso de tiempo dado y ${\textstyle tol}$ es la tolerancia de convergencia, fijada en ${\textstyle 10^{-5}}$ en todos los ejemplos. Aquí y en lo que sigue, las normas ${\textstyle L^{2}}$ se extienden a todo el dominio ${\textstyle \Omega {.}}$ En los casos del transporte de contaminantes estudiaremos el transporte de un solo contaminante, excepto en el ejemplo del modelo depredador-presa.

4.1 Pruebas de convergencia en malla

Las pruebas de convergencia en malla para el problema del transporte de un contaminante que hemos realizado consisten en calcular el error en norma ${\textstyle L^{2}}$ entre el valor exacto de una solución manufacturada y la aproximación numérica de una de las incógnitas del problema, el cual debe ser menor o igual que el valor teórico. Sabemos que si el error del método se comporta como el error de la interpolación, situación que puede considerarse óptima, este error debe cumplir con

(20)

donde ${\textstyle u_{h}}$ es la aproximación de la incógnita ${\textstyle u}$ en el espacio de elementos finitos, ${\textstyle C}$ es una constante positiva, ${\textstyle h}$ es el diámetro de la partición de elementos finitos (no confundir con la altura de la lámina de agua) y ${\textstyle d}$ el grado polinómico de la función de forma correspondiente a la interpolación de la aproximación ${\textstyle u_{h}}$ . Consideraremos pruebas de convergencia en malla para ${\textstyle 1\leq d\leq {4}}$ .

Si graficamos la ecuación (20) (en el caso de la igualdad) en el plano, con ${\textstyle \mathrm {log} \,e}$ en las ordenadas y ${\textstyle \mathrm {log} \,h}$ en las abscisas, tenemos una línea recta, donde ${\textstyle d+1}$ es la pendiente teórica óptima de convergencia. Las pruebas de convergencia en malla consisten entonces en verificar que la pendiente calculada sea precisamente ${\textstyle d+1}$ . Para ello, seleccionamos como solución exacta tanto para las velocidades ${\textstyle U_{1}}$ , ${\textstyle U_{2}}$ como para la elevación de la superficie libre del agua ${\textstyle \eta }$ y para la concentración de un contaminante ${\textstyle \varphi }$ , la función polinómica espacio-temporal

e imponemos el término fuente ${\textstyle \mathbf {F} }$ en la ecuación vectorial de CDR transitoria (7) mediante las siguientes expresiones tomadas de las ecuaciones (1), (2) y (3):

siendo ${\textstyle \mathbf {F} =[f_{1},f_{2},f_{3},f_{4}]^{T}}$ y ${\textstyle {\bar {\tau }}_{ji}}$ y ${\textstyle Q_{i}}$ dados por las ecuaciones (4) y (5), respectivamente. Hemos tomado batimetría constante con profundidad de 1 m, es decir, ${\textstyle H=1}$ m, aceleración de la gravedad ${\textstyle g=10}$ ${\textstyle \mathrm {m/s^{2}} }$ , viscosidad cinemática ${\textstyle \nu _{H}=10^{-3}}$ ${\textstyle \mathrm {m^{2}/s} }$ y los coeficientes empíricos de dispersión ${\textstyle k_{11}=k_{22}=10^{-3}}$ ${\textstyle \mathrm {m^{2}/s} ,}$ ${\textstyle k_{12}=k_{21}=0.}$ El parámetro de Coriolis ${\textstyle {\hat {f}},}$ las tensiones de superficie libre debidas al viento ${\textstyle \tau _{3i}^{s},}$ las variaciones de presión atmosférica ${\textstyle p_{a}}$ y el coeficiente de fricción en el fondo Ch no han sido tomados en cuenta para los ejemplos de las pruebas de convergencia en malla, como tampoco el término reactivo de la ecuación de transporte de contaminantes.

A continuación detallamos los resultados de los experimentos numéricos de las pruebas de convergencia en malla con solución analítica conocida, tanto para el método ASGS como para el OSS. En cada gráfica se puede apreciar con línea continua la pendiente teórica de convergencia y con línea con apéndices sobre ella la pendiente calculada. Presentamos los resultados para los distintos grados polinómicos de las funciones de forma estudiados, es decir, elementos lineales, cuadráticos, cúbicos y de cuarto orden triangulares ${\textstyle P_{1},\;P_{2},\;P_{3},\;P_{4}}$ y elementos cuadrangulares ${\textstyle Q_{1},\;Q_{2},\;Q_{3},\;Q_{4}}$ , respectivamente.

El dominio computacional es el cuadrado ${\textstyle \left[0,1\right]\times \left[0,1\right]}$ y el intervalo de tiempo es ${\textstyle \left[0,1\right]}$ , con incrementos para cada paso de tiempo ${\textstyle \delta t=0.2}$ . La malla de elementos finitos consiste en triángulos o cuadrados formando una malla regular. Para cada grado polinómico de las funciones de forma ${\textstyle P_{1},\;P_{2},\;P_{3},\;P_{4}}$ o ${\textstyle Q_{1},\;Q_{2},\;Q_{3},\;Q_{4},}$ hemos calculado los errores con norma ${\textstyle L^{2}}$ para ocho mallas de tamaño ${\textstyle h}$ , siendo ${\textstyle {\frac {1}{h}}=15,\;20,\;25,\;30,\;35,\;40,\;45,\;50,}$ el número de partes en que se ha dividido cada lado del cuadrado. En la Tabla 1 se muestran el número de elementos y el número de nodos para cada refinamiento, tanto para elementos triangulares como para elementos cuadrangulares. Los valores de las constantes algorítmcas ${\textstyle c_{i}}$ , ${\textstyle i=1,2,3,4,}$ se han calibrado a ${\textstyle c_{1}=12,}$ ${\textstyle c_{2}=2,}$ ${\textstyle c_{3}=1,}$ ${\textstyle c_{4}=1}$ y el integrador temporal utilizado es BDF3 (diferencias finitas hacia atrás de cuatro niveles de tiempo y tercer orden en ${\textstyle \delta t}$ ).

Tabla 1. Refinamiento para elementos triangulares y cuadrangulares

Triángulos		Número de nodos
${\textstyle 1/h}$	elem.	$P1$	$P2$	$P3$	$P4$
15	450	256	961	2116	3721
20	800	441	1681	3721	6561
25	1250	676	2601	5776	10201
30	1800	961	3721	8281	14641
35	2450	1296	5041	11236	19881
40	3200	1681	6561	14641	25921
45	4050	2116	8281	18496	32761
50	5000	2601	10201	22801	40401

Cuadrados		Número de nodos
${\textstyle 1/h}$	elem.	$Q1$	$Q2$	$Q3$	$Q4$
15	225	256	961	2116	3721
20	400	441	1681	3721	6561
25	625	676	2601	5776	10201
30	900	961	3721	8281	14641
35	1250	1296	5041	11236	19881
40	1600	1681	6561	14641	25921
45	2025	2116	8281	18496	32761
50	2500	2601	10201	22801	40401

En el encabezado de las gráficas de las pruebas de convergencia en malla de las Figuras 1 y 2 se adjuntan dos filas con los valores de las pendientes calculadas; la primera fila corresponde a los valores de las pendientes de las rectas que pasan por los primeros 5 puntos y la segunda fila son las pendientes de las rectas que pasan por los últimos 5 puntos de los 8 puntos del refinamiento de malla evaluados ${\textstyle \left({\frac {1}{h}}=15,\;20,\;25,\;30,\;35,\;40,\;45,\;50\right)}$ . De esta manera podemos visualizar numéricamente la tendencia de la convergencia de las pendientes calculadas hacia el valor teórico. Para el cálculo de dichas pendientes se ha utilizado el método de los mínimos cuadrados; la simbología ${\textstyle m(P1),}$ ${\textstyle m(P2),}$ ${\textstyle m(P3),}$ ${\textstyle m(P4),}$ corresponde a las pendientes para elementos triangulares: lineales, cuadráticos, cúbicos y de cuarto orden respectivamente, mientras que para elementos cuadrangulares lineales, cuadráticos, cúbicos y de cuarto orden la simbología es ${\textstyle m(Q1),}$ ${\textstyle m(Q2),}$ ${\textstyle m(Q3),}$ ${\textstyle m(Q4)}$ . Los valores de las pendientes teóricas para las variables de la velocidad y de concentración de un contaminante son ${\textstyle m(P1)=m(Q1)=2}$ , ${\textstyle m(P2)=m(Q2)=3}$ , ${\textstyle m(P3)=m(Q3)=4}$ , ${\textstyle m(P4)=m(Q4)=5}$ , y las pendientes teóricas para la variable ${\textstyle \eta }$ , que es la elevación de la superficie libre del agua, son ${\textstyle m(P1)=m(Q1)=1}$ , ${\textstyle m(P2)=m(Q2)=2}$ , ${\textstyle m(P3)=m(Q3)=3}$ , ${\textstyle m(P4)=m(Q4)=4}$ .

En la Figura 1, se muestra la convergencia en malla para la primera componente de la velocidad usando elementos triangulares y en la Figura 2 la convergencia en malla para la concentración de un contaminante. Las otras variables muestran igualmente convergencia óptima para todos los elementos considerados. Lo mismo sucede con los elementos cuadrangulares.


Figura 1. Convergencia en malla ASGS-OSS, elementos triangulares, velocidad $U_{1}$


Figura 2. Convergencia en malla ASGS-OSS, elementos triangulares, concentración $\varphi$

4.2 Transporte de un contaminante en una cavidad cuadrada

Este caso de referencia ha sido ampliamente examinado en la literatura con diferentes métodos, en [6] con el método de diferencias finitas, con el método de volúmenes finitos en [3,5], con el método SUPG en [4]; en todos los casos se ha considerado convección pura, es decir, el coeficiente empírico de dispersión del transporte de un contaminante se ha tomado nulo. En el presente trabajo vamos a reexaminar el caso de referencia con el método variacional ASGS y utilizando funciones de forma cuadráticas, cúbicas y de cuarto orden, en convección dominante, es decir, considerando el coeficiente empírico de dispersión ${\textstyle k_{11}=k_{22}=10^{-3}}$ ${\textstyle \mathrm {m^{2}/s} }$ y ${\textstyle k_{12}=k_{21}=0}$ ${\textstyle \mathrm {m^{2}/s} }$ .

El dominio computacional es un canal cuadrado de 9 km ${\textstyle \times 9}$ km. En cuanto al mallado, con el objeto de comparar los resultados entre elementos lineales y de alto orden, hemos considerado conveniente utilizar refinamientos de malla para elementos ${\textstyle Q_{1},}$ ${\textstyle Q_{2},}$ ${\textstyle Q_{3}}$ y ${\textstyle Q_{4}}$ que den el mismo número de nodos entre si. Tomamos como referencia el refinamiento de malla para ${\textstyle Q_{2}}$ , que consta de ${\textstyle 90\times 90}$ elementos cuadrados uniformes de 9 nodos, es decir, el espaciado en las direcciones ${\textstyle x}$ e ${\textstyle y}$ es ${\textstyle \delta x=\delta y=100}$ ${\textstyle \mathrm {m} }$ , con un total de 8100 elementos y 32761 nodos. De esta manera, la malla de referencia para elementos ${\textstyle Q_{1}}$ consta de 180 ${\textstyle \times }$ 180 elementos cuadrados uniformes de 4 nodos, es decir ${\textstyle \delta x=\delta y=50}$ ${\textstyle \mathrm {m} }$ , con un total de 32400 elementos y 32761 nodos. La malla de referencia para elementos ${\textstyle Q_{3}}$ consta de ${\textstyle 60\times 60}$ elementos cuadrados uniformes de 16 nodos, es decir ${\textstyle \delta x=\delta y=150}$ ${\textstyle \mathrm {m} }$ , con un total de 3600 elementos y 32761 nodos. Finalmente, la malla de referencia para elementos ${\textstyle Q_{4}}$ consta de 45 ${\textstyle \times }$ 45 elementos cuadrados uniformes de 25 nodos, es decir ${\textstyle \delta x=\delta y=200}$ ${\textstyle \mathrm {m} }$ , con un total de 2025 elementos y 32761 nodos.

El tamaño del paso de tiempo es uniforme, con ${\textstyle \delta t=20}$ ${\textstyle \mathrm {s} }$ , y el integrador temporal es BDF3 en todos los casos.

Un flujo uniforme con ${\textstyle u_{1}=u_{2}=0.5}$ ${\textstyle \mathrm {m/s} }$ y ${\textstyle H=0.2485}$ ${\textstyle \mathrm {m} }$ está impuesto en todo el dominio. La concentración inicial de contaminante está dada por la siguiente suposición de dos distribuciones Gausianas centradas respectivamente en ${\textstyle x_{1}=y_{1}=1400}$ ${\textstyle \mathrm {m} }$ y en ${\textstyle x_{2}=y_{2}=2400}$ ${\textstyle \mathrm {m} }$ , dada por:

\varphi (x,y,t=0)=\varphi _{1}\exp {\left({\frac {\left(x-x_{1}\right)^{2}+\left(y-y_{1}\right)^{2}}{\sigma _{1}^{2}}}\right)}+\varphi _{2}\exp {\left({\frac {\left(x-x_{2}\right)^{2}+\left(y-y_{2}\right)^{2}}{\sigma _{2}^{2}}}\right),}

donde ${\textstyle \varphi _{1}=10}$ , ${\textstyle \varphi _{2}=6.5}$ y ${\textstyle \sigma _{1}=\sigma _{2}=264}$ .

Dada la baja difusión del problema, la solución exacta es aquella concentración de contaminante que se mueve diagonalmente a través del dominio con velocidad constante, preservándose su forma durante todo el tiempo. La Figura 3 ilustra los resultados usando el elemento ${\textstyle Q_{2}}$ comparados con la solución teórica en diferentes tiempos, observándose en todos los casos resultados completamente satisfactorios comparados con la solución exacta.

En la Tabla 2, presentamos los valores máximos y mínimos de concentración de contaminante para las funciones de forma ${\textstyle Q_{1}}$ , ${\textstyle Q_{2}}$ , ${\textstyle Q_{3}}$ y ${\textstyle Q_{4}}$ , en donde observamos que a partir de ${\textstyle Q_{2}}$ el error de aproximación al valor exacto máximo es del orden de las centésimas, por lo cual en la Figura 3 hemos usado los resultados del elemento ${\textstyle Q_{2}}$ .




Figura 3. Contornos 3D (izquierda) y distribución del contaminante (derecha) a lo largo de la diagonal de la cavidad cuadrada, con elementos $Q_{2}$

Tabla 2. Valores máximos y mínimos de concentración de contaminante en la cavidad cuadrada para

s

	Solución exacta	$Q_{1}$	$Q_{2}$	$Q_{3}$	$Q_{4}$
Número de elementos	-	32400	8100	3600	2025
Número de nodos	-	32761	32761	32761	32761
Máximo valor de ${\textstyle \varphi }$	10	9.312	9.976	10.03	10.04
Mínimo valor de ${\textstyle \varphi }$	0	$-0.01032$	$-3.702\times {10}^{-10}$	$-2.986\times {10}^{-5}$	$-3.533\times {10}^{-6}$

4.3 Transporte de un contaminante en el golfo de Roses y en la desembocadura del río Guadalquivir

A continuación presentamos dos ejemplos numéricos del transporte de un contaminante. Ambos ejemplos corresponden a situaciones físicas reales, con geometrías complejas. El objetivo es mostrar que el modelo numérico propuesto puede ser utilizado en este tipo de situaciones, aunque no disponemos de valores experimentales con los cuales comparar los resultados y por consiguiente los mostraremos solamente de manera cualitativa, sin entrar en su discusión. Desde el punto de vista numérico, el valor de estos resultados es comprobar que no presentan ningún tipo de oscilación numérica, pese a utilizar igual interpolación para todas las variables y tratarse de situaciones de flujo con convección dominante.

Consideramos una concentración constante de un solo contaminante ${\textstyle \varphi =1}$ en el flujo de entrada y coeficientes de difusión constantes ${\textstyle k_{11}=k_{22}=10^{-3}}$ ${\textstyle \mathrm {m^{2}/s} }$ y ${\textstyle k_{12}=k_{21}=0}$ ${\textstyle \mathrm {m^{2}/s} }$ . La suposición de que los coeficientes de difusión son constantes es debido a que los efectos de turbulencia no han sido considerados en el modelo. También asumimos que la concentración del contaminante se anula en las fronteras sólidas, por lo que prescribimos a cero dicha concentración tanto en el perfil de las costas como en las riveras del río, en los dos problemas tratados. En el resto de los contornos consideramos la concentración de contaminante libre, con concentración inicial igual a cero. En ambos ejemplos numéricos hemos considerado batimetría constante con una profundidad de ${\textstyle H=20.0}$ m (no hemos considerado pues la disminución de la profundidad hasta ${\textstyle H=0}$ ), la aceleración de la gravedad ${\textstyle g=10}$ m/s ${\textstyle ^{2}}$ y la viscosidad cinemática ${\textstyle \nu _{H}=10^{-6}}$ m ${\textstyle ^{2}}$ /s. El parámetro de Coriolis ${\textstyle {\hat {f}}}$ , las tensiones de superficie libre debidas al viento ${\textstyle \tau _{3i}^{s}}$ , las variaciones de presión atmosférica ${\textstyle p_{a}}$ y el coeficiente de fricción en el fondo ${\textstyle {\mbox{Ch}}}$ no han sido tomados en cuenta en ninguno de los dos ejemplos.

El primer ejemplo es la simulación del transporte de un contaminante de concentración constante ${\textstyle \varphi =1}$ , debido a una corriente marina que viene desde el norte en el golfo de Roses sobre la costa catalana, al sur del Cabo de Creus. La geometría del problema y los resultados numéricos se muestran en las Figuras 4, 5 y 6. El dominio de simulación tiene dimensiones de aproximadamente 80 ${\textstyle \times }$ 40 km ${\textstyle ^{2}}$ . La malla de elementos finitos consiste en 6062 elementos triangulares lineales con 3217 nodos.

Una corriente marina de 0.1 m/s que llega desde el norte ha sido prescrita en el borde superior del dominio. Esto corresponde a una velocidad como flujo de entrada en la parte superior del dominio. En la costa la velocidad se prescribe a cero, mientras que en los contornos inferior y derecho del dominio se deja libre, donde la sobreelevación del agua se prescribe a cero.

Figura 4. Golfo de Roses. Geometría y malla de elementos finitos

$t$ =50000 s
$t$ =100000 s
$t$ =150000 s
$t$ =200000 s
	(a)	(b)

Figura 5. Flujo alrededor del golfo de Roses. (a) Velocidad (máx: 0.29 m/s). (b) Elevación (INCOG3), superficie libre del agua (máx: 2.52 mm, mín: -4.05 mm)

$t$ =50000 s
$t$ =100000 s
$t$ =150000 s
$t$ =200000 s
	(a)	(b)

Figura 6. Flujo alrededor del golfo de Roses. (a) Vectores de velocidad. (b) Concentración del transporte de un contaminante (INGOG 4)

En la segunda simulación numérica consideramos el transporte de un contaminante de concentración constante ${\textstyle \varphi =1}$ en la desembocadura del río Guadalquivir, en la costa del sur de España. La geometría del problema y los resultados numéricos se muestran en las Figuras 7, 8 y 9.

Figura 7. Desembocadura del río Guadalquivir. Geometría y malla de elementos finitos

$t$ =50000 s
$t$ =100000 s
$t$ =150000 s
$t$ =200000 s
	(a)	(b)

Figura 8. Flujo en la desembocadura del Guadalquivir. (a) Velocidad (máx: 0.83 m/s). (b) Elevación (INCOG3), superficie libre del agua (máx: 1.97 cm, mín: -1.42 cm)

$t$ =50000 s
$t$ =100000 s
$t$ =150000 s
$t$ =200000 s
	(a)	(b)

Figura 9. Flujo en la desembocadura del Guadalquivir. (a) Vector de velocidad. (b) Concentración del transporte de un contaminante (INGOG 4)

La malla de elementos finitos consiste en 4955 elementos triangulares lineales con 2645 nodos.

Un flujo de corriente marina de ${\textstyle 0.1}$ m/s está prescrito en el contorno inferior en dirección este-oeste y en la desembocadura del río se prescribe una corriente del río de 0.3 m/s, igualmente en dirección este-oeste. En el resto de la costa y en las riveras del río la velocidad se prescribe a cero, y se deja libre en los contornos izquierdo y superior, donde la sobreelevación del agua está prescrita a cero.

Los dos ejemplos numéricos considerados presentan flujos complejos, con importantes variaciones temporales y espaciales. Esto hace que los parámetros de estabilización sean muy variables elemento a elemento. Este hecho dificulta la convergencia de los esquemas iterativos. Puesto que ${\textstyle \tau _{1}}$ y ${\textstyle \tau _{3}}$ (en el caso de ${\textstyle N=1}$ contaminantes) son escalas de tiempo, a menudo se toman proporcionales a ${\textstyle \delta t}$ , y esto se puede justificar con argumentos diversos [24]. En los problemas en los que la convergencia del esquema iterativo es costosa, hemos tomado una cota máxima para ${\textstyle \tau _{1}}$ y ${\textstyle \tau _{3}}$ proporcional a ${\textstyle \delta t}$ , de manera que los hemos redefinido en cada elemento ${\textstyle e}$ como

(21)

donde ${\textstyle \alpha }$ y ${\textstyle \beta }$ son parámetros que hay que ajustar dependiendo del tamaño de la malla y del grado del polinomio de las funciones de forma. El valor de ${\textstyle \tau _{2}^{e}}$ se sigue calculando mediante la expresión dada en (18).

Mediante experimentación numérica hemos calibrado los parámetros ${\textstyle \alpha }$ y ${\textstyle \beta }$ . Para el primer ejemplo hemos considerado ${\textstyle \alpha =15}$ y ${\textstyle \beta =6}$ con un paso de tiempo uniforme ${\textstyle \delta t=}$ 100 s y para el segundo ejemplo ${\textstyle \alpha =3}$ y ${\textstyle \beta =5}$ con ${\textstyle \delta t=}$ 2500 s. En ambos casos hemos usado el integrador temporal BDF2 y presentamos resultados numéricos en diferentes instantes de tiempos de simulación, ${\textstyle t=50000,\;100000,\;150000}$ y ${\textstyle 200000}$ s.

4.4 Modelo depredador-presa

Vamos a considerar como problema a resolver el modelo depredador-presa mostrado en [33], que es un sistema acoplado de ecuaciones transitorias de convección-difusión-reacción no lineales en el término de reacción y que describen espacial y temporalmente la dinámica poblacional de una presa y un depredador en términos de funciones de densidad continuas. Usaremos este problema como ejemplo de sistema con términos reactivos no lineales, con dos concentraciones que compiten entre ellas (en este caso no necesariamente contaminantes). Nos centraremos en el problema propiamente del modelo depredador-presa, suponiendo que la velocidad es conocida. Si embargo, la extensión al caso en el que la velocidad resulta de resolver las ecuaciones del modelo de aguas someras es inmediata. En este problema además es frecuente considerar la posibilidad de que la velocidad de transporte del depredador sea distinta a la de la presa, por lo que cuando introduzcamos estas velocidades de transporte mantendremos esta posible generalización.

4.4.1 Formulación del modelo depredador-pressa

La forma general del modelo es como sigue:

donde ${\textstyle \varphi _{1}}$ y ${\textstyle \varphi _{2}}$ representan las densidades de población de la presa y del depredador, respectivamente. ${\textstyle k_{1}}$ y ${\textstyle k_{2}}$ son los coeficientes de difusión. ${\textstyle P(\varphi _{1})}$ es la función de crecimiento de la población de la presa. La función ${\textstyle E(\varphi _{1},\varphi _{2})}$ representa el acto de depredación resultante en una disminución de la población de la presa y en un incremento en la población del depredador. ${\textstyle \kappa }$ es la eficiencia de la depredación, la cual determina el incremento de la población del depredador, siendo ${\textstyle 0\leq \kappa \leq {1}}$ . Finalmente, ${\textstyle \mu (\varphi _{2})}$ es la función que determina la cantidad de mortalidad del depredador.

Una de las posibilidades para modelar las funciones ${\textstyle P}$ , ${\textstyle E}$ y ${\textstyle \mu }$ es:

donde la función ${\textstyle E(\varphi _{1},\varphi _{2})}$ está modelada por el llamado modelo de Hollinger tipo II, en el cual ${\textstyle B}$ representa la tasa de depredación y ${\textstyle H}$ es la densidad media de saturación de la presa. Para modelar el crecimiento de la población de la presa usamos el llamado modelo logístico, donde ${\textstyle C}$ representa la tasa de crecimiento de la población de la presa y ${\textstyle K}$ es la capacidad de carga del sistema, y denota la cantidad máxima de población de la presa que es soportada por el dominio; en este trabajo asumimos el valor de la unidad para este parámetro. La función de mortalidad ${\textstyle \mu }$ del depredador está dada por un término lineal, siendo ${\textstyle M}$ su coeficiente.

El modelo anterior considera que no hay transporte convectivo ni del depredador ni de la presa. Este transporte puede deberse a que depredador y presa se encuentren en un medio fluido en movimiento, que se correspondería con el modelo general que hemos tratado anteriormente, pero también puede ser necesario introducir la convección a través de un patrón de migración estacional o migración hacia regiones de disponibilidad de recursos. El término de advección resultante modela el movimiento masivo de poblaciones. Con esto se obtiene el sistema depredador-presa que finalmente consideraremos, dado por el sistema de ecuaciones:

al cual hay que añadir las condiciones de contorno y condiciones iniciales. En estas ecuaciones, ${\textstyle s_{11}=C,}$ ${\textstyle K=1,}$ ${\textstyle s_{12}={\frac {B}{H}},}$ ${\textstyle \alpha _{1}={\frac {1}{H}},}$ ${\textstyle s_{21}={\frac {\kappa B}{H}}}$ , ${\textstyle s_{22}=M.}$ El campo de velocidad para la presa es ${\textstyle {\boldsymbol {a}}_{1}=\left[a_{11},a_{12}\right]^{T}}$ y para el depredador ${\textstyle {\boldsymbol {a}}_{2}=\left[a_{21},a_{22}\right]^{T}}$ .

Escribiendo el sistema depredador-presa en el formato de la ecuación vectorial de CDR transitoria (7) tenemos:

{\begin{array}{c}{\begin{array}{c}\left[{\begin{array}{c}\partial _{t}\varphi _{1}\\\partial _{t}\varphi _{2}\end{array}}\right]-\left[{\begin{array}{cc}k_{1}&0\\0&k_{2}\end{array}}\right]\partial _{1}^{2}\left[{\begin{array}{c}\varphi _{1}\\\varphi _{2}\end{array}}\right]-\left[{\begin{array}{cc}k_{1}&0\\0&k_{2}\end{array}}\right]\partial _{2}^{2}\left[{\begin{array}{c}\varphi _{1}\\\varphi _{2}\end{array}}\right]+\left[{\begin{array}{cc}a_{11}&0\\0&a_{21}\end{array}}\right]\partial _{1}\left[{\begin{array}{c}\varphi _{1}\\\varphi _{2}\end{array}}\right]\end{array}}\\\\+\left[{\begin{array}{cc}a_{12}&0\\0&a_{22}\end{array}}\right]\partial _{2}\left[{\begin{array}{c}\varphi _{1}\\\varphi _{2}\end{array}}\right]+\left[{\begin{array}{cc}-s_{11}\left(1-\varphi _{1}\right)&s\left({\frac {\varphi _{1}}{1+\alpha _{1}\varphi _{1}}}\right)\\-s_{21}\left({\frac {\varphi _{2}}{1+\alpha _{2}\varphi _{1}}}\right)&s_{22}\end{array}}\right]\left[{\begin{array}{c}\varphi _{1}\\\varphi _{2}\end{array}}\right]=\left[{\begin{array}{c}0\\0\end{array}}\right].\end{array}}

(22)

La matriz de los parámetros de estabilización utilizada para este problema, que está descrita en [33], es ${\textstyle {\boldsymbol {\tau }}=\mathrm {diag} (\tau _{1},\tau _{2})}$ , con:

donde recordemos que ${\textstyle h}$ corresponde al diámetro del elemento y ${\textstyle d}$ es el orden polinomial de interpolación. Las variables ${\textstyle \varphi _{1}}$ y ${\textstyle \varphi _{2}}$ las tomamos calculadas en la iteración anterior del paso de tiempo actual.

4.4.2 Pruebas numéricas y resultados

Consideremos la ecuación (22) del modelo depredador-presa dentro de un cuadrado unitario ${\textstyle \Omega =\left[0,1\right]\times \left[0,1\right]}$ , con condiciones de contorno de Dirichlet homogéneas y condiciones iniciales de las densidades de población de la presa y del depredador dadas por las distribuciones normales (23) y (24), respectivamente, escritas a continuación y mostradas en la Figura 10:

$\varphi _{1}(x,y,t=0)=\exp \left[{-50\left(\left(x-0.25\right)^{2}+\left(y-0.25\right)^{2}\right)}\right],$	(23)
$\varphi _{2}(x,y,t=0)=\exp \left[{-50\left(\left(x-0.75\right)^{2}+\left(y-0.75\right)^{2}\right)}\right].$	(24)


Figura 10. Condición inicial de la densidad de población de la presa (izquierda) y del depredador (derecha)

En la Tabla 3 se encuentran los valores de los coeficientes de reacción ${\textstyle s_{11},}$ ${\textstyle s_{12},}$ ${\textstyle s_{21}}$ y ${\textstyle s_{22}}$ para varios casos de prueba que hemos considerado. Para todos estos casos de prueba los coeficientes de difusión se tomaron como ${\textstyle k_{11}=k_{22}=10^{-4}}$ m ${\textstyle ^{2}}$ /s. El campo de velocidades se mantiene constante tanto para el depredador como para la presa. Las componentes de velocidad de la presa son ${\textstyle a_{11}=0.5}$ m/s, ${\textstyle a_{12}=0.5}$ m/s y para el depredador son ${\textstyle a_{21}=-0.5}$ m/s, ${\textstyle a_{22}=-0.5}$ m/s, con lo cual las poblaciones del depredador y de la presa son conducidas en direcciones opuestas para encontrarse frente a frente una con la otra. No consideramos ningún término fuente y las constantes ${\textstyle \alpha _{1}}$ y ${\textstyle \alpha _{2}}$ se tomaron a la unidad. Con esto el sistema de ecuaciones del modelo depredador-presa se escribe:

{\begin{array}{c}{\begin{array}{c}\left[{\begin{array}{c}\partial _{t}\varphi _{1}\\\partial _{t}\varphi _{2}\end{array}}\right]-\left[{\begin{array}{cc}10^{-4}&0\\0&10^{-4}\end{array}}\right]\partial _{1}^{2}\left[{\begin{array}{c}\varphi _{1}\\\varphi _{2}\end{array}}\right]-\left[{\begin{array}{cc}10^{-4}&0\\0&10^{-4}\end{array}}\right]\partial _{2}^{2}\left[{\begin{array}{c}\varphi _{1}\\\varphi _{2}\end{array}}\right]+\left[{\begin{array}{cc}0.5&0\\0&-0.5\end{array}}\right]\partial _{1}\left[{\begin{array}{c}\varphi _{1}\\\varphi _{2}\end{array}}\right]\end{array}}\\\\+\left[{\begin{array}{cc}0.5&0\\0&-0.5\end{array}}\right]\partial _{2}\left[{\begin{array}{c}\varphi _{1}\\\varphi _{2}\end{array}}\right]+\left[{\begin{array}{cc}-s_{11}\left(1-\varphi _{1}\right)&\quad s_{12}\left({\frac {\varphi _{1}}{1+\varphi _{1}}}\right)\\-s_{21}\left({\frac {\varphi _{2}}{1+\varphi _{1}}}\right)&s_{22}\end{array}}\right]\left[{\begin{array}{c}\varphi _{1}\\\varphi _{2}\end{array}}\right]=\left[{\begin{array}{c}0\\0\end{array}}\right].\end{array}}

(25)

Tabla 3. Casos de prueba del modelo depredador-presa para diferentes coeficientes de reacción

	$s_{11}$	$s_{12}$	$s_{21}$	$s_{22}$
Caso 1	0	0	0	0
Caso 2	0	0	0	1
Caso 3	0	0	3	0.1
Caso 4	1	0	3	0.1
Caso 5	1	2	3	0.1

La malla de elementos finitos usada es regular y consta de ${\textstyle 50\times 50}$ elementos cuadrados de 4 nodos, es decir ${\textstyle \delta x=\delta y=0.02}$ , con un total de 2500 elementos y 2601 nodos. El intervalo de tiempo es ${\textstyle \left[0,1\right]}$ y el tamaño del paso de tiempo tomado es uniforme, con ${\textstyle \delta t=0.2}$ . Hemos usado el integrador temporal BDF2 y el método variacional multiescala ASGS descrito anteriormente. Presentamos resultados numéricos en diferentes instantes de tiempos de simulación, ${\textstyle t=0.2,}$ ${\textstyle 0.4}$ , ${\textstyle 0.6,}$ ${\textstyle 0.8}$ , ${\textstyle 1.0}$ , para todos los casos de prueba.

En las Figuras 11, 12, 13, 14 y 15 se encuentran los resultados de los casos 1, 2, 3, 4 y 5, respectivamente, considerados en la Tabla 3. Los resultados numéricos en todos los casos son los esperados y corresponden con la realidad física de cada caso, observando además que son coincidentes con los resultados mostrados en [33].

Las diferencias de los resultados de las densidades de población del depredador y de la presa con elementos de alto orden ${\textstyle Q_{2}}$ (con 10201 nodos), ${\textstyle Q_{3}}$ (con 22801 nodos) y ${\textstyle Q_{4}}$ (con 40401 nodos) están en el orden de las milésimas con respecto a los resultados con elementos lineales para ${\textstyle t=1.0}$ , tal como observamos en las Tablas 4 y 5, por lo que dependiendo de la escala del problema es suficiente la precisión con elementos lineales.

$t$ =0.2 s
$t$ =0.4 s
$t$ =0.6 s
$t$ =0.8 s
$t$ =1.0 s

Figura 11. Caso 1. Densidad de población de la presa (izquierda). Densidad de población del depredador (derecha)

$t$ =0.2 s
$t$ =0.4 s
$t$ =0.6 s
$t$ =0.8 s
$t$ =1.0 s

Figura 12. Caso 2. Densidad de población de la presa (izquierda). Densidad de población del depredador (derecha)

$t$ =0.2 s
$t$ =0.4 s
$t$ =0.6 s
$t$ =0.8 s
$t$ =1.0 s

Figura 13. Caso 3. Densidad de población de la presa (izquierda). Densidad de población del depredador (derecha)

$t$ =0.2 s
$t$ =0.4 s
$t$ =0.6 s
$t$ =0.8 s
$t$ =1.0 s

Figura 14. Caso 4. Densidad de población de la presa (izquierda). Densidad de población del depredador (derecha)

$t$ =0.2 s
$t$ =0.4 s
$t$ =0.6 s
$t$ =0.8 s
$t$ =1.0 s

Figura 15. Caso 5. Densidad de población de la presa (izquierda). Densidad de población del depredador (derecha)

Tabla 4. Densidad de población de la presa (

)

	Caso 1	Caso 2	Caso 3	Caso 4	Caso 5
${\textstyle Q_{1}}$	0.48614	0.48614	0.48614	0.69991	0.56545
${\textstyle Q_{2}}$	0.48847	0.48847	0.48847	0.70560	0.56853
${\textstyle Q_{3}}$	0.48847	0.48847	0.48847	0.70709	0.56928
${\textstyle Q_{4}}$	0.48848	0.48848	0.48848	0.70786	0.56982

Tabla 5. Densidad de población del depredador (

))

	Caso 1	Caso 2	Caso 3	Caso 4	Caso 5
${\textstyle Q_{1}}$	0.48614	0.22615	0.63921	0.70337	0.65541
${\textstyle Q_{2}}$	0.48847	0.22685	0.63872	0.71223	0.66228
${\textstyle Q_{3}}$	0.48847	0.22684	0.63954	0.71348	0.66352
${\textstyle Q_{4}}$	0.48848	0.22684	0.64002	0.71412	0.66422

5. Conclusiones

En este artículo se ha presentado la aproximación del modelo acoplado de las ecuaciones del movimiento de un fluido en aguas poco profundas junto con las ecuaciones de convección-difusión-reacción del transporte de contaminantes mediante formulaciones estabilizadas de elementos finitos de alto orden (hasta el cuarto orden).

Las formulaciones estabilizadas de elementos finitos para la ecuación vectorial de CDR transitoria, ASGS y OSS, nos permiten utilizar igual interpolación en todas las variables, que para el modelo del transporte de contaminantes en aguas poco profundas son las dos componentes de velocidad ${\textstyle U_{1}}$ , ${\textstyle U_{2}}$ , la elevación ${\textstyle \eta }$ de la superficie libre del agua y las concentraciones de contaminantes. Asimismo, las formulaciones ASGS y OSS para la ecuación vectorial de CDR transitoria nos permiten tratar con flujos de convección y reacción dominantes, incluyendo problemas no lineales en los términos convectivos y reactivos. Como ejemplo de problema de reacción no lineal hemos presentado el modelo depredador-presa, con dos concentraciones acopladas a través del término reactivo como incógnitas.

Para utilizar elementos finitos de alto orden, en la matriz de los parámetros de estabilización se ha incluido la dependencia con el orden polinomial de interpolación.

Los resultados de las gráficas de convergencia muestran claramente que las formulaciones para los métodos estabilizados ASGS y OSS son capaces de aproximar correctamente el problema vectorial transitorio de CDR con convección dominante, particularizado al sistema acoplado de cuatro ecuaciones correspondientes al flujo en aguas poco profundas y transporte de un contaminante.

Los resultados calculados del transporte de un contaminante en una cavidad cuadrada no muestran errores significativos de dispersión o disipación, por lo que concuerdan muy bien con los valores teóricos, y son mejores que los resultados presentados en la literatura, como en [3,4,5].

Asimismo, se observa que los resultados obtenidos del transporte de un contaminante en el golfo de Roses y en la desembocadura del río Guadalquivir son completamente satisfactorios, tanto para el campo de velocidades y el nivel de la superficie libre del agua como para la distribución de la concentración de un contaminante.

Finalmente los resultados del modelo depredador-presa son coherentes con la realidad física de cada caso, y muestran la robustez de las formulaciones estabilizadas ASGS y OSS para la ecuación vectorial de CDR transitoria en problemas de reacción dominante y de no linealidad en el término de reacción.

Agradecimientos

R. Codina agradece la ayuda recibida a través del programa ICREA Academia, del Gobierno de Catalunya.

Referencias

[1] Zlatev Z. The air pollution problem. In: Computer Treatment of Large Air Pollution Models, Environmental Science and Technology Library, vol 2, Springer, Dordrecht, 1995.

[2] Zienkiewicz O.C., Taylor R.L. The Finite element method. Fluid Mechanics. Butterworth- Heinemann, vol. 3, 5th edition, 2000.

[3] Li S., Duffy C.J. Fully-coupled modeling of shallow water flow and pollutant transport on unstructured grids. Procedia Environmental Sciences, 13:20982121, 2012.

[4] Behzadi F., Newman III J.C. A Semi-discrete SUPG method for contaminant transport in shallow water models. Procedia Computer Science, 80:1313-1323, 2016.

[5] Benkhaldoun F., Elmahi I., Seaïd M. Well-balanced finite volume schemes for pollutant transport by shallow water equations on unstructured meshes. J. Comput. Phys., 1:180-203, 2007.

[6] Komatsu T., Ohgushi A.K., Asai K. Refined numerical scheme for advective transport in diffusion simulation. J. Hydraul. Eng., 123:41-50, 1997.

[7] Begnudelli L., Sanders B.F. Unstructured grid finite-volume algorithm for shallow-water flow and scalar transport with wetting and drying. Journal of Hydraulic Engineering, 132(4):371-384, 2006.

[8] Cea L., Vázquez-Cedón M.E. Unstructured finite volume discretisation of bed friction and convective flux in solute transport models linked to the sallow water equations. Journal of Computational Physics, 231(8):3317-3330, 2012.

[9] Caleffi V., Valiani A. A 2D local discontinuous Galerkin method for contaminant transport in channel bends. Computers & Fluids, 88:629-642, 2013.

[10] Pavan S., Hervouet J.M., Ricchiuto M., Ata R. A second order residual based predictorcorrector approach for time dependent pollutant transport. Journal of Computational Physics, 318:122141, 2016.

[11] Cai L., Xie W.-X., Feng J.-H., Zhou J. Computations of transport of pollutant in shallow water. Applied Mathematical Modelling, 31:490498, 2007.

[12] Postma L., Hervouet J.-M. Compatibility between finite volumes and finite elements using solutions of shallow water equations for substance transport. Int. J. Numer. Methods Fluids, 53(9):14951507, 2007.

[13] Pavan S., Ata R., Hervouet J.-M. Finite volume schemes and residual distribution schemes for pollutant transport on unstructured grids. Environ. Earth Sci., 74:73377356, 2015.

[14] Liu H., Zhou J.G., Li M., Zhao Y. Multi-block lattice Boltzmann simulations of solute transport in shallow water flows. Advances in Water Resources, 58:24-40, 2013.

[15] Hammou H., Ginzburg I., Boulerhcha M. Two-relaxation-times lattice Boltzmann schemes for solute transport in unsaturated water flow with a focus on stability. Adv Wat Res., 4(6):779-793, 2011.

[16] Zhou J.G. A lattice Boltzmann method for solute transport. Int J Numer Methods Fluids, 61:848-63, 2009.

[17] Codina R., Badia S., Baiges J., Principe J. Variational multiscale methods in computational fluid dynamics. In: Encyclopedia of Computational Mechanics, E. Stein, R. Borst and T.J.R. Hughes (eds.), doi: 10.1002/9781119176817.ecm21172017, 2017.

[18] Rogers B., Fujihara M., Borthwick A.G.L. Adaptive Q-tree Godunov-type scheme for shallow water equations. Int J Numer Methods Fluids, 35:247-80, 2001.

[19] H. Yamamoto. Lecture 8: The shallow-water equations. June 2009, 9pp. https://www.coursehero.com/file/35480149/othersonline-HW4pdf/

[20] Hughes T.J.R. Multiscale phenomena: Green's functions, the Drichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Computer Methods in Applied Mechanics and Engineering, 127(1-4):387-401, 1995.

[21] Hughes T.J.R., Feijóo G.O., Mazzei L., Quincy J.B.. The variational multiscale method-a paradign for computational mechanics. Computer Methods in Applied Mechanics and Egineering, 166:3-24, 1998.

[22] Codina R. On stabilized finite element methods for linear systems of convection-diffusion-reaction equations. Computer Methods in Applied Mechanics and Engineering, 188:61-82, 2000.

[23] Codina R. A stabilized finite element method for generalized stationary incompressible flows. Computer Methods in Applied Mechanics and Enginineering, 190:2681-2706, 2001.

[24] Codina R. Comparison of some finite element methods for solving the diffusion-convection-reaction equation. Computer Methods in Applied Mechanics and Engineering, 156:185-210, 1998.

[25] Codina R., Blasco J. Analysis of a stabilized finite element approximation of the transient convection-diffusion-reaction equation using orthogonal subscales. Computing and Visualization in Science, 4:167-174, 2002.

[26] Codina R. Stabilized finite element approximation of transient incompressible flows using orthogonal subscales. Computer Methods in Applied Mechanics and Engineering, 191:4295-4321, 2002.

[27] Villota A., Codina R. Approximation of the shallow water equations with higher order finite elements and variational multiscale methods. Rev. int. métodos numér. cálc. diseño ing., 34(1), 28, 2018. URL https://www.scipedia.com/public/Villota_Codina_a.

[28] Codina R. Finite element approximation of the convection-diffusion equation: subgrid-scale spaces, local instabilities and anisotropic space-time discretizations. In Bail 2010 - Boundary and Interior Layers, Computational and Asymptotic Methods, Lecture Notes in Computational Science and Engineering, C. Clavero, J.L. Gracia and F.J. Lisbona (eds.), vol. 81, pages 85-97, Springer, 2011.

[29] Codina R., González-Ondina J.M., Díaz-Hernández G., Principe J. Finite element approximation of the modified Boussinesq equations using a stabilized formulation, International Journal for Numerical Methods in Fluids, 57:1249-1268, 2008.

[30] Villota A., Codina R. Approximation of the scalar convection-diffusion-reaction equation with stabilized finite element formulations of high order. Rev. int. métodos numér. cálc. diseño ing., 35(1), 6, 2019. URL https://www.scipedia.com/public/Villota_Codina_a

[31] Principe J., Codina R. On the stabilization parameter in the subgrid scale approximation of scalar convection-diffusion-reaction equations on distorted meshes. Computer Methods in Applied Mechanics and Engineering, 199:1386-1402, 2010.

[32] Codina R. On hp convergence of stabilized finite element approximations of the convection-diffusion equation. SeMA Journal, 75:591-606, 2018.

[33] Sheshachala S.K., Codina R. Finite element modeling of nonlinear reaction-diffusion-advection systems of equations. International Journal of Numerical Methods for Heat & Fluid Flow, 28(11):2688-2715, 2018. https://doi.org/10.1108/HFF-02-2018-0077