Approximation of the shallow water equations with higher order finite elements and variational multiscale methods

==Aproximación de las ecuaciones del movimiento de un fluido en aguas poco profundas con elementos finitos de alto orden mediante métodos variacionales de subescalas==

'''A. Villota<math>\mathrm{^{a,b}}</math> and R. Codina<math>\mathrm{^{a}}</math>'''

{|
|-
|<math>^{\mathrm{a}}</math>Universitat Politecnica de Catalunya, Jordi Girona 1-3, Edifici C1, 08034 Barcelona, Spain
|}

{|
|-
|<math>^{\mathrm{b}}</math>Escuela Politécnica Nacional, Ladrón de Guevara E 11-253, Quito, Ecuador
|}

==Abstract==

En este artículo presentamos aproximaciones de elementos finitos de alto orden usando métodos variacionales estabilizados de subescalas para aproximar las ecuaciones del movimiento de un fluido en aguas poco profundas (''shallow water equations'', en inglés). Escribimos estas ecuaciones como un sistema de ecuaciones de tipo convección-difusión-reacción (CDR) no lineal y transitorio, y planteamos los desarrollos en este marco general. Los métodos variacionales de multiescala (VMS) se basan en la descomposición de las incógnitas del problema continuo en una parte resoluble en el espacio de elementos finitos y otra que no puede ser capturada por la malla de elementos finitos y que la denominamos subescala. La subescala se aproxima en términos de la solución de elementos finitos, obteniéndose un esquema numérico robusto, que permite en particular igual interpolación para todas las incógnitas y la posibilidad de tratar con flujos de convección dominante (no consideraremos la posibilidad de tratar con choques). Los métodos VMS que consideraremos son los llamados de subscalas algebraicas (ASGS, por ''Algebraic SubGrid Scales'') y subescalas ortogonales (OSS, por ''Orthogonal Subgrid Scales'').

'''keywords'''

Aguas poco profundas, elementos finitos, estabilización, sistemas acoplados, alto orden.

==1 Introducción==

Las ecuaciones del movimiento de un fluido en aguas poco profundas se obtienen a partir de las ecuaciones Euler para fluidos invíscidos o a partir de las ecuaciones de Navier-Stokes para fluidos con viscosidad. Estas ecuaciones describen el movimiento en una capa delgada de fluido de densidad constante en equilibrio hidrostático, limitada en su parte inferior por la topografía del terreno y por arriba por una superficie libre, y promediando las variables a lo largo de la profundidad. Son una buena aproximación del comportamiento de un fluido cuando la profundidad del dominio estudiado es pequeña con respecto al tamaño global de este dominio. Es por esto que dichas ecuaciones describen bien el flujo en canales, ríos, lagos y lagunas, flujos de marea, corrientes marinas, avance de un frente de onda, arrastre de sedimentos, variación de concentración salina, flujo en la rotura de presas, flujos atmosféricos, transporte de contaminantes y tsunamis, entre otros; véase por ejemplo <span id='citeF-1'></span><span id='citeF-2'></span><span id='citeF-3'></span><span id='citeF-4'></span>[[#cite-1|[1,2,3,4]]].

Se ha realizado una extensa investigación numérica en el área de las ecuaciones del movimiento de un fluido en aguas poco profundas, aplicando  métodos de diferencias finitas <span id='citeF-5'></span>[[#cite-5|[5]]], métodos de volúmenes finitos <span id='citeF-6'></span><span id='citeF-7'></span><span id='citeF-8'></span>[[#cite-6|[6,7,8]]], métodos de elementos finitos <span id='citeF-9'></span><span id='citeF-10'></span><span id='citeF-11'></span><span id='citeF-12'></span><span id='citeF-13'></span>[[#cite-9|[9,10,11,12,13]]], métodos tipo Godunov <span id='citeF-14'></span><span id='citeF-15'></span>[[#cite-14|[14,15]]] y métodos de Boltzmann&nbsp;<span id='citeF-16'></span>[[#cite-16|[16]]].

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 <span id='citeF-17'></span>[[#cite-17|[17]]], entre otras referencias. El problema que se debe resolver puede escribirse en forma conservativa en términos de la profundidad del fluido y de las velocidades promediadas en la plano donde el flujo medio tiene lugar. La presencia del término fuente en la ecuación de momento, la consideración de batimetría variable, fricción en el fondo del lecho, tensión en la superficie libre del agua debida al viento y el efecto de Coriolis son aspectos que se pueden tener en cuenta en el modelo físico.

En este trabajo vamos ha realizar la aproximación de las ecuaciones del movimiento de un fluido en aguas poco profundas con elementos finitos de alto orden y usando los métodos VMS (''Variational Multi-Scale)'' llamados ASGS  (''Algebraic SubGrid Scale'') y OSS (''Orthogonal Subgrid Scale''). La contribución de este trabajo es tanto aplicar estas técnicas (por otro lado conocidas) al problema de flujo de aguas poco profundas como experimentar con elementos de orden elevado (hasta cuarto orden).

Veremos en primer lugar cómo a través de un cambio de variables adecuado y linearizando las ecuaciones de partida, cada iteración del proceso no lineal puede escribirse como un sistema acoplado de ecuaciones de convección-difusión-reacción (CDR). Plantearemos la aproximación de este sistema general de ecuaciones, particularizándolo después al problema de aguas poco profundas.

La estrategia de discretización que plantearemos consiste en dos pasos. Primero procederemos a discretizar en el tiempo mediante un esquema de integración temporal en diferencias finitas, y luego realizaremos una aproximación de elementos finitos en el espacio. Este procedimiento desacopla los errores provenientes de la discretización temporal de los correspondientes a la discretización espacial. Sin embargo, se tiene que señalar que 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. Para el método convencional de Galerkin ambos caminos dan lugar al mismo problema discreto, pero no es así para todos los métodos de estabilización basados en VMS. En lo concerniente a la discretización temporal, usaremos la regla trapezoidal generalizada, que es el método en diferencias finitas más sencillo; sin embargo, podríamos aplicar cualquier otro esquema en diferencias finitas.

En cuanto a la discretización espacial, es bien conocido que el método estándar de Galerkin puede fallar por dos razones: la presencia del término convectivo (no-lineal) dominante con respecto al término viscoso, y la necesidad de usar diferentes interpolaciones de elementos finitos para la profundidad y las velocidades promediadas, similarmente a lo que ocurre para la presión y la velocidad en el clásico problema de Stokes; dicha necesidad se plantea en forma de condición inf-sup. Ambas pueden superarse mediante una formulación estabilizada. Aquí vamos a usar las formulaciones basadas en subescalas y, en particular, el enfoque introducido por Hughes y otros en <span id='citeF-18'></span><span id='citeF-19'></span>[[#cite-18|[18,19]]] para la ecuación de CDR escalar  y generalizado a sistemas lineales en <span id='citeF-20'></span><span id='citeF-21'></span>[[#cite-20|[20,21]]] y que corresponde al método ASGS para resolver el problema de CDR vectorial estacionario. 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. Una importante característica introducida en <span id='citeF-22'></span><span id='citeF-23'></span>[[#cite-22|[22,23]]]  es que la componente no resuelta, denominada subescala, se considera ortogonal respecto del producto escalar de <math display="inline">L^2</math> al espacio de elementos finitos, lo cual se explicará más adelante. Esta idea fue primero introducida en <span id='citeF-23'></span>[[#cite-23|[23]]] como una extensión de un método de estabilización originalmente introducido para el problema de Stokes en <span id='citeF-24'></span>[[#cite-24|[24]]] y totalmente analizado para la ecuación estacionaria de Navier-Stokes en <span id='citeF-25'></span>[[#cite-25|[25]]]. Con estas ideas aplicadas al problema de CDR vectorial estacionario se obtiene el método OSS para resolver dicho problema. Finalmente, para el problema de CDR vectorial transitorio, que es el que usamos directamente para resolver las ecuaciones del movimiento de un fluido en aguas poco profundas, la formulación obtenida está basada en lo expuesto en <span id='citeF-26'></span><span id='citeF-27'></span><span id='citeF-22'></span>[[#cite-26|[26,27,22]]]. 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 su aproximación hay que recurrir a técnicas de captura de discontinuidades que no consideraremos en este trabajo.

El artículo está organizado como sigue. En la siguiente sección se describe la física del problema que queremos resolver y la obtención de la ecuación vectorial de CDR transitoria que representa el problema del movimiento de un fluido en aguas poco profundas. Las aproximaciones y formulaciones numéricas están descritas en la Sección&nbsp;3. En la Sección&nbsp;4 se presentan pruebas numéricas de convergencia en malla con soluciones analíticas conocidas y ejemplos numéricos de referencia en la literatura. Finalmente, en base a los resultados obtenidos se extraen algunas conclusiones en la Sección&nbsp;5.

==2 Planteamiento del problema, linealización y discretización temporal==

===2.1 Ecuaciones de gobierno===

Sea <math display="inline">\Omega </math> un dominio acotado de <math display="inline">\mathbb{R}^2</math> con coordenadas cartesianas <math display="inline">\mathbf{x} = (x_1,x_2)</math>. Consideremos un fluido que se mueve durante el intervalo de tiempo <math display="inline">[0,T]</math> en el dominio de <math display="inline">\mathbb{R}^3</math> dado por <math display="inline">\Omega \times (-H({\mathbf{x}}), \eta ({\mathbf{x}},t))</math>, donde <math display="inline">x_3 = -H({\mathbf{x}})</math> representa la batimetría del terreno y <math display="inline">x_3 = \eta ({\mathbf{x}},t)</math> es la superficie libre del fluido, <math display="inline">{\mathbf{x}}\in \Omega </math>, <math display="inline">t\in [0,T)</math>, que será incógnita del problema. Consideraremos que <math display="inline">h ({\mathbf{x}},t ) := H({\mathbf{x}}) + \eta ({\mathbf{x}},t)</math> es mucho menor que el diámetro de <math display="inline">\Omega </math>, por lo que supondremos válida la aproximación de las ecuaciones de movimento en aguas poco profundas, las cuales se obtienen por integración vertical de las ecuaciones generales de movimiento, bien sea de fluidos viscosos (ecuaciones de Navier-Stokes) o invíscidos (ecuaciones de Euler); véase por ejemplo <span id='citeF-17'></span><span id='citeF-28'></span>[[#cite-17|[17,28]]]. 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. Dichas ecuaciones pueden escribirse en forma indicial como

<span id="eq-1"></span>
<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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,\quad i = 1,2,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|-
| style="text-align: center;" | <math> \partial _{t}h+\partial _{i}\left(hU_{i}\right) =0,   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}
|}

con

<span id="eq-3"></span>
<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\overline{\tau }_{ij}  =\mu _{H}\left(\partial _{j}U_{i}+\partial _{i}U_{j}-\frac{2}{3}\delta _{ij}\partial _{k}U_{k}\right),\quad i,j = 1,2,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|-
| style="text-align: center;" | <math> 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}{\hbox{Ch}}^{-1}h^{-2},\quad i=1,2, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}
|}

y

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>U_{0}=\sqrt{U_{1}^{2}+U_{2}^{2}}, </math>
|}
|}

donde <math display="inline">U_{i}({\mathbf{x}},t)</math> <math display="inline">\left(i=1,2\right)</math> son las componentes de la velocidad promediadas en la profundidad, <math display="inline">\delta _{ij}</math> es la delta de Kronecker, <math display="inline">\partial _{t}</math> denota la derivada temporal y <math display="inline">\partial _{i}</math> la derivada con respecto a la <math display="inline">i</math>-ésima coordenada cartesiana <math display="inline">(i=1,2)</math> en <math display="inline">\Omega </math>.  En las ecuaciones anteriores y en lo que sigue, índices repetidos indican la notación de sumación de Einstein, y los subíndices <math display="inline">i,</math> <math display="inline">j,</math> <math display="inline">k</math> recorren desde 1 hasta 2. En las expresiones anteriores, <math display="inline">g</math> es la aceleración de la gravedad, <math display="inline">\rho </math> la densidad del fluido, <math display="inline">\bar{\tau }_{ij}</math> son las tensiones viscosas promediadas, <math display="inline">\mu _{H}</math> es el coeficiente de viscosidad dinámico, <math display="inline">\hat{f}</math> es el parámetro de Coriolis, <math display="inline">p_{a}</math> es la presión atmosférica, <math display="inline">\tau _{3i}^{s}</math> <math display="inline">\left(i=1,2\right)</math> son las tensiones en la superficie libre del agua debidas al viento, <math display="inline">{\hbox{Ch}}</math> es el coeficiente de Chézy para el rozamiento en el fondo y <math display="inline">U_{0}</math> es la norma euclídea de <math display="inline">\left(U_{1},U_{2}\right)</math>.

Las ecuaciones ([[#eq-1|1]]) y ([[#eq-2|2]]) representan un sistema de tres ecuaciones, cuyo vector de incógnitas es

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{U}^T =\left[U_{1}, U_{2},  h \right].</math>
|}
|}

Dichas ecuaciones constituyen las ecuaciones del movimiento de un fluido en aguas poco profundas. Se deben completar con condiciones iniciales y de contorno. Los valores iniciales para la elevación del agua y las velocidades promediadas 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 <math display="inline">\Omega </math>, mientras que en el resto del contorno se prescribe la componente normal de las tensiones viscosas y la elevación del fluido.

Las ecuaciones para aguas poco profundas descritas anteriormente, en su forma conservativa, se presentaron en esta forma en las referencias <span id='citeF-29'></span><span id='citeF-30'></span>[[#cite-29|[29,30]]]. Debido a su complejidad, en la literatura existen muchas variantes de estas ecuaciones con distintas simplificaciones, tanto para las ecuaciones que parten de las ecuaciones de Euler <span id='citeF-28'></span>[[#cite-28|[28]]], como para las ecuaciones obtenidas de las ecuaciones de Navier-Stokes&nbsp;<span id='citeF-17'></span>[[#cite-17|[17]]].

En nuestro trabajo utilizaremos las ecuaciones ([[#eq-1|1]]) y ([[#eq-2|2]]) 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 convección-difusión-reacción (CDR) transitoria de la forma

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial _{t}\mathbf{u}-\partial _{i}\left(\mathbf{K}_{ij}\partial _{j}\mathbf{u}\right)+\partial _{i}\left(\mathbf{A}_{i}\mathbf{u}\right)+\mathbf{Su}=\mathbf{F}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

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, <math display="inline">\mathbf{u}</math> es el vector de incógnitas, <math display="inline">\mathbf{F}</math> el vector de fuerzas conocido y <math display="inline">\mathbf{K}_{ij}</math>, <math display="inline">\mathbf{A}_{i}</math>, <math display="inline">\mathbf{S}</math> matrices dadas, que pueden depender de <math display="inline">\mathbf{u}</math> en problemas no lineales. Para una descripción general de cómo tratar el problema genérico descrito y sus aplicaciones en mecánica de fluidos, véase <span id='citeF-31'></span>[[#cite-31|[31]]].

===2.2 Cambio de variables===

Realizando el cambio de variables <math display="inline">u_{i}:=hU_{i}</math> y <math display="inline">P:=\frac{1}{2}g\left(h^{2}-H^{2}\right)</math>,  con <math display="inline">\partial _{t}P=hg\partial _{t}h</math>, las ecuaciones ([[#eq-1|1]]) y ([[#eq-2|2]]) se pueden escribir de la siguiente manera:

<span id="eq-6"></span>
<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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]</math>
|-
| style="text-align: center;" | <math>  \qquad +\hat{f}\left(-\delta _{i1}u_{2}+\delta _{i2}u_{1}\right)+gU_{0}{\hbox{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,\quad i=1,2,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|-
| style="text-align: center;" | <math>  \frac{1}{hg}\partial _{t}P+\partial _{i}u_{i}=0. </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}
|}

Las ecuaciones ([[#eq-6|6]]) y ([[#eq-7|7]]) se pueden escribir fácilmente en el formato de la ecuación ([[#eq-5|5]]), considerando ésta como una ecuación no lineal, con las matrices de coeficientes dependientes de la incógnita. Sin embargo, detallaremos esta identificación 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 de las ecuaciones de movimiento, consideremos una partición uniforme del intervalo de tiempo <math display="inline">0 = t^0 < t^1 <\dots < t^N = T</math>, con <math display="inline">\delta t=t^{n+1}-t^{n}</math>, donde <math display="inline">\delta t</math> es el tamaño del paso de tiempo, considerado constante. Aquí y en adelante, usaremos el superíndice <math display="inline">n</math> para indicar el nivel de tiempo&nbsp;<math display="inline">n</math>.

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 <math display="inline">f</math>, sea

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f^{n+\theta }=\theta f^{n+1}+\left(1-\theta \right)f^{n},\qquad 0\leq \theta \leq 1. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

Aproximaremos

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta _t^n  f :=  \frac{1}{\delta t} (f^{n+1} - f^n) = \frac{1}{\theta \delta t}(f^{n+\theta } - f^n) \approx (\partial _t f)^{n+ \theta } </math>
|}
|}

Esta aproximación es de segundo orden en <math display="inline">\delta t</math> en el caso en que <math display="inline">\theta = 1/2</math>  (método de Crank-Nicolson) y de primer orden si <math display="inline">\theta \not = 1/2</math>. Para <math display="inline">\theta = 1</math> corresponde al método de Euler implícito o el método de diferencias hacia atrás de primer orden, a menudo abreviado como 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

<span id="eq-9"></span>
<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta _{t}^{n}u_{i}+\partial _{j}\left(a_{j}u_{i}\right)+\partial _{i}p-\partial _{j}\hat{\tau }_{ji}+S_{ij}u_{j}+b_{i}  =0,\quad i=1,2,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|-
| style="text-align: center;" | <math> \frac{1}{h_{0}g}\delta _{t}^{n}p+\partial _{i}u_{i}  =0. </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}
|}

Usando un superíndice adicional para el contador de iteraciones, las variables introducidas en ([[#eq-9|9]])-([[#eq-10|10]]) para la iteración <math display="inline">i</math>-ésima en el paso de tiempo <math display="inline">n</math> son

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{3}  a_{j}=U_{j}^{n+\theta ,i-1}, \quad  u_{j}=\left(hU_{j}\right)^{n+\theta ,i},</math>
|-
| style="text-align: center;" | <math>  h_{0}=h^{n+\theta ,i-1},  \quad  p=P^{n+\theta ,i},</math>
|-
| style="text-align: center;" | <math>  \delta _{t}^{n}p=\frac{1}{\theta \delta t}\left(p-P^{n}\right), \quad  \delta _{t}^{n}u_{j}=\frac{1}{\theta \delta t}\left[u_{j}-\left(hU_{j}\right)^{n}\right],</math>
|-
| style="text-align: center;" | <math>  \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),</math>
|-
| style="text-align: center;" | <math>  S_{11}=S_{22}={\hbox{Ch}}^{-1}gU_{0}^{n+\theta ,i-1}h_{0}^{-3}, \quad      U_{0}^{n+\theta ,i-1}=\sqrt{\left(U_{1}^{n+\theta ,i-1}\right)^{2}+\left(U_{2}^{n+\theta ,i-1}\right)^{2}},</math>
|-
| style="text-align: center;" | <math>  S_{12}=-S_{21}=-\hat{f}, \quad  b_{i}=b_{i}^{*}+\partial _{j}\tau _{ji}^{*},</math>
|-
| style="text-align: center;" | <math>  b_{i}^{*}=-g\left(h_{0}-H\right)\partial _{i}H+\frac{1}{\rho }h_{0}\partial _{i}p_{a}-\frac{1}{\rho }\tau _{3i}^{S}. \quad   </math>
|}
|}

Obsérvese que la forma linealizada de las tensiones viscosas <math display="inline">\overline{\tau }_{ij}</math> dada en ([[#eq-3|3]]) corresponde a <math display="inline">\hat{\tau }_{ij}-\tau _{ij}^{*},</math> donde <math display="inline">\hat{\tau }_{ij}</math> está escrito en términos de las variables de la iteración actual, mientras que <math display="inline">\tau _{ij}^{*}</math> se calcula con los valores de las variables correspondientes a la iteración anterior.

===2.4 Forma vectorial===

Para escribir las ecuaciones ([[#eq-9|9]]) y ([[#eq-10|10]]) en la forma vectorial ([[#eq-5|5]]), aunque discretizada en el tiempo, escribimos las mismas como un sistema de tres ecuaciones con tres incógnitas:

<span id="eq-11"></span>
<span id="eq-12"></span>
<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta _{t}^{n}u_{1}-\frac{\mu _{H}}{\rho }\left(\frac{4}{3}\partial _{1}^{2}u_{1}+\frac{1}{3}\partial _{1}\partial _{2}u_{2}+\partial _{2}^{2}u_{1}\right)+a_{1}\partial _{1}u_{1}+a_{2}\partial _{2}u_{1}+\partial _{1}p+\left(\partial _{1}a_{1}+\partial _{2}a_{2}+S_{11}\right)u_{1}+S_{12}u_{2} </math>
|-
| style="text-align: center;" | <math>  \qquad +b_{1}^{*}+\partial _{1}\tau _{11}^{*}+\partial \tau _{21}^{*}=0,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|-
| style="text-align: center;" | <math>  \delta ^n_{t}u_{2}-\frac{\mu _{H}}{\rho }\left(\partial _{1}^{2}u_{2}+\frac{1}{3}\partial _{1}\partial _{2}u_{1}+\frac{4}{3}\partial _{2}^{2}u_{2}\right)+a_{1}\partial _{1}u_{2}+a_{2}\partial _{2}u_{2}+\partial _{2}p+\left(\partial _{1}a_{1}+\partial _{2}a_{2}+S_{22}\right)u_{2}+S_{21}u_{1}</math>
|-
| style="text-align: center;" | <math>  \qquad \,+b_{2}^{*}+\partial _{1}\tau _{12}^{*}+\partial _{2}\tau _{22}^{*}=0,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|-
| style="text-align: center;" | <math>  \frac{1}{h_{0}g}\delta _{t}^{n}p+\partial _{1}u_{1}+\partial _{2}u_{2}=0. </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}
|}

La forma vectorial discreta en el tiempo de ([[#eq-5|5]]) del sistema de ecuaciones ([[#eq-11|11]]), ([[#eq-12|12]]) y ([[#eq-13|13]]), linealizada con el método de Picard y escribiendo la divergencia del flujo convectivo como derivada convectiva, está entonces dada por

<span id="eq-46"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left[\begin{array}{c}\delta ^n_{t}u_{1}\\ \delta ^n_{t}u_{2}\\ \frac{1}{h_{0}g}\delta ^n_{t}p \end{array}\right] -\frac{\mu _{H}}{\rho }\left[\begin{array}{ccc}\frac{4}{3} & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 0 \end{array}\right]\partial _{1}^{2}\left[\begin{array}{c}u_{1}\\ u_{2}\\ p \end{array}\right]-\frac{\mu _{H}}{\rho }\left[\begin{array}{ccc}0 & \frac{1}{6} & 0\\ \frac{1}{6} & 0 & 0\\ 0 & 0 & 0 \end{array}\right]\partial _{1}\partial _{2}\left[\begin{array}{c}u_{1}\\ u_{2}\\ p \end{array}\right]</math>
|-
| style="text-align: center;" | <math>  -\frac{\mu _{H}}{\rho }\left[\begin{array}{ccc}1 & 0 & 0\\ 0 & \frac{4}{3} & 0\\ 0 & 0 & 0 \end{array}\right]\partial _{2}^{2}\left[\begin{array}{c}u_{1}\\ u_{2}\\ p \end{array}\right]-\frac{\mu _{H}}{\rho }\left[\begin{array}{ccc}0 & \frac{1}{6} & 0\\ \frac{1}{6} & 0 & 0\\ 0 & 0 & 0 \end{array}\right]\partial _{2}\partial _{1}\left[\begin{array}{c}u_{1}\\ u_{2}\\ p \end{array}\right]</math>
|-
| style="text-align: center;" | <math>  +\left[\begin{array}{ccc}a_{1} & 0 & 1\\ 0 & a_{1} & 0\\ 1 & 0 & 0 \end{array}\right]\partial _{1}\left[\begin{array}{c}u_{1}\\ u_{2}\\ p \end{array}\,\right]+\left[\begin{array}{ccc}a_{2} & 0 & 0\\ 0 & a_{2} & 1\\ 0 & 1 & 0 \end{array}\right]\partial _{2}\left[\begin{array}{c}u_{1}\\ u_{2}\\ p \end{array}\right]</math>
|-
| style="text-align: center;" | <math>  +\left[\begin{array}{ccc}\partial _{1}a_{1}+\partial _{2}a_{2}+S_{11} & S_{12} & 0\\ S_{21} & \partial _{1}a_{1}+\partial _{2}a_{2}+S_{22} & 0\\ 0 & 0 & 0 \end{array}\right]\left[\begin{array}{c}u_{1}\\ u_{2}\\ p \end{array}\right]=\left[\begin{array}{c}-Q_{1}\\ -Q_{2}\\ 0 \end{array}\right], </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}

con

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>Q_{1}  =b_{1}^{*}+\partial _{1}\tau _{11}^{*}+\partial _{2}\tau _{21}^{*},</math>
|-
| style="text-align: center;" | <math> Q_{2}  =b_{2}^{*}+\partial _{1}\tau _{12}^{*}+\partial _{2}\tau _{22}^{*}. </math>
|}
|}

La ecuación ([[#eq-46|46]]) corresponde a la ecuación vectorial de CDR transitoria ([[#eq-5|5]]) del movimiento de un fluido en aguas poco profundas discretizada en el tiempo y linealizada. Las matrices <math display="inline">\mathbf{K}_{ij}</math>, <math display="inline">\mathbf{A}_{i},</math> <math display="inline">\mathbf{S}</math> y <math display="inline">\mathbf{F}</math> son

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{K}_{11}=\frac{\mu _{H}}{\rho }\left[\begin{array}{ccc}\frac{4}{3} & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 0 \end{array}\right], \quad \mathbf{K}_{12}=\mathbf{K}_{21}=\frac{\mu _{H}}{\rho }\left[\begin{array}{ccc}0 & \frac{1}{6} & 0\\ \frac{1}{6} & 0 & 0\\ 0 & 0 & 0 \end{array}\right], \quad \mathbf{K}_{22}=\frac{\mu _{H}}{\rho }\left[\begin{array}{ccc}1 & 0 & 0\\ 0 & \frac{4}{3} & 0\\ 0 & 0 & 0 \end{array}\right],</math>
|-
| style="text-align: center;" | <math>  \mathbf{A}_{1}=\left[\begin{array}{ccc}a_{1} & 0 & 1\\ 0 & a_{1} & 0\\ 1 & 0 & 0 \end{array}\right], \quad \mathbf{A}_{2}=\left[\begin{array}{ccc}a_{2} & 0 & 0\\ 0 & a_{2} & 1\\ 0 & 1 & 0 \end{array}\right],</math>
|-
| style="text-align: center;" | <math>  \mathbf{S}=\left[\begin{array}{ccc}\partial _{1}a_{1}+\partial _{2}a_{2}+S_{11} & S_{12} & 0\\ S_{21} & \partial _{1}a_{1}+\partial _{2}a_{2}+S_{22} & 0\\ 0 & 0 & 0 \end{array}\right], \quad  \mathbf{F=}\left[\begin{array}{c}-Q_{1}\\ -Q_{2}\\ 0 \end{array}\right]. </math>
|}
|}

Obsérvese que las matrices <math display="inline">\mathbf{K}_{ij}</math> son simétricas. Los coeficientes de las matrices <math display="inline">\mathbf{A}_{i},</math> <math display="inline">\mathbf{S}</math> y <math display="inline">\mathbf{F}</math> 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. El vector de incógnitas es

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{u}^T =\left[u_{1}, u_{2}, p \right]. </math>
|}
|}

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

En esta sección vamos a describir la aproximación espacial mediante el método de los elementos finitos y estabilización numérica de la ecuación vectorial de CDR transitoria ([[#eq-5|5]]), considerando que las matrices de coeficientes están evaluadas con las incógnitas de una iteración anterior como resultado de la linealización. Esta aproximación la podremos aplicar directamente a ([[#eq-46|46]]), teniendo en cuenta que en este caso hemos discretizado el tiempo mediante un esquema de diferencias finitas y hemos escrito el término convectivo como <math display="inline">\mathbf{A}_{i}\partial _{i}\mathbf{u}</math>; las modificaciones necesarias por escribir así este término son obvias.

Es conocido que la aproximación de la ecuación estacionaria escalar de CDR por el método estándar de elementos finitos de Galerkin no es estable cuando el término convectivo es dominante con respecto al término difusivo. Para el caso vectorial sucede exactamente lo mismo, con el agravante de la no linealidad en los términos convectivos y de la compatibilidad en la interpolación de las incógnitas. Para resolver estos problemas, introduciremos los métodos de elementos finitos estabilizados ASGS y OSS basados en la formulación VMS.

===3.1 Forma diferencial del problema===

La ecuación diferencial que vamos a aproximar numéricamente es la ecuación vectorial transitoria de CDR en 2D, con condiciones de contorno de Dirichlet, Neumann y condiciones iniciales conocidas. Aunque todo lo desarrollado en esta sección es fácilmente aplicable a 3D, hemos preferido restringirnos a 2D por su directa aplicación a las ecuaciones del movimiento de un fluido en aguas poco profundas. El problema consiste en encontrar <math display="inline">\mathbf{u}:\Omega \times \gamma \rightarrow \mathbb{R}^{n_{inc}}</math>, <math display="inline">\Omega \subset \mathbb{R}^{2}</math>, <math display="inline">\gamma =\left[0,T\right]</math>, tal que

<span id="eq-61"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\partial _{t}\mathbf{u}+\mathcal{L\left(\mathbf{u}\right)} = \mathbf{F}, \qquad \left(\mathbf{x},t\right)\in \Omega \times \gamma , </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (61)
|}

con

<span id="eq-62"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{L\left(\mathbf{u}\right)}:=-\partial _{i}\left(\mathbf{K}_{ij}\partial _{j}\mathbf{u}\right)+\partial _{i}\left(\mathbf{A}_{i}\mathbf{u}\right)+\mathbf{Su}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (62)
|}

sujeta a condiciones de contorno de Dirichlet

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{u}=\mathbf{u}_{D},\qquad \left(\mathbf{x},t\right)\in \partial \Omega _{D}\times \mathbf{\gamma }, </math>
|}
|}

condiciones de contorno de Neumann

<span id="eq-63"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>n_{i}\mathbf{A}_{i}^{f}\mathbf{u}-n_{i}\mathbf{K}_{ij}\partial _{j}\mathbf{u}={\boldsymbol{}\varphi },\qquad \left(\mathbf{x},t\right)\in \partial \Omega _{N}\times \mathbf{\gamma }, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (63)
|}

y condiciones iniciales

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{u}=\mathbf{u}_{0},\qquad \mathbf{x}\in \Omega ,\quad t=0. </math>
|}
|}

Hemos considerado que la frontera de <math display="inline">\Omega </math> (denotada por <math display="inline">\partial \Omega </math>) está dada por <math display="inline">\partial \Omega =\partial \Omega _{D}+\partial \Omega _{N}</math> y las matrices de convección <math display="inline">\mathbf{A}_{i}</math> se pueden descomponer en dos términos, uno propiamente convectivo <math display="inline">\mathbf{A}_{i}^{c}</math> y otro que contribuye al flujo de la incógnita, <math display="inline">\mathbf{A}_{i}^{f}</math>, de manera que <math display="inline">\mathbf{A}_{i}=\mathbf{A}_{i}^{c}+\mathbf{A}_{i}^{f}</math>.

La ecuación ([[#eq-61|61]]) es un sistema de <math display="inline">n_{inc}</math> ecuaciones diferenciales en derivadas parciales, con <math display="inline">n_{inc} = 3</math> para las ecuaciones de aguas poco profundas. Las matrices de difusividad <math display="inline">\mathbf{K}_{ij}</math>, convección <math display="inline">\mathbf{A}_{i}</math> y de reacción <math display="inline">\mathbf{S}</math> son matrices cuadradas de orden <math display="inline">n_{inc}\times n_{inc}</math>, mientras que <math display="inline">\mathbf{u}</math> y <math display="inline">\mathbf{F}</math> son vectores de <math display="inline">n_{inc}</math> componentes.

Suponemos que las condiciones iniciales y de contorno de Dirichlet verifican la condición de compatibilidad <math display="inline"> \mathbf{u}_{0}\left(\mathbf{x}\right)|_{\partial \Omega _{D}}=\mathbf{u}_{D}\left(\mathbf{x},0\right)</math>.

===3.2 Forma débil del problema estacionario===

Para simplificar la exposición, empecemos considerando el problema estacionario de la ecuación vectorial de CDR, con condiciones de Dirichlet homogéneas en toda su frontera. Dicho problema consiste en hallar <math display="inline">\mathbf{u}:\Omega \rightarrow \mathbb{R}^{n_{inc}}</math> tal que

<span id="eq-64"></span>
<span id="eq-65"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{L\left(\mathbf{u}\right)}:=-\partial _{i}\left(\mathbf{K}_{ij}\partial _{j}\mathbf{u}\right)+\partial _{i}\left(\mathbf{A}_{i}\mathbf{u}\right)+\mathbf{Su}  =\mathbf{F},\qquad \mathbf{x}\in \Omega ,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (64)
|-
| style="text-align: center;" | <math> \qquad \qquad \qquad \qquad \qquad \mathbf{u}  =\mathbf{0},\mathbf{\qquad x}\in \partial \Omega{.}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
|}
|}

Para escribir la forma débil del problema ([[#eq-64|64]])-([[#eq-65|65]]), sea <math display="inline">\mathcal{W}</math> el espacio donde la ecuación diferencial tiene sentido distribucional. Por ejemplo, si <math display="inline">\mathbf{K}_{ij}</math> define una forma cuadrática definida positiva (que no es el caso en las ecuaciones de aguas poco profundas), dicho espacio es <math display="inline">\mathcal{W}:=\left(H_{0}^{1}(\Omega \right))^{n_{inc}}</math>, es decir,  el espacio de funciones vectoriales cuyas componentes son funciones cuadrado integrable, que tienen primera derivada cuadrado integrable y que se anulan en la frontera <math display="inline">\partial \Omega </math>. Suponemos también que <math display="inline">\mathbf{F}</math> multiplicado por funciones de <math display="inline">\mathcal{W}</math> es integrable, y que las matrices <math display="inline">\mathbf{K}_{ij},</math> <math display="inline">\mathbf{A}_{i}</math> y <math display="inline">\mathbf{S}</math> tienen coeficientes acotados.

Podemos ahora escribir la forma débil del problema ([[#eq-64|64]])-([[#eq-65|65]]) como sigue: hallar <math display="inline">\mathbf{u}\in \mathcal{W}</math> tal que

<span id="eq-66"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int _{\Omega }\partial _{i}\mathbf{v}^{T}\mathbf{K}_{ij}\partial _{j}\mathbf{u}{\,d}\Omega{+\int}_{\Omega }\mathbf{v}^{T}\partial _{i}\left(\mathbf{A}_{i}^{c}\mathbf{u}\right){\,d}\Omega{-\int}_{\Omega }\partial _{i}\mathbf{v}^{T}\mathbf{A}_{i}^{f}\mathbf{u}{\,d}\Omega{+\int}_{\Omega }\mathbf{v}^{T}\mathbf{Su}{\,d}\Omega </math>
|-
| style="text-align: center;" | <math>  =\int _{\Omega }\mathbf{v}^{T}\mathbf{F}{\,d}\Omega{+\int}_{\partial \Omega }(\mathbf{v}^{T}\mathbf{K}_{ij}\partial _{j}\mathbf{u})n_{i}{\,d}\Gamma{-\int}_{\partial \Omega }(\mathbf{v}^{T}\mathbf{A}_{i}^{f}\mathbf{u})n_{i}{\,d}\Gamma ,\qquad \forall \mathbf{v}\in \mathcal{W}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
|}

donde <math display="inline">n_i</math> es la componente <math display="inline">i</math>-ésima de la normal unitaria a <math display="inline">\partial \Omega </math>.

Los dos últimos términos del lado derecho de la ecuación ([[#eq-66|66]]) son valores conocidos, pues o bien la función de test <math display="inline">\mathbf{v}</math> es nula en el caso en el que <math display="inline">\partial \Omega _D = \partial \Omega </math> o bien usamos la condición de Neumann ([[#eq-63|63]]). Nos restringiremos al primer caso, como hemos indicado, con lo cual la forma variacional o débil del problema ([[#eq-64|64]])-([[#eq-65|65]])  consiste en encontrar <math display="inline">\mathbf{u}\in \mathcal{W}</math> tal que

<span id="eq-67"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>B(\mathbf{u},\mathbf{v})=\left\langle \mathbf{F},\mathbf{v}\right\rangle \qquad \forall \mathbf{v}\in \mathcal{W}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (67)
|}

donde la forma <math display="inline">B(\mathbf{u},\mathbf{v})</math> y la forma lineal <math display="inline">\left\langle \mathbf{F},\mathbf{v}\right\rangle </math> están dadas por

<span id="eq-68"></span>
<span id="eq-69"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>B(\mathbf{u},\mathbf{v})  =\int _{\Omega }\partial _{i}\mathbf{v}^{T}\mathbf{K}_{ij}\partial _{j}\mathbf{u}{\,d}\Omega{+\int}_{\Omega }\mathbf{v}^{T}\parti al _{i}\left(\mathbf{A}_{i}^{c}\mathbf{u}\right){\,d}\Omega{-\int}_{\Omega }\partial _{i}\mathbf{v}^{T}\mathbf{A}_{i}^{f}\mathbf{u}{\,d}\Omega{+\int}_{\Omega }\mathbf{v}^{T}\mathbf{Su}{\,d}\Omega , </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
|-
| style="text-align: center;" | <math>  \left\langle \mathbf{F},\mathbf{v}\right\rangle  =\int _{\Omega }\mathbf{v}^{T}\mathbf{F}{\,d}\Omega{.} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (69)
|}
|}

===3.3 Método estándar de Galerkin para el problema estacionario===

Para aproximar numéricamente la ecuación ([[#eq-67|67]]) en el espacio por el método de los elementos finitos, consideremos una partición <math display="inline">\{ \Omega ^e\} </math> de <math display="inline">\Omega </math>, con <math display="inline">e = 1,\dots , n_{el}</math>, siendo <math display="inline">n_{el}</math> el número de elementos de la partición. Sea <math display="inline">\mathcal{W}_{h}</math> un espacio conforme de elementos finitos para aproximar <math display="inline">\mathcal{W}</math> <math display="inline">\left(\mathcal{W}_{h}\subset \mathcal{W}\right)</math>. El problema discreto que se conoce como el método estándar de Galerkin, que aproxima la ecuación vectorial de CDR estacionaria ([[#eq-65|65]]), consiste en hallar <math display="inline">\mathbf{u}_{h}\in \mathcal{W}_{h}</math> tal que

<span id="eq-70"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>B(\mathbf{u}_{h},\mathbf{v}_{h})=\left\langle \mathbf{F},\mathbf{v}_{h}\right\rangle \qquad \forall \mathbf{v}_{h}\in \mathcal{W}_{h}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (70)
|}

En general, es necesario que la interpolación de las incógnitas satisfaga ciertas condiciones de compatibilidad (usualmente expresadas en forma de condiciones inf-sup). Sin embargo, con la estabilización que introduciremos más adelante esto no es necesario, por lo que presentamos a continuación algunos aspectos de implementación usando igual interpolación.

Observemos primero que podemos escribir

<span id="eq-71"></span>
<span id="eq-72"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>B(\mathbf{u}_{h},\mathbf{v}_{h}):=\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\left(\partial _{i}\mathbf{v}_{h}^{T}\mathbf{K}_{ij}\partial _{j}\mathbf{u}_{h}+\mathbf{v}_{h}^{T}\partial _{i}\left(\mathbf{A}_{i}^{c}\mathbf{u}_{h}\right)-\partial _{i}\mathbf{v}_{h}^{T}\mathbf{A}_{i}^{f}\mathbf{u}_{h}+\mathbf{v}_{h}^{T}\mathbf{Su}_{h}\right){\,d}\Omega ,</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (71)
|-
| style="text-align: center;" | <math>  \left\langle \mathbf{F},\mathbf{v}_{h}\right\rangle =\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\mathbf{v}_{h}^{T}\mathbf{F}{\,d}\Omega{.} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (72)
|}
|}

Escribamos <math display="inline">\mathbf{u}_{h}=\left(\hat{u}^{1},\hat{u}^{2},...,\hat{u}^{n_{inc}}\right)^{T}</math>. Si usamos las mismas funciones de interpolación de clase <math display="inline">C^0</math> (lagrangianas) para todos las componentes de la incógnita, podemos escribir

<span id="eq-73"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hat{u}^{q}({\mathbf{x}})=\overset{M}{\underset{m=1}{\sum }}u_{m}^{q}N_{m}\left(\mathbf{x}\right),\qquad q=1,...,n_{inc}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (73)
|}

donde <math display="inline">M</math> es el número de nodos de la malla de elementos finitos, <math display="inline">N_{m}\left(\mathbf{x}\right)</math> son las funciones de forma y <math display="inline">u_{m}^{q}</math> es el valor nodal de la componente <math display="inline">q</math> de la incógnita en el nodo <math display="inline">m</math> de la malla. Además, utilizando el método de Galerkin tomamos como funciones de prueba <math display="inline">\mathbf{v}_{h}</math> las dadas por:

<span id="eq-74"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{v}_{h}^{T}=\mathbf{v}_{h,p}^{T}=\left(\delta _{1p}N_{l},\delta _{2p}N_{l},...,\delta _{n_{inc}p}N_{l}\right),\quad p=1, \dots ,n_{inc},\quad l = 1,\dots , M, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (74)
|}

donde <math display="inline">\delta _{ip}</math> es la delta de Kronecker. Reemplazando ([[#eq-73|73]]) y ([[#eq-74|74]]) en ([[#eq-70|70]]), ([[#eq-71|71]]) y ([[#eq-72|72]]) se tiene:

<span id="eq-75"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\underset{e=1}{\overset{n_{el}}{\sum }} \sum _{q= 1}^{n_{inc}}\overset{M}{\underset{m=1}{\sum }} \left[ \int _{\Omega ^{e}} \partial _{i}N_{l}\left(K_{ij}\right)_{pq}\partial _{j} N_{m}{\,d}\Omega  +\int _{\Omega ^{e}} N_{l}\left(A_{i}^{c}\right)_{pq}\partial _{i} N_{m} {\,d}\Omega \right.</math>
|-
| style="text-align: center;" | <math>  \left.-\int _{\Omega ^{e}}\partial _{i}N_{l}\bigl(A_{i}^{f}\bigr)_{pq} N_{m}{\,d}\Omega  +\int _{\Omega ^{e}} N_{l}S_{pq}N_{m}{\,d}\Omega \right]u_{m}^{q} =\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}N_{l}f_{p}{\,d}\Omega ,</math>
|-
| style="text-align: center;" | <math>  p = 1,\dots ,n_{inc}, ~l = 1,\dots , M. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (75)
|}

En esta expresión, <math display="inline">(K_{ij})_{pq},</math> <math display="inline">(A_{i}^c)_{pq}</math>, <math display="inline">(A_{i}^{f})_{pq}</math> y <math display="inline">(S)_{pq}</math> son los coeficientes del elemento de fila <math display="inline">p</math>-ésima y columna <math display="inline">q</math>-ésima de las correspondientes matrices <math display="inline">\mathbf{K}_{ij}</math>, <math display="inline">\mathbf{A}_{i}^{c}</math>, <math display="inline">\mathbf{A}_{i}^{f}</math> y <math display="inline">\mathbf{S}</math>. El ensamblaje de la ecuación ([[#eq-75|75]]) da lugar, una vez impuestas las condiciones de contorno, a un sistema lineal de ecuaciones cuya solución proporciona los valores nodales de la aproximación numérica del problema estacionario ([[#eq-64|64]])-([[#eq-65|65]]).

===3.4 Método de las subescalas, problema estacionario===

Como ya dijimos antes, el método estándar de Galerkin presenta inestabilidades numéricas cuando el término convectivo es dominante con respecto al término difusivo y debido a las condiciones de compatibilidad entre la interpolación de las variables. El método que vamos a utilizar para la estabilización numérica es el método variacional de subescalas (VMS) introducido por Hughes y otros en <span id='citeF-18'></span><span id='citeF-19'></span>[[#cite-18|[18,19]]] para resolver las inestabilidades numéricas del método de Galerkin en el problema de CDR escalar estacionario y generalizado  en <span id='citeF-20'></span><span id='citeF-21'></span>[[#cite-20|[20,21]]] al problema de CDR vectorial estacionario.

La idea básica es descomponer la variable continua desconocida <math display="inline">\mathbf{u}\in \mathcal{W}</math> en una parte resoluble en el espacio de elementos finitos <math display="inline">\mathbf{u}_{h}\in \mathcal{W}_{h},</math> con <math display="inline">\mathcal{W}_{h}\subset \mathcal{W}</math>, y una parte en la escala de una submalla <math display="inline">\tilde{\mathbf{u}}\in \tilde{\mathcal{W}}</math>, <math display="inline">\tilde{\mathcal{W}}\subset \mathcal{W},</math> la cual no puede ser capturada por la malla de elementos finitos, siendo <math display="inline">\mathcal{W}=\mathcal{W}_{h}\oplus \tilde{\mathcal{W}}</math>, donde <math display="inline">\mathcal{\tilde{W}}</math> es cualquier espacio para completar <math display="inline">\mathcal{W}_{h}</math> en <math display="inline">\mathcal{W}</math>. Para evitar tecnicismos, podemos pensar de <math display="inline">\mathcal{W}</math> y <math display="inline">\mathcal{\tilde{W}}</math> como espacios dimensionalmente finitos, con una dimensión grande. Puesto que <math display="inline">\tilde{\mathcal{W}}</math> representa la componente de <math display="inline">\mathcal{W}</math> que no puede ser capturada por el espacio de elementos finitos, la llamamos espacio de las subescalas. De este modo, la variable continua desconocida podemos escribirla como <math display="inline">\mathbf{u}=\mathbf{u}_{h}+\tilde{\mathbf{u}}</math>, donde <math display="inline">\mathbf{u}_{h}</math> es la componente de <math display="inline">\mathbf{u}</math> en el espacio de elementos finitos y <math display="inline">\tilde{\mathbf{u}}</math> es su componente en el complemento (con respecto a un cierto producto interno) en <math display="inline">\mathcal{W}</math>.

Aplicando la descomposición mencionada a la ecuación continua ([[#eq-67|67]]) se obtiene:

<span id="eq-76"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>B(\mathbf{u}_{h}+\tilde{\mathbf{u}},\mathbf{v}_{h}+\tilde{\mathbf{v}})=\left\langle \mathbf{F},\mathbf{v}_{h}+\tilde{\mathbf{v}}\right\rangle \qquad \mathbf{\forall v}_{h}\in \mathcal{W}_{h},~\forall \tilde{\mathbf{v}}\in \tilde{\mathcal{W}}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (76)
|}

Hemos considerado el problema linealizado, por lo cual <math display="inline">B</math> es una forma bilineal (para el tratamiento generalizado de problemas no lineales, véase <span id='citeF-22'></span><span id='citeF-32'></span>[[#cite-22|[22,32]]]). Utilizando este hecho, la ecuación ([[#eq-76|76]]) se puede descomponer en el siguiente sistema de ecuaciones:

<span id="eq-77"></span>
<span id="eq-78"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>B\left(\mathbf{u}_{h},\mathbf{v}_{h}\right)+B\left(\tilde{\mathbf{u}},\mathbf{v}_{h}\right)=\left\langle \mathbf{F},\mathbf{v}_{h}\right\rangle \qquad \mathbf{\forall v}_{h}\in \mathcal{W}_{h},</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (77)
|-
| style="text-align: center;" | <math>  B\left(\mathbf{u}_{h},\tilde{\mathbf{v}}\right)+B\left(\tilde{\mathbf{u}},\tilde{\mathbf{v}}\right)=\left\langle \mathbf{F},\tilde{\mathbf{v}}\right\rangle \qquad \forall \tilde{\mathbf{v}}\in \tilde{\mathcal{W}}. </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (78)
|}
|}

La ecuación ([[#eq-77|77]]) corresponde a la escala resoluble <math display="inline">\mathbf{u}_{h}</math> en el espacio de elementos finitos y tiene dos términos a su lado izquierdo, el primero es la contribución del método de Galerkin y el segundo toma en cuenta la influencia de la subescala en <math display="inline">\mathbf{u}_{h}</math>. De la ecuación ([[#eq-78|78]]) debe obtenerse una aproximación a <math display="inline">\tilde{\mathbf{u}}</math>.

De forma análoga a la ecuación ([[#eq-71|71]]), para el segundo término del lado izquierdo de la ecuación ([[#eq-77|77]]) tenemos

<span id="eq-79"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>B\left(\tilde{\mathbf{u}},\mathbf{v}_{h}\right)=\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\left(\partial _{i}\mathbf{v}_{h}^{T}\mathbf{K}_{ij}\partial _{j}\mathbf{\tilde{\mathbf{u}}}+\mathbf{v}_{h}^{T}\partial _{i}\left(\mathbf{A}_{i}^{c}\mathbf{\tilde{\mathbf{u}}}\right)-\partial _{i}\mathbf{v}_{h}^{T}\mathbf{A}_{i}^{f}\mathbf{\tilde{\mathbf{u}}}+\mathbf{v}_{h}^{T}\mathbf{S\tilde{\mathbf{u}}}\right){\,d}\Omega , </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (79)
|}

e integrando por partes dentro de cada elemento y recordando que <math display="inline">\mathbf{A}_{i}=\mathbf{A}_{i}^{c}+\mathbf{A}_{i}^{f}</math> tenemos

<span id="eq-80"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>B\left(\tilde{\mathbf{u}},\mathbf{v}_{h}\right) =\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\left(-\partial _{j}\left(\partial _{i}\mathbf{v}_{h}^{T}\mathbf{K}_{ij}\right)\tilde{\mathbf{u}}-\partial _{i}\mathbf{v}_{h}^{T}\mathbf{A}_{i}\mathbf{\tilde{\mathbf{u}}}+\mathbf{v}_{h}^{T}\mathbf{S}\mathbf{\tilde{\mathbf{u}}}\right){\,d}\Omega </math>
|-
| style="text-align: center;" | <math>  +\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\partial \Omega ^{e}}\partial _{i}\mathbf{v}_{h}^{T}\mathbf{K}_{ij}n_{j}\tilde{\mathbf{u}}{\,d}\Gamma +\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\partial \Omega ^{e}}\mathbf{v}_{h}^{T}n_{i}\mathbf{A}_{i}^{c}\mathbf{\tilde{u}}{\,d}\Gamma{.} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (80)
|}

Sea <math display="inline">\mathcal{L}^{*}</math> el adjunto del operador <math display="inline">\mathcal{L}</math>, dado por

<span id="eq-81"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\mathcal{L}^{*}\left(\mathbf{v}_{h}\right)\right)^{T}=-\partial _{j}\left(\partial _{i}\mathbf{v}_{h}^{T}\mathbf{K}_{ij}\right)-\partial _{i}\mathbf{v}_{h}^{T}\mathbf{A}_{i}+\mathbf{v}_{h}^{T}\mathbf{S}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (81)
|}

donde hemos usado que en nuestro caso <math display="inline">\mathbf{K}^T_{ji} = \mathbf{K}_{ij}</math>. La ecuación ([[#eq-80|80]]) se puede escribir como

<span id="eq-82"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>B\left(\tilde{\mathbf{u}},\mathbf{v}_{h}\right)=\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\left(\mathcal{L}^{*}\left(\mathbf{v}_{h}\right)\right)^{T}\tilde{\mathbf{u}}{\,d}\Omega +\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\partial \Omega ^{e}}\partial _{i}\mathbf{v}_{h}^{T}\mathbf{K}_{ij}n_{j}\tilde{\mathbf{u}}{\,d}\Gamma +\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\partial \Omega ^{e}}\mathbf{v}_{h}^{T}n_{i}\mathbf{A}_{i}^{c}\mathbf{\tilde{u}}{\,d}\Gamma{.} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (82)
|}

Reemplazando ([[#eq-82|82]]) en la ecuación ([[#eq-77|77]]), ésta se puede escribir como

<span id="eq-83"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>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}\tilde{\mathbf{u}}{\,d}\Omega +\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\partial \Omega ^{e}}\left(\partial _{i}\mathbf{v}_{h}^{T}\mathbf{K}_{ij}n_{j}\tilde{\mathbf{u}}+\mathbf{v}_{h}^{T}n_{i}\mathbf{A}_{i}^{c}\mathbf{\tilde{u}}\right){\,d}\Gamma =\left\langle \mathbf{F},\mathbf{v}_{h}\right\rangle \quad \forall \mathbf{v}_{h}\in \mathcal{W}_{h}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (83)
|}

Ahora bien, integrando también por partes dentro de cada elemento los dos términos de la ecuación ([[#eq-78|78]]), ésta queda

<span id="eq-84"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\partial \Omega ^{e}}\tilde{\mathbf{v}}^{T}n_{i}\mathbf{K}_{ij}\partial _{j}\left(\mathbf{u}_{h}+\tilde{\mathbf{u}}\right){\,d}\Gamma +\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\tilde{\mathbf{v}}^{T}\mathcal{L}\left(\tilde{\mathbf{u}}\right){\,d}\Omega =\underset{e=1}{\overset{n_{el}}{\sum }}\int _{\Omega ^{e}}\tilde{\mathbf{v}}^{T}\left[\mathbf{F}-\mathcal{L}\left(\mathbf{u}_{h}\right)\right]{\,d}\Omega \qquad \forall \tilde{\mathbf{v}}\in \tilde{\mathcal{W}}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (84)
|}

Observemos que el primer término de la ecuación ([[#eq-84|84]]) se anula, ya que en la frontera de los elementos las componentes normales de los flujos de <math display="inline">\mathbf{u}</math> son iguales pero de signo contrario, debido a que los flujos difusivos deben ser continuos a través de los contornos inter-elementales. Como resultado de esto, la ecuación ([[#eq-84|84]]) es equivalente a encontrar <math display="inline">\tilde{\mathbf{u}}\in \tilde{\mathcal{W}}</math> tal que

<span id="eq-85"></span>
<span id="eq-86"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{L}\left(\tilde{\mathbf{u}}\right) =\mathbf{F}-\mathcal{L}\left(\mathbf{u}_{h}\right)+\mathbf{v}_{h,{ort}}~~\mathrm{en}\;\Omega ^{e},~~\forall \mathbf{v}_{h,{ort}}\in \tilde{\mathcal{W}}^{\perp }, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (85)
|-
| style="text-align: center;" | <math> \tilde{\mathbf{u}}  = \tilde{\mathbf{u}}_{esq}~~\mathrm{en}\;\partial \Omega ^{e}, </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (86)
|}
|}

para <math display="inline">e=1,...,n_{el}</math>, donde <math display="inline">\tilde{\mathbf{u}}_{esq}</math> es la subescala evaluada en los contornos interelementales. Es importante observar que la ecuación ([[#eq-85|85]]) se verifica para cualquier elemento <math display="inline">\mathbf{v}_{h,{ort}}</math> ortogonal a <math display="inline">\tilde{\mathcal{W}}.</math> La presencia del elemento <math display="inline">\mathbf{v}_{h,{ort}}</math> en la ecuación ([[#eq-85|85]]) garantiza que <math display="inline">\mathcal{L}\left(\tilde{\mathbf{u}}\right)-\left[\mathbf{F}-\mathcal{L}\left(\mathbf{u}_{h}\right)\right]</math> pertenece a <math display="inline">\tilde{\mathcal{W}}^{\perp }</math>, que es lo que implica la ecuación ([[#eq-84|84]]) sin el primer término. En otras palabras <math display="inline">\mathbf{v}_{h,{ort}}</math> garantiza que <math display="inline">\tilde{\mathbf{u}}</math> pertenezca al espacio de subescalas <math display="inline">\tilde{\mathcal{W}}</math> (y no a todo el espacio <math display="inline">\mathcal{W}</math>).

Las ecuaciones ([[#eq-83|83]]), ([[#eq-85|85]]) y ([[#eq-86|86]]) equivalen exactamente a los problemas ([[#eq-77|77]]) y ([[#eq-78|78]]), ya que hasta ahora no hemos hecho ninguna aproximación. Veremos ahora posibles aproximaciones para poder calcular de forma efectiva <math display="inline">\tilde{\mathbf{u}}</math>, las cuales darán lugar a los métodos ASGS y OSS. En esencia, lo que debemos hacer es seleccionar las funciones <math display="inline">\tilde{\mathbf{u}}_{esq}</math>, <math display="inline">\mathbf{v}_{h,{ort}}</math> y proponer una solución aproximada de la ecuación ([[#eq-85|85]]).

====3.4.1 El método ASGS para el problema de CDR vectorial estacionario====

Un caso particular del método de las subescalas descrito anteriormente es el método denominado ASGS. Para seleccionar <math display="inline">\tilde{\mathbf{u}}_{esq}</math> vamos a usar el siguiente argumento. Para tener una aproximación discreta óptima de <math display="inline">\mathbf{u}</math> que diera valores nodales exactos, se podría pedir que las subescalas se anularan en los nodos. En problemas unidimensionales esto da condiciones de contorno homogéneas en ([[#eq-86|86]]). Para los casos de más de una dimensión espacial, la elección de <math display="inline">\tilde{\mathbf{u}}_{esq}=\mathbf{0}</math> es solo una aproximación. En otras palabras, estamos tomando el espacio de las subescalas <math display="inline">\tilde{\mathcal{W}}</math> como un espacio de funciones burbuja, es decir, el espacio de funciones que se anulan en los contornos de los elementos (véase <span id='citeF-33'></span><span id='citeF-34'></span>[[#cite-33|[33,34]]]).  Esta aproximación, que adoptaremos en lo que sigue, no es sin embargo imprescindible (véase <span id='citeF-35'></span>[[#cite-35|[35]]]).

Con esta aproximación, el término de las integrales en los contornos de los elementos de la ecuación ([[#eq-83|83]]) se elimina, por lo cual el problema que debemos resolver es encontrar <math display="inline">\mathbf{u}_{h}\in \mathcal{W}_{h}</math> tal que

<span id="eq-87"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>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}\tilde{\mathbf{u}}{\,d}\Omega =\left\langle \mathbf{F},\mathbf{v}_{h}\right\rangle \qquad \forall \mathbf{v}_{h}\in \mathcal{W}_{h}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (87)
|}

en donde todavía hace falta determinar <math display="inline">\tilde{\mathbf{u}}</math>.

La siguiente aproximación del método ASGS es tomar <math display="inline">\mathbf{v}_{h,{ort}}=\mathbf{0}</math>, lo cual supone considerar que las subescalas están en el espacio de residuos, es decir, éste está compuesto por funciones de la forma <math display="inline">\mathbf{F}-\mathcal{L}\left(\mathbf{v}\right)</math>, con <math display="inline">\mathbf{v}\in \mathcal{W}</math>. Podemos entonces escribir simbólicament la solución del problema ([[#eq-85|85]])-([[#eq-86|86]]), asumiendo <math display="inline">\mathbf{\tilde{u}}_{esq}=\mathbf{0}</math> y considerando <math display="inline">\mathbf{v}_{h,{ort}}=\mathbf{0}</math>, como

<span id="eq-88"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tilde{\mathbf{u}}=\mathcal{L}^{-1}\left[\mathbf{F}-\mathcal{L}\left(\mathbf{u}_{h}\right)\right], </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (88)
|}

donde <math display="inline">\mathcal{L}^{-1}</math> es el operador inverso del operador <math display="inline">\mathcal{L}</math> con condiciones de Dirichlet homogéneas.

De ([[#eq-87|87]]) se observa que solo se necesita la componente de <math display="inline">\mathbf{\tilde{u}}</math> sobre el espacio <math display="inline">\mathcal{L}\left(\mathcal{W}_{h}\right)</math>, donde <math display="inline">\mathcal{L}\left(\mathcal{W}_{h}\right)</math> es el espacio de funciones de la forma <math display="inline">\mathcal{L}\left(\mathbf{v}_{h}\right)</math>, con <math display="inline">\mathbf{v}_{h}\in \mathcal{W}_{h}</math>; esto sugiere aproximar ([[#eq-88|88]]) como

<span id="eq-89"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tilde{\mathbf{u}}\approx {\tau }^{e}\left[\mathbf{F}-\mathcal{L}\left(\mathbf{u}_{h}\right)\right], </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (89)
|}

donde <math display="inline">{{\tau }}^{e}</math> es una matriz de <math display="inline">n_{inc}\times n_{inc}</math> componentes a la que denominamos matriz de parámetros de estabilización y que se calcula dentro del dominio de cada elemento. El diseño de la matriz <math display="inline">{{\tau }}^{e}</math> es una de las piedras angulares en el desarrollo de los métodos de elementos finitos estabilizados. Algunas formulaciones algebraicas para el diseño del parámetro de estabilización para la ecuación escalar de CDR en 1D, basadas en el principio del máximo discreto, pueden encontrarse en <span id='citeF-18'></span><span id='citeF-19'></span><span id='citeF-26'></span>[[#cite-18|[18,19,26]]]; la extensión a sistemas puede consultarse en <span id='citeF-20'></span>[[#cite-20|[20]]]. Asimismo, en <span id='citeF-27'></span><span id='citeF-22'></span>[[#cite-27|[27,22]]] se propone un diseño del parámetro de estabilización basado en un análisis de Fourier, considerando el problema de flujo en aguas poco profundas en <span id='citeF-12'></span>[[#cite-12|[12]]]. En este trabajo vamos a utilizar precisamente el diseño algorítmico del <math display="inline">{\tau }^e</math> desarrollado en <span id='citeF-12'></span>[[#cite-12|[12]]] y que está dado por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{{\tau }}^{e}=\mathrm{diag}(\tau ^e,\tau ^e,0), </math>
|}
|}

con el parámetro <math display="inline">\tau ^e</math> dado por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau ^e=\left[c_{1}\frac{\nu _{H}}{\left(h^{e}\right)^{2}}+c_{2}\frac{|\mathbf{a}|}{h^{e}}+c_{3}|S_{11}|+c_{4}|S_{12}|\right]^{-1}, </math>
|}
|}

para elementos finitos lineales, siendo <math display="inline">\nu _{H}=\frac{\mu _{H}}{\rho }</math> el coeficiente de viscosidad cinemática y <math display="inline">|\mathbf{a}|</math> la norma de la velocidad. Se entiende que las funciones que aparecen en esta expresión se calculan dentro de cada elemento <math display="inline">e</math>, cuyo diámetro es <math display="inline">h^e</math>, y <math display="inline">c_1</math>, <math display="inline">c_2</math>, <math display="inline">c_3</math> y <math display="inline">c_4</math> son constantes algorítmicas.

Para elementos de alto orden, la experimentación numérica y también el análisis teórico (que está fuera del objetivo de este artículo) nos han llevado a proponer como matriz de parámetros de estabilización:

<span id="eq-90"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>{{\tau }}^{e}=\mathrm{diag}(\tau ^e_{1},\tau ^e_{1},\tau ^e_{2}), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (90)
|}

con

<span id="eq-91"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau ^e_{1}=\left[c_{1}\frac{\nu _{H}}{\left(\frac{h^{e}}{d^{2}}\right)^{2}}+c_{2}\frac{|\mathbf{a}|}{\frac{h^{e}}{d}}+c_{3}|S_{11}|+c_{4}|S_{12}|\right]^{-1},\mbox{ }\tau ^e_{2}=\frac{\left(\frac{h^{e}}{d}\right)^{2}}{c_{1}\tau _{1}}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (91)
|}

donde <math display="inline">c_{1},</math> <math display="inline">c_{2},</math> <math display="inline">c_{3}</math> y <math display="inline">c_{4}</math> son los mismos coeficientes que antes, y que para los ejemplos numéricos hemos adoptado como <math display="inline">c_{1}=12,</math> <math display="inline">c_{2}=2</math> y <math display="inline">c_{3}=c_{4}=1</math> y <math display="inline">d</math> ''es el orden polinomial de interpolación''.

Reemplazando la ecuación ([[#eq-89|89]]) en ([[#eq-87|87]]), finalmente tenemos la formulación estabilizada ASGS para resolver la ecuación vectorial de CDR estacionaria, la cual consiste en encontrar <math display="inline">\mathbf{u}_{h}\in \mathcal{W}_{h}</math> tal que

<span id="eq-92"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>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 {\tau }^{e}\left[\mathcal{L}\left(\mathbf{u}_{h}\right)-\mathbf{F}\right]{\,d}\Omega =\left\langle \mathbf{F},\mathbf{v}_{h}\right\rangle \qquad \forall \mathbf{v}_{h}\in \mathcal{W}_{h}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (92)
|}

====3.4.2 El método OSS para el problema de CDR vectorial estacionario====

A continuación vamos a describir el método OSS propuesto en <span id='citeF-23'></span>[[#cite-23|[23]]]. El punto de partida es la elección del espacio de las subescalas, que ahora no será el espacio de residuos de elementos finitos como en el caso anterior, sino que se toma como espacio <math display="inline">\mathcal{\tilde{W}}</math> el espacio ortogonal en el sentido de <math display="inline">L^2</math> al espacio de elementos finitos, es decir,

<span id="eq-93"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{\tilde{W}=\mathcal{W}}_{h}^{\perp }, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (93)
|}

que es una legítima elección que cumple con <math display="inline">\mathcal{W}=\mathcal{W}_{h}\oplus \mathcal{\tilde{W}}.</math>

Tomando en cuenta ([[#eq-93|93]]) de ([[#eq-85|85]]) se sigue que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{v}_{h,{ort}}\in \mathcal{W}_{h},\quad  \mathbf{\tilde{u}}\in \mathcal{W}_{h}^{\perp },  </math>
|}
|}

lo que significa que <math display="inline">\mathbf{v}_{h,{ort}}</math> es una función del espacio de elementos finitos y por lo tanto numéricamente calculable, mientras que <math display="inline">\mathbf{\tilde{u}}</math> estará en el espacio ortogonal a <math display="inline">\mathcal{W}_{h}</math>. A este espacio <math display="inline">\mathcal{W}_{h}^{\perp }</math> escogido para las subescalas <math display="inline">\mathbf{\tilde{u}}</math> lo denominamos espacio de subescalas ortogonales.

Para resolver ([[#eq-85|85]]) y ([[#eq-86|86]]), al igual que en el método ASGS, asumimos <math display="inline">\mathbf{\tilde{u}}_{esq}=\mathbf{0}</math> y escribimos

<span id="eq-94"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{\tilde{u}}=\mathcal{L}^{-1}\left[\mathbf{F}-\mathcal{L}\left(\mathbf{u}_{h}\right)\right]+\mathcal{L}^{-1}\left(\mathbf{v}_{h,{ort}}\right), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (94)
|}

donde, como en el método ASGS, <math display="inline">\mathcal{L}^{-1}</math> es el inverso del operador <math display="inline">\mathcal{L}</math> con condiciones de Dirichlet homogéneas. De ([[#eq-87|87]]) se observa que solo se necesita la componente de <math display="inline">\mathbf{\tilde{u}}</math> sobre el espacio <math display="inline">\mathcal{L}\left(\mathcal{W}_{h}\right)</math> y, al igual que en el método ASGS, esto sugiere aproximar ([[#eq-94|94]]) como

<span id="eq-95"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{\tilde{u}}\approx {{\tau }}^{e}\left[\mathbf{F}-\mathcal{L}\left(\mathbf{u}_{h}\right)\right]+ {{\tau }}^{e}\mathbf{v}_{h,{ort}}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (95)
|}

Puesto que <math display="inline">\mathbf{\tilde{u}}\in \mathcal{W}_{h}^{\perp }</math>, de ([[#eq-95|95]]) se tiene

<span id="eq-96"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>0 = \int _\Omega \mathbf{\tilde{u}}^T \mathbf{v}_{h}{\,d}\Omega  \approx \overset{n_{el}}{\underset{e=1}{\sum }} \int _{\Omega ^e} \left({\tau }^{e}\left[\mathbf{F}-\mathcal{L}\left(\mathbf{u}_{h}\right)\right]\right)^T \mathbf{v}_{h}{\,d}\Omega  +\overset{n_{el}}{\underset{e=1}{\sum }}\int _{\Omega ^e}\left({\tau }^{e}\mathbf{v}_{h,{ort}}\right)^T \mathbf{v}_{h}{\,d}\Omega  \qquad \forall \mathbf{v}_{h}\in \mathcal{W}_{h}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (96)
|}

Si llamamos <math display="inline">\mathbf{\Pi }_{\tau }</math> la proyección <math display="inline">L^2</math> sobre el espacio de elementos finitos pesada con <math display="inline">{\tau }^e</math> dentro de cada elemento, de ([[#eq-96|96]]) se tiene que

<span id="eq-97"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{v}_{h,{ort}}=-\mathbf{\Pi }_{\tau }\left[\mathbf{F}-\mathcal{L}\left(\mathbf{u}_{h}\right)\right]. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (97)
|}

Reemplazando ([[#eq-97|97]]) en ([[#eq-95|95]]) tenemos

<span id="eq-98"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{\tilde{u}}\approx {\tau }^{e}\left[\mathbf{F}-\mathcal{L}\left(\mathbf{u}_{h}\right)-\mathbf{\Pi }_{\tau }\left[\mathbf{F}-\mathcal{L}\left(\mathbf{u}_{h}\right)\right]\right]. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (98)
|}

Sea ahora

<span id="eq-99"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{\Pi }_{\tau }^{\perp }=\mathbf{I}-\mathbf{\mathbf{\Pi }_{\tau }}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (99)
|}

donde <math display="inline">\mathbf{I}</math> es la identidad en <math display="inline">\mathcal{W}_{h}</math> y <math display="inline">\mathbf{\Pi }_{\tau }^{\perp }</math> es la proyección sobre el espacio ortogonal al espacio de elementos finitos. La ecuación ([[#eq-98|98]]) se puede escribir entonces como

<span id="eq-100"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{\tilde{u}}\approx {\tau }^{e}\mathbf{\Pi }_{\tau }^{\perp }\left[\mathbf{F}-\mathcal{L}\left(\mathbf{u}_{h}\right)\right]. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (100)
|}

Podemos reemplazar esta expresión en ([[#eq-83|83]]) y usar el hecho de que las integrales sobre los contornos de los elementos son cero, con lo cual obtenemos la formulación estabilizada de elementos finitos OSS que resuelve la ecuación vectorial de CDR estacionaria, la cual consiste en encontrar <math display="inline">\mathbf{u}_{h}\in \mathcal{W}_{h}</math> tal que

<span id="eq-101"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>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}{\tau }^{e}\mathbf{\Pi }_{\tau }^{\perp }\left[\mathcal{L}\left(\mathbf{u}_{h}\right)-\mathbf{F}\right]{\,d}\Omega =\left\langle \mathbf{F},\mathbf{v}_{h}\right\rangle \qquad \forall \mathbf{v}_{h}\in \mathcal{W}_{h}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (101)
|}

====3.4.3 Generalización====

Las ecuaciones ([[#eq-92|92]]) y ([[#eq-101|101]]), que definen respectivamente los métodos ASGS y OSS,  se pueden escribir en forma conjunta de la siguiente manera: encontrar <math display="inline">\mathbf{u}_{h}\in \mathcal{W}_{h}</math> tal que

<span id="eq-102"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>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}{\tau }^{e}\mathcal{\mathbf{P}}\left[\mathcal{L}\left(\mathbf{u}_{h}\right)-\mathbf{F}\right]{\,d}\Omega =\left\langle \mathbf{F},\mathbf{v}_{h}\right\rangle \qquad \forall \mathbf{v}_{h}\in \mathcal{W}_{h}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (102)
|}

donde

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathcal{\mathbf{P}}= \mathbf{I}  \hbox{ para el metodo ASGS},</math>
|-
| style="text-align: center;" | <math> \mathbf{\Pi }_{\tau }^{\perp }  \hbox{ para el metodo OSS}. </math>
|}
|}

===3.5 Ecuación de CDR vectorial transitoria===

Describimos a continuación la extensión de las formulaciones introducidas al problema transitorio. Como en el caso estacionario, consideramos para simplificar condiciones de contorno de Dirichlet homogéneas.

Las únicas modificacions que debemos tener en cuenta son las originadas por el término <math display="inline">\partial _t \mathbf{u}</math> en la ecuación&nbsp;([[#eq-61|61]]). En particular, la forma variacional continua será, en vez de ([[#eq-67|67]]): hallar <math display="inline">\mathbf{u} :\gamma \rightarrow \mathcal{W}</math> tal que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>(\partial _t \mathbf{u},  \mathbf{v}) + B ( \mathbf{u},  \mathbf{v}) = \langle  \mathbf{F},  \mathbf{v}\rangle \quad \forall \mathbf{v}\in \mathcal{W}, </math>
|}
|}

y cumpliendo la condición inicial. Aquí y en lo que sigue, <math display="inline">({f},{g}) =\int _\Omega {f}^T {g} {\,d}\Omega </math> para cualesquiera funciones vectoriales <math display="inline">{f}</math> y <math display="inline">{g}</math> de <math display="inline">L^2(\Omega )^2</math>.

Como hemos indicado anteriormente, la discretización en el tiempo la llevaremos a cabo usando la regla trapezoidal. Una vez linealizadas las ecuaciones mediante el método de Picard, el problema continuo en el espacio que habrá que resolver en cada paso de tiempo, escrito en forma variacional, será: dado <math display="inline">\mathbf{u}^n\in \mathcal{W}</math>, encontrar <math display="inline">\mathbf{u}^{n+1} \in \mathcal{W}</math> tal que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>(\delta _t^n \mathbf{u},  \mathbf{v}) + B ( \mathbf{u},  \mathbf{v}) = \langle  \mathbf{F},  \mathbf{v}\rangle \quad \forall \mathbf{v}\in \mathcal{W}, </math>
|}
|}

para <math display="inline">n = 1,2,\dots </math>, siendo <math display="inline">\mathbf{u}^0</math> la condición inicial. Recordemos de la notación introducida en ([[#eq-11|11]])-([[#eq-13|13]]) que se sobreentiende en esta ecuación que la incógnita <math display="inline">\mathbf{u}</math> (y el término de fuerzas) está aproximada en el instante de tiempo <math display="inline">t^{n+\theta } = t^n + \theta \delta t</math> (y en una iteración determinada).

Al plantear la descomposición <math display="inline">\mathbf{u} = \mathbf{u}_h + \tilde{\mathbf{u}}</math>, igual que hicimos en el caso estacionario, aparecerá el término <math display="inline">\partial _t\tilde{\mathbf{u}}</math>, que en principio no es nulo. Sin embargo, en todo lo que sigue haremos la suposición de que <math display="inline">\partial _{t}\tilde{\mathbf{u}}=\mathbf{0}</math>, lo cual implica que la variación temporal de las subescalas  es despreciable comparada con el resto de términos. En <span id='citeF-22'></span>[[#cite-22|[22]]] las subescalas que satisfacen esta condición son llamadas cuasi-estáticas, mientras que se denominan dinámicas si no se hace esta hipótesis (véase <span id='citeF-22'></span><span id='citeF-36'></span><span id='citeF-37'></span>[[#cite-22|[22,36,37]]]).

====3.5.1 Formulación ASGS para el problema de CDR vectorial transitorio====

Siguiendo los mismos pasos que para llegar a la aproximación ([[#eq-92|92]]) del problema estacionario, el caso transitorio aproximado mediante el método ASGS consiste en: dado <math display="inline">\mathbf{u}^n_h\in \mathcal{W}_h</math>, encontrar <math display="inline">\mathbf{u}^{n+1}_h \in \mathcal{W}_h</math> tal que

<span id="eq-103"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\delta ^n_{t}\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}{\tau }^{e}\left[\delta ^n_{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 \quad \forall \mathbf{v}_{h}\in \mathcal{W}_{h}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (103)
|}

Nótese que no aparece el término <math display="inline">\left(\delta ^n_{t}\tilde{\mathbf{u}},\mathbf{v}_{h}\right)</math>, puesto que hemos supuesto subescalas cuasi-estacionarias, pero sí que aparece <math display="inline">\delta ^n_{t}\mathbf{u}_h</math> en el término de estabilización.

La matriz de parámetros de estabilización <math display="inline">{\tau }^{e}</math> para elementos finitos de alto orden está dada por las ecuaciones ([[#eq-90|90]]) y ([[#eq-91|91]]), el operador <math display="inline">\mathcal{L}</math> por ([[#eq-62|62]]) y el operador <math display="inline">\mathcal{L}^{*}</math> por la ecuación ([[#eq-81|81]]).

====3.5.2 Formulación OSS para el problema de CDR vectorial transitorio====

También podemos seguir los mismos pasos que para el método ASGS en el caso del método OSS. La aproximación numérica a la que llegaremos será la siguiente: dado <math display="inline">\mathbf{u}^n_h\in \mathcal{W}_h</math>, encontrar <math display="inline">\mathbf{u}^{n+1}_h \in \mathcal{W}_h</math> tal que

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\delta ^n_{t}\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}{\tau }^{e}\mathbf{\Pi }_\tau ^\perp \left[\delta ^n_{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 \quad \forall \mathbf{v}_{h}\in \mathcal{W}_{h}. </math>
|}
|}

Podemos aplicar aquí los mismos comentarios que para ([[#eq-103|103]]). Sin embargo, hay una consideración importante a tener en cuenta, y es que puesto que <math display="inline">\delta _t^n \mathbf{u}_h</math> es una función que pertenece al espacio de elementos finitos <math display="inline">\mathcal{W}_h</math>, ''su componente en el espacio ortogonal es nula'', esto es, la ecuación variacional discreta anterior se reduce a

<span id="eq-104"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left(\delta ^n_{t}\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}{\tau }^{e}\mathbf{\Pi }_\tau ^\perp \left[\mathcal{L}\left(\mathbf{u}_{h}\right)-\mathbf{F}\right]{\,d}\Omega  =\left\langle \mathbf{F},\mathbf{v}_{h}\right\rangle \quad \forall \mathbf{v}_{h}\in \mathcal{W}_{h}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (104)
|}

==4 Resultados de los experimentos numéricos==

Con el fin de evaluar las formulaciones de los métodos variacionales de subescalas ASGS y OSS con elementos finitos de alto orden aplicados a las ecuaciones del movimiento de un fluido en aguas poco profundas, hemos realizado diferentes experimentos numéricos, tales como pruebas de convergencia en malla y la aproximación de varios problemas que se encuentran en la literatura. A continuación presentamos los resultados para un test de convergencia y para dos ejemplos que se han convertido en clásicos.

El método de linealización empleado es el de Picard. El criterio de convergencia en no linealidad lo hemos fijado en la norma <math display="inline">L^2</math>, de manera que el proceso iterativo lo hemos considerado convergido cuando

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \frac{\left\Vert \mathbf{u}_{h}^{i}-\mathbf{u}_{h}^{i-1}\right\Vert _{L^{2}}}{\left\Vert \mathbf{u}_{h}^{i}\right\Vert _{L^{2}}} < {tol}, </math>
|}
|}

donde <math display="inline">i</math> es el contador de iteraciones en un paso de tiempo dado y tol es la tolerancia de convergencia, fijada en <math display="inline">10^{-5}</math> en todos los ejemplos. Aquí y en lo que sigue, las normas <math display="inline">L^2</math> se extienden a todo el dominio <math display="inline">\Omega </math>. Asimismo, todas las unidades se consideran dadas en el SI.

===4.1 Pruebas de convergencia en malla===

Las pruebas de convergencia en malla que hemos llevado a cabo consisten en calcular el error en norma <math display="inline">L^{2}</math> entre el valor exacto y la aproximación numérica de una incógnita <math display="inline">u</math> (cualquiera de las del problema). Si el error del método se comporta como el error de interpolación, situación que puede considerarse óptima, este error debe comportarse como

<span id="eq-105"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>e=\left\Vert u-u_{h}\right\Vert _{L^{2}}\leqslant ch^{d+1}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (105)
|}

donde <math display="inline">u_{h}</math> es la aproximación de <math display="inline">u</math> en el espacio de elementos finitos, <math display="inline">c</math> es una constante positiva, <math display="inline">h</math> es aquí el diámetro de la partición de elementos finitos (no la altura de la lámina de agua) y <math display="inline">d</math> el grado polinómico de la función de forma correspondiente a la interpolación de la aproximación <math display="inline">u_{h}</math>. Aunque consideraremos valores de hasta <math display="inline">d = 4</math>, no consideraremos convergencia en <math display="inline">d</math>, por lo que la constante <math display="inline">c</math> la tomaremos fija (sabemos que decrece con <math display="inline">d</math>).

Si graficamos la ecuación ([[#eq-105|105]]) (en el caso de igualdad) en el plano <math display="inline">\mathbb{R}^{2}</math> con <math display="inline">\log e</math> en las ordenadas y <math display="inline">\log h</math> en las abscisas, tenemos una línea recta, donde <math display="inline">d+1</math> es la pendiente teórica óptima de convergencia. Las pruebas de convergencia en malla consisten entonces en verificar que la pendiente calculada sea precisamente <math display="inline">d+1</math>. Para ello, seleccionamos como solución exacta tanto para las velocidades <math display="inline">U_1</math>, <math display="inline">U_2</math> como para la elevación de la superficie libre del agua <math display="inline">\eta </math> la función polinómica espacio-temporal

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\eta \left(x,y,t\right)= U_1\left(x,y,t\right)= U_2\left(x,y,t\right)=x^{6}y^{6}\left(1-x\right)^{6}\left(1-y\right)^{6}t,  </math>
|}
|}

e imponemos el término fuente <math display="inline">\mathbf{F}</math> en la ecuación vectorial de CDR transitoria ([[#eq-5|5]]) mediante las siguientes expresiones tomadas de las ecuaciones ([[#eq-1|1]]), ([[#eq-2|2]]):

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f_{i}  =\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}\qquad i=1,2,</math>
|-
| style="text-align: center;" | <math> f_{3}  =\partial _{t}h+\partial _{i}\left(hU_{i}\right), </math>
|}
|}

siendo <math display="inline"> \mathbf{F}= [ f_{1},  f_{2},  f_{3}]^T</math> y <math display="inline">\bar{\tau }_{ji}</math> y <math display="inline">Q_{i}</math> dados por las ecuaciones ([[#eq-3|3]]) y ([[#eq-4|4]]), respectivamente. Hemos tomado <math display="inline">H=1</math> m, es decir, batimetría constante con profundidad de 1&nbsp;m, aceleración de la gravedad <math display="inline">g=10</math> m/s<math display="inline">^{2}</math> y viscosidad cinemática <math display="inline">\nu _{H}=\frac{\mu _{H}}{\rho }=10^{-3}</math> m<math display="inline">^{2}</math>/s. El parámetro de Coriolis <math display="inline">\hat{f}</math>, las tensiones de superficie libre debidas al viento <math display="inline">\tau _{3i}^{s}</math>, las variaciones de presión atmosférica <math display="inline">p_{a}</math> y el coeficiente de fricción en el fondo <math display="inline">{\hbox{Ch}}</math> no han sido tomados en cuenta para los ejemplos de las pruebas de convergencia en malla.

A continuación presentamos 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 <math display="inline">P_{1},\; P_{2},\; P_{3},\; P_{4}</math> y cuadrangulares <math display="inline">Q_{1},\; Q_{2},\; Q_{3},\; Q_{4}</math>, respectivamente.

El dominio computacional es el cuadrado <math display="inline">\left[0,1\right]\times \left[0,1\right]</math> y el intervalo de tiempo es <math display="inline">\left[0,1\right]</math>, con incrementos para cada paso de tiempo <math display="inline">\delta t= 0.2</math>. 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 <math display="inline">P_{1},\; P_{2},\; P_{3},\; P_{4}</math> o <math display="inline">Q_{1},\; Q_{2},\; Q_{3},\; Q_{4},</math> hemos calculado los errores con norma <math display="inline">L^{2}</math> para ocho mallas de tamaño <math display="inline">h</math>, siendo <math display="inline">\frac{1}{h}=15,\;,\;,\;,\;,\;,\;,\;,</math> el número de partes en las que se ha dividido cada lado del cuadrado. En el Cuadro&nbsp;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. El integrador temporal utilizado es BDF1.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |Table. 1 Refinamiento para elementos triangulares y cuadradangulares
|- style="border-top: 2px solid;"
| colspan='2' | Triángulos
|- style="border-top: 2px solid;"
| colspan='4' | Número de nodos
| elem. 
| <math>P_1</math>
| <math>P_2</math>
| <math>P_3</math>
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| <math>P_4</math>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 
|- style="border-top: 2px solid;"
| 40401  
|}
&nbsp;
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|- style="border-top: 2px solid;"
| colspan='2' | Cuadrados
|- style="border-top: 2px solid;"
| colspan='4' | Número de nodos
| elem. 
| <math>Q_1</math>
| <math>Q_2</math>
| <math>Q_3</math>
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| <math>Q_4</math>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 
|- style="border-top: 2px solid;"
| 40401  
|}

En el encabezado de las gráficas de las pruebas de convergencia en malla 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 <math display="inline">\left(\frac{1}{h}=15,\;,\;,\;,\;,\;,\;,\;\right)</math>. 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 <math display="inline">m(P1),</math> <math display="inline">m(P2),</math> <math display="inline">m(P3),</math> <math display="inline">m(P4),</math> 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 <math display="inline">m(Q1),</math> <math display="inline">m(Q2),</math> <math display="inline">m(Q3),</math> <math display="inline">m(Q4)</math>. Los valores de las pendientes teóricas para las variables de velocidad se corresponden con: <math display="inline">m(P1)=m(Q1)=2</math>, <math display="inline">m(P2)=m(Q2)=3</math>, <math display="inline">m(P3)=m(Q3)=4</math>, <math display="inline">m(P4)=m(Q4)=5</math>, y las pendientes teóricas para la variable <math display="inline">\eta </math> se corresponden con <math display="inline">m(P1)=m(Q1)=1</math>, <math display="inline">m(P2)=m(Q2)=2</math>, <math display="inline">m(P3)=m(Q3)=3</math>, <math display="inline">m(P4)=m(Q4)=4</math>.

En la Figura&nbsp;[[#img-1|1]] se muestra la convergencia en malla para la primera componente de la velocidad usando elementos triangulares. Las otras variables muestran igualmente convergencia óptima para todos los elementos considerados. Lo mismo sucede en el caso de elementos cuadrangulares.

Para los dos métodos de estabilización considerados, ASGS y OSS, los valores de las constantes algorítmicas <math display="inline">c_{i},</math> <math display="inline">i=1,2,3,4,</math> se han calibrado a <math display="inline">c_{1}=12,</math> <math display="inline">c_{2}=2,</math> <math display="inline">c_{3}=1,</math> <math display="inline">c_{4}=1</math>, como ya se ha mencionado. Para la convergencia del método OSS con elementos cuadráticos cuadrados <math display="inline">Q_2</math> fue necesario utilizar <math display="inline">c_1=15.</math> Aunque para los elementos lineales <math display="inline">Q_1</math> y de alto orden <math display="inline">Q_3,</math> <math display="inline">Q_4,</math> es suficiente utilizar <math display="inline">c_1 = 4,</math> hemos preferido usar <math display="inline">c_1=15</math> para todos los elementos con el objeto de uniformizar los datos en las gráficas de las pruebas de convergencia en malla. En las figuras de las gráficas de convergencia en malla se observa tanto gráfica como numéricamente que las pendientes calculadas son mayores o tienden al valor de la pendiente teórica <math display="inline">d+1</math>.

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Codina_771723027-CC_t_asgs_u1.png|213px|]]
|[[Image:draft_Codina_771723027-CC_t_oss_u1.png|213px|Convergencia ASGS-OSS, elementos triangulares, velocidad U₁.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 1:''' Convergencia ASGS-OSS, elementos triangulares, velocidad <math>U_{1}</math>.
|}

Los resultados de las gráficas de convergencia muestran claramente que las formulaciones generalizadas para los métodos estabilizados ASGS y OSS son capaces de aproximar correctamente el problema de CDR vectorial transitorio con convección dominante que estamos considerando, con igual interpolación para todas las variables e incluyendo no linealidad en los términos convectivo y de reacción.

===4.2 Flujo a través de un obstáculo elíptico===

Este caso de referencia ha sido ampliamente examinado en la literatura utilizando diferentes métodos, como volúmenes finitos <span id='citeF-6'></span><span id='citeF-7'></span>[[#cite-6|[6,7]]], método de Godnov <span id='citeF-14'></span>[[#cite-14|[14]]], método de los mínimos cuadrados con elementos finitos <span id='citeF-10'></span><span id='citeF-13'></span>[[#cite-10|[10,13]]], o esquemas WENO en diferencias finitas <span id='citeF-5'></span>[[#cite-5|[5]]].  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.

El dominio computacional es <math display="inline">\left[0,\right]\times \left[0,\right]</math> y la batimetría mostrada en la Figura&nbsp;[[#img-2|2]] está dada por

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> z\left(x,y\right)=0.8 \,{\exp }(-5(x-0.9)^{2}-50\left(y-0.5\right)^{2}). </math>
|}
|}

El movimiento se inicia con la condición inicial de la superficie libre del agua perturbada por un desplazamiento hacia arriba de <math display="inline">0.01</math>&nbsp;m  en la región <math display="inline">0.05\leq x\leq{0.15}</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \eta (x,y,0)=\left\{\begin{array}{cc} 0.01\; \hbox{m} & \qquad \hbox{ si}\quad{0.05}\leq x\leq{0.15},\\ 0\; \hbox{ m} & \qquad \hbox{ de otra manera}. \end{array}\right. </math>
|}
|}

El momento inicial en las direcciones <math display="inline">x,</math> <math display="inline">y</math> es nulo:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> U_{1}(x,y,0)=U_{2}(x,y,0)=0\; \hbox{m/s}. </math>
|}
|}

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 <math display="inline">Q_1,</math> <math display="inline">Q_2,</math> <math display="inline">Q_3,</math> y <math display="inline">Q_4</math> que den el mismo o aproximadamente el mismo número de nodos entre sí. Tomamos como referencia el refinamiento de malla para <math display="inline">Q_2</math>, que consta de <math display="inline">200\times{100}</math> elementos cuadrados uniformes de <math display="inline">9</math> nodos, es decir, el espaciado en las direcciones <math display="inline">x</math> e <math display="inline">y</math> es <math display="inline">\delta x=\delta y=0.01</math>&nbsp;m, con un total de <math display="inline">80601</math> nodos. De esta manera, la malla de referencia para elementos <math display="inline">Q_1</math> consta de <math display="inline">400\times{200}</math> elementos cuadrados uniformes de 4 nodos, es decir <math display="inline">\delta x=\delta y=0.005</math>&nbsp;m, con un total de <math display="inline">80601</math> nodos. La malla de referencia para elementos <math display="inline">Q_3</math> consta de <math display="inline">133\times{67}</math> elementos cuadrados uniformes de 16 nodos, es decir <math display="inline">\delta x=\delta y=0.015</math> m,  con un total de <math display="inline">80800</math> nodos. Finalmente, la malla de referencia para elementos <math display="inline">Q_4</math> consta de <math display="inline">100\times{50}</math> elementos cuadrados uniformes de 25 nodos, es decir <math display="inline">\delta x=\delta y=0.02</math>&nbsp;m, con un total de <math display="inline">80601</math> nodos.

El tamaño del paso de tiempo es uniforme, con <math display="inline">\delta t=0.001</math>&nbsp;s, y el integrador temporal es BDF1 en todos los casos.

La referencia del nivel del agua en reposo, es decir, con <math display="inline">\eta=0</math>&nbsp;m, es <math display="inline">z=H=1.0</math>&nbsp;m.

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Codina_771723027-OE_Batimetria.png|357px|Flujo a través de un obstáculo elíptico. Ilustración del dominio computacional y de la batimetría. ]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' Flujo a través de un obstáculo elíptico. Ilustración del dominio computacional y de la batimetría. 
|}

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Codina_771723027-OE_q3_2D_12.png|632px|]]
|[[Image:draft_Codina_771723027-OE_q4_2D_12.png|400px|]]
|-
|[[Image:draft_Codina_771723027-OE_q3_2D_24.png|632px|]]
|[[Image:draft_Codina_771723027-OE_q4_2D_24.png|400px|]]
|-
|[[Image:draft_Codina_771723027-OE_q3_2D_36.png|632px|]]
|[[Image:draft_Codina_771723027-OE_q4_2D_36.png|400px|]]
|-
|[[Image:draft_Codina_771723027-OE_q3_2D_48.png|632px|]]
|[[Image:draft_Codina_771723027-OE_q4_2D_48.png|400px|]]
|-
|[[Image:draft_Codina_771723027-OE_q3_2D_60.png|632px|]]
|[[Image:draft_Codina_771723027-OE_q4_2D_60.png|400px|Flujo a través de un obstáculo elíptico. Contornos 2D de la superficie libre del agua η.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 3:''' Flujo a través de un obstáculo elíptico. Contornos 2D de la superficie libre del agua <math>\eta{.}</math>
|}

La Figura&nbsp;[[#img-3|3]] muestra los contornos 2D de la superficie libre del agua <math display="inline">\eta </math> para distintos instantes de tiempo <math display="inline">t= 0.12, 0.24, 0.36, 0.48, 0.60</math> s. Los contornos 2D calculados se presentan para los elementos <math display="inline">Q_3</math> y <math display="inline">Q_4</math>, con similares resultados para elementos <math display="inline">Q_1</math> y <math display="inline">Q_2</math>. A pesar de ser calculados para todos los elementos con mallas con igual número total de nodos, los resultados conforme se aumenta el orden de los polinomios de interpolación presentan mayor detalle que los de menor orden; esto se aprecia mirando los gráficos horizontalmente para un instante de tiempo, comparando los resultados de elementos <math display="inline">Q_3</math> y <math display="inline">Q_4</math>.

<div id='img-4a'></div>
<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Codina_771723027-OE_elevacion_q2_0.png|600px|]]
|[[Image:draft_Codina_771723027-OE_3D_elevacion_0.png|400px|]]
|-
|[[Image:draft_Codina_771723027-OE_elevacion_q2_0.png|600px|]]
|[[Image:draft_Codina_771723027-OE_3D_elevacion_0.png|400px|]]
|-
|[[Image:draft_Codina_771723027-OE_elevacion_q2_0.png|600px|]]
|[[Image:draft_Codina_771723027-OE_3D_elevacion_0.png|400px|]]
|-
|[[Image:draft_Codina_771723027-OE_elevacion_q2_0.png|600px|]]
|[[Image:draft_Codina_771723027-OE_3D_elevacion_0.png|400px|]]
|-
|[[Image:draft_Codina_771723027-OE_elevacion_q2_0.png|600px|]]
|[[Image:draft_Codina_771723027-OE_3D_elevacion_0.png|400px|(a)]]
| (a) 
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 4:''' Flujo a través de un obstáculo elíptico. Contornos de la superficie libre del agua <math>\eta ,</math> con elementos <math>Q_2</math>. (a) Resultados del presente estudio, (b) resultados de Shin-Jye Liang y Ying-Chih Chen <span id='citeF-10'></span>[[#cite-10|[10]]].
|}

En la Figura&nbsp;[[#img-4|4]] mostramos una comparación de los resultados de los contornos 3D de la superficie libre del agua <math display="inline">\eta </math>. Los gráficos de la parte izquierda corresponden a los resultados de nuestro trabajo y los gráficos de la parte derecha son los resultados de Shin-Jye Liang y Ying-Chih Chen&nbsp;<span id='citeF-10'></span>[[#cite-10|[10]]], ambos estudios con iguales características de mallado y paso de tiempo (elementos <math display="inline">Q_2,</math> <math display="inline">\delta x=\delta y=0.01</math>&nbsp;m y <math display="inline">\delta t=0.001</math>&nbsp;s). Con respecto a los otros resultados 3D para elementos <math display="inline">Q_1,</math> <math display="inline">Q_3,</math> y <math display="inline">Q_4</math> de nuestro trabajo, las diferencias entre sí y con respecto al elemento <math display="inline">Q_2</math> son visualmente imperceptibles, por lo que en la Figura&nbsp;[[#img-4|4]](a) presentamos solo los resultados correspondientes al elemento <math display="inline">Q_2</math>, para poder así compararlos con los resultados de&nbsp;<span id='citeF-10'></span>[[#cite-10|[10]]]. En ambos trabajos el comportamiento de la onda es bastante similar, aunque en nuestro caso la evolución temporal es más disipativa por haber utilizado un integrador temporal de primer orden (recuérdese que nuestro propósito es el diseño de aproximaciones precisas y robustas en el espacio). Así, la perturbación que inicialmente es plana y se propaga hacia la derecha se ve afectada por el obstáculo elíptico en la parte inferior del lecho. El incremento de la amplitud de onda debido a la disminución de la profundidad del agua es obvio en <math display="inline">t=0.24</math>&nbsp;s, debido a la disminución de la velocidad de onda, siendo la más baja de todo el flujo en el momento en el cual la onda atraviesa el punto más alto del obstáculo elíptico y que aproximadamente se produce en el instante <math display="inline">t=0.24</math>&nbsp;s. Se observan ondas estacionarias de reflexión justo detrás del pico producido por el choque de la onda con el obstáculo elíptico, en <math display="inline">t= 0.24, 0.36, 0.48</math> y <math display="inline">0.60</math> s.


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-2'></span>Table. 2 Flujo a través de un obstáculo elíptico. Picos de elevación máxima
|- style="border-top: 2px solid;"
|   
| colspan='4' | Elevación de la superficie
|-
| colspan='4' | libre del agua, <math>\eta </math> (mm)
| <math>Q_1</math>
| <math>Q_2</math>
| <math>Q_3</math>
|- style="border-top: 2px solid;border-bottom: 2px solid;"
| <math>Q_4</math>0.12 
| 5.0975 
| 5.0939 
| 5.0716 
| 4.9786 0.24 
| 7.2415 
| 7.2120 
| 7.1718 
| 6.9990 0.36 
| 5.9430 
| 5.8503 
| 5.7924 
| 5.5978 0.48  
| 5.6073 
| 5.6526 
| 5.6161 
| 5.1742 0.60 
| 3.9634 
| 3.9812 
| 3.9523  
|- style="border-top: 2px solid;"
| 3.7503  
|}

En el Cuadro&nbsp;[[#table-2|2]] se encuentran tabuladas las elevaciones de la superficie libre del agua <math display="inline">\eta ,</math> correspondientes a los picos más altos de la onda para los instantes de tiempo <math display="inline">t= 0.12, 0.24, 0.36, 0.48</math> y <math display="inline">0.60</math> s. Observamos como en todos los casos la onda alcanza su pico máximo en <math display="inline">t= 0.24</math> s, debido precisamente a que alrededor de <math display="inline">t= 0.24</math> s la onda atraviesa el punto más alto del obstáculo elíptico. También observamos que para cada instante de tiempo la elevación de la superficie del agua <math display="inline">\eta </math> disminuye conforme se incrementa el orden del elemento, lo que es completamente coincidente con la mejor precisión en elementos de alto orden. Como podemos ver en el Cuadro&nbsp;[[#table-2|2]], para un instante de tiempo las variaciones de <math display="inline">\eta </math> conforme se aumenta el orden del elemento son del orden de las décimas, centésimas y en algunos casos milésimas de milímetro, por lo que dichas variaciones son visualmente imperceptibles, como se ha mencionado anteriormente, motivo por el cual en la Figura&nbsp;[[#img-4|4]](a) se ha presentado solo el contorno 3D correspondiente a elementos <math display="inline">Q_2</math>.

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Codina_771723027-OE_Contornos_2D.png|630px|Flujo a través de un obstáculo elíptico. Contornos 2D de la superficie libre del agua de: (a) resultados de Shin-Jye y Ying-Chih Chen <span id='citeF-10'></span>[[#cite-10|[10]]] con elementos Q₂ cuadrado de 9 nodos, (b) resultados de Ying-Chih Chen <span id='citeF-13'></span>[[#cite-13|[13]]] con elementos P₁ triángulos de 3 nodos, (c) resultados de Akoh et al.&nbsp;<span id='citeF-6'></span>[[#cite-6|[6]]] con elementos P₁.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 5:''' Flujo a través de un obstáculo elíptico. Contornos 2D de la superficie libre del agua de: (a) resultados de Shin-Jye y Ying-Chih Chen <span id='citeF-10'></span>[[#cite-10|[10]]] con elementos <math>Q_2</math> cuadrado de 9 nodos, (b) resultados de Ying-Chih Chen <span id='citeF-13'></span>[[#cite-13|[13]]] con elementos <math>P_1</math> triángulos de 3 nodos, (c) resultados de Akoh et al.&nbsp;<span id='citeF-6'></span>[[#cite-6|[6]]] con elementos <math>P_1</math>.
|}

Finalmente, en la Figura&nbsp;[[#img-5|5]] presentamos los resultados de los contornos 2D del nivel de superficie libre del agua <math display="inline">\eta </math> de: (a) Shin-Jye y Ying-Chih Chen&nbsp;<span id='citeF-10'></span>[[#cite-10|[10]]] (con elementos <math display="inline">Q_2,</math> <math display="inline">\delta x=\delta y=0.01</math> m y <math display="inline">\delta t=0.001</math> s, (b) Lian y Hsu&nbsp;<span id='citeF-13'></span>[[#cite-13|[13]]] (con elementos <math display="inline">P_1</math>, <math display="inline">\delta x=\delta y=0.01</math> m y <math display="inline">\delta t=0.0005</math> s) y (c) Akoh et al.&nbsp; <span id='citeF-6'></span>[[#cite-6|[6]]] (con elementos <math display="inline">P_1</math>, <math display="inline">\delta x=\delta y=0.01</math> m y <math display="inline">\delta t=0.0005</math> s), para ser comparados con nuestros resultados de la Figura&nbsp;[[#img-3|3]]. En general observamos que los resultados del presente trabajo dan estructuras de ondas y gradientes más nítidos y con mas detalle que los resultados de&nbsp;<span id='citeF-10'></span><span id='citeF-13'></span>[[#cite-10|[10,13]]] y&nbsp;<span id='citeF-6'></span>[[#cite-6|[6]]].

===4.3 Flujo de la rotura de una presa a través de una compuerta de esclusa===

A continuación consideramos el problema de la rotura de una presa, el cual se puede encontrar por ejemplo en&nbsp;<span id='citeF-38'></span><span id='citeF-39'></span>[[#cite-38|[38,39]]]. La ilustración del problema se muestra en la Figura&nbsp;[[#img-6|6]](a). La geometría del dominio computacional es una región cuadrada de 200 m <math display="inline">\times </math> 200 m con fondo plano sin fricción, y deslizamiento libre en todas las paredes del embalse (aguas arriba) y del desembalse (aguas abajo), excepto en las paredes interiores de separación del embalse y desembalse y de la compuerta de la esclusa. Dicha compuerta de esclusa de 75 m de ancho está localizada en la mitad del dominio. El nivel inicial del agua en el embalse de la presa (aguas arriba) es de 10 m, y el del desembalse (aguas abajo) es de 5 m de alto con el agua en reposo. En el instante en que falla la presa, el embalse de agua se libera a través de la compuerta de esclusa, creando ondas de propagación aguas abajo (desembalse). En cuanto a las condiciones de contorno, las velocidades en la dirección perpendicular a las paredes de la presa, tanto en el embalse como en el desembalse (partes derecha e izquierda correspondientemente de la Figura&nbsp;[[#img-6|6]](a)), es necesario prescribirlas a cero; además la velocidad tangencial en las paredes interiores de separación del embalse y desembalse, así como también en las paredes de la compuerta de la esclusa, también se debe prescribir a cero. El borde izquierdo del desembalse es libre, por lo tanto la velocidad perpendicular a este borde la dejamos sin prescribir a cero. La superficie libre del agua <math display="inline">\eta </math> no se prescribe, se deja libre en todo el contorno.

Debido a que la velocidad del flujo en la dirección <math display="inline">x</math> es muy variable, el comportamiento de los parámetros de estabilización también es muy variable, con valores  en nada comparables con el paso de tiempo, lo cual genera problemas de convergencia. Para evitar estos problemas, hemos redefinido <math display="inline">\tau _1</math> en cada elemento <math display="inline">e</math> como

<span id="eq-106"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tau ^e_{1} \gets \min (\tau ^e_{1},\alpha \delta t), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (106)
|}

donde <math display="inline">\delta t</math> es el paso de tiempo y <math display="inline">\alpha </math> es un parámetro que es necesario ajustar dependiendo del tamaño de la malla y del grado del polinomio de las funciones de forma. Para este caso estudiado hemos adoptado <math display="inline">\alpha=10</math> para elementos <math display="inline">P_1</math> y <math display="inline">P_2</math>, para elementos <math display="inline">P_3</math> se ha usado <math display="inline">\alpha=1</math> y para elementos <math display="inline">P_4</math> se ha tomado <math display="inline">\alpha=0.5</math>. El valor de <math display="inline">\tau ^e_{2}</math> se calcula con ([[#eq-91|91]]), pero en este caso con el <math display="inline">\tau ^e_{1}</math> calculado en ([[#eq-106|106]]).

<div id='img-6a'></div>
<div id='img-6b'></div>
<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Codina_771723027-RP_geometria_presa.png|205px|(a)]]
|[[Image:draft_Codina_771723027-RP_malla_3733.png|400px|(b)]]
|- style="text-align: center; font-size: 75%;"
| (a) 
| (b) 
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 6:''' Flujo de la rotura de una presa. (a) Geometría y condiciones iniciales del problema. (b) Malla con elementos triangulares <math>P_1</math> con 3733 nodos.
|}

En la Figura&nbsp;[[#img-6|6]](b) se encuentra la malla usada con 7164 elementos triangulares <math display="inline">P_1</math>, que dan un total de 3733 nodos.

<div id='img-7a'></div>
<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Codina_771723027-RP_7.png|600px|]]
|[[Image:draft_Codina_771723027-RP_7.png|400px|]]
|-
|[[Image:draft_Codina_771723027-RP_2D_h_7.png|600px|]]
|[[Image:draft_Codina_771723027-RP_corte_h_7.png|400px|Contorno 2D]]
| (a) Contorno 2D
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 7:''' Flujo de la rotura de una presa. Contornos 2D y corte longitudinal de la superficie libre del agua en 7.2 s. Comparación de: (a) Estudio presente, (b) ELRBFCM&nbsp;<span id='citeF-38'></span>[[#cite-38|[38]]].
|}

<div id='img-8a'></div>
<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Codina_771723027-RP_7.png|600px|]]
|[[Image:draft_Codina_771723027-RP_7.png|400px|]]
|-
|[[Image:draft_Codina_771723027-RP_2D_U1_7.png|600px|]]
|[[Image:draft_Codina_771723027-RP_corte_U1_7.png|400px|Contorno 2D]]
| (a) Contorno 2D
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 8:''' Flujo de la rotura de una presa. Contornos 2D y corte longitudinal de la velocidad en la dirección <math>x</math> en 7.2 s. Comparación de: (a) Estudio presente, (b) ELRBFCM&nbsp;<span id='citeF-38'></span>[[#cite-38|[38]]].
|}

<div id='img-9a'></div>
<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Codina_771723027-RP_3D_3.png|600px|]]
|[[Image:draft_Codina_771723027-RP_3D_h_3.png|400px|]]
|-
|[[Image:draft_Codina_771723027-RP_3D_RP_7.png|600px|]]
|[[Image:draft_Codina_771723027-RP_3D_h_7.png|400px|(a)]]
| (a) 
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 9:''' Flujo de la rotura de una presa. Contornos 3D de la superficie libre del agua, para los instantes de tiempo <math>t=3.5</math> s  y <math>t=7.2</math> s. Comparación de: (a) Estudio presente, (b) ELRBFCM&nbsp;<span id='citeF-38'></span>[[#cite-38|[38]]].
|}

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Codina_771723027-RP_p1p2p3p4h.png|333px|Flujo de la rotura de una presa. Corte longitudinal de la superficie libre del agua en t= 4.5 s. Comparación del refinamiento para elementos P₁, P₂, P₃ y P₄.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' Flujo de la rotura de una presa. Corte longitudinal de la superficie libre del agua en <math>t=</math> 4.5 s. Comparación del refinamiento para elementos <math>P_1</math>, <math>P_2</math>, <math>P_3</math> y <math>P_4</math>.
|}

En las Figuras&nbsp;[[#img-7|7]] y [[#img-8|8]] se presentan los contornos 2D y los cortes longitudinales indicados de la superficie libre del agua y de la velocidad en la dirección <math display="inline">x</math> en <math display="inline">t=</math>7.2 s, respectivamente. En las Figuras&nbsp;[[#img-7|7]](a) y [[#img-8|8]](a) se encuentran los resultados de nuestro trabajo, mientras que en las Figuras&nbsp;[[#img-7|7]](b) y [[#img-8|8]](b) están los resultados de C.K. Chou et al. con el método ELRBFCM&nbsp;<span id='citeF-38'></span>[[#cite-38|[38]]]. Observando los cortes longitudinales de las Figuras&nbsp;[[#img-7|7]] y [[#img-8|8]] vemos que nuestros resultados coinciden con los de C.K. Chou en el perfil de 25209 nodos, pero en nuestro trabajo hemos usado una malla con elementos triangulares <math display="inline">P_1</math> con 3733 nodos, lo que incide en un menor costo computacional. También los contornos 2D, tanto de la superficie libre del agua como de la velocidad en la dirección <math display="inline">x</math> en las Figuras&nbsp;[[#img-7|7]] y [[#img-8|8]], presentan mejores y más nítidos detalles que los resultados de&nbsp;<span id='citeF-38'></span>[[#cite-38|[38]]].

También presentamos la simulación del proceso de rotura de la presa en <math display="inline">t=3.5</math> s y <math display="inline">t=7.2</math> s, que se muestra en la Figura&nbsp;[[#img-9|9]](a) con 3733 nodos en el caso de nuestro trabajo y con 14254 nodos para los resultados de C.K. Chou et al.&nbsp;<span id='citeF-38'></span>[[#cite-38|[38]]] en la Figura&nbsp;[[#img-9|9]](b).

Finalmente presentamos una comparación de los resultados de elementos lineales <math display="inline">P_1</math> con elementos de alto orden <math display="inline">P_2</math>, <math display="inline">P_3</math> y <math display="inline">P_4</math> de los cortes longitudinales de la superficie libre del agua en el instante de tiempo <math display="inline">t=4.5</math> s en la Figura&nbsp;[[#img-10|10]]. Tomando como base la malla de 7164 elementos triangulares <math display="inline">P_1</math> con 3733 nodos y aumentado el orden de interpolación con el mismo número de elementos, obtenemos mallas con elementos triangulares <math display="inline">P_2</math> con un total de 14629 nodos, <math display="inline">P_3</math> con 32689 nodos y <math display="inline">P_4</math> con 57913 nodos. Observando las gráficas se concluye que al aumentar el orden de interpolación se obtienen mejores resultados tanto de la superficie libre del agua como de la velocidad en la dirección&nbsp;<math display="inline">x</math>, como era de esperar.

==5 Conclusiones==

En este artículo se ha presentado una revisión de los métodos estabilizados ASGS y OSS para una ecuación vectorial de tipo CDR general, estacionaria y transitoria, con aplicación a las ecuaciones del movimiento de un fluido en aguas poco profundas. Asimismo, se ha estudiado la utilización de elementos finitos de alto orden, cuadráticos, cúbicos y de cuarto orden, con elementos tanto triangulares como cuadrangulares.

Las formulaciones variacionales estabilizadas ASGS y OSS nos permiten utilizar igual interpolación en todas las variables, que en nuestro caso son tres, las dos componentes de la velocidad <math display="inline">U_{1},</math> <math display="inline">U_{2}</math> y la elevación <math display="inline">\eta </math> de la superficie libre del agua. Nos permiten también poder tratar con flujos de convección dominante. Estas dos dificultades se unen a la complejidad adicional de la no linealidad tanto en los términos convectivos como en la derivada temporal de la elevación de la superficie libre del agua. Para poder utilizar elementos finitos de alto orden fue necesario cambiar el cálculo del parámetro de estabilización, incluyendo en este cálculo el grado del polinomio.

De los resultados del test de convergencia realizados se puede concluir que puede ser ventajoso utilizar elementos de alto orden para el problema que hemos considerado. En términos generales, los errores para un número de nodos dado son menores cuanto mayor es el orden del polinomio utilizado. También se observó con experimentación numérica que la inclusión del término <math display="inline">\tau _{2}</math> en la matriz <math display="inline">{{\tau }}</math> es totalmente necesaria para el ejemplo de la rotura de la presa, ya que si se considera <math display="inline">\tau _{2}=0</math> el resultado es que con elementos <math display="inline">P_1</math> y <math display="inline">P_2</math> se presentan problemas de estabilidad alrededor de los 7 s, y con elementos <math display="inline">P_3</math> y <math display="inline">P_4</math> estos problemas se presentan desde <math display="inline">t=0.1</math> s. Sin embargo, para las pruebas de convergencia en malla y el ejemplo del flujo a través de un obstáculo elíptico, es indiferente el usar <math display="inline">\tau _{2}=0</math> o <math display="inline">\tau _{2}=\frac{\left(h/d\right)^{2}}{c_{1}\tau _{1}}</math>,  ya que las variaciones de los resultados en estos casos son mínimas.

En general, el comportamiento de los métodos estabilizado ASGS y OSS para el problema que hemos estudiado es completamente satisfactorio, mostrando que la comparación de nuestros resultados con los resultados de otros autores es favorable en todos los casos a nuestras formulaciones.

==Agradecimientos==

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

==Referencias==

===BIBLIOGRAPHY===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]''' P. Brufau and P. García-Navarro "Two-dimensional dam break flow simulation". IJNMF

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]''' J. C. Hsiek and T.Y. Yang. (2003) "Investigation on the suitability of two dimensional depth-averaged models for bend-flow simulation", Volume 129. ASCE Journal of Hydraulic Engineering 433-465

<div id="cite-3"></div>
'''[[#citeF-3|[3]]]''' E. Pena. (2002) "Estudio numérico y experimental del transporte de sedimentos en cauces aluviales"

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]''' E. Zahibo and N. Pelinovsky. (2006) "Analytical and numerical study of nonlinear effects at tsunami modeling", Volume 174. Applied Mathematics and Computation 795-809

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]''' Y. Xing and C.W. Shu "High order finite difference WENO schemes for a class of hyperbolic systems with source terms". JCP

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]''' R. Akoh and S. Li and F. Xiao "A CIP/multi-moment finite volume method for shallow water equations with source terms". IJNMF

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]''' Y. Xing  and C.W. Shu "High order well-balanced finite volume WENO schemes and discontinuos Galerkin methods for a class of hyperbolic systems with source terms". JCP

<div id="cite-8"></div>
'''[[#citeF-8|[8]]]''' V. Alessandro and V. Caleffi, V. and A. Zanni. (2002) "Case study: Malpasset dam-break simulation using a two-dimensional finite volume method", Volume 128. ASCE Journal of Hydraulic Engineering 460-472

<div id="cite-9"></div>
'''[[#citeF-9|[9]]]''' E. Hanert and V. Legat and E. Deleersnijder. (2002) "A comparison of three finite elements to solve the linear shallow water equations", Volume 5. Ocean Modelling 17&#8211;35

<div id="cite-10"></div>
'''[[#citeF-10|[10]]]''' S.-J. Liang and Y.-C. Chen. (2011) "Space-time least-squares finite-element method for shallow-water equations", Volume 19. Journal of Marine Science and Technology 571&#8211;578

<div id="cite-11"></div>
'''[[#citeF-11|[11]]]''' S. Takase and K. Kashiyama and S. Tanaka and Tayfun E. Tezduyar "Space-time SUPG formulation of the shallow-water equations". IJNMF

<div id="cite-12"></div>
'''[[#citeF-12|[12]]]''' R. Codina. (2001) "Finite element approximation of the shallow water equations using subgrid scale stabilization."

<div id="cite-13"></div>
'''[[#citeF-13|[13]]]''' S.J. Liang  and T.W. Hsu. (2009) "Least-squares finite-element method for shallow-water equations with source terms", Volume 25. Acta Mechanica Sinica 597&#8211;610

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]''' R. Le Veque "Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave-propagation algorithm". JCP

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]''' M. Fujihara and A.G.L. Borthwick. (2000) "Godunov-type solution of curvilinear shallow-water equations", Volume 126. ASCE Journal of Hydraulic Engineering 827-836

<div id="cite-16"></div>
'''[[#citeF-16|[16]]]''' J.G. Zhou "A lattice Boltzmann model for the shallow water equations". CMAME

<div id="cite-17"></div>
'''[[#citeF-17|[17]]]''' O.C. Zienkiewicz and R.L. Taylor. (2000) "The Finite Element Method. Volume 3: Fluid Mechanics", Volume 1. Butterworth-Heinemann, 5 Edition

<div id="cite-18"></div>
'''[[#citeF-18|[18]]]''' T.J.R. Hughes "Multiscale phenomena: Green's function, the  Dirichlet-to-Neumann formulation, subgrid scale models,  bubbles and the origins of stabilized formulations". CMAME

<div id="cite-19"></div>
'''[[#citeF-19|[19]]]''' T.J.R. Hughes and G.R. Feijóo and L. Mazzei and J.B. Quincy "The variational multiscale method&#8211;a paradigm for computational mechanics". CMAME

<div id="cite-20"></div>
'''[[#citeF-20|[20]]]''' R. Codina "On stabilized finite element methods for linear systems of convection-diffusion-reaction equations". CMAME

<div id="cite-21"></div>
'''[[#citeF-21|[21]]]''' R. Codina "A stabilized finite element method for generalized stationary incompressible flows". CMAME

<div id="cite-22"></div>
'''[[#citeF-22|[22]]]''' R. Codina "Stabilized finite element approximation of transient incompressible flows using orthogonal subscales". CMAME

<div id="cite-23"></div>
'''[[#citeF-23|[23]]]''' R. Codina "Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods". CMAME

<div id="cite-24"></div>
'''[[#citeF-24|[24]]]''' R. Codina and J. Blasco "A finite element formulation for the Stokes problem allowing equal velocity-pressure interpolation". CMAME

<div id="cite-25"></div>
'''[[#citeF-25|[25]]]''' R. Codina and J. Blasco. (2000) "Analysis of a pressure-stabilized finite element approximation of the stationary Navier-Stokes equations", Volume 87. Numerische Mathematik 59&#8211;81

<div id="cite-26"></div>
'''[[#citeF-26|[26]]]''' R. Codina "Comparison of some finite element methods for solving the diffusion-convection-reaction equation". CMAME

<div id="cite-27"></div>
'''[[#citeF-27|[27]]]''' R. Codina and J. Blasco. (2002) "Analysis of a stabilized finite element approximation of  the transient convection-diffusion-reaction equation using orthogonal subscales", Volume 4. Computing and Visualization in Science 167&#8211;174

<div id="cite-28"></div>
'''[[#citeF-28|[28]]]''' H. Yamamoto. (2009) "Lecture 8: The Shallow-Water Equations"

<div id="cite-29"></div>
'''[[#citeF-29|[29]]]''' J. Peraire and O.C. Zienkiewics and K. Morgan "Shallow water problems: a general explicit formulation". IJNME

<div id="cite-30"></div>
'''[[#citeF-30|[30]]]''' J. Peraire. (1986) "A finite element method for convection dominated flows"

<div id="cite-31"></div>
'''[[#citeF-31|[31]]]''' R. Codina and S. Badia and J. Baiges and J. Principe. (Por aparecer) "Variational Multiscale Methods in Computational Fluid Dynamics,  in ''Encyclopedia of Computational Mechanics''". John Wiley & Sons Ltd

<div id="cite-32"></div>
'''[[#citeF-32|[32]]]''' R. Codina and J.M. González-Ondina and G. Díaz-Hernández and J. Principe "Finite element approximation of the modified Boussinesq equations using a stabilized formulation". IJNMF

<div id="cite-33"></div>
'''[[#citeF-33|[33]]]''' C. Baiocchi and F. Brezzi and L.P. Franca "Virtual bubbles and Galerkin/least-squares type methods (Ga.L.S)". CMAME

<div id="cite-34"></div>
'''[[#citeF-34|[34]]]''' F. Brezzi and L.P. Franca and T.J.R. Hughes and A. Russo "<math>b = \int g</math>". CMAME

<div id="cite-35"></div>
'''[[#citeF-35|[35]]]''' R. Codina and J. Principe and J. Baiges "Subscales on the element boundaries in the variational two-scale finite element method". CMAME

<div id="cite-36"></div>
'''[[#citeF-36|[36]]]''' R. Codina and J. Principe and O. Guasch and S. Badia "Time dependent subscales in the stabilized finite element approximation of incompressible flow problems". CMAME

<div id="cite-37"></div>
'''[[#citeF-37|[37]]]''' R. Codina and J. Principe "Dynamic subscales in the finite element approximation of thermally coupled incompressible flows". IJNMF

<div id="cite-38"></div>
'''[[#citeF-38|[38]]]''' C.K. Chou and C.P. Sun and D.L. Young and V. Sladek. (2015) "Extrapolated local basis function collocation method for shallow water problems", Volume 50. Engineering Analysis with Boundary Elements 275&#8211;290

<div id="cite-39"></div>
'''[[#citeF-39|[39]]]''' S.J. Liang and J.H. Tang and M.S. Wu. (2008) "Solution of shallow-water equations using least-squares finite-element method", Volume 24. Acta Mechanica Sinica 523&#8211;532