Python and CFD --2d Diffusion

06-03

來自專欄 Fiction Studio

原文鏈接：

Jupyter Notebook Viewer?

nbviewer.jupyter.org

首先來看看二維的擴散方程。

$frac{partial u}{partial t} = u frac{partial ^2 u}{partial x^2} + u frac{partial ^2 u}{partial y^2}\ ag{1}$

其中擴散項採用中心差分格式：

$frac{u_{i,j}^{n+1} - u_{i,j}^n}{Delta t} = u frac{u_{i+1,j}^n - 2 u_{i,j}^n + u_{i-1,j}^n}{Delta x^2} + u frac{u_{i,j+1}^n-2 u_{i,j}^n + u_{i,j-1}^n}{Delta y^2}\ ag{2}$

整理得到：

$egin{align} u_{i,j}^{n+1} = u_{i,j}^n &+ frac{ u Delta t}{Delta x^2}(u_{i+1,j}^n - 2 u_{i,j}^n + u_{i-1,j}^n) + frac{ u Delta t}{Delta y^2}(u_{i,j+1}^n-2 u_{i,j}^n + u_{i,j-1}^n) end{align}\ ag{3}$

下面使用python實現：

import numpy as npfrom matplotlib import cmimport matplotlib.pyplot as pltfrom mpl_toolkits.mplot3d import Axes3D

開始定義需要使用的變數：

nx = 31ny = 31nt = 17nu = .05dx = 2 / (nx - 1)dy = 2 / (ny - 1)sigma = .25dt = sigma * dx * dy / nux = np.linspace(0, 2, nx)y = np.linspace(0, 2, ny)u = np.ones((ny, nx)) # create a 1xn vector of 1sun =np.ones((ny, nx))#非常常見的手法了u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2

先看看初始場：

fig = plt.figure()ax = fig.gca(projection=3d)X, Y = np.meshgrid(x, y)surf = ax.plot_surface(X, Y, u, rstride=1, cstride=1, cmap=cm.viridis, linewidth_=0, antialiased=False)ax.set_xlim(0, 2)ax.set_ylim(0, 2)ax.set_zlim(1, 2.5)ax.set_xlabel($x$)ax.set_ylabel($y$);

定義一個更新的函數：

def diffuse(nt): u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2 #初場定義 for n in range(nt + 1): un = u.copy() u[1:-1, 1:-1] = (un[1:-1,1:-1] + nu * dt / dx**2 * (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) + nu * dt / dy**2 * (un[2:,1: -1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1])) #更新場 u[0, :] = 1 u[-1, :] = 1 u[:, 0] = 1 u[:, -1] = 1 #邊界條件 fig = plt.figure() ax = fig.gca(projection=3d) surf = ax.plot_surface(X, Y, u[:], rstride=1, cstride=1, cmap=cm.viridis, linewidth_=0, antialiased=True) ax.set_zlim(1, 2.5) ax.set_xlabel($x$) ax.set_ylabel($y$);#計算完成後繪製u的分布

最後看看不同時刻的場分布：

diffuse(10)

diffuse(14)

diffuse(30)

diffuse(50)

很明顯的，隨著時間的推移，原本間斷分布的場開始向四周擴散，變得光滑起來。

來看看維基百科的定義：

Diffusionis the net movement of molecules or atoms from a region of high concentration (or high chemical potential) to a region of low concentration (or low chemical potential) as a result of random motion of the molecules or atoms. Diffusion is driven by a gradient in chemical potential of the diffusing species.

擴散現象的之所以存在是因為有濃度差。這是一種十分常見的物理現象，值得一提的是現在很多社會現象也引入類似於擴散的概念來建立微分方程（數學建模競賽了解一下）。

數值耗散

在對cfd方程離散後，離散方程與微分方程之間的誤差是由 $sum_{}^{}O(h^{n})$ 來表示，h為網格尺寸。如果n中有2，那麼認為離散的方程由耗散效應，如果n中有3，就認為離散方程由色散效應。

耗散效應就是把尖銳的解抹平，變得光滑，很多時候一些方程還會人為的引入二階導數項，雖然這樣增加了誤差，但是可以提高解的光滑性，避免奇點。