Julia CFD｜05 一维Burgers方程

本次利用Julia求解一维Burger方程。关于Burgers方程的具体描述，可以参阅维基百科。一维Burgers方程描述为：

Burger方程中同时包含了对流项与扩散项，式中u为速度，为介质粘度。

对时间项采用向前差分，对空间项采用向后差分，二阶导数项采用中心差分，可写成离散方程为：

将待求项提出来，可写成：

采用边界条件为：

这是一个周期边界。

本次采用的初始条件为：

其解析形式为：

可以利用Sympy包进行符号运算，类似Mathematica软件。

故采用Sympy进行简化。

using SymPy
 
t = Sym("t")
x, nu = symbols("x,nu", real=true)
# 定义phi的表达式
phi = exp(-(x - 4 * t)^2 / (4 * nu * (t + 1))) + exp(-(x - 4 * t - 2 * pi)^2 / (4 * nu * (t + 1)))
# 计算phi的偏导数
phiprime = diff(phi, x)
# 定义初始速度分布表达式
u = -2 * nu * phiprime / phi + 4
# 定义lambdify函数，将其改写为julia可以显式调用的函数
ufunc = lambdify(u, (t, x, nu))

下面定义初始条件。

using PyPlot
 
nx = 101  # 网格数量
nt = 100 # 时间步数
dx = 2 * pi / (nx - 1) # 网格尺寸
nu = 0.07  # 粘性系数
dt = dx * nu # 时间步长,这么算实际上是为了满足CFL条件
 
x = range(0, 2 * pi, length=nx) # 计算区域0~2 pi
un = zeros(nx)
t = 0
 
u = [ufunc(t, x0, nu) for x0 in x] # 定义初始值
PyPlot.plot(x,u,lw=2,color="red")
PyPlot.xlim([0,2 * pi])
PyPlot.xlabel("x(m)")
PyPlot.ylabel("u(m/s)")
PyPlot.ylim([0,10])
PyPlot.show()

初始速度分布如图所示。

回到Burgers方程上来。

for n in 1:nt
    un = copy(u)
    for i in 2:nx - 1
        u[i] = un[i] - un[i] * dt / dx * (un[i] - un[i - 1]) + nu * dt / dx^2 * (un[i + 1] - 2 * un[i] + un[i - 1]);
    end
    u[1] = un[1] - un[1] * dt / dx * (un[1] - un[end - 1]) + nu * dt / dx^2 * (un[2] - 2 * un[1] + un[end - 1]);
    u[end] = un[end] - un[end] * dt / dx * (un[end] - un[end - 1]) + nu * dt / dx^2 * (un[1] - 2 * un[end] + un[end - 1]);
end
 
u_analytical = [ufunc(nt * dt, xi, nu) for xi in x];
# figure(figsize=(11, 7), dpi=100)
PyPlot.plot(x,u, marker="o", lw=2, label="Computational")
PyPlot.plot(x, u_analytical, label="Analytical")
PyPlot.xlim(0,2 * pi)
PyPlot.ylim([0,10])
PyPlot.xlabel("x(m)")
PyPlot.ylabel("u(m/s)")
PyPlot.legend();
PyPlot.show()

计算结果如图所示。

可以看到，误差很大，我们在后续的文章中会提到如何改进算法以减小误差。

完整代码如下：

using SymPy
 
t = Sym("t")
x, nu = symbols("x,nu", real=true)
# 定义phi的表达式
phi = exp(-(x - 4 * t)^2 / (4 * nu * (t + 1))) + exp(-(x - 4 * t - 2 * pi)^2 / (4 * nu * (t + 1)))
# 计算phi的偏导数
phiprime = diff(phi, x)
# 定义初始速度分布表达式
u = -2 * nu * phiprime / phi + 4
# 定义lambdify函数，将其改写为julia可以显式调用的函数
ufunc = lambdify(u, (t, x, nu))
 
 
using PyPlot
nx = 101  # 网格数量
nt = 100 # 时间步数
dx = 2 * pi / (nx - 1) # 网格尺寸
nu = 0.07  # 粘性系数
dt = dx * nu # 时间步长,这么算实际上是为了满足CFL条件
 
x = range(0, 2 * pi, length=nx) # 计算区域0~2 pi
un = zeros(nx)
t = 0
 
u = [ufunc(t, x0, nu) for x0 in x] # 定义初始值
# 下面几行代码用于显示初始值
# PyPlot.plot(x,u,lw=2,color="red")
# PyPlot.xlim([0,2 * pi])
# PyPlot.xlabel("x(m)")
# PyPlot.ylabel("u(m/s)")
# PyPlot.ylim([0,10])
# PyPlot.show()
 
for n in 1:nt
    un = copy(u)
    for i in 2:nx - 1
        u[i] = un[i] - un[i] * dt / dx * (un[i] - un[i - 1]) + nu * dt / dx^2 * (un[i + 1] - 2 * un[i] + un[i - 1]);
    end
    # 下面两行用于处理边界值
    u[1] = un[1] - un[1] * dt / dx * (un[1] - un[end - 1]) + nu * dt / dx^2 * (un[2] - 2 * un[1] + un[end - 1]);
    u[end] = un[end] - un[end] * dt / dx * (un[end] - un[end - 1]) + nu * dt / dx^2 * (un[1] - 2 * un[end] + un[end - 1]);
end
 
u_analytical = [ufunc(nt * dt, xi, nu) for xi in x];
PyPlot.plot(x,u, marker="o", lw=2, label="Computational")
PyPlot.plot(x, u_analytical, label="Analytical")
PyPlot.xlim(0,2 * pi)
PyPlot.ylim([0,10])
PyPlot.xlabel("x(m)")
PyPlot.ylabel("u(m/s)")
PyPlot.legend();
PyPlot.show()