Skip to content

March 4, 2015

7

如何计算湍流能谱

作者: physixfan

最近在如何数值上计算湍流能谱(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的格点,这个时间复杂度是无论如何也算不出来了。。。

正确的算法是这样的:

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

\Phi_{ij}(\mathbf{k})=F\{R_{ij}(\mathbf{r})\}=F^*\{u_i(\mathbf{r})\}F\{u_j(\mathbf{r})\}=u^*_i(\mathbf{k})u_j(\mathbf{k})

其中^*表示复共轭。

因此正确的算法是,先计算速度场的傅里叶变换\mathbf{u}(\mathbf{k}),然后用上式来计算\Phi_{ij}(\mathbf{k}),然后就可以得到E_k(k)了。总时间复杂度的限制在傅里叶变换上,最后的总时间复杂度将是O(n^2\log n)

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

参考资料:

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

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

Read more from Magical Physics