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

ieic5tq8ymk e1668923225194

(本文經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方程,其中\vec{v}表示速度場:

\nabla \cdot\vec{v} = 0
\frac{\partial \vec{v}}{\partial t}+(\vec{v}\cdot\nabla)\vec{v} = -\frac{1}{\rho}\nabla p + \nu \nabla^2\vec{v}

第一個方程表示恆定密度下的質量守恆。第二個方程是動量守恆。但是出現了一個問題:不可壓縮流的連續方程沒有主要變數,也沒有明顯的方法來耦合速度和壓力。相比之下,可壓縮流的質量連續性為密度ρ提供一個演化方程,它耦合了關於ρ和p的狀態方程。

在不可壓縮流中,連續方程\nabla \cdot\vec{v} = 0提供一個動力上的約束,要求壓力場隨著膨脹率\nabla \cdot\vec{v} = 0消失而演變。為了耦合速度與壓力,解決方法是構建一個壓力場來滿足連續性,對動量方程求取散度,進而得到壓力的泊松方程(Poisson equation)!

步驟10:2D泊松方程(Poisson equation)

泊松方程是透過在Laplace方程右側加上一個源項而得到:

\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = b

因此,與Laplace方程不同,場內有一個有限值影響解。泊松方程起到”鬆弛”場內初始源的作用。而泊松方程的離散形式看起來幾乎與步驟9相同,只是多了源項:

\frac{p_{i+1,j}^{n}-2p_{i,j}^{n}+p_{i-1,j}^{n}}{\Delta x^2}+\frac{p_{i,j+1}^{n}-2 p_{i,j}^{n}+p_{i,j-1}^{n}}{\Delta y^2}=b_{i,j}^{n}

與先前一樣,重新排列這個方程,以便在點i,j處獲得p。因此得到:

p_{i,j}^{n}=\frac{(p_{i+1,j}^{n}+p_{i-1,j}^{n})\Delta y^2+(p_{i,j+1}^{n}+p_{i,j-1}^{n})\Delta x^2-b_{i,j}^{n}\Delta x^2\Delta y^2}{2(\Delta x^2+\Delta y^2)}

假設初始狀態所有地方的p=0,並使用邊界條件如下:

x=0,2和y=0,1的p = 0

源項包括域內的兩個初始尖峰,如下所示:

b_{i,j}=100i=\frac{1}{4}nx, j=\frac{1}{4}ny

b_{i,j}=-100i=\frac{3}{4}nx, j=\frac{3}{4}ny

b_{i,j}=0 在其他地方

迭代將在偽時間推進下來鬆弛初始尖峰,而泊松方程的鬆弛變得越來越慢。為什麼?

看看泊松方程可能的程式碼寫法。像往常一樣,載入我們最愛的Python函式庫。我們還想做一些可愛的3D圖,於是定義參數並完成初始化。你對下面的方法有什麼看法?

import numpy
from matplotlib import pyplot, cm
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
# Parameters
nx = 50
ny = 50
nt  = 100
xmin = 0
xmax = 2
ymin = 0
ymax = 1

dx = (xmax - xmin) / (nx - 1)
dy = (ymax - ymin) / (ny - 1)

# Initialization
p  = numpy.zeros((ny, nx))
pd = numpy.zeros((ny, nx))
b  = numpy.zeros((ny, nx))
x  = numpy.linspace(xmin, xmax, nx)
y  = numpy.linspace(xmin, xmax, ny)

# Source
b[int(ny / 4), int(nx / 4)]  = 100
b[int(3 * ny / 4), int(3 * nx / 4)] = -100

有了這個,我們打算透過偽時間來推進初始的猜測值。下方程式碼與步驟9用於求解Laplace方程的函式有何不同?

for it in range(nt):

    pd = p.copy()

    p[1:-1,1:-1] = (((pd[1:-1, 2:] + pd[1:-1, :-2]) * dy**2 +
                    (pd[2:, 1:-1] + pd[:-2, 1:-1]) * dx**2 -
                    b[1:-1, 1:-1] * dx**2 * dy**2) / 
                    (2 * (dx**2 + dy**2)))

    p[0, :] = 0
    p[ny-1, :] = 0
    p[:, 0] = 0
    p[:, nx-1] = 0

也許可以重用步驟9中的繪圖函式,你認為呢?

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.view_init(30, 225)
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
plot2D(x, y, p)

NS10 1

這是程式碼重複使用的驚奇!

現在,你可能會想:好吧,如果寫了這麼一個小巧的函式,它做了一些非常有用的事情,我想一遍又一遍使用它。我如何在不每次複製貼上的情況做到這一點?

如果你對此非常好奇,那麼你需要學習有關封裝(packaging)的知識。但這已經超出了我們CFD課程範疇。如果你真的想知道,就得自己去查找了。

參考資料

  1. 12 steps to Navier–Stokes: Step 10: 2D Poisson EquationLorena A. Barba

創用 CC 授權條款
本著作係採用創用 CC 姓名標示 4.0 國際 授權條款授權.

Leave a Comment

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