Hyperbolic Grid Generation

The motivation

The governing equations of aerodynamics are actually pretty intuitively simple, but actually solving practical engineering problems is way too fundamentally difficult for all but the simplest of problems. Engineering is about useful approximations though, and the first step of discretization, the process of turning an exact but unsolvable problem into a bunch of immensely easier approximate problems, is to generate a grid on which to solve the problem.

Oh, and of course I should point out that the primary motivation for me actually writing this up and putting it out here on the web is to try to figure out what all the reports I did years ago would have looked like if I’d had a set of tools on hand that allowed me to communicate things interactively. The trend of porting everything to JavaScript gets really old, but as rightfully maligned as it is, it does allow for a unique combination of portability and interactivity that I don’t think are very well utilized yet. So my goal here is to create marginally interesting things but also to refine the tools so that other people can create more interesting things.

So anyway, imagine trying to figure out how much lift and drag an airfoil produces. Unless airflow over the airfoil is laminar (smooth) and maybe even has a clever analytical solution, we’ll need a grid. Since we don’t really care where our grid ends as long as it’s far away, the obvious approach is to start at the airfoil and march a contour outward.

The example below marches a grid concentrically outward, offsetting the grid by a constant distance at each step. It has a big problem though. As the camber (curvature) of the airfoil increases, mesh eventually runs into itself, making fun patterns that unfortunately aren’t very useful.

The problem above is that nothing prevents concentric contours from forming cusps as the grid lines start to intersect. To fix this, we can try forcing areas at risk of collapsing to bulge outward. In other words, cells that get narrow in the around-direction can get longer in the outward-direction. Calling this area-preserving tendency $$\alpha$$, we confirm above that keeping cell area constant as we march outward fixes our grid.

A simple approach like this is what we’re aiming for, but I’ve had to gloss over some numerical surprises to get it to look so straightforward. Let’s write down the equations we need need to see how this all works.

The math

Two conditions define our strategy as we march outward:

1. Grid lines should stay perpendicular.
2. Cell area should stay the same.

Later on, the first condition will make the aerodynamics part much easier since lots of terms simply vanish when the grid is orthogonal, and the second condition prevents the problems we ran into above.

We’ll start by parameterizing the grid. We’ll call distance around the airfoil $$\eta$$ and the outward parameter $$\xi$$ (that’s the Greek letter “xi”, but we liked to call it “tornado” in school).

A single grid cell looks like this:

If we write the sides of this cell as vectors $$\left(\frac{\partial x}{\partial \eta}, \frac{\partial y}{\partial \eta}\right)$$ and $$\left(\frac{\partial x}{\partial \xi}, \frac{\partial y}{\partial \xi}\right)$$, then conditions (1) and (2) state identically that the dot product of the vectors representing the sides is zero (i.e. the sides are perpendicular) and that the cross product of the sides (i.e. the area) is a constant. In equation form,

$$\left\{ \begin{eqnarray} x_\xi x_\eta + y_\xi y_\eta &=& 0 \\ x_\xi y_\eta - x_\eta y_\xi &=& A(\xi,\eta) \end{eqnarray} \right.$$

where the dependence of $$A$$ on $$\xi$$ and $$\eta$$ just means we’ll leave room to take some liberties with our ‘constant’ and maybe tweak the area rule (while maintaining orthogonality, of course) so that we get the best grid we can. Some quick reorganization yields the equations

$$$$\left\{\begin{array}{c} x_\xi\\y_\xi\end{array}\right\} = -\frac{A(\xi,\eta)}{x_\eta^2+y_\eta^2} \left\{\begin{array}{c} -y_\eta\\x_\eta \end{array}\right\}.$$$$

Since we know $$x_\eta$$ and $$y_\eta$$, and $$A(\xi, \eta)$$ is just a function we can make up, this means we have everything we need to just march it outward. Piece of cake.

It wasn’t immediately clear to me why this method would be so unstable. The equations are very similar in form to the first order wave equation $$\frac{\partial f}{\partial t} = c\frac{\partial f}{\partial x}.$$ Lower order time integration and spatial differencing tends to be somewhat unstable for the first order wave equation, but it typically only takes a little bit of effort and some second order methods to at least acheive a stable if not highly accurate solution. It didn’t really dawn on me that fourth order Runge-Kutta integration in $$\xi$$ and sixth-order compact finite differences in $$\eta$$ might not work, but as it turns out, high order methods don’t really help at all.

In search of stability

To understand why the solution is so unstable, we turn to Fourier analysis. Problems that are difficult to analyze in space and time are often much simpler to analyze in the frequency domain. To that end, we’ll seek to represent the relative growth $$y(\xi+\Delta \xi, \eta) / y(\xi,\eta)$$ and $$x(\xi+\Delta \xi, \eta) / x(\xi,\eta)$$ over a single step. If the modulus of these growth factors is greater than one, then with repeated steps, our solution grows exponentially and we observe the instability we saw above.

Taking a Fourier transform in the $$\eta$$-direction, we write our state as $$\left\{\begin{eqnarray} x &=& \hat{x}(\xi) e^{i k_x \eta} \\ y &=& \hat{y}(\xi) e^{i k_y \eta}\end{eqnarray}\right.$$ Since the interesting cases are the ones where $$x$$ and $$y$$ interact, we’ll equate $$k_x$$ and $$k_y$$. And for simplicity, we’ll look at a first-order Euler step from step $$n \, (= \xi)$$ to $$n+1 \, (= \xi + \Delta \xi)$$ and an exact derivative in the $$\eta$$-direction so that the full update for each step outward becomes $$\left\{\begin{eqnarray}\frac{\hat{x}_{n+1} e^{ik\eta} - \hat{x}_n e^{ik\eta}}{\Delta \xi} &=& \phantom{-} \frac{A (i k \hat{y}_n e^{i k \eta})}{|ik\hat{x}_n|^2 + |ik \hat{y}_n|^2} \\ \frac{\hat{y}_{n+1} e^{ik\eta} - \hat{y}_n e^{ik\eta}}{\Delta \xi} &=& - \frac{A (i k \hat{x}_n e^{i k \eta})}{|ik\hat{x}_n|^2 + |ik\hat{y}_n|^2}\\ \end{eqnarray}\right..$$ Simplifying a bit yields $$\left\{ \begin{eqnarray} \hat{x}_{n+1} &=& \hat{x}_n - i \frac{ A \Delta \xi}{k} \frac{\hat{y}_n}{|\hat{x}_n|^2 + |\hat{y}_n|^2} \\ \hat{y}_{n+1} &=& \hat{y}_n + i \frac{ A \Delta \xi}{k} \frac{\hat{x}_n}{|\hat{x}_n|^2 + |\hat{y}_n|^2} \\ \end{eqnarray}\right.$$ To get the relative change over a single step, we divide through by $$\hat{x}_n$$ and $$\hat{y}_n$$ to get $$\left\{ \begin{eqnarray} \frac{\hat{x}_{n+1}}{\hat{x}_n} &=& 1 - i \left[\frac{ A \Delta \xi}{k(|\hat{x}_n|^2 + |\hat{y}_n|^2)} \right] \frac{\hat{y}_n}{\hat{x}_n} \\ \frac{\hat{y}_{n+1}}{\hat{y}_n} &=& 1 + i \left[\frac{ A \Delta \xi}{k (|\hat{x}_n|^2 + |\hat{y}_n|^2)}\right] \frac{\hat{x}_n}{\hat{y}_n} \\ \end{eqnarray}\right.$$ As a sanity check, note that if $$\Delta \xi = 0$$, i.e. if the step size outward is zero, then $$\hat{x}_{n+1}/\hat{x}_n$$ = $$\hat{y}_{n+1}/\hat{y}_n = 1$$ and the solution remains unchanged. Since the term in brackets is a real number that can be made arbitrarily small by taking a smaller step, we’ll just call it $$\epsilon$$ to get $$\left\{ \begin{eqnarray} \frac{\hat{x}_{n+1}}{\hat{x}_n} &=& 1 - i \epsilon \frac{\hat{y}_n}{\hat{x}_n} \\ \frac{\hat{y}_{n+1}}{\hat{y}_n} &=& 1 + i \epsilon \frac{\hat{x}_n}{\hat{y}_n} \\ \end{eqnarray}\right.$$ If we let $$\hat{y}_n/\hat{x}_n = re^{i\theta} = r\cos(\theta) + ir\sin(\theta)$$, then we see at last that $$\left\{ \begin{eqnarray} \frac{\hat{x}_{n+1}}{\hat{x}_n} &=& 1 - i r \epsilon \cos(\theta) + r \epsilon \sin(\theta) \\ \frac{\hat{y}_{n+1}}{\hat{y}_n} &=& 1 + \frac{i\epsilon}{r} \cos(\theta) + \frac{\epsilon}{r} \sin(\theta) \\ \end{eqnarray}\right.$$ defines our step. $$\epsilon$$ is a small number and $$r$$ and $$\theta$$ come from the relation of $$x$$ to $$y$$ so that we don’t really have any control over them. To make a long story short, plotting the growth factor for a single $$\xi$$-step on the unit circle below shows that there’s no $$\Delta \xi$$ small enough to keep things stable.

The best we can do and indeed what people seem to do is to add diffusion. In terms of stability, diffusion pulls things toward the origin of the complex plane and back inside the region of stability (the unit circle). The downside is that when we smooth out the bumps with diffusion, we also lose a bit of the orthogonality we were aiming for.

With some diffusion to keep things stable, the final result is below. In lieu of all of these caveats, it works reasonably well! In a future writeup, I’ll try to use this to do some actual fluid dynamics. You can find the code on github.