四层空间的交叉积，假想的人工制品出现

Question

我正在编写一个代码来模拟周期立方体中的湍流。 为此，我从形状为[64,64,64,3]的速度向量场开始，它意味着每个方向（x，y，z）上有64个网格点，在这个空间的每个点上有3个速度分量（ux，uy，uz）。 为了计算下一个时间步，必须从速度场计算涡度场，出现了这个问题。你知道吗

例如，让我们考虑一下表单的所谓ABC流

在这种情况下，速度场的旋转保持不变，即速度场u和涡度场omega相等。要计算出涡量场以外的速度场，可以使用如下公式

为了计算导数，我使用了伪谱方法，使得方程的读数是这样的

。你知道吗

其中uk和ωk是Fourier空间中的速度场和涡量场。你知道吗

我试着编写一个程序，从Fourier空间的ABC流场计算出涡度，然后从涡度场计算出速度场 . 在真实空间中，速度场和涡度场看起来应该是一样的，但正如你在这里看到的，它们不是

图中显示了立方体外单层的第一个矢量分量的实部，在实空间和傅里叶空间中，开始时的速度（pre）、涡度和由涡度计算出的速度（post）。 正如你所看到的，速度场看起来和预期的相似，而实际空间中的涡量场显示出条纹状的图案，尽管看起来应该是一样的。你知道吗

看一下实空间中涡量场的虚部，会发现一个条纹状的图案，它与实部的图案相吻合。你知道吗

我现在的问题是找到建立这种模式的原因。 涡度应该看起来像速度场，但它不是，我不知道为什么。 解决这个问题对于避免在计算下一个时间步时丢失信息非常重要。你知道吗

我的代码如下所示：

import numpy as np
import matplotlib.pyplot as plt

xx = np.arange(64)                 #creating 1D x coordinate
x,y,z = np.meshgrid(xx,xx,xx)      #calculating 3D x,y,z coordinates
a = np.zeros(shape=(64,64,64,3))   #creating 3D vector field 'a'
A, B, C = 3., 3., 3.               #setting values A,B,C
a[:,:,:,0] = A*np.cos(y/(2.*np.pi)) + B*np.sin(z/(2.*np.pi))  
a[:,:,:,1] = B*np.cos(z/(2.*np.pi)) + C*np.sin(x/(2.*np.pi))
a[:,:,:,2] = C*np.cos(x/(2.*np.pi)) + A*np.sin(y/(2.*np.pi))


k = np.zeros(shape=a.shape)                   # creating k array
k_sqr = np.zeros(shape=a.shape)               # creating k^2 array
k_lin = np.fft.fftfreq(64)*2.*np.pi           # calculating 1D k coordinate
k[:,:,:,0], k[:,:,:,1], k[:,:,:,2] = 
np.meshgrid(k_lin,k_lin,k_lin,indexing='xy')  # calculating 3D kx,ky,kz 
k_sqr[:,:,:,0] =  k[:,:,:,0]**2 + k[:,:,:,1]**2 + k[:,:,:,2]**2
k_sqr[:,:,:,1] =  k_sqr[:,:,:,0]
k_sqr[:,:,:,2] =  k_sqr[:,:,:,0]
k_sqr[0,0,0,:] += 1.      # set zero of k_sqr = 1. to avoide division by 0 

a_fft = np.fft.fftn(a,axes=(0,1,2))              # calculate FFT of a (=uk)

wk = np.cross(1j*k,a_fft)        # calculating vorticity wk by wk = ik x uk

w = np.fft.ifftn(wk,axes = (0,1,2))         # calculating iFFT von wk (w(x))

uk = np.cross(1j*k,wk)/k_sqr        # calculating uk by uk = (ik x wk)/k^2
uk[0,0,0,:] = a_fft[0,0,0,:]                                                        
# set value of position of 0-Division to the value it was before

u = np.fft.ifftn(uk,axes = (0,1,2))        # calculating iFFT of uk (u(x))

我还没有找到关于这个问题的任何信息。 如果这个问题已经存在，我很抱歉。 如果没有我很感激任何想法。你知道吗