Modelo de derrame de petróleo sobre terreno¶

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 coste 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 fondo y la fricción del fondo 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 términos de flujo y fuente en casos de agua en reposo, dando 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, pueden calcularse soluciones precisas sin imponer los 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 OilFlow2D; la información puede ampliarse en las referencias.

Supuestos del modelo de flujo viscoso¶

OilFlow2D utiliza las ecuaciones de aguas someras resultantes de la integración vertical de la ecuación de Navier-Stokes. Por 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 fondo 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 de Manning's n en las ecuaciones de cantidad de movimiento.
El modelo puede considerar la transferencia de calor para calcular la temperatura del petróleo 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 de 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, junto con todos los supuestos asociados. Ese 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 los componentes promediados 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 han obtenido suponiendo 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 fondo y de las fuerzas tangenciales generadas por el esfuerzo del fondo.

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

donde las pendientes 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 fondo 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 petróleo es prescrita por el usuario y no depende de los cambios ambientales durante la simulación.

Ecuaciones considerando la transferencia de calor¶

Esta sección presenta las ecuaciones que resuelve OilFlow2D para considerar la transferencia de calor del petróleo a la atmósfera y al suelo mientras fluye sobre el terreno. El modelo también determina cómo cambian la densidad, la viscosidad y el esfuerzo de fluencia del petróleo en función de la temperatura durante la simulación.

La ecuación de continuidad para la masa de petróleo se escribe como

\[\frac{\partial (\rho h)}{\partial t} + \frac{\partial}{\partial x} (\rho hu) + \frac{\partial}{\partial y} (\rho hv) = 0\]

y las leyes de conservación de la cantidad de movimiento lineal total a lo largo de los ejes $x$ e $y$ se expresan como

\[\begin{aligned} & \frac{\partial (\rho hu)}{\partial t} + \frac{\partial}{\partial x} (\rho hu^2 + \frac{1}{2} g_n \rho h^2 ) + \frac{\partial}{\partial y} (\rho huv) = -g_n \rho h \frac{\partial z_b}{\partial x} - \tau_{bx} \\ & \frac{\partial (\rho hv)}{\partial t} + \frac{\partial}{\partial x} (\rho huv) + \frac{\partial}{\partial y} (\rho hv^2 + \frac{1}{2} g_n \rho h^2) = -g_n \rho h \frac{\partial z_b}{\partial x} - \tau_{by} \end{aligned}\]

donde $\rho$ es la densidad volumétrica integrada en profundidad [kg/m$^3$ o lb/ft$^3$], $h$ es la profundidad vertical del flujo [m o ft] y ($u,\,v$) son los componentes del vector de velocidad del flujo integrado en profundidad $\mathbf{u}$ [m/s o ft/s], $z_b$ es la elevación del fondo [m o ft] y ($\tau_{bx},\,\tau_{by}$) son los componentes del vector de resistencia basal integrado en profundidad ${\tau_b}$. La ecuación de energía para la temperatura se escribe en forma no conservativa:

\[\frac{\partial T}{\partial t} + u\frac{\partial}{\partial x} T + v\frac{\partial}{\partial y} T = \mathbf{S}_{T}\]

donde el término fuente $\mathbf{S}_{T}$ considera la transferencia de calor del petróleo con el entorno en función de la profundidad del petróleo, $h$, y otros factores. Para cerrar el modelo, se utiliza la siguiente ecuación para determinar la densidad del fluido $\rho$ [kg/m$^3$ o lb/ft$^3$]

\[\rho = \rho_0 + \Lambda (T-T_0)\]

donde $\rho_0$ [kg/m$^3$ o lb/ft$^3$] es la densidad de referencia del petróleo a la temperatura $T_0$ [$^{\circ}$K] y $\Lambda$ es un parámetro experimental.

Propiedades del fluido y leyes de fricción¶

Además de la densidad, la temperatura del flujo también afecta a la viscosidad y al esfuerzo de fluencia, generando cambios en los esfuerzos de fricción entre el flujo y el terreno. La dependencia del esfuerzo de fluencia [Pa o lb/in$^2$] con la temperatura puede introducirse como una tabla de esfuerzos experimentales para cada temperatura o mediante fórmulas de regresión como la siguiente:

\[Ys = 10^{(Ays T^2+ Bys T - Cys)},\]

donde $Ays$, $Bys$ y $Cys$ son parámetros del flujo obtenidos de la caracterización del petróleo. La viscosidad del petróleo, $\mu$ [Pa$\cdot$s o lb$\cdot$s/in$^2$], se rige por la formulación de Andrade:

\[\mu = Av \ e^{\left( {Bv}/{T}\right)}\]

donde $Av$ y $Bv$ son parámetros del flujo obtenidos de la caracterización del petróleo y $T$ está en $^\circ$K. El modelo también ofrece la opción de introducir tablas que representan la variación de la viscosidad, el esfuerzo de fluencia y la densidad en función de la temperatura.

Término fuente de temperatura y mecanismos de transferencia de calor¶

El término fuente de temperatura se calcula como

\[\mathbf{S}_{T} = \dfrac{Q}{\rho C_p h}\]

donde $C_p$ es el calor específico [J/kg$^{\circ}$C o BTU/lb$^{\circ}$F], y $Q$ es el flujo de calor total [W/m$^2$ o BTU/ft$^2$] que cuantifica el intercambio de calor entre el petróleo y el entorno.

OilFlow2D considera la siguiente formulación para representar la transferencia de calor entre el petróleo y el entorno.

Radiación ambiental (petróleo$\rightarrow$entorno):
Radiación térmica solar recibida (sol$\rightarrow$petróleo): calculada en función de la radiación extraterrestre (R0), atenuada por la transmisión atmosférica (at), la nubosidad (ac), el sombreado (cs) y la reflexión (cr).
Transferencia por convección entre el petróleo y el suelo.
Transferencia por convección entre el petróleo y el aire.

El modelo de transferencia de calor considera la convección del aire, la radiación incidente y la radiación emitida. Los mecanismos de transferencia de calor actúan como intercambios de calor paralelos entre la superficie del fluido y el entorno, que pueden expresarse como

\[\mathcal{S_T}=\dfrac{\dot{Q}}{\rho h C_p} = \dfrac{\dot{Q}_{rad}+\dot{Q}_{conv}+\dot{Q}_{em}}{\rho h C_p}\]

La transferencia de calor por radiación entre el petróleo y el aire se calcula mediante la ecuación de Stefan-Boltzmann:

\[\dot{Q}_{em}=\varepsilon\sigma\left(T^4-T_{air}^4\right)\]

donde $\varepsilon$ es la emisividad de la superficie del petróleo (se supone 0.55), $\sigma = 5.67 \times 10^{-8}$ W/(m$^2$K$^{-4}$) es la constante de Stefan-Boltzmann y $T_{air}$ es la temperatura del aire. La convección del aire se calcula en función de la diferencia entre las temperaturas del fluido y del aire, $T$ y $T_{air}$ [K], y de un coeficiente de convección, $h_c$:

\[\dot{Q}_{conv}=h_c\left(T-T_{air}\right)\]

en términos del coeficiente convectivo de transferencia de calor $h_c$ [Wm$^{-2}$K], que introduce el efecto del movimiento del aire. Puede relacionarse con el número de Nusselt para establecer la importancia de la transferencia de calor por convección frente a la conducción:

\[Nu=\frac{h_c L}{k}\left(T-T_{air}\right)\]

donde $L$ [m] es una longitud característica y $k$ [W m/K] es la conductividad térmica.

En el modelo, $h_c$ se calcula en función de la velocidad del viento:

\[h_c = 8.55 + 2.56 \cdot U_w\]

Por último, la radiación incidente, $\dot{Q}_{rad}$, puede obtenerse a partir de datos medidos y se requiere en el modelo como parámetro de entrada.

Retención de petróleo¶

La retención de petróleo es un fenómeno en el que la velocidad de la fase de petróleo se reduce significativamente debido a su interacción con el suelo o la vegetación. Esta retención se produce por una combinación de adsorción, atrapamiento capilar en los espacios porosos y aumento de la resistencia viscosa, lo que reduce la movilidad del petróleo y prolonga su tiempo de permanencia en el subsuelo. En este contexto, la atención se centra en la resistencia viscosa entre el flujo de petróleo y el suelo. Para modelar este proceso se utiliza el enfoque Simple Shear Infinite Landslide ,. En este modelo, el esfuerzo cortante se equilibra con la componente de la fuerza gravitatoria en la dirección horizontal del flujo (véase la Figura ). Este esfuerzo cortante reduce la cantidad de movimiento en la dirección del flujo de petróleo. Sin embargo, este esfuerzo cortante solo afecta a la capa delgada más cercana al fondo. Por ello, debe definirse una profundidad del petróleo, $h_{\text{retention}}$, como el espesor de esta capa. Por tanto, el esfuerzo cortante causado por este fenómeno, $\tau_{\text{retention}}$, en la interfaz suelo-petróleo puede expresarse como:

\[\tau_{\text{retention}} = \rho g h_{\text{retention}} sin (\theta)\]

donde $\mu$ es la viscosidad del petróleo, $\rho$ es la densidad del petróleo, $g$ es la aceleración gravitatoria y $\theta$ es el ángulo de la pendiente del fondo (véase la Figura ). Con esta expresión pueden darse los siguientes casos:

Caso con pendiente ($S_0 = tg(\theta)$): En este caso, el esfuerzo cortante se define mediante la expresión , y la velocidad se considera cero cuando la profundidad del petróleo es $h\leq h_{\text{retention}}$.
Caso sin pendiente ($S_0=0, \; \theta=0$): En este caso, utilizando la expresión , $\tau_{\text{retention}}=0$. Por tanto, la retención de petróleo se produce imponiendo una velocidad de flujo de petróleo igual a cero cuando $h \leq h_{\text{retention}}$.

Es importante observar que la profundidad de retención del petróleo, $h_{\text{retention}}$, representa la profundidad mínima para el flujo de petróleo. En consecuencia, la velocidad de las celdas con una profundidad de petróleo menor que $h_{\text{retention}}$ será cero. Por esta razón, los valores de $h_{\text{retention}}$ deben ser del orden de milímetros.

Ecuaciones del hidrograma de derrame por rotura de tubería¶

El modelo Oil Pipeline Break de OilFlow2D para QGIS prepara un archivo de múltiples fuentes y un hidrograma de caudal para cada punto de rotura de tubería generado. El cálculo se realiza antes de la simulación de OilFlow2D; los archivos de hidrograma resultantes son leídos posteriormente por el modelo como series temporales de fuentes para la simulación del derrame de petróleo sobre terreno. Las ecuaciones siguientes resumen el cálculo utilizado por el módulo de tuberías.

Para un diámetro de tubería $D_p$ y un diámetro de fuga $D_l$, las áreas de la tubería y de la fuga son

\[A_p = {\pi D_p^2 \over 4}\]

\[A_l = {\pi D_l^2 \over 4}\]

La aceleración gravitatoria es $g = 9.80665$ m/s$^2$ cuando el proyecto utiliza unidades métricas, o $g = 32.174$ ft/s$^2$ cuando utiliza unidades inglesas. El coeficiente de pérdida de la fuga asociado al coeficiente de descarga introducido por el usuario $C_d$ es

\[K_l = {1 \over C_d^2} - 1\]

Altura inicial estacionaria de la tubería¶

El módulo calcula primero una altura de energía estacionaria a lo largo de la tubería a partir del caudal inicial $Q_0$ y de la altura de presión inicial introducida en el extremo de la tubería. Si $z_n$ es la elevación en el extremo aguas abajo y $H_{end}$ es la altura de presión introducida en ese extremo, la altura de energía aguas abajo es

\[h_n = z_n + H_{end}\]

Al desplazarse aguas arriba por segmentos, la altura de energía se actualiza como

\[h_i = h_{i+1} + {f_i \Delta x_i \over 2 g D_i A_i^2} Q_0^2\]

donde $h_i$ es la altura de energía en el nodo $i$, $\Delta x_i$ es la longitud del segmento, $D_i$ es el diámetro de la tubería, $A_i$ es el área de la tubería y $f_i$ es el factor de fricción de Darcy-Weisbach. El número de Reynolds utilizado para cada segmento es

\[Re_i = {|Q_0| D_i \over A_i \nu}\]

donde $\nu$ es la viscosidad cinemática. Para flujo laminar, el factor de fricción es

\[f_i = {64 \over Re_i}\]

Para flujo turbulento, el módulo resuelve la relación de Colebrook mediante iteración de Newton-Raphson:

\[ {1 \over \sqrt{f_i}} = -0.86 \ln \left({e_i \over 3.71 D_i} + {2.51 \over Re_i \sqrt{f_i}}\right) \]

donde $e_i$ es la rugosidad de la tubería. Si está habilitada la opción del diálogo para calcular el coeficiente de fricción a partir de la rugosidad, este cálculo de fricción también se utiliza en las ecuaciones de drenaje del derrame. De lo contrario, para el paso de drenaje del derrame se utiliza el coeficiente de fricción introducido por el usuario.

Flujo de fuga aguas arriba¶

Para cada rotura, la altura de presión utilizada en la rotura es

\[h_p = h_b - z_b\]

donde $h_b$ es la altura de energía calculada en la rotura y $z_b$ es la elevación de la rotura. El caudal de fuga aguas arriba en cada paso de tiempo se calcula mediante la ecuación del orificio:

\[Q_l = C_d A_l \sqrt{2 g h_p}\]

El intervalo de cálculo es $\Delta t$. La entrada desde la tubería aguas arriba de la rotura es inicialmente $Q_0$. Si $t_s$ es el tiempo de inicio del cierre de la válvula y $t_c$ es la duración del cierre, el módulo utiliza

\[Q_{in}(t) = \begin{cases} Q_0, & t \leq t_s \\ Q_0 \left(1 - {t - t_s \over t_c}\right), & t_s < t \leq t_s + t_c \\ 0, & t > t_s + t_c \end{cases}\]

En cada paso de tiempo, el volumen fugado, el volumen entrante y el volumen neto saliente son

\[V_l = Q_l \Delta t\]

\[V_{in} = Q_{in} \Delta t\]

\[V_{net} = V_l - V_{in}\]

El volumen disponible de tubería aguas arriba $V_a$ se reduce en función del volumen neto saliente. Mientras quede volumen disponible, la altura de presión se reduce en proporción al volumen restante:

\[V_a^{new} = V_a - V_{net}\]

\[h_p^{new} = h_p {V_a^{new} \over V_a}\]

Cuando se agota el volumen disponible, la altura de presión aguas arriba se establece en cero y termina la contribución aguas arriba.

Drenaje gravitatorio aguas abajo¶

Se supone que la parte aguas abajo de la tubería drena por gravedad desde el volumen de tubería situado aguas abajo de la rotura. Para la longitud aguas abajo contribuyente $L_a$, el volumen disponible es

\[V_a = L_a A_p\]

La pendiente local de la tubería utilizada para bajar la superficie del petróleo durante el drenaje es

\[m = {|z_s - z_b| \over L_a}\]

donde $z_s$ es la elevación al final del segmento contribuyente. En cada paso de tiempo el módulo resuelve la forma cuadrática

\[A Q_l^2 + B Q_l + C = 0\]

con

\[A = {8 \over g \pi^2} \left({1 + K_l \over D_l^4} - {1 \over D_p^4}\right)\]

\[B = f L_a\]

\[C = {P_l - P_s \over \rho g} + z_b - z_s\]

donde $\rho$ es la densidad del petróleo, $P_l$ es la presión en la fuga y $P_s$ es la presión en el extremo aguas arriba del segmento contribuyente. Para el drenaje gravitatorio aguas abajo, el módulo utiliza presión manométrica cero, por lo que el término de presión normalmente es cero. El caudal correspondiente a la raíz positiva es

\[Q_l = {-B + \sqrt{B^2 - 4 A C} \over 2 A}\]

El caudal de fuente escrito en el hidrograma del derrame es $Q_l$. Para la actualización interna del volumen de tubería aguas abajo, el módulo también calcula un término de flujo aguas abajo $Q_o$ cuando existe una contribución activa de un contorno aguas abajo o de una válvula; en caso contrario $Q_o = 0$. El volumen neto retirado del volumen de tubería contribuyente aguas abajo es

\[V_{net} = (Q_l - Q_o) \Delta t\]

y el volumen aguas abajo restante y la longitud contribuyente se actualizan como

\[V_a^{new} = V_a - V_{net}\]

\[L_a^{new} = {V_a^{new} \over A_p}\]

La elevación de la superficie del petróleo utilizada en el paso de tiempo siguiente se reduce en

\[z_s^{new} = z_s - m {V_{net} \over A_p}\]

El hidrograma aguas abajo termina cuando el volumen de tubería restante deja de cambiar más allá de la tolerancia del módulo.

Hidrograma de fuente combinado¶

Para cada punto de rotura, los hidrogramas aguas arriba y aguas abajo se combinan por paso de tiempo:

\[Q_{source}(t) = Q_{upstream}(t) + Q_{downstream}(t)\]

El hidrograma combinado se escribe como la serie temporal de fuente referenciada por el archivo de múltiples fuentes. OilFlow2D aplica cada fuente como una entrada puntual en las coordenadas de rotura generadas.

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

Para introducir el esquema de volúmenes finitos, se integra en un volumen o celda de cuadrícula $\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 la normal unitaria 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 (Figura ) 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 mediante 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}\]

donde $\widetilde{\boldsymbol{\Lambda}}_{k}$ es una matriz diagonal con los 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 del 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 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 produce 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:

\[\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 flujos que atraviesan el contorno del dominio, demostrando una 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 locales (RP) que forman 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 entropía en rarefacciones sónicas , que produce valores negativos de profundidad en las ecuaciones de aguas someras, incluso en ausencia de término fuente. La solución se restablece 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 aproximadas de los RP es clave para construir correcciones apropiadas que eviten resultados no físicos. En se mostró cómo pueden evitarse los errores de los enfoques integrales aplicados a los términos fuente imponiendo restricciones basadas en la física sobre la solución aproximada. Simplemente modificando los coeficientes de intensidad de fuente $\beta$, se restablecen 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 tener en cuenta el volumen de la celda y la longitud de los bordes compartidos $k$.

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

Considerando que cada RP $k$ se utiliza para transmitir 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 a partir de la aplicación directa de flujos unidimensionales conduce a intervalos de estabilidad reducidos.

El método de solución de OilFlow2D 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 producen 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 inferiores.

Condiciones de contorno abiertas¶

En OilFlow2D pueden utilizarse dos tipos principales de condiciones de contorno: contornos abiertos, donde 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 ). No hay restricción en el número de contornos de entrada o salida. Esta sección describe las condiciones de contorno abiertas.

OilFlow2D permite cualquier número de contornos de entrada y salida con diversas combinaciones de condiciones impuestas. El uso correcto de estas condiciones es un componente crítico de una simulación de OilFlow2D exitosa. La teoría de las ecuaciones de aguas someras indica que, para un flujo bidimensional subcrítico, es necesario proporcionar al menos una condición en los contornos de entrada y una en los de salida. Para un flujo supercrítico, deben imponerse todas las condiciones 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.

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 donde se prescriba la superficie del agua o una relación nivel-caudal (por ejemplo, Flujo uniforme). Tener solo una condición de caudal y ninguna condición de elevación de la superficie del agua 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 la elevación de la superficie del agua. Debe proporcionarse un archivo de condición de contorno asociado.\ - 5: Impone el caudal de agua y la elevación de la superficie del agua. - 6: Impone el caudal de entrada de agua. - 9: Impone una tabla de aforo nivel-caudal de valor único. - 10: Condición de entrada o salida libre\". Las velocidades y elevaciones de la superficie del agua son calculadas por el modelo. - 11: Condición de salida libre\". El modelo calcula las velocidades y elevaciones de la superficie del agua, pero solo se permite el flujo hacia fuera de la malla. - 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 condición 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 condición de contorno asociado.

Note

Si necesita imponer condiciones abiertas en segmentos de contorno adyacentes, hágalo de modo que cada segmento esté separado por un hueco de más de una celda (véase la Figura ). 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 temporal utilizado para 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 fondo.

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, asigna la velocidad a cada celda suponiendo una dirección perpendicular a la línea de contorno, como se muestra:

Tabla de aforo de caudal (BCTYPE 9)¶

Al utilizar una condición nivel-caudal de valor único, el modelo calcula primero el caudal en el contorno, interpola la elevación correspondiente de la superficie del agua a partir de la tabla de aforo e impone ese valor en el siguiente paso de tiempo. 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 las entradas de los valores de la tabla. Consulte la sección correspondiente para conocer el formato del archivo.

Dado que estas condiciones pueden 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. Desafortunadamente, no existe una forma general de seleccionar esa ubicación, por lo que será necesario experimentar numéricamente con el modelo real para obtener una ubicación razonable.

Note

En la mayoría de los ríos con poca pendiente, la relación nivel-caudal está afectada por la histéresis. En otras palabras, la curva nivel-caudal forma un bucle, con caudales mayores en la rama ascendente que en la rama de descenso 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 nivel-caudal en bucle no se consideran en esta versión de OilFlow2D.

Contornos abiertos de condición \"Free\" (BCTYPE 10, 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 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, 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 salida y la entrada de agua, mientras que BCTYPE 11 solo permite el flujo fuera de la malla.

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

Para aplicar esta condición de contorno, el usuario solo proporciona la pendiente del fondo $S_0$. El modelo utilizará $S_0$, Manning's n y el caudal para crear una tabla de aforo. 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 aforo. La tabla se calcula cada 0.05 m (0.16 ft.) desde la elevación mínima del fondo en la sección transversal de salida hasta 50 m (164 ft.) por encima de la elevación máxima del fondo en la sección. Si $S_0=-999$, el modelo calculará la pendiente media del fondo 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 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. 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 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 habitualmente 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 como 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 que 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 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\]

por consiguiente, no se requiere información adicional.

Cuando el estado del flujo es subcrítico tanto en la región de caudal de entrada como en la de salida, 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 el flujo no puede atravesar. Habitualmente, 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 igual a cero.

Determinar si una entrada o salida es supercrítica o subcrítica no es sencillo en una malla 2D. Una caracterización del régimen de flujo del contorno basada en celdas conduce a situaciones complicadas desde los puntos de vista físico y numérico. Por otra parte, las condiciones de contorno físicas o externas suelen referirse a magnitudes medias, como el nivel de la superficie del agua o el caudal total, que deben traducirse en profundidad o velocidad del agua en cada celda según el criterio del usuario. Para tratar 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 cuando 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 $w=Q/S_T$, y definiendo la sección transversal total mojada $S_T$ y la anchura 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 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 total de entrada $Q = Q(t)$ es la función habitual en las simulaciones de inundación, y es importante analizar la mejor forma de imponerlo, ya que afecta a toda la sección transversal de entrada y se trabaja con una representación discreta 2D en celdas computacionales. Pueden presentarse distintos casos.

Casos sencillos¶

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

El caudal total de entrada en el instante $t$, $Q_I(t)$, puede distribuirse a lo largo de la sección transversal de entrada utilizando un caudal constante por unidad de anchura, $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 tenerse en cuenta 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/llenado), y también cambia el número de celdas de contorno implicadas (Figura ).

Al trabajar con secciones de entrada como la de la Figura , un valor uniforme de $q_I$ como el indicado en produce 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 del agua 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 en cada celda de contorno $j$ es variable y se define en función del área total de la sección transversal, $S_T$, y del área transversal individual de la celda, $S_j$:

\[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 el caso general, a un conjunto de nuevas profundidades $h_j^{n+1}$ (Figura ) asociadas generalmente a diferentes niveles de 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 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 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 la anchura de la sección mojada $b_T$. Como muestra la Figura 3, se obtiene un nuevo nivel uniforme de agua en la sección, $d_S$:

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

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

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

se imponga al flujo de entrada un número de Froude máximo, $Fr_{s,max}$. Para ello, manteniendo la anchura de la sección $b_T$, se calcula una nueva área 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$ (Figura ). 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 de integración temporal discreta utilizado, este procedimiento no sigue el criterio de conservación de 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 consiste en imponer directamente el nivel global de la superficie del agua en la sección 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 y 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 de 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 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 OilFlow2D 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 ninguna condición externa. En caso contrario, todas las celdas están 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 cuando la condición de contorno es una curva de aforo de caudal.

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 aparece una subcapa viscosa muy fina que requeriría celdas extremadamente pequeñas para resolverse adecuadamente. OilFlow2D utiliza una condición de deslizamiento en los contornos cerrados y establece un flujo normal nulo a través del contorno, pero permite velocidades tangenciales. OilFlow2D 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 sobre 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 a la pared sólida (Figura ). En otras palabras, si el contorno está cerrado, el borde de contorno asociado $k_{\Gamma}$ es una pared sólida, con una componente de velocidad normal igual a cero. 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¶

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

En OilFlow2D, la malla de celdas triangulares puede cubrir áreas secas y mojadas, y el modelo tratará estas condiciones mediante dos algoritmos distintos según 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 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 la actualización de las variables en ambos lados. La situación denominada seca/mojada se produce 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 de agua en la celda mojada está por debajo del nivel del fondo en la celda seca.
El flujo es subcrítico.

El procedimiento que debe seguirse en ese caso se describe detalladamente en.

El algoritmo de secado y llenado de OilFlow2D 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 sus 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 de celdas secas/mojadas se calculan suponiendo que el borde es un contorno sólido y que las velocidades a ambos lados son cero.
El resto de las contribuciones de los bordes de celda se calcula según 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 mediante 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 próximos 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 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 realmente el contorno celda por celda:

\[\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 instante $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 instante $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:

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

Coeficientes de rugosidad de Manning's n¶

El valor de Manning's n, utilizado habitualmente para estimar las pérdidas de carga en canales y ríos, es una medida global que tiene en cuenta no solo los efectos de la rugosidad del fondo, sino también la fricción interna y 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, las estimaciones de Manning's n 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 término. Varios investigadores han comprobado 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 ello, la selección final de los coeficientes de Manning's n debe ser el resultado de un proceso de calibración en el que los resultados del modelo se ajusten a datos medidos.