Julia CFD｜04 一维扩散问题

本文描述利用Julia求解计算一维扩散方程。一维扩散方程为：

与前面方程不同的地方在于此方程包含二阶导数，因此首先对二阶导数进行离散。

采用中心差分格式对二阶导数进行离散。

考虑泰勒展开式：

两式相加，可得到二阶导数项：

改变排列顺序，可得到：

将二阶导数项代入到原方程中，可得到离散方程：

改变方程形式，可得到：

假设初始条件：, ； ，

注：这里定义一个sigma的目的是为了控制时间步长dt，以使其满足CFL条件。令：

得到：

这里的就是常说的Courant数。

可写代码如下：

using PyPlot
using Printf
 
nx = 41 # 网格节点数
dx = 2 / (nx - 1) # 网格尺寸
nt = 40       # 时间步数
nu = 0.3 # 扩散系数
sigma = 0.2  # 定义一个中间变量，用于确定时间步长dt
dt = sigma * (dx^2) / nu # 时间步长
 
# 定义初始条件
u = ones(nx)
u[Base.floor(Int, 0.5 / dx):Base.floor(Int, 1 / dx + 1)] .= 2
un = ones(nx) # 定义一个定义数组
 
 
for n in 1:nt # 时间循环，计算每一个时间点上网格点上的数据
    PyPlot.cla()
    un = copy(u)
 
    for i in 2:nx - 1 # 空间循环，计算u[i]
        u[i] = un[i] + nu * dt / dx^2 * (un[i + 1] - 2 * un[i] + un[i - 1])
    end
 
    # 下面为图形绘制代码
    PyPlot.plot(range(0, 2, length=nx), u, color="red", linewidth=3)  # 绘制速度随时间变化曲线
    PyPlot.ylim(0.9, 2.1)          # 确定y轴上下限
    PyPlot.title(@sprintf("times: %.4f s",(n * dt))) # 将时间打印到图形上
    PyPlot.pause(0.2) # 控制图形显示时间
end
PyPlot.show()

代码运行如图所示。

程序中利用参数sigma进行时间步长控制，下面来研究一下sigma取值对计算的影响。

将代码改造为下面的形式。

using PyPlot
using Printf
 
function diffusion(sigma=0.2)
    nx = 101 # 网格节点数
    dx = 2 / (nx - 1) # 网格尺寸
    nt = 40       # 时间步数
    nu = 0.3 # 扩散系数
    dt = sigma * (dx^2) / nu # 时间步长
 
    # 定义初始条件
    u = ones(nx)
    u[Base.floor(Int, 0.5 / dx):Base.floor(Int, 1 / dx + 1)] .= 2
 
    # 定义一个中间临时数组
    un = ones(nx)
 
    for n in 1:nt # 时间循环，计算每一个时间点上网格点上的数据
        PyPlot.cla()
        un = copy(u)
        for i in 2:nx - 1 # 空间循环，计算u[i]
            u[i] = un[i] + nu * dt / dx^2 * (un[i + 1] - 2 * un[i] + un[i - 1])
        end
 
        # 下面为图形绘制代码
        PyPlot.plot(range(0, 2, length=nx), u, color="red", linewidth=3)  # 绘制速度随时间变化曲线
        PyPlot.ylim(0.9, 2.1)          # 确定y轴上下限
        PyPlot.xlabel("x(m)")
        PyPlot.ylabel("u(m/s)")
        PyPlot.title(@sprintf("times: %.4f s",(n * dt))) # 将时间打印到图形上
        PyPlot.pause(0.2) # 控制图形显示时间
    end
    PyPlot.show()
end

增大sigma的值进行计算。

diffusion(0.3)

计算结果如下图所示。

除了时间步长有所增大外，计算过程是稳定的。继续增大sigma值。

diffusion(0.5)

计算结果如下图所示。

计算结果出现了振荡。减小sigma的值为0.49试试。

diffusion(0.49)

计算结果如下图所示。

起始有小幅振荡，后面趋于平稳。将sigma的值调整为0.51试试。

diffusion(0.51)

计算结果如下图所示。

崩了。

看来，sigma=0.5是临界值无疑了。