人工粘性CFD数值实验

Question

我想用matlab做下面的数值实验，但是发现在雷诺数Re=200，Δx=0.005和Re=2000,Δx=0.005时，也就是第一和第二个图像，仿真的效果相同，但是当Re=2000,Δx=0.0005的时候相差特别多，如下图

中间震荡非常大
下面是我的matlab程序：

%%  人工粘性法
%% Re=200 h=0.005;
clc,clear;
n_points=2001;% 网格数
dom_size=1;%区域大小
h=dom_size/(n_points-1);%网格大小
dt=0.0005;%时间步长
Re=2000;
%% Initialize 
u(n_points)=0;
for i=1:n_points
    if (i-1)*h<=0.5
        u(i)=1;
    end
    if (i-1)*h>0.5
        u(i)=0;
    end
end
unew=u;%时间迭代初始化
iterations=0;%迭代步数
%% Calculation
while iterations<=200
    for i=2:(n_points-1)
        unew(i)=u(i)+dt*(u(i+1)-2*(u(i))+u(i-1))/(Re*(h)^2)-dt*(u(i+1)-u(i-1))/(2*h);
    end
    u=unew;
    iterations=iterations+1;
end
x=0:h:1;
p=plot(x,u);
xlim([0 1]);
ylim([-0.1 1.1])
title(['t=' num2str((iterations-1)*dt) 's']);
xlabel('x');
ylabel('u');
p.LineWidth=1.5;
ax=gca;
ax.FontSize=18;
请各位大佬指正小弟，谢谢啦

Answer 1

引自免费微信小程序：皆我百晓生

在你的代码中，观察到的问题可能是因为时间步长（dt）相对于网格大小（h）来说过大，导致了数值不稳定。根据CFL条件（Courant-Friedrichs-Lewy condition），对于稳定性，时间步长应该满足以下关系：

[ \frac{dt}{h} \leq C ]

其中 ( C ) 是一个与方法有关的常数，对于二阶精度的中心差分通常建议 ( C \leq 0.5 )。

为了修正这个问题，你需要调整时间步长以确保满足CFL条件。在你的例子中，已经将网格大小设置为 h = 0.0005 和 h = 0.005，而时间步长设置为 dt = 0.0005。这个时间步长对于较小的网格大小可能是安全的，但对于较大的网格大小（例如 h = 0.005），它可能导致数值不稳定。

针对不同的Reynolds数和网格大小，你应该相应地调整时间步长。这里是一个修改后的版本，我会为每个情况提供一个合适的时间步长，并且会添加一些注释解释这些变化：

%%  人工粘性法
function plot_velocity profiles()
clc, clear;

% 定义变量
n_points = 2001; % 网格数
dom_size = 1; % 区域大小
Re = [200, 2000]; % 不同的Reynolds数
dx_values = [0.0005, 0.005]; % 不同的网格大小
time_end = 0.1; % 时间终点
max_iterations = 200; % 最大的迭代次数

% 遍历不同的Reynolds数和网格大小
for i = 1:length(Re)
    Re_current = Re(i);
    dx_current = dx_values(i);
    
    % 计算CFL限制下的最大允许时间步长
    dt_max = 0.5 * dx_current; % 假设CFL系数为0.5
    
    % 根据目标时间确定实际时间步长
    dt = min(dt_max, time_end / max_iterations);
    
    % 初始化
    u = zeros(n_points, 1);
    for j = 1:n_points
        if (j - 1) * dx_current <= 0.5
            u(j) = 1;
        else
            u(j) = 0;
        end
    end
    unew = u; % 时间迭代初始化
    iterations = 0; % 迭代步数
    
    % 计算
    while iterations < max_iterations
        for j = 2:(n_points - 1)
            unew(j) = u(j) + ...
                dt * (u(j + 1) - 2 * (u(j)) + u(j - 1)) / (Re_current * (dx_current)^2) - ...
                dt * (u(j + 1) - u(j - 1)) / (2 * dx_current);
        end
        u = unew;
        iterations = iterations + 1;
    end
    
    % 绘制结果
    x = linspace(0, dom_size, n_points);
    p = plot(x, u);
    xlim([0 1]);
    ylim([-0.1 1.1]);
    title(sprintf('t=%f s, Re=%d, dx=%f', (iterations - 1) * dt, Re_current, dx_current));
    xlabel('x');
    ylabel('u');
    p.LineWidth = 1.5;
    ax = gca;
    ax.FontSize = 18;
    
    % 显示图形
    pause(0.5); % 暂停一段时间以便查看每一步的结果
end
end

% 调用函数
plot_velocity_profiles();

这个版本的代码会分别对Reynolds数200和2000以及网格大小0.0005和0.005进行计算，并在每次迭代后暂停片刻以查看结果。通过这种方法，你可以观察到不同参数组合下的解的行为。如果你想要连续显示所有结果而不是逐个暂停，只需删除调用中的 pause(0.5) 行即可。