[系列]Python计算CFD问题<6>：2D线性对流问题
在第一个例子中，我们提到了解一维线性对流问题。在本例中，我们将利用python数值求解二维线性对流问题。
2D线性对流问题可以用此方程进行表达：
[attach]69963[/attach]
时间项采用向前差分，空间项采用先后差分。离散后的方程为：
[attach]69965[/attach]
换成迭代格式为：
[attach]69964[/attach]
初始条件为：
[attach]69966[/attach]
边界条件：
[attach]69967[/attach]
编制计算程序为：
- from mpl_toolkits.mplot3d import Axes3D ##New Library required for projected 3d plots
-
- import numpy as np
- import matplotlib.pyplot as plt
-
- ###variable declarations
- nx = 81
- ny = 81
- nt = 100
- c = 1
- dx = 2.0/(nx-1)
- dy = 2.0/(ny-1)
- sigma = .2
- dt = sigma*dx
-
- x = np.linspace(0,2,nx)
- y = np.linspace(0,2,ny)
-
- u = np.ones((ny,nx)) ##create a 1xn vector of 1's
- un = np.ones((ny,nx)) ##
-
- ###Assign initial conditions
- u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2
- ###Plot Initial Condition
- fig = plt.figure(figsize=(11,7), dpi=100)
- ax = fig.gca(projection='3d')
- X, Y = np.meshgrid(x,y)
- surf = ax.plot_surface(X,Y,u[:])
-
- u = np.ones((ny,nx))
- u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2
-
- for n in range(nt+1): ##loop across number of time steps
- un = u.copy()
- for i in range(1, len(u)):
- for j in range(1, len(u)):
- u[i,j] = un[i, j] - (c*dt/dx*(un[i,j] - un[i-1,j]))-(c*dt/dy*(un[i,j]-un[i,j-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')
- surf2 = ax.plot_surface(X,Y,u[:])
输出结果：
初始条件：
[attach]69968[/attach]
最终结果：
[attach]69969[/attach]
做CFD计算的童鞋可能更容易理解云图，这里以云图的方式显示：
[attach]69970[/attach]
[attach]69971[/attach]
只需修改图形显示代码：
- plt.figure()
- im=plt.imshow(u,interpolation='bilinear',origin='lower',cmap=cm.rainbow)
- plt.contour(X,Y,u)
- plt.colorbar(im,orientation='vertical')
【来源：http://nbviewer.ipython.org/gith ... ons/07_Step_5.ipynb】
|
