Skip to content

March 4, 2015

如何计算湍流能谱

作者: 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
  1. 横坐标标注什么?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));

    • 你好,非常感谢分享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的图)的物理意义又如何理解?切盼指教,非常感谢!

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

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

无觅相关文章插件,快速提升流量