# Why can't you simulate isotropic fluid flow on a square lattice?

There are easy methods for discrete simulations of gas dispersion in two dimensions. If you take a large square lattice, each cell of which is assumed to contain at most one gas molecule, and you move the molecules from cell to adjacent cell at random, the large-scale results are in many ways a good simulation of gas behavior. In particular, even though molecules are individually moving in only the four cardinal directions, at a large scale the simulation does not show any bias to these four directions. If you start with a large square lump of material, it will rapidly diffuse into a circle. This property of the simulation is called isotropy.

Long ago when I was involved in research in this area I was told that this worked fine for gas flow simulations, but not for incompressible fluids. If one wanted to make an isotropic two-dimensional simulation of incompressible fluid flow, I was told, one had to use a hexagonal lattice instead of a square lattice. I was told that one could deduce from the Navier-Stokes equations that that any simulation of incompressible fluid on a square lattice would necessarily be anisotropic.

Is this correct? If so, how does the proof go? If the argument is complex, is it possible to get an intuitive idea of why compressible and incompressible fluids are different in this regard? What is a reference that would include the full proof?

Addendum: Despite the answer I posted below, I still don't understand any of the details. I would be glad to award the bounty to someone who could explain it.

-
you mean that, in contrast to the square grid, simulation on a hexagonal grid would be perfectly isotropic? or only "less anisotropic"? –  begeistzwerst May 14 '12 at 15:24
I think the former, but I do not know if isotropy is an all-or-nothing property. –  MJD May 14 '12 at 15:27
i think it is, but still it might make sense to distinguish solutions by the extent to which they differ from an ideal isotropic one. anyway, my point is that i find it rather implausible that a switch to a hexagonal grid should render the simulation perfectly isotropic. but that's only a guess and i am far from being an expert. –  begeistzwerst May 14 '12 at 16:04
That's what I found remarkable also, but that's what I heard. I'd be just as happy to find out that it is not true. –  MJD May 14 '12 at 16:09
@begeistzwerst It seems that I understood correctly! Check out this answer. –  MJD May 15 '12 at 17:55

It seems that this might be the original source: Lattice gas automata for the Navier-Stokes equations, by Frisch et al. (PDF)

The [square lattice] HPP automaton is invariant under $\pi/2$ rotations. Such a lattice symmetry is insufficient to ensure the isotropy of the fourth degree tensor relating momentum flux to quadratic terms in the velocity. …
Observe that when the underlying microworld is two-dimensional and invariant under the hexagonal rotation group (multiples of $\pi/3$) the tensor $T$ is isotropic…

Unfortunately it will be a long time before I am able to understand this properly; I would still prefer to see an answer that explains what is going on.

-
Ah, this is not a limitation of general grid-based methods, but rather of the specific method of lattice gas automata which attempts to simulate fluids through cellular automata. As far as I am aware, the more commonly used grid-based fluid simulation methods that directly discretize the Navier-Stokes equations do not have this problem. –  Rahul May 15 '12 at 18:52
Thanks. I would be grateful for any elaboration or other explanation you could give. –  MJD May 15 '12 at 19:35
I'm afraid isotropy of fourth degree tensors is beyond my knowledge. I can make a conjecture: a second-order tensor being invariant under rotations of $\pi$ does not imply that it is necessarily a multiple of the identity, but invariance under rotations of $\pi/2$ does. Perhaps an analogous fact with $\pi/2$ and $\pi/3$ holds for fourth-order tensors. But you've probably thought of this already. –  Rahul May 18 '12 at 16:40
No, I know next to nothing about tensors of any order. –  MJD May 18 '12 at 17:51

As you say you are not very familiar with tensors (and I wouldn't call myself an expert either), I will try to explain everything in linear-algebraic language. To quote the paper,

...for the HPP model, the momentum flux tensor has the form $$P_{\alpha\beta} = p\delta_{ab} + T_{\alpha\beta\gamma\epsilon}u_\gamma u_\epsilon + O(u^4).$$ ...The tensor $T$ is, by construction, symmetric in both $(\alpha,\beta)$ and $(\gamma,\epsilon)$.

Here $T$ is "the fourth degree tensor relating momentum flux to quadratic terms in the velocity". We can think of the second-order tensors $P_{\alpha\beta}$ and $u_\gamma u_\epsilon$ as symmetric $2\times2$ matrices ($P_{\alpha\beta}$ must be symmetric to conserve angular momentum, while $u_\gamma u_\epsilon$ is $uu^T$ which is clearly symmetric), and $T$ as a linear transformation mapping the latter to the former.

A tensor is isotropic if it looks the same under any rotation of the coordinate frame. In this case, that means that if $T$ maps some matrix $X$ to a matrix $Y$, and you rotate the coordinate frame so that $X$ becomes $X'$ and $Y$ becomes $Y'$, then $T$ still maps $X'$ to $Y'$. Formally, if a rotation transforms vectors as $v' = Rv$, then matrices transform as $A' = RAR^T$, and we want $T(RXR^T)=RT(X)R^T$.

Things become clearer if we make a judicious choice of basis for the three-dimensional space of symmetric $2\times2$ matrices, namely $$E_1 = \begin{bmatrix}1 & 0\\0 & 1\end{bmatrix},\quad E_2 = \begin{bmatrix}1 & 0\\0 & -1\end{bmatrix},\quad E_3 = \begin{bmatrix}0 & 1\\1 & 0\end{bmatrix}.$$ Then it turns out that a rotation by an angle $\theta$ transforms this basis as \begin{align} E_1' &= E_1, \\ E_2' &= E_2\cos2\theta + E_3\sin2\theta, \\ E_3' &= E_3\cos2\theta - E_2\sin2\theta, \end{align} which one can also see by visualizing them as quadratic functions of the form $x \mapsto x^TEx$ over $\mathbb R^2$. In other words, a rotation transforms $E_2$ and $E_3$ into each other, so for $T$ to be invariant under such a transformation, it has to do "the same thing" to both $E_2$ and $E_3$, in some sense. But if you only look at rotations of $\pi/2$, you don't see this dependence, because then $E_2' = -E_2$ and $E_3' = -E_3$. So a tensor that treated $E_2$ and $E_3$ differently could be invariant under rotations of $\pi/2$, but it would not be invariant under arbitrary rotations. For example, we could choose $T$ such that $$T(E_1) = T(E_2) = 0,\quad T(E_3) = E_3.$$ You can verify that this tensor is invariant under a rotation of $\pi/2$, but not under a rotation of $\pi/4$ (where $E_2' = E_3$ and $E_3' = -E_2$), so it is not isotropic.

Anyway, since the title of the question still asks "Why can't you simulate isotropic fluid flow on a square lattice?", allow me to reiterate my comment that this is an issue with lattice gas automata in particular, not with fluid simulations on square grids in general. The problem with isotropy comes in when trying to show the convergence of the cellular automaton model to the Navier-Stokes equations. As far as I am aware, you can simulate isotropic fluid flow on a square grid through the usual methods of directly discretizing the Navier-Stokes equations and solving them numerically, without going through an intermediate cellular automata model.

-
Thanks very much. I will reread this carefully and digest it. What is a "quadratic height field"? I was not able to find this with a cursory web search. –  MJD Sep 2 '12 at 1:44
Oh, sorry, that's just my graphics background showing. I just meant that you think of a quadratic function $f(x,y) = a_{11}x^2 + a_{22}y^2 + 2a_{12}xy$ as the height of a surface $z=f(x,y)$ in $\mathbb R^3$. –  Rahul Sep 2 '12 at 1:49
I realized that my choice of basis was poor. I've updated the answer and it should be clearer now. –  Rahul Sep 2 '12 at 8:05