[问题讨论]使用Python学习CFD初级理论系列一一维扩散方程（5/10）

Oler 2019年4月27日浏览：2823

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

一维扩散方程为：

[问题讨论]使用Python学习CFD初级理论系列一一维扩散方程（5/10）的图1

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

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

考虑泰勒展开式：

[问题讨论]使用Python学习CFD初级理论系列一一维扩散方程（5/10）的图2

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

[问题讨论]使用Python学习CFD初级理论系列一一维扩散方程（5/10）的图3

改变排列顺序，可得到：

[问题讨论]使用Python学习CFD初级理论系列一一维扩散方程（5/10）的图4

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

[问题讨论]使用Python学习CFD初级理论系列一一维扩散方程（5/10）的图5

改变方程形式，可得到：

[问题讨论]使用Python学习CFD初级理论系列一一维扩散方程（5/10）的图6

假设初始条件：

[问题讨论]使用Python学习CFD初级理论系列一一维扩散方程（5/10）的图7

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

[问题讨论]使用Python学习CFD初级理论系列一一维扩散方程（5/10）的图8

得到：

[问题讨论]使用Python学习CFD初级理论系列一一维扩散方程（5/10）的图9

这里的 σ就是Courant数。

可写代码如下：

import numpy as np
import matplotlib.pyplot as plt
nx = 41 # 网格节点数
dx = 2/(nx - 1) # 网格尺寸
nt = 40       # 时间步数
nu = 0.3 # 扩散系数
sigma = 0.2  #定义一个中间变量，用于确定时间步长dt
dt = sigma * dx**2 / nu #时间步长
# 定义初始条件
u = np.ones(nx)
u[int(0.5 / dx):int(1 / dx + 1)] = 2
un = np.ones(nx) # 定义一个定义数组
plt.ion()
plt.show()
for n in range(nt): #时间循环，计算每一个时间点上网格点上的数据
   plt.cla()
   un = u.copy()
   for i in range(1,nx-1): #空间循环，计算u[i]
       u[i] = un[i]+ nu *dt /dx**2 * (un[i+1] -2 * un[i] + un[i-1])
   # 下面为图形绘制代码
   plt.plot(np.linspace(0,2,nx) , u , 'm' , lw =3 )  #绘制速度随时间变化曲线
   plt.ylim(0.9,2.1)          # 确定y轴上下限
   plt.title("times: %.4f s"% (n*dt)) # 将时间打印到图形上
   plt.pause(0.1) #控制图形显示时间
plt.ioff()
plt.show()