26

如何计算湍流能谱

 2 years ago
source link: https://www.physixfan.com/%e5%a6%82%e4%bd%95%e8%ae%a1%e7%ae%97%e6%b9%8d%e6%b5%81%e8%83%bd%e8%b0%b1/
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.

如何计算湍流能谱

作者: 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.


About Joyk


Aggregate valuable and interesting links.
Joyk means Joy of geeK