Skip to content

March 4, 2015

9

如何计算湍流能谱

作者: 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
9 Comments Post a comment
  1. Apr 28 2015

    。。。

    Reply
  2. May 30 2015

    横坐标标注什么?w纵坐标标注啥?E?湍流能谱代码%% understand the turbulent energy spectrum clear all U_max =10;N=100 [x,y] = meshgrid(linspace(0,1,N)); Vx = U_max*cos(2*pi*x+pi/2).*cos(2*pi*y); Vy = U_max*sin(2*pi*x+pi/2).*sin(2*pi*y); u=Vx;v=Vy; % 计算能量 Fu=fftshift(fft2(u)); Fv=fftshift(fft2(v)); phy11=conj(Fu).*Fu;phy22=conj(Fv).*Fv; phy=phy11+phy22; mesh(log(phy+1)); cx=(N+1)./2;cy=(N+1)./2; for i=1:N for j=1:N A(i,j)=round(sqrt((i-cx).^2+(j-cy).^2)); end end A(A==0)=1; for i=min(A(:)):max(A(:)) esp(i) = sum(phy(A==i));%统计在同一个波长尺度的能量 end plot(log(esp+1));

    Reply
    • 松鼠妹
      Oct 19 2016

      你好,非常感谢分享code!很有帮助。可否请教一下,如果要计算某个方向的Euu(kx), 这里的Euu是energy spectra,相当于你的Phy, 是不是不用A(i,j)=round(sqrt((i-cx).^2+(j-cy).^2)) 而只用A1(i)=round((i-cx),而esp11(i) = sum(phy11(A1==i))?尝试了下好像结果有点奇怪。另外,不知道对于2D u(x,y)的1D energy spectra (比如Euu(kx) 随kx的图)的物理意义又如何理解?切盼指教,非常感谢!

      Reply
  3. Sep 15 2015

    学长在读研吗,现在?

    Reply
  4. 冰冰
    Dec 24 2015

    通信里有这个类似的东西,自相关函数的傅里叶变换求出能谱密度,两个信号的互相关函数求出互能谱密度。

    Reply
  5. Jan 4 2016

    网站不错,过来看看!

    Reply
  6. Muyan
    Dec 21 2016

    太感谢了!查了好多教材都没有谈自相关函数的具体计算上,我也纳闷要是真按着原始定义计算这复杂度简直爆炸。

    Reply
  7. 白晓
    Oct 21 2018

    您好,打扰了。最近想探讨水翼流场的湍流级串方面的问题,但在计算流场的湍流能谱问题上一直一知半解,想请问您一下,一般课本里和您这里提出的方法是不是只适用于均匀各向同性湍流呢,对于非均匀湍流不同位置也就是X的影响应如何考虑呢?

    Reply
  8. Shanshan
    Feb 21 2019

    公式什么现在为什么都看不见啦?

    Reply

Share your thoughts, post a comment.

(required)
(required)

Note: HTML is allowed. Your email address will never be published.

Subscribe to comments