[Reseña] Cuando 30 K hacen mentir a Boussinesq — Demou (2018) NOB y la división de presión

Un salto de 30 K en agua es 24 veces el límite de 1.25 K de Gray–Giorgini, y la corrección de densidad se queda cerca del uno por ciento.

Gray y Giorgini dejaron dos números en 1976. En aire a temperatura ambiente, Oberbeck–Boussinesq se queda dentro de un error del 10 % solo hasta 28,6 K. En agua el límite es 1,25 K. El aire de una torre solar recorre cientos de kelvin. El agua de un PV/T se calienta 50 K. Demou, Frantzis y Grigoriadis (2018) sueltan la hipótesis. Las propiedades siguen a la temperatura. Solo el Poisson de presión se queda con coeficiente constante.

Al mover las paredes caliente y fría, con el mismo ΔT, la curva de densidad queda casi plana y la de viscosidad se inclina con fuerza.

T_hot: 333 K (60 °C)T_cold: 293 K (20 °C)

Dashed line marks the reference value at T₀ = 313 K. Δ% is the relative change between hot and cold. For ΔT ≈ 50 K, density moves about 2%, but viscosity changes by tens of percent — the OB approximation hides that.

En 50 K la densidad se mueve cerca de un 2 %. La viscosidad se mueve cerca de un 30 % en el mismo intervalo, y esa diferencia arma la capa límite asimétrica.

Ficha del artículo#

Autores: A. D. Demou, C. Frantzis, D. G. E. Grigoriadis
Revista: International Journal of Heat and Mass Transfer 121 (2018) 1148–1162
DOI: 10.1016/j.ijheatmasstransfer.2018.04.144
Aporte: en convección natural no-OB, las propiedades siguen dependiendo de la temperatura y el Poisson de presión de coeficiente variable pasa a uno de coeficiente constante más una corrección explícita, por la división de Dodd–Ferrante. El solver FFT se queda. La aceleración publicada es de 2,5–5×.

Dónde se rompe Boussinesq en realidad#

La Boussinesq estándar sostiene tres frases a la vez. La densidad varía solo en la gravedad. El resto de propiedades es constante. La disipación viscosa se deja fuera. Si ΔT crece, la viscosidad $\mu(T)$ , la conductividad $k(T)$ y el calor específico $C_p(T)$ se mueven juntos. En un líquido el cambio de densidad puede ser pequeño y el de viscosidad no.

Al conservar todas las propiedades, la ecuación de momento carga un coeficiente variable.

\frac{\partial u_i}{\partial t} + \frac{\partial(u_i u_j)}{\partial x_j} = -\frac{1}{\hat\rho}\frac{\partial P}{\partial x_i} + \frac{1}{\hat\rho}\frac{\Pr}{\sqrt{\mathrm{Ra}}}\frac{\partial}{\partial x_j}\!\left[\hat\mu\!\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right)\right] + \frac{1}{\mathrm{Fr}^2}\delta_{i3}

Los sombreros son propiedades normalizadas con el valor a la temperatura de referencia. La misma variabilidad entra en la energía.

La trampa del Poisson de coeficiente variable#

Un paso fraccionario plano pide a la corrección de presión un campo libre de divergencia. Ese paso es el que llega después de Chorin quita la presión de la ecuación de momento.

\frac{\partial}{\partial x_i}\!\left(\frac{1}{\hat\rho^{n+1}}\frac{\partial P^{n+1}}{\partial x_i}\right) = \frac{1}{\Delta t}\frac{\partial u_i^*}{\partial x_i}.

El factor $1/\hat\rho^{n+1}$ queda dentro de la divergencia. Si varía en espacio y tiempo, la matriz cambia en cada paso. Un solver directo por FFT solo acepta un laplaciano de coeficiente constante. La retirada es multigrid o un solver iterativo. Si el 60–80 % del tiempo de un paso de Navier–Stokes es la presión, esa retirada es la factura.

División Dodd–Ferrante — una línea#

El truco prestado viene del flujo bifásico (Dodd & Ferrante 2014). El gradiente de presión se corta en dos.

\frac{1}{\hat\rho^{n+1}}\frac{\partial P^{n+1}}{\partial x_i} \longrightarrow \frac{1}{\hat\rho_0}\frac{\partial P^{n+1}}{\partial x_i} + \left(\frac{1}{\hat\rho^{n+1}} - \frac{1}{\hat\rho_0}\right)\frac{\partial \widehat P}{\partial x_i}

El primer trozo multiplica la presión nueva por la constante $1/\hat\rho_0$ . El segundo multiplica la presión extrapolada $\widehat P = 2P^n - P^{n-1}$ . Al tomar la divergencia, el Poisson queda con coeficiente constante.

\frac{\partial^2 P^{n+1}}{\partial x_i \partial x_i} = \frac{\hat\rho_0}{\Delta t}\frac{\partial u_i^*}{\partial x_i} + \frac{\partial}{\partial x_i}\!\left[\left(1 - \frac{\hat\rho_0}{\hat\rho^{n+1}}\right)\frac{\partial \widehat P}{\partial x_i}\right]

El laplaciano de la izquierda es el mismo en cada paso. La FFT en dos direcciones deja una pila de Helmholtz en $z$ , resuelta de forma directa. Solo la corrección del lado derecho es explícita.

La magnitud que importa es $(1 - \hat\rho_0/\hat\rho)$ . En flujo bifásico se acerca a 1. En convección NOB de agua es muy pequeña.

Un ciclo en NumPy#

Las secciones 2.2–2.3 quedan en una losa 1D. No hay DNS 2D. Una corrección de presión.

import numpy as np
from numpy.fft import fft, ifft
 
# losa 1D: periódica, una pila de Helmholtz
Nz = 64
Lz = 1.0
dz = Lz / Nz
z = (np.arange(Nz) + 0.5) * dz
 
# densidad dependiente de la temperatura (polinomio del agua, simplificado)
def rho_hat(T_field, T0=313.0):
    dT = T_field - T0
    return 1.0 * (1 - 3.736e-4 * dT - 3.98e-6 * dT**2)
 
# perfil inicial de temperatura: lineal + pequeña perturbacion
T0 = 313.0
dT = 30.0
T_field = T0 - dT/2 + dT * z + 0.1 * np.sin(2*np.pi*z)
rho = rho_hat(T_field)
 
# densidad de referencia: la mas liviana (la mas caliente)
rho_ref = rho.min()
 
# divergencia ficticia (en un solver completo viene de ∇·u* tras el paso fraccional)
div_u_star = 0.01 * np.sin(2*np.pi*z)
dt = 1e-3
 
# presion extrapolada (Ec. 15 del articulo)
P_n_minus_1 = np.zeros(Nz)
P_n = 0.01 * np.cos(2*np.pi*z)
P_hat = 2 * P_n - P_n_minus_1
dP_hat = np.gradient(P_hat, dz)
 
# termino de correccion (segundo del miembro derecho de la Ec. 16)
correction = (1 - rho_ref / rho) * dP_hat
rhs = rho_ref / dt * div_u_star + np.gradient(correction, dz)
 
# Laplaciano de coeficiente constante (FFT, BC periodicas)
k = 2*np.pi * np.fft.fftfreq(Nz, dz)
laplacian_eigenvalues = -(k**2)
laplacian_eigenvalues[0] = 1.0  # quitar el espacio nulo
P_new = np.real(ifft(fft(rhs) / laplacian_eigenvalues))
P_new -= P_new.mean()
 
# tamano del termino de correccion
correction_magnitude = np.max(np.abs(1 - rho_ref / rho))
print(f"max |1 - rho_0/rho|  = {correction_magnitude:.4f}")
print(f"P_new range           = [{P_new.min():.3e}, {P_new.max():.3e}]")

La salida es:

max |1 - rho_0/rho|  = 0.0117
P_new range           = [-1.572e-06, 1.572e-06]

ρ_min / ρ_ref: 0.978ρ_max / ρ_ref: 1.000x-range:

The splitting scheme rewrites the variable Poisson coefficient as constant ρ₀ plus a residual proportional to |1 − ρ₀/ρ|. For water at ΔT = 50 K the correction stays under 3% — that is why a constant-coefficient FFT solver can replace an iterative variable solver without losing accuracy.

El control por defecto es agua a lo largo de 50 K. La corrección se queda dentro del 1 %. Al pasar al rango bifásico, el mismo término se acerca a 0,5. La salida 0,0117 del código con ΔT = 30 K es una ventana más estrecha que ese control de 50 K.

El límite de Gray–Giorgini para agua es 1,25 K. El ΔT = 30 K del código es $30/1.25 = 24$ veces ese límite. La línea del 10 % en aire es 28,6 K, así que 30 K también la pasa por poco. Un salto de 50 K en un PV/T es 40 veces el límite del agua. En esos 50 K la densidad se mueve cerca de un 2 % y la viscosidad cerca de un 30 %.

El coeficiente lineal del polinomio de densidad es $3.736\times 10^{-4}$ . Por 30 K da $3.736\times 10^{-4}\times 30 = 0.0112$ , al lado del $\max|1-\rho_0/\rho| = 0.0117$ impreso. El coeficiente cuadrático $3.98\times 10^{-6}$ por $(30/2)^2 = 225$ es cerca de $9.0\times 10^{-4}$ . Casi todo el 0,0117 es el término lineal. Con $T_0 = 313\,\mathrm{K}$ , $dt = 10^{-3}$ y una amplitud de divergencia de 0,01, el rango de $P$ es $\pm 1.572\times 10^{-6}$ .

El ejemplo bifásico es $\rho_{\min} = 1$ , $\rho_{\max} = 1000$ . Entonces $1 - 1/1000 = 0.999$ . Está a dos órdenes del 0,0117 del agua. Cuando la presión ocupa el 60–80 % de un paso, conservar la FFT es como aparece el 2,5–5× del artículo. Mirar solo $\mathrm{CFL}_{\max} = 0.2$ deja pasar $\mathrm{VSL} = \nu\Delta t/\Delta x^2$ , que cruza 0,03 antes. La teoría admite 0,3, y las estadísticas de orden alto ya se desvían un 3,5 %.

La frase de que la corrección a ΔT = 30 K se queda dentro del 1,2 % señala el mismo 0,0117 impreso. El mismo código con $\rho_{\min} = 1$ y $\rho_{\max} = 1000$ salta a 0,999. Por eso un solver bifásico con esta división sigue pidiendo un $\Delta t$ pequeño.

Con $N_z = 64$ y $L_z = 1$ , $dz = 1/64$ . La presión extrapolada es $\widehat P = 2P^n - P^{n-1}$ . En el listado la amplitud de $P^n$ es 0,01 y $P^{n-1}$ es 0, así que la amplitud de $\widehat P$ es 0,02. El 0,0117 impreso es el coeficiente sobre ese gradiente: la corrección explícita es el 1,17 % de él. La Table 2 recorre $Ra$ de $10^6$ a $10^8$ . Aflojar el techo de VSL de 0,03 a 0,3 en ese rango desplaza las estadísticas de orden alto un 3,5 %.

Tres tropiezos al validar#

La Table 2 (DHC de aire, $Ra = 10^6$ a $10^8$ ) tropieza en tres sitios.

Elegir $\hat\rho_0$ — Dodd–Ferrante usa la densidad mínima, para que $(1-\hat\rho_0/\hat\rho)$ quede al menos en 0 y por debajo de 1. El artículo dice que la densidad a la temperatura de referencia da casi la misma respuesta NOB, pero un ΔT grande gasta el margen de estabilidad.
CFL y VSL a la vez — con Adams–Bashforth explícito, $\mathrm{CFL}_{\max} = 0.2$ no es el primer límite. $\mathrm{VSL} = \nu\Delta t/\Delta x^2 < 0.03$ es el techo práctico. La teoría admite 0,3, y las estadísticas de orden alto se desvían un 3,5 %.
Forma convectiva de Arakawa — momento y energía cinética tienen que conservarse juntos en el discreto. Las diferencias centradas planas filtran energía, y los promedios largos se tuercen a Ra grande.

Fijar o no rho_0 en la densidad mínima#

Dodd–Ferrante fija $\hat\rho_0$ en la densidad mínima para que $(1-\hat\rho_0/\hat\rho)$ quede al menos en 0 y por debajo de 1. El artículo dice que la densidad a la temperatura de referencia basta en NOB. Pasado ΔT = 30 K esa diferencia gasta el margen. El 0,0117 impreso es una losa 1D con la densidad mínima. No es la factura de un DNS 2D.

El 2,5–5× supone que la presión ocupa el 60–80 % del paso y que la FFT tiene una dirección periódica que comer. Explicar la capa límite asimétrica con un 2 % de densidad deja fuera el 30 % de viscosidad. La temperatura en la que esa hipótesis se cierra no es la ventana OB de la convección natural en el número de Rayleigh 1708. En agua ya está cerrada a 1,25 K.