Modelo hidrodinámico¶

Los modelos hidráulicos unidimensionales no son adecuados para simular inundaciones cuando los flujos no están confinados o las velocidades cambian de dirección durante el transcurso del hidrograma. El costo de los modelos numéricos tridimensionales no simplificados puede evitarse utilizando las ecuaciones bidimensionales (2D) de aguas someras promediadas en profundidad.

Al trabajar con las ecuaciones de aguas someras, las aplicaciones realistas siempre incluyen términos fuente que describen la variación del nivel del lecho y la fricción del lecho que, si no se discretizan correctamente, pueden provocar inestabilidades numéricas. En la última década, el esfuerzo principal se ha centrado en mantener un equilibrio discreto entre los flujos y los términos fuente en casos de agua en reposo, lo que ha dado lugar a la noción de esquemas bien equilibrados o propiedad C [, , , ]. Recientemente, para incluir correctamente el efecto de los términos fuente en la solución débil, se han presentado solucionadores de Riemann aproximados aumentados [Rosatti et al. (2003), ]. De este modo, se pueden calcular soluciones precisas evitando la necesidad de imponer parámetros de ajuste dependientes del caso, que se utilizan con frecuencia para evitar valores negativos de la profundidad del agua y otras inestabilidades numéricas que aparecen al incluir términos fuente.

Esta sección presenta el sistema de ecuaciones, la formulación de las condiciones de contorno y el esquema de volúmenes finitos utilizado en RiverFlow2D; la información puede ampliarse en las referencias.

Supuestos del modelo hidrodinámico¶

RiverFlow2D utiliza las ecuaciones de aguas someras resultantes de la integración vertical de la ecuación de Navier-Stokes. Por lo tanto, el modelo no calcula aceleraciones verticales ni velocidades verticales y, en consecuencia, no puede resolver flujos secundarios.
Se supone que el esfuerzo cortante del lecho sigue las direcciones de la velocidad promediada en profundidad.
El modelo no incluye términos de dispersión ni de turbulencia. La disipación de la turbulencia y las pérdidas de energía se tienen en cuenta únicamente mediante el término n de Manning en las ecuaciones de cantidad de movimiento.
El modelo puede considerar la transferencia de calor para calcular la temperatura del aceite mientras fluye sobre el terreno, y considera la variación de la densidad, la viscosidad y el esfuerzo de fluencia en el tiempo y el espacio.

Ecuaciones del flujo considerando variaciones de temperatura prescritas¶

Los flujos de aguas someras pueden describirse matemáticamente mediante ecuaciones de conservación de masa y cantidad de movimiento promediadas en profundidad, con todos los supuestos asociados. Este sistema de ecuaciones diferenciales parciales se formulará aquí en forma conservativa de la siguiente manera:

\[\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{F(U)}}{\partial x}+\frac{\partial \mathbf{G(U)}}{\partial y}=\mathbf{S}(\mathbf{U},x,y)\]

donde $\mathbf{U}=\left( h , \; q_x , \; q_y \right)^{T}$ es el vector de variables conservadas, donde $h$ representa la profundidad del agua, $q_x=uh$ y $q_y=vh$ son los caudales unitarios, y $(u,v)$ son las componentes promediadas en profundidad del vector de velocidad $\mathbf{u}$ a lo largo de las coordenadas $(x,y)$, respectivamente. Los vectores de flujo vienen dados por:

\[\mathbf{F}=\left( q_x, \; \frac {q_y^2}{h} + \frac {1}{2} g h^2 , \; \frac {q_x q_y}{h} \right)^{T}, \qquad \mathbf{G}=\left( q_y, \; \frac{q_x q_y}{h}, \; \frac {q_y^2}{h} + \frac {1}{2} g h^2 \right)^{T}\]

donde $g$ es la aceleración de la gravedad. Los términos $\frac {1}{2}gh^2$ de los flujos se han obtenido suponiendo una distribución hidrostática de la presión en cada columna de agua, como se acepta habitualmente en los modelos de aguas someras. El vector de términos fuente incorpora el efecto de la fuerza de presión sobre el lecho y de las fuerzas tangenciales generadas por el esfuerzo del lecho

\[\mathbf{S}=\left( 0 , \; g h (S_{0x}-S_{fx}) , \; g h (S_{0y}-S_{fy}) \right)^{T}\]

donde las pendientes del lecho del nivel inferior $z_b$ son

\[S_{0x} = - \frac{\partial z_b}{\partial x} , \qquad S_{0y} = - \frac{\partial z_b}{\partial y}\]

y la contribución del esfuerzo del lecho se modela mediante la ley de fricción de Manning, de modo que:

\[S_{fx}= \frac{n^2u \sqrt{u^2+v^2}}{h^{4/3}}, \qquad S_{fy}= \frac{n^2v\sqrt{u^2+v^2}}{h^{4/3}}\]

con $n$ el coeficiente de rugosidad.

Mediante esta opción, el modelo puede considerar la variación de la temperatura y su efecto sobre la densidad y la viscosidad, pero se supone que la temperatura del aceite es prescrita por el usuario y no depende de los cambios ambientales durante la simulación.

//

Solución numérica de volúmenes finitos¶

Para introducir el esquema de volúmenes finitos, se integra en un volumen o celda de la cuadrícula $\Omega$ utilizando el teorema de Gauss:

\[\frac {\partial} {\partial t} \int_{\Omega} \mathbf{U}d\Omega + \oint_{ \partial \Omega } \mathbf{E} \mathbf{n} dl = \int_{\Omega} \mathbf{S}d \Omega\]

donde $\mathbf{E}=(\mathbf{F},\mathbf{G})$ y $\mathbf{n}=(n_x,n_y)$ es el vector normal unitario exterior al volumen $\Omega$. Para obtener una solución numérica del sistema, el dominio se divide en celdas computacionales, $\Omega_i$, utilizando una malla fija. Suponiendo una representación por tramos de las variables conservadas y una formulación ascendente y unificada de los flujos y términos fuente

\[\frac {\partial} {\partial t} \int_{\Omega_i} \mathbf{U}d\Omega + \sum_{k=1} ^{NE} (\mathbf{En}- \mathbf{\bar {S }})_{ k} l_k = 0\]

La solución aproximada puede definirse utilizando una matriz jacobiana aproximada $\widetilde{\mathbf{J}}_{\mathbf{n},k}$ del flujo normal no lineal $\mathbf{E_n}$ y dos matrices aproximadas $\widetilde{\mathbf{P}} = (\mathbf{\widetilde{e}}^1, \mathbf{\widetilde{e}}^2, \mathbf{\widetilde{e}}^3 )$, y $\widetilde{\mathbf{P}}^{-1}$, construidas utilizando los eigenvectores de la jacobiana, que hacen diagonal a $\widetilde{\mathbf{J}}_{\mathbf{n},k}$

\[\widetilde{\mathbf{P}}_{k}^{-1} \widetilde{\mathbf{J}}_{\mathbf{n},k} \widetilde{\mathbf{P}}_{k} =\widetilde{\boldsymbol{\Lambda}}_{k}\]

con $\widetilde{\boldsymbol{\Lambda}}_{k}$ una matriz diagonal con eigenvalores $\widetilde{\lambda}^{m}_{k}$ en la diagonal principal

\[\widetilde{\boldsymbol{\Lambda}}_{k} = \left( \begin{array}{ccccc} \widetilde{\lambda}^{ 1} & 0 & 0 \\ 0 & \widetilde{\lambda}^{ 2} & 0 \\ 0 & 0 & \widetilde{\lambda}^{3} \\ \end{array}\right)_{k}\]

Tanto la diferencia en el vector $\mathbf{U}$ a través del borde de la cuadrícula como el término fuente se proyectan sobre la base de eigenvectores de la matriz

\[\delta \mathbf{U}_{k}= \widetilde{\mathbf{P}}_{k} \mathbf{A}_{k} \quad (\mathbf{\bar {S }})_{k}= \widetilde{\mathbf{ P}}_{k} \mathbf{B}\]

donde $\mathbf{A}_{k} = ( \alpha ^1, \alpha ^2, \alpha ^3 )^{T}_{k}$ contiene el conjunto de intensidades de onda y $\mathbf{B} = ( \beta^1 , \beta^2 , \beta^3 )^{T}_{k}$ contiene las intensidades de los términos fuente. Los detalles se proporcionan en. La linealización completa de todos los términos, combinada con la técnica ascendente, permite definir la función de flujo numérico $(\mathbf{En}- \mathbf{\bar {S }})_{ k}$ como

\[(\mathbf{En}- \mathbf{\bar {S }})_{ k} = \mathbf{E}_i \mathbf{n}_k+ \sum^{3 }_{m=1} \left( \widetilde{\lambda}^-\; \theta \alpha \widetilde{\mathbf{ e}} \right)^{m}_{ k}\]

con $\widetilde{ {\lambda}}^{-} = \frac {1} {2}( \widetilde{ {\lambda}} -\vert \widetilde{ {\lambda}} \vert )$ y $\theta^{m}_{k} = \left(1- \frac {\beta}{\widetilde{\lambda}\alpha} \right)^{m}_{k}$ que, al insertarse en, dan un método de Godunov explícito de primer orden

\[\mathbf{U}^{n+1}_{ i} = \mathbf{U}^{n}_{ i} - \sum_{k=1} ^{NE} \left[ \mathbf{E}_i \mathbf{n}_k+ \sum^{3 }_{m=1} \left( \widetilde{\lambda}^-\; \theta \alpha \widetilde{\mathbf{ e}} \right)^{m}_{ k} \right] \frac { l_k } {A_i} \Delta t\]

Como la cantidad $\mathbf{E}_i$ es uniforme por celda $i$ y se cumple la siguiente propiedad geométrica en cualquier celda

\[\sum_{k=1}^{NE} \mathbf{n}_k l_k = 0\]

puede reescribirse como

\[\mathbf{U}^{n+1}_{ i} = \mathbf{U}^{n}_{ i} - \sum_{k=1} ^{NE} \left[ \sum^{3 }_{m=1} \left( \widetilde{\lambda}^-\; \theta \alpha \widetilde{\mathbf{ e}} \right)^{m}_{ k} \right] \frac { l_k \Delta t } {A_i}\]

El método de volúmenes finitos puede escribirse utilizando una formulación compacta de división de ondas de la siguiente manera:

\[\mathbf{U}^{n+1}_{ i} = \mathbf{U}^{n}_{ i} - \sum_{k=1} ^{NE} \left( \delta\mathbf{M}_{i,k}^{-} \right)^n \frac { l_k } {A_i} \Delta t\]

con

\[\delta\mathbf{M}_{i,k}^{-} = \sum^{3 }_{m=1} \left( \widetilde{\lambda}^-\; \theta \alpha \widetilde{\mathbf{ e}} \right)^{m}_{ k}\]

El uso de es eficiente al trabajar con condiciones de contorno. Al mismo tiempo, garantiza la conservación. En se demostró que, para un esquema numérico escrito en forma de división, la cantidad total de contribuciones calculadas dentro del dominio en cada borde de celda es igual al balance de los flujos que atraviesan el límite del dominio, demostrando la conservación exacta.

Optimizaciones numéricas¶

Una vez calculadas las propagaciones de ondas en $\delta\mathbf{M}_{i}^-$ en, el método de primer orden puede aplicarse promediando las contribuciones de los problemas de Riemann (RP) locales que conforman el contorno de la celda.

La solución aproximada siempre se construye como una suma de saltos o choques, incluso en casos que implican rarefacciones. Un problema ampliamente descrito de los solucionadores linealizados es la violación de la entropía en rarefacciones sónicas, que produce valores negativos de profundidad en las ecuaciones de aguas someras, incluso en ausencia de términos fuente. La solución se restablece mediante una redefinición adecuada de la solución aproximada por medio de correcciones de entropía.

La linealización temporal y espacial de los términos fuente en también puede tener consecuencias negativas, ya que pueden surgir inestabilidades numéricas al aproximar su valor. Su influencia sobre las soluciones aproximadas del RP es clave para construir correcciones apropiadas que eviten resultados no físicos. En se mostró cómo pueden evitarse los errores en los enfoques integrales realizados sobre los términos fuente si se imponen restricciones basadas en la física sobre la solución aproximada. Simplemente modificando los coeficientes de intensidad del término fuente $\beta$, se restauran las soluciones correctas cuando es necesario.

Región de estabilidad¶

Una vez aplicadas las correcciones numéricas, la región de estabilidad para el caso homogéneo puede utilizarse para calcular el tamaño del paso de tiempo. En el marco 2D, considerando mallas no estructuradas, la distancia relevante, que se denominará $\chi_i$ en cada celda $i$, debe considerar el volumen de la celda y la longitud de los $k$ bordes compartidos.

\[\chi_i = \frac {A_i}{\max_{k=1,NE} l_k}\]

Considerando que cada RP $k$ se utiliza para entregar información a un par de celdas vecinas de distinto tamaño, es relevante la distancia $\min(A_i,A_j)/l_k$. El paso de tiempo está limitado por

\[\Delta t \leq CFL \; \Delta t^{\widetilde{\lambda}} \qquad \Delta t^{\widetilde{\lambda}}= \frac { \min (\chi_i,\chi_j) }{ \max |\widetilde{\lambda}^{m} | }\]

con $CFL$=½, ya que la construcción de esquemas de volúmenes finitos mediante la aplicación directa de flujos unidimensionales conduce a rangos de estabilidad reducidos.

El método de solución de RiverFlow2D utiliza pasos de tiempo variables. El paso de tiempo máximo permitido está controlado por el número de Courant-Friederich-Lewy (CFL) establecido por el usuario, que es proporcional al tamaño local de la celda, pero también inversamente proporcional a la velocidad y la profundidad. Las celdas más pequeñas conducen a pasos de tiempo más pequeños. El valor teórico máximo de CFL es 1, pero en algunas ejecuciones puede ser necesario reducir este número a valores menores.

Condiciones de contorno abiertas¶

Hay dos tipos principales de condiciones de contorno que pueden utilizarse en RiverFlow2D: contornos abiertos, donde el flujo puede entrar o salir del área de modelización, y contornos cerrados, que son paredes sólidas sin flujo. No hay restricción en el número de contornos de entrada o salida. Esta sección describe las condiciones de contorno abiertas.

RiverFlow2D permite tener cualquier número de contornos de entrada y salida con diversas combinaciones de condiciones impuestas. El uso adecuado de estas condiciones es un componente crítico de una simulación exitosa de RiverFlow2D. La teoría de las ecuaciones de aguas someras indica que, para un flujo subcrítico bidimensional, se requiere proporcionar al menos una condición en los contornos de entrada y una en los contornos de salida. Para un flujo supercrítico, todas las condiciones deben imponerse en los contornos de entrada y no debe imponerse ninguna condición de contorno en los contornos de salida. La tabla siguiente ayuda a determinar qué condiciones utilizar en la mayoría de las aplicaciones.

Subcritical: Q or Velocity; Water Surface Elevation
Supercritical: Q and Water Surface Elevation; None

Note

Se recomienda tener al menos un contorno donde se prescriba Water Surface Elevation o stage-discharge (por ejemplo, Uniform Flow). Tener únicamente discharge y ninguna condición de Water Surface Elevation puede provocar inestabilidades debido a la violación de los requisitos teóricos de las condiciones de contorno de las ecuaciones de aguas someras.

Las opciones de condiciones de contorno abiertas se describen en la tabla siguiente.

& Impone Water Surface Elevation. Debe proporcionarse un archivo de condición de contorno asociado. - 5: Impone el caudal de agua y Water Surface Elevation. - 6: Impone el caudal de entrada de agua. - 9: Impone una tabla de clasificación stage-discharge de valor único. - 10: Condición de entrada o salida "Free\". El modelo calcula las velocidades y las elevaciones de la superficie del agua. - 11: Condición de salida "Free\". El modelo calcula las velocidades y las elevaciones de la superficie del agua, pero solo se permite el flujo hacia el exterior. - 12: Condición de salida Uniform Flow. & Impone Water Surface Elevation y fuerza direcciones de velocidad perpendiculares. Debe proporcionarse un archivo de condición de contorno asociado. - 18: Impone Water Surface Elevation y concentraciones de sedimentos o contaminantes. También fuerza direcciones de velocidad perpendiculares. Debe proporcionarse un archivo de condición de contorno asociado. - 26: Impone el caudal de entrada de agua y sedimentos. Debe proporcionarse un archivo de condición de contorno asociado.

Note

Si es necesario imponer condiciones abiertas en segmentos de contorno adyacentes, hágalo de modo que cada segmento esté separado por un espacio de más de una celda. Establecer dos o más condiciones abiertas sin esta separación provocará una detección incorrecta de los contornos abiertos.

Tipos de condiciones de contorno de variable única (BCTYPE 1 y 6)¶

Al imponer una sola variable (elevación de la superficie del agua o Q), el usuario debe proporcionar una serie temporal para la variable correspondiente. Para modelar un estado estacionario, la serie temporal debe contener valores constantes para todos los tiempos. No hay restricción sobre el intervalo de tiempo utilizado para la serie temporal. Al imponer la elevación de la superficie del agua, es importante comprobar que el valor impuesto sea superior a la elevación del lecho.

Caudal de agua convertido en velocidades (BCTYPE 6)¶

En esta condición de entrada, el programa calcula el área de flujo y la velocidad media del agua correspondiente al caudal impuesto, que puede variar con el tiempo. A continuación, la velocidad se asigna a cada celda suponiendo una dirección perpendicular a la línea de contorno, como se muestra:

Tabla de clasificación de caudales (BCTYPE 9)¶

Cuando se utiliza una condición stage-discharge de valor único, el modelo calcula primero el caudal en el contorno, luego interpola la elevación correspondiente de la superficie del agua a partir de la tabla de clasificación e impone ese valor para el paso de tiempo siguiente. Si el contorno está seco, funciona como un contorno con condición "Free\". Las elevaciones de la superficie del agua se imponen únicamente en los nodos mojados. Esta condición requiere proporcionar un archivo ASCII con las entradas de los valores de la tabla. Véase la sección para obtener detalles sobre el formato del archivo.

Dado que estas condiciones pueden generar reflexiones de ondas que se propaguen aguas arriba, es importante ubicar el contorno aguas abajo en un tramo suficientemente alejado del área de interés, minimizando así los efectos artificiales de remanso. Desafortunadamente, no existe una forma general de seleccionar dicho lugar, pero será necesario experimentar numéricamente con el modelo real para lograr una ubicación razonable.

Note

En la mayoría de los ríos de poca pendiente, la relación stage-discharge está afectada por histéresis. En otras palabras, la curva stage-discharge forma un bucle, con caudales mayores en la rama ascendente que en la rama de recesión del hidrograma. Esto se debe principalmente al gradiente de profundidad en la dirección del flujo, que cambia de signo a lo largo del hidrograma. En la práctica, esto implica que puede haber dos niveles posibles para el mismo caudal. Las relaciones stage-discharge en bucle no se consideran en esta versión de RiverFlow2D.

Contornos abiertos "Free\" (BCTYPE 10, 11)¶

En los contornos con condición Free, el modelo calcula las velocidades y las elevaciones de la superficie del agua aplicando las ecuaciones completas a partir de las celdas internas. En la práctica, esto equivale a suponer que las derivadas de las elevaciones de la superficie del agua y de las velocidades son 0. En situaciones de flujo subcrítico, es aconsejable utilizar estas condiciones únicamente cuando haya al menos otro contorno abierto donde se imponga la elevación de la superficie del agua o stage-discharge. BCTYPE 10 permite la salida y la entrada de agua, mientras que BCTYPE 11 solo permite que el flujo salga de la malla.

Condición de contorno Uniform Flow (BCTYPE 12)¶

Para aplicar esta condición de contorno, el usuario proporciona únicamente la pendiente del lecho $S_0$. El modelo utilizará $S_0$, el n de Manning y el caudal para crear una tabla de clasificación. Después, para cada intervalo de tiempo, el programa impondrá la elevación de la superficie del agua correspondiente al caudal del contorno, interpolando en la tabla de clasificación. La tabla de clasificación se calcula cada 0.05 m (0.16 ft.) comenzando desde la elevación mínima del lecho en la sección transversal de salida hasta 50 m (164 ft.) por encima de la elevación máxima del lecho en la sección. Si $S_0=-999$, el modelo calculará la pendiente media del lecho perpendicular a la línea de contorno.

Implementación numérica de los contornos abiertos¶

Muchos modelos de simulación se basan en esquemas numéricos fiables y conservativos. Al intentar ampliar su aplicación a problemas realistas que implican geometrías irregulares en los contornos, se debe prestar especial atención a preservar las propiedades del esquema original. La conservación, en particular, se deteriora si los contornos se discretizan sin cuidado.

En las celdas que forman la región de caudal de entrada, el flujo se caracteriza por el signo negativo del siguiente producto escalar en los $k_{\Gamma}$ bordes de contorno

\[\mathbf{q}_{i}\cdot\mathbf{n}_{i,k_{\Gamma}}= (h \mathbf{u})_{i}\cdot\mathbf{n}_{i,k_{\Gamma}} < 0\]

y por el estado del flujo, definido comúnmente mediante el número de Froude

\[Fr_{i}=\frac{\mathbf{u}_{i}\cdot\mathbf{n}_{i,k_{\Gamma}}}{c_i}\]

con $c_i=\sqrt{gh_i}$. Cuando el número de Froude definido en es mayor que uno, el flujo es supercrítico y todos los eigenvalores siguientes son negativos:

\[\lambda^{1}=\mathbf{u}_{i}\cdot\mathbf{n}_{i,k_{\Gamma}}+c_i<0 \qquad \lambda^{2}=\mathbf{u}_{i}\cdot\mathbf{n}_{i,k_{\Gamma}}<0 \qquad \lambda^{3}=\mathbf{u}_{i}\cdot\mathbf{n}_{i,k_{\Gamma}}-c_i<0\]

por lo tanto, deben imponerse los valores de h, u, v y $\phi$. La concentración de soluto del agua $\phi$ es independiente de los eigenvalores y, por lo tanto, debe proporcionarse en la región de entrada para todos los regímenes de flujo.\ Las celdas de la región de caudal de salida se definen mediante

\[\mathbf{q}_{i}\cdot\mathbf{n}_{i,k_{\Gamma}}= (h \mathbf{u})_{i}\cdot\mathbf{n}_{i,k_{\Gamma}} > 0\]

para flujo supercrítico, todos los eigenvalores siguientes son positivos:

\[\lambda^{1}=\mathbf{u}_{i}\cdot\mathbf{n}_{i,k_{\Gamma}}+c_i<0 \qquad \lambda^{2}=\mathbf{u}_{i}\cdot\mathbf{n}_{i,k_{\Gamma}}<0 \qquad \lambda^{3}=\mathbf{u}_{i}\cdot\mathbf{n}_{i,k_{\Gamma}}-c_i<0\]

en consecuencia, no se requiere información adicional.

Cuando en las regiones de caudal de entrada y salida el flujo es subcrítico, la información de actualización no está completa. Lo mismo ocurre en los bordes de celda que actúan como paredes sólidas, que no pueden ser atravesadas por el flujo. Habitualmente, la información adicional proporcionada aguas arriba y aguas abajo son funciones de caudal. Y, en los contornos sólidos, se define una función de caudal normal cero.

Decidir si una entrada o salida es supercrítica o subcrítica no es fácil en una malla 2D. Una caracterización basada en celdas del régimen de flujo en los contornos conduce a situaciones complicadas tanto desde el punto de vista físico como numérico. Por otro lado, las condiciones de contorno físicas o externas suelen referirse a cantidades medias, como el nivel de la superficie del agua o el caudal total, que deben traducirse a profundidad del agua o velocidad en cada celda, dependiendo del criterio del especialista. Para manejar estas situaciones, se requiere una conexión adecuada entre los modelos bidimensionales y unidimensionales en los contornos abiertos. El número de Froude de la sección se define una vez que la sección de contorno tiene un nivel de agua uniforme:

\[Fr_{s}=\frac{w}{\sqrt{g(S_T/b_T)}}\]

siendo la velocidad de la sección transversal $w=Q/S_T$ y definiendo la sección transversal húmeda total $S_T$ y el ancho total como:

\[S_{T}=\sum_{j=1}^{NB}S_j=\sum_{j=1}^{NB}h_jl_j \qquad, b_{T}=\sum_{j=1}^{NB}l_j\]

donde $NB$ es el número de celdas de contorno mojadas, $l_j$ es la longitud de cada borde que conforma el contorno mojado y $h_j$ es la profundidad del agua en cada celda de contorno.

Contorno de caudal de entrada¶

Esta es una de las condiciones de contorno que plantea mayores dificultades, porque debe definirse una representación correcta y conservativa del flujo entrante estacionario o no estacionario, y no existe una forma obvia única de implementarla. El hidrograma de caudal de entrada total $Q = Q(t)$ es la función habitual proporcionada en las simulaciones de inundación, y es importante analizar la mejor forma de imponerlo, ya que involucra toda la sección transversal de entrada y se trabaja con una representación 2D discreta en celdas computacionales. Pueden encontrarse distintos casos.

Casos simples¶

Cuando la sección transversal de entrada tiene forma rectangular, es decir, un fondo plano y limitada por paredes verticales, la sección transversal de entrada mojada es simplemente rectangular.

El caudal total de entrada en el tiempo $t$, $Q_I(t)$, puede distribuirse a lo largo de la sección transversal de entrada utilizando un caudal constante por unidad de ancho, $q_I (m^2s^{-1})$, que puede calcularse como

\[q_{j}=q_{I}=\frac{Q_I(t)}{b_T}\]

En este caso simple, $q_I$ es uniforme a lo largo del contorno de entrada, al igual que el módulo resultante de la velocidad, $w=q_I/h$, con $w=(u^2+v^2)^{1/2}$. Cabe señalar que la dirección del caudal entrante no es necesariamente la misma que la dirección normal al contorno de entrada. Sin embargo, esta dirección suele elegirse como información predeterminada.

Casos complejos¶

En problemas reales de geometría general, la sección transversal de entrada puede cambiar de forma a medida que cambia el nivel del agua (contorno de secado/mojado), y también cambia el número de celdas de contorno involucradas.

Al trabajar con secciones de entrada como la de la Figura , un valor uniforme de $q_I$ como el de conduce a un estado completamente irrealista, con agua más rápida en los bordes de la sección y más lenta en el centro de la sección transversal. Como las velocidades resultantes dependen del valor de la profundidad del agua $h$, aparecerán valores mayores en aquellas celdas donde la profundidad del agua sea menor.

Para buscar una distribución más adecuada, se impone un módulo uniforme de la velocidad del agua $w$ en toda la sección transversal del contorno de entrada. En este caso, el caudal unitario en cada celda de contorno $j$ es variable y se define dependiendo tanto del área total de la sección transversal, $S_T$, como del área transversal de cada celda individual, $S_j$, de la siguiente manera:

\[q_{j}=Q_I\frac{S_j}{S_T l_j}\]

Por otro lado, la actualización de los valores de profundidad del agua en las celdas de entrada proporcionada por el esquema numérico conduce, en el caso general, a un conjunto de nuevas profundidades del agua $h_j^{n+1}$ asociadas, en general, a diferentes niveles de la superficie del agua $d_j$ $d_j=h_j+z_j$.

$Evaluación de $d_{min}$.$

Para nuestros fines, se requiere un nivel horizontal de la superficie del agua en esa región, con el fin de facilitar la traducción entre los puntos de vista 2D y 1D en el contorno abierto. El valor de ese nivel uniforme de agua de la sección transversal se fija teniendo en cuenta la conservación de la masa, es decir, la redistribución conservativa del volumen de agua. Se encuentra el valor mínimo de los niveles de agua entre todas las celdas mojadas del contorno de entrada, $d_{min}$, y se evalúa el volumen de agua $V_S$ almacenado en la sección de entrada por encima de $d_{min}$ como

\[V_{S}=\sum_{j=1}^{NB}(d_j-d_{min})A_j|_{d_j>d_{min}}\]

y se define la superficie mojada por encima de ese nivel, $A_w$:

\[A_{w}=\sum_{j=1}^{NB}A_j|_{d_j>d_{min}}\]

Estos valores se utilizan para redistribuir el volumen sobre la sección de entrada, manteniendo constante el ancho de la sección mojada $b_T$. Como muestra la Figura 3, se obtiene un nuevo nivel uniforme del agua en la sección, $d_S$, mediante:

\[d_{s}=d_{min}+\frac{V_S}{A_w}\]

Además de ayudar a decidir el régimen de flujo en el contorno, las modificaciones descritas anteriormente facilitan el tratamiento de las condiciones de entrada supercríticas. Al modelar un flujo fluvial no estacionario, pueden encontrarse picos altos en el hidrograma. Si esos picos no se manejan correctamente desde el punto de vista numérico, pueden dar lugar a estados supercríticos locales y poco realistas en el contorno de entrada.

En ese caso de flujo de entrada supercrítico, es necesario especificar todas las variables en las celdas del contorno de entrada. Sin embargo, en muchos problemas prácticos solo se dispone del hidrograma de caudal como función del tiempo, sin datos, por lo general, sobre la distribución del nivel del agua o la dirección del caudal en el contorno de entrada.

La alternativa propuesta consiste, cuando el número de Froude de entrada es mayor que 1

\[Fr_{s}=\frac{w}{\sqrt{g(S_T/b_T)}}> 1\]

en imponer un número de Froude máximo, $Fr_{s,max}$, al flujo de entrada. Para ello, manteniendo el ancho de la sección $b_T$, se calcula una nueva área de la sección transversal mojada de entrada, $S_T^*$, a partir del $Fr_{s,max}$ impuesto:

\[S_T^*=(\frac{Q_I^2}{g Fr_{s,max}^2/b_T)})^{1/3}\]

Si $S_T^*$ es mayor que $S_T$, proporciona un nuevo nivel de la superficie del agua para la sección de entrada, $d^*$, también mayor que $d_s$. El incremento asociado en el volumen de agua se equilibra mediante una reducción del caudal impuesto $Q_I(t)$ en ese paso de tiempo.

En ocasiones, ambas condiciones, $Q_I(t)$ y $d(t)$, se conocen en entradas supercríticas. En esos casos, imponer ambos datos en el contorno de entrada es suficiente. Sin embargo, debido al método de integración temporal discreta utilizado, este procedimiento no sigue el criterio de conservación de la masa. Para garantizar que se conserve el balance de masa, se impone una de las condiciones y se modifica la otra, de modo que los flujos calculados en el paso siguiente conduzcan a la conservación de la masa. La mejor solución es imponer directamente el nivel global de la superficie del agua en la sección del contorno de entrada, $d(t)$, y adaptar el caudal de entrada discreto para garantizar que se conserve el volumen final. El valor impuesto de $d$ establece un volumen de entrada que puede transformarse en caudal dividiéndolo por el paso de tiempo. Este valor se suma al caudal, lo que conduce a un balance de masa correcto.

Cuando la celda de contorno pertenece a un contorno abierto donde se impone el caudal de entrada y el flujo es subcrítico, el caudal se calcula utilizando y se impone en la celda de contorno. Además, el nivel del agua se calcula como resultado de las contribuciones de los otros bordes de celda en al actualizar los valores conservados de la celda de contorno en el nivel de tiempo $n+1$, y se redistribuye cuidadosamente como se explicó anteriormente.

Contornos de salida¶

El análisis del flujo en el contorno de salida es más sencillo. Para una salida supercrítica no es necesario imponer condiciones externas. En RiverFlow2D, se realiza un barrido preliminar sobre las celdas mojadas del contorno de salida para evaluar el número de Froude de la celda. Si se encuentra una celda supercrítica, todo el flujo de la sección del contorno de salida se considera supercrítico y no es necesario imponer ninguna condición externa. De lo contrario, todas las celdas se encuentran en estado subcrítico y reciben un tratamiento análogo al del contorno de entrada descrito anteriormente. Como antes, se genera un nivel uniforme de agua en la sección transversal y se establece una distribución de velocidades en los casos en que una curva de clasificación de caudales es la condición de contorno que debe imponerse.

Contornos cerrados¶

Los contornos cerrados son paredes rígidas o sólidas que bloquean completamente el flujo, como las orillas de los ríos o las islas. Constituyen paredes verticales que el flujo nunca puede sobrepasar. Cerca de estos contornos existe una subcapa viscosa muy delgada que requeriría celdas extremadamente pequeñas para resolverse adecuadamente. RiverFlow2D utiliza una condición de deslizamiento en los contornos cerrados y el modelo establecerá un flujo normal cero a través del contorno, pero permitirá velocidades tangenciales. RiverFlow2D detecta automáticamente los contornos cerrados.

Este tipo de condición de contorno no requiere ningún tratamiento especial. Como ningún flujo debe atravesar el contorno, la condición física $\mathbf{u}\cdot\mathbf{n}=0$ se impone sobre la velocidad de celda $\mathbf{u}$ después de añadir todas las contribuciones de ondas del resto de los bordes de la celda, donde $\mathbf{n}$ es la normal a la pared sólida. En otras palabras, si el contorno está cerrado, el borde de celda asociado $k_{\Gamma}$ es una pared sólida, con una componente de velocidad normal cero. Como no hay contribuciones de ese borde, se establece $\delta\mathbf{M}_{i,k_{\Gamma}}^{-}=0$ en al actualizar los valores conservados de la celda de contorno en el nivel de tiempo $n+1$.

Modelado de celdas secas/mojadas¶

RiverFlow2D puede simular el secado y el mojado del lecho. Esta capacidad del modelo es importante al simular la propagación de una onda de inundación por un cauce inicialmente seco. En este caso, tanto el lecho del cauce como la llanura de inundación se inundarán. El lecho del cauce también puede volver a secarse cuando la onda de inundación retrocede.

En RiverFlow2D, la malla de celdas triangulares puede cubrir áreas secas y mojadas, y el modelo manejará estas condiciones utilizando dos algoritmos distintos y dependiendo de la siguiente clasificación de celdas.

Definiciones de celdas basadas en condiciones secas y mojadas¶

Una celda se considera seca si su profundidad de agua es inferior a una fracción de milímetro. No existe una situación de celda parcialmente seca. Un borde de celda se considera inactivo si separa dos celdas secas y se excluye del cálculo. De lo contrario, el borde de celda siempre contribuye a la actualización de las variables en ambos lados. La denominada situación seca/mojada tiene lugar en un borde de celda cuando se cumplen todas las condiciones siguientes:

Una de las celdas vecinas está mojada y la otra está seca.
El nivel del agua en la celda mojada está por debajo del nivel del lecho en la celda seca.
El flujo es subcrítico.

En ese caso, el procedimiento que debe seguirse está bien descrito en.

El algoritmo de secado y mojado de RiverFlow2D es una adaptación del propuesto originalmente por y mejorado posteriormente por y en el contexto de volúmenes finitos, y funciona de la siguiente manera:

Al comienzo de cada paso de tiempo, todas las celdas se clasifican como mojadas o secas de acuerdo con la definición.
Si una celda está seca y completamente rodeada de celdas secas, se elimina de los cálculos y las componentes de velocidad se establecen en cero durante el paso de tiempo en curso.
Todos los bordes de celda internos se clasifican como activos o inactivos de acuerdo con la definición.
Las contribuciones de los bordes de celdas secas/mojadas se calculan suponiendo que el borde es un contorno sólido y que las velocidades a ambos lados se establecen en cero.
El resto de las contribuciones de los bordes de celda se calcula de acuerdo con el esquema numérico descrito anteriormente.
Las celdas mojadas y las celdas secas rodeadas por al menos una celda mojada se mantienen en el cálculo y se resuelven con el esquema de actualización utilizando las contribuciones de los bordes de celda.

Este método genera soluciones numéricas estables sin velocidades espurias sobre áreas secas y ofrece errores de conservación de masa cercanos a la precisión de máquina, permitiendo utilizar la condición CFL clásica.

Conservación de volumen¶

La conservación de volumen o balance de volumen en el dominio de simulación puede definirse mediante una integral de contorno del caudal:

\[\Delta M(\Delta t)=\int_t^{t+\delta t} (\mathbf{Q}_{I}\cdot\mathbf{n}_{I}-\mathbf{Q}_{O}\cdot\mathbf{n}_{O})dt\]

donde $\mathbf{Q}_{I}$ y $\mathbf{Q}_{O}$ son las funciones de caudal total en los contornos de entrada y salida, respectivamente, y $\mathbf{n}_{I}$ y $\mathbf{n}_{O}$ son los vectores normales a los contornos. El caudal normal en las paredes sólidas es cero. Este balance se evalúa integrando en el contorno, celda por celda, de la siguiente manera

\[\Delta M(\Delta t)=\sum_{j=1}^{NB_I}q_{I,j}l_j (\mathbf{n}_{I}\cdot\mathbf{n}_{j}) \Delta t - \sum_{m=1}^{NB_O}q_{O,m}l_m (\mathbf{n}_{O}\cdot\mathbf{n}_{m}) \Delta t\]

donde $\mathbf{n}_{j}$ y $\mathbf{n}_{m}$ son las direcciones del flujo en las celdas de entrada y salida, respectivamente.

La variación de volumen en el dominio de cálculo solo puede deberse a

\[\Delta M(\Delta t) \neq 0\]

Por lo tanto, el error de masa de la solución numérica se mide comparando la cantidad total de agua calculada en el tiempo $t+\Delta t$

\[Vol(t+\Delta t)=\sum_{i=1}^{NCELLS} h_{i}^{n+1} S_i\]

con la cantidad total de agua existente en el tiempo $t$

\[Vol(t)=\sum_{i=1}^{NCELLS} h_{i}^{n} S_i\]

de la siguiente manera

\[Error = \left[Vol(t+\Delta t)- Vol(t)\right]- \Delta M(\Delta t)\]

Esto suele expresarse en términos relativos de la siguiente manera:

\[Relerror = \frac{\left[Vol(t+\Delta t)- Vol(t)\right]- \Delta M(\Delta t)}{Vol(t)+ \Delta M(\Delta t)}\]

Coeficientes de rugosidad n de Manning¶

El n de Manning que normalmente se estima para determinar las pérdidas de carga en el flujo de canales y ríos es una medida global que tiene en cuenta no solo los efectos de la rugosidad del lecho, sino también la fricción interna y las variaciones en la forma y el tamaño de la sección transversal del canal, las obstrucciones y la sinuosidad del río (Ven Te Chow, 1959). Por lo tanto, las estimaciones de n de Manning aplicables a modelos 1D deben ajustarse, porque las ecuaciones de los modelos 2D consideran el intercambio bidimensional de cantidad de movimiento dentro de la sección transversal, que en la simplificación 1D se concentra en un único valor. Varios investigadores han encontrado en aplicaciones prácticas de modelos 2D que los valores de n requeridos pueden ser un $30\%$ menores que los utilizados normalmente para modelos 1D en el mismo tramo de río (Belleudy, 2000). Sin embargo, los modelos 2D no tienen en cuenta la fricción lateral; por lo tanto, la selección final de los coeficientes n de Manning debe ser el resultado de un proceso de calibración en el que los resultados del modelo se ajusten a datos medidos.