Envolver un ala con una malla rectangular — transformación de coordenadas curvilíneas y métricas de malla

Un ingeniero quiere resolver el flujo alrededor de un ala. El solver que tiene a mano solo entiende mallas cartesianas. Pero el ala es curva. Por más fino que se piquen las celdas rectangulares, una superficie curva sigue saliendo en forma de escalera. La salida no es doblar la malla, sino doblar el espacio mismo. Se mapea la malla física deformada de vuelta a una malla lógica perfectamente recta, y allí se resuelven las ecuaciones. Este artículo recorre la transformación de coordenadas y las métricas de malla que hacen posible ese viaje de ida y vuelta. Al final quedará claro por qué el jacobiano detecta el plegado de la malla, y qué hace falta para que un flujo uniforme siga uniforme sobre una malla curva, verificado directamente en numpy.

Al subir shear y amplitude en el mapeador, las celdas se tuercen. El mapa de calor del jacobiano pinta cada celda según $J/J_0$.

Con amplitude cerca de 0.3 y wave number en 4 o más aparecen celdas rojas. En esas celdas $J$ es negativo. Si el generador cruza esa línea, un denominador métrico pasa por cero.

Dominio físico y dominio lógico#

Los métodos curvilíneos parten de dos dominios. El dominio físico es donde vive la malla deformada que envuelve el ala. Sus coordenadas son $(x, y)$. El dominio lógico es esa misma malla aplanada en un cuadrado unitario. Sus coordenadas son $(\xi, \eta)$.

El mapeo es la función que los enlaza.

Aquí $\xi$ corresponde al índice a lo largo de una dirección de la malla y $\eta$ a la otra. En el dominio lógico cada celda es un cuadrado de lado $\Delta\xi = \Delta\eta = 1/N$. El solver calcula diferencias solo sobre esta malla recta. Toda la deformación queda absorbida dentro de la función de mapeo.

Métricas: el tipo de cambio que construye la malla#

El problema es que la ecuación que se quiere resolver está escrita en derivadas respecto a $x, y$. La malla se mueve a lo largo de $\xi, \eta$. La tabla de tipos de cambio que conecta ambos son las métricas de malla.

Se parte de la regla de la cadena.

Términos como $\xi_x = \partial\xi/\partial x$ son los coeficientes métricos. Pero lo que realmente se tiene es el sentido inverso, es decir $x_\xi, x_\eta, y_\xi, y_\eta$, porque derivar el mapeo respecto a $\xi, \eta$ los da de forma directa. Ambos quedan ligados por la relación de matriz inversa.

Cada símbolo es una derivada primera del mapeo, y $J$ es el jacobiano de la sección siguiente. En resumen, una métrica es "la familia $x_\xi$, fácil de obtener, dividida por el jacobiano".

El jacobiano y el plegado de la malla#

El determinante jacobiano mide cuánto estira o encoge el área el mapeo.

Aquí $J$ es la razón entre un área unitaria del dominio lógico y el área que ocupa en el dominio físico. Si $J > 0$ el mapeo conserva la orientación. Si $J = 0$ una celda colapsa a un punto o a una línea. Si $J < 0$ la celda se da vuelta. Eso es el plegado de la malla (grid folding). En una malla plegada los denominadores métricos cruzan el cero y divergen, y el solver arroja NaN.

$J = x_\xi y_\eta - x_\eta y_\xi$ es el producto cruzado de los dos lados de la celda. Es el mismo término que el área de celda medida por un producto cruzado de diagonales. La celda que se puso roja en el mapa de arriba es donde $J < 0$.

La ley de conservación transformada#

Ahora se lleva la ecuación al dominio lógico. Tomemos una ley de conservación 2D.

$Q$ es la cantidad conservada y $F, G$ son los flujos en $x, y$. Al sustituir las métricas y multiplicar por $J$ se obtiene la forma de conservación fuerte (strong conservation form).

Los flujos transformados son

La clave es que productos como $J\xi_x = y_\eta$ y $J\xi_y = -x_\eta$ se reducen a simples derivadas del mapeo. Así el flujo transformado nunca tiene que dividir por una métrica para luego volver a multiplicarla. Si se respeta esta forma, la masa y el momento se conservan exactamente incluso sobre una malla curva.

Calcular y verificar las métricas en código#

Llevemos la teoría a numpy. Definimos una malla física ondulada, recuperamos las métricas y el jacobiano con diferencias centradas, y comprobamos si un flujo uniforme realmente se conserva.

import numpy as np
 
def body_fitted_map(xi, eta, amp=0.15, shear=0.3, wavek=2.0):
    """Mapea coords lógicas (xi, eta) en [0,1]^2 a coords físicas (x, y)."""
    x = xi + shear * eta + amp * np.sin(np.pi * wavek * eta)
    y = eta + amp * np.sin(np.pi * wavek * xi)
    return x, y
 
def compute_metrics(x, y, dxi, deta):
    """Recupera el jacobiano J y las métricas inversas por diferencias centradas."""
    x_xi  = np.gradient(x, dxi,  axis=1)   # dx/dxi
    x_eta = np.gradient(x, deta, axis=0)   # dx/deta
    y_xi  = np.gradient(y, dxi,  axis=1)
    y_eta = np.gradient(y, deta, axis=0)
    J = x_xi * y_eta - x_eta * y_xi        # determinante jacobiano
    xi_x, xi_y  =  y_eta / J, -x_eta / J
    eta_x, eta_y = -y_xi / J,  x_xi / J
    return J, (xi_x, xi_y, eta_x, eta_y)
 
# malla lógica
N = 41
xi  = np.linspace(0, 1, N)
eta = np.linspace(0, 1, N)
dxi, deta = xi[1] - xi[0], eta[1] - eta[0]
XI, ETA = np.meshgrid(xi, eta)             # axis0=eta, axis1=xi
 
X, Y = body_fitted_map(XI, ETA, amp=0.15, shear=0.3, wavek=2.0)
J, _ = compute_metrics(X, Y, dxi, deta)
 
print(f"J min/max = {J.min():.4f} / {J.max():.4f}")
print("malla:", "plegada (J<=0)" if J.min() <= 0 else "valida (J>0)")

Sigue la prueba de preservación del flujo libre (freestream preservation). Al introducir flujos constantes $F, G$ el residuo debe ser cero. Escribir el flujo transformado directamente como $J\xi_x = y_\eta$ hace que las derivadas mixtas centradas conmuten, así que el residuo cae al nivel del redondeo.

def freestream_residual(x, y, dxi, deta, F=1.0, G=0.5):
    """Residuo de flujo libre para flujo constante (F,G). Debe ser cero."""
    x_xi  = np.gradient(x, dxi,  axis=1)
    x_eta = np.gradient(x, deta, axis=0)
    y_xi  = np.gradient(y, dxi,  axis=1)
    y_eta = np.gradient(y, deta, axis=0)
    Fhat =  F * y_eta - G * x_eta          # = J(xi_x F + xi_y G)
    Ghat = -F * y_xi  + G * x_xi           # = J(eta_x F + eta_y G)
    dFhat = np.gradient(Fhat, dxi,  axis=1)
    dGhat = np.gradient(Ghat, deta, axis=0)
    return dFhat + dGhat
 
R = freestream_residual(X, Y, dxi, deta)
print(f"residuo de flujo libre max|R| (interior) = {np.abs(R[1:-1, 1:-1]).max():.2e}")
# ej.) J min/max = 0.62 / 1.38, residuo de flujo libre max|R| ~ 1e-16

Un error frecuente es resolver en la forma no conservativa, dividiendo la métrica por $J$ y luego volviéndola a multiplicar sobre el flujo. En tres dimensiones ese desajuste crece y siembra una fuente espuria en un flujo uniforme. Las métricas conservativas simétricas de Thomas y Lombard (1979) lo evitan. En 2D basta con escribir $\hat F$ directamente en términos de $y_\eta, x_\eta$ como arriba.

Los números salen del comentario del código. Con amp=0.15, shear=0.3, wavek=2 y N=41, $J$ va de 0.62 a 1.38. El máximo es $1.38/0.62$, cerca de 2.23 veces el mínimo. El área de celda oscila más del doble y todo $J$ sigue positivo. No hay celda plegada. En esa misma malla el residuo de flujo libre $\max|R|$ es del orden de 1e-16, con $F=1$ y $G=0.5$.

El área se mueve y el residuo se queda en el redondeo. La fórmula dice por qué. $\hat F = F y_\eta - G x_\eta$ y $\hat G = -F y_\xi + G x_\xi$. Al sumar la diferencia en $\xi$ y la diferencia en $\eta$, si el operador centrado intercambia las derivadas mixtas, la derivada en $\xi$ de $y_\eta$ cancela la derivada en $\eta$ de $y_\xi$. $F$ y $G$ son constantes y salen fuera. Queda la diferencia de dos derivadas mixtas del mapeo. El mismo operador deja esa diferencia en cero. Dividir $\xi_x$ por $J$ y volver a multiplicar por $J$ regresa, en el papel, a $y_\eta$. En la malla la división y el producto pueden tocar $J$ en sitios distintos. Esa rendija es una fuente falsa en un flujo uniforme. En 2D, escribir $\hat F$ directo desde $y_\eta$ y $x_\eta$ no la abre. El 1e-16 del comentario es la cancelación que quedó.

$\beta$ también da un número con la misma fórmula. En la pared, $\eta=0$, $dy/d\eta = \beta/(e^\beta-1)$. En $\eta=1$ es $\beta e^\beta/(e^\beta-1)$. El cociente de las dos pendientes es $e^\beta$. Con $\beta=4$ eso es cerca de 54.6. Las decenas de veces del tramo de pared son este cociente. El límite $\beta \to 0$ es una pendiente de 1, una malla uniforme.

Agrupar celdas cerca de la pared#

Para resolver una capa límite hay que apretar las celdas cerca de la pared. Eso es la misma transformación de coordenadas en una dimensión. Se estira un $\eta$ uniforme con una exponencial para construir la coordenada física $y$.

Aquí $\beta$ es el parámetro de agrupamiento. El jacobiano 1D es $dy/d\eta = \beta e^{\beta\eta}/(e^\beta - 1)$, y la malla se aprieta donde ese valor es pequeño.

Conviene jugar con $\beta$ en la simulación de abajo.

Al subir $\beta$ de 0 a 4, el primer espaciado cerca de la pared ($y=0$) se encoge de golpe, abriendo la razón frente al último espaciado a decenas de veces. El mismo número de celdas resuelve la capa límite mucho mejor.

Dividir por J y volver a multiplicar siembra una fuente en el flujo uniforme#

Una celda con $J$ negativo no recibe flujo. En la onda por defecto, amp=0.15, shear=0.3, wavek=2, N=41, $J$ queda entre 0.62 y 1.38. No hay plegado. Si el generador pinta aunque sea una celda roja, esa malla no entra al solver. El NaN sale de ese denominador.

La forma que se queda es una. $\hat F$ se escribe con $y_\eta$ y $x_\eta$. Dividir por $J$ y volver a multiplicar abre una rendija incluso en 2D, y en 3D es la fuente espuria que Thomas–Lombard cerraba. Girar la rigidez misma con la misma métrica es otro caso, en la transformación de coordenadas del tensor constitutivo.