问题讨论：使用Python学习CFD初级理论之一维Burgers方程（6/10）

import numpy as np
import sympy

x, nu, t = sympy.symbols('x nu 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 = phi.diff(x)

from sympy.utilities.lambdify import lambdify
# 得到u的表达式
u = -2 * nu * (phiprime / phi) +4
# 定义lambdify函数，将其写成python可计算的形式
ufunc = lambdify((t,x,nu),u)

from matplotlib import pyplot

nx = 101  # 网格数量
nt = 100 #数据步数
dx = 2 * numpy.pi / (nx - 1) # 网格尺寸
nu = 0.07  # 粘性系数
dt = dx * nu # 时间步长,这么算实际上是为了满足CFL条件

x = np.linspace(0, 2 * numpy.pi, nx) # 计算区域0~2 pi
un = np.empty(nx)
t = 0

u = np.asarray([ufunc(t, x0, nu) for x0 in x]) # 定义初始值

plt.plot(x,u,lw=3,color='red')
plt.xlim([0,2*np.pi])
plt.xlabel("x(m)")
plt.ylabel("time(s)")
plt.ylim([0,10])
plt.show()

for n in range(nt):
   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]

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

plt.figure(figsize=(8,4.8))
plt.plot(x,u, marker='o', lw=3, label='Computational')
plt.plot(x, u_analytical, label='Analytical',lw=3)
plt.xlim([0, 2 * np.pi])
plt.ylim([0, 10])
plt.xlabel("x(m)")
plt.ylabel("time(s)")
plt.legend()