CFD Python:求解Navier-Stokes方程的12個步驟(1)

ieic5tq8ymk e1668923225194

(本文經Prof. Lorena A. Barba授權轉載翻譯,原文鏈結請參照底下參考資料)

CFD Python:求解Navier-Stokes方程的12個步驟(1)
CFD Python:求解Navier-Stokes方程的12個步驟(2)
CFD Python:求解Navier-Stokes方程的12個步驟(3)
CFD Python:求解Navier-Stokes方程的12個步驟(4)

哈囉!歡迎來到求解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)
[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  1.  1.  1.  1.  1.]

現在用matplotlibplot功能來看看初始條件。我們已經引入matplotlib的繪圖庫pyplot且繪圖函式為plot,因此可以呼叫pyplot.plot。如果想要學習更多matplotlib,可參考有著範例繪圖的網站。

這裡使用簡單的二維繪圖語法: plot(x, y),其中x為均分網格節點:

pyplot.plot(numpy.linspace(0, 2, nx), u);

NS1

為什麼帽函式(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);
NS2
好的!帽函式(hat function)肯定移動至右邊了,但它不再是個hat,發生了什麼事情呢?

參考資料

  1. 12 steps to Navier–Stokes: Step 1: 1-D Linear ConvectionLorena A. Barba

Leave a Comment

發佈留言必須填寫的電子郵件地址不會公開。