基于MATLAB的涡流函数方法案例

分享一个基于MATLAB使用伪谱离散化 + 三阶龙格库塔时间积分实现涡流函数方法的CFD案例，并介绍一下相关理论。

1. 代码详解

这段 MATLAB 代码实现了一个二维湍流模拟，使用伪谱法在周期性边界条件下解算非线性涡度-流函数方程：

M = 256; % number of points
N = M;
Lx = 2*pi;
Ly = 2*pi;
nu = 5e-4; % kinematic viscosity
Sc = 0.7; % Schmidt number
beta = 0; % meridional gradient of Coriolis parameter
ar = 0.02; %random number amplitude
b = 1; % mean scalar gradient
CFLmax = 0.8;
tend = 200; % end time

x=linspace(0,Lx,M+1); x(end)=[];
dx = Lx/M;
kx=[0:M/2 -M/2+1:-1]*2*pi/Lx;

y=linspace(0,Ly,N+1); y(end)=[];
ky=[0:N/2 -N/2+1:-1]*2*pi/Ly;
dy = Ly/N;

time = 0;

index_kmax = ceil(M/3);
kmax = kx(index_kmax);
filter = ones(M,N);
filter(index_kmax+1:2*index_kmax+3,index_kmax+1:2*index_kmax+3)=0;

rng(64);
[u, v, omega, psi, ddx, ddy, idel2, kk, k2]=deal(zeros(M,N));

for j=1:N
    ddx(:,j)=1i*kx;
end

for i=1:M
    ddy(i,:)=1i*ky;
end

for i=1:M
    for j=1:N
        idel2(i,j)=-kx(i)^2-ky(j)^2;
    end
end
idel2=1./idel2;
idel2(1,1)=0;
for i=1:M
    for j=1:N
        kk(i,j)=kx(i)^2+ky(j)^2;
        k2(i,j)=kx(i)^2+ky(j)^2;

        if kk(i,j) >= 6^2 && kk(i,j) <= 7^2 % forcing
            kk(i,j) = -kk(i,j);
        end
        if kk(i,j) <= 2^2
            kk(i,j) = 8*kk(i,j); % large-scale dissipation
        end
    end
end
for i=1:M
    for j=1:N
        u(i,j) =  cos(2*x(i))*sin(2*y(j))+ar*rand;
        v(i,j) = -sin(2*x(i))*cos(2*y(j))+ar*rand;
    end
end
uhat = fft2(u);
vhat = fft2(v);
omegahat = ddx.*vhat - ddy.*uhat; % make vorticity
phi = rand(size(u));
phihat = fft2(phi);
dt = 0.5*min([dx dy]);
nstep = 1;
s2 = pcolor(x,y,phi'); 
title('Scalar');  
shading flat; 
axis equal tight; 
drawnow
while time < tend
    psihat = -idel2.*omegahat;
    uhat = ddy.*psihat;
    vhat = -ddx.*psihat;
    u = real(ifft2(uhat));

    v = real(ifft2(vhat));
    omegadx = real(ifft2(ddx.*omegahat));
    omegady = real(ifft2(ddy.*omegahat));
    facto = exp(-nu*8/15*dt*kk);
    factp = exp(-nu/Sc*8/15*dt*k2);
    r0o = -fft2(u.*omegadx+v.*omegady)+beta*vhat;

    r0p = -fft2(u.*real(ifft2(ddx.*phihat))+v.*real(ifft2(ddy.*phihat)))+b*vhat;
    omegahat = facto.*(omegahat + dt*8/15*r0o); % update omega
    phihat = factp.*(phihat + dt*8/15*r0p); % update phi

    %%%% Substep 2
    psihat = -idel2.*omegahat;

    uhat = ddy.*psihat;

    vhat = -ddx.*psihat;
    u = real(ifft2(uhat));

    v = real(ifft2(vhat));
    omegadx = real(ifft2(ddx.*omegahat));
    omegady = real(ifft2(ddy.*omegahat));
    r1o = -fft2(u.*omegadx+v.*omegady)+beta*vhat;
    r1p = -fft2(u.*real(ifft2(ddx.*phihat))+v.*real(ifft2(ddy.*phihat)))+b*vhat;
    omegahat = omegahat + dt*(-17/60*facto.*r0o + 5/12*r1o);
    phihat = phihat + dt*(-17/60*factp.*r0p + 5/12*r1p);
    facto = exp(-nu*(-17/60+5/12)*dt*kk);
    factp = exp(-nu/Sc*(-17/60+5/12)*dt*k2);
    omegahat = omegahat.*facto;
    phihat = phihat.*factp;
    %%%% Substep 3
    psihat = -idel2.*omegahat;
    uhat = ddy.*psihat;
    vhat = -ddx.*psihat;
    % max(max(abs(real(ifft2(1i*ddx.*uhat+1i*ddy.*vhat))))) % divergence
    u = real(ifft2(uhat));
    v = real(ifft2(vhat));
    omegadx = real(ifft2(ddx.*omegahat));
    omegady = real(ifft2(ddy.*omegahat));

    r2o = -fft2(u.*omegadx+v.*omegady)+beta*vhat;
    r2p = -fft2(u.*real(ifft2(ddx.*phihat))+v.*real(ifft2(ddy.*phihat)))+b*vhat;
    omegahat = omegahat + dt*(-5/12*facto.*r1o + 3/4*r2o);
    phihat = phihat + dt*(-5/12*factp.*r1p + 3/4*r2p);
    facto = exp(-nu*(-5/12+3/4)*dt*kk);
    factp = exp(-nu/Sc*(-5/12+3/4)*dt*kk);
    omegahat = omegahat.*facto;
    phihat = phihat.*factp;
    phihat = filter.*phihat;
    omegahat = filter.*omegahat;
    time = time + dt;
    nstep = nstep + 1;
    CFL = max(max(abs(u)))/dx*dt+max(max(abs(v)))/dy*dt;
    if mod(nstep,20)==0
        phi = real(ifft2(phihat));
        omega = real(ifft2(omegahat));
        dissipation = 2*nu*(real(ifft2(ddx.*uhat)).^2 + real(ifft2(ddy.*uhat)).^2 + real(ifft2(ddx.*vhat)).^2 + real(ifft2(ddy.*vhat)).^2);
        eta = (nu^3/mean(dissipation,'all'))^0.25;
        s2.CData = phi';
        drawnow
        fprintf(1,'step = %d    time = %g    dt = %g  CFL = %g    kmax*eta = %g %g\n', nstep, time, dt, CFL, eta.*kmax, eta.*kmax/sqrt(Sc));
    end
    dt = CFLmax/CFL*dt; %0.0005;
end

1.1初始化参数

M = 256; % 网格点数（沿x方向）
N = M; % 网格点数（沿y方向）
Lx = 2*pi; % x方向的域长度
Ly = 2*pi; % y方向的域长度
nu = 5e-4; % 动力粘性系数
Sc = 0.7; % 施密特数（分子扩散率与动量扩散率的比值）
beta = 0; % 科里奥利系数
ar = 0.02; % 随机数幅度
b = 1; % 标量梯度的均值
CFLmax = 0.8; % 最大CFL数
tend = 200; % 结束时间

• M 和 N 分别是 x 和 y 方向的网格点数。
• Lx 和 Ly 是计算域的长度。
• nu 是流体的运动粘度。
• Sc 是施密特数，用于模拟标量的扩散。
• beta 是科里奥利系数的经向梯度。
• ar 是用于初始化速度场的随机扰动的幅度。
• b 是标量梯度的均值，用于定义标量场的强度。

1.2 网格与频率

x = linspace(0,Lx,M+1); x(end)=[];
dx = Lx/M;
kx=[0:M/2 -M/2+1:-1]*2*pi/Lx;

y = linspace(0,Ly,N+1); y(end)=[];
ky=[0:N/2 -N/2+1:-1]*2*pi/Ly;
dy = Ly/N;

• x 和 y 是物理空间中的网格点。
• kx 和 ky 是对应的频率空间中的波数向量。它们用于傅里叶变换后的计算。

1.3 时间循环的初始化

time = 0;
index_kmax = ceil(M/3);
kmax = kx(index_kmax);
filter = ones(M,N);
filter(index_kmax+1:2*index_kmax+3,index_kmax+1:2*index_kmax+3)=0;

rng(64);
[u, v, omega, psi, ddx, ddy, idel2, kk, k2] = deal(zeros(M,N));

• time 初始化为 0。
• index_kmax 和 kmax 用于确定最大有效波数。
• filter 是一个频率滤波器，用于在频率空间中滤除高频噪声。
• u, v, omega, psi 等变量初始化为零矩阵。

1.4 初始化导数和拉普拉斯算子

for j=1:N
    ddx(:,j) = 1i*kx;
end

for i=1:M
    ddy(i,:) = 1i*ky;
end

for i=1:M
    for j=1:N
        idel2(i,j) = -kx(i)^2 - ky(j)^2;
    end
end
idel2 = 1./idel2;
idel2(1,1) = 0;

• ddx 和 ddy 分别是 x 和 y 方向的微分算子，使用频率空间中的傅里叶变换定义。
• idel2 是拉普拉斯算子的逆，用于求解流函数， idel2(1,1)=0 是为了避免数值发散。

1.5 初始条件

for i=1:M
    for j=1:N
        kk(i,j) = kx(i)^2 + ky(j)^2;
        k2(i,j) = kx(i)^2 + ky(j)^2;

        if kk(i,j) >= 6^2 && kk(i,j) <= 7^2 % 强制项
            kk(i,j) = -kk(i,j);
        end
        if kk(i,j) <= 2^2
            kk(i,j) = 8*kk(i,j); % 大尺度耗散
        end
    end
end

• kk 和 k2 是频率空间中的波数平方，用于计算非线性项。
• 对应不同的波数范围，对 kk 进行调整，用于引入大尺度耗散和强制项。

1.6 初始化速度场与标量场

for i=1:M
    for j=1:N
        u(i,j) = cos(2*x(i)) * sin(2*y(j)) + ar*rand;
        v(i,j) = -sin(2*x(i)) * cos(2*y(j)) + ar*rand;
    end
end

• u 和 v 是速度场的初始条件，由一个基态流（cos-sin 形式）和小幅度随机扰动组成。

1.7 时间推进循环

主循环依次计算涡度、流函数、速度场和标量场的时间推进。每个时间步分为三部分（子步骤），并且使用了一种三阶 Runge-Kutta 时间积分方法。每一步中，代码会更新各个场量并应用滤波器以防止数值发散。

最终，在循环结束后，代码会显示标量场的变化并打印出 CFL 数等关键信息。

2. 涉及的理论

该案例涉及多个计算流体力学的理论和数值方法，包括湍流建模、伪谱法、傅里叶变换、CFL条件、三阶Runge-Kutta方法，以及拉普拉斯方程的求解。

2.1 二维湍流

二维湍流在磁流体力学和地球物理流体力学（尤其是大气动力学）方向讨论的比较多

二维湍流在许多方面与三维湍流不同，其中最显著的区别是能量级联方向的差异。

在三维湍流中，能量从大尺度传递到小尺度（称为正能量级联），最终由于粘性作用而耗散。

而在二维湍流中，能量级联有两个方向：

• 正能量级联：从大尺度向小尺度传递，与三维湍流相似。
• 逆能量级联：从小尺度向大尺度传递，这在三维湍流中并不存在。

由于逆能量级联的存在，二维湍流通常会形成大尺度的涡旋结构，就如同本文案例中所展现的效果。

2.2 伪谱法（Pseudo-spectral Method）

伪谱法（拟谱法 ） 是一种数值方法，用于求解偏微分方程，特别是对流和扩散问题。在伪谱法中，非线性项在物理空间中计算，然后通过傅里叶变换转换到频率空间，以便利用快速傅里叶变换（FFT）高效地计算空间导数。

伪谱法的基本步骤如下：

1. 初始条件：定义物理空间中的场变量（如速度、标量）。
2. 傅里叶变换：将场变量从物理空间转换到频域空间。
3. 计算导数：在频率空间中，使用简单的乘法计算空间导数。
4. 傅里叶逆变换：将结果转换回物理空间，以更新场变量。
5. 时间推进：根据时间积分方案更新场变量。

伪谱法的优点是可以精确地计算导数（由于傅里叶变换的精度），但缺点是容易受到高频噪声的影响，因此通常需要使用滤波器来控制误差。

2.3 傅里叶变换（Fourier Transform）

傅里叶变换是一种将函数从时域或空间域转换到频域或波数域的方法。对于离散数据，使用的是离散傅里叶变换（DFT），其快速实现版本为快速傅里叶变换（FFT）。在湍流模拟中，傅里叶变换的作用包括：

• 频率空间的导数计算：在频率空间中，空间导数可以通过乘以虚数单位 1i 和波数 k 来实现。
• 解析偏微分方程：许多方程（如拉普拉斯方程）在频率空间中更容易解析。

在代码中，傅里叶变换用于计算速度场和涡度场的导数，以及解拉普拉斯方程来得到流函数。

2.4 CFL条件（CFL Condition）

CFL条件（Courant–Friedrichs–Lewy 条件）是数值模拟中的稳定性条件。它决定了时间步长 dt 必须足够小，以确保解的稳定性。

具体而言，CFL条件要求：

$$\text{CFL}=\frac{\text{最大速度}\times \Delta t}{\Delta x}<\text{CFLmax}$$

代码中每20步重新计算 CFL 数，用于调整时间步长 dt，以确保数值模拟的稳定性。

2.5 三阶Runge-Kutta方法

Runge-Kutta方法是一类用于求解常微分方程的数值积分方法。三阶Runge-Kutta方法是一种高级的时间积分方案，具有较高的精度和稳定性。

本文介绍的案例中使用的方案有以下三个子步骤：

1. 计算初始梯度，并更新场变量。
2. 基于第一次更新的结果，计算第二次更新。
3. 使用前两次更新的结果，计算最终的场变量。

这种方法可以有效地降低数值误差，并且适用于复杂的非线性系统。

2.6 拉普拉斯方程的求解

拉普拉斯方程是流体动力学中的常见方程，用于描述势场。

对于二维涡量方程，可以通过解拉普拉斯方程得到流函数 $(\psi )$：

$${\nabla }^2\psi =\omega$$

在频率空间中，拉普拉斯方程可以简化为代数方程：

$$\hat{\psi }=-\frac{\hat{\omega }}{k^2}$$

其中 ：k^2 是波数（Wavenumber）的平方，$\hat{\varphi }$ 表示流函数在傅里叶空间中的表示。

在代码中，通过 idel2（即 -1/k^2）计算 $\hat{\varphi }$。

2.7 能量和耗散

湍流模拟中，能量的传递和耗散是关键要素。通过波数 k 的调整，代码中引入了强制项和大尺度耗散项，以模拟真实湍流中能量的输入和耗散：

• 强制项：在中等尺度上，代码通过调整 kk 来引入能量（对应于物理上能量的输入）。
• 大尺度耗散：在大尺度上，通过增加 kk 的值来增强耗散，从而防止大尺度能量的积累。

这种方法在湍流研究中非常普遍，用于保持数值稳定性并模仿实际物理现象。