简介这是一份基于MATLAB的二维Navier-Stokes方程Fourier-Galerkin谱方法求解框架面向流体力学数值模拟方向的研究生、科研人员及对谱方法有基础的开发者用于求解不可压缩流体在二维空间中的速度场与压力场演化。压缩包共11个文件以10个M函数文件和1个txt说明文件为主整体仅6KBM脚本覆盖预处理、主循环、右端项计算与时间积分txt文件提供代码结构说明目录布局清晰便于按模块阅读。已有350人学习下载。通过研读代码可完整掌握谱方法空间离散、四阶Runge-Kutta时间推进以及周期边界条件的工程实现还可基于面向对象的类封装快速扩展自定义涡旋初始条件复现Taylor-Green涡、混合层等经典算例。整体代码精简而完整是理解谱方法实际编程的优质范例适合二次开发、论文复现与教学参考。1. 为什么是 Fourier-Galerkin 谱方法把 NS 方程推进到“每步只有三次 FFT”二维不可压 Navier-Stokes 方程的数值求解大多数工程师第一反应是有限差分或有限体积画网格、离散对流项、处理压力耦合、盯 CFL 条件。这套流程没有错但它把精度“锁”在了网格分辨率上——想要波数空间里多分辨一倍的结构网格点要翻四倍时间步也要跟着收紧。对于二维流动尤其是拟涡能enstrophy串级显著的涡街、混合层这类问题这种代价并不值得。Fourier-Galerkin 谱方法的思路完全不同解本身被展开成傅里叶基函数微分算子变成波数空间里的代数乘法Galerkin 投影保证展开系数满足的是“能量最优”的截断方程。只要流场足够光滑误差随截断波数增长是指数衰减的而不是有限差分里的代数收敛。配合 Matlab 内置的 FFT每一步的时间推进从“处理系数矩阵”压缩成“正变换 → 谱空间乘法 → 逆变换”这样三次 fft/ifft。本文用二维 NS 方程的流函数-涡量形式给出完整可运行的 Matlab 实现并解释去混淆、时间积分格式和分辨率设置这几个最容易出错的地方。适合已经写过 CFD 代码、想换谱方法提高精度的人也适合刚入门的同学对照公式看实现。2. 流函数-涡量形式的 Fourier 谱离散把压力项消掉再做 Galerkin2.1 为什么标量方程比原始变量形式更适合谱方法二维不可压 NS 方程写成原始变量速度-压力后连续方程是一个约束条件Galerkin 方法处理起来要引入压力 Poisson 方程或者投影算子。常见做法是改写为流函数-涡量ψ-ω形式定义涡量 ω ∂v/∂x − ∂u/∂y 和流函数 u ∂ψ/∂y、v −∂ψ/∂x连续方程自动满足。对动量方程取旋度后压力项消失剩下的是一个标量输运方程。这个处理对谱方法意义重大。Fourier-Galerkin 的核心操作是“把偏微分方程投影到截断的傅里叶模态上”压力项在这个形式里完全不存在规避了压力-速度耦合在谱框架下椭圆方程的复杂度。最终的方程组只有两个场∂ω/∂t J(ψ,ω) ν∇²ω 涡量输运方程 ∇²ψ −ω 运动学关系其中 J(ψ,ω) ∂ψ/∂x · ∂ω/∂y − ∂ψ/∂y · ∂ω/∂x 是雅可比行列式代表对流项的标量形式。第二个方程在每个时刻是一个 Poisson 方程在谱空间里直接除以 −k² 就解完了。2.2 Galerkin 投影下的离散方程每个波数都独立演化设计算域为 [0, Lx] × [0, Ly] 的周期域截断波数分别为 Nx/2 和 Ny/2。涡量场的谱展开写成ω(x,y,t) Σ ω̂(kx,ky,t) · exp(i kx x i ky y)这里的求和遍历所有满足 −Nx/2 ≤ kx Nx/2、−Ny/2 ≤ ky Ny/2 的离散波数。将展开代入涡量输运方程再乘以 exp(−i kx x − i ky y) 并在全域积分利用傅里叶基的正交性得到 Galerkin 截断后的常微分方程组dω̂/dt −Ĵ(ψ,ω) − ν (kx² ky²) ω̂其中 Ĵ(ψ,ω) 代表雅可比项在谱空间的投影系数。对流项的谱投影对应的正是物理空间乘积的 FFT——这成为下一步 Matlab 实现的关键接口。这里需要特别指出一个谱方法特有的问题雅可比项是 ψ 和 ω 的导数在物理空间相乘两个截断级数的乘积会产生高于截断波数的高频模态这些模态在没有被 Galerkin 截断的空间里发生“混叠”aliasing并折返回低频导致能量虚假堆积。避免混叠的常用处理是 2/3 规则在每次计算雅可比前把超过 2N/3 波数的谱系数直接清零。这个操作在后文实现中会体现。2.3 椭圆方程在谱空间的实现除以波数平方就够了每个时间步都要由 ω 恢复 ψ这一步叫“求逆 Laplace”。在谱空间这个操作简单到只是代数除法。但需要注意零波数模态的处理当 kx ky 0 时方程 ∇²ψ −ω 在谱空间的对应是 0 · ψ̂ −ω̂除数为零。这对应流函数的常数自由度与速度场无关直接置 ψ̂(0,0) 0 即可。这个细节单独写一段的原因在于很多第一次写谱方法代码的人会在这一步得到 NaN问题几乎都出在零频除法上。3. 用 Matlab 搭建谱方法求解器最小可运行实现3.1 计算域的网格与波数设置进入代码前先明确坐标系。设 Lx Ly 2π网格数 Nx Ny 128。物理空间用网格点均匀分布谱空间用整数波数。不同之处在于波数数组的排列顺序——Matlab 的 fft 输出从 0 频开始波数数组要按“负频在后半段”的规则构造% fft_output_grid.m clear; clc; % 计算域与分辨率 Lx 2*pi; Ly 2*pi; Nx 128; Ny 128; % 物理空间网格列向量和行向量 x (0:Nx-1) * Lx / Nx; % 列向量 y (0:Ny-1) * Ly / Ny; % 行向量 % 谱空间波数0,1,...,N/2-1,-N/2,...,-1 的排列 kx 2*pi/Lx * [0:Nx/2-1, -Nx/2:-1]; % 列向量 ky 2*pi/Ly * [0:Ny/2-1, -Ny/2:-1]; % 行向量 % 二维波数平方矩阵用于泊松求解和耗散项 [kx2, ky2] meshgrid(ky, kx); K2 kx2.^2 ky2.^2; K2(1,1) 1; % 零频占位后续除以 K2 后手动恢复 0这段代码几次用到了 Matlab 的隐式扩展机制x是列向量y是行向量x y自动生成二维场波数数组同样用列乘行的方式生成 K2 矩阵。K2(1,1) 1是零频保护在除以 K2 之后需要单独把结果的第一行第一列赋值为 0。这里列出常见误区波数数组的正确排列必须与 fft / ifft 的频点对齐错一个位置就会把高频模态当作低频处理。验证方法是用fft作用于任意光滑场后计算若干波数值的平方和与帕塞瓦尔定理对比。3.2 雅可比算子的 FFT 实现三步变换与 2/3 去混叠雅可比项 J(ψ,ω) 的谱系数按定义是物理空间两个乘积项的 FFT。具体做法是先在谱空间算出 ψx、ψy、ωx、ωy 的谱系数乘以 ikx 或 iky通过 ifft 回物理空间对应点相乘后做 fft。这是谱方法最典型的“伪谱”pseudospectral策略导数在谱空间算乘积在物理空间算一步 fft 在不同空间之间切换。function [Jhat] jacobian_psi_omega(psi_hat, omega_hat, kx, ky) % 输入: 谱系数 psi_hat, omega_hat已按 2/3 规则截断 % 输出: 雅可比 J(psi,omega) dpsi/dx * domega/dy - dpsi/dy * domega/dx 的谱系数 % 谱空间计算导数 dpsi_dx_hat 1i * kx .* psi_hat; dpsi_dy_hat 1i * ky .* psi_hat; domega_dx_hat 1i * kx .* omega_hat; domega_dy_hat 1i * ky .* omega_hat; % 变回物理空间做乘积 dpsi_dx real(ifft2(dpsi_dx_hat, symmetric)); dpsi_dy real(ifft2(dpsi_dy_hat, symmetric)); domega_dx real(ifft2(domega_dx_hat, symmetric)); domega_dy real(ifft2(domega_dy_hat, symmetric)); % 物理空间乘积 正变换得到谱系数 J_phys dpsi_dx .* domega_dy - dpsi_dy .* domega_dx; Jhat fft2(J_phys); endsymmetric参数通知 ifft2 输入共轭对称强制返回纯实数组可以避免浮点噪声在虚部累积。注意乘积 J_phys 包含的混叠频率还没有被剔除。下一步是在主程序的时间推进里每次调用完雅可比后立即执行去混淆操作。单独把去混淆写成独立函数是为了便于在不同分辨率下复用function [hat_filtered] dealias_2by3(hat, Nx, Ny) % 2/3 规则去混淆将谱系数中高于 2N/3 波数的分量置零 hat_filtered hat; kmax_x round(Nx/3); % 保留区间 [-N/3, N/3]相当于 2/3 规则 kmax_y round(Ny/3); % 用矩阵索引把四角高频区清零 hat_filtered(kmax_x1:end-kmax_x1, :) 0; hat_filtered(:, kmax_y1:end-kmax_y1) 0; end3.3 时间积分粘性项精确积分 对流项 RK4时间离散格式的选择直接影响稳定性。流场速度 U 下对流项的时间步限制约为 Δt CFL·Δx/U粘性项的限制约为 Δt Δx²/(4ν)。谱方法没有数值耗散如果粘性项走显式格式高波数模态会非常苛刻地压缩时间步。常见做法是采用积分因子integrating factor法把耗散项在时间推进中做精确积分剩余的对流项用显式格式。定义新的因变量ω̃ exp(ν k² t) · ω̂代入离散方程后粘性项被完全吸收方程退化为dω̃/dt −exp(ν k² t) · Ĵ这样处理之后每个子步先算雅可比乘上积分因子再用 RK4 更新。这个格式在 CFL 限制下可以工作得很好。下面给出 RK4 主循环的主体结构% ns2d_fourier.m 主时间循环 nu 1e-4; % 动力学粘性系数 dt 0.002; % 时间步长 T_end 5; % 总模拟时长 t 0; step 0; % 时间积分因子的谱空间矩阵 EF exp(nu * K2 * dt); % RK4 各子步使用 EF.^0.5 或 EF 的中间值 while t T_end % RK4 第一子步 t_sub t; psi_hat -omega_hat ./ K2; psi_hat(1,1) 0; % 由涡量恢复流函数谱 J1 jacobian_psi_omega(psi_hat, omega_hat, kx, ky); J1 dealias_2by3(J1, Nx, Ny); k1 -J1 .* EF.^0.5; % 半子步积分因子 % 其余三个子步按完全相同方式更新 t_sub, omega_hat, 计算 J2~J4 % 完整代码见文本附带的 ns2d_fourier.m % 组合子步 omega_hat (omega_hat (k1 2*k2 2*k3 k4)/6) ./ EF; t t dt; step step 1; if mod(step, 100) 0 % 计算能量谱并输出 end end参数说明需要特别注意以下三项。第一K2(1,1)必须保持为 1 再做指数运算否则零频增长到无穷。第二积分因子 EF 在 RK4 不同子步中的取法是关键——全子步用 EF半子步用 EF.^0.5这与 RK4 中间时刻的权重匹配。另一种更简单但精度略低的方案是完整步长一次显式欧拉但那会引入额外的粘性截断误差。第三恢复流函数psi_hat -omega_hat ./ K2这行代码最常见的错误是忘记处理零频更隐蔽的问题是如果 ω̂ 没有在当前时间步前做去混淆psi 的计算也会引入高频污染所以规范顺序永远是去混淆 → 解泊松 → 算雅可比 → 再更新。3.4 初始条件与能量谱诊断一个经典的验证算例是 Taylor-Green 涡TGV它有解析解适合检验谱方法的收敛性% 初始涡量Taylor-Green 涡的涡量场 kx0 1; ky0 1; V0 1; % 速度幅值 omega_hat zeros(Ny, Nx); % 解析涡量场 omega V0 * k0 * 2 * cos(kx0*x) .* cos(ky0*y) % 在谱空间直接构造避免物理空间插值 omega_phys 2 * V0 * cos(kx0 * (kx0 * x)) .* cos(ky0 * y); omega_hat fft2(omega_phys); omega_hat dealias_2by3(omega_hat, Nx, Ny);运行后每过一定步数计算能量谱 E(k)满足 ∑E(k) 总动能这是判断程序是否有 bug 的第一层检测。计算中务必保存能谱的时间序列用来判断计算是否发散如果能谱高频端出现能量回升spectral blocking说明去混淆失效或时间步过大。4. 边界处理、去混淆与多尺度工况的调参4.1 周期边界的适用性评估与网格伸缩Fourier 谱方法天然要求周期边界。实际工程里很多流动接近周期槽道湍流在流向和展向可以采用周期近似混合层、尾流的远场也常用周期域模拟。但对于壁面边界傅里叶基不满足无滑移条件这时需要换 Chebyshev-Galerkin 或切比雪夫配点法。以时间步长大小的粗调为例令 U 1、Δx 2π/128 ≈ 0.049CFL 数取 0.5则 Δt ≈ 0.024。对于粘性项ν1e-4 时显式限制约 Δt Δx²/(4ν)换算下来相对宽裕。但雷诺数 Re 1/ν 到达 1e5 量级时高波数模态仍然会因为非线性而产生真正的物理演化这时需要用自适应时间步。实现上可以在每个 RK4 子步后计算一次最大流速动态调整 dtUmax max(abs(u_phys(:))) max(abs(v_phys(:))); dt_new 0.3 * (2*pi/Nx) / Umax; % 取安全系数 0.3这里的逻辑是谱方法每步误差与最大流速和网格间距成正比谱方法没有数值扩散来“压制”误差所以时间步必须取保守值。实际经验是把安全系数设为 0.20.5比有限差分方法常用值低一个档次。4.2 去混淆的正确打开方式什么时候该做、什么时候可以不做初次接触谱方法的人容易陷入一个误区以为只要在初始时刻做一次去混淆后续就万事大吉。正确的理解是混叠发生在每一步的雅可比乘积中每个时间子步都要重新去混淆。多数教科书写“在物理空间乘满 N 个点后再 fft 作为 N 点频域卷积的近似”这个近似带来的误差在非线性项每次求值时都会重新发生。因此去混淆函数必须放在雅可比计算之后、RK4 组合之前。另一种去混淆是“零填充”法把谱系数扩展 3/2 倍比如 N128 扩到 192 点然后物理空间乘法再截回 N128。零填充法的好处是保留了完整的 −N/2 到 N/2 波数范围信息损失小于 2/3 规则代价是每个方向多付 1.5 倍的 FFT 开销。对二维问题这个开销是 2.25 倍但代码实现只多了几行。我一般建议在做高雷诺数涡量场模拟时用零填充做快速原型验证用 2/3 规则就够了。% 零填充去混淆实现替代之前 dealias_2by3 函数 function [Jhat] dealias_zeropad(Jhat_phys, Nx, Ny) % 输入物理空间乘积扩展到 3N/2 尺寸再截断 J_phys_pad zeros(round(3*Nx/2), round(3*Ny/2)); % 把原场放入左上角注意 fftshift 处理 J_phys_pad(1:Nx, 1:Ny) Jhat_phys; Jhat_pad fft2(J_phys_pad); % 截断波数回 Nx x Ny Jhat Jhat_pad(1:Nx/21, 1:Ny/21); end需要提醒的是这段代码省略了 fftshift 的配套操作直接照搬会出错。正确做法是先fftshift到中心化波数再做填充或者用 ifftshift 对齐。原因是 fft2 输出的频段排列是 0 正频、负频扩展矩阵的左上角并不是波数原点。建议初学者不要自己实现零填充用 Matlab 的 fftshift / ifftshift 验证对波数索引的影响后再写。4.3 高波数涡量场的初始化和静音区设置对湍流类问题初始涡量场往往用给定能谱 E(k) 的随机场经过“旋转”构造各向同性条件% 给定能量谱构造随机相位初始涡量 rng(0); rand_phase exp(1i * 2*pi * rand(Ny, Nx)); % 能量谱 E(k) A * k^alpha * exp(-beta*k^2) k_mag sqrt(K2); E_k 10 * k_mag.^2 .* exp(-k_mag.^2 / 1.5); omega_hat sqrt(E_k / (2*pi)) .* rand_phase; omega_hat dealias_2by3(omega_hat, Nx, Ny);设置随机初始场后常见的问题是能量在开头几百步内剧烈下降这不是数值错误随机相位场携带很强的压力和伪“声学”分量它们会在一开始的几个时间步内被快速衰减掉。物理上有意义的湍流统计要等初始过渡结束、能谱斜率接近 k^−3 之后再开始记录数据。如果能量下降幅度超过初始总量的 20%说明初始相位生成的涡量场频率占比过高典型原因是 E_k 的指数项参数太小高频能量过多。4.4 Matlab 语法的版本兼容与运行效率提示上面代码只用了 fft2、ifft2、meshgrid、列行向量扩展这些从 R2016b 开始就稳定的语法。Matlab 2026b 及 r2023b 都能直接运行。如果你只在 CPU 上做 128² 网格Fortran 和 Matlab 的性能差距不大——FFT 是库函数瓶颈在内存带宽。真正需要优化时把时间循环写成 parfor 困难很大因为每个步长依赖前一步更现实的优化是预先分配所有存储在循环外对 psi_hat、J1~J4 等做 zeros并关闭屏显。5. 验证与排错先诊断后调参的三个测试写完代码后不要急着跑大算例。第一步用 Taylor-Green 涡的解析解验证时间推进的正确性解析涡量场每个时刻的衰减是 ω(x,y,t) ω0 · exp(−2νk0²t)。跑 100 步后对比全域 L2 误差。如果误差随 Δt 呈四阶下降RK4 预期说明时间推进逻辑正确。第二步做网格收敛测试同一算例在 64²、128²、256² 运行相同物理时刻比较涡量场 L∞ 误差。谱方法的误差随 N 增长应比有限差分更快地下降。第三步验证守恒性检查总涡量零波数分量是否恒定、总动能单调衰减。最后一个实用的调试技巧是保存每个 RK4 子步的 K2 加权能量。若在子步中某一步出现 NaN打印该步的 K2 谱和最大涡量可以定位是 FFT 变换方向错误fft2 和 ifft2 配对反了还是 K2(1,1) 未被保护。把能量谱输出间隔设为 10 步绘制 log-log 图检查是否有高频翘尾是判断去混淆是否有效的终极手段。这一步做完再上大网格能省掉大量排错时间。本文还有配套的精品资源点击获取