Wrapping a Rectangular Grid Around a Wing — Curvilinear Coordinate Transforms and Grid Metrics

An engineer wants to solve the flow around a wing. The solver in hand only understands Cartesian grids. But the wing is curved. No matter how finely you chop the rectangular cells, a curved surface still comes out looking like a staircase. The fix is not to bend the grid but to bend space itself. You map the warped physical grid back onto a perfectly straight logical grid, then solve the equations there. This post walks through the coordinate transform and grid metrics that make that round trip possible. By the end, you will see why the Jacobian flags grid folding, and what it takes for a uniform flow to stay uniform on a curved grid — verified directly in numpy.

Raise shear and amplitude on the grid mapper and the cells twist. The Jacobian heatmap paints each cell by $J/J_0$.

Push amplitude near 0.3 and wave number to 4 or higher and red cells show up. $J$ in those cells is negative. Cross that line and a metric denominator passes through zero.

Physical and Logical Domains#

Curvilinear methods start with two domains. The physical domain is where the warped grid wrapping the wing lives. Its coordinates are $(x, y)$. The logical domain is that same grid flattened into a unit square. Its coordinates are $(\xi, \eta)$.

The mapping is the function that links them.

Here $\xi$ is the index along one grid direction and $\eta$ along the other. In the logical domain every cell is a square of side $\Delta\xi = \Delta\eta = 1/N$. The solver computes differences only on this straight grid. All of the warping is absorbed into the mapping function.

Metrics: The Exchange Rate the Grid Builds#

The trouble is that the equation we want to solve is written in derivatives with respect to $x, y$. The grid moves along $\xi, \eta$. The exchange-rate table that connects them is the set of grid metrics.

Start from the chain rule.

Terms like $\xi_x = \partial\xi/\partial x$ are the metric coefficients. What we actually hold is the reverse: $x_\xi, x_\eta, y_\xi, y_\eta$, because differentiating the mapping with respect to $\xi, \eta$ gives them directly. The two are tied through the inverse-matrix relation.

Each symbol is a first derivative of the mapping, and $J$ is the Jacobian from the next section. In short, a metric is "the easy-to-get $x_\xi$ family divided by the Jacobian."

The Jacobian and Grid Folding#

The Jacobian determinant measures how much the mapping stretches or shrinks area.

Here $J$ is the ratio between a unit area in the logical domain and the area it occupies in the physical domain. If $J > 0$ the mapping preserves orientation. If $J = 0$ a cell collapses to a point or a line. If $J < 0$ the cell flips over. That is grid folding. On a folded grid the metric denominators cross zero and blow up, and the solver spits out NaN.

$J = x_\xi y_\eta - x_\eta y_\xi$ is the cross product of the two cell edges. Same term as the cell area measured by a diagonal cross product. A cell that went red on the heatmap above is where $J < 0$.

The Transformed Conservation Law#

Now move the equation into the logical domain. Take a 2D conservation law.

$Q$ is the conserved quantity and $F, G$ are the fluxes in $x, y$. Substituting the metrics and multiplying by $J$ yields the strong conservation form.

The transformed fluxes are

The key point is that products like $J\xi_x = y_\eta$ and $J\xi_y = -x_\eta$ fold back into plain derivatives of the mapping. So the transformed flux never has to divide by a metric and then multiply it back. Keep this form and mass and momentum stay exactly conserved even on a curved grid.

Computing and Verifying Metrics in Code#

Port the theory to numpy. Define a wavy physical grid, recover the metrics and Jacobian with central differences, and check whether a uniform flux is really conserved.

import numpy as np
 
def body_fitted_map(xi, eta, amp=0.15, shear=0.3, wavek=2.0):
    """Map logical coords (xi, eta) in [0,1]^2 to physical coords (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):
    """Recover Jacobian J and inverse metrics by central differences."""
    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        # Jacobian determinant
    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)
 
# logical grid
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("grid:", "folded (J<=0)" if J.min() <= 0 else "valid (J>0)")

Next is the freestream-preservation test. Feed in constant fluxes $F, G$ and the residual should be zero. Writing the transformed flux directly as $J\xi_x = y_\eta$ lets the mixed central differences commute, so the residual drops to round-off.

def freestream_residual(x, y, dxi, deta, F=1.0, G=0.5):
    """Freestream residual for constant flux (F,G). Should be zero."""
    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"freestream residual max|R| (interior) = {np.abs(R[1:-1, 1:-1]).max():.2e}")
# e.g. J min/max = 0.62 / 1.38, freestream residual max|R| ~ 1e-16

A common mistake is solving in the non-conservative form, dividing the metric by $J$ and then multiplying it back onto the flux. In three dimensions that mismatch grows and seeds a spurious source in a uniform flow. The symmetric conservative metrics of Thomas and Lombard (1979) prevent it. In 2D, writing $\hat F$ directly in terms of $y_\eta, x_\eta$ as above is enough.

Start from the numbers in the code comment. With amp=0.15, shear=0.3, wavek=2, and N=41, $J$ runs from 0.62 to 1.38. The max is $1.38/0.62$, about 2.23 times the min. Cell area wobbles by more than a factor of two and every $J$ stays positive. Nothing has folded. On that same mesh the freestream residual $\max|R|$ is about 1e-16, for $F=1$ and $G=0.5$.

Follow the residual through the formula. $\hat F = F y_\eta - G x_\eta$ and $\hat G = -F y_\xi + G x_\xi$. Add the $\xi$ difference to the $\eta$ difference. As long as the central operator swaps mixed derivatives, the $\xi$ derivative of $y_\eta$ cancels the $\eta$ derivative of $y_\xi$. $F$ and $G$ are constant, so they factor out. What is left is the difference of two mixed derivatives of the map. The same operator makes that difference zero. Divide $\xi_x$ by $J$ and multiply by $J$ again and, on paper, you are back at $y_\eta$. On a mesh the division and the product can touch $J$ at different places. That gap is a fake source in a uniform stream. In 2D, writing $\hat F$ straight from $y_\eta$ and $x_\eta$ never opens it. The 1e-16 in the comment is the cancellation that survived.

$\beta$ produces a number from the same formula. At the wall, $\eta=0$, $dy/d\eta = \beta/(e^\beta-1)$. At $\eta=1$ it is $\beta e^\beta/(e^\beta-1)$. The ratio of the two slopes is $e^\beta$. At $\beta=4$ that is about 54.6. The "tens of times" in the wall section is this ratio. The limit $\beta \to 0$ is a slope of 1, a uniform mesh.

Clustering Cells Near the Wall#

To resolve a boundary layer you have to pack cells near the wall. That is the same coordinate transform in one dimension. Stretch a uniform $\eta$ with an exponential to build the physical coordinate $y$.

Here $\beta$ is the clustering parameter. The 1D Jacobian is $dy/d\eta = \beta e^{\beta\eta}/(e^\beta - 1)$, and the grid tightens wherever that value is small.

Play with $\beta$ in the simulation below.

Raise $\beta$ from 0 to 4 and the first spacing near the wall ($y=0$) shrinks sharply, opening the ratio against the last spacing to tens of times. The same cell count resolves the boundary layer far better.

Divide by J and multiply back, and a uniform flow grows a source#

A cell with negative $J$ does not get a flux. On the default wave, amp=0.15, shear=0.3, wavek=2, N=41, $J$ sits between 0.62 and 1.38. Nothing is folded. If the generator paints even one red cell, that mesh does not go into the solver. NaN comes out of that denominator.

One form stays. Write $\hat F$ from $y_\eta$ and $x_\eta$. Dividing by $J$ and multiplying back opens a gap even in 2D, and in 3D it is the spurious source Thomas–Lombard blocked. Turning the stiffness itself with the same metric is a separate case, in the coordinate transform of the constitutive tensor.