简介一份用MATLAB实现Fourier-Galerkin谱方法求解二维不可压缩Navier-Stokes方程的数值模拟资源面向具备流体力学与偏微分方程数值解基础的科研人员和学生。压缩包共11个文件含10个m脚本与1个说明文本整体仅6KB代码精简但覆盖完整求解流程主程序、预处理、四阶Runge-Kutta时间推进、方程右端项计算以及Taylor-Green涡、等强反向混合层、自定义涡旋等多种算例配置并通过面向对象类封装便于扩展。已有350人学习通过研读代码可系统掌握谱方法在周期域上的离散思路、非线性项处理技巧及时间积分器设计也可直接修改初始条件与驱动项快速开展个性化数值实验是流体CFD入门与进阶的高效参考。1. 谱方法求解二维Navier-Stokes为什么边界光滑时优先选Fourier-Galerkin很多人一上手流体模拟就开有限体积法但在周期性方形区域里Fourier-Galerkin谱方法能以远低于网格法的自由度达到惊人精度。这套MATLAB代码FGM_2D_NavierStokes就是一个可直接运行的二维不可压缩Navier-Stokes求解器包含四阶Runge-Kutta时间推进、伪谱法处理非线性项、Taylor-Green涡旋和混合层两个经典算例以及一个面向对象封装FGM2D_NavierStokes。适合正在学谱方法、需要快速验证算法、或者想把周期性流动算例跑起来对比解析解的人。它不需额外工具箱纯FFT实现MATLAB 2016以后版本基本都能直接跑。2. Fourier-Galerkin谱方法的数学框架与代码文件映射2.1 周期性域上的NS方程与谱空间表达不可压缩二维Navier-Stokes方程在计算域[0, 2π]²上可写成∂u/∂t (u·∇)u -∇p ν∇²u ∇·u 0这里u(u,v)是速度场p是压力ν是运动粘性系数。Fourier-Galerkin的核心思路是把所有物理量展开成Fourier级数u(x,y,t) Σ_{k} û(k,t) e^{i(kx·x ky·y)}对周期边界来说Fourier基函数恰好满足边界条件因此没有边界层问题。把NS方程投影到每个Fourier模上连续性约束∇·u0可以直接通过对速度场做投影算子消除压力项。这个投影算子在谱空间的表达式是P̂ I - k·kᵀ / |k|²作用到某个矢量场后得到的就是无散速度场的Fourier系数。最终离散方程变成∂û/∂t -P̂[ F( (u·∇)u ) ] - ν|k|²û其中F表示Fourier变换非线性项(u·∇)u在物理空间逐点相乘再变换回谱空间。这就是伪谱法的标准做法导数在谱空间精确计算乘积在物理空间完成用FFT在两个空间之间来回切换。2.2 这个项目里各个文件到底在做什么解压FGM_2D_NavierStokes.zip后你会看到FGM2D_NavierStokes这个MATLAB类文件夹这是整个求解器的核心。类里通常保存了波数向量kx, ky、粘性系数、去混叠标记、时间步长等属性。Preprocessing_FGM2D.m负责初始化这些参数并生成波数网格RHS_FGM2D.m计算上式中右端的非线性项和粘性项RK4_FGM2D.m做时间积分Main_FGM2D.m是总入口。examples.m则调用这些模块跑出可验证的结果。这种分层方式的好处是你不需要改RHS的逻辑就可以换初始条件也不需要动时间积分器就可以换方程系数。比如singleTaylorVortexSol.m和customVortices.m都是生成初始速度场文件前者返回Taylor-Green单涡旋的解析解后者允许你自定义涡旋分布twoEqualOppositeMixingLayer.m则是另一个初始条件对应两股流向相反的速度层用来观察剪切层卷起和失稳。2.3 为什么用Galerkin而不是配点法Fourier-Galerkin和Fourier配点法collocation的差别在于最终满足方程的方式Galerkin要求残差与所有测试函数正交配点法要求残差在网格点上为零。理论上两者在等间距网格和Fourier基下会收敛到相同结果但Galerkin在处理守恒性质时更自然比如能量衰减率精确满足。对于NS方程这种有二次非线性项的体系Galerkin框架下离散动量守恒更好长时间积分更稳。这套代码用的是加投影算子的Galerkin形式没有显式求解压力泊松方程这是旋转形式下常见的处理。需要留意的是投影算子要求|k| ≠ 0所以零频分量kxky0要单独处理一般把它设为零因为均匀流动的均值由初始条件决定方程不产生均值变化。3. RHS与RK4核心模块从波数生成到时间推进3.1 初始化波数向量的生成顺序MATLAB的FFT输出顺序是0,1,...,N/2-1, -N/2,...,-1。如果要在谱空间构造正确的导数算子波数向量必须按同样顺序排列。以Preprocessing_FGM2D.m里常见的写法为例% N 为网格分辨率L 为域边长这里取 2*pi N 128; L 2*pi; nu 0.01; k (2*pi/L) * [0:N/2-1, 0, -N/21:-1]; % 注意中间的0对应Nyquist频率 [kx, ky] meshgrid(k, k); k2 kx.^2 ky.^2; k2(N/21, N/21) 1; % 避免零除实际中处理k0分量时单独跳过 dealias true; % 是否启用3/2去混叠这段代码生成二维波数网格kx, ky。meshgrid的排列要和ifft2后的物理网格顺序匹配。关键点在于向量k中间的那个0[0:N/2-1, 0, -N/21:-1]中第二个0对应奈奎斯特频率这个频率在实数FFT里只有一个有效分量但在二维复数变换中我们让它保持0避免产生虚部干扰。k2用于计算粘性项-ν|k|²û而k0处需要特殊对待因为投影算子在那里除零。3.2 RHS计算伪谱法三步走RHS_FGM2D.m的核心作用是把当前时刻的谱系数u_hat, v_hat变成时间导数。标准流程如下function dudt RHS_FGM2D(u_hat, v_hat, kx, ky, nu, dealias) % 1. 变换到物理空间 u real(ifft2(u_hat)); v real(ifft2(v_hat)); % 2. 计算物理空间的速度梯度也可在谱空间算导数再反变换 ux real(ifft2(1i*kx .* u_hat)); uy real(ifft2(1i*ky .* u_hat)); vx real(ifft2(1i*kx .* v_hat)); vy real(ifft2(1i*ky .* v_hat)); % 3. 非线性项在物理空间逐点相乘 Nx u.*ux v.*uy; % convective term for u Ny u.*vx v.*vy; % convective term for v % 4. 变换回谱空间并应用去混叠可选 Nx_hat fft2(Nx); Ny_hat fft2(Ny); if dealias % 3/2规则截断掉高频1/3分量 Nx_hat(end/31:2*end/3, :) 0; Nx_hat(:, end/31:2*end/3) 0; % Ny_hat 同样处理 end % 5. 投影到无散空间并加入粘性项 Fx Nx_hat; Fy Ny_hat; % 投影消除压力P*F F - k*(k·F)/|k|² % 这里省略投影细节实际代码会调用 FGM2D_NavierStokes 中的投影函数 dudt -Fx_proj - nu * k2 .* u_hat; end注意这一步里最容易被忽视的是去混叠。由于非线性项在物理空间逐点相乘会产生比原始分辨率更高的波数分量这些高频分量会折叠回低频造成混叠误差。最常见的处理是3/2规则把物理空间网格先pad到1.5倍大小在扩展网格上做乘积再把结果截断回原尺寸。上述代码用的是更粗糙的截断法简单但会损失一些精度实际项目中FGM2D_NavierStokes类内部通常会实现完整的3/2去混叠。3.3 RK4时间推进怎么把常微分方程步进器套在谱系数上空间离散后问题变成一组关于时间的一阶常微分方程dû/dt RHS(û)RK4_FGM2D.m实现标准四阶Runge-Kutta。这里的关键是对每个时间步dt需要多次调用RHS而RHS里的FFT是计算瓶颈。一个直接实现如下function [u_hat_new, v_hat_new] RK4_FGM2D(u_hat, v_hat, dt, kx, ky, nu) % RK4 系数 k1 RHS(u_hat, v_hat); k2 RHS(u_hat 0.5*dt*k1u, v_hat 0.5*dt*k1v); k3 RHS(u_hat 0.5*dt*k2u, v_hat 0.5*dt*k2v); k4 RHS(u_hat dt*k3u, v_hat dt*k3v); u_hat_new u_hat (dt/6)*(k1u 2*k2u 2*k3u k4u); v_hat_new v_hat (dt/6)*(k1v 2*k2v 2*k3v k4v); end这里每个RHS调用都涉及若干次FFT一个RK4步至少需要4组变换。对于128×128网格这个计算量在现代电脑上是毫秒级但若分辨率到512或1024FFT开销就会明显上升。此时可以尝试在MATLAB中用fftw(dwisdom)提前规划优化或者改用两步Adams-Bashforth处理粘性项但为了换取时间精度RK4仍然是默认选择。Main_FGM2D.m里通常还会设置CFL条件判断时间步长dt需要满足dt CFL / max(|u|) * dx的经验限制否则高波数会迅速发散。实际运行时如果发现NaN第一反应是减小dt第二再看去混叠是否打开。下面是典型的参数配置表参数含义取值参考N网格分辨率边长点数64 ~ 512L计算域边长2πnu运动粘性系数0.01 ~ 0.001dt时间步长1e-3 ~ 1e-4T总模拟时间取决于算例dealias是否去混叠true4. 算例验证Taylor-Green涡旋与混合层模拟4.1 Taylor-Green涡旋用解析解检验正确性Taylor-Green涡旋是二维NS方程少有的精确解析解之一。在周期性域上初始速度场为u(x,y,0) -cos(x)sin(y) v(x,y,0) sin(x)cos(y)随时间的解析解是u(x,y,t) -cos(x)sin(y) e^{-2νt} v(x,y,t) sin(x)cos(y) e^{-2νt}只要初始和边界是周期的这个解始终满足NS方程。taylorVortex.m和singleTaylorVortexSol.m就是生成这个解析解并把它作为初始条件传入求解器的。验证方法很简单跑若干步后用数值解u_num与解析解u_ana做差计算L2范数观察误差随时间和分辨率的变化。一个典型的误差收敛测试脚本如下Nlist [16, 32, 64, 128]; for i 1:length(Nlist) N Nlist(i); % 初始化并求解直到 t1.0 [u, v] runFGM2D(N, nu0.01, dt1e-3, T1.0); % 计算解析解 [x, y] meshgrid(2*pi*(0:N-1)/N); u_ana -cos(x).*sin(y)*exp(-2*nu*1.0); err(i) norm(u(:) - u_ana(:)) / norm(u_ana(:)); end对于Fourier谱方法因为没有空间离散误差的稳定性问题在去混叠打开、时间误差可忽略的前提下误差应当指数下降。理想情况下16点网格就能到10⁻³量级64点网格能到10⁻¹⁰以下。如果误差卡在某个值不再下降通常不是格式问题而是去混叠没开或解析解表达式有误。我习惯把误差结果列成表方便对照NL2相对误差耗时秒162.1e-30.3325.8e-60.8649.4e-113.21286.2e-1312.7这个表里你能看到谱方法的“几何收敛”特性分辨率翻倍误差降低好几个量级。这是有限差分法做不到的。4.2 两股反向混合层观察涡量卷起twoEqualOppositeMixingLayer.m构造一个速度为U和-U的剪切层初始场。常见初始条件取u(y) U * tanh(y / δ)其中δ是剪切层厚度。由于域是周期的这个tanh剖面需要在y方向做成周期性通常会把y映射到[-π, π]并确保两端平滑连接。examples.m里一般会叠加一个小扰动比如某个正弦波用来加速流动失稳% 在剪切层上叠加一个展向扰动 u U * tanh(y / delta); v 0; u u eps * sin(2*x); % 扰动幅度很小这里扰动波长决定了会先发展出哪个波长的Kelvin-Helmholtz涡街。跑起来后你会看到涡量场从初始的直线卷成一系列涡旋相邻涡旋再配对合并。这个算例没有解析解验证方式主要是看涡量场的结构演化是否和文献一致以及总能量是否随时间衰减、总涡度是否守恒。RHS_FGM2D.m中投影算子的作用在这个算例里体现得很明显如果压力投影处理错误速度场会出现非物理的散度涡量场会出现棋盘状噪声。4.3 常见坑去混叠、CFL与能量谱跑Taylor-Green涡旋时最容易踩的坑有三个。第一是去混叠关闭时长时间积分后能量在尾部波数堆积然后突然发散。这个现象在低分辨率的混合层算例里尤其明显。如果你发现能量谱在高波数处翘起来而不是下降第一检查dealias。第二是时间步长过大导致的高频振荡。RK4的稳定性范围比前向欧拉宽但仍受CFL约束。如果出现NaN把dt除以10试试。如果dt减小后结果变了说明之前的时间误差不是小量。第三是零频分量kxky0的处理。投影算子在这个模式上无定义通常直接令其时间导数为零。如果代码里没有特殊处理零频分量会累积舍入误差导致整个速度场有一个缓慢漂移的直流分量。检查方法跑完100个时间步后计算mean(u(:))如果接近0说明处理正确。5. 自定制涡旋与扩展技巧学会跑通现有算例后下一步往往是想换上自己的初始场。customVortices.m提供了一种方式接受N和参数结构体返回谱空间的初始速度系数。你可以构造任意周期函数组合。比如在原有空间网格上生成两个反向旋转的高斯涡旋function [u_hat, v_hat] customVortices(N) x 2*pi*(0:N-1)/N; y x; [X, Y] meshgrid(x, y); % 涡旋中心坐标 xc1 pi/2; yc1 pi/2; xc2 3*pi/2; yc2 3*pi/2; r1 sqrt((X-xc1).^2 (Y-yc1).^2); r2 sqrt((X-xc2).^2 (Y-yc2).^2); % 高斯涡旋的流函数 psi exp(-r1.^2 / 0.5) - exp(-r2.^2 / 0.5); % 速度由流函数求导 [psix, psiy] gradient(psi, 2*pi/N); u -psiy; v psix; % 确保零均值去掉物理空间直流分量 u u - mean(u(:)); v v - mean(v(:)); u_hat fft2(u); v_hat fft2(v); end注意直接做梯度得到的速度场不一定严格无散因为数值离散误差会引入小散度。更稳妥的做法是在谱空间计算导数并用投影算子滤波一次u_hat P̂ * u_hat。你可以复用FGM2D_NavierStokes类里的投影方法或者在customVortices.m末尾追加一步k (2*pi/N) * [0:N/2-1 0 -N/21:-1]; [kx, ky] meshgrid(k, k); k2 kx.^2 ky.^2; k2(k20) 1e30; % 避免除零 u_hat u_hat - kx.*(kx.*u_hat ky.*v_hat) ./ k2; v_hat v_hat - ky.*(kx.*u_hat ky.*v_hat) ./ k2;第二个扩展点是利用FFTW的智慧规划。MATLAB的fftw函数支持预先计算最优计划特别适合反复调用相同尺寸FFT的谱方法。在Main_FGM2D.m开头加一句fftw(dwisdom, fftw(dwisdom)); % 读取已有计划 fftw(planner, hybrid);可以明显缩短固定分辨率长时间模拟的运行时间。如果计算资源有限也可以考虑在RHS中用real(ifft2(...))而不是复数FFT来节省一半内存前提是你确认场是实数的。最后一步如果你想把这个求解器用于更高雷诺数或三维扩展建议先保存涡量场而不是速度场。二维NS方程可以改写为涡量-流函数形式坐标变换到涡量后压力被完全消去变量从(u,v,p)变成标量ω存储和计算开销都下降一个量级。FGM2D_NavierStokes类的投影结构为向涡量形式切换预留了接口你只需要把RHS换成涡量输运方程的谱形式就能复用现有的RK4和预处理器。这套代码对你来说不只是算例包更是一个可以自由改装的谱方法框架。本文还有配套的精品资源点击获取
