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 hidrograma. El costo de los modelos numéricos tridimensionales no simplificados puede evitarse mediante las ecuaciones bidimensionales (2D) de aguas someras promediadas en profundidad.

En aplicaciones realistas de las ecuaciones de aguas someras siempre aparecen términos fuente que describen la variación del nivel del lecho y la fricción del lecho y que, si no se discretizan correctamente, pueden producir inestabilidades numéricas. Durante la última década, el esfuerzo principal se ha centrado en mantener un equilibrio discreto entre los términos de flujo y fuente en casos de agua en reposo, lo que condujo a la noción de esquemas bien equilibrados o propiedad C [, , , ]. Recientemente se han presentado solucionadores de Riemann aproximados aumentados para incluir correctamente el efecto de los términos fuente en la solución débil [Rosatti et al. (2003), ]. De este modo pueden calcularse soluciones precisas sin imponer parámetros de ajuste dependientes del caso, que se utilizan con frecuencia para evitar valores negativos de 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 HydroBID Flood; las referencias amplían esta información.

Supuestos del modelo hidrodinámico¶

HydroBID Flood utiliza las ecuaciones de aguas someras obtenidas mediante la integración vertical de la ecuación de Navier-Stokes. Por tanto, el modelo no calcula aceleraciones 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 turbulencia y las pérdidas de energía se contabilizan ú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 espacial y temporal de la densidad, la viscosidad y el esfuerzo de fluencia.

Ecuaciones de flujo con 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, junto 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, con $h$ como profundidad del agua, $q_x=uh$ y $q_y=vh$ como caudales unitarios, y $(u,v)$ como componentes promediadas en profundidad del vector de velocidad $\mathbf{u}$ a lo largo de las coordenadas $(x,y)$, respectivamente. Los vectores de flujo son:

\[\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 obtienen al suponer una distribución hidrostática de 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 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 nivel del fondo $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:

\[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}}\]

siendo $n$ el coeficiente de rugosidad.

Con 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 está 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 malla $\Omega$ mediante 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$ mediante una malla fija. Suponiendo una representación por tramos de las variables conservadas, como se muestra en la figura siguiente, 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 mediante 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 con los autovectores 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}\]

donde $\widetilde{\boldsymbol{\Lambda}}_{k}$ es una matriz diagonal con los autovalores $\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 del vector $\mathbf{U}$ a través del borde de la malla como el término fuente se proyectan sobre la base de autovectores 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 amplitudes de onda y $\mathbf{B} = ( \beta^1 , \beta^2 , \beta^3 )^{T}_{k}$ contiene las amplitudes de 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 introducirse en da 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 en cada celda $i$ y en cualquier celda se cumple la siguiente propiedad geométrica

\[\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 mediante una formulación compacta de separación de ondas:

\[\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 esta formulación es eficiente al tratar las condiciones de contorno y, al mismo tiempo, garantiza la conservación. En se demostró que, para un esquema numérico escrito en forma separada, la cantidad total de contribuciones calculadas dentro del dominio en cada borde de celda es igual al balance de los flujos que atraviesan el contorno del dominio, lo que demuestra la conservación exacta.

Optimizaciones numéricas¶

Una vez calculadas las propagaciones de onda 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 u ondas de choque, incluso en casos que incluyen rarefacciones. Un problema ampliamente descrito de los solucionadores linealizados es la violación de 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 recupera mediante una redefinición adecuada de la solución aproximada usando 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 RP aproximadas es la clave para construir correcciones adecuadas que eviten resultados no físicos. En se mostró que los errores de las aproximaciones integrales realizadas sobre los términos fuente pueden evitarse imponiendo restricciones físicas sobre la solución aproximada. Basta modificar los coeficientes de amplitud de fuente $\beta$ para recuperar las soluciones correctas cuando sea necesario.

Región de estabilidad¶

Una vez aplicadas las correcciones numéricas, la región de estabilidad del 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 bordes compartidos $k$.

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

Como cada RP $k$ se utiliza para transmitir información a un par de celdas vecinas de distinto tamaño, la distancia $\min(A_i,A_j)/l_k$ es relevante. 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 HydroBID Flood utiliza pasos de tiempo variables. El paso de tiempo máximo permitido está controlado por el número de Courant-Friederich-Lewy (CFL) definido 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 producen pasos de tiempo menores. El valor teórico máximo de CFL es 1, aunque en algunas ejecuciones puede ser necesario reducirlo.

Condiciones de contorno abiertas¶

En HydroBID Flood pueden utilizarse dos tipos principales de condiciones de contorno: contornos abiertos, por los que el flujo puede entrar o salir del área de modelación, y contornos cerrados, que son paredes sólidas sin flujo (véase la figura correspondiente). No hay restricción en el número de contornos de entrada o salida. Esta sección describe las condiciones de contorno abiertas.

HydroBID Flood permite cualquier número de contornos de entrada y salida con distintas combinaciones de condiciones impuestas. El uso correcto de estas condiciones es un componente crítico de una simulación exitosa de HydroBID Flood. La teoría de las ecuaciones de aguas someras indica que, para un flujo bidimensional subcrítico, se requiere proporcionar al menos una condición en los contornos de entrada y otra 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 en los contornos de salida. La tabla siguiente ayuda a determinar qué condiciones utilizar en la mayoría de las aplicaciones.

Subcrítico: Q o velocidad; elevación de la superficie del agua
Supercrítico: Q y elevación de la superficie del agua; ninguna

Note

Se recomienda tener al menos un contorno en el que se prescriba la superficie del agua o la relación nivel-caudal (por ejemplo, Uniform Flow). Es posible que disponer únicamente del caudal, sin una condición de elevación de la superficie del agua, produzca inestabilidades por incumplir 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 la elevación de la superficie del agua. Debe proporcionarse un archivo de condiciones de contorno asociado.\ - 5: Impone el caudal de agua y la elevación de la superficie del agua. - 6: Impone la entrada de caudal de agua. - 9: Impone una tabla de relación nivel-caudal de valor único. - 10: Condición de entrada o salida libre\". El modelo calcula las velocidades y las elevaciones de la superficie del agua. - 11: Condición de salida libre\". El modelo calcula las velocidades y las elevaciones de la superficie del agua, pero solo se permite el flujo hacia afuera. - 12: Condición de salida de flujo uniforme. & Impone la elevación de la superficie del agua y fuerza direcciones de velocidad perpendiculares. Debe proporcionarse un archivo de condiciones de contorno asociado.\ - 18: Impone la elevación de la superficie del agua y las concentraciones de sedimentos o contaminantes. También fuerza direcciones de velocidad perpendiculares. Debe proporcionarse un archivo de condiciones de contorno asociado. - 26: Impone la entrada de caudal de agua y sedimentos. Debe proporcionarse un archivo de condiciones de contorno asociado.

Note

Si necesita imponer condiciones abiertas en segmentos de contorno adyacentes, hágalo de forma que cada segmento esté separado por un hueco mayor que una celda (véase la figura correspondiente). 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 esa variable. 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 en la serie. Al imponer la elevación de la superficie del agua, es importante comprobar que el valor impuesto sea mayor que la elevación del lecho.

Conversión del caudal de agua en velocidades (BCTYPE 6)¶

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

Tabla de relación nivel-caudal (BCTYPE 9)¶

Cuando se utiliza una condición nivel-caudal de valor único, el modelo calcula primero el caudal en el contorno, luego interpola en la tabla de relación la elevación correspondiente de la superficie del agua e impone ese valor en el paso de tiempo siguiente. Si el contorno está seco, funciona como un contorno de condición "libre\". Las elevaciones de la superficie del agua solo se imponen en nodos mojados. Esta condición requiere proporcionar un archivo ASCII con los valores de la tabla. Consulte la sección correspondiente para obtener detalles sobre el formato del archivo.

Como esta condición puede generar reflexiones de onda 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. Lamentablemente, no existe una forma general de seleccionar ese lugar; 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 nivel-caudal está afectada por histéresis. En otras palabras, la curva nivel-caudal forma un bucle, con caudales mayores en la rama ascendente que en la rama de recesión del hidrograma. Esto se debe principalmente a que el gradiente de profundidad en la dirección del flujo cambia de signo a lo largo del hidrograma. En la práctica, esto implica que puede haber dos niveles posibles para el mismo caudal. Esta versión de HydroBID Flood no considera relaciones nivel-caudal en bucle.

Contornos abiertos "libres\" (BCTYPE 10 y 11)¶

En los contornos de condición libre, el modelo calcula las velocidades y las elevaciones de la superficie del agua aplicando las ecuaciones completas 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, se recomienda utilizar estas condiciones solo cuando exista al menos otro contorno abierto donde se imponga la elevación de la superficie del agua o la relación nivel-caudal. BCTYPE 10 permite la entrada y salida de agua, mientras que BCTYPE 11 solo permite el flujo hacia afuera de la malla.

Condición de contorno de flujo uniforme (BCTYPE 12)¶

Para aplicar esta condición de contorno, el usuario proporciona únicamente la pendiente del lecho $S_0$. El modelo utilizará $S_0$, n de Manning y el caudal para crear una tabla de relación. Después, en cada intervalo de tiempo, el programa impondrá la elevación de la superficie del agua correspondiente al caudal del contorno, interpolando en la tabla. La tabla se calcula cada 0.05 m (0.16 ft), desde la menor elevación del lecho en la sección transversal de salida hasta 50 m (164 ft) por encima de la mayor elevación del lecho de 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 extender su aplicación a problemas realistas con geometrías irregulares en los contornos, debe prestarse especial atención a conservar las propiedades del esquema original. En particular, la conservación se deteriora si los contornos se discretizan de forma descuidada.

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 bordes de contorno $k_{\Gamma}$

\[\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 normalmente 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 autovalores 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 autovalores y, por 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 autovalores 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 el estado del flujo es subcrítico tanto en la región de entrada como en la de salida, la información necesaria para la actualización no está completa. Lo mismo ocurre en los bordes de celda que actúan como paredes sólidas que el flujo no puede atravesar. Normalmente, la información adicional proporcionada aguas arriba y aguas abajo son funciones de caudal. En los contornos sólidos se define una función de caudal normal nulo.

En una malla 2D no es fácil decidir si una entrada o salida es supercrítica o subcrítica. 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 otra parte, 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 o velocidad en cada celda según el criterio del profesional. Para tratar estas situaciones, en los contornos abiertos se necesita una conexión adecuada entre los modelos bidimensional y unidimensional. 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 forma 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 presenta más dificultades, porque debe definirse una representación correcta y conservativa del flujo entrante estacionario o no estacionario y no existe una única forma evidente de implementarla. El hidrograma de caudal total de entrada $Q = Q(t)$ es la función habitual en una simulación 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 presentarse distintos casos.

Casos sencillos¶

Cuando la sección transversal de entrada es rectangular (como se muestra en la figura siguiente), es decir, tiene fondo plano y está limitada por paredes verticales, la sección transversal húmeda de entrada es simplemente rectangular.

El caudal total de entrada en el tiempo $t$, $Q_I(t)$, puede distribuirse a lo largo de la sección de entrada mediante 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 sencillo, $q_I$ es uniforme a lo largo del contorno de entrada y también lo es el módulo resultante de la velocidad, $w=q_I/h$, con $w=(u^2+v^2)^{1/2}$. Debe observarse que la dirección del caudal entrante no tiene por qué coincidir con 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 al cambiar el nivel del agua (contorno de secado/mojado), y también cambia el número de celdas de contorno involucradas (como se muestra en la figura siguiente).

En secciones de entrada como la mostrada en la figura anterior, un valor uniforme de $q_I$ 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. Como las velocidades resultantes dependen del valor de la profundidad del agua $h$, aparecerán valores mayores en las celdas donde la profundidad sea menor.

Para obtener 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 de cada celda de contorno $j$ es variable y se define en función tanto del área total de la sección transversal, $S_T$, como del área transversal individual de la celda, $S_j$, de la siguiente manera:

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

Por otra parte, la actualización de los valores de profundidad del agua en las celdas de entrada proporcionada por el esquema numérico conduce, en general, a un conjunto de nuevas profundidades $h_j^{n+1}$, asociadas a distintos niveles de superficie del agua $d_j$, con $d_j=h_j+z_j$, como se muestra en la figura siguiente.

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

Para nuestros fines, en esa región se requiere un nivel horizontal de la superficie del agua que ayude a traducir entre los puntos de vista 2D y 1D en el contorno abierto. El valor de ese nivel uniforme de la sección 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 el volumen de agua $V_S$ almacenado en la sección de entrada por encima de $d_{min}$ se evalúa como

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

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

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

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

\[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 facilitan el tratamiento de las condiciones de entrada supercríticas. Al modelar un flujo fluvial no estacionario pueden aparecer picos altos en el hidrograma. Si esos picos no se tratan correctamente desde el punto de vista numérico, pueden producir estados supercríticos locales e irreales en el contorno de entrada.

En el 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 y, por lo general, no hay datos sobre la distribución del nivel de agua ni sobre la dirección del caudal en el contorno de entrada.

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

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

consiste 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 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 superficie del agua para la sección de entrada, $d^*$, también mayor que $d_s$, como se muestra en la figura siguiente. El incremento asociado del volumen de agua se equilibra reduciendo el caudal impuesto $Q_I(t)$ en ese paso de tiempo.

En ocasiones se conocen ambas condiciones, $Q_I(t)$ y $d(t)$, en entradas supercríticas. En esos casos basta con imponer ambos datos en el contorno de entrada. Sin embargo, debido al método discreto de integración temporal utilizado, este procedimiento no cumple el criterio de conservación de la masa. Para garantizar 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 de entrada, $d(t)$, y adaptar el caudal discreto de entrada para garantizar la conservación del 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 y produce 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 mediante y se impone en la celda de contorno. Además, el nivel de agua se calcula como resultado de las contribuciones de los otros bordes de la celda en al actualizar los valores conservados de la celda de contorno en el nivel temporal $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 HydroBID Flood se realiza un barrido preliminar sobre las celdas mojadas del contorno de salida para evaluar el número de Froude de cada celda. Si se encuentra una celda supercrítica, todo el flujo de la sección de salida se considera supercrítico y no es necesario imponer una condición externa. De lo contrario, todas las celdas están en estado subcrítico y reciben un tratamiento análogo al descrito anteriormente para el contorno de entrada. Como antes, se genera un nivel uniforme de agua en la sección y se establece una distribución de velocidades cuando la condición de contorno que debe imponerse es una curva de relación de caudal.

Contornos cerrados¶

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

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

Modelación de celdas secas/mojadas¶

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

En HydroBID Flood, la malla de celdas triangulares puede cubrir áreas secas y mojadas, y el modelo tratará estas condiciones mediante dos algoritmos distintos según la clasificación siguiente de las celdas.

Definiciones de celdas según condiciones secas y mojadas¶

Una celda se considera seca si su profundidad de agua es menor que 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. En caso contrario, el borde de celda siempre contribuye a actualizar las variables de ambos lados. La denominada situación mojada/seca se produce en un borde de celda cuando se cumplen todas las condiciones siguientes:

Una de las celdas vecinas está mojada y la otra seca.
El nivel de 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 se describe detalladamente en.

El algoritmo de secado y mojado de HydroBID Flood 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 según 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 internos de las celdas se clasifican como activos o inactivos según la definición.
Las contribuciones de los bordes en situaciones mojadas/secas se calculan suponiendo que el borde es un contorno sólido y las velocidades de ambos lados se establecen en cero.
Las contribuciones del resto de los bordes se calculan conforme al 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.

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, lo que permite utilizar la condición CFL clásica.

Conservación del volumen¶

La conservación del volumen o balance volumétrico 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 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 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)\]

Normalmente se expresa 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¶

La 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, las variaciones de forma y tamaño de la sección transversal del canal, las obstrucciones y la sinuosidad del río (Ven Te Chow, 1959). Por tanto, deben ajustarse las estimaciones de n de Manning aplicables a modelos 1D, 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 término. Varios investigadores han encontrado, en aplicaciones prácticas de modelos 2D, que los valores de n necesarios 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 ello, 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.