[问题讨论]使用Python学习CFD初级理论系列一二维线性对流（8/10）

Oler 2019年4月28日浏览：2070

前面的案例大多数是一维的问题，从现在开始我们进入二维的世界。

事实上将一维问题扩展到二维甚至三维都是非常简单的，采用相同的思路。在2D空间中，结构网格可定义为：

[问题讨论]使用Python学习CFD初级理论系列一二维线性对流（8/10）的图1

注意这里所提到的结构网格，我们在后面还会详细介绍。

因此，可定义一阶差分格式：

[问题讨论]使用Python学习CFD初级理论系列一二维线性对流（8/10）的图2

下面来处理二维线性对流方程。

1
二维线性对流模型

二维线性对流控制方程为：

[问题讨论]使用Python学习CFD初级理论系列一二维线性对流（8/10）的图3

这里时间项采用向前差分，空间项采用向后差分，离散方程可写成以下格式：

[问题讨论]使用Python学习CFD初级理论系列一二维线性对流（8/10）的图4

式中，i为x方向角标，j为y方向角标，n为时间项角标。

可得待求项：

[问题讨论]使用Python学习CFD初级理论系列一二维线性对流（8/10）的图5

采用初始条件：

[问题讨论]使用Python学习CFD初级理论系列一二维线性对流（8/10）的图6

边界条件：

[问题讨论]使用Python学习CFD初级理论系列一二维线性对流（8/10）的图7

2
Python代码

先用代码将初始条件和边界条件表达出来。

from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

nx = 81 #x方向网格数量
ny = 81 #y方向网格数量
nt = 100 #时间步
c = 1 #常数
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
sigma  = 0.2
dt = sigma * dx

x = np.linspace(0,2,nx)
y  = np.linspace(0,2,ny)
u = np.ones((ny,nx))
un = np.ones((ny,nx))
u[int(0.5/dy):int(1/dy+1), int(0.5/dx):int(1/dx+1)] = 2

# 绘制初始条件
fig = plt.figure(figsize=(12,8))
ax = fig.gca(projection='3d')
x,y = np.meshgrid(x,y)
surf = ax.plot_surface(x,y,u[:],cmap = cm.viridis)
plt.show()

初始条件如下图所示。

[问题讨论]使用Python学习CFD初级理论系列一二维线性对流（8/10）的图8

下面开始迭代计算。我们可以分别采用循环和数组操作来实现，自己体会他们计算时间上的区别。

for n in range(nt + 1): ##loop across number of time steps
   un = u.copy()
   row, col = u.shape    for j in range(1, row):
       for i in range(1, col):
           u[j, i] = (un[j, i] - (c * dt / dx * (un[j, i] - un[j, i - 1])) -
                                 (c * dt / dy * (un[j, i] - un[j - 1, i])))
           u[0, :] = 1
           u[-1, :] = 1
           u[:, 0] = 1
           u[:, -1] = 1

fig = plt.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')
x,y = np.meshgrid(x,y)
surf2 = ax.plot_surface(x, y, u[:], cmap=cm.viridis)
plt.show()

计算结果如下图所示。

[问题讨论]使用Python学习CFD初级理论系列一二维线性对流（8/10）的图9

利用for循环计算速度很慢，下面改用数组运算试试。

for n in range(nt + 1): ##loop across number of time steps
   un = u.copy()
   u[1:, 1:] = (un[1:, 1:] - (c * dt / dx * (un[1:, 1:] - un[1:, :-1])) -
                             (c * dt / dy * (un[1:, 1:] - un[:-1, 1:])))
   u[0, :] = 1
   u[-1, :] = 1
   u[:, 0] = 1
   u[:, -1] = 1
fig = plt.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')
x,y = np.meshgrid(x,y)
surf2 = ax.plot_surface(x, y, u[:], cmap=cm.viridis)
plt.show()

相同的计算结果，如下图所示。

[问题讨论]使用Python学习CFD初级理论系列一二维线性对流（8/10）的图10