The Kim, Moin, and Moser (KMM) Velocity-Vorticity Spectral Formulation

A plane channel is the simplest wall-bounded turbulent flow we can simulate exactly. Two parallel walls, a fluid driven along the gap between them, and every eddy in that gap resolved on the grid rather than modelled. The two directions along the wall are statistically homogeneous, so we take them periodic: x downstream, and z across the span. The remaining direction, y, runs across the gap, from minus one at the lower wall to plus one at the upper. At each wall the fluid sticks: all three velocity components vanish there. That single condition, no-slip on a solid boundary, is the source of every difficulty in this lecture. Before deriving anything, let me put the destination on the board. Kim, Moin and Moser march two scalar fields. The first is the velocity component normal to the walls, and its equation is fourth order in space. The second is the vorticity component normal to the walls. Its own equation is a diffusion equation with a source term, and it is only second order in space. And look at what is missing. There is no pressure in either line. Two scalars have replaced three velocity components and a pressure, and every boundary condition the walls hand us will turn out to be a condition on something we are genuinely solving for. Here is the system we normally write down. Momentum for each of the three velocity components, with the pressure gradient sitting in it, and incompressible continuity beside it. Count them: four unknowns and four equations. But look at where the pressure sits. It enters only through its gradient in the momentum equations, and continuity does not contain it at all. Take the divergence of the momentum equations, use continuity, and the time derivative drops out. What is left is an elliptic problem: a Poisson equation whose right hand side is quadratic in the velocity. So the pressure is not marched at all. At every time step it is determined, everywhere at once, by the velocity field everywhere at once. And that is where the trouble starts. The first difficulty is the one we have just met. The pressure has no equation to march, so it has to be recovered from an elliptic solve at every step, over the whole channel at once. The second is that the velocity has to come out divergence free. Solve the momentum equations on their own and it will not. So velocity and pressure are two halves of one solve, not two solves in sequence. The third is the sharpest. The walls tell us the velocity there: it is zero. They tell us nothing at all about the pressure. And yet an elliptic problem demands a condition on every boundary, so whatever we supply is a numerical invention, sitting exactly where the interesting physics lives. Kim, Moin and Moser take the obvious way out. If the pressure is the problem, remove it from the equations altogether. That takes one vector identity and two derivatives, and it is the next thing we do.

Vorticity is the curl of the velocity, and the component we are going to need points along the wall-normal direction. Here is a wall-parallel plane, with the horizontal velocity drawn on it. This particular field circulates. It turns about the wall-normal axis, which is perpendicular to the plane you are looking at, and that turning is exactly what omega y measures. In components, only the two wall-parallel derivatives appear: the rate at which the downstream velocity varies across the span, minus the rate at which the spanwise velocity varies downstream. Now the derivation. Collect every nonlinear term into one symbol, N, so that the momentum equation reads: rate of change equals N, minus the pressure gradient, plus viscous diffusion. Take the curl of that equation. The pressure term is the gradient of a scalar, and the curl of a gradient is identically zero, so it simply disappears. That is the whole trick, and the rest of this lecture is bookkeeping around it. What comes back is the vorticity transport equation. Take its wall-normal component and the second of our two governing equations is already there: omega y diffuses, and it is stirred by a source built out of the nonlinear terms alone. The velocity equation takes more work. Start from the same momentum equation and take its divergence. Continuity kills the time derivative, and what is left is a Poisson equation for the pressure with a right hand side made only of velocities. Now take the wall-normal component of the momentum equation and hit the whole line with a Laplacian. Three of the four terms are harmless. The pressure term becomes the Laplacian of a y derivative, which is the y derivative of a Laplacian, and we have just worked out what the Laplacian of the pressure is. Substitute it in and the pressure is gone for good. Rearranged, this is the fourth-order equation: the Laplacian of v is what evolves, viscosity acts through a fourth derivative, and the source h collects the nonlinear terms. There is a neater way to say what we just did. Minus the Laplacian of a divergence-free velocity is the curl of its vorticity, so applying the Laplacian to the wall-normal momentum equation is taking the curl twice. The fourth-order operator is what a double curl looks like. So here they are together. A fourth-order equation for the wall-normal velocity, and a second-order equation for the wall-normal vorticity. They are coupled, but only through the source terms h, which are quadratic in the velocity and will be evaluated explicitly. The linear operators on the right are completely independent of each other. Fourth order sounds expensive, and it is not, because we are about to solve it in one direction only. Two questions are still open, though. Where did the other two velocity components go, and what boundary conditions does a fourth-order equation want? Both answers arrive as soon as we go to Fourier space.

Both wall-parallel directions are homogeneous and periodic, so both get Fourier series. Write the wall-normal velocity as a sum over modes: each mode carries a wavenumber in x, a wavenumber in z, and an amplitude that still depends on y and on time. Inside a single mode the wall-parallel derivatives are no longer derivatives. Differentiating in x multiplies by i k x, differentiating in z multiplies by i k z, and the Laplacian collapses to a second derivative in y minus k squared, where k squared is the sum of the two squares. That is the whole content of the transform, and it has already done something drastic. Neither of our two governing equations contains x or z at all any more. Put the two equations through that transform and what is left in each one is an ordinary differential equation in y alone, with time as a parameter. Fourth order for the velocity, second order for the vorticity. And here is the payoff the whole method is built on. Each pair of wavenumbers gets its own pair of equations, and those equations know nothing whatsoever about any other pair. The three-dimensional problem has come apart into a lattice of one-dimensional problems. Walk across that lattice. Here is one mode. Here is another, further out. And here is a third. Three completely independent boundary value problems, and the only thing that differs between them is the single number k squared sitting in the operator. Where have u and w gone? Nowhere. They have become algebra. Continuity, transformed, is one linear relation between the three transformed components. And the definition of the wall-normal vorticity, transformed, is a second linear relation between the same three. In both of them the y derivative of v hat is already known, because v is what we have just marched. So set them side by side as a two by two system for u hat and w hat. Multiply down one diagonal, multiply down the other, and subtract. The determinant is k x squared plus k z squared, which is k squared: nonzero for every mode but one, and we will come back to that exception. Solve it and the horizontal velocities fall out. u hat is one combination of the slope of v hat and the vorticity, w hat is another, and both are divided by k squared. No linear solve, no iteration, no elliptic problem: two multiplications and an addition, per mode, per step. And continuity is not being enforced here at all. It has been used, so the reconstructed field satisfies it identically, to machine precision. One thing is still missing, and it is not a small thing. This reconstruction divides by k squared, and there is exactly one mode in the lattice where k squared is zero. That mode, and the boundary conditions, come next.

A fourth-order equation wants four boundary conditions, two at each wall. A second-order equation wants two. So we need six, and the walls appear to offer only one thing: the velocity vanishes. Start with what is given. At either wall all three components are zero, for every x and every z, so mode by mode every one of the three amplitudes is zero as well. That is v hat equal to zero at both walls: two conditions for the fourth-order equation. For the other two, put the wall values into the transformed continuity relation. Since u hat and w hat vanish there, the y derivative of v hat has to vanish there too. So the picture at each wall is this: v hat pinned to zero, and pinned flat. Four conditions, and not one of them was invented for numerical convenience. The vorticity is easier still. Wall-normal vorticity is built from u and w alone, and both of those are zero at the wall, so omega y hat vanishes there. Two conditions for a second-order equation. Let me count it once more, because this is the part that makes the formulation feel inevitable. Each unknown, its order in y, and what the walls hand us for it. The wall-normal velocity obeys a fourth-order equation, and the walls give exactly four conditions: the value, and the slope, at each of the two walls. The wall-normal vorticity obeys a second-order equation, and the walls give exactly two. Nothing is left over, and nothing has to be guessed. Compare that with a boundary condition for the pressure, which does not exist. Now the exception. The reconstruction of the horizontal velocities divided by k squared, and there is exactly one mode where k squared is zero: no variation in x and no variation in z. That mode is the mean flow, averaged over wall-parallel planes, and everything we have built so far says nothing about it. One of its components is free. Continuity for that mode says the y derivative of v hat is zero, and v hat is zero at the wall, so v hat is identically zero right across the gap. There is no mean flow through the walls, which is exactly what we should expect. The two horizontal means are not determined algebraically, so they are marched directly. Average the x momentum equation over a wall-parallel plane and you get a one-dimensional diffusion equation for u bar, forced by the mean nonlinear term and by the mean pressure gradient. And there is the one piece of the pressure that survives the whole construction: a single number, the mean streamwise gradient, which is what drives the flow. Either you fix it and let the flow rate settle, or you fix the flow rate and adjust it every step. The spanwise mean gets the same treatment with no imposed gradient. Both are one-dimensional diffusion equations with zero at each wall, and the machinery of the next part solves them without noticing that anything special has happened.

The wall-normal direction is neither periodic nor homogeneous, so Fourier is the wrong basis there. Expand in Chebyshev polynomials instead, on the interval from minus one to plus one, which is exactly the gap. The first few look like this. T one is just y. T four already has four zeros inside the gap. T eight has eight, and notice where they crowd: the oscillations bunch up towards the two ends. So the expansion is a finite sum, N plus one coefficients for each wavenumber pair, and those coefficients are what the code actually stores. The collocation points come with the basis rather than being chosen. They are cosines of equally spaced angles, and here they are across the gap for N equal to eight. Look at the spacing. In the middle of the channel the points are as coarse as they ever get. Against the wall they are packed tight, with a spacing that shrinks like one over N squared. That is precisely the grid a wall-bounded flow wants, because the whole difficulty of this problem is a thin viscous layer against each wall. Now, what does the operator D squared minus k squared look like in this basis? Written for the coefficients, either as a tau formulation or as a Galerkin one, the second derivative connects a coefficient to only a few of its neighbours. The matrix is banded. For the second-order operator that means three diagonals. For the fourth-order operator, five. Either way it is a fixed number of entries per row, whatever N is, so one solve costs order N operations rather than order N cubed. Two further economies, and they are the ones that matter in practice. The matrix depends on the wavenumbers only through k squared, and on the time step, and neither of those changes as the simulation runs. So it is factored once, before the first step, and every step afterwards is a back substitution. And one more, which is free. The operator does not mix even polynomials with odd ones, so each system splits into two half-size systems, one per parity. All of that, for every wavenumber pair independently.

Here is one time step for the vorticity equation, written the way it is usually implemented. The viscous term is treated implicitly, by Crank-Nicolson: half of it at the old level, half at the new. The nonlinear source is treated explicitly, by a two-step Adams-Bashforth: three halves of the current value, minus a half of the previous one. No solve is involved, because both are already known. And look at what the left hand side is. It is the banded, prefactored operator from a moment ago, one per wavenumber pair, unchanged from step to step. The whole advance is a back substitution per mode, plus the work of forming h. Now the stability, which is the real reason for treating viscosity implicitly. The time step is no longer limited by the wall-normal spacing at all. That spacing goes like one over N squared, so an explicit viscous treatment would need a step going like one over N to the fourth. On a fine grid that is simply unusable. The velocity equation needs one more idea. It is fourth order, so rather than building a fourth-order operator we split it into two second-order problems. First solve for phi, the Laplacian of v hat. Then solve for v hat itself, with phi as the source. There is a catch here, and it is worth naming. All four boundary conditions are conditions on v hat: two values and two slopes. Neither of the two second-order problems has any natural condition on phi at the walls. The standard remedy is a Green's function argument. Solve the pair once with the real forcing and zero conditions on phi, then twice more with unit conditions and no forcing, and take the combination of the three that makes the slope of v hat vanish at both walls. That is a two by two solve per mode on quantities computed before the run. That leaves the source terms h, which are the only nonlinear work in the method, and they are evaluated pseudospectrally. Transform the modes back to a physical grid. Form the products there, where a product is just a multiplication. Then transform the result back to modes. Why not form the products mode by mode instead? Because a product of two Fourier series is a convolution, which costs order N squared per direction, while two transforms and a multiplication cost order N log N. That single choice is what makes spectral direct simulation affordable at all. The price is aliasing. A product of two modes inside the retained band generates content outside it, and on a finite grid that content folds back onto the modes you are keeping, as error. So in both wall-parallel directions the products are formed on an enlarged grid and everything above the retained band is thrown away rather than folded back. That is the three-halves rule, and a phase shift scheme does the same job by a different route. So let me put the efficiency in one place. First: there is no elliptic pressure solve, at any step, because there is no pressure. Second: the reconstructed velocity is divergence free by construction, so nothing has to be projected or corrected afterwards. Third: the Fourier transform in the two homogeneous directions uncouples the three-dimensional problem into an independent, banded, already-factored boundary value problem for every wavenumber pair. Fourth: every boundary condition the method needs is one the walls actually give, so there is nothing invented near the wall. What is left, per step, is dominated by the transforms: order N log N in the wall-parallel directions, times the number of points across the gap. The solves themselves are linear in the number of Chebyshev modes. That is the whole of it. Two curls removed the pressure. Fourier removed two of the three space dimensions from every solve. Chebyshev put the points where the walls need them. And no-slip turned out to supply exactly the six boundary conditions the two equations were waiting for.