(本文經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)
哈囉!歡迎來到求解Naiver-Stokes方程的教學,這是波士頓大學教授Lorena A. Barba從2009年春季於計算流體力學課程上教學的程式模組。這項課程假定您已具有基本程式能力(任何程式語言)以及偏微分方程與流體力學相關基礎知識。
此程式模組出自 Rio Yokota博士的主意(Barba研究室的後博士研究員),並由Barba教授與其學生在幾個學期的課程教學中不斷完善,此課程完全使用Pyhton教學,不懂Pyhton的學生也能透過此教程來學習。
這系列教學將會從零開始帶您從第一個步驟寫出自己的Navier-Stokes方程求解器,而我們即將開始。如果一開始什麼都不懂也不用擔心,我們將在接下來的過程中為您講解,而您也能透過Prof. Barba’s lectures on YouTube的影片來學習。
步驟1:一維線性對流方程
一維線性對流方程式是能用來學習CFD的基本模型。一個簡單的方程式能學到這麼多真令人驚喜,該方程如下:
\frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0
加上給定的初始條件(可理解為波),該方程式代表初始條件給定的波以速度c傳遞,且波形不變。假設初始條件為u(x,0)=u_0(x),則方程式精確解為u(x,t)=u_0(x-ct)。
我們使用前向差分來離散此方程的時間項,並用後向差分來離散空間項。離散的空間座標x用索引i = 0 至 N來表示,且離散後的時間步階以尺寸\Delta t表示
從微分定義來看可以得知:
\frac{\partial u}{\partial x}\approx \frac{u(x+\Delta x)-u(x)}{\Delta x}
因此,離散方程可寫成如下形式:
\frac{u_i^{n+1}-u_i^n}{\Delta t} + c \frac{u_i^n - u_{i-1}^n}{\Delta x} = 0
其中n與n+1代表時間上連續的兩個步階,而i – 1與i代表X方向上離散座標的兩個相鄰點。如果有給定的初始條件,那麼離散方程式的未知數就只有u_i^{n+1}。因此,能夠以一個隨著時間推進的方程式來求解未知數,如下:
u_i^{n+1} = u_i^n - c \frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n)
現在嘗試用Python來實踐它。
首先必須引入一些新的函式庫來幫助我們:
- numpy是一個提供許多矩陣操作的函式庫,類似MATLAB
- matplotlib是一個2D繪圖函式庫用來幫我們畫出結果
- time and sys提供基本的時間函式來幫助我們放慢動畫以便觀察結果
# Remember: comments in python are denoted by the pound sign
import numpy #here we load numpy
from matplotlib import pyplot #here we load matplotlib
import time, sys
#and load some utilities
#this makes matplotlib plots appear in the notebook (instead of a separate window)
%matplotlib inline
現在來定義一些新變數,我們想要在長度2的空間中定義等間距的網格點,因此x_i\in(0,2)。變數 nx表示網格點數量,且變數dx為任意一組兩個相鄰網格點的間距。
nx = 41 # try changing this number from 41 to 81 and Run All ... what happens?
dx = 2 / (nx-1)
nt = 25 #nt is the number of timesteps we want to calculate
dt = .025 #dt is the amount of time each timestep covers (delta t)
c = 1 #assume wavespeed of c = 1
此外,也需要設置初始條件。初始速度u_0在間距0.5 \leq x \leq 1給定u = 2,並在(0, 2)的任意某處皆為 u = 1(此為hat function)。
這邊使用函式ones()
來定義一個numpy的陣列,該陣列有著nx元素且各項元素值為1。
u = numpy.ones(nx) #numpy function ones()
u[int(.5 / dx):int(1 / dx + 1)] = 2 #setting u = 2 between 0.5 and 1 as per our I.C.s
print(u)
現在用matplotlib
的plot
功能來看看初始條件。我們已經引入matplotlib
的繪圖庫pyplot
且繪圖函式為plot
,因此可以呼叫pyplot.plot
。如果想要學習更多matplotlib
,可參考有著範例繪圖的網站。
這裡使用簡單的二維繪圖語法: plot(x, y)
,其中x為均分網格節點:
pyplot.plot(numpy.linspace(0, 2, nx), u);
為什麼帽函式(hat function)沒出現完美的直線?可以稍微思考看看。
現在是時候用有限差分來離散對流方程,對於陣列每個元素而言,我們需要執行下列操作
u_i^{n+1} = u_i^n - c \frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n)
把計算結果儲存在新(暫存)的陣列un,而un會做為下個時間步的u來計算。我們會不斷重複此操作,直到有足夠的時間步來讓我們觀察波動如何傳遞。
我們先用Numpy函式ones()
來初始化陣列un,陣列un用來保存第n+1個時間步的數值。
這邊可能有兩個迭代操作:一個是空間上的迭代,另一個是時間上的迭代,因此我們將會以巢狀迴圈的寫法來實踐。注意range()
函式的使用,當我們寫for i in range(1,nx)
,雖然會迭代整個陣列u,卻會略過第一個元素(編號0的元素),這是為什麼呢?
un = numpy.ones(nx) #initialize a temporary array
for n in range(nt): #loop for values of n from 0 to nt, so it will run nt times
un = u.copy() ##copy the existing values of u into un
for i in range(1, nx): ## you can try commenting this line and...
#for i in range(nx): ## ... uncommenting this line and see what happens!
u[i] = un[i] - c * dt / dx * (un[i] - un[i-1])
我們稍後會知道上述程式碼效率很差,在Python還有著更好的寫法,不過先繼續學習。
現在來試著畫出隨著時間推進後的陣列u。
pyplot.plot(numpy.linspace(0, 2, nx), u);
參考資料
- 12 steps to Navier–Stokes: Step 1: 1-D Linear Convection,Lorena A. Barba
本著作係採用創用 CC 姓名標示 4.0 國際 授權條款授權.