(本文經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 inline
def 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 國際 授權條款授權.