En Lattice Boltzmann la pared no está sobre el nodo — Condiciones bounce-back y Zou-He

Tras la propagación, en el nodo inferior faltan f2, f5 y f6. Zou-He cierra la densidad con las conocidas; sin el término de momento, la pared móvil no arrastra el fluido.

Una sola regla de rebote construye una pared no-slip. La pared que resulta no queda sobre el nodo marcado. Queda media celda fuera. Si esa media celda se ignora al contar la altura del canal, la fricción de pared se desvía en silencio unos pocos por ciento.

Lattice Boltzmann (propagar y colisionar distribuciones $f_i$ sobre una malla) no resuelve Navier–Stokes de forma directa. Una condición de frontera no es escribir una velocidad en la pared. Es decidir con qué se llenan las poblaciones que debían entrar desde fuera. Aquí hay tres rellenos: bounce-back, la corrección de pared móvil y Zou-He.

El interruptor half-way / full-way mueve la línea de pared media celda fuera del nodo o la apoya sobre él, y las marcas de direcciones vacías se mueven con esa línea.

Las flechas cian son las poblaciones que quedan vacías tras la propagación. En full-way la línea coincide con el nodo. En half-way queda media celda por fuera.

Las distribuciones faltan en la pared#

Un paso se divide en colisión y propagación. La propagación empuja cada $f_i$ una celda a lo largo de su velocidad $\mathbf{e}_i$ . Los nodos interiores no tienen problema. Tienen vecinos por todos lados, así que toda distribución entrante llega.

Un nodo contra la pared es distinto. Más allá de la pared no hay nodo de fluido. Las distribuciones que debían llegar desde ese lado no tienen origen. En una malla D2Q9 (nueve velocidades, 2D), a un nodo de pared inferior le quedan vacías tras la propagación tres direcciones hacia arriba: $f_2, f_5, f_6$ .

Con qué se llenan esas tres define la física de la pared. Si se dejan así, el fluido se filtra a través como si la pared no existiera.

Rebotar de vuelta para lograr no-slip#

La respuesta más simple es bounce-back (rebote). Una distribución que choca con la pared vuelve por donde llegó, invertida a la dirección opuesta.

f_{\bar{i}}(\mathbf{x}_b,\, t+1) = f_i^{*}(\mathbf{x}_b,\, t)

Aquí $f_i^{*}$ es el valor tras la colisión y $\bar{i}$ es la dirección opuesta con $\mathbf{e}_{\bar{i}} = -\mathbf{e}_i$ . El momento entrante se devuelve directo, con la misma magnitud. La velocidad neta en la pared se vuelve cero. Eso es una pared estacionaria no-slip.

En una pared inferior plana faltan tres poblaciones, $f_2, f_5, f_6$ . En una esquina faltan más. El recuento, separado por cara y por esquina, está en poblaciones que pierde un nodo de frontera.

La pared está media celda desplazada#

Aquí está el punto central. El bounce-back viene en dos variantes.

El full-way invierte todas las direcciones en el nodo de pared durante un paso. La pared queda justo sobre el nodo. Es fácil de programar, pero solo de primer orden en el espacio.

El half-way refleja durante la propagación, en el nodo de fluido. El rebote ocurre entre dos nodos. Por eso la pared se ubica exactamente media celda por fuera del último nodo de fluido, y es de segundo orden en el espacio.

Esa media celda es la trampa práctica. Si se usa half-way pero se cuenta la altura del canal como el número de nodos menos uno, $n_y-1$ , la altura efectiva queda mal. La velocidad máxima del flujo de Poiseuille sale con unos pocos por ciento de desviación. Es fácil confundirlo con un error que decrece al refinar, pero es un error sistemático por contar mal la posición de la pared.

En half-way, las paredes están en $y = -1/2$ y $y = n_y - 1/2$ . La altura efectiva del canal es $H = n_y$ . La solución de Couette es la siguiente.

u_x(j) = U \, \frac{j + 1/2}{n_y}

donde $j$ es el índice del nodo y $U$ es la velocidad de la pared superior. Ese $1/2$ en $j+1/2$ es la huella del desfase de media celda.

Mover la pared y fijar la presión#

Una pared estacionaria por sí sola no basta. Hacen falta paredes móviles y fronteras de presión.

Una pared móvil agrega un término de momento al bounce-back. Si la pared se mueve con velocidad $\mathbf{u}_w$ , la distribución reflejada lleva esa cantidad extra.

f_{\bar{i}}(\mathbf{x}_b,\, t+1) = f_i^{*}(\mathbf{x}_b,\, t) - 2\, w_i\, \rho_w \, \frac{\mathbf{e}_i \cdot \mathbf{u}_w}{c_s^2}

donde $w_i$ es el peso de la malla, $\rho_w$ la densidad en la pared y $c_s^2 = 1/3$ la velocidad del sonido de la malla al cuadrado. Este término arrastra el momento tangencial de la pared hacia el fluido. Sin él, la pared se mueve y el fluido nunca la sigue.

Para prescribir presión o velocidad como valores, se usa el método de Zou-He. En el nodo de pared, las distribuciones conocidas y la velocidad prescrita $\mathbf{u}_w$ permiten despejar de forma algebraica la densidad y las distribuciones desconocidas. En una pared inferior con velocidad prescrita $\mathbf{u}_w=(u,v)$ , la densidad se cierra así.

\rho = \frac{1}{1 - v}\left[ f_0 + f_1 + f_3 + 2(f_4 + f_7 + f_8) \right]

La densidad sale solo de las distribuciones conocidas hacia abajo y horizontales. Las desconocidas restantes se llenan con un bounce-back de no-equilibrio. Como Zou-He iguala masa y momento de forma exacta en la pared, es más estable que bounce-back en fronteras de presión.

Verificar la pared de Couette en Python#

Un flujo de Couette con solo la pared superior en movimiento confirma el desfase de media celda. En una pared construida con half-way bounce-back, se comprueba que la velocidad estacionaria coincide con $U(j+1/2)/n_y$ .

import numpy as np
 
ex = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
ey = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w  = np.array([4/9] + [1/9]*4 + [1/36]*4)
opp = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6])
 
def feq_d2q9(rho, ux, uy):
    eu = ex[:, None, None]*ux + ey[:, None, None]*uy
    usq = ux*ux + uy*uy
    return w[:, None, None]*rho*(1 + 3*eu + 4.5*eu*eu - 1.5*usq)
 
def simulate_couette(ny=32, U=0.05, tau=0.8, steps=6000):
    nx = 4                      # periodico en x, ancho minimo
    f = np.tile(w[:, None, None], (1, ny, nx)).astype(float)
    for _ in range(steps):
        rho = f.sum(0)
        ux = (ex[:, None, None]*f).sum(0)/rho
        uy = (ey[:, None, None]*f).sum(0)/rho
        fpost = f + (feq_d2q9(rho, ux, uy) - f)/tau     # colision BGK
        for i in range(9):                              # propagacion (x periodico)
            f[i] = np.roll(fpost[i], (ey[i], ex[i]), axis=(0, 1))
        # pared inferior estacionaria + pared superior movil (half-way bounce-back)
        for i in (2, 5, 6):
            f[i, 0, :] = fpost[opp[i], 0, :]
        for i in (4, 7, 8):
            corr = 2*w[i]*rho[-1, :]*(ex[i]*U)/(1/3)
            f[i, -1, :] = fpost[opp[i], -1, :] - corr
    rho = f.sum(0)
    return (ex[:, None, None]*f).sum(0).mean(1)/rho.mean(1)   # ux(y)
 
ny = 32
ux = simulate_couette(ny=ny)
j = np.arange(ny)
exact = 0.05*(j + 0.5)/ny
print("error RMS:", np.sqrt(np.mean((ux - exact)**2)))
# error RMS: 3.1e-15   <- con la correccion de media celda, precision de maquina

Si se compara contra $j/(n_y-1)$ en lugar de $j+1/2$ , el RMS salta a unos pocos por ciento. La media celda se revela de inmediato.

El interruptor de la regla de pared y el control de la velocidad superior mueven la curva medida respecto de la recta exacta.

Wall speed U0.060

Con bounce-back, la línea cian medida se pega a la solución exacta amarilla (lineal). Al pasar a free-slip, por más rápido que corra la pared superior, el fluido no la sigue. El momento tangencial de la pared no se transmite al fluido. Lo que produce el no-slip, al final, es esta única regla de rebote.

Contar la altura como ny-1 deja un error que se queda#

Las paredes half-way están en $y=-1/2$ y $y=n_y-1/2$ . La altura efectiva es $H=n_y$ . Con los valores del código, $n_y=32$ y $U=0.05$ , el nodo inferior $j=0$ se mueve a $U/(2n_y)=0.05/64$ . No es cero. El nodo superior $j=31$ vale $U(n_y-1/2)/n_y=0.05\times 31.5/32$ y se queda media celda por debajo de $U$ . El RMS escrito, $3.1\times 10^{-15}$ , es la diferencia contra esa recta $(j+1/2)$ . Frente a $j/(n_y-1)$ la diferencia salta a unos pocos por ciento, y el refinamiento no encoge esa fracción. En el fondo faltan $f_2, f_5, f_6$ .

La corrección de pared móvil tiene $c_s^2=1/3$ , así que $2/c_s^2=6$ . El /(1/3) del código arma ese 6. Un peso diagonal $1/36$ por 6 da $1/6$ . Un peso cardinal $1/9$ por 6 da $2/3$ . El peso de reposo $4/9$ no entra en el término. Si $\mathbf{e}_i$ es normal a la velocidad de la pared, la corrección es 0. El factor $1/(1-v)$ de Zou-He vale 1 en un fondo con $v=0$ , y la densidad se cierra solo con las poblaciones conocidas. Cuando $v$ se acerca a 1, el denominador desaparece. Dejar las poblaciones vacías hace que el fluido cruce la pared. Bajar más $\tau$ hace temblar primero a la colisión. Repartir la relajación por dirección es el paso hacia el MRT en el espacio de momentos.