(本文經Prof. Lorena A. Barba授權轉載翻譯自Prof. Lorena A. Barba的教材”12 Steps to Navier-Stokes”，原文鏈結請參照底下參考資料)
Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved BSD-3 license. (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
CFD Python：求解Navier-Stokes方程的12個步驟(1) CFD Python：求解Navier-Stokes方程的12個步驟(2) CFD Python：求解Navier-Stokes方程的12個步驟(3) CFD Python：求解Navier-Stokes方程的12個步驟(4) CFD Python：求解Navier-Stokes方程的12個步驟(5) CFD Python：求解Navier-Stokes方程的12個步驟(6) CFD Python：求解Navier-Stokes方程的12個步驟(7) CFD Python：求解Navier-Stokes方程的12個步驟(8) CFD Python：求解Navier-Stokes方程的12個步驟(9) CFD Python：求解Navier-Stokes方程的12個步驟(10) CFD Python：求解Navier-Stokes方程的12個步驟(11) CFD Python：求解Navier-Stokes方程的12個步驟(12)
上一個步驟解決了二維Burgers’ 方程：這是流體力學研究中重要的方程式，因為它包含完整的對流非線性。透過這個練習，我們也累積逐步撰寫Navier-Stokes求解器的經驗。
步驟9：2D Laplace方程
底下是二維Laplace方程：
\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = 0
我們知道如何離散二階導數。但是想想看，Laplace方程具有擴散現象的典型特徵。因此，它必須使用中心差分來離散，使離散與我們要模擬的物理現象一致。離散方程式為：
\frac{p_{i+1, j}^n - 2p_{i,j}^n + p_{i-1,j}^n}{\Delta x^2} + \frac{p_{i,j+1}^n - 2p_{i,j}^n + p_{i, j-1}^n}{\Delta y^2} = 0
注意Laplace方程與時間無關 – 即沒有p^{n+1}。與波動隨著時間的變化（如前幾個步驟）無關，Laplace方程計算的是給定邊界條件下的系統平衡狀態。
如果您修過熱傳導課程，您將認出Laplace方程就是穩態熱傳導方程。
我們不是計算系統在某個時間t的位置，而是迭代求解p_{i,j}^n，直到它滿足我們指定的條件。只有當迭代次數趨近於無窮大時，系統才能達到平衡，但我們可以透過迭代求解，直到相鄰兩次迭代之間的數值變化很小為止，此時近似平衡狀態。
重新整理離散化方程，求解p_{i,j}^n：
p_{i,j}^n = \frac{\Delta y^2(p_{i+1,j}^n+p_{i-1,j}^n)+\Delta x^2(p_{i,j+1}^n + p_{i,j-1}^n)}{2(\Delta x^2 + \Delta y^2)}
在兩個方向使用二階中心差分是求解Laplace算子最常用的方法。它也被稱為五點差分算子。
假設任一處的p =0來求解Laplace方程。然後，添加如下邊界條件：
在x=0處，p=0
在x=2處，p=y
在y=0,1處，\frac{\partial p}{\partial y}=0
在這些條件下，Laplace方程有一個解析解：
p(x,y)=\frac{x}{4}-4\sum_{n=1,odd}^{\infty}\frac{1}{(n\pi)^2\sinh2n\pi}\sinh n\pi x\cos n\pi y
練習
自己撰寫程式碼，使用迴圈來求解Poisson方程，並採用第一節課的寫法。然後，考慮如何使用函式編寫示例，並修改您的程式碼以套用新寫法。您能想到放棄舊寫法並採用模組化程式碼的原因嗎？
其他提示：
- 可視化每個迭代過程的步驟
- 思考邊界條件在做什麼
- 思考偏微分方程在做什麼
使用函數
記得使用Python編寫函式的課程嗎？我們將在此練習使用該作法。分別定義兩個函式：一個在3D投影圖中繪製數據，另一個用於迭代求解p，直到p的L1 Norm變化小於指定值。
import numpy
from matplotlib import pyplot, cm
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inlinedef plot2D(x, y, p):
fig = pyplot.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')
X, Y = numpy.meshgrid(x, y)
surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,
linewidth=0, antialiased=False)
ax.set_xlim(0, 2)
ax.set_ylim(0, 1)
ax.view_init(30, 225)
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')plot2D函數接受三個參數，分別是x向量，y向量和p矩陣。給定這三個值，它將產生一個3D投影圖，我們可以設置繪圖限制來提供一個漂亮的視角。
p_{i,j}^n = \frac{\Delta y^2(p_{i+1,j}^n+p_{i-1,j}^n)+\Delta x^2(p_{i,j+1}^n + p_{i,j-1}^n)}{2(\Delta x^2 + \Delta y^2)}
laplace2d函式接受五個參數，p矩陣，y向量，dx，dy和值l1norm_target。最後一個值定義兩個連續迭代之間的p矩陣應該多接近，以便迴圈中斷並返回計算結果的p值。
注意，您的程式碼執行上述單元將沒有輸出。這是因為雖然定義了函式，但尚未調用該函式。現在，它可以像numpy.linspace或命名空間中的其他函式一樣供您使用。
##variable declarations
nx = 31
ny = 31
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
##initial conditions
p = numpy.zeros((ny, nx)) # create a XxY vector of 0's
##plotting aids
x = numpy.linspace(0, 2, nx)
y = numpy.linspace(0, 1, ny)
##boundary conditions
p[:, 0] = 0 # p = 0 @ x = 0
p[:, -1] = y # p = y @ x = 2
p[0, :] = p[1, :] # dp/dy = 0 @ y = 0
p[-1, :] = p[-2, :] # dp/dy = 0 @ y = 1使用plot2D函式查看初始條件。如果函式已正確定義，則可以鍵入plot2D並按Tab鍵來自動完成選項。
成功了！這是問題的初始狀態，其中p的值除了沿著x=2的地方，p=y之外，其他地方都為零。現在，使用指定的L1目標值0.01運行laplace2d函式。
[提示：如果無法記住傳送到函式的變數順序，可以輸入 laplace2d()，iPython將彈出一個小提示框提醒您。]
p = laplace2d(p, y, dx, dy, 1e-4)現在，使用我們的繪圖函式來繪製新的p值。
plot2D(x, y, p)
參考資料
- 12 steps to Navier–Stokes: Step 9: 2D Laplace Equation，Lorena A. Barba
本著作係採用創用 CC 姓名標示 4.0 國際 授權條款授權.