如何计算湍流能谱

最近在如何数值上计算湍流能谱（turbulent energy spectrum）上面卡了好几天，现在终于问到了做过的人看了相应的书明白了正确的算法。写出来以供自己留个记录以及给搜索引擎喂食吧。

问题定义：现在已知速度场\(\mathbf{u}(\mathbf{x})\)，求能谱\(E_k(k)\)。（我研究的问题是2D的，所以以下都按照2D来写了，推广到3D也是非常容易的事情。）

能谱的含义是动能在k空间上的分布函数，其对k进行积分之后将得到总的动能。其定义为：\[E_k(k)=\int\frac{1}{2}\Phi_{ii}(\mathbf{k})\delta(|\mathbf{k}|-k)d\mathbf{k}\]

其中\(\Phi_{ii}=Tr(\Phi_{ij})\)（其中ij取12，如果是3D情况则ij取123）

其中\(\Phi_{ij}\)为关联函数（correlation function）的傅里叶变换：\[\Phi_{ij}(\mathbf{k})=F\{R_{ij}(\mathbf{r})\}\]

其中\(F\{.\}\)是傅里叶变换的记号，而\(R_{ij}\)为关联函数：\[R_{ij}(\mathbf{r})=\langle u_i(\mathbf{x})u_j(\mathbf{x}+\mathbf{r})\rangle\]

其中\(\langle . \rangle\)表示对位置\(\mathbf{x}\)的平均。

至此，根据定义，理论上如何用速度场\(\mathbf{u}(\mathbf{x})\)求能谱\(E_k(k)\)已经非常清楚了。但是问题是，根据这个定义进行数值计算的话，2D情况下单单是计算\(R_{ij}(\mathbf{r})\)的时间复杂度就已经是\(O(n^4)\)了！这是因为\(R_{ij}(\mathbf{r})\)本来就是一个二维的场量，而计算其中每个分量又需要对整个\(\mathbf{x}\)做平均。如果是1024*1024的格点，这个时间复杂度是无论如何也算不出来了。。。

正确的算法是这样的：…