Julia CFD｜03 CFL条件

在前面对流问题中，我们发现当网格网格数量增加时，数值计算会出现非物理解。

这里先写一个将网格数量作为参数的函数，以此来研究网格数量对计算结果的影响。

初始计算代码如下所示，这里采用固定的时间步长0.025 s。

using PyPlot
using Printf
 
function linearconv(nx)
    dx = 2 / (nx - 1)
    nt = 25
    dt = 0.025
    c = 1
 
    u = ones(nx)
    u[Base.floor(Int, 0.5 / dx):Base.floor(Int, 1 / dx)] .= 2
    un = ones(nx)
    PyPlot.plot(range(0, 2, length=nx), u, color="blue", linewidth=2.0, label="Time = 0.0s")
 
    for n in 1:nt
        un = copy(u)
        for i in 2:nx
            u[i] = un[i] - c * dt / dx * (un[i] - un[i - 1])
        end
    end   
 
    # 绘制结果数据
    PyPlot.plot(range(0, 2, length=nx), u, color="red", linewidth=2.0, label=@sprintf("Time = %.3fs",nt*dt))
    PyPlot.xlabel("x(m)")
    PyPlot.ylabel("velocity(m/s)")
    PyPlot.legend(loc=1)
    PyPlot.title(@sprintf("Number of Grid:%d", nx))
    PyPlot.show()
end
 
linearconv(82)

初始采用41个网格节点进行计算。

linearconv(41) #网格数为41

计算结果如下图所示。

这里用step1的参数计算结果作为参考，这里可以看到，波形存在较大的数值扩散，随时间推移，波形与初始形状偏离严重。理想状态是波形能够保持初始形状在物理空间上推进。

linearconv(61) #网格数为61

计算结果如修图所示。

随着网格数量增加，波形失真在减小。

linearconv(71)

计算结果如下图所示。

波形与初始形状越来越接近。

linearconv(81)

随着网格数量增加，波形失真越来越小，可以认为是计算越来越精确。

linearconv(91)

发生了什么情况，计算已经完全偏离了初始波形，出现了错误。这其中发生了什么呢？要回答这个问题，我们必须考虑一下代码中实际发生了什么。

在时间循环的每一次迭代中，我们使用当前时刻的波的数据来估计下一个时刻中的波的速度。最初，网格点数的增加返回了更准确的答案。数值扩散较少，波形看起来比第一个例子更像方波。

在每个时间步长内，波的传播距离必须小于网格尺寸dx。网格尺寸dx与总网格数量nx有关。因此，稳定计算条件需要满足下式：

其中是波速；称为Courant Number，该数值决定了计算是否稳定。当时，计算会发散。

按前面计算条件，m/s，s，可以计算得到临界m，则临界网格节点数量为，下面看看当网格节点数量为82时的计算结果。

linearconv(82)

计算结果如下图所示

与预料中的情况一直，可以看到计算已经出现了问题。在新版本的代码中，我们采用Courant数来根据dx来计算合适的dt。

下面改写代码：

using PyPlot
using Printf
 
function linearconv(nx, nt=25)
    c = 1
    dx = 2 / (nx - 1)
    # nt = 25
    sigma = 0.5
    dt = sigma * dx / c
 
    u = ones(nx)
    u[Base.floor(Int, 0.5 / dx):Base.floor(Int, 1 / dx)] .= 2
    un = ones(nx)
    PyPlot.plot(range(0, 2, length=nx), u, color="blue", linewidth=2.0, label="Time = 0.0s")
 
    for n in 1:nt
        un = copy(u)
        for i in 2:nx
            u[i] = un[i] - c * dt / dx * (un[i] - un[i - 1])
        end
    end   
 
    # 绘制结果数据
    PyPlot.plot(range(0, 2, length=nx), u, color="red", linewidth=2.0, label=@sprintf("Time = %.3fs",nt * dt))
    PyPlot.xlabel("x(m)")
    PyPlot.ylabel("velocity(m/s)")
    PyPlot.legend(loc=1)
    PyPlot.title(@sprintf("Number of Grid:%d", nx))
    PyPlot.show()
end

测试一下新的程序代码：

linearconv_m(82,25)

可以看到在网格数量为82时未出现问题，由于Courant数取为0.5，导致时间步长为只有原来时间步长的一半，在时间步长不变的情况下总时间也减半。这里将时间步数提高到50再重新计算。

linearconv(82,50)

结果如下图所示，得到了与前面网格数量41时相同的计算结果。

改造一下代码，将时间步数和网格数量均设置为1000。

linearconv(1000,1000)

嗯，满足了CFL条件后，不管网格怎么加密，计算结果依然不会崩溃。在Fluent中的Run Calculation面板中就有利用CFL进行时间步长调整的搞法。