如何计算湍流能谱

转自： http://www.physixfan.com/archives/2605

作者： physixfan

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

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

能谱的含义是动能在k空间上的分布函数，其对k进行积分之后将得到总的动能。其定义为：

Ek(k)=∫12Φii(k)δ(|k|−k)dk

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

其中 Φij 为关联函数（correlation function）的傅里叶变换：

Φij(k)=F{Rij(r)}

其中 F{.} 是傅里叶变换的记号，而 Rij 为关联函数：

Rij(r)=⟨ui(x)uj(x+r)⟩

其中 ⟨.⟩ 表示对位置 x 的平均。

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

正确的算法是这样的：

根据 Cross-Correlation Theorem，或者其特列 Wiener-Khinchin Theorem，cross-correlation 等价于傅里叶空间中的乘积！大一时候微积分应该是学过这个定理的，没想到这么多年后居然用到了...因此我们可以直接简便地计算关联函数的傅里叶变换：

Φij(k)=F{Rij(r)}=F∗{ui(r)}F{uj(r)}=u∗i(k)uj(k)

其中 ∗ 表示复共轭。

因此正确的算法是，先计算速度场的傅里叶变换 u(k) ，然后用上式来计算 Φij(k) ，然后就可以得到 Ek(k) 了。总时间复杂度的限制在傅里叶变换上，最后的总时间复杂度将是 O(n2logn) 。

注意：先计算动能场 E(x)=ui(x)ui(x) 然后直接对它进行傅里叶变换是得不到能谱的。

参考资料：

Pope, Turbulent Flows, 1st Edition, Page 219~221.

Press, Numerical Recipes in C, 2nd Edition, Page 545.