基于Python的CFD编程入门（4）伯格斯方程(Burgers’ Equation)

伯格斯方程(Burgers’ Equation)或称为贝特曼-伯格斯方程(Bateman–Burgers )是一个对流-扩散方程，出现在应用数学的各个领域，如流体力学、非线性声学、空气动力学等。

01

理论基础

一维伯格斯方程：

可以看到伯格斯方程由之前介绍过的一维非线性对流和扩散项所组成。

因此一维非线性对流方程也可称为无粘伯格斯方程。

因此可以使用前文同样的离散方法，对时间导数使用前向差分，对空间导数使用后向差分，对二阶导数使用中心差分。

调整格式：

令初始条件为：

在周期性边界条件：

伯格斯方程存在解析解：

02

代码实现

首先使用sympy库符号运算求解方程的解析解

import numpy as np
import sympy
from sympy.utilities.lambdify import lambdify
from matplotlib import pyplot as pltx = sympy.Symbol('x')
nu = sympy.Symbol('nu')
t = sympy.Symbol('t')
# 定义phi的表达式
phi = sympy.exp(-(x - 4 * t)**2 / (4 * nu * (t + 1))) + \
      sympy.exp(-(x - 4 * t - 2 * sympy.pi)**2 / (4 * nu * (t + 1)))
# 计算phi偏导数
phiprime = sympy.diff(phi,x)


u = -2 * nu / phi * phiprime + 4

ufunc = lambdify((t,x,nu),u)

# 定义初始条件
nx = 101 # 网格数量
nt = 100 # 时间步数
dx = 2 * np.pi / (nx -1) # 网格尺寸
nu = 0.04 # 粘性系数
dt = dx * nu # 时间步长

x = np.linspace(0,2*np.pi,nx)
un = np.empty(nx)
t = 0
u = np.asarray([ufunc(t, x0, nu) for x0 in x]) # 定义初始速度

u_analytical = np.asarray([ufunc(nt * dt, xi, nu) for xi in x]) # 准确值



得到解析解的图像：

接着我们求解数值解：

plt.ion() 
for n in range(nt):
    plt.cla()
    un = u.copy()
    for i in range(1, 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])
    u[0] = un[0] - un[0] * dt / dx * (un[0] - un[-2]) + nu * dt / dx ** 2 * (un[1] - 2 * un[0] + un[-2]) # 计算第一个点
    u[-1] = u[0]
    plt.plot(x, u_analytical,label='Analytical')
    plt.plot(x, u,label='Computational')
    plt.legend() 
    plt.title("time: %.4f s" % (n * dt)) 
    plt.pause(0.01) 

plt.ioff()
plt.show()

结果为：

迭代求解过程：

03

误差来源

伯格斯方程数值解与解析解的差异来源主要有以下几点：

• 两种情况下的粘度是一样的，但数值解显示出比解析解更多的时间和空间耗散（主要在时间尺度上，空间尺度相对较小）。

• 数值耗散是造成解析解和数值解之间差异的主要原因。当物理粘度降低，使物理耗散减小时，数值耗散的影响就会更加明显。
  
  
• 方程离散化过程中存在截断误差。

• 瞬态和对流项的一阶近似分别包含时间和空间上的数值扩散，因此需要对瞬态项使用高于一阶的方法。

—— end ——